[Python知识库]用Python求解数学规划问题

一、线性规划简介

线性规划(Linear Programming 简记为LP)是数学规划的一个重要分支。

规划问题分类

• 线性规划： 在一组线性约束条件的限制下，求一线性目标函数最大或最小的问题；
• 整数规划: 当约束条件加强，要求所有的自变量必须是整数时，成为整数规划（特别地，自变量只能为0或1时称为0-1规划）；
• 非线性规划: 无论是约束条件还是目标函数出现非线性项，那么规划问题就变成了非线性规划；
• 多目标规划: 在一组约束条件的限制下，求多个目标函数最大或最小的问题；
• 动态规划: 将优化目标函数分多阶段，利用阶段间的关系逐一进行求解的方法。

应用举例

旅行商问题、车辆路径规划问题、运输问题、最短路问题、最大流问题、中国邮递员问题等。

二、线性规划模型的三要素

线性规划模型主要包括三个部分：决策变量、目标函数、约束条件。

决策变量

决策变量是指问题中可以改变的量，例如生产多少货物，选择哪条路径等；线性规划的目标就是找到最优的决策变量。

在线性规划中决策变量包括实数变量，整数变量，0-1变量等。

目标函数

目标函数就是把问题中的决策目标量化，一般分为最大化目标函数和最小化目标函数。
在线性规划中，目标函数为一个包含决策变量的线性函数，约束条件是指问题中各种时间，空间，人力，物力等限制。

约束条件

在线性规划中约束条件一般表示为一组包含决策变量的不等式和等式。

三、线性规划模型的数学表示

一般形式

线性规划模型可以写成如下形式：
$$\begin{gathered} minz=x_1+x_2 \\ s.t.?\begin{array}{rl}{x}_1+2x_2& \le 1\\ 4x_1+3x_2& \le 2\\ x_1,x_2& \ge 0\end{array} \end{gathered}$$

矩阵形式

