python实现运筹学排序与统筹安排问题 python运筹学的包

1. cplex介绍

Cplex是IBM出的一款科学计算软件，从IBM官网可以下到最新版。最近惊奇的发现python竟然能直接安装使用cplex了！
安装方法如下：

pip install cplex
pip install docplex

这个版本是限制变量数量的，如果要使用无限制版，请购买正版cplex或申请学术版，然后将python/cplex文件夹整体复制到anaconda的site-package下。

2. 示例代码和说明

用一个简单案例尝试一下：
$\max x_1+2x_2+3x_3+x_4$
s.t. $-x_1+x_2+x_3+10x_4\le20$
$x_1-3x_2+x_3\le30$
$x_2-3.5x_4=0$
$x_1\le40$
$2\le x_4\le3$，$x_4\in Z$
$x_1,x_2,x_3\ge0$

import cplex
from cplex.exceptions import CplexError
my_obj = [1.0, 2.0, 3.0, 1.0]
my_ub = [40.0, cplex.infinity, cplex.infinity, 3.0]
my_lb = [0.0, 0.0, 0.0, 2.0]
my_ctype = "CCCI"
my_colnames = ["x1", "x2", "x3", "x4"]
my_rhs = [20.0, 30.0, 0.0]
my_rownames = ["r1", "r2", "r3"]
my_sense = "LLE"


def populatebyrow(prob):
    prob.objective.set_sense(prob.objective.sense.maximize)
    prob.variables.add(obj=my_obj, lb=my_lb, ub=my_ub, types=my_ctype,
                       names=my_colnames)
    rows = [[my_colnames, [-1.0, 1.0, 1.0, 10.0]],
            [my_colnames, [1.0, -3.0, 1.0, 0.0]],
            [my_colnames, [0.0, 1.0, -3.5, 0.0]]]
    prob.linear_constraints.add(lin_expr=rows, senses=my_sense,rhs=my_rhs, names=my_rownames)

try:
    my_prob = cplex.Cplex()
    handle = populatebyrow(my_prob)
    my_prob.solve()
    
except CplexError as exc:
    print(exc)

print("Solution status = ", my_prob.solution.status[my_prob.solution.get_status()])
print("Solution value  = ", my_prob.solution.get_objective_value())
x = my_prob.solution.get_values()
print('x: ',x)

输出为：

代码中：
my_ctype中，C表示连续变量，I表示整数型变量
my_sense，G表示大于等于，E表示等于，L表示小于等于
不过由于只用能矩阵和向量形式表示，因此对于某些复杂的模型，写起来还是比较麻烦的。

3. 使用建模语言

docplex（Decision Optimization CPLEX® Modeling for Python）是一个基于python的建模语言库，这里是官方文档。
关键语句包括：

• 定义模型，比如：m = Model(name=‘telephone_production’)
• 定义变量，比如：desk = m.continuous_var(name=‘desk’)；整数变量和01变量分别为integer_var和binary_var。
• 定义约束，比如：ct_assembly = m.add_constraint( 0.2 * desk + 0.4 * cell <= 400)
• 定义目标函数，比如：m.maximize(12 * desk + 20 * cell)
• 求解问题：m.solve()

下面是一个简单的例子：

# first import the Model class from docplex.mp
from docplex.mp.model import Model

# create one model instance, with a name
m = Model(name='telephone_production')

# by default, all variables in Docplex have a lower bound of 0 and infinite upper bound
desk = m.continuous_var(name='desk')
cell = m.continuous_var(name='cell')

# constraint #1: desk production is greater than 100
m.add_constraint(desk >= 100)

# constraint #2: cell production is greater than 100
m.add_constraint(cell >= 100)

# constraint #3: assembly time limit
ct_assembly = m.add_constraint( 0.2 * desk + 0.4 * cell <= 400)

# constraint #4: paiting time limit
ct_painting = m.add_constraint( 0.5 * desk + 0.4 * cell <= 490)

m.maximize(12 * desk + 20 * cell)
s = m.solve()
m.print_information()
m.print_solution()

进阶技巧包括：

