这段时间一直在复盘综合能源系统的优化调度正好前两天帮一个师弟调他的Matlab程序主题就是“计及新能源出力不确定性的电气设备综合能源系统协同优化”。说实话这个题目看起来像是把几个热门词堆在一起但真正动手做的时候每一步都有不少细节值得抠。这篇文章就以这个课题为例把我从建模到代码实现的完整思路、以及踩过的坑都整理出来。先说清楚这个课题到底在做什么构建一个包含电、热、气等多种能源形式耦合的综合能源系统在考虑风电、光伏出力不确定性的前提下通过协同优化各类设备的出力计划实现系统运行成本最低、弃风弃光最少、或者综合能效最高。适合谁来参考正在做综合能源优化调度相关毕设、论文的同学或者刚接触Matlabyalmip工具箱做优化建模的工程师。如果你已经能熟练跑通单能源系统的优化程序这篇文章能帮你把“协同”和“不确定性”这两块补上。1. 为什么“不确定性”和“协同”必须放在同一个框架里1.1 单看新能源不确定性是绕不开的坎风电、光伏的出力本质上是一个随机过程。风速受天气系统影响光照强度随云层变化剧烈这导致新能源机组的出力曲线没法像火电机组那样按计划精确执行。如果你在优化模型里直接把预测曲线当作确定值那实时的功率偏差就只能靠系统备用硬扛极端情况下会出现切负荷或者严重弃风。我在实际测试中发现当风电预测误差超过20%时确定性优化结果在实时运行中的失配成本会急剧上升。这个结论其实不难理解预测曲线给出的最优解在真实出力偏离之后可能根本不可行。1.2 协同优化为什么能部分对冲不确定性综合能源系统最大的优势在于多能互补。电、热、气不是三个孤立的系统而是通过热电联产机组、电锅炉、燃气锅炉、储能装置、甚至电转气设备耦合在一起。举个例子风电大发时电力系统可能面临消纳压力但此时如果电锅炉或者热泵能消纳多余电能转化为热能储存在蓄热罐里那风电的弃电率就能大幅下降。反过来当风电出力骤减时燃气锅炉和热电联产机组可以快速提升出力填补电力和热力的缺口。这种跨能源载体的互助能力实际上是给系统增加了额外的“柔性可调资源”。但前提是你的优化模型必须同时考虑多个能源子系统的运行约束而不是分开优化再拼凑结果。1.3 两者共同决定了优化结果的可用性把不确定性和协同分开处理效果都会打折扣。如果不考虑不确定性协同调度方案在实时运行中可能频繁失效如果只考虑不确定性但不做多能协同系统就少了很多调节手段优化结果会偏保守、经济性差。所以这个课题的价值就在于用协同优化的框架扩大可行域用不确定性建模提高解的鲁棒性两者结合才能在保证安全的前提下追求经济性。这也是为什么这个题目在最近几年的论文里出现频率这么高。2. 不确定性建模方式的选择场景法、鲁棒法与区间法的取舍2.1 三种主流方法的建模逻辑处理新能源出力不确定性学术界和工程上主要用三类方法场景法随机规划通过采样生成大量可能的新能源出力场景在目标函数中求期望成本。这种方法思路直观能直接给出概率意义下的最优解但计算量随场景数量线性增长。鲁棒优化只考虑不确定量落在某个集合内的情况优化方案需要在集合内最恶劣情形下仍然可行。这种方法求出的解偏保守但安全边际很高适合系统备用紧张、安全要求极高的场景。区间优化把不确定性表示成区间求解模型给出乐观解和悲观解。计算量介于前两者之间但结果是一个区间还需要决策者根据风险偏好再做权衡。在综合能源系统里我建议优先考虑场景法。原因很简单系统的设备约束本来就多非线性环节也多再用鲁棒优化去嵌套“最恶劣场景”模型容易变得又大又硬Matlab求解时间可能会从几分钟膨胀到几十分钟甚至不收敛。而场景法通过采样、削减、加权能落在混合整数线性规划的框架内用yalmip配合CPLEX或Gurobi能高效求解。2.2 场景生成与削减的实操思路场景生成最常见的做法是蒙特卡洛采样。假设风电预测误差服从正态分布预测值为P_forecast误差标准差取预测值的10%~15%那么每个场景的风电出力可以这样生成P_scenario P_forecast sigma * randn(N_scenarios, T);randn生成的随机数可能造成出力为负或者超过装机容量需要加边界约束P_scenario max(0, min(P_capacity, P_scenario));不过如果直接用几千个场景喂给优化模型变量和约束数量会爆炸。这时候要做场景削减。常见的削减算法是快速前向选择法和同步回代消除法核心思路是挑出概率大、代表性强的场景把相似度高的场景合并。我自己写代码时习惯用同步回代消除的简化版1. 计算所有场景两两之间的欧氏距离每个场景是一条时序曲线 2. 找到距离最近的一对场景合并其中概率较小者将其概率累加到另一个场景上 3. 重复以上步骤直到场景数量缩减到目标值这个算法代码量不大但效果很好。把1000个初始场景削减到30个左右优化结果和用全场景计算的期望成本偏差一般能控制在2%以内。2.3 如何选择场景数量场景数量不是越大越好。我做过一组对照实验10个场景下优化结果波动很大20个场景基本稳定30个场景再往上增加目标函数值的变化几乎可以忽略但求解时间明显变长。建议初学者先削减到20~30个场景把模型跑通后再逐步增加观察目标函数值的变化程度。如果20个到30个的结果变化很小说明场景收敛性已经满足要求。3. 综合能源系统协同优化的核心模型怎么搭3.1 确定系统架构与设备模型建模第一步是确定系统包含哪些设备。常见的架构是风电机组、光伏机组、热电联产机组CHP、燃气锅炉、电锅炉、储能电池、蓄热罐、电负荷、热负荷部分系统还有电转气设备。每个设备都要建立运行模型以典型的CHP机组为例电功率输出P_chp eta_e * F_gas热功率输出H_chp eta_h * F_gas其中F_gas为机组消耗的燃气功率eta_e为发电效率eta_h为供热效率。CHP的可行运行域通常用几个线性不等式约束来近似避免非线性。储能电池模型荷电状态递推SOC(t1) SOC(t) (P_charge(t) * eta_c - P_discharge(t) / eta_d) * delta_t充放电功率上下限约束、SOC上下限约束。同一个时刻不能同时充放电这是一组常见的整数约束。蓄热罐模型和电池类似只是能量载体是热能约束形式基本相同。3.2 目标函数设计与惩罚项目标函数以系统总运行成本最低为主流包括购电成本从上级电网购电购气成本天然气购买设备运行维护成本弃风弃光惩罚成本失负荷惩罚成本其中弃风弃光惩罚项设计很关键。如果你把弃风弃光成本设成0优化模型为了降低成本会优先让新能源出力压到最低表面上成本降低了实际弃风弃光惨不忍睹。我一般把弃风弃光的惩罚成本系数设得比常规机组发电成本高2~3倍迫使模型尽可能消纳新能源。目标函数统一写成min C_total sum_t( C_grid(t) * P_grid(t) ) sum_t( C_gas * F_gas_total(t) ) sum_t( sum_dev( c_om * P_dev(t) ) ) sum_t( C_waste * (P_wind_forecast(t) - P_wind_use(t)) ) sum_t( C_loss * L_loss(t) )3.3 耦合约束的数学表达协同优化的“协同”二字主要体现在耦合约束上。如果没有耦合约束各个能源系统就是独立的谈不上协同。以CHP为耦合点电力平衡方程需要包含CHP的发电项P_grid(t) P_wind(t) P_pv(t) P_chp_e(t) P_battery_discharge(t) P_load(t) P_eb_consume(t) P_battery_charge(t) P_waste(t)热力平衡方程H_chp(t) H_gas_boiler(t) H_tank_discharge(t) H_eb(t) H_load(t) H_tank_charge(t) H_waste_heat(t)其中H_eb(t)是电锅炉产热从而将电力系统和热力系统耦合起来。天然气平衡方程同理燃气量分配给CHP和燃气锅炉。电转气设备如果考虑进来可能在气网侧多一个耦合变量但会增加不少计算量初期不建议加。3.4 时间粒度和调度周期综合能源系统调度通常取24小时为周期每小时一个时段也就是T24。热力系统的惯性比较大小时级是合适的如果要考虑日内滚动修正可以取15分钟一个时段但规模会增大四倍求解压力也大不少。4. Matlab代码实现从数据准备到求解器调用的完整流程4.1 代码整体架构我习惯把程序拆成几个文件避免一个脚本堆到底main.m主函数控制整个流程data_input.m负荷曲线、新能源预测曲线、设备参数的输入scenario_generate.m生成新能源出力场景scenario_reduce.m场景削减build_model.m搭建优化模型run_optimization.m调用求解器求解plot_results.m可视化结果这样的结构好处很明显改参数、换场景、换求解器都只需要动对应的模块调bug时不用在几百行代码里翻来翻去。4.2 关键代码示例场景生成与削减% scenario_generate.m N_scenarios 1000; mu P_wind_forecast; % 预测值 sigma 0.12 * P_wind_forecast; % 标准差取12% scenarios nan(N_scenarios, T); for i 1:N_scenarios scenarios(i, :) mu sigma .* randn(1, T); scenarios(i, :) min(max(scenarios(i, :), 0), P_wind_capacity); end% scenario_reduce.m function [scen_reduced, prob_reduced] scenario_reduce(scenarios, N_target) % 同步回代消除法做场景削减 % scenarios: n*T 矩阵 % N_target: 目标场景数 n size(scenarios, 1); prob ones(n, 1) / n; while n N_target % 计算所有场景对儿的距离 dist nan(n, n); for i 1:n for j i1:n d sqrt(sum((scenarios(i, :) - scenarios(j, :)).^2)); dist(i, j) d; dist(j, i) d; end end % 找到距离最小的场景对 dist(dist 0) inf; flat_idx find(dist min(dist(:)), 1); [i, j] ind2sub([n, n], flat_idx); % 去掉概率小的场景 if prob(i) prob(j) prob(i) prob(i) prob(j); scenarios(j, :) []; prob(j) []; else prob(j) prob(j) prob(i); scenarios(i, :) []; prob(i) []; end n n - 1; dist []; end scen_reduced scenarios; prob_reduced prob; end4.3 yalmip建模与CPLEX/Gurobi求解Matlab里做优化建模强烈推荐用yalmip工具箱。它可以把优化问题描述得和数学表达式非常接近代码可读性好调试起来省心。% build_model.m yalmip(clear); % 决策变量 P_grid sdpvar(T, 1); % 向上级电网购电功率 P_chp sdpvar(T, 1); % CHP发电功率 H_chp sdpvar(T, 1); % CHP热功率 P_eb sdpvar(T, 1); % 电锅炉耗电功率 H_gb sdpvar(T, 1); % 燃气锅炉热功率 SOC_b sdpvar(T, 1); % 电池荷电状态 P_bc sdpvar(T, 1); % 电池充电功率 P_bd sdpvar(T, 1); % 电池放电功率 u_bc binvar(T, 1); % 电池充电状态 u_bd binvar(T, 1); % 电池放电状态 % ... 其余变量类似 % 目标函数 objective sum(C_grid_price .* P_grid) ... sum(C_gas * (P_chp H_chp / 0.8) / 0.35) ... sum(P_eb * 0.02 H_gb * 0.015) ... sum(C_waste * (P_wind_cut P_pv_cut)); % 约束条件 Constraints []; % 电力平衡约束 Constraints [Constraints, P_grid sum(wind_scen_temp, 3) P_pv P_chp P_bd ... P_load P_eb P_bc P_wind_cut]; % ... 其他约束同理 % 求解 options sdpsettings(solver, cplex, verbose, 2, savesolveroutput, 1); optimize(Constraints, objective, options);这里有个细节需要注意如果要在目标函数里考量多个场景的期望成本通常的做法是对每个场景分别写平衡约束目标函数中取各场景成本的期望值。变量会有场景编号维比如sdpvar(T, N_scen)。如果场景数一多变量规模就上来了这也是为什么场景削减很重要。4.4 没有商业求解器怎么办CPLEX和Gurobi都是商业软件学术用户可以通过学校申请免费许可证。如果你暂时拿不到可以用yalmip直接指定内置的cplex、gurobi如果都没有可以用sedumi或intlinprog。但对于混合整数规划问题intlinprog只能解决规模较小的问题。我的建议是先跑一个把小模型确保整体逻辑正确再上全套场景求解。不要一上来就怼1000个场景遇到数值问题会非常痛苦。5. 典型结果的解读不确定性带来的经济账5.1 对比确定性模型与不确定模型的成本差异我跑过的算例中确定性模型直接使用预测曲线算出的日运行成本往往看起来比场景法要低3%~5%。但这不是最优解的优势而是因为它“假装”新能源出力完全已知把所有偏差都用理想化手段消纳掉了。真实运行中误差一旦兑现系统需要额外调用备用资源这部分成本确定性模型根本没计入。场景法算出的成本里头本来就包含了误差引起的一部分预期调整成本两者的实际运行成本差距并没有表格上那么大。这一点在论文写作和项目汇报中需要特别说明否则容易被质疑“为什么不直接用确定性模型”。5.2 场景数量对成本和求解时间的影响还是用我前面提的对照实验来说明场景数量期望运行成本相对值求解时间51.082.3秒101.035.1秒201.0111.8秒301.00521.4秒501.00055.0秒可以看到10个场景是性价比最高的点30个场景之后收益很小但求解时间增长明显。这个规律在大多数综合能源系统模型中都成立可以作为初设参数时的参考。5.3 灵敏度分析误差标准差的影响风电预测误差标准差从5%调到30%优化结果会有什么变化理论上误差越大系统越需要依赖备用容量和储能调节运行成本越高。当标准差达到20%以上模型可能会主动减少风电的利用权重转而增加燃气轮机出力以保证功率平衡的可靠性。这种趋势用图形展示时X轴取预测误差标准差Y轴取系统总成本和弃风率曲线会呈现明显的非线性增长。我建议把所有灵敏度分析的结果都做成图表不仅方便写论文也有助于自己理解模型行为是否合理。5.4 收敛性与最优性的观察用CPLEX或者Gurobi求解时如果模型规模较大可以设置sdpsettings(cplex.mip.tolerances.mipgap, 1e-3)让求解器在相对最优间隙达到0.1%时停止。这样能在保证解质量的前提下显著缩短时间。判断求解是否正常要关注求解器输出的gap曲线如果gap下降平缓且最终到不了0附近多半是约束写错导致模型不可行或数值病态需要回头排查变量维度和单位一致性。6. 协同优化模型在调试中容易踩的坑6.1 时间索引错位Matlab数组索引从1开始而优化模型的时段从1到24这个没有问题。但要注意储能的递推方程里SOC(t1)在t24时会越界。我习惯把SOC的变量定义成T1维最后一个值作为终态SOC约束让它等于初态SOC。这样既满足循环约束又避免越界。6.2 电、热、气功率单位混用不同能源载体如果用不同的单位体系比如电力用MW天然气用m3/h热量用GJ/h目标函数里系数换算错的概率极高。我的建议是全部统一成MW或kW天然气热值按9.7 kWh/m3折算热量直接用MWh。这样所有变量的量级一致目标函数数值才稳定求解器也不容易出现数值困难。6.3 整数变量过多导致求解时间爆炸充放电状态、机组启停状态都是0-1变量。如果每个设备都加上启停状态24个时段、30个场景变量数量很容易突破上万。其中场景维和整数变量同时出现时求解时间指数级上升。我的解决思路是先别加启停约束只做连续功率优化。跑通后再逐步添加整数变量逐项检查求解时间变化。如果某个整数约束让求解时间从几秒变成几分钟就需要考虑它有没有更紧凑的线性化表达方式。6.4 电锅炉产热与热负荷平衡不同步电锅炉可以把电转成热但热负荷高峰期可能和风电大发期错开。这时候如果没有蓄热罐就会出现“电不够热、热不够电”的尴尬局面。我在调试初期忽略了蓄热罐的时序耦合结果优化结果里电锅炉几乎不启动系统全靠燃气锅炉供热风电弃电率居高不下。加蓄热罐之后电锅炉开始参与调节弃风率下降约15个百分点。这说明协同优化的价值不是设备越多越好而是设备之间的时间耦合和能量转换关系要匹配。蓄热罐、储电池这类时序转移设备在多能互补系统里承担的角色比单能源系统重要得多。6.5 边界约束写反风电出力场景生成之后记得要修正边界。randn可能出现负数如果不做max(0, ...)处理模型里就会出现“负的出力供应”求解器会用它来“吸收”功率导致结果完全失真。这种错误不容易发现因为目标函数值看起来正常但功率平衡约束实际上是靠负风电在满足的。每生成一批新场景先画几条典型曲线看看形状是否合理。肉眼扫一眼比任何测试代码都快。7. 进一步扩展的方向与我的个人体会这个项目跑通之后可以往几个方向继续深挖。一是把不确定性建模从场景法换成分布鲁棒优化用风电出力的矩信息和置信集合描述不确定性在解的保守性和计算复杂度之间找到新的平衡点。这个方向比较前沿适合论文创新点充实。二是把协同优化从日前调度扩展到日内滚动修正。不确定性是时变的提前24小时做好的计划需要在每个滚动时段基于最新预测结果重新优化形成一个反馈闭环。这种模型更贴近实际控制逻辑但程序结构和数据流复杂度会明显增加。三是考虑碳交易机制。在目标函数中引入碳配额和碳交易价格综合能源系统的低碳性和经济性就可能出现新矛盾优化结果也会随之改变。我在实际做这个项目的过程中最大的体会是模型的复杂度一定要匹配你的求解资源和对结果可靠性的要求。刚开始做的时候总觉得考虑的因素越多越好设备加了一堆约束写了上百条最后求解器算不动或者结果出来后发现大部分设备在最优解里根本没怎么用。反而是简化后的核心耦合结构把电、热、气三者的互动关系吃透在Matlab里更容易跑出有说服力的结果。还有一个小技巧在提交长任务求解前先用一个1~2天的简单算例做正确性校验再把维度扩到24小时、多场景。很多单位不一致、维度不匹配的问题在小模型中一眼就能看出来。如果你最近正在做类似的综合能源系统模型希望这篇文章能帮你少走几步弯路。尤其是场景法那套流程从场景生成到同步回代消除再配合yalmip调用CPLEX基本是当前学术界和工程界最通用也最容易上手的组合。照着这个思路一步步实现很快就能跑出自己的第一版协同优化结果。