简介基于二阶锥规划的主动配电网动态最优潮流代码包面向配电网优化调度研究人员和电力专业研究生重点解决含风电、并联电容器、静止无功发生器及有载调压变压器等多类设备的动态最优潮流建模问题。代码基于MATLAB的YALMIP与CPLEX工具实现通过二阶锥松弛将非凸潮流方程转换为可高效求解的SOCP模型能够覆盖多时段日前调度、无功电压优化等典型应用场景兼顾计算速度与精度。包内共23个文件约117.18MB以两套M源程序为核心并配有PDF与CAJ格式参考文献、PPT讲解课件与配套讲解视频、Word复现过程文档及运行日志便于读者对照源码理解建模细节、复现结果并做扩展研究。目前已有843人学习/下载适合需要掌握主动配电网优化调度算法实现、并希望快速上手的科研人员与高年级学生。1. 从配电网调度说起为什么这一版值得深挖光伏午间满发、储能和可调负荷在同一个时段抢出力配电网调度方案只要晚十分钟刷新就可能多一轮支路过载和电压越限。主动配电网优化调度要同时处理电压、潮流、储能SOC和离散挡位本质上是一个混合整数非线性规划二阶锥规划SOCP是这几年把这类代码在工程上跑通最直接的技术路线——先做变量替换把非线性项变成二次型再松弛成锥约束让商业求解器在分钟级内给出带最优性保证的解。标题里的“代码.rar”落地的就是这条链路。我按建模、代码、数据到避坑的顺序讲适合已经跑得通潮流计算、想自己做调度算例的读者。看完这套方案值不值得复现、坑在哪些地方你自己就能判断。2. 主动配电网调度的模型拆解从DistFlow到二阶锥约束2.1 调度问题到底在优化什么主动配电网与常规配电网最大的区别是“可控设备”从单一的无功补偿扩展到了有功与无功的联合调度。储能充放电功率、光伏逆变器无功出力、可中断负荷、电容器组投切、有载调压变压器OLTC挡位这些设备全部参与优化。对一次完整的日前调度通常取 24 个时段一个典型 33 节点系统就能产生上千个决策变量和约束时序上储能SOC还会把每个时段串成一条长链。把变量按工程角色分组最好理解。连续变量包括节点注入有功、无功储能充放电功率光伏无功离散变量包括储能充放状态、电容器组数、OLTC挡位状态变量则是节点电压幅值平方和支路电流平方。调度要做的事是在满足潮流方程、电压上下限、储能容量、支路载流量约束的前提下让网损、弃光、运行成本加权最小。这里有个值得反复说的工程判断如果只做单断面潮流非线性方程组用牛顿法直接解不需要优化一旦加上 24 小时时序耦合和离散决策问题就变成混合整数非线性规划MINLP商业求解器遇上一个中等规模配电网都会很吃力。二阶锥规划的切入点是保留离散变量的整数性把连续部分的潮流方程做凸松弛得到一个混合整数二阶锥规划MISOCP。理论上的收益是求解难度明显下降工程上的收益是Gurobi、Cplex、Mosek 都有成熟的分支切平面算法跑出来的结果有最优性间隙MIP Gap可以看。这也是为什么近年大量配电网调度代码都转向这个方向。2.2 DistFlow与变量替换松弛的不是精度是求解难度配电网辐射状结构让 DistFlow 支路潮流模型成为主流选择。对一条辐射支路 ij从节点 i 流向节点 j 的有功为 P_ij、无功为 Q_ij支路电阻为 r_ij电抗为 x_iju_i 表示节点电压幅值的平方l_ij 表示支路电流幅值的平方。DistFlow 的核心表达式是三组节点功率平衡方程描述每个节点注入功率与流出功率的关系电压降落方程是 u_j u_i − 2(r_ij × P_ij x_ij × Q_ij) (r_ij² x_ij²) × l_ij支路电流定义是 P_ij² Q_ij² u_i × l_ij。前两组是线性的第三组是二次等式正是非凸的来源。二阶锥松弛的关键是把第三组等式放宽为不等式 P_ij² Q_ij² ≤ u_i × l_ij再结合电压降落方程整理成标准二阶锥形式‖[2×P_ij; 2×Q_ij; u_i − u_j]‖ ≤ u_i u_j。这个写法在多数代码包里直接可见。它去掉电流平方变量 l_ij代价是把电压降落方程中的损耗项 (r_ij² x_ij²) × l_ij 当作小量忽略。对典型的 10kV 馈线这个近似带来的电压误差通常在 0.1% 以内工程上可以接受追求高精度时把 l_ij 作为变量保留再补一条锥约束完整描述损耗即可。为什么这个松弛在辐射网上“紧”直观说辐射网没有环流支路功率由负荷唯一决定目标函数中通常包含网损项而网损正比于 l_ij优化器会倾向把电流平方压到最小因此最优解处松弛不等式自然取等号。这不是经验玄学相关凸松弛紧性的证明已经有大量文献支持。但注意它成立有条件目标里完全不关心网损或网络含环时紧性可能被破坏后面避坑章会专门说。2.3 目标函数的三项网损、弃光与运行成本二阶锥调度代码的目标函数常见由三项组成。第一项是网损写成 sum(r_ij × l_ij)这也是让松弛取等的主要推手第二项是弃光惩罚表达式为 sum(光伏可用功率 − 实际出力)体现消纳意愿第三项是储能和离散设备的运行成本通常按充放电里程或开关次数计。这三项的量纲不统一。网损是功率弃光也是功率储能成本却往往折算成货币或等效里程。常见做法是把所有项都折算成等效电费用权重系数 ω₁、ω₂、ω₃ 控制优先级。有的代码还会把惩罚取平方项目的是让弃光尽量均匀地摊在各时段而不是集中在某一刻但平方项会破坏二阶锥结构把问题变成二次约束规划所以大多数复现版本宁愿保持线性惩罚也不轻易加平方。如果确实要加也建议分两段线性化近似避免求解器类型变化。到这里建模部分就齐了目标函数加潮流锥约束、节点功率平衡、电压限幅、支路载流、储能递推、离散设备约束就是一个完整的 SOCP 调度模型。下一步是把这些数学式变成能跑的代码。3. 把数学式写成可复现代码变量定义、约束拼装与求解器配置3.1 用YALMIP搭建最小SOCP骨架一份示例代码讲解最常见的代码载体是 MATLAB YALMIP因为 YALMIP 对锥约束的建模最直观求解器切换也方便。先按 IEEE 33 节点系统来定义变量节点数 nb33支路数 nl32时段 T24。下面这段是骨架把电压平方、支路功率、节点注入功率都声明为二维矩阵行号对应节点或支路列号对应时段。%% 基本参数 nb 33; nl 32; T 24; from lineData(:,1); to lineData(:,2); r lineData(:,3); x lineData(:,4); %% 决策变量 u sdpvar(nb, T, full); % 节点电压幅值平方 Pij sdpvar(nl, T, full); % 支路有功 Qij sdpvar(nl, T, full); % 支路无功 P sdpvar(nb, T, full); % 节点注入有功 Q sdpvar(nb, T, full); % 节点注入无功 %% 约束容器 Constraints []; %% 支路约束二阶锥 电压降落 for k 1:nl i from(k); j to(k); for t 1:T % 二阶锥松弛||[2P; 2Q; ui-uj]|| uiuj Constraints [Constraints, cone([2*Pij(k,t); 2*Qij(k,t); ... u(i,t)-u(j,t)], u(i,t)u(j,t))]; % 电压降落方程忽略损耗项 Constraints [Constraints, u(j,t) u(i,t) ... - 2*(r(k)*Pij(k,t) x(k)*Qij(k,t))]; end end逻辑说明cone(vector, radius)是 YALMIP 内置二阶锥函数第一参数是多维向量第二参数是一维半径展开后等价于约束向量的二范数不大于半径。把2*Pij、2*Qij和电压差放进同一个向量再把u(i,t)u(j,t)作为半径恰好对应上一章的旋转锥形式。电压降落方程必须保留否则锥约束能单独成立但电压之间失去物理联系。参数说明full声明是普通实矩阵而不是对称矩阵避免 sdpvar 默认的方形变量被当成对称量处理。from/to数组要和数据文件里的编号严格一致如果数据里节点号从 0 开始记得先在 MATLAB 里整体加 1否则索引越界或错位会在后面造成莫名其妙的结果。3.2 储能模型SOC递推与二进制互斥储能在调度代码里是最容易写错的部分。一个储能单元要有充、放电两个功率变量两个充放标志二进制变量一个SOC变量。互斥约束不能只写功率非负因为求解器为了同时满足目标会想办法让充放同时小额度出力所以必须用二进制变量直接限制。%% 储能参数 Pch_max 0.5; Pdis_max 0.5; % 最大充/放功率(MW) E_max 2.0; E_min 0.4; % 容量上下限(MWh) eta_ch 0.95; eta_dis 0.95; % 充/放电效率 E_init 1.0; % 初始SOC SOC_end_min 0.5; % 调度周期结束时的SOC下限 Pch sdpvar(1, T, full); Pdis sdpvar(1, T, full); E sdpvar(1, T1, full); bch binvar(1, T); % 充电状态 bdis binvar(1, T); % 放电状态 Constraints [Constraints, 0 Pch bch * Pch_max]; Constraints [Constraints, 0 Pdis bdis * Pdis_max]; Constraints [Constraints, bch bdis 1]; % 不能同时充放 Constraints [Constraints, E_min E E_max]; Constraints [Constraints, E(1) E_init]; for t 1:T % SOC递推E(t1) E(t) 充电入账 - 放电出账 Constraints [Constraints, E(t1) E(t) ... eta_ch * Pch(t) - Pdis(t) / eta_dis]; end Constraints [Constraints, E(T1) SOC_end_min]; % 末端SOC约束逻辑说明bch * Pch_max的作用是当二进制变量为 0 时强制充电功率为 0为 1 时把上限放开到额定功率。bch bdis 1是互斥约束同时防止充电和放电同时给值。SOC 递推式里充电效率在流入侧乘放电效率在流出侧除这是符合物理习惯的做法。末端约束E(T1) SOC_end_min相当于给储能留“后手”避免优化器把电全部放光。参数说明如果后续要研究储能套利建议把末端约束改成等式E(T1) E_init。这个改动会显著影响调度结果等式约束相当于强迫储能完成一次完整循环目标里的充放收益必须足以覆盖损耗才值得动作不等式约束则允许储能“只放不充”在某些风光场景下会得到看似更省但不可持续的策略。具体采用哪种取决于你要回答的问题是“今日最优”还是“循环寿命内最优”。3.3 求解器配置Gurobi的MIP Gap与数值处理模型拼装完成后求解器参数设置直接决定能否收敛以及结果可信度。SOCP 加二进制变量后求解器要跑分支定界每个节点都要解一个连续二阶锥规划耗时比纯LP高一个量级。常见的设置是给 Gurobi 或 Mosek 一个可以接受的 MIP Gap 和时间上限。%% 目标函数加权网损 弃光惩罚 储能成本 Objective sum(sum(r * ones(1,T) ... )) ; % 示意占位 % 实际更推荐用循环累加写清楚便于核对量纲 %% 求解器设置 ops sdpsettings(solver, gurobi, ... gurobi.mipgap, 0.001, ... gurobi.timelimit, 600, ... gurobi.numericfocus, 1); result optimize(Constraints, Objective, ops);逻辑说明mipgap设成 0.001意思是找到的整数解与下界相对差距在 0.1% 以内就停止。对于日前调度这个精度足够如果把 gap 压到 0求解时间可能从十几分钟变成数小时收益却微乎其微。timelimit是兜底保险避免某个算例卡死在分支树上。numericfocus设为 1让求解器多花一点时间处理数值病态这个参数在处理电压平方和功率量级相差较大的模型时很管用。参数说明Mosek 的优势是纯锥规划数值稳健但分支定界能力不如 Gurobi。如果模型中离散变量多优先 Gurobi如果几乎全是连续变量Mosek 更稳。Cplex 也可用语法上把gurobi.xxx换成cplex.xxx即可。YALMIP 里用sdpsettings(debug,1)可以在求解失败时打印详细诊断这一步比任何代码诊断插件都实际。4. 从rar到跑通解压结构、数据字段与对拍方法4.1 解压后的目录结构与运行环境拿到 rar 包后第一步先看目录结构而不是直接运行。规范代码通常会区分model/、data/、result/三个目录模型函数放一处算例数据放一处输出结果放一处。如果压缩包里只有一个主脚本加几个.mat或.xlsx文件说明作者把数据和模型耦合在一起复现时要小心路径依赖。解压时注意压缩包编码。在 Windows 上用 WinRAR 默认解压一般没问题但在某些环境里用命令行工具解压中文文件名可能出现乱码导致 MATLAB 读取路径失败。我一般会在压缩包上右键查看文件列表确认没有中文名或空格再解压到纯英文路径。如果压缩包有打开密码正确做法是先找发布者确认口令网络上那些所谓一键移除密码的工具既可能捆绑恶意程序也无法保证数据完整性不值得冒这个险。运行环境方面建议用 MATLAB R2020a 以上版本配 YALMIP 最新版和一个商业求解器。只有 YALMIP 没有求解器时模型能建模但优化会报 no solver 错误。YALMIP 自带的sedumi或sdpnal可以兜底跑小算例但 33 节点 24 时段这种规模最好还是上 Gurobi。4.2 数据文件字段到底在说什么调度代码的数据文件通常分四类节点数据、支路数据、分布式电源数据和储能数据。下面这张表是常见字段的对照不同代码命名略有差异但物理含义一致。数据文件典型字段含义在约束中的位置bus 节点表bus_id / Pd / Qd节点编号、有功负荷、无功负荷节点功率平衡方程右侧line 支路表from / to / r / x首端节点、末端节点、电阻、电抗电压降落与锥约束pv 电源表bus_id / Pmax / Qmax接入节点、最大有功、最大无功注入功率上限storage 储能表bus_id / E_max / P_max接入节点、容量、额定功率储能SOC与功率约束读数据时最容易犯的错是把负荷当成“正值”还是“负值”理解反。多数代码约定负荷是消耗功率注入节点的 DG 功率为正值净注入等于 DG 减负荷。如果数据表里 Pd 已经写成了负值再往功率平衡方程里加一个负号结果会是所有节点都在倒送功率电压全线偏高。4.3 跑通之后怎么确认结果没有错很多人第一次跑通代码看到目标函数是正数、有调度曲线就以为成功了。实际上最快的验证办法是对拍。第一步把 DG 出力和储能功率全部设零只留基础负荷跑一遍优化同时用同一个数据文件做一次普通潮流计算对比各节点电压。两者应该基本一致误差来自支路损耗项的近似处理。第二步检查全系统功率平衡每个时段所有节点注入的有功之和减去网损应该接近零。第三步看电压曲线33 节点系统正常运行电压平方在 0.95² 到 1.05² 之间也就是 0.9025 到 1.1025如果出现低于 0.85 或高于 1.15 的值大概率是约束漏写或数据单位错误。对拍结果和目标值也有参照意义。典型 33 节点系统 24 时段网损占总负荷的比例合理区间大致在 2% 到 5%如果算出来是 0.1% 或 20%说明某个约束没有正确生效不是模型“更优”而是模型“更错”。5. 复现避坑二阶锥调度代码最常见的五个坑5.1 现象电压曲线呈锯齿状松弛“不紧”现象是电压幅值在相邻时段之间跳来跳去变化幅度超过 0.03 p.u.而负荷曲线明明很平缓。原因目标函数里没有网损项或者网损项权重为 0锥约束在该时段没有动力取等号。松弛不紧时优化器给出的是一个可行但物理意义偏弱的解电压自由摆动。解决在目标函数里显式加入网损项 sum(r × l)并把权重调到一个可见的正值。如果商业场景里不想改变目标也可以加一个极小的正则项 1e-4 × u 的平方和。调完后重新求解再看电压曲线是否平滑。5.2 现象目标值低得离谱网损竟然为负现象目标函数里网损那部分是负数或者总成本小到不足以覆盖购电成本。原因网损项写成 r × P 而不是 r × l。P 是支路有功量纲不对而且功率在支路上可能是负方向流动累加后正负抵消甚至为负算出来的“网损”失去物理意义。解决检查代码里网损项的写法。正确形式是用支路电流平方 l 乘以电阻 r或者用支路功率平方除以电压平方近似。如果代码里没有定义 l 变量就用锥约束里隐含的关系转换不要直接用 Pij 乘 r。5.3 现象储能SOC在最后时段突然跳回初始值现象调度周期末尾储能SOC曲线出现一个陡峭的上升或下降接近初始值。原因末端约束写成 E(T1) E_init但初始约束或递推关系有一处时段偏移错误比如 E 数组长度是 T 而不是 T1最后一时刻的递推被漏掉求解器无从保证实际最后一个时段的SOC。解决按 3.2 节的写法把 E 声明为 T1 长度递推循环写到 T末端约束作用在 E(T1)。同时打印 E 曲线和 Pch、Pdis 曲线逐时段手算一遍最后一刻的平衡关系确认曲线的突变是递推覆盖产生的假象还是真实调度动作。5.4 现象Excel数据复制进MATLAB后维度错位约束数量翻倍现象约束矩阵维度报错或者求解时间比预期长一倍约束数量明显多于理论值。原因sdpvar 声明变量时用了矩阵形式但在循环里按单时段索引又用了:或者整列赋值产生维度扩展。另一个常见原因是最初把负荷数据做成列向量转置成行向量时漏掉节点顺序和支路顺序错位。解决固定一个约定节点、支路用行索引时段用列索引全程不改。数据读入后先执行assert(size(Pd,1)nb)之类的校验再进入建模循环。如果约束数量比预期大用length(Constraints)打印出来和理论值对比逐段二分找出是哪组约束被意外展开。5.5 现象求解器返回infeasible但模型检查了好几遍都“没问题”现象模型看起来逻辑完整总能找到一处局部调整让它松弛下去但整体一上求解器就报不可行。原因前三类原因占九成——第一某节点电压上限和支路载流量同时卡死储能容量太小没有可行域第二负荷与 DG 数据配错某时段净注入超过线路极限第三储能互斥约束与 SOC 末端等式约束冲突比如 SOC_end_min 设得过高而初始SOC和多时段充放电条件不足以达到。解决不要直接删约束用可行性问题定位。把目标换成常数 0只保留约束求解求解器会给出一个不可行削除集IIS。Gurobi 里有gurobi_iisYALMIP 里可以逐条注释约束二分定位。最快的一步是先把储能末端约束放宽到 0然后逐步收紧观察哪个值开始不可行。这个操作五分钟内能找到边界比反复读代码高效得多。6. 比跑通更进一步用对偶变量判断边界的紧张程度调度代码跑通只是第一步真正落到工程场景里要看哪些约束在限制系统。一个实用技巧是读约束的对偶变量——也就是拉格朗日乘子。在最优解处电压上限约束的对偶值越大说明放宽该节点电压上限对降低调度成本的边际价值越高反过来说这个节点正是电压支撑最薄弱的位置。YALMIP 里取对偶值很简单。先把电压约束单独存成一组再在 optimize 之后调用 dual。示例代码如下v_cons []; for t 1:T v_cons [v_cons, u(:,t) 1.05^2]; v_cons [v_cons, u(:,t) 0.95^2]; end Constraints [Constraints, v_cons]; ops sdpsettings(solver, gurobi, gurobi.mipgap, 0.001); result optimize(Constraints, Objective, ops); lambda dual(v_cons); % 与 v_cons 顺序一致 lam_ub reshape(lambda(1:nb*T), nb, T);逻辑说明dual返回与约束位置一一对应的乘子数组。把乘子按节点和时段画成热力图能一目了然看到哪个节点在哪个时段触碰了边界。对于电压下限约束乘子显著为正的节点就是最优位置安装储能或无功补偿的候选点。这个信息在普通潮流计算里看不出来是优化调度的附加价值。我个人的习惯是每次跑完算例都顺手把对偶值最大的三组约束打出来连同调度曲线放进结果报告里。这样做过几次之后对“约束突然失效”“目标值敏感”这些翻车场面会越来越有体感排查速度也会快不少。希望帮到你。本文还有配套的精品资源点击获取 主动配电网优化调度与二阶锥规划：从DistFlow建模到代码复现避坑
简介基于二阶锥规划的主动配电网动态最优潮流代码包面向配电网优化调度研究人员和电力专业研究生重点解决含风电、并联电容器、静止无功发生器及有载调压变压器等多类设备的动态最优潮流建模问题。代码基于MATLAB的YALMIP与CPLEX工具实现通过二阶锥松弛将非凸潮流方程转换为可高效求解的SOCP模型能够覆盖多时段日前调度、无功电压优化等典型应用场景兼顾计算速度与精度。包内共23个文件约117.18MB以两套M源程序为核心并配有PDF与CAJ格式参考文献、PPT讲解课件与配套讲解视频、Word复现过程文档及运行日志便于读者对照源码理解建模细节、复现结果并做扩展研究。目前已有843人学习/下载适合需要掌握主动配电网优化调度算法实现、并希望快速上手的科研人员与高年级学生。1. 从配电网调度说起为什么这一版值得深挖光伏午间满发、储能和可调负荷在同一个时段抢出力配电网调度方案只要晚十分钟刷新就可能多一轮支路过载和电压越限。主动配电网优化调度要同时处理电压、潮流、储能SOC和离散挡位本质上是一个混合整数非线性规划二阶锥规划SOCP是这几年把这类代码在工程上跑通最直接的技术路线——先做变量替换把非线性项变成二次型再松弛成锥约束让商业求解器在分钟级内给出带最优性保证的解。标题里的“代码.rar”落地的就是这条链路。我按建模、代码、数据到避坑的顺序讲适合已经跑得通潮流计算、想自己做调度算例的读者。看完这套方案值不值得复现、坑在哪些地方你自己就能判断。2. 主动配电网调度的模型拆解从DistFlow到二阶锥约束2.1 调度问题到底在优化什么主动配电网与常规配电网最大的区别是“可控设备”从单一的无功补偿扩展到了有功与无功的联合调度。储能充放电功率、光伏逆变器无功出力、可中断负荷、电容器组投切、有载调压变压器OLTC挡位这些设备全部参与优化。对一次完整的日前调度通常取 24 个时段一个典型 33 节点系统就能产生上千个决策变量和约束时序上储能SOC还会把每个时段串成一条长链。把变量按工程角色分组最好理解。连续变量包括节点注入有功、无功储能充放电功率光伏无功离散变量包括储能充放状态、电容器组数、OLTC挡位状态变量则是节点电压幅值平方和支路电流平方。调度要做的事是在满足潮流方程、电压上下限、储能容量、支路载流量约束的前提下让网损、弃光、运行成本加权最小。这里有个值得反复说的工程判断如果只做单断面潮流非线性方程组用牛顿法直接解不需要优化一旦加上 24 小时时序耦合和离散决策问题就变成混合整数非线性规划MINLP商业求解器遇上一个中等规模配电网都会很吃力。二阶锥规划的切入点是保留离散变量的整数性把连续部分的潮流方程做凸松弛得到一个混合整数二阶锥规划MISOCP。理论上的收益是求解难度明显下降工程上的收益是Gurobi、Cplex、Mosek 都有成熟的分支切平面算法跑出来的结果有最优性间隙MIP Gap可以看。这也是为什么近年大量配电网调度代码都转向这个方向。2.2 DistFlow与变量替换松弛的不是精度是求解难度配电网辐射状结构让 DistFlow 支路潮流模型成为主流选择。对一条辐射支路 ij从节点 i 流向节点 j 的有功为 P_ij、无功为 Q_ij支路电阻为 r_ij电抗为 x_iju_i 表示节点电压幅值的平方l_ij 表示支路电流幅值的平方。DistFlow 的核心表达式是三组节点功率平衡方程描述每个节点注入功率与流出功率的关系电压降落方程是 u_j u_i − 2(r_ij × P_ij x_ij × Q_ij) (r_ij² x_ij²) × l_ij支路电流定义是 P_ij² Q_ij² u_i × l_ij。前两组是线性的第三组是二次等式正是非凸的来源。二阶锥松弛的关键是把第三组等式放宽为不等式 P_ij² Q_ij² ≤ u_i × l_ij再结合电压降落方程整理成标准二阶锥形式‖[2×P_ij; 2×Q_ij; u_i − u_j]‖ ≤ u_i u_j。这个写法在多数代码包里直接可见。它去掉电流平方变量 l_ij代价是把电压降落方程中的损耗项 (r_ij² x_ij²) × l_ij 当作小量忽略。对典型的 10kV 馈线这个近似带来的电压误差通常在 0.1% 以内工程上可以接受追求高精度时把 l_ij 作为变量保留再补一条锥约束完整描述损耗即可。为什么这个松弛在辐射网上“紧”直观说辐射网没有环流支路功率由负荷唯一决定目标函数中通常包含网损项而网损正比于 l_ij优化器会倾向把电流平方压到最小因此最优解处松弛不等式自然取等号。这不是经验玄学相关凸松弛紧性的证明已经有大量文献支持。但注意它成立有条件目标里完全不关心网损或网络含环时紧性可能被破坏后面避坑章会专门说。2.3 目标函数的三项网损、弃光与运行成本二阶锥调度代码的目标函数常见由三项组成。第一项是网损写成 sum(r_ij × l_ij)这也是让松弛取等的主要推手第二项是弃光惩罚表达式为 sum(光伏可用功率 − 实际出力)体现消纳意愿第三项是储能和离散设备的运行成本通常按充放电里程或开关次数计。这三项的量纲不统一。网损是功率弃光也是功率储能成本却往往折算成货币或等效里程。常见做法是把所有项都折算成等效电费用权重系数 ω₁、ω₂、ω₃ 控制优先级。有的代码还会把惩罚取平方项目的是让弃光尽量均匀地摊在各时段而不是集中在某一刻但平方项会破坏二阶锥结构把问题变成二次约束规划所以大多数复现版本宁愿保持线性惩罚也不轻易加平方。如果确实要加也建议分两段线性化近似避免求解器类型变化。到这里建模部分就齐了目标函数加潮流锥约束、节点功率平衡、电压限幅、支路载流、储能递推、离散设备约束就是一个完整的 SOCP 调度模型。下一步是把这些数学式变成能跑的代码。3. 把数学式写成可复现代码变量定义、约束拼装与求解器配置3.1 用YALMIP搭建最小SOCP骨架一份示例代码讲解最常见的代码载体是 MATLAB YALMIP因为 YALMIP 对锥约束的建模最直观求解器切换也方便。先按 IEEE 33 节点系统来定义变量节点数 nb33支路数 nl32时段 T24。下面这段是骨架把电压平方、支路功率、节点注入功率都声明为二维矩阵行号对应节点或支路列号对应时段。%% 基本参数 nb 33; nl 32; T 24; from lineData(:,1); to lineData(:,2); r lineData(:,3); x lineData(:,4); %% 决策变量 u sdpvar(nb, T, full); % 节点电压幅值平方 Pij sdpvar(nl, T, full); % 支路有功 Qij sdpvar(nl, T, full); % 支路无功 P sdpvar(nb, T, full); % 节点注入有功 Q sdpvar(nb, T, full); % 节点注入无功 %% 约束容器 Constraints []; %% 支路约束二阶锥 电压降落 for k 1:nl i from(k); j to(k); for t 1:T % 二阶锥松弛||[2P; 2Q; ui-uj]|| uiuj Constraints [Constraints, cone([2*Pij(k,t); 2*Qij(k,t); ... u(i,t)-u(j,t)], u(i,t)u(j,t))]; % 电压降落方程忽略损耗项 Constraints [Constraints, u(j,t) u(i,t) ... - 2*(r(k)*Pij(k,t) x(k)*Qij(k,t))]; end end逻辑说明cone(vector, radius)是 YALMIP 内置二阶锥函数第一参数是多维向量第二参数是一维半径展开后等价于约束向量的二范数不大于半径。把2*Pij、2*Qij和电压差放进同一个向量再把u(i,t)u(j,t)作为半径恰好对应上一章的旋转锥形式。电压降落方程必须保留否则锥约束能单独成立但电压之间失去物理联系。参数说明full声明是普通实矩阵而不是对称矩阵避免 sdpvar 默认的方形变量被当成对称量处理。from/to数组要和数据文件里的编号严格一致如果数据里节点号从 0 开始记得先在 MATLAB 里整体加 1否则索引越界或错位会在后面造成莫名其妙的结果。3.2 储能模型SOC递推与二进制互斥储能在调度代码里是最容易写错的部分。一个储能单元要有充、放电两个功率变量两个充放标志二进制变量一个SOC变量。互斥约束不能只写功率非负因为求解器为了同时满足目标会想办法让充放同时小额度出力所以必须用二进制变量直接限制。%% 储能参数 Pch_max 0.5; Pdis_max 0.5; % 最大充/放功率(MW) E_max 2.0; E_min 0.4; % 容量上下限(MWh) eta_ch 0.95; eta_dis 0.95; % 充/放电效率 E_init 1.0; % 初始SOC SOC_end_min 0.5; % 调度周期结束时的SOC下限 Pch sdpvar(1, T, full); Pdis sdpvar(1, T, full); E sdpvar(1, T1, full); bch binvar(1, T); % 充电状态 bdis binvar(1, T); % 放电状态 Constraints [Constraints, 0 Pch bch * Pch_max]; Constraints [Constraints, 0 Pdis bdis * Pdis_max]; Constraints [Constraints, bch bdis 1]; % 不能同时充放 Constraints [Constraints, E_min E E_max]; Constraints [Constraints, E(1) E_init]; for t 1:T % SOC递推E(t1) E(t) 充电入账 - 放电出账 Constraints [Constraints, E(t1) E(t) ... eta_ch * Pch(t) - Pdis(t) / eta_dis]; end Constraints [Constraints, E(T1) SOC_end_min]; % 末端SOC约束逻辑说明bch * Pch_max的作用是当二进制变量为 0 时强制充电功率为 0为 1 时把上限放开到额定功率。bch bdis 1是互斥约束同时防止充电和放电同时给值。SOC 递推式里充电效率在流入侧乘放电效率在流出侧除这是符合物理习惯的做法。末端约束E(T1) SOC_end_min相当于给储能留“后手”避免优化器把电全部放光。参数说明如果后续要研究储能套利建议把末端约束改成等式E(T1) E_init。这个改动会显著影响调度结果等式约束相当于强迫储能完成一次完整循环目标里的充放收益必须足以覆盖损耗才值得动作不等式约束则允许储能“只放不充”在某些风光场景下会得到看似更省但不可持续的策略。具体采用哪种取决于你要回答的问题是“今日最优”还是“循环寿命内最优”。3.3 求解器配置Gurobi的MIP Gap与数值处理模型拼装完成后求解器参数设置直接决定能否收敛以及结果可信度。SOCP 加二进制变量后求解器要跑分支定界每个节点都要解一个连续二阶锥规划耗时比纯LP高一个量级。常见的设置是给 Gurobi 或 Mosek 一个可以接受的 MIP Gap 和时间上限。%% 目标函数加权网损 弃光惩罚 储能成本 Objective sum(sum(r * ones(1,T) ... )) ; % 示意占位 % 实际更推荐用循环累加写清楚便于核对量纲 %% 求解器设置 ops sdpsettings(solver, gurobi, ... gurobi.mipgap, 0.001, ... gurobi.timelimit, 600, ... gurobi.numericfocus, 1); result optimize(Constraints, Objective, ops);逻辑说明mipgap设成 0.001意思是找到的整数解与下界相对差距在 0.1% 以内就停止。对于日前调度这个精度足够如果把 gap 压到 0求解时间可能从十几分钟变成数小时收益却微乎其微。timelimit是兜底保险避免某个算例卡死在分支树上。numericfocus设为 1让求解器多花一点时间处理数值病态这个参数在处理电压平方和功率量级相差较大的模型时很管用。参数说明Mosek 的优势是纯锥规划数值稳健但分支定界能力不如 Gurobi。如果模型中离散变量多优先 Gurobi如果几乎全是连续变量Mosek 更稳。Cplex 也可用语法上把gurobi.xxx换成cplex.xxx即可。YALMIP 里用sdpsettings(debug,1)可以在求解失败时打印详细诊断这一步比任何代码诊断插件都实际。4. 从rar到跑通解压结构、数据字段与对拍方法4.1 解压后的目录结构与运行环境拿到 rar 包后第一步先看目录结构而不是直接运行。规范代码通常会区分model/、data/、result/三个目录模型函数放一处算例数据放一处输出结果放一处。如果压缩包里只有一个主脚本加几个.mat或.xlsx文件说明作者把数据和模型耦合在一起复现时要小心路径依赖。解压时注意压缩包编码。在 Windows 上用 WinRAR 默认解压一般没问题但在某些环境里用命令行工具解压中文文件名可能出现乱码导致 MATLAB 读取路径失败。我一般会在压缩包上右键查看文件列表确认没有中文名或空格再解压到纯英文路径。如果压缩包有打开密码正确做法是先找发布者确认口令网络上那些所谓一键移除密码的工具既可能捆绑恶意程序也无法保证数据完整性不值得冒这个险。运行环境方面建议用 MATLAB R2020a 以上版本配 YALMIP 最新版和一个商业求解器。只有 YALMIP 没有求解器时模型能建模但优化会报 no solver 错误。YALMIP 自带的sedumi或sdpnal可以兜底跑小算例但 33 节点 24 时段这种规模最好还是上 Gurobi。4.2 数据文件字段到底在说什么调度代码的数据文件通常分四类节点数据、支路数据、分布式电源数据和储能数据。下面这张表是常见字段的对照不同代码命名略有差异但物理含义一致。数据文件典型字段含义在约束中的位置bus 节点表bus_id / Pd / Qd节点编号、有功负荷、无功负荷节点功率平衡方程右侧line 支路表from / to / r / x首端节点、末端节点、电阻、电抗电压降落与锥约束pv 电源表bus_id / Pmax / Qmax接入节点、最大有功、最大无功注入功率上限storage 储能表bus_id / E_max / P_max接入节点、容量、额定功率储能SOC与功率约束读数据时最容易犯的错是把负荷当成“正值”还是“负值”理解反。多数代码约定负荷是消耗功率注入节点的 DG 功率为正值净注入等于 DG 减负荷。如果数据表里 Pd 已经写成了负值再往功率平衡方程里加一个负号结果会是所有节点都在倒送功率电压全线偏高。4.3 跑通之后怎么确认结果没有错很多人第一次跑通代码看到目标函数是正数、有调度曲线就以为成功了。实际上最快的验证办法是对拍。第一步把 DG 出力和储能功率全部设零只留基础负荷跑一遍优化同时用同一个数据文件做一次普通潮流计算对比各节点电压。两者应该基本一致误差来自支路损耗项的近似处理。第二步检查全系统功率平衡每个时段所有节点注入的有功之和减去网损应该接近零。第三步看电压曲线33 节点系统正常运行电压平方在 0.95² 到 1.05² 之间也就是 0.9025 到 1.1025如果出现低于 0.85 或高于 1.15 的值大概率是约束漏写或数据单位错误。对拍结果和目标值也有参照意义。典型 33 节点系统 24 时段网损占总负荷的比例合理区间大致在 2% 到 5%如果算出来是 0.1% 或 20%说明某个约束没有正确生效不是模型“更优”而是模型“更错”。5. 复现避坑二阶锥调度代码最常见的五个坑5.1 现象电压曲线呈锯齿状松弛“不紧”现象是电压幅值在相邻时段之间跳来跳去变化幅度超过 0.03 p.u.而负荷曲线明明很平缓。原因目标函数里没有网损项或者网损项权重为 0锥约束在该时段没有动力取等号。松弛不紧时优化器给出的是一个可行但物理意义偏弱的解电压自由摆动。解决在目标函数里显式加入网损项 sum(r × l)并把权重调到一个可见的正值。如果商业场景里不想改变目标也可以加一个极小的正则项 1e-4 × u 的平方和。调完后重新求解再看电压曲线是否平滑。5.2 现象目标值低得离谱网损竟然为负现象目标函数里网损那部分是负数或者总成本小到不足以覆盖购电成本。原因网损项写成 r × P 而不是 r × l。P 是支路有功量纲不对而且功率在支路上可能是负方向流动累加后正负抵消甚至为负算出来的“网损”失去物理意义。解决检查代码里网损项的写法。正确形式是用支路电流平方 l 乘以电阻 r或者用支路功率平方除以电压平方近似。如果代码里没有定义 l 变量就用锥约束里隐含的关系转换不要直接用 Pij 乘 r。5.3 现象储能SOC在最后时段突然跳回初始值现象调度周期末尾储能SOC曲线出现一个陡峭的上升或下降接近初始值。原因末端约束写成 E(T1) E_init但初始约束或递推关系有一处时段偏移错误比如 E 数组长度是 T 而不是 T1最后一时刻的递推被漏掉求解器无从保证实际最后一个时段的SOC。解决按 3.2 节的写法把 E 声明为 T1 长度递推循环写到 T末端约束作用在 E(T1)。同时打印 E 曲线和 Pch、Pdis 曲线逐时段手算一遍最后一刻的平衡关系确认曲线的突变是递推覆盖产生的假象还是真实调度动作。5.4 现象Excel数据复制进MATLAB后维度错位约束数量翻倍现象约束矩阵维度报错或者求解时间比预期长一倍约束数量明显多于理论值。原因sdpvar 声明变量时用了矩阵形式但在循环里按单时段索引又用了:或者整列赋值产生维度扩展。另一个常见原因是最初把负荷数据做成列向量转置成行向量时漏掉节点顺序和支路顺序错位。解决固定一个约定节点、支路用行索引时段用列索引全程不改。数据读入后先执行assert(size(Pd,1)nb)之类的校验再进入建模循环。如果约束数量比预期大用length(Constraints)打印出来和理论值对比逐段二分找出是哪组约束被意外展开。5.5 现象求解器返回infeasible但模型检查了好几遍都“没问题”现象模型看起来逻辑完整总能找到一处局部调整让它松弛下去但整体一上求解器就报不可行。原因前三类原因占九成——第一某节点电压上限和支路载流量同时卡死储能容量太小没有可行域第二负荷与 DG 数据配错某时段净注入超过线路极限第三储能互斥约束与 SOC 末端等式约束冲突比如 SOC_end_min 设得过高而初始SOC和多时段充放电条件不足以达到。解决不要直接删约束用可行性问题定位。把目标换成常数 0只保留约束求解求解器会给出一个不可行削除集IIS。Gurobi 里有gurobi_iisYALMIP 里可以逐条注释约束二分定位。最快的一步是先把储能末端约束放宽到 0然后逐步收紧观察哪个值开始不可行。这个操作五分钟内能找到边界比反复读代码高效得多。6. 比跑通更进一步用对偶变量判断边界的紧张程度调度代码跑通只是第一步真正落到工程场景里要看哪些约束在限制系统。一个实用技巧是读约束的对偶变量——也就是拉格朗日乘子。在最优解处电压上限约束的对偶值越大说明放宽该节点电压上限对降低调度成本的边际价值越高反过来说这个节点正是电压支撑最薄弱的位置。YALMIP 里取对偶值很简单。先把电压约束单独存成一组再在 optimize 之后调用 dual。示例代码如下v_cons []; for t 1:T v_cons [v_cons, u(:,t) 1.05^2]; v_cons [v_cons, u(:,t) 0.95^2]; end Constraints [Constraints, v_cons]; ops sdpsettings(solver, gurobi, gurobi.mipgap, 0.001); result optimize(Constraints, Objective, ops); lambda dual(v_cons); % 与 v_cons 顺序一致 lam_ub reshape(lambda(1:nb*T), nb, T);逻辑说明dual返回与约束位置一一对应的乘子数组。把乘子按节点和时段画成热力图能一目了然看到哪个节点在哪个时段触碰了边界。对于电压下限约束乘子显著为正的节点就是最优位置安装储能或无功补偿的候选点。这个信息在普通潮流计算里看不出来是优化调度的附加价值。我个人的习惯是每次跑完算例都顺手把对偶值最大的三组约束打出来连同调度曲线放进结果报告里。这样做过几次之后对“约束突然失效”“目标值敏感”这些翻车场面会越来越有体感排查速度也会快不少。希望帮到你。本文还有配套的精品资源点击获取