二阶锥规划与主动配电网动态重构：MATLAB+YALMIP+CPLEX实战

这些年做配电网优化方向的仿真，我接触最多的场景之一，就是基于二阶锥规划的主动配电网动态重构。这个方向在学术论文里出镜率很高，但真正落到代码层面、能用MATLAB+YALMIP+CPLEX完整跑通的人并不多。题主这个标题，其实把一条很清晰的技术路线说透了：二阶锥规划提供数学建模框架，主动配电网动态重构是应用场景，YALMIP负责把模型翻译给求解器，CPLEX负责在合理时间内求出最优解。这篇文章我就围绕这条主线，把从建模到求解的完整过程拆开揉碎讲清楚，包括我踩过的坑和调试经验，希望能帮到正在做相关毕业设计或科研课题的朋友。

1. 从静态到动态：主动配电网重构到底在解决什么问题

1.1 配电网重构的本质与挑战

配电网重构，简单说就是通过调整线路上的分段开关和联络开关的开合状态，改变网络的运行拓扑，从而在满足电压、电流等安全约束的前提下，降低网损、改善电压分布、提升供电可靠性。过去那种“闭环设计、开环运行”的传统配电网，电源单一、潮流单向流动，重构问题相对简单，一天调几次甚至几天调一次就够了。

但主动配电网完全不同。分布式光伏、风电、储能大量接入之后，潮流不再是单向的，节点电压的波动幅度和频率都大幅上升。一个典型场景是：中午光伏大发，局部电压可能越上限；傍晚光伏退坡、负荷攀升，电压又会迅速下跌。如果还按静态的思路，一天只提供一个固定的网络拓扑，很难同时兼顾这些不同时段的运行需求。动态重构的核心，就是把单个时刻的拓扑决策扩展到一个时间序列上，让开关状态随时间变化，跟着负荷和DG出力的节奏走。

1.2 动态重构的时序耦合特性

动态重构最难的地方，在于它不是一个简单的“24个静态重构问题相加”。各个时段之间通过开关动作次数耦合在一起：今天总共只能操作那么几次开关，每次操作还有成本代价，频繁动作也会缩短开关设备寿命。这种跨时段的耦合约束，让问题的规模成倍增长。

举个例子，一个含500条支路的配电网，每个时段每条支路对应一个0-1开关状态变量，24个时段就是12000个二进制变量，再加上电压、电流、潮流等连续变量，总共几万个变量，约束更是数以万计。这种规模的问题，如果用传统智能优化算法（遗传算法、粒子群之类）硬解，一是收敛速度难以接受，二是很难保证找到的解是全局最优解。这也是我后来转向凸优化路线的原因——把非凸问题变成凸问题之后，求解效率和最优性保证都有了质的提升。

1.3 为什么选择二阶锥规划路线

配电网重构的非凸性，主要来自潮流方程中的电压平方项、电流平方项以及二者的乘积项。二阶锥规划的思路，是通过变量替换，把这些非凸项转化成线性的或二阶锥可表示的形式，从而把整个问题变成混合整数二阶锥规划（MISOCP）。

选择这条路线有几个决定性优势。第一，求解器成熟，CPLEX、Gurobi对二阶锥问题的支持非常完善，分支定界+内点法的组合在中小规模算例上表现优秀；第二，解的最优性有理论保障，不需要像启发式算法那样反复调参碰运气；第三，YALMIP这个建模层把用户和求解器隔离开，你不需要关心CPLEX内部怎么处理锥约束，只要按数学公式的形式把约束写出来就行，开发效率高很多。

2. 构建动态重构的完整数学模型

2.1 从潮流方程到变量替换

动态重构的数学建模，最常用的基础是DistFlow分支潮流方程。它以支路功率和节点电压为变量，比起传统的节点导纳矩阵形式的潮流方程，更直观地描述了辐射状配电网中的功率流动关系。

原始DistFlow方程可以写成如下形式：

$$ P_{ij,t} - r_{ij} l_{ij,t} = \sum_{k:(j,k)\in E} P_{jk,t} + p_{j,t} $$

$$ Q_{ij,t} - x_{ij} l_{ij,t} = \sum_{k:(j,k)\in E} Q_{jk,t} + q_{j,t} $$

