数据驱动不确定集实战：基于 Python 与历史数据构建 3 种鲁棒优化模型

数据驱动不确定集实战基于 Python 与历史数据构建 3 种鲁棒优化模型当数据科学家面对现实世界中的决策问题时常常会遇到一个核心挑战如何在数据不确定性下做出最优决策传统优化方法假设输入参数是确定的但实际业务场景中需求预测、供应链波动、市场变化等关键参数往往存在显著波动。鲁棒优化提供了一种系统化的解决方案——通过构建合理的不确定集Uncertainty Set在最坏情况下仍能保证可行性和性能的决策方案。本文将聚焦数据驱动的不确定集构建方法通过Python实战演示如何利用历史数据构建三类典型不确定集盒式、椭球交集、基于矩信息并完整实现对应的鲁棒优化模型求解。不同于理论推导我们更关注从原始数据到可部署代码的完整链路每个模型均提供可直接运行的Jupyter Notebook示例。1. 环境准备与数据预处理在开始构建不确定集之前需要确保Python环境已安装必要的科学计算和优化库。推荐使用Anaconda创建专属环境conda create -n robust_opt python3.9 conda activate robust_opt pip install numpy pandas matplotlib cvxpy scipy statsmodels假设我们有一个制造业的生产计划优化场景历史需求数据存储在CSV文件中demand_history.csv包含过去12个月每日产品需求记录。首先进行数据探索和预处理import pandas as pd import numpy as np import matplotlib.pyplot as plt # 加载历史数据 demand pd.read_csv(demand_history.csv, parse_dates[date]) demand.set_index(date, inplaceTrue) # 可视化需求分布 plt.figure(figsize(12, 6)) demand[product_A].plot(titleHistorical Demand for Product A) plt.ylabel(Daily Demand) plt.grid(True) plt.show() # 计算统计特征 stats demand.describe().T[[mean, std, min, max]] stats[cov] stats[std] / stats[mean] # 变异系数 print(stats.round(2))关键预处理步骤包括异常值处理使用Tukey方法识别并修正异常值平稳性检验ADF检验确认时间序列平稳性相关性分析计算产品间的需求相关系数矩阵from statsmodels.tsa.stattools import adfuller # ADF平稳性检验 def check_stationarity(series): result adfuller(series.dropna()) print(fADF Statistic: {result[0]:.3f}) print(fp-value: {result[1]:.3f}) check_stationarity(demand[product_A]) # 多变量相关性矩阵 corr_matrix demand.corr() plt.matshow(corr_matrix, cmapcoolwarm) plt.colorbar() plt.title(Product Demand Correlation) plt.show()2. 盒式不确定集构建与线性鲁棒优化盒式不确定集Box Uncertainty Set是最基础的形式假设每个不确定参数在独立区间内波动。虽然结构简单但在许多场景下仍能提供有价值的鲁棒解。2.1 基于历史分位数的区间估计从历史数据中提取需求波动的合理边界# 设置置信水平 alpha 0.9 # 计算双边分位数 lower demand.quantile((1-alpha)/2) upper demand.quantile(1 - (1-alpha)/2) # 可视化产品A的区间估计 plt.figure(figsize(10, 5)) plt.fill_between(demand.index, lower[product_A], upper[product_A], alpha0.3, labelUncertainty Interval) demand[product_A].plot(labelActual Demand) plt.legend() plt.title(fDemand Range Estimation (Confidence: {alpha*100}%)) plt.show()2.2 生产计划鲁棒优化模型考虑一个多产品生产计划问题目标是最小化总成本生产成本库存成本同时满足需求约束import cvxpy as cp # 定义参数 n_products len(demand.columns) cost_prod np.array([2.5, 3.0, 4.2]) # 单位生产成本 cost_inv np.array([0.3, 0.4, 0.5]) # 单位库存成本 capacity np.array([200, 250, 180]) # 最大生产能力 # 决策变量 x cp.Variable(n_products) # 生产量 s cp.Variable(n_products) # 库存量 # 不确定参数需求 d cp.Parameter(n_products) d.value demand.mean().values # 初始设为均值 # 目标函数 objective cp.Minimize(cost_prod x cost_inv s) # 确定性约束 constraints [ x capacity, x 0, s 0 ] # 鲁棒约束满足所有可能需求 robust_constraints [ x s d (upper.values - demand.mean().values), # 上界偏离 x s d - (demand.mean().values - lower.values) # 下界偏离 ] # 合并约束 problem cp.Problem(objective, constraints robust_constraints) problem.solve() print(fOptimal Production: {x.value.round(1)}) print(fOptimal Inventory: {s.value.round(1)}) print(fTotal Cost: {problem.value:.2f})2.3 保守性分析与改进盒式不确定集的主要缺点是可能过于保守。通过计算实际需求被包含的百分比来评估保守程度# 计算历史需求被包含的比例 contained ((demand lower.values) (demand upper.values)).mean() print(fDemand contained ratio: {contained.mean():.1%}) # 调整alpha值平衡保守性与性能 alpha_values np.linspace(0.7, 0.95, 6) results [] for a in alpha_values: l demand.quantile((1-a)/2).values u demand.quantile(1 - (1-a)/2).values # 更新鲁棒约束 constraints[-2].args[1] d (u - demand.mean().values) constraints[-1].args[1] d - (demand.mean().values - l) problem.solve() results.append({ alpha: a, cost: problem.value, production: x.value.copy() })3. 椭球交集不确定集与二次锥优化椭球型不确定集能捕捉参数间的相关性提供更精确的不确定性描述。我们重点介绍椭球交集Ellipsoidal Intersection形式。3.1 构建数据驱动的椭球集基于历史数据的均值和协方差矩阵构建椭球from scipy.stats import chi2 # 计算统计量 mu demand.mean().values Sigma demand.cov().values n_samples len(demand) # 卡方分布临界值控制椭球大小 rho chi2.ppf(0.95, dfn_products) # 椭球集数学表达 # { d | (d - μ)^T Σ^{-1} (d - μ) ≤ ρ }3.2 鲁棒优化模型重构将椭球不确定集融入生产计划问题形成二次锥规划# 新增辅助变量 t cp.Variable() y cp.Variable(n_products) # 修改目标函数 objective cp.Minimize(cost_prod x cost_inv s t) # 椭球不确定集约束 constraints [ cp.SOC(t, y), x s mu y, cp.quad_form(y, np.linalg.inv(Sigma)) rho ] problem cp.Problem(objective, constraints) problem.solve(solvercp.ECOS) print(fRobust Solution (Ellipsoid):) print(fProduction: {x.value.round(1)}) print(fCost: {problem.value:.2f})3.3 与传统盒式集的对比实验通过蒙特卡洛模拟比较两种不确定集下的性能np.random.seed(42) n_sim 1000 # 生成测试场景多元正态分布 test_demand np.random.multivariate_normal(mu, Sigma, n_sim) test_demand np.maximum(test_demand, 0) # 需求非负 # 评估两种方案 def evaluate_solution(x_val, s_val, demands): shortage np.maximum(demands - x_val - s_val, 0) excess np.maximum(x_val s_val - demands, 0) cost cost_prod x_val cost_inv excess return cost.mean(), shortage.sum(axis1).mean() # 盒式集方案 box_cost, box_shortage evaluate_solution( results[-1][production], s.value, test_demand) # 椭球集方案 ellipsoid_cost, ellipsoid_shortage evaluate_solution( x.value, s.value, test_demand) print(fBox Set - Avg Cost: {box_cost:.2f}, Avg Shortage: {box_shortage:.2f}) print(fEllipsoid Set - Avg Cost: {ellipsoid_cost:.2f}, Avg Shortage: {ellipsoid_shortage:.2f})4. 基于矩信息的模糊集与分布鲁棒优化分布鲁棒优化DRO通过矩信息构建模糊集Ambiguity Set比传统鲁棒优化更充分利用数据统计特征。4.1 矩模糊集构建利用历史数据的一、二阶矩信息# 矩模糊集参数 epsilon 0.1 # 模糊度参数 support np.vstack([lower.values, upper.values]).T # 支撑集 # 模糊集数学表达 # { P ∈ P(Ξ) | E_P[ξ] μ, E_P[(ξ-μ)(ξ-μ)^T] ⪯ Σ εI }4.2 数据驱动的DRO模型将分布鲁棒约束转化为可求解形式# 引入对偶变量 lambda_ cp.Variable() M cp.Variable((n_products, n_products), PSDTrue) # 修改目标函数 objective cp.Minimize(cost_prod x cost_inv s lambda_ * epsilon cp.trace(Sigma M)) # 对偶约束 constraints [ M ⪰ 0, cp.quad_form(x s - mu - lambda_ * np.ones(n_products), np.linalg.inv(M)) 1 ] problem cp.Problem(objective, constraints) problem.solve(solvercp.MOSEK) print(fDRO Solution:) print(fProduction: {x.value.round(1)}) print(fProtection Cost: {lambda_.value:.2f})4.3 模糊度参数敏感性分析研究ε取值对解决方案的影响epsilons np.logspace(-2, 1, 10) dro_results [] for eps in epsilons: lambda_.value 0 M.value np.eye(n_products) # 更新目标函数 objective cp.Minimize(cost_prod x cost_inv s lambda_ * eps cp.trace(Sigma M)) problem.solve(solvercp.MOSEK) dro_results.append({ epsilon: eps, cost: problem.value, production: x.value.copy(), lambda: lambda_.value }) # 可视化分析 plt.figure(figsize(10, 5)) plt.semilogx(epsilons, [r[cost] for r in dro_results], o-) plt.xlabel(Ambiguity Parameter (ε)) plt.ylabel(Total Cost) plt.title(DRO Cost Sensitivity Analysis) plt.grid(True) plt.show()5. 模型选择与实施建议面对三种不同的不确定集建模方法实际应用中如何选择以下提供决策框架标准盒式集椭球交集矩模糊集计算复杂度低线性规划中二次锥规划高半定规划数据需求边缘统计量协方差矩阵一、二阶矩保守程度高中可调节参数相关性处理能力无强强适用场景快速原型开发中等规模问题高价值决策实施时的关键注意事项数据质量检查确保历史数据具有代表性非平稳数据需先进行差分或分解模型验证通过回测backtesting评估各模型在实际波动下的表现计算效率优化对大规模问题考虑稀疏矩阵表示使用分解算法加速求解鲁棒性-性能权衡通过调整参数α、ρ、ε控制系统保守程度# 示例稀疏矩阵处理 from scipy import sparse # 当产品数量大时使用稀疏协方差矩阵 if n_products 50: Sigma_sparse sparse.csc_matrix(Sigma) # 修改椭球约束为稀疏形式 # ...