电气互联系统有功-无功协同优化：碳约束下的建模与求解

1. “碳中和”目标下的电气互联系统为什么非做有功-无功协同不可先说一个我在做这类项目时最深的体会很多做电力系统优化的同学把重点全放在有功调度上无功这块要么忽略要么用固定的功率因数折算一下。但在“碳中和”目标下电气互联系统的优化边界已经被彻底拓宽了——你面对的不再是单一电网而是电力系统与天然气系统深度耦合的网络负荷侧多了电转气、氢能等元素电源侧多了大量分布式可再生能源。此时如果还把有功和无功拆开做结果不是次优而是可能直接得出不可行的运行方案。所谓电气互联系统通俗讲就是电网和气网通过两类关键元件“绑”在了一起一类是燃气机组把天然气转化为电力另一类是电转气设备Power to GasP2G把富余电能转化为天然气。前者是天然气系统向电力系统供电后者是电力系统反过来向天然气系统供气。这种双向耦合的直接后果就是电网的有功潮流变化会改变气网的供气需求气网的供气压降又反过来约束燃气机组的出力上限而机组无功出力调整又会改变线路无功潮流和电压分布——环环相扣谁也没法独立优化。为什么“碳中和”目标让这个问题更棘手因为高比例可再生能源接入后系统惯量降低、电压支撑能力变弱无功资源变得更加稀缺。光伏和风电在满发时可能挤占常规机组的出清空间但它们的无功调节能力有限甚至某些工况下还要吸收无功。另一个问题是碳排放约束把气网侧的成本结构改变了——碳价、碳配额都进入了目标函数这就让“电转气、气转电”的决策路径变得更加复杂。我做这个项目的核心思路可以概括成一句话在满足碳排放约束的前提下同时调节机组有功出力、无功出力、变压器分接头、无功补偿装置投切以及气源供气量让整个电气互联系统的综合运行成本最低。这个思路拆解下来涉及三个核心问题如何建立电网和气网统一的数学模型特别是描述耦合元件的变量关系。如何处理碳排放约束和目标函数中碳成本的量化。如何在Matlab中高效求解这样一个大规模非线性混合整数优化问题。下面我会沿着这三条线把从建模到求解再到代码实现的全过程完整分享出来。这篇文章适合正在做综合能源系统优化、无功优化或者碳约束调度的研究生和工程师参考也适合想用Matlab搭建一体化优化框架的读者当模板。2. 有功-无功协同的耦合机理数学上它们是怎么“纠缠”在一起的很多人觉得“有功-无功协同”就是把两个目标函数加到一起然后在约束里同时写上机组有功上下限和无功上下限。这是典型的误区。真正的耦合体现在等式约束和不等式约束交叉关联上下面我把数学模型层面的耦合点逐一拆开。2.1 电网侧数学模型潮流方程中的有功无功耦合电网侧采用交流潮流模型时节点i的功率平衡方程是P_i U_i * Σ(U_j * (G_ij * cos(θ_ij) B_ij * sin(θ_ij)))Q_i U_i * Σ(U_j * (G_ij * sin(θ_ij) - B_ij * cos(θ_ij)))从这个式子就能看出有功P_i和无功Q_i共享同一组电压幅值U和相角θ变量。改变任意一个节点的无功全网的电压分布都会变电压变了相角也跟着变相角变了有功潮流就变了。这就是电网侧最本质的耦合。在工程简化中有人用直流潮流模型只保留有功维度这在输电网经济调度里勉强可用但一旦涉及电压安全约束、无功补偿优化、变压器分接头调节直流潮流就完全失效了。本模型必须用交流潮流。对于输电线路还有一个关键约束是线路潮流极限。实际传输容量受热稳定极限、电压约束和稳定极限三重限制用视在功率S_ij的平方表示就是S_ij^2 P_ij^2 Q_ij^2 ≤ (S_ij^max)^2这个约束本身就是有功与无功的非线性耦合约束。如果只优化有功而忽略无功线路可能会在视在功率上越限而单纯看有功却一切正常。2.2 天然气网络模型气体流动方程与节点气压约束天然气网络的稳态模型常用节点气流平衡方程描述。对气网节点n气流平衡可以写成Σ(FLOW_k) G_source,m - G_load,m - G_p2g,m 0其中FLOW_k表示管道k的气流量正方向为流入节点反方向为流出G_source是气源供气量G_load是气负荷包括燃气机组耗气量和民用气负荷G_p2g是P2G设备消耗的电量转化成的天然气注入量。管道气流FLOW与两端气压的关系通常用Weymouth稳态方程FLOW_ij^2 K_ij^2 * (p_i^2 - p_j^2)这个式子是非线性的而且是本模型的第一个非线性难点来源。注意这里FLOW是有方向的实际工程中需要处理方向符号问题常见做法是引入双向流标志位或者直接假设气流方向已知——在做日前优化调度时气流方向通常可以从潮流结果预判但严谨起见最好在约束中保留方向变量。气压还有上下限约束p_min ≤ p_i ≤ p_max气源出力也有爬坡限制G_source,min ≤ G_source ≤ G_source,max这里有个值得注意的细节气网的惯性比电网大得多但稳态模型下管道存气linepack效应被忽略了。在日内滚动优化中linepack往往能提供宝贵的短时缓冲能力但本文的日前模型先不考虑保持模型相对简洁。2.3 耦合元件模型燃气机组与P2G的双向耦合燃气机组是气网到电网的耦合元件其模型包括两部分注入电网的有功出力P_g和旋转备用容量。从气网消耗的天然气量F_gas a * P_g^2 b * P_g c其中a、b、c是燃料消耗系数。这个二次函数也是典型的非线性环节。另一个方向——P2G设备的耦合模型为G_p2g η_p2g * P_p2g / GHV其中P_p2g是P2G设备消耗的有功功率η_p2g是转换效率GHV是天然气高位热值。P2G消耗的功率上限同时受电网上网容量和气网注入能力的双重约束P_p2g,min ≤ P_p2g ≤ min(P_p2g,max, G_p2g,max * GHV / η_p2g)到这里有功-无功耦合的数学画像就清晰了电网侧P、Q通过电压幅值和相角在潮流方程中耦合。气网侧气流量通过压力非线性方程耦合。系统级燃气机组和P2G在有功功率与气流量之间建立双向通道而这种耦合通过上述两个非线性环节传递。抽象成一句话有功决策变了无功优化空间也跟着变无功决策变了电压分布变了P2G消耗的有功上限和燃气机组的可用容量也变了。所谓“协同”就是在数学上同时处理这两组变量而不是把它们解耦成两个子问题分别求解后再拼起来。3. 协同优化模型搭建目标函数与关键约束的设计思路我在这类项目上走过的弯路不少最值钱的一条经验是目标函数决定了优化方向约束决定可行域边界但真正决定模型能否被高效求解的是约束的“数学形态”是线性、凸还是非凸。3.1 目标函数设计成本与碳排放的双重权衡本模型的目标函数我在实际项目中是这样设计的综合考虑了四种成本min F_total Σ(C_fuel) Σ(C_co2) Σ(C_p2g) Σ(C_grid)第一项是常规机组燃料成本是机组有功出力的二次函数C_fuel,i a_i * P_i^2 b_i * P_i c_i第二项是燃气机组的碳排放成本这里引入“碳中和”目标的关键我在实际建模中用了两种处理方式碳市场交易价格乘以排放量或者设定碳配额上限然后用阶梯碳价。模型里我建议使用线性碳价对单位碳排放量乘以碳价系数C_co2,i ε_i * e_i * P_i其中ε_i是碳价元/tCO2e_i是单位电量的碳排放强度tCO2/MWh。第三项是P2G运行成本主要包括设备运维成本和购电成本分摊C_p2g γ * P_p2g其中γ是P2G单位运行成本系数。第四项是从上级电网购电的成本只有当系统接入主网时才需要。这里我想特别说明为什么要这样设目标函数而不是简单加一个“碳排放最小”目标。因为纯减排目标在经济上不可行实际运行中调度员关心的是“在碳约束下的最经济运行”所以碳成本进入目标函数、碳约束进入不等式约束这两个维度分开处理才能保证结果既环保又可执行。3.2 等式约束详解电网潮流与气网平衡电场部分需要同时满足有功和无功平衡P_G,i - P_L,i - U_i * Σ(U_j * (G_ij * cos θ_ij B_ij * sin θ_ij)) - P_P2G,i 0Q_G,i - Q_L,i - U_i * Σ(U_j * (G_ij * sin θ_ij - B_ij * cos θ_ij)) - Q_P2G,i 0注意我在P2G节点的有功平衡里专门扣掉了P_P2G,i这是典型的“负荷型元件”处理方式P2G消耗有功同时如果是电解水制氢再合成甲烷的路径还需要向电网提供少量的无功给整流器等设备。气网侧的气流平衡在2.2节已经给出。这里补充说明一个细节对于含压缩机站的天然气网络压缩机的耗气量也是变量但在基础模型中往往先用固定比例折算C_comp β * FLOW_comp这种简化在工程上是可接受的因为压缩机耗气量占总输气量的比例通常只有1%到3%不会对优化结果产生方向性影响。3.3 不等式约束安全与设备极限的建模技巧不等式约束是我在实际代码中掉过坑最多的地方。常见的约束包括发电机组有功出力上下限P_min ≤ P ≤ P_max。无功出力上下限Q_min ≤ Q ≤ Q_max这里要注意无功上限与有功出力有关严格来说应该是Q_max(P)但很多模型直接给定常数上限会带来结果偏乐观的问题。节点电压幅值约束U_min ≤ U ≤ U_max在碳中和场景下电压约束的重要性会被放大新能源出力的随机性让电压越限风险显著增加。变压器分接头调节范围约束。气网节点气压约束和管道输气能力约束。线路传输容量约束。在众多不等式约束里最难处理的是线路容量约束。这个约束用线性化近似会丢失无功信息用完整交流潮流又带来大量非线性。我的做法是保留交流潮流形式的视在功率约束配合松弛技术确保可解性。还有个容易被忽略的约束是P2G设备的运行约束。P2G启动需要时间电解槽的爬坡速率有限-P_p2g,ramp ≤ P_p2g,t - P_p2g,t-1 ≤ P_p2g,ramp如果忽略这个约束优化器会让P2G在相邻时段随意跳变结果在实际中根本执行不了——这在电网侧调度中也是同样的道理。3.4 碳排放约束的建模与处理碳约束我倾向于直接写成总量约束的形式Σ(E_G,i) ≤ E_total,max其中E_G,i是第i台机组的碳排放量E_total,max是系统设定的碳排放上限。在“碳中和”目标下这个上限应该逐年收紧。如果进一步要求模型能反映不同碳价情景下的运行策略差异还可以把碳约束转化为拉格朗日惩罚项或者KKT乘子。但对优化模型来说显式碳约束更直观且方便做边际排放成本的影子价格分析——影子价格本身就是很有用的调度决策辅助信息。4. MatLab代码实现路线从数学模型到可运行程序模型搭完之后怎么把它在Matlab里落地是另一个大问题。我分享一下自己的完整技术路线包括求解工具选型、代码架构和核心函数的实现思路。4.1 求解工具与算法选型这个问题先说结论我选的是YALMIP建模工具箱 二阶锥松弛SOCP 商业求解器CPLEX或Gurobi的组合方案。为什么这么选电气互联系统优化问题本质上是MINLP混合整数非线性规划直接求解全局最优几乎不可能。工程上主要有四条路线内点法直接求解连续松弛的NLP速度快但容易陷入局部最优且很难处理整数变量。智能算法如粒子群和遗传算法不依赖梯度能处理非线性但每次求解的稳定性差而且没有最优性保证。Benders分解把原问题分解成主问题和子问题适合大规模系统但模型推导复杂代码实现成本高。二阶锥松弛把非凸的Weymouth方程和交流潮流约束松弛为凸约束配合Big-M法处理离散变量形成MISOCP问题商用求解器对这类问题的求解能力已经非常成熟。我实测下来二阶锥松弛商业求解器在几十个节点的系统上求解速度在几十秒到几分钟之间兼顾了精度和速度是最适合Matlab实现的技术路线。具体做法是把Weymouth方程做二阶锥转换定义变量α_i p_i^2任意变换后约束变成下式FLOW_ij^2 ≤ K_ij^2 * (α_i - α_j)同时把潮流方程中的二次项通过电压幅值平方替换处理。注意这里的松弛方向是有讲究的——必须松弛成小于等于的形式才能保持凸性。虽然松弛扩大了可行域但实际运行中优化结果一般会自动逼近等式边界。4.2 代码架构分模块设计的思路我的Matlab代码按模块组织便于调试和复用主文件main.m定义系统参数加载数据调用建模函数求解输出结果。数据定义文件case_data.m电力系统节点参数、气网管道参数、耦合元件参数等基础数据。建模函数build_opf_model.m用YALMIP定义优化变量、目标函数和约束。结果分析函数plot_results.m绘制电压分布、机组出力、气网压力、碳排量变化等曲线。核心的YALMIP建模框架大概是下面的思路以目标函数和关键约束为例伪代码% YALMIP变量定义 P_g sdpvar(n_gen, 1); % 发电机有功 Q_g sdpvar(n_gen, 1); % 发电机无功 U sdpvar(n_bus, 1); % 节点电压幅值 theta sdpvar(n_bus, 1); % 节点相角 P_p2g sdpvar(n_p2g, 1); % P2G有功消耗 G_source sdpvar(n_source, 1); % 气源供气量 FLOW sdpvar(n_pipe, 1); % 管道流量 % 目标函数 Objective sum(a .* P_g.^2 b .* P_g c) ... sum(carbon_price .* emission_coeff .* P_g) ... sum(gamma .* P_p2g); % 约束定义 Constraints []; % 节点有功平衡 Constraints [Constraints, A_incidence * P_g P_wind - P_load - P_p2g_map ... U .* (G_ij * cos(theta_ij) B_ij * sin(theta_ij))]; % 节点无功平衡 Constraints [Constraints, Q_g Q_comp - Q_load ... U .* (G_ij * sin(theta_ij) - B_ij * cos(theta_ij))]; % 燃气机组气负荷 Constraints [Constraints, F_gas a_fuel * P_g.^2 b_fuel .* P_g c_fuel]; % 气网节点平衡 Constraints [Constraints, A_gas * FLOW G_source - G_load - F_gas 0]; % 二阶梯松弛 for i 1:n_pipe Constraints [Constraints, FLOW(i)^2 K(i)^2 * (alpha_from(i) - alpha_to(i))]; end % 求解 ops sdpsettings(solver,gurobi,verbose,1); optimize(Constraints, Objective, ops);这段代码里最关键的几个细节值得说明第一潮流和网络拓扑相关的矩阵G_ij、B_ij不是手工输入的而是用Matpower的makeYbus函数自动生成的。这样即使换算例代码不用改结构。第二P2G映射关系P_p2g_map的处理是我调试时踩过的一个坑P2G设备接在哪个节点上必须在数据文件里明确指定并且要形成从变量索引到网络节点的映射矩阵否则气网和电网的平衡方程会错位。第三二阶梯约束的写法里我用的是而不是这关系到问题是否凸。如果用等式约束YALMIP会把它当作非线性等式来解速度慢几个数量级还可能因为初始点差解不出来。4.3 针对整数变量的处理变压器分接头与无功补偿投切变压器分接头和无功补偿电容/电抗器的投切是离散变量我一开始直接把分接头位置定义为整数变量binvar或intvar导致求解时间暴增。后来改成Big-M法处理先把分接头比转换为连续变量t并添加线性约束t_min Σ(step_i * z_i) ≤ t ≤ t_min Σ(step_i * z_i) (1 - z_i) * M其中z_i是0-1变量M是足够大的数。这种做法把离散决策从“选位置”转换成“选步进”在Gurobi里求解效率高多了。还有一个实用技巧如果只是做日前调度变压器分接头一天内不宜频繁动作可以额外加一条动作次数约束Σ|z_i,t - z_i,t-1| ≤ N_max这属于“运行可行性约束”本质上是从实际运维规则中提炼出来的不加的话结果容易在相邻时段之间抖动。5. 算例验证与结果分析看到数据背后的调度含义模型建好、代码跑通这只是第一步。下面我用一个中等规模算例走一遍从数据设置到结果分析的全部过程。5.1 算例设置改进的电气互联测试系统我采用的测试系统是IEEE 39节点电力系统与比利时20节点天然气系统的耦合版本这是综合能源系统研究中比较经典的组合。关键耦合元件设置为4台燃气机组接入气网节点承担系统约30%的负荷。2套P2G设备分别安装在风电富集节点额定功率分别为20MW和30MW转换效率取60%。碳价设为100元/tCO2。风电渗透率设为30%以典型日出力曲线输入。无功补偿装置配置在负荷密集区节点单组容量5Mvar共12组。5.2 模型收敛性与求解性能分析在Matlab R2022b环境下用YALMIPGurobi 10.0求解模型规模大致是连续变量280多个0-1变量36个约束方程620多个。实测求解时间在45秒左右。这个速度在科研和离线规划场景下完全够用。对比不做SOCP松弛、直接上NLP内点法的方案那个方案在部分初始点下会卡在局部最优求解时间反而更长且结果不稳定。SOCP方案的优势非常明显。需要说明的是SOCP松弛在某些极端算例下可能产生非物理的松弛解。我在实际种遇到过一种情况气网某些管道气流方向不确定时Weymouth方程的绝对值展开会让松弛不精确。解决方案是对每根管道引入方向变量d_i当d_i1时气流为正d_i0时气流为负这种处理让结果更接近真实物理过程。5.3 关键结果解读从数据到运行策略首先看优化出的机组出力方案。在30%风电渗透率下燃气机组的出力不是简单地按经济性排序的因为碳约束的存在部分高排放燃煤机组会被压缩出力转而由碳强度更低的燃气机组顶上。这就是碳价信号的传导效果。再看无功优化的效果。对比“有功-无功协同优化”与“仅有功优化最后按固定功率因数校核”两种策略协同优化下的系统网损降低了约8.3%最低节点电压从0.93p.u.提升到0.97p.u.以上且所有节点电压均在安全范围内。这个对比很好地说明了我前面强调的观点——只在末端校核无功是无法主动优化电压分布的。P2G设备的运行策略也很有意思。风电大发时段P2G满发消耗富余电力相当于给系统提供了一个可调节的“电负荷”缓解了弃风问题同时产出的天然气进入管网提高了气网供气充裕度。而在负荷高峰时段P2G自动降出力甚至停机把电力优先让给用户。这种“削峰填谷”双向调节作用正是电气互联系统相比单一电网的独特优势。5.4 碳价灵敏度分析碳中和政策的量化评估我还做了碳价的灵敏度分析观察碳价从0元/tCO2逐步提高到300元/tCO2的过程中系统运行策略的变化规律碳价为0时系统优先选择煤电碳排放总量最高。碳价为50元/tCO2时燃气机组开始逐步替代煤电。碳价达到200元/tCO2时P2G设备利用率显著上升系统碳排量在原有基础上下降了约25%。这个分析对实际政策制定和投资规划很有参考价值——它把碳价这一抽象的政策变量转化成了看得见的机组启停次序和碳排总量的变化曲线。作为工程师能给出这样的量化分析跟只说“提高碳价有利于减排”是两种完全不同的说服力。6. 工程化落地中的避坑清单与扩展思路6.1 五个最容易翻车的细节全部来自实际调试第一个坑大气网系统的Weymouth方程数值病态问题。气压值的量级通常在几兆帕到十兆帕之间平方后数值差异巨大导致矩阵条件数极差求解器经常报数值问题。我尝试过把气压单位换成bar、把方程两边同时除以基准值等方法最有效的还是把所有物理量都标幺化处理让变量落在0.1到2这个区间内。第二个坑数据不一致性。电网和气网的基准功率、基准电压如果不统一耦合元件功率和气网的对应关系会直接算错。比如燃气机组效率用默认高估时会出现气网气源供气量不足而电网出力却满发的矛盾结果。第三个坑SOCP松弛方向不正确。之前说过必须是我还见过有人写成的那是在人为缩小可行域有可能直接导致无解。第四个坑YALMIP和求解器的版本兼容性。不同版本对optimize返回状态的定义有细微差别。我在调试时曾遇到过的困惑包括明明有可行解但Gurobi报infeasible是求解器设置问题改成双精度求解精度后收敛性明显改善。第五个坑约束冗余导致求解变慢。初始建模阶段为了安全起见添加了大量冗余约束结果导致LP松弛的求解效率降低了接近一倍。后来专门做了约束削减把可被其他约束自动满足的冗余项删掉这个问题才缓解。6.2 模型向更高维度扩展的三个方向目前的模型是日前稳态优化实际工程应用有三个扩展方向。引入时序耦合约束加入机组启停成本、爬坡约束、储能充放电约束把模型扩展为多时段动态优化。这在Matlab中的实现并不复杂只需把变量从一维向量变成二维矩阵并把相邻时段的耦合约束加进去。进一步细化天然气系统的动态特性比如加入管道存气、暂态惯性等。日前调度用稳态模型是合理的但日内滚动修正时动态信息能显著提高调度精度。引入不确定性鲁棒优化或随机优化思路风光出力的预测误差需要在模型里显式建模用鲁棒边界优化的方法把最坏情况下的安全约束纳入考虑。另外如果系统的无功电压问题比较突出建议进一步加入变电站动态无功补偿设备的详细模型比如SVC和STATCOM的运行约束——它们的快速调节能力在应对新能源功率波动时效果非常明显。6.3 关于代码复用和工程部署的一个建议最后这点是给想把模型用于实际项目的读者的。不建议把这里说的所谓“Matlab代码实现”停留在一次性脚本的层面。把整个优化模型封装成一个可重复调用的黑盒对外只暴露系统数据和场景参数对内统一调取YALMIP建模和求解器接口。这样后续做多场景比对、方案筛选甚至对接上级调度系统都会轻松得多。我的惯用做法是封装成一个function [P_g_opt, Q_g_opt, U_opt, F_gas_opt] run_IES_optimization(caseData, carbonPrice, windScenario, options)函数所有情景参数都是输入参数内部自动完成建模、求解和结果整理。这样跑碳价灵敏度分析、不同的风电出力场景只需要在主脚本里写循环调用不用反复改模型代码。这套工作流在我后续的好几个项目里都在不断复用。回头再看电气互联系统的有功-无功协同优化表面上是一个数学优化问题实际上考验的是对两条物理网络的深刻理解和建模功力。把模型建立在清晰的物理认知之上用合适的松弛手段让求解器高效工作再通过严谨的结果分析反哺对系统运行规律的认识这套方法论才是最值得沉淀的东西。