企业官网建设流程全解析

做电力系统调度优化的朋友，应该都体会过那种“模型算出来挺漂亮，一到实际就跑偏”的无力感。风电预测给的是100MW，实际出力可能只有60MW；光伏中午预报大晴天，下午一片云飘过来出力直接腰斩；负荷就更不用说了，居民区和工业区的用电习惯根本不是一条平滑曲线。传统的确定性调度模型把这些预测值当成板上钉钉的常数，算出来的机组组合、备用容量、潮流分布，一旦遇到风光出力波动，轻则经济性大幅恶化，重则出现切负荷甚至潮流越限。这也是为什么近些年鲁棒优化在电力系统调度里越来越受重视——它不赌预测值的准确性，而是把不确定性“圈”在一个集合里，保证集合内所有可能场景下方案都可行。

这篇文章想分享的，就是一套我实际跑通的计及风、光、负荷不确定性的两阶段鲁棒优化完整实现，建模用YALMIP，求解器用Gurobi，核心算法是C&CG（列与约束生成），中间的双线性项用大M法做线性化。内容会从问题建模一路讲到Matlab代码细节和踩坑经验，适合正在做微电网调度、园区综合能源优化、电力系统日前计划这类课题的研究生和工程师参考。

1. 为什么确定性调度到了含新能源的系统里就不够用了

1.1 确定性模型的“赌徒心态”

传统经济调度模型里，风电、光伏、负荷都被写成了确定数值。模型长这样：

目标是最小化总运行成本，约束包括功率平衡、机组出力上下限、爬坡速率、网络安全等。求解器给出的结果是一个“在预测值完全准确时才最优”的方案。

问题在于，新能源的预测误差不是小概率事件，而是常态。风电短期预测误差的均方根误差普遍在10%~20%之间，极端天气下甚至能到30%以上。如果日前计划按预测值100MW安排火电出力，实际风电场发到60MW，那40MW的功率缺口就得靠实时市场的备用、储能或者紧急爬坡来补。要是备用容量不够，就得限负荷。

这种“赌预测准”的思路，在传统火电为主的系统里问题不大，因为负荷预测相对精准，机组调节能力也强。但风光占比提高之后，系统的不确定性来源从一个变成了三个，预测误差的叠加效应让确定性方案的实际运行成本远高于模型计算值。

1.2 鲁棒优化的核心思想：不赌，但做好准备

鲁棒优化的思路完全不同。它不依赖预测值的概率分布，而是给每个不确定参数划定一个变化范围，构造一个不确定集合（Uncertainty Set），然后寻找一个对所有集合内场景都可行的方案。

说句大白话：确定性优化是“我赌明天不下雨，所以只带一把伞”；鲁棒优化是“我不知道明天下不下雨，所以我包里常备雨衣，真下大了也不慌”。代价是雨衣占地方、有点重——对应到调度里，就是鲁棒方案的经济性会比确定性方案差一些，多出来的成本就是为了换取“不管风光怎么波动，系统都安全运行”的保障。

这个保守度是可以调节的。通过预算约束（Budget Constraint），我们可以限制“最坏情况下同时偏移到极端的不确定源数量”。比如三个不确定源，预算Γ取1，意味着优化只保证最多一个源同时取到极端值；Γ取3，则要求所有源同时取极端值时方案依然可行。Γ越大，方案越保守，成本越高。

1.3 两阶段结构：先定下来，再调整

两阶段鲁棒优化对应的是电力系统调度的实际时序逻辑：

• 第一阶段（日前决策）：机组启停、备用容量预留、储能充放电计划等必须提前确定下来。这些决策变量在不确定参数实现之前就要敲定，特点是“定了就改不了”。
• 第二阶段（实时调整）：等风光出力、负荷的实际值揭晓后，系统在给定第一阶段决策的基础上进行再调度，通过调整火电出力、储能充放电、切负荷等手段，以最小代价保证功率平衡和网络安全。

所以模型的目标函数由两部分构成：第一阶段的启动成本和基础运行成本，加上第二阶段在最坏场景下的调整成本。外层是第一阶段决策的min问题，内层是一个max-min嵌套——先枚举所有可能的不确定场景（max），再求每个场景下的最小再调度成本（min）。

