初识随机规划：用一个小例子理解随机规划

最新推荐文章于 2024-03-31 21:36:26 发布

刘兴禄

最新推荐文章于 2024-03-31 21:36:26 发布

阅读量9.9k

点赞数 59

分类专栏：运筹优化理论与实践最优化求解器(Gurobi+CPLEX+SCIP)

本文链接：https://blog.csdn.net/HsinglukLiu/article/details/112150067

版权

运筹优化理论与实践同时被 2 个专栏收录

47 篇文章 197 订阅

订阅专栏

最优化求解器(Gurobi+CPLEX+SCIP)

39 篇文章 138 订阅

订阅专栏

文章目录

初识随机规划：一个小小例子
参考文献
OlittleRer

初识随机规划：一个小小例子

本文中的确定性问题的例子来源于《运筹学》第四版，运筹学编写组编著。

随机规划和鲁棒优化都是运筹优化中比较高阶的内容。二者都是考虑不确定情形下的数学规划，但是两者又有不同。随机规划旨在优化不确定情形下的目标函数的期望等，比如总收益的期望值等。但是鲁棒优化致力于使得最坏的情况最好。所以相比而言，基本的鲁棒优化获得的结果会比较保守。

本文主要介绍随机规划。随机规划是一种非常常用的用来处理不确定参数下的优化方法。本节不涉及非常深入的部分，仅以一个非常小的例子，来跟大家介绍什么是随机规划，以及随机规划能干什么。

生产计划的例子

一个工厂生产2种产品 $A$ 和 $B$ ， $A$ 产品需要1个单位人工工时和4个单位的设备1的工时， $B$ 产品需要2个单位的人工工时和4个单位的设备2的工时。且该厂一共可提供人工工时、设备1工时和设备2的工时分别为8,16和12。且生产单位产品 $A$ 、 $B$ 的净利润分别为2和3。假设该工厂生产的产品是供不应求的，问，该工厂应该如何生产(可以是连续的量)，使得净利润最大化。

根据上述描述，我们引入下面的决策变量：
$x_1, x_2$ 分别代表生产产品 $A 、 B$ 的量，则该问题的数学模型为

$\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}$
我们用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$ 产品和2单位 $B$ 产品，可以获得14个单位的利润。

这是一个非常简单的例子，下面我们基于这个例子，来引入什么是随机规划(Stochastic Programming)。

参数的不确定性

我们将上述模型抽象称为一个参数形式的模型，其中 $c$ , $a$ , $b$ 均为参数。
$\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}$
上述模型中， $c_i$ 即为产品 $i$ 的净利润， $a_{ij}$ 即为资源消耗参数， $b_i$ 即为可用资源的数量。在确定性的问题中，我们假设所有参数，即 $c, a, b$ 都是可以提前给定的。比如商品的售价、可用的资源以及生产的时间等。

但是在实际的场景中，这些参数也许是不确定的，比如生产所需要的工时，可能会由于机器、人工的误差，导致一些波动(即 $a$ 出现波动)。一些产品的价格，也许会随着季节，供需关系的变动而变动(即 $c$ 出现波动)。生产资料也有可能出现同样的不确定性(即 $b$ 出现波动)。

如果我们按照确定性情形下的参数去做计划(即 $c_1=2, c_2=3, a_{11}=1, a_{12}=2, \cdots$ )，在实际生产活动中，由于随机的因素，可能导致真正的参数跟确定情形下考虑的参数不一致。比如员工们的状态不好，导致人工生产时间变长了一些(比如 $a_{11}$ 变成了1.05等)，或者设备启动等除了点小意外，使得需要的机器工时多了一些等。在这种情况下，我们原来考虑确定性参数的生产计划，在随机的情形下也许就会不可行，或者说受到较大的影响。

一种可行的解决方案就是：随机规划。

体现在上面的例子中，就是我们考虑模型的参数是不确定的。我们用下面的数学语言来描述
$\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}$
其中 $\tilde{c}, \tilde{a}, \tilde{b}$ 表示这些参数是不确定的。但是这些参数又不是完全没有已知信息。我们假设他们的一些概率分布信息是已知的。比如均值、方差等。

例如，我们假设已知
$\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}$

