初識隨機規劃:一個小小例子

劉興祿發表於2021-01-03

初識隨機規劃:一個小小例子

本文中的確定性問題的例子來源於《運籌學》第四版,運籌學編寫組編著。

生產計劃的例子

一個工廠生產2種產品 A A A B B B A A A產品需要1個單位人工工時和4個單位的裝置1的工時, B B B產品需要2個單位的人工工時和4個單位的裝置2的工時。且該廠一共可提供人工工時、裝置1工時和裝置2的工時分別為8,16和12。且生產單位產品 A A A B B B的淨利潤分別為2和3。假設該工廠生產的產品是供不應求的 ,問,該工廠應該如何生產(可以是連續的量),使得淨利潤最大化。

根據上述描述,我們引入下面的決策變數:
x 1 , x 2 x_1, x_2 x1,x2分別代表生產產品 A 、 B A、B AB的量,則該問題的數學模型為

max ⁡    2 x 1 + 3 x 2 s . t . x 1 + 2 x 2 ⩽ 8 4 x 1 ⩽ 16 4 x 2 ⩽ 12 x 1 , x 2 ⩾ 0 \begin{aligned} \max \,\,&2x_1 + 3x_2 \\ s.t. \quad &x_1 + 2x_2 \leqslant 8 \\ &4x_1 \qquad \leqslant 16 \\ &\qquad 4x_2 \leqslant 12 \\ &x_1, x_2 \geqslant 0 \end{aligned} maxs.t.2x1+3x2x1+2x284x1164x212x1,x20
我們用Pytho呼叫Gurobi對其進行求解,具體程式碼如下:

from gurobipy import * 
model = Model('production plan')
x = {} 
for i in range(1,3):
    x[i] = model.addVar(lb = 0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = 'x_' + str(i))

model.setObjective(2 * x[1] + 3 * x[2], GRB.MAXIMIZE)
model.addConstr(x[1] + 2 * x[2] <= 8, name = 'cons_1')
model.addConstr(4 * x[1] <= 16, name = 'cons_2')
model.addConstr(4 * x[2] <= 12, name = 'cons_2')

# solve 
model.optimize()

# print the results
print('Obj = ', model.ObjVal)
for key in x.keys():
    print(x[key].varName, ' = ', x[key].x) 

求解結果如下

Solved in 1 iterations and 0.02 seconds
Optimal objective  1.400000000e+01
Obj: 14.0
x_1  =  4.0
x_2  =  2.0

即生產4單位 A A A產品和2單位 B B B產品,可以獲得14個單位的利潤。

這是一個非常簡單的例子,下面我們基於這個例子,來引入什麼是隨機規劃(Stochastic Programming)。

引數的不確定性

我們將上述模型抽象稱為一個引數形式的模型,其中 c c c, a a a, b b b均為引數。
max ⁡    c 1 x 1 + c 2 x 2 s . t . a 11 x 1 + a 12 x 2 ⩽ b 1 a 21 x 1         ⩽ b 2            a 32 x 2 ⩽ b 3 x 1 , x 2 ⩾ 0 \begin{aligned} \max \,\,&c_1x_1 + c_2x_2 \\ s.t. \quad &a_{11}x_1 + a_{12}x_2 \leqslant b_1 \\ &a_{21}x_1 \qquad\,\,\,\,\,\,\, \leqslant b_2 \\ &\,\,\,\,\,\,\,\,\,\,\qquad a_{32}x_2 \leqslant b_3 \\ &x_1, x_2 \geqslant 0 \end{aligned} maxs.t.c1x1+c2x2a11x1+a12x2b1a21x1b2a32x2b3x1,x20
上述模型中, c i c_i ci即為產品 i i i的淨利潤, a i j a_{ij} aij即為資源消耗引數, b i b_i bi即為可用資源的數量。在確定性的問題中,我們假設所有引數,即 c , a , b c,a,b c,a,b都是可以提前給定的。比如商品的售價、可用的資源以及生產的時間等。

但是在實際的場景中,這些引數也許是不確定的,比如生產所需要的工時,可能會由於機器、人工的誤差,導致一些波動(即 a a a出現波動)。一些產品的價格,也許會隨著季節,供需關係的變動而變動(即 c c c出現波動)。生產資料也有可能出現同樣的不確定性(即 b b b出現波動)。

如果我們按照確定性情形下的引數去做計劃(即 c 1 = 2 , c 2 = 3 , a 11 = 1 , a 12 = 2 , ⋯ c_1=2, c_2=3, a_{11}=1, a_{12}=2, \cdots c1=2,c2=3,a11=1,a12=2,),在實際生產活動中,由於隨機的因素,可能導致真正的引數跟確定情形下考慮的引數不一致。比如員工們的狀態不好,導致人工生產時間變長了一些(比如 a 11 a_{11} a11變成了1.05等),或者裝置啟動等除了點小意外,使得需要的機器工時多了一些等。在這種情況下,我們原來考慮確定性引數的生產計劃,在隨機的情形下也許就會不可行,或者說受到較大的影響。

