企业官网建设流程全解析

晚上七点半，微网里的风电出力开始往下掉，负荷却往上蹿。我当时用的还是那套经典的确定性日前调度，24小时计划排得满满当当，结果一遇到这种预测偏差，只能临时手动切燃气轮机、拉储能，折腾下来成本超预算15%，还差点甩掉一片负荷。后来我把调度框架换成两阶段鲁棒微网优化调度，配合关键场景辨别算法去刻画不确定性，才把这类风险真正压住。这篇东西我不打算写成教材，就把我在Matlab里从建模、求解到调通的全过程，连同踩过的坑，一次性说清楚。文章面向正在研究微电网优化调度、两阶段鲁棒优化，或者拿Matlab做论文仿真但被C&CG迭代卡到怀疑人生的同学，也适合想弄清楚"鲁棒优化到底怎么落地"的工程师。

1. 为什么微网调度要搞"两阶段"还要"鲁棒"：从一次晚间负荷突变聊起的问题

1.1 确定性调度的短板：预测误差下为什么总在"救火"

先看传统做法。我们把微网里的风电、光伏、负荷都当成已知参数，建一个混合整数线性规划（MILP），求解让总运行成本最小的出力计划。这个模型本身没毛病，而且在预测完全准确的时候结果也确实漂亮。问题在于，真实世界里风电出力是波动的，负荷预测总有偏差，光伏更是每分钟都可能被云层遮住。一旦实际值和预测值对不上，日前排好的计划就变成了纸面方案。

我举一个实际案例。某天晚间19点到22点，天气预报的风速比实际偏高40%，光伏在18点半后已经归零，负荷却因为突然降温比预测高出8%。按照确定性模型，最优方案是让燃气轮机在低谷时段多发电、储能低谷充电高峰放电。但实际运行时，燃气轮机爬坡速率跟不上负荷突增，储能又在高峰段被计划排满，最后只能紧急调用柴油机补电——单小时成本翻了快三倍。这不是模型算法错，而是模型本身假设了一个不会出错的世界。

业内处理不确定性的主流路线有两条：随机优化（stochastic programming）和鲁棒优化（robust optimization）。随机优化要给不确定性变量配概率分布，然后搞场景树，场景一多计算量爆炸，而且微网里风电、光伏、负荷的联合分布你根本拿不到准确数据。鲁棒优化的思路更"蛮横"：我不赌概率，我把所有可能的不确定性范围圈起来，在这个范围里找出最坏的情况，然后保证这个最坏情况下系统也能安全运行且成本可控。这就是"两阶段鲁棒微网优化调度"的核心动机。

1.2 两阶段决策的真实含义：日前计划与实时调整的分工

"两阶段"不是玄学，它对应的是微网运行中真实存在的两个决策层次。

第一阶段（这里统一叫主问题决策）发生在日前，做的是那些"定了就很难改"的决定：哪些燃气轮机明天开机、哪些停机，储能的充放电状态怎么安排，与上级电网的购电协议要不要签。这些决策的特点是带二进制变量，要么0要么1，一旦定了之后日内一般不轻易翻转。第二阶段（子问题决策）发生在日内实时，做的是"在开机状态和储能状态都锁死的前提下，各台机组实际发多少电、储能实际充放多少功率、从电网买多少电"。这一层的变量是连续的，类似经济调度的出力分配。

这样的分层其实很符合微网的实际操作习惯：日前做机组组合，日内做经济调度。两阶段鲁棒模型把不确定性放进了第二阶段——假设第一阶段决策已经做出，不确定性在第二阶段显现，系统做最坏情况下的最低成本调整。我们优化的目标是让"第一阶段成本 + 第二阶段最坏场景下的最低运行成本"之和最小，这就是典型的min-max-min结构。

1.3 鲁棒优化和随机优化的差别：不需要精确概率分布

