用Python求解數學規劃問題

一、線性規劃簡介

線性規劃(Linear Programming 簡記為LP)是數學規劃的一個重要分支。

規劃問題分類

線性規劃：在一組線性限制條件的限制下，求一線性目标函數最大或最小的問題；
整數規劃: 當限制條件加強，要求所有的自變量必須是整數時，成為整數規劃（特别地，自變量隻能為0或1時稱為0-1規劃）；
非線性規劃: 無論是限制條件還是目标函數出現非線性項，那麼規劃問題就變成了非線性規劃；
多目标規劃: 在一組限制條件的限制下，求多個目标函數最大或最小的問題；
動态規劃: 将優化目标函數分多階段，利用階段間的關系逐一進行求解的方法。

應用舉例

旅行商問題、車輛路徑規劃問題、運輸問題、最短路問題、最大流問題、中國郵差問題等。

二、線性規劃模型的三要素

線性規劃模型主要包括三個部分：決策變量、目标函數、限制條件。

決策變量

決策變量是指問題中可以改變的量，例如生産多少貨物，選擇哪條路徑等；線性規劃的目标就是找到最優的決策變量。

線上性規劃中決策變量包括實數變量，整數變量，0-1變量等。

目标函數

目标函數就是把問題中的決策目标量化，一般分為最大化目标函數和最小化目标函數。

線上性規劃中，目标函數為一個包含決策變量的線性函數，限制條件是指問題中各種時間，空間，人力，物力等限制。

限制條件

線上性規劃中限制條件一般表示為一組包含決策變量的不等式和等式。

三、線性規劃模型的數學表示

一般形式

線性規劃模型可以寫成如下形式：

