先说一下背景。我最近在研究园区级综合能源系统的日前调度问题一开始天真地以为只要把光伏和风电的预测曲线当成确定值塞进优化模型算出来的就是能直接用的调度方案。结果一到实际运行就露馅了——某天光伏出力午间比我预测的高了30%电储能直接被充满还被迫弃光CHP机组照常满发多了的电又卖不掉整体运行成本比仿真结果高出快两成。后来我意识到综合能源系统协同优化真正难的地方不是把电、气、热多能流模型搭起来而是怎么处理新能源出力不确定性。这篇博客就围绕我用Matlab做的一套计及新能源出力不确定性的综合能源系统协同优化代码把建模思路、代码架构、求解过程和踩坑记录完整梳理一遍希望对正在做相关方向毕业设计或科研课题的朋友有参考价值。1. 为什么不确定性是综合能源系统协同优化的第一道坎1.1 一个实际才能体会到的痛点预测误差如何击穿调度方案光伏和风电的出力预测不管用数值天气预报还是统计学习方法误差都客观存在。我做过一段时间的实测数据统计在典型日尺度上光伏功率预测的平均绝对百分比误差能做到10%~15%就算不错遇到多云或突变天气单点误差飙升到30%以上很常见。风电更夸张风速的随机性导致出力波动范围很大尤其在阵风时段15分钟级出力可以从额定功率的20%跳到80%。如果优化模型把预测值当真值处理那么约束条件里的功率平衡方程在真实运行时大概率被打破。解决办法是实时调整一调整就牵一发动全身CHP机组爬坡跟不上、储能SOC被提前耗尽、电锅炉供热与热负荷错配。我印象最深的一次仿真确定性调度方案给出的购电计划是平滑的阶梯曲线但实际执行时几乎每隔一个时段就要重新修正最终经济性还不如一个保守的固定方案。这就是不确定性带来的隐形成本它不会直接出现在确定性模型的公式里却真实发生在每一次运行决策中。所以处理不确定性的本质不是追求一个精确的最优解而是追求一个在各种可能出力场景下都不太差、且经济性可接受的稳健解。1.2 不确定性建模的四种主流思路对比在综合能源系统优化领域处理新能源出力不确定性的主流方法大致有四类我用了很长一段时间逐一调研和测试这里直接给出对比结论。建模方法基本思想优点缺点典型应用场景随机规划场景法用离散场景集合近似连续概率分布模型直观能给出兼顾经济性的期望最优解场景数量大时求解规模膨胀缩减不当会丢失信息日前调度、机组组合鲁棒优化用不确定集合描述出力波动范围保证最坏情况下可行决策保守但可靠性极高求解效率高结果偏保守经济性偏差极端天气安全校核区间优化用区间数直接描述出力上下界只需上下界数据建模简单不能体现区间内概率分布信息数据匮乏场景模糊机会约束用模糊隶属度函数描述不确定性可以在一定置信水平下放松约束隶属度函数选取主观性强中长期规划我的实际选择是随机场景法为主结合同步回代缩减技术。原因是综合能源系统协同优化本身已经包含电、气、热多种异质能流的耦合约束鲁棒优化虽然求解快但所有不确定约束都按最坏情况去卡CHP机组出力和储能充放电会被严重限制经济性损失太大而场景法能比较精细地反映风电光伏出力的概率分布特性配合场景缩减技术可以用较少的典型场景逼近原始分布兼顾精度和计算效率。1.3 场景法落地时要回答的三个问题用场景法做不确定性优化代码实现前必须想清楚三件事场景从哪里来通常做法是用历史数据拟合风速的Weibull分布和光照强度的Beta分布再通过蒙特卡洛抽样生成大量初始场景。如果没有历史数据也可以用典型日出力曲线叠加正态分布噪声来构造。我在项目里采用的是后一种思路参数设置时参考了当地气象站多年平均数据保证抽样结果的统计特征符合常识。场景缩减到什么程度初始场景我取了1000个缩减到20个。缩减后的典型场景应该能保留原始场景集的均值、方差和相关性特征。这个数量是权衡求解时间和精度后定的如果你发现缩减后的场景间差异不明显可以把目标场景数提高到30或50。概率怎么赋场景缩减完成后每个典型场景需要分配一个概率值所有场景概率之和为1。目标函数中需要对各场景下的运行成本做期望运算也就是每个场景的成本乘以对应概率再累加。这一步如果写错求解结果会直接偏掉。2. 系统建模把电-气-热耦合关系写成数学语言2.1 能源集线器视角下的系统架构我构建的综合能源系统模型核心框架采用能源集线器Energy Hub思想。简单说变电站、燃气轮机组、电锅炉、燃气锅炉、储能装置这些电气设备和能源转换设备共同构成一个多输入多输出的能量转换单元。输入端是外部电网购入电力和天然气网购入燃气外加风光新能源出力输出端满足电负荷、热负荷两类刚性需求。能源集线器的输入输出关系可以写成矩阵形式[L_e] [η_ee η_ge] [P_e] [L_h] [η_eh η_gh] * [P_g]其中L_e和L_h分别是电、热负荷功率P_e和P_g分别是购入电能和购入燃气功率η矩阵中的元素代表各条能量转换路径的效率。这个模型的可贵之处在于它把复杂的多能耦合关系压缩成矩阵运算后续无论是写约束还是做灵敏度分析都非常方便。我在Matlab里就是按照这个思路先搭好能量转换系数矩阵再写各个设备的详细约束。2.2 关键设备建模与运行约束CHP热电联产机组是综合能源系统的核心设备它的特点是电出力P_chp和热出力Q_chp之间存在耦合关系。我用的是抽凝式机组线性化模型Q_chp(t) a * P_chp(t) b P_chp_min ≤ P_chp(t) ≤ P_chp_max -ΔP_chp ≤ P_chp(t1) - P_chp(t) ≤ ΔP_chpa是热电比斜率b是常数项。要注意的是CHP机组在低负荷时热电比会变化线性模型只能在某个工况区间内保持较高精度。如果实际项目中机组长期在部分负荷运行建议用分段线性化处理把可行域划分成多个线性段。电锅炉和燃气锅炉建模相对简单Q_eb(t) η_eb * P_eb(t), 0 ≤ P_eb(t) ≤ P_eb_max Q_gb(t) η_gb * F_gb(t), 0 ≤ F_gb(t) ≤ F_gb_max其中P_eb是电锅炉消耗的电功率F_gb是燃气锅炉消耗的燃气功率η是各自热效率。电锅炉的价值在于消纳富余风电——在风电大发时段如果系统向上级电网反送电不经济就把多余电能转化为热能储存或直接供热。储能系统我同时考虑了电储能和热储能模型统一写成SOC_e(t1) SOC_e(t) (η_ch * P_ch(t) - P_dis(t) / η_dis) * Δt SOC_min ≤ SOC_e(t) ≤ SOC_max 0 ≤ P_ch(t) ≤ P_ch_max * u_ch(t) 0 ≤ P_dis(t) ≤ P_dis_max * u_dis(t) u_ch(t) u_dis(t) ≤ 1u_ch和u_dis是充放电状态二进制变量这一步引入了整数变量模型变成MILP问题。这也是为什么求解器需要选择Cplex或Gurobi这类商业MILP求解器的主要原因。电储能主要用来削峰填谷和应对不确定性带来的实时波动热储能则配合CHP和电锅炉一起平滑热负荷曲线。电气设备层面的网络约束我没有把配电网三相潮流全部写进去因为综合考虑电-气-热耦合和不确定性场景后非线性潮流会严重拖慢求解速度。实际处理是在园区级尺度下采用功率平衡约束替代详细潮流并给变压器、线路设置了最大传输容量约束。如果你的课题更侧重配电网层面可以把DistFlow线性化潮流加进来Matlab里用Yalmip结合Cplex也能处理但求解时间会明显上升。2.3 目标函数与完整约束体系协同优化的目标是最小化系统总运行成本包括购电费用、购气费用、设备运维费用和弃风弃光惩罚min Σ_s ρ_s * { Σ_t [ c_e(t)·P_grid(t) c_g·F_gas(t) c_om·(P_chp(t)Q_chp(t)) λ·(P_w_avail(t)-P_w_use(t)) ] }其中ρ_s是场景s的概率c_e(t)是分时电价c_g是气价c_om是运维成本系数λ是弃风弃光惩罚系数P_w_avail是新能源可用出力P_w_use是实际消纳出力。等式约束主要是各时段的电功率平衡和热功率平衡电平衡P_grid(t) P_chp(t) P_w_use(t) P_dis_e(t) P_load(t) P_eb(t) P_ch_e(t) 热平衡Q_chp(t) Q_eb(t) Q_gb(t) P_dis_h(t) H_load(t) P_ch_h(t)不等式约束包括设备出力上下限、爬坡约束、储能SOC范围、充放电状态互斥约束、上级电网交互功率限制等。这些约束全部写成矩阵形式后用Yalmip的sdpvar和constraint组合非常方便这也是我推荐用Yalmip而不是纯Matlab手写求解器接口的原因——代码可读性和可修改性都好很多。3. Matlab代码架构从场景生成到求解器的一站式实现3.1 整体文件结构与调度流程收到不少同学的私信问这种代码从哪下手我给的统一建议是先画清楚数据流。整套Matlab程序我拆成了五个脚本文件职责边界非常清晰IES_Optimization/ ├── main_IES.m # 主程序串起整个求解流程 ├── params_IES.m # 系统参数与设备参数 ├── generate_scenarios.m # 新能源出力不确定性场景生成与缩减 ├── build_model.m # 用Yalmip构建优化模型 └── plot_results.m # 结果可视化与对比分析主程序的执行顺序是先运行params_IES.m加载参数再调用generate_scenarios.m生成典型场景及概率接着进入build_model.m构建决策变量、约束和目标函数调用求解器求解后把结果传给plot_results.m画图。有个容易忽视的点所有脚本都不建议设置成函数形式来调用参数结构体因为调试时你经常需要在工作区里查看中间变量。我直接把params_IES.m写成脚本里面的变量以params结构体形式存放例如params.chp.a、params.chp.p_max。这样在主程序里任何位置都能快速查看和修改参数迭代效率高很多。3.2 不确定性场景生成蒙特卡洛抽样与同步回代缩减场景生成是整套代码里最体现不确定性建模精髓的部分。我采用典型出力曲线 随机扰动的抽样方式先定义风电和光伏的预测典型曲线假设预测误差服从正态分布在每个时段叠加随机噪声生成1000个初始场景。% 参数初始场景数 N_scen 1000; T 24; % 风电典型出力曲线标幺值 wind_base [0.28 0.26 0.24 0.22 0.21 0.23 0.25 0.30 0.38 0.45 0.42 ... 0.36 0.33 0.35 0.30 0.26 0.24 0.28 0.35 0.42 0.48 0.45 0.38 0.30]; % 光伏典型出力曲线标幺值 pv_base [0 0 0 0 0 0.05 0.15 0.30 0.55 0.75 0.88 0.95 0.92 0.80 ... 0.62 0.40 0.20 0.08 0 0 0 0 0 0]; % 抽样生成初始场景矩阵每一行是一个完整日场景 wind_scen repmat(wind_base, N_scen, 1) .* (1 0.15 * randn(N_scen, T)); pv_scen repmat(pv_base, N_scen, 1) .* (1 0.10 * randn(N_scen, T)); % 修正越界值 wind_scen(wind_scen 0) 0; wind_scen(wind_scen 1) 1; pv_scen(pv_scen 0) 0; pv_scen(pv_scen 1) 1; % 组装场景矩阵把风、光场景拼接成1000×48的矩阵 scen_all [wind_scen, pv_scen]; prob_all ones(N_scen, 1) / N_scen;随机生成的场景数量太多直接带入优化模型会导致变量维度过大必须做场景缩减。我采用的是工程上最常用的同步回代缩减方法每次迭代找一对概率距离最小的场景删掉其中一个同时把被删场景的概率加到距离最近的那个场景上直到剩余场景数达到目标值。function [scen_red, prob_red] sbr_scenario_reduction(scen, prob, K) % scen: 原始场景矩阵, N×M % prob: 场景概率列向量, N×1 % K: 目标保留场景数 N size(scen, 1); while N K % 计算场景间欧氏距离矩阵 D pdist2(scen, scen); D(D 0) inf; % 自身距离设为无穷 % 每个场景找最近的相邻场景 [d_min, j_idx] min(D, [], 2); % 概率加权距离作为删除指标 cost prob .* d_min; % 找到删除代价最小的场景 [~, i_del] min(cost); % 将其概率累加到最近场景上并删除该场景 prob(j_idx(i_del)) prob(j_idx(i_del)) prob(i_del); scen(i_del, :) []; prob(i_del) []; N N - 1; end scen_red scen; prob_red prob; end这段代码虽然短但有几点需要额外留意。第一pdist2在场景数很大时内存占用可观如果你的机器配置一般建议用分块计算或直接写双重循环。第二每次删除场景后要重新计算距离矩阵所以缩减20个场景可能要迭代980次程序会跑一小段时间这是正常的。第三缩减结果的稳定性和初始场景有关建议固定随机种子rng(42)保证实验可复现。缩减之后我得到的典型场景可以直观理解为晴天高光伏阴天低光伏大风夜无风夜这几类典型天气模式每一类带上一个概率权重。这个概率分布会直接影响优化结果中储能充放电策略和CHP出力水平。3.3 Yalmip建模与求解器调用模型构建部分我完全依赖Yalmip工具箱它最大的好处是让建模语法和求解器解耦——你不用关心Cplex底层的C API怎么调用只需要用sdpvar定义决策变量、用constraint写约束、用optimize求解。下面给出核心建模代码的骨架。%% 决策变量定义 P_grid sdpvar(T, 1); % 购电功率 P_chp sdpvar(T, 1); % CHP电出力 Q_chp sdpvar(T, 1); % CHP热出力 P_eb sdpvar(T, 1); % 电锅炉电功率 P_w_use sdpvar(T, 1); % 实际消纳风电功率 P_ch_e sdpvar(T, 1); % 电储能充电功率 P_dis_e sdpvar(T, 1); % 电储能放电功率 SOC_e sdpvar(T 1, 1); % 电储能荷电状态 u_ch_e binvar(T, 1); % 充电状态指示 u_dis_e binvar(T, 1); % 放电状态指示 %% 约束集合 Constraints []; %% 功率平衡约束 for t 1:T Constraints [Constraints, P_grid(t) P_chp(t) P_w_use(t) P_dis_e(t) ... P_load(t) P_eb(t) P_ch_e(t)]; end %% 设备运行约束 Constraints [Constraints, params.chp.p_min P_chp params.chp.p_max, Q_chp params.chp.a * P_chp params.chp.b, P_eb_min P_eb P_eb_max, P_w_use 0, P_w_use sum(wind_PV_scenario, 2)]; % 注意上面的风光伏场景是缩减后的某个典型场景需要遍历所有场景实际代码里因为我们要对20个场景同时建模每个场景下都有一组决策变量和平衡约束目标函数里对场景做期望运算。我习惯用三维数组来存放场景相关的决策变量例如P_w_use(s, t)表示场景s下t时段的风电消纳功率P_grid(s, t)表示场景s下t时段购电功率。这样目标函数写成Objective 0; for s 1:N_red Obj_scen 0; for t 1:T Obj_scen Obj_scen ... c_e(t) * P_grid(s, t) ... c_gas * F_gas(s, t) ... c_om * (P_chp(s, t) Q_chp(s, t)) ... lambda * (P_w_avail(s, t) - P_w_use(s, t)); end Objective Objective prob_red(s) * Obj_scen; end求解调用ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); ops.cplex.mip.tolerances.mipgap 0.01; % 设置MIP间隙 result optimize(Constraints, Objective, ops);这里有个细节P_w_avail(s,t)是场景缩减后已知的风光最大可用功率它在优化前就是常数P_w_use(s,t)是决策变量。弃风弃光量的本质就是两者差值。惩罚系数λ设置很关键设置太小会导致模型主动弃风设置太大则牺牲经济性保消纳。我经过试算后把λ设为上网电价的两倍左右效果比较平衡。4. 算例验证与结果分析确定性调度和随机调度的差距有多大4.1 算例系统的参数设定我用一个典型园区综合能源系统做验证系统包含一台CHP机组额定电功率300kW、一台电锅炉200kW、一台燃气锅炉200kW、电储能200kWh/100kW、热储能100kWh/50kW配置风电150kW、光伏100kW。电负荷峰值约350kW热负荷峰值约180kW。参数类别具体参数数值分时电价峰时段 10:00-16:00, 19:00-22:000.85 元/kWh分时电价平时段 08:00-10:00, 16:00-19:00, 22:00-24:000.55 元/kWh分时电价谷时段 00:00-08:000.32 元/kWh天然气价格折算热值成本0.62 元/kWh弃风弃光惩罚单位惩罚成本1.60 元/kWhCHP机组热电比斜率 a / 截距 b1.2 / 20 kW电锅炉效率η_eb0.95电储能初始SOC / 容量0.5 / 200 kWh风光预测典型曲线就是上一节里展示的标幺值曲线乘以各自的额定容量即可得到功率量值。场景缩减后保留了20个典型场景每个场景都包含24个时段的风电、光伏可用出力。4.2 三种调度模式的结果对比我设计了三种模式做对比模式A确定性调度。把风光预测均值当成确定值完全不考虑不确定性。模式B随机场景调度本文方案。用缩减后的20个场景做期望优化。模式C鲁棒调度。用盒式不确定集合不确定区间取预测值的±15%最坏情况约束。调度模式总运行成本元/日弃风弃光率CHP日发电量kWh求解时间s确定性A87265.3%31206.8随机场景B91451.2%336038.5鲁棒C94020.4%348552.1表格里的数虽然只是我的算例结果但反映的规律具有普遍性确定性方案成本最低但这是建立在预测完全准确的假设上实际运行时这个成本根本守不住随机场景方案成本略高但优势在于它对各种可能场景都做了准备实际运行的可靠性大幅提升鲁棒方案最保守适用于极端天气频发、安全要求极高的场景。4.3 结果剖析不确定性成本与调度策略差异深入看CHP机组出力曲线可以发现模式A下CHP机组午间出力会压得比较低因为光伏出力预测较高系统认为不需要CHP多发电。但在随机场景调度中模型意识到午间光伏可能比预测低30%以上因此会让CHP保持一个相对高的出力水平同时让电储能提前预留一点容量——这就是协同优化体现的地方多能互补的灵活性和设备间的协调动作共同对冲新能源不确定性带来的风险。另外弃风弃光率从5.3%降到1.2%主要贡献来自电锅炉。在随机场景中部分场景下夜间风电出力远超预测此时电锅炉启动把多余风电转化为热能储存在热储能中替代一部分燃气锅炉出力。这种电气设备电锅炉、储能与热力负荷之间的协同调度正是综合能源系统相比单一电力系统调度的优势所在。我在结果分析时还计算了一个指标叫不确定性成本——模式B和模式A的成本差值除以模式A的成本约4.8%。这个数字如果换算到全年运行就是一笔相当可观的费用。做项目汇报时把这个指标讲清楚比单纯罗列成本数值更有说服力。5. 我在调试这套Matlab程序时踩过的坑5.1 场景缩减的尺度效应缩得太多反而丢信息一开始我把1000个场景缩减到5个求解确实飞快但优化结果和真实情况偏差很大。原因在于5个场景无法覆盖风光出力的相关性结构——比如光伏高风电低和光伏低风电高这两种关键模式缩减后可能只剩一种模型就会偏向某种特定天气失去代表性。后来我做了敏感性分析目标场景数从5、10、15、20、30一路往上试发现20个场景是一个明显的拐点再往上增加场景数成本改善幅度很小但求解时间增长明显。这个平衡点需要针对你的具体系统参数去试不能拍脑袋定。5.2 Yalmip求解时最容易翻车的三个细节第一二进制变量初始化问题。binvar定义的变量如果没给初值Cplex在某些版本下可能陷入很慢的branch and bound过程。建议在sdpsettings里设置好MIP gap比如0.01让求解器在精度达标后提前退出而不是硬求到全局最优。第二约束中的数值尺度。如果某个约束里的系数相差几个数量级比如电价是0.3而CHP容量是300000W求解器数值稳定性会变差。经验做法是统一量纲——功率用kW、能量用kWh、费用用元把数值控制在0.01~10000范围内Cplex处理起来会很舒服。第三等式约束里的温度、SOC传递关系要特别注意时段的错位。储能SOC约束如果写成SOC(1)和SOC(T1)的循环闭合容易出现索引不一致导致的约束错乱。我调试时有一个晚上所有方案都退化最后发现是SOC_e(t1)和SOC_e(t)的索引写反了这种低级错误排查起来特别费时间。5.3 CHP非线性特性的线性化处理前面提到CHP电热耦合约束Q a*P b这个线性化在小范围内还好但如果CHP机组的可行域是典型的三角形或四边形区域由最小电出力、最大电出力、最大热出力、最小热电比等围成只用一条直线约束会丢失可行域边界信息导致优化结果出现在物理上不可行的工况点。处理办法有两种。一种是把可行域拆成多个线性不等式组合画出四个顶点然后用凸包约束表达另一种是采用二进制变量选择工况区间做分段线性化。我最终采用的是第一种因为代码实现简单只需要加三四个线性不等式约束求解速度没有明显下降。如果你自己要在Matlab里验证CHP可行域建议先画出P-Q平面上的可行域多边形再反推线性约束而不是凭空写公式。5.4 求解时间与精度平衡的实战调参经验这套模型求解时间是38秒左右对于学术研究完全够用。但如果做多日滚动优化或需要嵌入实时控制这个时间就需要压缩。我的几条经验控制整数变量数量。储能充放状态两对二进制变量、CHP启停两对二进制变量24时段就是接近100个整数变量规模不大但会影响求解速度。如果你不需要模拟启停过程可以把启停二进制变量去掉只保留充放电互斥约束MILP退化成QP求解时间能缩短一半。设置合理的MIP gap。学术研究里gap设1%足够说明问题没必要追求0误差。减少场景数的同时用概率密度加权。有些场景虽然发生概率不高但对应极端出力直接删掉会影响鲁棒性可以考虑把低概率极端场景合并进相邻高概率场景而不是简单删掉。6. 这套方案还能往哪些方向扩展代码框架搭好之后后续扩展空间其实很大。我目前正在把电转气P2G设备和碳排放约束加进模型因为双碳目标下综合能源系统协同优化不仅要算经济账还要算碳账。添加P2G设备时只需要在能量耦合矩阵里增加一条电→天然气的转换路径在设备约束里加P2G的电功率输入范围和转换效率目标函数里加对应的运维成本项整个框架不需要大改。另外如果研究对象从园区级扩展到区域级需要考虑多能源网络潮流约束那么场景生成模块可以和拉丁超立方抽样或准蒙特卡洛方法结合提高采样效率求解器方面可以从Cplex切换到Gurobi它在大规模MILP问题上表现更好。最后再分享一个实际操作的体会这套Matlab代码在Windows和Linux环境下我都跑过最容易出问题的环节不是模型本身而是环境配置——Yalmip把Cplex、Gurobi这些求解器路径配好后建议先跑一个官方demo确认求解器可用再跑综合能源系统模型。否则一旦报错你很难分清是模型写错还是求解器接口没接好。别问我怎么知道的我在这上面浪费过整整一天时间。