双碳目标下很多做能源优化调度的同行都在往“综合能源系统”这个方向靠。所谓综合能源系统核心就是把电、气、热、冷这几类能源放在同一个优化框架里统一调度让不同能源之间通过燃气轮机、电锅炉、吸收式制冷机这些耦合设备相互转换最终达到降低运行成本、提升新能源消纳、压减碳排放的目的。要用Matlab实现这套“低碳运行优化调度”说白了就是三步把物理系统抽象成数学模型把模型写成代码然后求解并读懂结果。这篇就围绕这个过程把模型、代码、算例和踩过的坑一次讲透。1. 场景设定与建模思路先搞清系统里有哪些设备、能做什么1.1 典型综合能源系统的物理架构我不喜欢上来就抛公式先看一个最常用的园区级冷热电联供系统结构外部电网、天然气网分别向园区供电供气园区内部有光伏PV、风电WT这类新能源还有燃气轮机GTCHP联产、燃气锅炉GB、电锅炉EB、蓄电池ES、蓄热罐HS、吸收式制冷机AC和电制冷机EC。负荷侧分成电负荷、热负荷、冷负荷三类。这套系统的关键特征就是“耦合”。燃气轮机烧天然气发电子同时产生余热这个余热可以优先供给热负荷也可以拖着吸收式制冷机制冷电锅炉则把电能转成热能相当于给“低谷电”找了一条出路蓄电池和蓄热罐负责平抑供需波动。能量在这里不是单向流动的所以调度不能像传统电力系统那样只盯功率平衡得把电、气、热、冷四条能量链条同步考虑。建模前必须想清楚一个问题你想做“规划”还是“运行调度”这两者模型形态差异很大。规划问题偏向容量配置选型、定容决策变量包含设备是否建设、容量大小时间尺度以年为单位运行调度问题则是在设备容量已定的前提下优化未来24小时或更短周期的各设备出力计划时间步长通常取1小时。这篇讨论的是后者也就是“给定设备怎么开、开多少”。1.2 双碳目标在调度模型里的落点双碳目标落到运行调度层面主要做两件事一是把碳排放相关成本嵌入目标函数让系统主动选择低碳运行方式二是在约束条件里加入碳排放上限碳配额迫使系统压排放。这里有一个很关键的认知转变传统调度只考虑“经济性最优”但在双碳语境下碳排放不再是一个“事后统计量”而是直接影响决策的约束和成本项。我在实际建模中的常见做法是引入碳交易机制系统获得一个免费碳排放配额实际排放超过配额时需要在碳市场购买配额计入成本低于配额时则可以出售多余配额获得收益。碳交易成本通常用一个阶梯价格的凸函数表示即超额越多、边际碳价越高。用生活化的方式理解就像家里用水用电有阶梯气价、阶梯电价用的越多单价越贵。调度模型里多设置几个碳排放阶梯档位系统在优化时会自动控制高碳设备的使用时段和出力水平避免“超档”导致成本暴涨这就是“低碳”这个属性在模型里的实际体现。2. 优化调度模型构建把调度规则写成数学语言2.1 决策变量分类与定义先明确决策变量这是写代码前必须完成的一步。我的习惯是先把所有变量列一张表再写代码否则后面拼矩阵时容易乱。类别变量维度说明连续变量P_grid(t)24x1从电网购电功率连续变量V_gas(t)24x1购入天然气体积流量连续变量P_GT(t), H_GT(t)24x1燃气轮机发电功率与热出力连续变量H_GB(t)24x1燃气锅炉热出力连续变量H_EB(t)24x1电锅炉热出力耗电折算连续变量P_ch(t), P_dis(t)24x1蓄电池充、放电功率连续变量H_ch_hs(t), H_dis_hs(t)24x1蓄热罐蓄、放热功率连续变量H_AC(t), C_EC(t)24x1吸收式制冷机耗热、电制冷机耗电连续变量C_AC(t), C_EC_c(t)24x1制冷出力0-1变量u_ch(t), u_dis(t)24x1蓄电池充放电状态标识0-1变量u_hs_ch(t), u_hs_dis(t)24x1蓄热罐蓄放热状态标识需要特别说明的是蓄电池和蓄热罐这两个储能同一时刻只能处在一个状态充或放如果不加0-1变量去约束优化器会计算出“同时充电又放电”这种物理上不可能的方案而这个方案恰恰会让储能变成无限容量的“能量搬运工”使目标函数虚低。所以储能建模一定要带二进制变量这点后面会细讲。2.2 目标函数经济性 碳成本目标函数采用最小化综合运行成本包括五个部分[ \min \sum_{t1}^{24} \left( C_{ele}(t) C_{gas}(t) C_{om}(t) C_{co2}(t) \right) ]其中购电成本C_ele(t) 分时电价(t) × P_grid(t)购气成本C_gas(t) 天然气单价 × V_gas(t)运维成本C_om(t) Σ 各设备的单位运维系数 × 对应出力。这几个都是线性项容易处理。碳成本部分最常见的两种写法写法一碳排放配额线性奖惩[ C_{co2}(t) \lambda_c \cdot \left( E_{actual}(t) - E_{quota}(t) \right) ]E_actual为实际碳排放折算值E_quota为免费配额λc为碳交易价格。排放低于配额时括号内为负相当于系统通过“卖碳”获得收益但这在代码里要小心处理因为目标直接做成了可正可负的线性项intlinprog完全支持不需要额外处理符号问题。写法二阶梯碳价更贴近真实碳市场把实际排放与配额差值分成好几个区间每个区间对应不同碳价比如超排0~300kg按50元/tCO2计300kg以上按80元/tCO2计。这种分段线性函数同样可以转化为线性约束加二进制变量但如果只是做初步研究直接用写法一就够了后面再扩展阶梯模型也来得及。实际碳排放E_actual的来源外购电力的间接碳排放电网排放因子×购电量加上燃气设备燃烧天然气的直接碳排放天然气排放因子×天然气消耗量再扣除可再生出力对应的零碳部分。注意光伏和风电在运行阶段基本算零碳所以模型天然会倾向于优先用新能源。2.3 约束条件平衡、容量、爬坡、储能约束是调度模型的灵魂。约束不全结果不物理约束太紧模型无解。我给你梳理一套最常用、也最不容易出错的约束清单。电功率平衡约束[ P_{PV}(t) P_{WT}(t) P_{GT}(t) P_{grid}(t) P_{dis}(t) P_{load}(t) P_{ch}(t) H_{EB}(t)/\eta_{EB} P_{EC}(t) ]左边是供应侧右边是需求侧。注意电锅炉的输入是电功率热出力要除以电热转换效率折回电功率电制冷机同理。很多第一次搭模型的人容易在这里漏掉输入侧功率折算然后发现电平衡怎么都对不上。热功率平衡约束[ H_{GT}(t) H_{GB}(t) H_{EB}(t) H_{dis_hs}(t) H_{load}(t) H_{ch_hs}(t) H_{AC}(t) ]冷功率平衡约束[ COP_{AC} \cdot H_{AC}(t) COP_{EC} \cdot P_{EC}(t) C_{load}(t) ]制冷这里经常写错吸收式制冷机消耗的是热功率H_AC乘以制冷系数COP后得到冷出力电制冷机消耗电功率P_EC同样乘以制冷COP。COP取值一般在0.7~1.2吸收式、3~5电制冷单位是“冷量/输入能量”别把热和电的单位混了。燃气轮机联产约束[ H_{GT}(t) \alpha_{HR} \cdot P_{GT}(t) ]背压式机组的热电比固定为α_HR这是最简单的联产模型工程上用得多。抽凝式机组的热电比在一个区间内可调模型会复杂不少初学阶段先用固定热电比把系统跑通再扩展不迟。设备出力上下限每个设备的出力需要落在额定范围内形式都是P_min ≤ P(t) ≤ P_max。这里有个实操要点如果设备有最小技术出力比如燃气轮机P_min不为0那么在低负荷时段优化器会面临“要么不出力、要么至少开到某个功率”的选择此时小修小补无法通过连续变量描述又得靠0-1变量来建模机组启停模型复杂度会显著上升。初始版本可以先把最小技术出力设成0即允许完全停机简化处理。爬坡约束对于燃气轮机和燃气锅炉这类响应速度受限的设备还需要加爬坡约束[ -P_{ramp} \le P_{GT}(t1) - P_{GT}(t) \le P_{ramp} ]这里P_ramp是设备在相邻时段内允许的最大出力变化量。时间步长是1小时的话爬坡率一般取额定功率的20%~50%每小时具体要看设备手册。储能SOC动态约束蓄电池的荷电状态SOC按下式递推[ SOC(t1) SOC(t) \frac{\eta_{ch} \cdot P_{ch}(t)}{E_{cap}} - \frac{P_{dis}(t)}{\eta_{dis} \cdot E_{cap}} ]同一时刻不能同时充放[ u_{ch}(t) u_{dis}(t) \le 1 ]充电功率受到上限约束时写法要带上二进制变量[ 0 \le P_{ch}(t) \le P_{ch_max} \cdot u_{ch}(t) ]很多新手不理解为什么这里必须乘u_ch直接写0 ≤ P_ch ≤ P_ch_max不也一样吗区别在于若P_ch和P_dis都没有乘各自的0-1变量仅靠“u_chu_dis≤1”约束优化器完全可以让P_ch100、P_dis50同时成立u_ch和u_dis虽然不会同时为1但这个物理上荒谬的出力组合并没有被禁止。正确写法必须是“功率上限由状态变量控制”也就是充电功率只有在u_ch1时才能取正值。蓄热罐的建模也完全一样。2.4 双碳指标怎么算、怎么写进约束如果只是想分析“低碳”效果最简单的方式是不把碳排放作为硬性约束只在目标函数里加碳成本。但如果希望模型具备“有碳配额约束”的实验场景就需要把碳排放上限作为一个不等式约束加进去[ \sum_{t1}^{24} E_{actual}(t) \le E_{max} ]E_max可以取无碳约束场景下系统总排放的一定百分比比如80%、70%看排放压减效果。这样跑出来的调度方案与基准场景对比碳排下降量和成本上升量就都有了做敏感性分析也方便。我通常会把E_max写成脚本开头的参数而不是硬编码在约束矩阵里便于批量扫描。如果把碳排放写成表达式放进约束矩阵要特别注意碳排放系数矩阵是作用在哪些变量上的。比如外购电的碳排放系数乘以P_grid(t)天然气碳排放系数乘以V_gas(t)这部分对所有时段求和后要小于E_max在intlinprog里就是一行A_ineq约束。别把光伏、风电、储能的变量错误地写进碳排放表达式里那是零碳设备系数为零。3. Matlab实现与求解从YALMIP到intlinprog3.1 求解器选型YALMIPCPLEX还是纯Matlab建模方式上有两条主流路线一是用YALMIP这类高级建模语言把目标函数和约束写得像数学公式一样再交给CPLEX、Gurobi或内置求解器二是直接用intlinprog/lsqlin等底层接口手工拼装矩阵。网上很多论文代码是YALMIP风格直观但需要额外下载YALMIP工具包和求解器实际工程中直接写intlinprog的代码也很多因为不依赖外部商业软件部署方便。我的建议是如果你想快速把模型调通YALMIP的开发和调试效率高得多尤其适合改约束、增删变量如果模型规模不大几十个变量、一百多个约束MATLAB自带的intlinprog完全够用还省去安装求解器的麻烦。MATLAB 2020a之后的intlinprog性能已经不错单算24时段调度基本秒出结果。如果想做大规模随机优化再考虑CPLEX/Gurobi不迟。3.2 用YALMIP建模的核心代码YALMIP描述综合能源调度非常顺手这里给出目标函数和关键约束的示例骨架%% 决策变量定义 P_grid sdpvar(24, 1); % 购电功率 V_gas sdpvar(24, 1); % 购气量 P_GT sdpvar(24, 1); % 燃气轮机发电 H_GB sdpvar(24, 1); % 燃气锅炉热出力 H_EB sdpvar(24, 1); % 电锅炉热出力 P_ch sdpvar(24, 1); % 蓄电池充电 P_dis sdpvar(24, 1); % 蓄电池放电 H_AC sdpvar(24, 1); % 吸收式制冷机耗热 P_EC sdpvar(24, 1); % 电制冷机耗电 u_ch binvar(24, 1); % 蓄电充状态 u_dis binvar(24, 1); % 蓄电放状态 %% 目标函数 Cost_ele price_ele * P_grid; % 购电成本 Cost_gas price_gas * sum(V_gas); % 购气成本price_gas为单位热值价 Cost_om k_GT * sum(P_GT) k_GB * sum(H_GB) k_EB * sum(H_EB) ... k_ES * sum(P_ch P_dis); % 运维成本 E_actual k_grid_co2 * sum(P_grid) k_gas_co2 * sum(V_gas); % 碳排放量 Cost_co2 lambda_co2 * (E_actual - E_quota); % 碳交易成本 Objective Cost_ele Cost_gas Cost_om Cost_co2; %% 约束条件 Constraints []; % 电平衡 Constraints [Constraints, P_grid P_PV P_WT P_GT P_dis ... P_load P_ch H_EB / eta_EB P_EC]; % 热平衡 Constraints [Constraints, alpha_HR * P_GT H_GB H_EB H_load H_AC]; % 冷平衡 Constraints [Constraints, COP_AC * H_AC COP_EC * P_EC C_load]; % 燃气轮机上下限与联产 Constraints [Constraints, 0 P_GT P_GT_max]; % 储能约束 Constraints [Constraints, SOC(2:24) SOC(1:23) eta_ch * P_ch(1:23) / E_cap ... - P_dis(1:23) / (eta_dis * E_cap)]; Constraints [Constraints, SOC(1) SOC_0, SOC(end) SOC_0]; % 周期调度 Constraints [Constraints, u_ch u_dis 1]; Constraints [Constraints, 0 P_ch P_ch_max .* u_ch]; Constraints [Constraints, 0 P_dis P_dis_max .* u_dis]; % 碳排放上限可选 Constraints [Constraints, E_actual E_max]; %% 求解 options sdpsettings(solver, cplex, verbose, 1); % 或 gurobi optimize(Constraints, Objective, options);YALMIP有个好处SOC递推约束可以直接用向量切片写不用手动循环。SOC(2:24) SOC(1:23) ...这种写法非常简洁效率也高。注意这里SOC变量最好声明成连续sdpvar但初始值SOC_0在第一个时段通过等式约束固定最后一个时段也要回到SOC_0周期调度否则蓄电池等于在调度周期结束时剩余电量不计成本结果会失真会把储能当作免费能量源。价格参数怎么设置分时电价可以设三档峰时1.2元/kWh10:00-15:00、18:00-22:00、平时0.75元/kWh7:00-10:00、15:00-18:00、谷时0.35元/kWh23:00-次日7:00。天然气的价格普通工商业在3.5~4.2元/m³之间如果是学术演示可以折算成单位热值价格比如0.85元/kWh天然气低位热值约9.7kWh/m³按3.8元/m³折算约为0.39元/kWh各文献取值差异大你用统一口径并在论文里注明即可。3.3 不用YALMIP直接用intlinprog手拼矩阵有些场景不希望引入外部工具包或者需要部署成独立exe这时候用intlinprog更稳。手拼矩阵的核心思路是把所有决策变量排成一个长向量x目标函数变成fx约束变成A_ineqx ≤ b_ineq和A_eqx b_eq。变量排列顺序建议和上面YALMIP顺序保持一致比如x [P_grid; V_gas; P_GT; H_GB; H_EB; P_ch; P_dis; H_AC; P_EC; u_ch; u_dis; SOC(2:24)];然后逐条约束“翻译”成矩阵行。这个过程最繁琐但也最能加深理解。我写过一个辅助工具用循环生成目标函数系数向量f比如购电成本对应f(1:24)price_ele购气成本对应f(25:48)price_gas碳成本对购电变量和购气变量的系数分别是lambda_co2k_grid_co2和lambda_co2k_gas_co2。这样拼装后调用intlinprogintcon find(ismember(1:numel(x), [变量索引序列])); % 指定整数变量 [x_opt, fval] intlinprog(f, intcon, A_ineq, b_ineq, A_eq, b_eq, lb, ub);还需要注意intlinprog默认变量下界为0如果想允许某个变量为负一般不需要设备出力都是非负的必须显式设置lb和ub向量。另外intlinprog求解成功后返回的是最优目标函数值fval需要把x_opt切片还原成各设备出力序列再画图、统计、算碳排放这一步建议单独写一个结果解析函数。4. 算例结果解读与敏感性分析别只盯着最优值4.1 典型日调度结果应该长什么样这里给一个我实际调过的案例某园区综合能源系统光伏装机800kW风机300kW燃气轮机额定500kW锅炉600kW电锅炉400kW蓄电池容量600kWh最大充放功率150kW蓄热罐容量600kWh最大蓄放热100kW制冷由一台吸收式250kW和一台电制冷400kW构成。运行24小时时间步长1小时。跑完之后有几件必做的事第一步是看电平衡堆叠图横轴1-24时段纵轴功率把光伏、风电、购电、燃气轮机发电、蓄电池放电从下往上堆叠再叠加电负荷曲线。正常情况下应该看到白天光伏出力大购电量自动下降夜间谷时电锅炉和蓄电池充电负荷上升刻意把负荷从峰时往谷时搬峰时段购电价格高燃气轮机出力抬升蓄电池放电。第二步看蓄电池SOC曲线正常是一条在0.2到0.9之间波动的连续曲线充放切换点最早出现在谷转平电价节点最晚出现在峰时段开始前。如果SOC曲线出现锯齿状高频抖动大概率是充放二进制约束没写对或者是爬坡约束缺失。第三步看碳排放分时曲线与碳交易成本。我发现加了碳约束后一个比较明显的变化是中午光伏大发、负荷也不高的时候系统会把燃气轮机停机多用外部电网的绿色电力前提是你设置了电网排放因子按时段变化或者电网电的碳排相对燃气轮机更低晚间负荷高峰则优先燃气轮机满发因为碳排放配额约束把它限定在特定时段使用这样总排放最小。下面给一个简化的结果表格数值我做了脱敏处理但趋势准确场景总运行成本元总碳排放kg新能源消纳率%无碳约束18620542086.5碳配额基准90%19350487094.2碳配额基准80%20790433097.8碳配额基准70%22800378098.5注意一个趋势碳排放压得越低成本越高但新能源消纳率在提升这是多目标问题的经典矛盾形态在论文里可以说成“低碳边际成本递增”。4.2 参数敏感性分析怎么做综合能源调度代码最大的价值就是可以快速做“what-if”分析。我常用的敏感性分析套路如下碳价扫描让lambda_co2从0到300元/tCO2线性增长跑20次记录每次的总成本、总排放、燃气轮机总发电量。你会得到一条排放随碳价上升而下降的曲线曲线中可能出现一个拐点这个拐点对应“燃气轮机从满发过渡到减少出力的临界碳价”。论文里放这张图非常加分因为它直观说明碳价工具的有效性区间。配额比例扫描E_max从无约束排放的95%逐级降到70%观察调度退化的过程。配额从95%降到85%时往往碳排放下降很明显、成本增幅不大再往下压成本开始陡增说明系统已经接近减排瓶颈。储容与碳排的联动逐步增大蓄电池和蓄热罐容量观察碳排放变化。储能容量越大系统越有能力把高碳时段的负荷搬到低碳时段比如谷电时段电锅炉蓄热碳排放越低但边际效果递减。这类分析很能说明“为什么双碳目标下要配置储能”。所有这些敏感性分析都只需要在代码外层套一个for循环把要扫描的参数定义成数组每次只改目标函数或约束中的一个参数。我在工程里还会顺手把结果存成结构体或表格避免跑完一组数据后还要手动抄写。5. 常见问题与排查技巧实录5.1 求解结果不收敛或很慢怎么办28个连续变量加48个二进制变量的模型对intlinprog来说是小菜一碟但如果把时间粒度从1小时改成15分钟变量数变4倍求解时间可能不增反降这是因为约束矩阵的稀疏性变了。遇到求解慢的情况优先检查三点二进制变量数量是否过多。储能设备如果加一吨“启停状态”每台设备24小时就是24个0-1变量三台设备就是72个还好。但如果给燃气轮机也加了启停0-1变量再考虑启停费用变量数量直接翻倍。初版模型不要把机组的启停变量全加上默认机组始终在线、只在出力上下限之间波动。intlinprog的时间上限和MIPGap设置。设置options optimoptions(intlinprog, MaxTime, 120, RelativeGapTolerance, 0.01)让求解器在相对误差1%以内就接受当前解。对工程应用来说1%的优化误差完全可以接受换来的求解时间大幅下降。约束是否存在退化。如果设备出力上限设得太小或储能初始SOC给得太紧模型可能无解intlinprog会提示“No feasible solution found”这种时候先放松约束比如把储能初始SOC的可变范围扩大或者去掉某个设备的上下限定位无解原因。5.2 储能SOC递推约束里最常见的错误我见过大量初学者在写SOC递推时把时间索引搞错。SOC(t1)和SOC(t)的约束是“相邻时段”的关系如果写成SOC(2:24) SOC(1:24) ...维度就对不上MATLAB直接报错。正确写法是左边用SOC(2:24)23个值右边SOC(1:23)23个值最后单独处理SOC(1)的初值。另外如果目标是让储能参与连续多日调度周期性的SOC(24)SOC(0)条件必须保留否则模型会在最后一个时段把电池彻底放空反正剩余电量不算成本这个结果不符合工程预期。还有一种隐蔽错误把SOC的范围约束设置成0到1但SOC在递推公式中可能超过这个范围导致无解。建议SOC上下限留5%的安全余量比如0.1到0.9这样约束矩阵更不容易锁死。5.3 冷热电平衡约束的维度错误很多人在把设备出力和负荷直接相减时忘记设备效率。比如电锅炉输入功率100kW效率90%输出热功率是90kW不是100kW。在冷平衡里电制冷机消耗电功率100kWCOP等于3.5制冷出力是350kW单位是冷吨或kW冷量和输入电功率是不同量纲。所有耦合设备的输入输出必须用效率或COP桥接这一条吃透了建出来的模型才具备物理意义。5.4 让模型更贴近工程现实的进阶方向基础调度模型跑通以后有两类常见扩展值得做。第一类是考虑新能源出力不确定性的随机优化和鲁棒优化把光伏、风电预测误差看作随机变量用一个典型场景集描述目标函数从单一路径优化变成期望值优化。MATLAB里做随机优化需要多场景建模对同一个系统跑多组负荷和新能源数据求解规模成倍增长此时要考虑分解算法。第二类是引入需求响应把一部分可平移负荷、可削减负荷作为灵活变量写进模型让负荷曲线跟随系统信号变化这会显著提高优化空间但同时增加建模难度。进阶的方向还有碳捕集与封存CCUS设备建模、多园区互联协同调度以及电、气、热网络的潮流模型耦合每一步扩展都会让模型更接近工程实际。5.5 模型验证与结果合理性检查每次跑完优化我习惯性做几项“合理性体检”查看所有设备的出力序列是否落在上下限范围内检查能量平衡约束的残差是否接近机器精度观察储能SOC初末值是否一致对比分时电价低谷时段是否出现了明显的储能充电和电锅炉蓄热动作最后再把总成本拆回到各个子项确认没有出现某一个成本分量异常波动。如果以上检查都通过基本可以确认模型没有原则性错误。数值计算上别忘在写代码前把单位统一功率用kW、能量用kWh、价格用元/kWh天然气的单位热值价格也要折算一致否则结果会出现灾难性的量级错误。我在实际做这个项目时最大的体会是调度优化这个方向数学模型本身并不难难的是“把模型写对、把约束写全、把结果讲清楚”。先跑通一个最简单的电热耦合模型再一步步加入冷负荷、储能、碳交易、需求响应每加一块都做一次旧场景回归测试如果结果异常就用控制变量法把新增模块单独隔离测试。这种“逐步搭乐高”的做法比一上来就堆一个大而全的模型要靠谱得多。最后再分享一个小技巧把碳价设成一个可扫描的参数让代码自动跑一组不同碳价下的调度结果并输出排放-成本曲线这张曲线几乎成了我所有项目报告里的标准配图无论写论文还是做汇报都特别能说明问题。