线性规划

解决线性规划问题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

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))

PYTHON

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
TEXT

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

$$min(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*}$$

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)
PYTHON

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
TEXT

$$min(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*}$$

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)
PYTHON

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
TEXT

$$max \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*}$$

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)
PYTHON

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
TEXT

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

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)
PYTHON

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.
TEXT

​
​ 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
TEXT