企业官网建设流程全解析

搞电力系统调度方向的，应该都懂这种体验：花了好几周把确定性机组组合模型跑通，工况测试也很完美，结果把真实的风电预测数据接进去，第一轮仿真直接翻车——不是功率平衡被打破，就是备用容量见底。原因就一个：你把预测曲线当成真实出力来用了。风电实际出力和预测值差个20%是家常便饭，光伏遇到云层抖动更是断崖式变化，负荷侧也没好到哪里去。所以最近几年，把风、光、负荷不确定性直接建模进优化问题里的“两阶段鲁棒优化”越来越多，配合大M法做线性化、C&CG算法做迭代求解，几乎成了这个方向的标配组合。

这篇文章我从实际代码落地的角度，把两阶段鲁棒优化怎么建模、大M法到底用在哪里、C&CG迭代是个什么逻辑，以及Matlab里怎么一步步实现，完整拆一遍。最后会附上我在调参和排错过程中总结的几个坑。适合正在做电力系统优化方向研究生课题、或者刚接触鲁棒优化想快速上手的同学参考。

1. 风、光、负荷同时波动时，确定性优化为什么撑不住

1.1 预测曲线为什么不靠谱

电网调度的核心，是提前一天甚至提前数小时做出机组启停和出力计划。但这个计划的输入——风电出力、光伏出力、负荷预测——本质上都是随机量。

风电出力受风速影响，风速预测误差通常用正态分布近似，但实际尾部分布特别厚，极端场景下偏离预测值30%以上毫不稀奇。光伏出力和云层运动强相关，分钟级变化可能达到装机的50%。负荷预测虽然相对稳定，但遇到极端天气、节假日、突发工业负荷，系统性偏差也常有。

确定性优化的做法是把这些预测值当成已知参数代入模型，求出唯一最优解。看似合理，但问题是：预测误差根本没有被纳入约束校验。你的最优解可能刚好卡在一个窄通道里，真实场景稍微一偏，约束直接违反。

1.2 确定性模型在实际运行中的两种典型失效

第一种是备用不足。比如某一时段火电机组已经压到最低技术出力，光伏预测值是满发，但实际阴天，出力腰斩，系统频率就会往下掉。模型里没有场景校验，自然不会预留足够的上调空间。

第二种是经济性失真。预测误差较大时，确定性优化给出的机组组合在真实场景下根本不可行，需要实时市场高价购买调节能力，账面上算出来的“最优成本”和实际运行成本差一大截。

这本质上不是一个“数据质量”问题，而是“建模范式”问题。参数不确定的情况下，单点预测的最优解在完整场景空间里往往不是好解。

1.3 鲁棒优化的思路：我不要最优，我要“最坏情况下也还行”

鲁棒优化换了种思路：不追求在某个随机分布下期望成本最小，而是保证在最坏输入场景下方案依然可行，并且成本可控。

这轮到我当时的理解：它本质上是“先做事后复盘”的两阶段决策结构。第一阶段先把能定的决策定了，比如机组开不开机、启停成本花了多少；第二阶段等不确定参数揭晓后，再做经济调整，比如机组出力上调下调、切负荷量是多少、弃风弃光量是多少。

把这个两级结构写进优化目标，就成了一个套路非常固定的数学模型，也就是下一节要拆解的框架。

2. 两阶段鲁棒模型建模：把“先决策后调整”的调度逻辑落成数学表达

2.1 变量分类：哪些是第一阶段，哪些是第二阶段

第一阶段变量，实际中对应的是“必须在不确定性揭晓前作出的决策”。最常见的是机组组合中的0/1启停状态变量，也包括需要提前敲定的和上级电网的购电计划、储能日前充放电计划。这些变量一旦定了，当天没法改。

第二阶段变量，对应的是“不确定性揭晓后可以调整的变量”。比如火电机组的实时出力增量、切负荷量、弃风弃光量、储能实时出力。这些是在看到实际风光出力之后才做的决策。

这里有一个初学者容易搞混的地方：第一阶段变量里也有连续变量，比如机组的基础出力和计划出力，不是只有0/1变量。区别不在于变量类型，而在于决策时刻不同。

