调度员最怕的不是负荷突增而是头天晚上排好的计划第二天早上九点被光伏出力的暴涨彻底打乱。含分布式电源的配电网日前优化调度之所以难做难点不在优化算法本身而在于你做的每一个决定离实际执行的时刻还有24小时——这段时间里光照会变、风会停、负荷会偏离预测。这篇文章想聊的是怎么用Matlab实现一个两阶段框架下的日前调度模型第一阶段定“大方向”第二阶段算“细账”。适合正在做配电网调度、微电网优化或者毕业设计需要搭模型的研究生和工程师也适合那些已经跑通单阶段优化想进一步处理不确定性的朋友。我的习惯是不空谈模型直接上可复现的东西。下面会从数学建模、分布式电源建模、CCG求解、到MatlabYalmip代码架构完整过一遍包括我自己踩过的坑。1. 为什么日前调度要拆成两个阶段——先定大方向再算细账1.1 单阶段日前调度模型到底缺了什么单阶段日前调度模型是把24小时看作一个整体按照预测的负荷、光伏、风电出力一次性算出未来整天的机组出力、储能充放电计划和购电曲线。听起来没问题但工程上一算就知道不对劲。光伏出力和实际负荷不可能和预测完全一致。假设预测明天中午光伏出力是300kW结果实际来了500kW储能如果按预测结果安排好了放电计划配电网就会出现功率倒送、电压越限。反过来如果预测风电有200kW实际只有50kW那晚高峰的功率缺口就得临时从上级电网买高价电或者直接切负荷。单阶段模型把“预测误差”当作不存在来处理本质上是在赌。更麻烦的是有些决策一旦定了就不好改。储能是充电还是放电状态切换涉及电池寿命损耗不能让它每15分钟翻一次微燃机启动一次要烧不少燃料、还有最小运行时间限制不可能说停就停向上级电网的购电计划在日前市场里已经签了合约日内临时多买电要掏惩罚性电价。这些决策有一个共同特点执行周期长、切换代价高、受预测误差影响大。它们必须在日前阶段就定好不能等实际出力出来后再反应。而另一类决策则相反——微燃机的出力微调、储能功率小幅修正、部分可中断负荷的启停这些调整时间尺度短、代价相对可控可以等不确定性稍微明朗后再做。把这两类决策放在同一个单阶段模型里要么保守得离谱要么脆得不行。1.2 哪些决策必须留到第二阶段两阶段模型的分界线就是按“决策的柔韧性”来划的。第一阶段日前阶段处理的是储能的充放电状态与基准充放电功率、微燃机的启停状态、从上级电网购买的基准功率曲线。这些决策的特点是一旦确定后续很难低成本更改。第二阶段日内调整阶段处理的是分布式电源实际出力出来后对微燃机出力的调整量、储能功率的微调量、切负荷量、弃光弃风量。打个比方这就像公司做年度预算年初定下来今年要投哪些产线、招多少人第一阶段这个盘子不能随便动等季度财报出来了再根据实际经营情况调整供应商订单量、关停部分临时项目第二阶段。你不会等到年底再后悔年初投资决策错了但也不可能每个星期都把产线拆了重装。从数学上看第一阶段决策变量进了第二阶段的约束条件第二阶段能调多少微燃机出力取决于第一阶段微燃机是否启动储能第二阶段能充多少电取决于第一阶段定的充电状态和当前SOC。这就是两个阶段耦合的地方也是后面CCG算法要反复在主问题和子问题之间迭代的根本原因。如果把功率平衡写出来第一阶段是系统级的预测平衡第二阶段则要校核最坏情况下的平衡并给出调整方案。两阶段模型本质上是在问一个更现实的问题如果明天实际情况是“最坏的那一种”我这套日前计划还扛得住吗扛不住的话日内调整要花多少钱2. 两阶段数学模型的变量分层与成本函数设计2.1 目标函数第一阶段成本与第二阶段成本的拆解目标函数分两层这是整个模型的骨架。第一阶段的目标是确定性成本表达式概念如下min Σ_t [ C_ugrid(t) * P_grid(t) f_cost(MT(t)) SU_cost * y_start(t) C_battery * |P_b(t)| * Δt ]每一项的含义向配电网购电成本、微燃机燃料成本一般是非线性二次函数要按分段线性化处理、启动成本、储能充放电的寿命折算成本。购电成本用的是日前分时电价这对应前面说的日前市场交易。第二阶段的目标是在不确定性实现后的最坏情况调整成本Q(x) max_{u∈U} min_{y∈Ω(x,u)} Σ_t [ C_adj * y_adj(t) C_load_cut * P_cut(t) C_curtail * (P_pv_curt(t) P_wt_curt(t)) ]注意这里的结构外层的max是在找最坏的不确定性场景内层的min是在该场景下找最经济的调整策略。这就是两阶段鲁棒优化的标准形式。第二阶段成本里各惩罚项的系数设置直接影响结果合理性。切负荷惩罚一般设为峰值电价的5到10倍比如峰值电价1.2元/kWh切负荷惩罚取6到12元/kWh否则模型为了省钱会动不动就切负荷。弃光弃风惩罚反而要低取0.05到0.1元/kWh低于正常发电成本让模型只在真正过剩时才弃。惩罚系数的优先级排序是切负荷 购电调整 弃光弃风。这个次序在调试时特别重要如果设置的系数倒挂模型行为会非常反直觉。2.2 第一阶段决策变量与约束体系第一阶段决策变量清单如下变量类型说明P_grid(t)连续24小时向上级电网的购电功率z_ch(t)0-1储能充电状态1为充电P_b_ch(t), P_b_dis(t)连续储能充/放电功率SOC(t)连续储能荷电状态u_mt(t)0-1微燃机启停状态P_mt(t)连续微燃机基准出力y_start(t)0-1微燃机启动标志约束分几类。系统功率平衡是基础每个时段所有电源出力加购电功率等于负荷预测值加储能充电功率P_grid(t) P_pv_pred(t) P_wt_pred(t) P_mt(t) P_b_dis(t) P_load_pred(t) P_b_ch(t)储能自身的动态约束SOC递推方程、同一时刻不能同时充放电、充放电功率不能超过额定限值、SOC保持在上下限内。这里有个容易漏掉的约束——日末SOC回弹约束。如果模型只优化一天且不约束调度周期末的SOC那储能会在最后几个时段把电量全部放光等于偷了第二天的能量。一般加上SOC(24) SOC(0) 或者 SOC(24) SOC(0)。微燃机约束包括爬坡约束向上/向下爬坡率、最小运行时间和最小停机时间。启停状态的逻辑关系用大M法写y_start(t) 1 表示第t时段启动了机组需要满足 u_mt(t) - u_mt(t-1) y_start(t)。上级电网购电约束也很实际联络线容量有限P_grid(t) 不能超过变压器额定容量否则配电网节点电压会出问题。2.3 第二阶段决策变量与约束体系第二阶段决策变量是第一阶段计划在实际场景下的“修正量”微燃机出力调整量 ΔP_mt(t)可正可负受爬坡约束和出力上下限约束储能功率微调量 ΔP_b(t)同样受限幅约束但注意不能改变第一阶段定的充放电状态切负荷量 P_cut(t)上限是该时段的负荷值弃光弃风量 P_curt(t)上限是对应DG的实际出力第二阶段约束里最核心的是任意不确定性场景下的功率平衡P_grid(t) P_pv_act(t) P_wt_act(t) P_mt(t) ΔP_mt(t) P_b_dis(t) ΔP_b(t) P_load_act(t) - P_cut(t) P_b_ch(t) ΔP_b(t)注意这个式子里的左右两边都要处理ΔP_b(t) 在充电侧还是放电侧取决于第一阶段确定的储能状态。为了严谨一般直接把储能净功率 N_b(t) P_b_dis(t) - P_b_ch(t) ΔP_b(t)然后功率平衡写成一个等式。第二阶段的约束里有不确定性参数 u光伏、风电的实际出力、实际负荷都在不确定集合U中变动。对每个给定的u内层min问题是一个线性规划如果第一阶段变量固定这决定了后面可以用对偶理论处理。第二阶段的模型边界要说明清楚很多配电网两阶段模型会简化系统级功率平衡而不计潮流或者用DistFlow线性化处理电压约束。前者适合辐射状网络结构比较简单、电压问题不突出的场景后者需要额外引入节点注入功率和电压变量。如果做IEEE 33节点这类算例建议至少用DistFlow线性化否则调度结果在真实网络里根本转不起来。3. 分布式电源建模光伏、风电、储能、微燃机的核心约束3.1 光伏与风电的出力不确定模型光伏和风电是配电网不确定性的主要来源建模要把物理特性和数学表达结合好。物理上光伏出力可以写成P_PV(t) η * A * S(t) * (1 - 0.005 * (T_amb(t) - 25))η是光电转换效率A是光伏板面积S是辐照度T_amb是环境温度。这个式子的含义是辐照度越强出力越高但温度过高反而降低转换效率。工程上如果手头没有这么细的物理参数更常见的做法是用预测曲线乘不确定系数P_pv_act(t) P_pv_pred(t) * (1 u_pv(t)), u_pv(t) ∈ [-l_pv, u_pv]u_pv(t) 就是光伏出力偏差率上下限表示“光伏预测可能偏大或偏小多少”。这个形式做不确定集合非常方便后面CCG的子问题就是靠它构造出来的。风电的物理模型是用风速到功率的分段函数P_WT(V) 0, V V_cutin 或 V V_cutout P_WT(V) P_rated * (V - V_cutin) / (V_rated - V_cutin), V_cutin V V_rated P_WT(V) P_rated, V_rated V V_cutoutV_cutin是切入风速V_rated是额定风速V_cutout是切出风速。实际调度模型同样可以简化为预测值乘不确定系数。风机还有一个特点是响应速度很快逆变器控制下基本不设爬坡约束但如果有无功调节需求需要额外考虑无功能力边界本文先不展开。3.2 储能与微燃机连接两个阶段的调度枢纽储能是整个两阶段模型的“调度枢纽”因为它是唯一能在时间维度上转移能量的设备。SOC递推约束是储能建模的核心SOC(t1) SOC(t) (η_ch * P_b_ch(t) - P_b_dis(t) / η_dis) * Δt / Cap这里η_ch是充电效率η_dis是放电效率Cap是储能容量Δt是时段长度一般是1小时。充放电效率不对称会带来能量损耗体现在SOC变化量上。同一时刻不能同时充放电的约束用0-1变量表达% Yalmip写法 C [C, 0 P_ch(t) P_ch_max * z_ch(t)]; C [C, 0 P_dis(t) P_dis_max * (1 - z_ch(t))];z_ch(t) 1表示充电此时P_dis(t)被压到0反之亦然。储能寿命折算成本一般按充放电吞吐量计算也可以用一个很小的单位折旧成本乘功率绝对值。微燃机的成本函数是二次的f_cost(P_mt) a * P_mt^2 b * P_mt cGurobi、Cplex这些商业求解器默认不支持二次成本直接进MILP虽然支持二次目标但配合大量整数变量求解效率会大打折扣所以工程上要分段线性化。做法是取几个典型出力点把二次曲线近似成折线。爬坡约束写成-P_ramp_down P_mt(t) - P_mt(t-1) P_ramp_up最小启停时间约束稍微复杂但配电网算例里如果微燃机参数量小可以先跳过等模型跑通再加。4. 两阶段鲁棒优化求解主问题-子问题与CCG迭代4.1 max-min子问题的对偶变换两阶段鲁棒优化的模型写出来以后不能直接丢给求解器因为第二阶段是max-min嵌套结构商业求解器不认识。CCG算法Column-and-Constraint Generation列与约束生成是处理这个问题最主流的方法。CCG的核心思想是把问题拆成一个主问题MP和一个子问题SP然后迭代。主问题包含第一阶段变量、辅助变量η以及每次迭代时从子问题抓出来的“最坏场景约束”。初始时主问题里没有这些约束所以算出来的目标值偏乐观这就是下界LB。子问题是给定第一阶段决策x*后去求解最坏场景下的最小调整成本这个成本加上第一阶段成本构成上界UB。两个界不断逼近直到小于收敛阈值。子问题的核心难点是max-min的可解性。看这个形式Q(x*) max_{u∈U} min_{y∈Ω(x*,u)} c^T y内层min是一个线性规划LP对偶问题也是LP。根据强对偶定理如果内层原问题可行且目标值有界那么内层min的最优值等于其对偶max的最优值。这样把内层min替换成对偶问题的max整个子问题就变成Q(x*) max_{u∈U, λ∈Λ} f(x*, u, λ)两个max合并成一个max原问题变成单层优化。这里λ是对偶变量它对应的约束集合Λ是从内层原问题翻过来的。4.2 双线性项的线性化处理对偶变换后你会遇到一个麻烦目标里出现u和λ的乘积项比如u_pv(t) * λ_t这是一个双线性项而λ_t是有界变量、u_pv(t)在区间和预算约束下取值。双线性项让问题变成非凸的Gurobi处理不了除非用bilinear特殊模式但工程上不推荐。解决双线性项最常用的是大M线性化。思路是不确定变量u_pv(t)虽然有区间但最优解一定落在不确定集合的极点附近所以我们把u_pv(t)拆成两个0-1标志和偏差量的组合。或者更直接一点引入辅助变量w_t u_pv(t) * λ_t然后用一组带大M的约束线性化。-M * z_low(t) w_t - λ_t * u_low M * z_low(t) -M * z_up(t) w_t - λ_t * u_up M * z_up(t) z_low(t) z_up(t) 1这样做的含义是每个时段要么取u的下界要么取u的上界对应的w_t近等于λ_t乘对应界值。M取值不能太小要大于λ_t的物理边界也不能太大太大会导致数值病态。我一般取λ_t理论最大值的2到5倍。如果你发现子问题在大M线性化后求解时间无法接受还有一条路不确定集合的极点是“除Gamma个时段取上界外其余取下界”的组合理论上可以用多面体极点枚举加上约束生成来做但实际配电网24时段Gamma接近12的时候组合数量是天文数字工程上还是大M线性化最稳。4.3 CCG主循环流程与收敛判据整个CCG迭代流程如下初始化LB -infUB inf迭代计数k 1构造不含cut的主问题MP。求解MP得到第一阶段最优决策x和辅助变量η令LB MP目标值。固定x*求解子问题SP得到最坏场景u和最优调整成本Q。计算UB min(UB, 第一阶段成本(x*) Q*)。如果UB - LB ε比如0.01停止迭代输出方案。根据u*构造新的第二阶段变量y_new和对应约束加入MPk k1回到步骤2。注意CCG加cut的方式和Benders分解不一样。Benders加的是一张对偶信息组成的割平面CCG加的是“完整的最坏场景变量和约束”也就是把子问题在最坏场景u*下的原问题原封不动搬进主问题。这也是CCG比Benders收敛快的原因——每次迭代都引入了实际场景变量可行域的逼近更激进。关于收敛判据工程上gap阈值取0.011%通常够用严谨一点取0.001。调试的时候可以每一轮打印LB和UB正常情况下LB从低往高爬、UB从高往低降两条线慢慢合拢。如果LB和UB出现剧烈波动不收敛基本可以确定是子问题对偶变换或双线性项线性化出了问题。5. Matlab代码架构与实操要点从数据到结果5.1 求解器选型与Yalmip配置两阶段鲁棒优化最终落到代码端我的建议是用Matlab Yalmip 商业求解器Gurobi优先、Cplex次之。理由很简单Yalmip的建模语法能大幅减少写约束的时间和出错概率而且对第二阶段对偶变换后的双线性项线性化支持很友好。如果你的机器上没有Gurobi也可以先用HIGHS这个开源求解器Yalmip直接支持。但实测下来配电网两阶段MILP规模上来后HIGHS的求解速度会明显落后于Gurobi特别是有三四百个0-1变量的时候。学术界Gurobi可以提供免费license注册一个学术账号就能用。Yalmip配置求解器的核心代码ops sdpsettings(solver, gurobi, verbose, 2, gurobi.MIPGap, 0.001, gurobi.TimeLimit, 600);这里MIPGap设到0.001TimeLimit设600秒防止某个子问题卡死。我在实际调试里发现verbose开2能盯着迭代过程开0则啥都看不到、出了问题不好定位。5.2 输入数据准备与代码文件结构输入数据建议整理成Excel或CSV统一读入不要硬编码在脚本里。最少需要以下数据数据名称维度说明负荷预测曲线24x1单位kW光伏预测出力曲线24x1单位kW风电预测出力曲线24x1单位kW分时电价曲线24x1单位元/kWh微燃机参数1x1结构体容量、爬坡率、成本系数储能参数1x1结构体容量、功率上限、效率、SOC范围不确定集合参数结构体偏差上下界、预算Gamma代码文件结构建议分五个文件main_ccg_two_stage.m —— 主程序初始化数据、调CCG求解、输出结果data_loader.m —— 从Excel/CSV读数据build_mp.m —— 构建主问题Yalmipbuild_sp.m —— 构建子问题Yalmipplot_results.m —— 绘制功率平衡图、SOC曲线、迭代收敛图主问题构建的核心片段长这样% 第一阶段决策变量 z_ch binvar(24, 1); % 储能充电状态 P_ch sdpvar(24, 1); % 充电功率 P_dis sdpvar(24, 1); % 放电功率 SOC sdpvar(25, 1); % 荷电状态从0到24 P_mt sdpvar(24, 1); % 微燃机出力 u_mt binvar(24, 1); % 微燃机状态 eta sdpvar(1, 1); % 第二阶段成本代理变量 % 约束集合 MP_cons []; % 储能约束 MP_cons [MP_cons, 0 P_ch P_ch_max * z_ch]; MP_cons [MP_cons, 0 P_dis P_dis_max * (1 - z_ch)]; MP_cons [MP_cons, SOC(1) SOC_init]; MP_cons [MP_cons, SOC(2:end) SOC(1:end-1) (eta_ch * P_ch - P_dis / eta_dis) / Cap]; MP_cons [MP_cons, SOC_min SOC SOC_max]; % 功率平衡预测场景 MP_cons [MP_cons, P_grid P_pv_pred P_wt_pred P_mt P_dis P_load_pred P_ch]; % 第二阶段成本代理约束迭代时动态添加cut子问题构建要特别注意给第一阶段传入的x*赋值Yalmip可以用assign()固定变量值assign(P_mt, value(P_mt)); assign(P_ch, value(P_ch)); assign(P_dis, value(P_dis)); assign(z_ch, value(z_ch));然后重新定义第二阶段的调整变量构建对偶问题。代码里有个常见的坑如果只用assign而不把第一阶段变量从问题里去掉子问题的优化变量会包含第一阶段变量导致求解结果是“连第一阶段也可以改”的结果得到的Q*会偏小UB会失真。正确做法是在子问题里把第一阶段变量当作参数、重新定义独立的sdpvar。5.3 调试经验与常见错误处理两阶段模型调试比单阶段容易翻车我把实际踩过的坑列出来。第一先跑确定性版本。把不确定集合的预算Gamma设为0偏差上下界都设为0整个模型就退化成单阶段。先确认这个退化版本的结果合理再逐步加入不确定性。我见过不少同学直接上完整两阶段模型结果子问题一直不可行最后发现是第一阶段功率平衡本身写错了这种问题在退化版本里一眼就能看出来。第二子问题不可行的处理。第一阶段决策比较激进时比如把微燃机停了、储能放空了极端场景下可能无法满足功率平衡。解决办法是给子问题加松弛变量比如允许切负荷的上限放大或者在功率平衡等式里加一个人工松弛项并给极高惩罚。这样即使第一阶段决策很极端子问题也能给出一个有限的目标值CCG迭代能继续跑下去。第三大M值不要拍脑袋。双线性项线性化的M取值过大会让子问题求解出现数值问题Gurobi报告“numerical trouble”或返回错误的dual。我的经验是先跑一次不加线性化的近似版本观察λ_t的大致范围再设定M 2 * max(λ_t)。如果模型规模变化M也要对应调整。第四绘图验证结果合理性。收敛后至少画三张图24小时系统功率平衡堆叠图购电、各DG、储能出力、负荷、储能SOC曲线、CCG迭代的LB/UB收敛曲线。SOC曲线如果出现剧烈锯齿形振荡多半是储能状态切换频率约束漏了收敛曲线如果UB在中间某轮突然跳高往往是子问题里出现了违反不确定性约束的场景。我自己刚开始写这个模型的时候在子问题对偶变换和双线性项的线性化上卡了将近一周后来换了个思路——先完全不处理双线性项用枚举几个极端场景的方式验证对偶变换的方向是否正确确认对偶本身没问题后再回来怼大M线性化。把问题拆小、每一步先用肉眼可验证的简单算例校准比一头扎进完整模型里高效得多。这套模型跑通之后后续加需求响应、加节点电压约束、改成交替迭代的分布式求解都是在这个骨架上做增量开发性价比很高。