两阶段鲁棒优化在微网多电源容量配置中的MATLAB实现与CCG求解

简介：这是一份面向微电网规划与电力系统优化方向的两阶段鲁棒优化容量配置Matlab代码包，适合电气工程、自动化及数学相关专业学生用于课程设计、毕业设计或算法复现。代码基于参数化编程，提供2014/2019a/2021a版本兼容的M文件与运行结果，并附带案例数据，可快速调整参数并观察鲁棒优化在风电、光伏多电源容量配置中的实际效果。压缩包共427个文件，主要包含276个xls数据表、110个mat数据文件、5个m源码脚本以及csv、docx说明文档等，包体约91.62MB，目录层次清晰，方便按模块查阅。资源已有105人学习，作者为资深算法工程师，附赠可运行样例和结果，能帮助读者理解两阶段鲁棒优化的建模、求解与结果分析全过程。

1. 两阶段鲁棒优化在微网容量配置里到底在解决什么问题

风电和光伏装得越多，微网规划时那个经典矛盾就越明显：按确定性负荷曲线来配，投资方案很“划算”，可一旦遇上连续阴雨天叠加低风速，备用容量不足就只能切负荷；把容量往大了配，成本又压不住。两阶段鲁棒优化不猜风光的准确出力，而是把不确定量放进一个有边界的集合里，先用投资决策锁定容量，等恶劣场景“揭晓”后再做运行调度。这种 two-stage 结构给出的多电源容量配置，能在最恶劣的可预料场景下守得住可靠性底线，代价只是多付出10%左右的成本。下面把两阶段鲁棒优化算法的微网多电源容量配置全过程拆开讲，从不确定集、min-max-min 数学模型，到 MATLAB 代码实现、C&CG 迭代求解、运行结果指标解读，适合做微网规划、储能容量论证和相关课题复现的工程师直接对照调整。

2. 两阶段鲁棒优化模型：不确定集与 min-max-min 三层结构

2.1 微网容量配置里的不确定量怎么描述

微网多电源容量配置要回答的问题是：风电机组装多少、光伏板装多少、储能电池配多大、燃气轮机留多少余量，才能在经济性和可靠性之间取得平衡。负荷的不确定性、风电和光伏的出力波动，是这个问题最核心的输入。

这里把风电、光伏和负荷统一写成一个不确定参数向量u = [P_w, P_pv, P_L]^T，用盒式不确定集描述：

U = { u | u_L ≤ u ≤ u_U }

盒式集合的好处是线性约束，调 MATLAB 的linprog就能处理，不需要引入二阶锥或半定松弛，迭代效率高。更细一点的做法是加一个 1-范数预算约束（budget constraint），限制风光同时达到极端值的场景数量，避免“所有不确定量同时取最差值”这种在现实中几乎不出现的组合把配置推得过度保守。

预算约束的数学形式长这样：

∑ | u_i - u_i^c | / Δu_i ≤ B

其中u_i^c是预测中心值，Δu_i是偏离幅度，B是预算参数。B=0退化回确定性模型，B=n意味着允许所有维度的不确定量同时到达边界。工程上B取总维度数的三分之一到二分之一，能在鲁棒性和经济性之间取一个平衡点。很多刚开始做鲁棒优化的同学直接不写预算约束，把盒式集合当成全部，得到的容量配置常常保守到没法用，这个细节值得重视。

2.2 为什么要拆成“两阶段”而不是直接求一个超大优化

容量配置里存在两类决策变量：第一阶段是投资决策，包括风机、光伏、储能和燃气轮机的安装容量，记为x；第二阶段是给定不确定场景下的运行决策，包括各机组出力、储能充放电功率、切负荷量，记为y。

确定性模型写成：

min c^T x + d^T y s.t. Ax ≤ b Bx + Cy + Du ≤ e

确定性模型里的u是固定值，所以只有“一个场景”，求出来是一组确定的装机容量和一组确定的运行方案。两阶段鲁棒模型换成 min-max-min 结构：

min_x ( c^T x + max_{u∈U} min_{y∈F(x,u)} d^T y )

内层的 min 是给定投资方案和某个场景下，找使运行成本最小的调度方式；中间的 max 是在不确定集里找使这个最优运行成本最大的场景；外层的 min 是调整投资方案，使“投资成本 + 最恶劣场景运行成本”的总和最小。