2.2 目标函数：最坏场景下的调整成本也要算进来

一个典型的两阶段鲁棒机组组合目标函数长这样：

$$ \min_{x, p} \left[ \sum_{i} C_i^{start} x_i + F_i(p_i) \right] + \max_{u \in \mathcal{U}} \min_{y \in \mathcal{F}(x,p,u)} \left[ \sum_j c_j y_j \right] $$

拆开看就是三层：

• 第一层最小化：决策第一阶段变量x和p，包括启停成本加上运行成本。
• 第二层最大化：在不确定性集合中找一个“最坏场景”u。
• 第三层最小化：在给定第一阶段决策和最坏场景下，以最小代价做出实时调整y。

计算流程上，可能出现的成本形态是多样的。如果切负荷有惩罚，那就把切负荷量乘一个大惩罚系数。如果允许弃风弃光，也可以加对应的违约金。毕竟鲁棒优化不排斥第二阶段的灵活性，它追求的是“无论场景多坏，总有办法用最小的调整代价扛住”，扛不住才给高惩罚。

2.3 不确定性集合：盒式加预算Γ

不确定性参数u（风、光、负荷的一组向量）不能随便在±∞里跑，否则模型要么无解要么保守到没意义。工程上最常用的是带预算的盒式集合：

$$ \mathcal{U} = \left{ u ; \middle| ; u_i = u_i^f + \zeta_i \Delta u_i, ; |\zeta_i| \le 1, ; \sum_i |\zeta_i| \le \Gamma \right} $$

含义是：每个不确定参数在预测值附近波动，波动幅度上限是预测误差界的delta值；但同时限制所有参数同时偏离预测值的个数（或总偏离程度）不超过预算Γ。

Γ是个很有物理意义的参数：Γ=0就是把不确定集合压缩成一个点，退化成确定性优化；Γ取到全部不确定量的个数，则是完全保守，所有参数一起达到最坏情况；实际工程中Γ取3到8个通常就够覆盖风险了，这相当于承认“不可能所有风电场和负荷节点同时冲到极端”。

这里我补充一个建模细节：$|\zeta_i| \le 1$ 是绝对值形式。处理时通常引入辅助变量把绝对值展开成线性约束，或者直接用二进制变量表示ζ的正负偏移，这也是下一节大M法操作的对象之一。

3. 大M法真正的用武之地：对偶之后的双线性项与线性化实现

3.1 子问题为什么要做对偶变换

C&CG算法的子问题长这样：

$$ SP(x^) = \max_{u \in \mathcal{U}} \min_{y \in \mathcal{F}(x^, u)} b^T y $$

外层max和内层min套在一起，没法直接扔给求解器。标准的数学处理是把内层min问题写成线性规划，然后做对偶变换。

对偶是个很不直观的东西，但道理不复杂：对偶变换可以把“最小化”问题变成“最大化”问题，于是max-min就变成了max-max，直接合并成一个max问题。约束也同时从关于y的约束变换成关于对偶变量的约束。

具体操作上，对偶变换需要把原问题写成标准LP形式，逐条约束到对偶变量一一对应。这一段是整个建模最容易出错的地方，因为任何一条约束漏掉，后面的对偶约束就不对。

3.2 双线性项出现在哪里

对偶变换完成后，目标函数里多了一项对偶变量λ和不确定性变量u的乘积：

$$ \lambda^T (D u) $$

问题在于：λ是第二阶段问题的对偶变量，你需要求解才能知道；u是外层max的决策变量，也在优化里。两个都是变量，乘在一起就是双线性项，也就是非凸项。非凸意味着不能直接扔给MILP求解器。

这就是大M法要出场的地方。

3.3 大M线性化的具体写法

核心思路：把双线性项里的其中一个变量离散成有限个取值，然后用二进制变量去选择，再引入辅助变量，用一组带大M的约束把乘积项线性化。

举个具体例子。如果不确定参数u_i = u_i^f + ζ_i Δu_i，ζ_i ∈ [-1,1]，我们可以把ζ_i离散成有限取值，比如区间上下限和预测值对应的几个关键层级。引入二进制变量z_{i,k}，令：

