初識隨機規劃：一個小小例子
初識隨機規劃：一個小小例子
本文中的確定性問題的例子來源於《運籌學》第四版，運籌學編寫組編著。
生產計劃的例子
一個工廠生產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
A、B的量，則該問題的數學模型為
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+2x2⩽84x1⩽164x2⩽12x1,x2⩾0
我們用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+a12x2⩽b1a21x1⩽b2a32x2⩽b3x1,x2⩾0
上述模型中，
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~12x2⩽b~1a~21x1⩽b~2a~32x2⩽b~3x1,x2⩾0
其中
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.13x2⩽83.73x1⩽163.90x2⩽12x1,x2⩾0
但是也有可能是這樣的
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.11x2⩽83.78x1⩽164.22x2⩽12x1,x2⩾0
總之，就是有無數種可能。我們無法窮舉。
隨機規劃模型（Stochastic Programming）
為了能夠處理這種情況，我們生成
K
K
K個非常非常具有代表性的場景(Scenarios), 假設每種場景出現的可能性為
p
k
,
k
∈
K
p_k, k \in K
pk,k∈K，且
∑
k
∈
K
q
k
=
1
\sum_{k \in K}q_k = 1
∑k∈Kqk=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.k∈K∑pk(c1x1k+c2x2k)a11kx1k+a12kx2k⩽8,∀k∈Ka21kx1k⩽16,∀k∈Ka22kx2⩽12,∀k∈Kx1k,x2k⩾0,∀k∈K
即，在每一種場景
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 )
max∑k∈Kpk(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.k∈K∑pki∈I∑cixiki∈I∑aijkxik⩽bj,∀j∈J,∀k∈Kxik⩾0,∀i∈I,∀k∈K
其中
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
相關文章
- 演算法系列-動態規劃(1)：初識動態規劃演算法動態規劃
- [Vuex系列] - 初嘗Vuex第一個例子Vue
- 動態規劃中初識狀態壓縮（入門）動態規劃
- 動態規劃初級動態規劃
- expdp一個例子
- 從零搭建一個IdentityServer——初識OpenIDConnectIDEServer
- 初級演算法-動態規劃演算法動態規劃
- 如何得到一個隨機密碼隨機密碼
- 初識dagger(一)
- 初識 webpack (一)Web
- 初識Django（一）Django
- 一、初識NettyNetty
- RocketMq（一）初識MQ
- 一.初識JavaJava
- 初識 Three.js：第一個小 demoJS
- js隨堂初體驗（一）JS
- Python正規表示式初識（四）Python
- UA MATH567 高維統計II 隨機向量6 亞高斯隨機向量的應用: 半正定規劃H5隨機
- 計算機相關知識的小小科普回顧計算機
- 【2018.04.19 ROS機器人作業系統】機器人控制：運動規劃、路徑規劃及軌跡規劃簡介之一ROS機器人作業系統
- 發個好地方規劃法規
- Flutter | 通過一個小例子帶你認識動畫 AnimationFlutter動畫
- 初識ABP vNext（1）：開篇計劃&基礎知識
- 一個例子看懂call,applyAPP
- 一個隨機數的類c++隨機C++
- 初識機器學習機器學習
- 敏捷規劃，讓你做一個有計劃的開發人敏捷
- Docker初認識（一）Docker
- kafka初認識（一）Kafka
- Pythonrandom模組（獲取隨機數）常用方法和使用例子Pythonrandom隨機
- 大規模知識庫中的隨機遊走推理和學習隨機
- 尤拉計劃695：隨機長方形隨機
- 為初學者介紹的 Linux tee 命令（6 個例子）Linux
- 三個規則劃分企業商機管理系統
- 一個小小工程師魔幻的 10年工程師
- 反射機制 小小談反射
- 初識 SpringMVC，執行配置第一個Spring MVC 程式SpringMVC
- 一個複雜的json例子JSON