上述模型也可以写成如下的矩阵形式：
$$\begin{gathered} minz=c^{T}x \\ s.t.\{\begin{array}{rl}Ax& \le b\\ x& \ge 0\end{array} \end{gathered}$$
对于有 $n$ 个决策变量, $m$ 个约束的线性规划模型, $c,x$ 为 $n$ 维列向量，$b$ 为 $m$ 维列向量, $A$ 为 $m\times n$ 维矩阵。

标准形式

线性规划的目标函数可能是最大化，也可能是最小化，约束条件的符号可能是小于等于，也可能是大于等于。因此为了编程方便，一般统一为最小化目标函数，小于等于约束。

• 最大化目标函数可以添加负号变为最小化约束；
• 大于等于约束可以两边乘以负号变为小于等于约束；
• 等于约束可以变为一个大于等于约束和一个小于等于约束，但在编程中一般支持直接写等式约束,可以不进行转换。

四、线性规划模型的求解方法

方法总览

线性规划模型的求解方法主要有：图解法、单纯形法、椭球法、卡玛卡算法、内点法等。

其中内点法因为求解效率更高，在决策变量多，约束多的情况下能取得更好的效果，目前主流线性规划求解器都是使用的内点法。

使用scipy求解

对上面提到的线性优化问题
$$\begin{gathered} minz=x_1+x_2 \\ s.t.?\begin{array}{rl}{x}_1+2x_2& \le 1\\ 4x_1+3x_2& \le 2\\ x_1,x_2& \ge 0\end{array} \end{gathered}$$
使用scipy求解的步骤如下：

Step1 : 导入相关库

import numpy as np 
from scipy import optimize as op

Step2: 定义决策变量

# 给出变量取值范围
x1 = (0, None)
x2 = (0, None)

Step3: 将原问题化为标准形式

编程时默认为最小化目标函数，约束为小于等于约束。

Step4: 定义目标函数系数和约束条件系数

c = np.array([1, 1])  # 目标函数系数,2x1列向量

A = np.array([[1, 2], [4, 3]])  # 不等式约束系数A，2x2维矩阵
b = np.array([1, 2])  # 不等式约束系数b，2x1列向量

Step5: 求解

res = op.linprog(c, A, b, bounds=(x1, x2))  # 调用函数进行求解
res

输出结果：

    con: array([], dtype=float64)
     fun: 5.428987546487586e-11
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([1., 2.])
  status: 0
 success: True
       x: array([3.02272240e-11, 2.40626514e-11])

求解实例

例1：等式不等式约束

$$\begin{gathered} maxz=2x_1+3x_2?5x_3 \\ s.t.?\begin{array}{rl}{x}_1+x_2+x_3& =7\\ 2x_1?5x_2+x_3& \ge 10\\ x_1+3x_2+x_3& \le 12\\ x_1,x_2,x_3& \ge 0\end{array} \end{gathered}$$

转化为标准形式：
$$\begin{gathered} min?z=?2x_1?3x_2+5x_3 \\ s.t.?\begin{array}{rl}{x}_1+x_2+x_3& =7\\ ?2x_1+5x_2?x_3& \le ?10\\ x_1+3x_2+x_3& \le 12\\ x_1,x_2,x_3& \ge 0\end{array} \end{gathered}$$
对标准形式的问题编程求解

import numpy as np
from scipy import optimize as op

# 给出变量取值范围
x1 = (0, None)
x2 = (0, None)
x3 = (0, None)

c = np.array([-2, -3, 5])  # 目标函数系数,3x1列向量

A_ub = np.array([[-2, 5, -1], [1, 3, 1]])  # 不等式约束系数A，2x3维矩阵
B_ub = np.array([-10, 12])  # 不等式约束系数b, 2x1维列向量
A_eq = np.array([[1, 1, 1]])  # 等式约束系数Aeq，3x1维列向量
B_eq = np.array([7])  # 等式约束系数beq，1x1数值

res = op.linprog(c, A_ub, B_ub, A_eq, B_eq, bounds=(x1, x2, x3))  # 调用函数求解
res

输出结果：

     con: array([1.80714288e-09])
     fun: -14.57142856564504
 message: 'Optimization terminated successfully.'
     nit: 5
   slack: array([-2.24614993e-10,  3.85714286e+00])
  status: 0
 success: True
       x: array([6.42857143e+00, 5.71428571e-01, 2.35900788e-10])

即当 $x_1=6.43,x_2=5.71,x_3=0,$ 时，目标函数取得最小值 $?z=?14.57$，从而对于原问题而言，当 $x_1=6.42,x_2=0.57,x_3=0,$ 时，目标函数取得最大值 $z=14.57$。

例2：包含非线性项

$$\begin{gathered} minf(x)=x_1^2+x_2^2+x_3^2+8 \\ s.t.?\begin{array}{rl}{x}_1^2?x_2+x_3^2& \ge 0\\ x_1+x_2^2+x_3^2& \le 20\\ ?x_1?x_2^2+2& =0\\ x_2+2x_3^2& =3\\ x_1,x_2,x_3& \ge 0\end{array} \end{gathered}$$

由于存在非线性项，不能沿用例1中的 linprog 函数求解，这里使用自定义函数的方法编写目标函数和约束条件，并使用 scipy.optimize 中的 minimize 函数求解。

Step1:导入相关库

import numpy as np
from scipy.optimize import minimize

Step2:使用函数的形式表示目标和约束

# 定义目标函数
def objective(x):
    return x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + 8

# 定义约束条件
def constraint1(x):
    return x[0] ** 2 - x[1] + x[2] ** 2  # 不等约束

def constraint2(x):
    return -(x[0] + x[1] ** 2 + x[2] ** 2 - 20)  # 不等约束

def constraint3(x):
    return -x[0] - x[1] ** 2 + 2  # 等式约束

def constraint4(x):
    return x[1] + 2 * x[2] ** 2 - 3  # 等式约束

注：每一个函数的输入为一个 n 维列向量 x ，其中 x[0] 表示该列向量的第一个元素，即 $x_1$ 。

Step3:定义约束条件

con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}

# 4个约束条件
cons = ([con1, con2, con3, con4])

# 决策变量的符号约束
b = (0.0, None)  # 即决策变量的取值范围为大于等于0
bnds = (b, b, b)

注：每一个约束为一个字典，其中 type 表示约束类型：ineq 为大于等于，eq 为等于；fun 表示约束函数表达式，即 step2 中的自定义函数。

Step4:求解

x0 = np.array([0, 0, 0])  # 定义初始值
solution = minimize(objective, x0, method='SLSQP', \
                    bounds=bnds, constraints=cons)

注：minimize为最小化目标函数，且约束条件中默认为大于等于约束。

Step5：打印求解结果

x = solution.x
print('目标值: ' + str(objective(x)))
print('最优解为')
print('x1 = ' + str(round(x[0], 2)))
print('x2 = ' + str(round(x[1], 2)))
print('x3 = ' + str(round(x[2], 2)))
solution

输出结果为：

目标值: 10.651091840572583
最优解为
x1 = 0.55
x2 = 1.2
x3 = 0.95

     fun: 10.651091840572583
     jac: array([1.10433471, 2.40651834, 1.89564812])
 message: 'Optimization terminated successfully'
    nfev: 71
     nit: 15
  status: 0
 success: True
       x: array([0.55216734, 1.20325918, 0.94782404])

即当 $x_1=0.55,x_2=1.20,x_3=0.95,$ 时，目标函数取得最小值 $z=10.65$。

调用求解器求解

