企业官网建设流程全解析

简介：基于二阶锥规划的主动配电网最优潮流求解程序包，面向电力系统专业研究生、配电网规划与运行研究人员，以及具备一定MATLAB/CPLEX基础的学习者。资源以IEEE33节点配电网为算例，实现含风电（Wind）、并联电容器（CB）、静止无功发生器（SVG）、有载调压变压器（OLTC）及储能系统（ESS）的多时段24h协同优化，可帮助读者掌握主动配电网最优潮流的建模方法与二阶锥松弛求解技巧。程序包共17个文件，压缩包5.51MB。其中2个m文件为MATLAB主体程序，提供骨灰级注释便于逐行理解；12个log文件为CPLEX求解过程记录，便于对照分析迭代与收敛情况；png为配电网结构示意图，pptx为潮流计算原理说明，另有参考文献zip可供延伸学习。这套程序包目前已有835人学习下载，注释详尽、案例完整，尤其适合初学者循代码逐步搭建主动配电网优化框架，并快速迁移至自己的研究场景。

1. SOCP 这块硬骨头，CPLEX 怎么啃下来

分布式光伏、风电大量接入 10kV 馈线后，潮流方程从线性代数问题变成非凸优化问题，传统的牛拉法只能做潮流计算，没法直接做“运行点寻优”。二阶锥规划（SOCP）是当前主动配电网最优潮流里工程化程度最高的松弛手段。这套基于 MATLAB 的代码，把 24h 多时段、风电机组（Wind）、电容器组（CB）、静止无功发生器（SVG）、有载调压变压器（OLTC）和储能（ESS）全部拉进一个 CPLEX 可解的锥优化模型，在 IEEE33 节点系统上做完整的最优潮流分析，注释细到每个变量和每行约束都有说明。比较适合做配电网课题、正在入门凸优化应用、或想把非凸最优潮流模型换成可求解 SOCP 的工程师。

2. 从 DistFlow 到二阶锥松弛：非凸潮流怎么变成 CPLEX 能啃的模型

2.1 辐射网下的 DistFlow 方程

配电网最优潮流不像输电网那样适合用节点导纳矩阵直接写，因为 10kV 馈线大多是辐射状结构，支路功率方向明确，用支路潮流方程写起来更直观，这就是 DistFlow 模型。记支路ij首端有功功率为P_ij、无功功率为Q_ij，母线电压幅值平方为V_i，支路电流幅值平方为I_ij，那么 in 每个时段t，节点j的功率平衡可以写成：

P_ij(t) - sum(P_jk(t)) - R_ij * I_ij(t) = P_load_j(t) - P_gen_j(t) Q_ij(t) - sum(Q_jk(t)) - X_ij * I_ij(t) = Q_load_j(t) - Q_gen_j(t)

其中sum(P_jk(t))表示以j为首端的所有下游支路功率之和。电压递推关系则是：

V_j(t) = V_i(t) - 2 * (R_ij * P_ij(t) + X_ij * Q_ij(t)) + (R_ij^2 + X_ij^2) * I_ij(t)

这三个等式只刻画了潮流守恒，真正让模型变难的是最后一条：I_ij(t) * V_i(t) = P_ij(t)^2 + Q_ij(t)^2。这个约束里既有电压平方、电流平方，又有支路功率的二次项，整体是非凸的，CPLEX 这类求解器没法直接吃进去。项目里面对这个等式做了二阶锥松弛，把“等号”改成“不小于”，再变换成标准锥形式。

2.2 为什么要松弛成锥而不是线性化

把I_ij >= (P_ij^2 + Q_ij^2) / V_i展开时，可以写成矩阵范数形式：

|| [2*P_ij; 2*Q_ij; I_ij - V_i] ||_2 <= I_ij + V_i

这个锥约束在数学上是凸的。对辐射状配电网，只要电压幅值有合理上下界、目标函数是网损最小，松弛后的最优解往往会落在原非凸问题的可行域边界上，也就是说松弛是精确的。实际调试时你会发现，CPLEX 日志里每个支路电流约束基本都是紧的，几乎没有“锥内点”的浪费。

