Crossin的程式設計教室 2020-03-29
以下文章來源于邯鄲路220号子彬院 ，作者少年吉
邯鄲路220号子彬院
分享資料分析、Python、ML、DL，機率論，數理統計等相關知識。做一個有情懷的公衆号！
線性規劃簡介及數學模型表示線性規劃簡介一個典型的線性規劃問題線性規劃模型的三要素線性規劃模型的數學表示圖解法和單純形法圖解法單純形法使用python求解簡單線性規劃模型程式設計思路求解案例例1：使用scipy求解例2：包含非線性項的求解從整數規劃到0-1規劃整數規劃模型0-1規劃模型案例：投資的收益和風險問題描述與分析建立與簡化模型
線性規劃簡介及數學模型表示
線性規劃簡介
在人們的生産實踐中，經常會遇到如何利用現有資源來安排生産，以取得最大經濟效益的問題。此類問題構成了運籌學的一個重要分支—數學規劃，而線性規劃(Linear Programming 簡記LP)則是數學規劃的一個重要分支，也是一種十分常用的最優化模型。
而随着計算機的發展，線性規劃的方法被應用于廣泛的領域，已成為數學模組化裡最為經典，最為常用的模型之一。線性規劃模型可用于求解利潤最大，成本最小，路徑最短等最優化問題。
一個典型的線性規劃問題
某機床廠生産甲、乙兩種機床，每台銷售後的利潤分别為4000元與3000元。生産甲機床需用A、B機器加工，加工時間分别為每台2小時和1小時：生産乙機床需用A、B、C三種機器加工，加工時間為每台各一小時。若每天可用于加工的機器時數分别為機器10小時、機器8小時和C機器7小時，問該廠應生産甲、乙機床各幾台，才能使總利潤最大？
問題分析
這個問題是一個十分典型的線性規劃問題，首先對問題提取出關鍵資訊：
決策：生産幾台甲、乙機床
優化目标：總利潤最大
限制：生産機床的使用時間有限
将上訴三個要素寫成數學表達式，就是一個典型的線性規劃模型：
規劃問題的分類
- 線性規劃： 在一組線性限制條件的限制下，求一線性目标函數最大或最小的問題；
- 整數規劃：當限制條件加強，要求所有的自變量必須是整數時，成為整數規劃（特别地，自變量隻能為0或1時稱為0-1規劃）；
- 非線性規劃：無論是限制條件還是目标函數出現非線性項，那麼規劃問題就變成了非線性規劃；
- 多目标規劃：在一組限制條件的限制下，求多個目标函數最大或最小的問題；
- 動态規劃：将優化目标函數分多階段，利用階段間的關系逐一進行求解的方法；
應用舉例：旅行商問題、車輛路徑規劃問題、運輸問題、最短路問題、最大流問題、中國郵差問題
線性規劃模型的三要素
線性規劃模型主要包括三個部分：決策變量、目标函數、限制條件
決策變量
決策變量是指問題中可以改變的量，例如生産多少貨物，選擇哪條路徑等；線性規劃的目标就是找到最優的決策變量。
線上性規劃中決策變量包括實數變量，整數變量，0-1變量等。
目标函數
目标函數就是把問題中的決策目标量化，一般分為最大化目标函數和最小化目标函數
線上性規劃中，目标函數為一個包含決策變量的線性函數，例如
限制條件是指問題中各種時間，空間，人力，物力等限制。
線上性規劃中限制條件一般表示為一組包含決策變量的不等式，例如
此外，決策變量的取值範圍稱為符号限制，例如
線性規劃模型的數學表示
線性規劃模型可以寫成如下形式：
其中 為subject to的縮寫。
上訴模型也可以寫成如下的矩陣形式：
對于有 個決策變量， 個限制的線性規劃模型, 為 維列向量， 為 維列向量, 為 維矩陣。
線性規劃模型的标準形式
線性規劃的目标函數可能是最大化，也可能是最小化，限制條件的符号可能是小于等于，也可能是大于等于。
是以為了程式設計友善，一般統一為最小化目标函數，小于等于限制
最大化目标函數可以添加負号變為最小化限制：
大于等于限制可以兩邊乘以 變為小于等于限制：
等于限制可以變為一個大于等于限制和一個小于等于限制，但在程式設計中一般支援直接寫等式限制,可以不進行轉換：
圖解法和單純形法
圖解法
對于較為簡單且隻有兩個決策變量的線性規劃問題可以使用圖解法。
考慮如下線性規劃模型：
Image Name
從圖中可以看出，當紅線(即目标函數)經過多邊形的頂點P(即表示兩個限制條件的直線交點)時， 軸截距取得最大值。
即：當 時，目标函數取得最大值為
單純形法
對于決策變量比較多的線性規劃模型，圖解法不再适用。
單純形法是1947 年G. B. Dantzig提出的一種十分有效的求解方法，極大地推廣了線性規劃的應用，直到今日也在一些線性規劃的求解器中使用。
從圖解法的例子中，我們可以看出，限制條件所圍成的區域為一個凸多邊形，當決策變量多于兩個時，限制條件圍成的區域為一個凸多面體，稱之為可行域。其中每一個面（稱之為超平面）即代表一個限制條件。
可以證明：線性規劃的最優解一定在可行域的邊界上
單純形法的思路就是在可行域的一個頂點處找到一個初始可行解，判斷該解是不是最優，若不是，則疊代到下一個頂點處進行重複判斷。因為最優解的搜尋範圍從整個可行域縮小到了可行域的有限個頂點，算法的效率得到了極大的提升。
具體的找初始可行解的方法，判斷解是否最優的條件，如何進行疊代這裡不做詳細展開，有興趣可以查閱相關資料
此外，求解線性規劃的方法還有橢球法、卡瑪卡算法、内點法等。
其中内點法因為求解效率更高，在決策變量多，限制多的情況下能取得更好的效果，目前主流線性規劃求解器都是使用的内點法。
使用python求解簡單線性規劃模型
程式設計思路
1. 選擇适當的決策變量
在解決實際問題時，把問題歸結成一個線性規劃數學模型是很重要的一步，但往往也是困難的一步，模型建立得是否恰當，直接影響到求解。而選适當的決策變量，是我們建立有效模型的關鍵之一。
2.将求解目标簡化為求一個目标函數的最大/最小值
能把要求解的問題簡化為一個最值問題是能否使用線性規劃模型的關鍵，如果這一點不能達到，之後的工作都有沒有意義的。
3. 根據實際要求寫出限制條件（正負性，資源限制等）
線性規劃的限制條件針對不同的問題有不同的形式，總結來說有以下三種：等式限制、不等式限制、符号限制
考慮如下線性規劃問題
Step1: 導入相關庫
import numpy as np
from scipy import optimize as op
Step2: 定義決策變量
# 給出變量取值範圍
x1=(0,None)
x2=(0,None)
x3=(0,None)
Step3: 将原問題化為标準形式
注意：程式設計時預設為最小化目标函數，是以這裡改為 ；第二個限制為大于等于限制，這裡化為小于等于限制；
Step4: 定義目标函數系數和限制條件系數
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數值
Step5: 求解
res=op.linprog(c,A_ub,B_ub,A_eq,B_eq,bounds=(x1,x2,x3)) #調用函數進行求解
res
con: array([0.])
fun: -14.571428571428571
message: 'Optimization terminated successfully.'
nit: 3
slack: array([0. , 3.85714286])
status: 0
success: True
x: array([6.42857143, 0.57142857, 0. ])
即當 時，目标函數取得最大值
求解案例
例1：使用scipy求解
#導入相關庫
import numpy as np
from scipy import optimize as op
#定義決策變量範圍
x1=(0,None)
x2=(0,None)
x3=(0,None)
#定義目标函數系數
c=np.array([2,3,1])
#定義限制條件系數
A_ub=np.array([[-1,-4,-2],[-3,-2,0]])
B_ub=np.array([-8,-6])
#求解
res=op.linprog(c,A_ub,B_ub,bounds=(x1,x2,x3))
res
con: array([], dtype=float64)
fun: 7.0
message: 'Optimization terminated successfully.'
nit: 2
slack: array([0., 0.])
status: 0
success: True
x: array([0.8, 1.8, 0. ])
即當 時，目标函數取得最小值
例2：包含非線性項的求解
由于存在非線性項，不能沿用例一中的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 # 等式限制
注意：每一個函數的輸入為一個 維列向量 ，其中 表示該列向量的第一個元素，即 。
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: 86
nit: 15
njev: 15
status: 0
success: True
x: array([0.55216734, 1.20325918, 0.94782404])
即當 時，目标函數取得最小值z = 10.65
從整數規劃到0-1規劃
整數規劃模型
規劃中的變量（部分或全部）限制為整數時，稱為整數規劃。若線上性規劃模型中，變量限制為整數，則稱為整數線性規劃。
當決策變量均為整數時，稱純整數規劃；
當決策變量中部分為整數，部分為實數時，稱混合整數規劃；
一般用 表示 為整數
将第一節中的線性規劃圖解法的例子添加整數限制，則可行域變為了多邊形内的整點，如下圖所示：
可以看出，可行域變成了離散的點，這也使得整數規劃問題比線性規劃問題要更難求解，但現實中的許多決策變量都隻能取整數，是以混合整數規劃問題也成為了了研究最多的線性規劃問題。
注意：整數規劃最優解不能按照實數最優解簡單取整而獲得
整數規劃的兩個常用求解方法：分支定界算法、割平面法
分枝定界法
step1不考慮整數限制的情況下求解得到最優解 （一般不是整數）；
step2以該解的上下整數界限建立新的限制，将原整數規劃問題變為兩個問題（分枝）；
step3分别對兩個子問題求解（不考慮整數限制），若解剛好為整數解則結束；若不為整數解則繼續進行分枝；
step4以最開始的目标函數值作為上界，子問題求解中得到的任一整數解為下界（定界），對子問題進行剪枝，減小問題規模；
step5重複以上步驟直到得到最優解
割平面法
step2通過該解做一個割平面（二維情況下為一條直線），縮小可行域；
step3在縮小後的可行域中求最優解（不考慮整數限制）
step4重複步驟2和步驟3，直到最優解滿足整數限制
0-1規劃模型
當整數規劃問題中的整數型決策變量限制為隻能取0或1時，稱為0-1整數規劃，簡稱為0-1規劃。
因為0-1規劃問題的解空間比一般的整數規劃問題較少，求解起來較為容易，且所有的整數規劃問題都可以化為0-1規劃問題，是以在建立混合整數規劃模型求解實際問題時，應盡量使用0-1決策變量進行模組化。
例如：有十個工廠可供決策時，可以使用10個0-1變量，當取值為0時時代表不使用這個工廠，取值為1時使用該工廠。
0-1規劃的常用求解方法：分支定界算法、割平面法、隐枚舉法
案例：投資的收益和風險
問題描述與分析
Image NameImage Name
符号規定
模型假設
建立與簡化模型
根據模型假設和符号規定，我們可以寫出模型的第一個優化目标為總體風險盡可能小,而總體風險是所有投資中風險最大的一個
第二個優化目标為淨收益盡可能大。根據題意，交易費用為一個分段函數（非線性函數），是以需要進行簡化：由于題目給定的定值 相對于總投資額 很小，可以忽略不計，是以将交易費簡化為 ,是以目标函數為
對于一個多目标優化模型，常用的考慮方式為先固定其中一個目标，再優化另一個目标。
在本題中，可以給定一個投資者能夠承受的風險界限 ,使得最大投資風險下損失比例小于 ,即 ，将其作為新的限制，就可以把多目标優化轉化為單目标優化，即：
使用python scipy庫求解
反映了投資者對風險的偏好程度，從 開始，以步長為0.001進行循環搜尋，使用python編寫代碼如下
#導入相關庫
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as op
#定義a的取值
a = 0
profit_list = [] #記錄最大收益
a_list = [] #記錄a的取值
while a<0.05:
#定義決策變量取值範圍
x1=(0,None)
#定義目标函數系數
c=np.array([-0.05,-0.27,-0.19,-0.185,-0.185])
#定義不等式限制條件左邊系數
A = np.hstack((np.zeros((4,1)),np.diag([0.025,0.015,0.055,0.026])))
#定義不等式限制條件右邊系數
b=a*np.ones((4,1));
#定義等式限制條件左邊系數
Aeq=np.array([[1,1.01,1.02,1.045,1.065]])
#定義等式限制條件右邊系數
beq=np.array([1]);
#求解
res=op.linprog(c,A,b,Aeq,beq,bounds=(x1,x1,x1,x1,x1))
profit = -res.fun
profit_list.append(profit)
a_list.append(a)
a = a+0.001
#繪制風險偏好a與最大收益的曲線圖
plt.figure(figsize=(10,7))
plt.plot(a_list,profit_list)
plt.xlabel('a');plt.ylabel('Profit')
Text(0, 0.5, 'Profit')
從上圖中可以看出
1.風險越大，收益也就越大；
2，當投資越分散時，投資者承擔的風險越小，這與題意一緻。即：
冒險的投資者會出現集中投資的情況，保守的投資者則盡量分散投資。
3，在 附近有一個轉折點，在這點左邊，風險增加很少時，利潤增長很快。在這一點右邊，風險增加很大時，利潤增長很緩慢，是以對于風險和收益沒有特殊偏好的投資者來說，應該選擇曲線的拐點作為最優投資組合，大約是 ，總體收益為 ，所對應投資方案為：
風險度 ，收益 ， ，，，
模型求解的其他思路
在上面的例子中，我們使用固定風險水準來最大化收益的方法來将多目标轉化為單目标，也可考慮其他思路：
- 在總盈利在水準 以上的情況下，尋找風險最低的投資方案，即：
2.對風險和收益賦予權重 和 , 成為投資偏好系數，即
作者：少年吉