随机规划：求解报童问题期望值模型的算法方案

1 引言

上一篇关于不确定优化的文章（不确定优化入门：用简单实例讲明白随机规划、鲁棒优化和分布鲁棒优化 ）发表后，被部分大佬认为是一篇科普文，有点诚惶诚恐，毕竟我从今年才开始认真学习不确定优化，水平还很有限。

秉承着先了解大概框架再深入学习各类算法细节的心态，我在那篇文章中，用简化版的报童问题，把求解不确定优化的各类算法方案做了直观上的实践和比较，完成了算法框架的初步梳理。

现在，已经到学习算法方案详细设计的阶段了。

我的学习路径一如既往的是从简单到复杂，所以，本文要研究的问题是随机规划中最简单的一类：决策变量只有1个，不确定参数也只有1个且分布函数已知；解决问题的基本思路是最优化其目标函数的期望值，有些书中也称之为期望值模型。

正文见下。

2 数学模型

先看一下确定性优化模型
$$\begin{gathered} \text{min}f(x) \\ \text{s.t}{g}_j(x)\le 0,j=1,2,...,N \end{gathered}$$
其中$x$是决策变量，$f$是目标函数表达式，$g$是约束条件表达式。

增加随机变量$\xi$后，上述模型就变为不确定优化模型
$$\begin{gathered} \text{min}E[f(x,\xi )] \\ \text{s.t}E[g_j(x,\xi )]\le 0,j=1,2,...,N \end{gathered}$$
式中，$E$表示期望，假设$\xi$的概率密度函数为$\varphi (\xi )$，目标函数和约束条件的表达式可以表示为
$$\begin{gathered} E[f(x,\xi )]=\int f(x,\xi )\varphi (\xi )d\xi \\ E[g_j(x,\xi )]=\int g_j(x,\xi )\varphi (\xi )d\xi \end{gathered}$$
如果$\xi$是离散随机变量，且其分布函数为$P(\xi ={\xi }_i)={\theta }_i,i\in I$，目标函数和约束条件的表达式可以表示为
$$\begin{gathered} E[f(x,\xi )]=\sum _{i\in I}{\theta }_{i}f(x,{\xi }_i) \\ E[g_j(x,\xi )]=\sum _{i\in I}{\theta }_i{g}_j(x,{\xi }_i) \end{gathered}$$

3 报童问题

有了模型后，我们还需要知道用什么算法可以求解。

本节主要使用报童问题作为实例，来阐述在决策变量只有1个、不确定参数也只有1个的情况下，可以得到最优解的各类算法方案。

报童问题可以描述为：报童每天需要采购一定数量的报纸用于当天的销售。已知每份报纸的成本价$c=5$，销售价$p=8$，需求量$d$是个不确定参数，通过历史的数据可知其分布服从正态分布，均值是$\mu =100$，方差为$\sigma =20$，如果当天卖不完，会按回收价$s=4$将未卖完的报纸卖给回收站。

现在需要确定报童的最佳订购量$x$，使得报童的净收入$\theta (x)$最大化。$\theta (x)$的表达式为
$$\theta (x)=p\cdot E[\min (x,d)]+s\cdot E[\max (x-d,0)]-cx$$
第一项是售卖报纸的收益，第二项是回收报纸的收益，第三项是购买报纸的成本。

3.1 直接最优化

常规求解思路是：对$\theta (x)$求导使其梯度等于0，即可得到最佳$x$。

但由于公式中涉及到$x$和$d$的大小判断，$\theta (x)$的导数并不好算，为此，可以先把上面的$x$转换成一个新的变量$x-d$，接着转变$\theta (x)$的表达式。

1. 对于$\min (x,d)$项，将其做如下转化
   $$\min (x,d)=d-\max (d-x,0)$$
   上式的正确性可以分$x>d$和$x\le d$来依次验证：$x>d$时，左边为$d$，右边为$d-0=d$；$x\le d$时，左边为$x$，右边为$d-(d-x)=x$。