2. 两阶段鲁棒优化的数学模型拆解

2.1 不确定集合的构造与离散化

风、光、负荷的不确定集合，最常用的是盒式集合加预算约束。先把每个不确定源写成名义值加偏差的形式：

P_w ∈ [P_w^0 − ΔP_w, P_w^0 + ΔP_w] P_v ∈ [P_v^0 − ΔP_v, P_v^0 + ΔP_v] P_l ∈ [P_l^0 − ΔP_l, P_l^0 + ΔP_l]

其中P_w^0是风电预测值，ΔP_w是允许的最大偏差，三个不确定源的范围共同构成一个超立方体。

实际建模时，我不建议直接在这个连续区间上处理。原因有两个：一是连续区间与第二阶段对偶变量的乘积会构成复杂的双线性项，线性化非常麻烦；二是从工程角度看，真正让调度员头疼的是“某个源出力取到上界还是下界”，中间值的影响往往是单调的。所以实操中更常见的做法是引入0-1变量，把每个不确定源离散成两个极端状态：

u_i = u_i^0 + Δu_i · z_i, z_i ∈ {0, 1}

z_i = 0表示该不确定源取下界，z_i = 1表示取上界。再叠加预算约束：

Σ z_i ≤ Γ

Γ是预算参数，控制同时取极端的不确定源数量上限。离散化之后，不确定集合变成了一个包含有限个极端场景的组合集合，这正好为后续C&CG算法中“把最坏场景带回主问题”提供了天然的离散场景支持。

三个不确定源的具体参数化在算例里常这样设置（以下数值仅为示范，实际应基于预测数据统计得到）：

2.2 第一阶段决策模型

第一阶段决策变量主要包括：

• 机组启停状态u_g（0-1变量）
• 机组启动动作变量（用于计算启动成本）
• 预留的旋转备用容量
• 储能的充放电状态与SOC初值（如果有储能）

第一阶段模型可以写成：

min Σ (启动成本 + 空载成本) + η

s.t. 启停逻辑约束、最小启停时间约束、机组出力上下限约束、备用容量约束

这里的η是一个辅助变量，在C&CG迭代中用来逼近第二阶段的期望（最坏场景）成本，这也是两阶段鲁棒优化中“第二阶段成本被第一阶段决策耦合”的标准处理手法。

2.3 第二阶段再调度模型

给定第一阶段决策x（机组启停状态、备用等）之后，第二阶段模型针对某一具体的不确定场景u，求解最小再调度成本：

min Σ (机组出力调整成本 + 切负荷惩罚 + 弃风弃光惩罚)

s.t. 功率平衡约束 线路潮流约束（直流潮流或交流潮流线性化） 机组爬坡速率约束 切负荷量不超过该节点负荷需求 弃风弃光量不超过该时刻风光出力

写成紧凑的矩阵形式会更方便后续的对偶推导：

Q(x, u) = min_y b^T y s.t. D y ≥ h − A x − E u y ≥ 0

其中y代表第二阶段的连续决策变量（机组出力调整量、切负荷量、弃风弃光量等），h是常数向量，A、D、E是系数矩阵，u是不确定向量。

2.4 主问题与子问题的耦合结构

C&CG算法把完整的两阶段鲁棒问题拆成主问题（MP）和子问题（SP）交替求解：

主问题（在第k次迭代时）：

min c^T x + η s.t. 第一阶段约束 η ≥ b^T y_j, j = 1, 2, ..., k D y_j ≥ h − A x − E u_j^*, j = 1, 2, ..., k y_j ≥ 0

这里u_j^*是前k次迭代中，子问题找到的“最坏场景”。每迭代一次，主问题就会多一组变量y_j和对应的约束。这个机制通俗点理解就是：主问题一开始“乐观”地只考虑已知的场景，子问题每次迭代都会找出一个让主问题方案吃瘪的坏场景，把它“塞”回主问题，逼着主问题改进决策。

子问题（给定第一阶段决策x^*后）：