这里还有个常见误用：有人为了省事，直接把I_ij设成常数，或者把P_ij^2+Q_ij^2做一阶泰勒展开。前者适合配电网规划估算，但做 24h 运行优化时会低估网损；后者在运行点附近小扰动下勉强可用，遇到 OLTC 抽头切换或 ESS 大功率充放就容易失真。这也是为什么在 IEEE33 节点上做多时段潮流优化，SOCP 比线性潮流更稳妥。

2.3 在 YALMIP 中写下第一个锥约束

项目里大量使用 YALMIP 建模，锥约束不要手写norm，而是直接用cone函数，CPLEX 识别更高效，模型也更干净。以支路ij为例：

% 支路电流平方 I_ij，电压平方 V_i，支路功率 P_ij、Q_ij for ij = 1:n_branch for t = 1:24 idx_i = br_from(ij); % 首端节点编号 idx_j = br_to(ij); % 末端节点编号 V_i = V_sq(idx_i, t); % 首端电压幅值平方 L = I_sq(ij, t); % 支路电流幅值平方 P = P_branch(ij, t); Q = Q_branch(ij, t); % || [2P; 2Q; L - V] ||_2 <= L + V Constraints = [Constraints, cone([2*P; 2*Q; L - V_i], L + V_i)]; end end

cone的第一个入参是方向向量，第二个入参是范数的上界标量，整体表达的是||向量||_2 <= 标量。注意向量里第二项是L - V_i，不是L + V_i，这在抄模型时非常容易写反。如果写反，求解结果会变得很奇怪，电压曲线在某几个节点上突然偏低，但约束检查又提示 infeasible。参数br_from和br_to可以直接从 IEEE33 初始数据生成，也可以用[1:32;2:33]这种简单矩阵手工构造，本质上是告诉 YALMIP 哪些节点之间允许有潮流。

3. 多时段元件建模：Wind/CB/SVG/OLTC/ESS 的锥约束怎么进模型

3.1 Wind 的预测功率边界与弃风惩罚

风电机组在最优潮流里通常不是简单设成恒定有功注入，而是给一个预测出力上限，实际出力由优化器决定。这么做是为了兼顾弃风：当线路电压越过上限或储能 SOC 接近满时，调度可以主动压低风力出力。代码里一般写成：

% n_wind 台风电机组，T 为 24 小时 P_wind_max = forecast_wind; % 24h 预测曲线，维度 n_wind * T P_wind = sdpvar(n_wind, T); % 实际调度出力，连续变量 % 出力下限 0，上限不超过预测 Constraints = [Constraints, 0 <= P_wind <= P_wind_max];

变量定义要放在一个Constraints序列里不断追加。forecast_wind在代码里可能是从 Excel 或.mat文件读取，也可能在 m 文件里直接写成 24 个数的数组。注意风电场接入节点通常是无功支撑较弱的末端，如果只约束有功，不约束无功，模型会从线路末端吸取大量无功，导致 CPLEX 求解时间上升，最好给风机加Q_wind的容量约束，例如-0.2 * P_wind <= Q_wind <= 0.2 * P_wind。

3.2 CB 离散投切与 SVG 连续无功

电容器组和 SVG 都是无功补偿设备，区别在 CB 是离散投切，SVG 是连续调节。CB 建模时用binvar表示每组投切状态，再乘以单组无功容量，得到总无功注入：

n_cb_group = 5; % 5 组电容器 q_cb_single = 0.1; % 每组 0.1 Mvar，标幺化后处理 u_cb = binvar(n_cb_group, 24); % 投切状态，1 表示投入 Q_cb = q_cb_single * sum(u_cb, 1); % 24 个时段的 CB 总无功 % 防止频繁投切：一天最大动作次数限制 for g = 1:n_cb_group Constraints = [Constraints, sum(abs(diff(u_cb(g,:)))) <= 4]; end

这里的sum(abs(diff(...)))是在统计相邻时段投切状态变化次数。diff对二进制变量做差分，结果可能是-1、0、1，取绝对值再求和就是动作次数。如果不加这个约束，CPLEX 求解出来的 CB 策略看着很漂亮，但实际没法执行，因为每半小时切一次电容器会严重缩短开关寿命。SVG 相对简单，直接用连续变量并限制上下限：

Q_svg = sdpvar(n_svg, 24); Constraints = [Constraints, -Q_svg_max <= Q_svg <= Q_svg_max];