• 使用字典、元祖等进行定义，比如：
  costs = {(1,3): 2, (1,5):4, (2,4):5, (2,5):3}source = range(1, 3)tm.minimize(tm.sum(x[i,j]*costs.get((i,j), 0)))x = {(i,j): tm.continuous_var(name='x_{0}_{1}'.format(i,j)) for i in source for j in target}
• 使用列表推导式，比如：
  x = {(i,j): tm.continuous_var(name='x_{0}_{1}'.format(i,j)) for i in source for j in target}tm.minimize(tm.sum(x[i,j]*costs.get((i,j), 0) for i in source for j in target))
• 使用KPI，在report中查看。

下面是一个稍复杂的例子：

from collections import namedtuple

from docplex.mp.model import Model
from docplex.util.environment import get_environment


# ----------------------------------------------------------------------------
# Initialize the problem data
# ----------------------------------------------------------------------------

FOODS = [
    ("Roasted Chicken", 0.84, 0, 10),
    ("Spaghetti W/ Sauce", 0.78, 0, 10),
    ("Tomato,Red,Ripe,Raw", 0.27, 0, 10),
    ("Apple,Raw,W/Skin", .24, 0, 10),
    ("Grapes", 0.32, 0, 10),
    ("Chocolate Chip Cookies", 0.03, 0, 10),
    ("Lowfat Milk", 0.23, 0, 10),
    ("Raisin Brn", 0.34, 0, 10),
    ("Hotdog", 0.31, 0, 10)
]

NUTRIENTS = [
    ("Calories", 2000, 2500),
    ("Calcium", 800, 1600),
    ("Iron", 10, 30),
    ("Vit_A", 5000, 50000),
    ("Dietary_Fiber", 25, 100),
    ("Carbohydrates", 0, 300),
    ("Protein", 50, 100)
]

FOOD_NUTRIENTS = [
    ("Roasted Chicken", 277.4, 21.9, 1.8, 77.4, 0, 0, 42.2),
    ("Spaghetti W/ Sauce", 358.2, 80.2, 2.3, 3055.2, 11.6, 58.3, 8.2),
    ("Tomato,Red,Ripe,Raw", 25.8, 6.2, 0.6, 766.3, 1.4, 5.7, 1),
    ("Apple,Raw,W/Skin", 81.4, 9.7, 0.2, 73.1, 3.7, 21, 0.3),
    ("Grapes", 15.1, 3.4, 0.1, 24, 0.2, 4.1, 0.2),
    ("Chocolate Chip Cookies", 78.1, 6.2, 0.4, 101.8, 0, 9.3, 0.9),
    ("Lowfat Milk", 121.2, 296.7, 0.1, 500.2, 0, 11.7, 8.1),
    ("Raisin Brn", 115.1, 12.9, 16.8, 1250.2, 4, 27.9, 4),
    ("Hotdog", 242.1, 23.5, 2.3, 0, 0, 18, 10.4)
]

Food = namedtuple("Food", ["name", "unit_cost", "qmin", "qmax"])
Nutrient = namedtuple("Nutrient", ["name", "qmin", "qmax"])


# ----------------------------------------------------------------------------
# Build the model
# ----------------------------------------------------------------------------

def build_diet_model(**kwargs):
    # Create tuples with named fields for foods and nutrients

    foods = [Food(*f) for f in FOODS]
    nutrients = [Nutrient(*row) for row in NUTRIENTS]

    food_nutrients = {(fn[0], nutrients[n].name):
                      fn[1 + n] for fn in FOOD_NUTRIENTS for n in range(len(NUTRIENTS))}

    # Model
    mdl = Model(name='diet', **kwargs)

    # Decision variables, limited to be >= Food.qmin and <= Food.qmax
    qty = mdl.continuous_var_dict(foods, lb=lambda f: f.qmin, ub=lambda f: f.qmax, name=lambda f: "q_%s" % f.name)

    # Limit range of nutrients, and mark them as KPIs
    for n in nutrients:
        amount = mdl.sum(qty[f] * food_nutrients[f.name, n.name] for f in foods)
        mdl.add_range(n.qmin, amount, n.qmax)
        mdl.add_kpi(amount, publish_name="Total %s" % n.name)

    # Minimize cost
    mdl.minimize(mdl.sum(qty[f] * f.unit_cost for f in foods))

    return mdl

