做无功优化的同行应该都有体会传统教材里讲“无功优化”基本等价于调电压、降网损目标函数写来写去就那么几项算例也基本都在纯电网环境里折腾。但这两年做实际项目风向完全变了拿过来的系统往往是电网、天然气网甚至热网拧在一起的综合能源系统目标也从“网损最小”变成了“碳排放、购电购气成本、新能源弃电”揉在一起的多目标问题。我这个项目标题看起来不短——“‘碳中和’目标下电气互联系统有功-无功协同优化模型”但本质上要解决的问题非常聚焦把电气互联系统中原本分开做的有功调度和无功电压控制放到同一个优化框架里让电网和气网的耦合设备燃气轮机、P2G真正参与调节最终给出一个能直接在Matlab里跑起来的求解方案。这篇博文会把模型怎么建、决策变量怎么编排、约束条件怎么处理、求解器怎么配以及我实际跑代码时踩过的坑完整过一遍。适合准备做电气互联系统优化研究的同学也适合工程上需要把协同优化落到代码里的工程师。1. 项目整体思路与核心问题拆解1.1 碳中和目标到底给无功优化加了什么新约束过去做无功优化的核心假设是电网里有一堆火电机组既能发有功也能发无功电压稳定性和频率稳定性基本靠它们撑着。但碳中和目标下的电源结构完全变了风机、光伏大规模并网这类新能源机组的可调无功能力很弱而且出力是间歇性的与此同时煤电被压缩碳排放约束又让火电要么降出力、要么上碳捕集导致传统“靠同步电源扛电压”的体系越来越吃力。更关键的是碳排放本身成了一种有价格的资源。当碳价进入调度模型后一台机组发多少有功不再只取决于煤价和气价还要看它每发一度电产生多少碳、要付多少碳税。这个变化直接影响了发电组合而发电组合变了系统的无功分布和电压水平又会跟着变。也就是说在碳中和背景下有功优化和无功优化已经不能被拆开看了——你调有功改变了机组出力组合无功支撑能力随之改变你调无功改变了电压分布网损和有功需求也跟着变。1.2 电气互联系统一张电网和一张气网互相牵制电气互联系统说白了就是电力系统和天然气系统通过耦合元件绑在一起。两个最典型的耦合设备是燃气轮机和P2G电转气燃气轮机烧天然气发电是“气转电”P2G用富余的电电解水制氢或合成天然气是“电转气”。打个比方天然气网像一条地下的能源管道网电网像一条悬在空中的电力网燃气轮机和P2G就是两个网之间的“换乘站”。换乘站的存在让两边不再是独立系统气网供气的压力、流量会限制燃气轮机的发电出力电网低谷时段的富余电力又可以通过P2G转成天然气反过来影响气网的流量分布。这种强耦合关系意味着任何一边的边界条件变了另一边都会受到传导影响。所以我在这套模型里没有把气网简化为一个无限供气的“理想气源”而是把它建成了带节点压力、管道流量、气源出力上限的完整网络模型。1.3 有功-无功为什么要协同优化而不是解耦分步算行业内以前有一种简化做法先做有功经济调度确定机组出力再做无功优化确定电压和变压器分接头。这个流程在电源结构单一、系统裕度大的时候勉强能用但在电气互联系统里有两个硬伤。第一个硬伤是电流的物理特性决定的。交流潮流方程里有功功率不仅跟电压相角有关也跟电压幅值有关无功功率同样同时受相角和幅值影响。把两个变量硬拆成两个顺序问题本质上等于忽略了它们之间的耦合项。第二个硬伤是耦合设备的行为没法归到单一类别。燃气轮机的发电有功受气网供气约束但它本身又是同步电机它的无功容量极限会随有功出力变化P-Q容量曲线。P2G运行时消耗有功本身也需要无功支撑等于在电网里平白多了一个无功负荷。这些设备如果不在同一个优化模型里同时看到“有功侧”和“无功侧”就没法合理调度它们。所以这个项目的核心思路是把两个网络、两个优化维度、两类决策变量一次性塞进一个优化问题里求解器同时给出有功出力、无功出力、节点电压、气源供气量和P2G运行功率而不是分阶段迭代。协同优化的本质不是“分步走”加一个反馈修正而是“一起算”。2. 协同优化模型目标函数与约束条件2.1 目标函数怎么设计经济、碳排、网损一次收敛大多数碳中和背景下的优化模型都不会做真正意义上的多目标加权模糊处理而是把碳排放折算成成本放进一个单目标函数里。我也采用这个思路目标函数形式如下min F ΣC_G,i(P_G,i) ΣC_gas,j(f_g,j) τ_carbon × E_total λ_loss × P_loss第一项是常规发电机的燃料成本采用二次函数C_G,i a_i·P² b_i·P c_i其中a、b、c是燃料成本系数。第二项是天然气网络中气源节点的购气成本一般用线性或分段线性函数单位和电网侧一样要折算成元/小时。第三项是碳排放成本τ_carbon是碳价元/吨E_total是系统总碳排放量吨/小时包括煤电排放、气电排放。第四项是网损惩罚项P_loss是全网的支路有功损耗λ_loss取一个电价折算系数代表“少损耗一兆瓦时电量相当于多卖一兆瓦时电”的边际收益。实际编程时有一点要注意目标函数里各项的数量级差异可能极大。机组二次成本动辄几千上万元每小时碳排放成本在小碳价下可能只有几百元每小时网损惩罚项乘出来可能只有几十元。如果不做任何处理求解器会把注意力放在量级最大的项上小幅改善大头成本而完全忽略网损和碳排的优化。我的处理办法是给每一项做归一化用各自初始点的值做基准把四项都压到同一量级再放进目标函数。2.2 电网侧的等式与不等式约束电网部分的约束主要来自交流潮流方程。对每个节点i有功和无功注入必须满足P_inj,i V_i · ΣV_j · (G_ij·cosθ_ij B_ij·sinθ_ij) Q_inj,i V_i · ΣV_j · (G_ij·sinθ_ij − B_ij·cosθ_ij)其中P_inj,i是该节点的发电机有功注入减去负荷有功需求Q_inj,i同理。G和B是节点导纳矩阵的实部和虚部θ_ij是节点i和节点j的相角差。这套方程是典型的非线性等式约束必须在优化问题里作为ceq返回给求解器。不等式约束则包括节点电压幅值上下限常规取0.95pu到1.05pu新能源接入比例高的节点建议把上限稍微放低比如1.02pu防止电压偏高触发逆变器脱网。发电机有功和无功出力上下限有功上限来自机组额定容量无功上限来自P-Q容量曲线我为了方便把无功上限做成了随有功出力变化的函数而不是固定常数。支路潮流传输极限如果跑稳态潮流后发现有支路过载可以在不等式约束里加入支路电流或视在功率限制。不过我实测下来这个约束会让雅可比矩阵的计算量明显变大如果算例场景不涉及严重过载可以暂时不启用。电网侧还有一个容易忽略的细节基准功率。我在模型里所有功率量都转换到标幺值pu下计算基准功率取100MVA但天然气网络里的功率和流量用的是实际单位MW、kg/s天然气热值、转换效率这些参数全都涉及单位换算。编程时我在数据文件里把“换算因子”集中定义避免每个函数里各写一套。2.3 天然气网络侧的约束建模天然气网络的建模是很多人第一次写这类模型时最头疼的部分。我用的方法是节点流量平衡加管道流量方程。每个气网节点的净流入流量必须为零Σf_mn W_supply,i − L_load,i − F_gt,i F_p2g,i 0其中f_mn是通过管道mn的流量W_supply是气源供气量L_load是纯气负荷F_gt是燃气轮机耗气F_p2g是P2G产气。这个等式中F_gt和F_p2g就是电气互联的耦合项。管道流量采用经典的Weymouth稳态方程f_mn sign(p_m² − p_n²) · C_mn · sqrt(|p_m² − p_n²|)C_mn是管道常数与管径、长度、摩阻系数有关p_m和p_n是管道两端节点压力。这个方程有两个特点非线性很强且有方向性——天然气的流向取决于压力差的正负。直接写进约束时要注意在Matlab里用sign乘以sqrt绝对值避免出现负数开根号的复数问题。气源节点还有一个压力上下限约束。气源不是“要多少给多少”压力上限限制最大供气能力压力下限保证管网末端用户的用气压力。我把气源压力、节点压力全部加了下限0.7pu上限1.2pu以初始稳态解为基准做的标幺化。2.4 燃气轮机和P2G耦合元件的建模细节耦合元件怎么建直接决定了模型能不能真实反映电气互联效果。燃气轮机这里我用线性近似把有功出力映射为耗气量F_gt α_gt β_gt · P_gt其中α_gt是空载气耗β_gt是增量气耗系数。这个线性近似对中小型燃气轮机够用如果要精确可以考虑二次函数但会增加不少非线性求解负担。P2G则反过来F_p2g η_p2g · P_p2g / H_NH_N是天然气高位热值η_p2g是电转气的综合转换效率。我在算例里用的效率是0.6即1MWh电转成天然气的热值约为0.6MWh当量。P2G的运行功率有上下限下限一般取额定功率的10%到20%因为电解槽设备很少允许零负荷运行。这两个耦合设备建在哪个节点上也有讲究。燃气轮机通常建在电网负荷中心节点同时连着气网的某个中间节点P2G则更适合建在新能源富集节点附近这样才能把本来要弃掉的电转成气存储。我在IEEE30节点算例里把P2G放在了风电接入节点燃气轮机放在最大负荷节点这样能最直观地看到电气互联带来的时空互补效果。3. Matlab实现代码架构与求解配置3.1 代码文件的组织方式这个项目我按“数据、模型、求解、可视化”四层来组织文件避免把所有逻辑堆在一个脚本里。整个目录结构如下|-- main.m % 主程序读取数据、编排变量、调用求解器、输出结果 |-- dataIEEE30Gas.m % IEEE30节点电网数据 6节点天然气网数据 |-- objFun.m % 目标函数输入x输出F |-- nonlcon.m % 约束函数输入x输出c和ceq |-- extractX.m % 把决策向量x拆解成Pg、Qg、V、theta等子向量 |-- plotResult.m % 电压剖面、机组出力、气源流量绘图main.m里第一步调用dataIEEE30Gas.m加载数据第二步构建决策变量的上下界lb和ub第三步设置初始点x0第四步调用fmincon。这里强烈建议把“决策变量顺序”写成一个注释表放在main.m顶部比如“x(1:6)是6台发电机的有功x(7:12)是无功x(13:42)是30个节点的电压幅值……”不然调试时你会无数次忘记第37个变量到底代表什么。3.2 决策变量编排与索引封装决策变量向量的编排顺序看起来简单实际是这类代码最容易翻车的地方。我的编排如下% x [Pg(1:NG); Qg(1:NG); V(1:Nbus); theta(1:Nbus); ... % P_p2g; F_supply(1:NGasSource); p_node(1:NGasNode)]电力系统部分6台发电机的有功出力、6台发电机的无功出力30个节点的电压幅值标幺值、30个节点的相角弧度。电气互联系统部分1个P2G的运行功率2个气源节点的供气量6个气网节点的压力。总共是66303012681个决策变量。为了在目标函数和约束函数里不每次重复写索引我封装了一个extractX函数function [Pg, Qg, V, theta, P_p2g, F_sup, pNode] extractX(x) % 电力部分 Pg x(1:6); Qg x(7:12); V x(13:42); theta x(43:72); % 互联与气网部分 P_p2g x(73); F_sup x(74:75); pNode x(76:81); end这样做的好处是如果后续你想加储能、加需求响应、加氢能设备只需要扩展一个子向量和对应的索引注释目标函数和约束函数里不用大改。3.3 目标函数与约束函数的关键代码段目标函数里最核心的是碳排放计算和网损计算。碳排放我这里只算供给侧的直接排放煤电排放系数乘以煤电有功出力加上燃气轮机耗气量对应的排放系数。风电和光伏是零碳排。function F objFun(x) % 提取决策变量 [Pg, Qg, V, theta, P_p2g, F_sup, pNode] extractX(x); % 机组燃料成本(元/h)a*P^2 b*P c genCost sum(a .* Pg.^2 b .* Pg c); % 气源购气成本(元/h) gasCost sum(gasPrice .* F_sup); % 碳排放成本(元/h) E_coal EF_coal .* Pg(1:coalGen); % 煤电排放t/MWh E_gas EF_gas .* F_gt_cons; % 燃气耗气排放t/MWh当量 E_total sum(E_coal) sum(E_gas); carbonCost carbonPrice * E_total; % 网损惩罚项先算各支路损耗 lossMW calcLoss(V, theta, Ybus); lossPenalty lambda_loss * lossMW; % 汇总 F genCost gasCost carbonCost lossPenalty; end约束函数里电网功率平衡用2×Nbus个等式约束气网节点流量平衡用NGasNode个等式约束电压和压力上下限用不等式约束。核心代码结构如下function [c, ceq] nonlcon(x) [Pg, Qg, V, theta, P_p2g, F_sup, pNode] extractX(x); % 电网潮流等式节点注入功率等于潮流方程计算值 [P_calc, Q_calc] calcPowerFlow(V, theta, Ybus); ceq(1:Nbus) P_inj(x) - P_calc; ceq(Nbus1:2*Nbus) Q_inj(x) - Q_calc; % 气网流量平衡 ceq(2*Nbus1:2*NbusNGasNode) gasFlowBalance(x); % 不等式约束电压上下限、压力上下限、支路限流等 c [V - Vmax; Vmin - V; pNode - pMax; pMin - pNode; ...]; end这里我专门提一句ceq和c的顺序不能乱fmincon靠索引判断哪个约束不满足所以约束函数里每一行写清楚对应哪个节点哪个方程非常重要。调试时可以临时在nonlcon里设置断点打印max(abs(ceq))快速定位是哪一类方程不满足。3.4 求解器选型与参数配置对这类带非线性等式约束和不等式约束的优化问题Matlab内置的fmincon是最省事的选择。我用的求解配置如下options optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... MaxIterations, 3000, ... MaxFunctionEvaluations, 50000, ... OptimalityTolerance, 1e-6, ... StepTolerance, 1e-8, ... ConstraintTolerance, 1e-6);选内点法而不是SQP是因为这类模型约束规模不小我这里有60个潮流等式约束加6个气网等式约束内点法在处理大规模稀疏非线性规划时收敛性更稳定。SQP的优点是迭代步长比较直接但碰到潮流方程这种强非线性约束容易在迭代途中跑到无解区域反而不如内点法的“障碍项”来得平滑。求解器选型也可以做对比实验求解方案适用场景优缺点fmincon内点法中小规模、约束复杂收敛稳但易陷入局部最优fmincon SQP变量少、约束简单迭代快但对初值敏感YALMIP外部求解器需要形式化建模可读性好但安装配置复杂遗传/粒子群内点法大规模多局部极值全局搜索强但耗时高对于81个变量的模型内点法通常在几十秒内能收敛。我建议先用内点法跑通如果后续扩展成动态多时段模型、变量上千再考虑混合求解策略先用粒子群粗搜一组好初值再交给内点法精修。4. 算例设计与结果分析4.1 测试系统选择与数据修改算例系统我用了修改版IEEE30节点电网加一个6节点天然气网。IEEE30节点系统是电力系统优化的标准测试系统支路参数、机组参数都有公开数据方便复现。气网部分我是按参考文献里的6节点天然气系统改的两个气源分别放在气网节点1和节点6燃气轮机接在电网最大负荷节点对应的气网节点上P2G接在电网风电接入节点附近。修改的地方主要有三处一是加大了新能源渗透率在节点12和节点22各加了一个风电场装机容量分别是50MW和80MW二是把原有6台火电机组中的一台替换为燃气轮机让它通过气网供气三是把常规负荷调高到和新能源装机匹配的程度制造“低谷期风机出力接近负荷”的场景这样P2G才能发挥作用。所有计算都在标幺值下进行基准容量100MVA。天然气网络里的压力基准是5MPa流量基准按能量流换算成MW。这一步换算要特别小心我在代码里写了十几个“单位换算常量”比如1kg天然气约等于13.9kWh热值管道常数C_mn是直接抄文献还是自己算对结果影响很大。我建议文献里的管道参数能直接用的尽量用别自己推导否则调一个晚上都调不平。4.2 三种场景对比从解耦到协同的梯度验证为了说明协同优化的价值我设计了三个场景做递进式对比场景A独立无功优化。只优化无功固定有功调度结果为初始值目标只有网损和电压偏差不接气网、不算碳排。场景B有功经济调度碳税但气网固定。有功调度考虑煤电、气电和碳价但天然气网络只是简单一个固定供气成本的理想源不做压力约束。场景C本文完整模型。电气互联有功无功协同碳税P2G可调。三种场景的负荷曲线、风电出力完全一致碳价统一取100元/吨。跑完后的结果对比如下指标场景A独立无功场景B有碳税但气网理想化场景C完整协同模型系统总运行成本元/h入口条件固定约3.12万2.86万2.71万碳排放总量t/h约42.536.832.1全网有功损耗MW8.99.67.3最低节点电压pu0.980.940.97求解时间s82147这个结果很有代表性。场景A虽然网损控制得不错但因为固定了有功出力碳排放成本完全没被优化总成本是最高的。场景B引入了碳税发电组合向气电转移碳排明显下降但因为气网被理想化没有约束燃气轮机的实际供气能力同时也没靠P2G调整新能源出力导致最低电压差点越限。场景C把电气互联、P2G、天然气约束全部纳入后低谷时段富余风电被P2G消纳转化成气峰荷时段燃气轮机再把这部分气发成电碳排进一步下降网损反而降到最低。4.3 从结果倒推协同优化为什么能同时省钱又省碳看结果不能只看表还得解释机理。场景C的最优解里P2G的出力曲线和风电曲线高度一致夜间风电功率高、负荷低P2G启动吸收多余电力把电能变成天然气存储白天负荷上升燃气轮机用气发电补足缺口。等效于气网成了电网的“储能池”。无功侧的价值体现在电压水平上。场景C的电压最低点0.97pu明显高于场景B的0.94pu。原因有两方面一是P2G接入点比纯风电接入点多了无功支撑手段二是协同优化让燃气轮机在有功出力选择时兼顾了其无功能力没有单纯为了降煤耗把某台机组压到过低出力导致无功供给不足。另外我还做了一个敏感性实验把碳价从0元/吨摸到200元/吨观察碳排放总量的变化。结果是碳价越高P2G出力越大燃气轮机占比越高碳排放下降但气源购气成本上升。这个结果说明模型的行为符合预期碳价是驱动电气互联系统“电气转气再转电”空间套利的关键信号。5. 常见问题与调试实录5.1 潮流方程怎么调都收敛不了这是这类项目里最常遇到的坑我自己的经验是大概率不是模型错了而是初始点距离可行域太远。fmincon内点法虽然鲁棒但它在一个完全不可行的初始点附近经常陷入“约束矛盾导致步长缩小到接近零”的僵局。解决办法是“热启动”先用一个简化的牛顿法潮流计算程序算出常规电网的电压和相角解然后把潮流解作为优化问题的初始点。这个做法能减少95%以上的初值问题。具体来说先用mpower这类工具或自己写一个极坐标牛顿法求解一次纯潮流得到V0和theta0再把它塞进x0的对应位置。气网部分同理先把气网稳态流量算一遍得到压力初值。如果真的碰到“可行解不存在”先检查约束之间是否矛盾。常见原因是变压器变比上下限卡得太死或者气源最大供气量不够燃气轮机发电需求。调试时我用过一个土办法把全部不等式约束的上限放大10倍看等式约束能不能满足。如果放大后等式还是不平说明等式约束本身写错了如果等式平了但不等式调不回来说明可行域真的不存在。5.2 目标函数各项量级差太多导致优化方向被绑架前面说过成本项、碳排项、网损项量级差异大。如果你发现优化结果里网损没变小、碳排放也没变低只有发电燃料成本在降大概率就是量级问题。这一步没有捷径就是把每一项在初始点的值打印出来看看确认它们都在同一量级。我习惯把四项成本在初始点的值除以总成本换算成百分比权重然后手动调节归一化系数让四项占比大致在40%、30%、20%、10%这个区间再交给优化器。5.3 matlab安装与版本兼容性optimoptions参数报错也是坑这个项目依赖Matlab优化工具箱我建议用R2023b以上版本新版本的内点法实现更稳定而且optimoptions函数对各算法的参数校验更严格。严格是好事但也会带来一个常见问题你在老版本写的options参数名到新版本会直接报错比如“MaxFunEvals已更改为MaxFunctionEvaluations”。这种报错网上搜索一下就能解决但别以为是代码逻辑问题。还要注意一个正版环境的问题Matlab的许可证激活偶发报错尤其在学校或单位把实验室的旧license文件混用的情况。MathWorks官方下载安装和激活的工具链其实很成熟直接从官网或学校正版授权渠道走流程就好。我自己遇到过装完启动后闪退的情况最后定位是电脑上旧版本MATLAB的路径环境变量残留把注册表里老版本路径清理干净就正常了。正版授权机制虽然偶尔折腾但比用来路不明的所谓绿色版稳定得多排查起来也更有据可循。5.4 求解时间过长分阶段初值与变量松弛81个变量其实不大但如果你把天然气压力、所有节点电压、相角全部同时作为优化变量某些matlab版本的内点法会比较慢尤其是约束函数里潮流方程要重新组装一次的时候。实测来看61秒到2分钟是正常的但如果超过5分钟大概率是约束函数写得不高效。排查技巧约束函数里不要每次都重新算节点导纳矩阵Ybus把它设为全局变量或嵌套函数共享变量只算一次。还有就是如果发现某次求解耗时爆炸式增长可以尝试把部分等式约束转化成不等式约束的“松弛”形式比如把“节点注入功率等于潮流计算值”放宽为“注入等于计算值的误差小于1e-4”。这样性能快很多代价是结果精度略降但用于方案预研完全够。一点个人经验我实际做完这个项目的体会是协同优化模型的代码并没有比传统无功优化复杂多少真正的成本在建模细节和调试。数据单位、决策变量索引、初始点这三点控制好你就能从一个“看起来很高大上”的数学模型得到一个能稳定收敛的实用程序。如果准备自己复现我强烈建议先用一个小系统比如5节点或14节点电网加3节点气网把框架跑通再上IEEE30这个规模。这个模型后续还可以继续扩展把单时段变成24小时滚动优化加入储能和氢能设备碳价改成阶梯式碳税都是现成可以往下做的方向。电气互联系统的调度问题跟着碳中和的节奏走后面需要解决的场景会越来越多。
碳中和下电气互联系统有功-无功协同优化模型及Matlab实现
做无功优化的同行应该都有体会传统教材里讲“无功优化”基本等价于调电压、降网损目标函数写来写去就那么几项算例也基本都在纯电网环境里折腾。但这两年做实际项目风向完全变了拿过来的系统往往是电网、天然气网甚至热网拧在一起的综合能源系统目标也从“网损最小”变成了“碳排放、购电购气成本、新能源弃电”揉在一起的多目标问题。我这个项目标题看起来不短——“‘碳中和’目标下电气互联系统有功-无功协同优化模型”但本质上要解决的问题非常聚焦把电气互联系统中原本分开做的有功调度和无功电压控制放到同一个优化框架里让电网和气网的耦合设备燃气轮机、P2G真正参与调节最终给出一个能直接在Matlab里跑起来的求解方案。这篇博文会把模型怎么建、决策变量怎么编排、约束条件怎么处理、求解器怎么配以及我实际跑代码时踩过的坑完整过一遍。适合准备做电气互联系统优化研究的同学也适合工程上需要把协同优化落到代码里的工程师。1. 项目整体思路与核心问题拆解1.1 碳中和目标到底给无功优化加了什么新约束过去做无功优化的核心假设是电网里有一堆火电机组既能发有功也能发无功电压稳定性和频率稳定性基本靠它们撑着。但碳中和目标下的电源结构完全变了风机、光伏大规模并网这类新能源机组的可调无功能力很弱而且出力是间歇性的与此同时煤电被压缩碳排放约束又让火电要么降出力、要么上碳捕集导致传统“靠同步电源扛电压”的体系越来越吃力。更关键的是碳排放本身成了一种有价格的资源。当碳价进入调度模型后一台机组发多少有功不再只取决于煤价和气价还要看它每发一度电产生多少碳、要付多少碳税。这个变化直接影响了发电组合而发电组合变了系统的无功分布和电压水平又会跟着变。也就是说在碳中和背景下有功优化和无功优化已经不能被拆开看了——你调有功改变了机组出力组合无功支撑能力随之改变你调无功改变了电压分布网损和有功需求也跟着变。1.2 电气互联系统一张电网和一张气网互相牵制电气互联系统说白了就是电力系统和天然气系统通过耦合元件绑在一起。两个最典型的耦合设备是燃气轮机和P2G电转气燃气轮机烧天然气发电是“气转电”P2G用富余的电电解水制氢或合成天然气是“电转气”。打个比方天然气网像一条地下的能源管道网电网像一条悬在空中的电力网燃气轮机和P2G就是两个网之间的“换乘站”。换乘站的存在让两边不再是独立系统气网供气的压力、流量会限制燃气轮机的发电出力电网低谷时段的富余电力又可以通过P2G转成天然气反过来影响气网的流量分布。这种强耦合关系意味着任何一边的边界条件变了另一边都会受到传导影响。所以我在这套模型里没有把气网简化为一个无限供气的“理想气源”而是把它建成了带节点压力、管道流量、气源出力上限的完整网络模型。1.3 有功-无功为什么要协同优化而不是解耦分步算行业内以前有一种简化做法先做有功经济调度确定机组出力再做无功优化确定电压和变压器分接头。这个流程在电源结构单一、系统裕度大的时候勉强能用但在电气互联系统里有两个硬伤。第一个硬伤是电流的物理特性决定的。交流潮流方程里有功功率不仅跟电压相角有关也跟电压幅值有关无功功率同样同时受相角和幅值影响。把两个变量硬拆成两个顺序问题本质上等于忽略了它们之间的耦合项。第二个硬伤是耦合设备的行为没法归到单一类别。燃气轮机的发电有功受气网供气约束但它本身又是同步电机它的无功容量极限会随有功出力变化P-Q容量曲线。P2G运行时消耗有功本身也需要无功支撑等于在电网里平白多了一个无功负荷。这些设备如果不在同一个优化模型里同时看到“有功侧”和“无功侧”就没法合理调度它们。所以这个项目的核心思路是把两个网络、两个优化维度、两类决策变量一次性塞进一个优化问题里求解器同时给出有功出力、无功出力、节点电压、气源供气量和P2G运行功率而不是分阶段迭代。协同优化的本质不是“分步走”加一个反馈修正而是“一起算”。2. 协同优化模型目标函数与约束条件2.1 目标函数怎么设计经济、碳排、网损一次收敛大多数碳中和背景下的优化模型都不会做真正意义上的多目标加权模糊处理而是把碳排放折算成成本放进一个单目标函数里。我也采用这个思路目标函数形式如下min F ΣC_G,i(P_G,i) ΣC_gas,j(f_g,j) τ_carbon × E_total λ_loss × P_loss第一项是常规发电机的燃料成本采用二次函数C_G,i a_i·P² b_i·P c_i其中a、b、c是燃料成本系数。第二项是天然气网络中气源节点的购气成本一般用线性或分段线性函数单位和电网侧一样要折算成元/小时。第三项是碳排放成本τ_carbon是碳价元/吨E_total是系统总碳排放量吨/小时包括煤电排放、气电排放。第四项是网损惩罚项P_loss是全网的支路有功损耗λ_loss取一个电价折算系数代表“少损耗一兆瓦时电量相当于多卖一兆瓦时电”的边际收益。实际编程时有一点要注意目标函数里各项的数量级差异可能极大。机组二次成本动辄几千上万元每小时碳排放成本在小碳价下可能只有几百元每小时网损惩罚项乘出来可能只有几十元。如果不做任何处理求解器会把注意力放在量级最大的项上小幅改善大头成本而完全忽略网损和碳排的优化。我的处理办法是给每一项做归一化用各自初始点的值做基准把四项都压到同一量级再放进目标函数。2.2 电网侧的等式与不等式约束电网部分的约束主要来自交流潮流方程。对每个节点i有功和无功注入必须满足P_inj,i V_i · ΣV_j · (G_ij·cosθ_ij B_ij·sinθ_ij) Q_inj,i V_i · ΣV_j · (G_ij·sinθ_ij − B_ij·cosθ_ij)其中P_inj,i是该节点的发电机有功注入减去负荷有功需求Q_inj,i同理。G和B是节点导纳矩阵的实部和虚部θ_ij是节点i和节点j的相角差。这套方程是典型的非线性等式约束必须在优化问题里作为ceq返回给求解器。不等式约束则包括节点电压幅值上下限常规取0.95pu到1.05pu新能源接入比例高的节点建议把上限稍微放低比如1.02pu防止电压偏高触发逆变器脱网。发电机有功和无功出力上下限有功上限来自机组额定容量无功上限来自P-Q容量曲线我为了方便把无功上限做成了随有功出力变化的函数而不是固定常数。支路潮流传输极限如果跑稳态潮流后发现有支路过载可以在不等式约束里加入支路电流或视在功率限制。不过我实测下来这个约束会让雅可比矩阵的计算量明显变大如果算例场景不涉及严重过载可以暂时不启用。电网侧还有一个容易忽略的细节基准功率。我在模型里所有功率量都转换到标幺值pu下计算基准功率取100MVA但天然气网络里的功率和流量用的是实际单位MW、kg/s天然气热值、转换效率这些参数全都涉及单位换算。编程时我在数据文件里把“换算因子”集中定义避免每个函数里各写一套。2.3 天然气网络侧的约束建模天然气网络的建模是很多人第一次写这类模型时最头疼的部分。我用的方法是节点流量平衡加管道流量方程。每个气网节点的净流入流量必须为零Σf_mn W_supply,i − L_load,i − F_gt,i F_p2g,i 0其中f_mn是通过管道mn的流量W_supply是气源供气量L_load是纯气负荷F_gt是燃气轮机耗气F_p2g是P2G产气。这个等式中F_gt和F_p2g就是电气互联的耦合项。管道流量采用经典的Weymouth稳态方程f_mn sign(p_m² − p_n²) · C_mn · sqrt(|p_m² − p_n²|)C_mn是管道常数与管径、长度、摩阻系数有关p_m和p_n是管道两端节点压力。这个方程有两个特点非线性很强且有方向性——天然气的流向取决于压力差的正负。直接写进约束时要注意在Matlab里用sign乘以sqrt绝对值避免出现负数开根号的复数问题。气源节点还有一个压力上下限约束。气源不是“要多少给多少”压力上限限制最大供气能力压力下限保证管网末端用户的用气压力。我把气源压力、节点压力全部加了下限0.7pu上限1.2pu以初始稳态解为基准做的标幺化。2.4 燃气轮机和P2G耦合元件的建模细节耦合元件怎么建直接决定了模型能不能真实反映电气互联效果。燃气轮机这里我用线性近似把有功出力映射为耗气量F_gt α_gt β_gt · P_gt其中α_gt是空载气耗β_gt是增量气耗系数。这个线性近似对中小型燃气轮机够用如果要精确可以考虑二次函数但会增加不少非线性求解负担。P2G则反过来F_p2g η_p2g · P_p2g / H_NH_N是天然气高位热值η_p2g是电转气的综合转换效率。我在算例里用的效率是0.6即1MWh电转成天然气的热值约为0.6MWh当量。P2G的运行功率有上下限下限一般取额定功率的10%到20%因为电解槽设备很少允许零负荷运行。这两个耦合设备建在哪个节点上也有讲究。燃气轮机通常建在电网负荷中心节点同时连着气网的某个中间节点P2G则更适合建在新能源富集节点附近这样才能把本来要弃掉的电转成气存储。我在IEEE30节点算例里把P2G放在了风电接入节点燃气轮机放在最大负荷节点这样能最直观地看到电气互联带来的时空互补效果。3. Matlab实现代码架构与求解配置3.1 代码文件的组织方式这个项目我按“数据、模型、求解、可视化”四层来组织文件避免把所有逻辑堆在一个脚本里。整个目录结构如下|-- main.m % 主程序读取数据、编排变量、调用求解器、输出结果 |-- dataIEEE30Gas.m % IEEE30节点电网数据 6节点天然气网数据 |-- objFun.m % 目标函数输入x输出F |-- nonlcon.m % 约束函数输入x输出c和ceq |-- extractX.m % 把决策向量x拆解成Pg、Qg、V、theta等子向量 |-- plotResult.m % 电压剖面、机组出力、气源流量绘图main.m里第一步调用dataIEEE30Gas.m加载数据第二步构建决策变量的上下界lb和ub第三步设置初始点x0第四步调用fmincon。这里强烈建议把“决策变量顺序”写成一个注释表放在main.m顶部比如“x(1:6)是6台发电机的有功x(7:12)是无功x(13:42)是30个节点的电压幅值……”不然调试时你会无数次忘记第37个变量到底代表什么。3.2 决策变量编排与索引封装决策变量向量的编排顺序看起来简单实际是这类代码最容易翻车的地方。我的编排如下% x [Pg(1:NG); Qg(1:NG); V(1:Nbus); theta(1:Nbus); ... % P_p2g; F_supply(1:NGasSource); p_node(1:NGasNode)]电力系统部分6台发电机的有功出力、6台发电机的无功出力30个节点的电压幅值标幺值、30个节点的相角弧度。电气互联系统部分1个P2G的运行功率2个气源节点的供气量6个气网节点的压力。总共是66303012681个决策变量。为了在目标函数和约束函数里不每次重复写索引我封装了一个extractX函数function [Pg, Qg, V, theta, P_p2g, F_sup, pNode] extractX(x) % 电力部分 Pg x(1:6); Qg x(7:12); V x(13:42); theta x(43:72); % 互联与气网部分 P_p2g x(73); F_sup x(74:75); pNode x(76:81); end这样做的好处是如果后续你想加储能、加需求响应、加氢能设备只需要扩展一个子向量和对应的索引注释目标函数和约束函数里不用大改。3.3 目标函数与约束函数的关键代码段目标函数里最核心的是碳排放计算和网损计算。碳排放我这里只算供给侧的直接排放煤电排放系数乘以煤电有功出力加上燃气轮机耗气量对应的排放系数。风电和光伏是零碳排。function F objFun(x) % 提取决策变量 [Pg, Qg, V, theta, P_p2g, F_sup, pNode] extractX(x); % 机组燃料成本(元/h)a*P^2 b*P c genCost sum(a .* Pg.^2 b .* Pg c); % 气源购气成本(元/h) gasCost sum(gasPrice .* F_sup); % 碳排放成本(元/h) E_coal EF_coal .* Pg(1:coalGen); % 煤电排放t/MWh E_gas EF_gas .* F_gt_cons; % 燃气耗气排放t/MWh当量 E_total sum(E_coal) sum(E_gas); carbonCost carbonPrice * E_total; % 网损惩罚项先算各支路损耗 lossMW calcLoss(V, theta, Ybus); lossPenalty lambda_loss * lossMW; % 汇总 F genCost gasCost carbonCost lossPenalty; end约束函数里电网功率平衡用2×Nbus个等式约束气网节点流量平衡用NGasNode个等式约束电压和压力上下限用不等式约束。核心代码结构如下function [c, ceq] nonlcon(x) [Pg, Qg, V, theta, P_p2g, F_sup, pNode] extractX(x); % 电网潮流等式节点注入功率等于潮流方程计算值 [P_calc, Q_calc] calcPowerFlow(V, theta, Ybus); ceq(1:Nbus) P_inj(x) - P_calc; ceq(Nbus1:2*Nbus) Q_inj(x) - Q_calc; % 气网流量平衡 ceq(2*Nbus1:2*NbusNGasNode) gasFlowBalance(x); % 不等式约束电压上下限、压力上下限、支路限流等 c [V - Vmax; Vmin - V; pNode - pMax; pMin - pNode; ...]; end这里我专门提一句ceq和c的顺序不能乱fmincon靠索引判断哪个约束不满足所以约束函数里每一行写清楚对应哪个节点哪个方程非常重要。调试时可以临时在nonlcon里设置断点打印max(abs(ceq))快速定位是哪一类方程不满足。3.4 求解器选型与参数配置对这类带非线性等式约束和不等式约束的优化问题Matlab内置的fmincon是最省事的选择。我用的求解配置如下options optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... MaxIterations, 3000, ... MaxFunctionEvaluations, 50000, ... OptimalityTolerance, 1e-6, ... StepTolerance, 1e-8, ... ConstraintTolerance, 1e-6);选内点法而不是SQP是因为这类模型约束规模不小我这里有60个潮流等式约束加6个气网等式约束内点法在处理大规模稀疏非线性规划时收敛性更稳定。SQP的优点是迭代步长比较直接但碰到潮流方程这种强非线性约束容易在迭代途中跑到无解区域反而不如内点法的“障碍项”来得平滑。求解器选型也可以做对比实验求解方案适用场景优缺点fmincon内点法中小规模、约束复杂收敛稳但易陷入局部最优fmincon SQP变量少、约束简单迭代快但对初值敏感YALMIP外部求解器需要形式化建模可读性好但安装配置复杂遗传/粒子群内点法大规模多局部极值全局搜索强但耗时高对于81个变量的模型内点法通常在几十秒内能收敛。我建议先用内点法跑通如果后续扩展成动态多时段模型、变量上千再考虑混合求解策略先用粒子群粗搜一组好初值再交给内点法精修。4. 算例设计与结果分析4.1 测试系统选择与数据修改算例系统我用了修改版IEEE30节点电网加一个6节点天然气网。IEEE30节点系统是电力系统优化的标准测试系统支路参数、机组参数都有公开数据方便复现。气网部分我是按参考文献里的6节点天然气系统改的两个气源分别放在气网节点1和节点6燃气轮机接在电网最大负荷节点对应的气网节点上P2G接在电网风电接入节点附近。修改的地方主要有三处一是加大了新能源渗透率在节点12和节点22各加了一个风电场装机容量分别是50MW和80MW二是把原有6台火电机组中的一台替换为燃气轮机让它通过气网供气三是把常规负荷调高到和新能源装机匹配的程度制造“低谷期风机出力接近负荷”的场景这样P2G才能发挥作用。所有计算都在标幺值下进行基准容量100MVA。天然气网络里的压力基准是5MPa流量基准按能量流换算成MW。这一步换算要特别小心我在代码里写了十几个“单位换算常量”比如1kg天然气约等于13.9kWh热值管道常数C_mn是直接抄文献还是自己算对结果影响很大。我建议文献里的管道参数能直接用的尽量用别自己推导否则调一个晚上都调不平。4.2 三种场景对比从解耦到协同的梯度验证为了说明协同优化的价值我设计了三个场景做递进式对比场景A独立无功优化。只优化无功固定有功调度结果为初始值目标只有网损和电压偏差不接气网、不算碳排。场景B有功经济调度碳税但气网固定。有功调度考虑煤电、气电和碳价但天然气网络只是简单一个固定供气成本的理想源不做压力约束。场景C本文完整模型。电气互联有功无功协同碳税P2G可调。三种场景的负荷曲线、风电出力完全一致碳价统一取100元/吨。跑完后的结果对比如下指标场景A独立无功场景B有碳税但气网理想化场景C完整协同模型系统总运行成本元/h入口条件固定约3.12万2.86万2.71万碳排放总量t/h约42.536.832.1全网有功损耗MW8.99.67.3最低节点电压pu0.980.940.97求解时间s82147这个结果很有代表性。场景A虽然网损控制得不错但因为固定了有功出力碳排放成本完全没被优化总成本是最高的。场景B引入了碳税发电组合向气电转移碳排明显下降但因为气网被理想化没有约束燃气轮机的实际供气能力同时也没靠P2G调整新能源出力导致最低电压差点越限。场景C把电气互联、P2G、天然气约束全部纳入后低谷时段富余风电被P2G消纳转化成气峰荷时段燃气轮机再把这部分气发成电碳排进一步下降网损反而降到最低。4.3 从结果倒推协同优化为什么能同时省钱又省碳看结果不能只看表还得解释机理。场景C的最优解里P2G的出力曲线和风电曲线高度一致夜间风电功率高、负荷低P2G启动吸收多余电力把电能变成天然气存储白天负荷上升燃气轮机用气发电补足缺口。等效于气网成了电网的“储能池”。无功侧的价值体现在电压水平上。场景C的电压最低点0.97pu明显高于场景B的0.94pu。原因有两方面一是P2G接入点比纯风电接入点多了无功支撑手段二是协同优化让燃气轮机在有功出力选择时兼顾了其无功能力没有单纯为了降煤耗把某台机组压到过低出力导致无功供给不足。另外我还做了一个敏感性实验把碳价从0元/吨摸到200元/吨观察碳排放总量的变化。结果是碳价越高P2G出力越大燃气轮机占比越高碳排放下降但气源购气成本上升。这个结果说明模型的行为符合预期碳价是驱动电气互联系统“电气转气再转电”空间套利的关键信号。5. 常见问题与调试实录5.1 潮流方程怎么调都收敛不了这是这类项目里最常遇到的坑我自己的经验是大概率不是模型错了而是初始点距离可行域太远。fmincon内点法虽然鲁棒但它在一个完全不可行的初始点附近经常陷入“约束矛盾导致步长缩小到接近零”的僵局。解决办法是“热启动”先用一个简化的牛顿法潮流计算程序算出常规电网的电压和相角解然后把潮流解作为优化问题的初始点。这个做法能减少95%以上的初值问题。具体来说先用mpower这类工具或自己写一个极坐标牛顿法求解一次纯潮流得到V0和theta0再把它塞进x0的对应位置。气网部分同理先把气网稳态流量算一遍得到压力初值。如果真的碰到“可行解不存在”先检查约束之间是否矛盾。常见原因是变压器变比上下限卡得太死或者气源最大供气量不够燃气轮机发电需求。调试时我用过一个土办法把全部不等式约束的上限放大10倍看等式约束能不能满足。如果放大后等式还是不平说明等式约束本身写错了如果等式平了但不等式调不回来说明可行域真的不存在。5.2 目标函数各项量级差太多导致优化方向被绑架前面说过成本项、碳排项、网损项量级差异大。如果你发现优化结果里网损没变小、碳排放也没变低只有发电燃料成本在降大概率就是量级问题。这一步没有捷径就是把每一项在初始点的值打印出来看看确认它们都在同一量级。我习惯把四项成本在初始点的值除以总成本换算成百分比权重然后手动调节归一化系数让四项占比大致在40%、30%、20%、10%这个区间再交给优化器。5.3 matlab安装与版本兼容性optimoptions参数报错也是坑这个项目依赖Matlab优化工具箱我建议用R2023b以上版本新版本的内点法实现更稳定而且optimoptions函数对各算法的参数校验更严格。严格是好事但也会带来一个常见问题你在老版本写的options参数名到新版本会直接报错比如“MaxFunEvals已更改为MaxFunctionEvaluations”。这种报错网上搜索一下就能解决但别以为是代码逻辑问题。还要注意一个正版环境的问题Matlab的许可证激活偶发报错尤其在学校或单位把实验室的旧license文件混用的情况。MathWorks官方下载安装和激活的工具链其实很成熟直接从官网或学校正版授权渠道走流程就好。我自己遇到过装完启动后闪退的情况最后定位是电脑上旧版本MATLAB的路径环境变量残留把注册表里老版本路径清理干净就正常了。正版授权机制虽然偶尔折腾但比用来路不明的所谓绿色版稳定得多排查起来也更有据可循。5.4 求解时间过长分阶段初值与变量松弛81个变量其实不大但如果你把天然气压力、所有节点电压、相角全部同时作为优化变量某些matlab版本的内点法会比较慢尤其是约束函数里潮流方程要重新组装一次的时候。实测来看61秒到2分钟是正常的但如果超过5分钟大概率是约束函数写得不高效。排查技巧约束函数里不要每次都重新算节点导纳矩阵Ybus把它设为全局变量或嵌套函数共享变量只算一次。还有就是如果发现某次求解耗时爆炸式增长可以尝试把部分等式约束转化成不等式约束的“松弛”形式比如把“节点注入功率等于潮流计算值”放宽为“注入等于计算值的误差小于1e-4”。这样性能快很多代价是结果精度略降但用于方案预研完全够。一点个人经验我实际做完这个项目的体会是协同优化模型的代码并没有比传统无功优化复杂多少真正的成本在建模细节和调试。数据单位、决策变量索引、初始点这三点控制好你就能从一个“看起来很高大上”的数学模型得到一个能稳定收敛的实用程序。如果准备自己复现我强烈建议先用一个小系统比如5节点或14节点电网加3节点气网把框架跑通再上IEEE30这个规模。这个模型后续还可以继续扩展把单时段变成24小时滚动优化加入储能和氢能设备碳价改成阶梯式碳税都是现成可以往下做的方向。电气互联系统的调度问题跟着碳中和的节奏走后面需要解决的场景会越来越多。