对线性优化问题
$$\begin{gathered} minz=x_1+x_2 \\ s.t.?\begin{array}{rl}{x}_1+2x_2& \le 1\\ 4x_1+3x_2& \le 2\\ x_1,x_2& \ge 0\end{array} \end{gathered}$$
使用 python 调用求解器 cplex （该求解器一般用于求解线性问题）求解的步骤如下：

Step1 : 导入相关库

import numpy as np
from pyomo.environ import *
import pyutilib.subprocess.GlobalData
pyutilib.subprocess.GlobalData.DEFINE_SIGNAL_HANDLERS_DEFAULT = False

Step2: 定义目标函数及约束条件

# 定义目标函数
def objective(model):
    return model.x1 + model.x2

# 定义约束条件
def constraint(model):
    model.cons.add(expr=model.x1 + 2 * model.x2 <= 1)  # 必须为非严格不等式或等式
    model.cons.add(expr=4 * model.x1 + 3 * model.x2 <= 2)

Step3: 建立问题模型

model = ConcreteModel(name="optimize_prob")  # 这个名字随便起
# 决策变量定义
model.x1 = Var(bounds=(0,None),within=NonNegativeReals,initialize=0.1)  
model.x2 = Var(bounds=(0,None),within=NonNegativeReals,initialize=0.2)

model.cons = ConstraintList()  
constraint(model)
model.obj = Objective(rule=objective, sense=minimize)  # 目标函数

solver_path = "D:\\algorithm_tools\\solver\\cplex\\cplex"  # 求解器路径
opt = SolverFactory('cplex', executable=solver_path)       # 调用求解器求解

Step4: 求解问题

results = opt.solve(model, tee=True)  # tee=True 为打印求解过程
results.write()
print(model.x1.value)  # 打印最优解
print(model.x2.value)
print(model.obj())     # 打印最优值

最终的结果为：

# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: tmph0hv701k
  Lower bound: 0.0
  Upper bound: 0.0
  Number of objectives: 1
  Number of constraints: 3
  Number of variables: 3
  Number of nonzeros: 5
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: 0.0
  Termination condition: optimal
  Termination message: Dual simplex - Optimal\x3a Objective = 0.0000000000e+00
  Error rc: 0
  Time: 0.03191494941711426
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
0.0
0.0
0.0

求解实例

例1：等式不等式约束

$$\begin{gathered} maxz=2x_1+3x_2?5x_3 \\ s.t.?\begin{array}{rl}{x}_1+x_2+x_3& =7\\ 2x_1?5x_2+x_3& \ge 10\\ x_1+3x_2+x_3& \le 12\\ x_1,x_2,x_3& \ge 0\end{array} \end{gathered}$$

求解过程如下：

import numpy as np
from pyomo.environ import *
import pyutilib.subprocess.GlobalData
pyutilib.subprocess.GlobalData.DEFINE_SIGNAL_HANDLERS_DEFAULT = False

# 定义目标函数
def objective(model):
    return 2 * model.x1 + 3 * model.x2 - 5 * model.x3


# 定义约束条件
def constraint(model):
    model.cons.add(expr=model.x1 + model.x2 + model.x3 == 7)  # 必须为非严格不等式或等式
    model.cons.add(expr=2 * model.x1 - 5 * model.x2 + model.x3 >= 10)
    model.cons.add(expr=model.x1 + 3 * model.x2 + model.x3 <= 12)


model = ConcreteModel(name="optimize_prob")  # 这个名字随便起
model.x1 = Var(bounds=(0, None), within=NonNegativeReals, initialize=0.1)  
model.x2 = Var(bounds=(0, None), within=NonNegativeReals, initialize=0.2)
model.x3 = Var(bounds=(0, None), within=NonNegativeReals, initialize=0.5)

model.cons = ConstraintList()  
constraint(model)
model.obj = Objective(rule=objective, sense=maximize)  # 目标函数

solver_path = "D:\\solver\\cplex\\cplex"  # 求解器路径
opt = SolverFactory('cplex', executable=solver_path)  # 调用求解器求解

results = opt.solve(model, tee=True)  # tee=True 打印求解过程
results.write()
print("x1 = " + str(model.x1.value))  # 打印最优解
print("x2 = " + str(model.x2.value))
print("x3 = " + str(model.x3.value))
print("optimize value = " + str(model.obj()))  # 打印最优值

最终的输出结果为：

# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: tmp78fzob_g
  Lower bound: 14.571428571428575
  Upper bound: 14.571428571428575
  Number of objectives: 1
  Number of constraints: 4
  Number of variables: 4
  Number of nonzeros: 10
  Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: 0.0
  Termination condition: optimal
  Termination message: Dual simplex - Optimal\x3a Objective = 1.4571428571e+01
  Error rc: 0
  Time: 0.033006906509399414
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x1 = 6.428571428571429
x2 = 0.5714285714285716
x3 = 0.0
optimize value = 14.571428571428573