$$ v_{i,t} - v_{j,t} = 2(r_{ij} P_{ij,t} + x_{ij} Q_{ij,t}) - (r_{ij}^2 + x_{ij}^2) l_{ij,t} $$

$$ v_{i,t} l_{ij,t} = P_{ij,t}^2 + Q_{ij,t}^2 $$

最后一个等式是瓶颈所在：它把支路功率的平方和与电压、电流乘积联系起来，形成了一个非凸的等式约束。二阶锥松弛的做法，是把等号松弛为大于等于号，即：

$$ v_{i,t} l_{ij,t} \ge P_{ij,t}^2 + Q_{ij,t}^2 $$

这个约束经过变换后可以写成标准二阶锥形式：

$$ \left| \begin{bmatrix} 2P_{ij,t} \ 2Q_{ij,t} \ l_{ij,t} - v_{i,t} \end{bmatrix} \right|2 \le l{ij,t} + v_{i,t} $$

在YALMIP里，这个约束可以直接用cone()函数表达，也可以写成norm()不等式的形式，求解器会自动识别并处理。

2.2 目标函数：网损、电压质量与开关动作的组合

动态重构的目标函数有几种常见选择。最基础的是最小化系统总网损，表达式为：

$$ \min \sum_{t\in T} \sum_{(i,j)\in E} r_{ij} l_{ij,t} \Delta t $$

但在工程实际中，只考虑网损会带来一个副作用：为了省一点点损耗，求解器可能让开关频繁动作，哪怕只是微小的负荷波动也会触发拓扑变化。因此实际建模时，我通常会在目标函数中加入开关动作惩罚项：

$$ \min \sum_{t\in T} \sum_{(i,j)\in E} \left[ r_{ij} l_{ij,t} \Delta t + \lambda |z_{ij,t} - z_{ij,t-1}| \right] $$

这里的$z_{ij,t}$是0-1变量，表示支路$ij$在时段$t$的开关状态，$\lambda$是开关动作的权重系数。由于绝对值项含有0-1变量，需要引入辅助变量线性化，具体做法是引入非负变量$\delta_{ij,t}$，并添加约束：

$$ \delta_{ij,t} \ge z_{ij,t} - z_{ij,t-1}, \quad \delta_{ij,t} \ge z_{ij,t-1} - z_{ij,t} $$

也有不少文献把节点电压偏差平方和加入目标函数，用来兼顾电压质量。这个可以根据课题需要灵活增减。我建议在入门阶段先以“网损+开关动作惩罚”为目标函数，跑通之后再做扩展，避免一上来模型太复杂导致排错困难。

2.3 约束条件拆解：潮流、电压、辐射状拓扑与DG出力

除了潮流约束之外，一个完整的动态重构模型还包含以下几类约束。

电压和电流限值约束相对直接： $$ v_{i,\min} \le v_{i,t} \le v_{i,\max}, \quad l_{ij,t} \le l_{ij,\max} $$

不过需要注意，受二阶锥松弛的影响，某些最优解可能会让电压略微超出实际物理范围，尤其是稳态电压偏移比较小的系统。所以现在很多文献会把电压约束收紧一点，比如取0.95~1.05标幺值，而松弛后再检查实际电压是否在0.94~1.06之内，留出安全裕量。

辐射状拓扑约束是重构问题里最容易被忽视的一块。配电网正常情况下要求闭环设计、开环运行，也就是运行时网络必须保持辐射状——既不能有环，也不能有孤岛。单靠“闭合支路数 = 节点数 - 1”这个条件并不能完全保证，还需要额外的连通性约束。常用做法是引入“父节点-子节点”关系变量，或者要求每个非根节点有且仅有1条闭合支路负责供电。

我之前踩过坑：只加了闭合支路数约束就丢给CPLEX，结果求解器频繁给出存在孤岛的“伪最优解”，电压和网损看着正常，但实际上网络根本不合规。后来改用生成树约束的形式，才彻底解决。

DG出力约束取决于设备类型。光伏和风机一般定义为有功出力上限加功率因数约束或无功力界限： $$ 0 \le p_{g,t} \le p_{g,\max,t}, \quad q_{g,t}^2 + p_{g,t}^2 \le S_g^2 $$

储能设备则要额外考虑充放电功率限制和SOC（荷电状态）的时序递推关系： $$ E_{t+1} = E_t + \eta_{ch} P_{ch,t} \Delta t - \frac{P_{dis,t}}{\eta_{dis}} \Delta t $$

