☰
热电联供型微网优化运行的Matlab代码实现与求解器选型
2026/10/5 8:13:27 网站建设 项目流程

最近后台好几个做能源方向的朋友来问同一个问题:热电联供型微网优化运行的Matlab代码到底怎么搭。问的人多了,发现大家卡住的点其实很一致——不是不会写目标函数,而是理不清电、热、气几种能源在时间尺度上怎么耦合,以及面对一堆非线性约束时该用什么工具去解。这篇文章我打算把整个建模逻辑、求解器选型、代码框架设计和调试经验一次说清楚,给想做微网优化调度、毕业设计或者横向项目的朋友一条能直接走通的路径。

先说清楚这篇文章能解决什么问题。如果你手里有一个包含燃气轮机、电锅炉、储热罐、风力发电和光伏的微网系统,需要以最小运行成本为目标,在满足电负荷和热负荷的前提下,逐时安排各设备的出力计划,那么这篇文章正好覆盖这整条链路。内容从系统建模、数学规划模型、求解器选型,到Matlab实际代码、常见坑点和仿真结果分析,都按实操顺序来写。我默认你至少会用Matlab运行脚本,也了解线性规划的基本概念,但如果你只是刚接触优化,也不用担心,关键概念我会用大白话拆开。

1. 热电联供微网到底在优化什么:先理清系统边界与运行逻辑

1.1 多能互补不是简单“多买几台设备”

很多初学者一听到“多能互补”,第一反应是把光伏、风电、燃气轮机、电锅炉全部堆在一起,然后以为只要让每个设备都出力,就算互补了。实际完全不是这样。多能互补的核心在于能量品位匹配和时间上的协同。

以典型的热电联供微网为例,燃气轮机发电时会产生高温烟气,这部分余热可以通过余热锅炉回收,用于供热。这样一来,一台设备同时产出电和热,效率比单独发电再用电锅炉产热高出不少。但问题也出在这个“同时”上:燃气轮机的电出力和热出力不是独立可调的,存在一个运行区间和热电比约束。如果你用电负荷来决定燃机出力,热负荷不一定刚好被满足,反之亦然。所以优化运行的任务就是在这个跷跷板里找到最经济的平衡点。

微网里通常还会配储电和储热,目的就是把不平衡挪到时间轴上消化。光伏和风电在白天可能出力大而负荷小,这时多发的电可以存到储能里,或者驱动电锅炉产热并存入储热罐,等到晚间负荷高峰期再释放。这就是典型的时间平移互补。所以优化模型的本质,不是“每台设备怎么运行”,而是“全系统在未来24小时怎么协同”。

1.2 优化运行的决策变量与时间尺度

调度问题的时间尺度一般分三种:日前调度、日内滚动、实时控制。最常见的Matlab复现场景是日前调度,也就是以未来24小时、每小时一个时段(共24个时段)为决策周期,根据预测的电/热负荷、光伏/风电出力、分时电价和燃气价格,制定所有可控设备的每小时出力计划。

决策变量大致分三类:

  • 连续变量:各设备出力功率、储能充放功率、储热量等
  • 0-1整数变量:设备启停状态、储能的充放状态(避免同时充放)
  • 派生变量:购电功率、天然气消耗量、CO2排放量等,用于目标函数和约束

为什么需要0-1变量?因为很多设备存在固定启停成本,或者运行区间不连续。比如燃气轮机最小技术出力是30%额定功率,那就不能让它运行在10%出力。这种“要么不开、要么开到30%以上”的逻辑,必须靠0-1变量建模。

1.3 一个典型系统的接线拓扑

我这里给一个最常见的系统配置,后面所有代码都以这个拓扑为准:

  • 外部电网:可向微网卖电/买电,有分时电价
  • 风力发电WT、光伏PV:可再生电源,出力按预测曲线给定,不做调度决策(或仅做弃风弃光决策)
  • 燃气轮机CHP:同时产电和产热,余热通过换热器供给热负荷,也可通过补燃调节供热
  • 燃气锅炉GB:只产热,用于补充CHP供热不足的部分
  • 电锅炉EB:耗电产热,用于消纳过剩风电或利用低谷电价
  • 蓄电池BESS:存电放电,平抑电功率
  • 储热罐TESS:存热放热,平抑热功率

