數學建模——規划模型
數學規划是運籌學的一個重要分支、而線性規划又是數學規划的一部分主要內容，所有實際問題都可以歸總為“線性規划”問題。線性規划（linear programming，LP)有比較完善的理論基礎和有效的解決方法。在實際問題中有極其廣泛的應用。
一、線性規划模型
1、建立線性規划模型的步驟
規划問題的數學模型由三個要素組成
①決策變量：是問題中要確定的未知量，用於表明規划問題中用數量表示的方案、措施等，可由決策者決定和控制
②目標函數：是決策變量的函數，優化目標通常是求該函數的最大值或最小值。
③約束條件：是決策變量的取值所受到的約束和限制條件，通常用含有決策變量的等式或不等式表示。
建立線性規划模型通常需要以下三個步驟
（1）分析問題，找到決策變量
（2）根據問題所給的條件，找到決策變量必須滿足的一組線性等式或者不等式約束，即確定約束條件
（3）根據問題的目標、構造關於決策變量的一個線性函數，即目標函數
有了決策變量、約束條件和目標函數這三個要素之后，一個線性規划模型就建立起來了。
2、線性規划模型的形式
線性規划模型的一般形式為
max(或min) z=\(c1x1+c2x2+....+c_nx_n\)
\(s.t.\left\{\begin{aligned}a_{11}x_{1}+a_{12}x_{2}+...+a_{1n}x_{n}\leq(或=、\geq)b_{1}\\a_{11}x_{1}+a_{12}x_{2}+...+a_{1n}x_{n}\leq(或=、\geq)b_{2}\\a_{m1}x_{1}+a_{m2}x_{2}+...+a_{mn}x_{n}\leq(或=、\geq)b_{m}\\x{1},x{2},...x{n}\leq0\end{aligned}\right.\)
或者簡寫為
max（或min) z=\(\displaystyle \sum^{n}_{y=1}{C_{j}x_{j}}\)
\(s.t.\left\{\begin{aligned}\\\sum^{n}_{j=1}a_{ij}\leq(或=,\geq)b_{i} ,i=1,2,3,...m,\\x_{j}\geq0,j=1,2,3,...n\end{aligned}\right.\)
其向量表示形式為
max(或min) z=CTx
\(s.t.\left\{\begin{aligned}\\\sum^{n}_{j=1}P_{ij}\leq(或=,\geq)b\\x\geq0\end{aligned}\right.\)
其矩陣表示形式為
max(或min) z=CTx
\(s.t.\left\{\begin{aligned}\\Ax\leq(或=,\geq)b\\x\geq0\end{aligned}\right.\)X
其中c=\((c_1,c2,...c_n)^T\)為目標函數的系數向量；又稱為價值向量；x=\((x_1,x_2,...x_n)^T\)稱為決策向量；A=\((a_{ij})_{mxn}\)稱為約束方程組的系數矩陣；而P\(j\)=\((a_{1j},a_{2j},a_{3j},....a_{mj})^T\),j=1,2,3,4...n,n為A的列向量，又稱為約束方程組的系數向量；b=\((b_{1j},b_{2j}....b_{mj})^T\)稱為約束方程組的常數向量。
3、用python求解線性規划模型
（1）cvxpy庫，用於解決凸優化問題，即用於求解凸規划模型,具體用法，可參考cvxpy 常用功能匯總 - 嗶哩嗶哩 (bilibili.com)
（2）scipy庫中的optimize linprog模塊，可參考scipy.optimize.linprog函數參數最全詳解_佐佑思維的博客-CSDN博客_scipy.optimize.linprog
例如：max z=\(70_{x1}+50_{x2}+60_{x3}\)
\(s.t.\left\{\begin{aligned}\\2x_{1}+4x_{2}+3x_{3}\leq150,\\3x_{1}+x_{2}+53_{3}\leq160,\\7x_{1}+3x_{2}+5x_{3}\leq200,\\x_{i}\geq0,i=1,2,3\end{aligned}\right.\)
cvxpy的具體鏈接：cvxpy 常用功能匯總 - 嗶哩嗶哩 (bilibili.com)
# 使用cvxpy來實現線性規划問題
import cvxpy as cp
import numpy as np
c=np.array([70,50,60]) # 定義目的向量
a=np.array([[2,4,3],[3,1,5],[7,3,5]]) #定義約束矩陣
b=np.array([150,160,200]) # 定義約束條件的右邊向量
x=cp.Variable(3,pos=True) # 定義3個決策變量(x1,x2,x3)
obj=cp.Maximize(c@x) # 構造約束條件最大值(70x1+50x2+60x3)
cons=[a@x<=b] # 構造約束條件
pro=cp.Problem(obj,cons) # 求解問題
pro.solve(solver='GLPK_MI',verbose=True)
#用不同的求解器求出來的幾個是不一樣的，用GLPKMI解決問題，verbose=True表明顯示具體求解信息
print('最優解是',x.value)
print('最優值是',pro.value)
# Minimize 表示最小值，Maximize表示最大值，出現等於的約束條件可以將約束條件卸載cons里面
# 輸出的結果：
# 最優的解 [15.90909091 29.54545455 -0. ]
# 最優的值 2590.909090909091
cvxpy的求解器(solver)
LP:指線性規划 QP:指二次規划（二次函數） SOCP:指二次錐規划
SDP:半正定規划 EXP:指數規划 POW:冪規划 MIP:混合整數規划
scipy的具體鏈接：scipy.optimize.linprog函數參數最全詳解_佐佑思維的博客-CSDN博客_scipy.optimize.linprog
# 使用scipy庫實現線性規划
from scipy import optimize as op
import numpy as np
c=np.array([-70,50,60]) # 目標函數
A_ub=np.array([[2,4,3],[3,1,5],[7,3,5]]) # 約束條件的左邊
B_ub=np.array([[150,160,200]]) # 約束條件的右邊
x1=(0,None) # 三個變量及其范圍(0到正無窮）
x2=(0,None)
x3=(0,None)
res=op.linprog(-c,A_ub,B_ub,bounds=(x1,x2,x3),method='revised simplex') # 使用method='revised simplex'高精度運算
print(res)
# 輸出的結果：
# con: array([], dtype=float64)
# fun: -2590.909090909091
# message: 'Optimization terminated successfully.'
# nit: 2
# slack: array([ 0. , 82.72727273, 0. ])
# status: 0
# success: True
# x: array([15.90909091, 29.54545455, 0. ])
# 求解的值 2590.909090909091
# 求解的x [15.90909091 29.54545455 0. ]
二、整數線性規划
1、整數規划的類型
在一些問題（例如汽車企業在制定生產計划是，要求所生產的不同類型的汽車數量必須是整數）等在問題上，我們必須要求一部分或者全部決策變量必須是整數值的規划問題稱為整數規划（interger programming IP)。
從決策變量的取值范圍來看，整數規划通常分為以下幾種類型
（1）純整數規划：全部決策變量都必須整數值的整數規划模型
（2）混合整數規划：決策變量中優一部分必須去整數值，另一部分可以不取整數的整數規划模型
（3）0-1整數規划：決策變量只能取0或1的整數規划
2、整數線性規划模型一般形式為
max（或min) z=\(\displaystyle \sum^{n}_{y=1}{C_{j}x_{j}}\)
\(s.t.\left\{\begin{aligned}\sum^n_{j=1}a_{ij}x{j}\leq(或=,\geq)b_{i},i=1,2,....m\\x_{j}\geq0,j=1,2,...n\\x_1,x_2,....x_n 中部分或者全部取整數\end{aligned}\right.\)
3、python求解模型
整數處理方式：
1、純整數規划問題及混合整數規划問題
x=cp.Variable(3,pos=True,interger=True) # interger 化為整數形式
2、0-1整數規划
# 在約束條件中加入 0<=x,x<=1
x=cp.Variable(3,pos=True,interger=True)
cons=[a@x<=b,0<=x,x<=1]
例題：
min z=\(\displaystyle\sum^4_{i=1}\sum^5_{j=1}c_{ij}x{ij}\)
\(s.t.\left\{\begin{aligned}\sum^4_{i=1}a_{ij}\leq=1,i=1,2,....5\\\sum^5_{i=1}x_{ij}\leq0,j=1,2,3,4\\x_{ij}=0或1,i=1,2,3,4,j=1,2,...5\end{aligned}\right.\)
cvxpy的解法
import cvxpy as cp
import numpy as np
c=np.array([[15,13.8,12.5,11,14.3],[14.5,14,13.2,10.5,15],[13.8,13,12.8,11.3,14.6],[14.,13.6,13,11.6,14]])
x=cp.Variable((4,5),integer=True) # 定義目標變量，全為整數
obj=cp.Minimize(cp.sum(cp.multiply(c,x))) # 求解最小值 構造目標函數
con=[cp.sum(x,axis=1)<=2,cp.sum(x,axis=0)==1,0<=x,x<=1] # 構造約束條件
# cp.sum(x,axis=1)表示在每一列中不能大於2，0<=x,x<=1,表示0-1整數問題處理方法
pro=cp.Problem(obj,con)
pro.solve(solver='GLPK_MI')
print('最優解：',pro.value)
print('最優值',x.value)
三、非線性規划模型
我們遇到的規划模型中的目標函數和約束條件均為決策變量的線性函數，我們稱其為線性規划模型。如果在目標函數或者約束條件中含有決策變量的非線性函數，則該模型為非線性模型（nonlinear progamming,NIP)。
1、非線性規划模型形式
非線性規划模型的一般形式為
min f(x)
\(s.t.\left\{\begin{aligned}g_i(x)\leq0,i=1,2,...,p\\h_j(x)\geq0,i=1,2,...q\end{aligned}\right.\)
其中x=\((x_1,x_2,...x_n)^T\)是n維歐氏空間\(R^n\)中的一個n位向量，函數f（x),\(g_i(x)\),i=1,2,...,p和\(hi(x)\),j=1,2...q是定義在\(R^n\)上的實值函數，且其中至少一個是x的非線性函數。
特別，若在模型中，p=q=0，我們稱其為無約束非線性規划或無約束優化，即求多元函數的極值問題；若目標函數f(x)是一個二次函數，而\(g_i(x)\),i=1,2,...p和\(h_j(x)\),j=1,2,...q都是線性函數，我們稱其為二次歸規划。若目標函數f(x)為凸函數,則均為線性函數，我們稱其為凸規划。
2、用python求解非線性規划模型
對於一般的非線性問題，由於不是凸規划，就不能用cvxpy庫求解，求解一般的非線性規划問題使用scipy.optimize模塊的minimize函數求解。[minimize的使用](python 非線性規划方式(scipy.optimize.minimize)_python_腳本中心 - 編程客棧 (cppcns.com))
min T=\(\displaystyle\sum^3_{i_=1}t_i\)
\(s.t.\left\{\begin{aligned}(t1+t3)1.5cosa_1=400\\t_21.5a_2=760\\(1.47-1.5sina_1)(t_1+t_3)+(2.11-1.5sina_1)t2=1000\\t1=t3\\t_i\geq0,i=1,2,3\\o<a_i<Π/2,i=1,2\end{aligned}\right.\)
import numpy as np
from scipy.optimize import minimize
obj=lambda x:2*x[0]+x[1] # 定義目標的函數
con=[{'type':'eq','fun':lambda x:2*x[0]*1.5*np.cos(x[2]-400)},
{'type':'eq','fun':lambda x:x[1]*1.5*np.cos(x[3])-760},
{'type':'eq','fun':lambda x:(1.46-1.5*np.sin(x[2]))*2*x[0]+(2.11-1.5*np.sin(x[3]))*x[1]-1000}] # 定義約束條件
b1=[0]*4
b2=[np.inf,np.inf,np.pi/2,np.pi/2]
bd=list(zip(b1,b2)) # 將兩個綁定在一起，以元組的形式
x=minimize(obj,np.random.rand(4),constraints=con,bounds=bd)
# 根據minimize的格式，np.random.rand(4)表示初始化四個數，constraints約束條件，method表示方法
print(x)
print('角度',x.x[2:]*180/np.pi)