线性规划

解决线性规划问题Lingo用的比较多，但是Lingo存在很多缺点。

  1. Lingo是收费的，这样的小体量软件，大部分学校也不会专门采购。对于学生来说，只能用盗版，安装繁琐。
  2. Lingo只支持windows系统，对于使用linux系统的用户来说，不太方便。
  3. Lingo的UI非常难看，比Matlab还要远古。
  4. Lingo的语法错误提示非常诡异，有时候报错信息不太明确，甚至还有非常神奇的特性，不容易找到错误。上周听讲座时，老师写Lingo代码报错了，但是老师也不知道哪里错了。于是把代码剪下来，再重启Lingo，粘贴进去又可以运行了。在现场这让我大开眼界，而这种情况在Python和Matlab里是不会出现的，物理学在Lingo里是不存在的。
  5. Lingo的求解速度不太理想，对于大规模问题，求解速度很慢,而且不会显示进度条，有时候无法掌握求解状态。
  6. Lingo虽然简单，但是也要花时间学习，对于一些简单的线性规划问题，使用Lingo有点大材小用。

使用熟悉的工具解决这类问题显然更好，而Python的轮子是最多的，有许多的开源求解器可以使用，比如scipy.optimize.linprog，但是这个求解器只支持简单的线性规划问题，不支持整数规划问题。有一个更强大的工具：glpk。glpk开源，可以用来求解线性规划问题。glpk支持多种语言，这里使用python调用glpk求解线性规划问题。使用pulp库可以更方便的调用glpk求解线性规划问题。
pulp作为单词的意思是纸浆，但是这个库的名字是Python Ultimate LP

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import pulp as pl
# Define the model
lp_problem = pl.LpProblem("Max_z", pl.LpMaximize)
#定义决策变量
x1 = pl.LpVariable('x1', lowBound=0, cat='Continuous')
x2 = pl.LpVariable('x2', lowBound=0, cat='Continuous')
#定义目标函数
lp_problem += 3*x1 + 2*x2, "z"
#定义约束条件
lp_problem += x1 + x2 <= 4
lp_problem += 2*x1 + x2 <= 7
lp_problem += x1 <= 3
#求解
print(lp_problem)  # 打印问题
print(pl.LpStatus[lp_problem.status]) # 打印问题状态
lp_problem.solve()  # 求解问题
print("最优解为：", pl.value(lp_problem.objective))

Max_z:
MAXIMIZE
3*x1 + 2*x2 + 0
SUBJECT TO
_C1: x1 + x2 <= 4

_C2: 2 x1 + x2 <= 7

_C3: x1 <= 3

VARIABLES
x1 Continuous
x2 Continuous

Not Solved
最优解为： 11.0

但是后面我采用了Gurobi，是付费软件，但是在校学生有免费的许可证，非常容易申请，下面是一些题。