2. 对于$s\cdot E[\max (x-d,0)]$项，不需要调整。

3. 对于$cx$项，将其做如下转化
   $$cx=c(x-d)+cd=-c\max (d-x,0)+c\max (x-d,0)+cd$$
   上式的正确性也可以分$x>d$和$x\le d$来验证，这里就不赘述了。

4. 将上述三式重新组合，目标函数$\theta (x)$的表达式变为
   $$\theta =(p-c)\cdot E(d)-\{(p-c)\cdot E[\max (d-x,0)]+(c-s)\cdot E[\max (x-d,0)]\}$$
   定义$b=p-c,h=c-s,\alpha =b\cdot E[\max (d-x,0)]+h\cdot E[\max (x-d,0)]$，则最大化$\theta$等价于最小化$\alpha$。

事实上，$\alpha$中的两个表达式也可以分别理解为：因购买量过少导致的脱销损失和因购买量过多导致的滞销损失。此外，从上式中还可以看出，由于$\alpha \ge 0$，所以考虑不确定后，期望收益肯定不高于不考虑不确定时的收益值。

为了求得最小化的$\alpha$，将其表达式求导并令其等于0
$$b\cdot E\frac{\partial \max (d-x,0)}{\partial x}+h\cdot E\frac{\partial \max (x-d,0)}{\partial x}=0$$
需要注意的是，一维求导一般用$dx$就可以，但是本文已经使用$d$表示需求量，为了区分，此处使用了$\partial x$的形式。

先看第一项。$x\ge d$时，$\max (d-x,0)=0$，对$x$求导后是0；当$x<d$时，$\max (d-x,0)=d-x$，对$x$求导后是-1。设$x$的分布概率为$P(x)$，第一项可以表示为
$$b\cdot E\frac{\partial \max (d-x,0)}{\partial x}=b{\int }_{x_{\min }}^{d}P(x)\cdot (-1)dx=-bPr(x<d)=-b[1-Pr(d\le x)]$$
式中，$Pr(d\le x)$为累积分布函数。

同理，可以求得第二项为的结果为
$$h\cdot E\frac{\partial \max (x-d,0)}{\partial x}=hPr(d\le x)$$
将上述两式带回梯度值等于0的等式约束中，可以得到
$$Pr(d\le x)=\frac{b}{b+h}=\frac{p-c}{p-s}$$

这里还有个小的点需要注意，我们要计算的是$Pr(d\le x)$，而不是$Pr(x\le d)$。虽然$x$是决策变量，但是随机变量是$d$，后续还需要根据$d$的累积分布函数值反求$x$。

根据之前的定义，$b=1,h=3$，此时最优解满足
$$Pr(d\le x)=\frac{b}{b+h}=0.75$$

调用如下代码，可以反算出$x=113.49$。

IN [40]: import scipy
IN [41]: scipy.stats.norm.ppf(0.75, loc=100, scale=20)
Out[41]: 113.48979500392163

3.2 样本均值近似

用最优化的方法推导最优解的解析表达式，虽然优雅，但是对数学的要求比较高，而且如果问题复杂度提升了，能否推导出来都是一个问题。

既然如此，数值的方法就也值得一试。其中，最常见的方法，就是样本均值近似。在该方法中，通过抽样的方式把随机变量转变为一组离散参数，这样就可以把不确定优化问题转化为确定优化问题。

以下为使用样本均值近似实现不确定优化的具体代码。首先使用正态分布随机产生100000个$d$值，然后设置决策变量范围为80~120，依次计算不同决策变量值下的总收益，并将最优解保留下来。

import numpy as np

