做配电网优化调度的同学这两年应该没少听“主动配电网”和“源-荷-储协同互动”这两个词。说实话我刚接触这个方向的时候也花了不少时间才把论文里的公式搬到Matlab里跑通。IEEE33节点作为最经典的中压辐射状配网测试系统配上一套合理的源-荷-储协同调度策略几乎是每一位做配网优化的人绕不开的练功房。这篇博文就围绕“考虑源-荷-储协同互动的主动配电网优化调度研究【IEEE33节点】Matlab代码实现”这个项目把整个模型的构建思路、数学表达式、Matlab代码实现过程以及我在实际调试中踩过的坑完整梳理一遍。适合刚接触配电网优化的研究生、做园区或台区能量管理的工程师以及一切打算用Matlab复现同类论文的朋友。1. 项目背景与问题定位源-荷-储协同互动到底解决什么问题1.1 主动配电网与传统配电网的区别传统配电网本质上是一个被动网络结构上类似城市自来水管网电源在上游用户在下游能量单向流动调度员能做的只是被动地跟随负荷变化最多调调变压器分接头、切切电容。分布式光伏大规模接入之后情况变了配电网末端开始出现电源潮流不再全是单向的电压波动、过载、倒送电的问题接踵而至。如果还用老办法管理末端电压越限、馈线过流会频繁出现。主动配电网的“主动”二字核心是它把分布式电源、储能、可控负荷当成可管理的资源通过优化调度主动调节网络状态。传统配电网只能被动响应主动配电网则能主动调配源、荷、储资源让网络在经济性、安全性、消纳可再生能源能力上都上一个台阶。1.2 源、荷、储各自的角色与协同目标在IEEE33节点三个角色里源指的是光伏、风电这类分布式电源。它们的共同特点是出力波动大、难以精确预测调度时要结合预测曲线给一个合理的出力区间。荷除了固定负荷还应该包含可平移负荷、可控负荷比如电动汽车充电桩、空调负荷、可定时运行的家电。储指的是电池储能系统它是整个协同框架里的灵活性核心。三者协同的目标我一般归纳成三点经济性最优调度要让综合运行成本最小包括主网购电成本、分布式电源运维成本、储能充放电损耗成本和弃风弃光惩罚。安全性节点电压保持在0.95~1.05pu范围内支路电流不超过上限变压器不过载。消纳能力尽量让光伏、风电全量上网避免不必要的弃置。所以源-荷-储协同互动不是把每台设备单独控制而是整体协调储能低谷充电、高峰放电光伏和风电在电压安全范围内尽量多发可平移负荷避开负荷高峰。这样整个电网的运行曲线才能被“削峰填谷”地拉平。1.3 为什么用IEEE33节点做验证IEEE33节点测试系统是配电网优化领域最常用的标准算例。它共33个节点、32条支路基准电压12.66kV基准功率10MVA总有功负荷约3715kW无功负荷约2300kVar。系统虽然不大但包含了辐射状配网最典型的几个问题负荷分布不均、末端电压偏低、支路潮流差异大。很多初学者可能不知道这个系统的原始状态其实并不“健康”不加任何DG和储能时33节点的末端电压大约只有0.90~0.93pu左右低于电网运行规程要求的0.95pu。正因如此它特别适合用来验证优化调度的效果——接入光伏、储能并合理调度后末端电压能明显回到0.95以上效果直观可见。另外IEEE33节点数据公开、文献多、参数表容易核对是入行配网优化最合适的第一站。2. 优化调度模型的设计思路数学建模是重中之重2.1 目标函数怎么定我见过很多初学者的误区一上来就写代码结果模型边界不清、约束矛盾调半天不知道问题在哪。正确的思路是先写数学模型再翻译成代码。这个项目的场合设为日前24小时优化调度步长1h。假设在IEEE33节点上接入两个光伏电站节点17和节点32附近总容量约600kW、一个风电场节点21附近容量约400kW、两套电池储能节点7和节点24每套额定功率200kW容量1MWh。目标函数可以写成这样一个单目标加权和min F 购电成本 DG运行成本 储能损耗成本 弃风弃光惩罚 网损折算成本逐项拆开购电成本主网分时电价乘以每个时段的购电功率。峰时电价约1.2元/kWh平时约0.7元/kWh谷时约0.3元/kWh。这个电价结构是储能“低储高放”的经济驱动力没有峰谷电价储能的调度意义会小很多。DG运行成本光伏和风电的运维成本通常很小但要在目标函数里占一个位置否则优化器可能为了降网损随意压光伏出力。储能损耗成本充放电过程中电池的循环损耗用一个较小的成本系数折算避免储能频繁动作。弃风弃光惩罚如果预测的可发功率没有被消纳就要付出一个较高的惩罚成本比如2.0元/kWh。这样优化器会尽可能让DG满发。网损折算成本把网损电量乘以电价折算成钱让模型主动降低网络损耗。这种多目标转加权单目标的方法简单直接也是文献里最常见的处理方式。比“再跑一次多目标帕累托”省事很多而且工程上很实用。2.2 约束条件逐条拆解约束条件是整个模型的骨架我建议按下面几个模块来组织第一个模块潮流约束。这是配网优化和普通经济调度的最大区别。辐射状配电网通常使用DistFlow方程然后利用二阶锥松弛SOCP把非凸问题转成凸问题。对每条支路和每个时段方程可以写成有功平衡支路流入的有功减去下游支路流出有功和线损等于该节点净负荷原负荷减去DG出力再加储能充电功率减去放电功率。无功平衡形式与有功类似。电压关系支路末端电压的平方等于首端电压平方减去两倍的有功与电阻乘积加无功与电抗乘积再加上电流平方项。二阶锥约束关键的一步写成norm([2*P_ij; 2*Q_ij; l_ij - v_i], 2) l_ij v_i这种形式。这里有个重要的理论保证在辐射状配电网中二阶锥松弛是精确的也就是说松弛后的最优解就是原问题的最优解不需要担心松弛引入了误差。这也是为什么我强烈推荐用SOCP而不是用启发式算法的原因之一。第二个模块电压和电流安全约束。节点电压幅值平方的上下限0.95pu到1.05pu支路电流的平方也有上限通常按线路载流量设定。如果潮流解出来后某个节点电压偏低优化器会自动调整DG无功出力或储能充放电策略。第三个模块DG出力约束。光伏和风电的有功出力在0到预测可用功率之间无功出力在逆变器容量约束范围内可调。这里我建议给无功一个合理的调节范围因为分布式电源参与调压是主动配电网的重要功能。第四个模块储能系统约束。这是代码实现时的重点和难点主要包括SOC(t1) SOC(t) η_ch·P_ch(t)·Δt/E_cap - P_dis(t)·Δt/(η_dis·E_cap)其中η_ch和η_dis为充放电效率通常取0.95E_cap为电池容量。还需要限制同一时刻不能同时充放电用两个二进制变量u_ch和u_dis加约束u_ch u_dis 1完成充放电功率不超过额定值荷电状态保持在0.1到0.9之间调度周期初末SOC相等体现储能的长期可持续运行。第五个模块可平移负荷约束。为了让“荷”侧真正参与互动我假设总负荷中约10%是可平移的平移后一天内的总用电量不变每个时段可平移量有上限。这个约束在数学上很好写但实际效果很能说明需求响应如何削峰填谷。2.3 时序模拟与滚动优化思路日前24小时调度看起来只是一个静态优化问题但理解它的时间耦合关系很重要。储能SOC让不同时段之间产生耦合所以模型本质上是带时间耦合的优化问题不是24个独立的静态问题。很多新手没有意识到这一点如果把每个时段单独优化储能就会在每个时段只做对自己最有利的动作整体上完全失去“低储高放”的协调性。正确做法是把24个时段的变量全部放在一个模型里同时求解这样储能才能站在全天视角做决策。如果再往后走一步可以在日前调度的基础上加日内滚动修正比如每15分钟滚动一次用最新的光伏预测数据刷新模型。但第一步先把日前调度做扎实后面扩展会非常顺。3. Matlab代码实现全流程解析3.1 环境搭建与工具选型Matlab平台做配网优化我建议用YALMIP建模型再调用商业求解器Cplex或Gurobi。YALMIP是一个建模工具箱你只需要用sdpvar定义变量、用Constraints [Constraints, ...]堆约束、用optimize调用求解器剩下的模型转换和求解都由它代劳。为什么不直接用fmincon或粒子群因为这个问题本质是非凸的虽然二阶锥松弛后是凸的但用通用非线性求解器很容易陷入局部最优而且对初值极其敏感。粒子群一类的启发式算法虽然也能跑但每次运行结果都不一样论文复现时很容易被审稿人质疑。Cplex和Gurobi对SOCP和MISOCP都有原生支持求解速度快结果可复现。建议先安装YALMIP再去IBM官网申请Cplex学术版或Gurobi学术版然后在Matlab里setpath添加路径。安装完成后用yalmiptest命令测一下有没有正确识别求解器。下面是求解器对比表方便大家选型求解器类型学术免费SOCP支持备注Cplex商业是强IBM申请配网优化最常用Gurobi商业是更强申请流程略复杂性能顶尖SCS开源是支持适合SOCP混合整数能力弱ECOS开源是支持轻量适合中等规模fminconMatlab自带是不支持非凸不建议用在这个场景实际项目中我最常用的是Cplex稳定性高报错信息也友好。3.2 数据准备IEEE33节点参数与典型日曲线Matlab代码的第一步永远是准备数据。IEEE33节点的拓扑数据一般包括支路编号、首端节点、末端节点、电阻、电抗。网上常见的数据格式是% 支路编号 首端 末端 R(ohm) X(ohm) branch_data [ 1 1 2 0.0922 0.0470 2 2 3 0.4930 0.2511 3 3 4 0.3660 0.1864 4 4 5 0.3811 0.1941 5 5 6 0.8190 0.7070 % ... 共32条 ];节点负荷数据通常也按节点编号对应列出有功单位kW无功单位kVar。一个常见的坑是根节点的编号有些数据把变电站根节点记为0有些记为1。两种都能算但一定要统一不然支路首末端对应关系全错。我在代码里习惯让根节点是1变量编号从1开始。负荷曲线、光伏出力曲线和风电出力曲线可以手动录入典型日数据也可以用Beta分布、Weibull分布来随机生成。对于验证模型来说手动录一份平滑的典型日曲线就够了。要特别注意单位换算基准功率是10MVA负荷是kW换算到标幺值时需要除以10000算错一个单位电压结果会离谱到让你怀疑人生。3.3 核心代码框架从变量定义到求解调用我把核心代码框架整理成下面几段完整代码量比较大这里给出骨架和关键约束的写法大家拿到手之后代入完整数据就能跑。第一步是定义决策变量T 24; n_bus 33; n_branch 32; % 电压平方、电流平方、支路有功、支路无功 V2 sdpvar(n_bus, T, full); % 节点电压幅值平方 I2 sdpvar(n_branch, T, full); % 支路电流幅值平方 P_branch sdpvar(n_branch, T, full); Q_branch sdpvar(n_branch, T, full); % 主网购电功率 P_buy sdpvar(T, 1, full); % 分布式电源出力 P_pv sdpvar(npv, T, full); % npv为光伏电站数量 P_wt sdpvar(nwt, T, full); % 储能变量 P_ch sdpvar(n_ess, T, full); % 充电功率 P_dis sdpvar(n_ess, T, full); % 放电功率 SOC sdpvar(n_ess, T, full); % 荷电状态 u_ch binvar(n_ess, T); % 充电状态 u_dis binvar(n_ess, T); % 放电状态第二步是搭目标函数Objective sum(price .* P_buy) ... % 购电成本 sum(sum(c_pv * P_pv)) ... % 光伏运维成本 sum(sum(c_wt * P_wt)) ... % 风电运维成本 sum(sum(c_ess * (P_ch P_dis))) ... % 储能损耗成本 sum(sum(c_cur * (P_pv_max - P_pv))) ... % 弃光惩罚 sum(sum(c_cur * (P_wt_max - P_wt))); % 弃风惩罚第三步是堆约束。这里挑两个最有代表性的约束说明写法。首先是二阶锥潮流约束Constraints [Constraints, V2(:, 1) 1.0]; % 根节点电压 for k 1:n_branch i branch_data(k, 2); j branch_data(k, 3); for t 1:T % DistFlow 电压方程 Constraints [Constraints, V2(j, t) V2(i, t) ... - 2 * (R(k) * P_branch(k, t) X(k) * Q_branch(k, t)) / Vbase ... % 注意单位一致化 (R(k)^2 X(k)^2) * I2(k, t)]; % 二阶锥松弛 Constraints [Constraints, ... norm([2 * P_branch(k, t); 2 * Q_branch(k, t); ... I2(k, t) - V2(i, t)], 2) I2(k, t) V2(i, t)]; end end然后是储能SOC约束dt 1; % 时间步长单位h E_cap 1000; % 容量单位kWh % 初始SOC Constraints [Constraints, SOC(:, 1) 0.5]; for t 1:T-1 Constraints [Constraints, SOC(:, t1) SOC(:, t) ... eta_ch * P_ch(:, t) * dt / E_cap ... - P_dis(:, t) * dt / (eta_dis * E_cap)]; end % 充放电互斥与功率上限 Constraints [Constraints, P_ch u_ch * P_max]; Constraints [Constraints, P_dis u_dis * P_max]; Constraints [Constraints, u_ch u_dis 1]; % SOC上下限与周期末回到初始值 Constraints [Constraints, 0.1 SOC 0.9]; Constraints [Constraints, SOC(:, end) 0.5];最后调用求解器ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); optimize(Constraints, Objective, ops); % 提取结果 V2_opt value(V2); P_ess value(P_dis) - value(P_ch); P_buy_opt value(P_buy);我想强调一个细节约束里所有变量必须统一单位。很多人的程序跑出来电压全是0.9左右甚至负值十有八九是阻抗单位、功率单位或者基准值没统一。IEEE33节点里电阻电抗的单位是欧姆而电压用标幺值这时候要用基准阻抗Zbase换算不然潮流方程里的各项数量级对不上。3.4 结果输出与绘图让数据自己说话求解完成后必须把结果可视化否则你根本不知道调度策略有没有起作用。我一般至少画四张图第一张是24小时主网购电功率曲线对比优化前和优化后。优化后曲线的高峰会明显下降低谷会稍微抬升这就是储能削峰填谷的直接证据。第二张是各节点电压幅值对比图。我习惯取负荷最高的时刻比如晚上19点画一条随节点编号变化的电压曲线优化前末端电压可能只有0.92优化后应该能压在0.95以上。这张图是最能说明调度效果的核心图。第三张是储能SOC和充放电功率图。你会看到SOC在电价低谷爬升、高峰下降与分时电价高度对应说明储能的确在响应电价信号。第四张是DG出力曲线。光伏、风电基本满发说明弃风弃光惩罚起了作用。绘图代码很简单figure; plot(1:33, sqrt(V2_opt(:, 19)), b-o); hold on; plot(1:33, sqrt(V2_base(:, 19)), r--s); xlabel(节点编号); ylabel(电压(pu)); legend(优化后, 优化前); grid on;这张电压对比图建议做成论文的第一个结果图审稿人一眼就能看到方法的有效性。4. 常见问题与排查技巧实录4.1 求解器报错与安装问题“No suitable solver found”是最常见的报错。原因很简单YALMIP没找到可用的求解器。解决办法是确认Cplex或Gurobi已安装路径添加到了Matlab搜索路径中并且用yalmiptest测试通过。我遇到过一个很奇怪的情况Cplex装了但yalmiptest里Cplex显示红色不支持SOCP后来发现是Cplex版本太老换到12.10以上版本就好了。如果模型带二进制变量YALMIP会调用MILP或MISOCP求解器这时Cplex是首选。如果你用了开源求解器如SCS模型里一旦出现binvarYALMIP很可能直接报错因为SCS不支持整数变量。4.2 SOC越界或结果振荡储能SOC越界是另一个高发问题。我最初调试时遇到过SOC在某个时段突然跳到0或者超过1一查发现是SOC连续性约束里dt写错了。这里dt是小时数如果步长是1小时dt1没有问题如果步长是15分钟dt0.25就必须改。另外初始SOC和末尾SOC相等这个约束你要小心如果初始值设得过高比如0.9而SOC上限也是0.9那模型在第一个时段就只能放电不能充电整个调度策略就变形了。我一般把初始SOC设在0.5给储能留足双向调节空间。还有一个细节SOC上下限不能设得太紧。我见过有人把SOC上限定为0.8结果储能几乎没有调节能力电压也没抬起来。建议用0.1~0.9既有安全裕度又能发挥储能作用。4.3 单位制混乱的问题这个问题我要单独拿出来说因为它坑过太多人。IEEE33节点的基准功率是10MVA基准电压是12.66kV由这两个值算出的基准阻抗为Zbase Ubase^2 / Sbase 12.66^2 / 10 16.02 ohm支路阻抗的标幺值等于欧姆值除以16.02。如果你直接在DistFlow方程里用欧姆单位的电阻而电压用标幺值那各项数量级完全对不上电压解出来就不是0.95~1.05的范围。负荷数据同样如此一个节点100kW的有功在标幺值体系里是0.01pu100kW / 10000kVA。很多人在这一步忘记除以基准功率导致模型中负荷大了一百倍电压直接跌穿下限储能再怎么调都救不回来。我的建议是代码最前面直接做好全部数据标幺化并且加注释写明基准值。宁可多写几行转换代码也不要中途混着单位算。4.4 模型规模大、求解慢的处理33节点24时段的模型本身不大但如果后面扩展到几百个节点、多时间尺度或更多储能求解时间会明显增加。如果Cplex求解时间过长我一般按下面顺序处理先去掉储能互斥的二进制变量拉长求解时间限制观察结果差异把时间尺度从24小时缩短为代表性4个时段峰、平、谷、尖峰做初步验证设置sdpsettings(cplex.mip.tlimits, 300)限制求解时间得到次优可行解检查是不是约束写得过紧比如每个DG的无功上下限设得太窄导致可行域变小不要在模型里加不必要的非线性项比如效率用分段线性函数表示避免引入过多整数变量。下面是一个常见问题速查表方便翻查现象可能原因解决办法YALMIP找不到求解器未安装Cplex/Gurobi或路径未添加添加路径运行yalmiptest验证电压结果全部偏低或为负数单位制混乱基准阻抗/功率换算错误统一用标幺值检查ZbaseSOC越界或曲线突变dt步长错误或初始SOC设得不合理核对时间步长初始SOC设0.5储能不充不放电价差太小或充放电互斥约束卡死检查分时电价检查u_chu_dis约束求解时间过长二进制变量过多或约束过松缩短时间范围限制求解时间上限连续运行结果不同用了启发式算法或并行随机性改用Cplex/Gurobi确定性求解节点编号越界告警根节点编号0/1不一致统一节点编号体系并核对支路表4.5 一个关于二阶锥约束的独家经验二阶锥约束的写法很多人容易写错。YALMIP里用norm函数是最方便的但不能写成sqrt(...)加平方项YALMIP会自动识别norm并转换为锥约束。如果手动写成或一大串求解器可能识别不了凸性。还有一点二阶锥约束不可逆。也就是说不能反过来写成norm(...) ...那样是典型的非凸约束Cplex会直接拒绝求解。之前有同学问我为什么模型里加了一个看起来“更精确”的约束反而报错原因就是这个。5. 结果分析与扩展方向让研究价值最大化5.1 结果曲线应该怎么分析代码跑通只是第一步怎么解读结果才是论文或者项目报告的核心。拿到优化结果后推荐从四个维度分析第一是经济性维度优化后的购电成本比优化前降低了多少储能通过峰谷套利创造了多少收益这部分可以列一个成本对比表清晰直观。第二是安全性维度所有节点电压是否都进入0.95~1.05pu区间重载支路是否还有过载风险如果末端电压仍然偏低说明DG容量或储能安装位置需要调整可以尝试更换接入节点再跑一次。第三是消纳维度光伏和风电的消纳率是否接近100%如果存在弃风弃光要看是哪个时段、哪个节点再针对性调整惩罚系数。第四是储能运行维度SOC曲线是否平滑储能是否频繁切换充放电如果一天内充放电次数过多会增加电池寿命损耗可以在目标函数里加大储能损耗成本权重。5.2 模型可以往哪些方向扩展这个项目做好之后可以非常自然地往下面几个方向扩展考虑不确定性光伏出力预测不可能完全准确可以把确定性模型改成两阶段鲁棒优化或者用机会约束规划处理预测误差多时间尺度协同在日前调度的基础上加日内滚动修正用15分钟级别的刷新频率应对突发波动配电网重构联合优化把联络开关状态也作为决策变量问题从SOCP升级为MISOCP依然可以用Cplex求解多目标优化把碳排放、电压偏差、DG消纳率放入目标函数用加权求和或epsilon约束法求帕累托前沿分布式求解与多主体调度把储能、DG、负荷聚合商作为独立利益主体用ADMM算法做分布式协调这是当前园区级能量管理系统的热点扩展到更大系统IEEE33只是起点可以替换成IEEE123节点系统或者其他实际馈线数据模型框架完全不用改。我个人在做这类项目时最大的体会是数学模型一定要先于代码变量和约束在纸上列清楚了写代码就是翻译工作。另外跑通一个确定性模型之后别急着做复杂扩展先把结果图、参数敏感性分析补齐把“为什么储能装在这里”“为什么这个时段充电”这类问题解释清楚研究的价值自然就出来了。最后再分享一个小习惯我习惯把每次实验的求解器日志、目标函数值和关键结果都存成一个结构化表格换参数对比时会轻松很多。这个项目做完基本就掌握了配网优化调度的完整技能链后面再往综合能源、多时间尺度、不确定性方向走都会事半功倍。