随机存储模型（随机需求下的库存控制）

在确定性的模型中，我们把需求看成是固定不变的已知常量。但是，在现实世界中，更多的情况却是需求为一个随机变量。为此，在本节中我们将介绍需求是随机变量，特别是需求服从均匀分布和正态分布这两种简单情况的存贮模型。

一、报童模型（Newsvendor Model）

报童模型是运筹学和库存管理中的一个经典模型，用于解决单周期库存问题（需求为随机变量的单一周期存储模型）。该模型适用于那些需求具有不确定性且只能在销售周期开始前进行一次性订购的商品，如报纸、鲜花、季节性商品等。报童模型的目标是确定最佳订购量，以最大化预期利润或最小化预期成本。

报童每天销售报纸的数量是一个随机变量，根据以往的经验，每日售出$r$份报纸的概率$P(r)$是已知的。报童每售出一份报纸赚 $k$ 元，如果报纸未能售出，每份赔 $h$ 元，问报童每日最好准备多少报纸$Q$，才能使损失的期望值最小呢？

符号定义

• r：每日实际需求的随机变量，其概率分布 $P(r)$ 已知。
• k：每售出一份报纸的利润（单位：元）。
• h：每份未售出报纸的损失（单位：元）。

决策变量

• $Q$：报童每日准备的报纸数量。

1.1 基于期望损失最小化的最优订购量推导

数学模型
若供大于求（$Q>r$），则报纸没卖完的损失：
$$\sum _{r=0}^{Q}h(Q-r)p(r)$$
若供小于求（$Q\le r$），则因为缺货少赚钱的成本：
$$\sum _{r=Q+1}^{\infty }k(r-Q)p(r)$$
因此期望损失为：
$$\min C(Q)=\sum _{r=0}^{Q}h(Q-r)p(r)+\sum _{r=Q+1}^{\infty }k(r-Q)p(r)$$