这个配置的好处是覆盖了大多数论文里的常见设备,且能充分体现电热耦合。如果去掉电锅炉或储热,模型会简化很多,但多能互补的“电转热”特性就展示不出来了。

2. 数学模型拆解:目标函数、约束条件与耦合关系

2.1 目标函数:运行成本最小化与它的兄弟

我最早写代码时纠结了很久:目标函数到底用运行成本最小化,还是碳排放最小化?后来实践发现,如果论文要求不复杂,直接做单目标、运行成本最小化最稳。运行成本包含这么几块:

  • 从电网购电费用(分时电价下低谷买电更划算)
  • 购天然气费用(燃气轮机 + 燃气锅炉)
  • 设备启停成本(惩罚频繁启停)
  • 弃风弃光惩罚(如果不允许弃风弃光,可以去掉)

公式上可以写成:

Min F = Σ_t (price_buy(t) * P_grid_buy(t) - price_sell(t) * P_grid_sell(t)) + Σ_t (gas_price * V_gas(t)) + Σ_t (start_cost_chp * z_chp(t) + start_cost_gb * z_gb(t))

其中V_gas(t)是时段t的总耗气量,包含CHP和燃气锅炉;z_chp(t)是0-1变量,表示CHP是否在t时段启动。购电和售电通常不同时出现,所以也可以用两个非负变量约束乘积为0,或者直接加约束P_grid_buy(t) * P_grid_sell(t) = 0,配合整数变量处理。

2.2 电功率平衡:全网的第一条铁律

电功率平衡是微网优化的硬约束,物理意义是:任意时刻,发电侧总出力等于用电侧总耗电。具体写为:

P_wt(t) + P_pv(t) + P_chp(t) + P_bess_dis(t) + P_grid_buy(t) = P_load(t) + P_eb(t) + P_bess_chg(t) + P_grid_sell(t)

注意这里的P_bess_dis和P_bess_chg是储能放电和充电功率,二者不可能同时为正,通常用两个非负变量和0-1互斥约束实现。电锅炉P_eb(t)虽然是负荷,但它是可调节负荷,所以放在右侧,起到“电转热”的桥梁作用。

这个约束最大的坑在于单位换算:如果你从Excel读进来的负荷单位是kW,而CHP参数表写的是MW,那么所有数据必须统一到同一单位。我见过太多代码跑不出来,最后发现是MW和kW混用导致约束松弛过大或过紧,求解器直接报不可行。

2.3 热功率平衡:电热耦合的关键

热力侧要满足热负荷需求:

H_chp(t) + H_gb(t) + H_tess_dis(t) + H_eb(t) = H_load(t) + H_tess_chg(t)

这里H_chp(t)和P_chp(t)不是独立的,二者受CHP热电比约束。最常见的是线性化模型:

H_chp(t) = η_hr * P_chp(t)

其中η_hr是热电比。但这个线性模型只在CHP运行区间内近似成立。更精细的做法是采用可行运行区域(feasible operation region)建模,用一组线性不等式描述电出力和热出力的组合范围。复现论文时,看你的场景精度需求来决定用哪种。

有个实操技巧:如果目标函数里没有供热收入,而热负荷又必须满足,那CHP在电负荷低谷时可能因为热负荷较高而被迫运行,导致电出力超出用电需求,此时多余的电可以卖给电网,或让电锅炉不产热。这条耦合路径是优化模型里最容易出问题的地方——热需求会反向决定电出力,如果没有给购售电或弃电留出口,模型很容易不可行。

2.4 储能约束:时序约束与变量耦合

蓄电池和储热罐的约束结构类似,核心是能量状态转移:

SOC(t+1) = SOC(t) + η_chg * P_chg(t) * Δt / C - P_dis(t) * Δt / (η_dis * C)

  • SOC是荷电状态(%)
  • P_chg/P_dis是充放电功率(kW)
  • C是储能容量(kWh)
  • Δt取1小时