接着上面说的，随机优化的结果是一个期望成本，它隐含假设你知道每个场景的概率。鲁棒优化不需要这个假设，它只需要一个"不确定集"——就是不确定性变量可能的取值范围。比如风电出力的预测值是50kW，允许偏差±10%，那不确定集就是一个盒式区间。你在这个区间里找最坏情况，然后保证它发生系统也不出事。

代价是鲁棒优化通常比随机优化保守，毕竟它只盯着最坏场景。但如果不确定集构造合理，比如加一个预算约束Γ，限制"最坏情况不会在所有时段同时发生"，保守程度会被大幅压缩。这个预算约束我用下来是对工程实用性影响最大的参数之一，后面会专门讲怎么调。

2. 两阶段鲁棒优化的数学骨架：min-max-min结构、不确定集与对偶转换

2.1 目标函数和约束怎么拆成"主问题+子问题"

两阶段鲁棒的标准形式，我用代码块把它表达出来，因为纯文本写公式容易糊：

目标函数： min_x ( c^T * x + max_{u ∈ U} min_{y ∈ Ω(x, u)} d^T * y ) 第一阶段约束： A * x ≥ b x ∈ {0,1}^p × R^n 第二阶段约束（对给定x和u）： Ω(x, u) = { y | B * y ≥ h - C * x - D * u, y ≥ 0 }

这里x代表第一阶段决策，u是不确定性向量，y是第二阶段连续决策。整个问题的难点在于中间那个max-min：对于每个x，你得先找到让y最优目标值最大的那个u（最恶劣场景），再算这个场景下y的最小成本。这在计算上不是一个可以直接丢给求解器的问题，得拆。

最常用的拆法是列约束生成算法（C&CG，Column-and-Constraint Generation）。C&CG的套路是：先不管u，只针对一个初始场景u1求解第一阶段问题，得到一组x；把x代入第二阶段，求解max-min子问题，找到最恶劣场景u*；如果这个u对应的目标值还没达到收敛条件，就把关于u的约束和变量加到主问题里，重新求解x。如此反复迭代，直到上界和下界之差小于容忍度。

2.2 不确定集的选择：盒式、椭球式与预算约束

不确定集决定了鲁棒优化的"防守范围"。我常用的有三种。

最基础的是盒式不确定集，形式是u ∈ [u_mean - Δu, u_mean + Δu]，每个时刻的偏差独立，最大偏差Δu根据历史预测误差统计得到。优点是建模简单，缺点是太保守——如果风电和负荷在24个时段都取到最坏偏差，那运行成本会高得离谱，工程上根本接受不了。

改进方案是加预算约束Γ。还是盒式集合，但限制偏差的总量不超过一个阈值。比如Γ = 4，含义是"整个调度周期内，至多4个时段的不确定性同时取到边界值"。这样既保留了对极端情况的防护，又不会让所有时段一起背锅。Γ的值一般取调度时段数的15%到30%，我系统测试下来，Γ取5（24时段系统）时成本比Γ取24要低不少，且失负荷概率仍然很低。

椭球不确定集更平滑，但会引入二阶锥约束，求解负担变大。我的经验是：工程调度用盒式+预算约束性价比最高，除非你的微网里储能占比很大、需要精细刻画相关性，才值得上椭球。

2.3 子问题里max-min怎么通过强对偶变成可求解的单层问题

第二部分子问题是个双层优化：内层对y求min，外层对u求max。YALMIP不能直接处理这种嵌套结构，所以得把它转化成单层问题。

核心工具是线性规划的强对偶定理，或者更通用的KKT条件。因为内层对y的线性规划是凸的，满足强对偶条件，我可以把内层的min_y问题替换成它的对偶问题max_λ，于是max-min变成：

max_{u∈U} max_{λ ≥ 0} λ^T * (h - C*x - D*u)

把两个max合并，就变成一个在(u, λ)上联合求max的线性问题，可以直接丢给Gurobi或CPLEX求解。这里有三个关键点：