两者在目标网损中的权重完全不同：CB 基本是离散投切，纳入目标没有成本项，容易和 OLTC 一起形成“整点切一刀”的锯齿策略；SVG 可连续调节，通常也会给一个小权重避免高频抖动。

3.3 OLTC 抽头如何避免变比平方的非线性

有载调压变压器通过改变变比k_t来调节电压。直接写V_secondary = k_t^2 * V_primary会产生k_t^2和非线性乘积，破化 SOCP 结构。常见做法是把每个离散档位对应的k^2拆成固定数值，然后用二进制整数选择：

tap_pos = intvar(n_tap, 24); % 整数变量，档位位置 % 假设 9 档，tap_pos 范围 -4 到 4 Constraints = [Constraints, -4 <= tap_pos <= 4]; % 通过重复矩阵或者表查映射，把 tap_pos 换成变比平方 k2_table = [0.975^2 0.98^2 0.985^2 0.99^2 1^2 1.01^2 1.015^2 1.02^2 1.025^2]; k2 = sdpvar(n_tap, 24); % 辅助连续变量 for t = 1:24 for r = 1:n_tap % 用 implies 方式填值，实际工程中常用查找表 + 大 M end end

严格说这里需要把整数变量与连续变量耦合，代码里通常会用一个表格矩阵做索引，或者用一列二进制变量对每个档位独热编码。CPLEX 支持 MISOCP，因此 OLTC 离散档位不会破坏整体模型结构，但会显著增加分支定界节点数量。调参时优先改善的就是这个部件：如果 24h 模型求解太慢，先放宽抽头动作次数约束，比调求解器参数更有效。

3.4 ESS 的 SOC 递推是 24h 模型的时间耦合核心

ESS 是唯一让不同小时之间产生耦合的元件。其余 Wind、CB、SVG 在时间维度上都是独立断面，只是参数滚动更新；ESS 的荷电状态SOC(t+1)依赖SOC(t)，这一条约束让整个模型变成真正的多时段问题。代码里通常用以下方式建模：

% dim: n_ess * (T+1)，多出一列存初始 SOC E_soc = sdpvar(n_ess, 25); P_ch = sdpvar(n_ess, 24); % 充电功率，>=0 P_dis = sdpvar(n_ess, 24); % 放电功率，>=0 dt = 1; % 时段间隔 1h，也可以写成 1 的标幺值 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]; Constraints = [Constraints, 0 <= P_dis(:,t) <= P_dis_max]; Constraints = [Constraints, E_min <= E_soc(:,t) <= E_max]; end Constraints = [Constraints, E_soc(:,1) == E_soc(:,25)]; % 24h 循环调度

eta_ch与eta_dis通常取 0.9 到 0.95，不要用同一个混着算，否则充电到放电之间会产生虚拟能量增益，CPLEX 会利用这个漏洞“凭空发电”，得到很小的网损，但物理上不可能。E_soc(:,1) == E_soc(:,25)是日循环边界条件，如果做跨天调度，可以改成和前一天终端 SOC 绑定。功率上限建议用额定功率，但有些代码里会默认 0.5MW，标幺化后要仔细换算。

4. MATLAB+YALMIP+CPLEX 链路：从 IEEE33 数据到 24h 求解

4.1 IEEE33 数据组织和标幺值选择

项目里的IEEE33BW.m和IEEE33_2.m就是入口脚本。打开后最前面通常是一堆基础数据：33 个节点、32 条支路的电阻电抗、每个节点的有功无功负荷，以及各台设备的接入母线编号。这些数据必须转成标幺值，否则 SOCP 迭代时数值范围为差 6 个量级，CPLEX 还没开始分支定界就已经先报数值警告。

baseMVA = 10; % 功率基准 10MW，视具体系统调整 baseKV = 12.66; % 电压基准 12.66kV，IEEE33 基准值 % 负荷标幺化 load_p_pu = load_p_mw / baseMVA; load_q_pu = load_q_mvar / baseMVA;

baseMVA的选取和网络电压等级相关，IEEE33 基准负荷约 3.7MW，用baseMVA=10会让大部分变量落在 0.01 到 1 之间，数值条件较好。如果直接用 MW 做单位，支路功率和网损差 3 个数量级，CPLEX 求解器内部的尺度化步骤会花掉大量时间。项目日志里的clone0.log到clone11.log就是不同参数组合下的 CPLEX 运行日志，通过对比这些日志可以清楚看到标幺值对迭代次数和求解时间的影响。

