企业官网建设流程全解析

接手微网容量配置项目时，我一直被同一个问题反复拷问：业主问"你这套优化方案的底气从哪来"。我用两阶段鲁棒优化算法做微网多电源容量配置，在Matlab里搭了完整的代码框架，前后折腾了三周才把模型、主子问题迭代和数值稳定性全部跑顺。这篇文章把当时的完整思路和代码骨架整理出来——不只是贴代码，更重要的是讲清楚每一步为什么这样做。适合正在做微电网规划、储能配置，或者对鲁棒优化算法感兴趣的工程师和研究生参考。

1. 为什么微网容量配置必须拥抱两阶段鲁棒

1.1 确定性优化的阿喀琉斯之踵

先还原一下场景。常规的容量配置思路，是拿历史负荷和风光出力曲线，塞进一个确定性优化模型里，一次求解出"光伏装多少、风电装多少、柴油机和储能配多大"。这类模型本身没有错，问题出在输入数据的单点假设上——负荷取预测均值，光伏出力取典型曲线，风电出力也取期望值，然后求解器给出一个在"平均天气剧本"下成本最低的方案。

这个方案看起来完美，但只要你把时间尺度拉长到项目全生命周期，就会发现问题：微网容量一旦落地，基本20年不变，而期间的天气模式、负荷增长、电价波动都是不确定的。用一个单点预测值去对抗长期不确定性，本质上是在赌未来二十年的剧本不会偏离预测。我在实际项目中见过太多"确定性最优解"在连续阴雨天里被购电成本压垮的案例——那些看似省了百分之十几投资额的方案，在最不利场景下运行一年就能把节省的钱全部亏回去。

1.2 两阶段鲁棒的直观思想

两阶段鲁棒优化解决的就是这个"决策不可逆 + 环境不可知"的矛盾。它的核心结构是min-max-min：外层min是投资决策（第一阶段），内层是两个min和max的组合——在给定投资方案下，考虑所有可能的不确定场景（max），对每个场景都做出最优运行调度（min），最后找到"最坏场景下运行成本最低"的投资方案。

用开餐厅来类比最直观：你决定后厨设备买多少（第一阶段，花一大笔钱），之后每一天的备菜量、进货量则根据当天客流和菜价调整（第二阶段，属于日常运行）。你没法控制客流和菜价，但希望在最差的一天——客流爆满且菜价飙涨——餐厅依然不会亏损。鲁棒优化帮你把设备容量定在"最坏一天也扛得住"的水平，而不是"平均一天刚好够用"的水平。对应到微网，第一阶段变量是各电源的安装容量；第二阶段变量是每个时段里柴油机出力、储能充放电功率、向主网购电功率；不确定性则来自光伏出力、风电出力和负荷预测的偏差。

1.3 适用场景与前置条件

当然，不是所有项目都得用两阶段鲁棒。如果项目负荷曲线平滑、可再生能源占比低、弃光弃风风险小，那么确定性优化加灵敏度分析就够了。但以下情况我会强烈建议上两阶段鲁棒：可再生能源渗透率超过40%、电网购电容量受限或电价峰谷差极大、业主对停电风险容忍度很低、或者你手里有多套历史气象数据却无法确定未来真实分布。

前置条件有两个必须提前评估。第一，求解器支撑。虽然模型本身规模不大，但加上时间耦合（储能SOC）和迭代生成的场景后，需要Gurobi或CPLEX这种级别的商业求解器，开源的CBC和SCS在大型算例下会很吃力。第二，建模功底。两阶段鲁棒的难点从来不在算法本身，而在把物理系统准确翻译成数学约束，尤其是储能充放电的二进制变量、不确定性集合的构造、以及子问题对偶时双线性项的处理。这三块做扎实了，代码实现反而是水到渠成的事。

2. 数学模型：投资变量、运行变量与不确定性集合的三层设计

2.1 决策变量的归属

数学建模第一步，是明确哪些变量属于第一阶段、哪些属于第二阶段，这个归属不清，后面所有代码都会乱。

第一阶段变量（投资层）是容量决策，整个项目周期内不发生改变：

