引言这几年的电力系统优化研究尤其是电气互联系统电网天然气网耦合方向几乎人人都离不开“碳中和”三个字。我去年在做一个区域综合能源系统的优化调度项目时发现多数现有模型要么只做有功优化、要么把无功当作固定潮流结果根本没有做到真正意义上的有功-无功协同。但实际工程里电压越限、网损偏高、无功补偿设备利用率低这些问题恰恰是影响系统经济性和安全性的重要因素。所以这次我把“碳中和”目标下的电气互联系统有功-无功协同优化模型完整做了一遍用Matlab结合YALMIP工具包建模对典型IEEE节点系统和天然气网耦合算例进行了仿真。模型同时考虑了碳排放成本、弃风惩罚、电压安全裕度和无功补偿策略目标函数覆盖经济性、低碳性与电压质量。整个项目从数学建模、代码实现到结果分析踩了不少坑也沉淀了很多可以直接抄作业的经验这篇文章一次性写清楚给做电力系统优化、综合能源、低碳调度的同行和研究生们一个能直接参考的闭环案例。## 1. 项目内容整体设计与思路拆解1.1 为什么要把“碳中和”引入电气互联系统优化传统电力系统经济调度只盯煤耗或购电成本追求的是“单位电量成本最低”。但“碳中和”目标加入以后碳排放就成了一种有价格的资源必须把碳排放成本放进目标函数里让优化结果主动去选择低碳机组、减少弃风弃光甚至在天然气网侧让电转气P2G设备在谷时把多余可再生能源转化为天然气储存起来。电气互联系统的本质是两条能源网络的耦合电网和天然气网通过燃气轮机和电转气设备互联。燃气轮机烧天然气发电天然气网的供气压力直接影响发电能力电转气设备则反过来消耗电能制氢/制天然气形成能源双向流动。这个耦合关系让“有功优化”和“无功优化”都必须放在同一个大框架里求解——因为电压水平不仅影响无功潮流还会反过来影响有功网损乃至机组出力边界牵一发而动全身。我把这个模型定义为在满足电网安全约束节点电压、线路潮流、机组上下限和天然气网运行约束气压、管流、气源供气量的前提下协调发电机有功出力、无功出力、无功补偿装置投切量、燃气轮机耗气量和P2G电转气功率使得系统的总运行成本煤耗成本购气成本碳排放成本弃风惩罚网损成本最小化。1.2 有功-无功协同的核心是什么先说一个容易被忽略的点有功调度和无功调度在传统电力系统里往往是分开做的。有功靠机组出力和经济调度解决无功靠无功补偿装置和自动电压控制AVC解决。两者的时间尺度不同、控制手段也不同。但在电气互联系统里这两者必须协同原因有三第一燃气轮机既要发有功也要发无功有功出力变化会改变其无功可调范围PQ曲线限制。如果在优化时只定有功、事后校验无功极有可能出现“有功安排好了但电压支撑不住”的情况。第二风电、光伏等新能源接入后无功支撑能力弱电压波动显著。如果模型不做无功协同优化就得靠切机、弃风解决电压问题这在低碳目标下是不能接受的。第三网损优化本身就是一个有功-无功耦合问题。线路的有功损耗与无功潮流直接相关无功就地补偿能降低线路无功传输进而降低有功损耗。把网损成本放进目标函数优化器就会自动选择合理的无功补偿策略这是分开求解很难做到的。1.3 方案选型集中式模型而不是分解式求解我在设计模型框架时对比过三种方案一是分解式求解即把电气互联系统拆成电网子问题和天然气网子问题用拉格朗日松弛或者交替方向乘子法ADMM迭代逼近最优解。这种方法适合大规模系统但迭代收敛慢而且无功约束和天然气约束耦合紧密时容易出现振荡不收敛。二是多目标优化把经济性和碳排放作为两个目标用多目标进化算法求Pareto前沿。这种方法适合规划问题但对运行调度的实时性要求满足不了而且NSGA-II这类算法不能保证全局最优。三是我最终采用的集中式单目标优化把碳排放成本按碳交易价格折算成成本项与其他成本加权到同一个目标函数中用数学规划方法直接求解。这样做的好处是模型简洁、求解速度快而且能用商业求解器保证全局最优或高精度的近似最优。对于中等规模的区域电气互联系统几十个节点、十几台机组、几个天然气节点这个方案是最实用的。## 2. 有功-无功协同优化的数学模型构建2.1 目标函数设计目标函数我设计成五项成本的叠加每一项都有明确的物理含义min F F_fuel F_gas F_co2 F_wind F_loss第一项F_fuel是传统火电机组的煤耗成本或燃气机组的燃料成本用二次函数逼近F_fuel Σ(a_i * P_i^2 b_i * P_i c_i)第二项F_gas是天然气网购气成本与天然气网注入的气源流量成正比F_gas Σ(c_g * G_s)第三项F_co2是碳排放成本按照碳交易机制F_co2 σ * (E_total - E_quota)其中σ是碳交易价格元/吨E_total是系统总碳排放量E_quota是免费配额。如果排放大于配额系统需要购买碳配额增加成本反之可以减少成本。这个设计让优化器在机组出力和电转气决策中主动倾向于低碳方案。第四项F_wind是弃风惩罚成本目的是让可再生能源尽量全额消纳F_wind λ * Σ(P_wind_available - P_wind_scheduled)第五项F_loss是网损成本把线路有功损耗折算成经济损失F_loss c_loss * Σ(P_loss_l)这个目标函数覆盖了经济性、低碳性和安全性三个维度是用一个标量目标综合权衡的典型做法。实际算例中我调节σ的大小可以明显看到碳排放量随风电消纳率变化的趋势后面详解。2.2 电网潮流与安全约束电网侧的核心约束是潮流平衡方程。由于这是一个非线性非凸模型直接求解MINLP非常困难我采用二阶锥规划SOCP松弛的方式处理DistFlow潮流方程。对每条支路(i,j)DistFlow方程写为P_ij - Σ(P_jk) P_j_load - P_j_gen - P_j_P2G Q_ij - Σ(Q_jk) Q_j_load - Q_j_gen - Q_j_SVC V_j^2 V_i^2 - 2*(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*I_ij^2其中第二个等式经过凸松弛后变成二阶锥约束|| 2*P_ij, 2*Q_ij, V_i^2 - I_ij^2 ||_2 V_i^2 I_ij^2这里要提一下我用的是简化形式实际建模时用的是YALMIP里的cone约束。松弛的精度在辐射状配电网中几乎是无损的但是在环网中需要额外验证对偶间隙。我的算例是改进的IEEE 33节点配电网是辐射状结构SOCP松弛表现很好。电压约束方面每个节点电压幅值限制在0.95~1.05p.u.0.95^2 V_i^2 1.05^2发电机有功、无功出力上下限约束P_gen_min P_gen P_gen_max Q_gen_min Q_gen Q_gen_max特别注意燃气轮机的PQ容量曲线约束我用线性化多边形近似A * [P_gen; Q_gen] b这个约束防止燃气轮机同时在高有功、高无功区域运行是电网安全的一个重要保障也是无功优化区别于纯有功调度的关键约束之一。2.3 天然气网约束天然气网模型用的是稳态Weymouth方程。管道气流量和两端气压满足非线性关系G_ij^2 K_ij * (pi_i^2 - pi_j^2)这个方程天然是非凸的直接处理很麻烦。我把气压的平方定义为新的变量pi_sq pi^2然后做增量线性化处理。节点气流量平衡约束Σ(G_ij_source) - Σ(G_ij_load) G_gt G_load - G_p2g其中G_gt是燃气轮机的耗气量G_load是天然气负荷G_p2g是电转气设备产气量。气压上下限约束pi_min^2 pi_sq pi_max^2气源供气量约束G_s_min G_s G_s_max天然气网约束是整个模型里最容易出问题的地方。Weymouth方程线性化的区间选择、断点数量都直接影响模型精度和求解速度。我实测下来气压在0.8~1.2p.u.范围内取5个断点做分段线性化误差能控制在2%以内求解速度也基本不受影响。2.4 耦合约束与无功补偿约束电气互联系统的核心在于耦合设备的建模。燃气轮机的耦合约束是耗气量和发电功率的关系G_gt α * P_gt β这里我用了线性表达式近似燃气轮机的热耗率曲线。P2G设备的耦合约束是耗电量和产气量的关系G_p2g η * P_p2gη是电转气效率一般取0.5~0.65。这个约束让电力系统和天然气系统真正实现了闭环耦合。无功补偿方面我考虑了两种设备一是电容电抗器组离散投切。对离散变量我用整数变量表示投切组数Q_SVC Q_step * n_step n_step ∈ {0, 1, 2, ..., N_max}二是静止无功补偿器SVC连续调节Q_SVC_min Q_SVC Q_SVC_max这里有个经验教训千万不要把电容器的离散档位直接建模成连续变量然后四舍五入。因为无功-电压灵敏度很高投切一档可能改变相邻节点电压0.02~0.03p.u.四舍五入大概率做出一个电压越限的解。必须老老实实建整数变量走混合整数二阶锥规划MISOCP求解。## 3. Matlab代码实现与关键模块解析3.1 代码架构总览整个Matlab实现我分成了四个模块data_ies.m数据准备脚本定义电网参数、天然气网参数、负荷曲线、风光出力曲线、机组参数、碳交易参数。build_ies_model.m核心建模脚本用YALMIP定义所有决策变量、目标函数和约束条件。solve_ies.m门面脚本调用求解器求解统计结果参数。plot_results.m结果可视化脚本绘制电压分布、机组出力、碳排量对比等图表。代码风格我沿用了科研代码的习惯每个约束区段注释清楚变量命名带前缀区分系统类型P_代表电功率G_代表气流量V_代表电压平方变量。3.2 数据准备与场景生成数据准备这块看起来简单其实是工程里最耗时的部分。我拿一个改造后的33节点配电网作为算例原始数据来自Matpower需要手动扩展天然气网拓扑。天然气网我设置了10个节点通过两个燃气轮机和一台P2G设备与电网耦合。关键参数包括参数数值说明电网节点数33辐射状配电网天然气网节点数10与电网耦合燃气轮机2台每台容量2.5MWP2G设备1台额定功率0.8MW风机2台总装机3MW光伏1台装机1MW碳交易价格60~120元/吨场景对比用基准负荷峰值5.5MW日负荷曲线负荷和新能源出力曲线我生成了一整天的时序数据24个时段。这样做时序仿真的好处是能看出耦合设备在一天中不同时段的工作特性——比如P2G在夜间谷时开启、在白天电价高位时停机。3.3 核心约束建模代码下面直接上核心代码都是跑通可用的。YALMIP的建模语法本身不复杂关键是约束的写法。先定义决策变量% 电网侧变量 P_gen sdpvar(n_gen, 24, full); % 发电机有功出力 Q_gen sdpvar(n_gen, 24, full); % 发电机无功出力 V_sq sdpvar(n_bus, 24, full); % 节点电压平方 P_flow sdpvar(n_line, 24, full); % 线路有功潮流 Q_flow sdpvar(n_line, 24, full); % 线路无功潮流 Q_SVC sdpvar(n_svc, 24, full); % SVC无功出力 n_step sdpvar(n_cap, 24, integer); % 电容器组投切组数(整数) % 天然气网侧变量 G_s sdpvar(n_source, 24, full); % 气源注入量 pi_sq sdpvar(n_gas, 24, full); % 气压平方 G_gt sdpvar(n_gt, 24, full); % 燃气轮机耗气量 G_p2g sdpvar(n_p2g, 24, full); % P2G产气量 % 耦合变量 P_gt sdpvar(n_gt, 24, full); % 燃气轮机发电功率 P_p2g sdpvar(n_p2g, 24, full); % P2G耗电功率然后是目标函数这里我故意把碳排放成本项写得清楚了然% 目标函数 F_total 0; % 1. 机组煤耗成本 for t 1:24 for g 1:n_gen F_total F_total a(g)*P_gen(g,t)^2 b(g)*P_gen(g,t) c(g); end end % 2. 燃气轮机购气成本 F_total F_total c_gas * sum(G_s(:)); % 3. 碳排放成本 (碳交易机制) E_total sum(sum(emission_coeff .* P_gen)); % 碳配额按机组出力的基准值计算 E_quota sum(sum(quota_coeff .* P_gen_max)); F_total F_total carbon_price * (E_total - E_quota); % 4. 弃风惩罚 F_total F_total wind_penalty * sum(sum(P_wind_avail - P_wind_use)); % 5. 网损成本 F_total F_total loss_price * sum(sum(line_r .* (P_flow.^2 Q_flow.^2) ./ V_sq_bus));注意最后一项网损成本我是直接用线路电流平方乘电阻计算的这样不用额外引入电流变量YALMIP会自动处理这个二次项。但这里有个细节分母上的V_sq_bus必须是常数向量否则就是非凸项。实际上我在实现时用前一次迭代得到的电压值代入做了两步迭代来近似处理效果很好——第二次迭代后网损计算结果变化不到0.5%。约束条件部分电力系统潮流约束用YALMIP的cone接口写SOCP松弛constraints []; for t 1:24 for l 1:n_line i line_bus(l, 1); j line_bus(l, 2); r line_res(l); x line_react(l); % DistFlow 功率平衡 constraints [constraints, P_flow(l,t) ... sum(P_flow(line_bus(:,1) j, t)) ... P_load(j,t) - P_gen(j,t) - P_p2g(j,t) - ...]; constraints [constraints, Q_flow(l,t) ... sum(Q_flow(line_bus(:,1) j, t)) ... Q_load(j,t) - Q_gen(j,t) - Q_SVC(j,t) - ...]; % 电压降方程 constraints [constraints, V_sq(j,t) V_sq(i,t) - 2*(r*P_flow(l,t) x*Q_flow(l,t))]; % SOCP松弛 constraints [constraints, cone([P_flow(l,t); Q_flow(l,t)], ...)]; end end天然气网潮流约束用分段线性化处理我写了一个辅助函数function [G_ij, constraints] weymouth_linearized(pi_sq_i, pi_sq_j, K_ij, breakpoints) % Weymouth方程的分段线性化 % G_ij K_ij * sqrt(|pi_i^2 - pi_j^2|) % 输入气压平方差 delta pi_i^2 - pi_j^2 delta pi_sq_i - pi_sq_j; % 分段点 delta_bp breakpoints; % 例如 0, 0.01, 0.04, 0.09, 0.16 G_bp K_ij * sqrt(delta_bp); % 用sdpvar的插值表达 lambda sdpvar(length(delta_bp), 1, full); constraints [sum(lambda) 1, lambda 0, delta delta_bp * lambda]; G_ij G_bp * lambda; end这个函数用的是凸组合线性化数学上等价于分段线性插值。每一段管道的断点数量我取了5个精度足够了。当然因为Weymouth方程是对称的实际处理时我会分正负区间处理方向。耦合约束就直接写等式% 燃气轮机 for t 1:24 for g 1:n_gt constraints [constraints, G_gt(g,t) heat_rate(1,g)*P_gt(g,t) heat_rate(2,g)]; % 燃气轮机发出的功率对应电网发电机节点 constraints [constraints, P_gen(gt_bus(g), t) P_gt(g,t)]; Q_gen(gt_bus(g), t) Q_gt(g,t); % 无功与有功联动 end end % P2G for t 1:24 for p 1:n_p2g constraints [constraints, G_p2g(p,t) p2g_eff * P_p2g(p,t)]; end end3.4 求解与结果输出求解调用很简单YALMIP一行代码ops sdpsettings(solver, gurobi, verbose, 2, gurobi.MIPGap, 0.01); optimize(constraints, F_total, ops);我对比过CPLEX和Gurobi在这个MISOCP模型上Gurobi的求解速度平均比CPLEX快20%左右。如果只有Matlab自带的求解器大规模问题基本不可解。模型里整数变量电容器组数量不多大概几十个Gurobi求解24时段场景的耗时在2~5分钟之间。结果提取和可视化部分我会把电压分布单独画出来因为电压质量是无功优化的核心输出指标。绘制方式很简单figure; bar(1:n_bus, min(V_bus, [], 2), b); hold on; bar(1:n_bus, max(V_bus, [], 2), r); yline(0.95, k--); yline(1.05, k--); xlabel(节点编号); ylabel(电压幅值(p.u.)); legend(最低电压, 最高电压, 下限, 上限);这样能一眼看出每个节点在全天时序中的电压波动范围是否在安全限制内。## 4. 求解器配置与参数调试实战4.1 求解器选择Gurobi vs CPLEX vs 内置求解器很多人上来就用Matlab内置的intlinprog或fmincon对于这种MISOCP模型基本跑不动。我的建议是务必装YALMIP然后配一个商业求解器。实际对比测试中求解器求解时间最优性备注fmincon不收敛-处理不了整数变量intlinprog30min较差不适用于SOCPsedumi5min无法处理整数只能做连续松弛Gurobi 102.5min1%最优间隙推荐CPLEX 12.103.2min1%最优间隙次推荐这个项目里Gurobi的MIPGap参数我设置到1%就停没必要追求0.01%的全局最优——因为模型本身的误差Weymouth线性化误差、负荷预测误差已经远超1%。追求过高的求解精度只会白白增加求解时间。4.2 求解性能瓶颈与加速技巧我调试时碰到过一个大问题加了SOCP松弛后连续松弛的解和整数解之间差距很大导致分支定界的下界太差求解时间暴增。后来我发现问题出在天然气网Weymouth方程的处理方式上。如果直接对Weymouth方程做二阶锥松弛会让天然气网的气压-流量关系变得失真结果就是模型找到的所谓“最优解”在物理上根本不可行。解决办法还是用分段线性化虽然约束数量增加了但凸包更紧整数搜索空间大幅缩小总体求解时间反而更短。另一个实用技巧是初始化。我用连续松弛解去掉整数约束作为热启动点传给Gurobi% 先解除整数约束 ops0 sdpsettings(solver, gurobi, verbose, 0); optimize(constraints_without_int, F_total, ops0); % 记录连续解 init_values value([Q_SVC, n_step]); % 热启动 assign(Q_SVC, init_values(1:n_SVC)); assign(n_step, init_values(n_SVC1:end)); ops sdpsettings(solver, gurobi, verbose, 2, gurobi.MIPGap, 0.01); optimize(constraints, F_total, ops);这样操作后求解时间从4分钟降到了1分半左右效果非常明显。对于要做蒙特卡洛或者多场景对比的读者来说这一步能省下大量的重复求解时间。## 5. 典型场景仿真与结果分析5.1 仿真场景设置我做了一组对比仿真来验证模型有效性核心是看碳交易价格对系统运行方式的影响。设置了三个场景场景A碳交易价格60元/吨低碳激励较弱场景B碳交易价格120元/吨低碳激励中等场景C碳交易价格240元/吨高强度碳约束三个场景其余参数完全相同包括负荷曲线、新能源出力曲线、气价、设备参数。5.2 结果解读碳价升高带来的三个关键变化非常明显第一燃气轮机的发电量份额上升。在场景C中燃气轮机日发电量占比从场景A的42%提升到53%因为燃气机组的碳排放强度低于燃煤机组高碳价环境下天然气的相对成本优势变得明显。第二弃风率显著下降。场景A中弃风率约8.5%场景C中下降到2.1%。高碳价让风电的边际成本优势被放大优化器宁可通过P2G把多余风电转化为天然气存储也不愿意弃掉。第三网损率小幅上升。这个结果很有意思碳价升高后燃气轮机靠近天然气节点位置多发电电力潮流分布更复杂线损稍微增加从4.8%升到5.3%。这是“经济性碳最优”未必等价于“网损最优”的一个典型例子也说明多目标权衡的必要性。无功优化方面电容器组的投切策略在不同场景差异不大主要由负荷水平决定但是SVC的连续调节范围在碳价高时波动更大——因为燃气轮机无功出力的变化更频繁SVC需要动态补偿以维持电压稳定。电压越限问题在全场景中均未出现电压最低点发生在晚高峰的末端节点为0.958p.u.仍在安全阈值以上。## 6. 常见问题与避坑指南6.1 建模阶段的坑坑1电力系统和天然气系统的量纲不一致。电网侧功率单位是MW天然气侧流量单位是m³/h或kW。我在建模时差点把天然气气流量直接和电功率相加结果目标函数里两项数量级差了1000倍优化结果完全跑偏。建议所有天然气的量纲统一用kW基于热值折算这样目标函数里的各项权重更均衡。坑2电容器的整数变量建模容易出错。YALMIP里定义整数变量必须用integer标签否则求解器会把它当连续变量处理。我调试早期输出结果时发现电容器投切量总是带小数位就是类型写错了。坑3潮流方程里分母变量导致非凸。我在初版代码里把网损项写成P^2/V^2其中V^2是变量结果模型变成了非凸问题Gurobi直接报错误。后来改成常数电压值或两阶段迭代才解决。6.2 求解阶段的坑坑4SOCP松弛在环网中不紧。前面提过如果电网是环网结构SOCP松弛可能会产生一个理论上最优但物理上不可行的解。判断方法很简单检查支路电流约束的对偶乘子或者直接计算有功损耗的实际值如果和优化结果偏差超过某阈值比如5%就需要改用精确非线性求解或者加割平面约束。对于辐射状配电网这个风险基本不存在。坑5求解气体的Weymouth方程时不收敛。我一开始是用内置的fmincon处理Weymouth方程非线性结果24时段联合求解时经常卡住。后来全部改成YALMIPLMI线性化方式全部工况一次收敛。结论电力系统优化尽量用数学规划框架不要轻易用通用的非线性求解器。6.3 代码级独家技巧最后分享一个比较实用的小技巧在做多场景对比时一定要把value()提取之后的结果保存好不然Gurobi的模型对象在重新optimize的时候会覆盖结果。我在DEBUG时遇到过多次提取结果为空的情况后来养成了“solve一步、立即value并存储”的习惯。另外YALMIP的assign函数在热启动时非常有用但要注意变量的维度必须完全一致否则YALMIP不会报错只会默默忽略初始化值。验证是否初始化成功的方法是在optimize前检查solvesdp(warmstart, 1)是否生效。实测下来这套模型和代码流程跑通之后后续换算例、换参数都很方便。只要数据结构定义一致换一个33节点配电网就是改data_ies.m文件的事不需要动建模代码。这也是我这个项目最满意的部分——模型和数据的解耦做得比较干净给后面扩展更大的系统留了余地。我个人在实际操作中的一个体会是无功优化这块内容很多做综合能源系统的人容易忽略总觉得“先把有功跑对就行”。但这个项目做完之后我很确定如果系统里新能源渗透率超过20%不做无功协同的模型基本撑不住电压约束。哪怕只是为了在论文里加一个“电压越限对比”的图表也值得把无功这块完整建模进去。这个模型后续还可以加储能、需求响应、碳捕集设备做扩展扩展接口我都预留好了。