配电网重构与多时间断面潮流联合优化：从SOCP建模到YALMIP/CPLEX求解

配电网重构和多时间断面潮流放在一起基本就是拿着IEEE33和PG69这类经典算例去解决“含分布式电源的配电网日前运行优化”问题。核心诉求很直接负荷和光伏一天都在波动分段开关和联络开关怎么调才能让全天网损最小、电压不越限这类模型我习惯用二阶锥松弛把非凸潮流方程磨平成凸约束再用YALMIP建模、CPLEX求解整套流程跑下来就是标题里的“多时间断面潮流与配电网重构”联合优化。这篇东西适合正在做配电网重构、DG接入、储能日前调度方向的研究生也适合配网运行方式优化的工程师。下面从模型选型、算例数据、数学约束、YALMIP代码、多时段耦合、求解器调试一路讲过来重点写那些论文里不会写、但你实际跑模型一定会踩的坑。1. 问题本质与方案选型为什么是SOCP、YALMIP和CPLEX1.1 重构问题的本质在拓扑和潮流之间找平衡要理解配电网重构先想一个画面一条馈线一开始按辐射状运行联络开关全部断开。某个时刻光伏大发靠近馈线末端的电压被顶得很高而另一条支路负荷很重、线损很高。如果这时候能合上某个联络开关、断开某个分段开关让部分负荷由另一条路径供电电压和网损都可能明显改善。重构要做的事情就是从所有满足“辐射状连通性”的拓扑里选一个最优的拓扑。这里有两个硬性要求很容易在建模时忽略。一是任何时候都不能形成环网配电网绝大多数时间是开环运行闭环会导致保护装置不配合规划上不能接受二是不能甩负荷重构后所有节点必须还在同一棵树里。这两个要求落到模型上就是第3节要细说的辐射状约束。我自己初学的时候就因为在MATLAB里只加了“闭合支路数节点数-1”这个约束结果跑出来的拓扑又是闭环又是孤岛整个结果根本不能用。1.2 多时间断面的真需求时序耦合不能省单断面重构在文献里早就很成熟了为什么非要多时间断面因为光伏、负荷、储能全是时变的。如果只挑一个负荷最重的断面做重构很可能得到一个“中午光伏大发时最优、晚上负荷高峰时很糟糕”的拓扑。多时间断面把24小时的潮流串起来一起优化开关动作次数也作为约束或惩罚项进入模型这才和现场的智能开关操作次数限制对得上。加入储能之后多断面的必要性就更强了。储能的SOC是一个跨时段的连续状态变量单断面模型根本没法表达“白天多充、晚上多放”这种耦合关系。所以说多时间断面潮流不是简单的“把单断面算24遍”而是要让不同时段的潮流通过储能SOC、开关状态、切换次数真正耦合起来。这也是这类MISOCP模型比普通单断面重构难解的主要原因。2. 算例系统怎么选IEEE33和PG69的差别与数据准备2.1 系统参数对比与选择思路IEEE33和PG69是配电网重构和分布式电源研究里用得最多的两个公共算例。IEEE33系统是33节点、32条支路、5个联络开关基准电压12.66kV总负荷大约3715kW加2300kvar结构相对简单分支少适合用来验证算法正确性。PG69系统是69节点、68条支路、7个联络开关基准电压同样是12.66kV总负荷约3802kW加2694kvar馈线更长、分支更多更贴近真实配电网形态适合做规模测试和分布式电源接入研究。选哪个取决于你的目的是什么。如果是验证一个新想法、调试模型逻辑一定先用IEEE33跑通。等模型逻辑、数据索引、代码风格都稳定了再换PG69跑规模。两个系统的数据格式其实是通用的只要把bus数据和branch数据按同样的表格结构整理好同一套YALMIP代码直接改参数就能跑。我见过很多人一上来就用PG69调模型结果一个小时代码逻辑错误根本定位不了非常浪费时间。2.2 数据整理与时序曲线处理两个系统的数据在公开论文里都有我一般习惯整理成两个表bus_data包含节点编号、有功负荷、无功负荷branch_data包含首端节点、末端节点、电阻、电抗、初始开关状态。读取用readtable或者xlsread都行关键是索引编号要一致。MATLAB里所有数组从1开始而IEEE33原始数据里节点编号也是1到33天然对得上不需要额外偏移。PG69同样。时序数据方面负荷曲线我通常取典型日的24点负荷系数比如凌晨低谷0.5、晚高峰1.0这种然后让每个节点的负荷等于基准负荷乘该时刻系数。光伏出力曲线用归一化的典型日曲线中午1.0、晚上0。注意一个细节如果所有节点用同一个全局系数场景是简化版本更接近实际的做法是不同节点用不同的负荷形态但作为算法验证全局系数完全够了。功率基准值我强烈建议统一改成标幺值。比如基准容量取10MVA基准电压取12.66kV那么所有功率都在p.u.量级电阻电抗的p.u.值在0.001到0.1之间二阶锥约束右侧的V_i加I_k大约在1到2附近数值条件好很多。如果你直接用实际kW和欧姆建模大M要跟着跳好几个数量级求解器数值稳定性会非常差。3. 数学建模DistFlow、SOCP、重构约束3.1 DistFlow方程组的基础配电网的潮流方程最常用的是DistFlow形式。考虑一条支路k从节点i流向节点j定义P_k和Q_k为支路首端流向末端的有功和无功V_i为节点i电压幅值的平方I_k为支路电流幅值的平方。DistFlow的基本方程是节点功率平衡流入节点j的支路功率减去支路损耗等于节点j的负荷减去注入电源功率再减去所有子支路流出的功率。电压降方程V_j等于V_i减去支路上的电压降。具体展开就是V_j V_i - 2*(r_kP_k x_kQ_k) (r_k^2 x_k^2)*I_k。支路电流关系中包含一个关键的非凸等式I_k (P_k^2 Q_k^2) / V_i。这条约束把有功、无功、电压和电流的平方关系绑在一起是非凸的。如果不做任何松弛放进带0-1变量的重构模型里就是MINLP求解难度陡增。3.2 二阶锥松弛的核心原理二阶锥松弛的动机其实很朴素把上面那个非凸等式放宽成不等式。允许电流平方不小于物理需求即I_k (P_k^2 Q_k^2) / V_i。在电压非负的配电网里这个放宽意味着给潮流解留出了一个“电流比实际需要更大”的空间。那为什么松弛之后最优解还能回到等式因为重构问题最常见的优化目标是网损最小而网损正是r_k*I_k的累加目标函数会尽量压低I_k。既然I_k被约束抬高了下界最优解自然会把它压到下界上于是松弛后的解自动满足原等式。这就是“松弛是精确的”的工程直觉。严格的理论条件包括辐射状网络、合理的电压边界、目标函数单调递增于电流等实际判断还是靠跑完以后回代验证。把这个不等式整理成YALMIP可以直接吃进去的标准二阶锥形式是后面所有代码的核心。标准形式是|| [2P_k; 2Q_k; V_i - I_k] ||_2 V_i I_k注意这里的V是电压幅值的平方不是幅值本身。这套变量定义要贯穿始终一旦中途混入“真实电压幅值”锥约束的尺度立刻乱掉。3.3 重构开关的big-M建模当支路开关状态z_k1时支路满足完整的DistFlow当z_k0时支路断开P_k、Q_k、I_k全部应为0同时支路两端电压不再受电压降方程限制。如何把这种“或开或闭”的逻辑表达成线性约束常规做法是大M法。对于功率和电流可以写-Mpz_k P_k Mpz_k -Mpz_k Q_k Mpz_k 0 I_k Ml*z_k对于电压关系需要把电压降方程改写成大M形式。定义expr为V_j - V_i 2*(r_kP_k x_kQ_k) - (r_k^2 x_k^2)I_k则用Mv(1-z_k)把expr夹在上下界之间。z_k1时expr必须等于0恢复原始电压方程z_k0时expr被Mv放宽两端电压可以自由取值。大M的取值是个实操经验点。M太小会把可行域削掉导致原本应该可行的拓扑被误判不可行M太大会让分支定界变得非常慢数值上也可能出问题。我个人的习惯是Mp取全网总负荷的1.2到2倍Ml取最大允许电流平方的1.2倍Mv取Vmax^2 - Vmin^2的1.2倍左右。这样每个大M都有明确的物理尺度不会莫名其妙的病态。3.4 辐射状约束虚拟流法辐射状约束是重构建模里最容易出错的地方。只约束“闭合支路数节点数-1”是完全不够的可能出现闭环加孤岛的组合。一种经典且好实现的方法是“虚拟流法”也叫单商品流法。思路是这样的给每条支路定义两个方向的非负虚拟流变量f_pos和f_neg。f_pos表示从支路首端流向末端f_neg表示从末端流回首端。每条支路的两个方向虚拟流之和不能超过(n_bus-1)再乘以开关状态z_k也就是断开支路不能通过任何虚拟流。然后设置源和需求根节点变电站节点1向外净流出n_bus-1个单位其他每个非根节点净流入1个单位。这个方法的物理直觉很清晰。把每个非根节点想象成需要1单位“虚拟货物”根节点是唯一供应商。如果网络里有环环路会破坏每个节点的净流入守恒如果网络不连通孤岛里的节点拿不到那1单位虚拟流。所以只要虚拟流约束满足网络必然连通且无环也就是一棵生成树。3.5 目标函数的设计目标函数最常见的是全天网损最小即最小化sum_t sum_k r_k * I_k(t)。这也是最容易和二阶锥松弛配合的目标因为网损对I_k是单调递增的能保证SOCP松弛在最优解处紧起来。如果想同时限制开关操作次数可以在目标里给切换变量加权或者把总切换次数作为约束。实际配电系统里频繁操作开关既不经济也增加故障概率所以“开关动作次数限制”几乎是重构问题必带的条件。具体线性化方式我在第5节详细说。除了网损也可以用电压偏差最小、DG消纳最大等目标但要注意目标函数如果不是对I_k单调递增二阶锥松弛的精确性就需要额外验证不能默认成立。4. YALMIP建模实操从变量定义到CPLEX求解4.1 变量声明与数据索引先把参数读进来假设已经有branch_data和bus_data两个表格。关键变量是支路有功P、支路无功Q、节点电压平方V、支路电流平方I以及开关状态z。多时间断面下P、Q、I、V都带时间维度而z是否带时间维度取决于重构策略。如果全天只能重构一次z就是(n_branch,1)的二进制变量如果允许每个时段调整但限制总切换次数z就是(n_branch,n_time)的二进制变量。n_bus size(bus_data, 1); n_branch size(branch_data, 1); n_time 24; P sdpvar(n_branch, n_time, full); Q sdpvar(n_branch, n_time, full); V sdpvar(n_bus, n_time, full); I sdpvar(n_branch, n_time, full); z binvar(n_branch, 1); % 全天一个拓扑 % 如果每时段都可以调整开关用下面这行 % z binvar(n_branch, n_time); f_pos sdpvar(n_branch, 1, full); f_neg sdpvar(n_branch, 1, full);一个规模上的概念IEEE33加24个时段P、Q、I已经是3乘32乘24等于2304个连续变量再加上V的792个变量z的37个二进制变量总规模约3000个变量。这已经不能算小模型了PG69会更大。所以变量声明阶段就要对模型规模有心理准备。4.2 约束搭建的关键代码约束搭建是整个模型的核心。下面这一段是支路层面的约束包括开关状态与功率电流的关联、二阶锥约束、电压降的大M约束Constraints []; Mp 2 * sum(bus_data(:, 3)); % 有功大M Ml 1.2 * max_I_square; % 电流平方大M按数据设定 Mv 1.2 * (Vmax^2 - Vmin^2); % 电压平方差大M for t 1:n_time for k 1:n_branch i branch_data(k, 1); j branch_data(k, 2); r branch_data(k, 3); x branch_data(k, 4); Constraints [Constraints, -Mp*z(k) P(k,t) Mp*z(k)]; Constraints [Constraints, -Mp*z(k) Q(k,t) Mp*z(k)]; Constraints [Constraints, 0 I(k,t) Ml*z(k)]; Constraints [Constraints, ... norm([2*P(k,t); 2*Q(k,t); V(i,t) - I(k,t)], 2) V(i,t) I(k,t)]; expr V(j,t) - V(i,t) 2*(r*P(k,t) x*Q(k,t)) - (r^2 x^2)*I(k,t); Constraints [Constraints, -Mv*(1-z(k)) expr Mv*(1-z(k))]; end % 节点功率平衡 for m 1:n_bus to_idx find(branch_data(:,2) m); from_idx find(branch_data(:,1) m); Constraints [Constraints, ... sum(P(to_idx,t) - branch_data(to_idx,3).*I(to_idx,t)) sum(P(from_idx,t)) ... Pg(m,t) Pd(m,t)]; Constraints [Constraints, ... sum(Q(to_idx,t) - branch_data(to_idx,4).*I(to_idx,t)) sum(Q(from_idx,t)) ... Qg(m,t) Qd(m,t)]; end Constraints [Constraints, Vmin^2 V(:,t) Vmax^2]; end节点功率平衡方程的含义可以这样看流入节点m的所有支路功率经过支路损耗后到达m从m流出的所有支路功率从m出发。这两部分加起来再加上节点注入电源功率必须等于节点负荷。断开支路因为P、Q、I全为0自动退化不需要额外处理。4.3 虚拟流辐射状约束代码如果z是全天一个值虚拟流变量也只需要单断面。约束写起来比较直接Constraints [Constraints, 0 f_pos (n_bus - 1) * z]; Constraints [Constraints, 0 f_neg (n_bus - 1) * z]; root_out sum(f_pos(branch_data(:,1)1)) sum(f_neg(branch_data(:,2)1)); root_in sum(f_pos(branch_data(:,2)1)) sum(f_neg(branch_data(:,1)1)); Constraints [Constraints, root_out - root_in n_bus - 1]; for m 2:n_bus in_m sum(f_pos(branch_data(:,2)m)) sum(f_neg(branch_data(:,1)m)); out_m sum(f_pos(branch_data(:,1)m)) sum(f_neg(branch_data(:,2)m)); Constraints [Constraints, in_m - out_m 1]; end这段代码的索引细节很容易出错。f_pos(k)代表从branch_data(k,1)流向branch_data(k,2)的虚拟流f_neg(k)代表反方向。所以对节点m来说流入量是所有“首端为m、流向m之外的f_neg”加上所有“末端为m、从外流向m的f_pos”流出量正好相反。建议写完后单独测试一下这个子模型固定一个已知辐射状拓扑把z设成对应值检查虚拟流约束是否可行再故意设一个带环的拓扑看约束是否不可行。这个测试值得做一次。4.4 求解器调用与结果提取目标函数和求解设置如下Obj 0; for t 1:n_time Obj Obj sum(branch_data(:,3) .* I(:,t)); end ops sdpsettings(solver, cplex, verbose, 2, cplex.timelimit, 3600); sol optimize(Constraints, Obj, ops); if sol.problem 0 P_val value(P); Q_val value(Q); V_val value(V); I_val value(I); z_val value(z); fprintf(objective %.4f\n, value(Obj)); else disp(sol.info); endYALMIP的好处是建模和求解器解耦。写成solvercplex没问题如果机器上装的是Gurobi把cplex改成gurobi就能直接跑模型部分不用动。这也是我坚持用YALMIP而不是直接写CPLEX接口的原因——调试阶段省非常多事。4.5 对偶间隙与松弛精确性验证跑完MISOCP之后绝对不能直接拿结果去写论文必须先验证SOCP松弛是不是精确的。验证方法很直接回到原始非凸等式逐条支路检查电流平方是否真的等于功率平方除以电压平方for t 1:n_time for k 1:n_branch i branch_data(k,1); violation(k,t) abs(I_val(k,t) * V_val(i,t) - (P_val(k,t)^2 Q_val(k,t)^2)); end end max_violation max(violation(:));如果max_violation在1e-4以下标幺值体系下说明松弛紧重构结果可信。如果这个数很大优先怀疑三件事大M是不是取太大导致数值病态电压上下限是不是太宽给松弛留了多余空间目标函数里是不是加了会削弱“压缩电流”的东西。5. 多时间断面耦合与DG/储能处理5.1 储能SOC模型多时间断面里储能是最典型的跨时段耦合元件。设储能节点有功为P_ess(t)充电为正放电为负但实际建模时我习惯拆成两个非负变量P_ch和P_dis分别表示充电功率和放电功率SOC(t1) SOC(t) (P_ch(t) * eta_ch - P_dis(t) / eta_dis) * dt / E_rate这个方程把各个时段真正串起来了。SOC(1)通常设0.5末端要求SOC(241)不小于0.5保证一个调度周期内能量平衡。SOC本身要限制在0到1之间。如果不加任何约束优化问题可能出现同时充电和放电的荒谬解因为这样既能“吸收”多余光伏又能“支撑”电压代价却可能被网损目标掩盖。解决办法有两个一是加一个很小的惩罚项在目标里加epsilon乘以(P_chP_dis)让同时充放电变得不划算二是用二进制变量严格互斥P_ch受Mu限制P_dis受M(1-u)限制。对精度要求高的场景建议用后者。5.2 开关动作次数的线性化如果允许每个时段都调整开关但全天总切换次数不能超过K就需要把“是否发生切换”表达成线性约束。这里有个经典坑直接用abs(z(k,t1)-z(k,t))YALMIP虽然能处理但引入的模型结构往往很差求解速度明显下降。我建议手动引入二进制变量delta做精确线性化delta binvar(n_branch, n_time - 1); for k 1:n_branch for t 1:n_time - 1 Constraints [Constraints, delta(k,t) z(k,t1) - z(k,t)]; Constraints [Constraints, delta(k,t) z(k,t) - z(k,t1)]; Constraints [Constraints, delta(k,t) z(k,t) z(k,t1)]; Constraints [Constraints, delta(k,t) 2 - (z(k,t) z(k,t1))]; end end Constraints [Constraints, sum(sum(delta)) K];简单解释一下这组约束。前两行保证一旦开关状态发生变化delta至少要等于1。后两行保证开关状态没变化时delta被压成0两个状态都是0时z_t加z_t1等于0给出delta小于等于0两个状态都是1时2减去两者之和等于0同样给出delta小于等于0。这样delta就精确等于开关变化次数。这个线性化方式在整数规划里非常干净CPLEX处理起来速度很快。5.3 时序数据与断面数量的折中多时间断面不是断面越多越好断面数量对求解时间影响是指数级的。我调试的时候通常先用8个代表断面把流程跑通验证模型逻辑和结果合理性再逐步增加到24个断面跑正式结果。如果系统是PG69加上储能24断面的MISOCP在普通笔记本上可能要跑几十分钟甚至更久要有心理准备。必要时可以给CPLEX设置强时间限制先拿到一个可行解用于分析再慢慢往上提质量。6. 常见问题与调试实录6.1 CPLEX Community Edition的规模限制这是很多同学实际跑代码遇到的第一个坎。CPLEX Community Edition是免费版本但它对模型规模有严格限制我记得变量和约束都限制在1000量级具体数字建议以IBM官方文档为准。IEEE33单断面重构模型小CE能跑但是IEEE33加24时段、或者PG69加24时段变量数量直接冲到几千甚至上万CE会直接报规模超出限制根本进不了求解阶段。解决思路有几个。高校和科研单位可以申请IBM的学术版授权完全免费求解大规模MISOCP无压力。如果申请不到Gurobi也有学术版YALMIP切换求解器只需要改一个参数。开源方案可以试试SCIP它支持MISOCP速度比CPLEX和Gurobi慢但胜在没有任何授权限制用来验证算法逻辑完全够用。6.2 求解报错信息与状态码YALMIP返回的sol.problem值是判断求解状态的第一手信息。problem等于0表示成功等于1表示不可行等于2表示无界其他数字往往和数值问题有关。初学者看到不可行往往直接蒙了但实际上不可行最常见的原因是节点功率平衡写反了方向。我在代码里统一约定“流入节点m的支路功率加流出节点m的支路功率等于注入减去负荷”如果direction搞反模型一定不可行。数值警告也不少见尤其是你用了太大的M。一个M设为1e8的约束在CPLEX内部是把量级1的量放到1e8级别去衡量分支定界过程中数值退化非常严重。解决办法就是按第3节的经验针对不同约束分别取有物理含义的M值别一个M走天下。6.3 拓扑约束不满足怎么办跑完之后必须检查开关状态对应的网络拓扑。我习惯用MATLAB的graph和conncomp函数验证G graph(zeros(n_bus)); for k 1:n_branch if z_val(k) 0.5 G addedge(G, branch_data(k,1), branch_data(k,2)); end end cc conncomp(G); num_components max(cc); closed_branches sum(z_val 0.5); n_loops closed_branches - n_bus num_components;正常结果是num_components等于1n_loops等于0。如果连通分量大于1说明虚拟流约束写错了方向如果n_loops大于0说明辐射状约束没有完全生效。多数情况下问题都出在根节点净流出的符号上耐心推一遍就好。6.4 数值尺度问题配电网模型比输电网更容易出现数值问题因为支路阻抗量级差异大。一边是0.01欧姆的短线路一边是0.5欧姆的长线路直接在SI单位制下建模约束矩阵的条件数会非常难看。我强烈建议全部走标幺值。基准容量取10MVA基准电压取12.66kV所有功率和阻抗都落到0.001到1这个区间CPLEX的数值稳定性会好很多也能减少很多莫名其妙不可行的问题。6.5 对称性和求解加速技巧多时间断面模型有一个隐蔽的麻烦如果所有时段的结构和变量完全对称只是数据不同CPLEX在分支定界时会在对称区域里反复搜索浪费大量时间。一个很好用的加速技巧是给MIP提供一个初始可行解。我先跑一个单断面重构得到一组开关状态然后用assign函数把它赋给z把ops.usex0设为1。CPLEX会把这组解作为MIP start经常能把求解时间从半小时压到几分钟。assign(z, single_period_z); ops sdpsettings(solver, cplex, usex0, 1);另外对PG69这种规模较大的系统可以按支路负荷大小或历史开断频率给开关变量设置分支优先级CPLEX会优先分支更关键的变量剪枝效率往往有明显提升。不过这个属于偏高级的调参建议先跑通原始模型再加这些优化。7. 结果怎么看与后续扩展重构结果的一般形态是网损下降、电压曲线被拉平。以IEEE33为例基础拓扑的日网损通常在200kW量级重构之后能降到130到160kW左右具体取决于负荷曲线形态和是否允许多次开关。这个下降幅度可以作为结果的sanity check如果重构后网损反而比基础拓扑高很多那一定哪里出了问题。多时间断面带来的真正价值是它能给出一个“全天综合最优”的拓扑。我跑过不少案例中午光伏大发的单断面最优拓扑和晚上负荷高峰的单断面最优拓扑明显不同而多断面模型给出的折中拓扑虽然不一定在任何一个单断面达到最优但日积分网损比固定采用任意一个单断面最优方案都要低几个百分点。这就是多时间断面“耦合”的意义。后续扩展方向其实非常多。可以把确定性优化改成鲁棒优化或随机优化应对光伏预测误差可以在重构的同时做储能容量配置变成规划与运行联合优化可以把单相DistFlow换成三相不平衡模型也可以在中压配网里考虑N-1约束和故障重构。这些都是这个基础模型往上长出来的分支。最后分享一个我自己实际写代码的体会在动手写YALMIP之前先在纸上把每一个变量的下标、每一个约束的物理意义完整写一遍尤其把支路方向约定写清楚。这一步花不了半小时但能省掉后面整整一两天的调试时间。模型逻辑顺了代码只是翻译模型逻辑没顺代码写得再快也是白搭。