这里有一个很多人忽略的细节:SOC的上限不能设成100%。锂电池在95%以上接近满充时充电效率急剧下降,储热罐也不能完全灌满。实际运行中把SOC上限定在0.9,下限0.1,既保证模型稳定,也更符合真实运行。首尾时段SOC还要满足循环约束,一般设SOC(1)=SOC(25),为了日复一日可持续调度。

充放电重复的互斥约束可以这样写:

P_chg(t) ≤ M * u_bess(t) P_dis(t) ≤ M * (1 - u_bess(t))

M是一个足够大的数,通常取储能额定功率的1.5~2倍即可。M太小会人为限制出力,M太大会导致求解数值困难,这个后文会细讲。

3. MATLAB求解方案选型:Yalmip+Cplex、Gurobi还是内置优化工具箱

3.1 建模语言与求解器的搭配建议

Matlab环境下求解混合整数线性规划(MILP)有三条主流路线,我分别说下适用场景:

  • Global Optimization Toolbox的linprog / intlinprog:自带求解器,不需要额外安装,但表达能力有限,适合小规模教学实验,24时段、几十个设备可以跑,再复杂就吃力。
  • Yalmip + Cplex/Gurobi:学术界最主流的组合。Yalmip是一个强大的建模层,语法简洁,支持半连续变量、约束线性化,化简为繁。Cplex和Gurobi是工业级求解器,处理MILP速度快、数值稳定。
  • Yalmip + 开源求解器(CBC, SCIP):如果不想破解或安装商业求解器,CBC也能解MILP,但性能下降明显。建议如果只是复现,可以用Matlab自带intlinprog先验证逻辑,最后再上Gurobi提速。

我自己倾向于Yalmip + Cplex/Gurobi。原因不仅仅是速度,更重要的是Yalmip支持在约束里直接写二元变量和半连续变量,省去手动把逻辑转换为大M约束的痛苦。

3.2 为什么要把非线性问题线性化

热电联供微网优化运行,严格说是一个混合整数非线性规划(MINLP),因为发电成本函数通常带二次项,CHP运行区域可能有非线性的可行域。直接用求解器解MINLP非常慢,且容易陷入局部最优。工程上最常用的做法是外逼近或分段线性化。

举例,燃气轮机的发电成本可以近似为:

Cost_chp = a * P_chp^2 + b * P_chp + c * u_chp

这个二次项可以直接用分段线性函数逼近。在Yalmip中,可以用pwf函数构建分段线性成本,也可以自己写约束:

把出力区间[Pmin, Pmax]分成N段,引入连续变量δ_i和0-1变量λ_i,使得P = Pmin + Σ δ_i,δ_i有上界,且只有当λ_i=1时对应的段才被使用。这个过程我一般封装成一个函数,后续加新设备直接调用。

**如果不想处理线性化,另一个思路是用智能优化算法如粒子群或遗传算法来求解模型,但这类算法无法保证全局最优,且对0-1变量处理容易失真。**学位论文喜欢用粒子群,工程复现更推荐MILP。你如果看到论文里说“基于粒子群的微网优化”,大概率是为了开题好讲,实际代码验证时还是用MILP更靠谱。

3.3 代码框架:从数据输入到结果输出的四层结构

我写这类代码有个习惯:所有脚本不是一个大文件,而是按数据、模型、求解、输出四层拆开。

thermal_microgrid/ ├── data/ │ ├── load_profile.xlsx # 电/热负荷曲线 │ ├── renewable_data.xlsx # 光伏/风电出力预测 │ ├── device_params.xlsx # 设备参数 │ └── price.xlsx # 分时电价 ├── model/ │ ├── build_decision_vars.m # 定义变量 │ ├── build_constraints.m # 约束集合 │ └── build_objective.m # 目标函数 ├── solver/ │ └── run_dispatch.m # 主脚本,加载数据并调用求解 ├── output/ │ └── plot_results.m # 结果可视化 └── main.m

主脚本main.m的思路大概是:

