随机规划场景中的两类目标利润概率模型

1 引言

上个月因回老家结婚和过春节，持续学习的计划暂时被搁置一旁。这个月开始，要重回正轨了。

去年，我总结了几篇关于随机规划的文章：①随机规划：求解报童问题期望值模型的算法方案；②一文了解经典报童模型的扩展问题；③经典报童问题的2类扩展实例：带广告的报童问题和多产品报童问题。在这些文章中，目标函数都是围绕最大化利润的期望值。

本文将不再关注利润的期望值，而是研究与目标利润达成概率相关的两类目标函数：①最大化达成目标利润的概率；②约束目标利润的达成概率。

为更好地说明期望值模型与这两类新模型之间的差异，以下是一个不太严谨的例子：在股票市场中，期望值模型回答的是，投资方案A的平均利润为100万，但对于此次具体投资，利润并不保证；“最大化达成目标利润的概率”模型则表明，要达到100万以上的目标利润，投资方案B的实现概率最高；“约束目标利润的达成概率”模型则告诉我们，要达到100万以上的目标利润，且实现概率不低于90%，可以选择投资方案C。

通常情况下，决策者对不利情况的发生越敏感，越会关注目标利润的达成概率。接下来的内容，将主要介绍如何求解这两类问题。

正文如下。

2 最大化达成目标利润的概率

本节将继续使用报童模型来阐述该类问题的求解方案。

在报童模型中，我们定义如下变量：$ D$ 表示需求量， $g(D)$ 和 $G(D)$ 是对应的概率分布函数和累积分布函数； $Q$ 是订购量；$R$ 是售卖价格；$C$ 是成本；$V$是滞销损失；$S$ 是缺货损失。假设变量满足如下条件： $R>C>V$。

实际销量可以表示为：
$$A=\min (Q,D)$$

接下来，我们将分别讨论不考虑缺货损失和考虑缺货损失这两种情况。

2.1 不考虑缺货损失

此时$S=0$，报童模型的利润表达式可以写为
$$Z=RA+V\cdot \max (0,Q-D)-CQ$$
由于$\max (0,Q-D)=Q-\min (Q,D)$，所以上式等价于
$$Z=(R-V)A-(C-V)Q$$

假设目标利润为$B$，则最大化达成目标利润的概率可以表达为
$$\max \text{Pr}(Z\ge B)$$

为了求解该问题，假设$Q$为定值，此时$Z$的最大值为$Z_m=(R-C)Q$。因此，如果
$$B>(R-C)Q$$
那么
$$\text{Pr}(Z\ge B)=0$$
反之，令$Z=B$，可以得到对应的$D_B$为
$$D_B=\frac{B+(C-V)Q}{R-V}$$
因此
$$\text{Pr}(Z\ge B)=\text{Pr}(D\ge D_B)=1-G(D_B)$$
从上面几个式子可以看出，$Q$越小，$D_B$越小，$G(D_B)$越小，而$\text{Pr}(Z\ge B)$越大。因此，要最大化$\text{Pr}(Z\ge B)$，应该取$Q$的最小值。结合$B\le (R-C)Q$的约束，可得
$$Q^{\ast }=\frac{B}{R-C}$$
这里有个有意思的现象：$Q$的最优解，和$D$的分布没有任何关系。

2.2 考虑缺货损失

此时$S>0$，报童模型的利润表达式可以写为
$$Y=RA+V\cdot \max (0,Q-D)-S\max (0,D-Q)-CQ$$
由于$\max (0,D-Q)=S(D-A)$，所以上式等价于
$$Y=(R-V+S)A-SD-(C-V)Q$$

继续假设$Q$为定值，相比不考虑缺货损失时的$Z$，$Y$关于$D$不再是单调函数：当$D<Q$时，单调递增；当$D>Q$时，单调递减。即，给定目标利润$B$后，会同时存在两个恰好达成目标利润的$D$（不考虑目标利润达不成的情况）。如下图所示。

当$D<Q$时，令$Y=B$，可以得到
$$D_1=\frac{B+(C-V)Q}{R-V}$$

当$D>Q$时，仍令$Y=B$，可以得到
$$D_2=\frac{(R-C+S)Q-B}{S}$$
目标函数可以写为