一種可行的解決方案就是:隨機規劃。

體現在上面的例子中,就是我們考慮模型的引數是不確定的。我們用下面的數學語言來描述
max ⁡    c ~ 1 x 1 + c ~ 2 x 2 s . t . a ~ 11 x 1 + a ~ 12 x 2 ⩽ b ~ 1 a ~ 21 x 1         ⩽ b ~ 2            a ~ 32 x 2 ⩽ b ~ 3 x 1 , x 2 ⩾ 0 \begin{aligned} \max \,\,&\tilde{c}_1x_1 + \tilde{c}_2x_2 \\ s.t. \quad &\tilde{a}_{11}x_1 + \tilde{a}_{12}x_2 \leqslant \tilde{b}_1 \\ &\tilde{a}_{21}x_1 \qquad\,\,\,\,\,\,\, \leqslant \tilde{b}_2 \\ &\,\,\,\,\,\,\,\,\,\,\qquad \tilde{a}_{32}x_2 \leqslant \tilde{b}_3 \\ &x_1, x_2 \geqslant 0 \end{aligned} maxs.t.c~1x1+c~2x2a~11x1+a~12x2b~1a~21x1b~2a~32x2b~3x1,x20
其中 c ~ , a ~ , b ~ \tilde{c}, \tilde{a}, \tilde{b} c~,a~,b~表示這些引數是不確定的。但是這些引數又不是完全沒有已知資訊。我們假設他們的一些概率分佈資訊是已知的。比如均值、方差等。

例如,我們假設已知
E ( a ~ i j ) = μ i j , σ i j = 1 E ( b ~ i ) = μ i , σ i = 1 \begin{aligned} \mathbb{E}(\tilde{a}_{ij}) = \mu_{ij}, \quad \sigma_{ij} = 1 \\ \mathbb{E}(\tilde{b}_{i}) = \mu_{i}, \quad \sigma_{i} = 1 \end{aligned} E(a~ij)=μij,σij=1E(b~i)=μi,σi=1

但是有了這些分佈資訊,我們如何去處理呢?一個可選的方法就是,根據這些分佈資訊,生成若干個具有代表性的場景(Scenarios),每一個場景就對應一組引數,一個模型,並且這個場景有其對應的概率。我們可以用這些場景,去代表現實情況的所有可能。基於這些場景,我們可以給出一個綜合考慮了隨機因素的生產計劃。這種方法也叫基於場景的建模方法(Scenario-based)。

為了方便,我們只考慮 a a a不確定, c c c b b b的不確定在這裡先不做考慮。
首先,我們假設確定性模型中的引數 a a a的取值就等於每個引數的均值,且所有引數的方差均為0.2,且服從正態分佈,即
E ( a ~ 11 ) = μ 11 = 1 , σ 11 = 0.2 E ( a ~ 12 ) = μ 12 = 2 , σ 12 = 0.2 E ( a ~ 21 ) = μ 21 = 4 , σ 21 = 0.2 E ( a ~ 32 ) = μ 32 = 4 , σ 32 = 0.2 \begin{aligned} \mathbb{E}(\tilde{a}_{11}) = \mu_{11} = 1, \quad \sigma_{11} = 0.2 \\ \mathbb{E}(\tilde{a}_{12}) = \mu_{12} = 2, \quad \sigma_{12} = 0.2 \\ \mathbb{E}(\tilde{a}_{21}) = \mu_{21} = 4, \quad \sigma_{21} = 0.2 \\ \mathbb{E}(\tilde{a}_{32}) = \mu_{32} = 4, \quad \sigma_{32} = 0.2 \end{aligned} E(a~11)=μ11=1,σ11=0.2E(a~12)=μ12=2,σ12=0.2E(a~21)=μ21=4,σ21=0.2E(a~32)=μ32=4,σ32=0.2
比如,根據這個資訊,在實際的情況中,引數可能是這樣的
max ⁡    2 x 1 + 3 x 2 s . t . 1.10 x 1 + 2.13 x 2 ⩽ 8 3.73 x 1         ⩽ 16            3.90 x 2 ⩽ 12 x 1 , x 2 ⩾ 0 \begin{aligned} \max \,\,&2x_1 + 3x_2 \\ s.t. \quad &1.10x_1 + 2.13x_2 \leqslant 8 \\ &3.73x_1 \qquad\,\,\,\,\,\,\, \leqslant 16 \\ &\,\,\,\,\,\,\,\,\,\,\qquad 3.90x_2 \leqslant 12 \\ &x_1, x_2 \geqslant 0 \end{aligned} maxs.t.2x1+3x21.10x1+2.13x283.73x1163.90x212x1,x20
但是也有可能是這樣的
max ⁡    2 x 1 + 3 x 2 s . t . 0.95 x 1 + 2.11 x 2 ⩽ 8 3.78 x 1         ⩽ 16            4.22 x 2 ⩽ 12 x 1 , x 2 ⩾ 0 \begin{aligned} \max \,\,&2x_1 + 3x_2 \\ s.t. \quad &0.95x_1 + 2.11x_2 \leqslant 8 \\ &3.78x_1 \qquad\,\,\,\,\,\,\, \leqslant 16 \\ &\,\,\,\,\,\,\,\,\,\,\qquad 4.22x_2 \leqslant 12 \\ &x_1, x_2 \geqslant 0 \end{aligned} maxs.t.2x1+3x20.95x1+2.11x283.78x1164.22x212x1,x20
總之,就是有無數種可能。我們無法窮舉。

