机器学习结合运筹学，有钱途～

学习笔记6—报童模型（最详细最全总结含公式推导和应用举例）

在这篇文章中我们介绍了标准的报童模型，在标准模型中需求分布是已知的，然而在现实世界中，需求往往是不确定的，不容易得知需求的分布。这样的问题也被称为不确定优化问题，常采用随机规划、鲁棒优化、分布式鲁棒优化、模糊规划以及机会约束优化等方法来解决。

该问题本质上是一个预测加优化问题，预测未知参数和优化求解。这里介绍了两种方法：数据驱动方法SAA和大数据驱动方法ERM。

首先回顾一下标准的报童模型，若需求分布F已知，最优订货量x*可由c1/(c1+c2)分位数给出。

一、数据驱动方法SAA（用历史数据的均值近似未知参数再进行优化）

这里先明确一下数据驱动的概念，和后面的大数据驱动做以区分。

数据驱动：基于数据进行决策，很广泛的概念

大数据驱动：大数据驱动属于数据驱动，但更强调使用大规模，复杂和高维度的数据集。大数据通常具有4v特征，即体积大（Volume），速度快（Velocity），多样性大（Variety），准确性（Veracity），通常依赖于高级的数据分析技术如机器学习ML、深度学习DL等。

Sample average approximation（SAA）是一种非常简单的可以估计不确定性参数的方法，即用历史数据的均值作为真实的一种近似。关于SAA方法的进一步细节见参考文献[1]。

[1] Shapiro A, Dentcheva D, Ruszczy´nski AP (2009) Lectures on Stochastic Programming: Modeling and Theory, Vol. 9 (SIAM, Philadelphia).

则标准报童模型的优化问题可转化为如下优化问题：

缺点：SAA的缺点在于没有考虑特征的信息，例如需求可能受价格、天气等因素的影响。

对于上面这个模型，用SAA方法，使用gurobi求解代码如下：

from gurobipy import *
import numpy as np
np.random.seed(3)

cost    = 0.3  #成本
retail  = 0.5 #售价
recover = 0.2 #滞销的价格
samples = 1000 #样本数量
sigma   = 50
mu      = 300
demand = np.random.normal(mu, sigma, samples)

m = Model()
m.ModelSense = -1

order    = m.addVar(name='order')
profit   = m.addVars(samples,obj=1.0/samples,name='profit')
sales    = m.addVars(samples,ub=demand,name='sales')
discount = m.addVars(samples,name='discount')

m.addConstrs((profit[i] == sales[i] * retail - order * cost + recover * discount[i] for i in range(samples)),name='profit')
m.addConstrs((sales[i]+discount[i] == order for i in range(samples)),name='demand')
m.update()
m.optimize() 

print("进货量:",order.x)
print("总体利润:",m.objVal)

结果如下，可见取1000个样本的结果算出来跟最优结果322很接近，毕竟是近似算法嘛。

二、大数据驱动方法ERM（用机器学习或深度学习预测未知参数和优化）

大数据驱动方法考虑了特征的信息（例如影响报纸需求的时间、天气、地点以及其他经济指标，统称为特征feature），则优化问题变为了如下所示(论文中符号和我们前面不一样哈)：

采用ERM算法解决上述优化问题，ERM算法说简单点就是训练一个预测模型使其误差函数最小（比如常用的回归评价指标MAE，MSE，RMSE），则问题就变成了：

具体内容见参考文献[2]。

[2] Ban, G.Y., Rudin, C., 2018. The big data newsvendor: Practical insights from machine learning. Operations Research 67, 90–108.

进一步的，根据模型的误差函数形式可将大数据驱动方法ERM分为两阶段模型（predict, then optimize）和一阶段模型（smart predict, then optimize），统称为端对端预测后优化，即输入端到输出端。两个模型的区别就在于两阶段模型预测模型采用的误差函数是预测误差（例如MAE，MSE），预测和优化是分开的，先预测，后优化，这属于机器学习和运筹学的简单结合；而一阶段模型采用的是决策误差（就是求出来的解对应的目标函数值和最优解对应的目标函数值之间的误差），这里预测和优化是结合在一起的，就属于机器学习和运筹学的深度结合。很显然，预测好不一定决策好，一阶段模型的结果往往会更精确一点。那如果模型的预测精度能达到100%呢，两阶段模型是不是就好了，显然不太可能。

两阶段模型（predict, then optimize）：这里的特征是随机生成的，我们可以将自己的数据代入进去。

from sklearn.datasets import make_regression
from sklearn.neural_network import MLPRegressor
import numpy as np
import gurobipy as gp
from gurobi_ml import add_predictor_constr
from gurobipy import GRB

np.random.seed(3)

cost    = 0.3  #成本
retail  = 0.5 #售价
recover = 0.2 #滞销的价格
samples = 1000 #样本数量

n_features, n_targets = 10, 1 # 假设需求有10个特征
X, y = make_regression(n_features=n_features, n_samples = samples, n_targets = n_targets, noise=1.0) # 随便生成一点学习训练所需数据
y = np.random.normal(300, 50, samples)  # 生成的y里有负值，这里替换一下，x则保留

pred_model = MLPRegressor([30, 20], max_iter=100, random_state=1) # 建立神经网络
pred_model.fit(X, y) # 训练神经网络

# 随机抽取一个X和y
index = np.random.choice(X.shape[0], 1, replace=False)
X_examples = X[index, :]
y_examples = y[index]

m = gp.Model()
m.ModelSense = -1
order    = m.addVar(name='order')
profit   = m.addVars(samples,obj=1.0/samples,name='profit')
sales    = m.addVars(samples,ub=y,name='sales')
discount = m.addVars(samples,name='discount')
input_vars = m.addMVar(X_examples.shape, lb=X_examples, ub=X_examples) # 预测模型输入变量
output_vars = m.addMVar(y_examples.shape, lb=-gp.GRB.INFINITY, ub=gp.GRB.INFINITY)# 预测模型输出变量
pred_constr = add_predictor_constr(m, pred_model, input_vars, output_vars) # 添加模型约束
pred_constr.print_stats()
m.addConstrs((profit[i] == sales[i] * retail - order * cost + recover * discount[i] for i in range(samples)),name='profit')
m.addConstrs((sales[i]+discount[i] == order for i in range(samples)),name='demand')
m.update()
m.optimize() 

print("进货量:",order.x)
print("总体利润:",m.objVal)

一阶段模型（smart predict, then optimize）：一阶段模型在两阶段模型的基础上自定义预测模型的误差函数即可，需要注意的是MLP不支持自定义误差函数，可以用pytorch构建神经网络，然后训练模型更新参数。

思路很清晰，但是有一个问题是线性规划、整数线性规划等具有线性目标函数的问题，其最优解关于需求的函数一阶导数要么为0，要么不存在。那么神经网络参数更新就继续不下去。面对这一挑战，有三种解决方法：一种是寻找替代的梯度信息以更新模型参数；一种是重新设计一个替代损失函数；另一种是选择基于树的预测模型[3]。

[3] Elmachtoub, A. N., Liang, J. C. N., & McNellis, R. (2020, November). Decision trees for decision-making under the predict-then-optimize framework. In International Conference on Machine Learning (pp. 2858-2867). PMLR.

关于更多的信息请参考下面这篇文章。作者还开发了一个pyepo工具将这些方法都写在了里面，真的是大树底下好乘凉呀，在这里感谢这些大佬的辛勤付出！

端对端预测后优化 优化｜当机器学习上运筹学：PyEPO与端对端预测后优化 - 知乎