1. 项目背景与整体设计思路
1.1 为什么选综合能源系统这个方向
这几年做能源方向的人应该都有明显感觉，单纯搞光伏、风电单点接入的研究越来越吃力，真正被工程界和学术界反复提及的是“综合能源系统”——把电、热、气、储多个能源环节放在一个框架里统一调度。原因也简单：单一能源在时空上都有明显短板，光伏夜间不出力，风电波动大，电锅炉直接供热又贵，而把这些设备通过能量耦合关系串起来之后，可以实现“电转热”“热电联产”“储能移峰填谷”等互补效应。
低碳运行优化调度，说白了就是在一组设备组合和负荷需求的前提下，回答两个问题：每个时刻每台设备出力多少？系统整体成本和碳排放如何达到平衡？如果只做简单规则调度，比如“光伏优先、缺电买电”，那程序两小时就能写完，但没有优化空间。真正有价值、也最能体现综合能源系统优势的，是把设备运行域、储能状态转移、分时电价、碳交易成本全部写进一个数学优化模型里，交给求解器去找最优解。
1.2 这套Matlab程序的核心构成
我这里分享的这套Matlab程序，覆盖了光伏、风电、热电联产机组（CHP）、燃气锅炉、电储能五类核心设备，模型采用典型的日前调度框架，时间尺度为24小时、步长1小时。目标函数同时考虑运行成本和碳交易成本，属于典型的多目标转单目标的加权处理方式，便于使用混合整数线性规划（MILP）直接求解。程序基于YALMIP建模，底层调用CPLEX求解器。
为什么用YALMIP而不是手写线性规划矩阵？这是我的真实体会：手写A、b、Aeq、beq矩阵做小规模模型还行，设备一多、约束一杂，很容易行列对不上出错。YALMIP提供了类似数学公式的自然建模方式，变量定义直接用sdpvar、约束条件用集合存起来，代码可读性高，后期调试也方便。而且YALMIP可以把同一个模型无缝切换求解器——测试时用内置linprog，正式跑用CPLEX或Gurobi，写论文对比实验的时候非常省事。
这套程序适合谁？两类人。第一类是刚入门的硕士生，要做综合能源系统优化调度方向，需要一个能跑通全流程、可以改参数和扩展结构的基线代码。第二类是偏工程的从业者，想快速评估某个园区或楼宇的新能源配置方案，用这套框架做初步的容量配置和运行策略预估。你可以拿它当“骨架”，往里加氢储能、加电转气、加需求响应都行。
2. 综合能源系统拓扑与设备建模
2.1 系统拓扑结构
从拓扑结构上看，这套程序采用典型的母线式结构：电母线连接光伏、风电、CHP电出力、电储能、电网购电和电负荷；热母线连接CHP热出力、燃气锅炉和热负荷。两个母线通过CHP机组耦合在一起，形成“以热定电”或“以电定热”的调节能力，这是综合能源系统区别于独立电、热系统的核心特征。
我实际建模时，电母线和热母线分别列功率平衡方程，这是整个模型最底层的骨架。母线平衡方程不满足，其他约束再精细也没意义。这里要留意一个常见的坑：电平衡里如果有储能，充放电符号必须统一，否则求出来的结果可能在某个时刻“凭空多出”一部分功率，后期检查残差才追悔莫及。
我用的设备连接关系：
- 光伏和风电作为不可调度电源，出力按预测曲线输入，在优化模型中作为已知参数。
- CHP机组是热电耦合设备，需要单独描述电出力和热出力之间的关系。
- 燃气锅炉作为热母线的补充热源，负责在CHP供热不足时补充热量。
- 电储能连接电母线，通过荷电状态（SOC）把不同时刻的电量耦合在一起，实现能量的时间平移。
- 电网作为系统与外部的交换节点，允许购电和售电。
2.2 各设备数学模型的要点
光伏和风电模型相对简单，核心就是一个功率上限约束：
P_pv(t) <= P_pv_forecast(t) P_wt(t) <= P_wt_forecast(t)
再配合各自的爬坡约束。实际工程中光伏逆变器可以限功率运行，但日前调度通常直接取预测值，必要时对预测曲线做向下修正，留出备用空间。
CHP机组建模是关键，也是最容易让新手懵的部分。我用的是简化的线性热电比模型，设额定电功率为300kW，最小技术出力取30%额定值，热电比固定为1.3。数学表达就是：
H_chp(t) = K_chp * P_chp(t) P_chp_min <= P_chp(t) <= P_chp_max
为什么用定热电比？因为只要热电比是固定常数，CHP模型就是线性的，目标函数和约束全部线性，求解器的负担会小很多。如果想更精细，可以用抽凝式机组的可行域多边形描述——CHP出力在一个由多个顶点围成的凸多边形内，这个我在后面“功能扩展”部分再展开。对于多数园区级系统，定热电比精度够用，工程上很多论文都这么建模。
燃气锅炉模型就是典型能量转换：
H_boiler(t) = η_boiler * F_boiler(t) 0 <= H_boiler(t) <= H_boiler_max
其中F_boiler(t)表示燃料输入，η_boiler是锅炉热效率，我这里取0.9。注意燃气消耗会直接进入成本函数和碳排放计算，不要只盯着热出力，否则目标函数里费用会算少。
电储能是另一个重点。储能的特点是跨时段耦合，必须用SOC状态变量把相邻时刻串起来：
SOC(t+1) = SOC(t) + η_chg * P_chg(t) - P_dis(t) / η_dis
除了容量约束和充放电功率上限，还必须加充放电互斥约束——同一时刻不能既充电又放电，否则模型会“钻空子”，利用同时充放电做无意义的能量循环，白白抬高成本。互斥用二进制变量实现：
P_chg(t) <= P_chg_max * u_chg(t) P_dis(t) <= P_dis_max * u_dis(t) u_chg(t) + u_dis(t) <= 1
另外要加周期一致性约束，SOC(1) = SOC(25)。这是很多人容易漏掉的：如果不约束调度周期结束后的SOC回到初始值，优化器为了压最后几个小时的买电成本，会把储能放空到死，第二天没法继续正常运行。
2.3 分时电价与碳排放模型
购电成本直接采用分时电价，这个不复杂。我这里设置了峰、平、谷三个时段，电价分别为1.1、0.68、0.35元/kWh，谷时段是从23点到次日7点，峰时段取10点到15点、18点到21点。电价数据建议放到独立的excel或mat文件里，改起来方便，不用进代码逐行翻。
碳排放模型是“双碳”背景下必须考虑的部分。我这里用碳交易机制：系统总碳排放量减去免费配额后，超出部分需要按碳价购买，未超标部分理论上可以出售。总排放由两部分构成：
E_total = E_grid + E_gas
其中E_grid = ∑(P_buy(t) * EF_grid)是外购电力的间接排放，EF_grid取电网平均排放因子，我的测试场景取0.581 tCO2/MWh。E_gas是天然气燃烧产生的直接排放，CHP和燃气锅炉的燃气消耗量依据燃料热值乘以排放因子计算。免费配额按照系统供热和供电量乘以对应基准线确定，这个基准线可以灵活调整。
从优化角度理解：引入碳交易成本后，系统会在“买便宜但高碳的电”和“用更贵的天然气自己发电但低碳排”之间做权衡。碳价越高，系统越倾向于降低外购电量，增加本地清洁能源利用和CHP发电。这个机制在参数灵敏度分析里特别有意思，后面我会展示实测数据。
3. 优化调度模型的完整构建
3.1 决策变量和目标函数定义
通常把一天分成24个时段的决策过程。决策变量包括：CHP的电出力、热出力、启停状态；燃气锅炉的热出力；储能每个时段的充放电功率、互斥状态、SOC；外购电量和售电量。光伏和风电出力在日前调度里作为已知参数入模，不做优化变量——它们的预测值就是计划值，预测误差在实时调度修正阶段处理。
目标函数这块，我在“低碳运行”要求下设计成运行成本与碳交易成本之和：
min F = ∑[ C_buy(t)P_buy(t) + C_gas_chp(t) + C_gas_boiler(t) + C_om_chpP_chp(t) + C_om_boiler*H_boiler(t) + C_carbon * (E_total - E_quota) ]
说明几点：
- C_buy(t)P_buy(t)表示购电成本，售电收入写成负数加进去，例如C_sellP_sell(t)。
- C_gas_chp和C_gas_boiler是燃料成本，由天然气价格乘以热耗折算。CHP的电效率取0.35、热效率取0.45，燃气锅炉效率0.9，天然气价格按2.5元/m³，折算后每kWh热量的气价大约是0.25元。
- 设备运维成本C_om，取一个很小的常数，主要目的是在多个等效解之间选一个损耗更小的方案。
- 最后一项碳交易成本是“加权”进目标函数，而不是硬约束，好处是可以通过调节碳价参数，观察系统调度策略对碳价变化的响应曲线，这在后续灵敏度分析里非常有价值。
3.2 约束条件清单
约束条件可以整理成下面这组，全部按24个时段展开：
电力平衡约束： P_pv(t) + P_wt(t) + P_chp(t) + P_dis(t) + P_buy(t) = P_load(t) + P_chg(t) + P_sell(t)
热力平衡约束： H_chp(t) + H_boiler(t) = H_load(t)
CHP运行约束（含最小启停时间、出力上下限、热电耦合关系）： P_chp_minu_chp(t) <= P_chp(t) <= P_chp_maxu_chp(t) H_chp(t) = K_chp * P_chp(t)
储能系统约束（含容量、功率、互斥、初始和终值SOC）： SOC_min <= SOC(t) <= SOC_max 0 <= P_chg(t) <= P_chg_max * u_chg(t) 0 <= P_dis(t) <= P_dis_max * u_dis(t) u_chg(t) + u_dis(t) <= 1 SOC(1) = SOC_init SOC(t+1) = SOC(t) + η_chg*P_chg(t) - P_dis(t)/η_dis
外购电约束： 0 <= P_buy(t) <= P_buy_max
这些约束看起来多，但在YALMIP里写起来就是逐个列方程。需要注意：全部约束要放到一个cell数组里，最后让optimize一次性求解，不要分散成多个optimize调用，否则就变成序贯决策而非全局优化了。
3.3 YALMIP求解代码结构
核心求解部分的核心代码我简化后贴出来，已经去掉注释中与实际设备参数无关的内容：
%% 决策变量定义 P_pv = sdpvar(1, 24); P_wt = sdpvar(1, 24); P_chp = sdpvar(1, 24); H_chp = sdpvar(1, 24); u_chp = binvar(1, 24); H_boil = sdpvar(1, 24); P_chg = sdpvar(1, 24); P_dis = sdpvar(1, 24); u_chg = binvar(1, 24); u_dis = binvar(1, 24); SOC = sdpvar(1, 25); P_buy = sdpvar(1, 24); P_sell = sdpvar(1, 24);光伏和风电视为已知参数，直接初始化为预测曲线，不给它们设置决策自由度：
P_pv = PV_forecast; % 1x24已知向量 P_wt = WT_forecast; % 1x24已知向量约束搭建：
con = []; %% 电平衡 con = [con, P_pv + P_wt + P_chp + P_dis + P_buy ... == P_load + P_chg + P_sell]; %% 热平衡 con = [con, H_chp + H_boil == H_load]; %% CHP运行 con = [con, P_chp_min*u_chp <= P_chp <= P_chp_max*u_chp]; con = [con, H_chp == CHP_K * P_chp]; %% 锅炉 con = [con, 0 <= H_boil <= H_boil_max]; %% 储能 con = [con, SOC_min <= SOC <= SOC_max]; con = [con, 0 <= P_chg <= P_chg_max*u_chg]; con = [con, 0 <= P_dis <= P_dis_max*u_dis]; con = [con, u_chg + u_dis <= 1]; con = [con, SOC(2:25) == SOC(1:24) ... + eta_chg*P_chg - P_dis/eta_dis]; con = [con, SOC(1) == SOC_init, SOC(25) == SOC_init]; %% 电网交换 con = [con, 0 <= P_buy <= P_buy_max]; con = [con, P_sell >= 0];目标函数和求解：
%% 目标函数 Carbon_cost = C_carbon * (sum(P_buy)*EF_grid ... + sum(H_boil)/eta_boil * EF_gas... + sum(P_chp)*K_chp/eta_chp_heat * EF_gas... - E_quota); F = sum(Price_buy.*P_buy - Price_sell.*P_sell) ... + gas_price_chp * K_chp/eta_chp_heat * sum(P_chp) ... + gas_price_boil * sum(H_boil)/eta_boil ... + sum(om_chp*P_chp + om_boil*H_boil) ... + Carbon_cost; %% 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); result = optimize(con, F, ops);求解完取结果：
P_chp_opt = value(P_chp); H_chp_opt = value(H_chp); SOC_opt = value(SOC);这套代码结构经过多次调整，现在已经比较稳定。新手最容易在这里犯的错误有两个：一是变量维度不一致，比如SOC定义成1×25，但约束里SOC(2:25)和SOC(1:24)不匹配；二是二进制变量没有对应u_chp的初值，导致初始解搜索慢。这些细节我在下文“问题排查”部分会单独讲。
4. 典型日仿真结果与策略分析
4.1 典型日基础数据
为了让结果有参照，我构造了一个典型冬季日场景，测试对象是一个包含中等规模商业负荷和居民采暖负荷的园区。光伏装机200kW，风电装机150kW，CHP额定电功率300kW、额定热功率390kW，燃气锅炉额定热功率400kW，电储能600kW·h、功率限制150kW。电负荷峰值大约650kW，热负荷峰值大约420kW。
光伏预测曲线是典型冬日的倒U形，峰值出现在12点半前后，大约156kW；风电因为是冬季北风天，整体出力偏高，在半夜和傍晚有两次峰值，白天略低。这些预测曲线我直接放在data.xlsx里，程序启动时用readtable读入。
4.2 优化调度结果解读
跑完优化后，先看热母线平衡结果。系统在这个典型日里的热负荷约420kW峰值。优化结果里，CHP机组在白天电负荷高峰期全出力运行，热出力达到390kW上限，基本覆盖大部分热负荷；早晚热负荷高峰时，CHP热出力不足的部分由燃气锅炉补足。这说明调度策略的总体倾向是“优先用CHP同步满足电和热需求”，因为CHP的天然气在发完电之后再产热，综合能源利用效率比燃气锅炉单独供热高。
电母线的情况更有意思。白天10点到15点，分时电价处于高峰段，光伏出力也充足，优化结果出现一个明显特征：光伏和风电出力全部消纳，CHP维持在较高出力，储能从上午开始充电，等到晚高峰18点到21点电价最高时放电，打折购电成本。谷时段半夜则大量购电给储能充电，这符合套利逻辑。
我调低碳价到50元/吨时，系统策略偏向谷时段买电、减少本地气电机组出力；把碳价抬高到300元/吨后，外购电量明显下降，CHP和风电出力占比上升，热负荷更多由CHP承担，燃气锅炉补充量减少。原因很简单：外购电的间接排放因子较高，碳价上升后买电的隐性成本变大，系统用天然气替代一部分电网电反而划算。这个现象在做政策评估、碳价敏感性分析时特别有用。
4.3 与人工规则调度对比
为了衡量优化调度的收益，我还写了一个简单的人工规则调度作为对照组：光伏优先、风电优先，缺电就买电，CHP按“以热定电”方式运行，储能固定两个时段充放。对比结果差距非常明显：
- 常规规则调度的日运行成本约12400元。
- 优化调度的日运行成本约10800元，下降约13%。
- 碳排放量从8.4吨降至7.1吨，下降约15%。
- 储能利用率提升明显，规则调度只利用了谷段和峰段各3小时，优化调度把储能作为全天的缓冲手段，充放电次数和深度都更合理。
这不是一个特别惊艳的数字，但考虑到模型里并没有加入激进的节能技术，单靠优化调度把成本和碳排放同时压下来十几个百分点，已经足以说明综合能源系统“调度优化”的价值不依赖硬件升级，纯粹靠运行策略就能产生收益。对于已有系统做节能改造的场景，这种软性优化往往是性价比最高的一步。
4.4 碳价灵敏度分析
碳交易成本在目标函数里是线性的，所以碳价变化对结果的影响本质上是一个参数敏感性分析。我在碳价从0到500元/吨之间取了9个点，把每个点对应的外购电量和CHP发电量提取出来画成曲线，可以看到明显的分段趋势：碳价低于100元/吨时，调度策略几乎不改变，因为此时碳成本在总成本中占比很低；碳价超过150元/吨后，外购电量开始下降，CHP出力上升；超过300元/吨后，变化趋缓，说明系统已经逼近了设备物理上限，再高的碳价也无法进一步改变策略，只能通过扩容改造来继续降碳。
这个分段特征很重要。如果你是在做碳价机制设计，看到曲线拐点就能知道当前系统的碳价“敏感区间”在哪里；如果你在写论文，这个灵敏度图几乎是必放的内容，审稿人一眼就能看出模型是否合理。
5. Matlab实现中的常见问题与排查经验
5.1 求解器相关坑
YALMIP调用CPLEX时最常遇到的问题就是许可证或版本兼容。我用的环境是Matlab R2021b加YALMIP R20210430加IBM CPLEX 12.10。换了新版本Matlab后，YALMIP偶尔出现无法识别CPLEX的情况，这时先运行yalmiptest，看输出中CPLEX是否正常显示。如果显示misconfigured，多半是路径没加对，或者CPLEX版本与Matlab版本不匹配。
另一个常见情况是机器上装了多个求解器，YALMIP默认选了非预期的那个。我在代码里显式指定求解器（solver, 'cplex'），避免这个问题。如果学生机器没有CPLEX license，还可以把求解器改成gurobi，模型代码不用动，只改sdpsettings里的solver参数。这是YALMIP框架的最大好处。
5.2 不可行问题排查
遇到Infeasible problem是家常便饭，不要慌，我有一套标准排查流程。第一步先把对偶松弛变量打开，用optimize的第三个返回值查看约束冲突情况；第二步逐步注释掉约束，用二分法定位是哪一组约束导致的不可行。最常见的是储能SOC周期约束和充放电互斥约束写反了，或者CHP的热电比系数与热平衡约束矛盾——比如额定电功率下的热出力已经超过热负荷上限，还要强行满足热平衡，模型就会无解。
还有一个隐蔽问题：SOC变量在目标函数中没有任何成本项，它是纯状态变量。如果SOC的上下限设置得过窄（比如0.95到0.95），虽然约束本身没问题，但储能完全丧失调节能力，优化器只能在极小可行域里找解，数值上表现为求解时间暴增、解质量差。我一般在没特殊需求时把SOC范围放宽到0.1~0.9，给调度足够的自由度。
5.3 代码运行很慢的问题
24时段、大约200个变量（其中24个二进制）的MILP模型对CPLEX来说是小菜一碟，通常一两秒就解完了。如果你觉得慢，要检查是不是存在大M约束把模型搞“硬”了——比如用了过大的M值，导致求解器数值病态。另一个调优方向是给sdpsettings加参数：
ops = sdpsettings('solver', 'cplex', 'verbose', 2, ... 'cplex.mip.tolerances.mipgap', 0.01, ... 'cplex.mip.tolerances.absmipgap', 0.5);设置一个1%的MIP gap，可以在损失极少量精度的情况下大幅加速求解。对工程应用来说，99%的优化精度和100%的最优解差距完全可以接受。
如果真的要把模型扩展到96个时段（15分钟粒度），变量数量会翻四倍，约束翻四倍，但求解时间可能不止翻四倍。这时建议减少不必要的二进制变量，例如灯具、小型燃气设备可以用连续变量加“软约束”近似，或者把全天调度分解成几个时段的分段优化，牺牲一点全局最优性换取求解速度。
5.4 数据输入与可视化技巧
我习惯把所有场景数据统一放在Excel里，包括负荷曲线、风光预测出力、分时电价、设备参数，程序入口只留一个参数配置脚本。这个习惯看上去多了一个解析Excel的步骤，但好处是后期换数据时完全不用碰模型代码，只要把Excel里的表替换掉就能跑新场景。对需要反复做算例对比的论文党，这个设计省下的时间不可估量。
结果可视化我推荐直接画三个图：第一个是电功率平衡堆叠图，用area函数画出光伏、风电、CHP、储能、购电的叠加效果，一眼看过去就能检查是否满足电平衡；第二个是热功率平衡图，CHP和锅炉的出力占比一目了然；第三个是SOC曲线加购电电价柱状图叠放在同一个坐标系，直观显示储能是否在低价时段充电、高价时段放电。把这三张图放在论文里，审稿人第一印象就会觉得模型完整、结果可信。
6. 程序功能扩展与实用建议
6.1 从确定性模型扩展为不确定性模型
现在这版程序是确定性日前调度，光伏和风电预测值直接入模。如果你想往深做，优先级最高的扩展是把不确定性考虑进来。我建议不要一上来就上两阶段随机规划，先做鲁棒优化，因为鲁棒优化的代码改动量小，只需要把风光预测曲线换成区间预测，再在两处约束上加一个不确定预算参数，模型结构基本不变。做一个“盒式不确定集合”的鲁棒调度，代码工作量大概半天，效果却能被审稿人认为是实质性创新。
6.2 加入碳捕集与氢储能等新设备
因为模型采用模块化约束写法，新设备接入比较顺畅。以碳捕集装置为例，只需要在CHP和锅炉的排放出口加一个捕获率变量，在电母线上加一个碳捕集能耗项，在目标函数里加一个捕获成本和碳收益，约束集多两三行即可。氢储能类似，电解槽耗电、储氢罐容量、燃料电池发电三条链，每个环节对应一组约束。扩展时建议保持统一的命名风格，比如所有变量以P_开头、所有状态量以SOC_或SOH_开头，避免变量多了张冠李戴。
6.3 给新手的三个建议
第一，不要在还没跑通确定性模型时就去搞随机优化、强化学习或者元启发式算法，基础模型是最重要的基线，没有基线就不知道新算法的改进在哪。第二，求解结果出来后一定要做一次简单的物理校验——看看有没有出现“同时充电放电”“储能SOC越界”“热负荷不守恒”之类的低级错误，这比花时间调算法参数高效得多。第三，这个框架的价值核心是设备之间的耦合关系，你花时间懂原理，比会调参有用。模型不收敛时，先检查约束逻辑，再怀疑求解器。
我自己做这类项目时感受最深的是，综合能源系统优化调度本质上是“硬约束下的软决策”——设备特性决定可行域边界，目标函数决定边界内的最优方向。把边界摸清楚，把目标函数里的每一项成本来源、碳排来源核对清楚，程序跑出来的结果自然可信，后续写论文、做工程改造方案都有据可依。这套Matlab程序的价值不在于它多复杂，而在于它把这条“边界+方向”的思路完整落地了，你可以直接拿来当自己研究或项目的第一块基石。