外层x的决策要先于u揭晓，所以叫“这里现在”阶段；内层y的调度在u已知后执行，所以叫“等待并观察”阶段。这种拆法的工程含义很直观：先花钱建电站、把储能容量定下来，之后无论老天爷给出多差的天气，微网都能通过调度兜住负荷。如果直接对所有可能场景做确定性优化，每一个连续不确定变量都会让问题规模指数膨胀，求解器会直接被压垮。两阶段结构把不确定量从决策层中剥离，交给列与约束生成算法（C&CG）在有限次迭代里求解，计算量大幅缩小，这也是这个模型能落到 MATLAB 上的根本原因。

2.3 与单阶段鲁棒、随机规划的差别

随机规划要求事先给出不确定量的概率分布，对分布误差很敏感。微网里的风电、光伏历史数据往往只有两三年的覆盖，分布估计得并不准，你用的正态分布或贝塔分布可能跟实际情况差距很大。两阶段鲁棒优化不需要分布，只需要知道取值区间。在实际规划中，这个区间可以从历史最小值/最大值或者分位数来定，数据需求少、边界清晰，比基于分布假设的方案更稳妥。

单阶段鲁棒把所有变量捆绑在一起做最坏情况优化，结果是“投资决策和运行决策都按最坏场景算”，过度保守。两阶段只要求存在一个可行的运行方案能应对最坏场景，不要求某一组运行决策在所有场景下都可行，经济性明显更好。举例：单阶段模型会把储能功率严格按“最坏时刻需要最大充放电”来设计，两阶段模型则允许在坏场景下额外切一点可中断负荷，储能的装机就能降下来，这就是两阶段这个设计最核心的工程价值。

3. MATLAB 实现：从参数定义到 C&CG 迭代求解

3.1 主问题 MP：把运行成本看作一个待逼近的变量

两阶段鲁棒不能用常规求解器直接求解 min-max-min，因为中间那层 max 没有一个可微的目标让你写梯度。C&CG 算法的思路是迭代逼近：先给一个初始场景，求解主问题 MP，再把主问题的x拿出来给子问题 SP 找最恶劣场景，找到后把场景加入主问题再解，反复进行。

主问题 MP 写成：

min_{x, η, y_k} c^T x + η s.t. Ax ≤ b Bx + C y_k + D u_k ≤ e, k=1..K d^T y_k ≤ η, k=1..K

这里η是对最恶劣场景运行成本的逼近值。每轮迭代识别出的场景u_k是已知常数，因此 MP 是个带参数K个运行块的标准线性规划。K越大，η被约束得越紧，MP 的目标值越接近真实的两阶段鲁棒最优值。

% 主问题 MP 建模（YALMIP 语法） % x: 第一阶段投资决策（连续变量） % y{k}: 第 k 个已识别场景下的运行决策 % eta: 标量辅助变量，逼近最恶劣运行成本 K = size(u_list, 2); % 已识别场景个数 Constraints = [A * x <= b]; % 投资约束 Cost_MP = c' * x + eta; % 总成本 = 投资成本 + 运行逼近值 for k = 1:K Constraints = [Constraints, B * x + C * y{k} + D * u_list(:, k) <= e]; Constraints = [Constraints, d' * y{k} <= eta]; end ops = sdpsettings('solver', 'gurobi', 'verbose', 0); % 也可用 cplex 或 linprog optimize(Constraints, Cost_MP, ops); x_mp = value(x);

这段代码的关键逻辑是：u_list里存的是所有已被识别为“恶劣”的场景，每增加一列，就多一组运行变量y{k}和一条约束d'*y{k} <= eta。eta被迫取这些场景运行成本的最大值，所以 MP 的目标值是从下往上逼近真实最优值的。如果子问题找出的最恶劣成本大于当前eta，说明 MP 低估了风险，必须把这个场景加入u_list再解。

3.2 子问题 SP：固定 x 之后找最恶劣场景

给定x = x_mp后，子问题是内层 min 和外层 max 的嵌套。处理方法是对内层运行问题取对偶，把 min 换成 max，从而把整个子问题化简成单层最大化问题。内层运行问题是一个线性规划，强对偶条件成立，所以对偶转换没有损失精度。

对偶后的子问题里会出现双线性项——对偶变量乘不确定量，这使 SP 在一般形式下是 NP-hard 的。但有一个关键简化：在盒式不确定集下，最大值一定出现在某个顶点上。这意味着不需要搜索连续空间，只需枚举2^n个顶点（n是不确定量个数）就能找到最恶劣场景。