第一，这要求第二阶段问题里y的系数矩阵不包含二进制变量，否则强对偶不成立。如果模型里第二阶段也有启停变量或者储能状态变量，你得先用大M法或者把它提升到第一阶段加以线性化，否则对偶间隙会直接弄崩收敛精度。这点我在第5章会重点展开。

第二，对偶变量λ本身是无约束的，但要保证对偶问题有界，需要加上内层原始问题约束对应的对偶可行域约束，写代码时别忘了把对偶变量符号约束对应好。

第三，子问题求出的最优u*就是当前x下的最恶劣场景，它的目标函数值加上第一阶段成本，构成原始问题的上界（UB）。主问题每次迭代求解得到的x和对应目标值，因为约束少于原始问题，给出的则是下界（LB）。迭代就是不断缩小UB-LB的间隙。

3. 关键场景辨别算法：它不是在"挑场景"，而是在给C&CG做瘦身

3.1 C&CG列约束生成的标准流程

我先描述一下C&CG的标准循环，大家都熟悉：

1. 初始化：把不确定变量u固定到某个初始场景，比如预测平均值，求解主问题MP，得到x和LB的候选值。
2. 把x代入子问题SP，求解max-min，得到最恶劣场景u和SP的目标值。此时UB = c^Tx + SP目标值。
3. 若UB - LB < ε，停止；否则把u*对应的第二阶段变量及其约束加入MP，回到第1步。

逻辑上很简单。它最大的优势是每次迭代加入的是"整段约束+变量"，比经典的Benders分解只加一条割线要紧得多，所以收敛次数通常远少于Benders。我用下来，大多数24时段系统，C&CG在5到15次迭代内能收敛，而Benders往往需要几十上百次迭代。

3.2 为什么切割会越加越多、主问题越解越慢

C&CG有个很实际的问题：每轮子问题都会产出一个最恶劣场景，这些场景会被全部加进主问题，主问题的规模会随迭代次数线性膨胀。你听上去每次只多一个场景，好像无所谓，但主问题是含有二进制变量和几百条约束的MILP，每多一个场景，就是几十个变量和几十条约束。迭代到第八轮、第十轮的时候，MP的求解时间会从几秒变成几分钟，整个算法在算力上就"僵住"了。

更麻烦的是，很多场景其实长得差不多。比如风电出力最恶劣的情况，可能连续几轮迭代解出的u*都是"第7时段和第13时段取上界"，只是其他时段的微调不同。这些场景对约束的贡献高度重叠，把它们全部加进MP，等于把大量冗余约束堆进模型，白花钱。

关键场景辨别算法要解决的就是这个问题。它的思路很直接：在把子问题解出的场景加入MP之前，先做一次"筛选体检"，只放那些对解有实质影响的新鲜场景进去。

3.3 关键场景辨别的判据设计与场景库管理

我用的判据分为三个层面：

第一层是场景重复度。维护一个场景库S，里面存着所有已经加入MP的场景。每轮SP解出新的u*，先计算它与场景库里每个场景的归一化欧氏距离：

dist = ||u* - s_i|| / ||s_i||

如果最小距离小于一个阈值ε_s（我通常取0.05），就认为这个场景和库里的某个场景在结构上几乎相同，直接丢弃，不更新MP。这能杜绝大量"换汤不换药"的重复切割。

第二层是目标值增益。即使场景不重复，也要看它对这个x到底"狠不狠"。计算该场景下SP的目标值（也就是如果加入该场景会给MP带来的约束强度）与当前LB之间的相对提升幅度。如果提升幅度小于0.5%，说明这个场景虽然新，但对优化方向没有实质影响，可以延迟或者丢弃。

第三层是紧约束检测。这一步更精细：SP求解完会给出对偶乘子，乘子里数值较大的那部分意味着当前x在该场景下"卡得很紧"。如果一个场景对应的最大对偶乘子都很小，说明即使这个场景发生，第一阶段决策也毫不费力，那它就不构成有效约束。用这个指标可以进一步过滤。

