尧图网页设计 
最近后台好几个做能源方向的朋友来问同一个问题热电联供型微网优化运行的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(t1) SOC(t) η_chg * P_chg(t) * Δt / C - P_dis(t) * Δt / (η_dis * C)SOC是荷电状态%P_chg/P_dis是充放电功率kWC是储能容量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求解方案选型YalmipCplex、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有上界且只有当λ_i1时对应的段才被使用。这个过程我一般封装成一个函数后续加新设备直接调用。**如果不想处理线性化另一个思路是用智能优化算法如粒子群或遗传算法来求解模型但这类算法无法保证全局最优且对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_totalEmission_total包括购电对应的间接碳排放和天然气燃烧的直接碳排放。通过改变λ的取值可以画出一条帕累托前沿展示“多花多少钱才能减多少碳”。这个做法比直接跑NSGA-II稳定得多因为模型仍然是MILP只不过每次求一个加权解。我实际做下来λ从0到0.5变化时系统碳排放可以下降约18%而运行成本只上升5%左右。原因很简单低碳调度会让CHP多发电、替代购电因为燃气发电的碳排放因子通常低于燃煤电网同时电锅炉消纳风电也减少了弃风。这说明多能互补本身就有环境红利多目标优化只不过把这个红利显性化了。另外如果你未来要把这套代码扩展到园区级综合能源系统只需要增加氢能、冰蓄冷这类设备模型骨架完全不用变相当于在现有约束集合上增加“能量母线和转换设备模块”。这个扩展方向值得提前留好代码接口比如把设备参数结构体按统一字段定义后续加设备时多填一行参数即可。我自己在做过几个园区综合能源项目后最大的体会是这个领域的技术门槛不在单一设备建模而在如何把“源网荷储”用一个统一的优化框架串起来。这篇文章从系统边界、数学模型、求解器选型、代码框架、调试技巧讲到仿真结果分析和进阶方向基本覆盖了我做热电联供微网优化运行项目时完整的思考链路。最后说两条实操建议第一跑通模型后马上做参数错误注入测试比如故意把光伏出力翻倍看模型能否正确处理弃光这是检验模型鲁棒性的好办法第二代码注释不要只写“这是电平衡”要写“为什么这一项在等式右边”因为三个月后你回头看代码只会记得当时是怎么想的不会记得公式长什么样。希望这篇文章能让你少踩几个我踩过的坑顺利把算例跑通。 热电联供型微网优化运行的Matlab代码实现与求解器选型
最近后台好几个做能源方向的朋友来问同一个问题热电联供型微网优化运行的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(t1) SOC(t) η_chg * P_chg(t) * Δt / C - P_dis(t) * Δt / (η_dis * C)SOC是荷电状态%P_chg/P_dis是充放电功率kWC是储能容量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求解方案选型YalmipCplex、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有上界且只有当λ_i1时对应的段才被使用。这个过程我一般封装成一个函数后续加新设备直接调用。**如果不想处理线性化另一个思路是用智能优化算法如粒子群或遗传算法来求解模型但这类算法无法保证全局最优且对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_totalEmission_total包括购电对应的间接碳排放和天然气燃烧的直接碳排放。通过改变λ的取值可以画出一条帕累托前沿展示“多花多少钱才能减多少碳”。这个做法比直接跑NSGA-II稳定得多因为模型仍然是MILP只不过每次求一个加权解。我实际做下来λ从0到0.5变化时系统碳排放可以下降约18%而运行成本只上升5%左右。原因很简单低碳调度会让CHP多发电、替代购电因为燃气发电的碳排放因子通常低于燃煤电网同时电锅炉消纳风电也减少了弃风。这说明多能互补本身就有环境红利多目标优化只不过把这个红利显性化了。另外如果你未来要把这套代码扩展到园区级综合能源系统只需要增加氢能、冰蓄冷这类设备模型骨架完全不用变相当于在现有约束集合上增加“能量母线和转换设备模块”。这个扩展方向值得提前留好代码接口比如把设备参数结构体按统一字段定义后续加设备时多填一行参数即可。我自己在做过几个园区综合能源项目后最大的体会是这个领域的技术门槛不在单一设备建模而在如何把“源网荷储”用一个统一的优化框架串起来。这篇文章从系统边界、数学模型、求解器选型、代码框架、调试技巧讲到仿真结果分析和进阶方向基本覆盖了我做热电联供微网优化运行项目时完整的思考链路。最后说两条实操建议第一跑通模型后马上做参数错误注入测试比如故意把光伏出力翻倍看模型能否正确处理弃光这是检验模型鲁棒性的好办法第二代码注释不要只写“这是电平衡”要写“为什么这一项在等式右边”因为三个月后你回头看代码只会记得当时是怎么想的不会记得公式长什么样。希望这篇文章能让你少踩几个我踩过的坑顺利把算例跑通。