最近有个学弟拿着一个标题来找我:“硕士论文复现,计及P2G厂站的电-气综合能源系统规划研究,附Matlab代码。”他想让我帮忙看看这个课题怎么入手。说实话,这类标题在电气工程和能源系统方向的毕业论文里相当常见:P2G(Power-to-Gas)厂站作为电力系统和天然气系统的耦合环节,在规划阶段同时考虑投资决策和运行模拟,最后用Matlab把模型跑通。这篇文章我就把整个复现过程中的思路、模型、代码结构和踩过的坑一起整理出来,给同样要做论文复现的人一条可以直接参考的路径。
先说适用范围:如果你手头也有类似方向的硕士论文需要复现,或者正在做电-气综合能源系统规划选题,又或者单纯想搞懂Yalmip怎么搭MILP模型,这篇文章应该能帮你省下不少时间。我默认你至少会用Matlab的基本语法,但不要求你已经有商业求解器,因为后面我会专门说明如何配置求解环境。
1. 复现之前:先搞懂这套系统到底在做什么
1.1 论文复现不是“翻译代码”,而是还原思路
很多同学拿到一篇硕士论文,第一反应是找“源代码”或者“现成包”。实际上,绝大多数毕业论文并不会附完整代码,标题里出现“附Matlab代码”往往是作者自己在复现时写的版本,或者培训机构打包出来的材料。真正的复现工作,应该是把论文里的数学模型用代码重新表达一遍,并让它在标准算例上得出可解释的结果。
我的习惯是分四步走:第一步,通读摘要和结论,弄清这篇论文解决了什么问题,创新点到底是P2G厂站建模、规划模型改进,还是求解算法改良。第二步,看核心章节的数学模型,把目标函数、约束条件一条条摘出来,理清变量下标是节点、时段还是场景。第三步,根据模型找算例数据,硕士论文通常会用IEEE节点系统或者某个区域系统的简化数据,自己编数据时也要和论文口径尽量一致。第四步才是动手写代码,写的时候不是逐行翻译公式,而是要让每个约束都能在代码里找到对应的可计算表达式。
这一篇要复现的论文,核心就是“电-气综合能源系统规划”。所谓规划,就是在已知负荷、风电出力等条件下,回答“P2G厂站应该建在哪、建多大容量、电网气网需不需要扩容、系统总成本最低是多少”这几个问题。理解了这个问题,代码的逻辑就清晰了。
1.2 电-气综合能源系统和P2G厂站的定位
先把这个系统拆开看。电力系统一侧,有常规机组、风电场、电网线路和电负荷;天然气系统一侧,有气源、输气管道、压缩机和气负荷。传统规划是把这两个系统分开做的,互不影响。但P2G厂站一出现,电网和气网就被“挂钩”了:P2G消耗电能,通过电解水制氢,再把氢进一步甲烷化后注入天然气管道,或者直接以氢能形式供给其他用户。
这样一来,原本可能被弃掉的风电,可以转化成天然气被储存利用,电气系统之间的能量流动变成一个闭环。这也是论文标题里“计及P2G厂站”的意义所在:规划时如果不考虑这个耦合环节,就可能低估风电消纳能力,也会错过跨系统优化带来的成本下降空间。复现的时候,你要时刻记住,P2G不是一个简单的负荷,也不是简单的气源,而是一个可以双向调节的柔性耦合设备。它在数学上表现为“消耗多少电功率,就相应产生多少天然气流量”,两者之间靠效率和运行约束连接。
我复现时用的算例规模不大:电网取33节点配电网,气网取20节点的简化输气网,两个系统通过若干个候选P2G节点耦合。规模小不是问题,关键是模型要完整,求解器能算出结果,曲线趋势能解释得通,这就达到复现的目的了。
2. 核心建模:P2G厂站怎么“吃电产气”
2.1 P2G过程拆分:从电解水到燃气
P2G的技术路线主要有两种:一种是电解水制氢,氢气直接利用或储存;另一种是电解水后再进行甲烷化反应,把氢气和二氧化碳合成为甲烷,这样可以完全兼容现有天然气管道和用气设备。论文里如果想突出“注入气网”这个环节,大多采用甲烷化路线。
从能量转换角度看,P2G整体的效率通常在55%到75%之间。简化建模时,不需要把电解槽、甲烷化反应器分别建模,只需要定义:
[ G_{P2G}(t) = \eta_{P2G} \cdot P_{P2G}(t) ]
其中,(P_{P2G}(t)) 是第t个时段P2G消耗的电功率,(G_{P2G}(t)) 是注入气网的天然气流量(通常换算成热值或等效功率),(\eta_{P2G}) 是综合转换效率。这里有个容易踩的坑:如果不把单位统一,算出来的结果会非常离谱。比如电功率用MW,天然气流量用立方米每小时,就需要乘上天然气的热值系数再转换为MW或MWh。
P2G厂站还要考虑运行范围。实际设备不会在零到额定容量之间任意运行,通常会有一个最小技术出力比例,比如额定值的10%或20%。此外,升降负荷速率也有限制。因此在模型里,除了容量约束,还要加上爬坡约束和启停状态变量。不过,如果论文主要研究规划,运行部分往往只做典型日简化,爬坡约束可以按整个时段步长设置,不需要太复杂。
2.2 耦合约束与能量平衡
P2G在系统里同时出现在两个网络中:电能侧,它是一个“电负荷”;天然气侧,它是一个“气源”。所以建模时,需要把P2G的功率变量同时写到电网节点功率平衡方程和气网节点流量平衡方程里。
电网节点功率平衡是经典形式:
[ \sum_{g \in i} P_{g,t} + \sum_{w \in i} P_{w,t} + P_{line,t}^{in} - P_{line,t}^{out} - P_{P2G,i,t} = L_{i,t} ]
天然气网络的节点平衡则要复杂一些。简单起见,很多硕士论文把天然气管网看成稳态模型,只考虑管道流量约束,忽略储气动态。那么对每个气网节点有:
[ \sum_{s \in i} S_{s,t} + \sum_{p2g \in i} G_{P2G,i,t} - \sum_{l \in i} G_{load,l,t} = \sum_{k} Q_{k,t}^{out} - \sum_{k} Q_{k,t}^{in} ]
天然气管道的流量由管道两端压力决定,通常用Weymouth方程描述,这是一个非线性约束。实际代码里必须把它线性化或松弛化,否则Matlab里的求解器很难直接处理。关于这一点,我在第4章会专门展开。
2.3 规划模型的目标函数与决策变量
规划模型的目标函数一般是最小化总成本。总成本分成两大部分:投资成本与运行成本。投资成本是待建P2G厂站的建设费用,可能还要包括电网扩容、气管网扩容的费用;运行成本包括机组燃料成本、购气成本、弃风惩罚等。
用公式表达就是:
[ \min \ C^{inv} + \sum_{t \in T} \left( C^{gen}_t + C^{gas}_t + C^{curtail}_t \right) ]
决策变量包括两类。投资决策变量通常是0-1变量,比如某个候选节点是否建设P2G;连续变量则是P2G建设容量。运行决策变量各时段都有:机组出力、气源产气量、P2G消耗电功率、P2G产气流量、管道流量、节点相角或压力等。
这里有个关键点需要提醒:规划模型里的投资变量和运行变量是耦合的。比如P2G建设容量是(C_{p2g}),运行功率就不能超过这个容量,还要受“是否建设”这个0-1变量的限制。写成约束就是:
[ 0 \le P_{P2G,i,t} \le C_{p2g,i} \cdot x_{p2g,i} ]
如果(x_{p2g,i}=0),那么这个节点就不能消耗电功率;如果为1,容量上限就是(C_{p2g,i})。这类约束是混合整数规划的基本写法,代码里必须保证所有运行变量都通过类似方式与投资变量挂钩,否则就会出现“没建厂站却在运行”的荒谬结果。
3. Matlab代码实现与复现步骤
3.1 整体代码框架与文件目录
我复现时用的是Matlab + Yalmip + Gurobi这套组合。Yalmip是一个建模工具箱,它最大的好处是不需要自己手写大规模稀疏矩阵,可以直接用符号变量写约束,然后调任意底层求解器。Gurobi负责求解混合整数线性规划(MILP),性能远好于Matlab自带的intlinprog,尤其是节点规模稍大时差距明显。
代码建议按下面的目录组织:
main.m % 主程序,跑整个流程 load_data.m % 读取系统参数和负荷数据 build_network.m % 构建电网/气网节点导纳矩阵、管道参数 define_variables.m % 定义所有决策变量 constraints.m % 添加所有约束 objective.m % 添加目标函数 solve_and_post.m % 求解、提取结果、画图这样拆的好处是,论文里每修改一个假设,只需改动对应的函数。比如想改P2G效率,不用在几百行代码里找参数,直接改load_data里的eta_p2g。
主程序的运行逻辑很简单:
clear; clc; close all; load_data; build_network; define_variables; constraints; objective; solve_and_post;当然,实际过程中我不会把所有变量都塞到全局空间,但作为复现脚本,这种线性流程最容易调试。如果想做成更正式的项目,可以用结构体把数据包起来,比如data.bus、data.branch、data.gas,避免命名冲突。
3.2 用Yalmip搭建优化模型的关键代码段
这里我直接贴一段简化版的建模代码,对应第2章里的P2G耦合约束:
% 候选P2G节点数量 n_bus = 33; T = 24; % 投资变量 C_p2g = sdpvar(n_bus,1); % P2G建设容量 x_p2g = binvar(n_bus,1); % 是否建设 % 运行变量 P_p2g = sdpvar(n_bus,T); % P2G消耗电功率 G_p2g = sdpvar(n_bus,T); % P2G注入气网功率 % 容量上限约束 Constraints = [Constraints, 0 <= C_p2g <= C_max .* x_p2g]; Constraints = [Constraints, sum(C_p2g) <= C_total_max]; % P2G运行约束:功率不超过容量 Constraints = [Constraints, 0 <= P_p2g <= C_p2g * ones(1,T)]; % P2G能量转换 Constraints = [Constraints, G_p2g == eta_p2g * P_p2g]; % 电网节点平衡约束(简化) for i = 1:n_bus Constraints = [Constraints, ... sum(P_gen(i,:),2) + ... % 机组出力 P_wind(i,:) + ... P_line_in(i,:) - P_line_out(i,:) - P_p2g(i,:) == P_load(i,:)]; end有几个地方需要特别说明。
C_p2g * ones(1,T)这一步很关键。C_p2g是33×1的列向量,P_p2g是33×24的矩阵,直接写P_p2g <= C_p2g会报维度错误。乘上ones(1,T)是为了把列向量扩展成矩阵,让每列都共享同一组容量上限。
binvar定义了0-1变量,这是Yalmip的语法。如果你用Gurobi,Yalmip会自动转换成MIP格式,不需要自己声明整数类型。sdpvar则用来定义连续变量,虽然名字里带“SDP”,其实主要用于线性、二次、二阶锥等各种问题。
目标函数可以按这样加:
% 投资成本:单位容量投资 * 容量 inv_cost = cost_per_MW * sum(C_p2g); % 运行成本:逐时段求和 op_cost = sum(sum( ... c_gen .* P_gen + ... % 机组燃料成本 c_gas .* G_source + ... % 购气成本 c_curtail .* P_wind_curtail)); % 弃风惩罚 Objective = inv_cost + op_cost; % 求解 optimize(Constraints, Objective, sdpsettings('solver','gurobi','verbose',1));求解完毕后,用value()提取变量结果:
C_p2g_opt = value(C_p2g); x_p2g_opt = value(x_p2g); P_p2g_opt = value(P_p2g);只要Yalmip安装正确,求解器配置没问题,这一段代码就能跑通。当然,真实论文里约束远不止这些,还会有机组出力上下限、爬坡约束、气源容量约束、节点压力上下限等,但建模思路是完全一样的。
3.3 算例数据准备与参数设置
没有原始数据是复现时最大的痛点。硕士论文通常只给部分参数,我复现时以IEEE 33节点配电网为基础,气网参照比利时20节点天然气管网做简化,再通过5个候选节点把两个网络耦合起来。
常用参数可以按下面这个表设置:
| 参数名称 | 数值 | 说明 |
|---|---|---|
| P2G综合效率 | 0.65 | 电解水+甲烷化整体效率 |
| P2G单位投资成本 | 5000元/kW | 不同论文差异较大,按目标年份调整 |
| P2G最大单点容量 | 10 MW | 候选点容量上限 |
| 风电装机 | 28 MW | 分布在多个节点 |
| 电负荷峰值 | 60 MW | 典型日负荷曲线缩放 |
| 气负荷峰值 | 30 MW | 折算为等效功率 |
| 天然气热值 | 36 MJ/m³ | 用于单位转换 |
负荷曲线我一般取24个时段的典型日数据,分成春夏秋冬四个代表日,模型里用多场景方式处理。多场景会显著增加变量数量,但在单台电脑上求解几个典型日不算难事。
有个单位换算的经验:如果电力侧用MW,天然气侧也用MW(即把流量乘以热值转换为功率),那么P2G的效率就可以直接用无量纲系数。这是最稳妥的做法,可以避免在后面结果分析时被一堆系数绕晕。等模型跑通后,再按论文要求把天然气流量换算回立方米每小时用于画图。
3.4 结果后处理:把指标变成图表
跑出最优解只是第一步,论文里需要看的是结果图表。我的后处理脚本会一次性画出下面几种图:
- 系统拓扑图:在电网单线图上标出P2G建设位置和容量,用不同颜色表示容量大小。
- P2G各时段运行功率曲线:可以看到风电出力大的夜间,P2G消耗功率明显上升;风电出力小的时段,P2G基本停机。
- 弃风率对比柱状图:不装P2G和安装P2G后的弃风量对比,这是体现P2G价值最直观的图。
- 气网节点压力分布:检查天然气管道压力是否越限。
画图代码我习惯用figure+subplot组合,先画出来再统一调整样式。导出图片时用exportgraphics(gcf, 'result.png', 'Resolution', 300),这样论文插图直接够用。如果遇到新版Matlab的exportgraphics在某些旧版本不可用,也可以用print(gcf, '-dpng', '-r300', 'result.png')。
还有个细节:做结果对比时,最好把“不含P2G”和“含P2G”两种场景都跑一遍。很多论文的价值就体现在这个对比里,比如总成本下降多少、弃风率降低多少、P2G建在哪几个节点。复现时如果只跑一个场景,很容易漏掉这个关键结论。
4. 复现过程中最常见的六类坑(附排查方法)
4.1 求解器安装与许可证问题
Yalmip本身只需要下载解压,然后把文件夹加入Matlab路径。真正麻烦的是底层求解器。Gurobi和Cplex都提供学术许可证,用学校邮箱申请很方便,安装时注意版本要和Matlab系统兼容。安装完之后,在Matlab里运行yalmiptest,看到输出里Gurobi显示available就说明配置成功。
如果拿不到商业求解器,可以先用Matlab自带的intlinprog试试。在Yalmip里只需要把solver改成intlinprog:
optimize(Constraints, Objective, sdpsettings('solver','intlinprog'));不过intlinprog对大规模MILP问题会明显吃力,特别是有上千个0-1变量时,求解时间可能从几分钟变成几个小时。我的建议是:初学阶段用intlinprog验证模型正确性,正式跑算例时再切回Gurobi。
这里必须多说一句:网上流传的各种所谓“离线包”和“密钥文件”不建议碰,一方面有法律风险,另一方面容易带恶意脚本。学校能提供学术许可就要用学术许可,没有也没关系,换开源求解器SCIP也完全可以跑通论文算例。
4.2 维度不匹配和稀疏矩阵构造错误
这是新手最容易卡住的地方,也是我帮人调试时见到最多的问题。Yalmip虽然比手写矩阵友好,但变量维度不匹配照样会报错或者产生错误模型。特别是C_p2g * ones(1,T)这类扩展写法,稍不留神就会变成“隐式扩大约束”的错误逻辑。
排查维度问题有几个技巧。第一,定义变量之后就打印size(),确认每个变量的行列数。第二,写约束时尽量保持同一个物理量使用同一维度,比如所有节点变量用n_bus × T,所有时段变量用1 × T。第三,遇到Yalmip报“Unable to perform assignment because size of left side is X and right side is Y”时,不要急着堆repmat,先想清楚这个约束数学上到底是逐点约束还是矩阵约束。
有时候模型不报错但结果异常,也可能是维度扩展写错了。比如我想让每个节点的P2G容量不超过该节点上限,写成P_p2g <= C_p2g就不会触发维度错误,因为Yalmip会把列向量和矩阵做广播运算,但这个广播不一定是你要的。最稳妥的写法是明确扩展成P_p2g <= repmat(C_p2g, 1, T),肉眼一看就明白。
4.3 管道非线性约束的处理不当
天然气管道流量与节点压力的关系是论文模型里最大的坑。Weymouth方程是:
[ Q_{ij}^2 = K_{ij}^2 (p_i^2 - p_j^2) ]
这个约束里有平方项,直接放到MILP模型里是没法求解的。常见的处理方式有两种。
第一种是增量分段线性化。把管道流量和节点压力差关系拆成多段直线,用一组连续变量和二进制变量表示强制落在某一段。这个方法精度高,但变量数量会随分段数增加。
第二种是二阶锥松弛。将(p_i^2 - p_j^2)替换成中间变量,并把等式写成不等式[ Q_{ij}^2 \le K_{ij}^2 (p_i^2 - p_j^2) ]的形式,这样问题就变成混合整数二阶锥规划(MISOCP),Yalmip可以直接用optimize求解,Gurobi从9.0开始也原生支持二阶锥约束,不需要额外处理。
复现论文时,我建议先看原文用的是什么方法。如果原文没说,就先用分段线性化,因为它在MILP框架内实现起来更直觉,后处理也容易画图。如果节点数很多导致计算太慢,再换成二阶锥松弛,求解时间通常能降一个量级。
4.4 MILP求解太慢、收敛性差
规划模型里如果候选P2G节点有10个,每个节点有0-1变量,再加上机组启停变量,MILP规模很容易膨胀。Gurobi求解器默认的MIPGap是1e-4,对论文复现来说没有必要这么严格。可以在求解设置里放宽一点:
op = sdpsettings('solver','gurobi','gurobi.MIPGap',0.01); optimize(Constraints, Objective, op);百分之1的间隙对规划结果影响不大,但求解时间可能从半小时降到两分钟。另外,给所有变量设置合理的上下界也很重要。Yalmip默认变量范围是正负无穷,这会让分支定界过程搜索空间巨大。即使模型里没有显式约束,也应该给关键变量加上边界,比如:
C_p2g = sdpvar(n_bus,1); Constraints = [Constraints, C_p2g >= 0, C_p2g <= 50];设置初始可行解也很有帮助。可以先固定投资变量为0(即不建P2G),求解一次得到运行成本,再把投资变量设为1,得到一个粗略的可行解,然后用Yalmip的assign赋值给变量,再调用optimize,求解器会用这个初始点开始搜索,收敛会快很多。
4.5 结果数值不合理但代码能跑
这种情况最让人头大。代码没报错,求解状态是“solved”,但结果明显不对劲,比如P2G建设容量极小,弃风率反而更高,或者气网流量为负。
我总结下来,最常见的原因是单位不一致。电源侧用kW,负荷侧用MW,天然气侧再用m³/h,这些单位混在一起,模型还能解,但解读全乱了。建议整个项目统一用标幺值,或统一用MW和MWh。如果论文给了基础功率,就在load_data里先把所有数据折算到同一基准。
第二个原因是热值系数错了。天然气的热值按36 MJ/m³算,1 m³/h约等于0.01 MW,如果漏乘这个系数,P2G产气量就会被低估或高估一个数量级。检查办法很简单:单独设置一个只有一台P2G、无其他约束的小测试模型,输入1 MW电功率,看输出是不是0.65 MW等热值,如果不是,就说明单位换算出错了。
第三是目标函数中某个成本项权重过大,导致求解器通过降低这项成本来“优化”。比如弃风惩罚设得特别高,模型可能倾向于建设极贵的储能或P2G来消除弃风,结果总成本反而更高。看到这类结果时,要把目标函数拆分打印出来,看投资成本、燃料成本、惩罚成本各自占比多少,问题往往一目了然。
4.6 版权、引用与代码分享规范
复现论文不是为了抄袭,而是为了把方法跑通并验证可用性。如果你准备把复现代码放到GitHub或者自己的博客,一定要在README里注明原始论文标题、作者、年份和DOI,同时写明这份代码是基于论文模型自己的实现。如果参考了别人的开源代码,还必须遵守对应的开源协议,比如MIT、GPL等。
Matlab代码中如果引用了第三方工具箱,也要注意许可证兼容问题。Yalmip是BSD协议,可以放心用;Gurobi虽然免费给学术使用,但开源项目分发时不能捆绑Gurobi的安装包,只能让用户自己申请。这些看起来都是小事,但真到分享阶段都是必须处理的雷区。
5. 从复现到迁移:还能怎么扩展这套代码
5.1 加入储氢罐,打破“即产即用”假设
基础的P2G模型默认产气后立刻注入气网,不允许存储。实际系统中加一个储氢罐可以显著提升灵活性:风电大发时,可以多产氢存起来,等气价高或者气负荷高峰时再释放。
这段代码的改动并不复杂,只需要增加一个状态变量表示储氢量:
E_h2 = sdpvar(1,T); % 储氢罐能量状态 Constraints = [Constraints, E_h2(:,1) == E_h2_init]; Constraints = [Constraints, E_h2(:,t+1) == E_h2(:,t) + ... G_p2g_partial(:,t) - G_release(:,t)]; Constraints = [Constraints, 0 <= E_h2 <= E_h2_max];有了储氢环节,P2G就不需要严格满足“产气量=注入气网量”,而是可以用额外的变量表示氢气流向储罐或燃料电池/燃气轮机。这类改动适合作为论文第4章的扩展场景。
5.2 改成多目标规划或考虑碳交易
原始论文如果只做单目标成本最小,你可以把碳排放量作为第二个目标。最常用的方法是epsilon约束法:把碳排放设成一个约束,比如总碳排放不超过某个阈值,然后观察总成本如何随阈值变化,画出帕累托前沿。
Matlab里实现epsilon约束法很方便。外层循环用for epsilon = [0.9, 0.8, ...],内层在约束中加入total_emission <= epsilon * emission_base,依次求解,把成本记录到一个数组里即可。这个结果放到论文里可以写:“随着碳排放约束收紧,系统总成本上升至xxx,P2G配置容量增加,说明P2G在低碳转型中起关键作用。”逻辑很顺。
如果论文涉及碳交易机制,也可以在目标函数中加入碳价乘以碳排放量的项。这样P2G的价值就能直接反映在成本上,比单纯看弃风率更有说服力。
5.3 从气网平移到热网/电热耦合
P2G的思路稍作修改就是P2H(Power-to-Heat)。电转热设备比如电锅炉、热泵,与P2G一样都是消耗电能、产出另一种能量。区别在于热网通常不需要Weymouth方程,而是用热力管道传输延迟和温度混合方程建模。如果你能跑通电气系统,再换成热网时只需要把网络约束替换成热网节点功率平衡,模型框架不用变。
这类“换汤不换药”的扩展最适合在毕设里做不同场景对比:同一个规划模型,分别考虑P2G、P2H、P2G+P2H,看哪种技术路线经济性最好、对风电消纳贡献最大。代码上的改动集中在耦合设备参数和网络约束部分,其他都不动,非常能体现工作量。
5.4 用AI工具辅助Matlab代码生成
近几年AI代码辅助工具进步很快,也有不少人在问“Codex能不能像执行Python一样直接操作Matlab任务”。我的实测感受是:AI可以用来生成一段模型约束代码,但它不会帮你理解论文里的物理建模逻辑。比如你让它写Weymouth线性化,它写得像模像样,可参数设置、分段数选择、求解器兼容性这些细节仍然要自己把关。
我自己的做法是,先把论文中的公式逐条写在注释里,再让AI工具按注释生成初步代码,然后逐段检查约束是不是和公式一致。这样既省时间,又保留了核心的建模控制权。说到底,论文复现的本质是验证你对模型的理解,而不是生成一段能跑的代码。
最后再分享一点个人体会。我在复现这类论文时,最大的收获不是得到了一堆可用的Matlab代码,而是真正理解了规划模型里“投资决策”和“运行模拟”之间怎么互相作用。P2G厂站的位置和容量不是拍脑袋定的,而是由风电出力、电网阻塞、气网压力、设备效率和经济性共同决定的结果。你把这个过程亲手用代码实现一遍,才算是把“计及P2G厂站的电-气综合能源系统规划”这个课题真正吃透了。后续如果你想在这个方向深入,建议把代码里每个约束对应的物理含义都标注清楚,然后慢慢把单目标扩展成多目标,把典型日扩展成全年8760小时场景,再到加入不确定性鲁棒优化。这条路走通之后,再做其他综合能源系统规划论文,基本就是改网络数据和设备参数的事。