Q(x^) = max_{u ∈ U} min_y { b^T y | D y ≥ h − A x^− E u, y ≥ 0 }

这就是经典的max-min双层结构。外层的max寻找最坏场景，内层的min求解该场景下的最小再调度成本。这个双层结构不能直接求解，需要借助强对偶定理把它压平成单层优化问题。

3. 大M法：把双层嵌套压平的关键线性化

3.1 强对偶转换

子问题中内层的min是一个线性规划，满足强对偶条件（原问题可行且有界时，对偶问题最优值相等）。我们对内层min取对偶，得到：

Q(x^) = max_{u, λ} λ^T (h − A x^− E u) s.t. D^T λ ≤ b, λ ≥ 0 u ∈ U

这是一个纯max问题，变量是对偶变量λ和不确定变量u。到这里，max-min的双层嵌套变成了单层max，可以直接扔给求解器。

3.2 双线性项的出现与大M线性化

但注意目标函数λ^T (h − A x^* − E u) 里藏着麻烦：展开后会有λ^T E u这样的项，即对偶变量与不确定变量的乘积。λ是连续变量，u经过离散化后是0-1变量，λ_i × z_j这类项是双线性项，非线性非凸，Gurobi和CPLEX都直接处理不了。

处理思路是利用大M法把“连续变量 × 0-1变量”的乘积线性化。引入辅助变量v_ij = λ_i · z_j，用以下一组线性不等式刻画v_ij的行为：

0 ≤ v_ij ≤ M · z_j λ_i − M(1 − z_j) ≤ v_ij ≤ λ_i + M(1 − z_j)

这组约束的含义是：当z_j = 0时，第一个约束把v_ij强制压成0；当z_j = 1时，第二个约束把v_ij限制在λ_i的邻域内，当M足够大时v_ij就等于λ_i。可以把这个机制理解成“开关”——z_j是开关，打开时λ_i通过，关闭时乘积归零。

项目实践中，我曾经绕过离散化、直接用连续u变量做线性化，需要引入更复杂的多面体近似，实现难度高而且不稳定。相比之下，0-1离散化+大M法是目前最成熟、最稳的方案，特别是与预算约束天然搭配，所以我一直推荐这条路。

3.3 大M取值的实操经验

大M法看起来简单，真正翻车最多的地方就是M的取值。

M取得太小，会错误地切掉可行解，导致子问题求出来的目标值偏小，C&CG迭代出来的结果不对；M取得太大，会造成数值病态，Gurobi求解时会出现“Model may be unreliable”的警告，矩阵条件数恶化，求解速度大幅下降，甚至出现错误的整数解。

根据我的实操经验，M的取值可以参考以下逻辑：

• 先不加入大M线性化约束，求解一个松弛的子问题，观察对偶变量λ的范围
• M取λ_i可能最大值的2到5倍
• 一般情况下，电力系统调度模型里M取1e3到1e4是比较安全的区间
• 尽量不要超过1e6，超过1e6之后Gurobi的数值稳定性会明显变差

如果模型做完标幺化处理（所有功率除以基准值），λ的量级通常在0到几十之间，M取500到1000就够用了。这个细节直接决定求解稳定性，值得多花几分钟做一次敏感性测试。

4. C&CG算法：主问题和子问题怎么互相“喂料”

4.1 为什么选C&CG而不是Benders分解

两阶段鲁棒优化早期常用Benders分解，每次迭代只给主问题加一条割平面约束。C&CG的区别在于，它不仅加约束，还向主问题中直接添加一组完整的第二阶段决策变量和对应的场景约束。每轮迭代主问题的规模会增加，但换来的是割平面质量显著提升，迭代次数大幅减少。

打个比方：Benders分解像是一个学生每次考试只被告诉“总成绩不及格”，需要自己猜哪道题错了；C&CG则直接告诉他“这次考试哪道题错了、错在哪、正确答案应该是什么”。虽然每次批改的工作量更大，但学生进步快得多。对于两阶段鲁棒优化这种场景，C&CG已经是事实上的标准算法，收敛速度比Benders快一个数量级。

