做综合能源系统协同优化的那段时间我每天跟两样东西较劲一个是电、热、气三条能量流之间的耦合关系另一个就是新能源出力的不确定性。这两个问题叠在一起才是综合能源系统调度真正麻烦的地方。刚开始我图简单直接用风电光伏出力预测的确定值建模代码跑得飞快调度方案看起来也漂亮但一到实际落地就露馅了——光伏午后实际出力比预测低了近两成储能被迫提前放空夜间只能高价买电补缺口。所以这次我把计及新能源出力不确定性的电气设备综合能源系统协同优化完整做了一遍用Matlab实现了场景法随机优化的全流程包括不确定性场景生成、场景削减、设备协同建模、Yalmip建模求解和结果分析整套代码可以直接改参数复用到不同规模的园区综合能源系统上。这篇文章我尽量把建模思路、代码架构和踩坑经历都讲透方便你自己复现和改造。1. 先想清楚一件事不确定性为什么是个真问题1.1 综合能源系统的耦合结构决定了它牵一发动全身综合能源系统之所以叫综合是因为电、热、气三类能源之间不是各走各的而是通过设备实现强耦合。以最常见的园区级系统为例燃气轮机或内燃机组成的CHP机组是核心枢纽它烧天然气同时发电和产热电锅炉和热泵把电能转化为热能燃气锅炉是热力侧的备选电储能和热储能分别在两个网络里做缓冲。这样一来电力侧的波动会通过CHP机组直接传导到热力侧——风机出力减了电功率缺额由CHP多发电补上CHP多发电又带来更多余热热力侧就得相应调整锅炉出力或储能充放。这种耦合关系放大了不确定性的影响范围。风光出力不确定表面上只影响电力平衡实际上一路传导最终影响的是整个系统的运行成本和安全裕度。传统电力系统调度只需考虑单一能量流的平衡综合能源系统则要同时盯住多个能量流和多个时间尺度的耦合约束这也是协同优化比单能系统调度难一个量级的原因。1.2 风光的出力误差到底有多大先说一个比较扎心的事实风电短期预测误差普遍在10%到20%光伏由于受云层影响在阴晴交替天气下的误差可能更高。这个误差不是随机小扰动而是有时间结构的——午间光伏高发时段误差大夜间小风时段的风电预测偏差常常超过均值的一半。预测误差还会跟负荷峰谷叠加傍晚电负荷爬坡时如果恰好赶上风电预测偏高实际出力不足就会造成严重的功率缺额。我在实际算例中做过统计一个100MW光伏电站15分钟级预测的平均绝对误差约8%但极端时段误差可以到25%以上。用确定性模型处理这些数据等于假设误差恒为零调度结果自然经不起实际运行检验。更重要的是综合能源系统里有储能和CHP这些可调设备它们的调节能力是有限的、有成本的。确定性模型会充分利用预测值把CHP安排在最优工况点储能安排在最优充放曲线但实际风光一偏离这些最优工况可能瞬间变得不可行。1.3 确定性模型的三个典型翻车场景我见过的确定性模型翻车案例总结下来有三个典型套路。第一个是光伏高估场景。确定性模型看到午间光伏预测出力充足就让燃气锅炉保持低载、储能多充电备用。结果当天云层偏厚光伏实际出力只有预测的80%午间热负荷却因为气温骤降升高了这时候热力侧供给不足要么紧急启动燃气锅炉满发要么从电网买高价电来电锅炉补热运行成本直接超预算。第二个是风电低估场景。夜间风电预测出力偏低确定性模型安排CHP高发。实际风速上来了风电大发电力侧供大于求但CHP的热出力又没法快速降下来——CHP的热电比运行域约束摆在那里降电功率必须同步降热功率可夜间热负荷还居高不下。最后只能弃风弃风损失记在账上。第三个是爬坡失配场景。确定性模型按预测曲线安排储能充放和CHP爬坡没有预留调节裕量。实际风光在时段间剧烈波动时CHP爬坡速率不够储能容量也不够系统只能切负荷或者甩负荷。这三个场景我都在算例里复现过结论是一致的不考虑不确定性优化结果的经济性指标很好看但实际可执行性很差。协同优化四个字重心不在协同而在优化时必须把不确定性当成一等公民对待。2. 不确定性建模场景法和鲁棒区间两条路怎么选2.1 场景法把不确定变成多个确定的未来应对不确定性最直观的思路就是把一个不确定的未来变成一堆可能的确定性未来每个未来带一个概率权重。这就是场景法Stochastic Programming的核心思想。具体到风电光伏出力不确定性的场景生成常规做法分三步。第一步是获取预测值和误差分布。预测值可以从历史数据、数值天气预报或简单的时间序列模型得到误差分布则用预测残差的统计特性来拟合。风电残差通常用正态分布近似光伏残差因为受云层遮挡影响用Beta分布拟合效果更好。如果想更精细可以直接用历史残差做核密度估计不预设分布形态但计算量会大一些。第二步是采样生成场景。用蒙特卡洛方法从误差分布中抽样叠加到预测曲线上就得到一组出力场景。这里有个细节风光出力在时间上是相关的不能逐时段独立采样否则生成的场景时序上剧烈跳变明显不符合实际。我的做法是使用多变量正态分布利用历史数据的协方差矩阵刻画时段间相关性或者用一阶自回归模型来生成时序场景。Matlab里用mvnrnd函数就可以生成多变量正态随机数方便。第三步是场景削减。直接采样1000个场景每个场景都带24个时段的变量这个随机优化模型规模会非常大求解时间难以接受。场景削减的思路是保留概率分布特征的前提下用较少的场景代替原始场景集。常见方法有同步回代消除算法SBR和基于概率距离的快速前向选择Fast Forward Selection。核心逻辑每次都删掉一个对概率分布影响最小的场景并把它的概率累加到距离它最近的场景上直到场景数达到预设值。我在实际项目里一般把1000个初始场景削减到20到50个计算精度损失不到3%求解时间却能降低一个数量级。2.2 鲁棒优化不赌概率只保底线跟场景法不同鲁棒优化不问每个场景的概率是多少而是问最坏情况是什么。它把不确定性描述成一个集合比如风电出力在预测值加减一定偏差的区间内波动优化目标是在这个集合中所有可能实现下都可行的方案里找一个成本最低的。具体到综合能源系统常用的是带预算约束的盒式不确定集每个时段风光出力可以在预测值上下偏离一定比例但所有时段的总偏离量被一个参数Γ限制住。Γ就是鲁棒调节参数Γ0时退化为确定性模型Γ越大保护的极端情况越多优化结果越保守。这个方法的数学处理通常需要把约束条件中的不确定参数通过对偶理论转化为确定性约束推导过程稍复杂但一旦转化完成还是一个规模可控的线性规划或混合整数线性规划问题。2.3 我的选择建议两条路线没有绝对优劣取决于工程场景。场景法的优点是可解释性强能看到每个典型场景下的运行策略还能算出期望成本、最坏成本等统计指标缺点是模型规模随场景数增大且场景削减过程本身有信息损失。鲁棒优化的优点是计算效率高不需要生成和削减场景缺点是结果偏保守——毕竟它把所有偏差都当成同概率对待而实际上极端天气出现的概率很低。我在工程中更倾向于混搭策略主体用场景法做随机优化得到经济性最优的调度方案然后取一个较小的鲁棒预算Γ或者对关键时段比如负荷峰值时段加上鲁棒边界检查方案在最坏情况下的安全性。如果边界检查不过再针对性增加备用或调整储能策略。这种随机优化为主、鲁棒校验兜底的做法兼顾了经济性和可靠性也是实际工程里比较容易落地的思路。3. 设备建模和目标函数协同优化到底在优化什么3.1 关键设备模型怎么建综合能源系统的设备模型精度太低没意义精度太高又非线性严重没法用商用求解器处理。我的经验是用混合整数线性规划MILP能表达的程度去建模把非线性关系做分段线性化或逻辑约束处理。CHP机组是最关键的设备。它有两个输出电功率和热功率这两个输出之间不是固定的比例关系而是落在一个可行运行域里。常见做法是把可行域近似为一个凸多边形用一组线性不等式约束来描述。比如某型号内燃机CHP电功率范围20到100kW热功率范围30到120kW电热运行域的顶点坐标通过厂家数据拟合得到使用时直接把这些顶点代入形成凸组合约束。CHP还需要建模启停状态、最小启停时间和爬坡约束这些都需要引入二进制变量。储能设备用荷电状态SOC递推方程描述SOC(t1) SOC(t) η_c·P_c(t)·Δt - P_d(t)·Δt/η_d。其中η_c和η_d分别是充放电效率P_c和P_d是充放电功率。这里要特别注意的是SOC递推方程在整个调度周期上形成强耦合是模型中最容易出数值问题的地方。还有燃气锅炉的出力上下限和爬坡约束购售电联络线的功率约束以及电网交互电价的分时变化处理这些相对简单直接写线性约束就行。3.2 目标函数不止是省钱协同优化的目标函数最基础的是运行成本最小化。运行成本包含几块天然气购气成本CHP和燃气锅炉共用、从电网购电的成本分时电价、设备运行维护成本按出力线性折算、启停成本固定成本项。如果系统包含弃风弃光和失负荷的情况还要加上相应的惩罚项。工程上还有一个现实约束目标函数里要不要加碳排放现在很多项目确实有碳约束需求最简单的方式是在目标函数中增加一个碳排放成本项排放系数乘以天然气消耗量和购电量再乘以碳价。这样做的好处是不用额外加约束求解难度不变还能调节碳价来观察不同政策下的设备运行策略变化。我实际建模时还会加一个小的松弛技巧在目标函数中加入对储能充放次数的惩罚项或对CHP出力波动的惩罚项防止优化结果出现频繁切换的锯齿型出力曲线这种曲线理论上最优实际运行根本执行不了。惩罚项的系数要调到一个合适区间太小不起作用太大又会让经济性变差。3.3 约束体系平衡约束和运行约束缺一不可约束体系分三层。第一层是能量平衡约束每个时段电功率平衡——CHP发电、光伏出力、风电出力、储能放电加上购电等于电负荷加储能充电加售电热功率平衡——CHP产热、燃气锅炉产热、热储能放热等于热负荷加热储能充热。这一层约束在场景法里要遍历每个场景逐时段写。第二层是设备运行约束前面已经描述包括出力上下限、爬坡约束、储能SOC递推、启停逻辑等。第三层是系统安全约束比如联络线功率限值、关键设备备用容量、失负荷比例上限等。还有一个经常被忽略的约束是热电比耦合约束——CHP的电出力和热出力必须落在可行运行域内这个约束在协同优化里是承上启下的关键把它写错了整个模型的物理意义就不对了。4. Matlab实现求解器选择和代码架构4.1 为什么选Yalmip加商用求解器Matlab上做优化建模最省力的路径是Yalmip工具箱加Gurobi或CPLEX求解器。Yalmip是一个建模语言层它把优化问题的变量定义、目标函数和约束表达得跟数学公式非常接近然后自动转换成求解器需要的标准形式调用底层求解器解算。用Yalmip写优化模型比直接调求解器API或手写单纯形法快一个数量级。求解器方面Gurobi和CPLEX都能处理MILP问题对中小规模的综合能源系统模型来说两者性能差距在实际使用中基本感觉不到。选型主要看许可证获取方便程度Gurobi有学术授权个人学习和科研用途申请很方便CPLEX在IBM收购后授权渠道有变化但也是成熟稳定的选择。如果你的模型规模很大比如场景数上千或者系统节点数上百建议两个求解器都装上实测对比再决定用哪个。求解器相关参数里最值得关注的是MIP Gap容忍度我一般设为1%到2%没必要追求0.1%以下——为了那点精度多花几倍求解时间工程上不划算。4.2 代码架构按数据、模型、求解三块拆写Matlab优化代码最容易犯的错误是把所有东西塞在一个脚本里变量铺满整个文件改参数要翻半天。我建议按三个模块组织参数数据模块、模型构建模块、求解与后处理模块。参数数据模块用一个函数或脚本集中定义所有系统参数负荷曲线、风光预测值、设备参数、电价气价、场景数据。我习惯用struct保存系统参数比如sys.chp.Pmax、sys.ess.eta_c这样模型构建模块引用起来清晰而且参数修改不影响模型逻辑。场景数据单独用一个三维数组保存scen(n, t, type)第一维是场景编号第二维是时段第三维是类型风电、光伏或温度相关负荷。模型构建模块是核心它接收参数数据内部定义Yalmip变量、写约束和目标函数最后调用optimize求解。这个模块最好写成一个独立的函数输入参数、返回结果结构体。后处理模块负责从求解结果中提取各设备出力曲线计算成本构成画图。三个模块分开以后换一组数据只改参数模块换求解器只动求解部分改约束条件只动模型构建模块整体可维护性好很多。4.3 变量定义和约束骨架的写法下面给一段变量定义和约束搭建的示意代码这是基于常规实现思路的通用写法你可以直接参考。% 基础参数 T 24; % 调度时段 N_s 20; % 场景数量 dt 1; % 时间间隔单位小时 % Yalmip变量所有决策变量都是24维或N_s x 24维 P_chp sdpvar(1, T, full); % CHP电出力 H_chp sdpvar(1, T, full); % CHP热出力 H_boiler sdpvar(1, T, full); % 燃气锅炉热出力 P_es_c sdpvar(1, T, full); % 储能充电功率 P_es_d sdpvar(1, T, full); % 储能放电功率 SOC sdpvar(1, T1, full); % 储能荷电状态 P_buy sdpvar(1, T, full); % 购电功率 P_sell sdpvar(1, T, full); % 售电功率 z_chp binvar(1, T, full); % CHP启停状态 % 约束集合 Constraints []; for t 1:T % 电功率平衡含场景期望 for s 1:N_s Constraints [Constraints, ... P_chp(t) P_pv_scen(s,t) P_wt_scen(s,t) P_es_d(t) P_buy(t) ... P_load(t) P_es_c(t) P_sell(t)]; end % 热功率平衡 Constraints [Constraints, ... H_chp(t) H_boiler(t) H_hs_d(t) H_load(t) H_hs_c(t)]; % CHP可行域约束多边形逼近 % 这里用顶点凸组合的方式表达 end % 目标函数期望运行成本最小 Objective gas_cost grid_cost om_cost penalty_cost;变量定义里有一个细节值得说如果决策变量在全部时段和全部场景下是一致的也就是单阶段随机优化先决策、后观察场景那么决策变量只需要定义T维不需要定义N_s乘T维。场景只在约束的右侧参数风光出力中体现约束需要遍历所有场景但决策变量不随场景变化。这一点很多初学者会搞混把决策变量也定义成场景相关的——那样模型就变成每个场景独立决策完全失去了随机优化的意义结果也完全没有鲁棒性。5. 代码模块拆解从场景生成到求解后处理5.1 场景生成与削减的代码逻辑场景生成部分我用的是多元正态分布加协方差矩阵的方式打一个简单的示例% 假设风电预测值 wt_forecast (1 x T) % 历史残差协方差矩阵 Sigma (T x T) n_scen_raw 1000; % 初始场景数 rng(42); residuals mvnrnd(zeros(1, T), Sigma, n_scen_raw); wt_scen_raw repmat(wt_forecast, n_scen_raw, 1) residuals; % 限幅处理保证出力在[0, 额定容量]范围内 wt_scen_raw min(max(wt_scen_raw, 0), wt_capacity);这里有一个工程细节直接对预测残差采样然后叠加生成的场景里个别点可能超出物理允许范围必须做限幅处理。如果不限幅后面优化时会出现无解或者解明显不合理的情况。光伏场景同理但要注意光伏夜间出力一定为0所以残差采样只适用于白天时段夜间时段直接置0。场景削减我常用同步回代消除的简化实现。核心思路是迭代计算场景两两之间的概率距离找到被删除后分布变化最小的场景删除它并把概率合并到最近邻场景。概率距离一般用欧氏距离乘以两者概率计算公式不复杂但循环次数多用向量化写法能明显加速。削减后的场景保存在scen_prob每个场景的概率和scen_data每个场景的风光出力曲线里后续模型构建直接引用这两个变量。5.2 目标函数和约束写入的完整套路写目标函数时我有一个建议把成本拆成几项分别累加不要混在一起写一行。这样做的好处是后处理阶段可以直接分析各项成本占比定位系统运行的瓶颈。比如燃气成本项gas_cost gas_price * sum(gas_consum_chp gas_consum_boiler); % gas_consum_chp 可以由CHP电出力和热出力通过效率曲线换算这里CHP的燃料消耗量计算要特别留意它不是一个单纯的线性函数而是电出力和热出力共同决定的。工程上常用简化方式把燃料消耗量近似为电出力的一次函数加热出力的一次函数即gas_cons a*P_chp b*H_chp c*z_chp系数通过厂家数据拟合。这个线性近似在MILP框架里非常有效既保证了精度又不引入非线性。约束写入时场景相关的约束必须用循环遍历。这里有一个性能优化技巧如果能用Yalmip的矩阵运算批量生成约束尽量多用矩阵形式少用for循环。比如风光场景数据是一个N_s x T的矩阵对应的平衡约束可以用repmat配合矢量化操作一次生成全部约束。场景数较大时这种写法和for循环相比求解时间能差好几倍因为约束生成阶段的Matlab开销从O(N_s)循环降到了矩阵运算级别。5.3 求解设置与结果后处理求解调用很简单ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.01); result optimize(Constraints, Objective, ops);verbose建议设成1或2能实时看到求解器日志。如果模型无解日志里会显示约束冲突信息这在调试时非常有用。求解返回的result结构体里有solveTime、problem等字段检查problem是否为0非0则对应不同的错误类型。求解完成后用value()函数提取各决策变量数值P_chp_opt value(P_chp); SOC_opt value(SOC);提取完画图时我习惯把风光场景的包络线同时画出来——就是把所有削减后场景的最大值和最小值画成阴影区域然后把优化出的CHP、储能曲线画在图上。这样一眼就能看出调度方案在不确定性影响下的稳健程度比只画期望场景曲线有说服力得多。6. 跑通代码之后结果分析和参数敏感性6.1 确定性模型和随机优化模型的差距在哪同样一套系统参数分别用确定性模型只用预测曲线和场景法随机优化模型跑一遍结果对比非常有意思。确定性模型的总成本通常略低大约低3%到8%但它是虚假的优——它假设风光完全按预测出力不需要额外的备用调节。随机优化模型因为要考虑多种场景会主动预留一部分调节能力比如储能不会在某个时段全部充满而是保留一部分容量应对风光低出力的场景CHP也不会盯死在经济最优工况点而是稍微偏离一点以保留爬坡能力。这些预留体现在成本上就是随机优化比确定性模型多出来的那部分。工程上应该这样看待这笔钱它是为了换取方案在实际运行中可执行性所支付的保险费用。我在汇报时通常会展示一个对比表确定性方案在预测场景下成本低但把它放到真实场景集合中测试失负荷概率可能达到10%以上随机优化方案在所有测试场景下都能满足约束虽然期望成本略高但没有严重的失负荷风险。把这个对比做出来决策者很容易理解为什么必须考虑不确定性。6.2 鲁棒调节参数和场景数怎么影响结果如果你做了鲁棒校验或者直接用鲁棒优化会观察到Γ参数对成本和保守性的敏感度曲线。Γ从0增大时总成本单调上升但失负荷或缺电风险单调下降。曲线通常在中间区域有一个性价比最高的拐点再增加Γ成本上升明显但风险下降趋缓。找出这个拐点的方法是扫描Γ从0到最大预算值记录每个取值下的成本然后对比边际成本增量。我习惯做5到10个样本点画出一条平滑曲线拐点位置一目了然。场景数的影响则体现在两方面一是求解时间二是解的质量。场景数为5时求解极快但结果可能只覆盖了少数典型运行状态某个特殊场景下约束会失效场景数增加到30后期望成本和约束满足率基本稳定再增加场景数求解时间暴涨结果几乎不变。建议不同项目用一组场景数做一次收敛性测试找到自己模型的甜点场景数标准做法是测试5、10、20、50、100几个档次画出成本和求解时间随场景数的变化曲线。6.3 结果展示和报告呈现的技巧最后汇报时最有力的展示是三件套第一张图是系统电功率平衡图把CHP、风电、光伏、储能、购电和负荷画成堆叠面积图直观展示各设备协同出力关系第二张图是热功率平衡图展示CHP产热、锅炉产热、热储能和热负荷的匹配关系第三张图是储能SOC曲线看它的充放策略是否合理——一个明显特征是储能倾向于在低价时段充电、高价时段放电但如果SOC曲线频繁满充满放说明目标函数里缺少对电池寿命的考虑需要加惩罚项。文字结论部分我会突出三个数字期望运行成本、最坏场景下的成本、失负荷风险概率。这三个数字分别对应经济性、稳健性和安全性决策者看这三个数字就能快速判断方案是否可接受。7. 调试踩坑记录这五个问题我花了一周7.1 量纲混用让目标函数失真第一次跑模型时我发现在不同设备参数里混用了kW和MW两种单位——光伏和风机容量用的是MWCHP和负荷用的则是kW。目标函数里所有成本项都按不同量纲累加结果燃气成本比购电成本低了三四个数量级优化结果直接让CHP满发、完全不买电明显不合理。排查半天才意识到是量纲问题。解决方法是写一个单位换算函数所有数据进出模型前统一换算成同一单位制并且在参数定义处加注释说明。7.2 场景削减后概率没有归一化场景削减算法本身实现没问题但我在做概率合并时忘了对削减后的场景概率集合做归一化处理。结果所有场景概率之和是0.96而不是1期望成本被系统性低估了4%。这个偏差看起来不大但在做方案对比时足够造成误判——两个方案的期望成本差距本来就只有几个百分点归一化误差直接掩盖了真实差异。从那以后我在每个处理阶段都加一个校验assert(abs(sum(scen_prob) - 1) 1e-8)。7.3 CHP可行域多边形出现不符合实际的顶点组合CHP可行域用多边形顶点组合约束表达时我从厂家样本数据里直接取了几个工况点作为顶点结果这些点构成的多边形包含了物理上不可能达到的工况——低电功率配高热量输出的区域。模型确实有解但调度方案在某些时段要求CHP发出超出能力范围的组合功率实际根本做不到。这个问题的根源是厂家样本点覆盖不全多边形凸包范围超出了真实可行域。修正方法是把样本点画出来人工剔除极端点或者用凹包算法紧贴实际数据边界。只要这一步出了问题后面所有优化结果都是空中楼阁。7.4 大M法的M取值过大导致数值问题处理CHP启停逻辑和储能充放互斥时我用了大M法线性化逻辑约束。第一个版本把M设成10000想着肯定够大结果求解器频繁报数值警告有些约束在容差范围内处于松弛和满足的边界解的质量很不稳定。后来把M改成跟该约束相关的最大物理量纲——比如电功率相关的约束M取系统最大供电能力的1.1倍热功率约束同理——数值问题立刻消失。大M不是越大越好取一个物理上有意义的紧上界才是正确的做法。7.5 热负荷需求曲线的口径没对齐最后这个坑藏在数据准备阶段热负荷数据来自供热公司的日累计热量表单位是GJ而我模型中热功率单位是kW时间尺度是小时换算关系是1kW·h等于0.0036GJ。我第一版直接把GJ数值当kW用导致热负荷被放大近三百倍燃气锅炉永远满发CHP热功率约束全部顶到上限系统完全丧失灵活性。这类单位问题在数据交换时特别容易发生强烈建议每个数据源接入模型前写一个明确的单位映射表并且在模型构建后用物理常识做一次校验算出来的CHP热出力是否在设备容量范围内、储能SOC初末值是否接近等。这些检查几分钟就能做完能省掉几天的排查时间。整套流程做下来我的体会是计及不确定性的协同优化难点不在算法本身有多高深而在于每个环节都要想清楚物理意义和工程约束。场景生成、削减、设备建模、MILP求解、结果分析每一环都有它容易出错的地方把这一条链路走通并形成自己的代码模板后面应对不同规模和配置的综合能源系统项目基本就是换参数、调约束的事很快就能出结果。 计及新能源出力不确定性的综合能源系统协同优化与Matlab实现
做综合能源系统协同优化的那段时间我每天跟两样东西较劲一个是电、热、气三条能量流之间的耦合关系另一个就是新能源出力的不确定性。这两个问题叠在一起才是综合能源系统调度真正麻烦的地方。刚开始我图简单直接用风电光伏出力预测的确定值建模代码跑得飞快调度方案看起来也漂亮但一到实际落地就露馅了——光伏午后实际出力比预测低了近两成储能被迫提前放空夜间只能高价买电补缺口。所以这次我把计及新能源出力不确定性的电气设备综合能源系统协同优化完整做了一遍用Matlab实现了场景法随机优化的全流程包括不确定性场景生成、场景削减、设备协同建模、Yalmip建模求解和结果分析整套代码可以直接改参数复用到不同规模的园区综合能源系统上。这篇文章我尽量把建模思路、代码架构和踩坑经历都讲透方便你自己复现和改造。1. 先想清楚一件事不确定性为什么是个真问题1.1 综合能源系统的耦合结构决定了它牵一发动全身综合能源系统之所以叫综合是因为电、热、气三类能源之间不是各走各的而是通过设备实现强耦合。以最常见的园区级系统为例燃气轮机或内燃机组成的CHP机组是核心枢纽它烧天然气同时发电和产热电锅炉和热泵把电能转化为热能燃气锅炉是热力侧的备选电储能和热储能分别在两个网络里做缓冲。这样一来电力侧的波动会通过CHP机组直接传导到热力侧——风机出力减了电功率缺额由CHP多发电补上CHP多发电又带来更多余热热力侧就得相应调整锅炉出力或储能充放。这种耦合关系放大了不确定性的影响范围。风光出力不确定表面上只影响电力平衡实际上一路传导最终影响的是整个系统的运行成本和安全裕度。传统电力系统调度只需考虑单一能量流的平衡综合能源系统则要同时盯住多个能量流和多个时间尺度的耦合约束这也是协同优化比单能系统调度难一个量级的原因。1.2 风光的出力误差到底有多大先说一个比较扎心的事实风电短期预测误差普遍在10%到20%光伏由于受云层影响在阴晴交替天气下的误差可能更高。这个误差不是随机小扰动而是有时间结构的——午间光伏高发时段误差大夜间小风时段的风电预测偏差常常超过均值的一半。预测误差还会跟负荷峰谷叠加傍晚电负荷爬坡时如果恰好赶上风电预测偏高实际出力不足就会造成严重的功率缺额。我在实际算例中做过统计一个100MW光伏电站15分钟级预测的平均绝对误差约8%但极端时段误差可以到25%以上。用确定性模型处理这些数据等于假设误差恒为零调度结果自然经不起实际运行检验。更重要的是综合能源系统里有储能和CHP这些可调设备它们的调节能力是有限的、有成本的。确定性模型会充分利用预测值把CHP安排在最优工况点储能安排在最优充放曲线但实际风光一偏离这些最优工况可能瞬间变得不可行。1.3 确定性模型的三个典型翻车场景我见过的确定性模型翻车案例总结下来有三个典型套路。第一个是光伏高估场景。确定性模型看到午间光伏预测出力充足就让燃气锅炉保持低载、储能多充电备用。结果当天云层偏厚光伏实际出力只有预测的80%午间热负荷却因为气温骤降升高了这时候热力侧供给不足要么紧急启动燃气锅炉满发要么从电网买高价电来电锅炉补热运行成本直接超预算。第二个是风电低估场景。夜间风电预测出力偏低确定性模型安排CHP高发。实际风速上来了风电大发电力侧供大于求但CHP的热出力又没法快速降下来——CHP的热电比运行域约束摆在那里降电功率必须同步降热功率可夜间热负荷还居高不下。最后只能弃风弃风损失记在账上。第三个是爬坡失配场景。确定性模型按预测曲线安排储能充放和CHP爬坡没有预留调节裕量。实际风光在时段间剧烈波动时CHP爬坡速率不够储能容量也不够系统只能切负荷或者甩负荷。这三个场景我都在算例里复现过结论是一致的不考虑不确定性优化结果的经济性指标很好看但实际可执行性很差。协同优化四个字重心不在协同而在优化时必须把不确定性当成一等公民对待。2. 不确定性建模场景法和鲁棒区间两条路怎么选2.1 场景法把不确定变成多个确定的未来应对不确定性最直观的思路就是把一个不确定的未来变成一堆可能的确定性未来每个未来带一个概率权重。这就是场景法Stochastic Programming的核心思想。具体到风电光伏出力不确定性的场景生成常规做法分三步。第一步是获取预测值和误差分布。预测值可以从历史数据、数值天气预报或简单的时间序列模型得到误差分布则用预测残差的统计特性来拟合。风电残差通常用正态分布近似光伏残差因为受云层遮挡影响用Beta分布拟合效果更好。如果想更精细可以直接用历史残差做核密度估计不预设分布形态但计算量会大一些。第二步是采样生成场景。用蒙特卡洛方法从误差分布中抽样叠加到预测曲线上就得到一组出力场景。这里有个细节风光出力在时间上是相关的不能逐时段独立采样否则生成的场景时序上剧烈跳变明显不符合实际。我的做法是使用多变量正态分布利用历史数据的协方差矩阵刻画时段间相关性或者用一阶自回归模型来生成时序场景。Matlab里用mvnrnd函数就可以生成多变量正态随机数方便。第三步是场景削减。直接采样1000个场景每个场景都带24个时段的变量这个随机优化模型规模会非常大求解时间难以接受。场景削减的思路是保留概率分布特征的前提下用较少的场景代替原始场景集。常见方法有同步回代消除算法SBR和基于概率距离的快速前向选择Fast Forward Selection。核心逻辑每次都删掉一个对概率分布影响最小的场景并把它的概率累加到距离它最近的场景上直到场景数达到预设值。我在实际项目里一般把1000个初始场景削减到20到50个计算精度损失不到3%求解时间却能降低一个数量级。2.2 鲁棒优化不赌概率只保底线跟场景法不同鲁棒优化不问每个场景的概率是多少而是问最坏情况是什么。它把不确定性描述成一个集合比如风电出力在预测值加减一定偏差的区间内波动优化目标是在这个集合中所有可能实现下都可行的方案里找一个成本最低的。具体到综合能源系统常用的是带预算约束的盒式不确定集每个时段风光出力可以在预测值上下偏离一定比例但所有时段的总偏离量被一个参数Γ限制住。Γ就是鲁棒调节参数Γ0时退化为确定性模型Γ越大保护的极端情况越多优化结果越保守。这个方法的数学处理通常需要把约束条件中的不确定参数通过对偶理论转化为确定性约束推导过程稍复杂但一旦转化完成还是一个规模可控的线性规划或混合整数线性规划问题。2.3 我的选择建议两条路线没有绝对优劣取决于工程场景。场景法的优点是可解释性强能看到每个典型场景下的运行策略还能算出期望成本、最坏成本等统计指标缺点是模型规模随场景数增大且场景削减过程本身有信息损失。鲁棒优化的优点是计算效率高不需要生成和削减场景缺点是结果偏保守——毕竟它把所有偏差都当成同概率对待而实际上极端天气出现的概率很低。我在工程中更倾向于混搭策略主体用场景法做随机优化得到经济性最优的调度方案然后取一个较小的鲁棒预算Γ或者对关键时段比如负荷峰值时段加上鲁棒边界检查方案在最坏情况下的安全性。如果边界检查不过再针对性增加备用或调整储能策略。这种随机优化为主、鲁棒校验兜底的做法兼顾了经济性和可靠性也是实际工程里比较容易落地的思路。3. 设备建模和目标函数协同优化到底在优化什么3.1 关键设备模型怎么建综合能源系统的设备模型精度太低没意义精度太高又非线性严重没法用商用求解器处理。我的经验是用混合整数线性规划MILP能表达的程度去建模把非线性关系做分段线性化或逻辑约束处理。CHP机组是最关键的设备。它有两个输出电功率和热功率这两个输出之间不是固定的比例关系而是落在一个可行运行域里。常见做法是把可行域近似为一个凸多边形用一组线性不等式约束来描述。比如某型号内燃机CHP电功率范围20到100kW热功率范围30到120kW电热运行域的顶点坐标通过厂家数据拟合得到使用时直接把这些顶点代入形成凸组合约束。CHP还需要建模启停状态、最小启停时间和爬坡约束这些都需要引入二进制变量。储能设备用荷电状态SOC递推方程描述SOC(t1) SOC(t) η_c·P_c(t)·Δt - P_d(t)·Δt/η_d。其中η_c和η_d分别是充放电效率P_c和P_d是充放电功率。这里要特别注意的是SOC递推方程在整个调度周期上形成强耦合是模型中最容易出数值问题的地方。还有燃气锅炉的出力上下限和爬坡约束购售电联络线的功率约束以及电网交互电价的分时变化处理这些相对简单直接写线性约束就行。3.2 目标函数不止是省钱协同优化的目标函数最基础的是运行成本最小化。运行成本包含几块天然气购气成本CHP和燃气锅炉共用、从电网购电的成本分时电价、设备运行维护成本按出力线性折算、启停成本固定成本项。如果系统包含弃风弃光和失负荷的情况还要加上相应的惩罚项。工程上还有一个现实约束目标函数里要不要加碳排放现在很多项目确实有碳约束需求最简单的方式是在目标函数中增加一个碳排放成本项排放系数乘以天然气消耗量和购电量再乘以碳价。这样做的好处是不用额外加约束求解难度不变还能调节碳价来观察不同政策下的设备运行策略变化。我实际建模时还会加一个小的松弛技巧在目标函数中加入对储能充放次数的惩罚项或对CHP出力波动的惩罚项防止优化结果出现频繁切换的锯齿型出力曲线这种曲线理论上最优实际运行根本执行不了。惩罚项的系数要调到一个合适区间太小不起作用太大又会让经济性变差。3.3 约束体系平衡约束和运行约束缺一不可约束体系分三层。第一层是能量平衡约束每个时段电功率平衡——CHP发电、光伏出力、风电出力、储能放电加上购电等于电负荷加储能充电加售电热功率平衡——CHP产热、燃气锅炉产热、热储能放热等于热负荷加热储能充热。这一层约束在场景法里要遍历每个场景逐时段写。第二层是设备运行约束前面已经描述包括出力上下限、爬坡约束、储能SOC递推、启停逻辑等。第三层是系统安全约束比如联络线功率限值、关键设备备用容量、失负荷比例上限等。还有一个经常被忽略的约束是热电比耦合约束——CHP的电出力和热出力必须落在可行运行域内这个约束在协同优化里是承上启下的关键把它写错了整个模型的物理意义就不对了。4. Matlab实现求解器选择和代码架构4.1 为什么选Yalmip加商用求解器Matlab上做优化建模最省力的路径是Yalmip工具箱加Gurobi或CPLEX求解器。Yalmip是一个建模语言层它把优化问题的变量定义、目标函数和约束表达得跟数学公式非常接近然后自动转换成求解器需要的标准形式调用底层求解器解算。用Yalmip写优化模型比直接调求解器API或手写单纯形法快一个数量级。求解器方面Gurobi和CPLEX都能处理MILP问题对中小规模的综合能源系统模型来说两者性能差距在实际使用中基本感觉不到。选型主要看许可证获取方便程度Gurobi有学术授权个人学习和科研用途申请很方便CPLEX在IBM收购后授权渠道有变化但也是成熟稳定的选择。如果你的模型规模很大比如场景数上千或者系统节点数上百建议两个求解器都装上实测对比再决定用哪个。求解器相关参数里最值得关注的是MIP Gap容忍度我一般设为1%到2%没必要追求0.1%以下——为了那点精度多花几倍求解时间工程上不划算。4.2 代码架构按数据、模型、求解三块拆写Matlab优化代码最容易犯的错误是把所有东西塞在一个脚本里变量铺满整个文件改参数要翻半天。我建议按三个模块组织参数数据模块、模型构建模块、求解与后处理模块。参数数据模块用一个函数或脚本集中定义所有系统参数负荷曲线、风光预测值、设备参数、电价气价、场景数据。我习惯用struct保存系统参数比如sys.chp.Pmax、sys.ess.eta_c这样模型构建模块引用起来清晰而且参数修改不影响模型逻辑。场景数据单独用一个三维数组保存scen(n, t, type)第一维是场景编号第二维是时段第三维是类型风电、光伏或温度相关负荷。模型构建模块是核心它接收参数数据内部定义Yalmip变量、写约束和目标函数最后调用optimize求解。这个模块最好写成一个独立的函数输入参数、返回结果结构体。后处理模块负责从求解结果中提取各设备出力曲线计算成本构成画图。三个模块分开以后换一组数据只改参数模块换求解器只动求解部分改约束条件只动模型构建模块整体可维护性好很多。4.3 变量定义和约束骨架的写法下面给一段变量定义和约束搭建的示意代码这是基于常规实现思路的通用写法你可以直接参考。% 基础参数 T 24; % 调度时段 N_s 20; % 场景数量 dt 1; % 时间间隔单位小时 % Yalmip变量所有决策变量都是24维或N_s x 24维 P_chp sdpvar(1, T, full); % CHP电出力 H_chp sdpvar(1, T, full); % CHP热出力 H_boiler sdpvar(1, T, full); % 燃气锅炉热出力 P_es_c sdpvar(1, T, full); % 储能充电功率 P_es_d sdpvar(1, T, full); % 储能放电功率 SOC sdpvar(1, T1, full); % 储能荷电状态 P_buy sdpvar(1, T, full); % 购电功率 P_sell sdpvar(1, T, full); % 售电功率 z_chp binvar(1, T, full); % CHP启停状态 % 约束集合 Constraints []; for t 1:T % 电功率平衡含场景期望 for s 1:N_s Constraints [Constraints, ... P_chp(t) P_pv_scen(s,t) P_wt_scen(s,t) P_es_d(t) P_buy(t) ... P_load(t) P_es_c(t) P_sell(t)]; end % 热功率平衡 Constraints [Constraints, ... H_chp(t) H_boiler(t) H_hs_d(t) H_load(t) H_hs_c(t)]; % CHP可行域约束多边形逼近 % 这里用顶点凸组合的方式表达 end % 目标函数期望运行成本最小 Objective gas_cost grid_cost om_cost penalty_cost;变量定义里有一个细节值得说如果决策变量在全部时段和全部场景下是一致的也就是单阶段随机优化先决策、后观察场景那么决策变量只需要定义T维不需要定义N_s乘T维。场景只在约束的右侧参数风光出力中体现约束需要遍历所有场景但决策变量不随场景变化。这一点很多初学者会搞混把决策变量也定义成场景相关的——那样模型就变成每个场景独立决策完全失去了随机优化的意义结果也完全没有鲁棒性。5. 代码模块拆解从场景生成到求解后处理5.1 场景生成与削减的代码逻辑场景生成部分我用的是多元正态分布加协方差矩阵的方式打一个简单的示例% 假设风电预测值 wt_forecast (1 x T) % 历史残差协方差矩阵 Sigma (T x T) n_scen_raw 1000; % 初始场景数 rng(42); residuals mvnrnd(zeros(1, T), Sigma, n_scen_raw); wt_scen_raw repmat(wt_forecast, n_scen_raw, 1) residuals; % 限幅处理保证出力在[0, 额定容量]范围内 wt_scen_raw min(max(wt_scen_raw, 0), wt_capacity);这里有一个工程细节直接对预测残差采样然后叠加生成的场景里个别点可能超出物理允许范围必须做限幅处理。如果不限幅后面优化时会出现无解或者解明显不合理的情况。光伏场景同理但要注意光伏夜间出力一定为0所以残差采样只适用于白天时段夜间时段直接置0。场景削减我常用同步回代消除的简化实现。核心思路是迭代计算场景两两之间的概率距离找到被删除后分布变化最小的场景删除它并把概率合并到最近邻场景。概率距离一般用欧氏距离乘以两者概率计算公式不复杂但循环次数多用向量化写法能明显加速。削减后的场景保存在scen_prob每个场景的概率和scen_data每个场景的风光出力曲线里后续模型构建直接引用这两个变量。5.2 目标函数和约束写入的完整套路写目标函数时我有一个建议把成本拆成几项分别累加不要混在一起写一行。这样做的好处是后处理阶段可以直接分析各项成本占比定位系统运行的瓶颈。比如燃气成本项gas_cost gas_price * sum(gas_consum_chp gas_consum_boiler); % gas_consum_chp 可以由CHP电出力和热出力通过效率曲线换算这里CHP的燃料消耗量计算要特别留意它不是一个单纯的线性函数而是电出力和热出力共同决定的。工程上常用简化方式把燃料消耗量近似为电出力的一次函数加热出力的一次函数即gas_cons a*P_chp b*H_chp c*z_chp系数通过厂家数据拟合。这个线性近似在MILP框架里非常有效既保证了精度又不引入非线性。约束写入时场景相关的约束必须用循环遍历。这里有一个性能优化技巧如果能用Yalmip的矩阵运算批量生成约束尽量多用矩阵形式少用for循环。比如风光场景数据是一个N_s x T的矩阵对应的平衡约束可以用repmat配合矢量化操作一次生成全部约束。场景数较大时这种写法和for循环相比求解时间能差好几倍因为约束生成阶段的Matlab开销从O(N_s)循环降到了矩阵运算级别。5.3 求解设置与结果后处理求解调用很简单ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.01); result optimize(Constraints, Objective, ops);verbose建议设成1或2能实时看到求解器日志。如果模型无解日志里会显示约束冲突信息这在调试时非常有用。求解返回的result结构体里有solveTime、problem等字段检查problem是否为0非0则对应不同的错误类型。求解完成后用value()函数提取各决策变量数值P_chp_opt value(P_chp); SOC_opt value(SOC);提取完画图时我习惯把风光场景的包络线同时画出来——就是把所有削减后场景的最大值和最小值画成阴影区域然后把优化出的CHP、储能曲线画在图上。这样一眼就能看出调度方案在不确定性影响下的稳健程度比只画期望场景曲线有说服力得多。6. 跑通代码之后结果分析和参数敏感性6.1 确定性模型和随机优化模型的差距在哪同样一套系统参数分别用确定性模型只用预测曲线和场景法随机优化模型跑一遍结果对比非常有意思。确定性模型的总成本通常略低大约低3%到8%但它是虚假的优——它假设风光完全按预测出力不需要额外的备用调节。随机优化模型因为要考虑多种场景会主动预留一部分调节能力比如储能不会在某个时段全部充满而是保留一部分容量应对风光低出力的场景CHP也不会盯死在经济最优工况点而是稍微偏离一点以保留爬坡能力。这些预留体现在成本上就是随机优化比确定性模型多出来的那部分。工程上应该这样看待这笔钱它是为了换取方案在实际运行中可执行性所支付的保险费用。我在汇报时通常会展示一个对比表确定性方案在预测场景下成本低但把它放到真实场景集合中测试失负荷概率可能达到10%以上随机优化方案在所有测试场景下都能满足约束虽然期望成本略高但没有严重的失负荷风险。把这个对比做出来决策者很容易理解为什么必须考虑不确定性。6.2 鲁棒调节参数和场景数怎么影响结果如果你做了鲁棒校验或者直接用鲁棒优化会观察到Γ参数对成本和保守性的敏感度曲线。Γ从0增大时总成本单调上升但失负荷或缺电风险单调下降。曲线通常在中间区域有一个性价比最高的拐点再增加Γ成本上升明显但风险下降趋缓。找出这个拐点的方法是扫描Γ从0到最大预算值记录每个取值下的成本然后对比边际成本增量。我习惯做5到10个样本点画出一条平滑曲线拐点位置一目了然。场景数的影响则体现在两方面一是求解时间二是解的质量。场景数为5时求解极快但结果可能只覆盖了少数典型运行状态某个特殊场景下约束会失效场景数增加到30后期望成本和约束满足率基本稳定再增加场景数求解时间暴涨结果几乎不变。建议不同项目用一组场景数做一次收敛性测试找到自己模型的甜点场景数标准做法是测试5、10、20、50、100几个档次画出成本和求解时间随场景数的变化曲线。6.3 结果展示和报告呈现的技巧最后汇报时最有力的展示是三件套第一张图是系统电功率平衡图把CHP、风电、光伏、储能、购电和负荷画成堆叠面积图直观展示各设备协同出力关系第二张图是热功率平衡图展示CHP产热、锅炉产热、热储能和热负荷的匹配关系第三张图是储能SOC曲线看它的充放策略是否合理——一个明显特征是储能倾向于在低价时段充电、高价时段放电但如果SOC曲线频繁满充满放说明目标函数里缺少对电池寿命的考虑需要加惩罚项。文字结论部分我会突出三个数字期望运行成本、最坏场景下的成本、失负荷风险概率。这三个数字分别对应经济性、稳健性和安全性决策者看这三个数字就能快速判断方案是否可接受。7. 调试踩坑记录这五个问题我花了一周7.1 量纲混用让目标函数失真第一次跑模型时我发现在不同设备参数里混用了kW和MW两种单位——光伏和风机容量用的是MWCHP和负荷用的则是kW。目标函数里所有成本项都按不同量纲累加结果燃气成本比购电成本低了三四个数量级优化结果直接让CHP满发、完全不买电明显不合理。排查半天才意识到是量纲问题。解决方法是写一个单位换算函数所有数据进出模型前统一换算成同一单位制并且在参数定义处加注释说明。7.2 场景削减后概率没有归一化场景削减算法本身实现没问题但我在做概率合并时忘了对削减后的场景概率集合做归一化处理。结果所有场景概率之和是0.96而不是1期望成本被系统性低估了4%。这个偏差看起来不大但在做方案对比时足够造成误判——两个方案的期望成本差距本来就只有几个百分点归一化误差直接掩盖了真实差异。从那以后我在每个处理阶段都加一个校验assert(abs(sum(scen_prob) - 1) 1e-8)。7.3 CHP可行域多边形出现不符合实际的顶点组合CHP可行域用多边形顶点组合约束表达时我从厂家样本数据里直接取了几个工况点作为顶点结果这些点构成的多边形包含了物理上不可能达到的工况——低电功率配高热量输出的区域。模型确实有解但调度方案在某些时段要求CHP发出超出能力范围的组合功率实际根本做不到。这个问题的根源是厂家样本点覆盖不全多边形凸包范围超出了真实可行域。修正方法是把样本点画出来人工剔除极端点或者用凹包算法紧贴实际数据边界。只要这一步出了问题后面所有优化结果都是空中楼阁。7.4 大M法的M取值过大导致数值问题处理CHP启停逻辑和储能充放互斥时我用了大M法线性化逻辑约束。第一个版本把M设成10000想着肯定够大结果求解器频繁报数值警告有些约束在容差范围内处于松弛和满足的边界解的质量很不稳定。后来把M改成跟该约束相关的最大物理量纲——比如电功率相关的约束M取系统最大供电能力的1.1倍热功率约束同理——数值问题立刻消失。大M不是越大越好取一个物理上有意义的紧上界才是正确的做法。7.5 热负荷需求曲线的口径没对齐最后这个坑藏在数据准备阶段热负荷数据来自供热公司的日累计热量表单位是GJ而我模型中热功率单位是kW时间尺度是小时换算关系是1kW·h等于0.0036GJ。我第一版直接把GJ数值当kW用导致热负荷被放大近三百倍燃气锅炉永远满发CHP热功率约束全部顶到上限系统完全丧失灵活性。这类单位问题在数据交换时特别容易发生强烈建议每个数据源接入模型前写一个明确的单位映射表并且在模型构建后用物理常识做一次校验算出来的CHP热出力是否在设备容量范围内、储能SOC初末值是否接近等。这些检查几分钟就能做完能省掉几天的排查时间。整套流程做下来我的体会是计及不确定性的协同优化难点不在算法本身有多高深而在于每个环节都要想清楚物理意义和工程约束。场景生成、削减、设备建模、MILP求解、结果分析每一环都有它容易出错的地方把这一条链路走通并形成自己的代码模板后面应对不同规模和配置的综合能源系统项目基本就是换参数、调约束的事很快就能出结果。