$$P=\text{Pr}(Y\ge B)=G(D_2)-G(D_1)$$
此时，$Q$同时影响$D_1$和$D_2$的值，其和$P$的单调性关系不再容易判断。为了求得$P$的极值，直接求导
$$\frac{dP}{dQ}=g(D_2)\frac{d{D}_2}{dQ}-g(D_1)\frac{d{D}_1}{dQ}$$
将$D_1$和$D_2$的表达式带入上式，并令其值为0，得到
$$g(\frac{(R-C+S)Q-B}{S})\frac{R-C+S}{S}-g(\frac{B+(C-V)Q}{R-V})\frac{C-V}{R-V}=0$$
该等式的根$Q^{\ast }$即为最优解。

举个例子：假设$D\sim N(\mu ,{\sigma }^2)$，且$\mu =20$,$\sigma =4$,$R=6$,$C=4$,$S=4$,$V=1$,$B=20$。
运行以下代码，便可以得到，最优解$Q^{\ast }$的值为$20.94$，对应的最大概率为$0.7503$。

import math
from scipy.optimize import fsolve
import numpy as np
from scipy.stats import norm


# 迭代计算模型
def f1(x, mu, sigma):
    return norm.pdf((6 * x - 20) / 4, mu, sigma) * 3 / 2 - norm.pdf((3 * x + 20) / 5, mu, sigma) * 3 / 5


if __name__ == '__main__':

    mu = 20
    sigma = 4
    sol = fsolve(f1, mu, args=(mu, sigma, ))  # 使用mu值为初值
    max_prob = norm.cdf((6 * sol - 20) / 4, mu, sigma) - norm.cdf((3 * sol + 20) / 5, mu, sigma)
    print('sol: {}, max_prob: {}'.format(sol, max_prob))

3 约束目标利润的达成概率

为了求解该类问题，本节考虑如下的机会约束规划模型（CCP）
$$\max \overline{f}$$
KaTeX parse error: No such environment: eqnarray at position 7: \begin{̲e̲q̲n̲a̲r̲r̲a̲y̲}̲\nonumber \text…
式中，$\alpha$和$\beta$是约束条件和目标函数的达成概率，这意味着同时限制了目标利润的达成概率和约束条件的满足概率；如果不考虑约束条件，则退化为仅限制目标利润的达成概率。由于上述表达式可以统一为
Pr { g ( x , ξ ) ≤ 0 } ≥ α \text{Pr}\{g(\pmb x,\pmb \xi) ≤ 0 \} ≥ \alpha Pr{g(x,ξ)≤0}≥α
因此，在问题建模时，不必特别关注限制的是目标函数，还是约束条件的达成概率。

CCP的经典求解方法是将其转化为一个等价类，这可以理解为：将原包含随机变量的表达式转化为一个不含随机变量的新表达式，并保证新表达式的解与原公式的解保持一致。

接下来，我们将根据随机变量是一维还是多维分别进行分析。

3.1 一维变量

此时，$\xi$是一维变量，假设$g(x,\xi )=h(x)-\xi$，则上式转化为
Pr { g ( x ) ≤ ξ } ≥ α \text{Pr}\{g(\pmb x) ≤ \xi \} ≥ \alpha Pr{g(x)≤ξ}≥α
其中，$h(x)$可以为$x$的线性或非线性函数，$\xi$是一维随机变量，累积分布函数为$G(\cdot )$。

由于 Pr { g ( x ) ≤ ξ } = 1 − G ( g ( x ) ) \text{Pr}\{g(\pmb x) ≤ \xi \}=1-G(g(\pmb x)) Pr{g(x)≤ξ}=1−G(g(x))，带入上式，可以得到
$$g(x)\le 1-\alpha$$

举个实例：
Pr { x 1 2 − x 2 3 ≤ ξ } ≥ 0.90 \text{Pr}\{x_1^2 - x_2^3 ≤ \xi \} ≥ 0.90 Pr{x12​−x23​≤ξ}≥0.90
该约束条件就等价于
$$x_1^2-x_2^3\le G^{-1}(1-0.90)=0.7184$$

3.2 多维变量