即当 $x_1=6.43,x_2=0.57,x_3=0$ 时，目标函数取得最大值 $14.57$。

例2：包含非线性项

$$\begin{gathered} minf(x)=x_1^2+x_2^2+x_3^2+8 \\ s.t.?\begin{array}{rl}{x}_1^2?x_2+x_3^2& \ge 0\\ x_1+x_2^2+x_3^2& \le 20\\ ?x_1?x_2^2+2& =0\\ x_2+2x_3^2& =3\\ x_1,x_2,x_3& \ge 0\end{array} \end{gathered}$$

由于存在非线性项，不能沿用例1中的 cplex 求解器求解，这里使用可求解非线性规划问题的求解器 bonmin 对该问题进行求解。代码如下

from pyomo.environ import *
import pyutilib.subprocess.GlobalData
pyutilib.subprocess.GlobalData.DEFINE_SIGNAL_HANDLERS_DEFAULT = False

# 定义目标函数
def objective(model):
    return model.x1 ** 2 + model.x2 ** 2 + model.x3 ** 2 + 8

# 定义约束条件
def constraint(model):
    model.cons.add(expr=model.x1 ** 2 - model.x2 + model.x3 ** 2 >= 0)  # 必须为非严格不等式或等式
    model.cons.add(expr=model.x1 + model.x2 ** 2 + model.x3 ** 2 <= 20)  
    model.cons.add(expr=-model.x1 - model.x2 ** 2 + 2 == 0) 
    model.cons.add(expr=model.x2 + 2 * model.x3 ** 2 == 3)  

model = ConcreteModel(name="optimize_prob") 
model.x1 = Var(bounds=(0, None), within=NonNegativeReals, initialize=0.1)  # 决策变量定义
model.x2 = Var(bounds=(0, None), within=NonNegativeReals, initialize=0.2)
model.x3 = Var(bounds=(0, None), within=NonNegativeReals, initialize=0.5)

model.cons = ConstraintList() 
constraint(model)
model.obj = Objective(rule=objective, sense=minimize)  # 目标函数

solver_path = "D:\\solver\\bonmin\\bonmin"  # 求解器路径
opt = SolverFactory('bonmin', executable=solver_path)  # 调用求解器进行求解

results = opt.solve(model, tee=True)  # tee=True 为打印求解过程
results.write()
print("x1 = " + str(model.x1.value))  # 打印最优解
print("x2 = " + str(model.x2.value))
print("x3 = " + str(model.x3.value))
print("optimize value = " + str(model.obj()))  # 打印最优值

最终的输出结果为：

x1 = 0.5521673357043836
x2 = 1.2032591841796418
x3 = 0.9478240384753155
optimize value = 10.651091838843191

即当 $x_1=0.55,x_2=1.20,x_3=0.95$ 时，目标函数取得最大值 $10.65$。

注：在工程上，对于大规模非线性问题，可通过先对其线性化处理转化为线性规划问题，再调用cplex求解器的方法求解该问题。

五、常用求解器简介

cplex

官网：https://www.ibm.com/cn-zh/analytics/cplex-optimizer

Cplex是IBM公司开发的一款商业版的优化引擎，也有免费版，但免费版的有规模限制，不能求解规模过大的问题。

Cplex专门用于求解大规模的线性规划（LP）、二次规划（QP）、带约束的二次规划（QCQP）、二阶锥规划（SOCP）等四类基本问题，以及相应的混合整数规划（MIP）问题。有以下优势：

• 能解决一些非常困难的行业问题；
• 求解速度非常快；
• 提供超线性加速功能的优势。

gurobi

官网：https://www.gurobi.com/ 中文官网：http://www.gurobi.cn/

Gurobi是由美国Gurobi公司开发的新一代大规模数学规划优化器，适用于LP、QP等场景，提供了C，C++，java，python, MATLAB, R语言等多种语言的接口。

bonmin

官网：https://www.coin-or.org/Bonmin/index.html

BONMIN （基本开源非线性混合整数规划）是用于解决一般 MINLP（混合整数非线性规划）问题的开源代码BO，NMIN 是 开源软件。

scip

官网：https://www.scipopt.org/

SCIP 是目前用于混合整数规划 (MIP) 和混合整数非线性规划 (MINLP) 的最快的非商业求解器之一。它也是一个用于约束整数规划和分支切割和定价的框架，允许完全控制求解过程并访问详细信息，直到求解器的内部。

六、总结

本文在线性规划问题的基本概念之上，简单介绍了利用 python求解线性规划问题以及简单的非线性规划问题的两种方法，即使用 scipy 模块和调用合适的求解器，并附有详细的操作步骤，最后简单介绍了一些常用的求解器。