m i n z = x 1 + x 2 s . t . { x 1 + 2 x 2 ≤ 1 4 x 1 + 3 x 2 ≤ 2 x 1 , x 2 ≥ 0 min\quad\quad\quad z=x_1+x_2 \\ s.t. \quad \left\{ \begin{aligned} x_1+2x_2&\le1\\ 4x_1+3x_2&\le2\\ x_1,x_2&\ge0\\ \end{aligned} \right. minz=x1+x2s.t.⎩⎪⎨⎪⎧x1+2x24x1+3x2x1,x2≤1≤2≥0

矩陣形式

上述模型也可以寫成如下的矩陣形式：

m i n z = c T x s . t . { A x ≤ b x ≥ 0 min\quad\quad z=c^Tx \\ s.t. \quad \left\{ \begin{aligned} Ax&\le b\\ x&\ge0\\ \end{aligned} \right. minz=cTxs.t.{Axx≤b≥0

對于有 n n n 個決策變量, m m m 個限制的線性規劃模型, c , x c, x c,x 為 n n n 維列向量， b b b 為 m m m 維列向量, A A A 為 m × n m \times n m×n 維矩陣。

标準形式

線性規劃的目标函數可能是最大化，也可能是最小化，限制條件的符号可能是小于等于，也可能是大于等于。是以為了程式設計友善，一般統一為最小化目标函數，小于等于限制。

最大化目标函數可以添加負号變為最小化限制；
大于等于限制可以兩邊乘以負号變為小于等于限制；
等于限制可以變為一個大于等于限制和一個小于等于限制，但在程式設計中一般支援直接寫等式限制,可以不進行轉換。

四、線性規劃模型的求解方法

方法總覽

線性規劃模型的求解方法主要有：圖解法、單純形法、橢球法、卡瑪卡算法、内點法等。

其中内點法因為求解效率更高，在決策變量多，限制多的情況下能取得更好的效果，目前主流線性規劃求解器都是使用的内點法。

使用scipy求解

對上面提到的線性優化問題

使用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：等式不等式限制

m a x z = 2 x 1 + 3 x 2 − 5 x 3 s . t . { x 1 + x 2 + x 3 = 7 2 x 1 − 5 x 2 + x 3 ≥ 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 , x 2 , x 3 ≥ 0 max\quad\quad z=2x_1+3x_2-5x_3 \\ s.t. \quad \left\{ \begin{aligned} x_1+x_2+x_3&=7\\ 2x_1-5x_2+x_3&\ge10\\ x_1+3x_2+x_3&\le12\\ x_1,x_2,x_3&\ge0\\ \end{aligned} \right. maxz=2x1+3x2−5x3s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1+x2+x32x1−5x2+x3x1+3x2+x3x1,x2,x3=7≥10≤12≥0

轉化為标準形式：

m i n − z = − 2 x 1 − 3 x 2 + 5 x 3 s . t . { x 1 + x 2 + x 3 = 7 − 2 x 1 + 5 x 2 − x 3 ≤ − 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 , x 2 , x 3 ≥ 0 min\quad\quad -z=-2x_1-3x_2+5x_3 \\ s.t. \quad \left\{ \begin{aligned} x_1+x_2+x_3&=7\\ -2x_1+5x_2-x_3&\le-10\\ x_1+3x_2+x_3&\le12\\ x_1,x_2,x_3&\ge0\\ \end{aligned} \right. min−z=−2x1−3x2+5x3s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1+x2+x3−2x1+5x2−x3x1+3x2+x3x1,x2,x3=7≤−10≤12≥0

對标準形式的問題程式設計求解

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 , x_1=6.43, x_2=5.71, x_3=0, x1=6.43,x2=5.71,x3=0, 時，目标函數取得最小值 − z = − 14.57 -z=-14.57 −z=−14.57，進而對于原問題而言，當 x 1 = 6.42 , x 2 = 0.57 , x 3 = 0 , x_1=6.42, x_2=0.57, x_3=0, x1=6.42,x2=0.57,x3=0, 時，目标函數取得最大值 z = 14.57 z=14.57 z=14.57。

例2：包含非線性項

m i n f ( x ) = x 1 2 + x 2 2 + x 3 2 + 8 s . t . { x 1 2 − x 2 + x 3 2 ≥ 0 x 1 + x 2 2 + x 3 2 ≤ 20 − x 1 − x 2 2 + 2 = 0 x 2 + 2 x 3 2 = 3 x 1 , x 2 , x 3 ≥ 0 min\quad f(x)=x_1^2+x_2^2+x_3^2+8 \\ s.t. \quad \left\{ \begin{aligned} x_1^2-x_2+x_3^2&\ge0\\ x_1+x_2^2+x_3^2&\le20\\ -x_1-x_2^2+2&=0\\ x_2+2x_3^2&=3\\ x_1,x_2,x_3&\ge0\\ \end{aligned} \right. minf(x)=x12+x22+x32+8s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧x12−x2+x32x1+x22+x32−x1−x22+2x2+2x32x1,x2,x3≥0≤20=0=3≥0

由于存在非線性項，不能沿用例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 x_1 x1 。

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 , x_1=0.55, \; x_2=1.20, \; x_3=0.95, x1=0.55,x2=1.20,x3=0.95, 時，目标函數取得最小值 z = 10.65 z=10.65 z=10.65。

調用求解器求解

對線性優化問題

使用 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：等式不等式限制

求解過程如下：

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 x_1=6.43, \; x_2=0.57,\; x_3=0 x1=6.43,x2=0.57,x3=0 時，目标函數取得最大值 14.57 14.57 14.57。

例2：包含非線性項

由于存在非線性項，不能沿用例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 x_1=0.55, \; x_2=1.20,\; x_3=0.95 x1=0.55,x2=1.20,x3=0.95 時，目标函數取得最大值 10.65 10.65 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 子產品和調用合适的求解器，并附有詳細的操作步驟，最後簡單介紹了一些常用的求解器。

用Python求解數學規劃問題

一、線性規劃簡介

規劃問題分類

應用舉例

二、線性規劃模型的三要素

決策變量

目标函數

限制條件

三、線性規劃模型的數學表示

一般形式

矩陣形式

标準形式

四、線性規劃模型的求解方法

方法總覽

使用scipy求解

求解執行個體

例1：等式不等式限制

例2：包含非線性項

調用求解器求解

求解執行個體

例1：等式不等式限制

例2：包含非線性項

五、常用求解器簡介

cplex

gurobi

bonmin

scip

六、總結

繼續閱讀