做配电网优化的人手里迟早得有一套动态最优潮流DOPF程序。我最近被问到最多的就是光伏、储能、调压器都塞进模型之后原来能解的静态 OPF 开始疯狂报 infeasible好不容易把二阶锥松弛写进去又不知道怎么调求解器参数。这个项目标题里的两个词其实已经把路线点明白了——二阶锥规划、主动配电网动态最优潮流。翻译成大白话就是在含分布式电源、储能、有载调压变压器等主动设备的辐射状配电网里把未来 24 小时每个时段的潮流调度一次算出来目标函数可以是网损最小、电压偏差最小或者 DG 消纳最大而约束项里既有非凸的交流潮流方程又有储能 SOC 这类跨时段耦合项。这篇文章我会从一个 IEEE 33 节点主动配电网案例出发完整走一遍“建模 - 二阶锥松弛 - YALMIP 建模 - CPLEX 求解 - 结果校验”的流程。重点不是甩一个跑完就完事的代码包而是讲清楚每个关键写法背后的原因为什么 DistFlow 要那样改锥、为什么标幺化会坑掉一大半新手、为什么求解器报 infeasible 时先别怀疑求解器。适合正在做配电网优化、分布式资源调度、微电网能量管理方向的同学也适合刚接触优化求解、想快速搭建一个可扩展框架的人。1. 动态最优潮流到底难在哪1.1 “动态”两个字让变量直接翻倍静态最优潮流只算一个截面比如某时刻网损最小。动态最优潮流则要把 T 个时段全部串起来典型是 24 个整点或 96 个 15 分钟断面。最直接的冲击就是变量数量膨胀假设一个 33 节点的网络支路 32 条每个时段光支路有功、无功、电流平方、节点电压平方这几个核心变量就接近 200 个24 个时段直接奔着四五千个连续变量去再加入储能 SOC、OLTC 档位这类状态量规模轻松过万。这对求解器来说还不是致命问题真正麻烦的是变量之间有了时间耦合。储能是这种耦合的典型。某时段充电会影响后续所有时段的 SOC 可行域OLTC 一天调节次数有限属于跨时段的计数约束分布式电源出力受气象和日前预测影响本质上也要按时间序列描。方程之间环环相扣你已经不是在求“某一个时刻”的最优解而是在求“一整条可行轨迹”的最优解。这就像排班单独安排一个人某天休息很容易要让一个班组连续一周每天都有人值班且总工时最小就是完全不同的复杂度。如果只算单时段很多配置可以拍脑袋给一旦放到 24 时段任何一段不行整体就不可行。很多初学者第一次跑动态 OPF 时程序报 infeasible其实不是潮流方程写错而是储能首末 SOC 约束太紧、或者负荷高峰时段和光伏高峰时段错配导致整个调度周期根本串不出一条可行轨迹。1.2 潮流约束的非凸性才是求解器的拦路虎维度增加还只是“量变”真正让最优潮流变得棘手的是潮流方程本身的非凸性。在 DistFlow 支路方程里节点电压幅值平方 U、支路电流平方 L、有功和无功 P、Q 之间存在关系L (P² Q²) / U。这个式子把三个变量揉在一起而且除以 U 意味着这是一个非凸区域。非凸优化模型的麻烦在于局部极小值不一定是全局极小值通用求解器完全没有能力判断自己停下来的点到底是不是全局最优。你用内点法直接怼这个非线性问题大概率得到一堆局部最优甚至不收敛的结果而且没人敢说这个结果行不行。于是从业者们想到了二阶锥松弛SOCP Relaxation。思路不复杂先把非凸等式 L (P²Q²)/U 放宽成不等式 L ≥ (P²Q²)/U也就是允许电流平方比物理值更大。这个不等式经过代数变形后能写成标准二阶锥约束。形成凸可行域后再用内点法或分支割平面法求全局最优。这个操作相当于把一个坑坑洼洼的山地变成了一大片规整的下坡路求解器终于能大步往前走了。再叠加上 OLTC 档位、电容器组投切、储能充放电逻辑这些离散变量问题最终会变成混合整数二阶锥规划MISOCP。这也是这个项目标题里为什么点名 CPLEX 的原因——CPLEX 是处理混合整数规划和二阶锥规划都成熟的老牌商用求解器YALMIP 负责建模转换CPLEX 负责在锥约束和整数约束之间来回切效率比纯开源求解器高不少。2. 二阶锥松弛怎么把非凸问题变成凸问题2.1 DistFlow 方程和它的等价变形二阶锥松弛不能套在任意潮流方程上它有一套完整的适用前提——辐射状配电网的支路潮流方程。这套方程在配电网优化里被叫作 DistFlow本质上是从基尔霍夫定律推出来的只保留了有功、无功、电压幅值平方和电流幅值平方四个量。对一个从节点 i 到节点 j 的支路标准形式是三条等式首先是节点功率平衡进入节点 j 的有功功率加上本地发电等于本地负荷和流出功率之和无功也同理。然后是电压降方程末端电压平方 U_j 等于首端电压平方 U_i 减去线路上的压降再补上电流平方引起的修正。最后是电流定义式L_ij (P_ij² Q_ij²) / U_i。前两条等式是线性的问题全部集中在第三条上。P²、Q²、除以 U三者放在一起直接破坏了凸性。松弛操作就是把这第三条等式改成不等式L_ij ≥ (P_ij² Q_ij²) / U_i。把它变个形左右两边同时乘 U_i得到 U_i · L_ij ≥ P_ij² Q_ij²。这是一个旋转二阶锥rotated quadratic cone的标准形式。再经过一次变量替换可以写成更利于代码实现的 2-范数形式|| [2P_ij; 2Q_ij; U_i - L_ij] ||₂ ≤ U_i L_ij。这个形式在 YALMIP 里可以直接照抄成一条约束CPLEX 内部会把它当作原生锥约束处理不会粗暴地展开成一堆二次项。实际写代码的时候我见过不少人手动展开写成平方和不等式不仅费劲而且容易触发 CPLEX 的二次约束非凸检查。2.2 松弛“精不精确”怎么判断和处理有人会问你把等式放宽成不等式那求出来的还是原问题的最优解吗答案要看情况。对于辐射状配电网如果目标函数是网损最小这类“压着线走”的目标松弛后得到的解通常刚好落在锥边界上也就是说 L_ij 和 (P_ij² Q_ij²)/U_i 相等这时候二阶锥松弛是精确的求出的结果就是原非凸问题的全局最优解。但事情没这么绝对。当网里有大量 DG 注入尤其是无功功率边界被压得很紧的时候锥松弛可能会被“撑开”约束没有全激活解出来的 L_ij 明显大于物理上需要的最小电流平方。此时回代校验会发现误差这个情况在专业文献里叫“松弛间隙”duality gap不为零。我自己排查这类问题时最常用的方法是一段十行的回代脚本把求解得到的 U、P、Q 重新代入 L (P²Q²)/U算一个绝对误差超过 1e-4 甚至 1e-3 就要警惕了。如果确实有间隙最常见的修补方法是给目标函数加一个关于 L_ij 的小惩罚项比如把目标改成网损 ε · Σ L让求解器有动力把锥压回边界。ε 通常取 1e-6 到 1e-4具体数值靠你实测去调。加完惩罚之后要重新检查间隙如果还不行就得检查是不是无功边界给得太别扭、或者目标函数的各项权重尺度差太多。2.3 主动配电网的各类设备约束怎么塞进 SOCP把二阶锥当作“平台”之后剩下来的工作就是往这个框架里放设备模型。最基础的是分布式电源光伏和风机的有功出力通常有个上限无功功率容量可以近似成一个圆形或者梯形约束写成线性不等式就够用。微型燃气轮机和燃料电池也类似有功无功都在一个凸多边形里。储能是动态问题里最核心的一段。它的 SOC 递推是线性方程SOC(t1) SOC(t) dt · (η_ch · P_ch(t) - P_dis(t)/η_dis) / E。这里 P_ch 和 P_dis 分别代表充电、放电功率上界是储能额定功率SOC 本身有上下限。如果想让模型更严谨还需要禁止同一时段既充电又放电这就要引入二进制变量模型从 SOCP 变成 MISOCP。OLTC 一类带档位的设备天然是整数变量。YALMIP 里用 intvar 声明CPLEX 负责把离散档位和连续潮流耦合在一起。可投切电容器组也类似用整数变量表示投入组数。可以说模型的最终形态基本就是“线性约束 二阶锥约束 离散整数变量”三件套这也是为什么 YALMIP CPLEX 这套组合在这个场景下这么流行它俩对这三类约束都支持得最完整。3. YALMIPCPLEX完整实现与参数坑3.1 环境准备和测试求解器先把环境装利索。MATLAB 方面我用的版本是 R2023a 以上YALMIP 直接从 GitHub 拉最新版放到 MATLAB 搜索路径里即可。CPLEX 需要装对应平台的版本并且确保 MATLAB 能通过cplex命令访问到求解器。这个环节最常见的坑不是安装而是 MATLAB 和 CPLEX 的位数/版本不匹配导致 yalmiptest 里 CPLEX 那一栏显示红色。所以我每次搭好环境第一件事就是跑一遍yalmiptest。这个函数会把 YALMIP 支持的所有求解器挨个测一遍能识别出来的会显示 OK。但凡 CPLEX 不 OK先别急着改代码回到环境配置上排查。路径配好了、版本对得上了再跑一个极小的二阶锥算例比如二维的 min ||x||₂验证整条链路通不通。这一步花不了两分钟但能省掉后面一晚上的调试时间。3.2 数据准备、标幺化和数据结构数据准备看着简单实际是最容易出错的一环。IEEE 33 节点系统的基准电压是 12.66 kV基准功率常取 10 MVA。线路参数如果原始数据给的是欧姆和千瓦必须先换算成标幺值再进模型。基准阻抗 Z_base V_base² / S_base 12.66² / 10 ≈ 16.02 Ω。很多新手把线路 R、X 零点几个欧姆直接丢进模型算出来的结果偏向离谱根源全在单位混用。我自己的习惯是把负荷、光伏出力、线路参数全部放在一个结构体里按 t 时刻排列。比如 P_load 是 24 × 33 矩阵第 t 行第 j 列表示 t 时段节点 j 的有功负荷P_pv 是 24 × 33 的 PV 最大可用出力矩阵R、X 和 bus_from、bus_to 是支路拓扑数据。这样做的好处是所有约束组装都能严格对齐时间维度和节点维度不会出现维度错位后求解器神秘报错的情况。3.3 变量定义与核心约束代码下面这段代码是整套模型的骨架。我故意省略了太多细节只保留最容易写歪的地方方便你对照自己手里的工程。决策变量都用 sdpvar 声明维度跟着时间、支路、节点走T 24; n_bus 33; n_branch 32; n_storage 2; % 假设节点18和节点25接入储能 % 连续变量 P_br sdpvar(T, n_branch, full); % 支路首端有功 Q_br sdpvar(T, n_branch, full); % 支路首端无功 Lij sdpvar(T, n_branch, full); % 支路电流平方 U sdpvar(T, n_bus, full); % 节点电压平方 % 储能变量 P_ch sdpvar(T, n_storage, full); P_dis sdpvar(T, n_storage, full); SOC sdpvar(T, n_storage, full); % OLTC档位引入整数变量模型变为MISOCP tap intvar(1, T, full);核心约束组装分成三块。第一块是支路方程和二阶锥松弛第二块是节点功率平衡第三块是储能递推。DistFlow 和锥约束的写法如下Constraints []; for t 1:T for k 1:n_branch i bus_from(k); j bus_to(k); % DistFlow 电压降方程 Constraints [Constraints, ... U(t,i) - U(t,j) 2*(R(k)*P_br(t,k) X(k)*Q_br(t,k)) ... - (R(k)^2 X(k)^2)*Lij(t,k)]; % 二阶锥松弛核心中的核心 Constraints [Constraints, ... norm([2*P_br(t,k); 2*Q_br(t,k); U(t,i)-Lij(t,k)]) ... U(t,i) Lij(t,k)]; end end注意这里norm([...]) ...就是标准二阶锥约束。如果你在旧资料里看到cone函数那是 YALMIP 的另一套写法本质一样。我更喜欢 norm 写法的原因很简单一眼能看懂它在表达什么调试时不容易把维度搞错。节点功率平衡部分先按拓扑判断每个节点的父支路和子支路然后列平衡方程。下面是只含负荷和储能的示意光伏或 DG 再加对应的注入项即可for t 1:T for j 1:n_bus inflow 0; outflow 0; for k 1:n_branch if bus_to(k) j % 支路末端注入节点j的功率 首端功率 - 线路损耗 inflow inflow P_br(t,k) - R(k)*Lij(t,k); end if bus_from(k) j outflow outflow P_br(t,k); end end % 本地注入项储能充电是负荷放电是发电 P_gen 0; for s 1:n_storage if storage_bus(s) j P_gen P_gen P_dis(t,s) - P_ch(t,s); end end Constraints [Constraints, ... inflow - outflow P_load(t,j) - P_gen - P_pv(t,j)]; end end无功功率平衡的写法一模一样只需要把有功换成无功、把 R(k)*Lij 换成 X(k)*Lij再补上无功负荷、光伏无功上限、储能无功为 0 这几个条件。我说“一模一样”是真的一模一样复制粘贴再改变量名就行。储能递推约束注意时刻对齐和效率系数E 1.0; % 储能容量MWh eta_ch 0.95; % 充电效率 eta_dis 0.95; % 放电效率 Pmax 0.2; % 储能额定功率MW for s 1:n_storage Constraints [Constraints, SOC(1,s) 0.5]; for t 2:T Constraints [Constraints, ... SOC(t,s) SOC(t-1,s) ... (eta_ch*P_ch(t-1,s) - P_dis(t-1,s)/eta_dis) / E]; end Constraints [Constraints, SOC(:,s) 1, SOC(:,s) 0.2]; Constraints [Constraints, SOC(T,s) 0.5]; % 充放电功率限幅 Constraints [Constraints, 0 P_ch(:,s) Pmax, 0 P_dis(:,s) Pmax]; % 若加二进制变量阻断同时充放电 % u binvar(T,1,full); % Constraints [Constraints, P_ch(:,s) Pmax*u, P_dis(:,s) Pmax*(1-u)]; end这段代码里我注释掉了同时充放电的二进制变量。实际建模时如果项目要求必须符合物理时序建议把那段注释取消如果只是先验证算法框架SOCP 连续松弛跑得快但解里很可能出现同一时段又充又放的理论解回代工程现场会被质疑。这个取舍要看你的场景是学术计算还是工程可实施。3.4 求解器参数设置和常见状态码模型组装完目标函数和求解器设置直接决定最终能不能跑出好结果。以网损最小为例目标函数很简单Objective sum(sum(Lij .* R)); % 各支路电流平方乘电阻再累加所有时段YALMIP 里点乘是逐元素乘Lij 和 R 的维度需要匹配。R 如果只有 1×32要用 repmat 或广播展开成 24×32。这里不展开但代码里容易因为维度问题报错。求解器设置我一般写成options sdpsettings(solver, cplex, verbose, 2); options.cplex.mip.tolerances.mipgap 0.01; options.cplex.optimalitytarget 3; options.cplex.timelimit 600; sol optimize(Constraints, Objective, options);mipgap 0.01的意思是最优性间隙允许 1%对工程场景通常够了能大幅缩短求解时间。optimalitytarget 3这条设置很关键尤其是在旧版 CPLEX 里遇到锥约束和整数变量混合的模型不设置这个值可能触发“quadratic constraint non-convex”的误报。timelimit防止模型卡死超过十分钟直接停方便你回头优化模型而不是干等。求解结束之后一定看一眼sol.problemif sol.problem 0 disp(求解成功); else disp(sol.info); % 返回具体错误信息 end0 表示成功非 0 就要结合sol.info排查。常见的问题码包括不可行、无界、数值问题等。我很少背数字表直接把sol.info打出来YALMIP 会给出人话解释。3.5 结果提取与合理性检查求解完别忘了把变量还原成物理量否则一堆标幺值账根本对不上。电压幅值就是sqrt(value(U))电流平方直接value(Lij)储能 SOC 曲线value(SOC)。这些量建议一次性存成结构体方便后续画图。看结果的时候我有一套固定动作先画 24 时段全网电压曲线看所有节点电压是否落在 0.95~1.05 之间再画储能 SOC 曲线看首末值是否满足约束最后检查支路电流是否越限。如果前两步没问题但最后电流越了限说明支路容量约束没加或者加错了。另外要习惯用“状态信息”而不是“目标函数值”来判断求解质量。目标函数 0.05 和 0.052 看起来差不多但如果求解状态在 MIP 节点上挣扎半天、gap 一直降不下去你就要考虑是不是整数变量和连续变量的量级差太多导致数值病态。此时通常是调单位或调参数而不是盲目加约束。4. 常见问题与排查技巧实录4.1 求解器返回 infeasible问题基本出在数据第一次跑动态 OPF十个人有八个会碰到 infeasible。这时候先别怀疑二阶锥写错因为我见过的案例里九成以上都是数据问题一个是单位没换算。线路阻抗用欧姆、负荷用千瓦但基准容量取的 MVA电流平方和电压平方的公式全乱套。检查方法很简单取一个纯负荷且无 DG 的时段手算一下根节点注入功率和平衡节点的电压至少保证量纲一致。一个是拓扑方向搞反。DistFlow 是沿“根节点 - 叶子节点”的方向写的如果某条支路的 bus_from 和 bus_to 反了电压降方程符号错功率平衡也跟着错连锁反应就是整个约束系统自相矛盾。还有一个是储能约束太紧。初始 SOC 明明只有 0.3又要求 24 小时后回到 0.5中间还有一大段负荷高峰根本没有充电窗口整段时间轨迹就断了。解决办法是先放宽 SOC 上下限、去掉末状态约束等模型能解之后再逐步收紧逐步逼近真实要求。排 infeasible 最有效的手段是“二分法删约束”。把储能、OLTC、DG 一段段注释掉直到模型能解再放回来缩小问题范围。这个手段比盯着报错日志大海捞针高效得多。我自己的习惯是先把动态问题退化成单时段静态问题单时段能解再慢慢加时间耦合项加一个检查一次基本能定位到是哪条约束把模型锁死了。4.2 锥松弛间隙大电压曲线失真怎么办有时候模型能解目标函数也收敛了但电压曲线在重负荷时段出现反常的剧烈凹陷或者支路电流平方和 P²Q²/U 对不上。这是二阶锥松弛间隙大的典型征兆。我自己的检查脚本长这样P_ value(P_br); Q_ value(Q_br); U_ value(U); L_ value(Lij); gap max(max(abs(L_ - (P_.^2 Q_.^2) ./ U_))); disp([最大锥松弛间隙 , num2str(gap)]);间隙如果超过 1e-4就要处理。最常见的修复方案是往目标函数里加一个关于 Lij 的小系数惩罚项把锥“拽”回边界。不过这个惩罚项的系数不能太大太大相当于把总损耗目标都压掉了解出来的调度方案对网损不再敏感。我从 1e-6 开始试逐步放大每次重新求解并检查 gap直到 gap 小于阈值。另一个有效技巧是给电压变量 U 提供一个高质量的初始值。先用普通潮流工具算出 24 个时段的基准电压在 YALMIP 求解前用assign赋予初值再让 CPLEX 从这组初值开始求解。对锥松弛类问题热启动效果很明显求解器不需要在寻找初始可行解上花太多功夫锥约束也更容易切到边界。4.3 求解速度慢多半不是求解器的锅如果你发现 33 节点的 24 时段模型跑了五分钟以上还不出结果先不要觉得是 CPLEX 不行大概率是模型本身太“胖”。最典型的原因是把约束一股脑写进一个巨大的循环里YALMIP 需要为每条约束生成内部表达式几万条约束堆下来建模时间占比比求解时间还大。优化方向有几个。第一把循环里的重复表达式抽出来比如支路阻抗平方 R²X² 提前算好。第二尽量避免在约束里写sum(sum(...))这种维度广播过大的表达式改成矩阵运算代替循环。第三合理使用 YALMIP 的批量约束写法比如Constraints [Constraints, U 1.05^2, U 0.95^2]这种整矩阵约束比双层循环逐一 append 快得多。如果模型里带了 OLTC 和储能充放电二进制变量那问题就变成 MISOCP分支定界过程本身就可能比较慢。这时可以合理放宽 mipgap。我在 33 节点算例上把 mipgap 从 1e-4 放宽到 1e-3速度提升好几倍目标函数只多了不到 0.2%工程上完全可接受。4.4 用回代测试给结果“验真”最后一步也是我最坚持的一步用二阶锥解算出来的调度方案重新跑一遍精确的交流潮流回代验证电压和支路潮流是否真的满足物理方程。这不是多此一举因为“二阶锥松弛精确”是一种理想情况实际模型里加了很多设备近似可能继承了一点误差。回代的方式有现成工具MATPOWER 的runpf就能给一组 PQ 和 PV 节点算精确潮流。你也可以自己写一个简单的牛顿法去核对关键节点电压工作量不大。回代后如果最大电压偏差超过 0.005 p.u.我就要回头检查锥松弛间隙和各个设备约束的近似程度。这步做扎实了后续论文里的解、工程可实施性都有底。我曾经踩着个坑某次程序跑完目标函数很好看电压也在限值内但回代后某些节点电压直接掉到 0.91。查了一天发现问题出在光伏无功约束我用了线性近似而实际逆变器的无功能力曲线是圆的模型松弛掉了一块导致最优解刚好落在被松弛掉的那块区域。从那以后回代校验成了我所有优化代码出结果后的固定动作。在我个人实际跑了大量配电网优化算例之后有一个体会特别深这套 YALMIP CPLEX 二阶锥松弛的组合最大价值不是“能算”而是可扩展性强。从 33 节点换到 123 节点从单目标网损改成多目标加权电压偏差从确定性负荷换成随机场景基本只需要改数据、改约束、改目标函数不用推翻整体架构。代码里的求解参数、数据单位、校验脚本最好从一开始就用结构体和脚本固定下来后续折腾新模型时能省掉大量返工。最后再分享一个小技巧如果你发现二阶锥松弛在某个边界算例上老是翻车不妨把锥约束从norm([2P;2Q;U-L]) UL这种显式写法换成旋转锥写法再用 YALMIP 的cone函数有时候不同写法会让内部预处理的路径完全不同莫名就稳了。别每次都只背一种写法多试几种你对自己模型的理解会深很多。希望这篇内容能帮你在动态最优潮流这条路上少走几周弯路。 
长河君 · 2026/10/10 15:04:12 · 阅读 1.2k 二阶锥松弛与YALMIP+CPLEX实现主动配电网动态最优潮流
做配电网优化的人手里迟早得有一套动态最优潮流DOPF程序。我最近被问到最多的就是光伏、储能、调压器都塞进模型之后原来能解的静态 OPF 开始疯狂报 infeasible好不容易把二阶锥松弛写进去又不知道怎么调求解器参数。这个项目标题里的两个词其实已经把路线点明白了——二阶锥规划、主动配电网动态最优潮流。翻译成大白话就是在含分布式电源、储能、有载调压变压器等主动设备的辐射状配电网里把未来 24 小时每个时段的潮流调度一次算出来目标函数可以是网损最小、电压偏差最小或者 DG 消纳最大而约束项里既有非凸的交流潮流方程又有储能 SOC 这类跨时段耦合项。这篇文章我会从一个 IEEE 33 节点主动配电网案例出发完整走一遍“建模 - 二阶锥松弛 - YALMIP 建模 - CPLEX 求解 - 结果校验”的流程。重点不是甩一个跑完就完事的代码包而是讲清楚每个关键写法背后的原因为什么 DistFlow 要那样改锥、为什么标幺化会坑掉一大半新手、为什么求解器报 infeasible 时先别怀疑求解器。适合正在做配电网优化、分布式资源调度、微电网能量管理方向的同学也适合刚接触优化求解、想快速搭建一个可扩展框架的人。1. 动态最优潮流到底难在哪1.1 “动态”两个字让变量直接翻倍静态最优潮流只算一个截面比如某时刻网损最小。动态最优潮流则要把 T 个时段全部串起来典型是 24 个整点或 96 个 15 分钟断面。最直接的冲击就是变量数量膨胀假设一个 33 节点的网络支路 32 条每个时段光支路有功、无功、电流平方、节点电压平方这几个核心变量就接近 200 个24 个时段直接奔着四五千个连续变量去再加入储能 SOC、OLTC 档位这类状态量规模轻松过万。这对求解器来说还不是致命问题真正麻烦的是变量之间有了时间耦合。储能是这种耦合的典型。某时段充电会影响后续所有时段的 SOC 可行域OLTC 一天调节次数有限属于跨时段的计数约束分布式电源出力受气象和日前预测影响本质上也要按时间序列描。方程之间环环相扣你已经不是在求“某一个时刻”的最优解而是在求“一整条可行轨迹”的最优解。这就像排班单独安排一个人某天休息很容易要让一个班组连续一周每天都有人值班且总工时最小就是完全不同的复杂度。如果只算单时段很多配置可以拍脑袋给一旦放到 24 时段任何一段不行整体就不可行。很多初学者第一次跑动态 OPF 时程序报 infeasible其实不是潮流方程写错而是储能首末 SOC 约束太紧、或者负荷高峰时段和光伏高峰时段错配导致整个调度周期根本串不出一条可行轨迹。1.2 潮流约束的非凸性才是求解器的拦路虎维度增加还只是“量变”真正让最优潮流变得棘手的是潮流方程本身的非凸性。在 DistFlow 支路方程里节点电压幅值平方 U、支路电流平方 L、有功和无功 P、Q 之间存在关系L (P² Q²) / U。这个式子把三个变量揉在一起而且除以 U 意味着这是一个非凸区域。非凸优化模型的麻烦在于局部极小值不一定是全局极小值通用求解器完全没有能力判断自己停下来的点到底是不是全局最优。你用内点法直接怼这个非线性问题大概率得到一堆局部最优甚至不收敛的结果而且没人敢说这个结果行不行。于是从业者们想到了二阶锥松弛SOCP Relaxation。思路不复杂先把非凸等式 L (P²Q²)/U 放宽成不等式 L ≥ (P²Q²)/U也就是允许电流平方比物理值更大。这个不等式经过代数变形后能写成标准二阶锥约束。形成凸可行域后再用内点法或分支割平面法求全局最优。这个操作相当于把一个坑坑洼洼的山地变成了一大片规整的下坡路求解器终于能大步往前走了。再叠加上 OLTC 档位、电容器组投切、储能充放电逻辑这些离散变量问题最终会变成混合整数二阶锥规划MISOCP。这也是这个项目标题里为什么点名 CPLEX 的原因——CPLEX 是处理混合整数规划和二阶锥规划都成熟的老牌商用求解器YALMIP 负责建模转换CPLEX 负责在锥约束和整数约束之间来回切效率比纯开源求解器高不少。2. 二阶锥松弛怎么把非凸问题变成凸问题2.1 DistFlow 方程和它的等价变形二阶锥松弛不能套在任意潮流方程上它有一套完整的适用前提——辐射状配电网的支路潮流方程。这套方程在配电网优化里被叫作 DistFlow本质上是从基尔霍夫定律推出来的只保留了有功、无功、电压幅值平方和电流幅值平方四个量。对一个从节点 i 到节点 j 的支路标准形式是三条等式首先是节点功率平衡进入节点 j 的有功功率加上本地发电等于本地负荷和流出功率之和无功也同理。然后是电压降方程末端电压平方 U_j 等于首端电压平方 U_i 减去线路上的压降再补上电流平方引起的修正。最后是电流定义式L_ij (P_ij² Q_ij²) / U_i。前两条等式是线性的问题全部集中在第三条上。P²、Q²、除以 U三者放在一起直接破坏了凸性。松弛操作就是把这第三条等式改成不等式L_ij ≥ (P_ij² Q_ij²) / U_i。把它变个形左右两边同时乘 U_i得到 U_i · L_ij ≥ P_ij² Q_ij²。这是一个旋转二阶锥rotated quadratic cone的标准形式。再经过一次变量替换可以写成更利于代码实现的 2-范数形式|| [2P_ij; 2Q_ij; U_i - L_ij] ||₂ ≤ U_i L_ij。这个形式在 YALMIP 里可以直接照抄成一条约束CPLEX 内部会把它当作原生锥约束处理不会粗暴地展开成一堆二次项。实际写代码的时候我见过不少人手动展开写成平方和不等式不仅费劲而且容易触发 CPLEX 的二次约束非凸检查。2.2 松弛“精不精确”怎么判断和处理有人会问你把等式放宽成不等式那求出来的还是原问题的最优解吗答案要看情况。对于辐射状配电网如果目标函数是网损最小这类“压着线走”的目标松弛后得到的解通常刚好落在锥边界上也就是说 L_ij 和 (P_ij² Q_ij²)/U_i 相等这时候二阶锥松弛是精确的求出的结果就是原非凸问题的全局最优解。但事情没这么绝对。当网里有大量 DG 注入尤其是无功功率边界被压得很紧的时候锥松弛可能会被“撑开”约束没有全激活解出来的 L_ij 明显大于物理上需要的最小电流平方。此时回代校验会发现误差这个情况在专业文献里叫“松弛间隙”duality gap不为零。我自己排查这类问题时最常用的方法是一段十行的回代脚本把求解得到的 U、P、Q 重新代入 L (P²Q²)/U算一个绝对误差超过 1e-4 甚至 1e-3 就要警惕了。如果确实有间隙最常见的修补方法是给目标函数加一个关于 L_ij 的小惩罚项比如把目标改成网损 ε · Σ L让求解器有动力把锥压回边界。ε 通常取 1e-6 到 1e-4具体数值靠你实测去调。加完惩罚之后要重新检查间隙如果还不行就得检查是不是无功边界给得太别扭、或者目标函数的各项权重尺度差太多。2.3 主动配电网的各类设备约束怎么塞进 SOCP把二阶锥当作“平台”之后剩下来的工作就是往这个框架里放设备模型。最基础的是分布式电源光伏和风机的有功出力通常有个上限无功功率容量可以近似成一个圆形或者梯形约束写成线性不等式就够用。微型燃气轮机和燃料电池也类似有功无功都在一个凸多边形里。储能是动态问题里最核心的一段。它的 SOC 递推是线性方程SOC(t1) SOC(t) dt · (η_ch · P_ch(t) - P_dis(t)/η_dis) / E。这里 P_ch 和 P_dis 分别代表充电、放电功率上界是储能额定功率SOC 本身有上下限。如果想让模型更严谨还需要禁止同一时段既充电又放电这就要引入二进制变量模型从 SOCP 变成 MISOCP。OLTC 一类带档位的设备天然是整数变量。YALMIP 里用 intvar 声明CPLEX 负责把离散档位和连续潮流耦合在一起。可投切电容器组也类似用整数变量表示投入组数。可以说模型的最终形态基本就是“线性约束 二阶锥约束 离散整数变量”三件套这也是为什么 YALMIP CPLEX 这套组合在这个场景下这么流行它俩对这三类约束都支持得最完整。3. YALMIPCPLEX完整实现与参数坑3.1 环境准备和测试求解器先把环境装利索。MATLAB 方面我用的版本是 R2023a 以上YALMIP 直接从 GitHub 拉最新版放到 MATLAB 搜索路径里即可。CPLEX 需要装对应平台的版本并且确保 MATLAB 能通过cplex命令访问到求解器。这个环节最常见的坑不是安装而是 MATLAB 和 CPLEX 的位数/版本不匹配导致 yalmiptest 里 CPLEX 那一栏显示红色。所以我每次搭好环境第一件事就是跑一遍yalmiptest。这个函数会把 YALMIP 支持的所有求解器挨个测一遍能识别出来的会显示 OK。但凡 CPLEX 不 OK先别急着改代码回到环境配置上排查。路径配好了、版本对得上了再跑一个极小的二阶锥算例比如二维的 min ||x||₂验证整条链路通不通。这一步花不了两分钟但能省掉后面一晚上的调试时间。3.2 数据准备、标幺化和数据结构数据准备看着简单实际是最容易出错的一环。IEEE 33 节点系统的基准电压是 12.66 kV基准功率常取 10 MVA。线路参数如果原始数据给的是欧姆和千瓦必须先换算成标幺值再进模型。基准阻抗 Z_base V_base² / S_base 12.66² / 10 ≈ 16.02 Ω。很多新手把线路 R、X 零点几个欧姆直接丢进模型算出来的结果偏向离谱根源全在单位混用。我自己的习惯是把负荷、光伏出力、线路参数全部放在一个结构体里按 t 时刻排列。比如 P_load 是 24 × 33 矩阵第 t 行第 j 列表示 t 时段节点 j 的有功负荷P_pv 是 24 × 33 的 PV 最大可用出力矩阵R、X 和 bus_from、bus_to 是支路拓扑数据。这样做的好处是所有约束组装都能严格对齐时间维度和节点维度不会出现维度错位后求解器神秘报错的情况。3.3 变量定义与核心约束代码下面这段代码是整套模型的骨架。我故意省略了太多细节只保留最容易写歪的地方方便你对照自己手里的工程。决策变量都用 sdpvar 声明维度跟着时间、支路、节点走T 24; n_bus 33; n_branch 32; n_storage 2; % 假设节点18和节点25接入储能 % 连续变量 P_br sdpvar(T, n_branch, full); % 支路首端有功 Q_br sdpvar(T, n_branch, full); % 支路首端无功 Lij sdpvar(T, n_branch, full); % 支路电流平方 U sdpvar(T, n_bus, full); % 节点电压平方 % 储能变量 P_ch sdpvar(T, n_storage, full); P_dis sdpvar(T, n_storage, full); SOC sdpvar(T, n_storage, full); % OLTC档位引入整数变量模型变为MISOCP tap intvar(1, T, full);核心约束组装分成三块。第一块是支路方程和二阶锥松弛第二块是节点功率平衡第三块是储能递推。DistFlow 和锥约束的写法如下Constraints []; for t 1:T for k 1:n_branch i bus_from(k); j bus_to(k); % DistFlow 电压降方程 Constraints [Constraints, ... U(t,i) - U(t,j) 2*(R(k)*P_br(t,k) X(k)*Q_br(t,k)) ... - (R(k)^2 X(k)^2)*Lij(t,k)]; % 二阶锥松弛核心中的核心 Constraints [Constraints, ... norm([2*P_br(t,k); 2*Q_br(t,k); U(t,i)-Lij(t,k)]) ... U(t,i) Lij(t,k)]; end end注意这里norm([...]) ...就是标准二阶锥约束。如果你在旧资料里看到cone函数那是 YALMIP 的另一套写法本质一样。我更喜欢 norm 写法的原因很简单一眼能看懂它在表达什么调试时不容易把维度搞错。节点功率平衡部分先按拓扑判断每个节点的父支路和子支路然后列平衡方程。下面是只含负荷和储能的示意光伏或 DG 再加对应的注入项即可for t 1:T for j 1:n_bus inflow 0; outflow 0; for k 1:n_branch if bus_to(k) j % 支路末端注入节点j的功率 首端功率 - 线路损耗 inflow inflow P_br(t,k) - R(k)*Lij(t,k); end if bus_from(k) j outflow outflow P_br(t,k); end end % 本地注入项储能充电是负荷放电是发电 P_gen 0; for s 1:n_storage if storage_bus(s) j P_gen P_gen P_dis(t,s) - P_ch(t,s); end end Constraints [Constraints, ... inflow - outflow P_load(t,j) - P_gen - P_pv(t,j)]; end end无功功率平衡的写法一模一样只需要把有功换成无功、把 R(k)*Lij 换成 X(k)*Lij再补上无功负荷、光伏无功上限、储能无功为 0 这几个条件。我说“一模一样”是真的一模一样复制粘贴再改变量名就行。储能递推约束注意时刻对齐和效率系数E 1.0; % 储能容量MWh eta_ch 0.95; % 充电效率 eta_dis 0.95; % 放电效率 Pmax 0.2; % 储能额定功率MW for s 1:n_storage Constraints [Constraints, SOC(1,s) 0.5]; for t 2:T Constraints [Constraints, ... SOC(t,s) SOC(t-1,s) ... (eta_ch*P_ch(t-1,s) - P_dis(t-1,s)/eta_dis) / E]; end Constraints [Constraints, SOC(:,s) 1, SOC(:,s) 0.2]; Constraints [Constraints, SOC(T,s) 0.5]; % 充放电功率限幅 Constraints [Constraints, 0 P_ch(:,s) Pmax, 0 P_dis(:,s) Pmax]; % 若加二进制变量阻断同时充放电 % u binvar(T,1,full); % Constraints [Constraints, P_ch(:,s) Pmax*u, P_dis(:,s) Pmax*(1-u)]; end这段代码里我注释掉了同时充放电的二进制变量。实际建模时如果项目要求必须符合物理时序建议把那段注释取消如果只是先验证算法框架SOCP 连续松弛跑得快但解里很可能出现同一时段又充又放的理论解回代工程现场会被质疑。这个取舍要看你的场景是学术计算还是工程可实施。3.4 求解器参数设置和常见状态码模型组装完目标函数和求解器设置直接决定最终能不能跑出好结果。以网损最小为例目标函数很简单Objective sum(sum(Lij .* R)); % 各支路电流平方乘电阻再累加所有时段YALMIP 里点乘是逐元素乘Lij 和 R 的维度需要匹配。R 如果只有 1×32要用 repmat 或广播展开成 24×32。这里不展开但代码里容易因为维度问题报错。求解器设置我一般写成options sdpsettings(solver, cplex, verbose, 2); options.cplex.mip.tolerances.mipgap 0.01; options.cplex.optimalitytarget 3; options.cplex.timelimit 600; sol optimize(Constraints, Objective, options);mipgap 0.01的意思是最优性间隙允许 1%对工程场景通常够了能大幅缩短求解时间。optimalitytarget 3这条设置很关键尤其是在旧版 CPLEX 里遇到锥约束和整数变量混合的模型不设置这个值可能触发“quadratic constraint non-convex”的误报。timelimit防止模型卡死超过十分钟直接停方便你回头优化模型而不是干等。求解结束之后一定看一眼sol.problemif sol.problem 0 disp(求解成功); else disp(sol.info); % 返回具体错误信息 end0 表示成功非 0 就要结合sol.info排查。常见的问题码包括不可行、无界、数值问题等。我很少背数字表直接把sol.info打出来YALMIP 会给出人话解释。3.5 结果提取与合理性检查求解完别忘了把变量还原成物理量否则一堆标幺值账根本对不上。电压幅值就是sqrt(value(U))电流平方直接value(Lij)储能 SOC 曲线value(SOC)。这些量建议一次性存成结构体方便后续画图。看结果的时候我有一套固定动作先画 24 时段全网电压曲线看所有节点电压是否落在 0.95~1.05 之间再画储能 SOC 曲线看首末值是否满足约束最后检查支路电流是否越限。如果前两步没问题但最后电流越了限说明支路容量约束没加或者加错了。另外要习惯用“状态信息”而不是“目标函数值”来判断求解质量。目标函数 0.05 和 0.052 看起来差不多但如果求解状态在 MIP 节点上挣扎半天、gap 一直降不下去你就要考虑是不是整数变量和连续变量的量级差太多导致数值病态。此时通常是调单位或调参数而不是盲目加约束。4. 常见问题与排查技巧实录4.1 求解器返回 infeasible问题基本出在数据第一次跑动态 OPF十个人有八个会碰到 infeasible。这时候先别怀疑二阶锥写错因为我见过的案例里九成以上都是数据问题一个是单位没换算。线路阻抗用欧姆、负荷用千瓦但基准容量取的 MVA电流平方和电压平方的公式全乱套。检查方法很简单取一个纯负荷且无 DG 的时段手算一下根节点注入功率和平衡节点的电压至少保证量纲一致。一个是拓扑方向搞反。DistFlow 是沿“根节点 - 叶子节点”的方向写的如果某条支路的 bus_from 和 bus_to 反了电压降方程符号错功率平衡也跟着错连锁反应就是整个约束系统自相矛盾。还有一个是储能约束太紧。初始 SOC 明明只有 0.3又要求 24 小时后回到 0.5中间还有一大段负荷高峰根本没有充电窗口整段时间轨迹就断了。解决办法是先放宽 SOC 上下限、去掉末状态约束等模型能解之后再逐步收紧逐步逼近真实要求。排 infeasible 最有效的手段是“二分法删约束”。把储能、OLTC、DG 一段段注释掉直到模型能解再放回来缩小问题范围。这个手段比盯着报错日志大海捞针高效得多。我自己的习惯是先把动态问题退化成单时段静态问题单时段能解再慢慢加时间耦合项加一个检查一次基本能定位到是哪条约束把模型锁死了。4.2 锥松弛间隙大电压曲线失真怎么办有时候模型能解目标函数也收敛了但电压曲线在重负荷时段出现反常的剧烈凹陷或者支路电流平方和 P²Q²/U 对不上。这是二阶锥松弛间隙大的典型征兆。我自己的检查脚本长这样P_ value(P_br); Q_ value(Q_br); U_ value(U); L_ value(Lij); gap max(max(abs(L_ - (P_.^2 Q_.^2) ./ U_))); disp([最大锥松弛间隙 , num2str(gap)]);间隙如果超过 1e-4就要处理。最常见的修复方案是往目标函数里加一个关于 Lij 的小系数惩罚项把锥“拽”回边界。不过这个惩罚项的系数不能太大太大相当于把总损耗目标都压掉了解出来的调度方案对网损不再敏感。我从 1e-6 开始试逐步放大每次重新求解并检查 gap直到 gap 小于阈值。另一个有效技巧是给电压变量 U 提供一个高质量的初始值。先用普通潮流工具算出 24 个时段的基准电压在 YALMIP 求解前用assign赋予初值再让 CPLEX 从这组初值开始求解。对锥松弛类问题热启动效果很明显求解器不需要在寻找初始可行解上花太多功夫锥约束也更容易切到边界。4.3 求解速度慢多半不是求解器的锅如果你发现 33 节点的 24 时段模型跑了五分钟以上还不出结果先不要觉得是 CPLEX 不行大概率是模型本身太“胖”。最典型的原因是把约束一股脑写进一个巨大的循环里YALMIP 需要为每条约束生成内部表达式几万条约束堆下来建模时间占比比求解时间还大。优化方向有几个。第一把循环里的重复表达式抽出来比如支路阻抗平方 R²X² 提前算好。第二尽量避免在约束里写sum(sum(...))这种维度广播过大的表达式改成矩阵运算代替循环。第三合理使用 YALMIP 的批量约束写法比如Constraints [Constraints, U 1.05^2, U 0.95^2]这种整矩阵约束比双层循环逐一 append 快得多。如果模型里带了 OLTC 和储能充放电二进制变量那问题就变成 MISOCP分支定界过程本身就可能比较慢。这时可以合理放宽 mipgap。我在 33 节点算例上把 mipgap 从 1e-4 放宽到 1e-3速度提升好几倍目标函数只多了不到 0.2%工程上完全可接受。4.4 用回代测试给结果“验真”最后一步也是我最坚持的一步用二阶锥解算出来的调度方案重新跑一遍精确的交流潮流回代验证电压和支路潮流是否真的满足物理方程。这不是多此一举因为“二阶锥松弛精确”是一种理想情况实际模型里加了很多设备近似可能继承了一点误差。回代的方式有现成工具MATPOWER 的runpf就能给一组 PQ 和 PV 节点算精确潮流。你也可以自己写一个简单的牛顿法去核对关键节点电压工作量不大。回代后如果最大电压偏差超过 0.005 p.u.我就要回头检查锥松弛间隙和各个设备约束的近似程度。这步做扎实了后续论文里的解、工程可实施性都有底。我曾经踩着个坑某次程序跑完目标函数很好看电压也在限值内但回代后某些节点电压直接掉到 0.91。查了一天发现问题出在光伏无功约束我用了线性近似而实际逆变器的无功能力曲线是圆的模型松弛掉了一块导致最优解刚好落在被松弛掉的那块区域。从那以后回代校验成了我所有优化代码出结果后的固定动作。在我个人实际跑了大量配电网优化算例之后有一个体会特别深这套 YALMIP CPLEX 二阶锥松弛的组合最大价值不是“能算”而是可扩展性强。从 33 节点换到 123 节点从单目标网损改成多目标加权电压偏差从确定性负荷换成随机场景基本只需要改数据、改约束、改目标函数不用推翻整体架构。代码里的求解参数、数据单位、校验脚本最好从一开始就用结构体和脚本固定下来后续折腾新模型时能省掉大量返工。最后再分享一个小技巧如果你发现二阶锥松弛在某个边界算例上老是翻车不妨把锥约束从norm([2P;2Q;U-L]) UL这种显式写法换成旋转锥写法再用 YALMIP 的cone函数有时候不同写法会让内部预处理的路径完全不同莫名就稳了。别每次都只背一种写法多试几种你对自己模型的理解会深很多。希望这篇内容能帮你在动态最优潮流这条路上少走几周弯路。