1. 先搞清楚这道题到底在解什么做配电网多时段优化的人几乎都绕不开“YALMIP CPLEX 二阶锥松弛”这一套组合。我在IEEE33和PG69两个经典算例上把这套流程完整落过地最终算的是多时间断面的配电网最优潮流既包含时序负荷、储能的充放电调度也考虑了电压约束和网损优化。如果你是刚接触这方面的人我建议不要一上来就啃数学论文先把“这道题在解决什么”想明白后面建模会顺利很多。1.1 为什么偏偏是IEEE33和PG69IEEE33节点系统是配电网文献里出现频率最高的算例之一33个节点、32条支路、额定电压12.66kV总负荷大约3.715MW加2.3Mvar。规模不上不下既不会让初学的人看到一大堆节点发懵又能体现电压损耗、线路重载、末端压降这些真实配电网问题。更重要的是它的数据在网上很容易找到几乎所有配电网重构、分布式电源接入、储能调度的文章都用过它你拿自己的模型去跑一遍跟论文里的结果一对就知道有没有做错。PG69这个系统更刺激一点69个节点比IEEE33多了一倍多总负荷大概3.8MW加2.7Mvar。它的干线更长越靠近末端电压下降越明显所以对“潮流计算电压约束”来说是个更严格的测试场景。很多人在IEEE33上跑得通换到PG69就出现收敛问题或者锥松弛不精确恰恰就是因为系统规模变大、末端节点电压过低让模型对约束的真实压力显现出来。我把两个系统的特点拉成一张表方便对比算例节点数支路数额定电压典型总负荷特点IEEE33333212.66kV3.715MW 2.3Mvar规模适中资料丰富适合入门验证PG69696812.66kV3.8MW 2.7Mvar干线路长末端压降明显适合压力测试1.2 “多时间断面”到底多在哪单时段潮流优化其实很简单就是给定某一时刻的负荷和DG出力算一个最优运行点。但现实里光伏出力从早到晚是条曲线负荷也是波动的储能更是晚上充电、白天放电这天然是一个跨时段耦合的问题。所谓多时间断面就是你把一天按小时切成24个断面甚至按15分钟一个断面切成96个断面所有断面共享同一个网络拓扑每个断面有自己的负荷、新能源出力和电压分布但断面与断面之间又通过储能SOC、机组爬坡、需求响应等约束串在一起。普通的潮流计算不用管这些但最优潮流如果忽略时间耦合它给出的储能策略一定是错的。比如储能单时段模型里它只需要满足“这一小时充的电等于这一小时放的电”这明显不符合实际。多时段模型里你要约束每两个相邻时段的电量递推关系还要限制充放电状态不能同时为1这些都是时间维带来的麻烦。这也是我在IEEE33和PG69上做多时段建模时最深的体会空间维度靠DistFlow方程时间维度靠储能和爬坡两者加在一起问题才完整。1.3 这套工具组合的逻辑YALMIP守前场CPLEX收尾工具选型很多人问过为什么不用Gurobi不用Mosek偏偏用CPLEX我的理由很简单首先很多高校和研究所本来就有CPLEX的学术许可证装起来不费劲其次CPLEX对混合整数二阶锥规划也就是MISOCP支持非常成熟第三YALMIP作为MATLAB里的建模层让你用近乎数学公式的方式写出锥约束和整数变量避免了手动把问题转成求解器接口的麻烦。流程是这样的YALMIP负责把二阶锥约束、储能SOC约束、整数变量这些建模逻辑转化成CPLEX能识别的标准形式CPLEX内部用分支切割和内点法把MISOCP解出来。整个过程你不需要手写任何内点法代码也基本不用关心算法的迭代细节。这也是工程上最舒服的地方建模思路集中在“问题长什么样”而不是“求解器内部怎么做”。2. 二阶锥松弛建模中最关键的一跳2.1 DistFlow方程的非线性卡点配电网潮流计算如果从零开始写最经典的做法是用DistFlow方程。以一条支路ij为例从节点i流向节点j的有功功率记作P_ij无功功率记作Q_ij节点电压幅值的平方记作v_i支路电流幅值的平方记作l_ij。DistFlow的核心方程可以写成这几条P_ij p_j^load - p_j^dg Σ P_jk r_ij * l_ijQ_ij q_j^load - q_j^dg Σ Q_jk x_ij * l_ijv_j v_i - 2 (r_ij P_ij x_ij Q_ij) (r_ij² x_ij²) l_ijl_ij (P_ij² Q_ij²) / v_i前三条都是线性约束问题就出在最后一条。l_ij等于一个分式分母是v_i分子是P_ij和Q_ij的平方和这是一个非凸等式。你去做优化的时候这个等式会让整个可行域变得极度扭曲普通的线性规划和二次规划都没法直接处理。很多初学的人会问那直接用牛顿法跑潮流不就行了吗可以但你要搞清楚潮流计算和最优潮流是两回事。潮流计算是给定运行点求解状态量你需要的是唯一可行解而最优潮流是要在约束里找使目标最小化的那个点如果模型非凸你得到的很可能只是局部最优甚至根本不收敛。要做全局最优的分析就得想办法把这团非凸的东西“掰”成凸问题。2.2 从非凸等式到旋转锥约束二阶锥松弛的思路很直接。原来那个等式约束太苛刻了我把它放宽成一个不等式l_ij ≥ (P_ij² Q_ij²) / v_i这个不等式的意思是支路电流的平方至少要达到欧姆定律要求的最小值。因为目标函数里通常会包含网损项求解器为了压低损耗会尽量让l_ij变小最后它会自动被压到这个不等式边界上也就是说最优解仍然满足等式。于是我们把“直接求等式”变成了“求不等式的最优解”。这个不等式还带一个分式直接交给CPLEX也不方便所以要再转写成标准的二阶锥形式。经过数学变形可以得到下面这个约束|| [2P_ij; 2Q_ij; v_i - l_ij] ||₂ ≤ v_i l_ij这看起来唬人其实很简单。把两边同时平方展开左边是4P_ij²加4Q_ij²再加(v_i - l_ij)²右边是(v_i l_ij)²化简之后得到的就是4P_ij² 4Q_ij² ≤ 4 v_i l_ij等价于v_i l_ij ≥ P_ij² Q_ij²正好是我们想要的松弛不等式。在YALMIP里写这个约束不需要你自己去展开矩阵直接用cone函数就行Cons [Cons, cone([2*Pij(b,t); 2*Qij(b,t); V2(i,t)-Iij(b,t)], V2(i,t)Iij(b,t))];cone函数第一个参数是一个向量第二个参数是标量r它建模的数学含义就是||向量||₂ ≤ r。YALMIP会自动判断这是一个二阶锥约束并在传给CPLEX之前把它转换成标准SOCP格式。这一步是整套模型里最关键也最容易被忽略的地方。2.3 什么时候二阶锥松弛是精确的松弛之后可行域变大了不等于最优解一定落在原可行域上。如果你做完优化发现l_ij明显大于(P_ij²Q_ij²)/v_i说明锥松弛“不紧”解出来的结果并不对应一个真实物理潮流。根据理论和实际经验在IEEE33和PG69这类辐射状配电网里只要满足两个条件锥松弛基本是精确的第一每条支路的电阻和电抗都大于0第二目标函数对支路电流平方l_ij是严格单调递增或者至少不是反向激励。典型的目标比如最小化网损、最小化购电成本这些目标都会促使求解器把电流压到最低最后锥约束自然收紧。实际操作中我推荐在目标里显式加入网损项哪怕你的核心目标是调度储能或削减峰值也放一个很小的网损惩罚系数。这样既不会明显扭曲经济性目标又能保证锥松弛的精确性。跑完之后检查一下锥松弛间隙也是一个必要的自检步骤。2.4 加入离散动作后问题升级为MISOCP多时段时间断面如果只做连续变量纯SOCP就够了。但只要你想考虑有载调压变压器分接头、电容器组投切、储能充放电状态问题就变成混合整数二阶锥规划这些离散量必须用整数变量或者0-1变量建模。CPLEX处理MISOCP的能力相当强这也是我坚持用它的原因之一。在YALMIP里定义一个储能充放电状态变量只需要写u binvar(ns, T); % ns个储能T个时段后面再配合两个不等式保证充电功率大于0时放电功率为0Pch Pmax * uPdis Pmax * (1 - u)如果没有这个0-1变量模型很可能出现同一时段又充电又放电的荒谬结果。这个细节虽然简单但很多人第一次建模时都会漏掉导致结果完全没法看。3. 多时间断面的完整数学建模和变量组织3.1 时间耦合约束到底从哪来多时段问题里最容易写错的就是时间耦合约束其中最典型的是储能SOC递推关系。以一个时间段间隔为Δt的模型为例假设Δt以小时为单位储能荷电状态E_t的计算公式是E_{t1} E_t η_ch * Pch_t * Δt - Pdis_t * Δt / η_dis其中η_ch是充电效率η_dis是放电效率。这个递推式把相邻两个时段绑在一起不能独立求解。与此同时还要让充放电功率不同时为正所以引入了前面说的0-1变量u_t。除了储能分布式电源出力也有爬坡速率限制-ΔP_ramp ≤ Pg_{t1} - Pg_t ≤ ΔP_ramp这个约束本质上也是时间耦合只不过比储能递推稍微温和一点它不限制电量只限制每一段之间的变化幅度。如果模型还包含有载调压变压器、电容器组那么这些设备的分接头档位和投切动作也应该跨时段限制。最简单的做法是两个相邻时段最多动作一次或者限制一天内动作总次数。这些约束不加求解器就会给出每15分钟疯狂切换一次的理想方案工程上根本没法执行。3.2 每个时段都要满足的DistFlow约束时间维度增加了但每个时段内部仍然要满足配电网潮流方程。对于每个断面t我要对每一条支路b写一组DistFlow约束。先给支路定义from和to两个方向数组from br(:, 1); % 支路起点编号 to br(:, 2); % 支路终点编号 R_ohm br(:, 3); % 电阻单位欧姆 X_ohm br(:, 4); % 电抗单位欧姆然后利用前面说过的SOCP锥约束对每个时间断面建立支路电压降和锥约束。代码如下% 电压降线性约束 Cons [Cons, V2(to(b), t) V2(from(b), t) ... - 2*(R_pu(b)*Pij(b,t) X_pu(b)*Qij(b,t)) ... (R_pu(b)^2 X_pu(b)^2)*Iij(b,t)]; % 二阶锥约束 Cons [Cons, cone([2*Pij(b,t); 2*Qij(b,t); ... V2(from(b),t) - Iij(b,t)], V2(from(b),t) Iij(b,t))];节点功率平衡也要逐时段写。对于每个节点i所有流出支路功率减去所有流入支路功率再加上线路损耗必须等于该节点净注入功率。这里的净注入等于DG出力减去负荷公式为Σ Pij_{out} - Σ Pij_{in} Σ r_b * l_ij_b Pg_i - Pload_i实际写代码的时候用find函数找每个节点的出线和进线支路索引就行。3.3 目标函数怎么设计才算合理目标函数是整个模型战略性的部分。我做IEEE33和PG69多时段算例时最常使用的是三部分相加购电成本、网损、储能成本。购电成本从上级电网买电的费用通常用分时电价乘以根节点注入功率网损所有支路电流平方乘以支路电阻再求和这个也是让SOCP约束保持紧凑的关键项储能成本可以理解成电池循环损耗也可以加一个对充电次数的软约束。YALMIP里目标函数可以直接写成Objective sum(price .* Pg) ... sum(sum(R_pu .* Iij)) * baseMVA ... 0.001 * sum(sum(Pch Pdis));有人会问为什么网损要用p.u.值再乘baseMVA因为配电网里潮流功率的真实量级是MW而线路电阻在p.u.下通常只有0.005左右两者乘出来的网损可能在1e-4这种量级直接作为目标会被求解器当成噪声忽略反而丢失了锥松弛的收紧作用。乘回baseMVA后网损就恢复成几十千瓦到几百千瓦的真实量级和购电成本保持在同一个数量级CPLEX在数值上才会认真优化它。3.4 变量的维度设计是少走弯路的重点多时段模型的变量维度设计我建议一开始就按“支路数×时段数”或“节点数×时段数”来建矩阵而不是每个时段单独设一套变量。这样YALMIP内部矩阵维数小求解速度也快。以IEEE33为例T取24那么Pij sdpvar(nb, T); % 支路有功nb行T列 Qij sdpvar(nb, T); % 支路无功 Iij sdpvar(nb, T); % 支路电流平方 V2 sdpvar(n, T); % 节点电压平方 Pg sdpvar(ng, T); % 上级电网注入有功 Qg sdpvar(ng, T); % 上级电网注入无功后面在写约束时只要把t从1到T循环一遍用Pij(b,t)索引某一个具体时间段即可。这种二维变量矩阵既直观又能利用YALMIP对结构化变量的优化比把所有变量拉成一维长向量再慢慢拼约束要舒服得多。4. CPLEXYALMIP实战落地从环境配置到代码骨架4.1 安装环节最容易翻车先讲清楚很多人的模型本身没有错但卡在第一步CPLEX装好了YALMIP却找不到求解器。YALMIP不是自带求解器的它只是一个建模层它把所有约束翻译成标准模型文件后需要调用外部求解器来算。所以必须保证CPLEX的MATLAB接口路径被正确添加。CPLEX安装好后注意找到它的MATLAB接口目录一般是安装路径下的cplex/matlab文件夹。在MATLAB里执行addpath(genpath(D:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab)); savepathYALMIP的安装也是类似把下载下来的文件夹整个加入路径addpath(genpath(D:\yalmip)); savepath装完之后运行yalmiptest你会看到YALMIP检查所有已安装求解器的报告。里面会有一行显示CPLEX相关的测试是否通过。如果显示missing或者error基本就是路径没配对。关于CPLEX获取方式高校用户建议走学术计划申请教育版证书个人使用也可以看看官方社区里免费的社区版不过社区版对问题规模有限制IEEE33的24时段问题勉强能顶PG69的多时段MISOCP很可能超限长期做研究还是用完整许可证靠谱。4.2 数据准备和标幺化IEEE33和PG69的原始数据里面线路参数通常是以欧姆为单位给出的负荷以kW为单位而优化求解器在处理SOCP这种含有二次约束的问题时数值范围太离谱会导致收敛困难。我强烈建议计算前先把所有量统一到标幺值下。一个实用的标幺选择是基准电压Vb取12.66kV基准功率Sb取10MVA。这样IEEE33的总负荷3.715MW就变成0.3715p.u.数值合理。线路阻抗的转换公式是Z_b Vb² / Sbr_pu R_ohm / Z_bx_pu X_ohm / Z_b打个具体比方IEEE33第一条支路的电阻是0.0922ΩVb12.66kVSb10MVA时Z_b等于16.02Ω换算出来的标幺电阻大约0.00576。这个数值在二阶锥约束中不会过小目标函数中的网损项也不会被淹没。如果不做这一步直接用欧姆值和kW去建模数学上虽可行但数值条件数会很差CPLEX求解效率和精度都会下降。4.3 一个可直接套用的求解代码骨架我习惯把建模代码拆成几个逻辑块数据定义、变量声明、约束循环、目标函数、求解、结果提取。下面是一个基于IEEE33、24时段、含储能的多时段SOCP最简骨架%% 1. 基础数据 T 24; n 33; nb 32; baseMVA 10; % branch数据假设已经整理成from, to, R_ohm, X_ohm Zb 12.66^2 / baseMVA; R_pu R_ohm / Zb; X_pu X_ohm / Zb; %% 2. 定义变量 Pij sdpvar(nb, T); Qij sdpvar(nb, T); Iij sdpvar(nb, T); V2 sdpvar(n, T); Pg sdpvar(1, T); % 根节点注入有功 Qg sdpvar(1, T); %% 3. 储能变量 ns 2; Pch sdpvar(ns, T); Pdis sdpvar(ns, T); E sdpvar(ns, T); u binvar(ns, T); %% 4. 约束 Cons []; for t 1:T % 根节点电压设为1.0 Cons [Cons, V2(1,t) 1.0]; % 支路DistFlow与锥约束 for b 1:nb i from(b); j to(b); Cons [Cons, V2(j,t) V2(i,t) ... - 2*(R_pu(b)*Pij(b,t) X_pu(b)*Qij(b,t)) ... (R_pu(b)^2 X_pu(b)^2)*Iij(b,t)]; Cons [Cons, cone([2*Pij(b,t); 2*Qij(b,t); V2(i,t)-Iij(b,t)], ... V2(i,t)Iij(b,t))]; end % 节点有功平衡 for i 1:n outb find(from i); inb find(to i); Cons [Cons, sum(Pij(outb,t)) - sum(Pij(inb,t)) ... sum(R_pu(outb).*Iij(outb,t)) Pg(1,t) - Pload(i,t)]; end % 节点无功平衡格式类似这里省略 end %% 5. 储能时间耦合 for s 1:ns for t 1:T-1 Cons [Cons, E(s,t1) E(s,t) 0.9*Pch(s,t) - Pdis(s,t)/0.9]; Cons [Cons, Pch(s,t) 0.5*u(s,t) ... , Pdis(s,t) 0.5*(1-u(s,t))]; end end %% 6. 目标函数 Objective sum(Pg(1,:) .* price) ... sum(sum(R_pu .* Iij)) * baseMVA ... 0.001*sum(sum(Pch Pdis)); %% 7. 求解 ops sdpsettings(solver,cplex,verbose,2,showprogress,1); ops.cplex.mip.tolerances.mipgap 1e-4; optimize(Cons, Objective, ops);这段代码的核心思路就是先把网络约束逐时段写入再补时间耦合约束。实际运行时要根据你自己数据里的Pload(i,t)和price(t)去完善我这里为了让骨架更干净省略了无功平衡的对应写法你做的时候别漏。4.4 求解器参数千万不要全默认很多人习惯直接optimize(Cons, Objective)在简单小问题上没问题但多时段的MISOCP尤其PG69这种上百个支路的算例不做参数调整会让求解时间从几十秒变成几小时。我常用的CPLEX关键参数有三个ops.cplex.mip.tolerances.mipgap相对MIP间隙。默认可能是1e-4甚至更小对于工程分析我一般放1e-4。如果你只需要评估策略趋势放到1e-3能大幅度加速。ops.cplex.timelimit求解时间上限。实际工程中设一个时间限制比如600秒超时后直接取当前最好可行解避免无休止地卡在分支定界上。ops.cplex.threads并行线程数。CPLEX默认会用满所有核心但有时候并行反而导致内存爆炸手动调到物理核心数会稳一点。另外verbose开到2可以在命令行看到实时间隙变化对判断求解是否卡住很有帮助。我就是靠这个判断是加约束还是调参数不用等半天最后看一个冷冰冰的infeasible。5. 我踩过的坑和排查手册5.1 求解器报“infeasible”怎么办新模型的第一个报错大多是不可行。我见过太多人直接开始改动约束但我习惯反过来做减法检查。第一步把储能时间耦合约束全部注释掉只算每个时段独立的潮流如果这时候还是不可行说明问题出在网络约束本身跟时间维度无关。第二步检查负荷是不是有个别节点没接上导致功率不平衡第三步检查电压上下限约束是不是设得太紧比如把0.93到1.07改成0.95到1.05可能某个时段本身就解不出来。还有一种隐蔽情况我一开始用SB1MVA做标幺IEEE33的负荷达到3.7p.u.储能的充电功率0.5p.u.变成了0.5MW看似合理但电压降方程的数值尺度还是有点别扭。后来把SB改成10MVA问题立刻顺了很多。所以infeasible排查的第一步永远是回头看标幺基准别急着怀疑模型逻辑。5.2 怎么验证锥松弛到底紧不紧SOCP结果出来之后不要只盯着目标函数值就完事。我会单独算一下每一条支路的锥松弛间隙powerij value(Pij); powerqj value(Qij); current value(Iij); volt2 value(V2); % 对每条支路的每个时段计算松弛间隙 for t 1:T for b 1:nb v_i volt2(from(b), t); gap(b,t) v_i * current(b,t) - powerij(b,t)^2 - powerqj(b,t)^2; end end max_gap max(gap(:));如果max_gap在1e-6这个量级甚至更低说明锥松弛是紧的解出来的结果可以直接当真实潮流用。如果gap跑到1e-3或者更大说明松弛不紧你得回去检查目标函数里是不是漏了网损项。这个方法花不了几行代码但能让结果可信度提高一大截我在IEEE33上跑多时段时一般gap都能压到1e-6以下。5.3 求解时间太长怎么办第一次面对PG69的96时段、带储能的模型时我差点怀疑电脑坏了连续跑了两个小时都没出结果。后来发现问题是三个叠加整数变量太多、MIP gap设太小、verbose都没开看不出来进展。解决手段有三个方向。第一个方向是变量瘦身储能充放电状态虽然有0-1变量但如果允许储能长时间不动作就不要给每个存储电池都加状态变量合理简化能砍掉一半整数变量。第二个方向是目标函数里的惩罚系数不要设置太高否则会造成数值上的瓶颈求解器会在同一批次上反复切割。第三个方向是时间步长策略先跑T6验证模型再逐步放大T12、T24这样能快速定位是哪类约束导致复杂度暴涨。等T24模型能在几十秒内跑完再考虑加更多时段。5.4 数值警告满天飞的排查思路CPLEX有时候会给出类似Overflow或者Numerical difficulties的警告。我在多时段建模中遇到过的常见诱因是约束两边数值差距过大。比如线路电阻标幺值是1e-5而功率变量是几十乘完之后约束右侧的量级跨度介于几个数量级之间自然容易出问题。解决思路还是回到标幺化和约束重缩放上。你可以把基准功率调大比如1MVA换成10MVA或者把目标函数中的权重缩放一下让所有约束等式右侧都在0.001到1000之间数值困难就很少出现了。6. 最后分享一个我坚持至今的调试习惯这套多时段二阶锥松弛模型我前前后后写过不下十个版本每一次改动完后我都会先固定T2或者T3只测极短时间断面。这样做的好处显而易见两个小时跑不出来的问题压缩成两个时段可能十秒就能出结果我可以迅速验证约束是否写错、锥约束是否被正确识别、数值是否有警告。短时段的解虽然经济性意义不大但它的结构足以暴露建模错误比一上来就跑96时段然后对着日志发呆要高效得多。另外一个习惯是结果可视化。多时段优化跑完把电压分布画成二维热力图或者各时段曲线很多时候一眼就能看出问题。比如某条支路电流阶段处出现明显尖峰很可能就是约束里哪个节点索引对不上。数值上看着正常的解画出来未必正常这个经验我在IEEE33和PG69上都验证过无数次。你如果也在这两个算例上做多时段建模我建议早点养成这两个习惯能帮你少熬好几个通宵。