% 子问题 SP：顶点枚举法找最恶劣场景 % x_mp: 由 MP 求得的投资容量 % u_L, u_U: 不确定量下界/上界向量 n_u = length(u_L); worst_cost = -Inf; u_worst = u_L; % 遍历所有顶点组合，2^n_u 个候选 for idx = 0:(2^n_u - 1) u_try = u_L; for j = 1:n_u if bitget(idx, j) % 第 j 维取上界 u_try(j) = u_U(j); end end % 固定 u_try，解内层运行 LP，得到该场景下最小运行成本 [cost_val, ~] = solve_operation(x_mp, u_try); if cost_val > worst_cost worst_cost = cost_val; u_worst = u_try; end end

参数说明：bitget(idx, j)用于把循环索引idx的二进制位拆出来，决定当前顶点在第j维取上界还是下界。solve_operation是内层运行优化函数，输入投资容量和不确定场景，输出该场景最小的运行成本。这个函数通常只包含功率平衡、储能递推、机组出力范围等运行约束，用linprog或quadprog就能解。

顶点枚举的前提是目标函数关于u的极值落在盒子顶点上。对偶问题里双线性函数关于u是线性的，确实在顶点取极值，但当不确定量维度超过 12 时，2^n的数量级会非常可观，枚举法就不适用了。那种情况建议改用 big-M 线性化把顶点选择写成混合整数线性规划，或者把不确定集里的连续变量保留，用 Benders 风格的对偶切割逼近。

3.3 C&CG 整体迭代框架

C&CG 的完整流程：用不确定集中心u_c = (u_L + u_U)/2初始化u_list，求解 MP 得到(x_mp, eta_val)，固定x_mp求解 SP 得到(cost_worst, u_new)，计算上下界间隙，如果不满足收敛条件就把u_new加入u_list，重复迭代。

数值上，MP 的目标值是下界LB = c^T x_mp + eta_val，SP 给出的c^T x_mp + cost_worst是上界UB，两者之差就是间隙。用相对间隙做停止准则比绝对阈值更可靠，因为容量配置问题的目标值跨度很大，几百万和几千万的数量级差异下，绝对阈值没有可比性。

% C&CG 主循环 tol = 1e-3; % 相对间隙阈值 max_iter = 30; K = 1; u_list(:, 1) = (u_L + u_U) / 2; % 初始场景取中心 for iter = 1:max_iter % 1) 求解主问题 [x_mp, eta_val] = solve_mp(u_list); % 2) 固定 x_mp，求解子问题，得到最恶劣场景 [cost_worst, u_new] = solve_sp(x_mp); % 3) 计算上下界与相对间隙 LB = c' * x_mp + eta_val; UB = c' * x_mp + cost_worst; rel_gap = abs(UB - LB) / max(abs(UB), 1); % 4) 收敛则退出，否则把新场景加入 u_list if rel_gap <= tol fprintf('收敛于迭代 %d 次, 目标值 %.2f 万元\n', iter, UB); break; end K = K + 1; u_list(:, K) = u_new; end

这个循环的关键在收敛判断这一步。UB - LB趋向零说明 MP 对最恶劣场景成本的逼近已经足够准，再加新场景也不会让目标值明显变化。实践中如果迭代 20 次 gap 还不降，检查两个地方：一是子问题里运行可行性约束是否完备，比如储能 SOC 递推有没有加错时间下标；二是u_new是否真的比已有场景更恶劣，如果cost_worst小于等于当前eta，说明主问题已经识别了所有关键场景，gap 数值上的微小抖动是求解器数值误差造成的，不是算法问题。

3.4 求解器配置与关键参数表

MATLAB 里做两阶段鲁棒优化，最常用的组合是 YALMIP + Gurobi 或 YALMIP + Cplex。YALMIP 负责把模型翻译成求解器需要的标准形式，Gurobi 负责高性能求解线性规划和混合整数规划。如果只是课程作业级的小算例，用 MATLAB 自带的linprog也完全可以，但注意linprog默认对偶单纯形法在重复迭代场景下比较慢。

B这个参数影响最明显：B越小鲁棒性越弱，B越保守，备用的储能容量与燃气轮机功率跟着涨。边界u_L / u_U的取法也很讲究。直接用历史最大最小值，边界会被极端尖峰拉得很宽，最恶劣场景在数学上合法但现实中基本不会出现；用 5% 和 95% 分位数更符合工程实践。

