简介基于二阶锥规划的主动配电网最优潮流求解程序包面向电力系统专业研究生、配电网规划与运行研究人员以及具备一定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^2Q_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 endcone的第一个入参是方向向量第二个入参是范数的上界标量整体表达的是||向量||_2 标量。注意向量里第二项是L - V_i不是L V_i这在抄模型时非常容易写反。如果写反求解结果会变得很奇怪电压曲线在某几个节点上突然偏低但约束检查又提示 infeasible。参数br_from和br_to可以直接从 IEEE33 初始数据生成也可以用[1:32;2:33]这种简单矩阵手工构造本质上是告诉 YALMIP 哪些节点之间允许有潮流。对比项DistFlow 原始形式SOCP 松弛形式适用网络辐射状配电网辐射状配电网核心变量支路功率、电压、电流平方同样变量多一个锥约束约束性质非线性等式非凸凸锥约束求解器支持需要 IPOPT 等非线性求解器CPLEX、MOSEK、Gurobi 原生支持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(t1)依赖SOC(t)这一条约束让整个模型变成真正的多时段问题。代码里通常用以下方式建模% dim: n_ess * (T1)多出一列存初始 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(:,t1) 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标幺化后要仔细换算。元件决策变量类型典型约束多时段耦合Wind连续有功/无功0 P 预测上限无CB二进制投切动作次数限制较弱SVG连续无功-Qmax Q Qmax无OLTC整数抽头档位范围、动作次数较弱ESS连续充放电功率SOC 更新、功率上下限强4. MATLABYALMIPCPLEX 链路从 IEEE33 数据到 24h 求解4.1 IEEE33 数据组织和标幺值选择项目里的IEEE33BW.m和IEEE33_2.m就是入口脚本。打开后最前面通常是一堆基础数据33 个节点、32 条支路的电阻电抗、每个节点的有功无功负荷以及各台设备的接入母线编号。这些数据必须转成标幺值否则 SOCP 迭代时数值范围为差 6 个量级CPLEX 还没开始分支定界就已经先报数值警告。baseMVA 10; % 功率基准 10MW视具体系统调整 baseKV 12.66; % 电压基准 12.66kVIEEE33 基准值 % 负荷标幺化 load_p_pu load_p_mw / baseMVA; load_q_pu load_q_mvar / baseMVA;baseMVA的选取和网络电压等级相关IEEE33 基准负荷约 3.7MW用baseMVA10会让大部分变量落在 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-6CPLEX 会陷入长时间分支定界而收益只是网损多准确几个小数点。timelimit是硬性保护多时段含 OLTC、CB 的 MISOCP 模型很容易超过 30 分钟。savesolveroutput为 1 时会保留 CPLEX 原始日志到 YALMIP 结果对象里排错时可以直接看sol.solveroutput.info。4.3 求解完成后的结果对象怎么读取optimize返回后不要立刻value()所有变量先看sol.info和sol.solvertimeif sol.problem 0 fprintf(求解成功用时 %.2f s\n, sol.solvertime); else disp(sol.info); endsol.problem 0表示求解正常完成1表示 infeasible2表示 unbounded9表示 NaN 值。如果sol.problem是 0 但后面value(V_sq)里出现 NaN通常不是求解器问题而是 YALMIP 变量没有全部进入约束集合某个孤立变量没有被任何表达式引用value之后是空。这时候去检查代码里是否有某个sdpvar变量定义了却没进入Constraints。CPLEX 参数作用项目中的建议值cplex.mip.tolerances.mipgap分支定界最优化间隙1e-4cplex.timelimit求解时间上限3600 秒cplex.threads并行线程数4 或 8cplex.mip.display求解日志显示频率2cplex.mip.limits.nodes最大节点数1e6 或留空5. 求解失败时CPLEX 日志和 check 函数怎么定位问题5.1 先看sol.info再看求解状态遇到Infeasible时很多人第一反应是扩大 CPLEX 容差这是错误方向。infeasible只说明约束集合本身没有交点和数值容差关系不大。项目日志里如果出现infeasible优先检查 OLTC 的抽头变量和 CB 的投切变量是否作用到母线电压上。比如 CB 只是定义出了Q_cb但没有写进节点无功平衡方程那这个变量孤立存在不会导致 infeasible反而是写进平衡方程但上下限冲突更容易触发。在实际项目中的排查顺序是先注释掉 ESS 的 SOC 递推约束改成单时段断面分别求解看每个时段是否可行。如果单时段可行、24h 不可行说明 SOC 初值或容量边界设置不合理。再把 OLTC 抽头范围从-4:4扩大到-8:8看是否依然 infeasible。若可行说明电压约束和变比范围冲突。最后检查潮流方程里的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 随机选了一个最优顶点。这时候调整目标权重比调求解参数重要。可以给动作惩罚设一个很小的正数只要大于数值噪声就能稳定输出平滑策略。日志特征真实含义项目里优先动作Infeasible约束无交集检查 ESS SOC 和 OLTC 档位边界Unbounded目标无下界检查风机无功和网损表达式方向Numerical difficulties数值刚度大调整lpmethod和标幺值Node limit exceeded分支定界节点爆炸放宽 mipgap 或抽头档位Integer optimal求到最优整数解读value前先看剩余间隙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 经常给出同时充放的病态解。本文还有配套的精品资源点击获取 配电网多时段SOCP最优潮流：基于MATLAB+YALMIP+CPLEX建模
简介基于二阶锥规划的主动配电网最优潮流求解程序包面向电力系统专业研究生、配电网规划与运行研究人员以及具备一定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^2Q_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 endcone的第一个入参是方向向量第二个入参是范数的上界标量整体表达的是||向量||_2 标量。注意向量里第二项是L - V_i不是L V_i这在抄模型时非常容易写反。如果写反求解结果会变得很奇怪电压曲线在某几个节点上突然偏低但约束检查又提示 infeasible。参数br_from和br_to可以直接从 IEEE33 初始数据生成也可以用[1:32;2:33]这种简单矩阵手工构造本质上是告诉 YALMIP 哪些节点之间允许有潮流。对比项DistFlow 原始形式SOCP 松弛形式适用网络辐射状配电网辐射状配电网核心变量支路功率、电压、电流平方同样变量多一个锥约束约束性质非线性等式非凸凸锥约束求解器支持需要 IPOPT 等非线性求解器CPLEX、MOSEK、Gurobi 原生支持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(t1)依赖SOC(t)这一条约束让整个模型变成真正的多时段问题。代码里通常用以下方式建模% dim: n_ess * (T1)多出一列存初始 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(:,t1) 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标幺化后要仔细换算。元件决策变量类型典型约束多时段耦合Wind连续有功/无功0 P 预测上限无CB二进制投切动作次数限制较弱SVG连续无功-Qmax Q Qmax无OLTC整数抽头档位范围、动作次数较弱ESS连续充放电功率SOC 更新、功率上下限强4. MATLABYALMIPCPLEX 链路从 IEEE33 数据到 24h 求解4.1 IEEE33 数据组织和标幺值选择项目里的IEEE33BW.m和IEEE33_2.m就是入口脚本。打开后最前面通常是一堆基础数据33 个节点、32 条支路的电阻电抗、每个节点的有功无功负荷以及各台设备的接入母线编号。这些数据必须转成标幺值否则 SOCP 迭代时数值范围为差 6 个量级CPLEX 还没开始分支定界就已经先报数值警告。baseMVA 10; % 功率基准 10MW视具体系统调整 baseKV 12.66; % 电压基准 12.66kVIEEE33 基准值 % 负荷标幺化 load_p_pu load_p_mw / baseMVA; load_q_pu load_q_mvar / baseMVA;baseMVA的选取和网络电压等级相关IEEE33 基准负荷约 3.7MW用baseMVA10会让大部分变量落在 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-6CPLEX 会陷入长时间分支定界而收益只是网损多准确几个小数点。timelimit是硬性保护多时段含 OLTC、CB 的 MISOCP 模型很容易超过 30 分钟。savesolveroutput为 1 时会保留 CPLEX 原始日志到 YALMIP 结果对象里排错时可以直接看sol.solveroutput.info。4.3 求解完成后的结果对象怎么读取optimize返回后不要立刻value()所有变量先看sol.info和sol.solvertimeif sol.problem 0 fprintf(求解成功用时 %.2f s\n, sol.solvertime); else disp(sol.info); endsol.problem 0表示求解正常完成1表示 infeasible2表示 unbounded9表示 NaN 值。如果sol.problem是 0 但后面value(V_sq)里出现 NaN通常不是求解器问题而是 YALMIP 变量没有全部进入约束集合某个孤立变量没有被任何表达式引用value之后是空。这时候去检查代码里是否有某个sdpvar变量定义了却没进入Constraints。CPLEX 参数作用项目中的建议值cplex.mip.tolerances.mipgap分支定界最优化间隙1e-4cplex.timelimit求解时间上限3600 秒cplex.threads并行线程数4 或 8cplex.mip.display求解日志显示频率2cplex.mip.limits.nodes最大节点数1e6 或留空5. 求解失败时CPLEX 日志和 check 函数怎么定位问题5.1 先看sol.info再看求解状态遇到Infeasible时很多人第一反应是扩大 CPLEX 容差这是错误方向。infeasible只说明约束集合本身没有交点和数值容差关系不大。项目日志里如果出现infeasible优先检查 OLTC 的抽头变量和 CB 的投切变量是否作用到母线电压上。比如 CB 只是定义出了Q_cb但没有写进节点无功平衡方程那这个变量孤立存在不会导致 infeasible反而是写进平衡方程但上下限冲突更容易触发。在实际项目中的排查顺序是先注释掉 ESS 的 SOC 递推约束改成单时段断面分别求解看每个时段是否可行。如果单时段可行、24h 不可行说明 SOC 初值或容量边界设置不合理。再把 OLTC 抽头范围从-4:4扩大到-8:8看是否依然 infeasible。若可行说明电压约束和变比范围冲突。最后检查潮流方程里的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 随机选了一个最优顶点。这时候调整目标权重比调求解参数重要。可以给动作惩罚设一个很小的正数只要大于数值噪声就能稳定输出平滑策略。日志特征真实含义项目里优先动作Infeasible约束无交集检查 ESS SOC 和 OLTC 档位边界Unbounded目标无下界检查风机无功和网损表达式方向Numerical difficulties数值刚度大调整lpmethod和标幺值Node limit exceeded分支定界节点爆炸放宽 mipgap 或抽头档位Integer optimal求到最优整数解读value前先看剩余间隙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 经常给出同时充放的病态解。本文还有配套的精品资源点击获取