if __name__ == '__main__':

    # 报童模型参数
    c = 5
    s = 4
    p = 8

    # 需求分布参数
    mu = 100
    sigma = 20
    d = np.random.normal(mu, sigma, 100000)

    # 最优解
    best_f = 0
    best_x = 0

    # 遍历所有决策变量x：范围80~120
    for cur_x in range(80, 120):

        cur_f = 0
        # 需求参数
        for cur_d_index in range(len(d)):
            cur_d = d[cur_d_index]
            # 计算当前决策变量和当前需求值时的值
            cur_f += p * min(cur_x, cur_d) + s * max(cur_x - cur_d, 0) - c * cur_x

        # 更新最优解
        if cur_f > best_f:
            best_f = cur_f
            best_x = cur_x

    print('best_x: {}, best_f: {}.'.format(best_x, best_f / len(d)))

运行上述代码后，可以得到最佳订购量是113，该值和解析推导的结果是吻合的。

best_x: 113, best_f: 274.4170922340697.

3.3 两阶段规划

相比直接求解原模型，两阶段规划是更常见的求解随机规划问题的方法。该方法在建模时，会把原问题划分为两个阶段，通常第一阶段不包含随机变量，第二阶段包含随机变量。

在报童问题中，第一阶段是：在观测到实际需求$d$之前，就需要确定订货量 $x$，目标函数表达式为
$$-cx+E[Q(x,d)]$$
此处，$E[Q(x,d)]$是第二阶段的收益，目前未知，暂时用一个表达式代替。

到了第二阶段，需求量变为已知，此时$x$也已经确定，$E[Q(x,d)]$的表达式就明确了
$$E[Q(x,d)]=p\cdot E[\min (x,d)]+s\cdot E[\max (x-d,0)]$$

把上式带入上上式后，我们发现，和直接建模的目标函数表达式是一模一样的，看起来好像没啥区别。这主要是因为本文所用实例的问题维度比较低，如果问题复杂度变高，两阶段规划模型和原模型就会不一样。后续文章肯定还会再讨论该方案，此处暂时就不再深入了。

3.4 结果分析

前面三节给出了求解期望值模型的三种算法方案，本节针对优化结果再做一些简单分析。

如果不考虑随机性，最优订购量是100，收益是300；考虑随机性后，最优订购量是113，预期收益是274。

预期收益降低，在“解析推导”章节中已经说明了原因，此处不再赘述；本节着重理解为什么考虑随机性后最优订购量增会高于100。

事实上，当$x<d$时，脱销损失的系数值$h=3$，而$x>d$时，滞销损失的系数值$b=1$。由于目标函数是最小化脱销损失+滞销损失，所以最优解会高于需求量均值$d=100$，可以理解为：宁愿滞销、不能脱销，即宁愿多赔，不能少赚。

如果调整模型参数为$p=6,s=3$，此时滞销损失系数值高于脱销损失系数值，最优订购量小于100，变成91，即宁愿少赚，不能多赔。

4 在线求教

关于文章主题的内容，到这里已经结束了。本节主要是针对我在学习不确定优化过程中，遇到的一些困扰，进行在线求教，有愿意帮忙的大佬可以直接后台联系我。

(1) 需要帮忙勘误。国内关于不确定优化的经典教材或视频教程偏少，我自己也是边学边总结，很担心文章中的一些认知和理解有问题，会误导大家，所以迫切需要这个方向的大佬在我发文前帮忙做一下免费勘误。

(2) 需要RSOME指导。之前有个RSOME的讲座，主要讲的是如何利用RSOME求解鲁棒优化和分布鲁棒优化问题，当时主讲人的意思是RSOME也能用于求解随机规划问题。我目前还没用过RSOME，但很想实践一下并将其优化结果和代码添加至文章中，所以需要有擅长RSOME的大佬帮忙指导一下内容或方向。

在此感谢！

5 相关阅读

报童问题：https://www.jianshu.com/p/dfac740faa34

报童模型两阶段表示：https://blog.csdn.net/robert_chen1988/article/details/79350624

报童问题详细推导及利用Python的SAA方法求解：https://blog.csdn.net/zzzssszzzzsz/article/details/119875340

报童模型——如何做一次性采购的供给决策？：https://www.jianshu.com/p/6c7858159898