• x_pv：光伏安装容量（kW）
• x_wt：风电安装容量（kW）
• x_dg：柴油发电机安装容量（kW）
• x_bat：储能安装容量（kWh）

第二阶段变量（运行层）随调度时段变化，在本文采用24小时典型日的滚动调度框架：

• P_dg(t)：柴油机在时段t的出力
• P_pv(t)、P_wt(t)：光伏、风电在时段t的实际出力（受不确定性和安装容量共同约束）
• P_bat_ch(t)、P_bat_dch(t)：储能充放电功率
• SOC(t)：储能荷电状态，连接相邻时段
• P_grid(t)：时段t向主网的购电功率

这里有一个容易踩的坑：储能容量和功率是两个维度，一个决定能量上限，一个决定功率上限。简化处理时我通常把储能容量x_bat当作能量上限，而充放电功率上限按容量乘以一个倍率系数来换算，例如0.5C，也就是功率上限 = 0.5 × x_bat。这个系数要让业主确认，因为不同储能厂商的倍率特性差别很大。

2.2 不确定性集合到底怎么选

两阶段鲁棒的性能高度依赖不确定性集合的形状，集合太小则鲁棒效果不足，集合太大会让解过度保守。

实际工程中我推荐盒式集合加预算约束的组合： U = { u | u_nom - u_dev ≤ u ≤ u_nom + u_dev, Σ |u_i - u_nom,i| / u_dev,i ≤ Γ }

这个设计的妙处在于预算参数Γ。它限制的是"所有时段同时偏离预测值的总量"，而不是逐点都取极端值。否则如果允许24个时段的光伏同时比预测低20%，模型求出来的解只会是"光伏全装最大"，没有任何参考价值。Γ取值的经验规律是：Γ越小，决策越接近确定性优化；Γ越大，决策越保守。具体怎么定量选，我在第6章专门讲。

除了盒式集合，还有椭球集合（适合不确定量服从正态分布的场合）和场景集合（适合有典型历史场景的场合）。椭球集合在求解时需要二阶锥约束，对求解器要求更高，我在微网规划里用得少，主要因为业主更喜欢用"偏差百分比"这种直白的语言沟通，而盒式集合的参数含义一目了然。

2.3 完整模型

下面给出一个能够直接移植到代码的紧凑模型。目标函数写作：

min_x c_inv^T x + max_{u∈U} min_{y∈Ω(x,u)} c_run^T y

其中c_inv是单位容量投资成本向量，c_run是运行成本向量（包括燃料成本、购电成本、运维成本）。第一阶段的约束包括：

• 容量上下限：x_min ≤ x ≤ x_max
• 投资预算约束：c_inv^T x ≤ B_inv

第二阶段的运行约束Ω(x,u)包括：

• 功率平衡：P_pv(t) + P_wt(t) + P_dg(t) + P_bat_dch(t) + P_grid(t) = P_load(t) + P_bat_ch(t)
• 柴油机出力范围与爬坡约束：P_dg_min ≤ P_dg(t) ≤ P_dg_max，以及相邻时段爬坡限制
• 光伏/风电出力受容量和不确定性共同约束：P_pv(t) ≤ u_pv(t) × x_pv
• 储能SOC递推：SOC(t+1) = SOC(t) + η_ch × P_bat_ch(t) - P_bat_dch(t) / η_dch
• 储能容量约束：SOC_min × x_bat ≤ SOC(t) ≤ SOC_max × x_bat
• 联络线功率上限：0 ≤ P_grid(t) ≤ P_grid_max

这个模型的阶段耦合点在于：第一阶段的容量x通过约束"P_dg(t) ≤ x_dg"、"P_pv(t) ≤ u_pv(t)×x_pv"等进入第二阶段，从而把投资和运行绑定在一起。鲁棒性就体现在：不确定参数u不是固定的，而是在集合U内由max层主动寻找最坏取值。

3. C&CG迭代：主问题、子问题与对偶处理的完整逻辑

3.1 从min-max-min到主子问题迭代

两阶段鲁棒问题不能直接用常规优化器求解，因为max-min嵌套结构让模型非凸。工程中最常用的求解框架是列与约束生成算法（C&CG），也叫Zeng-Zhao算法。它的核心思路是：不一次性枚举所有不确定性场景，而是让子问题主动找出"最坏场景"，把这个场景作为新的约束补充到主问题中，反复迭代直到收敛。

