企业官网建设流程全解析

引言

这几年的电力系统优化研究，尤其是电气互联系统（电网+天然气网耦合）方向，几乎人人都离不开“碳中和”三个字。我去年在做一个区域综合能源系统的优化调度项目时，发现多数现有模型要么只做有功优化、要么把无功当作固定潮流结果，根本没有做到真正意义上的有功-无功协同。但实际工程里，电压越限、网损偏高、无功补偿设备利用率低这些问题，恰恰是影响系统经济性和安全性的重要因素。

所以这次我把“碳中和”目标下的电气互联系统有功-无功协同优化模型完整做了一遍，用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设备与电网耦合。

关键参数包括：

负荷和新能源出力曲线我生成了一整天的时序数据，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 end

3.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，然后配一个商业求解器。

实际对比测试中：

这个项目里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_SVC+1: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 求解阶段的坑

坑4：SOCP松弛在环网中不紧。

前面提过，如果电网是环网结构，SOCP松弛可能会产生一个理论上最优但物理上不可行的解。判断方法很简单：检查支路电流约束的对偶乘子或者直接计算有功损耗的实际值，如果和优化结果偏差超过某阈值（比如5%），就需要改用精确非线性求解或者加割平面约束。对于辐射状配电网，这个风险基本不存在。

坑5：求解气体的Weymouth方程时不收敛。

我一开始是用内置的fmincon处理Weymouth方程非线性，结果24时段联合求解时经常卡住。后来全部改成YALMIP+LMI线性化方式，全部工况一次收敛。结论：电力系统优化尽量用数学规划框架，不要轻易用通用的非线性求解器。

6.3 代码级独家技巧

最后分享一个比较实用的小技巧：在做多场景对比时，一定要把value()提取之后的结果保存好，不然Gurobi的模型对象在重新optimize的时候会覆盖结果。我在DEBUG时遇到过多次提取结果为空的情况，后来养成了“solve一步、立即value并存储”的习惯。

另外，YALMIP的assign函数在热启动时非常有用，但要注意变量的维度必须完全一致，否则YALMIP不会报错，只会默默忽略初始化值。验证是否初始化成功的方法是在optimize前检查solvesdp('warmstart', 1)是否生效。

实测下来，这套模型和代码流程跑通之后，后续换算例、换参数都很方便。只要数据结构定义一致，换一个33节点配电网就是改data_ies.m文件的事，不需要动建模代码。这也是我这个项目最满意的部分——模型和数据的解耦做得比较干净，给后面扩展更大的系统留了余地。

我个人在实际操作中的一个体会是：无功优化这块内容，很多做综合能源系统的人容易忽略，总觉得“先把有功跑对就行”。但这个项目做完之后我很确定，如果系统里新能源渗透率超过20%，不做无功协同的模型基本撑不住电压约束。哪怕只是为了在论文里加一个“电压越限对比”的图表，也值得把无功这块完整建模进去。这个模型后续还可以加储能、需求响应、碳捕集设备做扩展，扩展接口我都预留好了。