隨機規劃模型(Stochastic Programming)

為了能夠處理這種情況,我們生成 K K K個非常非常具有代表性的場景(Scenarios), 假設每種場景出現的可能性為 p k , k ∈ K p_k, k \in K pk,kK,且 ∑ k ∈ K q k = 1 \sum_{k \in K}q_k = 1 kKqk=1。則生產計劃的隨機規劃模型就可以寫成
max ⁡    ∑ k ∈ K p k ( c 1 x 1 k + c 2 x 2 k ) s . t . a 11 k x 1 k + a 12 k x 2 k ⩽ 8 , ∀ k ∈ K a 21 k x 1 k         ⩽ 16 , ∀ k ∈ K            a 22 k x 2 ⩽ 12 , ∀ k ∈ K x 1 k , x 2 k ⩾ 0 , ∀ k ∈ K \begin{aligned} \max \,\,& \sum_{k \in K}p_k (c_1 x_1^k + c_2 x_2^k ) \\ s.t. \quad & a_{11}^k x_1^k + a_{12}^kx_2^k \leqslant 8, \quad \forall k \in K \\ &a_{21}^kx_1^k \qquad\,\,\,\,\,\,\, \leqslant 16, \quad \forall k \in K \\ &\,\,\,\,\,\,\,\,\,\,\qquad a_{22}^kx_2 \leqslant 12, \quad \forall k \in K \\ &x_1^k, x_2^k \geqslant 0, \quad \forall k \in K \end{aligned} maxs.t.kKpk(c1x1k+c2x2k)a11kx1k+a12kx2k8,kKa21kx1k16,kKa22kx212,kKx1k,x2k0,kK
即,在每一種場景 k k k下,我們都需要決策出那種場景 k k k下的生產計劃 x 1 k , x 2 k x_1^k, x_2^k x1k,x2k,我們最終的計劃,要使得我們考慮的 K K K個場景下的總期望收益最大化,即 max ⁡   ∑ k ∈ K p k ( c 1 x 1 k + c 2 x 2 k ) \max \, \sum_{k \in K}p_k (c_1 x_1^k + c_2 x_2^k ) maxkKpk(c1x1k+c2x2k)

我們將上述模型再寫地緊湊一些,即為
max ⁡    ∑ k ∈ K p k ∑ i ∈ I c i x i k s . t . ∑ i ∈ I a i j k x i k ⩽ b j , ∀ j ∈ J , ∀ k ∈ K x i k ⩾ 0 , ∀ i ∈ I , ∀ k ∈ K \begin{aligned} \max \,\,& \sum_{k \in K}p_k \sum_{i \in I}c_i x_i^k \\ s.t. &\quad \sum_{i \in I}{a_{ij}^k x_i^k} \leqslant b_j, \forall j \in J, \forall k \in K \\ &x_i^k \geqslant 0, \forall i \in I, \forall k \in K \end{aligned} maxs.t.kKpkiIcixikiIaijkxikbj,jJ,kKxik0,iI,kK
其中 I = { 1 , 2 } I = \{1, 2\} I={1,2}表示產品的集合, J = { 1 , 2 , 3 } J = \{1, 2, 3\} J={1,2,3}表示資源的集合。

這就是非常常見的隨機規劃模型。一些相關論文的模型本質上跟上面的形式是類似的。

Python呼叫Gurobi求解隨機規劃模型

接下來我們用Python呼叫Gurobi對上述隨機規劃進行求解。

我們考慮5種場景,即 K = 5 K=5 K=5,且 p i = 0.2 , i = 1 , 2 , ⋯   , 5 p_i = 0.2, i = 1, 2, \cdots, 5 pi=0.2,i=1,2,,5
結果如下