4.2 目标函数和 CPLEX 参数设置

目标函数常见组合是网损最小、弃风惩罚、抽头动作惩罚三部分加权：

% 网损：所有支路电阻 * 电流平方 * dt loss = sum(sum(R_branch .* I_sq)) * dt; % 弃风惩罚：预测值减实际出力 curtailment = sum(P_wind_max - P_wind); Objective = loss + penalty_curtail * curtailment;

实际代码里可能只有前两项，也可能额外加上 CB 动作次数惩罚。penalty_curtail不宜设得太大，否则会把网损优化的空间完全压掉，导致储能不充不放一直维持 SOC 边界。比较合理的做法是先单独跑一次不弃风情况下的网损下限，再把惩罚系数按网损的 3 到 5 倍设置。

求解器参数在sdpsettings里传入：

ops = sdpsettings('solver', 'cplex', ... 'verbose', 2, ... 'cplex.mip.tolerances.mipgap', 1e-4, ... 'cplex.timelimit', 3600, ... 'savesolveroutput', 1); sol = optimize(Constraints, Objective, ops);

cplex.mip.tolerances.mipgap控制最优性间隙，10kV 配电网下网损通常在 0.05 到 0.2MW 之间，设成1e-4已经足够；如果继续设成1e-6，CPLEX 会陷入长时间分支定界，而收益只是网损多准确几个小数点。timelimit是硬性保护，多时段含 OLTC、CB 的 MISOCP 模型很容易超过 30 分钟。savesolveroutput为 1 时会保留 CPLEX 原始日志到 YALMIP 结果对象里，排错时可以直接看sol.solveroutput.info。

4.3 求解完成后的结果对象怎么读取

optimize返回后，不要立刻value()所有变量，先看sol.info和sol.solvertime：

if sol.problem == 0 fprintf('求解成功，用时 %.2f s\n', sol.solvertime); else disp(sol.info); end

sol.problem == 0表示求解正常完成，1表示 infeasible，2表示 unbounded，9表示 NaN 值。如果sol.problem是 0 但后面value(V_sq)里出现 NaN，通常不是求解器问题，而是 YALMIP 变量没有全部进入约束集合，某个孤立变量没有被任何表达式引用，value之后是空。这时候去检查代码里是否有某个sdpvar变量定义了却没进入Constraints。

5. 求解失败时，CPLEX 日志和 check 函数怎么定位问题

5.1 先看sol.info，再看求解状态

遇到Infeasible时，很多人第一反应是扩大 CPLEX 容差，这是错误方向。infeasible只说明约束集合本身没有交点，和数值容差关系不大。项目日志里如果出现infeasible，优先检查 OLTC 的抽头变量和 CB 的投切变量是否作用到母线电压上。比如 CB 只是定义出了Q_cb，但没有写进节点无功平衡方程，那这个变量孤立存在不会导致 infeasible，反而是写进平衡方程但上下限冲突更容易触发。

在实际项目中的排查顺序是：

1. 先注释掉 ESS 的 SOC 递推约束，改成单时段断面分别求解，看每个时段是否可行。如果单时段可行、24h 不可行，说明 SOC 初值或容量边界设置不合理。
2. 再把 OLTC 抽头范围从-4:4扩大到-8:8，看是否依然 infeasible。若可行，说明电压约束和变比范围冲突。
3. 最后检查潮流方程里的I_ij松弛是否写漏了节点编号，特别是br_from和br_to从 1 开始索引还是从 0 开始。MATLAB 的索引从 1 开始，如果原始数据来自 Python 习惯，很容易错位。

5.2 用check(Constraints)检查约束紧度

YALMIP 的check函数会返回每条约束的残差，这是验证松弛平坦度和模型准确性最快的方法。代码片段如下：

% 求解结束后检查所有约束 residual = check(Constraints); [max_res, idx] = max(abs(residual)); if max_res > 1e-5 fprintf('最大误差出现在第 %d 条约束\n', idx); Constraints(idx) end

