主动配电网最优潮流：二阶锥规划建模与IEEE33节点求解实践

简介：面向配电网优化调度研究者的二阶锥规划（SOCP）最优潮流求解代码包，在动态最优潮流中充分计及风电、CB、SVG、OLTC等有源与无功设备，基于MATLAB+YALMIP+CPLEX构建并求解多时段SOCP模型，兼顾精度与效率。资源共23个文件，以M源码、log运算跟踪、PDF文献与CAJ论文为主，辅以PPT讲解、PNG拓扑图和DOC复现笔记，压缩包约117.18MB，内容涵盖IEEE33节点配电网结构、24小时动态优化实现及店主编写的复现全过程文档。目前已有843人学习浏览，适合具备一定最优潮流基础、希望深入研究主动配电网协调优化与二阶锥方法的读者。通过源码与资料互补，可快速掌握SOCP建模技巧、多设备协同策略及YALMIP+CPLEX调参排错思路。

1. 为什么主动配电网最优潮流用二阶锥规划而不是传统牛拉法

配电网优化调度里有个绕不开的分歧：潮流计算用牛拉法、前推回代都行，为什么基于二阶锥规划的主动配电网最优潮流偏偏要建SOCP模型？牛拉法只能求给定运行点的潮流解，而主动配电网最优潮流要的是一组设备决策——风电、CB、SVG、OLTC、储能怎么协同，才能让24小时运行成本最低、电压不越限。这是优化问题，不是单纯的计算问题。

二阶锥规划把非凸潮流方程松弛为凸锥约束，整套模型交给CPLEX一次求解，IEEE33节点24时段用例分钟内出全局最优解。这套代码覆盖从模型推导、参数设置到求解验证的全过程，附带讲解视频，适合做配电网优化调度、微网优化调度研究和工程复现的人直接上手。

2. SOCP松弛推导与YALMIP建模：从非凸潮流方程到锥约束

2.1 DistFlow方程的非凸性来源与二阶锥松弛

主动配电网的辐射状结构决定了分支潮流方程（DistFlow）是描述功率流动最自然的形式。以支路ij为例，受端净注入功率与送端流出的关系为Pj_net = Pij - rij * Iij²，其中支路电流幅值平方Iij² = (Pij² + Qij²) / Vi²，把电流写成支路功率与节点电压的函数。Q侧同理，Qj_net = Qij - xij * Iij²。

这里藏着两个非凸点：一是Iij²被电压平方除，形成分式非线性；二是Vi²乘以Iij²等于Pij² + Qij²，这是一个二次等式约束。单个节点看似乎只是几个方程，但放到32条支路、24个时段里，每一个等式约束都向外拓展出一个非凸曲面，外层再套整数变量搜索时，解空间形状会变得非常可怕。

SOCP松弛的做法是引入两个辅助变量去掉分式：vi = Vi²，lij = Iij²，再把二次等式Vi² * lij = Pij² + Qij²松弛成不等式Vi² * lij ≥ Pij² + Qij²。这个不等式经过换元恰好等价于旋转锥：

|| [2Pij; 2Qij; lij - vi] ||₂ ≤ lij + vi

展开后左边四项平方和小于等于右边平方，整理正是lij * vi ≥ Pij² + Qij²。这个锥约束是凸的，CPLEX这类求解器可以高效处理。

松弛的代价是什么？最优解可能越出原问题的等式曲面。但理论保证在辐射状配电网中，只要目标函数是网损的增函数（网损最小、购电成本最小都满足），松弛解会自动回到等式上，这就是SOCR精确性的含义。这个性质是整个二阶锥规划能用于主动配电网最优潮流求解的基石，也是第5章验证环节要检查的对象。

2.2 YALMIP的cone与旋转锥写法

YALMIP里声明二阶锥约束用cone函数，第一个参数是锥体向量，第二个参数是锥半径。对应旋转锥约束，写出来是这样：

% YALMIP下SOCP模型的变量声明与锥约束写法 % nb为节点数，nl为支路数 Pij = sdpvar(nb, nb, 'full'); % 支路有功，声明为全矩阵便于下标访问 Qij = sdpvar(nb, nb, 'full'); % 支路无功 vi = sdpvar(nb, 1); % 节点电压幅值平方 lij = sdpvar(nb, nb, 'full'); % 支路电流幅值平方 % 旋转锥约束：||[2Pij; 2Qij; lij - vi]|| <= lij + vi Constraints = [Constraints, cone([2*Pij(i,j); 2*Qij(i,j); ... lij(i,j) - vi(i)], lij(i,j) + vi(i))];

