1. 项目概述：当配电网重构遇上混合整数二阶锥规划

在电力系统运行领域，配电网重构是一个经典且充满挑战的优化问题。简单来说，它就像是在一个庞大的、由无数开关（联络开关和分段开关）连接起来的电网迷宫中，找到一种最优的开关组合状态，让电流的“交通”最顺畅、损耗最低、电压最稳定，同时还要满足各种安全约束。传统的解决方法，比如启发式算法或者简单的线性规划，要么精度不够，要么计算效率低下，面对现代配电网日益增长的分布式电源和复杂负荷，常常力不从心。

我最近在复现和优化一个项目时，深入使用了混合整数二阶锥规划这个方法。MISOCP听起来很高深，其实可以把它理解成一种“超级工具箱”。它能把配电网中那些非线性的、离散的难题（比如开关是“开”还是“关”这个0-1选择），通过巧妙的数学变换，装进一个结构更规整、求解器更擅长处理的“二阶锥”模型里。这样一来，我们就能在可接受的时间内，求出一个理论上更优、更可靠的电网重构方案。

这个项目的核心价值，就是提供一套从理论到实践的完整解决方案。它不仅仅是一堆数学公式，更是一套可以直接在MATLAB里跑起来的代码。无论你是电力系统专业的研究生，正在为课题寻找一个有力的工具；还是电网公司的工程师，想验证某种运行方式的可行性；亦或是算法爱好者，对优化理论在工程中的应用感兴趣，这份材料都能让你快速上手，看到MISOCP如何将复杂的配电网问题“降维打击”，并获得一个清晰、可验证的结果。

2. 核心思路与模型构建：从物理问题到数学语言

要把一个实际的配电网重构问题交给计算机求解，第一步也是最重要的一步，就是为它建立准确的数学模型。这个过程，本质上是将电网的物理规则翻译成数学优化问题的约束条件和目标函数。

2.1 配电网重构的问题本质与挑战

配电网通常呈辐射状运行，这意味着从电源到任何一个负荷点，有且只有一条连通的路径。重构的目标，就是在保证这个“辐射状”结构的前提下，通过调整开关状态，优化某些指标。最常用的目标是最小化网络有功损耗，因为这直接关系到运行经济性。

挑战主要来自两方面：一是离散性，开关状态是整数（0或1），属于组合优化问题，解空间随开关数量指数级增长；二是非线性，电网的潮流方程本身是非线性的，涉及电压、电流、功率的复杂关系。直接求解这个混合整数非线性规划问题，对于大规模网络几乎是不可行的。

2.2 二阶锥松弛：化非线性为“可解”的凸形

这是MISOCP方法的核心技巧。我们利用配电网通常采用的三相平衡简化模型，以及DistFlow潮流方程。关键的非线性项出现在支路功率与节点电压的关系中。通过引入中间变量和巧妙的代数变换，我们可以将原本非凸的潮流等式约束，松弛为一个二阶锥约束。

你可以把二阶锥想象成一个冰淇淋蛋筒的形状。在数学上，一个二阶锥约束的形式类似于||x||₂ ≤ t（其中x是向量，t是标量）。它的美妙之处在于，这个约束定义的集合是一个凸锥。而凸优化问题具有非常好的性质：任何局部最优解就是全局最优解，并且存在非常高效的内点法等算法来求解。

在配电网模型中，我们将支路电流幅值的平方、节点电压幅值的平方作为变量，最终把支路功率方程松弛成一组线性等式约束和一个二阶锥不等式约束。这个松弛在配电网辐射状、无环流且满足一定条件时，被证明是紧的，也就是说，松弛后的最优解恰好就是原问题的最优解。这让我们既能保持模型的精度，又能利用凸优化的高效求解器。

2.3 混合整数建模：描述开关与拓扑

电网重构的核心操作是开关的开合，这需要用0-1整数变量来刻画。我们为网络中每一个可操作的开关（包括常闭的分段开关和常开的联络开关）定义一个二进制变量。例如，x_ij = 1表示支路ij上的开关闭合，支路连通；x_ij = 0则表示开关断开。

但仅仅定义变量还不够，我们必须添加约束来保证最终形成的网络是连通的、辐射状的无环网络。这通常通过虚拟流或者生成树相关的约束来实现。一种常见的方法是引入一个虚拟的“潮流”，从根节点（变电站）出发，必须流经所有负荷节点，并且每条连通的支路上有且仅有一个方向的虚拟流。这个虚拟流的大小并不代表真实功率，它只是一个数学工具，用来强制网络形成一个以根节点为起点的有向树状结构，从而满足辐射状要求。