注意 YALMIP 的check返回值在约束满足时为负数或 0，这个负数值是松弛不等式左侧减右侧的差。通常max_res在1e-7到1e-5之间。如果某个 SOC 约束的残差达到1e-2，问题大概率不是求解精度，而是cone的参数顺序写反。比如把cone([P;Q], L)写成cone(L, [P;Q])，YALMIP 会把它解释成L <= norm([P;Q])，这个约束完全变了形。建议在代码里加一条调试命令：

% 输出锥约束的数量，对比与 n_branch*24 是否一致 fprintf('锥约束数量: %d\n', n_branch * 24);

5.3 CPLEX 参数里的几个“救命开关”

如果模型可行，但求解时间极长，或者日志里大量出现Numerical difficulties，需要调整 CPLEX 的线性求解器行为。常见做法是：

ops = sdpsettings(ops, 'cplex.lpmethod', 4); % barrier 方法 ops = sdpsettings(ops, 'cplex.mip.strategy.search', 1); % 深度优先 ops = sdpsettings(ops, 'cplex.mip.stagnation.nodes', 50000);

cplex.lpmethod设为 4 是 use barrier，对 SOCP 的根节点求解通常比默认单纯形更好。mip.strategy.search设为 1 是深度优先，适合中等规模但整数变量集中的模型。stagnation.nodes用来检测目标值长时间不变的停滞，达到阈值后 CPLEX 会自动加全局切开平面。

还有一类问题表现为求解结果可用，但某两个时段间 CB 投切完全相反，这通常是目标函数缺少对动作次数的惩罚，导致所有无功方案网损相同，CPLEX 随机选了一个最优顶点。这时候调整目标权重比调求解参数重要。可以给动作惩罚设一个很小的正数，只要大于数值噪声，就能稳定输出平滑策略。

6. 24h 结果的工程化校验：从 v²、I² 再到物理可信度

6.1 从优化变量回到传统潮流结果

SOCP 求解得到的是电压平方、电流平方和支路功率，直接画电压曲线时很多人会犯一个错误：用sqrt(value(V_sq))得到幅值，却没乘以母线基准电压。IEEE33 的基准电压是 12.66kV，标幺化后电压平方落在 0.950. 到 1.05、之间，换算成实际电压要乘回baseKV。代码里可以这样写：

V_pu = sqrt(value(V_sq)); % 标幺电压，维度 n_bus*24 V_kv = V_pu .* 12.66; % 实际电压 kV bus_10_charge = V_pu(10,:); % 挑一个末端节点观察 figure stairs(1:24, bus_10_charge)

观察第 10 号节点在夜间和午间的电压波动。风电机组高发时段如果电压升高超过 1.03pu，说明无功支撑或变压器档位调整不足以完全抑制电压抬升。这时候 CB 和 SVG 的出力会体现为无功注入，而 OLTC 的档位变化可以在value(tap_pos)上看到明显的整点跳变。要注意，SOCP 模型对电压平方做松弛，输出的V_sq理论上是满足网损最小的最优值，但不一定满足该节点单相电压的实际谐波和三相平衡约束，做工程推广时还需要回到三相潮流进一步校验。

6.2 一个马上可以验证的扩展：储能 SOC 下界灵敏度

这套代码最值得动手改的地方是 ESS 的E_min。把 SOC 下界从 0.1 改成 0.3，重新求解，然后对比网损和弃风量：

E_min_list = [0.1 0.2 0.3 0.4]; for k = 1:length(E_min_list) E_min = E_min_list(k); % 重复上一轮建模和求解 optimize(Constraints, Objective, ops); loss_result(k) = value(loss); wind_curtail_result(k) = value(curtailment); end

这个灵敏度分析可以很容易判断储能容量是否还有扩容空间。如果 SOC 下界从 0.1 提到 0.3 后网损只增加 0.5%，说明储能多数时段处于高 SOC 区域，容量冗余较充足；如果网损增加超过 5%，则说明储能深度充放对削峰填谷影响很大，此时加强对充放电功率的时序约束比单纯增加容量更划算。另有一个小技巧：把P_ch和P_dis的上下限从固定值改成P_ch_max * u_ess和P_dis_max * (1-u_ess)，引入一个充电放电互斥的二进制变量u_ess，可以避免求解结果出现同一储能同时充电和放电的“套利”假象。虽然目标函数本身会抑制这种浪费，但在无网损惩罚或双峰电价场景下，不加互斥约束时 CPLEX 经常给出同时充放的病态解。

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