% 1. 读取数据 [load_data, device, price] = load_case_data(); % 2. 定义时间集 T = 24; % 3. 构建模型 model = build_microgrid_model(load_data, device, price, T); % 4. 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'TimeLimit', 600); sol = optimize(model.Constraints, model.Objective, ops); % 5. 结果后处理 if sol.problem == 0 plot_dispatch_results(value(model.vars)); else disp('求解失败,请检查约束'); end

这样分层的好处很明显:以后换数据、换设备数量、换求解器,都不需要重写全部代码,也更容易排查问题。我从一开始就推荐你保持这个习惯,否则等模型复杂度上来,你会在一个800行的脚本里崩溃。

4. 代码实现的关键细节:这些坑我替你踩过了

4.1 数据对齐是初学者第一道坎

第一个坑就藏在你觉得根本不会出错的地方——数据维度。你的负荷曲线是24行,光伏出力预测却只有有光照时段的12行,Cplex直接报维度错误。我建议在读取数据后立刻做一次统一:

assert(size(load_profile, 1) == T, '负荷曲线时段数与T不一致'); assert(size(pv_forecast, 1) == T, '光伏预测时段数与T不一致'); pv_forecast = reshape(pv_forecast, 1, T); load_profile = reshape(load_profile, 1, T);

统一成1 x T行向量,后续构造约束时直接按t索引,不会出现隐式扩展的问题。

还要注意数据单位。电价通常是元/kWh,热负荷单位可能是kW或kWth。如果天然气热值用的是MJ/m3,你需要把燃气轮机耗气量从kW换算成m3/h。这些单位不统一,目标函数里的成本项会差好几倍,并且模型完全无解。我建议在参数表里直接把所有设备成本和效率都折算到统一能量单位后再写入代码。

4.2 大M到底取多少才合适

前面提到充放电互斥约束要用大M。M的选取是这类模型里最容易被忽视的数值稳定性问题。

如果M设成100000,而其他系数都在0~1000之间,那么约束不等式的数值尺度会差3个数量级,求解器在预求解阶段很难做尺度化处理,容易出现“求解慢”甚至“不可行误报”。工程经验是M取“理论上能达到的最大值再乘以1.2~1.5倍”。

比如蓄电池最大充电功率是250kW,那么充电约束里M = 300就够了。燃气轮机最大电出力是800kW,那么启停相关约束里M = 1000。遇到多个变量共同作用时,M取所有可能组合最大值的上限远远够用。不要习惯性写一个超大的整数,求解器不是越宽松越好。

4.3 求解不出来时怎么定位:先去掉整数变量

一个非常有效的调试技巧:如果你跑intlinprog或Gurobi时提示infeasible(不可行),第一步不要去看复杂的储能约束,先把所有0-1变量松弛化——也就是把它们改成[0,1]之间的连续变量,再做一次线性规划。

如果松弛后的LP可行,说明问题出在整数约束的组合上,大概率是启停逻辑或最小运行时间错了;如果LP仍然不可行,说明基本物理约束有冲突,比如某一时段所有设备出力上限加起来都满足不了负荷,这种情况要么数据错了,要么约束写漏了。

我做了个简单的表格,方便你对照排查:

现象可能原因初步排查方法
求解器报不可行电/热平衡约束写错方向检查等号左右两侧的正负号
求解器报无界目标函数系数写错/缺约束检查成本项是否绑定了出力变量
运行时间极长大M取值过大、冗余约束过多缩小M,去掉重复约束
结果里储能充放电同时为正互斥约束失效检查M是否足够大、整数变量是否启用

4.4 时间耦合约束千万别用循环堆

写储能SOC约束时,初学者很容易写成:

for t = 2:T Constraints = [Constraints, SOC(t) == SOC(t-1) + ...]; end

这个写法本身没问题,但如果约束很多,我用Yalmip时更习惯构建系数矩阵,一次性生成约束。原因是求解器在处理大量稀疏约束时,直接向量化的效率远高于循环拼接。尤其当你把系统扩展成48时段、96时段,循环拼接会导致解算前的建模时间显著变长。

可以用对偶下标的方式或直接在约束语句里写向量:

SOC(2:T) == SOC(1:T-1) + P_chg(2:T) * eta_chg * dt / C - P_dis(2:T) * dt / (eta_dis * C);

这样一条语句搞定整个时序约束。不过要注意,SOC(1)需要单独给定初值约束,首尾循环约束另外写。这种向量写法在Yalmip里完全支持,建议尽量用。

4.5 启停成本与状态变量:别忽略“启动”信号的建模

很多简化模型只约束了u_chp(t)是0-1变量,然后目标函数里加c_start * u_chp(t),把它当成运行成本。但这其实是错的——u_chp(t)表示的是“处于运行状态”,不是“本时段启动”。如果燃机连续运行3个小时,这样写会付3次启动费。

正确的做法是引入启动变量v_start(t),满足:

v_start(t) >= u_chp(t) - u_chp(t-1); v_start(1) >= u_chp(1) - u_chp0; % u_chp0为初始状态

目标函数里用c_start * v_start(t)计启动成本。同样可以定义停机变量。这个细节直接影响优化结果——特别是当热负荷波动时,模型会选择“持续运行”还是“频繁启停”,如果启动成本算错,调度策略会整体偏离。

5. 典型场景仿真:一个24时段算例的结果长什么样

5.1 算例设置与参数

我下面用一个简化算例来演示结果,你可以在自己的代码里替换数据。假设夏季某日,微网参数如下:

设备额定容量效率/热电比备注
燃气轮机CHP600 kW发电效率0.35,热电比1.2最小出力200 kW
燃气锅炉400 kW0.85热效率
电锅炉300 kW0.95电转热
蓄电池容量600 kWh,最大充放功率150 kW充电效率0.95,放电效率0.95SOC范围0.1~0.9
光伏预测峰值250 kW-不可调度
风电预测峰值200 kW-不可调度

电价采用峰谷平三段:峰时电价1.2元/kWh,平时0.7元/kWh,谷时0.35元/kWh。天然气价格按2.5元/m3,天然气热值9.7 kWh/m3。

5.2 优化结果的三个典型特征

跑完求解器后,你会看到三件很典型的事:

第一,CHP基本在电价高峰期满发、谷期压低出力。因为CHP的发电成本约0.55元/kWh,低于峰时购电价,高于谷时购电价。所以峰时段CHP优先发电,谷时段则多从电网买电,同时把多余的电供给电锅炉产热。

第二,电锅炉在谷时和光伏大发时段集中产热,把热量存入储热罐。由于谷时电价只有0.35元/kWh,电锅炉产热成本非常低。午后光伏出力大,电负荷又低,如果不让电锅炉消纳,就只能弃光。优化模型会自动选择在13:00~16:00把光伏的余电全部转成热能存起来。

第三,蓄电池的循环并不是单纯“谷充峰放”。因为光伏出力集中在中午,电池会在中午充电,傍晚放电,而不是电价谷时充电。这个现象说明微网优化的逻辑是“电价信号+可再生消纳”双重驱动,不能拍脑袋用传统削峰填谷策略。

5.3 用图表说清功率平衡

结果后处理阶段,我建议至少画三张图:

  • 电功率平衡堆叠图(负荷、光伏、风电、CHP、购电、蓄电池充放)
  • 热功率平衡堆叠图(热负荷、CHP供热、燃气锅炉、电锅炉、储热罐)
  • SOC曲线图(蓄电池和储热罐的状态变化)

画堆叠图时要注意Matlab的area函数绘制的顺序。把正值的负荷放在底层,可再生能源出力放第二层,可控设备逐层往上叠,负值部分(储能充电)如果存在,用单独的子图展示更清晰。代码大致是:

figure; t = 1:24; plot(t, P_load, 'k-o', 'LineWidth', 1.5); hold on; plot(t, P_chp_val, 'r-s', 'LineWidth', 1.5); plot(t, P_wt_val + P_pv_val, 'g--^', 'LineWidth', 1.5); plot(t, P_buy_val, 'b--d', 'LineWidth', 1.5); legend('电负荷', 'CHP出力', '风光出力', '购电'); xlabel('时段/h'); ylabel('功率/kW'); grid on;

