这两年做电-气互联系统优化调度的人明显多了原因大家都清楚碳中和从一句口号变成硬约束之后电网调度里再也不能把碳排放当空气。而我这个课题最核心的一句话是有功和无功必须放在同一个优化问题里一起算。单纯做经济调度或者单纯做无功优化在新能源渗透率不高的年代勉强够用可一旦风电、光伏大量接入系统电压问题跟有功潮流强耦合分开算出来的结果经常是次优的甚至在边界工况下直接凑不出可行解。这篇博客把我在碳中和目标下搭建电气互联系统有功-无功协同优化模型的全过程拆开讲模型怎么建、约束怎么处理、Matlab里用YALMIP加Cplex怎么调通以及我实际跑代码时踩过的一堆坑。适合正在做综合能源系统优化、无功优化方向的研究生或者想把手里的经济调度代码升级成电-气协同碳成本版本的工程师。1. 项目概述为什么非要做这个模型1.1 电气互联系统到底指什么先把这个名词掰清楚。这里的电气互联系统不是单纯指电网而是电力系统与天然气系统的耦合系统英文文献里常写成 electricity-gas interconnected system。常见的耦合设备就两类燃气轮机gas turbine把天然气转成电电转气设备P2Gpower-to-gas把电变成氢或天然气。前者在电网里当灵活调节电源后者在新能源大发时把多余电量转成气相当于给电网加了一个储能池。为什么这个耦合在碳中和背景下越来越重要因为风电光伏出力不受调度端完全控制系统需要大量灵活性资源来平衡波动。燃气机组启动快、调节范围宽是天然的好帮手但天然气也排碳如果直接当成零成本清洁电用碳排放总量会很难看。P2G可以消纳弃风弃光再把产出的气送回气网形成闭环但P2G本身耗电而且能量转换效率只有六成左右怎么用、用多少必须靠优化模型算不能拍脑袋。把两个网络放一起优化跟分开优化最大的区别在于电网缺电时可以从气网调气发电气网气量不足时可以反过来让电网多出力P2G还能反过来补气。这种双向耦合如果不在同一个模型里算很容易出现电网说我有电气网说我缺气的尴尬局面。1.2 有功-无功协同解决的是传统调度的盲区传统经济调度ED只看有功平衡把各机组有功出力求出来之后再交给无功优化模块去调电压和无功补偿两步分开做。这种串行方案在系统裕度大、新能源占比低的时候问题不大但现在是高比例新能源接入问题就来了。举一个我实际遇到过的场景某个风电场集中接入的节点离主网电气距离远有功出力一大节点电压就往下掉需要大量无功支撑。但传统有功调度把风电场出力已经定死了无功优化模块这时候就算把附近的电容器全投上、机组无功顶到上限电压还是压不住只能被迫弃风。其实如果初始就把这个节点的有功出力降一点、把另一个节点的无功补偿配合好是可以既保住电压又少弃风的。这就是协同的价值——有功和无功本来就是一个问题的两面电压约束会限制有功传输能力有功潮流又会改变无功损耗分布拆开算必然顾此失彼。另一方面无功优化本身也需要有功配合。减少网损可以通过调整无功补偿实现但机组无功出力变化会引起有功损耗变化这是一条耦合路径反过来调节变压器分接头会改变无功潮流分布也会改变有功网损。协同模型一次性把这些耦合关系放进一个优化问题里得到的解才是真正全局最优的。2. 模型构建思路与关键原理2.1 目标函数经济性和低碳性怎么揉在一起这个模型的目标函数我最终写成四项相加全部以运行成本形式表达量纲统一方便求解器处理。第一项是火电机组发电成本用二次函数近似[ C_{G}\sum_{i \in G} (a_i P_i^2 b_i P_i c_i) ](a_i、b_i、c_i) 是成本系数(P_i) 是有功出力。无功出力对机组成本的影响相对小模型里没有直接进成本项但通过网损间接体现。第二项是天然气供应成本也就是气源产气费用。气网和电网一样有源节点每个气源有单位成本总和就是购气成本。第三项是碳排放成本。碳中和目标下不能只算经济账我引入的是碳交易机制给系统一个碳排放配额实际排放超过配额的部分按碳价购买低于配额可以出售。为简化我按实际排放量乘一个碳价系数 (\lambda_{CO_2}) 加入目标函数(\lambda) 取不同值就能画灵敏度曲线观察碳价对调度结果的影响。排放来源主要是火电燃烧和燃气轮机耗气折算出二氧化碳排放量。第四项是弃风弃光惩罚。用可再生能源的理论可发功率减去实际出力差值乘以一个较高的惩罚系数。这一步非常关键否则优化求解器为了省钱会天然的弃风弃光因为新能源边际成本是零但也要有出力成本不惩罚它就全弃了模型结果也就没意义了。这样四项加起来就是一个单目标优化问题避免了多目标加权的权重选取麻烦也更容易解释结果。如果想做成真正的多目标Pareto前沿可以用NSGA-II之类但作为Matlab工程实现加权单目标更适合配合商业求解器高效求解。2.2 电力网络约束潮流方程与二阶锥处理电力侧约束是模型里最重的部分。我用的是支路潮流模型Branch Flow ModelBFM也叫DistFlow特别适合配网和辐射型网络在一般输电网也能用。对每条支路 ((j,k))定义 (P_{jk})、(Q_{jk}) 为首端有功、无功功率(V_j) 为节点电压幅值约束如下节点功率平衡 [ \sum_{k:(j,k)} P_{jk} - \sum_{i:(i,j)} (P_{ij} - I_{ij}^2 r_{ij}) P_j^{gen} - P_j^{load} ] 无功同理把有功换成无功、电阻换成电抗。这个式子的物理含义很直白每个节点的净注入功率等于流入该节点的所有支路功率之和。支路潮流本身是非线性的因为 (P^2Q^2) 和 (V^2 I^2) 的关系。工程上最常见的处理是二阶锥松弛SOCP relaxation把非凸的等式松弛成不等式[ P_{jk}^2 Q_{jk}^2 \le S_{jk}^{2} ]实际上更标准的写法是引入辅助变量 (l_{jk}I_{jk}^2)、(v_jV_j^2)把约束转成二阶锥形式[ \left| \begin{matrix} 2P_{jk} \ 2Q_{jk} \ l_{jk}-v_j \end{matrix} \right|2 \le l{jk}v_j ]在YALMIP里这个约束一行就能写出来但理解它为什么是松弛很重要。严格说原问题里等号要求功率和电流的关系完全吻合松弛成不等式后相当于放大了可行域如果最优解上等式恰好取得称松弛紧那这个解就是原问题的真实最优解。对辐射型网络在满足一些条件下松弛是紧的这已经有理论保证。实际工程中我习惯算完后检查每条支路的松弛间隙万一有大于1e-4的支路就意味着这个解在物理上可能不存在需要处理。电压约束相对简单 [ \underline{V}_j \le V_j \le \overline{V}_j ] 通常取0.95到1.05 pu。机组出力约束、爬坡约束都是有上下界的线性约束这类约束直接作为sdpvar变量的边界条件加入即可。2.3 天然气网络约束Weymouth方程的线性化天然气网络稳态模型比电网更糙因为天然气的流动规律不是简单的基尔霍夫定律而是受管道压差驱动。对每条管道 (m-n)流量 (q_{mn}) 与两端压力平方差之间的关系服从Weymouth方程[ q_{mn}^2 K_{mn}^2 (p_m^2 - p_n^2) ]其中 (K_{mn}) 是管道常数与管径、长度、温度、摩擦系数有关。这个方程是非线性的而且涉及压力的平方处理起来比潮流方程更棘手。学术界处理它的办法大致有三类。第一类是直接保留非线性用IPOPT这类非线性求解器硬解小系统没问题系统稍微大一点就容易陷入局部最优或收敛困难。第二类是分段线性化把 (p_m^2 - p_n^2) 作为变量将 (q_{mn}^2) 用若干段线性函数逼近转化为混合整数线性规划MILP精度可控但要引入大量0-1变量。第三类是二阶锥松弛跟电网侧的思想一样先做变量替换再把等式松弛成锥约束。我实际采用的是分段线性化少量整数变量的方案原因是气网规模不大管道数在十条左右MILP规模完全能接受且求解稳定不会像CONOPT那样动不动就报奇点。天然气节点平衡约束和气源约束都是线性的[ \sum_{m} q_{mn} q_n^{source} - q_n^{load} - q_n^{P2G} q_n^{GT} 0 ]燃气轮机的耗气量与其有功出力有关我用了二次近似 [ q_n^{GT} \alpha_n P_n^{GT} \beta_n ] 这个线性表达式抓住主要特征即可太精细的耗气曲线反而会让模型变得不稳定。2.4 耦合设备与碳成本怎么进模型电-气耦合设备是模型里最出彩的部分两种设备方向相反燃气轮机消耗天然气输出电功率增加碳排放。在模型中它的耗气量进入气网节点平衡约束有功出力进入电网节点平衡约束碳排放按燃料热值换算后进入目标函数。P2G消耗电功率输出天然气本身不排碳但消耗的电在电网侧可能是火电发出来的其间接排放已经体现在电网的购电成本里不需要重复计算。P2G的效率我取了0.6也就是1 MW电进去0.6 MW气出来实际工程中电解水制氢再甲烷化的全链条效率也就在这个量级。碳排放量怎么折算我在代码里做一个简单的换算模块火电机组按天然气占比和煤耗估算单位发电量排放因子气电按单位耗气量的排放因子算最后统一转成吨二氧化碳乘以碳价加入目标函数。这样碳价从20元/吨扫到200元/吨就能清楚看到调度策略如何从煤电硬扛转向气电新能源P2G消纳。3. Matlab代码实现与实操细节3.1 环境配置与数据准备这套代码我跑在Matlab R2022b上建模用YALMIP工具箱求解外部调用Cplex。如果手头没有Cplex学术许可也可以用Gurobi或者退一步用开源的SCS求解二阶锥问题但整数变量一大SCS就会很吃力所以最好还是上商业求解器。数据方面电网部分我用的IEEE 30节点系统但要先做改造把部分火电节点改成燃气轮机节点在特定节点加风电和光伏再在合适的节点配置P2G设备与气网相连。气网部分没有标准算例我参考了经典的电-气互联6节点气网数据管道常数、气源上限都做了标幺化处理。这里必须提一个容易翻车的点单位制。电力系统习惯标幺值基准功率100 MVA、基准电压按各电压等级取但天然气网络的流量单位通常是kcf/h或者m³/h压力是bar直接混用会导致系数尺度差十几个数量级求解器直接罢工。我的做法是把气网也做成标幺制取一个基准流量和一个基准压力平方让K系数在0.1到10之间这样目标函数里各项成本的量级也能对上。3.2 模型核心代码变量定义与约束组装用YALMIP建模的关键就是把模型从数学公式翻译成约束和变量。下面这个代码块是模型的主骨架我尽量贴了实际能跑的写法%% 初始化 mpc loadcase(case30); % IEEE 30节点数据 baseMVA 100; % 基准功率 nB size(mpc.bus, 1); % 节点数 nG size(mpc.gen, 1); % 机组数 %% 定义决策变量 P_g sdpvar(nG, 1); % 机组有功出力 (pu) Q_g sdpvar(nG, 1); % 机组无功出力 (pu) V sdpvar(nB, 1); % 节点电压幅值 (pu) P_w sdpvar(nW, 1); % 风电场实际出力 Q_w sdpvar(nW, 1); % 风电场无功出力 P_p2g sdpvar(nP2G, 1); % P2G消耗有功 Q_svc sdpvar(nSVC, 1); % SVC无功补偿量 q_src sdpvar(nGasSrc, 1); % 气源产气量 %% 目标函数 fuelCost sum(a .* P_g.^2 b .* P_g c); gasCost sum(gasPrice .* q_src); carbonCost lambda_carbon * (sum(emit_g .* P_g) sum(emit_gt .* P_gt)); curtailPenalty M_pen * sum(P_w_max - P_w); objective fuelCost gasCost carbonCost curtailPenalty; %% 约束 C []; % 节点有功平衡 C [C, A_P * [P_g; P_w; -P_load] 0]; % 节点无功平衡 C [C, A_Q * [Q_g; Q_w; Q_svc; -Q_load] 0]; % 电压上下限 C [C, 0.95 V 1.05]; % 机组出力上下限 C [C, P_g_min P_g P_g_max, Q_g_min Q_g Q_g_max]; % P2G约束 C [C, 0 P_p2g P_p2g_max]; C [C, q_p2g_out eta_p2g * P_p2g]; % 电转气线性关系 % 气网平衡 C [C, A_gas * [q_src; q_p2g_out; -q_gt] 0]; % Weymouth分段线性化约束 C [C, weymouth_linearized(q_flow, p_sq, K, seg_breakpoints)];需要注意上面代码里A_P、A_Q、A_gas是我根据网络拓扑生成的关联矩阵这是电力和天然气网络的接线表必须保证每个节点的注入和流出严格相等这也是最容易写错的地方。我每次都写一个小函数把mpc.branch转成这个矩阵并打印一下节点注入总和验证平衡。3.3 关键难点Weymouth分段线性化代码气网管道约束是我调试时间最长的地方。分段线性化思路不复杂把压力平方差 (\Delta p^2) 的范围分成N段对每段用不同的斜率近似 (q^2)但要保证分段函数连续。实现时我用增量模型incremental model引入连续分段变量和0-1选择变量function C weymouth_linearized(q, dp2, K, seg, nSeg) % q: 管道流量变量 % dp2: 管道两端压力平方差变量 % seg: 分段断点位置 % nSeg: 分段数 C []; for k 1:nSeg % 每段的斜率和截距 slope(k) (Q_cal(seg(k1)) - Q_cal(seg(k))) / (seg(k1) - seg(k)); intercept(k) Q_cal(seg(k)) - slope(k) * seg(k); end % 引入0-1变量和连续分段权重变量 % 注意这里要求 dp2 的取值单调落在某个分段区间 % 每个管道的端点平方压力差要与网络其他节点耦合 C [C, q sum(lambda_w .* (slope .* dp2 intercept))]; C [C, sum(lambda_w) 1]; C [C, sum(alpha) 1]; end实际上完整的增量线性化要引入 SOS2 约束或0-1逻辑YALMIP里可以用sos2直接声明。但那需要单独的处理技巧代码容易膨胀。我的建议是如果气网管道少于15条直接用YALMIP的implies逻辑或者binvar来做分支选择都行重点是别一开始就追求理论上最优的建模方式先让模型跑通再逐步优化。3.4 求解设置与结果输出求解设置建议开verbose2这样能实时看到求解器迭代过程对定位问题非常有帮助options sdpsettings(solver,cplex, verbose, 2, cplex.mip.tolerances.mipgap, 0.001); diagnostics optimize(C, objective, options); if diagnostics.problem 0 % 提取结果 P_g_opt value(P_g); V_opt value(V); fprintf(目标值: %.2f\n, value(objective)); % 输出到图表 bar(P_g_opt); end结果分析环节我最看重三个东西一是各节点电压断面看有没有节点被压在边界上二是发电机无功出力分布判断无功源是否配置合理三是碳排放总量随碳价的变化曲线这个能直接说明碳中和目标如何引导调度结果。4. 常见问题与调试经验实录4.1 求解器报错排查速查跑这套模型半年多我把最常见的报错整理成一张表照着查能省大量时间报错现象常见原因解决办法No suitable solverYALMIP没识别到Cplex/Gurobi重新运行yalmiptest检查求解器路径是否加入MatlabInfeasible problem约束互相矛盾比如功率平衡要求负荷等于出力但机组上限太低先用只含线性约束的简化模型验证网络数据再加非线性约束Numeric issues变量量级差异太大尤其气网参数没标幺化检查所有K系数、流量、压力是否在同一标幺基下mipgap一直不收敛整数变量过多或分段数太多减少分段数或把一些整数变量改成正实数变量Memory limit exceededMILP规模扩大增大Cplex内存限制或减小测试系统规模先验证方法4.2 模型不可行的定位技巧模型报不可行是最让人头疼的尤其是约束一多根本不知道是哪一个约束在打架。我的定位办法是二分法松弛变量。先给所有等式约束加一个很小的正松弛变量 (s)并给目标函数加上对 (s) 的大惩罚。如果加了松弛之后问题可解那value(s)大于0的地方就是矛盾源。这个方法虽然土但非常管用。我就在调试时发现气网和电网的基准功率不一致导致P2G转换等式系数差了一个数量级模型怎么都不可行加了松弛变量一查就定位到了。另一种情况是电压约束0.95到1.05之间本身不可行。电网数据在基准情况下正常但加上P2G大功率消耗后节点电压直接被拉低。这种情况可以先放宽电压范围到0.9到1.1确认其他约束没问题再逐步收缩范围观察是哪个节点在制造瓶颈。4.3 二阶锥松弛不紧怎么办运行SOCP模型后我习惯检查每条支路的松弛间隙。画出来如果有个别点离0很远说明那个支路的松弛是松的解在物理上不严格成立。常见补救手段有两个一是在目标函数里加一项对松弛量的微小惩罚逼迫解向等式靠拢代价是目标结果轻微失真二是改用精确的潮流方程把问题交给非线性求解器做局部优化。我实测的经验是辐射型为主的系统SOCP松弛几乎都是紧的但一旦出现环网或者系统重负荷接近极限松弛容易出现间隙。对于我这个电-气互联系统因为气网侧做了分段线性化本身就不是完全非线性精确解所以轻微的松弛间隙反而没必要过于纠结——整个模型本来就是近似框架下的最优解一致性比数学上的严格紧更重要。4.4 参数标幺化最容易被忽视的坑最后专门说说标幺化。很多第一次做电-气耦合模型的同学电网用100 MVA标幺气网直接用工程单位两边变量一相加Cplex直接提示数值问题。我自己的处理方式是统一到一个基准下气网的流量基准取每小时10,000 m³压力平方基准取60 bar²这样管道K系数就变成0.1到2之间的小数和目标函数里电网成本项能对上。另外电压和无功的基准也要统一。无功功率在IEEE数据里通常以Mvar为单位在YALMIP建模时如果用标幺值要注意节点无功负荷、无功补偿上限都要除以基准功率。这些小细节错一个结果就会出现无功缺额巨大但显示正常的诡异现象。5. 我实际做完这套模型之后的体会这套模型跑通之后最直观的感受是碳中和目标真的能改变调度结果。碳价低的时候系统倾向于用便宜的火电硬扛弃风率也高碳价抬到一定水平燃气机组和P2G开始进场弃风率下降但购气成本上升两条曲线有一个明确的交叉点那个交叉点就是既要低碳又要经济的平衡点。这就是为什么我觉得做优化的人不能只盯目标函数值灵敏度分析才是模型真正有价值的部分。最后分享一个小建议如果你也是刚开始搭这类模型千万不要一上来就试图把电网、气网、碳交易、新能源、无功补偿全部塞进一个大模型里。先拿一个小系统比如9节点电网加3节点气网把电-气耦合跑通再加入无功和碳成本最后再换大系统。我就是在IEEE 30节点上直接改结果每次报错都要从头排查效率低得令人崩溃。而当我退回去先用小系统验证建模逻辑后再换大系统其实非常顺利一小时就能跑完所有算例。这套思路也推荐给你。