设最优订货量为$Q^{\ast }$，则有
$$\{\begin{array}{l}C(Q^{\ast })\le C(Q^{\ast }+1)\\ C(Q^{\ast })\le C(Q^{\ast }-1)\end{array}$$

$$\begin{array}{rl}& C(Q^{\ast })\le C(Q^{\ast }+1)\\ ⟺& \sum _{r=0}^{Q^{\ast }}h(Q^{\ast }-r)p(r)+\sum _{r=Q^{\ast }+1}^{\infty }k(r-Q^{\ast })p(r)\le \sum _{r=0}^{Q^{\ast }+1}h(Q^{\ast }+1-r)p(r)+\sum _{r=Q^{\ast }+2}^{\infty }k(r-Q^{\ast }-1)p(r)\\ ⟺& \sum _{r=0}^{Q^{\ast }}p(r)\ge \frac{k}{k+h}\\ & \\ & C(Q^{\ast })\le C(Q^{\ast }-1)\\ ⟺& \sum _{r=0}^{Q^{\ast }-1}p(r)\le \frac{k}{k+h}\end{array}$$

故最佳订货量$Q^{\ast }$满足不等式：
$$\sum _{r=0}^{Q^{\ast }-1}p(r)\le \frac{k}{k+h}\le \sum _{r=0}^{Q^{\ast }}p(r)$$
可以把上式改为：
$$P(r<Q^{\ast })\le \frac{k}{k+h}\le P(r\le Q^{\ast })$$
如果只考虑连续型随机变量，上式又可以写为：
$$P(r\le Q^{\ast })=\frac{k}{k+h}$$

1.2 基于期望利润最大化的最优订购量推到

在报童问题中我们的目标是使订购量过剩和过小两种损失最小从而推出公式$P(X<Q^{\ast }-1)\le \frac{k}{k+h}\le P(X\le Q^{\ast })$，能否改变一下优化目标使期望利润最大，推出该公式?

报童的利润取决于实际需求 r 和准备数量 Q 的关系。分为两种情况：

1. 供小于求（$Q<r$）​售出数量为$Q$，利润为$k\cdot Q$。
2. 供大于求（$Q\ge r$）​售出数量为$r$。未售出数量为 $Q-r$，损失为 $h\cdot (Q-r)$，利润为 $k\cdot r-h\cdot (Q-r)$。

将上述两种情况统一，利润函数可以表示为：
$$G(Q,r)=\{\begin{array}{ll}k\cdot Q& \text{if}Q<r\\ k\cdot r-h\cdot (Q-r)& \text{if}Q\ge r\end{array}$$
综合两种情况，期望利润为：
$$G(Q)=\sum _{r=0}^Q\left[k\cdot r-h\cdot (Q-r)\right]P(r)+\sum _{r=Q+1}^{\infty }k\cdot Q\cdot P(r)$$

​潜在最大利润：当订购量 Q 完全匹配需求 r 时，无任何过剩或缺货损失，此时利润最大：
$$G_{\text{max}}=\sum _{r=0}^{\infty }k\cdot r\cdot P(r)$$
实际利润$G(Q)$等于潜在最大利润减去总损失$C(Q)$：
$$G(Q)=G_{\text{max}}-C(Q)$$故
$$\max G(Q)\iff \min C(Q)$$两者的最优解$Q^{\ast }$必然相同，因为最大化$G(Q)$等价于最小化$C(Q)$。

二、算例

2.1 算例一

某报亭出售某种报纸，每售出一百张可获利15元，如果当天不能售出，每一百张赔20元。每日售出该报纸份数的概率P(d )根据以往经验如下表所示。试问报亭每日订购多少张该种报纸能使其赚钱的期望值最大。

销售量（百张）	5	6	7	8	9	10	11
概率P(x)	0.05	0.10	0.20	0.20	0.25	0.15	0.05

解：要使其赚钱的期望值最大，也就是使其因售不出报纸的损失和因缺货失去销售机会的损失的期望值之和为最小。已知 k = 15，h = 20，则有
$$\frac{k}{k+h}=\frac{15}{15+20}=0.4286$$
另有
$$\sum _{x=0}^{7}P(x)=P(5)+P(6)+P(7)=0.35$$
$$\sum _{x=0}^{8}P(x)=P(5)+P(6)+P(7)+P(8)=0.55$$
故$Q=8$时，不等式
$$\sum _{x=0}^{7}p(x)\le \frac{k}{k+h}\le \sum _{r=0}^{8}p(x)$$
成立，因此最优订货量为800张。

---

2.2 算例二

某书店拟在年前出售一批新年挂历。每售出一本可盈利20元，如果年前不能售出，必须削价处理。由于削价，一定可以售完，此时每本挂历要赔16元。根据以往的经验，市场的需求量近似服从均匀分布，其最低需求为550本，最高需求为1100本，该书店应订购多少新年挂历，使其损失期望值为最小？
$$f(x)=\{\begin{array}{ll}\frac{1}{b-a}& \text{if }a\le x\le b\\ 0& \text{otherwise}\end{array}$$
解：由题意知挂历的需求量是服从区间[550，1100]上的均匀分布的随机变量， k = 20，h = 16，则其需求量小于$Q^{\ast }$的概率为
$$P(x\le Q^{\ast })=\frac{Q^{\ast }-550}{1100-550}=\frac{20}{20+16}$$
由此得$Q^{\ast }=856$（本）

---

2.3 算例三

某化工公司与一客户签订了一项供应一种独特的液体化工产品的合同。客户每隔六个月来购买一次，每次购买的数量是一个随机变量，通过对客户以往需求的统计分析，知道这个随机变量服从以均值$\mu =1000$ （公斤），方差$\sigma =100$ （公斤）的正态分布。化工公司生产一公斤此种产品的成本为15元，根据合同固定售价为20元。合同要求化工公司必须按时提供客户的需求。一旦化工公司由于低估了需求产量不能满足需要，那么化工公司就到别的公司以每公斤19元的价格购买更高质量的替代品来满足客户的需要。一旦化工公司由于高估了需求，供大于求，由于这种产品在两个月内要老化，不能存贮至六个月后再供应给客户，只能以每公斤5元的价格处理掉。化工公司应该每次生产多少公斤的产品才使该公司获利的期望值最大呢？

根据题意得 k =5－1= 4，h = 15－5= 10

$$P(x\le Q^{\ast })=\frac{k}{k+h}=\frac{4}{4+10}=0.29$$

需求为随机变量服从以均值$\mu =1000$ （公斤），方差$\sigma =100$ （公斤）的正态分布，即
$$\Phi (\frac{Q^{\ast }-\mu }{\sigma })=0.29$$

通过查阅标准正态表，得
$$\frac{Q^{\ast }-\mu }{\sigma }=0.55$$
将$\mu =1000$，$\sigma =100$代入得$Q^{\ast }=-0.55\times 100+1000=945$公斤

2.3.1 使用蒙特卡罗方法进行求解

import numpy as np
import matplotlib.pyplot as plt

# 参数定义
mu = 1000  # 需求均值
sigma = 100  # 需求标准差
N = 10000  # 样本数量
k = 4 
h = 10

# 生成需求样本
demand_samples = np.random.normal(mu, sigma, N)

# 可能的订货量范围
Q_range = np.arange(500, 1500, 1)  # 假设订货量在500到1500之间

# 存储每个订货量的损失成本
average_cost_list = []

for Q in Q_range:
    cost_list = []
    for D in demand_samples:
        if D <= Q:
            cost = h * (Q - D)
        else:
            cost = k * (D - Q)
        cost_list.append(cost)
    average_cost_list.append(np.mean(cost_list))

# 找到最小平均成本对应的订货量
optimal_Q = Q_range[np.argmin(average_cost_list)]
print(f"Optimal order quantity: {optimal_Q}")

# 可视化结果
plt.plot(Q_range, average_cost_list)
plt.xlabel("Order Quantity (Q)")
plt.ylabel("Average cost")
plt.title("Order Quantity vs. Average cost")
plt.axvline(optimal_Q, color="r", linestyle="--", label=f"Optimal Q: {optimal_Q}")
plt.legend()
plt.show()

多次运行接近解析解945，运行结果如下。

PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
Optimal order quantity: 942
PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
Optimal order quantity: 945
PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
Optimal order quantity: 945
PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
Optimal order quantity: 944

2.3.2 使用样本均值近似算法求解

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# 参数定义
mu = 1000  # 需求均值
sigma = 100  # 需求标准差
N = 10000  # 样本数量
k = 4
h = 10

# 生成需求样本
demand_samples = np.random.normal(mu, sigma, N)

# 计算临界比率
critical_ratio = k / (k + h)

# 计算标准正态分布的分位数
z = norm.ppf(critical_ratio)

# 计算最优订货量
optimal_Q = mu + z * sigma
print(f"分位数{z}")
print(f"Optimal order quantity: {optimal_Q}")

# 可视化结果
plt.hist(demand_samples, bins=50, alpha=0.75, label='Demand Samples')
plt.axvline(optimal_Q, color='r', linestyle='--', label=f'Optimal Q: {optimal_Q:.2f}')
plt.xlabel('Demand')
plt.ylabel('Frequency')
plt.title('Demand Distribution and Optimal Order Quantity')
plt.legend()
plt.show()

PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
分位数-0.5659488219328631
Optimal order quantity: 943.4051178067136
PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
分位数-0.5659488219328631
Optimal order quantity: 943.4051178067136
PS C:\Users\pengkangzhen\PycharmProjects\operation-research> & C:/Users/pengkangzhen/anaconda3/envs/py3.11_ml/python.exe c:/Users/pengkangzhen/PycharmProjects/operation-research/蒙特卡罗/报童模型.py
分位数-0.5659488219328631
Optimal order quantity: 943.4051178067136