这段时间在做10kV配电网的无功优化分布式光伏渗透率接近35%正午时段末端电压经常冲到1.07 p.u.原先靠并联电容器和OLTC档位组合调节响应慢、动作次数又老被考核。后来我把整套模型换成了基于二阶锥优化SOCP的配电网无功优化并把燃气轮机、储能、P2G这些电气综合能源设备一起放进MATLAB里做多目标优化问题从“能不能收敛”变成了“怎么把目标权重调合理”。这套方法的核心思路很简单把配电网潮流方程里的二次项做凸松弛得到可以全局求解的二阶锥规划再用MATLABYALMIPCPLEX这套组合把模型写出来。这篇文章写给正在做配电网优化、或者想把手头非线性无功优化模型升级成凸模型的工程师和学生。我会讲清楚SOCP背后的原理给出能跑通的代码骨架也会把我在实际算例里踩到的坑一条条列出来。1. 为什么配电网无功优化需要把二阶锥松弛当成正经方案1.1 传统非线性模型的两个难处配电网无功优化最底层的约束是潮流方程。无论是节点注入功率方程还是支路潮流方程电压和功率之间都是二次乘积关系整个可行域是非凸的。所谓非凸用大白话说就是你站在一个山谷里的某个位置往下看没法确定最低点到底在山谷深处还是更远的山背后局部最优和全局最优之间隔着算不完的山脊。我最早用传统内点法跑这个模型IEEE 33节点还好换到69节点就开始看运气了。初始点给得好十几秒出结果初始点稍微离谱一点海森矩阵奇异、迭代发散、报错信息千奇百怪。后来试过粒子群全局搜索能力确实强但每次运行结果都不一样同一份数据跑三次网损能差2%到5%写报告的时候很难交代为什么“最优值”是波动的。这不是算法不够努力而是非凸模型本身就没有稳定的求解路径。1.2 SOCP松弛为什么能“占便宜”二阶锥规划属于凸优化家族可行域是锥的凸集合目标函数只要是凸的局部最优一定等于全局最优。更重要的是内点法求解凸优化的收敛性有理论保证几十个变量的模型跑一遍结果就是确定的同一个算例换台电脑跑数值结果也可以复现。这在工程上太重要了。打个比方原来的配电网潮流约束像一片山坳SOCP松弛像把山坳填平成漏斗形状。填平过程中某些关键的山谷线被保留在漏斗壁上你顺着漏斗壁往底部走走到的就是原问题的最低点附近。这就是“精确松弛”的含义我们丢掉了等号把约束放宽但目标函数会把解拉回到边界上最后得到的解依然满足原问题。这里有个基本前提配电网通常是辐射状结构树状拓扑、单一平衡节点。对于辐射状配电网只要目标函数是电压、网损和电源成本的单调组合DistFlow的SOCP松弛在大多数工程场景下是精确的这也是近年大量配电网优化文献都采用SOCP的原因。1.3 什么时候该怀疑松弛结果SOCP松弛不等于万能。如果网架是环网或者目标函数里加了奇怪的常量约束松弛后的解可能落在“漏斗内部”也就是原问题不可行那就必须回到非线性模型或者先做网架解环处理。所以任何SOCP结果都要做一步松弛间隙验证后面我会专门讲怎么看这个间隙。拿到结果直接信是做优化的人最容易犯的错误。2. DistFlow潮流方程与SOCP约束的落地细节2.1 变量选错后面所有事都白做在MATLAB里写SOCP第一步不是往代码里堆约束而是选对变量。配电网无功优化我建议直接用DistFlow支路潮流模型变量取四组节点电压幅值平方 V_i单位 p.u.支路电流幅值平方 I_k单位 p.u.支路首端有功 P_k单位 p.u.或MW支路首端无功 Q_k单位 p.u.或Mvar。为什么不用电压相角和节点注入复功率直接建模因为一旦引入相角潮流方程里就出现cos和sin凸性彻底没了。V 和 I 取平方后约束里最多是二次乘积而且可以二次锥化。这是SOCP建模最关键的直觉所有讨厌的东西都藏在平方变量里先平方再锥化。2.2 支路潮流等式与节点功率平衡对任意支路 k首端节点为 i末端节点为 jDistFlow的三条主约束是有功平衡P_k - r_k * I_k 流向 j 的所有子支路有功之和 节点 j 的净有功负荷无功平衡Q_k - x_k * I_k 流向 j 的所有子支路无功之和 节点 j 的净无功负荷电压降方程V_j V_i - 2 * (r_k * P_k x_k * Q_k) (r_k^2 x_k^2) * I_k。其中 r_k、x_k 是支路电阻和电抗P_k、Q_k 取支路首端流向末端的方向。节点净负荷是负荷减去发电发电包括光伏、储能放电、燃气轮机P2G耗电则作为正负荷叠加进去。很多初学者把节点功率平衡写成“注入功率减去流出功率等于负荷”方向搞反后整个系统要么无解要么出现负网损。我写代码时习惯先画一张辐射状树的父子节点表再用find(from j)找子支路索引避免手推每条支路的连接关系。2.3 电流平方的二次等式怎么变成二阶锥DistFlow里原本有一个等式约束I_k (P_k^2 Q_k^2) / V_i这个式子是二次的而且带除法非凸。SOCP的标准做法是把它松弛成不等式|| [2 * P_k; 2 * Q_k; I_k - V_i] ||_2 I_k V_i这个形式的二阶锥约束在YALMIP里直接写Con [Con, cone([2*P(k,t); 2*Q(k,t); I(k,t)-V(i,t)], I(k,t)V(i,t))];为什么要写成这个形式因为把两侧乘开正好等价于 I_k * V_i P_k^2 Q_k^2。也就是说松弛后的可行域比原问题大多出来的东西是“电流平方大于等于真实值”的假解。由于目标函数里包含网损项 r_k * I_k优化器会本能地把 I_k 往下压最终 I_k 会贴到等号上假解被排除。这就是整个方法最精妙的地方目标函数充当了恢复等号的裁判。2.4 旋转锥转标准锥光伏逆变器和储能配电网里的分布式光伏通常要求逆变器具备无功调节能力而逆变器容量约束是S_dg^2 P_dg^2 Q_dg^2这同样是一个非凸约束但可以转成标准二阶锥。设 S_ref 为逆变器额定视在功率则写成cone([2P_dg; 2Q_dg; S_dg - S_ref], S_dg S_ref)实际代码里我会把 S_dg 当成一个辅助变量S_ref 取常数然后放一堆这个锥约束。储能充电功率和放电功率的平方约束处理方式完全相同。只要看到“某容量值大于等于两个分量平方和”这类结构就往标准锥上套。3. 多目标无功优化的目标设计与权重求解3.1 四个常用目标怎么取舍无功优化的目标往往不止一个。我常用的有四个配电网有功网损最小反映运行经济性节点电压偏移最小反映电压质量综合运行成本最小包括从上级网购电成本、燃气轮机燃料成本、储能退化成本弃光惩罚最小让分布式光伏尽量多发电。理论上可以把四个目标叠成一个加权和但工程上我不建议上来就堆四个目标。目标越多权重越难解释做敏感性分析时你会被各种反直觉的结果折腾到怀疑人生。我一般只保留三个网损、电压偏移、运行成本弃光惩罚可以放在电网侧购电成本里体现光伏出力被削减等同于收入损失。3.2 加权和法的量纲统一操作加权和法是把多目标问题变成单目标问题最简单的方法但直接写成 w1Ploss w2Vdev w3*Cost 是新手最爱踩的坑。网损一般在小数点后几位运行成本可能是几千上万直接加权后网损和电压偏移基本被成本项吞掉权重怎么调都感觉不灵敏。我处理这个问题比较粗暴先分别单独最小化每个目标得到各个单目标最优值然后每个目标项都除以自己的单目标最优值。这样三项都在1附近起步权重才有公平可言。%% 先求两根基准线 opts sdpsettings(solver,cplex,verbose,0); optimize(Con, Ploss, opts); base_Ploss value(Ploss); optimize(Con, Vdev, opts); base_Vdev value(Vdev); optimize(Con, Cost, opts); base_Cost value(Cost); %% 三目标加权优化 w [0.4 0.3 0.3]; obj w(1)*Ploss/base_Ploss w(2)*Vdev/base_Vdev w(3)*Cost/base_Cost; optimize(Con, obj, opts);如果某个单目标最优值是0比如理想情况下电压偏移可以做到0那就不能直接除。我的处理办法是不用绝对零目标而是给电压偏移加一个很小的参考基准或者用 (Vdev - Vdev_min) 作为实际目标项总之不要让分母出现0。3.3 帕累托前沿扫描与可视化想观察目标之间的冲突关系可以固定两个目标扫描权重画帕累托前沿。以网损和电压偏移为例权重从0到1每间隔0.05取一个点wList 0:0.05:1; Pareto zeros(length(wList), 2); for i 1:length(wList) obj wList(i)*Ploss/base_Ploss (1-wList(i))*Vdev/base_Vdev; optimize(Con, obj, opts); Pareto(i,1) value(Ploss); Pareto(i,2) value(Vdev); end plot(Pareto(:,1)*1000, Pareto(:,2)*100, o);从帕累托前沿上能直观看到“电压调得越狠网损跟着涨”这样的冲突关系这比单点优化结果强得多。需要注意的是加权和法在帕累托前沿是凹形状时取不到中间点如果发现曲线中间明显断开就改用电约束法一次固定一个目标阈值去最小化另一个目标。配电网无功优化里前沿多数情况比较平滑加权和够用。4. 电气综合能源系统的设备级建模扩展4.1 燃气轮机和热电联产机组电气综合能源系统里和配电网关系最紧密的设备是燃气轮机和热电联产机组。燃气轮机输出有功功率 P_gt无功输出受功率因数约束-abs(P_gt) * tan(acos(pf_min)) Q_gt abs(P_gt) * tan(acos(pf_min))由于 P_gt 是非负变量这条约束是线性的不会破坏凸性。燃料消耗按有功出力线性化写成F_fuel a_gt * P_gt b_gt热电联产机组可以把发电余热用于热负荷热出力和电出力近似线性关系再配上电锅炉或者热泵就能在模型里加入热平衡约束。注意我这里做的是“设备级耦合”而不是把整个天然气管网和热力管网都建模进来。在无功优化这个主题下管网方程会让问题规模急剧膨胀而且会让SOCP结构变得不干净。工程上先把电气耦合设备当成可调电源和可调负荷接入电力节点是性价比最高的做法。4.2 储能SOC状态方程和同时充放电问题储能是时序优化里最麻烦也最有价值的一环。SOC递推约束必须写soc(t1) soc(t) eta_ch * Pch(t) - Pdis(t) / eta_dis其中 eta_ch、eta_dis 是充放电效率Pch、Pdis 都是非负变量。为了限制同一时刻不能既充电又放电严格做法是引入0-1整数变量变成混合整数二阶锥规划CPLEX和Gurobi都能解但求解时间会明显上涨。我在工程里通常用连续松弛完全去掉二进制变量只保留 Pch 0、Pdis 0然后目标函数里给充放电加一个很小的成本系数。因为同时充放电只会让SOC白白损耗效率带成本目标的情况下优化器天然会避免实际解里同时充放电的时段基本不会出现。用这个方法可以把整个模型维持成纯SOCP24小时时段的求解速度会快很多。如果后续要跟调度部门对接采纳严格约束再升级成MISOCP也不难。储能接入节点后节点净有功负荷里要减掉 Pdis加上 Pch。电压质量紧约束的末端节点挂储能效果最明显这是无功电压控制里“有功调节也能辅助电压”的典型案例。4.3 P2G与电锅炉如何进入节点功率平衡P2G设备消耗有功功率 P_p2g产生天然气可以用效率 eta_p2g 折算成天然气产出。在无功优化模型里不需要把天然气管网潮流建模出来只需要把 P_p2g 视作可调电负荷放进对应节点的有功平衡方程。类似地电锅炉消耗有功功率产生热能热出力满足H_eb eta_eb * P_eb热负荷平衡约束写成H_eb H_chp H_load这里的 H_chp 是热电联产供热部分由电出力换算得到。这套设备级建模的核心思想是电气综合能源系统的耦合点最终都落在电力节点的有功注入上所有设备的无功能力要么来自功率因数范围要么来自逆变器或机组的无功限额全部可以写成线性或二阶锥约束不会破坏SOCP的整体结构。5. MATLAB代码骨架与松弛间隙验证5.1 程序主结构与数据组织我用YALMIP作为建模层求解器用CPLEX。整个程序不建议写成一个两千行的单文件而是拆成四段数据准备、变量声明、约束组装、求解后处理。数据准备部分包括节点表、支路表、负荷曲线、光伏曲线、储能参数、燃气轮机参数。节点表和支路表是核心我建议把from、to、r、x四个向量严格按支路索引对齐别用手填矩阵的方式存拓扑后面所有循环都要靠索引向量取数。5.2 核心建模片段YALMIP下面这段是核心建模骨架已经能体现DistFlow和SOCP松弛的全部逻辑。跑通之前需要把from、to、r、x、Pload、Qload等数据准备好并定义光伏、储能、燃气轮机的注入变量。%% distflow_socp_core.m nb 33; % 节点数 nt 24; % 调度时段 nl 32; % 支路数 V sdpvar(nb, nt); % 节点电压幅值平方 I sdpvar(nl, nt); % 支路电流幅值平方 P sdpvar(nl, nt); % 支路首端有功 Q sdpvar(nl, nt); % 支路首端无功 Ppv sdpvar(nb, nt); % 光伏有功注入 Qpv sdpvar(nb, nt); % 光伏无功注入 Pgt sdpvar(nb, nt); % 燃气轮机有功出力 Qgt sdpvar(nb, nt); % 燃气轮机无功出力 Pdis sdpvar(nb, nt); % 储能放电 Pch sdpvar(nb, nt); % 储能充电 soc sdpvar(nb, nt1); % 荷电状态 Con []; for t 1:nt for k 1:nl i from(k); j to(k); kids find(from j); % 以j为首端的子支路 if isempty(kids) child_P 0; child_Q 0; else child_P sum(P(kids,t)); child_Q sum(Q(kids,t)); end Pgen Ppv(j,t) Pdis(j,t) - Pch(j,t) Pgt(j,t); Qgen Qpv(j,t) Qgt(j,t); Con [Con, P(k,t) - r(k)*I(k,t) child_P Pload(j,t) - Pgen]; Con [Con, Q(k,t) - x(k)*I(k,t) child_Q Qload(j,t) - Qgen]; Con [Con, V(j,t) V(i,t) - 2*(r(k)*P(k,t) x(k)*Q(k,t)) ... (r(k)^2 x(k)^2)*I(k,t)]; Con [Con, cone([2*P(k,t); 2*Q(k,t); I(k,t)-V(i,t)], ... I(k,t)V(i,t))]; end Con [Con, V(1,t) 1.0]; % 根节点电压固定为1.0 p.u. end %% 目标函数网损 电压偏移 运行成本 Ploss sum(sum(r .* I)); Vdev sum(sum((V - 1).^2)); Cost sum(sum(cgrid .* Pgrid cgt .* Pgt cbat .* (Pch Pdis))); obj w1*Ploss w2*Vdev w3*Cost; opts sdpsettings(solver,cplex,verbose,1); optimize(Con, obj, opts);代码里的Pgrid是从根节点出线功率之和构造方式可以直接用sum(P(find(from1),t))。节点功率平衡里对根节点没有做强制等式因为变电站母线等同于上层电网的松弛母线注入量由优化自由决定成本项会自然约束它不要太大。5.3 松弛间隙到底该看什么指标SOCP模型求出来的解第一件事不是看目标函数而是看松弛到底紧不紧。最简单实用的指标是每条支路的SOCP松弛残差res_k || [2P_k; 2Q_k; I_k - V_i] ||_2 / (I_k V_i)理论上这个值接近0代表 I_k * V_i 约等于 P_k^2 Q_k^2也就是松弛解正好落回原潮流等式上。计算方式tol_soc 0; for k 1:nl res norm([2*value(P(k,end)); 2*value(Q(k,end)); ... value(I(k,end))-value(V(from(k),end))]) ... / (value(I(k,end)) value(V(from(k),end))); tol_soc max(tol_soc, res); end fprintf(SOCP最大松弛残差: %.3e\n, tol_soc);残差小于1e-5基本可以认为松弛是精确的小于1e-3勉强可用超过1e-2就要警惕。出现大的残差优先检查是不是目标函数里没有网损项或者电压调节项导致优化器没有动力把 I_k 压回边界。5.4 一组实测量级我在一台i5-9400、16GB内存的机器上用MATLAB R2019b加CPLEX 12.10跑过的量级大致如下算例节点数时段数决策变量规模求解时间最大SOCP松弛残差15节点1524约13001s到2s1e-6量级33节点3324约450010s到20s1e-6量级69节点6924约1000090s到180s1e-5量级这些时间是把YALMIP建模时间也算进去的。如果加了储能二进制变量和燃气轮机爬坡约束33节点问题运行时间可能翻到两三分钟做权重扫描时长会很难受所以要提前规划好算例规模。6. 求解器选型、参数缩放与排查经验6.1 求解器怎么选CPLEX、Gurobi、SeDuMi纯SOCP问题在CPLEX、Gurobi、MOSEK、SeDuMi、SDPT3里都能解但性能差异很大。我的经验是CPLEX和Gurobi对大算例的数值稳定性明显好SeDuMi在5000个变量以上的模型里经常出现收敛慢或数值警告。如果自己是学生并且拿不到商业求解器授权先用YALMIP加SeDuMi把15节点的模型跑通理解建模流程然后申请Gurobi学术授权跑大算例这是最务实的路线。混合整数SOCP就完全不同SeDuMi和SDPT3不支持整数变量必须用CPLEX、Gurobi或者MOSEK。所以只要模型里加了储能同时充放电限制的0-1变量求解器选型基本就锁死了商业求解器。6.2 单位与缩放问题这个坑非常隐蔽。我在第一次跑69节点模型时目标函数里网损用kW成本用元约束里电压用kV结果CPLEX一直警告“problem is badly scaled”求解时间暴涨。后来把所有功率统一换算成p.u.电压统一用p.u.平方成本系数按基准功率折算后同样的问题求解时间降了一个数量级。建议每写完一个约束用YALMIP的check(Con)函数看一眼各约束残差量级。约束残差分布如果在1e-6到1e-4之间说明缩放基本健康如果出现1e2甚至1e3量级的残差一定有你忘记归一化的物理量。6.3 不可行和收敛慢的排查顺序模型无解的时候不要急着往约束里加惩罚项先按顺序排查检查根节点电压约束。V(1,t)1.0 一般没问题但如果你同时给根节点电压设置了上下限1.05和0.95就会互相打架。检查储能SOC末态约束。很多时序优化要求末时段SOC回到初值这在负荷和光伏曲线波动剧烈时可能无解改成软约束加惩罚项让末态尽量靠近初值。检查节点功率平衡方向。Ppv、Pdis、Pch、Pgt的符号和Pload之间的正负关系必须一致这是最常见的低阶错误。检查光伏无功范围。光伏逆变器无功上限和当前有功出力有关不能给一个固定的大常数否则可能超出逆变器视在功率锥约束。6.4 时序优化的提速建议做24小时甚至96时段优化时建议先用15节点或33节点算例把参数调通再上大网络。我在做69节点24时段储能燃气轮机的模型时一开始权重扫描设了21个点点与点之间还有耦合变量跑了两个多小时才出完帕累托前沿。后来改成三步走固定一组权重先用33节点跑通全时段把约束容差和最优性容差保持默认不追求1e-9的精度权重扫描时每次重新求解前删除上一次的YALMIP变量缓存并在sdpsettings里设置cplex.timelimit, 300避免单个点卡死。最后再分享一条经验我吃过最大的亏不是模型写不对而是拿一个没做松弛间隙检验的SOCP结果直接汇报。配电网优化里的SOCP只是一个工具箱松弛是否精确必须检查否则你拿到的是一个放缩解物理上可能根本不可行。先跑通15节点做残差检验再上33节点和电气综合能源设备这条路最稳。