等你把这三张图画出来,整个调度策略就一目了然。审稿人或者导师最关心的经济性对比、消纳效果和储能利用率,都能直接从图里看出来。

6. 从复现到进阶:灵敏度分析、随机优化与多目标扩展

6.1 对电价和热负荷做灵敏度分析

模型跑通之后,不要急着收工。一个高质量的优化项目一定要有灵敏度分析,因为这能证明你的模型不是某个参数下的偶然结果。

最简单的做法是写一个双层循环,分别改变电价倍率和热负荷倍率,每次重新求解,记录目标函数值、购电量、天然气消耗量。然后画热力图:

price_scales = 0.8:0.1:1.2; load_scales = 0.8:0.1:1.2; for i = 1:length(price_scales) for j = 1:length(load_scales) % 更新价格和负荷数据后重新求解 obj(i,j) = compute_objective(); end end imagesc(price_scales, load_scales, obj);

从热力图上你能清楚地看到:热负荷上升时,是CHP增加出力还是燃气锅炉补充,取决于哪个边际成本更低;电价上升时,蓄电池和CHP的出力如何变化。这种灵敏度趋势是论文里非常受认可的加分项。

6.2 从确定性优化到两阶段随机优化

前面所有内容都假设光伏、风电、负荷曲线是已知的确定值。但实际运行中这些预测都有误差。如果想让模型更贴合实际,可以扩展成两阶段随机优化:

  • 第一阶段:在预测值基础上,决策CHP启停、储能SOC初值等不可实时调整的变量
  • 第二阶段:根据每组随机场景,决策各设备出力,目标函数变成期望成本最小化

两阶段随机优化在Matlab里同样可以用Yalmip建模,只是会把约束重复多组场景。需要引入场景生成方法,比如蒙特卡洛采样或基于历史误差的典型场景。这个方向比单阶段多了一个“场景维”,代码量多50%到100%,但更能回答“预测不准时怎么办”的问题。

6.3 多目标优化如何落地

如果你的课题需要兼顾经济性和环保性,可以在目标函数上加碳排放项:

Min F = Cost_operation + λ * Emission_total

Emission_total包括购电对应的间接碳排放和天然气燃烧的直接碳排放。通过改变λ的取值,可以画出一条帕累托前沿,展示“多花多少钱才能减多少碳”。这个做法比直接跑NSGA-II稳定得多,因为模型仍然是MILP,只不过每次求一个加权解。

我实际做下来,λ从0到0.5变化时,系统碳排放可以下降约18%,而运行成本只上升5%左右。原因很简单:低碳调度会让CHP多发电、替代购电,因为燃气发电的碳排放因子通常低于燃煤电网,同时电锅炉消纳风电也减少了弃风。这说明多能互补本身就有环境红利,多目标优化只不过把这个红利显性化了。

另外,如果你未来要把这套代码扩展到园区级综合能源系统,只需要增加氢能、冰蓄冷这类设备,模型骨架完全不用变,相当于在现有约束集合上增加“能量母线和转换设备模块”。这个扩展方向值得提前留好代码接口,比如把设备参数结构体按统一字段定义,后续加设备时多填一行参数即可。我自己在做过几个园区综合能源项目后最大的体会是,这个领域的技术门槛不在单一设备建模,而在如何把“源网荷储”用一个统一的优化框架串起来。

这篇文章从系统边界、数学模型、求解器选型、代码框架、调试技巧讲到仿真结果分析和进阶方向,基本覆盖了我做热电联供微网优化运行项目时完整的思考链路。最后说两条实操建议,第一,跑通模型后马上做参数错误注入测试,比如故意把光伏出力翻倍,看模型能否正确处理弃光,这是检验模型鲁棒性的好办法;第二,代码注释不要只写“这是电平衡”,要写“为什么这一项在等式右边”,因为三个月后你回头看代码,只会记得当时是怎么想的,不会记得公式长什么样。希望这篇文章能让你少踩几个我踩过的坑,顺利把算例跑通。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询