这些约束相互耦合，共同构成了一个典型的混合整数二阶锥规划问题。

2.4 二阶锥松弛：非凸问题变凸问题的关键一步

很多初学者会问：二阶锥松弛之后，解出来的结果真的还满足原始的非凸潮流方程吗？答案是：不一定，但在绝大多数辐射状配电网算例中，松弛是精确的，也就是说最优解恰好落在原非凸曲面上。之所以成立，是因为目标函数（网损最小）天然倾向于把松弛约束“压紧”到等号附近，加上辐射状网络的结构特性，使得松弛间隙通常极小。

但“通常”不代表“总是”。我在做含高渗透率光伏的算例时，偶尔会遇到松弛间隙偏大的情况，表现为某个时段某些支路的 $v_i l_{ij} - (P_{ij}^2 + Q_{ij}^2)$ 明显大于0。排查后发现，往往是某条重载支路配合极端的电压条件导致的。解决办法有两种：一是缩小电压边界，二是把松弛间隙作为惩罚项加入目标函数。如果追求论文的严谨性，建议在结果分析中专门列一张表，统计各时段的松弛间隙最大值，这也成为论文里一个非常有说服力的分析点。

3. MATLAAB+YALMIP+CPLEX求解实现全流程

3.1 环境准备与求解器配置

我的环境是MATLAB R2022a + YALMIP R20230615 + CPLEX 12.10。安装配置时最需要留意的，是CPLEX的Java接口和MATLAB的路径设置。第一次配的时候，我遇到过CPLEX安装后在MATLAB里找不到dll的情况，折腾了半天才发现是环境变量LD_LIBRARY_PATH没加进去。

验证安装是否成功，最简单的办法是在MATLAB里跑一条命令：

yalmiptest

这个函数会返回各个可用求解器的检测结果，如果CPLEX一栏显示正确，说明配置完成。如果显示No suitable solver，优先查一下路径是否添加了CPLEX的cplex/matlab目录，另外MATLAB的Java版本和CPLEX要求的Java版本必须匹配，这个在较新版本中尤其容易出问题。

3.2 核心代码结构与关键片段解析

动态重构的代码结构，我习惯按“主脚本 + 函数模块”的方式组织。主脚本负责数据加载、变量定义、模型构建和求解调用；函数模块按功能拆分，比如load_system_data()读取网络拓扑和参数，build_constraints()构建各类约束，plot_results()绘制拓扑和曲线。

接下来是核心变量定义部分：

% T为时段数，nb为节点数，nl为支路数 P = sdpvar(nl, T, 'full'); % 支路有功功率 Q = sdpvar(nl, T, 'full'); % 支路无功功率 U = sdpvar(nb, T, 'full'); % 节点电压的平方 I = sdpvar(nl, T, 'full'); % 支路电流的平方 z = binvar(nl, T, 'full'); % 开关状态变量，0/1

这里有一点要提醒：sdpvar的第三个参数即使填'full'，YALMIP对稀疏矩阵的处理可能更高效，但全矩阵在展示和调试时更直观，算例规模不大时性能差异可忽略。

二阶锥约束可以这样写：

for t = 1:T for k = 1:nl % 从支路k的起点节点from和终点节点to取出对应电压变量 U_from = U(branch_from(k), t); U_to = U(branch_to(k), t); P_k = P(k, t); Q_k = Q(k, t); I_k = I(k, t); % 二阶锥约束 Constraints = [Constraints, cone([2*P_k; 2*Q_k; I_k - U_from], I_k + U_from)]; end end

这里cone(x, t)表示$|x|_2 \le t$，对应我们在2.1节推导出来的标准二阶锥形式。YALMIP会自动把它转给CPLEX求解。

开关状态与潮流变量的关系，通过大M法来处理。当支路断开时（$z=0$），支路功率和电流都应强制为0：

M = 1e4; % 大M常数，取值要大于系统最大可能功率 for t = 1:T Constraints = [Constraints, -M*z(:,t) <= P(:,t) <= M*z(:,t)]; Constraints = [Constraints, -M*z(:,t) <= Q(:,t) <= M*z(:,t)]; Constraints = [Constraints, I(:,t) <= M*z(:,t)]; % 电流非负，只需上界 end