变量说明：Pij、Qij声明成全矩阵是为了下标运算方便，但代价是变量数量按节点对增长。IEEE33节点如果声明成33×33，会产生大量不存在的支路变量，求解器要用等于0的约束把它们固定住，这会把模型搞大好几倍。我一般会先用line矩阵读入支路首末端编号，把变量声明成nl×1，再通过映射矩阵组装节点功率平衡约束。压缩包里IEEE33_2.m和IEEE33BW.m两个版本，差异点就在这些变量怎么组织，clone日志显示店主在这一步反复对比过。

cone函数的参数不用手动展开旋转锥，YALMIP识别到第二个参数是线性表达式且第一个参数里出现配对项时，会自动生成对应旋转锥的底层表达式。实际调试时如果发现YALMIP报“The solver is not applicable”，多数情况是有二次等式混进了模型，趁早用sdisplay检查有没有漏松弛的项。

2.3 目标函数与松弛紧性的配合

SOCP松弛的精确性与目标函数形式强相关。常见三种目标：网损最小、购电成本最小、弃风量最小，都是随电流或注入功率单调变化的函数，SOCP松弛不会产生伪解。但若目标函数改为max(节点电压)或max(变压器负载率)，松弛后可能出现不可行却“最优”的解。因此建模时目标函数优先保持对电流和损耗单调增的性质，这是复现任何SOCP代码时最先检查的一行。

3. IEEE33节点24h算例：Wind/CB/SVG/OLTC/ESS的SOCP实现

3.1 IEEE33节点结构与基准参数

IEEE33是标准辐射状配网算例：基准电压12.66kV，基准功率10MVA，33个节点、32条支路，节点1经根变压器接上级网络。这套代码在IEEE33基础上接入风电、CB、SVG、OLTC和ESS，基本覆盖主动配电网最常见的可控资源。压缩包里的潮流计算.pptx和IEEE33节点配电网结构.png把网架和潮流分布讲得很清楚，建议读代码前先看这两份材料。

设备接入位置不同，控制量和变量类型也不同，我整理参数时习惯用一张表收着：

这个表格有一个工程含义：连续变量占比越高，求解越接近纯SOCP；整数变量占比上升，问题快速向MI-SOCP滑落。CB分5组还是10组、OLTC取5档还是9档，直接影响CPLEX分支定界的搜索树规模，后面第4章会看到这个影响在实际求解时长上的体现。

3.2 风电出力曲线与CB的分组整数建模

风电在日前优化调度里通常按预测出力曲线处理：24个时段各给一个有功注入值，作为参数写进节点功率平衡方程。乔珊的论文《主动配电网多源协同运行优化研究》对风电场景概率建模做了更细的讨论，工程复现先把确定性预测曲线跑通，再扩充随机场景更稳。

CB的分组投切是典型的离散决策。每组电容器容量固定，投切的是一组还是两组，对应整数的0/1组合。建模时用整数变量表示组数，再乘以单组容量：

% CB分组投切的整数建模 Q_CB_on = integer(sdpvar(N_cb_groups, 1)); % 每组投切状态，0或1 Constraints = [Constraints, Q_CB_on >= 0, Q_CB_on <= 1]; Q_CB = Q_CB_on' * Q_CB_step; % 总补偿无功，Q_CB_step为每组容量向量 Q_SVG = sdpvar(1); Constraints = [Constraints, Q_SVG_min <= Q_SVG <= Q_SVG_max]; % 节点无功平衡方程：负荷 - 风电 - CB - SVG Constraints = [Constraints, Q_inj(node) == Q_load(node) - Q_Wind(node) ... - Q_CB - Q_SVG];

参数说明：Q_CB_step是行向量，元素为每组电容器容量，单位Mvar，典型单组容量在0.05~0.3Mvar之间；Q_SVG_min/Q_SVG_max由SVG容量决定，常用对称区间。把CB和SVG放在同一个无功平衡方程里，恰好体现“离散粗调+连续细调”的配合逻辑——白天电压偏高时切CB、调SVG进相运行，夜间负荷轻时反之，这是主动配电网优化调度区别于普通无功优化的核心场景。

3.3 OLTC变比与ESS储能的时序约束

OLTC建模有两个关键点：变比离散化和动作次数限制。变比范围一般取0.95~1.05，档距0.025，共5个档位；用整数变量选档位后，二次侧电压等于一次侧电压乘变比。动作次数限制把不同时段耦合起来，需要引入二进制切换变量：

