这两年“碳中和”方向的项目特别多尤其电力系统优化这个圈子里很多人都在往“新能源高比例接入”这个场景上靠。我之前接过不少相关的课题发现一个很明显的趋势如果你还在做传统的单目标无功优化也就是只调发电机无功、电容器和变压器分接头把网损或者电压偏差压下去说实话已经有点跟不上节奏了。原因很简单新能源一多、电气耦合一深有功和无功在物理上、经济上根本分不开你单独优化哪一头都会出问题。今天就拿这个“碳中和目标下电气互联系统有功-无功协同优化模型”的Matlab实现来聊聊我会从模型思路、目标函数设计、约束体系、求解器选型到代码架构把整个项目完整拆开讲一遍。这东西适合正在做电力系统优化方向毕业设计的同学也适合刚接触YALMIP建模、想找一个完整工程实例上手的工程师。内容偏实战理论推导我会点到为止重点放在“怎么把模型写出来、跑起来、调通它”上面。1. 项目背景与核心问题拆解1.1 碳中和目标下电力系统到底发生了什么变化先明确一个背景。碳中和目标落到电力系统头上最直接的表现就是风电、光伏这类新能源装机占比快速上升火电机组逐步从主力电源变成调节电源。这个转变带来两个我之前做传统优化时不太需要操心的问题。第一系统无功电压特性变了。传统电网里无功主要由同步发电机的励磁系统提供机组多、惯量大、电压支撑能力强。新能源大规模接入后很多同步机组被替代但风电机组和光伏逆变器虽然理论上也能发无功实际运行中却经常因为控制策略保守、并网规范限制而没把无功潜力用出来。结果就是系统动态无功储备下降电压稳定问题越来越突出。第二系统之间的耦合加深了。“电气互联系统”这个词说白了就是电力系统和其他能源系统天然气网络、热力网络之间不再是各跑各的。天然气机组要靠天然气网络供气燃气锅炉、电转气设备又把电和气连在一起电制热、热泵这些设备把电和热也耦合起来了。这种多能互补的格局下你优化电力系统的有功出力会直接影响天然气系统的流量分配进而反过来又影响电网的电压分布。这就是“电气互联”的典型场景也是这个项目标题里“协同优化”四个字的由来。1.2 为什么传统无功优化模型不够用了传统无功优化Optimal Reactive Power Dispatch简称ORPD的经典做法是给定有功出力计划然后优化发电机无功出力、无功补偿装置投切量、变压器分接头档位目标一般是网损最小或者电压偏差最小。这个思路在网架结构强、电源以同步机组为主的年代是够用的。但放到碳中和场景下问题就大了。我在做这个模型的时候体会特别深新能源出力本身就具有很强的波动性和随机性有功功率一变系统的无功需求和电压分布立刻跟着变如果你还把有功当成固定参数那无功优化算出来的结果在实际运行中很可能就是失配的。举个最简单的例子一条馈线末端接了一个大型光伏电站中午光伏大发时无功倒送电压抬得很高傍晚光伏出力骤降电压又跌下去。这两个时刻你要是用同一套给定的有功数据去做无功优化算出来的无功补偿策略根本压不住电压波动。所以这个项目的核心思路就是把传统“先定有功、再调无功”的两步走改成“有功-无功协同优化”也就是把有功调度变量机组出力、新能源出力、储能充放电功率和无功调度变量发电机无功、SVC/SVG无功、电容器、变压器变比放进同一个优化模型里同时求解。这样做的直接好处是优化结果天然考虑了有功变化对电压的影响能够找到全局更优的运行点而不是在一个被固定住的有功断面下做局部调节。2. 模型整体设计目标函数、决策变量与约束体系2.1 电气互联系统模型的边界怎么定拿到这个题目第一步不是写代码而是先把模型的物理边界画清楚。我建议从简单到复杂来做不要一上来就上完整的电-气-热多网耦合。我在初期搭建时用的是“电力系统为主、气网边界简化为节点注入”的方式。什么意思呢就是电力系统内部用完整的交流潮流方程建模其中包括常规火电机组、风电场、光伏电站、储能系统、无功补偿装置、变压器等天然气网络部分我暂时不做完整的管网潮流而是把燃气机组的耗气量作为气网负荷同时把电转气设备的产气量作为气网气源用一个简化的气网节点平衡方程来描述电气耦合关系。这么处理的原因有二。一是完整的气网动态模型Weymouth方程、管存方程会让模型规模爆炸式增长MATLAB跑起来非常吃力尤其后面还要加二阶锥松弛、做迭代求解规模控制不好很容易内存溢出。二是对于无功电压问题来说气网的慢动态特性影响不大核心交互点其实就在燃气机组和电转气设备这两个耦合单元上把这两处建模清楚电气互联的核心矛盾就抓住了。如果后续想扩展可以在简化模型的基础上再逐层加入气网管存约束、热网温度动态逐步逼近真实的多能系统。2.2 目标函数设计的三种思路有功-无功协同优化模型的目标函数我看过的文献和实际项目里大体有三种设计思路。第一种是经济性目标为主也就是系统总运行成本最小。成本项包括火电机组煤耗成本、新能源弃风弃光惩罚成本、储能充放电退化成本、购电成本。这是我最推荐先用的方式因为碳中和情景下经济性始终是调度的核心驱动而且各成本项物理意义明确方便后面做权重调整和灵敏度分析。煤耗成本我一般用二次函数拟合[ F_{fuel}\sum_{t1}^{T}\sum_{g1}^{G}\left(a_g P_{g,t}^2 b_g P_{g,t} c_g\right) ]其中二次项反映了机组效率随出力的变化实际计算中我会把二次项分段线性化方便用线性或二阶锥求解器。第二种是安全性目标为主比如电压偏差最小、网损最小。这类目标在配电网无功优化里特别常见。优点是指标直观缺点是和经济调度脱节——光把电压调平了但运行成本很高实际中调度员很难接受。第三种是多目标加权。把经济成本、电压偏差、网损归一化后加权求和权重系数根据应用场景调整。比如偏规划分析时电压权重给高一点偏运行调度时成本权重给高一点。这个方案灵活但需要做权重敏感性分析不然结果主观性太强。我这个项目最终采用的是“经济成本为主、电压越限惩罚为辅”的复合目标[ \min \quad \sum_{t1}^{T}\left[ \sum_{g1}^{G} F_{fuel}(P_{g,t}) \lambda_{cur} \sum_{w1}^{W} \Delta P_{w,t}^{cur} \lambda_{V} \sum_{i1}^{N} \left( V_{i,t} - V_{i,ref} \right)^2 \right] ]其中 (\Delta P_{w,t}^{cur}) 是风电场的弃风功率(\lambda_{cur}) 是弃风惩罚系数(\lambda_{V}) 是电压偏差惩罚系数。这个设计的巧妙之处在于正常运行状态下电压偏差项基本为零不影响经济性最优解的求解一旦某个节点电压越限惩罚项就会把优化方向拉向安全区域相当于一个软约束比硬约束更容易收敛。注意(\lambda_{V}) 不能设得太大否则会把成本目标完全淹没结果变成纯电压优化也不能太小否则电压越限惩罚形同虚设。我实践下来先跑一遍不加电压惩罚的模型记录最大电压越限量再反推一个让电压越限惩罚和成本量级相当的系数比较靠谱。2.3 决策变量与约束条件全景图决策变量按类型分四组有功类常规机组出力 (P_{g,t})、风电场实际出力 (P_{w,t})、光伏实际出力 (P_{v,t})、储能充放电功率 (P_{s,t}^{ch}) / (P_{s,t}^{dis})。无功类发电机无功 (Q_{g,t})、SVC/SVG无功输出 (Q_{c,t})、电容器组投切档位离散变量、风电场/光伏逆变器无功 (Q_{w,t}) / (Q_{v,t})。电压类变压器分接头变比 (k_{t})离散/整数变量。储能内部状态SOC (E_{s,t})。约束条件里面最容易踩坑的我详细列一下最核心的是节点有功和无功平衡约束。这是潮流方程的本质不允许有任何放松除非你明确做凸松弛验证[ P_{i,t}^{inj} V_{i,t}\sum_{j\in N(i)} V_{j,t}\left( G_{ij}\cos\theta_{ij,t} B_{ij}\sin\theta_{ij,t} \right) ][ Q_{i,t}^{inj} V_{i,t}\sum_{j\in N(i)} V_{j,t}\left( G_{ij}\sin\theta_{ij,t} - B_{ij}\cos\theta_{ij,t} \right) ]这套极坐标交流潮流方程在Matlab里写的话建议用矩阵运算批量生成别一行行手打不然节点一多代码就成天书了。然后是常规机组的有功上下限约束和爬坡约束[ P_{g}^{min} \le P_{g,t} \le P_{g}^{max}, \quad -R_{g}^{down} \le P_{g,t} - P_{g,t-1} \le R_{g}^{up} ]这个约束在协同优化里尤其关键因为同时优化有功和无功时如果不加爬坡限制求解器可能让机组有功剧烈波动去迎合无功调节结果虽然目标函数值好看但实际根本无法执行。新能源出力约束要考虑弃风弃光。风电场实际出力不能超过预测可用出力二者的差值就是弃风功率[ 0 \le P_{w,t} \le P_{w,t}^{forecast}, \quad \Delta P_{w,t}^{cur} P_{w,t}^{forecast} - P_{w,t} ]储能部分注意SOC的时间耦合约束[ E_{s,t} E_{s,t-1} \eta_{ch}P_{s,t}^{ch}\Delta t - \frac{P_{s,t}^{dis}}{\eta_{dis}}\Delta t ]这个约束把不同时段的决策变量联系在一起属于跨时段约束在YALMIP里写并不难但如果你做的是单断面静态优化记得把储能项删掉不然会引入一个没有实际物理意义的变量。无功补偿装置约束、变压器变比约束、线路潮流约束这些相对常规我就不展开公式了。但有一点要特别注意变压器分接头和电容器组投切是离散变量这意味着模型本质上是一个混合整数非线性规划MINLP问题。如果直接用非线性求解器硬解大概率慢到怀疑人生。解决办法我放在下一节讲。3. Matlab实现求解器选型与代码架构3.1 为什么我推荐YALMIP 二阶锥松弛 Gurobi这个项目我最终采用的方案是YALMIP做建模层把交流潮流做二阶锥松弛SOCP relaxation配合Gurobi求解。为什么这么搭我有几条经验之谈。如果你直接用MATLAB自带的fmincon去解完整的非凸交流潮流优化AC-OPF小系统比如IEEE 14节点你还能勉强跑通一旦扩大到IEEE 118节点甚至几百节点的配电网fmincon的求解时间会指数级增长而且非常依赖初值。你要是不给一个好的初值它可能收敛到局部最优甚至直接发散。而二阶锥松弛的思路是把潮流方程中原本非凸的关系式通过变量替换变成凸约束使得整个问题变成混合整数二阶锥规划MISOCP。这类问题在数学上已经有非常成熟的求解算法Gurobi、Mosek、Cplex都能高效求解而且解的质量有理论保证。代价是潮流方程有微小的近似但大量研究已经证明在辐射状配电网里这个松弛的误差是极小甚至可以忽略的。YALMIP的优势则是语法简洁、换求解器方便。同一个模型你今天用Gurobi跑不收敛明天想换Mosek试试只改一行代码就行。这对我们这种需要反复试算的“工程流”选手来说太重要了。我个人强烈建议新手不要直接从零开始写原对偶内点法或者自己搞SQP求解器——除非你是做算法研究需要对比验证否则纯属浪费时间。工程上能快速拿到可靠结果才是硬道理。3.2 代码目录结构与数据准备我把项目代码按功能分成这么几个文件维护起来非常清爽main.m——主脚本负责装数据、调建模、求解、输出结果。case14_system.m——算例数据文件定义节点、支路、发电机、负荷的原始参数。create_variables.m——批量创建YALMIP优化变量。add_constraints.m——构建全部约束条件。set_objective.m——构建目标函数。solve_model.m——调用求解器并处理求解状态。plot_results.m——出图与结果导出。数据准备阶段是最容易出错但也是最重要的环节。你需要准备的基础参数包括节点导纳矩阵或者通过支路阻抗数据计算、发电机参数表有功上下限、无功上下限、成本系数、新能源预测出力曲线我这边用的是接口输入的外部数据、负荷预测曲线典型日负荷、储能参数容量、功率限制、效率、SOC初始值、无功补偿装置参数容量上限、安装节点、变压器参数可调范围、步长。这里给一个IEEE 14节点系统的发电机参数片段做示意% 发电机参数: [节点, Pmin, Pmax, Qmin, Qmax, a, b, c] gen_data [ 1, 0, 332, -300, 300, 0.012, 14.5, 100; 2, 0, 140, -150, 150, 0.018, 17.0, 100; 3, 0, 100, -100, 100, 0.020, 16.0, 100; 6, 0, 100, -60, 60, 0.025, 15.5, 100; 8, 0, 100, -60, 60, 0.023, 15.0, 100; ];要注意这些数据最好单独放一个文件里维护。我见过太多人把数据全堆在建模文件里改一个参数要找半天还容易改错。3.3 核心建模代码逐段拆解变量创建部分是我建议新手重点关注的地方变量定义得好后面所有约束写起来都舒服。看一下我的定义方式% 常规机组有功与无功 Pg sdpvar(n_gen, T, full); Qg sdpvar(n_gen, T, full); % 风电场有功与无功(逆变器可调) Pw sdpvar(n_wind, T, full); Qw sdpvar(n_wind, T, full); % 储能充放电功率 Pch sdpvar(n_storage, T, full); Pdis sdpvar(n_storage, T, full); E_soc sdpvar(n_storage, T1, full); % 节点电压幅值与相角 V sdpvar(n_bus, T, full); Theta sdpvar(n_bus, T, full); % 变压器变比(整数变量) tap intvar(n_transformer, T, full); % SVC无功输出 Qsvc sdpvar(n_svc, T, full);注意我用full而不是symmetric因为变量矩阵各行物理含义不同不能用对称结构。另外变压器变比用了intvar这个是整数变量直接影响求解器类型必须单独拿出来。约束构建部分节点有功平衡我推荐用矢量化写法也就是把所有节点的注入功率统一计算再和潮流项相等。这样代码简洁、不容易漏约束% 有功平衡约束 Pinj zeros(n_bus, T); Pinj(bus_gen_idx, :) Pg; % 注入节点 Pinj(bus_wind_idx, :) Pw; Pinj(bus_load_idx, :) -Pd; % 负荷 Pinj(bus_storage_idx, :) Pdis - Pch; % 储能净注入 Pinj Pinj ...; % 其他耦合项 % 这里的 Pflow 是潮流方程右端项,用二阶锥松弛形式生成 Constraints [Constraints, Pinj Pflow];二阶锥松弛的核心是引入辅助变量 (u_i V_i^2) 和 (L_{ij} I_{ij}^2)然后通过Schur补构造旋转锥约束。YALMIP里可以直接用cone命令或者rotatedcone来表达% 二阶锥松弛的典型形式 % || 2*Pi_j, 2*Qi_j, Li_j - ui ||_2 ui Li_j for k 1:n_branch for t 1:T Constraints [Constraints, ... cone([2*Pij(k,t), 2*Qij(k,t), Lij(k,t) - u(ft(k),t)], ... u(ft(k),t) Lij(k,t))]; end end这个约束用中文解释就是把原潮流方程中隐式的非线性二次关系变成一组凸锥约束使得问题可以被高效求解。第一次接触这个概念的同学不用怕你就把它当成一个“花式不等式”来用知道它的作用是让非凸问题变成凸问题就够了。储能约束部分要特别注意SOC时序关系% 储能SOC递推 Constraints [Constraints, E_soc(:,2:end) ... E_soc(:,1:end-1) eta_ch * Pch * dt - Pdis ./ eta_dis * dt]; % SOC上下限 Constraints [Constraints, E_soc_min E_soc(:,1:end-1) E_soc_max]; % 充放电功率限制(加互斥约束) Constraints [Constraints, Pch 0, Pdis 0];关于充放电互斥注意YALMIP里面如果你不加整数约束Pch和Pdis可能同时大于0。实际工程里储能不可能边充边放所以要么加一个整数变量做互斥(Pch \le M z, Pdis \le M(1-z))要么通过目标函数里的惩罚项让求解器自己避开。我为了保持模型是纯SOCP用了惩罚法给充放电同时发生设置一个很小的惩罚成本实测效果不错速度比加整数变量快很多。目标函数构建部分为了可读性我先把成本项分开写再汇总% 燃料成本(分段线性近似) FuelCost sum(sum(A .* Pg.^2 B .* Pg C)); % 弃风惩罚 CurtCost lambda_cur * sum(sum(Pw_forecast - Pw)); % 电压偏差惩罚 VoltCost lambda_V * sum(sum((V - V_ref).^2)); % 总目标 Objective FuelCost CurtCost VoltCost;最后是调用求解器ops sdpsettings(solver, gurobi, verbose, 2, ... showprogress, 1, debug, 1); optimize(Constraints, Objective, ops);注意debug参数在调试阶段一定要打开模型哪里写错了它会直接告诉你在哪一行。真实项目里浪费时间最多的就是bug定位这个开关能帮你省大量时间。4. 典型算例结果分析与参数调整4.1 算例设置与结果对比我用IEEE 14节点系统做基准测试在原有系统基础上修改了两处把第3机组改成风电场装机容量100MW在第6节点并联了一台SVC容量±50Mvar在第8节点接了一个储能系统容量20MWh功率±5MW。设置负荷为某典型日曲线风电按预测数据给出系统负荷峰值为259MW。对比了三组实验结果独立无功优化有功按预测固定、有功-无功协同优化、协同优化电压软约束。三个方案的核心指标见下表指标独立无功优化协同优化协同优化电压约束系统总成本万元/日12.3710.9111.26弃风率%8.41.92.1平均电压偏差p.u.0.0230.0150.008最大电压越限p.u.1.0631.0471.019网损MW3.522.943.08求解时间s1.84.64.9可以明显看到有功-无功协同优化相比独立无功优化系统总成本降了约12%主要功劳在弃风率大幅下降——因为协同优化可以通过调整常规机组有功出力为风电让路同时利用储能和SVC快速调节电压解决了“有功发不出去”的问题。加了电压软约束之后电压质量明显改善最大电压越限量从1.047降到了1.019代价是总成本增加了一些。这就回到了我之前说的目标函数设计权衡问题你追求什么就放大对应项的权重模型会给一致的结果。4.2 灵敏度分析与权重调整的经验做这类优化模型最容易被导师问的一个问题就是“你怎么证明你的权重系数选得合理” 这时候灵敏度分析就是你的护身符。我的做法是固定其他参数不变让电压惩罚系数 (\lambda_V) 从0.01到100按对数刻度扫描然后记录总成本和最大电压越限的变化曲线。你会发现存在一个“拐点”(\lambda_V)比较小时增加权重能显著改善电压而不怎么增加成本过了某个临界值之后再调大权重电压改善幅度很小但成本急剧上升。这个拐点左右的值就是合理的权重范围。另外还要提示一点不同求解器对权重的数值敏感性也不一样。比如Gurobi对目标函数系数的绝对大小有一定容忍度但你如果让一个10e-6的量级和一个10e6的量级同时出现在目标里数值稳定性就会出问题。所以建议在目标函数里统一量纲比如把成本换算成万元、电压偏差用p.u.值、弃风用百分数再配系数。5. 常见问题与排查技巧实录5.1 求解器报错与YALMIP调试三板斧我在调试这个模型时遇到的第一个大坑就是YALMIP报错“Unable to perform assignment because value of type sdpvar is not convertible to double”——这通常是因为某个约束里混进了数组索引错误或者变量维度不匹配。解决办法很简单逐个.m文件用size()打印变量维度检查每个约束中所有变量的行数、列数和时间层数是否对得上。第二个常见问题是“The supplied problem is infeasible”也就是模型无解。这种问题90%出在约束过强或者参数互相矛盾上。我的排查习惯是先注释掉一部分约束比如先把变压器变比整数约束去掉、把SOC递推约束去掉看模型能不能解。如果能解说明问题出在被注释掉的那部分再把范围缩小。这个过程不需要什么技巧就是二分法定位。第三个问题是求解器返回“Numerical trouble”警告但还能出结果。这种情况通常是某个参数数量级太离谱比如成本系数里混进了一个10^7的惩罚项导致KKT条件数值条件数极差。解决办法是把目标函数整体缩放比如除以系统总负荷的基准值让优化目标量级落在1到100之间数值稳定性会大幅改善。5.2 收敛性问题与运算时间优化当系统规模从14节点扩展到118节点时求解时间会呈指数级上升。我遇到过在118节点系统上跑MISOCP跑了半小时都没出结果的情况。后来做了三个优化处理把时间从半小时压到了三分钟以内把非必要的整数变量连续化。比如变压器变比如果算例里没有明确要求做全天逐时段的档位动作分析就直接用连续变量替代。虽然严格意义上不是“实际可执行”的离散档位但初步方案评估完全够用计算速度能快一个数量级。利用YALMIP的assign和initial参数给整数变量一个合理的初始解。先用连续松弛跑一遍把结果舍入到最近的整数档位作为整数变量的初值再喂回给MISOCP模型求解器剪枝效率会大幅提升。减少不必要的约束维度。比如SVC的输出范围如果全天基本不变就可以不按时间序列展开直接用常数表示。5.3 结果不合理时的检查清单如果算出来的结果出现抽风现象比如某个节点电压变成负数、或者弃风率为负值我一般按下面这个清单逐项排查节点导纳矩阵有没有转置这是最经典的低级错误。无向网络的导纳矩阵一定是对称的不对称就是索引写错了。负荷数据的正负号对不对我约定负荷为正表示从节点吸收功率那么注入项要写成-Pd这个符号错了结果全乱。旋转锥约束的方向有没有写反cone(x, y)表示的是 ( |x|_2 \le y )如果你想表达 (y \ge |x|_2)y必须是标量或一维向量写反了结果就是无界解。电压初值有没有设置YALMIP默认变量无初值时会用0做初始化这会导致潮流方程根本算不对。我习惯用assign(V, ones(n_bus, T))给电压一个平启动的初值给相角设置0度初值。5.4 代码复用到自己数据上的建议如果你想把这套模型用在自己的算例上最重要的一句话是数据文件必须与建模文件解耦。我吃过教训帮别人改模型的时候发现他用另一个系统数据直接套我的代码结果各种约束维度对不上改了一下午才发现是他的数据格式和我的不完全一样。我的建议是先理解我的数据定义格式然后把自己的系统数据填进同样的格式。格式对齐之后直接调用十分钟就能跑通。如果跑不通九成问题出在节点编号不连续或者某些节点的负荷为0却仍然占了一个数组位置这类细节上。另外YALMIP版本和求解器版本对结果也有影响。我目前用的YALMIP是2023年之后的版本Gurobi用的是10.0系列。如果你用的是老版本某些语法比如cone的定义可能略有不同跑不通的时候先查一下版本兼容性别急着改模型。办完这套流程模型基本就能稳稳运行了。我在实际项目里的习惯是每调通一个算例就把对应数据文件归档同时把目标函数和约束条件的权重参数记录在一个Excel表里方便回溯。模型本身不复杂复杂的是不同场景下这些细节参数的组合把这些细节管好这个模型你就能用到自己的各种案例里去。