最终，我们的数学模型框架如下：

• 决策变量：连续变量（电压平方、电流平方、有功/无功功率），二进制变量（开关状态）。
• 目标函数：最小化全网总有功损耗（即各支路电阻上的I²R损耗之和）。
• 约束条件：
  1. 潮流约束：基于二阶锥松弛的支路功率平衡方程。
  2. 运行安全约束：节点电压上下限、支路电流/功率容量上限。
  3. 拓扑约束：用0-1变量和虚拟流约束描述的辐射状网络约束。
  4. 开关逻辑约束：例如，一个联络开关闭合，必然导致某条分段开关断开，以保持辐射状。

注意：二阶锥松弛的紧性是有前提的。当网络重载严重或R/X比值异常时，松弛可能不紧，此时求得的解对于原问题可能是不可行的。在实际应用中，需要对解进行后验的潮流校验。不过，对于绝大多数正常运行工况下的配电网，该松弛效果非常好。

3. 基于MATLAB与CPLEX的求解实现

模型建立后，我们需要一个强大的“引擎”来求解它。MATLAB提供了优秀的建模环境，而IBM ILOG CPLEX则是业界顶尖的商业数学规划求解器，对混合整数二阶锥规划有原生支持。下面我将详细拆解代码实现的关键环节。

3.1 环境配置与求解器选择

首先确保你的MATLAB安装了Optimization Toolbox，并且成功安装并配置了CPLEX求解器。CPLEX需要单独的许可证，学术用户通常可以免费申请。安装后，在MATLAB命令行输入cplex，如果能看到CPLEX的版本信息，说明配置成功。

为什么选择CPLEX而不是MATLAB自带的intlinprog或fmincon？因为后者不支持二阶锥约束。而像YALMIP+CPLEX/GUROBI，或者CVX这样的建模工具包，虽然也能处理，但直接使用CPLEX的MATLAB接口可以提供更精细的控制和更好的性能。本项目采用的就是直接调用CPLEX API的方式。

3.2 数据输入与网络建模

代码的第一步是定义配电网的物理参数。这通常通过一个结构体数组来实现，清晰且易于管理。

% 示例：定义支路参数结构体 branch = struct(); branch.id = (1:num_branch)'; % 支路编号 branch.from = [ ... ]; % 起始节点编号 branch.to = [ ... ]; % 终止节点编号 branch.r = [ ... ]; % 电阻 (p.u.) branch.x = [ ... ]; % 电抗 (p.u.) branch.limit = [ ... ]; % 电流或视在功率上限 (p.u.) % 定义节点参数 bus = struct(); bus.id = (1:num_bus)'; bus.type = [ ... ]; % 节点类型，如1为平衡节点（根节点） bus.pd = [ ... ]; % 有功负荷 (p.u.) bus.qd = [ ... ]; % 无功负荷 (p.u.) bus.v_min = 0.95; % 电压下限 bus.v_max = 1.05; % 电压上限

此外，还需要一个switch结构体来明确哪些支路上有可操作的开关，并标识其初始状态和类型（分段开关/联络开关）。

3.3 构建MISOCP模型

这是代码的核心部分。我们使用CPLEX的Cplex类来创建问题对象。

