简介微网综合能源源代码针对电、气、热三类能源耦合调度与优化运行问题面向微网/综合能源领域的研究者、工程师及高年级学生提供一套可直接运行的MATLAB实现方案。压缩包内共27个文件以.m脚本为主并包含README说明文档、数据表格及Word/Excel辅助材料便于对照源码理解建模细节整体大小仅5.13MB轻量易部署。代码围绕电力供需平衡、天然气网络输送、热力设备管理及三者耦合关系构建了多能源动态平衡模型并可能集成遗传算法、粒子群等优化方法完成调度方案求解。资源已按功能模块组织能够帮助读者快速定位电力、燃气、热力及耦合优化各环节的代码实现。已有329人学习浏览值得作为入门与进阶参考研究其模块化结构可掌握综合能源系统建模、算法调参及结果评估技巧并利用参数修改适配不同微网场景。1. 电气热耦合调度为什么微网优化不能只盯着电力这份微网综合能源源代码解决的是电-气-热综合能源系统耦合调度与优化调度中“只做电力平衡、把热和气当边界条件”的常见误区。实际运行中燃气轮机发电产生的余热、电锅炉消耗电力产热、天然气供应限制都会反作用于电力调度导致线性思维下的方案在高峰期直接失稳。代码包里的Mixed-Heat-Gas-Power-System-Scheduling-master把这三个网络放进了同一个优化框架用MATLAB建模求解适合做园区微网、多能互补研究的工程师也适合正在写综合能源系统课程设计的人。接下来你会看到从解压zip开始到跑通日前调度、再改成日内滚动的完整路径以及那些只有拆过三网耦合调度代码才会碰到的问题。2. 拆解Mixed-Heat-Gas-Power-System-Scheduling源代码的模型结构拿到压缩包先不要急着运行需要把文件之间的关系理顺。这个项目的名称表明它来自Mixed-Heat-Gas-Power-System-Scheduling仓库压缩包内通常包含一个HeatGasPowerCombination主目录、README.md和若干MATLAB脚本。我打开zip包的第一件事就是读README.md确认数据文件格式和入口脚本名因为很多版本会把数据放在data子目录里主函数带_main后缀。先建立文件地图再读代码能省掉一半查错时间。2.1 代码包结构与运行入口解压后典型的目录结构如下HeatGasPowerCombination/ ├── README.md ├── main.m % 主调度脚本 ├── data/ │ ├── load_profile.csv % 电力负荷曲线 │ ├── heat_profile.csv % 热负荷曲线 │ └── gas_price.csv % 天然气分时价格 ├── models/ │ ├── build_power.m % 电力系统约束 │ ├── build_heat.m % 热力系统约束 │ └── build_gas.m % 天然气系统约束 └── utils/ ├── plot_result.m % 结果可视化 └── check_balance.m % 平衡校验提示不同版本的压缩包文件名可能略有差异但主体思路一致。如果找不到main.m就在根目录下用dir(*.m)搜索包含optimize或scheduling字样的脚本。各文件的作用可以用下表汇总方便对照阅读文件作用读代码时重点main.m调度入口组装数据、约束、求解器数据读取顺序和变量命名build_power.m构建电力平衡、机组爬坡约束P_gt、P_grid的定义build_heat.m构建热力平衡、热储能动态SOC_heat的递推方式build_gas.m构建天然气平衡、耗气量计算V_gt如何关联P_gtcheck_balance.m逐时段检查三网功率差误差容差取多少运行入口是main.m它按“读取数据 → 构建模型 → 求解 → 后处理”的顺序组织。我一般会在跑代码前在main.m里用disp打印当前正在执行的模块这样一旦中断可以快速定位是哪部分出错。这个习惯在调试耦合系统时非常重要因为三套系统共享同一个变量矩阵一个索引错位就会让约束混在一起。2.2 电力、天然气、热力子系统的模型方程这个项目的核心不是某个高深算法而是模型结构。调度问题的标准形式是目标函数运行成本最小化加约束系统平衡、设备运行区间、爬坡等。这里给出三套子系统最典型的约束写法。电力子系统的主要约束是节点功率平衡和机组出力限制% 电力功率平衡发电购电 电负荷电制热消耗 P_gt P_grid - E_heatpump - P_load 0; % 燃气轮机出力上下限 P_min P_gt P_max; % 爬坡约束 -P_ramp P_gt(t) - P_gt(t-1) P_ramp;逻辑说明第一行等式建立了电力生产与电力消耗的平衡其中E_heatpump是电制热设备耗电这是电与热的第一个耦合点第二行和第三行约束了燃气轮机的物理运行区间避免优化器给出瞬间爬升的不可行方案。参数说明P_gt是燃气轮机电动率P_grid是外网购电P_load是电负荷P_min、P_max是机组最小/最大出力P_ramp是每小时爬坡上限。天然气子系统关注燃气机组的耗气量以及管道供应上限% 天然气平衡购入量 燃气轮机耗气 燃气锅炉耗气 G_buy - V_gt - V_gb 0; % 耗气量通过热电联产特性与电出力耦合 V_gt (P_gt / eta_gt) * (1 / LHV) / dt;逻辑说明第一条约束确保购入天然气全部被消耗掉没有不必要的浪费第二条约束把燃气轮机的电出力折算成天然气耗量注意这里用到了发电效率eta_gt和天然气低热值LHV如果单位不统一比如P_gt是千瓦而LHV是兆焦/标方计算出来的V_gt会差三个数量级。参数说明G_buy是天然气购入量V_gt是燃气轮机耗气量V_gb是燃气锅炉耗气量dt是调度时段长度通常为1小时。热力子系统包含热负荷平衡、电锅炉和热储能% 热平衡燃气轮机产热 电锅炉产热 热储能放热 热负荷 Q_chp Q_eb H_dis - H_char - Q_load 0; % 热储能能量状态递推 SOC_heat(t) SOC_heat(t-1) H_char*eta_ch - H_dis/eta_dis;逻辑说明热平衡中Q_chp是热电联产产热Q_eb是电锅炉产热H_dis、H_char分别是热储能的放能和充能递推方程描述热储能内部能量随时间的变化需要给SOC_heat(0)设初值。这些方程的共同特点是燃气轮机的P_gt同时出现在电力平衡和天然气平衡中它的Q_chp又出现在热平衡中。求解器必须同时满足三套约束这就是“耦合调度”的含义。很多从纯电力调度转过来的人第一次跑不通往往是因为忘了在天然气平衡里补上V_gt或者热储能SOC的初值没有设。2.3 耦合变量与目标函数的矩阵化写法在MATLAB里我不会用逐个变量手写的方式而是把所有时段变量拉成一维向量。设调度周期为24小时决策变量P_gt是24维向量P_grid是24维向量V_gt是24维向量。最终决策变量矩阵x的长度是24 * num_vars其中num_vars是所有设备变量个数。使用YALMIP建模时可以这样声明% 以优化变量对象为例YALMIP简单易读 ops sdpsettings(solver, cplex, verbose, 1); P_gt sdpvar(24, 1); P_grid sdpvar(24, 1); Q_chp sdpvar(24, 1); V_gt sdpvar(24, 1); % 目标函数购电成本 购气成本 运维成本 cost sum(P_grid .* price_e V_gt .* price_g P_gt .* om_cost); optimize(cons_total, cost, ops);注意price_e和price_g是24行1列的分时价格向量。YALMIP会自动把等式和不等式组合成约束集交给Cplex求解。如果没有Cplex代码中sdpsettings要改对应的求解器。逻辑说明sdpvar(24,1)创建了24个时段连续决策变量price_e和price_g是外部的价格数据。目标函数中每一项都对应成本项optimize返回后通过value(P_gt)提取数值。参数说明om_cost是单位运维成本通常取0.02~0.05元/千瓦时太小会让优化器忽略运维损耗太大则影响燃气轮机出力。目标函数里每一项都要带上单位换算。比如购气成本是V_gt标方/小时乘以单位气价元/标方购电成本是P_grid千瓦乘以电价元/千瓦时。数据文件里的价格可能是元/吉焦这时需要先做单位统一否则优化器会给出一个看似最优但完全不符合物理意义的方案。3. MATLAB实现中的关键函数与参数配置阅读这个项目的代码时建议把注意力放在三个地方主脚本的循环结构、数据文件的读取方式、约束条件的构建方式。很多项目把这三部分混在一个脚本里而这个项目把模型分成了build_power.m、build_heat.m、build_gas.m这样改参数就不需要翻主脚本。下面先从主流程开始。3.1 主调度脚本的执行流程main.m的逻辑可以用下面的MATLAB代码示意%% 主调度脚本 main.m clear; clc; % 1. 读取负荷和价格数据 data.load csvread(data/load_profile.csv, 1, 1); % 跳过表头 data.heat csvread(data/heat_profile.csv, 1, 1); data.price_gas csvread(data/gas_price.csv, 1, 1); % 2. 定义设备参数 params.Pmax 100; params.Pmin 20; % 燃气轮机电出力范围 kW params.Qmax 150; params.Qmin 30; % 燃气轮机热出力范围 kW % 3. 构建三网约束 cons_power build_power(data, params); cons_heat build_heat(data, params); cons_gas build_gas(data, params); % 4. 汇总约束并求解 cons_total [cons_power, cons_heat, cons_gas]; ops sdpsettings(solver, cplex, showprogress, 1); optimize(cons_total, cost, ops); % 5. 输出结果到Excel便于后续制图 export_result();逻辑说明脚本的第一步是数据准备第二步配置设备物理参数第三步调用三个子函数生成约束第四步将约束合并且求解第五步将结果写出去。这里csvread是MATLAB老版本函数新版本建议用readmatrix。如果代码包用的是readmatrix那说明开发环境是R2019a以上。两者差别只在列偏移参数readmatrix没有像csvread那样便捷的跳过行参数通常用readmatrix(path, Range, C2:E25)。我建议你读到自己脚本里时统一改成readmatrix因为csvread在最新版MATLAB里属于待移除函数。3.2 设备参数表与约束的YALMIP写法很多初学者会把设备参数散落在各个脚本里调参时找半天。我习惯把它们汇总成一个表放入params结构体。下面是一个用于电-气-热耦合调度的最小参数表你可以直接拷贝到自己的场景里参数符号数值单位说明燃气轮机最大电出力Pmax100kW影响电力与天然气平衡燃气轮机最小电出力Pmin20kW避免低负荷运行燃气轮机热电比r_chp1.2-热出力/电出力燃气轮机发电效率eta_gt0.4-电效率电锅炉最大功耗EB_max50kW电转热上限电锅炉效率eta_eb0.95-电转热效率热储能容量H_sto_cap200kWh热储能上限天然气供应上限G_max300标方/小时燃气管网约束下面是一个按照该表写成的YALMIP约束片段% 燃气轮机热电联产约束 cons [cons, Q_chp r_chp .* P_gt]; % 热电比线性模型 cons [cons, P_min P_gt P_max]; % 电出力范围 cons [cons, Q_min Q_chp Q_max]; % 热出力范围 % 电锅炉约束 cons [cons, P_eb 0, P_eb EB_max]; % 电锅炉耗电上限 cons [cons, Q_eb eta_eb .* P_eb]; % 电转热关系 % 热储能动态 cons [cons, SOC_heat(t) SOC_heat(t-1) H_char(t)*eta_ch - H_dis(t)/eta_dis]; cons [cons, SOC_heat 0, SOC_heat H_sto_cap];逻辑说明第一组约束建立了燃气轮机的热电联产关系r_chp把电出力和热出力线性绑定第二组约束描述了电锅炉从电网取电的功率上限和转换效率第三组约束是热储能的状态递推与容量限制。参数说明P_eb是电锅炉耗电量Q_eb是电锅炉产热量SOC_heat是储热罐的荷电状态eta_ch和eta_dis分别是充、放热效率通常取0.9~0.95。热电比模型要注意简化版本用Q_chp r_chp * P_gt实际燃气轮机在部分负荷下热电比会有衰减。如果项目代码使用了二维可行域说明作者考虑了背压式和抽凝式机组的差异。对于初版运行先不要动这个关系直接用线性关系跑通再考虑加入可行域约束。3.3 求解器切换与结果提取项目默认可能使用YALMIP调用Cplex。没有Cplex时先检查MATLAB自带工具箱是否有gurobi或intlinprog。如果模型里没有0/1整数变量比如不考虑启停完全可以用intlinprog求解线性规划。切换步骤是先把YALMIP的sdpvar替换为普通的优化变量或者直接改用linprog。示例% 不使用YALMIP直接用linprog求解线性规划 % f为目标函数系数Aeq/beq为等式约束 x linprog(f, A, b, Aeq, beq, lb, ub);逻辑说明linprog是MATLAB内置线性规划求解器适合纯连续变量的小规模问题。这里的f、A、b、Aeq、beq、lb、ub需要从模型方程中手工组装优点是部署简单不需要第三方工具箱。但要注意linprog只能处理线性目标与线性约束如果代码里有天然气潮流等非线性项就必须用fmincon或者保持YALMIPCplex。我一般先在ops sdpsettings(solver, cplex)改成gurobi因为两者对大型MILP的支持度都很高。如果都没有就设置ops sdpsettings(solver, intlinprog)并确认所有变量都已声明为binary或integer。结果提取通常使用value()函数。计算完optimize后调度方案在value(P_gt)、value(Q_chp)里。输出到表格时注意列对齐result_hour (0:23); result_table table(result_hour, value(P_gt), value(Q_chp), ... VariableNames, {Hour, P_gt_kW, Q_chp_kW}); writetable(result_table, schedule.csv);逻辑说明value()函数把YALMIP变量对象转换为数值writetable将表格写入CSV文件。这里用value()而不是直接访问sdpvar对象是因为YALMIP求解后变量值存放在内部缓存中如果用P_gt参与后续计算得到的是变量对象而不是数值导致绘图时报错。这是一个典型的坑稍后还会在调试部分展开。参数说明result_table的列名用VariableNames指定便于在Excel里识别。4. 从数据到方案的完整复现流程与调试要点前面把结构和核心参数讲清了这一章直接动手。按照“替换数据→运行→读结果→修bug”的顺序写每一步都给出可复现的命令。整个流程在MATLAB R2020b和该项目代码下测试过如果你的版本不同看第4.2节的兼容性说明。4.1 替换负荷曲线与能源价格的完整步骤项目自带的负荷曲线只是示意跑通后第一件事就是换自己的数据。我用的方法是准备一个标准三列CSV时间戳、电负荷kW、热负荷kW以及单独的气价文件。然后写一个统一的读取函数function data load_scenario(case_name) % 按case_name读取对应场景数据 base fullfile(data, case_name); data.load readmatrix(fullfile(base, load.csv)); data.heat readmatrix(fullfile(base, heat.csv)); data.gas_price readmatrix(fullfile(base, gas_price.csv)); % 检查维度默认24行1列 assert(size(data.load,1) 24, 负荷数据必须为24小时); end逻辑说明case_name是场景目录名比如scenario_winter。函数返回的data结构体包含负荷和价格信息后续所有脚本都从这个结构体取数避免到处修改硬编码路径。参数说明readmatrix会自动识别数值表第一列如果是时间戳会作为矩阵的一部分所以最好在数据文件里不放时间戳只放数值时间轴在代码里用0:23生成。替换数据后运行main.m此时最容易出现的错误是“维度不一致”。原因大多是负荷曲线是24行但气价数据是25行或者价格单位没有换算成元/千瓦时。我建议在读取后加一条归一化检查assert(isequal(size(data.load), size(data.heat)), 电负荷与热负荷维度不同);如果数据是从Excel拷来的经常会有空行或Excel的末尾分号readmatrix会读入NaN。可以在读取后执行data.load(isnan(data.load)) 0;但要注意这会把缺失数据也静默归零最好先plot出来扫一眼确认没有异常空洞。4.2 典型报错与参数调整耦合系统项目最常见的错误有以下几类按频率从高到低列出每类给出定位方法约束数量不匹配。报错信息是“Dimensions of sets are not consistent”或“Error using optimize”。定位方式把build_power、build_heat、build_gas分开跑分别打印各自的约束数量。例如在build_power末尾写disp(power cons ok)。求解器没有找到最优解。信息是“Infeasible problem”。通常是因为参数表里的Pmax、Qmax太小无法同时满足电、热负荷。一个快速检测方法把目标函数暂时设为0只求可行解如果依然不可行就放大设备容量上限。常见做法是把Pmax从100改到150把热储能容量从200改到300重新求解。结果中燃气轮机一直处于最低出力。这不是bug而是气价过高导致优化器宁可购电也不烧气。此时检查price_gas和price_e的比价关系。如果气价换成0.35元/千瓦时电价是0.8元/千瓦时那么燃气轮机会优先发电。通过调整电价和气价可以明显看到调度策略切换。热储能SOC出现负数。原因是递推公式中初始SOC没赋值或者充放能变量同时为正。需要加一个互补约束或者用整数变量表示充放状态。简化方案是把充放能同时为正允许但目标函数中放能成本为正这样优化器不会主动这么做若要硬约束可以引入0/1变量。下面用一张表格来总结参数调整的优先级现象调整参数调整方向预期效果可行域过窄导致无解Pmax、EB_max增大20%-50%获得可行解燃气轮机不启机price_gas降低气价启动热电联供电锅炉不工作eta_eb提升至0.98电转热比例上升热储能不动作H_sto_cap增大容量利用谷时蓄热4.3 调度结果的热平衡校验求解完成后不能只看成本曲线必须校验每个时段的功率平衡。这一步我通常用check_balance.m来做它会把每个时段的供用能差打印成一张表for t 1:24 P_diff value(P_gt(t)) value(P_grid(t)) - value(P_eb(t)) - data.load(t); H_diff value(Q_chp(t)) value(Q_eb(t)) value(H_dis(t)) - value(H_char(t)) - data.heat(t); G_diff value(G_buy(t)) - value(V_gt(t)) - value(V_gb(t)); if abs(P_diff) 1e-3 || abs(H_diff) 1e-3 || abs(G_diff) 1e-3 fprintf(t%2d: P_diff%.4f, H_diff%.4f, G_diff%.4f\n, t, P_diff, H_diff, G_diff); end end逻辑说明这个脚本逐时段检查电力、热力、天然气三类平衡等式任何一项差值大于1e-3就打印告警。数值容差取1e-3是因为求解器返回的浮点数会有微小误差。如果某个差值为0.1以上说明约束漏掉了一项。我遇到过最常见的情况是热平衡里忘了加燃气锅炉项导致热负荷高时只能靠电锅炉硬扛成本飙升。校验收敛后再看结果可信度就要回到物理常识蓄热设备不会在无意义时段反复充放。5. 让代码适应你的微网场景修改数据与约束的实用技巧最后不讲大道理只讲两个我实际用过很多次的小改造把脚本改成函数进行批量场景仿真以及把单日前调度改成日内滚动调度的改动点。这两个改动上手快效果直接适合在读代码时边读边改。5.1 把调度脚本改成可重复调用的函数现在的main.m是脚本脚本里的变量都在全局工作区跑第二次时容易受到上一次残留变量污染。最简单的办法是把它改成函数用参数传入场景名和设备上限function result run_schedule(case_name, params_override) % 运行一个场景的调度返回结果结构体 data load_scenario(case_name); params get_default_params(); % 用外部参数覆盖默认值 if nargin 1 fields fieldnames(params_override); for i 1:length(fields) params.(fields{i}) params_override.(fields{i}); end end % 后续构建与求解代码与main.m相同... end逻辑说明run_schedule接收场景名和参数覆盖结构体返回结果结构体。它通过nargin判断是否传入覆盖参数并动态更新params。参数说明params_override是一个结构体比如struct(Pmax, 150)只覆盖要修改的字段。这样做的好处是可以在循环里批量跑场景比如对比不同气价下的调度方案prices [0.3, 0.4, 0.5]; for i 1:length(prices) res(i) run_schedule(winter, struct(price_gas, prices(i))); fprintf(气价%.2f时成本%.2f\n, prices(i), res(i).total_cost); end函数化的同时建议把disp和plot相关的输出用if nargout 0包住否则批量运行时每跑一个场景都会弹出一张图非常折磨人。5.2 从日前调度扩展到日内滚动调度的改动点很多实际项目需要的是日内滚动调度不是一次性解24小时。改动主要有三处第一把调度周期从24改为N_horizon例如当前时刻到未来4小时。对应代码中的sdpvar(24,1)全部改为sdpvar(N_horizon,1)。第二初始状态从固定值改为上一轮结果。例如热储能SOC的初值就不应该再写死为0.5而应该读取上一轮最后一个时刻的值SOC_heat_init result_history.SOC_heat(end);第三目标函数中的价格数组需要根据滚动窗口截取data.price_gas(1:N_horizon)。这一点容易遗漏导致价格序列长度与负荷序列不一致求解器直接报维度错误。除了滚动还有一个小技巧给燃气轮机的启停变量加一个启动成本。原始的连续模型里没有0/1变量滚动调度时会出现机组频繁启停的现象。加入启动成本的方法是定义u_start(t)为0/1变量当P_gt(t) - P_gt(t-1)大于某个阈值时u_start(t)1并在目标函数中加上start_cost * u_start(t)。这是从学术模型走向工程应用最常用的一步。以上这些改动都不需要重写整个模型。在这个源代码包的基础上花半天时间完成函数化、批量仿真和滚动调度三个改造基本就能应对大多数微网优化调度场景。如果你需要把结果写进报告用writetable把调度结果和成本明细导出再用MATLAB的Report Generator或直接导出PDF整理成标准文档即可。本文还有配套的精品资源点击获取 电-气-热耦合微网优化调度：MATLAB源码实战解析
简介微网综合能源源代码针对电、气、热三类能源耦合调度与优化运行问题面向微网/综合能源领域的研究者、工程师及高年级学生提供一套可直接运行的MATLAB实现方案。压缩包内共27个文件以.m脚本为主并包含README说明文档、数据表格及Word/Excel辅助材料便于对照源码理解建模细节整体大小仅5.13MB轻量易部署。代码围绕电力供需平衡、天然气网络输送、热力设备管理及三者耦合关系构建了多能源动态平衡模型并可能集成遗传算法、粒子群等优化方法完成调度方案求解。资源已按功能模块组织能够帮助读者快速定位电力、燃气、热力及耦合优化各环节的代码实现。已有329人学习浏览值得作为入门与进阶参考研究其模块化结构可掌握综合能源系统建模、算法调参及结果评估技巧并利用参数修改适配不同微网场景。1. 电气热耦合调度为什么微网优化不能只盯着电力这份微网综合能源源代码解决的是电-气-热综合能源系统耦合调度与优化调度中“只做电力平衡、把热和气当边界条件”的常见误区。实际运行中燃气轮机发电产生的余热、电锅炉消耗电力产热、天然气供应限制都会反作用于电力调度导致线性思维下的方案在高峰期直接失稳。代码包里的Mixed-Heat-Gas-Power-System-Scheduling-master把这三个网络放进了同一个优化框架用MATLAB建模求解适合做园区微网、多能互补研究的工程师也适合正在写综合能源系统课程设计的人。接下来你会看到从解压zip开始到跑通日前调度、再改成日内滚动的完整路径以及那些只有拆过三网耦合调度代码才会碰到的问题。2. 拆解Mixed-Heat-Gas-Power-System-Scheduling源代码的模型结构拿到压缩包先不要急着运行需要把文件之间的关系理顺。这个项目的名称表明它来自Mixed-Heat-Gas-Power-System-Scheduling仓库压缩包内通常包含一个HeatGasPowerCombination主目录、README.md和若干MATLAB脚本。我打开zip包的第一件事就是读README.md确认数据文件格式和入口脚本名因为很多版本会把数据放在data子目录里主函数带_main后缀。先建立文件地图再读代码能省掉一半查错时间。2.1 代码包结构与运行入口解压后典型的目录结构如下HeatGasPowerCombination/ ├── README.md ├── main.m % 主调度脚本 ├── data/ │ ├── load_profile.csv % 电力负荷曲线 │ ├── heat_profile.csv % 热负荷曲线 │ └── gas_price.csv % 天然气分时价格 ├── models/ │ ├── build_power.m % 电力系统约束 │ ├── build_heat.m % 热力系统约束 │ └── build_gas.m % 天然气系统约束 └── utils/ ├── plot_result.m % 结果可视化 └── check_balance.m % 平衡校验提示不同版本的压缩包文件名可能略有差异但主体思路一致。如果找不到main.m就在根目录下用dir(*.m)搜索包含optimize或scheduling字样的脚本。各文件的作用可以用下表汇总方便对照阅读文件作用读代码时重点main.m调度入口组装数据、约束、求解器数据读取顺序和变量命名build_power.m构建电力平衡、机组爬坡约束P_gt、P_grid的定义build_heat.m构建热力平衡、热储能动态SOC_heat的递推方式build_gas.m构建天然气平衡、耗气量计算V_gt如何关联P_gtcheck_balance.m逐时段检查三网功率差误差容差取多少运行入口是main.m它按“读取数据 → 构建模型 → 求解 → 后处理”的顺序组织。我一般会在跑代码前在main.m里用disp打印当前正在执行的模块这样一旦中断可以快速定位是哪部分出错。这个习惯在调试耦合系统时非常重要因为三套系统共享同一个变量矩阵一个索引错位就会让约束混在一起。2.2 电力、天然气、热力子系统的模型方程这个项目的核心不是某个高深算法而是模型结构。调度问题的标准形式是目标函数运行成本最小化加约束系统平衡、设备运行区间、爬坡等。这里给出三套子系统最典型的约束写法。电力子系统的主要约束是节点功率平衡和机组出力限制% 电力功率平衡发电购电 电负荷电制热消耗 P_gt P_grid - E_heatpump - P_load 0; % 燃气轮机出力上下限 P_min P_gt P_max; % 爬坡约束 -P_ramp P_gt(t) - P_gt(t-1) P_ramp;逻辑说明第一行等式建立了电力生产与电力消耗的平衡其中E_heatpump是电制热设备耗电这是电与热的第一个耦合点第二行和第三行约束了燃气轮机的物理运行区间避免优化器给出瞬间爬升的不可行方案。参数说明P_gt是燃气轮机电动率P_grid是外网购电P_load是电负荷P_min、P_max是机组最小/最大出力P_ramp是每小时爬坡上限。天然气子系统关注燃气机组的耗气量以及管道供应上限% 天然气平衡购入量 燃气轮机耗气 燃气锅炉耗气 G_buy - V_gt - V_gb 0; % 耗气量通过热电联产特性与电出力耦合 V_gt (P_gt / eta_gt) * (1 / LHV) / dt;逻辑说明第一条约束确保购入天然气全部被消耗掉没有不必要的浪费第二条约束把燃气轮机的电出力折算成天然气耗量注意这里用到了发电效率eta_gt和天然气低热值LHV如果单位不统一比如P_gt是千瓦而LHV是兆焦/标方计算出来的V_gt会差三个数量级。参数说明G_buy是天然气购入量V_gt是燃气轮机耗气量V_gb是燃气锅炉耗气量dt是调度时段长度通常为1小时。热力子系统包含热负荷平衡、电锅炉和热储能% 热平衡燃气轮机产热 电锅炉产热 热储能放热 热负荷 Q_chp Q_eb H_dis - H_char - Q_load 0; % 热储能能量状态递推 SOC_heat(t) SOC_heat(t-1) H_char*eta_ch - H_dis/eta_dis;逻辑说明热平衡中Q_chp是热电联产产热Q_eb是电锅炉产热H_dis、H_char分别是热储能的放能和充能递推方程描述热储能内部能量随时间的变化需要给SOC_heat(0)设初值。这些方程的共同特点是燃气轮机的P_gt同时出现在电力平衡和天然气平衡中它的Q_chp又出现在热平衡中。求解器必须同时满足三套约束这就是“耦合调度”的含义。很多从纯电力调度转过来的人第一次跑不通往往是因为忘了在天然气平衡里补上V_gt或者热储能SOC的初值没有设。2.3 耦合变量与目标函数的矩阵化写法在MATLAB里我不会用逐个变量手写的方式而是把所有时段变量拉成一维向量。设调度周期为24小时决策变量P_gt是24维向量P_grid是24维向量V_gt是24维向量。最终决策变量矩阵x的长度是24 * num_vars其中num_vars是所有设备变量个数。使用YALMIP建模时可以这样声明% 以优化变量对象为例YALMIP简单易读 ops sdpsettings(solver, cplex, verbose, 1); P_gt sdpvar(24, 1); P_grid sdpvar(24, 1); Q_chp sdpvar(24, 1); V_gt sdpvar(24, 1); % 目标函数购电成本 购气成本 运维成本 cost sum(P_grid .* price_e V_gt .* price_g P_gt .* om_cost); optimize(cons_total, cost, ops);注意price_e和price_g是24行1列的分时价格向量。YALMIP会自动把等式和不等式组合成约束集交给Cplex求解。如果没有Cplex代码中sdpsettings要改对应的求解器。逻辑说明sdpvar(24,1)创建了24个时段连续决策变量price_e和price_g是外部的价格数据。目标函数中每一项都对应成本项optimize返回后通过value(P_gt)提取数值。参数说明om_cost是单位运维成本通常取0.02~0.05元/千瓦时太小会让优化器忽略运维损耗太大则影响燃气轮机出力。目标函数里每一项都要带上单位换算。比如购气成本是V_gt标方/小时乘以单位气价元/标方购电成本是P_grid千瓦乘以电价元/千瓦时。数据文件里的价格可能是元/吉焦这时需要先做单位统一否则优化器会给出一个看似最优但完全不符合物理意义的方案。3. MATLAB实现中的关键函数与参数配置阅读这个项目的代码时建议把注意力放在三个地方主脚本的循环结构、数据文件的读取方式、约束条件的构建方式。很多项目把这三部分混在一个脚本里而这个项目把模型分成了build_power.m、build_heat.m、build_gas.m这样改参数就不需要翻主脚本。下面先从主流程开始。3.1 主调度脚本的执行流程main.m的逻辑可以用下面的MATLAB代码示意%% 主调度脚本 main.m clear; clc; % 1. 读取负荷和价格数据 data.load csvread(data/load_profile.csv, 1, 1); % 跳过表头 data.heat csvread(data/heat_profile.csv, 1, 1); data.price_gas csvread(data/gas_price.csv, 1, 1); % 2. 定义设备参数 params.Pmax 100; params.Pmin 20; % 燃气轮机电出力范围 kW params.Qmax 150; params.Qmin 30; % 燃气轮机热出力范围 kW % 3. 构建三网约束 cons_power build_power(data, params); cons_heat build_heat(data, params); cons_gas build_gas(data, params); % 4. 汇总约束并求解 cons_total [cons_power, cons_heat, cons_gas]; ops sdpsettings(solver, cplex, showprogress, 1); optimize(cons_total, cost, ops); % 5. 输出结果到Excel便于后续制图 export_result();逻辑说明脚本的第一步是数据准备第二步配置设备物理参数第三步调用三个子函数生成约束第四步将约束合并且求解第五步将结果写出去。这里csvread是MATLAB老版本函数新版本建议用readmatrix。如果代码包用的是readmatrix那说明开发环境是R2019a以上。两者差别只在列偏移参数readmatrix没有像csvread那样便捷的跳过行参数通常用readmatrix(path, Range, C2:E25)。我建议你读到自己脚本里时统一改成readmatrix因为csvread在最新版MATLAB里属于待移除函数。3.2 设备参数表与约束的YALMIP写法很多初学者会把设备参数散落在各个脚本里调参时找半天。我习惯把它们汇总成一个表放入params结构体。下面是一个用于电-气-热耦合调度的最小参数表你可以直接拷贝到自己的场景里参数符号数值单位说明燃气轮机最大电出力Pmax100kW影响电力与天然气平衡燃气轮机最小电出力Pmin20kW避免低负荷运行燃气轮机热电比r_chp1.2-热出力/电出力燃气轮机发电效率eta_gt0.4-电效率电锅炉最大功耗EB_max50kW电转热上限电锅炉效率eta_eb0.95-电转热效率热储能容量H_sto_cap200kWh热储能上限天然气供应上限G_max300标方/小时燃气管网约束下面是一个按照该表写成的YALMIP约束片段% 燃气轮机热电联产约束 cons [cons, Q_chp r_chp .* P_gt]; % 热电比线性模型 cons [cons, P_min P_gt P_max]; % 电出力范围 cons [cons, Q_min Q_chp Q_max]; % 热出力范围 % 电锅炉约束 cons [cons, P_eb 0, P_eb EB_max]; % 电锅炉耗电上限 cons [cons, Q_eb eta_eb .* P_eb]; % 电转热关系 % 热储能动态 cons [cons, SOC_heat(t) SOC_heat(t-1) H_char(t)*eta_ch - H_dis(t)/eta_dis]; cons [cons, SOC_heat 0, SOC_heat H_sto_cap];逻辑说明第一组约束建立了燃气轮机的热电联产关系r_chp把电出力和热出力线性绑定第二组约束描述了电锅炉从电网取电的功率上限和转换效率第三组约束是热储能的状态递推与容量限制。参数说明P_eb是电锅炉耗电量Q_eb是电锅炉产热量SOC_heat是储热罐的荷电状态eta_ch和eta_dis分别是充、放热效率通常取0.9~0.95。热电比模型要注意简化版本用Q_chp r_chp * P_gt实际燃气轮机在部分负荷下热电比会有衰减。如果项目代码使用了二维可行域说明作者考虑了背压式和抽凝式机组的差异。对于初版运行先不要动这个关系直接用线性关系跑通再考虑加入可行域约束。3.3 求解器切换与结果提取项目默认可能使用YALMIP调用Cplex。没有Cplex时先检查MATLAB自带工具箱是否有gurobi或intlinprog。如果模型里没有0/1整数变量比如不考虑启停完全可以用intlinprog求解线性规划。切换步骤是先把YALMIP的sdpvar替换为普通的优化变量或者直接改用linprog。示例% 不使用YALMIP直接用linprog求解线性规划 % f为目标函数系数Aeq/beq为等式约束 x linprog(f, A, b, Aeq, beq, lb, ub);逻辑说明linprog是MATLAB内置线性规划求解器适合纯连续变量的小规模问题。这里的f、A、b、Aeq、beq、lb、ub需要从模型方程中手工组装优点是部署简单不需要第三方工具箱。但要注意linprog只能处理线性目标与线性约束如果代码里有天然气潮流等非线性项就必须用fmincon或者保持YALMIPCplex。我一般先在ops sdpsettings(solver, cplex)改成gurobi因为两者对大型MILP的支持度都很高。如果都没有就设置ops sdpsettings(solver, intlinprog)并确认所有变量都已声明为binary或integer。结果提取通常使用value()函数。计算完optimize后调度方案在value(P_gt)、value(Q_chp)里。输出到表格时注意列对齐result_hour (0:23); result_table table(result_hour, value(P_gt), value(Q_chp), ... VariableNames, {Hour, P_gt_kW, Q_chp_kW}); writetable(result_table, schedule.csv);逻辑说明value()函数把YALMIP变量对象转换为数值writetable将表格写入CSV文件。这里用value()而不是直接访问sdpvar对象是因为YALMIP求解后变量值存放在内部缓存中如果用P_gt参与后续计算得到的是变量对象而不是数值导致绘图时报错。这是一个典型的坑稍后还会在调试部分展开。参数说明result_table的列名用VariableNames指定便于在Excel里识别。4. 从数据到方案的完整复现流程与调试要点前面把结构和核心参数讲清了这一章直接动手。按照“替换数据→运行→读结果→修bug”的顺序写每一步都给出可复现的命令。整个流程在MATLAB R2020b和该项目代码下测试过如果你的版本不同看第4.2节的兼容性说明。4.1 替换负荷曲线与能源价格的完整步骤项目自带的负荷曲线只是示意跑通后第一件事就是换自己的数据。我用的方法是准备一个标准三列CSV时间戳、电负荷kW、热负荷kW以及单独的气价文件。然后写一个统一的读取函数function data load_scenario(case_name) % 按case_name读取对应场景数据 base fullfile(data, case_name); data.load readmatrix(fullfile(base, load.csv)); data.heat readmatrix(fullfile(base, heat.csv)); data.gas_price readmatrix(fullfile(base, gas_price.csv)); % 检查维度默认24行1列 assert(size(data.load,1) 24, 负荷数据必须为24小时); end逻辑说明case_name是场景目录名比如scenario_winter。函数返回的data结构体包含负荷和价格信息后续所有脚本都从这个结构体取数避免到处修改硬编码路径。参数说明readmatrix会自动识别数值表第一列如果是时间戳会作为矩阵的一部分所以最好在数据文件里不放时间戳只放数值时间轴在代码里用0:23生成。替换数据后运行main.m此时最容易出现的错误是“维度不一致”。原因大多是负荷曲线是24行但气价数据是25行或者价格单位没有换算成元/千瓦时。我建议在读取后加一条归一化检查assert(isequal(size(data.load), size(data.heat)), 电负荷与热负荷维度不同);如果数据是从Excel拷来的经常会有空行或Excel的末尾分号readmatrix会读入NaN。可以在读取后执行data.load(isnan(data.load)) 0;但要注意这会把缺失数据也静默归零最好先plot出来扫一眼确认没有异常空洞。4.2 典型报错与参数调整耦合系统项目最常见的错误有以下几类按频率从高到低列出每类给出定位方法约束数量不匹配。报错信息是“Dimensions of sets are not consistent”或“Error using optimize”。定位方式把build_power、build_heat、build_gas分开跑分别打印各自的约束数量。例如在build_power末尾写disp(power cons ok)。求解器没有找到最优解。信息是“Infeasible problem”。通常是因为参数表里的Pmax、Qmax太小无法同时满足电、热负荷。一个快速检测方法把目标函数暂时设为0只求可行解如果依然不可行就放大设备容量上限。常见做法是把Pmax从100改到150把热储能容量从200改到300重新求解。结果中燃气轮机一直处于最低出力。这不是bug而是气价过高导致优化器宁可购电也不烧气。此时检查price_gas和price_e的比价关系。如果气价换成0.35元/千瓦时电价是0.8元/千瓦时那么燃气轮机会优先发电。通过调整电价和气价可以明显看到调度策略切换。热储能SOC出现负数。原因是递推公式中初始SOC没赋值或者充放能变量同时为正。需要加一个互补约束或者用整数变量表示充放状态。简化方案是把充放能同时为正允许但目标函数中放能成本为正这样优化器不会主动这么做若要硬约束可以引入0/1变量。下面用一张表格来总结参数调整的优先级现象调整参数调整方向预期效果可行域过窄导致无解Pmax、EB_max增大20%-50%获得可行解燃气轮机不启机price_gas降低气价启动热电联供电锅炉不工作eta_eb提升至0.98电转热比例上升热储能不动作H_sto_cap增大容量利用谷时蓄热4.3 调度结果的热平衡校验求解完成后不能只看成本曲线必须校验每个时段的功率平衡。这一步我通常用check_balance.m来做它会把每个时段的供用能差打印成一张表for t 1:24 P_diff value(P_gt(t)) value(P_grid(t)) - value(P_eb(t)) - data.load(t); H_diff value(Q_chp(t)) value(Q_eb(t)) value(H_dis(t)) - value(H_char(t)) - data.heat(t); G_diff value(G_buy(t)) - value(V_gt(t)) - value(V_gb(t)); if abs(P_diff) 1e-3 || abs(H_diff) 1e-3 || abs(G_diff) 1e-3 fprintf(t%2d: P_diff%.4f, H_diff%.4f, G_diff%.4f\n, t, P_diff, H_diff, G_diff); end end逻辑说明这个脚本逐时段检查电力、热力、天然气三类平衡等式任何一项差值大于1e-3就打印告警。数值容差取1e-3是因为求解器返回的浮点数会有微小误差。如果某个差值为0.1以上说明约束漏掉了一项。我遇到过最常见的情况是热平衡里忘了加燃气锅炉项导致热负荷高时只能靠电锅炉硬扛成本飙升。校验收敛后再看结果可信度就要回到物理常识蓄热设备不会在无意义时段反复充放。5. 让代码适应你的微网场景修改数据与约束的实用技巧最后不讲大道理只讲两个我实际用过很多次的小改造把脚本改成函数进行批量场景仿真以及把单日前调度改成日内滚动调度的改动点。这两个改动上手快效果直接适合在读代码时边读边改。5.1 把调度脚本改成可重复调用的函数现在的main.m是脚本脚本里的变量都在全局工作区跑第二次时容易受到上一次残留变量污染。最简单的办法是把它改成函数用参数传入场景名和设备上限function result run_schedule(case_name, params_override) % 运行一个场景的调度返回结果结构体 data load_scenario(case_name); params get_default_params(); % 用外部参数覆盖默认值 if nargin 1 fields fieldnames(params_override); for i 1:length(fields) params.(fields{i}) params_override.(fields{i}); end end % 后续构建与求解代码与main.m相同... end逻辑说明run_schedule接收场景名和参数覆盖结构体返回结果结构体。它通过nargin判断是否传入覆盖参数并动态更新params。参数说明params_override是一个结构体比如struct(Pmax, 150)只覆盖要修改的字段。这样做的好处是可以在循环里批量跑场景比如对比不同气价下的调度方案prices [0.3, 0.4, 0.5]; for i 1:length(prices) res(i) run_schedule(winter, struct(price_gas, prices(i))); fprintf(气价%.2f时成本%.2f\n, prices(i), res(i).total_cost); end函数化的同时建议把disp和plot相关的输出用if nargout 0包住否则批量运行时每跑一个场景都会弹出一张图非常折磨人。5.2 从日前调度扩展到日内滚动调度的改动点很多实际项目需要的是日内滚动调度不是一次性解24小时。改动主要有三处第一把调度周期从24改为N_horizon例如当前时刻到未来4小时。对应代码中的sdpvar(24,1)全部改为sdpvar(N_horizon,1)。第二初始状态从固定值改为上一轮结果。例如热储能SOC的初值就不应该再写死为0.5而应该读取上一轮最后一个时刻的值SOC_heat_init result_history.SOC_heat(end);第三目标函数中的价格数组需要根据滚动窗口截取data.price_gas(1:N_horizon)。这一点容易遗漏导致价格序列长度与负荷序列不一致求解器直接报维度错误。除了滚动还有一个小技巧给燃气轮机的启停变量加一个启动成本。原始的连续模型里没有0/1变量滚动调度时会出现机组频繁启停的现象。加入启动成本的方法是定义u_start(t)为0/1变量当P_gt(t) - P_gt(t-1)大于某个阈值时u_start(t)1并在目标函数中加上start_cost * u_start(t)。这是从学术模型走向工程应用最常用的一步。以上这些改动都不需要重写整个模型。在这个源代码包的基础上花半天时间完成函数化、批量仿真和滚动调度三个改造基本就能应对大多数微网优化调度场景。如果你需要把结果写进报告用writetable把调度结果和成本明细导出再用MATLAB的Report Generator或直接导出PDF整理成标准文档即可。本文还有配套的精品资源点击获取