4.2 算法伪代码

完整流程如下：

1. 初始化：UB = +∞，LB = −∞，k = 1
2. 选取初始场景u_1^*（通常取名义值，即所有不确定源取预测值）
3. 求解主问题，得到最优解(x_k^, η_k^)，更新下界：LB = max(LB, c^T x_k^* + η_k^*)
4. 固定x_k^，求解子问题，得到最坏场景u_{k+1}^和最优值Q(x_k^)，更新上界：UB = min(UB, c^T x_k^+ Q(x_k^*))
5. 如果(UB − LB) / UB ≤ ε（ε通常取0.01或更小），则停止迭代，输出当前解
6. 否则，将新场景u_{k+1}^*对应的变量y_{k+1}和约束η ≥ b^T y_{k+1}、D y_{k+1} ≥ h − A x − E u_{k+1}^*加入主问题，k = k + 1，回到第3步

这个流程里有个微小但很重要的细节：更新UB时，每一轮都要取历史所有子问题值中的最小值，而不能只取当前值。原因是子问题在给定x_k^*下求出的最坏场景成本本来就是一个可行方案的真实成本上界，历史上任何一轮的x_k^*都是一个可行的第一阶段决策，对应一个上界。取最小值能得到更紧的上界，有助于提前收敛。

4.3 收敛判据的工程处理

理论上UB和LB会随着迭代逐渐逼近，实际代码里我还会加两个辅助收敛条件：

• 如果连续三轮迭代中新场景u_k^*与已有场景的差异小于某个阈值（比如所有场景变量的汉明距离为0），说明子问题已经找不到新信息了，提前终止，避免无谓的迭代开销
• 设置最大迭代次数（比如30次），防止数值问题导致的死循环

子问题求解时有一个常见的坑：如果对偶问题不可行或有界性异常，多半是原问题第二阶段约束A、D矩阵写错了，或者是M取值过小导致线性化约束过度收紧。排查时建议先把预算Γ设为0，退化为确定性场景，如果此时子问题正常，再逐步增加Γ，缩小问题定位范围。

5. Matlab + YALMIP代码实现与踩坑记录

5.1 环境与求解器配置

我的环境是Matlab R2022b + YALMIP + Gurobi 10.0。Gurobi求解MILP的能力业界领先，CPLEX也可以，但Gurobi在鲁棒优化这类含大量0-1变量和辅助变量的问题上表现更稳。

YALMIP的安装不用多说，重点说Gurobi参数配置。在C&CG迭代中，主问题和子问题都是MILP，有些参数值得重点关注：

ops = sdpsettings('solver', 'gurobi', ... 'verbose', 0, ... 'gurobi.mipgap', 1e-4, ... 'gurobi.numericsfocus', 1, ... 'gurobi.timelimit', 300);

numericsfocus设为1或2，可以在出现数值警告时显著提升求解稳定性，代价是求解时间略微增加，在鲁棒优化这种数值敏感性较高的问题里非常值得。timelimit一定要设，否则遇到病态模型Gurobi可能跑几小时不出结果。

5.2 主问题代码骨架

主问题用YALMIP建模的核心逻辑如下：

