最近有个师弟拿仿真结果来找我说他在做计及新能源出力不确定性的电气设备综合能源系统协同优化时把储能和电转气都加进模型后系统总成本反而比不加还高怀疑自己代码写错了。我让他把风电、光伏的预测曲线翻出来看了一眼问题马上就定位到了——他把预测值当成真实值用做的是确定性优化目标函数里又没加弃风弃光惩罚项。这正是综合能源系统协同优化里最常见的认知盲区新能源出力不确定性不处理多能耦合再复杂也是白搭结果不仅失真甚至会给出反直觉的错误结论。这篇文章就是围绕这个方向做的完整复盘从不确定性建模、多能系统设备建模、优化模型构建到Matlab代码实现和结果分析全部串起来讲一遍。内容适合正在做综合能源系统、微电网、电-气耦合方向课题的研究生也适合刚接触Yalmip建模、想跑通一个完整算例的工程师。我会把每一步的为什么也讲清楚——为什么这样建模、为什么选这个不确定性处理方式、为什么这里必须加某个约束尽量让你拿到之后不仅会跑还能改、能调、能讲清楚。1. 难点拆解新能源出力不确定性与多能耦合为什么是组合难题1.1 风光预测误差确定性优化在这里为什么必然翻车先说一个基本事实风电和光伏出力本质上是一个随机过程。风速受天气系统影响光照受云层和气溶胶影响这些物理过程决定了预测模型无论怎么改进都不可能给出精确的出力曲线只能给出预测值加误差分布这种形式。实际工程和学术研究中常用的做法是把光伏出力用Beta分布来描述风电出力用Weibull分布或者正态误差来描述然后在这个基础上生成多个可能出现的出力场景。为什么不能直接用预测值做优化我拿一个最简单的电力平衡约束举例。假设某个时段系统里风机预测出力是100MW但实际可能只有70MW如果调度计划严格按照100MW来安排那这一时段就有30MW的缺口需要从电网购电或者切负荷。更麻烦的是实时市场中购电价格往往比日前调度高切负荷更是有惩罚成本。也就是说预测误差会在真实运行时转化为经济损失确定性优化完全看不到这一层。所以计及不确定性的本质不是让优化算法更聪明而是让优化模型具备预见偏差的能力——在做调度决策时就把可能出现的出力偏差考虑进去让决策方案在大部分场景下都能安全、经济运行。这就是所谓鲁棒性和经济性的权衡也是这个题目最核心的科学问题之一。1.2 电气-气-热多网络耦合协同优化的复杂度来源综合能源系统的综合两个字不是白叫的。我做的系统里同时包含电力网络、天然气网络和热力网络三者通过设备耦合在一起CHP机组消耗天然气发电、同时回收余热供热P2G装置用电制氢或制天然气燃气锅炉直接气转热电锅炉和热泵则是电转热。每台设备都是两个网络的连接点一个网络里的决策会直接影响另一个网络的供需平衡。这就带来两个层面的复杂度。第一是空间耦合电力平衡约束、天然气平衡约束和热力平衡约束必须同时满足单独优化某一个网络往往会在其他网络上造成不可行。第二是时间耦合储能电池的荷电状态、储热罐的储热量、机组爬坡约束都跨越多个时段上一时段的决策直接限制下一时段的可调空间。尤其是储能设备它本质上是在时间轴上搬运能量如果建模时忽略SOC的递推关系优化结果必然失真。再加上前面说的新能源不确定性整个问题就变成一个高维度、强耦合的随机优化问题。坦白讲这个问题没有解析解工程上可行的路线只有两条一条是场景法把随机问题转化为多个确定性场景下的期望值优化另一条是鲁棒优化直接在不确定集合内寻找最坏情况下的最优解。这两种路线我在第3节会展开对比这里先让大家对问题的复杂度有个整体认知。2. 系统建模先搞清设备模型再谈优化算法2.1 典型园区综合能源系统的设备拓扑做这个方向的第一步不是写代码而是把系统拓扑画清楚。我这里用的拓扑是一个典型的园区级电气-气-热综合能源系统包含的设备和能量流关系如下表设备输入能量输出能量在网络耦合中的角色风电机组风能电能电源出力不确定光伏阵列太阳能电能电源出力不确定CHP机组天然气电能热能电气热三网耦合核心燃气锅炉天然气热能气-热耦合P2G装置电能天然气电-气耦合电锅炉电能热能电-热耦合储能电池电能电能时间尺度上的电力平移储热罐热能热能时间尺度上的热力平移为什么把CHP和P2G同时放进来这就要说到协同优化的核心意义了。CHP机组把天然气变成电和热但它有一个先天性问题——热电产出比例受限。当系统里风电大发、电负荷又低的时候如果CHP为了供热必须开机多余的电就只能弃掉。这时候P2G就有了价值它能把多余的电变成天然气存起来既能消纳新能源又能补充气网相当于给系统加了一个电力-天然气双向流动的阀门。很多论文里管这种运行模式叫热电联产与电转气的协调运行本质上就是让电、气、热三个网络在时间上互补。2.2 关键设备的数学模型与参数选取设备模型是优化问题的骨架模型写得准不准直接决定优化结果的工程有效性。我逐个说一下我实际用的模型以及参数选择时的注意事项。光伏出力模型。工程上最常用的是标准工况折算模型P_pv(t) P_STC × (G(t)/G_STC) × [1 - β(T_c(t) - T_STC)]其中G为光照强度T_c为电池板温度。在优化调度里更常见的简化做法是直接给出预测出力曲线和误差分布把光伏当成一个负的、带不确定性的负荷。两种模型我都试过前者适合做容量规划后者适合做运行调度。运行调度里我推荐直接对预测出力乘以一个服从Beta分布的随机系数来生成场景这样处理起来最简单物理含义也清晰。风电出力模型。风电出力和风速之间是三次方关系P_w 0风速小于切入风速或大于切出风速、P_w P_r × (v - v_in)/(v_r - v_in)切入和额定之间、P_w P_r额定到切出之间。在不确定性建模里我一般直接用预测出力加上正态分布误差来生成场景误差标准差取预测出力的10%-20%这个区间是行业里比较常见的设定。CHP机组模型。CHP有两种建模方式定热电比和可变热电比。定热电比模型简单认为热出力P_h c_m × P_e其中c_m是常数适合粗略估算可变热电比模型则考虑了CHP的可运行域需要引入二进制变量来表示开停机状态热出力在一个区间内可以灵活调节。我的建议是如果你的研究重点是协同优化策略一定要用可变热电比模型因为热电比可调是CHP灵活性的核心来源定热电比会严重低估系统的调度空间。P2G装置模型。P2G的输入是电输出是天然气模型相对简单G_p2g(t) η_p2g × P_p2g(t)效率η_p2g一般在55%-70%之间。但P2G建模有一个容易被忽略的点它消耗的电能是可弃属性的——也就是说只有在系统存在新能源富余电量的时段P2G运行才有经济价值。所以目标函数里要配套设置弃风弃光惩罚项否则优化器不会主动让P2G消纳新能源而是直接在低谷电价时段买电制气这不符合物理直觉。储能电池模型。储能电池的SOC递推方程是SOC(t1) SOC(t) η_ch × P_ch(t) - (1/η_dis) × P_dis(t)同时要满足SOC上下限约束、充放电功率约束以及一个常被忽略的约束——调度周期始末SOC要相等否则优化结果等于在透支电池能量没有实际可操作性。2.3 目标函数与约束体系的搭建逻辑目标函数我选的是系统运行总成本最小化主要包括四部分向电网购电成本、购买天然气成本、设备运行维护成本、弃风弃光惩罚成本。这里有一个关键设计为什么必须加弃风弃光惩罚因为如果不加惩罚优化器在要不要消纳新能源这件事上是无差别的它只会盯着购电成本和购气成本的最小化可能出现宁可弃掉风电也不调整CHP出力的错误决策。加入惩罚项之后目标函数才真正体现优先消纳新能源的调度逻辑。约束体系方面核心约束有五大类电力平衡约束电网购电 风电 光伏 CHP发电 储能放电 电负荷 P2G耗电 电锅炉耗电 储能充电。天然气平衡约束气网购气 P2G产气 CHP耗气 燃气锅炉耗气。热力平衡约束CHP余热 燃气锅炉产热 电锅炉产热 储热罐放热 热负荷 储热罐储热。设备出力上下限约束每台设备的出力不能超过其额定容量。爬坡约束与储能SOC递推约束用于处理时间耦合。这套约束体系写完之后问题规模大概是多少以一个24时段的调度周期为例如果场景数是20个那么决策变量的数量级在几千到上万之间属于典型的混合整数线性规划MILP问题用Yalmip工具箱加CPLEX或Gurobi求解器可以高效求解。这也是为什么这个方向几乎都用Matlab代码实现的原因之一——Yalmip在学术界的生态太成熟了模型写起来直观调试也方便。3. 不确定性建模场景法与鲁棒优化的选型对比3.1 场景法蒙特卡洛采样与同步回代削减场景法是我在这个项目里的主力方法思路非常直白既然风光出力是随机的那我们就生成多组可能的出力曲线称为场景然后用这些场景下的加权平均成本作为优化目标得到一个在各种可能出现的情况下平均表现最好的调度方案。具体操作分三步。第一步确定误差分布。光伏场景我用的做法是对每个时段的预测值P_pv_pre(t)生成随机系数k(t) ~ Beta(α, β)使该时段场景出力为k(t) × P_pv_pre(t)。为什么用Beta分布因为Beta分布定义在[0,1]区间天然适合描述实际出力占预测出力比例这种变量。风电场景则用正态误差P_w_scene(t) P_w_pre(t) ε(t)其中ε(t) ~ N(0, σ_t²)σ_t取预测值的15%。第二步蒙特卡洛采样。对每个时序节点采样就得到一条完整的出力曲线也就是一个场景。我一般先采800到1000个初始场景这个数量足够覆盖分布的主要特征。第三步场景削减。1000个场景直接放进优化模型里计算量会很大而且很多场景之间差异很小对优化结果几乎没有独立贡献。所以要用削减算法把它们压缩到20个左右同时保留原始概率分布的主要信息。这里最常用的是同步回代削减法核心逻辑是反复寻找距离最近的两个场景删掉其中一个把它的概率加到另一个上面直到场景数量达到目标值。距离通常用欧氏距离也就是两个场景出力曲线在每个时段差值的平方和再开方。同步回代的关键实现细节有三个一是场景概率的归一化删完最后一批场景后所有场景概率之和必须重新归一化为1否则目标函数的期望成本算出来会偏小二是距离矩阵的计算要注意内存和速度1000个场景的距离矩阵是1000×1000每次迭代都要更新用Matlab向量化写法比for循环快一个数量级三是削减后的场景集合应该尽量覆盖分布的两端不能只留下接近均值的场景否则会低估不确定性带来的风险这一点可以通过对比削减前后场景总体的方差来验证。3.2 鲁棒优化盒式不确定集与不确定性预算场景法之外另一种主流路线是鲁棒优化。鲁棒优化的核心思想是我不去枚举具体场景而是把风光出力限定在一个不确定集合里然后在最坏情况下做优化决策——也就是说求得的方案保证在不确定集合内任何出力组合下都是可行的。最基础的是盒式不确定集每个时段的风光出力可以在预测值加减一个偏差范围内波动比如P_w(t) ∈ [P_w_pre(t) - Δ_w(t), P_w_pre(t) Δ_w(t)]。这种情况下问题会变成min-max形式求解难度比场景法高不少。为了控制保守程度工程上会引入不确定性预算Γ用来限制各个时段同时达到最坏情况的个数。比如Γ8表示在24个时段里最多允许8个时段的风电同时取到区间边界的偏差其余时段取预测值。Γ设得越大方案越保守运行成本越高Γ0退化为确定性优化。我再解释得直白一点。盒式不确定集本身太悲观了——如果允许每个时段风电都取最小值优化器为了避免缺电风险就会把燃气机组出力调得很高成本飙升。但实际情况是不可能每个时段都同时出现最坏偏差。不确定性预算Γ这个参数就是在告诉优化器别把最坏情况想得那么极端同时掉链子的时段是有限制的。这样一来鲁棒方案的成本会比纯盒式低同时应对不确定性能力依然比确定性方案强。3.3 怎么选才合适精度、保守性与求解规模的权衡很多做这个方向的同学上来就问场景法和鲁棒优化哪个好说实话没有绝对答案要看你的研究目的和算力条件。我把两种方法的对比整理成一张表对比维度场景法鲁棒优化建模难度较低关键是误差分布设对较高需要设计不确定集求解规模场景数×问题规模呈线性增长问题规模本身不大但约束形式复杂保守程度适中偏向期望值最优偏高偏向最坏情况最优结果解释性强每个场景物理含义清晰较弱最优方案对应的人为构造场景适用场景论文研究、经济调度安全要求高的系统、保守调峰决策如果只是跑通代码、验证协同优化的效果我建议你优先用场景法尤其是蒙特卡洛加同步回代这个组合。原因很实在代码逻辑直观结果可解释性强调试时哪里出问题一眼就能看到。鲁棒优化适合用在后续深入研究里比如探讨系统面对极端天气的抗风险能力时它的优势才会完全显现。我个人的做法是先把场景法跑通作为基线再用鲁棒优化做对比一篇文章里两种方法都用上从稳健性和经济性两个维度把问题聊透。4. Matlab代码实现从场景生成到Yalmip求解4.1 代码整体结构与数据准备整个Matlab工程我建议按模块拆成五个脚本文件而不是写一个大而全的main.m。这样做的理由是这个问题的代码量通常在几百行上下混在一起调试时会让人崩溃拆开后每个模块可以独立验证。我的工程结构是这样的data_definition.m—— 定义系统参数、负荷曲线、预测出力曲线、设备容量、能源价格。scenario_generation.m—— 生成风光出力场景并执行同步回代削减。model_build.m—— 定义决策变量、目标函数和约束条件构建Yalmip模型。solve_and_postprocess.m—— 调用求解器求解提取结果并计算各项成本指标。plot_results.m—— 绘制电功率平衡图、热功率平衡图、SOC曲线、各机组出力曲线。数据准备阶段最容易出错的地方是量纲统一。举个例子天然气热值单位通常用kW·h/m³但购气价格可能按元/m³给CHP的耗气量是立方米发电量和发热量是kW·h。千万要先统一成同一个能量单位体系我习惯全部换算成kW·h天然气价格也换算成元/kW·h后再写进目标函数。这一步错了后面的所有结果都是错的而且数值上还不会差得离谱特别隐蔽。4.2 场景生成与削减核心代码场景生成部分我直接给出可以用的代码框架。光伏场景生成和风电场景生成的逻辑稍有不同光伏用Beta系数缩放预测曲线风电用正态误差叠加。%% 光伏场景生成对预测出力应用Beta分布随机系数 % P_pv_pre: 1x24 光伏预测出力向量 % n_scenes: 需要生成的场景数量 n_scenes 1000; pv_scenes zeros(n_scenes, 24); for s 1:n_scenes % Beta分布的参数决定出力波动幅度这里选alpha9, beta7 % 对应的分布峰值为0.6附近表示多数情况下实际出力约为预测值的60%左右 k betarnd(9, 7, 1, 24); pv_scenes(s, :) k .* P_pv_pre; end %% 风电场景生成预测值叠加正态误差 % P_w_pre: 1x24 风电预测出力向量 % sigma_percent: 误差标准差占预测值的比例一般取0.15 sigma_percent 0.15; wind_scenes zeros(n_scenes, 24); for s 1:n_scenes epsilon sigma_percent .* P_w_pre .* randn(1, 24); wind_scenes(s, :) P_w_pre epsilon; % 去掉负值或超出额定容量的无效场景点 wind_scenes(s, wind_scenes(s, :) 0) 0; wind_scenes(s, wind_scenes(s, :) P_w_rated) P_w_rated; end场景削减最关键的就是同步回代法。完整代码如下这段代码我直接用了向量化写法1000个场景削减到20个用时在几秒量级不需要额外优化。function [reduced_scenes, reduced_prob] backward_reduction(scenes, prob, target_num) % scenes: n x 24 的原始场景矩阵 % prob: n x 1 的原始场景概率向量 % target_num: 目标场景数量 n size(scenes, 1); if n target_num reduced_scenes scenes; reduced_prob prob; return; end % 计算所有场景两两之间的距离 % 这里用向量化计算距离定义为24个时段出力差的平方和 dist_matrix zeros(n, n); for i 1:n diff scenes(i, :) - scenes; % 广播到所有场景 dist_matrix(i, :) sum(diff.^2, 2); end remaining true(n, 1); current_num n; while current_num target_num % 找出当前所有剩余场景中距离最近的一对 min_dist Inf; del_idx 0; keep_idx 0; for i 1:n if ~remaining(i) continue; end for j i1:n if ~remaining(j) continue; end if dist_matrix(i, j) min_dist min_dist dist_matrix(i, j); del_idx j; keep_idx i; end end end % 删除del_idx将概率加到keep_idx上 prob(keep_idx) prob(keep_idx) prob(del_idx); remaining(del_idx) false; current_num current_num - 1; % 更新与del_idx相关的距离为Inf避免后续再选到它 dist_matrix(:, del_idx) Inf; dist_matrix(del_idx, :) Inf; end reduced_scenes scenes(remaining, :); reduced_prob prob(remaining); % 概率归一化 reduced_prob reduced_prob / sum(reduced_prob); end这里有一个我在实际调试中深刻体会过的坑删除场景后如果不对所有距离做失效标记下次循环还可能选中这个已经被删除的场景。上面的代码用把距离矩阵对应行列置为Inf的方式来处理简单有效。另外场景削减的目标数量没有绝对标准我试过10个、20个、50个几组20个在这个系统规模下是性价比比较高的选择计算时间和结果精度平衡得比较好。4.3 Yalmip建模核心代码逐段解读模型构建我用的是Yalmip工具箱决策变量用sdpvar定义连续变量、binvar定义二进制变量。求解24时段调度问题时所有变量都定义成1×24的行向量便于后续直接拼接约束矩阵。%% 决策变量定义 % 电力相关 P_buy sdpvar(1, 24); % 从电网购电功率 P_chp sdpvar(1, 24); % CHP发电功率 P_pv sdpvar(1, 24); % 光伏实际使用出力 P_wind sdpvar(1, 24); % 风电实际使用出力 P_p2g sdpvar(1, 24); % P2G耗电 P_eb sdpvar(1, 24); % 电锅炉耗电 P_ess_ch sdpvar(1, 24); % 储能充电功率 P_ess_dis sdpvar(1, 24); % 储能放电功率 u_ess binvar(1, 24); % 储能充放电状态1表示充电 % 天然气与热力相关 G_buy sdpvar(1, 24); % 购买天然气 G_p2g sdpvar(1, 24); % P2G产气量 H_gb sdpvar(1, 24); % 燃气锅炉产热 H_chp sdpvar(1, 24); % CHP余热回收 H_eb sdpvar(1, 24); % 电锅炉产热 H_hs_ch sdpvar(1, 24); % 储热罐充热 H_hs_dis sdpvar(1, 24); % 储热罐放热 SOC_ess sdpvar(1, 24); % 储能电池荷电状态 SOC_hs sdpvar(1, 24); % 储热罐储热量 %% 目标函数运行总成本最小化 % 购电成本 购气成本 运维成本 弃风弃光惩罚 % 电价和天然气价格按实际时段数据给定 objective sum(P_buy .* Price_electric) ... sum(G_buy .* Price_gas) ... sum(P_chp .* c_om_chp) sum(H_gb .* c_om_gb) ... sum(P_p2g .* c_om_p2g) sum(P_eb .* c_om_eb) ... sum(P_ess_ch .* c_om_ess) sum(P_ess_dis .* c_om_ess) ... sum(P_pv_pre - P_pv) * punish_pv ... % 弃光惩罚 sum(P_wind_pre - P_wind) * punish_wind; % 弃风惩罚这里要特别说明一个建模技巧sum(P_pv_pre - P_pv) * punish_pv这个惩罚项要求P_pv的取值尽量接近预测值P_pv_pre从而实现优先消纳新能源的调度逻辑。注意P_pv被定义为sdpvar连续变量它不等于预测值而是由优化器在约束条件下决策的变量。这个设计是整个目标函数最关键的地方——它把新能源消纳变成了一个可量化、可优化的经济指标。约束条件部分我挑几个典型的列出来看。电力平衡约束是第一个要写的Constraints []; % 电力平衡购电风光CHP储能放电 负荷P2G电锅炉储能充电 for t 1:24 Constraints [Constraints, P_buy(t) P_pv(t) P_wind(t) P_chp(t) P_ess_dis(t) ... P_load(t) P_p2g(t) P_eb(t) P_ess_ch(t)]; end储能SOC递推约束和充放电互斥约束是这个模型里最容易写错的地方。充放电互斥约束用的是二进制变量加一个大M不等式的方式——当u_ess(t)1时充电功率可以大于0放电功率必须为0反之亦然。我之前见过很多初学代码的同学不写互斥约束结果求解出来的方案是储能一边充电一边放电白白损耗能量这个约束不能省。% SOC递推约束第一个时段用初始SOC之后的时段用上一时段SOC SOC_ess(1) SOC_ess_init; for t 2:24 Constraints [Constraints, SOC_ess(t) SOC_ess(t-1) ... eta_ess_ch * P_ess_ch(t) - (1/eta_ess_dis) * P_ess_dis(t)]; end % 调度周期末SOC回到初始值防止透支电池 Constraints [Constraints, SOC_ess(24) SOC_ess_init]; % 充放电互斥约束用二进制变量限制 Constraints [Constraints, P_ess_ch M * u_ess, P_ess_dis M * (1 - u_ess)];天然气平衡约束是综合能源系统的特色环节% 天然气平衡购气P2G产气 CHP耗气燃气锅炉耗气 % CHP耗气量 P_chp / eta_chp_e / Gas_heat_value % 燃气锅炉耗气量 H_gb / eta_gb / Gas_heat_value for t 1:24 Constraints [Constraints, G_buy(t) G_p2g(t) ... P_chp(t) / eta_chp_e / Gas_heat_value H_gb(t) / eta_gb / Gas_heat_value]; end什么时候用for循环写约束、什么时候可以用向量化一次写完我的经验是涉及多时段递推关系的约束必须用for循环因为每个时段都依赖前一个时段但上下限约束、平衡约束这类时段独立的约束可以直接向量化写比如Constraints [Constraints, 0 P_chp P_chp_max]Yalmip会自动扩展成24个约束。这个习惯能显著减少代码量但前提是你对Yalmip的向量约束展开机制足够熟悉否则建议全部用for循环稳一点。4.4 求解器配置与结果提取模型构建完成后求解部分很简单但有几个配置细节会明显影响求解效率。我用的求解器是CPLEX配置代码如下options sdpsettings(solver, cplex, verbose, 2); % 对于MILP问题设置允许的MIP gap为0.1%不需要调到0因为这样会显著拖慢求解 options.cplex.mip.tolerances.mipgap 0.001; sol optimize(Constraints, objective, options); if sol.problem ~ 0 disp(求解失败); disp(sol.info); else disp(求解成功); end求解器选择上有个小知识点如果模型里不含二进制变量问题退化为LP或QP用内点法求解器如SDPT3、SeDuMi也行但如果模型里含二进制变量比如CHP开停机、储能充放电状态就一定要用MILP求解器也就是CPLEX或Gurobi。我实际测试过同样的模型在CPLEX下求解时间在几秒到几十秒之间换成内点法求解器可能几分钟都跑不完而且大概率会报错。跑这个方向CPLEX或Gurobi是必需品直接在Yalmip里调用即可。结果提取也有一个容易踩的坑用value()函数提取解之后变量是double型向量可以直接用来计算各项指标。但要注意如果模型不可行或求解器在迭代中途退出value()返回的结果是不可信的一定要先判断sol.problem是否为0。另外场景法求出来的期望运行成本应该等于所有场景下成本的平均值严格说是概率加权平均这个值才是论文里要报的期望总成本不是某个单一场景下的成本数值。5. 结果分析与复现避坑指南5.1 协同优化效果怎么评估模型求解完成之后不能只盯着总成本这一个数。我做结果分析时习惯从四个维度看第一总成本和成本构成。把购电成本、购气成本、运维成本和弃风弃光惩罚分开汇总看成本大头在哪里。一般来说购电和购气成本占比在80%以上如果弃风弃光惩罚占比明显偏高说明不确定性处理能力不足需要调整场景数量或惩罚系数。第二新能源消纳率。消纳率等于实际使用的风光电量除以理论可发电量。这个指标是判断协同优化效果的核心——如果系统的储能、P2G、电锅炉等灵活性设备都发挥了作用消纳率应该明显高于不协同的情况。我建议把消纳率画成时序曲线配合弃风弃光功率来看能直观看出哪个时段弃风最严重、对应设备当时在干什么。第三设备利用情况。看CHP的发电出力和余热回收出力曲线看P2G的产气量分布在哪些时段看储能SOC曲线的变化形态。SOC曲线尤其值得研究如果SOC曲线一直在上限附近来回撞说明储能容量相对于系统波动来说太大了有点浪费如果一天下来SOC变化幅度很小说明储能没出力或者模型约束写得太紧。第四与确定性优化的对比。这是验证计及不确定性价值的关键实验。跑两个模型一个用预测值直接做确定性优化另一个用场景法做不确定性优化然后对比两个方案在各个场景下的真实成本。确定性方案在预测场景下成本可能很低但在其他场景下可能会有较高的缺电惩罚成本期望成本反而更高。这个对比是论文里很有说服力的图表之一。5.2 复现这个项目时我踩过的五个典型坑第一个坑是量纲不统一。这个在第4.1已经讲过但值得再强调一次。天然气热值、效率百分比、价格单位、功率单位任何一个不一致都会让结果数值对不上甚至让约束关系错误。我的建议是建一个单位换算参数区把所有物理量先换算成统一的kW·h再写进约束和目标函数。第二个坑是SOC初值设置不当。如果初始SOC设得太高调度周期内储能的放电空间很小灵活性被浪费如果设得太低前几个时段储能可能强制充电导致系统运行成本升高。我测试下来初始SOC设在0.2到0.3之间配合周期末回到初始值的约束储能利用效果最好。另外末端SOC等于初始SOC这个约束绝对不能去掉否则优化结果在工程上是不可执行的。第三个坑是惩罚系数太小。弃风弃光惩罚系数如果设置过低优化器宁可弃掉新能源也不愿意让燃气机组调低出力因为调低出力会损失供热。但系数太高也不好会让优化器不计代价地消纳新能源导致P2G和储能在低效区间运行整体经济性反而下降。我常用的标定方法是先用确定性优化跑一遍看看系统在完全没有不确定性约束时会有多少弃风弃光然后按弃风弃光电量的成本损失来反推惩罚系数的量级。第四个坑是CHP的热电比选择。如果直接用定热电比模型CHP的热电产出被锁死储能、P2G这些设备就失去了调节空间协同优化的价值大打折扣。我在第2.2里说过一定要用可变热电比模型让CHP在可运行域内灵活调节热电比这样才能体现出综合能源系统多能互补的调度灵活性。第五个坑是场景削减目标数量不校验。很多人削减完场景就直接丢进模型也不检查削减后的场景集是否保留了原始分布的关键特征。我建议削减前后做一次对比算一下原始场景和削减后场景的风光出力均值与方差如果方差明显变小说明削减把很多波动性场景削掉了模型会低估不确定性得到的结果偏乐观。这时候要增加目标场景数量或调整削减算法。5.3 参数敏感性分析的正确打开方式最后聊一下参数敏感性分析。很多人在论文里做敏感性分析就是把某个参数改几个值看看总成本怎么变然后画一条曲线了事。这种做法的信息量很低审稿人也一眼能看出来是凑数。我自己的经验是抓住三个最有工程意义的参数做分析不确定性预算Γ、储能容量、P2G容量。不确定性预算Γ的敏感性分析最有意思。把Γ从0确定性优化逐渐增加到24最保守工况画出运行成本和新能源消纳率的变化曲线。Γ小的时候系统运行成本低但消纳率也低Γ增大会让成本上升但消纳率也会随之改善。曲线中会出现一个拐点过了这个拐点后继续增大Γ成本上升很快但消纳率几乎不再提升说明不确定性预算已经超出了实际需要的保守程度这个拐点附近就是调度决策的最佳保守区间。储能容量和P2G容量的敏感性分析则要成对做。把储能容量从基准值的0.5倍增加到2倍同时记录系统总成本和消纳率的变化你会看到典型的边际递减规律容量增大初期效果明显但超过某个值之后收益趋缓经济上的最优容量就在这个拐点附近。P2G容量也是类似的逻辑。这两个参数的联合敏感性分析可以直接为系统规划和设备扩容提供参考比单参数扫描有意义得多。从我自己复现这个课题的经历来看最值得投入精力的不是把模型写得更复杂而是把不确定性的场景质量做好。场景质量直接决定了后续所有优化结果的可靠性一个分布建模失真的场景集无论求解器多强大、优化算法多先进最后产出的都只是看起来合理的数字。用场景法把基线跑通再往鲁棒优化方向深挖这条路走下来无论是做课题还是发论文都能拿得出手。 计及新能源出力不确定性的综合能源系统协同优化建模与求解
最近有个师弟拿仿真结果来找我说他在做计及新能源出力不确定性的电气设备综合能源系统协同优化时把储能和电转气都加进模型后系统总成本反而比不加还高怀疑自己代码写错了。我让他把风电、光伏的预测曲线翻出来看了一眼问题马上就定位到了——他把预测值当成真实值用做的是确定性优化目标函数里又没加弃风弃光惩罚项。这正是综合能源系统协同优化里最常见的认知盲区新能源出力不确定性不处理多能耦合再复杂也是白搭结果不仅失真甚至会给出反直觉的错误结论。这篇文章就是围绕这个方向做的完整复盘从不确定性建模、多能系统设备建模、优化模型构建到Matlab代码实现和结果分析全部串起来讲一遍。内容适合正在做综合能源系统、微电网、电-气耦合方向课题的研究生也适合刚接触Yalmip建模、想跑通一个完整算例的工程师。我会把每一步的为什么也讲清楚——为什么这样建模、为什么选这个不确定性处理方式、为什么这里必须加某个约束尽量让你拿到之后不仅会跑还能改、能调、能讲清楚。1. 难点拆解新能源出力不确定性与多能耦合为什么是组合难题1.1 风光预测误差确定性优化在这里为什么必然翻车先说一个基本事实风电和光伏出力本质上是一个随机过程。风速受天气系统影响光照受云层和气溶胶影响这些物理过程决定了预测模型无论怎么改进都不可能给出精确的出力曲线只能给出预测值加误差分布这种形式。实际工程和学术研究中常用的做法是把光伏出力用Beta分布来描述风电出力用Weibull分布或者正态误差来描述然后在这个基础上生成多个可能出现的出力场景。为什么不能直接用预测值做优化我拿一个最简单的电力平衡约束举例。假设某个时段系统里风机预测出力是100MW但实际可能只有70MW如果调度计划严格按照100MW来安排那这一时段就有30MW的缺口需要从电网购电或者切负荷。更麻烦的是实时市场中购电价格往往比日前调度高切负荷更是有惩罚成本。也就是说预测误差会在真实运行时转化为经济损失确定性优化完全看不到这一层。所以计及不确定性的本质不是让优化算法更聪明而是让优化模型具备预见偏差的能力——在做调度决策时就把可能出现的出力偏差考虑进去让决策方案在大部分场景下都能安全、经济运行。这就是所谓鲁棒性和经济性的权衡也是这个题目最核心的科学问题之一。1.2 电气-气-热多网络耦合协同优化的复杂度来源综合能源系统的综合两个字不是白叫的。我做的系统里同时包含电力网络、天然气网络和热力网络三者通过设备耦合在一起CHP机组消耗天然气发电、同时回收余热供热P2G装置用电制氢或制天然气燃气锅炉直接气转热电锅炉和热泵则是电转热。每台设备都是两个网络的连接点一个网络里的决策会直接影响另一个网络的供需平衡。这就带来两个层面的复杂度。第一是空间耦合电力平衡约束、天然气平衡约束和热力平衡约束必须同时满足单独优化某一个网络往往会在其他网络上造成不可行。第二是时间耦合储能电池的荷电状态、储热罐的储热量、机组爬坡约束都跨越多个时段上一时段的决策直接限制下一时段的可调空间。尤其是储能设备它本质上是在时间轴上搬运能量如果建模时忽略SOC的递推关系优化结果必然失真。再加上前面说的新能源不确定性整个问题就变成一个高维度、强耦合的随机优化问题。坦白讲这个问题没有解析解工程上可行的路线只有两条一条是场景法把随机问题转化为多个确定性场景下的期望值优化另一条是鲁棒优化直接在不确定集合内寻找最坏情况下的最优解。这两种路线我在第3节会展开对比这里先让大家对问题的复杂度有个整体认知。2. 系统建模先搞清设备模型再谈优化算法2.1 典型园区综合能源系统的设备拓扑做这个方向的第一步不是写代码而是把系统拓扑画清楚。我这里用的拓扑是一个典型的园区级电气-气-热综合能源系统包含的设备和能量流关系如下表设备输入能量输出能量在网络耦合中的角色风电机组风能电能电源出力不确定光伏阵列太阳能电能电源出力不确定CHP机组天然气电能热能电气热三网耦合核心燃气锅炉天然气热能气-热耦合P2G装置电能天然气电-气耦合电锅炉电能热能电-热耦合储能电池电能电能时间尺度上的电力平移储热罐热能热能时间尺度上的热力平移为什么把CHP和P2G同时放进来这就要说到协同优化的核心意义了。CHP机组把天然气变成电和热但它有一个先天性问题——热电产出比例受限。当系统里风电大发、电负荷又低的时候如果CHP为了供热必须开机多余的电就只能弃掉。这时候P2G就有了价值它能把多余的电变成天然气存起来既能消纳新能源又能补充气网相当于给系统加了一个电力-天然气双向流动的阀门。很多论文里管这种运行模式叫热电联产与电转气的协调运行本质上就是让电、气、热三个网络在时间上互补。2.2 关键设备的数学模型与参数选取设备模型是优化问题的骨架模型写得准不准直接决定优化结果的工程有效性。我逐个说一下我实际用的模型以及参数选择时的注意事项。光伏出力模型。工程上最常用的是标准工况折算模型P_pv(t) P_STC × (G(t)/G_STC) × [1 - β(T_c(t) - T_STC)]其中G为光照强度T_c为电池板温度。在优化调度里更常见的简化做法是直接给出预测出力曲线和误差分布把光伏当成一个负的、带不确定性的负荷。两种模型我都试过前者适合做容量规划后者适合做运行调度。运行调度里我推荐直接对预测出力乘以一个服从Beta分布的随机系数来生成场景这样处理起来最简单物理含义也清晰。风电出力模型。风电出力和风速之间是三次方关系P_w 0风速小于切入风速或大于切出风速、P_w P_r × (v - v_in)/(v_r - v_in)切入和额定之间、P_w P_r额定到切出之间。在不确定性建模里我一般直接用预测出力加上正态分布误差来生成场景误差标准差取预测出力的10%-20%这个区间是行业里比较常见的设定。CHP机组模型。CHP有两种建模方式定热电比和可变热电比。定热电比模型简单认为热出力P_h c_m × P_e其中c_m是常数适合粗略估算可变热电比模型则考虑了CHP的可运行域需要引入二进制变量来表示开停机状态热出力在一个区间内可以灵活调节。我的建议是如果你的研究重点是协同优化策略一定要用可变热电比模型因为热电比可调是CHP灵活性的核心来源定热电比会严重低估系统的调度空间。P2G装置模型。P2G的输入是电输出是天然气模型相对简单G_p2g(t) η_p2g × P_p2g(t)效率η_p2g一般在55%-70%之间。但P2G建模有一个容易被忽略的点它消耗的电能是可弃属性的——也就是说只有在系统存在新能源富余电量的时段P2G运行才有经济价值。所以目标函数里要配套设置弃风弃光惩罚项否则优化器不会主动让P2G消纳新能源而是直接在低谷电价时段买电制气这不符合物理直觉。储能电池模型。储能电池的SOC递推方程是SOC(t1) SOC(t) η_ch × P_ch(t) - (1/η_dis) × P_dis(t)同时要满足SOC上下限约束、充放电功率约束以及一个常被忽略的约束——调度周期始末SOC要相等否则优化结果等于在透支电池能量没有实际可操作性。2.3 目标函数与约束体系的搭建逻辑目标函数我选的是系统运行总成本最小化主要包括四部分向电网购电成本、购买天然气成本、设备运行维护成本、弃风弃光惩罚成本。这里有一个关键设计为什么必须加弃风弃光惩罚因为如果不加惩罚优化器在要不要消纳新能源这件事上是无差别的它只会盯着购电成本和购气成本的最小化可能出现宁可弃掉风电也不调整CHP出力的错误决策。加入惩罚项之后目标函数才真正体现优先消纳新能源的调度逻辑。约束体系方面核心约束有五大类电力平衡约束电网购电 风电 光伏 CHP发电 储能放电 电负荷 P2G耗电 电锅炉耗电 储能充电。天然气平衡约束气网购气 P2G产气 CHP耗气 燃气锅炉耗气。热力平衡约束CHP余热 燃气锅炉产热 电锅炉产热 储热罐放热 热负荷 储热罐储热。设备出力上下限约束每台设备的出力不能超过其额定容量。爬坡约束与储能SOC递推约束用于处理时间耦合。这套约束体系写完之后问题规模大概是多少以一个24时段的调度周期为例如果场景数是20个那么决策变量的数量级在几千到上万之间属于典型的混合整数线性规划MILP问题用Yalmip工具箱加CPLEX或Gurobi求解器可以高效求解。这也是为什么这个方向几乎都用Matlab代码实现的原因之一——Yalmip在学术界的生态太成熟了模型写起来直观调试也方便。3. 不确定性建模场景法与鲁棒优化的选型对比3.1 场景法蒙特卡洛采样与同步回代削减场景法是我在这个项目里的主力方法思路非常直白既然风光出力是随机的那我们就生成多组可能的出力曲线称为场景然后用这些场景下的加权平均成本作为优化目标得到一个在各种可能出现的情况下平均表现最好的调度方案。具体操作分三步。第一步确定误差分布。光伏场景我用的做法是对每个时段的预测值P_pv_pre(t)生成随机系数k(t) ~ Beta(α, β)使该时段场景出力为k(t) × P_pv_pre(t)。为什么用Beta分布因为Beta分布定义在[0,1]区间天然适合描述实际出力占预测出力比例这种变量。风电场景则用正态误差P_w_scene(t) P_w_pre(t) ε(t)其中ε(t) ~ N(0, σ_t²)σ_t取预测值的15%。第二步蒙特卡洛采样。对每个时序节点采样就得到一条完整的出力曲线也就是一个场景。我一般先采800到1000个初始场景这个数量足够覆盖分布的主要特征。第三步场景削减。1000个场景直接放进优化模型里计算量会很大而且很多场景之间差异很小对优化结果几乎没有独立贡献。所以要用削减算法把它们压缩到20个左右同时保留原始概率分布的主要信息。这里最常用的是同步回代削减法核心逻辑是反复寻找距离最近的两个场景删掉其中一个把它的概率加到另一个上面直到场景数量达到目标值。距离通常用欧氏距离也就是两个场景出力曲线在每个时段差值的平方和再开方。同步回代的关键实现细节有三个一是场景概率的归一化删完最后一批场景后所有场景概率之和必须重新归一化为1否则目标函数的期望成本算出来会偏小二是距离矩阵的计算要注意内存和速度1000个场景的距离矩阵是1000×1000每次迭代都要更新用Matlab向量化写法比for循环快一个数量级三是削减后的场景集合应该尽量覆盖分布的两端不能只留下接近均值的场景否则会低估不确定性带来的风险这一点可以通过对比削减前后场景总体的方差来验证。3.2 鲁棒优化盒式不确定集与不确定性预算场景法之外另一种主流路线是鲁棒优化。鲁棒优化的核心思想是我不去枚举具体场景而是把风光出力限定在一个不确定集合里然后在最坏情况下做优化决策——也就是说求得的方案保证在不确定集合内任何出力组合下都是可行的。最基础的是盒式不确定集每个时段的风光出力可以在预测值加减一个偏差范围内波动比如P_w(t) ∈ [P_w_pre(t) - Δ_w(t), P_w_pre(t) Δ_w(t)]。这种情况下问题会变成min-max形式求解难度比场景法高不少。为了控制保守程度工程上会引入不确定性预算Γ用来限制各个时段同时达到最坏情况的个数。比如Γ8表示在24个时段里最多允许8个时段的风电同时取到区间边界的偏差其余时段取预测值。Γ设得越大方案越保守运行成本越高Γ0退化为确定性优化。我再解释得直白一点。盒式不确定集本身太悲观了——如果允许每个时段风电都取最小值优化器为了避免缺电风险就会把燃气机组出力调得很高成本飙升。但实际情况是不可能每个时段都同时出现最坏偏差。不确定性预算Γ这个参数就是在告诉优化器别把最坏情况想得那么极端同时掉链子的时段是有限制的。这样一来鲁棒方案的成本会比纯盒式低同时应对不确定性能力依然比确定性方案强。3.3 怎么选才合适精度、保守性与求解规模的权衡很多做这个方向的同学上来就问场景法和鲁棒优化哪个好说实话没有绝对答案要看你的研究目的和算力条件。我把两种方法的对比整理成一张表对比维度场景法鲁棒优化建模难度较低关键是误差分布设对较高需要设计不确定集求解规模场景数×问题规模呈线性增长问题规模本身不大但约束形式复杂保守程度适中偏向期望值最优偏高偏向最坏情况最优结果解释性强每个场景物理含义清晰较弱最优方案对应的人为构造场景适用场景论文研究、经济调度安全要求高的系统、保守调峰决策如果只是跑通代码、验证协同优化的效果我建议你优先用场景法尤其是蒙特卡洛加同步回代这个组合。原因很实在代码逻辑直观结果可解释性强调试时哪里出问题一眼就能看到。鲁棒优化适合用在后续深入研究里比如探讨系统面对极端天气的抗风险能力时它的优势才会完全显现。我个人的做法是先把场景法跑通作为基线再用鲁棒优化做对比一篇文章里两种方法都用上从稳健性和经济性两个维度把问题聊透。4. Matlab代码实现从场景生成到Yalmip求解4.1 代码整体结构与数据准备整个Matlab工程我建议按模块拆成五个脚本文件而不是写一个大而全的main.m。这样做的理由是这个问题的代码量通常在几百行上下混在一起调试时会让人崩溃拆开后每个模块可以独立验证。我的工程结构是这样的data_definition.m—— 定义系统参数、负荷曲线、预测出力曲线、设备容量、能源价格。scenario_generation.m—— 生成风光出力场景并执行同步回代削减。model_build.m—— 定义决策变量、目标函数和约束条件构建Yalmip模型。solve_and_postprocess.m—— 调用求解器求解提取结果并计算各项成本指标。plot_results.m—— 绘制电功率平衡图、热功率平衡图、SOC曲线、各机组出力曲线。数据准备阶段最容易出错的地方是量纲统一。举个例子天然气热值单位通常用kW·h/m³但购气价格可能按元/m³给CHP的耗气量是立方米发电量和发热量是kW·h。千万要先统一成同一个能量单位体系我习惯全部换算成kW·h天然气价格也换算成元/kW·h后再写进目标函数。这一步错了后面的所有结果都是错的而且数值上还不会差得离谱特别隐蔽。4.2 场景生成与削减核心代码场景生成部分我直接给出可以用的代码框架。光伏场景生成和风电场景生成的逻辑稍有不同光伏用Beta系数缩放预测曲线风电用正态误差叠加。%% 光伏场景生成对预测出力应用Beta分布随机系数 % P_pv_pre: 1x24 光伏预测出力向量 % n_scenes: 需要生成的场景数量 n_scenes 1000; pv_scenes zeros(n_scenes, 24); for s 1:n_scenes % Beta分布的参数决定出力波动幅度这里选alpha9, beta7 % 对应的分布峰值为0.6附近表示多数情况下实际出力约为预测值的60%左右 k betarnd(9, 7, 1, 24); pv_scenes(s, :) k .* P_pv_pre; end %% 风电场景生成预测值叠加正态误差 % P_w_pre: 1x24 风电预测出力向量 % sigma_percent: 误差标准差占预测值的比例一般取0.15 sigma_percent 0.15; wind_scenes zeros(n_scenes, 24); for s 1:n_scenes epsilon sigma_percent .* P_w_pre .* randn(1, 24); wind_scenes(s, :) P_w_pre epsilon; % 去掉负值或超出额定容量的无效场景点 wind_scenes(s, wind_scenes(s, :) 0) 0; wind_scenes(s, wind_scenes(s, :) P_w_rated) P_w_rated; end场景削减最关键的就是同步回代法。完整代码如下这段代码我直接用了向量化写法1000个场景削减到20个用时在几秒量级不需要额外优化。function [reduced_scenes, reduced_prob] backward_reduction(scenes, prob, target_num) % scenes: n x 24 的原始场景矩阵 % prob: n x 1 的原始场景概率向量 % target_num: 目标场景数量 n size(scenes, 1); if n target_num reduced_scenes scenes; reduced_prob prob; return; end % 计算所有场景两两之间的距离 % 这里用向量化计算距离定义为24个时段出力差的平方和 dist_matrix zeros(n, n); for i 1:n diff scenes(i, :) - scenes; % 广播到所有场景 dist_matrix(i, :) sum(diff.^2, 2); end remaining true(n, 1); current_num n; while current_num target_num % 找出当前所有剩余场景中距离最近的一对 min_dist Inf; del_idx 0; keep_idx 0; for i 1:n if ~remaining(i) continue; end for j i1:n if ~remaining(j) continue; end if dist_matrix(i, j) min_dist min_dist dist_matrix(i, j); del_idx j; keep_idx i; end end end % 删除del_idx将概率加到keep_idx上 prob(keep_idx) prob(keep_idx) prob(del_idx); remaining(del_idx) false; current_num current_num - 1; % 更新与del_idx相关的距离为Inf避免后续再选到它 dist_matrix(:, del_idx) Inf; dist_matrix(del_idx, :) Inf; end reduced_scenes scenes(remaining, :); reduced_prob prob(remaining); % 概率归一化 reduced_prob reduced_prob / sum(reduced_prob); end这里有一个我在实际调试中深刻体会过的坑删除场景后如果不对所有距离做失效标记下次循环还可能选中这个已经被删除的场景。上面的代码用把距离矩阵对应行列置为Inf的方式来处理简单有效。另外场景削减的目标数量没有绝对标准我试过10个、20个、50个几组20个在这个系统规模下是性价比比较高的选择计算时间和结果精度平衡得比较好。4.3 Yalmip建模核心代码逐段解读模型构建我用的是Yalmip工具箱决策变量用sdpvar定义连续变量、binvar定义二进制变量。求解24时段调度问题时所有变量都定义成1×24的行向量便于后续直接拼接约束矩阵。%% 决策变量定义 % 电力相关 P_buy sdpvar(1, 24); % 从电网购电功率 P_chp sdpvar(1, 24); % CHP发电功率 P_pv sdpvar(1, 24); % 光伏实际使用出力 P_wind sdpvar(1, 24); % 风电实际使用出力 P_p2g sdpvar(1, 24); % P2G耗电 P_eb sdpvar(1, 24); % 电锅炉耗电 P_ess_ch sdpvar(1, 24); % 储能充电功率 P_ess_dis sdpvar(1, 24); % 储能放电功率 u_ess binvar(1, 24); % 储能充放电状态1表示充电 % 天然气与热力相关 G_buy sdpvar(1, 24); % 购买天然气 G_p2g sdpvar(1, 24); % P2G产气量 H_gb sdpvar(1, 24); % 燃气锅炉产热 H_chp sdpvar(1, 24); % CHP余热回收 H_eb sdpvar(1, 24); % 电锅炉产热 H_hs_ch sdpvar(1, 24); % 储热罐充热 H_hs_dis sdpvar(1, 24); % 储热罐放热 SOC_ess sdpvar(1, 24); % 储能电池荷电状态 SOC_hs sdpvar(1, 24); % 储热罐储热量 %% 目标函数运行总成本最小化 % 购电成本 购气成本 运维成本 弃风弃光惩罚 % 电价和天然气价格按实际时段数据给定 objective sum(P_buy .* Price_electric) ... sum(G_buy .* Price_gas) ... sum(P_chp .* c_om_chp) sum(H_gb .* c_om_gb) ... sum(P_p2g .* c_om_p2g) sum(P_eb .* c_om_eb) ... sum(P_ess_ch .* c_om_ess) sum(P_ess_dis .* c_om_ess) ... sum(P_pv_pre - P_pv) * punish_pv ... % 弃光惩罚 sum(P_wind_pre - P_wind) * punish_wind; % 弃风惩罚这里要特别说明一个建模技巧sum(P_pv_pre - P_pv) * punish_pv这个惩罚项要求P_pv的取值尽量接近预测值P_pv_pre从而实现优先消纳新能源的调度逻辑。注意P_pv被定义为sdpvar连续变量它不等于预测值而是由优化器在约束条件下决策的变量。这个设计是整个目标函数最关键的地方——它把新能源消纳变成了一个可量化、可优化的经济指标。约束条件部分我挑几个典型的列出来看。电力平衡约束是第一个要写的Constraints []; % 电力平衡购电风光CHP储能放电 负荷P2G电锅炉储能充电 for t 1:24 Constraints [Constraints, P_buy(t) P_pv(t) P_wind(t) P_chp(t) P_ess_dis(t) ... P_load(t) P_p2g(t) P_eb(t) P_ess_ch(t)]; end储能SOC递推约束和充放电互斥约束是这个模型里最容易写错的地方。充放电互斥约束用的是二进制变量加一个大M不等式的方式——当u_ess(t)1时充电功率可以大于0放电功率必须为0反之亦然。我之前见过很多初学代码的同学不写互斥约束结果求解出来的方案是储能一边充电一边放电白白损耗能量这个约束不能省。% SOC递推约束第一个时段用初始SOC之后的时段用上一时段SOC SOC_ess(1) SOC_ess_init; for t 2:24 Constraints [Constraints, SOC_ess(t) SOC_ess(t-1) ... eta_ess_ch * P_ess_ch(t) - (1/eta_ess_dis) * P_ess_dis(t)]; end % 调度周期末SOC回到初始值防止透支电池 Constraints [Constraints, SOC_ess(24) SOC_ess_init]; % 充放电互斥约束用二进制变量限制 Constraints [Constraints, P_ess_ch M * u_ess, P_ess_dis M * (1 - u_ess)];天然气平衡约束是综合能源系统的特色环节% 天然气平衡购气P2G产气 CHP耗气燃气锅炉耗气 % CHP耗气量 P_chp / eta_chp_e / Gas_heat_value % 燃气锅炉耗气量 H_gb / eta_gb / Gas_heat_value for t 1:24 Constraints [Constraints, G_buy(t) G_p2g(t) ... P_chp(t) / eta_chp_e / Gas_heat_value H_gb(t) / eta_gb / Gas_heat_value]; end什么时候用for循环写约束、什么时候可以用向量化一次写完我的经验是涉及多时段递推关系的约束必须用for循环因为每个时段都依赖前一个时段但上下限约束、平衡约束这类时段独立的约束可以直接向量化写比如Constraints [Constraints, 0 P_chp P_chp_max]Yalmip会自动扩展成24个约束。这个习惯能显著减少代码量但前提是你对Yalmip的向量约束展开机制足够熟悉否则建议全部用for循环稳一点。4.4 求解器配置与结果提取模型构建完成后求解部分很简单但有几个配置细节会明显影响求解效率。我用的求解器是CPLEX配置代码如下options sdpsettings(solver, cplex, verbose, 2); % 对于MILP问题设置允许的MIP gap为0.1%不需要调到0因为这样会显著拖慢求解 options.cplex.mip.tolerances.mipgap 0.001; sol optimize(Constraints, objective, options); if sol.problem ~ 0 disp(求解失败); disp(sol.info); else disp(求解成功); end求解器选择上有个小知识点如果模型里不含二进制变量问题退化为LP或QP用内点法求解器如SDPT3、SeDuMi也行但如果模型里含二进制变量比如CHP开停机、储能充放电状态就一定要用MILP求解器也就是CPLEX或Gurobi。我实际测试过同样的模型在CPLEX下求解时间在几秒到几十秒之间换成内点法求解器可能几分钟都跑不完而且大概率会报错。跑这个方向CPLEX或Gurobi是必需品直接在Yalmip里调用即可。结果提取也有一个容易踩的坑用value()函数提取解之后变量是double型向量可以直接用来计算各项指标。但要注意如果模型不可行或求解器在迭代中途退出value()返回的结果是不可信的一定要先判断sol.problem是否为0。另外场景法求出来的期望运行成本应该等于所有场景下成本的平均值严格说是概率加权平均这个值才是论文里要报的期望总成本不是某个单一场景下的成本数值。5. 结果分析与复现避坑指南5.1 协同优化效果怎么评估模型求解完成之后不能只盯着总成本这一个数。我做结果分析时习惯从四个维度看第一总成本和成本构成。把购电成本、购气成本、运维成本和弃风弃光惩罚分开汇总看成本大头在哪里。一般来说购电和购气成本占比在80%以上如果弃风弃光惩罚占比明显偏高说明不确定性处理能力不足需要调整场景数量或惩罚系数。第二新能源消纳率。消纳率等于实际使用的风光电量除以理论可发电量。这个指标是判断协同优化效果的核心——如果系统的储能、P2G、电锅炉等灵活性设备都发挥了作用消纳率应该明显高于不协同的情况。我建议把消纳率画成时序曲线配合弃风弃光功率来看能直观看出哪个时段弃风最严重、对应设备当时在干什么。第三设备利用情况。看CHP的发电出力和余热回收出力曲线看P2G的产气量分布在哪些时段看储能SOC曲线的变化形态。SOC曲线尤其值得研究如果SOC曲线一直在上限附近来回撞说明储能容量相对于系统波动来说太大了有点浪费如果一天下来SOC变化幅度很小说明储能没出力或者模型约束写得太紧。第四与确定性优化的对比。这是验证计及不确定性价值的关键实验。跑两个模型一个用预测值直接做确定性优化另一个用场景法做不确定性优化然后对比两个方案在各个场景下的真实成本。确定性方案在预测场景下成本可能很低但在其他场景下可能会有较高的缺电惩罚成本期望成本反而更高。这个对比是论文里很有说服力的图表之一。5.2 复现这个项目时我踩过的五个典型坑第一个坑是量纲不统一。这个在第4.1已经讲过但值得再强调一次。天然气热值、效率百分比、价格单位、功率单位任何一个不一致都会让结果数值对不上甚至让约束关系错误。我的建议是建一个单位换算参数区把所有物理量先换算成统一的kW·h再写进约束和目标函数。第二个坑是SOC初值设置不当。如果初始SOC设得太高调度周期内储能的放电空间很小灵活性被浪费如果设得太低前几个时段储能可能强制充电导致系统运行成本升高。我测试下来初始SOC设在0.2到0.3之间配合周期末回到初始值的约束储能利用效果最好。另外末端SOC等于初始SOC这个约束绝对不能去掉否则优化结果在工程上是不可执行的。第三个坑是惩罚系数太小。弃风弃光惩罚系数如果设置过低优化器宁可弃掉新能源也不愿意让燃气机组调低出力因为调低出力会损失供热。但系数太高也不好会让优化器不计代价地消纳新能源导致P2G和储能在低效区间运行整体经济性反而下降。我常用的标定方法是先用确定性优化跑一遍看看系统在完全没有不确定性约束时会有多少弃风弃光然后按弃风弃光电量的成本损失来反推惩罚系数的量级。第四个坑是CHP的热电比选择。如果直接用定热电比模型CHP的热电产出被锁死储能、P2G这些设备就失去了调节空间协同优化的价值大打折扣。我在第2.2里说过一定要用可变热电比模型让CHP在可运行域内灵活调节热电比这样才能体现出综合能源系统多能互补的调度灵活性。第五个坑是场景削减目标数量不校验。很多人削减完场景就直接丢进模型也不检查削减后的场景集是否保留了原始分布的关键特征。我建议削减前后做一次对比算一下原始场景和削减后场景的风光出力均值与方差如果方差明显变小说明削减把很多波动性场景削掉了模型会低估不确定性得到的结果偏乐观。这时候要增加目标场景数量或调整削减算法。5.3 参数敏感性分析的正确打开方式最后聊一下参数敏感性分析。很多人在论文里做敏感性分析就是把某个参数改几个值看看总成本怎么变然后画一条曲线了事。这种做法的信息量很低审稿人也一眼能看出来是凑数。我自己的经验是抓住三个最有工程意义的参数做分析不确定性预算Γ、储能容量、P2G容量。不确定性预算Γ的敏感性分析最有意思。把Γ从0确定性优化逐渐增加到24最保守工况画出运行成本和新能源消纳率的变化曲线。Γ小的时候系统运行成本低但消纳率也低Γ增大会让成本上升但消纳率也会随之改善。曲线中会出现一个拐点过了这个拐点后继续增大Γ成本上升很快但消纳率几乎不再提升说明不确定性预算已经超出了实际需要的保守程度这个拐点附近就是调度决策的最佳保守区间。储能容量和P2G容量的敏感性分析则要成对做。把储能容量从基准值的0.5倍增加到2倍同时记录系统总成本和消纳率的变化你会看到典型的边际递减规律容量增大初期效果明显但超过某个值之后收益趋缓经济上的最优容量就在这个拐点附近。P2G容量也是类似的逻辑。这两个参数的联合敏感性分析可以直接为系统规划和设备扩容提供参考比单参数扫描有意义得多。从我自己复现这个课题的经历来看最值得投入精力的不是把模型写得更复杂而是把不确定性的场景质量做好。场景质量直接决定了后续所有优化结果的可靠性一个分布建模失真的场景集无论求解器多强大、优化算法多先进最后产出的都只是看起来合理的数字。用场景法把基线跑通再往鲁棒优化方向深挖这条路走下来无论是做课题还是发论文都能拿得出手。