把一套含电、气、热三类能源的综合能源优化程序从“能跑”调到“跑得稳”我前后折腾了小半年。最典型的教训是同样的园区数据第一版用非线性求解器直接算潮流方程初值稍微给偏一点CHP出力的结果就能差出15%第二版把潮流方程线性化速度快了但最优解代入实际潮流一验算线路过载都发现不了第三版终于走到正路上来——用MATLAB加Yalmip建模把电网DistFlow方程里的电压-电流二次项做二阶锥松弛天然气Weymouth方程做凸松弛热网在质调节假设下线性化整套模型交给Cplex做混合整数二阶锥规划MISOCP不猜初值、收敛稳定、全局最优有保障。这篇文章就把第三版的路子从头到尾拆开物理建模、凸松弛原理、Yalmip代码、求解器调参、松弛紧性校验按你拿到一个问题后实际操作顺序来写适合正在做园区综合能源、微电网调度和能源互联网方向课题的同学也适合想从传统非线性求解转凸优化的工程师。1. 综合能源优化问题的骨架三条能源网络的耦合与调度边界1.1 优化目标不是单一的电费而是整个能源站的总账很多第一次接触综合能源优化的人会把目标函数直接写成“购电费用最小”。这其实只对了一半。含电气热的系统里天然气的购买费用、CHP和锅炉的启停损耗、设备的运行维护成本甚至弃风弃光惩罚全部都要进目标函数。不然很容易出现一种奇怪结果电费降下来了气费涨得更多总账反而难看。以最典型的日前经济调度为例目标函数长这样Objective sum(price_e .* P_buy * dt) ... % 购电费用 sum(price_g .* F_gas * dt) ... % 购气费用 sum(C_start .* u_start) ... % 启停成本 sum(C_om .* P_device); % 运维成本其中dt是各时段长度单位小时price_e是分时电价price_g是气价。注意这里有个隐蔽细节如果电价单位是“元/MWh”电功率是MW时间长度是h乘出来直接就是元但如果电价表给的是“元/kWh”很多人直接代进去目标函数数值一下子就差了1000倍求解器大概率直接报数值问题。1.2 CHP、锅炉、电转气物理设备如何把三条网“缝”起来电气热三个网络不是各算各的它们靠耦合设备连在一起。最常见的三样CHP热电联产天然气进去电和热同时出来。电出力与热出力之间有一个可行的配比范围不是随便组合。燃气锅炉天然气进去只出热。结构简单效率相对固定。电锅炉/热泵电进去出热。在弃风弃光时段特别好用相当于把多余的电转化成热储存或直接供热。如果系统里还有P2G电转气那就是第四条耦合链路电变成天然气回灌到气网。做这类课题时我把这些设备想象成“配货中转站”——电、气、热三条传送带各自有容量限制CHP就是那个从气传送带取货、同时往电和热传送带放货的机器但放货比例有约束。二阶锥规划要干的事就是在所有传送带容量和转运规则都满足的前提下找到总成本最低的配货方案。1.3 为什么只做能量平衡不够还要带上网络约束这是很多简化模型最容易被诟病的地方。如果只写“电功率平衡 热功率平衡 气源平衡”相当于默认所有线路和管道都是无限容量、没有压降、没有电压限制。这种模型求出来的解拿到真实物理网络里很可能是不可行的CHP在某时段满出力局部节点电压可能越上限热网末端用户可能因为管道压降不够而供热不足。所以完整的综合能源优化必须把网络约束写进优化模型。电网要有电压约束、支路容量约束气网要有节点压力约束、管道流量与压力关系热网要有供回水温度和节点热功率约束。但这恰恰引出了一个问题这些物理约束几乎都是非线性的直接扔给求解器就会回到我最开始说的初值敏感、局部最优的泥潭。这也是二阶锥规划在这类问题里成为主流选择的核心原因。2. 为什么偏偏是二阶锥凸松弛的数学依据与适用边界2.1 DistFlow的非凸等式与“松弛为锥”的完整推导配电网最常用的潮流模型是DistFlow适用于辐射状网络。对一条支路(i, j)用P_ij、Q_ij表示支路首端有功、无功功率v_i是节点i电压幅值的平方l_ij是支路电流幅值的平方r_ij、x_ij是支路电阻、电抗那么电压降落方程v_j v_i - 2(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*l_ij电流定义方程l_ij (P_ij^2 Q_ij^2) / v_i问题出在最后这个等式。它不是凸的直接放进优化模型会得到一个非凸可行域求解时对初值极其敏感。二阶锥松弛的做法就是把这个等式“放松”成不等式l_ij (P_ij^2 Q_ij^2) / v_i这个不等式等价于标准二阶锥|| [2*P_ij; 2*Q_ij; v_i - l_ij] ||_2 v_i l_ij不信可以自己推一下两边平方展开以后就是P_ij^2 Q_ij^2 v_i * l_ij正好是上面的式子。这里的直觉可以这样理解原约束是一条弯曲的非凸曲线松弛成锥以后可行域变大了。如果求出来的最优解恰好落在锥的“边界”上也就是不等式取等号那么这个解和原非凸问题完全一致我们说松弛是“紧”的。文献里的精确性定理告诉我们对辐射状配电网、目标函数合理、没有特殊逆向潮流场景的绝大多数情况SOCP松弛都是紧的。工程项目的实践也基本符合这个结论。2.2 天然气Weymouth方程和热网约束的凸化处理天然气管道稳态流量和节点压力之间最常用的是Weymouth方程f_pipe^2 C^2 * (alpha_i - alpha_j)其中alpha_i是节点i压力的平方C是管道常数。这是一个非凸二次等式。工程做法同样是松弛成不等式f_pipe^2 C^2 * (alpha_i - alpha_j)再借助一个漂亮的变换写成标准二阶锥d alpha(i) - alpha(j); Constraints [Constraints, d 0]; Constraints [Constraints, cone([2*f_pipe(k); C2 - d], C2 d)];C2是C^2。这个写法的好处是无论Yalmip还是Cplex/Gurobi都把它当成真正的二阶锥约束来处理而不是回调非线性求解器。需要说明的是气网流向如果可能反转上述写法就不够用一般要引入二进制变量刻画方向模型复杂度会明显上升。实际项目里通常可以根据气源和负荷位置预判流向或者只在少数管道上保留方向变量。热网是三者里最麻烦的。热功率本身是流量和温差的乘积H c_p * m * (T_s - T_r)。如果同时把流量m、供水温度T_s、回水温度T_r都当变量这就是双线性项直接突破SOCP框架。工程上常用两条路绕开一是质调节把流量固定在设计工况值只优化温度二是量调节温度固定只优化流量。两种做法都能让热网约束保持在线性体系内。如果确实需要两者同时可变我建议用“固定流量解一次SOCP再把解出来的流量代回去更新温降系数重新求解”的外层迭代通常两三次就能收敛。这一点后面建模部分会细说。2.3 SOCP与LP、NLP方案的取舍不同建模方式没有绝对优劣只有适合和不适合。我自己的使用感受可以整理成下面这个表建模方式网络精度问题是否凸初值依赖全局最优求解速度线性化/直流近似低是无有秒级二阶锥松弛高是无有秒到分钟直接求解非线性等式高否强不保证分钟到不收敛所以对含电气热的综合能源优化我的默认选择是连续变量走SOCP需要机组启停就加0-1变量变成MISOCP。只有在SOCP实在没法表达某些强非线性设备特性时才考虑对那一小块做分段线性化或外层迭代。3. 电-气-热联合建模从物理方程到代数约束的完整推演3.1 电网侧建模与标幺化电网建模第一步是统一单位。我习惯把功率统一成MW电压统一成标幺值pu。电压幅值的平方v一般限制在0.95到1.05的平方之间也就是0.9025到1.1025这个范围比直接用电压百分数要直观得多。支路数据和DistFlow约束的代码骨架是这样的nb 5; nl 4; branch [1 2; 2 3; 2 4; 4 5]; % from, to r [0.01; 0.02; 0.015; 0.01]; x [0.03; 0.05; 0.04; 0.03]; v sdpvar(nb, 1); l sdpvar(nl, 1); Pij sdpvar(nl, 1); Qij sdpvar(nl, 1); Constraints []; for k 1:nl i branch(k,1); j branch(k,2); Constraints [Constraints, ... v(j) v(i) - 2*(r(k)*Pij(k) x(k)*Qij(k)) ... (r(k)^2 x(k)^2)*l(k)]; Constraints [Constraints, ... cone([2*Pij(k); 2*Qij(k); v(i)-l(k)], v(i)l(k))]; end Constraints [Constraints, v 0.95^2, v 1.05^2, l 0];节点功率平衡我建议用关联矩阵一次性写避免手写循环把方向搞混。核心逻辑是流入节点的功率加上该节点注入购电、分布式电源、CHP电出力等于流出节点功率加该节点负荷。这里不展开完整代码但提醒一句branch矩阵里的方向定义必须和功率平衡公式保持一致否则算出来的潮流符号会全反。3.2 气网侧建模与压降方程气网变量主要是节点压力平方alpha、管道流量f_pipe、气源注入量F_source和节点气负荷F_load。节点气负荷来自CHP的气耗和燃气锅炉的气耗。约束就那么几类气源点注入有上限节点压力有上下限管道Weymouth约束用上一章的SOCP写法治压缩机如果存在耗气量常简化成流量的一定比例或者干脆用分段线性近似压缩机的升压比约束一般也是线性的。这样做完整个气网约束保持在SOCP体系内不会引入非线性求解需求。3.3 热网侧建模双线性项的两个务实降阶思路热网建模比电网、气网更“工程”一些。完整的动态热网模型涉及管道热延迟稳态优化一般不考虑延迟只考虑节点热功率平衡和温度混合关系。节点热功率约束是H_node c_p * m * (T_s - T_r)c_p取4.182 kJ/(kg·℃)流量m单位kg/s温差单位℃算出来热功率单位是kW再除以1000才是MW。这是单位最容易翻车的地方之一。如果流量和温度都是变量这就不再是SOCP。我在工程里最常用的做法是质调节流量按设计工况固定下来那么热功率约束变成线性约束T_s、T_r是决策变量H_node也是变量。温度还有上下限约束比如供水温度不超过120℃回水温度不低于40℃。管道温降的指数项我会做成分段线性近似节点温度混合约束本质是“流入该节点所有管道热量之和等于流出热量之和”在流量固定以后同样线性化。如果你不想固定流量那就用外层迭代。第一次按设计流量解SOCP拿出各支路流量后重新计算温降系数和混合系数再解第二次。实测下来两三次迭代足够稳定而且每轮都是凸问题不会出现NLP那种中途发散的情况。3.4 耦合设备可行域与目标函数CHP是最需要认真建模的设备它的电出力和热出力之间存在一个二维可行域。精确模型是复杂的非线性区域工程上常用多边形近似用一组线性不等式描述Constraints [Constraints, P_chp Pmin, P_chp Pmax]; Constraints [Constraints, H_chp Hmin, H_chp Hmax]; Constraints [Constraints, P_chp alpha_chp.*H_chp beta_chp];alpha_chp、beta_chp来自可行域多边形顶点的拟合。加上启停变量后变成u_chp binvar(n_chp, 1); Constraints [Constraints, Pmin*u_chp P_chp Pmax*u_chp]; Constraints [Constraints, Hmin*u_chp H_chp Hmax*u_chp];燃气锅炉的模型就简单很多H_gb eta_gb * F_gb效率常数。电锅炉是H_eb eta_eb * P_eb。P2G是F_p2g eta_p2g * P_p2g。这些设备方程全部线性。目标函数方面除了前面说的购电购气成本还要考虑设备启停和运维。这里提一个常见问题目标函数里如果没有网损项或惩罚项松弛锥可能“偷懒”导致SOCP解虽然可行但不紧。解决办法后面专门讲建模阶段先记住目标函数最好加一个极小系数的网损惩罚项比如1e-5 * sum(r.*l)能显著改善松弛紧性。4. MATLABYalmip程序实现变量、约束与求解接口4.1 数据组织和决策变量声明写这类程序我强烈建议把所有网络参数装进struct或者独立脚本不要在代码里到处写魔法数字。下面是一个单时段的示例多时段只需把变量扩展成24列约束按t1:T循环即可。% 统一单位电功率MW电压pu温度℃压力MPa % 电网参数 grid.nb 5; grid.branch [1 2; 2 3; 2 4; 4 5]; grid.r [0.01; 0.02; 0.015; 0.01]; grid.x [0.03; 0.05; 0.04; 0.03]; grid.P_load [0; 0.2; 0.15; 0.25; 0.1]; % 决策变量 v sdpvar(grid.nb, 1); l sdpvar(size(grid.branch,1), 1); Pij sdpvar(size(grid.branch,1), 1); Qij sdpvar(size(grid.branch,1), 1); P_buy sdpvar(1, 1); P_chp sdpvar(1, 1); H_chp sdpvar(1, 1); H_gb sdpvar(1, 1); F_gas sdpvar(1, 1);变量名用Pij这种带方向的东西没问题但要注意Yalmip里sdpvar默认就是全变量不需要加full参数。有些老版本写法加了也不报错但新版本可能把你当成矩阵对称去处理坑过不少人。4.2 电网DistFlow与SOCP约束组装电网约束组装我在3.1已经给出了代码这里补充节点功率平衡的写法。推荐用关联矩阵一次性构造% 关联矩阵 incidence(j,k)支路k从节点j流出取1流入节点j取-1 incidence zeros(grid.nb, size(grid.branch,1)); for k 1:size(grid.branch,1) incidence(grid.branch(k,1), k) 1; % from incidence(grid.branch(k,2), k) -1; % to end % 节点注入 P_buy P_chp按实际接入位置映射 P_inj zeros(grid.nb, 1); P_inj(1) P_buy; % 节点1是上级电网连接点 P_inj(3) P_chp; % CHP接在节点3 Constraints [Constraints, incidence*Pij P_inj grid.P_load];这里的符号约定是incidence(j,k)1表示支路功率从j流出-1表示流入j。节点功率平衡写成“流出减流入加注入等于负荷”。实际项目里还要加无功平衡方法和有功一致把Qij和Q_load放进去。如果只做有功调度无功可省略但电压约束会失真所以一般至少保留无功潮流。4.3 气网、热网和设备约束的写法气网的SOCP约束写法前面已经给了。热网在质调节假设下的约束写法如下% 质调节流量固定为设计工况 m_fixed T_s sdpvar(nh, 1); T_r sdpvar(nh, 1); H_node sdpvar(nh, 1); % 节点热功率 for n 1:nh Constraints [Constraints, ... H_node(n) 4.182e-3 * m_fixed(n) * (T_s(n) - T_r(n))]; Constraints [Constraints, T_s 70, T_s 120]; Constraints [Constraints, T_r 40, T_r 70]; endCHP可行域和启停约束、燃气锅炉和电锅炉的效率方程按3.4的写法直接拼进Constraints。这里特别提醒如果H_node是变量并且它同时出现在热网节点平衡和CHP热出力约束里不要重复声明两个不同变量去“相等”直接在热网节点平衡里使用H_chp和H_gb本身少一层中间变量就少一批数值问题。% 热网节点热平衡源节点注入等于负荷 Constraints [Constraints, H_chp H_gb sum(H_load)];4.4 求解设置与结果提取模型组装完成后求解设置是我的固定套路ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.param.mip.tolerances.mipgap 1e-3; ops.cplex.param.timelimit 3600; sol optimize(Constraints, Objective, ops); if sol.problem 0 fprintf(求解成功总成本 %.2f 元\n, value(Objective)); Pij_opt value(Pij); v_opt value(v); l_opt value(l); else disp(sol.info); endvalue()是Yalmip提取变量数值的标准接口。多时段程序里P_buy是24维向量value(P_buy)直接就得到24个时段的购电功率曲线画图非常方便。5. 求解器选型与参数调优从“能算”到“算得快”5.1 为什么默认Cplex/Gurobi而不是SedumiYalmip背后支持很多求解器但如果综合能源模型里出现了binvar问题就变成MISOCP。SDPT3和Sedumi这类纯内点法求解器只能处理连续SOCP不认二进制变量。Cplex和Gurobi对MISOCP的支持非常成熟分支定界、割平面、锥检测一条龙基本不用你操心。如果你的问题没有整数变量那么Mosek也很快但从生态和接口稳定性来看我还是优先推荐Cplex或Gurobi。需要提醒的是Yalmip默认会自动挑一个它能探测到的求解器。如果装了解算器但Yalmip没识别出来先运行yalmiptest检查接口是否配好。很多“Solver not found”其实是mex接口没安装不是模型问题。5.2 数值缩放是最容易被低估的一步这条经验值钱。电气热综合能源系统里电压平方可能只有0.95热网温度可能是一百多天然气压力如果用了帕斯卡单位就是几百万这几个数量级差到天上去了。Cplex内点法做锥约束检测的时候对数值尺度极其敏感。变量量级差距超过1e6常见的表现就是报Numerical problems或者明明有可行解却判成不可行。我踩过最惨的一次就是天然气节点压力用了Pa做单位模型怎么调都不可行后来统一改成MPa同一个模型秒解。从那以后我的固定习惯是电功率全部MW电压全部pu温度全部℃压力全部MPa成本全部万元或元任何变量尽量压在1e-3到1e3这个区间。5.3 常用求解器参数与中等规模算例表现Cplex的常用参数ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.param.mip.tolerances.mipgap 1e-3; % MIP间隙 ops.cplex.param.mip.tolerances.integrality 1e-5; % 整数容差 ops.cplex.param.threads 8; % 并行线程 ops.cplex.param.timelimit 3600; % 时间上限Gurobi对应的是ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 1e-3; ops.gurobi.TimeLimit 3600; ops.gurobi.NumericFocus 2;NumericFocus是我个人很喜欢的一个参数数值有问题时调到2或3求解器会更谨慎地处理矩阵计算代价是速度慢一点但总比不可行强。中等规模算例电网100多个节点、气网30个节点、热网20个节点、24时段含80个左右0-1变量在我自己的测试环境下纯连续SOCP一般十几秒能出结果加整数变量后往往要到几分钟。如果你遇到几小时跑不完的情况先别改参数回头检查是不是给太多设备加了0-1变量——很多设备在调度周期内根本没有启停需求固定成1就行能砍掉大量分支。6. 松弛紧性检验与排错实战6.1 残差检查判断二阶锥松弛是否“紧”SOCP解出来后不能直接拿去汇报。你要先判断松弛是不是紧的。判断方法很简单检查每个支路是不是满足原来的等式约束。res_max 0; for k 1:nl i branch(k,1); j branch(k,2); res_k abs(v_opt(i)*l_opt(k) - (Pij_opt(k)^2 Qij_opt(k)^2)); res_max max(res_max, res_k); end fprintf(最大锥残差: %.3e\n, res_max);在标幺制下残差如果小于1e-5基本可以认为松弛是紧的解可以直接用。如果残差到了1e-2甚至更大说明某些支路被松弛得很松最优解虽然SOCP可行但放到原物理方程里根本不对。最有效的补救办法在目标函数里加一个很小系数的网损惩罚项1e-5 * sum(r.*l)。原因很朴素目标函数里没有网损项时潮流变量对模型来说是“免费的”求解器怎么松快怎么来加上网损惩罚后它就有动力把l压在真实物理值附近。气网也要做同样检查看f_pipe^2和C^2*(alpha_i - alpha_j)的偏差。我见过不少项目电网上来就检查气网完全忽略最后结果里气网管道流量和压力的关系对不上整个调度方案直接被否掉。6.2 Infeasible和Numerical problems的排查链路遇到不可行我有一套固定的排查顺序比盲目调参有效率得多第一步把约束按模块拆开只留电网约束求解再分别只留气网、热网约束确认每个模块单独都是可行的。这个步骤能快速定位问题是出在哪一圈。第二步检查单位。这是最频繁的幕后黑手。尤其是热功率方程里那个4.182的系数用kJ/(kg·℃)还是kW·h/(kg·℃)结果差3600倍。压力用Pa还是MPa又差1e6倍。一次不可行排查先花十分钟把所有单位统一往往就解决了。第三步检查负荷有没有“无源节点”。比如某个热负荷节点边上没有热源、没有管道连接热量进不去模型当然不可行。电网里如果出现孤立负荷节点同理。第四步如果模型是MISOCP且报Numerical problems把整数约束临时改成连续变量看SOCP能不能解。连续能解、加整数就出问题多半是数值缩放不良参考5.2去处理。6.3 常见报错速查表下面这张表是从我自己的踩坑记录里整理出来的按出现频率排序现象常见原因处理办法Solver not found求解器mex接口没配置运行yalmiptest检查Yalmip接口Numerical problems变量量级差距大、单位不统一统一为MW/pu/℃/MPa缩放模型Infeasible problem约束冲突、单位换算错、孤立节点按模块拆解逐层定位求解时间过长0-1变量太多、约束过紧检查启停变量必要性调大MIPGap松弛残差偏大目标函数缺少潮流项惩罚加网损惩罚项如1e-5*sum(r.*l)这里再补充一个Yalmip特有的坑如果把sdpvar变量直接赋值给普通double矩阵MATLAB会报类型错。遇到这种情况不是模型错是数据结构设计问题把存储容器改成cell或者直接重声明变量就好。6.4 结果验证与非线性模型对照的实用技巧SOCP结果出来以后我习惯做一次“物理后验”而不是直接信求解器。做法是把求出来的P、Q、v代入原始非凸DistFlow方程看电压降落方程是不是成立潮流等式是不是吻合。这不是重新求解非凸问题只是验证速度快也不涉及初值。气网和热网同理。如果项目允许我会再用非线性求解器跑一个不要求全局最优的可行解来对照看成本和SOCP结果差多少。如果SOCP成本远低于NLP可行解大概率是SOCP松弛不紧或者某些设备可行域被过度放大了。反过来如果两者成本接近那说明松弛质量很好结果可以放心用。还有一个容易被忽略的工程检查看热网温度结果是否落在物理合理区间。比如某个时段热负荷为0但求解器给出的CHP仍在供热那基本可以断定热网模型少加了CHP最小热出力约束或者热负荷平衡那里漏了用能下限。这种“数学可行、工程怪异”的解靠残差是看不出来的必须结合物理直觉去审视出力曲线。最后说点个人习惯。我现在接到电气热综合能源优化程序第一件事不是写代码而是先画三张网络图把所有设备的耦合关系列成一张表统一好单位再动手。很多“程序跑不起来”的坑其实都是物理关系没理顺、单位不一致导致的。另外一个非常实用的小技巧先跑纯SOCP不带整数的连续模型确认收敛和残差指标合格之后再加启停二进制变量。这样能很快区分是建模错误还是求解器数值问题。如果你也卡在“同样的模型人家收敛我不收敛”这种问题上先从这一步开始排查大概率有奇效。