% 第一阶段变量 u_start = binvar(n_g, 1); % 机组启停 start_action = binvar(n_g, 1); % 启动动作 eta = sdpvar(1, 1); % 第二阶段成本逼近变量 % 第一阶段约束：启停逻辑、最小启停时间、备用容量等 Constraints = []; % 历史割集合（每轮迭代追加） % u_history 是 n_unc × k 的矩阵，存储历史最坏场景 for j = 1:k y_j = sdpvar(n_y, 1); % 第j轮添加的第二阶段变量 Constraints = [Constraints, ... eta >= b' * y_j, ... D * y_j >= h - A * u_start - E * u_history(:, j), ... y_j >= 0]; end objective = c' * u_start + eta; optimize(Constraints, objective, ops);

这段代码的关键在于：u_history(:, j)是已知常数，是子问题在历史迭代中求出来的最坏场景。主问题的决策变量是u_start、eta和所有历史添加的y_j。每轮迭代y_j会多一组，主问题规模线性增长。

5.3 子问题代码骨架

子问题经过强对偶+大M线性化之后，YALMIP代码结构如下：

% 输入：当前第一阶段解 u_start_star lambda = sdpvar(n_con, 1); % 对偶变量 z = binvar(n_unc, 1); % 不确定变量的0-1表示 % 对偶可行域 Constraints = [D' * lambda <= b, lambda >= 0]; % 不确定集合与预算约束 Constraints = [Constraints, sum(z) <= Gamma]; % 大M线性化：v = lambda * z' v = sdpvar(n_con, n_unc); M = 1000; % 根据λ的量级调整 for i = 1:n_con for j = 1:n_unc Constraints = [Constraints, ... 0 <= v(i, j) <= M * z(j), ... lambda(i) - M * (1 - z(j)) <= v(i, j) <= lambda(i) + M * (1 - z(j))]; end end % 目标函数：λ^T(h - A*x*) - sum(sum(E .* v)) % YALMIP默认求解min，因此max问题取负号 Objective = lambda' * (h - A * u_start_star) - sum(sum(E .* v)); optimize(Constraints, -Objective, ops); % 提取最坏场景 u_star = u_nominal + delta_u .* value(z); Q_star = value(Objective);

这代码里有几个关键点值得多说一嘴。大M法的核心是v(i, j)这个辅助变量，n_con代表对偶约束数量，n_unc代表不确定源数量（通常等于3乘以时段数，如果按每个时段的风、光、负荷都独立不确定来建模的话）。注意这里我没有用双循环以外的向量化写法来简化展示，实际项目里建议用repmat或kron做向量化，当n_con和n_unc变大时双循环会非常慢。

另一个容易踩的坑是目标函数里E .* v的符号。E矩阵是第二阶段约束中u的系数矩阵，它可能是负的也可能有正有负，取决于原问题的结构。目标函数整理成λ^T(h − A x^* − E u)之后，展开项−λ^T E u对应到辅助变量就是−sum(sum(E .* v))。符号反了，子问题求出来的“最坏场景”就变成了“最好场景”，C&CG会完全失效。我第一次实现时就在这里卡了一整天，最后是通过对比确定性场景（Γ=0）下的子问题和直接优化结果才发现的。

5.4 迭代主循环

主循环的写法相对直接：

Gamma = 2; % 预算参数 max_iter = 30; epsilon = 1e-3; UB = inf; LB = -inf; k = 0; % 初始场景：名义值 u_his = u_nominal; while k < max_iter k = k + 1; % 求解主问题 optimize(mp_constraints, mp_objective, ops); x_k = value(u_start); eta_k = value(eta); LB = max(LB, value(mp_objective)); % 求解子问题 [u_new, Q_k] = solve_subproblem(x_k, Gamma, ops); UB = min(UB, c' * x_k + Q_k); % 收敛判断 if (UB - LB) / abs(UB) < epsilon break; end % 向主问题添加新场景 u_his = [u_his, u_new]; add_new_cut_to_mp(u_new); % 追加 y_k 变量和相关约束 end

注意一个容易忽略的细节：LB初始化用−inf，但如果主问题在第一轮就不可行，说明第一阶段约束写得有问题，需要回头检查是不是某个约束把所有机组都禁用了或者备用容量要求过严。实际项目中我遇到过主问题第一轮就不可行的情况，排查下来是备用容量约束里留的余量太大，超出了所有机组的调节能力总和。

5.5 高频问题排查表

每次在社区里看到有人问“为什么我的C&CG不收敛”或者“为什么结果明显不对”，基本都是下面这几个原因之一：

5.6 初始化技巧

一个容易被人忽略但非常有效的技巧是用确定性场景的结果初始化C&CG。具体做法是：先令所有不确定源取名义值，求解一个确定性两阶段问题，把得到的机组启停方案作为C&CG主问题第一轮迭代的初始可行解。这样第一轮子问题就能给出一个质量相当不错的上界，能明显减少迭代次数。

如果模型中包含储能系统，建议把储能的初始SOC也纳入第一阶段决策变量处理。因为SOC的初值直接影响后续所有时段的充放电可行域，如果放到第二阶段再决定，会破坏“日前决策+实时调整”的时序逻辑。

6. 从能跑到能用的进阶建议

6.1 预算Γ的调参逻辑

预算参数Γ是鲁棒优化模型提供的一个“保守度旋钮”。Γ=0时模型退化为确定性优化，成本最低但对不确定性没有任何抵御能力；Γ等于不确定源总数时模型最保守，要求所有源同时取极端值时方案依然可行，成本最高但安全性最强。

实际工程中，我会这样给决策者分析：跑一组Γ从0到最大值的结果，画出成本-保守度曲线。曲线上每个点的斜率可以理解为“再买一单位安全性的边际成本”。大多数情况下，Γ从小往大增加时，成本会先陡增后趋缓。拐点位置的Γ值往往就是推荐的折中选择——既能覆盖绝大部分风险，又不会让成本高到离谱。

6.2 与大M法相关的几个数值稳定性经验

大M法是最容易引起数值隐患的环节。除了前面说的M取值，还有几个经验值得分享：

• 所有变量尽量统一量纲，功率用MW，成本用元，不要混着用kW和万元
• 如果模型中同时存在数量级相差1e6以上的系数（比如M取1e5，而其他约束系数在个位数），Gurobi的内置预求解器有时会自行简化模型，反而导致线性化约束被错误消去。这种情况可以在Gurobi参数里关掉Presolve的部分功能，或者把M缩小
• 大M法的M值即使是同一个模型，不同约束上也可以不一样。比如某个对偶变量λ_i理论上限是100，那么对应的M取300就足够；另一个λ_j上限是1000，M可以取3000。分别给值比统一一个大M要稳得多

6.3 算例验证的推荐做法

验证代码时，我建议按以下顺序推进：

先在修改的6节点系统上测试，三个不确定源分别对应风机节点、光伏节点和负荷节点。预测值随便设一组合理的数，偏差范围按10%~20%取。先跑Γ=0的确定性情形，确认主问题和子问题在名义场景下求解正常，再逐步增加Γ。

然后对比确定性方案和鲁棒方案在极端场景下的表现：把确定性方案代入子问题，得到它在最坏场景下需要付出的高额再调度成本（可能伴随切负荷）；再把鲁棒方案代入，成本会低得多。这个对比是论文里最常用的结果展示，也是验证鲁棒优化价值最直观的方式。

等6节点算例跑通后，再移植到IEEE 14节点或118节点系统。移植时唯一要注意的是网络约束矩阵的维度和稀疏性，YALMIP里用稀疏矩阵构造约束比用全矩阵快很多，内存占用也小很多。

套用这套框架，实际做过的项目里，6节点系统一般5轮以内收敛，118节点系统大概需要10到20轮，单轮求解时间视Gurobi性能从几秒到几分钟不等。如果迭代超过30轮还没收敛，几乎可以肯定是代码里有bug，别怀疑是算法本身的问题。

6.4 扩展方向

这套两阶段鲁棒优化的框架本身是非常通用的。换一个不确定集合的定义方式（比如用数据驱动的多面体集合代替盒式集合），换一个目标函数（比如加入碳排放约束），或者把第二阶段的线性规划换成混合整数规划（比如考虑储能充放电状态的实时调整），大框架都不需要动，只需要改对应的约束和变量定义。

更进一步，如果想考虑不确定性之间的时空相关性，可以考虑用椭球不确定集合或者数据驱动的不确定集合，这时大M法的角色会弱化，但C&CG主-子问题迭代的框架依然成立。

最后说一点个人体会。这套方法表面上看是数学建模和代码实现的功夫，真正吃透之后，最有价值的其实是它对“不确定性成本”的量化能力。当决策者问“鲁棒方案为什么比确定性方案贵这么多”，用一条Γ-成本曲线就能讲清楚：你愿意花多少钱，换来对多少不确定性的保障。这种把工程决策翻译成优化语言的能力，是做电力系统优化工作的人最值钱的本事。