提示：不要为了“看起来更鲁棒”而把不确定集边界设到历史极值。鲁棒优化的价值是控制尾部风险，不是消除一切风险。边界过宽导致投资成本虚高，评审或决策层看到成本翻倍后，往往会直接否定这个技术路线。

4. 多电源容量配置建模与运行结果分析

4.1 风机、光伏、储能、燃气轮机的决策变量与约束

多电源容量配置的“多”一般落在四类电源上：风电装机（风机台数乘单机容量）、光伏装机（峰值功率）、储能（额定容量与最大充放电功率）、燃气轮机（额定功率）。第一阶段决策变量写为：

x = [ N_w, P_pv, E_bat, P_ch, P_dis, P_gt, P_grid ]

N_w是风机台数，P_pv是光伏峰值功率，E_bat是储能容量，P_ch / P_dis是储能充放电功率上限，P_gt是燃气轮机额定功率，P_grid是微网与外部电网的联络线功率上限。

约束分成两类。投资约束：

N_w ≤ N_w_max P_pv ≤ P_pv_max E_bat ≤ E_bat_max

这些上限来自可用土地、屋顶面积和资金上限。运行约束是每个时段的功率平衡：

P_w(t) + P_pv(t) + P_gt(t) + P_dis(t) + P_grid(t) = P_L(t) + P_ch(t)

储能 SOC 递推是多电源系统里最容易写错的地方：

SOC(t+1) = SOC(t) + η_ch * P_ch(t) * Δt - P_dis(t) / η_dis * Δt 0 ≤ SOC(t) ≤ E_bat

充放电功率同时受上下限约束，且同一时刻只能处在一种状态。实际建模中常常用二进制变量来防“既充电又放电”，但这会把运行子问题从 LP 变成 MILP，增加求解负担。一个折中是忽略同充同放约束，让线性规划天然在最优解里避开这种浪费行为——功率平衡约束加上成本为正的前提，使它不会同时充电和放电。

燃气轮机约束包括出力上下限和爬坡率：

0 ≤ P_gt(t) ≤ P_gt | P_gt(t+1) - P_gt(t) | ≤ r_gt * Δt

联络线约束考虑微网与主网交换功率的上限，以及关口功率不反向倒送的场景约束。

4.2 典型日与逐时出力曲线怎么嵌入两阶段框架

容量配置的工程实践里，不能拿全年 8760 小时逐时数据去套两阶段模型，那样不确定集的维度巨大，C&CG 迭代几十轮也收敛不了。常见做法是选典型日：春夏秋冬各挑一个典型日，再加上从历史数据中提取的极端日，比如持续无风加阴雨 48 小时的场景，把这些天的逐时曲线折合成 24 小时序列作为场景基架。

% 典型日 24 小时负荷/光伏/风电出力率（标幺值） load_curve = [0.72 0.68 0.65 0.63 0.64 0.68 0.74 0.85 ... 0.95 1.00 0.97 0.90 0.86 0.88 0.92 0.96 ... 1.05 1.10 1.08 0.98 0.88 0.80 0.75 0.70]; pv_curve = [0 0 0 0.02 0.08 0.16 0.30 0.45 ... 0.60 0.70 0.72 0.68 0.62 0.55 0.42 0.28 ... 0.15 0.05 0 0 0 0 0 0]; wind_curve = [0.35 0.38 0.40 0.36 0.32 0.30 0.34 0.40 ... 0.44 0.38 0.36 0.40 0.42 0.39 0.35 0.33 ... 0.36 0.44 0.50 0.48 0.45 0.42 0.38 0.35]; % 各时段风电/光伏实际出力 = 装机容量 * 出力率 P_w(:) = N_w * wind_curve(:); P_pv(:) = P_pv_cap * pv_curve(:);

把典型日的曲线作为运行层输入的确定性部分，不确定集再对曲线做整体缩放，就能保留“一天内出力形状固定、但整体水平不确定”的物理特征。实际项目中常见做法是给典型日曲线乘一个[0.85, 1.15]区间的缩放系数，并用预算约束限制同时达到高比例波动的时段数。这样运行层仍保持线性，子问题顶点枚举的维度等于不确定场景数，计算量完全可控。

4.3 运行结果指标：装机容量、弃风弃光率与鲁棒代价

运行结果通常输出三类指标：各类电源最优装机容量、最恶劣场景下的运行成本、C&CG 迭代收敛曲线。

