★ FEATURED ARTICLE 做配电网优化调度的朋友应该都有体会开源代码不少但一套能直接跑、还带分布式电源和两阶段调度的完整Matlab代码真不好找。最近我把手头的日前两阶段优化调度模型整理了一遍基于IEEE 33节点配电网加入了分布式光伏、风电和储能用Yalmip建模、Cplex求解把第一阶段的日前计划与第二阶段的实时修正串了起来。这篇文章就把模型原理、数学公式、代码结构和调试时踩过的坑全部摊开讲适合正在做分布式电源接入、配电网经济运行或者准备投文章的同学参考。1. 项目背景与模型思路1.1 为什么配电网调度要采用“日前两阶段”配电网里的分布式电源一多原来的“负荷预测定出力”玩法就不好使了。光伏和风电的出力受天气影响特别大早上还阳光明媚中午一片云飘过来出力能瞬间掉一半。你如果只用一组预测曲线去做24小时调度一旦实际出力偏差大电压越限、线路过载、弃光弃风这些问题马上就会冒出来。这时候就需要“两阶段”思路。第一阶段叫“日前计划”在当天零点之前基于明天的负荷、新能源预测曲线制定24小时各时段的可控DG出力、储能充放电计划、向上级电网的购电计划。第二阶段叫“日内修正”在实时运行中把新能源和负荷的实际值一点一点露出来再对日前计划做最小幅度的调整让系统既满足安全约束又能把成本控制在可接受范围。这个逻辑很像我们出门旅行。查天气预报会先定一个大致的行程表这就是日前计划到了当天遇到临时下雨或者堵车你再微调景点顺序、替换交通工具这就是日内修正。两阶段的好处在于不是把宝全押在预测上而是给不确定性留了缓冲。1.2 两阶段模型整体架构这套模型的整体结构不复杂但每个模块都得扣细节。第一阶段做的是“预决策”变量包括分布式光伏、风电的日前计划出力储能每个时段的充放电功率和SOC轨迹以及配电网根节点向上级电网的购电功率。目标函数是最小化总运行成本约束条件主要是各时段的潮流方程、节点电压上下限、支路电流容量、DG出力上下限、储能SOC递推关系。第二阶段是在第一阶段确定的“基准点”上做修正决策。常见的做法有两种一种是用多场景随机规划生成若干组光伏、风电、负荷场景每个场景下都能调整出力目标函数变成“期望成本最小”另一种是鲁棒优化考虑最坏场景下的可行性和成本。我这版代码用的是场景法好处是物理意义直观Matlab里面用Yalmip写起来也顺手。每个场景下第二阶段决策变量可以向量化表达求解规模可控。1.3 分布式电源与配电网的建模要点先讲分布式电源。光伏和风电在优化里通常当成“负的负荷”或者可控出力电源处理。如果是“不可控”的新能源其实更准确的说法是“可弃电”也就是允许在一定惩罚成本下削减出力。这样模型里就要加弃光弃风变量约束是实际出力不超过预测出力。储能模型则要处理充放电状态互斥、功率上下限、SOC递推以及避免同时充放电的约束。再看配电网。配电网和输电网不一样电阻和电抗比值比较大不能忽略有功损耗潮流计算也更讲究。这里我用了DistFlow支路潮流模型加上二阶锥松弛把非凸潮流约束变成可高效求解的锥约束。Yalmip里可以直接用cone定义锥约束Cplex能原生识别求解速度很快。IEEE 33节点配电网是经典测试算例单辐射状网络带联络开关但我在基础版里先固定开环运行避免整数变量一下子太多先把两阶段调度逻辑跑通再说。2. 数学模型拆解2.1 目标函数第一阶段成本最小化先给第一阶段目标函数。[ \min \sum_{t1}^{24} \left( c_t^{buy} P_{t}^{buy} \sum_{g1}^{n_g} c_g P_{g,t} \sum_{d1}^{n_d} c^{cur} P_{d,t}^{cur} \sum_{b1}^{n_b} c^{bat}\left(P_{b,t}^{dis}P_{b,t}^{ch}\right) \right) ]其中第一项是向上级电网购电成本(c_t^{buy})是分时电价(P_{t}^{buy})是根节点购电功率。第二项是可控DG运行成本通常是燃气轮机或者柴油机发电成本一般建模成线性或分段线性。第三项是弃光弃风惩罚这个系数不能设得太小否则模型会为了省钱疯狂弃掉新能源也不能设得太大否则数值求解容易出问题我一般取500~1000元/MWh具体看你研究场景。第四项是储能充放电成本。严格来说储能本身不“烧钱”但每充放一次电池寿命都有损耗所以我会在目标函数里加一个很小的单位退化成本。注意这里用的是(P_{ch}P_{dis})也就是不管充电还是放电只要动作就有成本这样才能避免模型为了凑约束让储能白白空转。第二阶段的目标函数是在第一阶段基础上对每个随机场景(s)求最小调整成本[ \min \sum_{s} \pi_s \sum_{t1}^{24} \left( c^{adj,} \Delta_{s,t}^{} c^{adj,-} \Delta_{s,t}^{-} \right) ](\Delta^{})和(\Delta^{-})表示实际出力相比日前计划的向上、向下调整量目标就是让实际运行尽量贴着计划走。2.2 约束条件潮流、电压、DG出力、储能SOCDistFlow潮流方程是这套代码的核心对每个节点(j)、每个时段(t)满足[ P_{j,t} P_{i,t} - \sum_{k: j \to k} P_{k,t} - R_{ij} l_{ij,t} - P_{load,j,t} P_{dg,j,t} ][ Q_{j,t} Q_{i,t} - \sum_{k: j \to k} Q_{k,t} - X_{ij} l_{ij,t} - Q_{load,j,t} Q_{dg,j,t} ]这里(i)是父节点(k)是子节点(R_{ij}, X_{ij})是支路阻抗(l_{ij,t})是支路电流幅值平方。节点电压的平方(U_{j,t})通过下面的方程耦合[ U_{j,t} U_{i,t} - 2(R_{ij}P_{ij,t} X_{ij}Q_{ij,t}) \left(R_{ij}^2X_{ij}^2\right) l_{ij,t} ]再加上二阶锥约束[ \left|\begin{bmatrix} 2P_{ij,t} \ 2Q_{ij,t} \ l_{ij,t}-U_{i,t} \end{bmatrix}\right|2 \leq l{ij,t} U_{i,t} ]这个锥约束的作用是把非凸的潮流方程松弛成凸问题。只要目标函数有促使网损变小的项松弛通常都是紧的结果可信。DG约束方面光伏和风电出力不能超过预测值[ 0 \leq P_{dg,d,t} \leq P_{dg,d,t}^{forecast} ]可控DG出力在上下限之间并且爬坡率限制也要加上。储能约束是最容易写错的[ SOC_{b,t1} SOC_{b,t} \eta_{ch} P_{b,t}^{ch} - \frac{P_{b,t}^{dis}}{\eta_{dis}} ][ 0 \leq SOC_{b,t} \leq SOC_{b}^{max} ]这里我额外加了一个“充放电互斥”约束用二进制变量(u_{b,t})表示状态虽然会让模型变成混合整数二阶锥规划但求解器比如Cplex和Gurobi都能搞定。你要是完全不用二进制变量也可以用一个“充电和放电功率乘积为0”的约束但是那样非线性太强不建议。2.3 第二阶段修正与场景生成第二阶段最关键的输入是随机场景。我这里用预测误差模型生成假设光伏和风电的实际出力等于预测值加一个服从正态分布的误差项负荷也类似。然后对每个时段独立抽样再对海量样本做场景削减保留典型场景。我用的场景削减方法是基于概率距离的快速前向选择法从1000个场景里挑出10个代表性场景让它们的概率分布和原始样本的Wasserstein距离最小。Matlab里可以用自带的kmeans聚类近似也可以用scenario工具箱。我代码里用的是自己写的简化版聚类200行左右效果够用。场景数量是关键。太少了模型结果偏乐观太多了求解时间指数上涨。我实测IEEE 33节点配电网5个场景就已经能覆盖大部分不确定性10个场景跑出来的结果和5个差别不大但求解时间翻了一倍以上。所以我默认设成5个场景你们可以按自己的算力调整。3. Matlab代码实现解析3.1 代码总体结构完整代码不是單个脚本而是一个工程文件夹我按职责拆成了下面几部分Case_33bus/ ├── main_optimize.m # 主程序入口 ├── data/ │ ├── load_profile.m # 负荷数据 │ ├── pv_wind_profile.m # 新能源出力预测 │ ├── system_data.m # 线路、节点、DG参数 │ └── price_profile.m # 分时电价 ├── model/ │ ├── build_distflow.m # DistFlow约束 │ ├── build_storage.m # 储能约束 │ ├── build_stage1.m # 第一阶段建模 │ └── build_stage2.m # 第二阶段建模 ├── solve/ │ ├── solve_optimizer.m # 调用YalmipCplex │ └── scenario_reduce.m # 场景削减 ├── result/ │ └── plot_result.m # 绘图与输出主程序就几行把数据加载、建模、求解、结果展示串起来。这种分文件结构的好处是改数据不用翻代码做二次开发也方便。你们拿到代码后最先要改的就是data文件夹里的system_data.m把33节点拓扑改成自己系统。3.2 数据准备与参数设置系统参数和数据不是随便填的里面有不少坑。IEEE 33节点的线路参数我建议统一用有名值基准容量取1 MVA基准电压取12.66 kV这样潮流约束里的电阻电抗数值差别不会太大。你要是不统一单位Cplex解完可能会因为数值病态给你个“infeasible”排查半天发现只是阻抗单位用错了。分时电价我这里设成三个时段峰时1.2元/kWh平时0.7元/kWh谷时0.35元/kWh。分布式光伏预测曲线用了一个夏天晴天出力的典型形状早上6点开始上升中午12点达到峰值下午5点降下来。风电则用一个平稳但有点波动的曲线。负荷数据用IEEE 33节点标准日负荷曲线peak负荷大约5.6 MVA。数据定义用Matlab结构体params.baseMVA 1; params.baseKV 12.66; params.branch [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; ... ]; params.load load_profile(); params.pv pv_wind_profile().pv; params.wind pv_wind_profile().wind; params.price price_profile(); params.horizon 24; params.scenarioNum 5;3.3 基于Yalmip的建模核心代码建模部分我用Yalmip因为它可以用很接近数学表达式的语法写约束维护性比手写大矩阵好太多。下面这段是第一阶段建模的精华。% 定义变量 P_buy sdpvar(24,1); P_pv sdpvar(24,length(pvBus)); P_wind sdpvar(24,length(windBus)); P_ch sdpvar(24,length(batteryBus)); P_dis sdpvar(24,length(batteryBus)); SOC sdpvar(25,length(batteryBus)); u_bat binvar(24,length(batteryBus)); % 充放电状态 % 目标函数 objective sum(price.*P_buy) ... sum(sum(dgCost .* P_g)) ... curCost * (sum(sum(P_pvForecast - P_pv)) sum(sum(P_windForecast - P_wind))) ... batCost * (sum(sum(P_ch)) sum(sum(P_dis))); % 储能SOC递推 for t 1:24 for b 1:nBattery constraints [constraints, ... SOC(t1,b) SOC(t,b) eta_ch*P_ch(t,b) - P_dis(t,b)/eta_dis]; constraints [constraints, ... 0 P_ch(t,b) u_bat(t,b)*P_ch_max(t,b)]; constraints [constraints, ... 0 P_dis(t,b) (1-u_bat(t,b))*P_dis_max(t,b)]; end end这里最容易被忽略的是SOC下标。我用了SOC(25,1)因为24个时段有25个状态点初始SOC是第1个结束SOC是第25个。很多新手写成SOC(24)最后一天的状态递推就会越界。第二阶段建模我采用“场景数组化”的方式把所有场景的变量一次性展开用三维数组存。Yalmip数组索引写起来会麻烦一点但求解时效率高关键是避免for循环里反复调用optimize。你要是一个场景一个场景地调用求解器不仅慢还失去了两阶段模型整体优化的意义。3.4 求解器配置与结果输出求解器我用的是Cplex 12.10通过Yalmip接口调用。核心配置就三行options sdpsettings(... solver,cplex,... verbose,2,... savesolveroutput,1,... cplex.mip.tolerances.mipgap,1e-4);mipgap设到1e-4既保证精度又不让求解器死在整数变量上。如果你们用Gurobi可以把solver改成gurobiYalmip会自动适配。结果输出我主要画四张图24小时购电功率、DG出力曲线、储能SOC曲线、节点电压分布。还有一个表格打印总成本、购电成本、DG成本、弃电惩罚成本。跑完main_optimize.m后工作区里会生成result结构体里面存了所有变量的值方便后续写论文或者做参数分析。4. 运行效果与算例验证4.1 IEEE 33节点算例结果我在一台i5-12400、16GB内存的电脑上跑默认算例5个随机场景第一阶段加第二阶段总共约5000个连续变量、240个整数变量Cplex求解时间约45秒。每次跑完总成本在6500元左右其中购电成本占大头约4800元储能单位退化成本约300元没有发生弃电因为在晴天场景下光伏出力被完整消纳了。节点电压方面未接入DG时33节点配电网末端节点电压偏低大约0.92 p.u.。接了分布式光伏和风电之后末端电压抬升到0.97 p.u.附近个别中午光伏出力高峰时段节点18电压接近1.03 p.u.但没超过上限。这说明分布式电源对电压支撑有明显作用但也带来倒送功率和电压偏高的风险。两阶段模型的意义在这里就体现了——日前计划会提前协调DG出力和储能充电避免中午光伏大发时电压越上限。4.2 两阶段对比分析为了看两阶段到底“值不值”我做了三组对比实验。第一组是纯日前确定性调度不考虑任何不确定性全天使用预测曲线作为真实值。结果总成本最低约6100元但这个成本是“事后诸葛亮”成本因为实际光伏出力不可能完全等于预测值。第二组是日前不考虑不确定性但日内强制按实际场景运行如果不做修正电压越限、功率不平衡会发生必须在第二阶段额外购买调整功率。我加了一个很大的外购惩罚成本模拟紧急调整最终总成本飙到7900元。第三组就是我用的两阶段随机优化模型总成本约6500元。它比第二组少了1400元紧急调整成本只比第一组多了400元“保险成本”。这就是两阶段的优势你多花一点准备成本却避免了不确定性带来的巨额惩罚。如果你写论文这张三组成本对比表基本就是核心结果了。方案日前成本/元日内调整成本/元总成本/元确定性日前610006100日前固定紧急调整610018007900两阶段随机优化580070065004.3 灵敏度分析与扩展场景这套代码还能直接做灵敏度分析。我把光伏渗透率从0.5倍逐步提高到2倍发现总成本先降后升。光伏多了购电成本下降但弃光和电压越限风险上升最终导致惩罚成本增加。渗透率1.2倍左右是当前网络条件下的经济最优值。储能容量也值得测。把储能容量从500 kWh加到2000 kWh总成本下降约6%但继续往上加收益就不明显了因为储能容量受限于充放电功率和配电网络承载能力。这个结论写论文的时候很有用可以在结论里说“适度配置储能最优过度配置边际收益递减”。我还尝试过把第二阶段从场景法改成鲁棒优化用“盒式预算约束”描述不确定性模型会保守一点总成本大概比场景法高10%但结果更鲁棒。如果你想融合这两个方向代码里可以在model/build_stage2.m里替换约束模块。5. 常见问题与调试心得5.1 求解器报错排查清单周围朋友跑这套代码遇到最多的问题我整理成了一张速查表。报错现象常见原因解决办法Yalmip提示没有求解器Cplex未正确安装或路径未添加运行yalmiptest把Cplex目录加入Matlab路径Infeasible problem参数单位不统一或约束过强检查线路阻抗、功率基准值放宽DG出力上下限求解时间过长场景数太多或整数变量爆炸减少场景数到5个把不必要的二进制约束改成连续约束二阶锥约束报错用了而不是cone用cone([2P;2Q;l-U], lU)定义SOC结果不连续充放电效率来回乘除导致数值误差在SOC递推式中统一用pu值避免量纲混用变量名冲突工作区里残留旧变量主程序开头加clear; clc; close all5.2 收敛性调优技巧有两类问题最让人头大一类是模型很长但解不出来一类是解出来了但结果明显不对。对于“解不出来”首先检查二阶锥松弛。DistFlow的二阶锥约束在Yalmip里要写成二阶锥标准形式而不是单纯的不等式。其次I think there is no need to mention that simplified, but its good.第二个经验是Big-M参数不要取得太大。写“充放电互斥”时如果用Big-M表达状态和功率的关系M值设小一点比如功率上限的1.2倍就够了。设成1e6的话Cplex的数值病态会让你怀疑人生。第三如果模型实在收敛慢可以先固定储能SOC初值和终值比如SOC(1)0.2,SOC(25)0.2这样能省去一大部分可行域搜索速度能快30%左右。代价是储能无法参与跨日套利但很多论文本来就会假设调度周期内SOC首末相等所以这样处理是合理的。5.3 代码二次开发建议很多同学拿代码不是为了做复现而是想改造成自己的模型。我的建议是先跑通原版再逐步替换模块。如果你想加需求响应可以在目标函数里增加一个可削减负荷变量约束是削减量不能超过用户合同上限同时给削减成本加一个阶梯价格。这样配电网就不是单纯“源随荷动”而是“源荷互动”。如果你想做三相不平衡配电网需要把DistFlow改成三相解耦形式在IEEE 33节点基础上加变压器中性点模型。这个改动比较大建议至少看懂原版代码里build_distflow.m的每一行再说。如果你只是想换算例网络比如换成IEEE 123节点直接把system_data.m里的线路和负荷矩阵替换就行。Yalmip建模部分用的是节点编号数组只要网络拓扑数据格式一致代码基本不用动。但要留意123节点网络需要配网重构也就是联络开关控制如果不加整数变量结果会有偏差。我个人在实际操作中的体会是两阶段调度模型最花时间的不是建模而是调试随机场景和参数。你花一晚上把约束写对了第二天可能又因为场景削减后概率不为1而掉坑。建议你每写完一个模块都先把对应约束的维度打印出来检查一遍。下面这个小技巧是我一直在用的在optimize之前插入这一行能瞬间定位变量维度和约束数量是否正确。fprintf(Variables: %d, Constraints: %d\n, length(recover(depends(objective))), length(constraints));如果变量数量和约束数量对不上优先查for循环里的下标八成是某一行把t写成了t-1。这个代码里的所有索引我都检查过你们在用的时候重点检查把节点数从33改成其他网络时pvBus、windBus、batteryBus这些下标数组是否越界。最后再分享一个扩展技巧这套两阶段模型已经预留了“滚动时域”接口。你可以把horizon24改成horizon4每个小时重新跑一次只执行第一个时段的决策这就变成分布式电源参与实时调度的MPC框架了。想从“日前规划”升级成“日内滚动优化”的话这是最省事的路径。
阅读完成 · 觉得有帮助？
