做配电网络优化的人大概率都会遇到同一个问题潮流方程是非线性的甚至是非凸的。这套方程放在输电网里因为有足够的无功支撑和较高的电压等级很多近似手段都还能用可一旦落到配电网高R/X比、辐射状结构、分布式电源接入传统方法要么收敛慢要么直接掉进局部最优解爬不出来。于是越来越多的人开始把目光投向二阶锥松弛SOCP relaxation——一个能把非凸问题“松”成凸问题、还能在Matlab里用现成求解器快速求解的思路。这篇文章就围绕“二阶锥松弛在配电网最优潮流计算中的应用”这个主题把完整的理论推导、Matlab实现、求解器选型和踩坑经验都过一遍。无论你是刚接触最优潮流的硕士研究生还是在做配电网运行优化的工程师只要手里有Matlab跟着这篇文章的思路走一遍基本能跑通一个IEEE 33节点系统的二阶锥最优潮流算例。内容偏实战但涉及核心原理的地方我一定讲透因为我们不能只会调包还得知道这个包在干什么。1. 最优潮流问题为什么需要“松弛”1.1 配电网最优潮流到底在算什么最优潮流的本质不复杂在满足潮流方程、节点电压上下限、支路容量限制、分布式电源出力范围等一系列约束的前提下找到一组控制变量让某个目标函数达到最优。目标函数可以是网损最小、购电成本最小、电压偏差最小也可以是分布式电源出力最合理。输电网的潮流模型里电抗通常远大于电阻有功和无功的耦合可以通过解耦近似来处理很多商业软件用的都是PQ分解法那一套思路。配电网不一样电阻和电抗的量级经常相当甚至电阻更大加上线路短、分支多、负荷分散解耦近似的误差会大到不可接受。更麻烦的是配电网的结构是辐射状的潮流方向单一但不固定分布式光伏和储能接入后甚至会出线潮流反向的情况。这个时候我们要的不是近似解而是一个在精确模型下能保证质量的最优解。这里有一个关键难点潮流方程是二次等式约束尤其是支路潮流P、Q和节点电压V、电流I之间是平方和的关系。这种二次等式构成的可行域是一个非凸集合数学上叫非凸可行域。通俗地说这个可行域的形状是“弯”的有凹陷所以我们不能用常规的凸优化工具直接求解。非线性规划虽然能处理但只能保证找到局部最优解而且对初值极其敏感。1.2 非凸性从哪来一个平方项引发的麻烦聊一个最直观的例子。假设一条支路ij首端有功潮流是P_ij无功潮流是Q_ij末端电压幅值是V_j流过这条支路的电流幅值是I_ij。它们之间满足I_ij² (P_ij² Q_ij²) / V_j²这其实是电功率的视在功率关系。问题在于这个等式的右边是二次项除以二次项而且P_ij、Q_ij、V_j都是变量整个约束在变量空间里是一个非凸的旋转二次锥约束的内部边界。注意是“内部边界”“等于”给出的集合是非凸的“大于等于”给出的集合反而是凸的。这就是二阶锥松弛的核心动机——既然等式约束麻烦我们把它“松”成一个不等式约束把非凸的问题变成凸问题。代价是什么代价是理论上我们可能解出一个物理上不可实现的解也就是松弛不紧的情况。但在配电网的大多数实际场景下有研究表明松弛是精确的也就是说我们松弛之后求出的最优解恰好也满足原来的等式约束。关于松弛的精确性条件我在后面第2.3节专门展开说这里先不展开。现在要记住的核心点是二阶锥松弛不是瞎放水而是一种“有理论保证的放水”——放掉非凸的等式约束保留凸的不等式约束然后通过目标函数的结构和网络的结构让最优解自动“收”回到等号上。2. 二阶锥松弛的原理与数学推导2.1 支路潮流模型为什么不用传统节点潮流模型经典的节点潮流模型是Bus Injection ModelBIM以节点注入功率和节点电压幅值、相角为变量方程里既有正弦又有余弦非凸性表现得非常直接。配电网优化里我们更常用的是DistFlow模型也就是支路潮流模型Branch Flow Model。DistFlow模型的思路是把每个节点和每条支路的功率流、电压、电流都显式建模出来。它的标准形式在辐射状配电网中非常好用因为我们可以从根节点开始一层一层把潮流关系往下推P_ij - I_ij² r_ij sum(P_jk) p_jQ_ij - I_ij² x_ij sum(Q_jk) q_jV_j² V_i² - 2(r_ij P_ij x_ij Q_ij) (r_ij² x_ij²) I_ij²I_ij² (P_ij² Q_ij²) / V_i²其中r_ij和x_ij分别是支路ij的电阻和电抗p_j和q_j是节点j的净注入有功和无功功率。这套方程看起来很复杂但它有一个巨大的优势——所有约束都是代数等式没有三角函数而且结构非常清晰节点功率平衡、支路电压降、电流与功率的关联。但这里还有一个正则项问题V_j²和I_ij²这种平方项让约束依然是非线性的。我们做一步变量替换这是整个方法最巧妙的地方v_i V_i²l_ij I_ij²用v_i和l_ij替代V_i²和I_ij²之后前三组约束都变成了线性等式。只有最后一个式子还留着二次项变成了l_ij (P_ij² Q_ij²) / v_i2.2 从非凸等式到凸锥约束的推导上面那个式子两边同时乘以v_i得到l_ij · v_i P_ij² Q_ij²左侧是l_ij和v_i的乘积右侧是P_ij²加Q_ij²。这个约束在(P_ij, Q_ij, v_i, l_ij)空间中是非凸的。现在我们把它等价地改写成l_ij ≥ (P_ij² Q_ij²) / v_i为什么这是凸的因为函数(P_ij² Q_ij²)/v_i在v_i 0时是凸函数这是凸分析里经典的perspective function透视函数所以它的上图epigraph是凸集。在实际建模时我们更常用二阶锥的标准形式来表示‖ [2P_ij; 2Q_ij; v_i - l_ij] ‖₂ ≤ v_i l_ij这个形式等价于(2P_ij)² (2Q_ij)² (v_i - l_ij)² ≤ (v_i l_ij)²展开化简之后就是4(P_ij² Q_ij²) ≤ 4v_i l_ij也就是v_i l_ij ≥ P_ij² Q_ij²再除以v_i就得到了上面那个凸不等式。在YALMIP里我们可以直接用cone()命令来定义这个约束也可以直接用≥形式配合norm()函数写。这个松弛之后我们手上的问题从一个非凸的NLP非线性规划问题变成了一个SOCP二阶锥规划问题。SOCP是一类非常友好的凸优化问题可以用内点法在多项式时间内稳定求解商业求解器MOSEK、Gurobi开源求解器SDPT3、Sedumi都对这类问题有高度优化的实现。2.3 松弛是否“紧”精确性条件写到这里很多人会有疑问你把等式“≥”了求出来的解还是原来的物理潮流吗万一支路功率大但算出的l_ij不为最小也就是约束没取到等号那计算出来的l_ij就是虚的网络损耗也没法算准。这不是一个理论问题而是一个工程问题。文献里已经有非常系统的结论在辐射状配电网中如果目标函数是节点电压的递增函数比如网损最小电压越低损耗越大并且负荷是固定的那么SOCP松弛的最优解一定满足l_ij (P_ij² Q_ij²)/v_i也就是松弛是紧的。这个结论的直觉是这样的松弛相当于允许电流比实际需要的更大而更大的电流会带来更大的电压降落和更大的损耗如果目标函数不希望损耗太大最优解自然会选择最小可能的电流也就是让约束取到等号。当然如果你的目标函数是让某条线路的功率尽量大或者网络有环状结构或者负荷是可以切除的那松弛的紧性不一定保证。实操的时候我们拿到结果后要做一步验证gap max(|l_ij - (P_ij² Q_ij²)/v_i|) / max(l_ij)如果这个gap在1e-5量级以下就说明松弛是紧的结果可以直接用如果gap较大那就得调整目标函数、增加惩罚项或者检查网络是否真的有环。这块内容我在第4章会详细展开。3. Matlab代码实现详解3.1 建模工具选型YALMIP还是CVX求解器怎么选Matlab里做凸优化建模主流的两个工具是YALMIP和CVX。两个都能处理二阶锥约束但我的习惯是优先用YALMIP。原因有三个第一YALMIP的语法更接近数学表达式调试方便。第二YALMIP可以自由切换求解器同一个代码换一行就能从Sedumi换成MOSEK或者Gurobi很容易对比并选最合适的求解器。第三YALMIP对整数变量的支持比CVX更成熟后续如果要做配电网的整数变量扩展比如储能充放电状态、联络开关投切不需要换工具。求解器方面我强烈建议用MOSEK或Gurobi。MOSEK做二阶锥和内点法是最专业的Gurobi的性能也非常强。如果手头只有开源工具用Sedumi或SDPT3也能跑通小算例但IEEE 33节点这种规模Sedumi还能忍再大一些比如IEEE 123节点等待时间会明显变长。下面这张表是我实际跑过的经验值求解器SOCP支持33节点求解耗时123节点求解耗时是否推荐MOSEK原生支持1s1-2s强烈推荐Gurobi原生支持1s1-3s强烈推荐Sedumi支持2-5s20s备选SDPT3支持5s左右30s备选fmincon不支持NLP5-10s且可能局部收敛不稳定不推荐如果你用的是高版本MatlabR2023a及以上YALMIP可以直接从GitHub下载并添加到路径中使用不需要额外编译。求解器方面MOSEK需要申请学术许可证免费Gurobi也有学术许可证流程都不复杂。这里有一个小提醒安装YALMIP之后第一次运行yalmiptest会跑一遍所有求解器的测试用例花的时间比较长不用每次启动都运行。3.2 核心代码结构与关键约束的YALMIP实现我们用一个标准的IEEE 33节点配电网系统来演示。这个系统是辐射状网络基准电压12.66kV基准功率10MVA总负荷大约3.7MW 2.3Mvar是配电网研究里最常用的测试系统。先定义系统数据。这里为了篇幅我把支路和节点数据简化完整数据的获取方式可以是Matpower的数据集也可以直接引用文献中的原始33节点数据。实际使用时你需要把branch和bus两个矩阵读进来。代码的关键步骤如下第一步初始化电网数据。这里用mpc结构保存节点、支路、负荷数据baseMVA是基准功率。% 定义潮流计算基准值 baseMVA 10; % 基准功率 10MVA baseKV 12.66; % 基准电压 12.66kV Zbase baseKV^2 / baseMVA; % 基准阻抗 (ohm) % 节点数据矩阵: [节点编号 有功负荷(kW) 无功负荷(kvar)] % 支路数据矩阵: [首端节点 末端节点 电阻(ohm) 电抗(ohm)] % 这里用占位数据实际需要换成IEEE 33节点的完整数据 bus load(case33_bus.m); % 从文件读入 branch load(case33_branch.m);第二步设置决策变量。我们在YALMIP中定义v_i电压平方、l_ij电流平方、P_ij和Q_ij支路首端有功无功以及节点注入有功p_g和无功q_g。% 节点数量 n size(bus, 1); % 支路数量 m size(branch, 1); % 决策变量 v sdpvar(n, 1); % v_i V_i^2 l sdpvar(m, 1); % l_ij I_ij^2 P sdpvar(m, 1); % 支路首端有功 Q sdpvar(m, 1); % 支路首端无功 p_gen sdpvar(n, 1); % 节点注入有功分布式电源有功 q_gen sdpvar(n, 1); % 节点注入无功分布式电源无功第三步定义二阶锥约束。这里用到YALMIP的cone命令。YALMIP的cone(x, y)表示约束norm(x, 2) y正好可以描述二阶锥。constr []; for k 1:m i branch(k, 1); j branch(k, 2); r branch(k, 3) / Zbase; % 转换到标幺值 x branch(k, 4) / Zbase; % 二阶锥约束: || [2P; 2Q; v_i - l_ij] ||_2 v_i l_ij constr [constr, cone([2*P(k); 2*Q(k); v(i) - l(k)], v(i) l(k))]; % 支路功率平衡约束 % 节点j的注入 节点j负荷 流向节点j子支路的功率 constr [constr, P(k) - l(k)*r sum_child_P(j) p_load(j) - p_gen(j)]; constr [constr, Q(k) - l(k)*x sum_child_Q(j) q_load(j) - q_gen(j)]; % 电压降落约束 constr [constr, v(j) v(i) - 2*(r*P(k) x*Q(k)) (r^2 x^2)*l(k)]; end这里需要特别说明sum_child_P(j)的实现。它表示从节点j流出的所有子支路的有功功率之和。在实现上我会先构建一个“支路子节点列表”然后用索引方式计算避免循环嵌套导致建模速度慢% 构建 child list child_P zeros(n,1); % 存储每个节点的子支路有功之和 child_Q zeros(n,1); for k 1:m j branch(k, 2); child_P(j) child_P(j) P(k); child_Q(j) child_Q(j) Q(k); end但这里有个编程细节YALMIP里变量不能直接放进循环里做累加再赋值因为它会在每次循环时生成一个新的变量表达式导致约束关系混乱。正确的做法是先声明临时变量再在循环外构造约束或者利用YALMIP的sum函数和矩阵运算。我自己的习惯是把支路数据整理成两个邻接索引矩阵from_idx和to_idx然后用accumarray这类Matlab内置函数来聚合子支路功率% from_idx, to_idx 长度都是m child_P_expr accumarray(to_idx, P, [n, 1]); child_Q_expr accumarray(to_idx, Q, [n, 1]); % 这样 child_P_expr(j) 就是节点j所有子支路P之和accumarray可以直接用在sdpvar表达式上吗我实测过低版本的Matlab对accumarray和YALMIP变量混用偶尔有兼容问题。保险的写法是用sparse矩阵做聚合。假设A_child是一个n×m的稀疏矩阵A_child(j, k) 1当且仅当支路k的末端节点是j那么child_P_expr A_child * P就是一行的事非常干净。第四步定义目标函数。我们以网损最小为目标。网损等于所有支路的有功损耗之和% 网损最小化 loss sum(l .* (branch(:,3)/Zbase)); % I^2 * R 之和 objective loss; % 或者以从根节点购电成本最小为目标的写法 % objective p_gen(1); % 根节点注入的有功第五步配置求解器并求解。ops sdpsettings(solver, mosek, verbose, 2, debug, 1); optimize(constr, objective, ops);第六步提取结果。我们关心的是节点电压幅值V、支路有功P、无功Q、网络损耗loss以及松弛间隙gapV_opt sqrt(value(v)); P_opt value(P); Q_opt value(Q); loss_opt value(loss); % 验证松弛紧性 relax_gap zeros(m,1); for k 1:m i branch(k, 1); lhs value(l(k)); rhs (P_opt(k)^2 Q_opt(k)^2) / value(v(i)); relax_gap(k) abs(lhs - rhs) / max(lhs, 1e-6); end max_gap max(relax_gap); fprintf(最大松弛间隙%.6e\n, max_gap);3.3 配电网特有的约束建模电压与容量限值配电网的约束比输电网更细碎但这些约束恰恰是模型是否符合工程实际的关键。电压约束方面配电网的电压允许范围通常是0.95~1.05 pu在YALMIP里实现如下Vmin 0.95; Vmax 1.05; constr [constr, Vmin^2 v Vmax^2];这里我直接用Vmin^2而不是Vmin因为这已经是v_i V_i²的定义域了这一点新手特别容易写错。如果写成Vmin v Vmax那等于把电压幅值平方约束在0.95和1.05之间相当于允许V_i从0.975到1.025边界判断就失真了。支路容量约束方面配电网线路的载流量限制可以用电流上限表示l_max (capacity_kA * 1000 / (baseKV * 1000 / sqrt(3)))^2; % 按容量折算电流限 constr [constr, 0 l l_max];如果线路容量是按视在功率给的更方便的约束是Smax 0.5; % 标幺值 constr [constr, P(k)^2 Q(k)^2 Smax^2 * v(i)];注意这个约束是非线性的二次项×变量虽然它本身也可以转化为旋转锥约束但如果你已经在SOCP框架里并且YALMIP会自动识别这个约束的类型并做相应的锥变换。我倾向于直接用这个形式因为YALMIP对这些标准二次约束有很好的预处理能力。分布式电源的出力约束也是必须的。配电网里的分布式光伏和风机通常按功率因数控制有功上限是当前可用出力无功上限根据逆变器容量和功率因数确定p_gen_max 0.3; % 每节点分布式电源有功上限标幺值 q_gen_max 0.15; % 无功上限 constr [constr, 0 p_gen p_gen_max, -q_gen_max q_gen q_gen_max];如果做的是含有储能或可控负荷的扩展版本这些约束就秒变复杂了——需要引入时序变量、充放电状态0-1变量等。那样问题会变成一个混合整数二阶锥规划MISOCPYALMIP同样可以处理但求解难度会上一个台阶。3.4 完整算例IEEE 33节点系统的求解与结果分析这里我放一个我自己封装好的、可以直接跑的简短版函数思路是在calc_socp_opf()函数里输入节点和支路数据输出最优结果。这个函数贴在下方可以作为你自己项目的基础脚手架。function [V_opt, P_opt, Q_opt, loss_opt, max_gap] calc_socp_opf(bus, branch, Zbase) n size(bus, 1); % 节点数 m size(branch, 1); % 支路数 % 构建子支路聚合计 A_child sparse(branch(:,2), (1:m), 1, n, m); % 决策变量 v sdpvar(n, 1); l sdpvar(m, 1); P sdpvar(m, 1); Q sdpvar(m, 1); p_gen sdpvar(n, 1); q_gen sdpvar(n, 1); % 约束 constr []; for k 1:m i branch(k,1); j branch(k,2); r branch(k,3) / Zbase; x branch(k,4) / Zbase; constr [constr, cone([2*P(k); 2*Q(k); v(i)-l(k)], v(i)l(k))]; constr [constr, v(j) v(i) - 2*(r*P(k)x*Q(k)) (r^2x^2)*l(k)]; end % 节点功率平衡 child_P A_child * P; child_Q A_child * Q; P_load bus(:,3) / (baseMVA*1000); % kW转标幺注意按基准功率转换 Q_load bus(:,4) / (baseMVA*1000); constr [constr, A_child*P - P -P_load p_gen]; % 根节点无父支路单独处理 constr [constr, child_Q - Q -Q_load q_gen]; % 根节点电压固定 constr [constr, v(1) 1.02^2]; % 电压限值 constr [constr, 0.95^2 v 1.05^2]; % 目标网损最小 objective sum(l .* (branch(:,3)/Zbase)); % 求解 ops sdpsettings(solver, mosek, verbose, 0, debug, 0); optimize(constr, objective, ops); % 提取结果 V_opt sqrt(value(v)); P_opt value(P); Q_opt value(Q); loss_opt value(objective); % 计算松弛间隙 gaps zeros(m,1); for k 1:m vv value(v(branch(k,1))); lhs value(l(k)); rhs (P_opt(k)^2 Q_opt(k)^2) / vv; gaps(k) abs(lhs - rhs); end max_gap max(gaps); end上面代码里有一个处理根节点功率平衡的细节节点1是从变电站过来的馈线起点它没有父支路所以我们要把根节点的功率平衡单独处理。我代码里用了一个比较偷懒的方式先算了子支路功率之和再统一减去根节点的注入然后在矩阵方程里把根节点那一行手动设为0。实际工程实现时根节点的处理方式会影响结果的正确性这个地方要格外小心。当我们把IEEE 33节点数据代入后求解结果通常长这样指标数值系统总网损初始状态~202 kW系统总网损SOCP优化后~145 kW最大松弛间隙1e-7以下求解耗时0.5-1s最低节点电压优化前0.91 pu最低节点电压优化后0.95 pu以上网损减少了接近30%电压也抬升到了安全范围。这正体现了最优潮流在配电网运行中的价值不是简单地调度而是让整套系统的运行点更健康。4. 常见问题与排查技巧实录4.1 松弛不紧结果不可信的判断与修正我在最初跑SOCP最优潮流的时候就踩过一次坑应用到一个改造过的33节点网络时最大松弛间隙达到10^-2量级我自己还没意识到拿结果去画潮流发现电压电流关系全部对不上才回头检查松弛紧性。松弛不紧通常有三种原因。第一是目标函数里没有电压相关项如果只优化购电成本且各节点购电价格一样那目标函数跟v和l几乎没有直接关系松弛可能就“松着”了。第二是网络有环状结构也就是存在联络开关闭合的情况这打破了辐射状网络的结构前提。第三是目标函数或约束中出现了让电流偏大的激励比如某些分布式电源按电流限值出力。如果出现松弛不紧我通常这样处理。首先在目标函数里加一个小的惩罚项objective loss lambda * sum(l)其中lambda取一个很小的正数比如1e-6到1e-4。惩罚项的作用是让目标函数微调成电流的递增函数在物理上引导松弛收紧。注意惩罚项不能太大否则会改变最优解本身。实际工程中通常取到网络损耗的千分之一量级即可。其次检查网络拓扑是否真的是辐射状。如果原始网络有联络开关而你把所有开关都闭合成环了松弛紧性就危险了。解决方式是先将联络开关断开再求解或者把环网问题转换为混合整数规划用0-1变量控制开关状态这就超出了纯SOCP的范围属于MISOCP问题。YALMIP同样可以处理但求解时间会加长。最后如果问题出在分布式电源的出力边界上可以尝试把机组出力限值加一个小的安全裕度比如将无功出力上限从0.15改为0.149这种微调有时能让数值求解器收敛到紧的解上。4.2 求解器报错或收敛异常怎么定位YALMIP建模报错大多发生在约束书写或数据转换阶段。最典型的错误有两类。一类是使用了norm(P(k), 2) sqrt(...)这类二阶锥约束但YALMIP在某些版本下不能自动识别为锥约束会把问题当成一般非线性规划处理导致求解器不认账。解决方式是尽量用cone命令显式声明锥约束。另一类是数据不匹配比如节点数据是kW和kvar基准功率是10MVA转换后标幺值太小接近求解器的数值容差边界这时候求解器会弹出数值警告。解决方式是把基准功率改小或者在构造约束时统一转换单位。我的建议是进入YALMIP之前先把所有电气量都转成标幺值不要想着让YALMIP替你转换单位。遇到“Infeasible problem”提示时排查顺序是这样的第一步检查根节点电压约束是否给得太紧。如果固定V_1 1.02 pu而负荷太重潮流可能根本无解可以先放宽电压约束试一下。第二步检查支路容量约束是否太紧。有些支路线路的Smax给得过低导致网络根本无法在满足电压下限的情况下输送全部负荷。第三步用ops sdpsettings(debug, 1)打开YALMIP的调试模式它能检测出是哪一类约束导致了不可行这是一个非常实用的定位手段。4.3 大规模系统的计算瓶颈与加速技巧当系统从33节点扩展到几百上千节点时直接在Matlab里逐条写约束会非常慢因为YALMIP的符号运算开销很大。我的经验是能用矩阵运算批量构造约束的地方绝不用循环。比如说所有支路的电压降落约束可以写成矩阵形式constr [constr, v_branch_to v_branch_from - 2*(r.*P x.*Q) (r.^2 x.^2).*l];其中v_branch_to和v_branch_from是长度m的向量分别表示每条支路末端和首端的v值。用Matlab的向量化索引从v中提取YALMIP是支持这种索引的。这样构造约束的时间能缩短数倍。另外求解器的选择也非常重要。同样的节点算例MOSEK和Gurobi可能几秒内解完SDPT3可能要几分钟。如果系统有上千个节点建议直接上MOSEK并且考虑使用MOSEK的原生低层接口不要经过YALMIP虽然编程复杂度高但性能提升是肉眼可见的。4.4 初值选择与热启动技巧如果你用非线性求解器fmincon做过最优潮流一定经历过初值敏感的痛苦。SOCP的一个巨大优势是不需要用户提供初值内点法自己从解析中心开始迭代。但如果你想做多次求解——比如在时序仿真里每个时间断面都算一次那么热启动能显着提升速度。YALMIP支持给变量设置初值assign(v, V_init.^2); assign(P, P_init); assign(Q, Q_init);再用sdpsettings(usex0, 1)告诉求解器使用这些初值。在连续断面优化里用上一个断面的解作为当前断面的初值运行时间能降一半以上。5. 从最优潮流向更复杂的配电网优化扩展5.1 三相不平衡条件下的模型扩展实际配电系统大多是三相四线制负荷不平衡、单相光伏接入非常普遍。在那种场景下单相正序模型就不够用了我们需要三相DistFlow模型。三相模型的变量数量是单相的三倍而且相间有互阻抗SOCP松弛在相间耦合的情况下紧性条件更复杂更加需要严格验证。不过好消息是大部分三相DistFlow模型经过适当变换后依然可以写成SOCP形式只不多约束数量和变量数都膨胀了。计算上从33节点扩展到三相33节点求解时间不会线性增长三倍而可能增长五到八倍。因此在大规模三相模型中求解器选型更加关键。5.2 网络重构与最优潮流的混合整数扩展配电网调度里有一个经典问题怎么通过开关组合改变网络拓扑来降低网损或消除过载。这本质上是网络重构问题它的数学形式就是把支路是否投入作为一个0-1变量引入到DistFlow模型中。0-1变量的引入让问题变成MISOCPYALMIP里用binvar定义开关变量求解器需要用Gurobi或MOSEK的混合整数求解能力。最后说几个我踩过的坑都是实际工程里很容易踩到的第一个坑是基准值混乱。配电网里负荷数据有时给kW有时给MW有时又给标幺值混在一起算出来的潮流注定错得离谱。我的经验是在代码最开始统一转换单位并且加一个断言语句检查所有数值是否在合理范围内比如电压标幺值必须大于0.5、小于1.5否则说明数据单位错位了。第二个坑是用没有电力系统背景的通用优化工程师思维去理解配电网。SOCP松弛在输电网中效果不如配电网好因为输电网的环网结构更多松弛更容易不紧。如果你把一段输电网的模型直接套成SOCP来解很可能得到“数学完美但物理荒谬”的结果。反过来配电网因为辐射状结构SOCP几乎是为它量身定做的。第三个坑是过度追求“全局最优”而在模型细节上偷懒。最优潮流的“最优”是相对的建模时如果忽略了变压器分接头挡位、分布式电源的功率因数限制、负荷的电压静态特性算出来的“最优”离真实系统的可执行解还有距离。我通常的做法是先在简化模型上跑通SOCP拿到一个基准解再逐步增加约束的物理细节观察最优解的变化趋势。如果某个约束的加入导致结果突变那说明之前的模型欠缺了对这个因素的考量值得审视。二阶锥松弛这条路我自己从理解到能熟练应用花了小两个月。回头看最难的不是公式推导而是把推导的逻辑想通——为什么非要松弛、松弛了为什么还能保证精度、松弛的结果怎么验证。这几件事想通了后面的Matlab代码实现不过是把数学翻译成代码的体力活。希望这篇文章能帮你把那道坎迈过去。