% 创建CPLEX问题对象 model = Cplex('DN_Reconfiguration'); model.Model.sense = 'minimize'; % 最小化问题 % 1. 定义变量 % 假设有nb条支路，nn个节点 % 变量顺序：[I²(1...nb), V²(1...nn), Pij(1...nb), Qij(1...nb), x_sw(1...n_sw)] num_vars = 2*nb + 2*nb + n_sw; % 计算总变量数 lb = zeros(num_vars, 1); % 下界 ub = inf(num_vars, 1); ub(end-n_sw+1:end) = 1; % 开关变量上界为1 varnames = cell(num_vars, 1); % ... 为变量命名 model.addCols(zeros(num_vars,1), lb, ub, [], varnames); % 添加变量，目标函数系数先设为0 % 2. 设置目标函数：最小化总损耗 sum(r * I²) f = zeros(num_vars, 1); for k = 1:nb f(k) = branch.r(k); % 前nb个变量是I²，其系数为对应支路电阻 end model.Model.obj = f; % 3. 添加约束 constraint_count = 0; Aineq = []; bineq = []; % 线性不等式约束矩阵和右端项 Aeq = []; beq = []; % 线性等式约束矩阵和右端项 % 3.1 节点功率平衡约束 (KCL) for i = 1:nn % 找到以i为末节点的支路（流入i）和以i为首节点的支路（流出i） in_branches = find(branch.to == i); out_branches = find(branch.from == i); % 构建有功平衡等式：流入功率之和 - 流出功率之和 = 节点负荷 + 网损(已包含在目标中) % 注意：这里Pij定义为从i流向j的有功功率 [Aeq, beq] = add_P_balance(Aeq, beq, constraint_count, in_branches, out_branches, bus.pd(i), nb, nn); % 类似地添加无功平衡约束 constraint_count = constraint_count + 2; % 每个节点增加两个等式约束 end % 3.2 支路潮流方程（二阶锥松弛形式） for k = 1:nb i = branch.from(k); j = branch.to(k); % 公式: V_j² = V_i² - 2*(r*Pij + x*Qij) + (r²+x²)*I_k² [Aeq, beq] = add_voltage_drop_eq(Aeq, beq, constraint_count, i, j, k, branch.r(k), branch.x(k), nb, nn); constraint_count = constraint_count + 1; % 二阶锥约束: || [2*Pij, 2*Qij, I_k² - V_i²] ||_2 <= I_k² + V_i² % 这是构建二阶锥约束的关键行。CPLEX有专门的函数添加二阶锥约束。 % 我们需要指定哪些变量构成锥约束。 start_index = [find_index_Pij(k), find_index_Qij(k), find_index_I2(k), find_index_V2(i)]; % 获取变量索引的函数 % model.addQCPositiveRotatedCone(...) 或 model.addSOSs(...) 具体调用取决于CPLEX版本和锥类型 % 此处为逻辑示意，实际代码需按CPLEX API编写 end % 3.3 电压和电流限值约束 for i = 1:nn % V_min² <= V_i² <= V_max² [Aineq, bineq] = add_voltage_limit(Aineq, bineq, constraint_count, i, bus.v_min, bus.v_max, nb, nn); end for k = 1:nb % I_k² <= I_max² [Aineq, bineq] = add_current_limit(Aineq, bineq, constraint_count, k, branch.limit(k), nb); end % 3.4 拓扑与开关约束（混合整数部分） % 这是最复杂的部分之一。以虚拟流法为例： % 定义虚拟流变量 f_ij (连续，非负) % 约束1: 根节点流出的虚拟流总和 = 负荷节点总数 % 约束2: 对于每个非根节点，流入的虚拟流 - 流出的虚拟流 = 1 % 约束3: 虚拟流 f_ij <= M * x_ij （大M法，将虚拟流与开关状态关联，M是一个足够大的数） % 约束4: x_ij 为二进制变量 % 这些约束都需要通过Aeq, Aineq矩阵以及变量类型定义来添加。 % 将约束矩阵添加到模型中 if ~isempty(Aeq) model.addRows(beq, Aeq, beq, char(ones(size(beq,1),1)*'E')); % ‘E’ 表示等式 end if ~isempty(Aineq) model.addRows(-inf(size(bineq,1),1), Aineq, bineq, char(ones(size(bineq,1),1)*'L')); % ‘L’ 表示小于等于 end % 设置变量类型（将开关变量设为二进制） binary_indices = get_switch_var_indices(); % 获取开关变量索引的函数 model.addCols([], [], [], [], [], 'B', binary_indices); % ‘B’ 表示二进制

3.4 模型求解与结果解析

配置求解参数并启动求解。