主问题是最小化投资成本加一个辅助变量η，η代表"当前已发现的所有场景中，最大的运行成本"。主问题形式如下： MP: min c_inv^T x + η s.t. (第一阶段约束) 对每个已发现的场景u_i： (第二阶段约束，针对场景i)，且 η ≥ c_run^T y_i

子问题则是在给定主问题解x后，寻找使运行成本最大的不确定性场景： SP: Q(x) = max_{u∈U} min_{y∈Ω(x*,u)} c_run^T y

C&CG的收敛性已经在理论上被证明：当不确定性集合U是多面体时，算法有限步内收敛到全局最优。这对我们来说是个好消息，意味着不需要担心迭代过程陷入局部最优。

3.2 子问题的对偶变换与双线性项处理

实现子问题的关键难点在于嵌套结构。对于固定的x*，内层是一个以y为决策变量的线性规划（LP）。LP满足强对偶条件，因此可以把内层min问题替换为其对偶问题，从而把max-min结构变成一个等价的max问题：

内层min的对偶形式为： max_{λ,μ} (b + C×u + D×x*)^T λ + e^T μ s.t. A^T λ + E^T μ = d λ ≥ 0

这里λ和μ是对偶乘子，分别对应第二阶段的不等式约束和等式约束。合并到外层max后，子问题变成： max_{u∈U, λ, μ} (b + C×u + D×x*)^T λ + e^T μ s.t. A^T λ + E^T μ = d，λ ≥ 0，u ∈ U

问题来了：目标函数中出现了λ^T C×u这种双线性项，两个变量相乘，子问题不再是一个线性规划。处理双线性项的常规手段是McCormick包络或大M法。具体做法是：先估计λ的上界λMax（根据对偶可行域有限性，λ确实存在上界），然后引入辅助变量Z_ij = λ_i × u_j，用一组线性约束包络住这个乘积：

Z_ij ≥ 0 Z_ij ≤ u_j^max × λ_i Z_ij ≤ λ_i^max × u_j Z_ij ≥ u_j^max × λ_i + λ_i^max × u_j - λ_i^max × u_j^max

这组约束的逻辑很简单：当u_j和λ_i都在各自的上下界内时，这四个不等式恰好把乘积的取值框在可行范围内。需要注意的是，这个线性化只有在λ和u都有界时才成立，因此需要对λ_i的上界做预估。我在实际项目中通常先用一个较小的上界试算，再检查求解结果中是否有对偶变量接近边界，若有则放大上界重新求解。

3.3 C&CG与Benders分解的取舍

很多教材会拿C&CG和Benders分解做对比。直观说，Benders分解是把子问题的对偶信息（割平面）逐步加入主问题，而C&CG是直接把最坏场景对应的完整变量和约束加入主问题。C&CG的优势在于收敛快很多，尤其是子问题含有整数变量时更是如此，因为C&CG不需要引入复杂的组合割。

我在微网容量配置中优先用C&CG的另一个原因是代码结构更清晰：主问题就是一个参数化的场景集合，子问题就是一个独立求解的max问题，两个部分各自好调试，出了问题也容易定位。如果模型规模很大（比如源网荷储一体化规划），我才会考虑Benders，因为这时主问题每轮迭代求解负担更重，Benders的轻量割平面可能整体更快。

4. Matlab代码落地：YALMIP建模与主循环实现

4.1 代码架构与文件规划

代码我建议按功能拆成五个文件，不要把所有逻辑堆在一个脚本里，否则调试时你会疯掉：

• main_ccg.m：C&CG主循环，负责迭代控制与结果输出
• load_system_data.m：读取系统参数，输出结构体par
• solve_MP.m：求解主问题
• solve_SP.m：求解子问题，返回最坏场景和运行成本
• build_run_constraints.m：构建第二阶段的运行约束（被MP和SP共用）

建模工具我用YALMIP，因为它的语法最接近数学表达式，省去大量矩阵拼接工作。求解器接口通过sdpsettings统一设置，方便在Gurobi和CPLEX之间切换对比。