scenario_num = 5 
sto_model = Model('Stochastic Programming')
prob = [0.2, 0.2, 0.2, 0.2, 0.2] 
# prob = [0.1, 0.2, 0.3, 0.15, 0.25]
profit = [2, 3] 
b = [8, 16, 12] 

# generate parameters under scenarios 
a_11 = np.random.normal(1, 0.2, scenario_num)  # mu, sigma, sample_number 
a_12 = np.random.normal(2, 0.2, scenario_num)
a_21 = np.random.normal(4, 0.2, scenario_num)
a_32 = np.random.normal(4, 0.2, scenario_num) 
    
x_sto = {}
for i in range(2): 
    # creat variables 
    for s in range(scenario_num):
        x_sto[i, s] = sto_model.addVar(lb = 0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = 'x_' + str(i) + '_' + str(s))
    
# set objective 
obj = LinExpr(0)
for key in x_sto.keys():
    product_ID = key[0]
    scenario_ID = key[1] 
    obj.addTerms(prob[scenario_ID] * profit[product_ID], x_sto[key])
sto_model.setObjective(obj, GRB.MAXIMIZE)

# add constraints:
# constraints 1 
for s in range(scenario_num):
    lhs = LinExpr(0)
    lhs.addTerms(a_11[s], x_sto[0, s]) 
    lhs.addTerms(a_12[s], x_sto[1, s])  
    sto_model.addConstr(lhs <= b[0], name = 'cons_' + str(0))  
    
# constraints 2      
for s in range(scenario_num):
    lhs = LinExpr(0)
    lhs.addTerms(a_21[s], x_sto[0, s]) 
    sto_model.addConstr(lhs <= b[1], name = 'cons_' + str(1))   
    
# constraints 3      
for s in range(scenario_num):
    lhs = LinExpr(0)
    lhs.addTerms(a_32[s], x_sto[1, s])  
    sto_model.addConstr(lhs <= b[2], name = 'cons_' + str(2))     

# solve 
sto_model.optimize()

print('Obj = ', sto_model.ObjVal)
# for key in x_sto.keys():
#     print(x_sto[key].varName, ' = ', x_sto[key].x)   
for s in range(scenario_num): 
    print(' ------  scenario   ', s, '  ------- ') 
    for i in range(2):
        print(x_sto[i,s].varName, ' = ', x_sto[i,s].x)  
    print('\n\n') 

求解結果如下

Solved in 4 iterations and 0.01 seconds
Optimal objective  1.488085782e+01
Obj =  14.880857821250329
 ------  scenario    0   ------- 
x_0_0  =  3.8650223745517267
x_1_0  =  2.4927517036690205



 ------  scenario    1   ------- 
x_0_1  =  3.999245964351986
x_1_1  =  3.0288732634647424



 ------  scenario    2   ------- 
x_0_2  =  1.706111126833073
x_1_2  =  3.0604202728741234



 ------  scenario    3   ------- 
x_0_3  =  3.9331434475778084
x_1_3  =  2.4207092574101026



 ------  scenario    4   ------- 
x_0_4  =  4.003138883536607
x_1_4  =  2.127567340098422

改變場景的出現概率,比如考慮

prob = [0.1, 0.2, 0.3, 0.15, 0.25]

則求解結果如下:

Solved in 5 iterations and 0.01 seconds
Optimal objective  1.406161189e+01
Obj =  14.061611891889607
 ------  scenario    0   ------- 
x_0_0  =  4.371132256077677
x_1_0  =  1.6122311882908937



 ------  scenario    1   ------- 
x_0_1  =  3.9240152749431267
x_1_1  =  2.280958241360689



 ------  scenario    2   ------- 
x_0_2  =  2.1745705126924615
x_1_2  =  3.2383131316117613



 ------  scenario    3   ------- 
x_0_3  =  4.215519897219448
x_1_3  =  2.0395998835922193



 ------  scenario    4   ------- 
x_0_4  =  3.992288784058322
x_1_4  =  1.82358745935411

可以看到,隨機規劃時給出了一系列的解決方案。在每一種場景下,都會有一組解。

也就是說,隨機規劃實際上是給出了你考慮所有情形下的解。這些解的平均目標函式是達到了最大。在實際中,發生了哪一種對應的情況,我們就去執行那種場景下對應的解即可。

上面的例子給大家直觀的展示了什麼是隨機規劃。想要繼續深入的讀者,可以自行閱讀相關書籍。


作者:劉興祿,清華大學,清華伯克利深圳學院,博士在讀
聯絡方式:hsinglul@163.com
部落格主頁:https://blog.csdn.net/HsinglukLiu/article/details/107848461

相關文章