我最后实现的流程是：SP解出u* → 重复度检查 → 目标增益检查 → 紧约束检查 → 决定是否加入MP和场景库。三重门槛全过，才更新MP。实测在24时段微网算例中，标准C&GG需要12轮、往MP加了12个场景；加了辨别后只加了5个场景就收敛，MP单轮求解时间从峰值17秒降到6秒，总计算时间大概缩短了40%。

3.4 算法收敛与计算时间对比：一个简单算例的直观感受

简单看一组来自我算例测试的对比数据（微网含1台燃气轮机、1套储能、风电场、负荷，24时段，Gurobi 9.5作为求解器）：

两条路线最终目标值几乎一样，说明辨别算法没有牺牲精度，但把计算时间砍掉了一半多。对小规模系统你体感不强，一旦微网节点数上百、或者把时序运行耦合约束加进来，这个差距会被放大到小时级别。

4. Matlab代码实现的关键环节：YALMIP建模、MP/SP迭代与场景库动态更新

4.1 工作环境与求解器配置

我的环境是Matlab R2022b + YALMIP + Gurobi。YALMIP负责把优化问题建模成高层的数学表达式，Gurobi负责实际求解MILP和LP。如果你手头没有Gurobi，CPLEX也行，甚至Mosek都能跑大部分步骤；不过子问题SP里那个大M法线性化部分，Mosek处理起来会比Gurobi略慢，我测过，不算大问题。

安装之后先跑一下自检：

yalmiptest;

确保Gurobi和YALMIP之间通信正常。一个容易踩的坑是Gurobi版本和YALMIP不兼容，表现是solve函数报Output too long或者Unknown solver。我建议Gurobi 9.x配合YALMIP R20200516以上的版本，别用太旧的YALMIP。

4.2 主问题MP的建模与切割追加方式

主问题MP本质上是一个带离散变量的大MILP。建模代码如下（省略具体参数赋值，思路是结构性的）：

% —— 主问题MP建模 —— x_s = binvar(n_gen, T); % 机组开停状态 p_g = sdpvar(n_gen, T); % 各机组出力 p_es = sdpvar(2, T); % 储能充放功率 cost_second = sdpvar(1); % 第二阶段成本变量 % 目标：第一阶段成本 + 第二阶段成本变量 objective = sum(sum(cost_start * max(0, diff(x_s,1,2)))) + ... sum(sum(c_gen .* p_g)) + sum(p_es(1,:).*c_ch + p_es(2,:).*c_dis) + cost_second; Constraints = []; Constraints = [Constraints, % 机组出力上下限 p_g <= repmat(P_G_MAX',1,T) .* x_s, ... p_g >= repmat(P_G_MIN',1,T) .* x_s]; % 储能约束 Constraints = [Constraints, SOC(:,2:end) == SOC(:,1:end-1) + p_es(1,:)*eta_ch - p_es(2,:)/eta_dis]; Constraints = [Constraints, SOC(:,1) == SOC_init, SOC >= SOC_min, SOC <= SOC_max];

关键的切割追加逻辑在迭代循环里。每次SP解出一个新场景u_k，我就往MP里加一组对应这个场景的等式约束、不等式约束和变量：

p_g_k = sdpvar(n_gen, T); % 该场景下的机组出力 p_es_k = sdpvar(2, T); % 该场景下的储能功率 p_buy_k = sdpvar(1, T); % 该场景下的购电功率 Constraints = [Constraints, % 功率平衡：负荷 = 各类电源出力 + 储能放电 - 充电 + 购电 p_load - u_k .* p_wind == sum(p_g_k,1) + p_es_k(2,:) - p_es_k(1,:) + p_buy_k, p_g_k <= repmat(P_G_MAX',1,T) .* x_s, p_g_k >= repmat(P_G_MIN',1,T) .* x_s, % …… 其余约束类似 ];