但是有了这些分布信息，我们如何去处理呢？一个可选的方法就是，根据这些分布信息，生成若干个具有代表性的场景(Scenarios)，每一个场景就对应一组参数，一个模型，并且这个场景有其对应的概率。我们可以用这些场景，去代表现实情况的所有可能。基于这些场景，我们可以给出一个综合考虑了随机因素的生产计划。这种方法也叫基于场景的建模方法(Scenario-based)。至于如何使得生成的场景比较具有代表性，那就又是一大部分内容了，这里我们不做具体展开。(不过可以先提一下，一般比较常用的就是Sample Average Approximation, SAA)

为了方便，我们只考虑 $a$ 不确定， $c$ 和 $b$ 的不确定在这里先不做考虑。
首先，我们假设确定性模型中的参数 $a$ 的取值就等于每个参数的均值，且所有参数的方差均为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}$
比如，根据这个信息，在实际的情况中，参数可能是这样的
$\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}$
但是也有可能是这样的
$\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}$
总之，就是有无数种可能。我们无法穷举。

随机规划模型（Stochastic Programming）

为了能够处理这种情况，我们生成 $K$ 个非常非常具有代表性的场景(Scenarios), 假设每种场景出现的可能性为 $p_k, k \in K$ ，且 $\sum_{k \in K}q_k = 1$ 。则生产计划的随机规划模型就可以写成
$\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}$
即，在每一种场景 $k$ 下，我们都需要决策出那种场景 $k$ 下的生产计划 $x_1^k, x_2^k$ ，我们最终的计划，要使得我们考虑的 $K$ 个场景下的总期望收益最大化，即 $\max \, \sum_{k \in K}p_k (c_1 x_1^k + c_2 x_2^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}$
其中 $I = \{1, 2\}$ 表示产品的集合， $J = \{1, 2, 3\}$ 表示资源的集合。

这就是非常常见的随机规划模型。一些相关论文的模型本质上跟上面的形式是类似的。

Python调用Gurobi求解随机规划模型

接下来我们用Python调用Gurobi对上述随机规划进行求解。代码中的Scenarios不是用SAA做的，而是我根据分布信息随机生成的，只是为了展示一下这个过程。

我们考虑5种场景，即 $K = 5$ ，且 $p_i = 0.2, i = 1, 2, \cdots, 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

可以看到，随机规划实际上是给出了一系列的解决方案。在每一种场景下，都会有一组解。

可以这么去理解，一个Scenario就是一种实际情况的可能。比如说，今天做好了生产计划，可是到了明天，设备1真的出问题了，其加工时长变长了，那就对应一个场景。如果加工时长没变，但是有员工1状态非常好，速度比平均状态快了，那就对应另外一种场景……所有的这些可能出现的情况，都可以用我们生成的这5种Scenario去代表，也就是说，到了第二天，我们发现真实情况是情况1，我们拿情况1与我们考虑的5种Scenario对比，发现情况1和Scenario 2非常相似，那么此时，Scenario 2就可以代表情况1。同样的，其他情况也可以通过这种方式去对应。也就是说，5种Scenario就近似的代表了所有可能发生的情况，并且我们的随机规划模型提供了每种Scenario下的相应的解，只需要对应去执行就可以了。
而我们的目标函数，是这些所有情况下的总收益的期望值最大。即，不管明天发生什么样的随机情况，我们期望能获得的收益就是那么多。至于具体获得多少，那就看具体实际情况是哪种Scenario了。

总结一下，随机规划实际上是给出了你考虑所有情形下的解。这些解的平均目标函数是达到了最大。在实际中，发生了哪一种对应的情况，我们就去执行那种场景下对应的解即可。

上面的例子给大家直观的展示了什么是随机规划。想要继续深入的读者，可以自行阅读相关书籍(参考文献中给出了几个参考书籍)。

参考文献

[1]: Shapiro, A and Ruszczynski, A, Handbooks in operations research and management science: Vol. 10. Stochastic programming, 2003
[2]: Birge, John R., and Francois Louveaux. Introduction to stochastic programming. Springer Science & Business Media, 2011.

作者：刘兴禄，清华大学，清华伯克利深圳学院，博士在读
联系方式：hsinglul@163.com
博客主页：https://blog.csdn.net/HsinglukLiu/article/details/107848461