min(2x+3y+5z)s.t.x+2y+2z≥303x+y+2z≥202x+y+10z≥40z={−3,if x≥16,if x<1x,y≥0min(2x+3y+5z)\\ s.t.\\ \begin{align*} x + 2y + 2z &\geq 30 \\ 3x + y + 2z &\geq 20 \\ 2x + y + 10z &\geq 40 \\ z &= \begin{cases} -3, & \text{if } x \geq 1 \\ 6, & \text{if } x < 1 \\ \end{cases} \\ x, y &\geq 0 \end{align*} 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import gurobipy 
from gurobipy import GRB

        
model = gurobipy.Model('q1')
# add variables
x = model.addVar(vtype=GRB.CONTINUOUS, name='x')
y = model.addVar(vtype=GRB.CONTINUOUS, name='y')
z = model.addVar(vtype=GRB.CONTINUOUS, name='z')
model.update() # if not update, the variables will not be added to the model
# set objective
model.setObjective(2*x + 3*y + 5*z, sense=GRB.MINIMIZE)
# add constraints
model.addConstr(x + 2*y + 2*z >= 30)
model.addConstr(3*x + y +2*z >= 20)
model.addConstr(2*x + y + 10*z <=50)
model.addConstr(2*x + y + 10*z >= 40)
model.addConstrs((i>=0 for i in [x,y,z]), 'nonneg') # non-negative constraints，但是没有必要，因为默认就是非负的
# solve
model.optimize()
# print solution

print(model.getVars()) # get all variables,return a list of variables

# use a function to display the output
def display_out(model):
    if model.status == gurobipy.GRB.OPTIMAL:
        for var in model.getVars():
            print(f'{var.varName} = {var.x}')
        # 输出目标值
        print(f'Objective value: {model.ObjVal}')
    else:
        print("No optimal solution found.")

display_out(model)
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows 11+.0 (26100.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 7 rows, 3 columns and 15 nonzeros
Model fingerprint: 0xb2aec4d2
Coefficient statistics:
  Matrix range     [1e+00, 1e+01]
  Objective range  [2e+00, 5e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [2e+01, 5e+01]
Presolve removed 4 rows and 0 columns
Presolve time: 0.01s
Presolved: 3 rows, 4 columns, 10 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   8.750000e+00   0.000000e+00      0s
       3    5.0714286e+01   0.000000e+00   0.000000e+00      0s

Solved in 3 iterations and 0.01 seconds (0.00 work units)
Optimal objective  5.071428571e+01
[<gurobi.Var x (value 0.9523809523809523)>, <gurobi.Var y (value 11.904761904761903)>, <gurobi.Var z (value 2.619047619047619)>]
x = 0.9523809523809523
y = 11.904761904761903
z = 2.619047619047619
Objective value: 50.71428571428571

min(2x+3y+5z)s.t.x+2y+2z≥303x+y+2z≥202x+y+10z≥40z={−3,if x≥16,if x<1x,y≥0min(2x+3y+5z)\\ s.t.\\ \begin{align*} x + 2y + 2z & \geq 30 \\ 3x + y + 2z & \geq 20 \\ 2x + y + 10z & \geq 40 \\ z & = \begin{cases} -3, & \text{if } x \geq 1 \\ 6, & \text{if } x < 1 \\ \end{cases} \\ x, y & \geq 0 \end{align*} 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import gurobipy as gp
from gurobipy import GRB

model2 = gurobipy.Model('q2')
# add variables
x = model2.addVar(vtype=GRB.CONTINUOUS,lb=0, name='x')
y = model2.addVar(vtype=GRB.CONTINUOUS,lb=0,name='y')
z = model2.addVar(vtype=GRB.CONTINUOUS, name='z')
b1 = model2.addVar(vtype=GRB.BINARY, name='b1') #add binary variable
model2.update() #用了addConstraint之后，需要update一下，否则会报错
#target function
model2.setObjective(2*x + 3*y + 5*z, sense=GRB.MINIMIZE)
# add constraints
model2.addConstr(x +2*y + 2*z >= 30)
model2.addConstr(3*x + y + 2*z >= 20)
model2.addConstr(2*x + y + 10*z >= 40)

'''
可以省去这种写法，直接在addVar的时候就加上约束条件  
model2.addConstrs((i>=0 for i in [x,y]))
'''

# 这里采用gurobipy的方法来添加条件约束,后面需要用到大M法，以后再说
model2.addGenConstrIndicator(b1,True,x>=1,name='x_cond')
model2.addGenConstrIndicator(b1,True,z == -3,name='x_cond1')
model2.addGenConstrIndicator(b1,False,z == 6,name='x_cond2')
#optimize
model2.optimize()
display_out(model2)
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows 11+.0 (26100.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 3 rows, 4 columns and 9 nonzeros
Model fingerprint: 0x74cceb2c
Model has 3 general constraints
Variable types: 3 continuous, 1 integer (1 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+01]
  Objective range  [2e+00, 5e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+01, 4e+01]
  GenCon rhs range [1e+00, 6e+00]
  GenCon coe range [1e+00, 1e+00]
Presolve removed 3 rows and 4 columns
Presolve time: 0.01s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.03 seconds (0.00 work units)
Thread count was 1 (of 8 available processors)

Solution count 1: 57 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.700000000000e+01, best bound 5.700000000000e+01, gap 0.0000%
x = -0.0
y = 9.0
z = 6.0
b1 = 0.0
Objective value: 57.0

max∑i=110∑j=15aij(xi+yj)(aij)=(12⋯567⋯10⋮⋮⋱⋮4647⋯50)约束条件∑i=110xi≤100∑j=15yj≤200∑i=15(xi+yj)≥150xi,yj≥0max \sum_{i=1}^{10}\sum_{j=1}^{5} a_{ij}(x_i+y_j) \\ (a_{ij}) = \begin{pmatrix} 1 & 2 & \cdots & 5 \\ 6 & 7 & \cdots & 10 \\ \vdots & \vdots & \ddots & \vdots \\ 46 & 47 & \cdots & 50 \\ \end{pmatrix} \\ 约束条件 \begin{align*} \sum_{i=1}^{10} x_i & \leq 100 \\ \sum_{j=1}^{5} y_j & \leq 200 \\ \sum_{i=1}^{5} (x_i+y_j) & \geq 150 \\ x_i, y_j & \geq 0 \end{align*} 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import gurobipy as gp
from gurobipy import GRB

model3 = gp.Model('q3')
# add variables
x = model3.addVars(10,lb=0,name='x')
y = model3.addVars(5,lb=0,name='y')
a = [[1+5*i +j for j in range(5)] for i in range(10)]
#add target function
model3.setObjective(gurobipy.quicksum(a[i][j] for i in range(10) for j in range(5)), sense=GRB.MAXIMIZE)
#add constraints
model3.addConstr(gp.quicksum(x[i] for i in range(10)) <= 100)
model3.addConstr(gp.quicksum(y[j] for j in range(5)) <= 200)
model3.addConstr(gp.quicksum(x[i] + y[j] for i in range(10) for j in range(5)) >= 150)

'''
可选，添加非负约束，但是之前在addVars的时候已经添加了
model3.addConstrs((x[i] >= 0 for i in range(10)), 'nonneg')
model3.addConstrs((y[j] >= 0 for j in range(5)), 'nonneg')
'''

#optimize
model3.optimize()
display_out(model3)
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows 11+.0 (26100.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 3 rows, 15 columns and 30 nonzeros
Model fingerprint: 0x53991e8c
Coefficient statistics:
  Matrix range     [1e+00, 1e+01]
  Objective range  [0e+00, 0e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+02, 2e+02]
Presolve removed 3 rows and 15 columns
Presolve time: 0.04s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.2750000e+03   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.05 seconds (0.00 work units)
Optimal objective  1.275000000e+03
x[0] = 0.0
x[1] = 100.0
x[2] = 0.0
x[3] = 0.0
x[4] = 0.0
x[5] = 0.0
x[6] = 0.0
x[7] = 0.0
x[8] = 0.0
x[9] = 0.0
y[0] = 0.0
y[1] = 0.0
y[2] = 0.0
y[3] = 0.0
y[4] = 200.0
Objective value: 1275.0

找满足下式的整数解minf=(x2−4y)2+(1−2x)2找满足下式的整数解 \\ min f = (x^2-4y)^2+(1-2x)^2 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import gurobipy as gp   
from gurobipy import GRB    
#
model4 = gp.Model('q4')
# add variables
x = model4.addVar(vtype=GRB.INTEGER,lb = -GRB.INFINITY ,name='x')
y = model4.addVar(vtype=GRB.INTEGER,lb= -GRB.INFINITY,name='y')
t1 = model4.addVar(vtype=GRB.CONTINUOUS, name='t1')
t2 = model4.addVar(vtype=GRB.CONTINUOUS, name='t2')
# 添加非线性约束,这里需要特殊的转化，将非线性约束转化为线性约束
model4.addConstr(t1 == x*x - 4*y, name="t1_constraint")
model4.addConstr(t2 == 1 - 2*x, name="t2_constraint")
# 设置目标函数
model4.setObjective(t1**2 + t2**2, GRB.MINIMIZE)
#solve
model4.optimize()
display_out(model4)
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows 11+.0 (26100.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 1 rows, 4 columns and 2 nonzeros
Model fingerprint: 0x422167fa
Model has 2 quadratic objective terms
Model has 1 quadratic constraint
Variable types: 2 continuous, 2 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 2e+00]
  QMatrix range    [1e+00, 1e+00]
  QLMatrix range   [1e+00, 4e+00]
  Objective range  [0e+00, 0e+00]
  QObjective range [2e+00, 2e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 1e+00]
Presolve time: 0.00s
Presolved: 4 rows, 5 columns, 10 nonzeros
Presolved model has 2 quadratic objective terms
Presolved model has 1 bilinear constraint(s)
Warning: Model contains variables with very large bounds participating
         in product terms.
         Presolve was not able to compute smaller bounds for these variables.
         Consider bounding these variables or reformulating the model.

​
​ Solving non-convex MIQCP
​
​ Variable types: 2 continuous, 3 integer (0 binary)
​
Root relaxation: objective 1.000000e+00, 0 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       1.0000000    1.00000  0.00%     -    0s

Explored 1 nodes (0 simplex iterations) in 0.05 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)

Solution count 1: 1 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.000000000000e+00, best bound 1.000000000000e+00, gap 0.0000%
x = 0.0
y = 0.0
t1 = 0.0
t2 = 1.0
Objective value: 1.0
#math
线性规划
https://silenzio111.github.io/2024/05/20/lp/
作者
silenzio
发布于
2024年5月20日
许可协议