4.2 主问题MP的建模实现

主问题MP的完整YALMIP代码如下。核心逻辑是：对每一个已发现的场景u_cell{i}，复制一份运行变量，构建对应的运行约束，并通过η约束保证主问题的目标值不小于所有场景的运行成本。

function [x_opt, eta, LB] = solve_MP(u_cell, par) K = length(u_cell); % 第一阶段容量变量 x = sdpvar(par.nGen, 1, 'full'); % 运行成本辅助变量 eta = sdpvar(1, 1); % 各场景运行变量（cell数组） y = cell(1, K); CostRun = zeros(1, K); Constraints = []; % 第一阶段约束：容量上下限与投资预算 Constraints = [Constraints, par.xMin <= x <= par.xMax]; Constraints = [Constraints, par.invCost' * x <= par.budget]; % 对每个已发现场景构建运行约束 for i = 1:K y{i} = sdpvar(par.nRun, 1, 'full'); [cons_i, cost_i] = build_run_constraints(x, y{i}, u_cell{i}, par); Constraints = [Constraints, cons_i]; Constraints = [Constraints, eta >= cost_i]; end % 目标：投资成本 + 最大运行成本（由eta代理） Objective = par.invCost' * x + eta; ops = sdpsettings('solver', par.mpSolver, 'verbose', 0); sol = optimize(Constraints, Objective, ops); x_opt = value(x); eta = value(eta); LB = value(Objective); if sol.problem ~= 0 warning('MP求解异常: %s', sol.info); end end

这里有一个细节值得强调：每轮迭代主问题都会新增一组变量和约束（对应新的最坏场景），导致主问题规模逐渐增大。如果迭代到七八轮之后求解开始变慢，可以考虑上一轮主问题的老场景约束其实已经不再起作用，通过解耦或热启动来加速，但这属于进阶优化，新手先确保逻辑正确再考虑提速。

4.3 子问题SP的对偶实现

子问题的完整实现比主问题复杂，因为需要处理对偶乘子和双线性项。我直接给一个经过验证的YALMIP代码框架，其中Z矩阵就是McCormick线性化引入的辅助变量。注意代码中的par.lambdaMax需要根据对偶可行域估计，这是整个实现中最容易出错的环节。

function [u_worst, obj_sp] = solve_SP(x_opt, par) % 不确定性变量 u = sdpvar(par.nUnc, 1, 'full'); % 对偶乘子：不等式约束对应lambda（非负），等式约束对应mu（自由） lambda = sdpvar(par.nIneq, 1, 'full'); mu = sdpvar(par.nEq, 1, 'full'); % 双线性项辅助变量 Z = sdpvar(par.nIneq, par.nUnc, 'full'); Constraints = []; % 对偶可行性约束 Constraints = [Constraints, par.A' * lambda + par.E' * mu == par.d]; Constraints = [Constraints, lambda >= 0]; % 不确定性集合：盒式+预算约束 Constraints = [Constraints, par.uMin <= u <= par.uMax]; Constraints = [Constraints, sum(abs(u - par.uNom) ./ par.uDev) <= par.Gamma]; % McCormick线性化：Z(i,j) = lambda(i) * u(j) M = par.lambdaMax; for i = 1:par.nIneq for j = 1:par.nUnc Constraints = [Constraints, ... Z(i,j) >= 0, ... Z(i,j) <= par.uMax(j) * lambda(i), ... Z(i,j) <= M(i) * u(j), ... Z(i,j) >= par.uMax(j) * lambda(i) + M(i) * u(j) - M(i) * par.uMax(j)]; end end % 对偶目标：线性部分 + 双线性部分 ObjLinear = par.b' * lambda + x_opt' * par.D' * lambda + par.e' * mu; ObjBilinear = sum(sum(par.C' .* Z)); Obj = ObjLinear + ObjBilinear; % YALMIP默认最小化，取相反数实现最大化 ops = sdpsettings('solver', par.spSolver, 'verbose', 0); sol = optimize(Constraints, -Obj, ops); u_worst = value(u); obj_sp = value(Obj); if sol.problem ~= 0 warning('SP求解异常: %s', sol.info); end end

代码里的par.C矩阵严格来说应该是一个与约束系数对应的张量结构，具体维度取决于你的模型定义。这里为了展示核心逻辑，把它简化为一个nIneq × nUnc的矩阵，实际工程中需要根据对偶约束的系数结构仔细拼装。我建议你在实现时把小规模算例对偶目标函数的数值和手算结果做对比，确认符号和转置方向没有错。

4.4 C&CG主循环与容差判断

主循环代码如下。初始化时我把u_cell的第一个元素设为预测值场景，这样第一次主问题求解就等价于一个确定性容量配置问题，给后续迭代提供一个可行的起点。

%% 两阶段鲁棒优化C&CG主循环 clear; clc; close all; addpath(genpath(pwd)); par = load_system_data(); par.Gamma = 2; par.tol = 1e-3; par.maxIter = 20; % 初始化：以预测场景作为第一个场景 u_cell{1} = par.uNom; LB = -inf; UB = inf; gap = inf; k = 1; while gap > par.tol && k <= par.maxIter fprintf('========== 第 %d 次迭代 ==========\n', k); % 求解主问题 [x_opt, eta, LB] = solve_MP(u_cell, par); % 求解子问题 [u_worst, obj_sp] = solve_SP(x_opt, par); % 更新上界：当前投资成本 + 最坏运行成本 UB_candidate = par.invCost' * x_opt + obj_sp; UB = min(UB, UB_candidate); % 计算相对间隙 gap = (UB - LB) / abs(UB); fprintf('LB = %.2f, UB = %.2f, gap = %.4f\n', LB, UB, gap); % 未收敛则把最坏场景加入主问题 if gap > par.tol k = k + 1; u_cell{k} = u_worst; end end %% 输出结果 disp('========== 鲁棒优化结果 =========='); disp('最优容量配置 = '); disp(x_opt); disp('总成本下界LB = '); disp(LB); disp('总成本上界UB = '); disp(UB);

注意一个数值细节：LB和UB的更新方向容易搞混。主问题是原问题的松弛问题（因为η只约束了已发现场景），所以主问题目标值给出的永远是一个下界；而"投资成本+子问题最优值"对应一个可行方案，给出的是上界。程序里我让LB直接取主问题目标值，UB取最小候选值，这个逻辑和C&CG理论完全一致。

5. 3机微网算例：收敛过程与鲁棒解对比分析

5.1 测试系统与参数设置

为了验证代码正确性，我用一个简化的3机微网系统做测试。系统包含光伏、风电、柴油发电机、储能电池以及与主网的联络线。典型日按24个时段建模，负荷曲线和风光出力预测曲线采用某园区夏季典型日数据，但不透露出处以免卷入数据合规问题。核心参数如下：

购电价格为分时电价，高峰时段1.2元/kWh，低谷时段0.4元/kWh。不确定性偏差设置为：光伏预测上下浮动20%，风电预测上下浮动30%，负荷预测上下浮动10%。不确定性预算参数Γ取2，表示所有时段偏差加权后总量被限制在两倍单时段偏差以内。

5.2 迭代收敛过程详解

运行主循环后，C&CG经历了5次迭代收敛，相对间隙降到0.3%以下。每次迭代的上下界变化如下表：

这里的成本口径是"全生命周期折算值"，即投资成本等年值加典型日运行成本折算到年。第一次迭代LB较低是因为主问题里只有一个预测场景，约束很少，松弛程度大。随着最坏场景不断加入，主问题可行域被逐步收紧，LB向上攀升，UB则保持在上界候选的最小值，两条曲线最终靠拢。这个单调收敛的模式是C&CG的正常表现，看到LB单调上升、UB基本持平，说明代码逻辑正确。

5.3 鲁棒决策与确定性决策的对比

为了体现鲁棒优化的价值，我把同样数据跑了一遍确定性优化（即固定不确定性为预测值），两个方案对比如下：

这个对比清晰地展示了两阶段鲁棒的代价与收益：投资多花约19.5%，但最坏情况下的运行成本下降31.3%，代价是平均运行成本小幅上升6.5%。业主可以根据风险偏好选择"中庸"的Γ值——例如Γ=1时投资约增加10%，最坏日运行成本下降18%左右。这套模型跑出来的是一整条"投资-保守度"帕累托曲线，而不是只有一个孤零零的最优点，这恰恰是两阶段鲁棒相对确定性优化最有说服力的地方。

6. 实战中踩过的坑：Big-M、求解器与预算参数调优

6.1 Big-M参数：并不是越大越好

子问题中McCormick线性化需要预先估计λ的上界，也就是代码里的par.lambdaMax。我最初偷懒给了一个全局统一的M=1e5，结果子问题求解时间飙升，而且几次迭代后数值明显不稳定。排查下来发现，过大的M会严重破坏约束矩阵的条件数，Gurobi在双重里耗费大量迭代。

我的经验做法分三步：第一步，先解一个去掉双线性项、只保留线性部分的子问题，记录λ的数值范围，通常就能获得一个合理的上界；第二步，将这个范围的10倍设为lambdaMax；第三步，求解完整子问题后检查λ值是否接近lambdaMax，如果触界，则逐步放大到50倍再试，直到没有变量触界为止。这套流程虽然多花一点时间，但能避免用拍脑袋的M值去赌数值稳定性。

提示：McCormick线性化只有在λ和u都有界时才严格成立，所以lambdaMax不能省。如果你发现即使把lambdaMax放大到很大，某些约束仍然被激活，那多半不是M的问题，而是对偶模型本身写错了。

6.2 求解时间爆炸的常见原因与优化

三机微网规模不大，单次MP和SP都只需秒级求解，所以迭代时间完全可控。但当你把系统扩展到十几台机组、采用多典型日联合规划时，主问题每轮会复制成千上万个变量和约束，求解时间会指数级恶化。我总结出三个实用的优化手段：

第一个是场景压缩。粗糙的历史数据往往包含上千个时段，我会先做场景聚类，把相似时段合并为不超过5到7个典型日，每个典型日再降采样为24个时段，这样能在保精度的前提下大幅压缩模型规模。第二个是热启动。上一轮主问题的最优解可以作为下一轮主问题的初始可行解传入求解器，显著减少分支定界的探索量。第三个是子问题并行——多典型日时各日子的最坏场景互相独立，可以并行求解，把total wall time压到原来的十分之一。

6.3 不确定性预算参数的经验选择

Γ是这套方法里最需要跟业主沟通的参数，因为它直接决定方案的保守程度。我给项目做参数标定时，一般流程是：先从历史数据中统计每个时段预测误差的分布，取85%分位数作为u_dev；然后让Γ从0到T（T是时段数）以0.5为步长扫描，画出"总投资-最坏运行成本"的帕累托曲线；最后和业主确认风险承受能力，选择曲线上拐点位置对应的Γ值。

集群算例的经验表明，Γ取T的15%到25%通常是一个合理区间。太小起不到鲁棒作用，太大则方案成本激增且几乎不再有意义地降低最坏运行成本。还有一个容易忽视的点：Γ的取值要结合实际情况审核。比如光伏在夜间出力本来就应该为0，你不能在夜间时段也给它设置20%的波动偏差，否则子问题会生成一个物理上不成立的"夜间光伏满发"场景，导致容量配置虚高。解决办法是在不确定性集合里给每个时段乘一个出力上限系数，光伏夜间强制为0。

我在实际运行中还发现，两阶段鲁棒优化的结果对不确定性集合的形状很敏感，但对求解器选择的敏感性相对较低。Gurobi和CPLEX在同一个模型上给出的最优解完全一致，差异主要体现在求解时间上。真正影响结果质量的，永远是模型的物理准确性和不确定性参数的取值依据。在这两点上多花时间，比反复调求解器参数要有价值得多。

最后再分享一个小技巧：所有C&CG代码里，我强烈建议在第一轮迭代时手动输入一个"已知的极端场景"作为初始u_cell元素，而不只是用预测场景。比如历史记录里真实发生过的"连续阴雨+负荷尖峰"那天的出力曲线。这样做能让第一次主问题的解更接近最终最优解，减少两三轮迭代，更重要的是能帮你快速验证子问题找出的最坏场景是否合理——如果这么明显的极端场景子问题都没找到，说明你的对偶实现大概率有bug。