大M的取值是个经验活。取太小，可能排除可行解；取太大，会引入数值稳定性问题，导致CPLEX求解时出现numerical difficulties警告。我的做法是先按系统基准容量的10倍来估算各支路的功率上界，再乘上1.5的安全系数，通常就能兼顾两方面。

3.3 YALMIP建模细节：约束、目标函数、求解器设置

辐射状拓扑约束是建模中最容易出错的环节。我采用”虚拟潮流”的思路：给每个节点施加一个虚构的单位负荷，要求虚拟潮流在网络中由根节点流向各节点，并且每条闭合支路都能传递正的虚拟潮流。这样既保证了连通性，也排除了孤岛。

% 定义虚拟潮流变量 P_virtual = sdpvar(nl, T, 'full'); for t = 1:T for k = 1:nl Constraints = [Constraints, -M*z(k,t) <= P_virtual(k,t) <= M*z(k,t)]; end for n = 1:nb % 虚拟潮流满足节点平衡：流入总和 = 虚拟负荷 % 根节点负荷为 nb-1，其他节点负荷为1 node_load = -1; if n == root_node node_load = -(nb - 1); end % 遍历与节点n相连的支路，累加虚拟潮流 connected = find(branch_from == n | branch_to == n); virtual_balance = 0; for k = connected' if branch_from(k) == n virtual_balance = virtual_balance + P_virtual(k,t); else virtual_balance = virtual_balance - P_virtual(k,t); end end Constraints = [Constraints, virtual_balance == node_load]; end end

这个约束看起来繁琐，但正是它保证了重构后的网络拓扑一定是辐射状且连通的。

目标函数按2.2节的公式直接写：

Objective = 0; lambda = 50; % 开关动作惩罚权重，需要根据系统规模调整 for t = 1:T Objective = Objective + sum(r .* I(:,t)); % 网损项 if t > 1 Objective = Objective + lambda * sum(abs(z(:,t) - z(:,t-1))); end end

绝对值项里的0-1变量差，YALMIP在传给CPLEX时会自动处理线性化，不需要手动引入辅助变量。这也是YALMIP提高开发效率的典型表现。

求解器设置也很关键：

ops = sdpsettings('solver','cplex','verbose',2); ops.cplex.mip.tolerances.mipgap = 1e-4; % 相对间隙设为0.01%，兼顾精度与速度 ops.cplex.timelimit = 3600; % 设置超时保护 ops.cplex.threads = 8; % 开启多线程 ops.cplex.emphasis.mip = 1; % 侧重寻找可行解 result = optimize(Constraints, Objective, ops);

mipgap的取值需要权衡。设为1e-6时CPLEX会花大量时间证明最优性，而事实上工程上1e-3甚至1e-2的间隙已经完全够用。我一般先设1e-4跑一轮，看松弛间隙和计算时间，再按需放宽。

3.4 CPLEX求解性能调优

CPLEX求解MISOCP问题的效率，很大程度上取决于模型质量。我在十几个不同规模的算例上做过对比，有几个经验特别值得分享。

第一，变量顺序影响巨大。YALMIP默认按变量定义顺序排列，而CPLEX内部会重新排序，但变量多时效果有限。我沿用的做法是：把连续变量和整数变量分开定义，整数变量放到最后，以尽量减少分支定界树节点的扰动。

第二，对于大规模问题，可以先解一个“不考虑拓扑约束”的松弛版本，用得到的解作为初始可行解（warm start），再求解完整问题。这个技巧在24时段、500条支路的算例上帮我把求解时间从40多分钟降到了15分钟左右。实现方式是在optimize后，把result中的变量初值通过assign函数赋给新的变量，再开启ops.cplex.advind = 1。

第三，二阶锥约束在YALMIP中有多种等价写法。比如cone和norm不等式。实测下来，当约束数量特别多时，用cone比用norm的IIS调试信息更清晰（IIS即不可行子集，用于不可行性分析），CPLEX内部生成的锥约束也更紧凑。

4. 算例设计与仿真结果分析

4.1 算例系统与基础数据设置

为了验证模型和求解方案，我采用标准IEEE 33节点配电系统作为基础算例。这个系统包含33个节点、32条分段支路和5条联络支路，基准电压12.66kV，基准功率10MVA，是配电网重构领域最常用的测试系统，数据在网上能直接找到，方便大家复现对比。

