做微网调度这些年我最深的一个体会是调度方案的“可用性”比“好看”重要得多。风光预测曲线画得再漂亮实际运行的时候风速一变、云层一压原计划当场失效。这也是为什么两阶段鲁棒优化在微网优化调度里越来越受重视——它不赌预测准不准而是先把不确定性的边界画出来然后在最坏的情况下也能给你一套兜底可用的调度策略。这篇文章就围绕“基于关键场景辨别算法的两阶段鲁棒微网优化调度Matlab代码实现”这个具体的案例来拆。核心关键词拆开是三块微网系统建模、两阶段鲁棒优化框架、用来识别恶劣场景的关键场景辨别算法。我会把数学模型怎么搭、主问题子问题怎么迭代、Matlab里怎么一步步实现、以及调试过程中那些容易翻车的细节全部讲透给想自己动手复现的同行一条能走通的路。1. 从“平均最优”到“最坏情况最优”两阶段鲁棒优化的设计逻辑1.1 为什么单阶段确定性调度在微网面前撑不住传统的确定性调度本质上是在做一件事把预测值当成真实值求解一个单层优化问题。你早上根据负荷预测曲线、光伏出力预测曲线、风电出力预测曲线再加上电价信息求解出未来24小时各机组出力和储能充放电计划这个计划在预测完全准确的前提下当然是最优的。问题是微网里的不确定性源太多了。光伏出力受云层遮挡影响分钟级波动可以超过装机容量的50%风电出力更是看天吃饭风速预测偏差10%功率偏差可以到30%以上负荷侧的大型冲击性负载一旦启动预测曲线直接失真。实际运行的时候你按预测值定好的储能充放电计划很可能让系统出现功率失衡这时候要么弃风弃光要么切负荷要么高价从上级电网买电经济性和可靠性双双失守。两阶段鲁棒优化的思路完全不一样。它把决策拆成两段第一阶段是在不确定性还没暴露之前就得拍板的变量比如机组启停状态、储能充放电基准计划第二阶段是等到不确定性实现之后系统还能做的“矫正动作”比如调整燃气轮机出力、调用储能余量、调节与主网的交换功率。优化目标不再追求预测场景下的期望最优而是确保在最恶劣的不确定性实现下第二阶段的调整成本依然可控、系统依然可行。用大白话讲确定性调度像是你根据天气预报决定明天穿什么赌的是预报准两阶段鲁棒像是你既准备了一件薄外套又包里塞了一件厚冲锋衣——天气真变了你还有第二层调整手段。代价是方案会稍微“保守”一点但换来的是极强的抗风险能力。1.2 两阶段鲁棒优化的数学骨架两阶段鲁棒优化的一般形式长这样min c1^T x max u∈U min y c2^T y这个式子要拆开看。x 是第一阶段决策变量对应那些一旦确定就很难改变的物理量比如机组启停启动后不能立刻停机、储能的日前充放电计划。这批变量的特点是“现在定后面不能改”术语叫 here-and-now 决策。u 是不确定变量也就是我们关心的光伏出力、风电出力、负荷功率。它落在不确定集 U 里面U 刻画了我们对“不确定性到底有多坏”的认知边界。max 这一层的意思是既然你不知道真实场景是哪个那就假设命运会跟你对着干在所有可能的 u 里挑一个让你的成本最高的。这不是悲观主义而是鲁棒优化的基本假设——我做的方案必须经得起这种最坏情况的考验。y 是第二阶段决策变量术语叫 wait-and-see。等你已经知道了 u 的真实取值系统就能根据实际情况做调整比如让储能多放一点电、让燃气轮机多带一点负荷。min y c2^T y 就是在最坏 u 下第二阶段的调整成本能做到多小。换个角度理解这个三层结构是对鲁棒性的一个“压力测试”。先做一个决策 x然后让不确定性尽量刁难你最后再看你还有多少补救能力。如果一个调度计划在这种连环压力测试下依然可行且成本可接受那它在真实运行中基本上就稳了。综合起来整个优化问题在做的是找一个第一阶段决策 x使得“第一阶段的固定成本 最坏情况下的第二阶段最小调整成本”达到最小。1.3 不确定集给“最坏情况”划好边界两阶段鲁棒里最关键的一个工程设计是 U 的选取。U 太小鲁棒性形同虚设U 太大调度过于保守经济性崩盘。实际中最常用的是盒式不确定集加预算约束。用风电出力来举例。设风电场预测出力是 P_w^forecast实际出力可以写为P_w,actual P_w,forecast ΔP_w其中 ΔP_w 是预测误差它被限制在一个区间内ΔP_w ∈ [-ΔP_w^max, ΔP_w^max]这就是“盒式”的意思——把每个不确定变量的偏差都框在一个矩形区间里。如果只做到这一步最坏情况就是所有风电机组同时偏到区间下限这在实际中几乎不可能发生方案会极度保守。预算约束的作用就是控制“同时在边界上的不确定变量数量”。比如微网里有10台分布式风机预算 Γ3那最坏情况只允许至多3台风机同时取到区间端点的偏差值剩余7台的偏差只能落在区间内部。这样 U 的形状从一个“满盒子”被切削成了“带棱角的超多面体”对应的最坏场景更贴近物理现实解的经济性也明显改善。用公式写预算约束就是∑ |ΔP_w| / ΔP_w^max ≤ Γ直观理解Γ 越大系统越保守Γ0 时退化为确定性模型Γ 取到不确定变量总数时就是纯盒式。实际调试中Γ 是调节方案“胆量”的核心旋钮。微网里风光渗透率越高Γ 通常要取得越大因为不确定性的相关性更强极端情况的破坏力更大。2. 关键场景辨别算法如何从连续不确定集里精准锁定“最狠”的场景2.1 问题的麻烦不确定集是连续的没法直接枚举把两阶段鲁棒模型摆出来之后第一个让人头疼的问题就是max u∈U 这一层怎么算。如果 U 是一个盒式连续区间这就意味着有无限多个可能的 u。一个由有限离散场景构成的集合可以用枚举来求解但连续集合没法一轮轮全部列出来。更麻烦的是子问题是一个 max-min 结构外层让 u 尽量坏内层让 y 尽量补救两股力量互相拉扯直接求解本身就有难度。这时候就得用到“生成”的思想我不需要把 U 里所有场景都拿出来参与计算我只关心那些对当前解打击最大的关键场景。每一轮迭代我固定住第一阶段的 x去求解子问题找到一个让第二阶段成本最大化的 u——这个 u 就是当前决策下的“最恶劣场景”。然后把 u 作为新场景加入主问题的场景集合重新求解主问题更新 x。这样做的好处非常明显每次迭代都在“定向挖掘”最坏场景而不是靠随机抽样去碰运气。它本质上是把连续优化问题转化成一个不断扩大的有限场景集合上的离散优化问题。场景集合不断扩大鲁棒解就不断逼近真正的两阶段鲁棒最优解。这个思想在学术界叫列与约束生成算法而“关键场景辨别”这个词把它的核心动作描述得更直白——在每一轮里通过子问题的求解从不确定集里辨别出当前威胁最大的那个点。2.2 关键场景辨别与传统场景削减法不是一回事很多做微网的人第一反应是不确定性处理不是可以先蒙特卡洛抽样再用K-means或者快速前向选择做场景削减吗这两条技术路线在目标上是有本质区别的。场景削减法backward reduction、K-means聚类等处理的是“概率意义上有代表性的场景”。它立足于历史数据把大量随机场景聚成几个典型簇然后用簇心代表整类场景。这套方法用于随机优化stochastic programming非常合适因为随机优化算的是期望成本需要用概率分布把场景加权起来它关注的是“平均怎么样”。但鲁棒优化关心的是“最坏怎么样”。假如你按概率聚类挑了5个高频出现的温和场景忽略掉了2个虽然概率低但一旦出现就代价巨大的极端场景那么你的调度方案在实际遇到极端天气时可能直接就不可行了。鲁棒优化的核心诉求恰恰是要覆盖这种极端情况所以场景选取标准必须围绕“对目标函数的破坏力”来设计而不是“出现的频次”。关键场景辨别算法在每一轮迭代时都是严格对着当前调度决策去找最恶劣场景的这个场景可能概率不高但一定是对解构成真实威胁的那个点。用行话讲它有“最优性导向”场景的产生不是独立于优化模型而是由优化模型的子问题内生出来的。2.3 子问题是怎么一步步把最恶劣场景“逼”出来的子问题的形式是一个 max-min 优化问题。内层 min y 是在给定 u 之后找到一个最经济的补救方案外层 max u 是在 U 里折腾所有可能的 u看哪个让我内层补救也救不回来、成本崩到最高。直接求双层优化很难但可以利用线性规划的对偶理论把内层 min 问题转成对偶的 max 问题这样结构就变成“max 对偶问题”和“max 不确定变量”挤在了一层不再有嵌套。内层原始问题是min c2^T y s.t. D y ≤ e − E x_current − F u它的对偶问题写出来是max λ^T (e − E x_current − F u) s.t. D^T λ c2 λ ≥ 0把对偶问题代回原结构后子问题变成max λ^T (e − E x_current − F u) s.t. λ ≥ 0, D^T λ c2, u ∈ U这里注意看目标函数里 λ 和 u 相乘出现了一个双线性项。处理这个双线性项通常有两种路径一是由于 u 的可行域是多面体在强对偶条件下最优解会出现在 U 的顶点可以枚举顶点二是把 u 离散化引入0-1辅助变量和大M法做线性化。第二种方法在 Matlab YALMIP 框架下非常常见虽然增加了整数变量但现代求解器处理几百个二进制变量轻轻松松而且换来的是模型结构清晰、不依赖人工推导顶点表。一旦这一层算完子问题给出的最优 u* 就是当前决策下的关键场景。把它加入主问题下一轮主问题里就会新增一组约束强迫第一阶段决策 x 在这个场景下也必须有可行的第二阶段响应。如此往复场景集合越来越大解也一步步“服软”直到上下界差满足收敛条件。3. Matlab代码实现从数学模型到可运行的迭代程序3.1 先明确微网系统长什么样动手写代码前先得把系统结构定下来。我用的标准配置是一个低压交流微网包含三类分布式电源光伏阵列装机容量 400 kW预测误差取 15%风电机组装机容量 300 kW预测误差取 20%燃气轮机最大出力 250 kW带斜坡约束和最小启停时间约束储能系统容量 600 kWh最大充放电功率 150 kW充放电效率 0.95SOC范围 20%~90%与上级电网的联络线最大交互功率 300 kW调度周期取24小时时间分辨率1小时。第一阶段变量包括燃气轮机启停状态、储能各时段充电/放电状态第二阶段变量包括各机组实际出力调整量、储能调整量、联络线功率调整量等。这个配置不算大但已经能完整体现两阶段鲁棒调度的所有核心细节求解规模也适合在普通PC上跑。如果读者手头有真实的微网数据替换成本非常低只需要把预测曲线和装置参数换成自己的即可。3.2 主问题建模场景集合不断扩大的日调度计划主问题负责在给定关键场景集合的情况下求解第一阶段决策和近似成本。它长这样min c1^T x η s.t. A x ≤ b 对于每一个已经辨别的场景 u_k都有 η ≥ c2^T y_k D y_k ≤ e − E x − F u_k这里的 η 是一个辅助变量用来替代原问题中 max u∈U min y 那部分成本。每加入一个关键场景 u_k就新增一组第二阶段变量 y_k 和一约束 η ≥ c2^T y_k。这样做的好处是所有场景对应的调度问题在同一个优化模型里求解主问题的规模会随迭代次数线性增长但每一轮也更接近真实解。在 Matlab 里我用 YALMIP 建模核心代码结构如下% 定义第一阶段变量 x_bin binvar(gen_n, 24); % 燃气轮机启停 x_pg sdpvar(gen_n, 24); % 燃气轮机基准出力 x_pc sdpvar(1, 24); % 储能充电基准 x_pd sdpvar(1, 24); % 储能放电基准 eta sdpvar(1); % 定义每个阶段的场景变量和约束伪代码示意 Constraints []; Objective c1 * [x_pg(:); x_pc(:); x_pd(:)] eta; for k 1:length(scenario_set) % 取出场景k下的不确定变量取值 u_k scenario_set{k}; % 定义第二阶段的调整变量 y_delta_pg sdpvar(gen_n, 24); y_delta_soc sdpvar(1, 24); y_pgrid sdpvar(1, 24); % 添加第二阶段约束和η约束 Constraints [Constraints, eta c2 * y_delta_pg(:)]; Constraints [Constraints, ...]; % 场景k下的运行约束 end ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints, Objective, ops);这里有一个工程细节值得强调每一轮主问题里场景对应的第二阶段变量是独立的但它们共享同一个第一阶段决策 x。所以主问题的结构天然是“一个大矩阵带块对角结构”这也是为什么现代求解器处理起来并不会太吃力。随着迭代轮次增加主问题规模变大单轮求解时间也会拉长这是正常的不用慌。3.3 子问题求解怎么找到当前决策下的最恶劣场景子问题的输入是第一阶段的决策值 x_current输出是 u_worst 和最坏情况下的第二阶段成本 Q(x_current)。子问题是把内层 min 对偶之后形成的单层 max 问题。我的实际做法是用 YALMIP 直接建模对偶变量不需要手写推导让求解器的预处理器自己去处理。如果你用的是纯 Matlab 单独求解器不经过 YALMIP那就需要自己对每个约束手工对偶那个工作量相当痛苦而且容易出错。子问题里参与优化的变量包括不确定变量 u 和对偶变量 λ。需要注意原问题如果原本是 LP取对偶后的变量类型是连续非负的可以直接用 sdpvar 声明如果内层问题里出现绝对值项比如成本项里包含功率偏差绝对值那就需要先引入辅助变量 t 并添加两条不等式约束 t ≥ Δ 和 t ≥ −Δ把它改写为线性形式再交给 YALMIP 处理。我在代码里是这样组织的function [u_worst, Q_val] solve_subproblem(x_current, system_data) % 不确定变量 u_wind sdpvar(length(wind_bus), 24); u_pv sdpvar(length(pv_bus), 24); % 对偶变量按原问题约束逐条引入 lambda_balance sdpvar(1, 24); lambda_gen sdpvar(length(generators), 24); lambda_soc sdpvar(1, 24); % 不确定集约束盒式 预算 Constraints [ub_lower u_wind ub_upper]; Constraints [Constraints, sum(abs(u_wind - forecast_wind) ./ max_deviation) Gamma]; % 目标函数是已经处理好双线性项的对偶目标 Objective -lambda_balance * (x_current - system_data.forecast_load) ... lambda_gen * (system_data.capacity - x_current) ...; % 因为有双线性项需要 u 离散化或引入辅助变量处理 % 实际中我倾向于把 u 的连续盒式区间离散成有限顶点集 % 再用 big-M 线性化λ与u的乘积 ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints, -Objective, ops); u_worst value(u_wind); Q_val value(Objective); end这段代码里最值得留意的就是双线性项的线性化处理。我这边的经验是与其硬着头皮保留 λ·u 的产品项然后让非线性求解器一头撞死不如直接把 u 的可行域顶点枚举出来、配合 big-M 辅助变量做线性化。虽然问题会多出一些整数变量但换来的是求解过程的绝对稳定性和可复现性这在工程项目里比节省几十毫秒求解时间重要得多。3.4 主循环迭代直到上下界收敛整个算法的主循环不复杂按“主问题求解 → 子问题求解 → 加入场景 → 再回主问题”的顺序转圈。收敛判据是上界 UB 和下界 LB 的间隙不超过设定阈值。LB 来自主问题的最优目标值 c1^T x ηUB 来自在当前 x 下子问题求出的真实第二阶段成本 c1^T x Q(x)。两者的差反映了“场景集合不完整”造成的近似误差。当初始终止条件% 两阶段鲁棒主循环 UB inf; LB -inf; tol 1e-3; max_iter 15; scenario_set {}; for iter 1:max_iter % 1. 解主问题得到第一阶段决策和LB [x_opt, LB_new] solve_master_problem(scenario_set); LB max(LB, LB_new); % 2. 固定x解子问题得到关键场景和UB [u_worst, Q_val] solve_subproblem(x_opt, system_data); UB_new ...; % 第一阶段成本 Q_val UB min(UB, UB_new); % 3. 判断收敛 if (UB - LB) / abs(UB) tol break; end % 4. 把新的最关键场景加入集合 scenario_set [scenario_set, u_worst]; end我实测下来常规的24节点微网案例这个循环通常跑5到10轮就能收敛到1e-3的相对间隙。如果10轮还不收敛优先检查是不是子问题目标函数里少了某块成本或者不确定集模型写错导致子问题每次都在造出不可能的极端场景。4. 算例调试与参数敏感性分析4.1 基础数据配置预测曲线和不确定集参数为了做对照实验我用了一组典型日的预测数据光伏出力曲线是典型的“倒钟形”中午12点到14点达到峰值 350 kW风速预测全天在 6~9 m/s 之间波动折算成风机出力是一条不太规则的曲线负荷曲线呈现早、晚双峰特征早峰在9点前后约 500 kW晚峰在19点前后约 550 kW。电价取分时电价峰段8:00-11:00、18:00-21:001.2 元/kWh平段 0.8 元/kWh谷段 0.4 元/kWh。燃气轮机的发电成本系数取 0.45 元/kWh 加上启停成本。储能不计算折旧成本但把充放电次数作为软约束处理。不确定集参数我设成光伏出力最大偏差取预测值的 15%风电取 20%负荷偏差取 5%。预算 Γ 先取 6相当于不确定变量总数的一半左右后面再专门做敏感性分析。4.2 预算 Γ 对调度经济性和保守性的影响把 Γ 从 0 逐步调到 12其他参数保持不变跑完整个两阶段鲁棒流程后调度总成本的变化趋势非常能说明问题Γ 值总成本元与确定性模型成本差最坏场景辨识迭代次数0确定性432601244893.8%4446126.6%6647209.1%88481111.2%1012全盒式495614.6%12从这个表里能明显看到一个规律成本随着 Γ 增加而上升但边际增幅在递减。Γ 从0到4成本增加6.6%从4到8增加4.6%从8到12只增加3.4%。这说明预算约束的设计价值就在这里——用比较小的经济性牺牲覆盖掉绝大部分极端场景风险。在真实项目中我不会盲目把 Γ 拉到满格而是结合历史极端天气出现的概率来选一般取不确定变量总数的1/3到1/2是性价比最高的区间。另外我注意到 Γ 越大迭代收敛需要的场景数量越多。原因也很好理解Γ 增大了 U 的“搜索空间”子问题需要更多轮次才能逼近全局最坏场景。这一点在算法性能调试时值得提前有心理预期。4.3 不确定区间宽度的影响警惕过度保守除了 Γ不确定区间的宽度同样关键。我把光伏、风电的预测误差上限从5%逐步加到30%固定 Γ6结果如下误差上限总成本元最恶劣场景下弃风弃光率5%45121.2%10%45983.8%15%47206.7%20%48909.5%30%524615.2%可以看到当误差上限超过某个阈值之后成本会快速抬升。误差上限取30%的时候成本比5%时高了整整16%而且最恶劣场景下的弃风弃光率也到了15%左右——说明这个方案已经开始靠牺牲清洁能源消纳来保安全。实际操作中我建议不确定区间的宽度要从历史预测误差数据里统计出来而不是拍脑袋定一个数。对每一预测时段统计过去N天同一时刻的预测偏差取95%分位数作为该时段的上下界比统一按百分比更贴近真实波动规律而且调度结果往往更经济。这些统计代码在 Matlab 里无非就是 quantile 函数几行的事性价比极高。5. 常见问题与排查实录5.1 子问题是双层结构直接求解失败怎么办我早期写过一版代码子问题直接用 YALMIP 的 max-min 嵌套去建模结果求解器一直报“无法处理嵌套优化问题”。这是因为标准求解器不支持一个优化问题里再套一个优化问题。解决办法就是把内层 min 对偶成 max。如果你不想手推对偶YALMIP 里可以用 duality 相关的工具函数辅助但更稳妥的是针对自己的原问题手推一遍。我把这个步骤遇到过的问题汇总了一下对偶变量非负性漏写很多约束对应到对偶变量时是“无符号限制”的比如等式约束对应的对偶变量是自由变量unrestricted in sign写错会让对偶问题的最优值跟原问题完全对不上。原问题如果是极小化对偶问题的目标函数符号容易搞反。检查对偶间隙是否为零的有效方法把原问题的原始最优解代入对偶目标函数验证是否同值。绝对值项忘记线性化成本中有 |y| 形式的项时直接进对偶会得到非线性项必须先用辅助变量改写。建议做法是先拿一个小规模的算例手算一遍确认对偶变换完全正确之后再上整套砂型代码。正式解决前先用尽量能手工验算的算例调试对偶逻辑这步省下来的排查时间通常非常多。5.2 优化结果不可行多半是第二阶段子问题缺了松弛变量两阶段鲁棒的子问题有时候会遇到不可行的情况。原因可能不只是模型写错——搜索到的那个极端场景在当前 x 下确实让系统没有任何可行的 y 可以满足所有约束。这在实际物理系统里意味着你第一阶段定的基准计划太“僵”遇到了远超预期的不确定性实现。标准做法是在第二阶段约束里引入松弛变量给关键约束比如功率平衡、SOC上下限配上很大的惩罚系数。相当于给紧急情况留了一扇“紧急通风口”并且用高惩罚成本逼着优化器尽量避免去开它。这样即使极端场景下系统不能严格满足所有硬约束模型也不会直接崩掉而是反应成一个可负担的代价。惩罚系数的选择也要注意。我习惯先把惩罚成本系数设成正常成本的 100 倍以上然后跑一遍看它是否被激活。如果发现最恶劣场景下松弛变量显著非零那说明方案本身已经有了物理上的不可行风险。这时候一定要回头检查第一阶段决策而不是一味调大惩罚系数把问题掩盖掉。5.3 Matlab 环境配置与求解器选型Gurobi 还是 Cplex这套代码依赖 YALMIP 工具箱和一个高性能混合整数线性规划求解器不是 Matlab 自带的 solve 函数能搞定的。我的实际配置是Matlab R2023bYALMIP 最新版GitHub 上直接 clone 下来路径加进 Matlab 搜索目录Gurobi 求解器学生或学术版免费性能和数值稳定性都比 Cplex 好尤其是处理 big-M 线性化带来的整数变量时用 sdpsettings(solver, gurobi) 指定求解器如果你用的是 Cplex代码结构基本不用改YALMIP 后端自动识别。不过我个人在双线性项线性化后的问题里Gurobi 的整数预处理能力明显更强求解速度能差到一倍以上。如果没有 Gurobi 授权Cplex 的学术版也完全可以跑通整个流程只是大算例下会慢一些。另外如果系统里没装求解器YALMIP 内置的免费求解器比如 sedumi、linprog也能跑但只能处理连续线性规划遇到 big-M 引入的二进制变量会直接报错。所以一开始装环境的时候把 Gurobi 或 Cplex 装好是最重要的一步。5.4 我踩过的一些坑一些程序运行层面的实际问题第一个坑是 SOC 约束跨时段耦合。储能SOC 是逐时累积的变量它在第二阶段场景里也要跟着调整。如果每个场景都独立定义一套 SOC 变量而场景之间的 SOC 没有共享或衔接关系那么即便每个场景单独都可行整体能量约束也可能是虚的。我的处理方式是把 SOC 初值设为固定参数并且保证主问题里的储能基准计划在所有场景下都不越界。这样做虽然会损失一点储能灵活性但换来的是模型在数学上更干净、收敛更稳。第二个坑是 big-M 值的选择。线性化 λ·u 时M 值取得太大会导致数值条件数恶化求解器给出的解精度漂移取得太小又可能把可行域错误地切掉。通用的做法是先用一个不包含二元变量的 LP 求解一次估算 λ 的合理范围然后让 M 比这个范围上限大一个数量级左右不要无脑写 1e6。第三个坑是收敛判据里的绝对值。当两阶段成本里有大量固定成本时UB 和 LB 的差距本身就有一个“平台期”直接看绝对间隙可能永远达不到 1e-3需要改成相对间隙除以UB作为收敛条件。我在前面主循环代码里写的就是相对间隙判据这个细节对实际跑通流程很关键。最后再做一次实验记录上面的例子在 i5 16G 内存的普通办公笔记本上单轮主问题 子问题的求解时间在 7 到 15 秒之间完整迭代 8 轮收敛到 0.1% 相对间隙总耗时约 2 分钟。这个计算成本在微网日前调度场景下完全能接受甚至可以在滚动更新的场景下多跑几轮。如果你在自己机器上跑出了远超这个数量级的耗时优先检查一下是不是场景数增长导致主问题规模爆炸或者是否不小心把第二阶段变量不该有的耦合写进去了。做两阶段鲁棒优化我最想跟大家强调的一点是不要一上来就求“完美复现”。先把确定性模型跑通再把不确定性集合加进去最后才上主问题-子问题的迭代框架。每一步拆开都能验证最后合起来就不会出现“跑不通但不知道哪里错了”的焦灼。这个从简到繁的过程也是我重写了三版代码后最深刻的体会。