★ FEATURED ARTICLE 拿到这个复现任务很多人第一反应是把论文里的模型公式抄下来然后直接上 MATLAB 调 YALMIP 写求解。但真正动起手来你会发现问题远不止“敲代码”这么简单——鲁棒优化里的 max-min 结构不是随手就能表达出来的两阶段问题的收敛判断、不确定集里那个预算参数的调节、储能 SOC 在 24 个小时段间的耦合约束每一个环节都可能卡住你半天。这篇文章把我从零开始复现“含可再生能源和储能的区域微电网最优运行多阶段鲁棒调度模型”的完整过程、踩坑点和可复用的代码思路都梳理出来。内容会覆盖为什么确定性调度不够用、多阶段鲁棒模型到底在“鲁棒”什么、CCG列与约束生成算法怎么一步步落地以及 MATLAB 里整套程序怎么组织。适合正在做微电网优化调度、综合能源系统研究或者毕业论文里需要把“鲁棒优化”这个抽象概念真正落到代码上的同学参考。1. 项目整体设计与思路拆解1.1 为什么“确定性最优调度”在微电网里不够用先看一个最简单的例子。某区域微电网里有一台燃气机组、一组光伏、一组储能通过一条联络线与上级电网相连。未来 24 小时的负荷预测曲线、光伏出力预测曲线都已经拿到按照确定性模型求解目标函数是运行成本最小约束是每个时段的功率平衡、机组出力上下限、储能容量限制。这个模型用 MATLAB 几分钟就能解出来结果也“合理”白天光伏多的时候少买电夜间光伏为零就多买电机组只在电价高或负荷大的时段出力。但问题是预测值只是“预测”。光伏预测出力是 100 kW实际可能因为云层遮挡只有 60 kW。如果调度计划是按照 100 kW 的出力去制定的那实际运行时系统就会出现功率缺额。缺额从哪里补只能由储能多放电、机组多出力、或者向主网多买电。可机组和联络线都有容量上限储能也有 SOC 约束如果日前计划已经把这些资源压得很满实际执行时根本调整不过来最终只能切负荷。切负荷的代价无论从失负荷惩罚成本还是从供电可靠性角度都比多烧一点燃气大得多。所以一个真正可执行的微电网调度方案必须把“不确定性来源”直接放进模型里而不是先在确定性世界里求出一个完美计划再指望现场靠运气扛过去。1.2 随机规划与鲁棒优化的取舍处理不确定性有两类主流范式随机规划和鲁棒优化。随机规划的思路是为每个不确定参数构造若干场景并为每个场景赋一个概率模型在期望意义下做决策。它的好处是结果不过分保守而且能利用历史数据的统计信息。但它有个绕不开的痛点你很难准确给出光伏出力的概率分布。气象预测误差分布在不同地区、不同季节、不同天气类型下完全不一样拿一个假设好的正态分布去套算出来的方案往往“看着科学用着心虚”。鲁棒优化的思路完全不同。它对不确定参数的建模只需要一个集合比如“光伏出力在预测值上下 20% 范围内波动”“所有时段总波动幅度不超过某个预算值”。模型不关心这些波动发生的概率是多少只要求在集合内的任何波动之下系统都能通过第二阶段调整维持安全运行。它问的是“最坏情况下会怎样”并且为这个最坏情况提前留好调整空间。这两种方式用一句话区分随机规划管住“大概率的风险”鲁棒优化兜住“最坏的情况”。对微电网这种量级小、调节手段有限、对可靠性要求高的系统鲁棒优化显然更贴合工程直觉。1.3 多阶段模型的真实含义标题里的“多阶段”需要说清楚。它并不是指把 24 小时切成很多段然后依次滚动求解这里的“多阶段”在结构上更接近一个大闭环第一阶段做一个“事前决策”比如机组启停、与主网的购售电计划、储能的基线出力这个决策必须在不确定性实现之前定下来第二阶段是“事后调整”比如燃气机组实际出力、储能实际充放电功率、不得已情况下的弃风和切负荷这个决策可以在看到真实风光出力之后做。工程上理解第一阶段相当于前一天晚上做好的“计划”第二阶段相当于当天根据实际来风、来光进行的“实时修正”。储能作为微电网里最灵活的调节资源往往横跨两个阶段它的日前出力计划是第一阶段的决策但实际 SOC 轨迹会根据实时情况微调。这个结构决定了储能建模不能只写一条 SOC 递推公式还要区分“计划值”和“实际值”。这个设计的好处是避免调度方案“过度保守”如果所有决策都必须在风光出力划定之前拍板那就等于把全系统所有调节能力都绑死在最坏风光场景上经济性会很差。两阶段结构让机组启停和购电计划保持相对稳定同时把出力调整级、储能微调级留给日内实时平衡既保证鲁棒性又不牺牲太多经济性。2. 核心数学模型搭建与细节处理2.1 优化目标从运行成本到惩罚项任何调度模型的出发点都是目标函数。这个复现项目里目标由四类成本组成第一是微型燃气轮机的燃料成本通常用出力的二次函数近似在 YALMIP 里用二次约束可以直接表达但如果想让模型严格保持混合整数线性规划MILP就把二次函数做分段线性近似用 SOS2 约束或 big-M 逻辑实现。第二是向上级电网购电的成本按分时电价计如果允许向主网倒送电还可以考虑售电收益项但通常反向潮流的价格远低于正向购电写不写都不影响主要结论。第三是弃风弃光惩罚成本。这个惩罚项的存在很重要鲁棒优化的调度结果里为了应对最坏场景可能会出现“把上网功率压得很低”的保守策略如果没有惩罚模型会倾向于大量弃掉便宜的可再生能源虽然安全但经济性很差。给弃电设一个适中的惩罚单价比如可再生能源上网电价的一半左右模型才会在安全性和利用率之间做出合理权衡。第四是切负荷惩罚成本。这个惩罚单价需要设得比购电价、燃料成本高一个数量级体现供电可靠性的优先级。我复现时把切负荷惩罚设为每千瓦时 5 到 10 元正常不会走到切负荷这个选项但模型需要保留这条路避免因为约束太紧直接无解。两阶段鲁棒目标的标准写法是min_x c^T x max_{u in U} min_{y in F(x, u)} d^T y外层 min 对应第一阶段成本中间 max 是在不确定集 U 中寻找最坏场景内层 min 是在看到该场景后选择最优调整方案。这个嵌套结构是整个项目的核心难点也是 CCG 算法要解决的东西。2.2 决策变量分类清楚划分两阶段决策变量是后续建模不混乱的前提。我复现时按下面这张表划分变量类别所属阶段包含内容第一阶段变量 x事前决策燃气机组启停状态、出力基线、向主网购电/售电功率、储能预调度充放电功率、储能充放电状态 0-1 变量第二阶段变量 y事后调整实际出力偏差、实际储能充放电调整量、弃风电功率、弃光电功率、切负荷功率、功率不平衡松弛量不确定参数 u外部输入风电出力偏差、光伏出力偏差特别提醒一点第二阶段变量如果包含 0-1 整数变量求解复杂度会急剧上升因为子问题里 max-min 结构无法直接通过强对偶定理转成单层优化。复现这个项目时我建议把第二阶段变量全部设计成连续变量第一阶段变量才允许含整数。这也是很多经典论文里两阶段鲁棒模型的通用做法。2.3 不确定集构建用预算参数控制保守度不确定集的形状决定了鲁棒模型的保守程度是整套模型里最“艺术”的地方。我复现时采用最常见的盒式加预算约束的不确定集U { u_w u_pv : -delta_w u_w delta_w -delta_pv u_pv delta_pv sum( |u_w| / delta_w |u_pv| / delta_pv ) Gamma }前两行把每个时段的风电、光伏出力偏差限制在预测值附近的区间内第三行是整个调度周期内的总偏差“预算”。Gamma 就是那个控制保守度的旋钮Gamma 越大要求系统应对的总偏差越大方案越保守Gamma 0 时退化为确定性模型。从工程经验看Gamma 的取值不需要从最大的开始试。我通常先跑 Gamma 0 的确定性方案作为下界然后从 Gamma 1 开始逐步增加到 6 或 8。当成本随 Gamma 变化的斜率突然变大的时候就说明再增加保守度要多付出很高的经济代价这个拐点附近的 Gamma 就是比较合适的运行保守水平。还有一个细节绝对值约束在 YALMIP 中要用 epigraph 技巧展开引入一个非负辅助变量 z把 |u_w| / delta_w 转成 z_w 再加两个不等式约束。YALMIP 里可以直接写 norm 或者 abs但求解器底层会做同样的展开。2.4 储能与功率平衡约束储能是模型中约束最多的元件每一个误区都会导致无解或者结果毫无物理意义。先看 SOC 递推约束SOC(t1) SOC(t) eta_ch * P_ch(t) * dt - P_dis(t) * dt / eta_dis其中 eta_ch 和 eta_dis 分别是充、放电效率dt 是时段长度。这里要注意充放电效率不要用一个笼统的 eta 去替代否则储能“存一度电再放一度电”的效率损失会算错。充放电互斥约束是另一个高频坑点。储能不能同时充电和放电需要引入 0-1 变量 z_ch0 P_ch(t) P_ch_max * z_ch(t) 0 P_dis(t) P_dis_max * (1 - z_ch(t))这套 big-M 写法把互斥关系用变量上界做了绑定。第一阶段的 0-1 变量数量是 24 个对 MILP 而言不算多。功率平衡约束必须体现“不确定参数”——这就是两阶段模型和确定性模型的核心差异P_g(t) P_w(t) P_pv(t) P_dis(t) P_buy(t) P_load(t) P_ch(t) P_sell(t)式中P_g、P_dis、P_buy、P_ch、P_sell 是待求变量P_w 和 P_pv 中含不确定偏差 u_w、u_pv。因为 P_load 是预测值也被外界数据误差影响但复现时通常只把可再生能源出力作为不确定来源负荷偏差可以并入切负荷项处理。功率平衡约束里的每一项都有物理边界这些边界不能拍脑袋写。联络线购电功率上限由变压器容量决定机组出力上限看铭牌储能充放电功率上限是由双向变流器决定的跟电池容量不是同一个概念。3. 求解算法CCG 分解与子问题转化3.1 两阶段问题的整体结构两阶段鲁棒模型天然不适合“整个丢给求解器”。求解器能直接处理的是一个单纯的 MILP而这里出现的是“min-max-min”三层嵌套结构。要解它必须用分解算法把它拆成一个大问题加一个小问题交替迭代求解。主问题Master Problem是带有若干“已知坏场景”的第一阶段优化问题。它的想法是我暂时不知道完整的不确定集合但我把若干最坏场景直接作为约束放进第一阶段模型中要求方案在所有这些场景下都可行。场景越多主问题越接近真实鲁棒模型的可行域同时目标函数值也会不断上升。子问题Subproblem是“给定一个第一阶段方案找出让系统最难受的那个不确定场景”并同时求出该场景下的最小运行调整成本。子问题输出的场景会被加入主问题形成迭代。这个思路跟 Benders 分解很像但 CCG 每一步是往主问题里添加“新的场景变量和对应约束”而不是添加割平面。很多资料里把 CCG 译作“列与约束生成”就是这个原因。3.2 CCG 迭代流程以 24 时段调度为例CCG 的总体框架是第一步初始化。给不确定集参数赋好值设置迭代次数 k 1设置上界 UB inf下界 LB -inf设置收敛精度 epsilon。第二步求解主问题。主问题里包含了前 k-1 轮迭代加入的所有场景。解出第一阶段决策变量值 x_*和目标函数值 obj_MP。因为主问题只是原问题的松弛版本所以 LB max(LB, obj_MP)。第三步固定 x_* 到子问题。子问题要在不确定集内搜索最坏场景并给出最坏场景下的最小调整成本。子问题目标函数值就是原问题里“max-min”那部分的值它加上第一阶段成本 c^T x_就构成原问题的一个上界候选值。更新 UB min(UB, c^T x_ obj_SP)。第四步收敛判断。如果 (UB - LB) / abs(UB) epsilon迭代结束否则把子问题中找到的最坏场景 u_* 添加为主问题的一个新场景即新增一组第二阶段变量 y_k并把它和第一阶段变量 x 之间的耦合约束加入主问题k k 1回到第二步。整个迭代过程在 MATLAB 里就是一个 while 循环循环条件很直接while (UB - LB) / abs(UB) tol iter MaxIter实操时我会把 tol 设为 1e-3 或 1e-4。按我的经验这个模型通常迭代 3 到 8 轮就能收敛不会出现跑几十轮还不收敛的情况。如果出现很可能模型有 bug而不是算法本身的问题。3.3 子问题的对偶转换技巧子问题本身也是“max-min”结构不能直接交给求解器。但如果第二阶段的决策变量全是连续变量且所有约束都是线性的那么内层的 min 问题就是一个线性规划。根据强对偶定理线性规划的最优值等于其对偶问题的最优值。这一步是整个复现工作中技术含量最高的地方。具体操作是把内层 min 先写成标准形式然后按规则写出对偶问题等式约束对偶变量自由不等式约束对偶变量非负原问题求 min对偶问题求 max。对偶之后内层的 min 变成了 max外面本来就是一个 max两个 max 合并成一个单层 max 问题。合并后的单层问题目标函数里会出现不确定变量和双层变量的乘积项也就是 u 和拉格朗日乘子 λ 相乘。这种双线性项是非凸的不能直接交给 CPLEX 或 Gurobi 求解。标准处理手段是用大 M 法做线性化引入辅助变量替代 u * λ 的乘积并通过 M 值和逻辑约束将这个辅助变量钉在正确的范围内。M 的取值需要足够大但不能大到引发数值病态我通常取相应约束右侧系数的 10 到 100 倍跑出来是稳定的。这里有一个取舍如果你觉得大 M 线性化很烦也可以用对偶 KKT 条件直接推导出等价的混合整数单层问题但那样推导量更大而且对初学者不那么友好。复现这个项目我用的是“对偶 大 M”路线已经足够解决 24 时段的小规模算例。3.4 收敛判断与实际运行时的心得收敛判据用相对间隙比绝对间隙更合理。成本量级在几千到几万元之间绝对间隙设成一个固定值在不同量级下意义完全不同。我用相对间隙 0.1% 作为停止条件效果很好。关于初始上界的取值有一个提速技巧。第一轮迭代前 UB 不知道可以先跑一个确定性模型Gamma 0把它的最优解当成一个可行方案代入子问题求一个初始 UB。这样第一轮迭代的区间就从“inf 到 -inf”变成了一个具体的数值区间收敛速度会明显加快。此外主问题和子问题每次迭代都要重新构建优化模型。如果在 MATLAB 循环里反复调用 sdpvar 定义全新变量会非常耗时。我的建议是把主问题模型作为基础模型只搭建一次在循环里只修改需要增加的约束和变量使用 YALMIP 的 assign 和 optimize 复用已有模型结构。4. MATLAB 实现从数据准备到结果输出4.1 算例数据怎么构造完全复现的工程第一步不是写代码而是准备好一套结构完整的算例数据。我自己用的是一套 24 时段、总负荷峰值约 1.5 MW 的区域微电网数据包含一台 800 kW 燃气机组、容量 500 kW / 1000 kWh 的储能、300 kW 光伏和 200 kW 风电联络线功率上限 800 kW。光伏预测曲线可以取一个类似“晴天”的钟形曲线早晨 6 点开始上升中午 12 点达到峰值19 点后归零。风电预测曲线取一条随机波动但整体平稳的曲线。负荷曲线按“两峰一谷”的工业与居民混合模式构造上午 9 点一个高峰晚上 20 点一个高峰凌晨低谷。这些曲线不要求来自真实数据但形状必须合理。为了让对比结果更清晰我把光伏峰值时刻设在与负荷高峰重叠的时间段附近这样不确定性对系统的影响会更明显。电价采用分时电价峰时 1.2 元/kWh平时 0.8 元/kWh谷时 0.4 元/kWh。这个价格差会让储能“低充高放”的价值充分体现出来。如果你的目标是把论文结果复现出数值量级一致的效果最好的办法是把算例参数做成独立的数据文件比如 Excel 表格统一管理。我曾经见过不少同学把 24 个时段的参数一条一条手写在脚本里一旦要调整数据改起来就是噩梦。独立数据文件 主脚本读取的方式后续跑参数敏感性分析会省很多时间。4.2 主循环代码框架下面给一个可以照着写的骨架代码。注意这里展示的是思路级别的框架不是可以直接复制就能跑的完整程序但你把这套骨架填进 YALMIP 语法后就能得到一版可用的复现程序% 初始化边界 LB -inf; UB inf; iter 1; Gap_list []; % 初始化主问题不带任何场景仅第一阶段变量与约束 [MP, x] create_master_problem(data); while (UB - LB) / abs(UB) tol iter MaxIter % 1. 求解主问题 optimize(MP.Constraints, MP.obj, options); x_star value(x); LB max(LB, value(MP.obj)); % 2. 固定 x_star求解子问题 [SP, u, y] create_subproblem(data, x_star); optimize(SP.Constraints, -SP.obj, options); u_star value(u); obj_SP -value(SP.obj); UB min(UB, c_phase1(x_star) obj_SP); % 3. 记录间隙 gap (UB - LB) / abs(UB); Gap_list(end1) gap; % 4. 将最坏场景加入主问题 MP add_scene_to_master(MP, u_star, data); iter iter 1; end这里有个细节值得说明子问题求的是 max-min而求解器默认求 min所以子问题的目标在代码里要加负号解完后再取负回来。第一次写这个逻辑的时候很容易搞混正负号建议每一步都用数值做一次合理性校验。4.3 储能约束的 YALMIP 写法储能 SOC 递推约束在 YALMIP 里可以直接写成向量形式速度比 for 循环快不少% 变量定义T 24 P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); SOC sdpvar(T1, 1); z_ch binvar(T, 1); % 充放电状态 % SOC 递推 SOC(2:end) SOC(1:end-1) eta_ch * P_ch * dt - P_dis * dt / eta_dis; % 充放电互斥 P_ch 0; P_ch P_ch_max * z_ch; P_dis 0; P_dis P_dis_max * (1 - z_ch); % SOC 边界与初末值 SOC_lb SOC SOC_ub; SOC(1) SOC_initial; SOC(end) SOC_initial; % 是否要求周期末回到初始值按工况决定周期性约束SOC 末端等于初始值要不要加我建议加。如果不加模型会倾向于在调度周期结束时把 SOC 放到最低把储能的剩余价值“吃干榨净”虽然对单次调度目标有利但不符合运行实际。加上周期性约束后系统要为下一个调度周期留好“余粮”这本身就是鲁棒性的一部分。4.4 结果输出与关键指标分析程序跑完后至少输出以下几张图和一组对比数据第一调度结果曲线图。把 24 时段的燃气机组出力、购电功率、储能充放电功率、光伏和风电预测出力画在同一张图里检查功率平衡关系是不是在每个时段都闭合。第二SOC 曲线图。SOC 必须在上下边界之间平滑变化如果出现跳跃多半是递推公式写错了。第三收敛曲线图。横轴是迭代次数纵轴是 LB 和 UB。正常形状是 LB 逐步上升、UB 逐步下降最后闭合到一个点。如果 LB 和 UB 出现“顶牛”不收敛优先检查主问题里是不是漏加了某个场景约束。第四确定性模型和鲁棒模型的结果对比表。重点对比三个指标总运行成本、最坏场景下的切负荷量、弃风弃光量。鲁棒模型的成本会比确定性模型略高高出来的部分就是“为不确定性付出的保险费用”也叫鲁棒代价。如果鲁棒代价只有百分之几而切负荷风险大幅下降这套方案的工程价值就体现出来了。5. 常见问题与排查技巧实录5.1 求解器报错“无可行解”从哪来复现过程中最常见的错误就是模型无解。我遇到过几种情况可以给到排查方向。储能约束自相矛盾是最隐蔽的一种。比如把 SOC 上下限设成 10%~90%但容量只有 1000 kWh某时段又要求 SOC 在一个小时内从 90% 充到 100%给定充放电功率上限根本做不到系统就无解。排查办法是先把储能约束单独提出来跑一个只含储能和功率平衡的极小模型看看还无解。用这种“拆零件”的方式定位问题比看一堆报错日志快得多。另一个典型的无解原因是功率平衡约束写成了严格的等号而不确定场景下系统又没有足够的调节余量。解决思路是在平衡约束里显式引入“失负荷”和“弃电”两个非负松弛变量配合大惩罚系数。这样代码几乎不会因数学上不可行而中断而且优化结果会告诉你“最坏场景下风险有多大”——这正是鲁棒分析的目标。5.2 求解结果里储能同时充放电如果发现储能某时段既充电又放电先看充放电互斥约束是否真的加入了主问题。我踩过的一个坑是主问题里加了互斥约束但加入新场景时忘了对场景对应的第二阶段变量复制一套充放电状态 0-1 变量导致新场景下的储能变量没有互斥限制优化器就“钻空子”同时充放电来满足平衡约束。检查方法很简单把 P_ch 和 P_dis 的曲线叠加在一张图里如果同一时段两根柱子都不为 0那就是约束漏了。修正后任何时段最多只有一个方向有值。5.3 收敛慢问题出在初值和 M 值收敛慢的绝大部分原因是 M 值取得不合适。M 过小辅助变量的取值域被错误截断子问题找不到真正的最坏场景主问题就会反复产生相同的场景迭代一直在原地打转。M 过大数值病态导致求解器精度下降间隙到不了收敛阈值。一个实用经验是M 值不要拍脑袋统一设成 1e6。按每个约束的实际物理边界来定比如某条约束右侧与功率上限相关M 就取功率上限的 10 倍。这样既不会截断也不会引入严重数值问题。下表是我复现过程中整理的几个高频问题和对应处理建议现象常见原因处理建议子问题返回的目标明显偏小对偶变量符号正负搞反用一个小算例手算验证对偶推导循环迭代超过 20 轮仍不收敛主问题新场景加入逻辑错误检查新场景的约束是否确实与 x 形成耦合Gamma 调大但成本不变不确定集合名义偏差过小检查 delta 设置确认预测值 ± 偏差是否覆盖实际波动范围SOC 曲线出现锯齿充放电效率写法不一致统一使用 eta_ch 和 eta_dis 两个独立效率结果与确定性模型完全一样子问题 max 部分总取到零场景检查不确定变量是否真正出现在子问题约束中5.4 鲁棒度调节的经验最后聊一下保守度调节。Gamma 从 0 到最大值的全过程我建议每个值各跑一遍把总成本画成一条“Gamma-成本”曲线。这条曲线的形状一般先平缓上升到某个值后开始陡增。那个“膝点”就是工程上最值得采用的保守度。不要盲目追求“完全对抗所有可能波动”。真实系统里风光出力同时出现最大正负偏差的概率极低Gamma 设得过大只会让日前计划变得非常僵硬储能永远留着一半电量不敢用可再生消纳率反而下降。我在实际调试中Gamma 取 4 到 6 之间时成本增加约 3%-5%但切负荷量降到了接近 0这个性价比是比较舒服的。复现这类模型时间久了我有一个很深的体会代码本身只是最后一步真正费时间的是把子问题对偶推导、不确定集设计、储能约束耦合这三件事想明白。MATLAB 里的每一行代码本质都是在把模型公式“翻译”给求解器听。翻译错了求解器再强也只会给你一个精美的错误结果。如果你刚接触这个项目我的建议是先手推一遍两阶段模型把主问题和子问题之间的变量传递关系画在纸上再去写代码。等你把这个框架跑顺了往里面加需求响应、电动汽车、多微网互联都只是在现有骨架上添砖加瓦的事。
阅读完成 · 觉得有帮助？