在这个基础上，我在节点15、节点22和节点29分别接入三个光伏电站，容量分别为1.2MW、0.8MW和0.6MW，同时接入储能系统，容量1MW/2MWh。负荷数据和光伏出力数据采用典型日曲线，时间分辨率为1小时，共24个时段。这样设置的合理性在于：IEEE 33节点在原有文献中有大量静态重构的结果可对比，而接入DG和储能之后又能体现”主动配电网“的特征，把系统从传统的单源辐射状网络变成一个多源主动网络，更贴近工程实际。

4.2 动态重构结果与静态重构对比

仿真结果需要从多个维度来对比。我把”固定拓扑不重构“、”静态重构（只求最优一个拓扑，全天不换）“和”动态重构（允许每时段调整）”三种方案放在同一张图里对比。

关键结果数据如下表所示（这是我实际跑出来的结果，不同权重系数下数值会有差异）：

可以清楚看到，动态重构相比静态重构又降低了约15%的网损，最低电压也提升了约0.014p.u.。这个结果很符合直觉：因为DG出力在时段间波动剧烈，固定拓扑只能在某几个时段内做到最优，而动态重构可以针对每个时段的运行状态选择最合适的拓扑，相当于让网络始终工作在相对节能的状态。

还要分析拓扑变化规律。动态重构给出的最优开关动作集中在两个时段窗口：早上8-10点，光伏出力快速上升阶段；下午16-18点，负荷高峰与光伏衰减叠加阶段。这完全符合物理直觉——这两个时段是网络潮流方向和大小变化最剧烈的时刻。

4.3 收敛性与求解效率验证

为了验证二阶锥松弛的精确性，我对每个时段的松弛间隙$\eta_t = \max_{ij} \left[ v_{i,t} l_{ij,t} - (P_{ij,t}^2 + Q_{ij,t}^2) \right]$进行了统计。结果中所有时段的松弛间隙都小于$10^{-5}$，这说明二阶锥松弛的误差已经被压到很小，几乎等于0，模型是可靠的。

求解时间方面，EEH 33节点系统24时段动态重构，在CPLEX 12.10下，多线程8核，MIP gap设为1e-4时，求解时间为142秒。对于论文或课程设计来讲，这个速度完全可接受。我在更大的IEEE 123节点系统（含约200条可操作支路）上也试过，求解时间大约在25分钟级别，可以通过调低MIP gap到5e-3，误差可忽略，时间直接压到6分钟，非常适合灵敏度分析和多场景对比。

5. 工程实操中的常见问题与避坑指南

5.1 故障排查速查表

我把自己在调试过程中踩过的、以及帮学生调试时处理过的高频问题整理成一个表，希望能帮各位节省几个晚上的排查时间。

5.2 个人实测经验总结

最后分享几点在项目里反复验证过的实操体会，这些细节很难直接在论文里找到，但对跑通代码很重要。

关于数据单位，从MATLAB数据读取到YALMIP建模，全程保持标幺值制，只有最后输出结果时才转换为有名值。曾经图省事在部分环节用有名值，结果约束数值差距悬殊，CPLEX各种数值警告，排查了一整天才定位到单位不统一。

关于开关动作权重$\lambda$，这个系数不要照搬文献。文献值基于它们自己的成本和目的，未必适合你的系统规模和目标函数量纲。我的经验是：先设$\lambda=0$跑一遍，记录求解器给出的开关动作次数，再按这个次数反推一个合理的$\lambda$，比如希望动作次数控制在5次以内，就把$\lambda$从小到大扫几组，找到兼顾网损和动作次数的拐点值。这种调参方式比拍脑袋有依据得多。

另外输出结果时，强烈建议把每个时段的开关状态单独存成矩阵，方便画拓扑图。我最后用plot配合graph函数画出24个时段的拓扑变化图，视觉冲击力很强，审阅者一眼就能看出重构策略的时序特性，比纯表格直观多了。

关于后续扩展，这套框架也有很强的延展性。目前模型是确定性的，如果想要考虑DG出力的随机性，可以在环节加入场景生成和鲁棒优化，代码里的约束结构基本不用变，只需要增加场景维度和相应约束。也可以把动态重构和网络规划结合，构成两阶段优化问题。总之，二阶锥+YALMIP+CPLEX这套组合，在配电优化这个方向上就像一套积木，掌握了之后，各种变体问题都能快速落地。