% OLTC档位选择与动作次数约束（示意） tap = integer(sdpvar(24, 1)); % 每个时段的档位 Constraints = [Constraints, 0 <= tap <= 4]; % 5个档位，编号0~4 delta = binvar(24, 1); % delta(t)=1表示变比发生变化 M = 10; % 大M常数，大于档位差上限 for t = 2:24 Constraints = [Constraints, tap(t) - tap(t-1) <= M * delta(t)]; Constraints = [Constraints, tap(t-1) - tap(t) <= M * delta(t)]; end Constraints = [Constraints, sum(delta) <= 4]; % 全天动作不超过4次

大M法的逻辑：档位差非零时delta必须为1，否则不等式无解；档位相同时delta可自由取0，配合sum(delta)<=4限制总动作次数。M只要大于最大档位差即可，取10是工程习惯，取得过大会引发数值病态，比如CPLEX在分支界探索时出现大范围数值溢出。

ESS是另一个多时段耦合点。储能SOC逐时递推，同时要避免充电、放电同时发生，用二进制变量互斥：

% ESS储能约束：SOC递推 + 充放电互斥 E_soc = sdpvar(25, 1); % 每时段末SOC，含t=0初值 P_ch = sdpvar(24, 1); % 充电功率 P_dis = sdpvar(24, 1); % 放电功率 u_ch = binvar(24, 1); % 1为充电时段 Constraints = [Constraints, E_soc(1) == E_init]; % 初始电量 for t = 1:24 Constraints = [Constraints, E_soc(t+1) == E_soc(t) ... + eta_ch * P_ch(t) - P_dis(t) / eta_dis]; % 效率不对称 Constraints = [Constraints, 0 <= P_ch(t) <= P_ch_max * u_ch(t)]; Constraints = [Constraints, 0 <= P_dis(t) <= P_dis_max * (1 - u_ch(t))]; Constraints = [Constraints, E_soc_min <= E_soc(t) <= E_soc_max]; end

eta_ch和eta_dis分别是充、放电效率，两个方向损耗不同所以必须分开写；功率上限乘二进制变量的写法把“充电时才能有充电功率”翻译成线性约束，这是标准做法。跑完模型后看储能出力曲线，可以和微电网优化调度里光伏+储能+分时电价的经典结果对照，验证动态过程是否符合预期。

注意：OLTC档位如果错声明成连续变量，求解器会给出“变比落在两档之间”的不可执行解。这是检查别人代码时最常发现的建模错误之一。

3.4 24h多时段模型的时间尺度与数据组织

24时段、每时段1小时是日前优化调度的标准粒度。数据组织上注意三点：负荷和风电曲线是24×1向量，节点负荷用基准值乘曲线系数；OLTC和ESS的跨时段约束要求变量维度从一开始就是24×1，不能只建单时段再循环；求解后用value()按时段取出结果时，注意SOC和OLTC变量的时段索引对齐。压缩包clone*.log记录了店主复现中每次求解的日志，能看到每个时段的约束数和变量数变化，排查维度问题时很有参考价值。

4. CPLEX求解配置与多时段动态调度的收敛性排查

4.1 YALMIP+CPLEX环境设置与参数说明

MATLAB里跑这套代码，前提是装好YALMIP和CPLEX，CPLEX 12系列与近几年发布的MATLAB版本兼容性都不错。求解选项在YALMIP里通过sdpsettings统一传入：

ops = sdpsettings('solver', 'cplex', 'verbose', 2, ... 'cplex.mip.tolerances.mipgap', 1e-4, ... 'cplex.mip.tolerances.absmipgap', 1e-3, ... 'cplex.output.clonelog', 1); sol = optimize(Constraints, Objective, ops);

参数说明：solver指定求解器为cplex；verbose控制求解器输出详细程度，2级能看到分支定界进度但不会刷屏，调试时够用；mipgap是相对最优间隙，混合整数问题几乎不可能收敛到严格0，设1e-4表示找到与下界相差0.01%以内的整数解就提前停止；absmipgap是绝对间隙，目标量级较大时防止相对间隙卡得太死；output.clonelog即CPLEX的clone日志开关，压缩包里大量clone*.log就是这里生成的，排查分支过程时很有用。

提示：verbose等级不要开到4以上，变量数千级的MI-SOCP模型日志本身会成为瓶颈，而且真正有用的分支定界信息在clone日志里都有。

4.2 变量规模与求解时间预估

24时段的MI-SOCP规模取决于设备数量和离散变量数。IEEE33全设备模型大约有几千个连续变量加一两百个整数/二进制变量，CPLEX默认参数下通常几十秒到几分钟收敛。判断模型轻重有个经验规律：整数变量每增加一档，分支定界节点数近似翻倍，CB分5组和分10组，求解时间不是2倍而是5到10倍，所以很多二阶锥规划代码里CB组数和OLTC档位都设置得比较粗。

