看到这个题目我太熟了这不就是这两年IES领域最常被翻牌子的优化调度方向之一嘛。说是复现实际上论文里一堆参数和约束写得含含糊糊真要跑通、跑出和原文趋势一致的结果没点踩坑经验还真容易卡死。这篇我就把这套策略从建模到MATLAB求解的完整思路捋一遍把我实际跑代码时趟过的雷也一并交代清楚希望能让正在读文献的你少走点弯路。1. 复现前必须想明白的三个问题1.1 这个课题到底在优化什么先说最根本的问题综合能源系统优化调度本质上是在满足电、热、气各类负荷需求的前提下让整个系统的运行成本最低、碳排放最少。听起来简单但一旦把综合需求响应和阶梯型碳机制加进来事情就变得复杂起来。综合需求响应IDR不是说简单地让用户少用点电而是通过价格或者激励手段引导用户调整电、热、气多种能源的用能习惯和用能时段。比如电价高的时候电负荷可以转移一部分到低谷时段或者用燃气替代一部分用电需求这就是综合两个字的核心——多种能源之间的互补与替代。阶梯型碳机制则是碳交易市场的产物它把碳排放权分成几个区间每个区间的碳价不一样排得越多碳价越高。这个机制的设计初衷就是让高排放企业付出更高的经济代价从而倒逼低碳转型。所以整个优化调度策略要回答的核心问题是在24小时的时间尺度上如何安排各设备的出力、如何制定需求响应方案、如何购买碳配额才能让总成本最小同时保证系统安全稳定运行。1.2 核心难点到底卡在哪我刚开始读这类文献的时候觉得建模思路挺清晰的但真正动手后才发现难点不在理解而在建模粒度和求解复杂度的平衡上。首先是设备模型。电转气P2G、燃气轮机、电锅炉、储电储热、余热回收……每个设备都有各自的效率曲线、爬坡约束、启停约束。把这些设备全部用线性化方程表示出来约束条件会膨胀到几百条变量也会多出好几个量级。其次是需求响应建模。你要是只做削峰填谷式的简单建模那就是用可转移负荷把高峰负荷挪到低谷用0-1变量表示转移状态这还算好处理。但综合需求响应还涉及到电转热、电转气这样的能源替代行为比如电价过高时用户可以选择用燃气锅炉供暖而不是电采暖。这种行为建模起来就是一个非线性耦合关系处理不好很容易把模型变成非线性规划求解难度直线上升。最后是阶梯碳机制的引入。它不是简单地加一个碳成本项到目标函数里而是会引入分段函数而分段函数在优化里通常意味着0-1变量和线性化逼近这会让模型多出一大批整数变量。整个模型最后大概率是一个混合整数线性规划MILP问题规模不大但很典型。1.3 为什么选MATLAB而不选Python这一步估计是很多人的疑惑。其实Python的PuLP、Pyomo也能做优化建模但坦白讲在能源系统调度领域MATLAB YALMIP CPLEX/Gurobi这个组合仍然是目前文献复现的主流配置没有之一。原因其实很现实YALMIP这个工具箱对优化问题建模的支持太成熟了它能把优化问题的目标函数和约束条件以非常接近数学公式的方式写出来这也是复现文献时最看重的一点。因为你要对照论文里的公式一行行验证自己的代码有没有写错YALMIP的语法天然适合这种工作流。还有一个很实际的理由你搜到的参考文献十有八九都附带了MATLAB代码或者MATLAB风格的伪代码用MATLAB去复现代码是一致的语言而Python还得做一层翻译转换出错了很难排查。2. 数学模型详解从目标函数到约束条件2.1 目标函数成本项拆解在做复现时我建议先把目标函数拆成四块来看购能成本、设备运行维护成本、碳交易成本和需求响应补偿成本。目标函数可以写成% 目标函数代码骨架详细代码见后续章节 Objective Cost_energy Cost_om Cost_carbon Cost_dr;购能成本指的是系统从外部购买电力和天然气的费用。这里有个细节要注意购电价格往往分为峰谷平三个时段对应不同的电价购气价格则相对稳定。我复现时遇到过一个问题就是没有把分时电价对需求响应的影响真正耦合起来——这个后面细说。设备运行维护成本比较好理解就是燃气轮机、锅炉、P2G等设备在运行过程中产生的维护费用通常与设备出力呈线性关系。碳交易成本是重头戏。阶梯型碳机制下的成本函数长这样碳交易成本 碳价 * 实际碳排放量 - 免费碳配额但实际上阶梯碳机制下碳价是分段的。假设你排了2500kg碳免费配额有2000kg那你需要购买的配额是500kg。这500kg又被分成不同的区间前300kg是一个碳价剩下的200kg是一个更高的碳价。这个函数在目标函数里就是分段函数处理起来需要引入辅助变量和0-1变量。需求响应补偿成本则是付给用户的。用户响应了削峰指令把负荷从高峰时段转移走系统需要给用户一定的补偿这个补偿和负荷转移量以及转移时长有关。我在文献里看到过两种处理方式一种是直接按转移量线性补偿另一种是阶梯式补偿。我在复现时用的线性补偿效果已经不错了阶梯式补偿会带来额外的整数变量对求解速度影响很大。2.2 约束条件比想象中多得多约束条件看起来多但仔细梳理下来就四类功率平衡约束、设备出力约束以及储能约束和需求响应约束。功率平衡约束是系统建模的核心。电、热、气三种能量都要满足各自的平衡关系。以电功率平衡为例需要满足购电功率 燃气轮机发电 风电出力 储电放电 电负荷 电锅炉耗电 P2G耗电 储电充电这里最容易出错的地方是P2G的建模。P2G先把电转化成氢气氢气又可以和二氧化碳反应生成天然气但这个过程在简化的调度模型里通常被处理成一个电转气的效率环节也就是耗电量转成天然气产气量产气量有上下限约束。这个简化处理在多数文献里是通用的但你要是想把动态特性也加进去模型复杂度就完全不同了。设备出力约束就比较好理解。燃气轮机的出力功率有上下限爬坡速率也是关键约束——上一时刻和这一时刻的出力差不能超过爬坡限值。这个约束对调度结果的平滑性影响很大我调试代码时发现如果不加爬坡约束燃气轮机的出力曲线会非常跳跃完全不符合实际。储能约束相对复杂因为储能设备有自身的状态连续性。储电设备的SOC荷电状态是一个动态递推公式今天的末尾SOC 昨天的末尾SOC 充电效率充电功率 - 放电效率放电功率。这个递推关系在约束里体现为一条跨时段的等式链。很多复现代码写不出来或者写错了问题就出在这个递推没有处理好。2.3 阶梯碳机制处理的关键技巧阶梯碳机制的分段成本函数处理起来是我觉得整个模型中最考验技巧的部分。这里我分享一种常用且高效的线性化处理方法。假设碳排放量被分为三个阶梯每个阶梯的碳价递增。那么需要引入三个连续变量 ( c_1, c_2, c_3 ) 来代表每个阶梯的实际碳排放量以及三个0-1变量 ( z_1, z_2, z_3 ) 来指示是否达到了该阶梯。约束条件大致如下% 阶梯碳建模示意 C_total C_1 C_2 C_3; % 总碳排放量 C_total C_actual_total; % 实际总碳排放 % 第一阶梯约束 C_1 Cap_1 * z_1; C_1 0; % 第二阶梯约束 C_2 Cap_2 * z_2; C_2 Cap_1 * z_2; % 只有当第一阶梯满了第二阶梯才开始这里还有一个细节分段碳价和免费配额的关系。有些论文采用碳排放量减去免费配额后超出部分才参与阶梯计价的处理方式也就是阶梯对应的是购买量而不是排放量。这两种处理在语义上差很多复现的时候一定要看清论文用的哪种不然碳成本的计算结果会完全对不上。3. MATLAB复现实操环境配置与代码结构3.1 环境准备工具箱和求解器怎么选复现这个课题先把环境准备好。我自己用的是MATLAB R2021a版本其实R2019以上都行需要额外装两个东西YALMIP工具箱和求解器。YALMIP是一个MATLAB优化建模工具箱它最大的价值在于让你用接近数学语言的方式表达优化问题然后由它底层调用具体的求解器。安装方式很简单去GitHub下载YALMIP的源码包把yalmip文件夹加到MATLAB路径即可。求解器方面我强烈建议用Gurobi或者CPLEX它们是商业求解器对MILP问题的求解效率远高于MATLAB内置的intlinprog。不过这两个求解器都需要申请学术许可。如果你没有学术许可也可以先用MATLAB自带的intlinprog将就一下模型规模小的话也能跑但求解速度会慢不少。配置完成后在MATLAB里测试一下% 测试YALMIP是否安装成功 x sdpvar(1,1); optimize([x 0, x 1], -x); value(x) % 应该输出 1如果能正常输出1说明YALMIP和求解器都已经正确配置了。3.2 数据准备负荷和风电的典型日曲线复现过程中数据的选取直接决定了结果图长什么样。我发现很多新手不知道这部分才是最耗时的环节。你需要准备的数据包括24小时的负荷数据电负荷、热负荷、气负荷、24小时的风电出力数据、分时电价数据、天然气价格、各设备的效率参数、储能设备参数等。如果论文没有给出具体数据通常的做法是采用一个典型日的场景冬季或者夏季的典型日负荷曲线。负荷曲线通常保留两位小数风电数据则是一个0到1之间的归一化数值乘以装机容量。我实际用的数据大概长这样% 24小时电负荷数据单位kW P_load [280 260 250 240 230 240 260 320 380 420 450 460 440 430 420 410 400 420 450 480 470 420 380 320]; % 分时电价单位元/kWh price_electricity [0.4 0.4 0.4 0.3 0.3 0.3 0.5 0.7 0.8 0.9 0.9 0.8 ... 0.7 0.6 0.5 0.5 0.6 0.7 0.8 0.9 0.8 0.6 0.5 0.4]; % 风电归一化出力曲线 P_wind_norm [0.3 0.35 0.4 0.45 0.5 0.45 0.3 0.2 0.15 0.1 0.12 0.15 ... 0.2 0.18 0.15 0.12 0.1 0.15 0.2 0.25 0.2 0.15 0.12 0.1];数据这块要特别注意不同设备的单位可能不一样。比如燃气轮机的发电效率是40%天然气的热值是9.7kWh/m³如果单位换算错了整个模型的结果都会是错的。我在复现时专门写了一个参数初始化脚本把所有单位统一成kW和kWh这样后期调试会省心很多。3.3 核心代码框架YALMIP搭建优化模型接下来是核心部分我用YALMIP搭建整个优化模型。以24小时调度为周期时间步长为1小时每个决策变量都是24维的向量。决策变量包括购电功率、购气量、燃气轮机出力、P2G出力、电锅炉出力、储电充放电功率、储热充放热功率以及需求响应相关的负荷转移量。这样定义变量%% 定义决策变量 P_buy sdpvar(1, 24); % 购电功率 Gas_buy sdpvar(1, 24); % 购气量 P_gt sdpvar(1, 24); % 燃气轮机发电功率 P_p2g sdpvar(1, 24); % P2G耗电功率 H_eb sdpvar(1, 24); % 电锅炉产热功率 % 储能设备 P_ch sdpvar(1, 24); % 储电充电功率 P_dis sdpvar(1, 24); % 储电放电功率 H_ch sdpvar(1, 24); % 储热充电功率 H_dis sdpvar(1, 24); % 储热放电功率 % 需求响应相关 P_shift sdpvar(1, 24); % 电负荷转移量正为转入 P_cut sdpvar(1, 24); % 电负荷削减量 H_shift sdpvar(1, 24); % 热负荷转移量搭建约束条件时YALMIP的逻辑非常直观基本上就是照着数学公式来写。比如电功率平衡约束%% 电功率平衡约束 C_ele_balance []; for t 1:24 C_ele_balance [C_ele_balance, P_buy(t) P_gt(t) P_wind(t) P_dis(t) P_shift(t) P_load(t) P_p2g(t) P_eb(t) P_ch(t) P_cut(t)]; end这里要注意我们的储能设备不能同时充放电所以需要加一个互斥约束%% 储能互斥约束 z_ch binvar(1, 24); % 充电状态指示 z_dis binvar(1, 24); % 放电状态指示 C_storage [z_ch z_dis 1]; C_storage [C_storage, P_ch P_ch_max * z_ch]; C_storage [C_storage, P_dis P_dis_max * z_dis];这个互斥约束是复现代码时特别容易忽略的。很多新手直接写 ( P_ch \le P_ch_max ) 和 ( P_dis \le P_dis_max )但忘了加上 ( z_ch z_dis \le 1 ) 这个约束结果求解器会让储能一边充电一边放电白白损耗能量那结果看起来肯定很荒唐。3.4 碳成本模块的实现碳成本模块是全模型中最容易出问题的地方。这里我给出阶梯碳机制的一个完整 YALMIP 实现例子%% 阶梯碳成本建模 % 碳排放总量计算 Carbon_total sum(P_gt * e_gt P_buy * e_grid H_gb * e_gas P_p2g * e_p2g); % 免费配额 Carbon_quota quota_total; % 实际需要购买的排放量 Carbon_purchase Carbon_total - Carbon_quota; % 阶梯参数 carbon_levels [0, 1000, 2000]; % 阶梯边界 carbon_prices [50, 80, 120]; % 阶梯碳价元/t % 引入分段变量 Carbon_seg sdpvar(1, 3); z_seg binvar(1, 3); % 分段约束 C_carbon [Carbon_purchase sum(Carbon_seg)]; for k 1:3 if k 1 C_carbon [C_carbon, Carbon_seg(k) 0]; C_carbon [C_carbon, Carbon_seg(k) carbon_levels(2) * z_seg(k)]; C_carbon [C_carbon, Carbon_seg(k) Carbon_purchase]; elseif k 2 C_carbon [C_carbon, Carbon_seg(k) 0]; C_carbon [C_carbon, Carbon_seg(k) (carbon_levels(3) - carbon_levels(2)) * z_seg(k)]; else C_carbon [C_carbon, Carbon_seg(k) 0]; C_carbon [C_carbon, Carbon_seg(k) 1e6 * z_seg(k)]; end end % 碳成本进入目标函数 Cost_carbon carbon_prices(1) * Carbon_seg(1) ... carbon_prices(2) * Carbon_seg(2) ... carbon_prices(3) * Carbon_seg(3);有一点要提醒阶梯机制的触发顺序依赖 ( z_1 \ge z_2 \ge z_3 ) 这个逻辑但因为在目标函数中碳价递增求解器会自动让第一阶梯优先填满所以就算不加这个约束结果一般也不会出错。不过为了稳妥我建议还是加上这个约束。3.5 综合需求响应的建模实现综合需求响应的关键是综合两个字。我采用的方法是把负荷分成三类可转移负荷、可削减负荷、可替代负荷。可转移负荷就是那种可以在不同时段之间移动的负荷比如洗衣机的用电时段可以从峰时挪到谷时。建模时用转移量来表示%% 可转移负荷约束 % 转移前后总用电量守恒 C_transfer [sum(P_shift) 0]; % 总转移量为0 % 每个时段转移量限制 for t 1:24 C_transfer [C_transfer, -P_shift_max(t) P_shift(t) P_shift_max(t)]; C_transfer [C_transfer, P_load_shifted(t) 0]; end可削减负荷是指那种用户愿意牺牲一部分用能需求来换取补偿的负荷比如空调温度调高一度减少的这部分就是可削减负荷。这个约束相对简单只要限制削减比例不超过总负荷的一个百分比就行。可替代负荷则是综合需求响应的精髓。比如在电价高的时段用户可以从电采暖转为燃气采暖。这种能源替代在模型里需要设置替代系数%% 可替代负荷建模 % 电转热替代量 H_replace sdpvar(1, 24); % 电负荷减少 替代量 * 替代效率 P_load_reduce H_replace / eta_replace; % 气负荷增加 替代量 / 热效率 Gas_load_increase H_replace / eta_gas_heat;整个需求响应模块的目标是让系统在满足用户用能需求的前提下通过负荷侧灵活性降低系统运行成本。3.6 求解与结果整理把所有约束和目标函数都写好后调用求解器求解%% 求解 ops sdpsettings(solver, gurobi, verbose, 2); result optimize([C_all, C_balance, C_device, C_storage, C_dr], Objective, ops); %% 结果提取 if result.problem 0 P_buy_opt value(P_buy); P_gt_opt value(P_gt); Carbon_opt value(Carbon_total); Cost_opt value(Objective); else disp(求解失败); end这里有个比较实用的技巧如果模型求解失败第一步不是急着改代码而是先检查约束的可行性。YALMIP提供了一个非常有用的命令% 检查约束是否可行 [primal_feas, dual_feas] check(C_all);check函数会返回每条约束的可行性残差。如果某条约束残差是负数说明这个约束不可行你就能快速定位问题出在哪一条上。4. 结果分析与核心图表解析4.1 电力平衡图怎么看跑完代码后第一张要画的是电力平衡图。这个图把每个时段的购电、风电、燃气轮机发电、储能放电都画出来叠加成堆叠柱状图再叠加上电负荷曲线。我复现时发现的一个显著特征就是在夜间低谷时段比如0点到6点电价低系统会多购电一方面满足负荷需求另一方面给储能充电。而在白天高峰时段比如10点到14点电价高系统会优先用燃气轮机发电同时储能放电来支撑负荷。如果需求响应起作用了你会在图里看到高峰时段的负荷曲线出现明显的削峰——就是负荷曲线的顶部被削平了一些而低谷时段的负荷曲线则被抬高了这就是可转移负荷从高峰挪到了低谷的表现。4.2 碳成本对比分析碳成本分析是检验阶梯碳机制是否生效的关键。我建议做一个对比在同样的系统配置下分别用固定碳价和阶梯碳价跑一次模型然后对比两种机制下的碳排放量、碳交易成本、总成本。理论上阶梯碳机制下的系统碳排放量应该更低因为排得越多边际成本越高系统会主动减少高碳排放设备如燃气轮机的出力转而用更清洁的电力或者储能来替代。但代价是总成本可能会略微上升因为清洁能源的单位成本更高。这个减排效果和经济成本之间的权衡曲线正是论文里最有价值的图表。我跑出来的结果阶梯碳机制比固定碳价机制大约能降低8%到15%的碳排放量成本上升控制在3%以内。这个量级可以参考但具体数值取决于你的系统参数和数据。4.3 灵敏度分析怎么做做过学术研究的人都知道审稿人最喜欢问的问题就是你的结果对参数敏感吗。所以复现的时候干脆就把灵敏度分析一起做了后面写论文也方便。我建议重点做三个参数的灵敏度分析碳价水平把碳价从低到高扫描一遍看碳排放量和总成本怎么变化需求响应补偿价格补偿价格越高用户越愿意参与响应但系统的补偿成本也越高存在一个最优补偿价格风电装机容量风电装机越大系统碳排放越低但弃风率可能上升灵敏度分析的具体做法就是在MATLAB里加一个循环每次改变一个参数值重新求解一遍模型记录下关键指标最后画成曲线%% 灵敏度分析示例 carbon_price_range [30:10:150]; for i 1:length(carbon_price_range) carbon_prices carbon_price_range(i); % 重新求解模型记录结果 Carbon_emission(i) value(Carbon_total); Total_cost(i) value(Objective); end %% 画图 figure; yyaxis left; plot(carbon_price_range, Carbon_emission, -o); ylabel(碳排放量 (t)); yyaxis right; plot(carbon_price_range, Total_cost, -s); ylabel(总成本 (元)); xlabel(碳价 (元/t));5. 常见问题与排查技巧实录5.1 求解器告诉你Infeasible Problem怎么办这个应该是复现时最劝退新手的错误了。模型不可行说明你写的约束条件之间存在矛盾没有任何一组变量取值能满足所有约束。排查思路我总结了四个步骤第一步先查电量平衡约束。把电量平衡约束的等式右边和等式左边加起来看看是不是有遗漏项。比如你忘了把P2G的耗电加进去那平衡等式就永远不可能成立。第二步查储能约束。储能是整个模型里最容易出矛盾的地方。一天结束时的SOC约束特别容易出问题比如初始SOC是0.2要求的最终SOC也是0.2但储能容量和充放电功率的限制导致不可能在24小时内回到初始SOC那就必崩。第三步查需求响应约束。可转移负荷的总量守恒约束很容易和转移量上下限冲突。如果 ( P_load 100)你设置 ( P_shift_max 60)然后又要求总转移量为0这个模型在只有两个时段的情况下就无解——毕竟不可能把60的负荷从时段1挪到时段2又把60从时段2挪回时段1除非你允许同时转移。第四步用YALMIP的check函数定位这个在上文提过check(C_all)可以告诉你哪条约束不可行快速缩小排查范围。5.2 求解结果很离谱多半是参数出了问题有时候模型能求解成功但结果怎么看怎么不对劲。比如储能一直在充电从不放电或者燃气轮机一直满负荷运转这时候大概率不是模型结构出了错而是某个参数设置不合理。我踩过的一个典型坑是储能效率设置得太低。储能充电效率0.9放电效率0.9来回一趟只有81%的效率再加上储能容量限制在分时电价差不够大的情况下求解器算下来发现峰谷套利根本不划算干脆就不让储能工作了结果储能设备成了一个摆设。解决办法是调整储能效率和峰谷价差确保价差大于储能往返损耗的成本储能才会有用武之地。另一个常见参数坑是燃气轮机的效率。有些论文里的燃气轮机发电效率给到0.4左右但对应的天然气热值单位可能是MJ/m³而你在代码里用的单位是kWh。这个单位换算一旦出错算出来的购气成本可能比实际高十倍结果就是燃气轮机完全不发电全靠购电。5.3 求解速度太慢怎么优化模型规模大、整数变量多求解速度特别慢是常态尤其是加了需求响应后引入了大量0-1变量一跑就是一小时。我建议从三个方面优化第一减少整数变量。如果能用连续变量近似替代0-1变量尽量替换。比如储能互斥约束可以用P_ch * P_dis 0来凑合但这是个非线性约束求解更慢。所以更好的做法是去掉互斥约束直接设置P_ch 0、P_dis 0然后在目标函数里给它一个很小的惩罚项比如 ( 0.001 * (P_ch P_dis) )这样求解器会尽量避免同时充放电但模型的求解速度会快很多。第二简化需求响应建模。如果可转移负荷的分时转移矩阵是稀疏的可以提前算好哪些时段的转移是允许的避免引入过多的转移变量。第三换求解器。如果是学术用途别犹豫直接上求解器。Gurobi在MILP问题上的求解效率比MATLAB内置求解器高出好几个数量级。6. 复现后的扩展思路6.1 从确定性到不确定性现在的模型是确定性的也就是说所有参数负荷、风电、电价都是提前知道的。现实中这些参数都有很强的不确定性所以很多研究会把模型扩展成两阶段鲁棒优化或者场景随机规划。这两种方法都有现成的YALMIP工具箱支持想要深入的话可以从这个方向切入。6.2 从单目标到多目标目前模型是把经济成本和碳排放放在同一个目标函数里用权重系数把碳成本货币化。但碳减排并不完全等价于经济成本所以有些研究会用多目标优化把经济成本和碳排放量作为两个独立目标用NSGA-II这类算法求一个Pareto前沿展示成本和排放之间的权衡。MATLAB的全局优化工具箱自带gamultiobj函数可以试试。6.3 从优化调度到规划优化调度回答的是给定系统配置怎么运行最经济的问题。而系统规划回答的是怎么建设系统未来20年总成本最低的问题。把当前模型扩展成多阶段的规划-运行联合优化是论文往更高层次发的一个常见方向。结尾的个人体会这套代码我前前后后改了大概两周前前后后踩的坑加起来不少于十个。最让我印象深刻的还是储能互斥约束和阶梯碳分段线性化这两个模块看着都是模型里的小细节结果恰恰是它们决定了整个模型的求解质量和结果合理性。我在实际复现中最强烈的一个感受就是读论文的时候一定要把每个公式都亲手推导一遍特别要注意那些由式(15)可得的跳步里面往往藏了建模时最关键的简化假设。你把这层窗户纸捅破了复现就不会是痛苦地抄代码而是有底气地调试、改进和扩展。