$$ \zeta_i = \sum_k z_{i,k} \cdot \zeta_{i,k}, \quad \sum_k z_{i,k} = 1, \quad z_{i,k} \in {0,1} $$

然后产品λ_i * ζ_i就变成λ_i和z_{i,k}的乘积。再引入辅助变量$w_{i,k} = \lambda_i \cdot z_{i,k}$，用如下大M约束来约束w：

$$ w_{i,k} \ge \lambda_i - M(1 - z_{i,k}) $$

$$ w_{i,k} \le \lambda_i + M(1 - z_{i,k}) $$

$$ w_{i,k} \ge -M z_{i,k} $$

$$ w_{i,k} \le M z_{i,k} $$

这组约束的逻辑是：当$z_{i,k}=1$时，w被强制等于λ_i；当$z_{i,k}=0$时，w被压到0。这样乘积项就被精确表示出来了，而整个问题变成了一个混合整数线性规划(MILP)。

3.4 大M取值的工程经验

大M的取值是个非常经典的坑。

M太小，比如小于λ_i可能达到的最大绝对值，约束就会过度束缚，把可行域错误地掐掉，得到的“最优解”其实是假解。M太大，比如取1e10，求解器在数值上会出现病态，矩阵条件数恶化，求解时间飙涨，甚至收敛到错误的整数解。

我的经验做法分两步：

1. 先把子问题里的大M项全部去掉（相当于把双线性项忽略），做一个LP松弛求解，记录λ_i在这个松弛过程中的最大绝对值量级。
2. 取这个量级的10到100倍作为M，再正式求解。

这个做法在大多数场景下能避开数值问题，又不至于过紧。

4. C&CG迭代框架拆解：主问题给乐观解，子问题找最坏场景

4.1 C&CG和Benders分解的本质区别

很多教程会把C&CG和Benders放在一起讲，因为两者的迭代框架很像，都是主问题-子问题交替求解。

但核心区别要搞清楚：Benders分解在每次迭代后给主问题添加的是对偶信息构成的割平面约束，也就是“你的决策在这个方向上不够好”这种信息；C&CG则更“暴力”，它是把第二阶段的决策变量和对应的约束原封不动地复制一份，加到主问题里。

C&CG的优势在于收敛速度，因为它加入了具体场景下的完整可行域信息，而不是只给一个割平面。劣势则是主问题规模会随着迭代次数增大——每轮都加一组变量和约束。所以实际使用中C&CG的迭代轮数通常不多，两阶段问题往往几十轮内就收敛到1e-3，主问题规模还是可控的。

4.2 主问题的结构：每轮迭代加一组“列+约束”

记第k轮已经识别出的最坏场景为u^(1), u^(2), ..., u^(k)。主问题MP如下：

$$ \min_{x, \theta, y^{(l)}} c^T x + F(x) + \theta $$

$$ \text{s.t.} \quad \theta \ge b^T y^{(l)}, \quad l = 1,\dots,k $$

$$ B x + C y^{(l)} \le g - D u^{(l)}, \quad l = 1,\dots,k $$

注意，θ是一个标量辅助变量，用来逼近“最坏场景下的第二层成本”。每轮迭代新发现的最坏场景u^(l)会带入一组新的第二阶段变量y^(l)和相应约束。这组新约束就是“列和约束”这个名称的来源。

因为主问题只考虑了截至目前发现的最坏场景，那些还没被发现的场景约束还没进来，所以主问题解出来的目标值是真实最优解的一个下界（LB）。

4.3 子问题的“压力测试”角色与上下界更新逻辑

子问题输入一个固定的第一阶段解x^(*)，然后在整个不确定集合里搜索，看能不能找到让第二阶段调整成本更高的场景。

$$ SP(x^) = \max_{u \in \mathcal{U}} \min_{y \in \mathcal{F}(x^,u)} b^T y $$

如果子问题返回的成本很高，说明主问题低估了风险，要拿这个场景回去“惩罚”主问题。如果子问题返回的成本和主问题的θ差不多，说明主问题考虑的场景已经足够全面，算法收敛。

上下界的更新逻辑也很直觉：