这里有个坑：MP里新加的第二阶段变量p_g_k等，必须和已有的x_s等第一阶段变量通过约束关联起来，否则它只是被"塞进"模型，对x_s产生不了约束力，切割等于白加。

4.3 子问题SP的KKT/强对偶转换在Matlab里的落地

子问题SP求解是整篇代码里最容易写错的地方。我的做法是：先写出内层规划的对偶问题，然后用对偶变量作为新的优化变量，和u一起做联合求解。

如果内层问题是纯LP，对偶问题可以手推。比如内层问题（给定x和u）是：

min_y d^T y s.t. B y ≥ h - C x - D u , y ≥ 0

其对偶是：

max_λ λ^T (h - C x - D u) s.t. B^T λ ≤ d λ ≥ 0

于是SP变成一个在(u, λ)上的联合LP：

max_{u, λ} λ^T (h - C x - D u) s.t. B^T λ ≤ d λ ≥ 0 u ∈ U

这个形式在Matlab里直接用YALMIP写就行：

% —— 子问题SP（强对偶转换后） —— lambda = sdpvar(n_con, 1); % 对偶变量 u_var = sdpvar(n_unc, 1); % 不确定变量 % 不确定集：盒式+预算约束 Constraints_SP = [Constraints_SP, u_mean - u_dev <= u_var <= u_mean + u_dev, sum(abs(u_var - u_mean) ./ u_dev) <= Gamma]; % 对偶可行域 Constraints_SP = [Constraints_SP, B' * lambda <= d, lambda >= 0]; % 目标：max λ^T (h - C*x_fixed - D*u_var) Objective_SP = -lambda' * (h - C*x_fixed - D*u_var); % 注意YALMIP是求min，所以取负号 optimize(Constraints_SP, Objective_SP, sdpsettings('solver','gurobi'));

如果你不想手推对偶，也可以用KKT条件把内层问题作为约束放进SP，配合互补松弛的线性化。这个做法不用求对偶表达式，但会引入更多变量和非线性项，我一般只在内层问题比较复杂、手推对偶容易出错时才用。

4.4 收敛判据、初始场景选择与关键代码框架

C&CG的收敛判据是UB-LB小于容忍度ε。UB是SP目标值加上第一阶段成本，LB是MP目标值。迭代主体我用伪代码形式放出来：

% 初始化 x_current = []; UB = inf; LB = -inf; k = 0; quiesce = false; while (UB - LB) > epsilon && ~quiesce k = k + 1; % 求解MP optimize(MP_Constraints, MP_Objective, MP_options); LB = value(MP_Objective); x_current = value([x_g, x_es]); % 取出第一阶段决策 % 求解SP optimize(SP_Constraints, SP_Objective, SP_options); SP_obj = -value(SP_Objective); % 还原符号 UB_temp = value(objective_first) + SP_obj; if UB_temp < UB UB = UB_temp; end u_star = value(u_var); % 最恶劣场景 % —— 关键场景辨别 —— if ~isKeyScene(u_star, scene_set, epsilon_s) || ... (UB_temp - LB) / LB < epsilon_gain % 不加入MP，继续下一轮 quiesce = true; % 或者进入更精细的局部搜索 continue; end % 加入新场景和对应变量约束到MP scene_set = [scene_set, u_star]; end

初始场景怎么选对收敛速度影响很大。我的经验是用预测均值场景，也就是u = u_mean，先让系统在"最乐观"的情况下给出一个初始x。这个阶段MP求解很快，得到的x比较平滑。之后SP再去找这个x下的最恶劣场景，逐轮逼近。比随机初始化收敛稳定得多。

还有一个细节：SP问题规模不大，但有时会因为对偶变量过多而变慢。可以预先在MP里把所有对偶变量的索引整理成一个稀疏矩阵传入Gurobi，别让YALMIP动态扩维。

5. 我在调试中踩过的坑：二进制变量、对偶间隙与"伪收敛"

5.1 子问题里的二进制变量怎么处理：大M法与二进制展开

这是新手最容易卡死的地方。如果你的微网模型里，储能充放电状态、机组启停状态这些二进制决策出现在第二阶段子问题中——比如实时调度还要决定储能是充还是放——那么内层min问题就不是一个LP，而是MILP，强对偶定理不成立，你没法直接转换成单层max-λ形式。

我踩过一次很惨的坑：直接把MILP子问题扔进C&CG，用Gurobi的默认解算，UB和LB怎么都收敛不到同一个值，差了大概2.7%。排查后发现是SP里含二进制变量，对偶转换出来的对偶间隙根本绕不过去。

解决办法有两种：

第一种是大M法。把第二阶段里的二进制变量全部线性化掉，比如储能充放电的互斥约束p_ch + p_dis ≤ M * z_state（z_state是0/1变量），然后再做对偶转换时把这部分约束换成它的对偶约束。这个大M要选得足够大但不至于破坏数值稳定性，M取该机组最大出力或储能最大功率的3到5倍，我一般这样定。

第二种是二进制展开法。把第二阶段二进制变量从子问题中"提出来"，放到第一阶段去，或者在每次迭代求解子问题时，把x_fixed当作参数、枚举部分典型组合来生成场景。这种方法在二进制变量少的时候可行，变量一多就会让场景数量爆炸。

我的最终建议是：尽量重新审视模型划分，把可能改变系统拓扑结构的离散决策（机组启停、储能状态）统一放到第一阶段，第二阶段只保留连续变量的经济分配。这样两阶段鲁棒模型才"干净"，也最符合微网实际运行逻辑——日前确定运行方式，日内只做连续功率调整。

5.2 为什么UB-LB卡住不动：切割失效的典型原因

我调试C&CG时遇到过最诡异的现象是：迭代三轮之后UB和LB完全卡住，两个数字在小数点后两位都不变化了。表面上看是收敛，但实际上UB比真实最优高了一大截，是典型的"伪收敛"。

后来排查发现，问题出在MP里追加的第二阶段变量没有和第一阶段变量形成足够强的耦合约束。具体来说，我在第4.2节写的那个功率平衡等式里，把u_k参与的地方写错成常量了——风电出力直接用了预测值，导致场景变化根本不影响MP约束结构。SP每次返回的u*虽然不同，但MP根本不"感知"这些场景，所以LB一直停在同一个值上。

这是个非常隐蔽的bug。排查方法很简单：看每个场景加入后LB有没有实质变化。如果连续两个不同的极端场景加进去、LB纹丝不动，那几乎可以确定切割约束是失效的，赶紧检查追加的约束和变量是否正确关联了x_s。

另一个常见的原因是场景重复太严重。SP每次返回的u*只在两个数值上变动，MP的信息量早就饱和了。这种就要靠关键场景辨别里的重复度检查来提前拦截，否则白白多迭代好几轮。

5.3 关键场景辨别阈值怎么调：太小没效果、太大丢精度

阈值ε_s和ε_gain过小，辨别算法形同虚设，场景还是会大量涌入MP；过大又会把真正紧约束的场景过滤掉，导致UB-LB虽然收敛快，但收敛到的解严重偏离真实最优。

我的标定经验是：用历史几轮SP的场景轨迹做统计。先跑一轮完整的标准C&CG，把每轮场景的归一化欧氏距离和目标值增益记录下来。取中位数和分位数作为初始阈值。对我来说，ε_s取0.03到0.08之间比较合适，ε_gain取0.002到0.01之间（0.5%左右）。然后对比同一算例在不同阈值下的最终UB，挑UB相对变化不超过0.1%的最大阈值。这样既保证精度又保证速度。

这里放一个测试小表：

第三组参数明显丢精度了，所以阈值不是越大越好，得卡在中间区间。

5.4 求解器数值问题：YALMIP+Gurobi的精度参数注意点

鲁棒调度模型里经常出现量级差异巨大的系数，比如成本系数是几百，设备容量是几万千瓦，对偶变量可能小到10^-4，这种数值差异会让Gurobi内部单纯形法出现数值警告。

我建议在sdpsettings里显式设置几个参数：

opts = sdpsettings('solver','gurobi', ... 'gurobi.MIPGap', 1e-4, ... 'gurobi.OptimalityTol', 1e-6, ... 'gurobi.FeasibilityTol', 1e-6, ... 'gurobi.NumericFocus', 2);

更重要的预处理是：把成本除以100或1000、功率以MW为单位、成本系数单位统一到元/MW，避免出现10^6级别的巨大数值。我遇到过一次MP求解器返回Success但结果明显不对，打开Gurobi日志发现模型有数值警告，就是因为风电和负荷功率用了W为单位，数值太大导致Presolve直接把一堆约束误删了。这种错误极难发现，只能靠这个经验提前规避。

6. 拿到调度方案之后：怎么用蒙特卡洛验证它真的"鲁棒"

6.1 后验评估流程：别让鲁棒优化变成"纸上盾牌"

优化求解收敛不等于方案可信。我每次跑完两阶段鲁棒，一定会再做一次蒙特卡洛后验评估：根据风电和负荷预测误差的历史分布，随机生成2000组不确定性场景（注意这些场景可以超出你设定的不确定集范围，用来测试方案的"安全感"），然后把鲁棒调度方案的第一阶段决策x固定下来，逐个场景求解第二阶段经济调度问题。

评估指标我常用三个：系统失负荷概率（LOLP）、缺供电量期望值（EENS）、平均运行成本。鲁棒方案和确定性方案后验对比通常长这样：

这个表很有说服力。确定性方案平均成本低，但在极端场景下会付出高额代价，甚至失负荷；鲁棒方案平均成本高出5%，但代价的"尾部"被压住了。实际工程里，尾部风险往往比平均值更致命——一次失负荷事故的惩罚成本可能顶得上几十次正常运行的利润。

6.2 蒙特卡洛实现的几个注意点

做后验评估的时候，我直接把不确定性场景当作参数传入第二阶段问题（不经过SP那个max-min结构），所以它就是一个标准的LP或QP，求解极快。2000个场景并行跑，Matlab里用parfor可以轻松在三分钟内完成。

这里也有个容易踩的坑：第二阶段问题本身可能无解。如果某个极端场景下燃气轮机爬坡不够、储能电量不足，系统会直接不可行。在确定性优化和鲁棒优化的主问题里，我们一般允许切负荷并设定很高的罚函数成本，让模型有"最后保底手段"。后验评估时保留同样的罚函数，才能让每个场景都有解，LOLP和EENS才有意义。

6.3 我在实际运行中积累的几点体会

这套流程跑了快一年，我个人最深的感触是：两阶段鲁棒优化的价值不在它算得有多快多准，而在于它逼着你在建模阶段就把"哪些决策是日前定的、哪些是日内调的、不确定性怎么圈"这三件事想清楚。这个梳理过程本身，比跑通代码更值钱。

另外，预算约束Γ的取值千万别照搬文献。不同微网的负荷曲线、风电置信度、储能容量差异很大，我建议把它当成一个可调参数做敏感性分析，画一条"Γ vs 成本/失负荷概率"曲线，再让运行人员根据风险偏好去选。成本端和风险端的平衡，只有结合具体场景才有意义。

最后再分享一个小技巧：C&CG迭代过程中，把每轮SP求出的最恶劣场景u*存下来做个热力图，能非常直观地看出哪些时段对不确定性最敏感。我接手过一套系统，迭代里反复出现的就是傍晚骤风时刻和晚间负荷尖峰，后来根据这个观察，专门在傍晚时段提高了燃气轮机的热备用，效果比单纯调大Γ还要好。这个套路大家以后做鲁棒优化分析时可以试试。