做综合能源系统优化调度最头疼的不是把目标函数写出来而是怎么在MATLAB里把所有物理约束、市场机制、储能时序都搭起来还能在合理时间内求出稳定结果。我前前后后用MATLAB做过好几版IES调度程序发现大家普遍卡在两个地方一是阶梯型碳交易机制怎么写成可线性化的成本函数二是氢能子系统电解槽、储氢罐、燃料电池怎么和电网络、热网络平滑耦合。这篇博文把我最近跑通的一套程序完整拆开讲程序考虑了阶梯型碳交易机制和氢能做24小时经济优化调度目标是最小化系统总运行成本。你只要照着模型改改负荷曲线和价格参数就能出结果特别适合写论文、做课程设计或者给园区能源调度做前期预研。1. 项目整体设计为什么这样组合1.1 典型综合能源系统的设备拓扑先交代一下我假设的系统结构这样后面所有公式都有落脚点。这个系统是一个典型的园区级综合能源系统和电网、天然气网双向连接。内部的主要设备有风电机组WT、光伏机组PV、燃气轮机GT、燃气锅炉GB、电解槽EL、储氢罐HS、氢燃料电池FC另外还加了一组电储能ESS用来削峰填谷。能量流动的路径是这样的风电和光伏发出来的电一部分直接供给电负荷一部分给电储能充电一部分送给电解槽制氢燃气轮机发电同时通过余热回收给热负荷供热燃气锅炉作为热负荷的补充热源燃气轮机发的电和燃料电池发的电都汇入电母线储氢罐承担氢的时序平移功能电解槽产出的氢先存入储氢罐再由储氢罐供给燃料电池和氢负荷。这套拓扑不算特别复杂但已经覆盖了电-热-氢三种能量形式的双向耦合。很多期刊论文里的IES模型都是这个骨架只是设备参数和算例数据不同。选这个拓扑的原因也很直接它既能体现多能互补的核心思想又不至于让模型规模失控等于是用最精简的设备集把综合能源系统的特点全表达出来了。1.2 阶梯型碳交易比固定碳价强在哪我在第一版程序里用的是固定碳价也就是每吨碳排放给一个固定价格比如50元/吨超配额就按这个价格买盈余就按这个价格卖。结果做出来的调度策略非常温和燃气轮机的出力不会明显受到碳价的约束因为固定碳价本质上是给碳排放加了一个常数项斜率优化器会在燃料成本和碳成本之间做一个线性的折中。后来我换了阶梯型碳交易机制效果立刻不一样了。阶梯碳价的核心逻辑是超出免费配额的部分分档设置碳价超得越多单价越贵。这和高收入人群个税累进税率的逻辑是一样的目的是让高排放行为付出递增的代价。把这个机制放进优化模型里系统就会更主动地控制碳排放。举个例子如果某时段燃气轮机出力很高、碳排放已经进入第二阶梯那么优化器会倾向于降低燃气轮机出力改用储氢罐里的氢通过燃料电池发电或者从电网低谷时段买电存到电储能里。这种策略性的转变是固定碳价模型里很难看到的。所以阶梯型碳交易不只是一个成本计算问题它会实质性地改变调度策略这也是这篇程序的一个核心创新点。1.3 氢能在这套模型里扮演什么角色氢能子系统是我认为整套模型最有意思的部分。电解槽本质上是把电能转化成氢能的一种时间平移装置光伏大发的中午时段如果电负荷用不完与其低价卖电给电网不如让电解槽制氢存起来等晚上用电高峰再让燃料电池发电。这样氢就成了一种跨时段的储能介质。和电储能相比氢能的优势主要体现在两个维度。第一个维度是能量密度储氢罐的单位能量容量比同等体积的电池高不少而且不会因为长期充放循环而明显衰减。第二个维度是灵活性氢既可以发电也可以直接满足氢负荷比如工业用氢、氢燃料汽车加注还可以通过燃料电池的余热回收参与供热。也就是说氢能是少数能同时接入电、热、氢三个网络的设备具有天然的多能耦合优势。我把氢能引入模型之后发现它对风光消纳的贡献非常明显。原来光伏大发时段只能靠电储能吸收现在电解槽可以直接把多余电能变成氢相当于给可再生能源多开了一条出口。这个行为在结果里能直观看到白天光伏峰值时段电解槽输入功率会明显增加储氢罐也开始积累能量。1.4 求解方案选型为什么用MILP模型搭建之前我先定了一个基调用混合整数线性规划MILP求解。原因有三点第一阶梯碳价本质上是一个分段线性函数通过引入0-1整数变量和大M法就能精确线性化不需要用到非线性的启发式算法。第二MILP可以保证全局最优解这一点在写论文对比算例时特别重要——审稿人只要看到你用的是启发式算法多半会问一句全局最优性怎么保证。第三MATLAB自带的intlinprog可以直接求解中小规模的MILP不需要额外的商业求解器也能跑通。当然如果你有YALMIP和Gurobi/CPLEX建模体验会更友好后面我会给出两种写法的要点。2. 碳交易与氢能的数学化建模2.1 碳排放核算与免费配额计算碳交易模型的第一步是核算系统一天的实际碳排放量。我这里的系统碳排放来源主要有两块一块是从电网购电对应的间接排放另一块是燃气轮机和燃气锅炉燃烧天然气产生的直接排放。计算式可以写成E_actual EF_grid * sum(P_buy) EF_gas * sum(P_gt H_gb)其中E_actual是全天实际碳排放量单位kgEF_grid是电网购电的排放因子我这里取0.85 kg/kWh这个值接近区域电网平均排放水平EF_gas是天然气燃烧的排放强度按输出能量折算取0.2 kg/kWh也就是每输出1 kWh的电或热功率对应0.2 kg的二氧化碳排放。P_buy是向电网购电的功率P_gt和H_gb分别是燃气轮机出力和燃气锅炉出力。这里有个细节容易混淆燃气轮机的排放应该按输入天然气的热值算还是按输出的电功率算两种写法都能建模关键是排放因子要和功率基准匹配。我偷了个懒按输出功率折算这样等式就统一成了功率平衡的形式省去了引入天然气热值转换的麻烦。免费配额E_quota我设成固定的1800 kg/天你也可以按负荷比例或者设备容量计算出动态配额。这个参数对结果影响很大我建议做敏感性分析时把它当做一个变量去扫描。配额给得越紧碳交易成本对调度的约束力越强氢能和储能的价值就越能体现出来。2.2 阶梯型碳成本的分段线性化碳交易的差额是E_diff E_actual - E_quota。如果E_diff 0说明实际排放超过了免费配额需要购买碳排放权如果E_diff 0说明配额有盈余可以把盈余卖出获利。我用的阶梯碳价结构如下表所示区间条件碳交易价格第一阶梯0 E_diff ≤ 500 kg50元/吨0.05元/kg第二阶梯500 E_diff ≤ 1000 kg60元/吨0.06元/kg第三阶梯E_diff 1000 kg70元/吨0.07元/kg盈余出售E_diff 040元/吨0.04元/kg注意盈余出售价格比购买价格低这个设计是有讲究的如果买卖价格相同优化器会想办法多买碳配额再卖出去套利虽然模型中配额是给定的但这种无风险套利路径在真实市场里也不存在所以我把出售价压低让模型更老实。现在说线性化的实现思路。阶梯碳成本是一个分段线性函数直接写进目标函数是凹函数优化器会找漏洞。标准做法是引入0-1整数变量z1、z2表示是否进入第二、第三阶梯再把E_diff拆成三段非负变量X1、X2、X3。数学表达是这样的E_diff X1 X2 X3 - Xs 0 ≤ X1 ≤ L1 0 ≤ X2 ≤ (L2 - L1) * z1 0 ≤ X3 ≤ M * z2 0 ≤ Xs ≤ M * (1 - z1)其中L1500L21000M是一个足够大的正数我取5000。Xs是盈余出售量只有当E_diff 0时才有值。还要加两个大M约束确保z1和z2与E_diff的数值是对应的E_diff ≥ L1 - M * (1 - z1) E_diff ≥ L2 - M * (1 - z2)为什么还要这两个约束因为如果没有它们优化器可能让z11但X20相当于白白多占了一个整数变量的自由度。加上之后只要实际排放没有超过对应阈值整数变量就强制为0模型逻辑就严密了。最终碳交易成本是C_co2 0.05 * X1 0.06 * X2 0.07 * X3 - 0.04 * Xs你可能已经注意到这种线性化方法本质上是把阶梯变成了三段不同价格的虚拟排放量之和。由于高价段的排放量只能在低价段填满之后才出现且低价段的单价更低所以优化器会自动先填满第一段、再填第二段、最后填第三段等效结果和真正的阶梯函数完全一致。2.3 电解槽-储氢罐-燃料电池协同模型氢能子系统的建模我全部采用功率平衡的能量等效方式不在模型中纠结于氢气质量单位这样和电功率、热功率的计量方式天然统一。电解槽的模型是一个线性转换关系P_h2_prod eta_el * P_el其中P_el是电解槽输入电功率eta_el取0.7P_h2_prod就是产氢功率以氢气热值计。这个效率水平对碱性电解槽来说是合理的质子交换膜电解槽可以到0.75以上但0.7稳妥一点。燃料电池的模型正好反过来P_fc eta_fc * P_fc_in其中P_fc是燃料电池输出的电功率P_fc_in是消耗的氢功率eta_fc取0.45。燃料电池发电时还有余热回收我设热回收系数为0.35所以燃料电池同时提供的热功率是H_fc 0.35 * P_fc_in。储氢罐的动态约束采用能量状态SOC形式H_s(t1) H_s(t) eta_hc * P_h2_prod(t) - P_h2_con(t) / eta_hdeta_hc和eta_hd分别是充氢效率和放氢效率都取0.95表示储氢过程中有一定损失。P_h2_con是储氢罐向外输出的氢功率主要用于燃料电池和直接供应氢负荷P_h2_con(t) P_fc_in(t) H_load(t)储氢罐容量我设成1000 kWh按氢气能量计SOC上下限是10%和90%。此外还要在约束里加上储氢罐的充放功率上限防止模型在一小时内把整个储氢罐灌满或放空这个约束很多人会漏掉漏掉之后调度结果里会出现很多不合理的脉冲式充放。3. 运行优化模型目标函数与约束体系3.1 目标函数构建整个优化调度的目标函数是系统日运行总成本最小化包括六项购电成本、售电收入、天然气成本、设备运维成本、碳交易成本和储能退化成本。数学形式如下min C_total sum(Price_buy .* P_buy) % 购电成本 - sum(Price_sell .* P_sell) % 售电收入负向计入 C_gas % 天然气成本 C_om % 运维成本 C_co2 % 碳交易成本 C_battery_degradation % 电储能退化成本购电采用分时电价我设计的电价结构是谷段0-6时0.3元/kWh平段7-17时0.5元/kWh峰段18-23时0.9元/kWh。售电上网电价固定为0.35元/kWh。天然气成本的计算C_gas gas_unit_price * (sum(P_gt / eta_gt) sum(H_gb / eta_gb))这里eta_gt是燃气轮机发电效率取0.3eta_gb是燃气锅炉热效率取0.9gas_unit_price是天然气折算到每kWh热量的价格取0.25元/kWh。这个折算方法用的是按能量输入计费的思路比直接按立方米计费更便于和功率模型集成。运维成本按设备的单位出力维护费率累加。我用的费率是燃气轮机0.03元/kWh燃气锅炉0.02元/kWh电解槽0.01元/kWh按输入电功率算燃料电池0.02元/kWh按输出电功率算电储能0.02元/kWh按充放电功率之和算风电和光伏的运维成本很低各取0.01元/kWh。电储能退化成本是我后加的一项按充放功率的0.02元/kWh计。为什么要加因为如果不加优化器会让电储能高频充放一天下来充放电循环次数可能达到十几个实际电池早废了。加上这个小惩罚结果会平稳很多。3.2 系统平衡约束与设备约束电功率平衡约束。这个约束是整套模型的核心它把所有发用电设备串在一条母线上P_buy(t) P_sell_pv(t) P_wt(t) P_gt(t) P_fc(t) P_dis(t) P_load(t) P_el(t) P_ch(t) P_sell(t)这里P_buy是购电功率P_wt和P_pv是可调度的风光出力0到预测上限之间P_dis和P_ch分别是电储能放电和充电功率。注意P_sell和P_buy在同一个小时不能同时为正所以我加了互斥约束或者用符号约束限制。热功率平衡约束H_gt(t) H_gb(t) H_fc(t) H_load(t)其中H_gt是燃气轮机余热回收功率按热电比1.2折算H_gt 1.2 * P_gt。这个热电比意味着燃气轮机每发1 kWh电同时产出1.2 kWh热符合典型的热电联产机组参数范围。氢功率平衡约束在2.3节已经给出这里再强调一下它是把电解槽产氢、储氢罐充放、燃料电池耗氢、直接氢负荷四者联系在一起的纽带。设备出力上下限约束P_gt_min ≤ P_gt(t) ≤ P_gt_max 0 ≤ H_gb(t) ≤ H_gb_max 0 ≤ P_el(t) ≤ P_el_max 0 ≤ P_fc(t) ≤ P_fc_max同时还要给燃气轮机和燃气锅炉爬坡约束一小时内的出力变化不能超过额定出力的30%。这个约束在纯线性的静态优化中很容易被忽略但真实机组必须满足。电储能约束E_s(t1) E_s(t) eta_ch * P_ch(t) - P_dis(t) / eta_dis E_s_min ≤ E_s(t) ≤ E_s_max E_s(1) E_s(25) % 周期平衡避免首日SOC流浪周期平衡约束是必须的。如果不加优化器会在第一天把储能耗尽第二天再从电网买电充满导致整个24小时结果失真。加上这个约束后储能效应就局限在当天内部的时间平移更符合日前调度的设定。设备参数汇总如下表设备参数数值燃气轮机出力范围60~200 kW燃气轮机电效率0.3燃气轮机热电比1.2燃气锅炉出力范围0~500 kW燃气锅炉热效率0.9电解槽输入范围0~200 kW电解槽制氢效率0.7燃料电池输出范围0~150 kW燃料电池电效率0.45储氢罐容量1000 kWh储氢罐SOC范围0.1~0.9电储能容量300 kWh电储能最大充放功率60 kW4. MATLAB程序实现从数据到结果4.1 24小时负荷与新能源出力输入优化调度需要一整天的时序数据。我这里的算例数据如下表表中风电和光伏给出的是可用出力上限实际运行中系统可以选择是否部分弃用时段电负荷热负荷氢负荷风电上限光伏上限1802003080027519025950370180201100470170151300568165151500680160201200712015030140108180140501004092301207080901026011070601501127010060501901226595508021013250905070200142409560601801523510070901601623811060120130172601305011010018300150409080193201804060502030020050401521260210603002222021570500231602106080024110205501000注意时段1对应凌晨0点。这套数据的特征是电负荷有两个峰中午时段受生产活动影响有一个小峰晚上18-20点达到全天最高热负荷在深夜和傍晚较高白天略降氢负荷相对平稳白天略高。风电和光伏都具有明显的时序特性特别是光伏在中午时段出力达到210 kW而同时电负荷只有250-270 kW多余电量正好可以被电解槽吸收。在MATLAB里我把这些数据统一存在data_IES.m文件中定义为结构体数组方便主程序调用。4.2 核心求解代码阶梯碳成本与氢能约束下面这段代码展示的是YALMIP框架下的核心建模逻辑。我用YALMIP是因为它的约束表达更直观适合快速验证模型逻辑。如果你没有安装YALMIP也可以把同样的约束整理成矩阵形式用MATLAB自带的intlinprog求解。%% 优化变量定义 P_buy sdpvar(1, 24, full); % 购电功率 kW P_sell sdpvar(1, 24, full); % 售电功率 kW P_gt sdpvar(1, 24, full); % 燃气轮机发电功率 kW H_gb sdpvar(1, 24, full); % 燃气锅炉热功率 kW P_el sdpvar(1, 24, full); % 电解槽输入电功率 kW P_fc sdpvar(1, 24, full); % 燃料电池输出电功率 kW P_ch sdpvar(1, 24, full); % 电储能充电功率 kW P_dis sdpvar(1, 24, full); % 电储能放电功率 kW E_s sdpvar(1, 25, full); % 电储能 SOC 状态 kWh H_s sdpvar(1, 25, full); % 储氢罐 SOC 状态 kWh %% 阶梯碳交易区间变量 X1 sdpvar(1, 1); % 第一阶梯排放量 kg X2 sdpvar(1, 1); % 第二阶梯排放量 kg X3 sdpvar(1, 1); % 第三阶梯排放量 kg Xs sdpvar(1, 1); % 盈余出售量 kg z1 binvar(1, 1); % 是否进入第二阶梯 z2 binvar(1, 1); % 是否进入第三阶梯 %% 目标函数 Cost_buy sum(Price_buy .* P_buy); Income_sell sum(Price_sell .* P_sell); Cost_gas gas_price * (sum(P_gt ./ eta_gt) sum(H_gb ./ eta_gb)); Cost_om sum(om_gt .* P_gt) sum(om_gb .* H_gb) ... sum(om_el .* P_el) sum(om_fc .* P_fc) ... sum(om_ess .* (P_ch P_dis)); Cost_co2 0.05 * X1 0.06 * X2 0.07 * X3 - 0.04 * Xs; Objective Cost_buy - Income_sell Cost_gas Cost_om Cost_co2; %% 碳排放与阶梯碳约束 E_actual EF_grid * sum(P_buy) EF_gas * sum(P_gt H_gb); E_diff E_actual - E_quota; E_diff X1 X2 X3 - Xs; M 5000; L1 500; L2 1000; Constraints []; Constraints [Constraints, 0 X1 L1]; Constraints [Constraints, 0 X2 (L2 - L1) * z1]; Constraints [Constraints, 0 X3 M * z2]; Constraints [Constraints, 0 Xs M * (1 - z1)]; Constraints [Constraints, E_diff L1 - M * (1 - z1)]; Constraints [Constraints, E_diff L2 - M * (1 - z2)]; Constraints [Constraints, z2 z1]; %% 电功率平衡 Constraints [Constraints, ... P_buy P_wt P_pv P_gt P_fc P_dis ... P_load P_el P_ch P_sell]; %% 热功率平衡 H_gt 1.2 * P_gt; H_fc 0.35 * (P_fc / eta_fc); Constraints [Constraints, H_gt H_gb H_fc H_load]; %% 氢平衡与储氢罐动态 P_h2_prod 0.7 * P_el; % 电解槽产氢功率 P_fc_in P_fc / 0.45; % 燃料电池耗氢功率 P_h2_con P_fc_in H_load; % 储氢罐输出氢功率 for t 1:24 Constraints [Constraints, H_s(t1) H_s(t) 0.95 * P_h2_prod(t) - P_h2_con(t) / 0.95]; end Constraints [Constraints, H_s_min H_s(1:24) H_s_max, H_s(1) H_s(25)]; %% 电储能相关约束 for t 1:24 Constraints [Constraints, E_s(t1) E_s(t) 0.95 * P_ch(t) - P_dis(t) / 0.95]; Constraints [Constraints, P_ch(t) 60, P_dis(t) 60]; Constraints [Constraints, E_s_min E_s(t) E_s_max]; end Constraints [Constraints, E_s(1) E_s(25)]; %% 求解 ops sdpsettings(solver, gurobi, verbose, 0); result optimize(Constraints, Objective, ops); %% 结果提取 P_buy_opt value(P_buy); P_fc_opt value(P_fc); E_s_opt value(E_s);这段代码基本是完整可运行的骨架。有几个需要解释的细节第一P_wt和P_pv是可调度变量但我在代码里直接用预测上限赋值如果你想开放弃风弃光就把P_wt和P_pv改成sdpvar变量加一个0 P_wt P_wt_forecast的约束就行。第二H_s(1) H_s(25)和E_s(1) E_s(25)是周期平衡约束保证储能SOC首尾一致。如果你不用YALMIP用intlinprog的话关键就是把上面的等式和不等式整理成标准矩阵形式。整数变量只有z1和z2其余全是连续变量所以规模并不大大约有150个连续变量和50个约束intlinprog几秒钟就能解完。4.3 结果输出与合理性校核程序跑完之后我习惯先画三张图第一张是电功率平衡图展示各设备的出力时序第二张是热功率平衡图第三张是储氢罐和电储能的SOC曲线。这三张图能快速判断结果是否物理上合理。有一次我跑出来的结果里储氢罐SOC曲线出现了锯齿状跳变白天充满、傍晚放空、晚上又充满。后来一查是因为我忘了加充放功率上限约束模型在边界点疯狂充放。加上P_h2_in_max和P_h2_out_max之后曲线就平滑了。还有一个很实用的校核手段把碳交易成本拆出来单独看。如果系统实际排放低于配额碳成本是负数盈余收益说明调度策略倾向于清洁供能如果碳成本占了总成本的很大比例说明系统对天然气和购电的依赖太高可以考虑增加氢能设备的容量。5. 调试、排错与扩展心得5.1 常见报错与对应排查方法我在调试这套程序的过程中踩了不少坑挑几个有代表性的分享出来。第一个坑是求解器报infeasible problem。出现不可行问题首先检查电功率平衡约束里有没有漏掉某个负荷项。我踩过一次是因为氢负荷放在了热平衡约束里导致氢平衡约束无解。排查方法是用check(Constraints)逐个约束看残差YALMIP会在约束不满足时标出哪个约束的残差最大照着找就行。第二个坑是阶梯碳价的整数变量没有绑定好导致碳成本变成负的巨额收益。这是因为X1、X2、X3三个变量之间缺少递进约束优化器会想办法让高价段变量取负值来套利。解决办法就是把X1 0、X2 0、X3 0的非负约束全部写死并加上z2 z1这个逻辑约束。第三个坑是储氢罐SOC的初值和末值不一致。如果把H_s(1) H_s(25)去掉结果就是储氢罐在调度周期内累积了大量氢能相当于白赚了一个初始状态。加上周期平衡约束后结果才是真正的经济调度。第四个坑是电价和碳排放因子的单位混淆。我在第一版里把电价单位写成了元/MWh结果成本爆炸。后来统一用元/kWh所有功率单位统一用kW时间尺度统一为1小时才把单位理顺。下面整理一份常见问题速查表问题现象可能原因排查方法求解器报infeasible电/热/氢平衡约束漏项或某个设备出力范围过小用check(Constraints)查看残差最大的约束碳交易成本异常为负阶梯变量X1 X2 X3缺少非负约束或递进逻辑检查z2z1和分段边界约束结果中储能SOC首尾不一致缺少周期平衡约束增加SOC(1)SOC(25)购电和售电同时为正缺少互斥约束在P_buy和P_sell之间加P_buy.*P_sell0或整数调整碳价对调度结果无影响免费配额设得过大实际排放远低于配额逐步减小E_quota做敏感性分析风光大量弃用但燃料电池满发氢能子系统容量太小增大电解槽/储氢罐容量重新求解5.2 模型扩展方向这套模型虽然已经覆盖了电-热-氢三种能量流但扩展空间还是很大的。我列几个我打算后续做的方向。第一个方向是引入需求响应。目前电负荷是刚性给定的如果把一部分可平移负荷变成决策变量比如电动汽车充电负荷可以让负荷曲线主动匹配新能源出力曲线进一步提升系统经济性。第二个方向是加入碳捕集设备。碳捕集和电转气在功能上是互补的碳捕集减少排放电转气把多余电能转化为天然气或氢能。两者协同起来系统对电网和天然气的依赖会进一步降低。第三个方向是改成多园区协同调度。单个园区的规模有限风光互补性不强。如果把多个园区放在一起让园区之间通过公共母线交换电力局部的弃风弃光很可能变成其他园区的可用电源整体效益更好。第四个方向是在求解层面加入滚动时域优化MPC。日前调度只能依赖预测数据当天执行时如果风光出力偏差大就需要滚动修正。修改思路其实不复杂就是把时间窗从24小时缩短到4-6小时每隔1小时滚动求解一次。这套程序我实际跑下来最深的体会是模型的价值不仅在于算得准更在于能让你清晰看到不同机制对调度结果的影响。固定碳价版本和阶梯碳价版本对比系统对燃气轮机的使用策略差异非常明显不装电解槽和装电解槽对比中午光伏时段的弃光率变化也一目了然。这种可解释性正是做优化调度最重要的东西也是这套带氢能和阶梯碳交易的MATLAB程序最值得借鉴的地方。你在自己的项目里完全可以拿这份模型当底子把设备、价格、碳配额参数换成自己的数据稍作调整就是一篇扎实的算例分析。