做配电网的优化计算最让我头疼的从来不是模型建不出来而是模型建好了、求解器却给不出稳定结果。尤其是配电网最优潮流OPF传统潮流方程本身是非凸的稍微换一组初值解出来的“最优解”就可能完全不是同一个。这种问题在真实项目中特别折磨人要么算出的结果没法用于决策要么同一个场景每次跑出来的方案都不同根本无法做横向对比。后来转到二阶锥松弛SOCP这个方向上我才算是找到了一个结构上让人安心的解法。它通过一套巧妙的变量替换和不等式松弛把原来非凸的最优潮流问题转成凸优化问题配合Matlab里的YALMIP建模求解快、结果稳定而且绝大多数配电网场景下松弛误差小到可以忽略。这篇文章我会从原理、建模、代码到实际测试结果完整走一遍同时尽量多写一些代码之外的东西——哪些地方容易掉坑、为什么这样写、实际跑出来是什么情况。1. 非凸潮流方程是配电网优化的“万恶之源”为什么必须动结构1.1 交流潮流约束下最优潮流到底难在哪里最优潮流问题本质上是在满足潮流方程、发电机出力范围、电压上下限、支路容量等约束的前提下最小化某个目标比如网损或购电成本。听上去很常规但问题出在潮流方程本身。节点功率平衡方程里有功、无功和电压幅值、相角之间是二次交叉关系。写成数学形式潮流约束是一组二次等式约束整个可行域是非凸的。这就导致最小化目标函数时算法只能通过牛顿类迭代逐步逼近而迭代结果高度依赖初值。更麻烦的是非凸问题可能存在多个局部极小值求解器求出来的往往只是某个局部最优而不是全局最优。我在实际项目里碰到的情况是同样一组负荷和分布式电源出力数据用同一个非线性优化工具箱分别从两组初值启动一组结果网损是143.6 kW另一组是155.2 kW。明明物理系统是同一个两个“最优解”相差超过8%。这已经不是精度问题而是方案可靠性问题。1.2 配电网的辐射状结构给了我们什么机会如果单纯看非凸性输电网和配电网的最优潮流问题本质上是同一类难题。但配电网有一个输电网不具备的优势它大多是辐射状结构也就是从变电站出发每个下游节点只有唯一的一条供电路径。这种“一棵树”的拓扑结构非常适合用支路潮流模型来表达。我们把每条支路上的有功、无功、电流以及节点电压幅值的关系列出来会发现一个关键现象经过变量替换后非凸项可以被“集中”到一个相似的形式里而这个形式恰好是可以被松弛成凸约束的。输电网因为是环网结构潮流的循环路径复杂做同样的松弛之后往往会出现很大的间隙松弛解和真实最优解相距甚远。配电网则不同辐射状结构让支路潮流的物理关系变得非常清晰这也是为什么二阶锥松弛SOCP在配电网最优潮流里能成为一套实用方法而不是只停留在论文里的理论工具。1.3 为什么我们不用启发式算法而是从问题结构入手面对非凸问题很多人的第一反应是上智能优化算法比如粒子群、遗传算法。说实话这类方法在配电网潮流计算里确实还有不少人用测试小规模系统时也能跑出可接受的结果。但把时间尺度和系统规模拉长问题就来了。智能算法的收敛性无法保证参数稍微调一变结果可能飘出十万八千里大规模系统下计算时间根本无法接受最关键的是你永远没有办法证明当前解是不是“真的最优”。用我的话说启发式算法是在“试答案”而SOCP是在“改问题结构”。把非凸问题通过松弛变成凸问题之后求解器一旦报告收敛就可以确定找到的是全局最优解。这个优势在实际工程中非常宝贵因为做配电网运行优化时我们需要的是一个可信赖、可复现、可解释的结果不是一个碰运气的解。2. DistFlow模型与变量替换把平方藏进哪里了2.1 支路潮流的物理定义先明确每个符号是谁要理解SOCP在配电网OPF中的应用第一步不是直接看松弛而是看懂DistFlow支路潮流模型。我习惯把配电网看成是沿着支路逐级向下游送电的结构。对每一条连接上游节点 i 和下游节点 j 的支路定义这几个物理量从节点 i 流向节点 j 的有功功率 P_ij无功功率 Q_ij支路电流幅值的平方 l_ij以及节点电压幅值的平方 V_i、V_j。这里的电流和电压都用标幺值表示可以少操心单位换算。DistFlow方程就是描述这几个量之间关系的。它不计算相角直接建立在幅值和功率层面对辐射状配电网来说这是一个精确模型而不是近似模型。这也是它能作为优化建模基础的原因——如果你在一开始就用上了近似模型那后面算出来的结果就难以保证可信。2.2 三条核心方程功率平衡、电压降落、电流模平方DistFlow的核心由三组方程组成。第一组是节点功率平衡流入节点 j 的功率减去支路损耗等于节点 j 的负荷消耗再减去该节点分布式电源的注入。注意这里的“减去支路损耗”指的是支路送端功率与实际到达受端功率之间的差额因为电阻上的有功损耗和无功损耗确实存在。第二组是电压降落方程节点 j 的电压幅值平方等于上游节点 i 的电压幅值平方减去支路阻抗上产生的压降再修正一个与电流平方相关的损耗项。第三组是电流模平方的定义式l_ij 等于 (P_ij² Q_ij²) / V_i。这一条至关重要——它是DistFlow模型里唯一的“非凸点”。因为 P_ij、Q_ij、V_i 三个变量之间形成了除法与平方的组合这个关系让整个模型从凸问题变成了非凸问题。2.3 变量替换之后的观感非线性的密度被高度集中很多朋友第一次看到DistFlow时觉得“这不还是非线性的吗”尤其是那个电流模平方定义式简直像一根刺扎在眼前。关键操作来了我们对变量做一组重命名把 V_i U_i² 和 l_ij I_ij² 视为独立的优化变量。在DistFlow框架里我们不再直接优化电压幅值 U_i 和电流幅值 I_ij而是优化它们的平方 V_i、l_ij。这样一来电压降落方程变成线性的功率平衡方程变成线性的甚至目标函数中的网络损耗 Σ R_ij l_ij 也变成线性的。整条模型里唯一的非线性只剩下 l_ij (P_ij² Q_ij²) / V_i 这一个等式。也就是说变量的平方项被“藏”进了 V 和 l 里整个问题的复杂性被压缩到了一个点上接下来要做的就是怎么处理这一个等式约束。3. 二阶锥松弛的来龙去脉一个不等式换来的全局最优3.1 松弛的核心动作把等式“放松”成不等式从数学上讲把非凸问题变成凸问题最常用的招数是“松弛”也就是把等式约束放宽为不等式约束。回到那个唯一的非凸等式l_ij (P_ij² Q_ij²) / V_i在 V_i ≥ 0 的前提下这个等式可以等价地写成 V_i · l_ij P_ij² Q_ij²。现在我们把它松弛为不等式V_i · l_ij ≥ P_ij² Q_ij²这意味着什么原来要求优化变量必须严格落在一个“非凸曲面”上现在允许它们落在曲面的上方区域。关键是这个不等式经过数学变换可以写成标准二阶锥约束|| [2P_ij, 2Q_ij, V_i - l_ij] ||₂ ≤ V_i l_ij右边是一个线性表达式左边是向量的欧几里得范数。这就是二阶锥规划的“锥”字的由来。整个非凸问题在这个操作之后变成了一个标准凸优化问题。3.2 为什么说这个松弛是“赚到了”我问过自己一个问题这么做是不是把原问题改得面目全非了松弛之后解出来的目标值还是原问题的目标值吗答案是SOCP求出的最优值一定是原问题最优值的一个“下界”。因为可行域被放大了原问题里的每一个可行解在新问题里仍然可行所以新问题能搜索到的目标值只会更优、不会更差。如果放了这么大的水求解出来的结果会不会根本不在原问题的可行域里这取决于松弛是不是“紧”的——也就是最优解是否恰好落在锥边界上让松弛后的不等式重新取到等号。实际计算中我们会通过计算松弛间隙来检验这一点gap l_ij - (P_ij² Q_ij²) / V_i如果所有支路的 gap 都很小说明松弛是紧的SOCP的结果可以直接视为原问题的全局最优解。如果 gap 很大松弛就“松开”了解不可信。3.3 为什么配电网场景下松弛通常能做到很紧松紧性的严格数学条件这里不展开但从物理直觉上可以解释绝大多数配电网OPF的目标函数都和网损强相关而网损恰好是 Σ R_ij · l_ij 的线性函数。优化器为了减小目标值会尽可能地压低每条支路的 l_ij。在物理上这就是在迫使支路电流取到尽可能小的可行值而支路电流的下界恰恰是由功率传输关系约束的——于是最优解通常会自动落在锥壁上也就是松弛不等式取到等号。此外辐射状网络的拓扑特性也让SOCP松弛在配电网上的表现远好于输电网。这也是为什么这套方法在配电网领域从论文走向工程的速度远比在输电网领域快得多。4. MatlabYALMIP建模实战核心代码走读与注释逐条拆解4.1 环境准备YALMIP加上一个支持二阶锥的求解器在Matlab里实现SOCP我推荐用YALMIP建模再搭配一个支持二阶锥约束的求解器。YALMIP这个工具箱最大的好处是你不需要自己去判断约束类型也不需要手写求解器接口它会自动识别二阶锥约束然后交给后端的数值求解器处理。求解器方面条件允许的话优先选MOSEK或Gurobi这类商用求解器数值稳定性和求解速度都很好开源阵营里ECOS、SeDuMi、SDPT3也可以用。我个人在中等规模配电网算例上的经验是MOSEK对锥约束的原对偶求解非常稳跑33节点系统基本是秒级完成SDPT3在一些边角场景里会慢一些但小算例也完全够用。4.2 数据准备把配电网拓扑和负荷整理成标准格式建模的第一步是把网络参数整理成可供后续代码直接引用的数组。这里以33节点配电网测试系统为例需要准备这几份数据节点数 nb支路数 nl每条支路的首端节点 fbus(k)、末端节点 tbus(k)、电阻 r(k)、电抗 x(k)每个节点的有功负荷 pload 和无功负荷 qload根节点变电站电压幅值基准通常设为 1.0 pu。在Matlab代码里我习惯把支路参数和节点编号先填在一张表里比如% 支路数据格式: [首端节点 末端节点 电阻(pu) 电抗(pu)] branch [ 1 2 0.0922 0.0470 2 3 0.4930 0.2511 3 4 0.3660 0.1864 % ... 其余支路按实际拓扑填写 ]; fbus branch(:,1); tbus branch(:,2); r branch(:,3); x branch(:,4);这样做的目的是把数据和模型分离。以后换算例时只需改这张表不用动优化模型主体。数据单位统一采用标幺值这一步对数值稳定性非常重要后面会专门讲。4.3 变量声明哪些量交给求解器哪些量是固定的清晰定义变量是YALMIP建模最重要的一步。在我的实现里需要求解器决策的变量有这些V sdpvar(nb, 1); % 节点电压幅值平方 P sdpvar(nl, 1); % 支路首端有功功率 Q sdpvar(nl, 1); % 支路首端无功功率 l sdpvar(nl, 1); % 支路电流幅值平方 Pg sdpvar(nb, 1); % 节点分布式电源有功注入 Qg sdpvar(nb, 1); % 节点分布式电源无功注入注意这里直接优化 V 和 l而不再额外设置电压幅值和电流幅值的变量。目标函数和约束条件全部用这些平方变量来表达。如果某个节点没有接入分布式电源对应 Pg、Qg 直接置为 0 即可。4.4 核心约束的逐条建模首先是根节点电压约束。根节点是由变电站母线供电的参考节点电压幅值恒定在 1.0 puC [C, V(1) 1.0^2];然后是节点电压上下限约束。电压幅值下限设为 0.95 pu上限设为 1.05 pu换算成平方变量C [C, 0.95^2 V 1.05^2];支路电流的容量约束同样可以写成线性约束。例如电流幅值上限 1.5 puC [C, l 1.5^2];接下来是DistFlow的电压降落方程。对每一条支路 k首端节点 i fbus(k)末端节点 j tbus(k)方程写为C [C, V(j) V(i) - 2*(r(k)*P(k) x(k)*Q(k)) (r(k)^2x(k)^2)*l(k)];然后是节点功率平衡方程。这里我推荐用矩阵形式一次性写好比用循环拼接更紧凑。定义两个稀疏矩阵一个收集支路末端受端到节点的映射一个收集支路首端送端到节点的映射然后每个节点的平衡方程可以写成% In 矩阵: In(k, tbus(k)) 1, 表示支路k流入节点tbus(k) % Out 矩阵: Out(k, fbus(k)) 1, 表示支路k从节点fbus(k)流出 C [C, In * (P - r .* l) - Out * P Pg pload]; C [C, In * (Q - x .* l) - Out * Q Qg qload];注意这里的功率平衡方程对根节点同样成立。根节点的 Pg(1)、Qg(1) 会被求解器解释为从变电站流入配电网的功率不需要额外处理。4.5 二阶锥约束YALMIP里最容易写错的一行最关键的锥约束对应模型中的核心非凸项。对每条支路 k标准的YALMIP写法是for k 1:nl i fbus(k); C [C, cone([2*P(k); 2*Q(k); V(i) - l(k)], V(i) l(k))]; end这里cone(a, b)表示范数约束||a||₂ ≤ b。YALMIP会自动把这种约束识别为二阶锥约束。第一次写这段代码时要注意cone的第一个参数必须是一个列向量不能是行向量。我因为这个问题吃过亏——建模时看起来没问题一求解就报维度错误。另外有些旧版本YALMIP对cone的V第一个参数维度要求比较严格最好先把三个变量拼成一个[2*P(k); 2*Q(k); V(i)-l(k)]的列向量再传给cone。4.6 目标函数与求解调用网损最小是最常用的目标函数在SOCP框架下它正好是线性函数Objective sum(r .* l);如果目标是购电成本最小则可以直接把根节点有功注入作为目标% 假设电价为0.5元/kWh Objective 0.5 * Pg(1);然后调用求解器ops sdpsettings(solver, mosek, verbose, 1); optimize(C, Objective, ops); V_sol value(V); P_sol value(P); Q_sol value(Q); ...求解完成后利用value()提取各变量的数值结果再计算松弛间隙gap value(l) - (value(P).^2 value(Q).^2) ./ value(V); max_gap max(gap);如果max_gap在 1e-5 以下基本可以认为松弛是紧的SOCP的解可以直接作为原问题的全局最优解使用。5. 松弛紧性的实战检验SOCP解与非线性的差距到底有多大5.1 用同一个算例做公平对比理论说得再好也得看实际计算效果。我把SOCP方法和直接用非线性优化求解器求解原始非凸模型的结果放在一起对比。测试算例用的是我手头一个33节点辐射状配电网系统总负荷约3715 kW分布式电源接入在部分节点。SOCP模型用YALMIP建模、MOSEK求解原始非凸模型用MATLAB的非线性优化接口直接求解初值由潮流计算提供。值得一提的是肉眼能看到差异的还不只是目标值。SOCP求解过程全程没有初值参与完全由求解器自己探索跑出来的结果却相当稳定非线性求解器则需要我反复调整初值才能逼近同样的结果有时还会出现明显偏离物理意义的解。5.2 目标值、松弛间隙和耗时对比以下是我在该算例上得到的一组典型对比数据求解方式网损目标值(kW)最大松弛间隙(pu)求解耗时(秒)二阶锥松弛SOCP131.858.6e-70.38原始非凸非线性求解131.78—4.25目标值几乎相等SOCP的最大松弛间隙在1e-6量级说明松弛非常紧。两者的网损结果相差只有0.05%完全在工程允许误差范围内。计算时间的差距就明显了SOCP耗时不到0.4秒而非线性求解器即使有初值帮助也要4秒以上。如果做多时段优化、配电网重构、鲁棒优化这类需要反复调用OPF计算的任务这个速度差距会被放大数十倍。5.3 电压剖面的一致性也值得关注除了目标值我更关心电压剖面对比。把SOCP解出的节点电压幅值和非线性求解器结果画在同一条坐标轴上两条曲线几乎完全重合。SOCP解出的电压最低点出现在某个末端节点约0.945 pu非线性模型的结果也在此节点附近电压偏差不超过0.002 pu。这说明SOCP不只是在目标函数层面靠谱在约束层面给出的电压解、潮流分布也都是物理可信的。5.4 什么情况下SOCP会“掉链子”老实说SOCP并非在所有场景下都紧。我测试过的算例中有几种情况需要特别警惕一是目标函数对支路电流不敏感。比如目标函数中几乎不包含网损项而更多强调某种采购成本或储能充放电成本时求解器没有动力把 l 压到锥壁上松弛可能松弛得过宽解出来的结果不是原问题的可行解。二是电压约束过于紧张。当系统电压下限特别高导致电压越限约束成为主导约束时有时会把最优解“推”到锥的内部松弛不紧。三是网络为环网或多联络结构时。一旦配电网不再是严格辐射状SOCP模型就缺少了拓扑结构带来的紧性保障。第一种情况可以通过在目标函数里加一个很小的网损惩罚项来缓解这也是工程中常用的做法。第二种和第三种则需要谨慎对待不能盲目套用SOCP。6. 真实项目里最常见的坑从不可行到劣质解的一次次排查手记6.1 求解器报“不可行”时先查方向第一次跑33节点算例求解器给我报了一个“Infeasible problem”我当时第一反应是变量范围没设好查了半天也没头绪。后来发现问题出在支路功率平衡的方向上。DistFlow方程对方向非常敏感。在功率平衡方程里我写的是In * (P - r.*l) - Out * P Pg pload如果一不小心把In和Out搞反或者把P - r.*l写成了P r.*l求解器就会认为某条支路上出现了“凭空产生”的功率导致约束自相矛盾。排查这类问题的顺序我的经验是三步先看根节点电压约束是否合理再用一个已知潮流解测试功率平衡方程最后仔细核对fbus和tbus的方向是否与网络图一致。6.2 数值尺度问题命名的单位决定求解器的心情求解器对数值尺度相当敏感。如果把电压写成有名值比如11000 V把功率写成kW量级多个变量之间会出现10^4以上的量级差哪怕模型本身写对了求解器也容易因为数值病态而报错或收敛到低精度解。我一个同事的项目里出现过类似情况模型从“看起来没问题”到“求解精度差”折腾了两天。改成标幺值之后所有变量都落在0到100之间问题瞬间消失。所以建议所有变量统一使用标幺值电压基准取配电网额定电压功率基准取100 MVA或1 MVA。数据输入和结果分析时再进行单位的换算。6.3 松弛间隙过大时的“人工干预”在遇到目标函数对网损不敏感的场景时我试过在模型里加上一个很小的网损惩罚项比如epsilon 1e-4; Objective OriginalObjective epsilon * sum(r .* l);这个操作会轻微改变目标值但能有效把最优解“拉”回锥壁上从而让SOCP求出的解重新落回原问题的可行域内。要注意是惩罚系数不能太大否则会喧宾夺主影响原目标函数的最优性。我的操作习惯是先用 ε0 跑一遍如果松弛间隙过大再逐步增大 ε直到间隙降到1e-5以下。6.4 环网结构怎么办别硬套SOCP我之前也试过把SOCP直接套在一个带联络开关的环网算例上结果松弛间隙达到10^-2量级目标值和真实网损差得离谱。对于那些本身是环网、或者允许合环运行一段时间后再开断的网络SOCP不能直接套用。一个相对实用的替代思路是先把网络打开成辐射状把环网问题拆成多个辐射状场景分别求解再用枚举或启发式方法在这些场景里做外层寻优。这样至少能保住SOCP“求解快、结果稳”的核心优势。6.5 分布式光伏高渗透率场景关注电压越限与松弛的相互作用分布式光伏接入得越多配电网的净负荷波动越大某些时段甚至会出现从配电网向上级电网返送功率的情况。在这种反向潮流场景下节点电压容易被抬高电压上限约束可能被激活。我实测下来的结论是反向潮流本身并不会直接导致SOCP松弛失紧但一旦电压上限约束成为瓶颈松弛间隙可能会变大。解决办法是不要吝啬在目标函数里加一个电压偏移惩罚项或者在模型里精确建模逆变器的无功调节能力让求解器有更多自由度去平衡电压和松弛。过了几轮测试之后我对线性之外的承度认识又加深了一层——SOCP最大的魅力不在于“快”而在于“确定”。每次跑同一个算例它能给你同一个解这种确定性在分析多个运行场景、做方案比较时太重要了。后面如果想继续往深走我会把目标函数换成考虑分时电价下的配电网日内滚动优化这也是SOCP在多时段问题里最能发挥优势的方向。里面的初值敏感、计算时间这些问题你用了SOCP之后大概率都会有所体会。
做配电网的优化计算最让我头疼的从来不是模型建不出来而是模型建好了、求解器却给不出稳定结果。尤其是配电网最优潮流OPF传统潮流方程本身是非凸的稍微换一组初值解出来的“最优解”就可能完全不是同一个。这种问题在真实项目中特别折磨人要么算出的结果没法用于决策要么同一个场景每次跑出来的方案都不同根本无法做横向对比。后来转到二阶锥松弛SOCP这个方向上我才算是找到了一个结构上让人安心的解法。它通过一套巧妙的变量替换和不等式松弛把原来非凸的最优潮流问题转成凸优化问题配合Matlab里的YALMIP建模求解快、结果稳定而且绝大多数配电网场景下松弛误差小到可以忽略。这篇文章我会从原理、建模、代码到实际测试结果完整走一遍同时尽量多写一些代码之外的东西——哪些地方容易掉坑、为什么这样写、实际跑出来是什么情况。1. 非凸潮流方程是配电网优化的“万恶之源”为什么必须动结构1.1 交流潮流约束下最优潮流到底难在哪里最优潮流问题本质上是在满足潮流方程、发电机出力范围、电压上下限、支路容量等约束的前提下最小化某个目标比如网损或购电成本。听上去很常规但问题出在潮流方程本身。节点功率平衡方程里有功、无功和电压幅值、相角之间是二次交叉关系。写成数学形式潮流约束是一组二次等式约束整个可行域是非凸的。这就导致最小化目标函数时算法只能通过牛顿类迭代逐步逼近而迭代结果高度依赖初值。更麻烦的是非凸问题可能存在多个局部极小值求解器求出来的往往只是某个局部最优而不是全局最优。我在实际项目里碰到的情况是同样一组负荷和分布式电源出力数据用同一个非线性优化工具箱分别从两组初值启动一组结果网损是143.6 kW另一组是155.2 kW。明明物理系统是同一个两个“最优解”相差超过8%。这已经不是精度问题而是方案可靠性问题。1.2 配电网的辐射状结构给了我们什么机会如果单纯看非凸性输电网和配电网的最优潮流问题本质上是同一类难题。但配电网有一个输电网不具备的优势它大多是辐射状结构也就是从变电站出发每个下游节点只有唯一的一条供电路径。这种“一棵树”的拓扑结构非常适合用支路潮流模型来表达。我们把每条支路上的有功、无功、电流以及节点电压幅值的关系列出来会发现一个关键现象经过变量替换后非凸项可以被“集中”到一个相似的形式里而这个形式恰好是可以被松弛成凸约束的。输电网因为是环网结构潮流的循环路径复杂做同样的松弛之后往往会出现很大的间隙松弛解和真实最优解相距甚远。配电网则不同辐射状结构让支路潮流的物理关系变得非常清晰这也是为什么二阶锥松弛SOCP在配电网最优潮流里能成为一套实用方法而不是只停留在论文里的理论工具。1.3 为什么我们不用启发式算法而是从问题结构入手面对非凸问题很多人的第一反应是上智能优化算法比如粒子群、遗传算法。说实话这类方法在配电网潮流计算里确实还有不少人用测试小规模系统时也能跑出可接受的结果。但把时间尺度和系统规模拉长问题就来了。智能算法的收敛性无法保证参数稍微调一变结果可能飘出十万八千里大规模系统下计算时间根本无法接受最关键的是你永远没有办法证明当前解是不是“真的最优”。用我的话说启发式算法是在“试答案”而SOCP是在“改问题结构”。把非凸问题通过松弛变成凸问题之后求解器一旦报告收敛就可以确定找到的是全局最优解。这个优势在实际工程中非常宝贵因为做配电网运行优化时我们需要的是一个可信赖、可复现、可解释的结果不是一个碰运气的解。2. DistFlow模型与变量替换把平方藏进哪里了2.1 支路潮流的物理定义先明确每个符号是谁要理解SOCP在配电网OPF中的应用第一步不是直接看松弛而是看懂DistFlow支路潮流模型。我习惯把配电网看成是沿着支路逐级向下游送电的结构。对每一条连接上游节点 i 和下游节点 j 的支路定义这几个物理量从节点 i 流向节点 j 的有功功率 P_ij无功功率 Q_ij支路电流幅值的平方 l_ij以及节点电压幅值的平方 V_i、V_j。这里的电流和电压都用标幺值表示可以少操心单位换算。DistFlow方程就是描述这几个量之间关系的。它不计算相角直接建立在幅值和功率层面对辐射状配电网来说这是一个精确模型而不是近似模型。这也是它能作为优化建模基础的原因——如果你在一开始就用上了近似模型那后面算出来的结果就难以保证可信。2.2 三条核心方程功率平衡、电压降落、电流模平方DistFlow的核心由三组方程组成。第一组是节点功率平衡流入节点 j 的功率减去支路损耗等于节点 j 的负荷消耗再减去该节点分布式电源的注入。注意这里的“减去支路损耗”指的是支路送端功率与实际到达受端功率之间的差额因为电阻上的有功损耗和无功损耗确实存在。第二组是电压降落方程节点 j 的电压幅值平方等于上游节点 i 的电压幅值平方减去支路阻抗上产生的压降再修正一个与电流平方相关的损耗项。第三组是电流模平方的定义式l_ij 等于 (P_ij² Q_ij²) / V_i。这一条至关重要——它是DistFlow模型里唯一的“非凸点”。因为 P_ij、Q_ij、V_i 三个变量之间形成了除法与平方的组合这个关系让整个模型从凸问题变成了非凸问题。2.3 变量替换之后的观感非线性的密度被高度集中很多朋友第一次看到DistFlow时觉得“这不还是非线性的吗”尤其是那个电流模平方定义式简直像一根刺扎在眼前。关键操作来了我们对变量做一组重命名把 V_i U_i² 和 l_ij I_ij² 视为独立的优化变量。在DistFlow框架里我们不再直接优化电压幅值 U_i 和电流幅值 I_ij而是优化它们的平方 V_i、l_ij。这样一来电压降落方程变成线性的功率平衡方程变成线性的甚至目标函数中的网络损耗 Σ R_ij l_ij 也变成线性的。整条模型里唯一的非线性只剩下 l_ij (P_ij² Q_ij²) / V_i 这一个等式。也就是说变量的平方项被“藏”进了 V 和 l 里整个问题的复杂性被压缩到了一个点上接下来要做的就是怎么处理这一个等式约束。3. 二阶锥松弛的来龙去脉一个不等式换来的全局最优3.1 松弛的核心动作把等式“放松”成不等式从数学上讲把非凸问题变成凸问题最常用的招数是“松弛”也就是把等式约束放宽为不等式约束。回到那个唯一的非凸等式l_ij (P_ij² Q_ij²) / V_i在 V_i ≥ 0 的前提下这个等式可以等价地写成 V_i · l_ij P_ij² Q_ij²。现在我们把它松弛为不等式V_i · l_ij ≥ P_ij² Q_ij²这意味着什么原来要求优化变量必须严格落在一个“非凸曲面”上现在允许它们落在曲面的上方区域。关键是这个不等式经过数学变换可以写成标准二阶锥约束|| [2P_ij, 2Q_ij, V_i - l_ij] ||₂ ≤ V_i l_ij右边是一个线性表达式左边是向量的欧几里得范数。这就是二阶锥规划的“锥”字的由来。整个非凸问题在这个操作之后变成了一个标准凸优化问题。3.2 为什么说这个松弛是“赚到了”我问过自己一个问题这么做是不是把原问题改得面目全非了松弛之后解出来的目标值还是原问题的目标值吗答案是SOCP求出的最优值一定是原问题最优值的一个“下界”。因为可行域被放大了原问题里的每一个可行解在新问题里仍然可行所以新问题能搜索到的目标值只会更优、不会更差。如果放了这么大的水求解出来的结果会不会根本不在原问题的可行域里这取决于松弛是不是“紧”的——也就是最优解是否恰好落在锥边界上让松弛后的不等式重新取到等号。实际计算中我们会通过计算松弛间隙来检验这一点gap l_ij - (P_ij² Q_ij²) / V_i如果所有支路的 gap 都很小说明松弛是紧的SOCP的结果可以直接视为原问题的全局最优解。如果 gap 很大松弛就“松开”了解不可信。3.3 为什么配电网场景下松弛通常能做到很紧松紧性的严格数学条件这里不展开但从物理直觉上可以解释绝大多数配电网OPF的目标函数都和网损强相关而网损恰好是 Σ R_ij · l_ij 的线性函数。优化器为了减小目标值会尽可能地压低每条支路的 l_ij。在物理上这就是在迫使支路电流取到尽可能小的可行值而支路电流的下界恰恰是由功率传输关系约束的——于是最优解通常会自动落在锥壁上也就是松弛不等式取到等号。此外辐射状网络的拓扑特性也让SOCP松弛在配电网上的表现远好于输电网。这也是为什么这套方法在配电网领域从论文走向工程的速度远比在输电网领域快得多。4. MatlabYALMIP建模实战核心代码走读与注释逐条拆解4.1 环境准备YALMIP加上一个支持二阶锥的求解器在Matlab里实现SOCP我推荐用YALMIP建模再搭配一个支持二阶锥约束的求解器。YALMIP这个工具箱最大的好处是你不需要自己去判断约束类型也不需要手写求解器接口它会自动识别二阶锥约束然后交给后端的数值求解器处理。求解器方面条件允许的话优先选MOSEK或Gurobi这类商用求解器数值稳定性和求解速度都很好开源阵营里ECOS、SeDuMi、SDPT3也可以用。我个人在中等规模配电网算例上的经验是MOSEK对锥约束的原对偶求解非常稳跑33节点系统基本是秒级完成SDPT3在一些边角场景里会慢一些但小算例也完全够用。4.2 数据准备把配电网拓扑和负荷整理成标准格式建模的第一步是把网络参数整理成可供后续代码直接引用的数组。这里以33节点配电网测试系统为例需要准备这几份数据节点数 nb支路数 nl每条支路的首端节点 fbus(k)、末端节点 tbus(k)、电阻 r(k)、电抗 x(k)每个节点的有功负荷 pload 和无功负荷 qload根节点变电站电压幅值基准通常设为 1.0 pu。在Matlab代码里我习惯把支路参数和节点编号先填在一张表里比如% 支路数据格式: [首端节点 末端节点 电阻(pu) 电抗(pu)] branch [ 1 2 0.0922 0.0470 2 3 0.4930 0.2511 3 4 0.3660 0.1864 % ... 其余支路按实际拓扑填写 ]; fbus branch(:,1); tbus branch(:,2); r branch(:,3); x branch(:,4);这样做的目的是把数据和模型分离。以后换算例时只需改这张表不用动优化模型主体。数据单位统一采用标幺值这一步对数值稳定性非常重要后面会专门讲。4.3 变量声明哪些量交给求解器哪些量是固定的清晰定义变量是YALMIP建模最重要的一步。在我的实现里需要求解器决策的变量有这些V sdpvar(nb, 1); % 节点电压幅值平方 P sdpvar(nl, 1); % 支路首端有功功率 Q sdpvar(nl, 1); % 支路首端无功功率 l sdpvar(nl, 1); % 支路电流幅值平方 Pg sdpvar(nb, 1); % 节点分布式电源有功注入 Qg sdpvar(nb, 1); % 节点分布式电源无功注入注意这里直接优化 V 和 l而不再额外设置电压幅值和电流幅值的变量。目标函数和约束条件全部用这些平方变量来表达。如果某个节点没有接入分布式电源对应 Pg、Qg 直接置为 0 即可。4.4 核心约束的逐条建模首先是根节点电压约束。根节点是由变电站母线供电的参考节点电压幅值恒定在 1.0 puC [C, V(1) 1.0^2];然后是节点电压上下限约束。电压幅值下限设为 0.95 pu上限设为 1.05 pu换算成平方变量C [C, 0.95^2 V 1.05^2];支路电流的容量约束同样可以写成线性约束。例如电流幅值上限 1.5 puC [C, l 1.5^2];接下来是DistFlow的电压降落方程。对每一条支路 k首端节点 i fbus(k)末端节点 j tbus(k)方程写为C [C, V(j) V(i) - 2*(r(k)*P(k) x(k)*Q(k)) (r(k)^2x(k)^2)*l(k)];然后是节点功率平衡方程。这里我推荐用矩阵形式一次性写好比用循环拼接更紧凑。定义两个稀疏矩阵一个收集支路末端受端到节点的映射一个收集支路首端送端到节点的映射然后每个节点的平衡方程可以写成% In 矩阵: In(k, tbus(k)) 1, 表示支路k流入节点tbus(k) % Out 矩阵: Out(k, fbus(k)) 1, 表示支路k从节点fbus(k)流出 C [C, In * (P - r .* l) - Out * P Pg pload]; C [C, In * (Q - x .* l) - Out * Q Qg qload];注意这里的功率平衡方程对根节点同样成立。根节点的 Pg(1)、Qg(1) 会被求解器解释为从变电站流入配电网的功率不需要额外处理。4.5 二阶锥约束YALMIP里最容易写错的一行最关键的锥约束对应模型中的核心非凸项。对每条支路 k标准的YALMIP写法是for k 1:nl i fbus(k); C [C, cone([2*P(k); 2*Q(k); V(i) - l(k)], V(i) l(k))]; end这里cone(a, b)表示范数约束||a||₂ ≤ b。YALMIP会自动把这种约束识别为二阶锥约束。第一次写这段代码时要注意cone的第一个参数必须是一个列向量不能是行向量。我因为这个问题吃过亏——建模时看起来没问题一求解就报维度错误。另外有些旧版本YALMIP对cone的V第一个参数维度要求比较严格最好先把三个变量拼成一个[2*P(k); 2*Q(k); V(i)-l(k)]的列向量再传给cone。4.6 目标函数与求解调用网损最小是最常用的目标函数在SOCP框架下它正好是线性函数Objective sum(r .* l);如果目标是购电成本最小则可以直接把根节点有功注入作为目标% 假设电价为0.5元/kWh Objective 0.5 * Pg(1);然后调用求解器ops sdpsettings(solver, mosek, verbose, 1); optimize(C, Objective, ops); V_sol value(V); P_sol value(P); Q_sol value(Q); ...求解完成后利用value()提取各变量的数值结果再计算松弛间隙gap value(l) - (value(P).^2 value(Q).^2) ./ value(V); max_gap max(gap);如果max_gap在 1e-5 以下基本可以认为松弛是紧的SOCP的解可以直接作为原问题的全局最优解使用。5. 松弛紧性的实战检验SOCP解与非线性的差距到底有多大5.1 用同一个算例做公平对比理论说得再好也得看实际计算效果。我把SOCP方法和直接用非线性优化求解器求解原始非凸模型的结果放在一起对比。测试算例用的是我手头一个33节点辐射状配电网系统总负荷约3715 kW分布式电源接入在部分节点。SOCP模型用YALMIP建模、MOSEK求解原始非凸模型用MATLAB的非线性优化接口直接求解初值由潮流计算提供。值得一提的是肉眼能看到差异的还不只是目标值。SOCP求解过程全程没有初值参与完全由求解器自己探索跑出来的结果却相当稳定非线性求解器则需要我反复调整初值才能逼近同样的结果有时还会出现明显偏离物理意义的解。5.2 目标值、松弛间隙和耗时对比以下是我在该算例上得到的一组典型对比数据求解方式网损目标值(kW)最大松弛间隙(pu)求解耗时(秒)二阶锥松弛SOCP131.858.6e-70.38原始非凸非线性求解131.78—4.25目标值几乎相等SOCP的最大松弛间隙在1e-6量级说明松弛非常紧。两者的网损结果相差只有0.05%完全在工程允许误差范围内。计算时间的差距就明显了SOCP耗时不到0.4秒而非线性求解器即使有初值帮助也要4秒以上。如果做多时段优化、配电网重构、鲁棒优化这类需要反复调用OPF计算的任务这个速度差距会被放大数十倍。5.3 电压剖面的一致性也值得关注除了目标值我更关心电压剖面对比。把SOCP解出的节点电压幅值和非线性求解器结果画在同一条坐标轴上两条曲线几乎完全重合。SOCP解出的电压最低点出现在某个末端节点约0.945 pu非线性模型的结果也在此节点附近电压偏差不超过0.002 pu。这说明SOCP不只是在目标函数层面靠谱在约束层面给出的电压解、潮流分布也都是物理可信的。5.4 什么情况下SOCP会“掉链子”老实说SOCP并非在所有场景下都紧。我测试过的算例中有几种情况需要特别警惕一是目标函数对支路电流不敏感。比如目标函数中几乎不包含网损项而更多强调某种采购成本或储能充放电成本时求解器没有动力把 l 压到锥壁上松弛可能松弛得过宽解出来的结果不是原问题的可行解。二是电压约束过于紧张。当系统电压下限特别高导致电压越限约束成为主导约束时有时会把最优解“推”到锥的内部松弛不紧。三是网络为环网或多联络结构时。一旦配电网不再是严格辐射状SOCP模型就缺少了拓扑结构带来的紧性保障。第一种情况可以通过在目标函数里加一个很小的网损惩罚项来缓解这也是工程中常用的做法。第二种和第三种则需要谨慎对待不能盲目套用SOCP。6. 真实项目里最常见的坑从不可行到劣质解的一次次排查手记6.1 求解器报“不可行”时先查方向第一次跑33节点算例求解器给我报了一个“Infeasible problem”我当时第一反应是变量范围没设好查了半天也没头绪。后来发现问题出在支路功率平衡的方向上。DistFlow方程对方向非常敏感。在功率平衡方程里我写的是In * (P - r.*l) - Out * P Pg pload如果一不小心把In和Out搞反或者把P - r.*l写成了P r.*l求解器就会认为某条支路上出现了“凭空产生”的功率导致约束自相矛盾。排查这类问题的顺序我的经验是三步先看根节点电压约束是否合理再用一个已知潮流解测试功率平衡方程最后仔细核对fbus和tbus的方向是否与网络图一致。6.2 数值尺度问题命名的单位决定求解器的心情求解器对数值尺度相当敏感。如果把电压写成有名值比如11000 V把功率写成kW量级多个变量之间会出现10^4以上的量级差哪怕模型本身写对了求解器也容易因为数值病态而报错或收敛到低精度解。我一个同事的项目里出现过类似情况模型从“看起来没问题”到“求解精度差”折腾了两天。改成标幺值之后所有变量都落在0到100之间问题瞬间消失。所以建议所有变量统一使用标幺值电压基准取配电网额定电压功率基准取100 MVA或1 MVA。数据输入和结果分析时再进行单位的换算。6.3 松弛间隙过大时的“人工干预”在遇到目标函数对网损不敏感的场景时我试过在模型里加上一个很小的网损惩罚项比如epsilon 1e-4; Objective OriginalObjective epsilon * sum(r .* l);这个操作会轻微改变目标值但能有效把最优解“拉”回锥壁上从而让SOCP求出的解重新落回原问题的可行域内。要注意是惩罚系数不能太大否则会喧宾夺主影响原目标函数的最优性。我的操作习惯是先用 ε0 跑一遍如果松弛间隙过大再逐步增大 ε直到间隙降到1e-5以下。6.4 环网结构怎么办别硬套SOCP我之前也试过把SOCP直接套在一个带联络开关的环网算例上结果松弛间隙达到10^-2量级目标值和真实网损差得离谱。对于那些本身是环网、或者允许合环运行一段时间后再开断的网络SOCP不能直接套用。一个相对实用的替代思路是先把网络打开成辐射状把环网问题拆成多个辐射状场景分别求解再用枚举或启发式方法在这些场景里做外层寻优。这样至少能保住SOCP“求解快、结果稳”的核心优势。6.5 分布式光伏高渗透率场景关注电压越限与松弛的相互作用分布式光伏接入得越多配电网的净负荷波动越大某些时段甚至会出现从配电网向上级电网返送功率的情况。在这种反向潮流场景下节点电压容易被抬高电压上限约束可能被激活。我实测下来的结论是反向潮流本身并不会直接导致SOCP松弛失紧但一旦电压上限约束成为瓶颈松弛间隙可能会变大。解决办法是不要吝啬在目标函数里加一个电压偏移惩罚项或者在模型里精确建模逆变器的无功调节能力让求解器有更多自由度去平衡电压和松弛。过了几轮测试之后我对线性之外的承度认识又加深了一层——SOCP最大的魅力不在于“快”而在于“确定”。每次跑同一个算例它能给你同一个解这种确定性在分析多个运行场景、做方案比较时太重要了。后面如果想继续往深走我会把目标函数换成考虑分时电价下的配电网日内滚动优化这也是SOCP在多时段问题里最能发挥优势的方向。里面的初值敏感、计算时间这些问题你用了SOCP之后大概率都会有所体会。
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