以日峰值负荷 1 MW 的中等规模微网为例，两阶段鲁棒配置的典型结果是：风电和光伏装机容量略低于确定性优化结果，储能容量明显偏高，燃气轮机保持一个兜底容量不变。原因是风光在恶劣场景下贡献比例低，多装只会推高投资成本；储能虽然也贵，但连续阴天场景里确能承担调峰作用，可靠性和经济性的边际平衡点在容量配置里非常清晰。

运行结果中还有一个指标值得关注：鲁棒代价（price of robustness），定义为两阶段鲁棒解的目标值与确定性优化目标值的比值。它通常落在 1.0~1.25 之间，意思是：为了守住最恶劣场景的可靠性，需要在正常场景下多付出 0~25% 的成本。

这张表反映的规律很典型：风电光伏因为最恶劣场景下出力低，鲁棒解选择少装，把省下的投资额度补给储能和燃气轮机。弃风弃光率在两类方案里差异并不大，真正的差异体现在切负荷概率——确定性方案在极端场景下可能切掉 8% 的负荷，而鲁棒方案把切负荷率压到了 0。这就是两阶段鲁棒优化在微网容量配置里最核心的价值：它不是让方案更省钱，而是让方案在极端天气下仍然撑得住。

5. 用蒙特卡洛回代检验验证鲁棒解的质量

5.1 回代检验步骤与 MATLAB 实现

两阶段鲁棒优化在理论上保证盒式不确定集内最恶劣场景可行，但盒式集始终是对真实不确定性的近似，边界参数设得合不合理，需要用历史分布回代检验。做法是：从历史风光负荷数据中估计每时段的均值与标准差，用蒙特卡洛生成大量随机场景，把已求得的装机容量代入运行仿真，统计切负荷率和弃风弃光率。

% 蒙特卡洛回代：统计切负荷概率 LOLP 与期望缺供电量 EENS n_sample = 1000; lolp_count = 0; eens_sum = 0; for s = 1:n_sample % 从历史分布抽样生成 24 小时曲线 load_s = load_mean + randn(24,1) .* load_std; pv_s = max(0, pv_mean + randn(24,1) .* pv_std); wind_s = max(0, wind_mean + randn(24,1) .* wind_std); % 代入容量配置，运行日调度 loss = run_day_schedule(load_s, pv_s, wind_s, config); if loss > 1e-6 lolp_count = lolp_count + 1; end eens_sum = eens_sum + loss; end fprintf('LOLP = %.2f%%, EENS = %.2f kWh\n', ... lolp_count / n_sample * 100, eens_sum / n_sample * 24);

回代结果直接影响不确定集参数的设计：如果 LOLP 高于可靠性目标，说明不确定集取窄了，应放大u_L / u_U区间或增加预算B；如果 LOLP 远低于目标且投资成本虚高，则应收窄边界。回代检验是针对结果的验证，同时也是校准参数的手段。

5.2 三个容易出错但必查的边界条件

第一个是储能 SOC 终端约束。两阶段鲁棒模型若只约束每个时段的 SOC 在 0 和容量之间，不要求日末回到初值，优化结果可能在 24 小时结束时留下大量剩余电量，导致储能容量配置偏小。修正方法是在约束里加SOC(24) ≥ SOC(0)，等价于“储能不能消耗初始电量”。

第二个是切负荷罚系数。子问题目标里的切负荷成本必须远高于正常供电成本，常见做法是取电价的 20~50 倍，否则优化会故意切负荷来降低总投资。但罚系数也不宜超过 1e6，过大会让求解器矩阵条件数恶化，Gurobi 报数值奇异，目标值出现异常抖动。

第三个是收敛判据的数值噪声。MP 每加一个场景，eta只会上升不会下降，但 SP 的相对间隙在迭代后期可能在 1e-4 附近波动，原因是单纯形法在不同基解之间跳变，切负荷量少数值上的差异导致目标出现微小振荡。这类波动不是不收敛，只要相对间隙低于预设阈值即可停止迭代。

判断一个两阶段鲁棒解是否合格，最终看回代仿真里的切负荷率，而不只是看 C&CG 是否收敛。容量配置的问题本质上是给不确定性定价，两阶段鲁棒的价值是算出“最坏情况下需要多少备用”，而这个备用是否值得投资，要用蒙特卡洛回代结合全年运行经济性来最终回答。

本文还有配套的精品资源，点击获取