给一个IEEE33典型配置的规模参考，具体数值随接入节点不同会漂移：

如果跑半小时不收敛，优先检查：整数变量是否声明正确；mipgap是否设得过小；OLTC大M常数是否过大导致数值病态。这三个原因占了八成以上。

4.3 模型不可行与求解器报错的定位顺序

MI-SOCP最典型的报错是infeasible和unbounded。前者是约束冲突，后者通常是某个变量缺边界。排查顺序我从调试经验整理如下：

1. 先用optimize(Constraints, [], ops)纯可行性求解，不带目标函数，区分约束冲突和目标函数问题。
2. 看diagnostics.problem字段，等于1是infeasible，等于2是unbounded，注意0才是求解成功。
3. 对可疑约束组逐个用check(Constraints)看残差，负残差超过1e-5的约束就是违约点。
4. 检查ESS末时段SOC是否要求回流到初值，储能模型最常见的死因是E_soc(25)=E_init与24时段充放电约束冲突。
5. 检查OLTC档位范围是否覆盖实际电压变化区间，档位给窄了会直接制造无解区域。

% 可行性诊断示例 ops0 = sdpsettings('solver', 'cplex', 'verbose', 0); diagnostics = optimize(Constraints, [], ops0); % 不带目标函数 if diagnostics.problem == 1 [res, ~] = check(Constraints); idx = find(res < -1e-4); % 找违约约束 fprintf('违约约束 %d 个，首个索引: %d\n', length(idx), idx(1)); end

check(Constraints)返回的向量里，正值表示满足且有余量，负值表示违约，绝对值越大违约越严重。这个技巧在处理多时段耦合模型时特别实用，能直接定位到具体时段、具体支路，而不是靠猜。

4.4 求解卡住时的常见对策

MI-SOCP比普通MIP慢，因为锥约束的分支定价过程每次下界计算都要做锥规划内点迭代。经验上的对策：CB和OLTC的档位参数预先算好离散点，不要写成求解器每节点现场计算的表达式；ESS的SOC上下限各留2%~5%裕度，避免数值抖动把求解器拖进病态边界；CPLEX parallel模式默认开启，多核机器上求解时间能压到单核三分之一左右。如果目标只是验证算法，还可以把mipgap放宽到1e-3，一天24时段的方案在10分钟内完成求解是完全可以实现的。

5. 用对偶间隙和潮流回代验证SOCP解的可执行性

5.1 验证松弛紧性的两条路径

SOCP模型求解完成后，第一件事不是看目标函数值，而是验证松弛是不是紧的。紧性验证有两条互补路径：约束层面的残差检查，直接看原始等式lij * vi与Pij² + Qij²的偏差；结果层面的潮流回代，把优化得到的注入功率代回原始潮流方程重新算一遍，对比节点电压。两者结合才能确定最优解能不能下发给调度执行。

5.2 约束残差与回代验证代码

% 松弛紧性检查：计算原始等式残差 V_opt = sqrt(value(vi)); % 由辅助变量还原电压幅值 L_opt = value(lij); P_opt = value(Pij); Q_opt = value(Qij); eq_residual = abs(L_opt .* V_opt.^2 - (P_opt.^2 + Q_opt.^2)); fprintf('最大等式残差: %.2e\n', max(eq_residual(:)));

max(eq_residual(:))在1e-4以下，说明SOCP松弛基本紧，最优解可直接作为调度方案。若某个时段残差到1e-2量级，说明该时段目标函数没有充分驱动松弛回到等式，常见处理是给目标函数加一个小正则项，比如0.001 * sum(lij)，用可忽略的目标偏移把电流变量压回物理值。

潮流回代更直接：把24个时段的节点注入功率逐个代回标准潮流计算程序，得到校验电压，与SOCP输出的电压比较，最大偏差小于0.001p.u.即认为解可执行。这一步对配电网优化调度特别重要——SOCP松弛解在凸意义上是数学最优，但调度员要的是真实物理网络能复现的控制指令。

5.3 判断结果可执行性的几条实用标准

最后给几条判别经验：电压幅值落在0.95~1.05p.u.内且没有恰好压在边界上；OLTC全天动作次数严格小于等于约束值，说明动作限制没被逼到饱和；ESS的SOC轨迹平滑，不出现锯齿形充放；CB投切次数合理，与SVG出力没有互相抵消耗散。四条全满足，SOCP求出的主动配电网优化调度方案基本可以直接用于日前调度计划编制。

本文还有配套的精品资源，点击获取