简介这是一份面向配电网重构与潮流计算研究的MATLAB源码针对支路潮流方程先采用二阶锥松弛再进一步使用多面体逼近从而将问题转化为混合整数线性规划MILP求解适合电气工程、自动化等相关专业学生开展毕设、课设或复现Jabr等经典文献方法。资源共3个文件以MATLAB脚本.m为主包含IEEE 33节点系统的MISOCP与MILP求解实现另有说明文档便于快速了解模型和运行方式整个压缩包仅8KB结构精简、便于直接阅读。该资源已有105人浏览学习。代码经运行验证可作为配电网最优重构、二阶锥松弛与多面体逼近对比研究的实用起点也可在现有基础上修改扩展实现其他功能下载后支持远程教学与答疑适合想尽快入门或深挖建模细节的读者。1. 配电网重构里为何要把潮流、SOCR、多面体逼近和 MILP 串成一条链路拿到配电网重构问题绝大多数人第一反应是启发式算法遗传、粒子群、模拟退火跑几百次也能给出一组开关组合。但这类方法有两个硬伤没有最优性下界无法证明当前解距离全局最优有多远遇到 N-1 校验、多时段储能、分布式电源随机性这类扩展时改造成本极高。于是工程上更常用的路径是把重构问题写成混合整数规划用 0-1 变量表示开关状态用潮流方程描述电气量关系再用大 M 或线性化技巧把非线性部分处理成可求解形式。这条链路里二阶锥松弛SOCR解决的是“潮流方程非凸”的问题多面体逼近解决的是“锥约束不是线性约束”的问题两步走完整个模型退化成纯 MILPMATLAB 自带的 intlinprog 就能求解。这篇文章就顺着这条链路把每一步的数学动机、约束写法、MATLAB 代码和参数调试顺序整理出来面向正在做主动配电网规划、重构或运行优化的工程师和研究生。2. 配电网重构的潮流基础DistFlow 方程为什么非凸二阶锥松弛怎样让它凸2.1 先选定潮流模型DistFlow 与 Branch Flow 的选择配电网重构里要反复修改拓扑每次修改都要重算潮流所以不能用需要节点导纳矩阵、依赖网络连通结构的传统牛顿法。业界和学术界更常用的是 DistFlow放射状支路潮流模型它把每条支路的功率、电压、电流分别设为变量只依赖支路两端状态。对辐射状网络的每条支路(i,j)定义变量P_ij为从 i 流向 j 的有功功率Q_ij为无功功率v_i为节点 i 的电压幅值平方l_ij为支路电流幅值平方。方程如下v_j v_i - 2(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*l_ij其中r_ij、x_ij是支路阻抗。看到没有这个方程本身是线性的。真正的非线性来自电流的定义式l_ij (P_ij^2 Q_ij^2) / v_i问题就在这里。l_ij 与 P、Q、v 是二次分式关系引入 0-1 开关变量 z_ij 之后整个可行域是非凸混合整数非线性规划MINLP普通求解器既没有全局最优性保证计算量也迅速膨胀。常见做法是用 Branch Flow 模型中的二阶锥松弛来处理这个等式把等号改成不等式(P_ij^2 Q_ij^2) v_i * l_ij再配合 v_i 0、l_ij 0 的约束它就等价于标准二阶锥形式|| [2*P_ij, 2*Q_ij, v_i - l_ij] ||_2 v_i l_ij这样处理之后原本非凸的潮流可行域被“撑”成了一个凸锥分布式电源接入、电容投切、储能充放电等约束都能直接往这个模型里加。2.2 松弛的紧性什么条件下 SOCR 不改变原问题的解二阶锥松弛扩大了可行域松弛解有可能落在原等式之外那就不是真实潮流。工程上要关心松弛是否“紧”exact如果求得的最优解满足 l_ij (P_ij^2 Q_ij^2)/v_i松弛就是精确的重构结果可以直接采用。根据常见的文献结论在主网提供足够电压支撑、网损占比不大、无逆向潮流等条件下配电网重构的 SOCR 解通常紧。但实际算例未必满足这些理论前提尤其是分布式光伏高渗透率时可能出现逆向潮流这时松弛质量会下降。我一般会在求解完成后做一步校验把 MILP 解代回原始潮流方程计算残差与目标函数偏差后面第 5 章会给出这段校验代码。2.3 MATLAB 里先跑一遍连续 SOCP得到最优性下界在嵌入 0-1 变量之前先把所有开关状态固定为闭合z1求解连续 SOCP。这一步有两个意义一是验证模型正确二是给后续 MILP 提供一个目标函数下界。MATLAB 中可以直接用 YALMIP 建模% 定义变量P、Q 为支路功率V 为节点电压平方L 为支路电流平方 P sdpvar(nBranch, 1); Q sdpvar(nBranch, 1); V sdpvar(nNode, 1); L sdpvar(nBranch, 1); % 约束DistFlow 电压方程 Constraints [V(fBus) - V(tBus) 2*(r.*P x.*Q) - (r.^2 x.^2).*L]; % 约束二阶锥松弛 for k 1:nBranch Constraints [Constraints, cone([2*P(k); 2*Q(k); V(fBus(k)) - L(k)], ... V(fBus(k)) L(k))]; end % 目标线损最小化 objective sum(r .* L); % 求解 ops sdpsettings(solver, gurobi, verbose, 1); optimize(Constraints, objective, ops);代码逻辑说明fBus、tBus是支路首末端节点编号r、x是支路阻抗向量cone(y, t)在 YALMIP 里表示||y||_2 t。上面这段求解出来的value(objective)是所有后续含开关方案的下界。这里有个建模框架的选型问题整理成表求解方案是否支持 0-1 变量是否支持 SOC 约束典型场景intlinprog支持不支持需多面体逼近纯 MILP完全依赖 MATLABconeprog不支持支持连续 SOCP 下界计算YALMIP Gurobi支持支持MI-SOCP效率最高YALMIP Mosek支持支持大算例且对数值稳定性要求高注意Gurobi、Mosek 能直接求解带锥约束的混合整数问题但这两个是商业求解器需要额外授权。如果项目环境只有 MATLAB 自带工具箱就必须走多面体逼近这也是标题里“进一步采用多面体逼近”的真正需求来源。3. 多面体逼近把二阶锥约束换成线性不等式组模型才进得了 MILP 框架3.1 为什么要多一步MI-SOCP 与 MILP 的生态差异第二章得到了一个带 0-1 变量和锥约束的混合整数二阶锥规划MI-SOCP。在学术论文里这已经可以求解但在实际工程环境里会遇到几个很现实的问题intlinprog 不支持锥约束商业求解器授权受限非线性锥约束在数值病态时容易产生不可行解或伪最优解。把锥用多面体逼近成线性约束之后整个模型就变成了标准 MILP求解器选型、许可证、调参经验全部复用成熟路线。多面体逼近的核心思想是用有限个超平面围出一个凸多面锥去逼近原来的旋转锥。多边形的边数 N 越大逼近精度越高但约束数量线性增加求解时间也随之上涨。实际算例里N 取 8 到 24 之间是性价比比较高的区间。3.2 内逼近与外逼近保守性与可行性的取舍对锥约束||[U; V]||_2 T用 N 个支撑半平面来近似U * cos(theta_k) V * sin(theta_k) c * T, k 1,2,...,N其中theta_k在 [-pi, pi] 上均匀取 N 个角度。当c 1时多面体是锥的外逼近得到的可行域比真实锥域大松弛解可能在原锥外界但不是真实可行解。当c cos(pi/N)时多面体位于原锥内部得到的可行域是真实可行域的子集解一定满足原锥约束但可能遗漏靠近锥边界的真实最优解。工程上两种做法都在用。求可行解优先的场合比如开关组合要直接下发给调度我倾向于用内逼近保证解的物理可行性求最优性下界或做松弛质量分析的场合可以用外逼近。还有一种两阶段法先用外逼近算出解再验证它是否满足原锥不满足就沿该点法向量加一条割平面重新求解迭代几次后逼近精度会非常高。下面给一个生成逼近系数的 MATLAB 函数function [A, b] polyConeApprox(N, mode) % 生成锥约束 ||[U; V]||_2 T 的多面体近似 % 输出形式: A * [U; V; T] b % N: 多边形边数mode: inner 内逼近 / outer 外逼近 theta linspace(0, 2*pi, N 1); theta theta(1:end-1); A zeros(N, 3); c 1; % 外逼近系数 if strcmp(mode, inner) c cos(pi / N); % 内逼近系数 end for k 1:N A(k, :) [cos(theta(k)), sin(theta(k)), -c]; end b zeros(N, 1); end参数说明每行不等式对应多边形的一条边。c 1时顶点恰好落在单位圆的切线方向上形成外切多边形c cos(pi/N)时把每条边向内平移使原锥的顶点到边的距离与圆心到边的距离对齐得到内接多边形。调用时把第 2.3 节里的cone([2*P; 2*Q; V_f - L], V_f L)替换为[A_p, b_p] polyConeApprox(16, inner); % 对每条支路加入 A_p * [2*P(k); 2*Q(k); V(fBus(k)) - L(k)] b_p3.3 逼近误差与边数 N 的关系内逼近的半径误差为rho 1 - cos(pi/N)外逼近的最大角度误差为pi/N。用一张表直观展示常见取值N内逼近半径误差约束条数/锥适用场景429.3%4快速可行性测试不推荐用于精度要求高的场景613.4%6拓扑调试阶段的粗略模型87.6%8重构方案的初步筛选123.4%12常规规划算例161.9%16推荐默认值误差与算力平衡较好240.85%24需要高精度目标函数值时才用需要说明的是半径误差并不直接等于目标函数误差因为目标函数对锥边界的敏感度取决于网损占系统总负荷的比例。但低网损场景下半径误差会近似传导向目标误差所以 N 直接决定你能把最优性 gap 压到多低。3.4 同时处理 0-1 变量与连续变量的乘积项配电网重构里还存在大量z_ij * x形式的乘积项比如支路功率在开关断开时应强制为 0用乘积约束直接表达就是非线性项。常见做法是引入中间变量 y 和 Big-M 约束做线性化y M * z y x - M * (1 - z) y x y -M * (1 - z)其中 M 取变量 x 实际物理范围的宽松上界。支路功率的 M 可以取该支路所带最大负荷的根号二倍电压平方变量取1.1^2 - 0.9^2再乘节点编号偏移量也可以关键是 M 不能取 1e9 这种无脑大数否则线性规划系数矩阵条件数恶化整数变量的分支定界过程会接近退化。4. 用 MATLAB 实现配电网重构 MILP从变量设计到求解器参数4.1 变量设计与数据组织配电网重构 MILP 的决策变量分成三类0-1 开关变量 z连续电气量变量P、Q、V、L以及线性化引入的辅助变量y_p、y_q。数据组织方面我通常用一个branch结构体存支路首末端、阻抗、初始状态一个load向量存各节点有功、无功负荷根节点单独编号。整模型的变量数量大约为6 * nBranch nNode大部分算例几千个变量intlinprog 都能在几十秒内完成。4.2 完整 MILP 模型用 4.3 的多面体约束替代锥约束这里给出一个可运行的 YALMIP 建模骨架核心差异是把第二章的cone替换成多面体线性约束并加上拓扑约束% 变量声明 z binvar(nBranch, 1); % 开关状态1为闭合 P sdpvar(nBranch, 1); % 支路有功 Q sdpvar(nBranch, 1); % 支路无功 V sdpvar(nNode, 1); % 节点电压平方 L sdpvar(nBranch, 1); % 支路电流平方 yP sdpvar(nBranch, 1); % 线性化辅助变量 z.*P yQ sdpvar(nBranch, 1); % 线性化辅助变量 z.*Q % 电压幅值约束 Constraints [0.95^2 V 1.05^2]; % DistFlow 电压方程与节点功率平衡 Constraints [Constraints, V(fBus) - V(tBus) ... 2*(r.*P x.*Q) - (r.^2 x.^2).*L]; % 多面体逼近替代二阶锥每条支路生成 N 条线性约束 [A_p, ~] polyConeApprox(16, inner); for k 1:nBranch U 2*yP(k); Vtmp 2*yQ(k); T V(fBus(k)) L(k); Constraints [Constraints, A_p * [U; Vtmp; T] zeros(16,1)]; end % 开关断开的支路功率为 0使用 Big-M 线性化 M 2 * sum(load) / nBranch; % 按单条支路最大可能载流量估算 Constraints [Constraints, yP M*z, yP P, yP P - M*(1-z), yP -M*(1-z)]; Constraints [Constraints, yQ M*z, yQ Q, yQ Q - M*(1-z), yQ -M*(1-z)]; % 辐射状拓扑单商品流约束见下方说明 f sdpvar(nBranch, 1); % 辅助流量变量 Constraints [Constraints, 0 f (nNode-1)*z]; nodeInc incidence(branch_table); % 节点支路关联矩阵 Constraints [Constraints, nodeInc(2:end,:) * f ones(nNode-1, 1)]; Constraints [Constraints, sum(z) nNode - 1]; % 目标网损最小 objective sum(r .* L); ops sdpsettings(solver, gurobi, verbose, 1); optimize(Constraints, objective, ops);代码逻辑说明incidence(branch_table)构造的是带方向的节点-支路关联矩阵乘上流量 f 后每个非根节点流出减流入等于 1这种单商品流约束配合sum(z) nNode - 1能同时保证网络连通且无环正好对应配电网的辐射状运行要求。第 6 行到第 8 行的 Big-M 线性化把“开关断开时功率为 0”和“开关闭合时功率等于实际值”两个语义合并进线性约束。4.3 intlinprog 与 YALMIP 的参数配置YALMIP 是建模层不是求解器。sdpsettings(solver, gurobi)会调用外部求解器如果不想装商业求解器可以把跨度放宽到intlinprog写法是sdpsettings(solver, intlinprog)。YALMIP 会自动把上面的线性约束和整数变量翻译成 intlinprog 需要的矩阵。MIPGap用sdpsettings(mipgap, 1e-4)控制分支定界提前终止条件。目标函数单位是 kW 到 MW 量级1e-4 的绝对 gap 对绝大多数规划问题已经足够。MaxTime重构算例规模不大时一般 300 秒足够。超过 600 秒还在不停分支通常说明 Big-M 选得过大或拓扑约束写错导致可行域过松。NodeLimitintlinprog 默认没有硬性限制但低配机器上建议设 5 万左右防止内存被分支树撑爆。常用参数整理如下参数建议值说明mipgap1e-4相对最优性间隙过小会拖长求解时间maxtime600超过后取当前最优可行解solverintlinprog / gurobi无商用许可可用前者多面体边数 N16误差 1.9%速度和精度平衡点Big-M支路容量 2~3 倍不要用 1e6 以上顺带提一句现在有 Codex 这类工具能像执行 Python 一样操作 MATLAB 任务把重构模型封装成函数后批量扫描 N-1 场景或负荷曲线是可行的。但注意模型本身的数值行为、松弛紧性验证这类判断还是得人来把关自动化工具更适合用于跑参数扫描而不是设计模型。4.4 求解流程建议连续下界 → MILP 可行解 → 校验我的习惯是分三步走先跑第 2.3 节的连续 SOCP 拿到目标下界再跑多面体 MILP 得到开关组合和网损最后把开关组合固定回代 DistFlow 方程校验电压和功率误差。这样做的好处是连续下界与 MILP 可行解之间的 gap 可以直接作为模型质量的量化指标。如果 gap 超过 5%优先怀疑多面体边数 N 太小其次检查径向拓扑约束是否写全。5. 验证松弛紧性与排错从目标函数回代到拓扑约束的全套检查方法5.1 校验 SOCR 松弛紧性的回代代码无论使用内逼近还是外逼近MILP 求出的解都需要回代原始潮流方程校验。这里的“原始方程”指的不是锥松弛后的不等式而是未松弛的等式约束% 取 MILP 求解结果 P_opt value(P); Q_opt value(Q); V_opt value(V); L_opt value(L); z_opt value(z); % 对闭合支路计算原始电流平方 l_orig (P_opt(active).^2 Q_opt(active).^2) ./ V_opt(fBus(active)); l_relax L_opt(active); % 偏差指标 maxErr max(abs(l_orig - l_relax) ./ l_orig); voltErr max(abs(V_opt - V_check)); % V_check 由 DistFlow 逐支路前推执行逻辑active是开关闭合的支路索引l_orig用功率和电压反推真实电流l_relax是优化模型的松弛变量。maxErr超过 1% 就说明二阶锥松弛在该算例上不紧这时要检查是不是出现了逆向潮流、电压越限或负荷过重的情况。对这类场景修正方案是把内逼近边数加大到 24或者改用两阶段割平面法在最优解处反复添加线性割平面逼近原锥。5.2 三个高频排错点在重构算例里的具体表现出现孤立节点目标函数能算出但有节点没有任何闭合支路连接。根因通常是单商品流约束里的注入量写错。检查incidence矩阵方向是否与fBus、tBus顺序一致以及根节点所在行有没有被排除在约束外。开关组合变化但网损几乎不变这种情况多半是 Big-M 取值远大于实际功率导致分支定界过程中功率可以“假装”在断开支路上流动目标函数对开关状态不敏感。把 M 压缩到支路载流量的 2 倍再跑。多面体逼近后原问题无可行解内逼近会缩小可行域这是正常代价。如果 N8 时无解而 N16 有解说明最优解落在锥边界附近提升边数即可如果 N24 仍无解优先检查拓扑约束而不是继续加大边数。调试顺序也有讲究先把 N 设成 4 跑通模型流程确认拓扑约束和 Big-M 线性化没有语法和维度问题再把边数提升到 8、16、24 各跑一遍画出网损随边数的收敛曲线这样最终上报的精度指标才有依据。本文还有配套的精品资源点击获取 配电网重构：DistFlow、二阶锥松弛、多面体逼近与MILP求解链路
简介这是一份面向配电网重构与潮流计算研究的MATLAB源码针对支路潮流方程先采用二阶锥松弛再进一步使用多面体逼近从而将问题转化为混合整数线性规划MILP求解适合电气工程、自动化等相关专业学生开展毕设、课设或复现Jabr等经典文献方法。资源共3个文件以MATLAB脚本.m为主包含IEEE 33节点系统的MISOCP与MILP求解实现另有说明文档便于快速了解模型和运行方式整个压缩包仅8KB结构精简、便于直接阅读。该资源已有105人浏览学习。代码经运行验证可作为配电网最优重构、二阶锥松弛与多面体逼近对比研究的实用起点也可在现有基础上修改扩展实现其他功能下载后支持远程教学与答疑适合想尽快入门或深挖建模细节的读者。1. 配电网重构里为何要把潮流、SOCR、多面体逼近和 MILP 串成一条链路拿到配电网重构问题绝大多数人第一反应是启发式算法遗传、粒子群、模拟退火跑几百次也能给出一组开关组合。但这类方法有两个硬伤没有最优性下界无法证明当前解距离全局最优有多远遇到 N-1 校验、多时段储能、分布式电源随机性这类扩展时改造成本极高。于是工程上更常用的路径是把重构问题写成混合整数规划用 0-1 变量表示开关状态用潮流方程描述电气量关系再用大 M 或线性化技巧把非线性部分处理成可求解形式。这条链路里二阶锥松弛SOCR解决的是“潮流方程非凸”的问题多面体逼近解决的是“锥约束不是线性约束”的问题两步走完整个模型退化成纯 MILPMATLAB 自带的 intlinprog 就能求解。这篇文章就顺着这条链路把每一步的数学动机、约束写法、MATLAB 代码和参数调试顺序整理出来面向正在做主动配电网规划、重构或运行优化的工程师和研究生。2. 配电网重构的潮流基础DistFlow 方程为什么非凸二阶锥松弛怎样让它凸2.1 先选定潮流模型DistFlow 与 Branch Flow 的选择配电网重构里要反复修改拓扑每次修改都要重算潮流所以不能用需要节点导纳矩阵、依赖网络连通结构的传统牛顿法。业界和学术界更常用的是 DistFlow放射状支路潮流模型它把每条支路的功率、电压、电流分别设为变量只依赖支路两端状态。对辐射状网络的每条支路(i,j)定义变量P_ij为从 i 流向 j 的有功功率Q_ij为无功功率v_i为节点 i 的电压幅值平方l_ij为支路电流幅值平方。方程如下v_j v_i - 2(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*l_ij其中r_ij、x_ij是支路阻抗。看到没有这个方程本身是线性的。真正的非线性来自电流的定义式l_ij (P_ij^2 Q_ij^2) / v_i问题就在这里。l_ij 与 P、Q、v 是二次分式关系引入 0-1 开关变量 z_ij 之后整个可行域是非凸混合整数非线性规划MINLP普通求解器既没有全局最优性保证计算量也迅速膨胀。常见做法是用 Branch Flow 模型中的二阶锥松弛来处理这个等式把等号改成不等式(P_ij^2 Q_ij^2) v_i * l_ij再配合 v_i 0、l_ij 0 的约束它就等价于标准二阶锥形式|| [2*P_ij, 2*Q_ij, v_i - l_ij] ||_2 v_i l_ij这样处理之后原本非凸的潮流可行域被“撑”成了一个凸锥分布式电源接入、电容投切、储能充放电等约束都能直接往这个模型里加。2.2 松弛的紧性什么条件下 SOCR 不改变原问题的解二阶锥松弛扩大了可行域松弛解有可能落在原等式之外那就不是真实潮流。工程上要关心松弛是否“紧”exact如果求得的最优解满足 l_ij (P_ij^2 Q_ij^2)/v_i松弛就是精确的重构结果可以直接采用。根据常见的文献结论在主网提供足够电压支撑、网损占比不大、无逆向潮流等条件下配电网重构的 SOCR 解通常紧。但实际算例未必满足这些理论前提尤其是分布式光伏高渗透率时可能出现逆向潮流这时松弛质量会下降。我一般会在求解完成后做一步校验把 MILP 解代回原始潮流方程计算残差与目标函数偏差后面第 5 章会给出这段校验代码。2.3 MATLAB 里先跑一遍连续 SOCP得到最优性下界在嵌入 0-1 变量之前先把所有开关状态固定为闭合z1求解连续 SOCP。这一步有两个意义一是验证模型正确二是给后续 MILP 提供一个目标函数下界。MATLAB 中可以直接用 YALMIP 建模% 定义变量P、Q 为支路功率V 为节点电压平方L 为支路电流平方 P sdpvar(nBranch, 1); Q sdpvar(nBranch, 1); V sdpvar(nNode, 1); L sdpvar(nBranch, 1); % 约束DistFlow 电压方程 Constraints [V(fBus) - V(tBus) 2*(r.*P x.*Q) - (r.^2 x.^2).*L]; % 约束二阶锥松弛 for k 1:nBranch Constraints [Constraints, cone([2*P(k); 2*Q(k); V(fBus(k)) - L(k)], ... V(fBus(k)) L(k))]; end % 目标线损最小化 objective sum(r .* L); % 求解 ops sdpsettings(solver, gurobi, verbose, 1); optimize(Constraints, objective, ops);代码逻辑说明fBus、tBus是支路首末端节点编号r、x是支路阻抗向量cone(y, t)在 YALMIP 里表示||y||_2 t。上面这段求解出来的value(objective)是所有后续含开关方案的下界。这里有个建模框架的选型问题整理成表求解方案是否支持 0-1 变量是否支持 SOC 约束典型场景intlinprog支持不支持需多面体逼近纯 MILP完全依赖 MATLABconeprog不支持支持连续 SOCP 下界计算YALMIP Gurobi支持支持MI-SOCP效率最高YALMIP Mosek支持支持大算例且对数值稳定性要求高注意Gurobi、Mosek 能直接求解带锥约束的混合整数问题但这两个是商业求解器需要额外授权。如果项目环境只有 MATLAB 自带工具箱就必须走多面体逼近这也是标题里“进一步采用多面体逼近”的真正需求来源。3. 多面体逼近把二阶锥约束换成线性不等式组模型才进得了 MILP 框架3.1 为什么要多一步MI-SOCP 与 MILP 的生态差异第二章得到了一个带 0-1 变量和锥约束的混合整数二阶锥规划MI-SOCP。在学术论文里这已经可以求解但在实际工程环境里会遇到几个很现实的问题intlinprog 不支持锥约束商业求解器授权受限非线性锥约束在数值病态时容易产生不可行解或伪最优解。把锥用多面体逼近成线性约束之后整个模型就变成了标准 MILP求解器选型、许可证、调参经验全部复用成熟路线。多面体逼近的核心思想是用有限个超平面围出一个凸多面锥去逼近原来的旋转锥。多边形的边数 N 越大逼近精度越高但约束数量线性增加求解时间也随之上涨。实际算例里N 取 8 到 24 之间是性价比比较高的区间。3.2 内逼近与外逼近保守性与可行性的取舍对锥约束||[U; V]||_2 T用 N 个支撑半平面来近似U * cos(theta_k) V * sin(theta_k) c * T, k 1,2,...,N其中theta_k在 [-pi, pi] 上均匀取 N 个角度。当c 1时多面体是锥的外逼近得到的可行域比真实锥域大松弛解可能在原锥外界但不是真实可行解。当c cos(pi/N)时多面体位于原锥内部得到的可行域是真实可行域的子集解一定满足原锥约束但可能遗漏靠近锥边界的真实最优解。工程上两种做法都在用。求可行解优先的场合比如开关组合要直接下发给调度我倾向于用内逼近保证解的物理可行性求最优性下界或做松弛质量分析的场合可以用外逼近。还有一种两阶段法先用外逼近算出解再验证它是否满足原锥不满足就沿该点法向量加一条割平面重新求解迭代几次后逼近精度会非常高。下面给一个生成逼近系数的 MATLAB 函数function [A, b] polyConeApprox(N, mode) % 生成锥约束 ||[U; V]||_2 T 的多面体近似 % 输出形式: A * [U; V; T] b % N: 多边形边数mode: inner 内逼近 / outer 外逼近 theta linspace(0, 2*pi, N 1); theta theta(1:end-1); A zeros(N, 3); c 1; % 外逼近系数 if strcmp(mode, inner) c cos(pi / N); % 内逼近系数 end for k 1:N A(k, :) [cos(theta(k)), sin(theta(k)), -c]; end b zeros(N, 1); end参数说明每行不等式对应多边形的一条边。c 1时顶点恰好落在单位圆的切线方向上形成外切多边形c cos(pi/N)时把每条边向内平移使原锥的顶点到边的距离与圆心到边的距离对齐得到内接多边形。调用时把第 2.3 节里的cone([2*P; 2*Q; V_f - L], V_f L)替换为[A_p, b_p] polyConeApprox(16, inner); % 对每条支路加入 A_p * [2*P(k); 2*Q(k); V(fBus(k)) - L(k)] b_p3.3 逼近误差与边数 N 的关系内逼近的半径误差为rho 1 - cos(pi/N)外逼近的最大角度误差为pi/N。用一张表直观展示常见取值N内逼近半径误差约束条数/锥适用场景429.3%4快速可行性测试不推荐用于精度要求高的场景613.4%6拓扑调试阶段的粗略模型87.6%8重构方案的初步筛选123.4%12常规规划算例161.9%16推荐默认值误差与算力平衡较好240.85%24需要高精度目标函数值时才用需要说明的是半径误差并不直接等于目标函数误差因为目标函数对锥边界的敏感度取决于网损占系统总负荷的比例。但低网损场景下半径误差会近似传导向目标误差所以 N 直接决定你能把最优性 gap 压到多低。3.4 同时处理 0-1 变量与连续变量的乘积项配电网重构里还存在大量z_ij * x形式的乘积项比如支路功率在开关断开时应强制为 0用乘积约束直接表达就是非线性项。常见做法是引入中间变量 y 和 Big-M 约束做线性化y M * z y x - M * (1 - z) y x y -M * (1 - z)其中 M 取变量 x 实际物理范围的宽松上界。支路功率的 M 可以取该支路所带最大负荷的根号二倍电压平方变量取1.1^2 - 0.9^2再乘节点编号偏移量也可以关键是 M 不能取 1e9 这种无脑大数否则线性规划系数矩阵条件数恶化整数变量的分支定界过程会接近退化。4. 用 MATLAB 实现配电网重构 MILP从变量设计到求解器参数4.1 变量设计与数据组织配电网重构 MILP 的决策变量分成三类0-1 开关变量 z连续电气量变量P、Q、V、L以及线性化引入的辅助变量y_p、y_q。数据组织方面我通常用一个branch结构体存支路首末端、阻抗、初始状态一个load向量存各节点有功、无功负荷根节点单独编号。整模型的变量数量大约为6 * nBranch nNode大部分算例几千个变量intlinprog 都能在几十秒内完成。4.2 完整 MILP 模型用 4.3 的多面体约束替代锥约束这里给出一个可运行的 YALMIP 建模骨架核心差异是把第二章的cone替换成多面体线性约束并加上拓扑约束% 变量声明 z binvar(nBranch, 1); % 开关状态1为闭合 P sdpvar(nBranch, 1); % 支路有功 Q sdpvar(nBranch, 1); % 支路无功 V sdpvar(nNode, 1); % 节点电压平方 L sdpvar(nBranch, 1); % 支路电流平方 yP sdpvar(nBranch, 1); % 线性化辅助变量 z.*P yQ sdpvar(nBranch, 1); % 线性化辅助变量 z.*Q % 电压幅值约束 Constraints [0.95^2 V 1.05^2]; % DistFlow 电压方程与节点功率平衡 Constraints [Constraints, V(fBus) - V(tBus) ... 2*(r.*P x.*Q) - (r.^2 x.^2).*L]; % 多面体逼近替代二阶锥每条支路生成 N 条线性约束 [A_p, ~] polyConeApprox(16, inner); for k 1:nBranch U 2*yP(k); Vtmp 2*yQ(k); T V(fBus(k)) L(k); Constraints [Constraints, A_p * [U; Vtmp; T] zeros(16,1)]; end % 开关断开的支路功率为 0使用 Big-M 线性化 M 2 * sum(load) / nBranch; % 按单条支路最大可能载流量估算 Constraints [Constraints, yP M*z, yP P, yP P - M*(1-z), yP -M*(1-z)]; Constraints [Constraints, yQ M*z, yQ Q, yQ Q - M*(1-z), yQ -M*(1-z)]; % 辐射状拓扑单商品流约束见下方说明 f sdpvar(nBranch, 1); % 辅助流量变量 Constraints [Constraints, 0 f (nNode-1)*z]; nodeInc incidence(branch_table); % 节点支路关联矩阵 Constraints [Constraints, nodeInc(2:end,:) * f ones(nNode-1, 1)]; Constraints [Constraints, sum(z) nNode - 1]; % 目标网损最小 objective sum(r .* L); ops sdpsettings(solver, gurobi, verbose, 1); optimize(Constraints, objective, ops);代码逻辑说明incidence(branch_table)构造的是带方向的节点-支路关联矩阵乘上流量 f 后每个非根节点流出减流入等于 1这种单商品流约束配合sum(z) nNode - 1能同时保证网络连通且无环正好对应配电网的辐射状运行要求。第 6 行到第 8 行的 Big-M 线性化把“开关断开时功率为 0”和“开关闭合时功率等于实际值”两个语义合并进线性约束。4.3 intlinprog 与 YALMIP 的参数配置YALMIP 是建模层不是求解器。sdpsettings(solver, gurobi)会调用外部求解器如果不想装商业求解器可以把跨度放宽到intlinprog写法是sdpsettings(solver, intlinprog)。YALMIP 会自动把上面的线性约束和整数变量翻译成 intlinprog 需要的矩阵。MIPGap用sdpsettings(mipgap, 1e-4)控制分支定界提前终止条件。目标函数单位是 kW 到 MW 量级1e-4 的绝对 gap 对绝大多数规划问题已经足够。MaxTime重构算例规模不大时一般 300 秒足够。超过 600 秒还在不停分支通常说明 Big-M 选得过大或拓扑约束写错导致可行域过松。NodeLimitintlinprog 默认没有硬性限制但低配机器上建议设 5 万左右防止内存被分支树撑爆。常用参数整理如下参数建议值说明mipgap1e-4相对最优性间隙过小会拖长求解时间maxtime600超过后取当前最优可行解solverintlinprog / gurobi无商用许可可用前者多面体边数 N16误差 1.9%速度和精度平衡点Big-M支路容量 2~3 倍不要用 1e6 以上顺带提一句现在有 Codex 这类工具能像执行 Python 一样操作 MATLAB 任务把重构模型封装成函数后批量扫描 N-1 场景或负荷曲线是可行的。但注意模型本身的数值行为、松弛紧性验证这类判断还是得人来把关自动化工具更适合用于跑参数扫描而不是设计模型。4.4 求解流程建议连续下界 → MILP 可行解 → 校验我的习惯是分三步走先跑第 2.3 节的连续 SOCP 拿到目标下界再跑多面体 MILP 得到开关组合和网损最后把开关组合固定回代 DistFlow 方程校验电压和功率误差。这样做的好处是连续下界与 MILP 可行解之间的 gap 可以直接作为模型质量的量化指标。如果 gap 超过 5%优先怀疑多面体边数 N 太小其次检查径向拓扑约束是否写全。5. 验证松弛紧性与排错从目标函数回代到拓扑约束的全套检查方法5.1 校验 SOCR 松弛紧性的回代代码无论使用内逼近还是外逼近MILP 求出的解都需要回代原始潮流方程校验。这里的“原始方程”指的不是锥松弛后的不等式而是未松弛的等式约束% 取 MILP 求解结果 P_opt value(P); Q_opt value(Q); V_opt value(V); L_opt value(L); z_opt value(z); % 对闭合支路计算原始电流平方 l_orig (P_opt(active).^2 Q_opt(active).^2) ./ V_opt(fBus(active)); l_relax L_opt(active); % 偏差指标 maxErr max(abs(l_orig - l_relax) ./ l_orig); voltErr max(abs(V_opt - V_check)); % V_check 由 DistFlow 逐支路前推执行逻辑active是开关闭合的支路索引l_orig用功率和电压反推真实电流l_relax是优化模型的松弛变量。maxErr超过 1% 就说明二阶锥松弛在该算例上不紧这时要检查是不是出现了逆向潮流、电压越限或负荷过重的情况。对这类场景修正方案是把内逼近边数加大到 24或者改用两阶段割平面法在最优解处反复添加线性割平面逼近原锥。5.2 三个高频排错点在重构算例里的具体表现出现孤立节点目标函数能算出但有节点没有任何闭合支路连接。根因通常是单商品流约束里的注入量写错。检查incidence矩阵方向是否与fBus、tBus顺序一致以及根节点所在行有没有被排除在约束外。开关组合变化但网损几乎不变这种情况多半是 Big-M 取值远大于实际功率导致分支定界过程中功率可以“假装”在断开支路上流动目标函数对开关状态不敏感。把 M 压缩到支路载流量的 2 倍再跑。多面体逼近后原问题无可行解内逼近会缩小可行域这是正常代价。如果 N8 时无解而 N16 有解说明最优解落在锥边界附近提升边数即可如果 N24 仍无解优先检查拓扑约束而不是继续加大边数。调试顺序也有讲究先把 N 设成 4 跑通模型流程确认拓扑约束和 Big-M 线性化没有语法和维度问题再把边数提升到 8、16、24 各跑一遍画出网损随边数的收敛曲线这样最终上报的精度指标才有依据。本文还有配套的精品资源点击获取