含分布式电源的配电网两阶段日前优化调度Matlab建模详解

前段时间我帮某课题组把一套含分布式电源的配电网日前两阶段优化调度模型完整落地成了Matlab代码整个过程比想象中曲折也踩了不少意料之外的坑。这个项目的核心一句话就能讲清楚在已知光伏、负荷等预测曲线的前提下把未来24小时的调度决策拆成“第一阶段定计划、第二阶段纠偏差”两步最终得到一套结构清晰、可以直接跑出结果的Matlab程序。它能处理分布式电源出力波动、储能充放电安排、与主网的交换功率以及最让人头疼的预测误差问题。这套模型适合正在做配电网优化、分布式能源接入、储能调度相关课题的研究生也适合有线性规划和Matlab基础、想快速上手两阶段调度代码的工程师。代码的价值不只是“能跑”而是把不确定性建模、场景生成、Yalmip建模、商用求解器调用这一整条链路打通了。下面我按自己的实际调试顺序把这套模型从数学到代码完整拆开讲一遍。1. 项目到底是什么解决什么问题1.1 为什么配电网调度需要“两阶段”配电网的日前调度本质上是“明天还不知道准确情况就得先拍板”的问题。光伏出力和负荷预测永远是带误差的如果把预测数据当成真实值只做一个单阶段优化得到的调度计划往往在真实场景下并不可行可能出现功率越限、电压越限、储能SOC对不上甚至需要频繁调整机组出力。两阶段调度就是把“计划”和“调整”分开。第一阶段属于“现在决策”依赖预测信息制定未来24小时的基础计划例如每台可控DG出多少电、储能什么时候充什么时候放、跟主网买多少电。第二阶段属于“看到场景后决策”在真实或虚拟的不确定场景下以最小调整成本去修正第一阶段计划从而保证各场景下系统都能安全运行。打个比方这很像出门前看天气预报。决定带不带伞是第一阶段的“日前计划”真下雨了选择打伞、躲雨还是打车是第二阶段的“实时调整”。模型里不会真的等到明天才调整而是用大量可能的风光、负荷场景把“明天可能发生的各种情况”提前模拟一遍。这样算出来的日前计划才是一个在所有可能场景下都不会离谱的结果。1.2 这套代码能做什么、适合谁从功能上看这套代码输入是各节点负荷预测曲线、光伏逆变器可发功率预测、发电机组参数、储能参数和网络拓扑输出则是一整套可落地的调度方案包括各时段各节点可控分布式电源的有功出力、储能充放电功率与荷电状态、主网购电功率、弃光量和切负荷量。从研究角度来说它把两阶段随机优化这个抽象概念做成了看得见摸得着的程序。你可以在同一套框架下修改预测误差标准差、场景数量、惩罚因子、网络参数观察目标成本和调度曲线怎么变。也可以把第二阶段从“随机场景期望成本”扩成“最坏场景下的鲁棒优化”这在思路上是完全兼容的。适合参考的人群主要有三类一是在做含高比例分布式电源配电网调度课题的学生二是想用Yalmip建模但不知道怎么把两阶段模型落到代码里的开发者三是已经有一些优化调度基础想快速搭建一套基准模型做对比验证的工程师。需要的前提是对Matlab基本编程有概念知道变量、约束、目标函数在优化模型里的关系。2. 两阶段日前优化调度模型的数学拆解2.1 目标函数第一阶段成本与第二阶段期望成本两阶段随机规划的目标函数最容易理解的形式是这样min F Phase1_cost Σ_s p_s × Phase2_cost_s其中Phase1_cost是基础计划对应的确定性成本Phase2_cost_s是第S个场景发生后的调整成本p_s是场景概率。把期望调整成本加进去之后第一阶段决策就不会为了降低眼前成本而把风险全甩给未来。第一阶段的成本通常包括主网购电费用可控分布式电源的发电成本可用分段线性函数近似可控机组的启停成本储能充放电折算的损耗成本第二阶段的成本通常带有惩罚性质对第一阶段计划的正负调整量成本弃风弃光惩罚切负荷惩罚这里有个关键点第二阶段成本里的惩罚系数一定要比正常购电、发电费用大上一截否则求解器会“很聪明”地用切负荷来逃避调整成本。我常用的做法是把切负荷惩罚设成正常电价的5到10倍弃光惩罚也类似。惩罚系数不是越大越好过大会让模型过度保守储能提前充满、DG不敢优化最终整个调度方案失去经济性。2.2 决策变量与关键约束设计两阶段模型里变量分成两类。第一阶段变量是公共的对所有场景都一样例如P_dg(g,t)第g台可控DG在t时段出力z_dg(g,t)启停状态0-1整数变量P_buy(t)主网购电功率P_ch(t)、P_dis(t)储能充电、放电功率SOC(t)储能荷电状态第二阶段变量按场景独立例如ΔP_dg(g,t,s)第S个场景下DG出力调整量ΔP_buy(t,s)联络线交换功率调整量ΔP_ch(t,s)、ΔP_dis(t,s)储能功率修正量P_curt(t,s)弃光功率P_shed(t,s)切负荷功率约束条件按功能分五组约束类型内容作用功率平衡约束每个节点/每个时段注入与流出功率平衡保证物理可行网络潮流约束支路有功、无功与节点电压关系保证辐射状配电网的运行安全机组约束DG出力上下限、爬坡率、启停逻辑保证发电设备可执行储能约束SOC递推、充放电功率限值、容量上下限保证储能连续运行第二阶段调整边界ΔP必须在第一阶段计划允许的范围内体现“纠偏”而非“推翻”储能SOC递推公式是SOC(t) SOC(t-1) η_ch*P_ch(t)*Δt - P_dis(t)*Δt/η_dis。很多第一次写两阶段代码的人会在第二阶段把储能SOC也单独设一套变量结果每个场景的SOC曲线互不相关导致计划看起来每个场景都对但实际储能容量根本对不上。我的建议是第一阶段SOC作为公共变量第二阶段的储能功率修正量也需要同步更新SOC且SOC仍然是一个跨所有场景公共的时序变量只是功率允许按场景小幅调整。这样物理含义才正确。2.3 不确定场景怎么建模两阶段模型最核心的部分不是数学公式而是不确定性场景怎么生成。我一般这样做先假设光伏出力和负荷预测误差服从正态分布然后用蒙特卡洛抽样生成大量原始场景再用同步回代削减法挑出几个代表性场景。场景生成公式很简单。预测值是P_hat误差比例是ε场景值就是P_scene P_hat × (1 ε)其中ε服从均值为0、标准差为δ的正态分布。δ不建议拍脑袋最好用历史预测与实际出力的统计数据算出来。没有历史数据时光伏预测误差标准差通常取10%到20%负荷预测误差取2%到5%。太小会让两阶段退化成单阶段太大又会过度保守储能几乎提前放空。场景数量方面我的经验是先抽200个原始场景再用同步回代削减到10到20个。10个场景已经能够覆盖绝大多数分布特征20个场景结果更精细但求解时间会明显增加。目标函数里用到的p_s就是削减后每个场景保留的概率权重所有场景权重之和应为1。这一步如果做错两阶段的“期望成本”就会失真算出来的计划会偏向大概率场景忽略小概率高风险场景。3. 基于Matlab的完整代码实现流程3.1 开发环境与工具箱准备我在复现这套代码时用的环境是Matlab R2020a以上版本配上Yalmip工具箱底层求解器用的是一套主流的商用MILP/SOCP求解器。Yalmip是建模语言它本身不求解任何问题只负责把你写的约束和目标转成求解器能识别的标准形式。安装好之后务必在Matlab命令行执行一次yalmiptest它会列出当前机器上Yalmip检测到了哪些求解器。如果显示Gurobi或CPLEX没有被识别多半是路径没加进Matlab。可以用addpath(genpath(求解器安装目录))加进去然后保存路径。这一步经常被忽略但90%的“No suitable solver”报错都出在这。如果不想用商业求解器Yalmip也支持开源求解器例如CBC、SCIP等。不过两阶段问题里有MILP凸约束建议还是优先用商业求解器求解速度和数值稳定性差距非常大。尤其是第二阶段场景多的时候开源求解器可能跑到几十分钟都出不来结果。3.2 工程文件结构与主程序框架代码不是把所有内容堆在一个脚本里就完事。我建议按下面结构组织main.m 主程序参数初始化、建模、求解、出图 load_case_data.m 读取负荷、光伏、网络和机组参数 scenario_generate.m 生成和削减不确定场景 build_two_stage_model.m 构建两阶段优化模型 plot_result.m 绘制调度结果曲线主程序框架大致是这样%% 初始化 clear; clc; close all; load_case_data; %% 场景生成 [load_scene, pv_scene, p_s] scenario_generate(load_hat, pv_hat, ns); %% 建模与求解 [model, solution] build_two_stage_model(...); plot_result(solution);这样拆的好处是后面想换场景生成方法、魔鬼网网约束、改目标函数都只需要动对应的函数不需要重写整个主程序。两阶段优化尤其依赖这种模块化结构因为变量多、约束多所有逻辑堆在一个脚本里很难定位问题。3.3 核心代码片段场景生成、Yalmip建模与求解先看场景生成函数。核心逻辑就是对每个时段、每个场景生成一个乘性误差function [load_scene, pv_scene, p_s] scenario_generate(load_hat, pv_hat, ns) T length(load_hat); raw_n 200; load_raw zeros(T, raw_n); pv_raw zeros(T, raw_n); sigma_load 0.03; sigma_pv 0.15; for s 1:raw_n load_raw(:,s) load_hat .* (1 sigma_load * randn(T,1)); pv_raw(:,s) pv_hat .* (1 sigma_pv * randn(T,1)); end % 这里可调用同步回代削减函数将200个场景削减为ns个 % load_scene: T×ns, pv_scene: T×ns, p_s: ns×1 end同步回代削减是一个经典算法每次删掉一个场景使剩余场景与删除场景的场景距离最小然后把被删场景的概率叠加到离它最近的场景上。这个算法用Matlab实现并不复杂也可以用kmeans聚类近似替代。聚类更快但场景概率不好直接得到回代削减物理含义更准确我推荐优先用回代。Yalmip建模部分变量声明是关键。第一阶段变量这样写P_buy sdpvar(24, 1, full); P_dg sdpvar(24, ng, full); P_ch sdpvar(24, 1, full); P_dis sdpvar(24, 1, full); SOC sdpvar(24, 1, full); z_dg binvar(24, ng, full);第二阶段变量要用到场景维度。如果第二阶段只有聚合的联络线调整量、储能调整量和切负荷量可以这样声明dP_buy sdpvar(24, ns, full); dP_ch sdpvar(24, ns, full); dP_dis sdpvar(24, ns, full); P_curt sdpvar(24, ns, full); P_shed sdpvar(24, ns, full);注意sdpvar(24, ns)生成矩阵时行是时段列是场景索引P_curt(t,s)是有效的。千万注意维度方向我见过很多人把24和ns写反后面加约束时出现维度不匹配却找不到原因。约束的写法我建议按模块加。先写第一阶段功率平衡和机组约束再写第二阶段每个场景的功率平衡。第二阶段场景平衡式可以写成constraints []; for t 1:24 % 第一阶段平衡 constraints [constraints, ... P_buy(t) sum(P_dg(t,:)) P_dis(t) - P_ch(t) ... load_hat(t)]; end for s 1:ns for t 1:24 constraints [constraints, ... P_buy(t) dP_buy(t,s) ... sum(P_dg(t,:)) sum(dP_dg(t,:,s)) ... P_dis(t) dP_dis(t,s) - P_ch(t) - dP_ch(t,s) ... load_scene(t,s) - P_shed(t,s)]; end end目标函数按之前公式实现objective sum(C_buy * P_buy) sum(sum(C_dg .* P_dg)); % 加启停成本 objective objective sum(sum(C_start .* max(0, diff(z_dg)))); % 加第二阶段期望惩罚 for s 1:ns objective objective p_s(s) * ( ... C_adj_buy * sum(abs(dP_buy(:,s))) ... C_curt * sum(P_curt(:,s)) ... C_shed * sum(P_shed(:,s)) ); end如果调整量带绝对值会引入非线性。Yalmip对绝对值有自动处理它会把它线性化但更好的做法是引入非负辅助变量把dP_buy拆成正调整和负调整。这样目标函数和约束都保持线性求解更稳。求解代码非常简单options sdpsettings(solver, gurobi, verbose, 2); sol optimize(constraints, objective, options); if sol.problem 0 P_buy_opt value(P_buy); else disp(求解失败); endsol.problem 0表示求解成功。这个判断一定要写不然遇到不可行问题时value(P_buy)返回的可能是NaN后面画图全乱。3.4 典型运行结果与调参经验一个调好的模型画出来的24小时曲线通常有三个特征主网购电功率会避开负荷高峰时段储能往往在谷时充电、峰时放电光伏大发时段可控DG出力会明显下调。这些现象说明模型把经济性和消纳放在一起权衡了。我在调试时最常用的手段是“变量归一化验证”。先只启用第一阶段成本把第二阶段惩罚系数全部设成0看结果是不是一个普通单阶段模型。如果单阶段结果正确再逐步加入第二阶段惩罚成本。这样能快速定位是第一阶段建错了还是第二阶段耦合出了问题。另一个经验是第二阶段调整量的边界一定要给。比如让dP_buy不超过联络线容量与第一阶段计划剩余容量的差让储能调整功率不超过充放电功率上限。如果不给边界第二阶段可以无限调整成本却只按惩罚系数线性增加求解器会宁可大调特调也不优化第一阶段两阶段就失去意义。4. 常见报错、调试方法与排查速查表4.1 求解器调用问题最常见的报错是No suitable solver for MILP。出现这个错误首先执行yalmiptest看一下求解器列表里有没有支持MILP的求解器。如果有但Yalmip没识别多半是安装路径不对如果压根没安装商业求解器建议直接换用支持MILP的开源求解器或者安装商业求解器并注册好许可证。Solver not found这类报错还有一个隐蔽原因用sdpsettings指定了求解器名称但名称拼写和Yalmip内部名称不一致。比如写了solver,gurbi少了一个字母Matlab不会报变量错误只会提示找不到求解器。遇到这类问题先删掉solver选项让Yalmip自动选择。4.2 模型表达与约束错误约束维度不匹配是新手最常踩的坑。sdpvar(24, ns)和P_load(:,s)相加没问题但如果你定义的是sdpvar(1, 24)再用(t)索引就会越界。我的调试习惯是每加一个变量就size(value(var))看一眼不要在写完整个模型才去排查否则会非常痛苦。还有一个容易错的是SOC变量下标。SOC(1)和SOC(0)在Matlab里是index超出问题很多人会把递推式里的SOC(t-1)写成了SOC(t-1)却没给初始值SOC0。需要在约束里先固定SOC(1) SOC0 eta_ch * P_ch(1) * dt - P_dis(1) * dt / eta_dis同时把SOC(1)初始值作为已知参数传入。如果模型里加入了线性化DistFlow潮流约束还要注意节点电压平方变量V_sq不能用V^2否则会产生非线性项导致求解器报Nonlinear constraints错误。正确做法是直接定义V_sq为变量并让V_sq通过潮流方程与支路功率耦合。这是一个非常核心的配网建模细节写错会直接让优化变成NLP求解难度指数级上升。4.3 不确定场景处理与计算性能问题场景数越多求解越慢。如果20个场景的小型网络模型已经要跑几分钟不要急着加场景数量而是先检查是否存在重复的场景。重复场景会增加矩阵规模但不会提供任何新信息。我实际测试过200个原始场景削减到10个和削减到20个目标函数值可能只差个位百分数但求解时间可能差5到10倍。所以场景削减务必保留概率权重不要简单取平均。还有一个隐蔽问题随机场景生成时可能产生负的光伏出力或负的负荷。光伏场景如果乘性误差导致pv_scene出现负值物理上完全不合理。生成场景后必须做钳位处理把负值改成0或者直接丢弃重抽。否则那些离谱场景会把调度结果带偏。4.4 排查速查表现象可能原因处理办法No suitable solver求解器路径没配置运行yalmiptest检查求解成功但结果NaN模型不可行sol.problem非0注释部分约束逐项排查维度不匹配sdpvar行列写反用size()检查变量维度二阶锥报错用V^2建模电压改用V_sq变量场景负值随机抽样越界生成后max(scene, 0)修正第二阶段调整量过大缺少调整上限约束增加dP ub - P约束目标值异常低漏掉第二阶段惩罚项检查目标函数是否包含所有场景5. 从“两阶段随机”到更多实际扩展5.1 转向两阶段鲁棒优化与CCG求解思路随机优化求期望成本鲁棒优化则求最坏场景下的最小成本。如果分布式电源历史数据很少或者你需要调度方案在极端天气下依然可靠可以把现在的随机两阶段模型改造成两阶段鲁棒优化。鲁棒模型通常写成min_x Phase1_cost(x) max_{u∈U} min_y Phase2_cost(x,u,y)。这个三层结构用CCG算法求解主问题先固定一个不确定集子问题寻找当前决策下的最坏场景再把坏场景对应的约束作为割平面加回主问题循环迭代直到收敛。从代码角度现有模型里第二阶段变量其实可以很方便地改成“对每个最坏场景求解”。你只需要把场景生成函数换成一个不确定集描述函数把期望目标换成max-min目标。有了现在的完整代码框架改造的工程量比从零写要小很多。5.2 日内滚动调度与MPC扩展日前计划做完之后现实中还需要日内修正。把代码中的时间窗从24小时缩短到滚动窗口比如每15分钟执行一次、预测未来4小时就是模型预测控制MPC。这种扩展对储能系统尤其重要因为日前计划无法准确预测午后的强光伏波动日内滚动可以及时修正储能出力。一个可行的做法是把当前模型的两阶段“职责”重新分配。日前第一阶段仍然输出基础计划日内MPC调用第二阶段模型但只执行第一个时段的决策执行后滚动到下一个时段重新求解。代码里需要加入状态更新函数把实测SOC和实测负荷传给下一轮优化。5.3 从单目标到多目标加约束、加设备、加时间尺度这套两阶段框架天然支持扩展。想加入需求响应只要在负荷侧增加可削减负荷变量和补偿成本项想加入网络重构增加联络线开关状态0-1变量和对应的拓扑约束想加入碳约束在目标函数里增加碳排放成本或碳配额约束即可。分布式光伏多的时候还可以在第二阶段加入无功调节和电压控制变量把配电网的电压安全抬高到与有功经济性同等重要的位置。我个人在实际调试中体会最深的一点是两阶段模型最怕的不是数学推不出来而是场景和变量之间的“耦合关系”被写错。一旦第二阶段变量没有与第一阶段变量形成合理的调整关系模型就会变成各自独立的24小时调度表失去两阶段的意义。所以拿到任何一套现成代码我都建议先把第二阶段调整量的边界约束找出来看它是不是真的“在计划附近小范围修正”而不是“另起炉灶”。这一步确认了代码基本就吃透了一大半。如果你正在做含分布式电源的调度方向我强烈建议手边保留这样一套两阶段随机规划的完整代码作为基准。后续无论是换求解器、改鲁棒约束、接日内滚动都会有一条清晰的参照线。最简单有效的练法是先把场景数从10改成50跑一次再想想为什么结果变化不大——原因想通了两阶段调度的核心也就通了。