• LB是主问题的最优目标值，因为主问题只考虑了已发现场景，必然 ≤ 真实鲁棒最优值。
• UB是x^()的真实成本c^T x^() + SP(x^())，因为SP是在整个不确定集合中的最坏成本，所以c^T x^() + SP(x^(*))是某个可行方案的真实成本，取所有迭代里的最小值就是真实最优值的上界。

C&CG的完整迭代流程：

1. 初始化：设LB=-∞，UB=+∞，迭代轮数k=1。给出初始最坏场景u^(1)，通常取预测值中点或区间中点。
2. 求解主问题MP，得到x^()、θ^()和LB。
3. 固定x^()，求解子问题SP，得到最坏场景u^(k+1)和子问题目标值，更新UB = min(UB, c^T x^() + SP(x^(*))）。
4. 若(UB-LB)/|UB| < ε，停止；否则k=k+1，把u^(k+1)加入主问题，回到步骤2。

5. Matlab代码骨架：从YALMIP变量定义到主-子循环完整实现

5.1 环境准备与工具箱配置

Matlab代码实现的方案里，YALMIP是我推荐的最省事工具箱，它把变量定义、约束拼接、求解器调用统一抽象，写MILP比手写矩阵方便太多。子问题对偶后的双线性项也能用YALMIP的binvar和约束定义直接拼出来。

求解器方面，MILP建议用Gurobi或者CPLEX，免费替代方案可以考虑SCIP。纯lp问题用Matlab自带linprog够用，但MILP性能会差不少。

% 检查并初始化YALMIP yalmip('clear'); % 如果没安装，可以用 addpath(genpath('D:\yalmip')) 类似方式加载 ops = sdpsettings('solver', 'gurobi', 'verbose', 0); ops.mip.tolerancemip = 1e-6; ops.mip.tolerancegap = 1e-4;

5.2 主问题MP代码骨架

主问题的构建需要每轮迭代往里加场景约束。核心数据结构是：所有阶段的变量重新定义在同一组变量上，通过场景索引区分。

% 基础参数 nG = 5; % 火电机组数 nW = 2; % 风电场数 nPV = 1; % 光伏电站数 nL = 3; % 负荷节点数 % 第一阶段变量 x_uc = binvar(nG, 1); % 机组启停0/1变量 p_g = sdpvar(nG, 1); % 机组计划出力 theta = sdpvar(1, 1); % 鲁棒辅助变量，表示最坏场景第二阶段成本 % 基础约束 Constraints = []; Constraints = [Constraints, p_g >= p_g_min .* x_uc]; Constraints = [Constraints, p_g <= p_g_max .* x_uc]; Constraints = [Constraints, sum(p_g) == D_forecast - sum(Pw_forecast) - sum(Ppv_forecast)]; Objective = sum(C_start .* x_uc) + sum(a_g .* p_g) + theta;

每轮迭代发现新的最坏场景u_k后，往主问题里追加一组第二阶段变量和约束：

% 第k次迭代：新增场景对应的一组决策变量 delta_p = sdpvar(nG, 1); % 机组出力修正量 lshed = sdpvar(nL, 1); % 切负荷量 wcurt = sdpvar(nW + nPV, 1); % 弃风弃光量 % 场景约束：实时功率平衡 Constraints = [Constraints, ... sum(delta_p) + sum(Pw_scene_k + wcurt(1:nW)) + sum(Ppv_scene_k + wcurt(nW+1:end)) ... == D_scene_k - lshed]; % 机组调节范围 Constraints = [Constraints, delta_p >= -ramp_down .* x_uc]; Constraints = [Constraints, delta_p <= ramp_up .* x_uc]; % 切负荷不能超过负荷值 Constraints = [Constraints, 0 <= lshed <= D_scene_k]; % 弃风弃光不能超过实际出力 Constraints = [Constraints, 0 <= wcurt <= [Pw_scene_k; Ppv_scene_k]];

注意这里的Pw_scene_k、D_scene_k是第k轮子问题返回的最坏场景值，不是预测值。

主问题的目标函数和约束拼接好以后直接optimize：

diagnose_MP = optimize(Constraints, Objective, ops); x_uc_val = value(x_uc); p_g_val = value(p_g); LB = value(Objective);

5.3 子问题SP代码骨架

子问题是整个实现的难点，因为要对偶、要引入大M线性化双线性项。这里给出一个更工程化的处理方式：将不确定量离散成若干个关键层级，用二进制变量选择层级，然后大M线性化。

% 不确定性离散点：比如每个风电场取{下界, 预测值, 上界}三档 % 这里假设每个不确定参数离散成3档 z_w = binvar(nW, 3, 'full'); % 风电不确定性档位选择 z_pv = binvar(nPV, 3, 'full'); % 光伏不确定性档位选择 z_l = binvar(nL, 3, 'full'); % 负荷不确定性档位选择 % 每个不确定量只能选择一档 Constraints_SP = []; for i = 1:nW Constraints_SP = [Constraints_SP, sum(z_w(i,:)) == 1]; end % 类似处理z_pv和z_l % 不确定量取值为档位值与二进制变量的线性组合 Pw_actual = sdpvar(nW, 1); for i = 1:nW Pw_actual(i) = z_w(i,:) * [Pw_min(i); Pw_forecast(i); Pw_max(i)]; end

然后是预算约束Γ。如果每个不确定量的极端偏离都对应|ζ|=1，那么预算约束可以线性表示为所有选中“上界或下界档位”的个数不超过Γ：

% 预算约束：所有偏离预测值档的个数不超过Gamma zeta_indicator = [z_w(:,1) + z_w(:,3); z_pv(:,1) + z_pv(:,3); z_l(:,1) + z_l(:,3)]; Constraints_SP = [Constraints_SP, sum(zeta_indicator) <= Gamma];

双线性项的处理是子问题里最核心的代码。子问题的目标函数对偶化后，通常会出现λ * Pw_actual这样的项。我们把Pw_actual换成档位的线性组合，乘积项变成λ * z_w(i,k)。引入辅助变量w_w(i,k)，用大M约束线性化：

M_big = 1000; % 根据对偶变量λ的量级调整，经验值见下节 lambda = sdpvar(nBus, 1); % 对偶变量，维度与功率平衡约束条数一致 w_w = sdpvar(nW, 3, 'full'); for i = 1:nW for k = 1:3 Constraints_SP = [Constraints_SP, w_w(i,k) <= M_big * z_w(i,k)]; Constraints_SP = [Constraints_SP, w_w(i,k) >= -M_big * z_w(i,k)]; Constraints_SP = [Constraints_SP, w_w(i,k) <= lambda(i) + M_big * (1 - z_w(i,k))]; Constraints_SP = [Constraints_SP, w_w(i,k) >= lambda(i) - M_big * (1 - z_w(i,k))]; end end

子问题目标函数写成：

% 对偶变换后的目标函数，双线性项被w_w代替 Objective_SP = sum(dual_coeff .* lambda) + sum(w_w .* Pw_levels); Optimize_SP = optimize(Constraints_SP, -Objective_SP, ops); % 求max，YALMIP里写成负号 u_next = value([Pw_actual; Ppv_actual; D_actual]); % 本轮识别的最坏场景 UB_k = value(c' * x_uc_val + Objective_SP);

5.4 主循环与收敛判断

LB = -inf; UB = inf; k = 0; u_hist = {}; Gap = inf; % 初始场景：取预测值中点 u_k = [Pw_forecast, Ppv_forecast, D_forecast]; while Gap > 1e-3 k = k + 1; u_hist{k} = u_k; % 1. 构建并求解主问题MP（把u_hist中所有场景的约束都加进去） [Constraints_MP, Objective_MP] = build_MP(u_hist); optimize(Constraints_MP, Objective_MP, ops); LB = value(Objective_MP); x_k = value(x_uc); p_k = value(p_g); % 2. 固定x_k，构建并求解子问题SP [Constraints_SP, Objective_SP] = build_SP(x_k, p_k); optimize(Constraints_SP, -Objective_SP, ops); SP_val = value(Objective_SP); UB = min(UB, value(c' * x_k + p_k' * a_g) + SP_val); % 3. 提取新的最坏场景 u_k = get_worst_scene(); % 从子问题变量中提取Pw/Ppv/D % 4. 更新gap Gap = (UB - LB) / abs(UB); fprintf('Iter %d: LB=%.4f, UB=%.4f, Gap=%.6f\n', k, LB, UB, Gap); end

这个骨架已经是能跑起来的主逻辑了，剩下就是根据你自己的具体系统替换约束和数据。

6. 调参三个关键坑：大M取值、预算Γ扫描、收敛判定

6.1 大M取值：调小了假可行，调大了数值病态

我在第3、4节都提过大M，这里集中说一次真正的操作经验。

一个很典型的翻车现象是：主问题LB明明还在上升，但子问题SP突然返回一个异常大的值，然后UB跳到一个天文数字，gap飙到百分之几千。这种时候先怀疑的不是模型约束，而是大M的取值。

大M过小时，线性化约束会把本应可行的对偶变量组合限制住，子问题找不到真正的最坏场景，SP值虚低或虚高都不一定，但LB和UB的gap就会不稳定。

处理建议是这样的：

1. 先不要一次性把整个鲁棒模型调通，而是把子问题单独拿出来，在某个固定x下，去掉大M约束，用非线性求解器或枚举法求一个参考解，确定λ的实际范围。
2. 再用这个范围定M，比如M = 100 * max|λ|。
3. 如果系统规模特别大，不确定量十几甚至几十个，建议考虑把所有不确定量的档位选择集中成一个二进制变量组，避免对偶变量的维度扩张导致M估计不准。

还有些求解器对大M很敏感，Gurobi算MILP时数值鲁棒性相对好一些，但M超过1e7之后依然会出现数值警告。这时候可以考虑把目标函数和约束里的量纲统一缩放，比如功率用p.u.而不是MW。

6.2 不确定预算Γ怎么选：画成本-保守度曲线

Γ的取值没有理论最优，除非你有明确的可靠性指标。

我常用方法：固定其他参数，把Γ从0扫到最大可能值（即所有不确定参数个数之和），每取一个Γ跑一遍完整的C&CG，记录目标函数值。画出来的曲线通常是一个先快速增长、然后趋于平缓的S形。

拐点处的Γ就有参考价值。比如风电、光伏、负荷不确定参数加起来有8个，曲线拐点在Γ=3附近，那就说明系统主要风险来自3个左右同时极端偏离的参数组合，再增加Γ只会徒增成本，对可靠性几乎没有额外提升。

这比拍脑袋取一个Γ要靠谱得多，审稿人问起来也更有理有据。

6.3 判断收敛别只看绝对间隙

很多初版代码用UB - LB < epsilon判断收敛。如果成本量级是1e5，epsilon取1e-4，循环可能永远跑不完；如果epsilon取太大（比如1），精度又难保证。

推荐写法是相对间隙：

$$ \text{Gap} = \frac{|UB - LB|}{\max(|UB|, |LB|, 1)} $$

一般取1e-3到1e-4。另外迭代后期gap下降会非常慢，如果你只需要方案可行性而不是严格最优，gap到5e-3也可以接受，实际工程中这个误差远小于预测误差本身。

6.4 迭代初始化：初始场景选得好，迭代次数减一半

主问题第一次迭代需要一个初始最坏场景。很多人直接用区间中点，也就是预测值本身。这确实是最简单的方式，但不一定最优。

如果C&CG前几次迭代发现每次子问题返回的最坏场景变化很大，gap下降慢，可以考虑的初始化方式有两种：

1. 把预测值场景先跑一遍确定性优化，用得到的x作为初始解传给C&CG。
2. 分别对每个不确定参数单独取自己的边界（其他参数取预测值），跑一遍子问题，选其中成本最高的场景作为初始场景。

第二种方法等于告诉算法“系统最怕哪个参数单独出问题”，虽然多花一次子问题求解的时间，但往往能省下好几轮外循环。

我在实际项目里的体会是，两阶段鲁棒优化这套东西，数学上并不复杂，真正的门槛在建模转换和工程细节。对偶不能漏约束、大M不能乱取、主问题场景约束不能加错，任何一个点出错，结果都会偏离真实解。前期耐心点，先把小规模案例跑透，再往上加复杂度和不确定性参数数量，会比一上来就套大模型顺畅得多。