此时，$\xi =(a_1,a_2,\cdot \cdot \cdot ,a_n,b)$，假设$g(x,\xi )$的表达式为
$$g(x,\xi )=a_1{x}_1+a_2{x}_2+\cdot \cdot \cdot +a_n{x}_n-b$$
上式中，$a_i$和$b$相互独立，且均为正态随机变量。

机会约束可以写为
Pr { ∑ i = 1 n a i x i ≤ b } ≥ α \text{Pr}\{\sum_{i=1}^na_ix_i ≤ b\} ≥ \alpha Pr{i=1∑n​ai​xi​≤b}≥α

令函数
$$y(x)=\sum _{i=1}^n{a}_i{x}_i-b$$
根据正态分布的叠加性，$y$也服从正态分布，且
$$\mu =\sum _{i=1}^{n}E(a_i)x_i-E(b)$$
$$\sigma =\sum _{i=1}^{n}V(a_i)x_i^2+V(b)$$
式中，$E(\cdot )$和$V(\cdot )$分别为期望值和方差。
再定义新的变量$X$
$$X=\frac{{\sum }_{i=1}^n{a}_i{x}_i-b-({\sum }_{i=1}^{n}E(a_i)x_i-E(b))}{\sqrt{{\sum }_{i=1}^{n}V(a_i)x_i^2+V(b)}}$$
显然，$X$服从标准正态分布。

针对不等式${\sum }_{i=1}^n{a}_i{x}_i\le b$，等价于
$$\frac{{\sum }_{i=1}^n{a}_i{x}_i-b-({\sum }_{i=1}^{n}E(a_i)x_i-E(b))}{\sqrt{{\sum }_{i=1}^{n}V(a_i)x_i^2+V(b)}}\le -\frac{{\sum }_{i=1}^{n}E(a_i)x_i-E(b)}{\sqrt{{\sum }_{i=1}^{n}V(a_i)x_i^2+V(b)}}$$

所以原约束条件等价于

Pr { X ≤ − ∑ i = 1 n E ( a i ) x i − E ( b ) ∑ i = 1 n V ( a i ) x i 2 + V ( b ) } ≥ α \text{Pr}\{X≤-\frac{\sum_{i=1}^nE(a_i)x_i - E(b)}{\sqrt{\sum_{i=1}^nV(a_i)x_i^2 + V(b)}}\} ≥ \alpha Pr{X≤−∑i=1n​V(ai​)xi2​+V(b) ​∑i=1n​E(ai​)xi​−E(b)​}≥α

到这里，就和一维变量的过程很相似了，此处直接给出上式的化简等价类如下
$$D^{-1}(\alpha )\le -\frac{{\sum }_{i=1}^{n}E(a_i)x_i-E(b)}{\sqrt{{\sum }_{i=1}^{n}V(a_i)x_i^2+V(b)}}$$
即
$$\sum _{i=1}^{n}E(a_i)x_i+D^{-1}(\alpha )\sqrt{\sum _{i=1}^{n}V(a_i)x_i^2+V(b)}\le E(b)$$

再举个实例：
Pr { a 1 x 1 + a 2 x 2 + a 3 x 3 ≤ b } ≥ 95 % \text{Pr}\{a_1x_1+a_2x_2+a_3x_3≤b\}≥95\% Pr{a1​x1​+a2​x2​+a3​x3​≤b}≥95%
其中$a_1,a_2,a_3$和$b$分别服从正态分布$N(1,1)$,$N(2,1)$,$N(3,1)$和$N(4,1)$。则该式的等价类为
$$x_1+2x_2+3x_3+1.645\sqrt{x_1^2+x_2^2+x_3^2+1}\le 4$$

4 总结

正文到此结束，以下是核心内容的总结：

(1) 最大化达成目标利润的概率：以报童问题为例，若不考虑缺货损失，最优解可以通过一个简单的公式得到； 而考虑缺货损失后，这个问题则可以转化为求解根的问题。

(2) 约束目标利润的达成概率：使用机会约束规划模型进行建模，根据随机变量是一维还是多维进行区分，当目标函数为特定类型时，可以将其转化为等价类进行求解。

5 相关阅读

The Newsboy Problem under Alternative Optimization Objectives：https://link.springer.com/article/10.1057/jors.1980.96

随机规划与模糊规划：https://book.douban.com/subject/1508496/