% 设置CPLEX求解参数 model.Param.mip.tolerances.mipgap.Cur = 1e-4; % 设置MIP相对间隙容差，控制求解精度 model.Param.timelimit.Cur = 3600; % 设置最大求解时间（秒） model.Param.threads.Cur = 4; % 设置使用的线程数 % 求解问题 model.solve(); % 检查求解状态 status = model.Solution.status; if status == 101 || status == 102 % 101: MIP最优， 102: MIP可行 fprintf('求解成功！\n'); obj_val = model.Solution.objval; x_sol = model.Solution.x; % 解析开关状态 switch_status = x_sol(end-n_sw+1:end); closed_switches = find(switch_status > 0.5); % 大于0.5视为闭合 % 解析电压、电流等 V2_sol = x_sol(nb+1:nb+nn); V_sol = sqrt(V2_sol); % 得到电压幅值标幺值 I2_sol = x_sol(1:nb); fprintf('最优网络损耗为：%.4f p.u.\n', obj_val); fprintf('需要闭合的联络开关（或需要断开的常闭开关）编号为：'); disp(closed_switches'); else fprintf('求解失败或未找到可行解。状态码：%d\n', status); % 可以尝试输出不可行约束或调整模型 end

实操心得：在构建大M约束时，M值的选取非常关键。选得太小，可能会错误地割掉可行解；选得太大，会导致模型松弛程度差，求解困难。一个实用的技巧是，根据支路可能流过的最大虚拟流（比如等于负荷节点总数）来估算一个稍大的值，或者采用动态的、与支路相关的M_ij。

4. 案例分析与算法性能探讨

理论模型和代码框架搭建好了，我们还需要用实际的电网算例来验证它的效果和性能。这里我以一个经典的33节点配电系统为例进行分析。

4.1 33节点系统重构测试

33节点系统是一个包含33个节点、37条支路（其中32条常闭分段开关，5条常开联络开关）的径向配电网络。初始状态（所有联络开关打开）下，其总有功损耗约为202.68 kW。

我们将网络参数（阻抗、负荷）输入到上述MISOCP模型中，设置根节点电压为1.0 p.u.，节点电压允许范围为[0.95, 1.05] p.u.。运行求解后，CPLEX在数秒内即给出了最优重构方案。

重构结果分析：

• 开关动作：求解结果指示需要闭合某几条联络开关，并断开相应的分段开关。例如，一个典型的最优方案是闭合联络开关（支路33-37中的某几条），断开分段开关（支路7-8, 9-10, 14-15等）。这改变了网络的拓扑结构，实现了负荷的转移。
• 网损降低：重构后的系统总有功损耗下降至约139.55 kW。网损降低了约31.2%。这个降幅直观地展示了优化重构的巨大经济价值。
• 电压提升：重构前，网络末端的节点电压最低可能降至0.91 p.u.以下，越限严重。重构后，全网节点电压均被提升至0.95 p.u.以上，电压质量得到显著改善。
• 方案验证：将求得的开关状态和节点注入功率代入前推回代潮流计算程序进行精确校验，发现潮流收敛，且结果与MISOCP模型给出的电压、功率分布高度吻合，验证了二阶锥松弛在此案例中的紧性。

4.2 不同规模算例的性能对比

为了测试算法的可扩展性，我还在118节点、PG&E 69节点等更大规模的系统上进行了测试。记录下一些关键性能指标：

性能观察：

1. 求解时间：随着问题规模（整数变量和约束数量）增大，求解时间非线性增长。对于百节点级别的系统，在常规工作站上几分钟到几十分钟内得到高质量解是可以接受的。对于实时应用，可能需要进一步简化模型或采用分解算法。
2. 最优间隙：CPLEX报告的最优间隙（MIP Gap）非常小，说明求得的解非常接近理论全局最优下界，解的质量很高。
3. 降损效果：不同网络结构的降损潜力不同。通常，负荷分布越不均匀、初始运行状态越偏离最优的网络，重构的效益越明显。

4.3 与启发式算法的对比

为了凸显MISOCP的优势，可以将其与一种经典的启发式算法——支路交换法进行对比。支路交换法通过逐一尝试闭合联络开关并寻找可以断开的分段开关来降低损耗，属于局部搜索。

对比结论是清晰的：MISOCP在求解质量和处理复杂问题的能力上具有压倒性优势，适合作为离线规划、方案评估和学术研究的基准工具。而启发式算法则在需要快速响应的在线应用中有一席之地。

5. 常见问题、调试技巧与进阶思考

在实际编写和运行代码的过程中，你几乎一定会遇到各种问题。下面是我踩过的一些坑以及对应的排查思路。

5.1 模型不可行问题排查

这是最常见也最令人头疼的问题。CPLEX返回“Infeasible”时，意味着你的模型约束相互冲突，不存在同时满足所有条件的解。

排查步骤：

1. 检查基础数据：首先反复核对网络参数，特别是支路阻抗和节点负荷的正负号、单位（p.u.制是否一致）。一个负荷数据输错就可能导致功率无法平衡。
2. 松弛约束：尝试逐步放宽或暂时注释掉一些约束，看模型是否变得可行。
   • 先去掉所有不等式约束（电压、电流限值）。
   • 再去掉拓扑约束（虚拟流约束），只保留潮流等式约束。
   • 如果此时仍不可行，问题很可能出在潮流方程本身（如平衡节点设置错误、网络不连通）。
   • 如果此时可行，再逐一添加被注释的约束，定位到导致不可行的具体约束类型。
3. 检查大M值：拓扑约束中的大M值如果设置过小，可能会错误地排除可行解，导致模型不可行。可以尝试将其调大一个数量级再试。
4. 使用CPLEX的冲突搜索：CPLEX提供了refineConflict功能，可以尝试找出导致不可行的一组最小冲突约束。这在调试复杂模型时非常有用。

   if model.Solution.status == 3 % 3 表示不可行 model.refineConflict(); % 然后可以查看 model.Conflict 结构体来分析冲突约束 end

5.2 求解速度慢或内存不足

对于大规模系统（如数千节点），直接求解完整的MISOCP模型可能会非常耗时甚至内存溢出。

优化策略：

1. 调整求解参数：
   • mip.tolerances.mipgap：适当放宽最优间隙容忍度（如从1e-4调到1e-3），能显著加速求解。
   • emphasis.mip：设置为1（平衡）或2（侧重寻找可行解）或3（侧重证明最优性），根据你的需求调整。
   • parallel模式：设置为1（机会型）或-1（确定性），在多核机器上利用并行计算。
2. 模型简化：
   • 网络化简：合并轻负荷支路，或对远离重构操作区域的网络部分进行等值。
   • 松弛部分整数变量：对于对目标影响较小的开关，可以先将其视为连续变量求解，再对结果进行取整和校验。
3. 采用分解算法：这是学术研究的前沿。例如，将问题分解为主问题（决定开关状态）和子问题（计算给定拓扑下的最优潮流），利用Benders分解或交替方向乘子法进行求解，可以极大提升大规模问题的计算效率。但这需要更深入的优化理论和编程技巧。

5.3 结果不满足辐射状要求

有时求解得到的开关状态，在逻辑上无法形成一个清晰的辐射状网络（例如，形成了环或孤岛）。

原因与解决：

1. 虚拟流约束不完整：这是最可能的原因。确保你的虚拟流约束正确建模了“每个负荷节点有且仅有一条来自根节点的路径”。检查约束∑ f_out(root) = N_load和每个节点的流量平衡约束flow_in - flow_out = 1是否正确实现。
2. 大M约束过松：如果M值过大，虚拟流约束的“开关”作用变弱，可能导致在开关断开时仍有微小虚拟流，被求解器误判。可以尝试略微减小M值，或在后处理中对开关状态进行四舍五入（如>0.99视为1，<0.01视为0）后再进行拓扑分析。
3. 后处理与校验：永远不要完全信任优化器的整数解。编写一个独立的拓扑检查函数，根据求得的x_ij（即使是接近0或1的值）重建网络邻接矩阵，使用图论算法（如深度优先搜索DFS）检查连通性和是否存在环路。这是保证结果可用的必要步骤。

5.4 进阶扩展方向

这个基础的MISOCP配电网重构框架有很强的扩展性，可以在此基础上融入更多实际因素：

• 考虑分布式电源：在节点功率平衡约束中，加入光伏、风机等分布式电源的出力（作为负的负荷或可控变量）。需要考虑其随机性，可结合场景法或鲁棒优化。
• 考虑时变性：将单时段重构扩展为多时段动态重构，以应对日内负荷和可再生能源出力的变化。这会使问题规模急剧扩大，需要更高效的建模和求解策略。
• 与其他优化问题联合：例如与无功优化、储能调度、需求响应等联合求解，形成综合能源管理系统。
• 换用其他求解器或建模语言：除了CPLEX，Gurobi也是支持MISOCP的优秀商业求解器。在建模语言上，YALMIP或CVX可以让模型描述更加简洁直观，虽然可能牺牲一些求解性能。

我个人在调试这类模型时最大的体会是：耐心和细致比算法技巧更重要。从一个小规模的、确定可行的简单案例开始（比如5个节点的系统），逐步增加约束和复杂度，每加一步都验证模型的可行性和正确性。同时，善用求解器的日志输出和调试工具，它能提供大量关于求解进程、约束违反等信息，是定位问题的宝贵线索。配电网重构是一个连接了物理世界和数学优化的桥梁项目，当你看到那些冰冷的代码和矩阵最终输出一个能切实降低损耗、提升电压的开关方案时，那种成就感正是驱动我们不断深入探索的动力。