# ----------------------------------------------------------------------------
# Solve the model and display the result
# ----------------------------------------------------------------------------


if __name__ == '__main__':
    mdl = build_diet_model()
    mdl.print_information()
    mdl.export_as_lp()
    if mdl.solve():
        mdl.float_precision = 3
        print("* model solved as function:")
        mdl.print_solution()
        mdl.report_kpis()
        # Save the CPLEX solution as "solution.json" program output
        with get_environment().get_output_stream("solution.json") as fp:
            mdl.solution.export(fp, "json")
    else:
        print("* model has no solution")

输出为：

Model: diet
 - number of variables: 9
   - binary=0, integer=0, continuous=9
 - number of constraints: 7
   - linear=7
 - parameters: defaults
* model solved as function:
objective: 2.690
  "q_Spaghetti W/ Sauce"=2.155
  "q_Chocolate Chip Cookies"=10.000
  "q_Lowfat Milk"=1.831
  "q_Hotdog"=0.930
*  KPI: Total Calories      = 2000.000
*  KPI: Total Calcium       = 800.000
*  KPI: Total Iron          = 11.278
*  KPI: Total Vit_A         = 8518.433
*  KPI: Total Dietary_Fiber = 25.000
*  KPI: Total Carbohydrates = 256.806
*  KPI: Total Protein       = 51.174

4. MIP log

打开log，在 Model.solve()中设置log_output=True。MIP的监控参数可以设置为如下：

我们看一个官方例子：

Tried aggregator 1 time. 
Presolve time =    0.00 sec. (0.00 ticks)
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: none, using 1 thread.
Root relaxation solution time = 0.00 sec (0.00 ticks)

         Nodes                                 Cuts/   
   Node  Left     Objective  IInf  Best Integer     Best Bound    ItCnt    Gap
*     0+    0                            0.0000     3261.8212        8     ---
*     0+    0                         3148.0000     3261.8212        8    3.62%
      0     0     3254.5370     7     3148.0000       Cuts: 5       14    3.38%
      0     0     3246.0185     7     3148.0000       Cuts: 3       24    3.11%
*     0+    0                         3158.0000     3246.0185       24    2.79%
      0     0     3245.3465     9     3158.0000       Cuts: 5       27    2.77%
      0     0     3243.4477     9     3158.0000       Cuts: 5       32    2.71%
      0     0     3242.9809    10     3158.0000     Covers: 3       36    2.69%
      0     0     3242.8397    11     3158.0000     Covers: 1       37    2.69%
      0     0     3242.7428    11     3158.0000       Cuts: 3       39    2.68%
      0     2     3242.7428    11     3158.0000     3242.7428       39    2.68%
     10    11     3199.1875     2     3158.0000     3215.1261       73    1.81%
*    20+   11                         3168.0000     3215.1261       89    1.49%
     20    13     3179.0028     5     3168.0000     3215.1261       89    1.49%
     30    15     3179.9091     3     3168.0000     3197.5227      113    0.93%
*    39     3      integral     0     3186.0000     3197.3990      126    0.36%
     40     3     3193.7500     1     3186.0000     3197.3990      128    0.36% 

Cover cuts applied:  9
Zero-half cuts applied:  2
Gomory fractional cuts applied:  1

Solution pool: 5 solutions saved. 
MIP-Integer optimal solution:  Objective =   3.1860000000e+03
Solution time =    0.01 sec.  (0.00 ticks)  Iterations = 131   Nodes = 44

log中，Node数表示生成了44个节点的树，共执行了131次（对偶单纯形法）迭代
第一列*表示找到了一个新的整数解；第二列表示当前节点，+表示用启发式方法找到了一个解，对应第四列为空；第三列表示剩余需要探索的节点数；第四列表示当前节点的目标函数值，或者显示停止分支的原因；第五列表示不满足整数约束的变量个数；第六列是当前找到的最优整数解；第七列是松弛问题的最优解，和第六列构成问题的上下界；第八列是总的迭代次数；第九列是上下界之间的差距。