考虑源荷随机特征的热电联供微网优化：场景法与Matlab实现

1. 项目概述与整体设计思路1.1 核心需求解析前阵子帮一个做综合能源方向的同学看代码发现他卡在“源荷随机性”这个概念上很久。他手里有一份热电厂出力的确定性调度代码跑起来没问题但只要一涉及风光预测误差、负荷波动整个模型就不知道该怎么改。这个项目标题“考虑源荷随机特征的热电联供微网优化研究Matlab代码实现”其实点得很明白不是单纯做个微网调度而是要把不确定性“请进”优化模型里用随机优化的思路去处理源侧和荷侧的波动。先交代一下背景。热电联供Combined Heat and PowerCHP微网是目前园区级能源系统的主流形态核心设备是燃气轮机或内燃机它在发电的同时回收余热供给热负荷能效比单纯的电热分产高出不少。微网里还会配风电、光伏这类可再生能源以及电储能、蓄热罐、燃气锅炉等辅助设备。整个系统的运行目标是在满足电负荷和热负荷的前提下让总运行成本最低。但问题来了——风电和光伏出力是随机波动的负荷也有预测误差如果优化调度时把这些都当成固定值那调度方案在实际执行时可能严重偏离预期极端情况下会导致切负荷或者弃风弃光。所以就有了“考虑源荷随机特征”这个研究方向。这个项目适合谁看两类人。一类是正在做微电网/综合能源系统方向毕业论文的硕士生另一类是刚接触随机优化、想知道怎么用Matlab落地实现的工程师。我会把随机建模的原理、场景生成与削减、优化模型的数学表达、以及YalmipCplex的完整实现逻辑全部讲透最后附上我在调试过程中踩过的坑。1.2 为什么不能只用确定性优化很多人一开始不理解既然调度问题本质是个优化问题我直接把风电出力取预测值、负荷取预测值解出来的结果不也能用吗为什么要折腾随机性我拿一个生活场景类比。你每天早上出门前决定要不要带伞依据是天气预报。如果天气预报说“明天降水概率30%”你大概率不带伞但如果预报说“降水概率80%”你肯定带。确定性优化相当于把天气预报当成“明天一定下雨”或“明天一定不下雨”来处理而随机优化是拿到“降水概率分布”后在所有可能天气下做权衡找到一个平均意义上最优的决策。放到微网调度里道理完全一样。燃气轮机的启停状态、电储能的充放电计划、与上级电网的交互功率这些决策都需要提前制定。如果你只按风电出力的期望值来安排计划实际风电比预期低20%时可能就得被迫高价购电或者切负荷。随机优化会预先把这种偏差考虑进去生成的调度方案在各种可能场景下都留有余地综合表现更稳健。具体来说确定性优化和随机优化的本质区别有三个维度输入数据确定性优化用单点预测值随机优化用概率分布或多场景集求解过程确定性优化只解一次随机优化需要处理多个场景的联合优化结果形态确定性优化给出一个固定调度计划随机优化给出的计划对所有场景期望意义上最优这个项目的核心就是用场景法Scenario-Based来处理随机性先通过概率分布抽样生成大量场景再用场景削减技术挑出有代表性的少数场景最后把所有场景作为一个整体放进优化模型里求解。下面我会一步步拆开讲。2. 源荷随机特征建模的原理与方法2.1 风电出力不确定性的数学描述风电出力的随机性主要来自风速的波动。学术研究和工程实践中普遍用两参数威布尔分布Weibull Distribution来描述风速的概率特性其概率密度函数为f(v) (k / c) * (v / c)^(k-1) * exp(-(v / c)^k)其中v是风速k是形状参数一般取2附近c是尺度参数。得到风速分布后通过风机出力特性曲线将风速映射到出力。标准的映射关系是一个分段函数风速小于切入风速或大于切出风速时出力为0风速在切入风速和额定风速之间时出力近似按线性或三次方关系上升风速在额定风速和切出风速之间时出力为额定功率这里给出一组我在测试中常用的典型参数切入风速3 m/s额定风速12 m/s切出风速25 m/s额定功率1.5 MW。用makedist和random函数在Matlab里抽样每次能得到一组风速样本再转换成功率样本就得到了风电出力场景。2.2 光伏出力不确定性的数学描述光伏出力的随机性主要来自光照强度波动。工程上通常假设光照强度服从贝塔分布Beta Distribution概率密度函数为f(r) Γ(αβ) / (Γ(α) * Γ(β)) * (r / r_max)^(α-1) * (1 - r / r_max)^(β-1)其中r是实际光照强度r_max是最大光照强度α和β是形状参数。光伏出力与光照强度近似成正比再考虑温度对组件效率的影响可以写成功率输出的表达式。在Matlab里面betarnd函数可以直接生成服从贝塔分布的随机数非常方便。2.3 负荷预测误差建模电负荷和热负荷的不确定性不像风光那样来自自然条件而是来自预测模型的误差。工程上一般假设预测误差服从正态分布即P_load P_forecast ε, ε ~ N(0, σ²)σ通常取预测值的3%到10%具体看历史数据的统计结果。这里有个容易忽略的点电负荷误差和热负荷误差不一定独立。在一些场景里电热负荷同时受天气影响比如温度既影响采暖负荷也影响空调负荷所以建模时可以设置一定的相关系数。但为了让代码简洁、便于复现我下面的实现先假设两者独立。2.4 场景生成与削减技术蒙特卡洛抽样与同步回代消除有了概率分布模型下一步就是生成场景。最直接的方法是蒙特卡洛抽样对每个随机变量根据它的分布独立抽样然后把风电、光伏、电负荷、热负荷的抽样结果组合成一个完整的场景向量。抽几千个场景就得到一个覆盖各种可能性的场景集。但场景太多会带来计算问题。假设你有2000个原始场景每个场景引入一组变量和约束优化问题的规模会被撑大几十倍求解时间可能从几秒钟变成几十分钟。而且很多场景之间高度相似——比如两个场景的风电出力偏差只有0.5%把它们同时放进模型里纯属浪费计算资源。这时候就需要场景削减。最常用的是同步回代消除法Simultaneous Backward ReductionSBR。核心思想很简单每次找出一对距离最近的场景删掉其中一个同时把被删场景的概率累加到保留场景上让总概率分布尽量不变。具体算法流程是这样计算所有场景对之间的概率距离常用的是欧氏距离乘以场景概率的加权值找到距离最小的一对场景假设场景i和场景j删掉其中概率较小的那个场景保留另一个把被删场景的概率加到保留下来的场景上重复步骤1到4直到场景数降到预设值比如20个或50个我用Matlab写过这个算法直接循环实现核心代码不超过30行。对于能接受复杂依赖的读者也可以用概率工具箱里的SceneReduction函数。但说实话自己写一遍更能理解它的含义。这里放一组我在代码中实测过的典型参数原始抽样场景数取2000削减后的场景数取20。2000到20看起来削减幅度很大但对比测试显示优化结果与用更多场景比如50个时的结果已经相当接近运行时间却从几分钟降到了几秒钟。实际使用中建议做一次灵敏度分析后面我会详细讲。3. 热电联供微网优化模型的完整构建3.1 系统架构与设备模型先说清楚这个微网系统里都有什么。我的代码实现采用的是一个典型的热电联供微网结构风电WT和光伏PV可再生能源出力具有随机性燃气轮机CHP核心设备同时产电和产热是决定系统经济性的关键燃气锅炉GB补充供热当CHP余热不足时启动电储能ESS存储低价电或风光富余电力高峰时放电蓄热罐HS存储CHP的余热实现热电解耦上级电网允许购电和售电购电价格采用分时电价值得注意的一个细节是热电联产机组的运行方式。CHP机组通常工作在“以热定电”或“以电定热”两种模式之一。本项目采用的是“以热定电”模式先根据热负荷需求确定CHP的产热量再通过热电比约束联动确定发电量。这种方式在实际园区中比较常见因为它能优先保障供热可靠性。代码中我通过热电比参数η_chp把电出力和热出力耦合起来。各设备的关键参数如下表所示设备类型容量效率/关键参数爬坡率限制燃气轮机3 MW电/ 3.6 MW热发电效率40%热电比1.20.5 MW/h燃气锅炉5 MW热效率85%0.8 MW/h电储能2 MWh充放电效率95%SOC范围0.1-0.90.5 MW蓄热罐4 MWh储放热效率90%容量范围0.15-0.950.6 MW风电2 MW--光伏1 MW--3.2 目标函数最低运行成本到底在优化什么优化目标是最小化系统在一个调度周期通常是24小时内的总运行成本。总成本的组成项比较多但每一项都有明确的工程意义我列个清单第一项是购电成本或售电收益。微网与上级电网存在功率交换用分时电价结算。峰时段电价高谷时段电价低这直接影响储能的充放电策略。第二项是燃料成本。燃气轮机和燃气锅炉都烧天然气成本按单位热值价格乘以消耗量计算而消耗量与设备的出力存在明确的效率换算关系。第三项是运维成本。我按设备出力的一定比例折算比如燃气轮机的运维费率取0.05元/kWh风光的运维费率低一些取0.01元/kWh。第四项是弃风弃光惩罚。当系统无法消纳全部可再生能源时会以惩罚系数的形式计入成本这样优化器会尽量避免弃风弃光。最后一项是切负荷惩罚。如果系统功率不足导致必须削减电负荷或热负荷会产生很高的惩罚成本我把系数设置成购电价的几十倍确保优化器只在极端情况下才舍得切负荷。把这几项加总写成一个目标函数就得到了我们的优化目标。特别提醒一下在第二阶段的期望成本计算里每个场景都会算一次从上级电网购电、弃风弃光、切负荷的成本再按场景概率加权求和。这种“第一阶段决策 第二阶段期望调整”的建模方式正是随机规划中典型的两阶段随机优化的结构。3.3 约束条件的体系与细节约束条件是优化模型的主体我拆成几类来说。第一类是功率平衡约束。对于每个时段和每个场景系统的电功率必须满足平衡关系风电出力加光伏出力加CHP发电加电储能放电加购电等于电负荷加电储能充电加售电。热功率同理CHP产热加燃气锅炉产热加蓄热罐放热等于热负荷加蓄热罐充热。平衡约束是所有调度问题的骨架表达式本身不复杂但容易在编写代码时搞混正负号的方向。第二类是设备出力上下限约束。每台设备都有最小技术出力和最大技术出力限制。特别要说的是CHP机组还有最小开机电出力要求低于这个值时机组无法稳定运行反映出燃气轮机不能频繁启停、不能在极低负荷下运行的工程实际。第三类是爬坡约束。燃气轮机和燃气锅炉在相邻时段的出力变化不能超过规定速率。这一点在随机优化里尤其重要因为场景中可能出现相邻时段风电出力大幅跳变的情况如果没有爬坡约束调度方案在物理上根本执行不了。第四类是储能系统约束。电储能的荷电状态State of ChargeSOC在每个时段都有动态递推关系并且SOC必须维持在上下限之间。蓄热罐的SOC递推关系与之类似但放热可以带一定的热损失系数。储能约束的实现是整个代码里最容易出错的地方因为SOC是连续耦合的状态变量需要在所有时段上递推。第五类是备用约束。考虑到源荷随机性系统需要在常规调度基础上留出一定备用容量。我的实现里考虑了两种备用旋转备用应对负荷或风电的突然波动和热备用保障供热可靠性。备用约束是随机优化区别于确定性优化最直观的体现它相当于给所有设备再罩上一层安全网。3.4 求解器选型为什么是Yalmip加Cplex模型建好之后需要一个高效求解器。很多初学者习惯用fminunc或者fmincon直接解非线性问题但这类方法在处理混合整数问题时往往效率很低而且容易陷入局部最优。热电联供微网优化中燃气轮机和锅炉的启停变量是0-1整数变量所以整个问题本质上是一个混合整数线性规划MILP问题。如果目标函数或约束中存在非线性项还需要做线性化处理。在Matlab生态里最顺手的组合是Yalmip工具箱加Cplex求解器或者Gurobi看你的授权情况。Yalmip扮演的是“建模语言”的角色提供了非常直观的语法来表达优化问题Cplex负责在底层做数学求解。我在这个项目里默认使用Yalmip加Cplex用起来大概就是ops sdpsettings(solver, cplex, verbose, 2); optimize(Constraints, Objective, ops);这两行代码看起来轻巧但背后完成了从建模到求解的全部工作。如果你的机器上没有安装Cplex也可以用intlinprog替代Yalmip会自动检测可用的求解器只不过大规模场景下收敛速度会慢一点。4. Matlab代码实现与实操详解4.1 代码整体架构模块化设计整个Matlab工程我按模块拆分成五个主要脚本每个脚本职责单一方便调试和维护data_input.m输入所有基础数据包括设备参数、分时电价、预测的负荷曲线、风光预测出力曲线scenario_generation.m蒙特卡洛抽样生成原始场景然后做场景削减输出代表性场景及概率chp_model.m构建CHP微网的优化模型定义决策变量、目标函数和约束条件solve_optimization.m调用Yalmip和Cplex求解提取结果plot_results.m绘制调度结果图对比不同场景下的运行状态这种模块化拆分最大的好处是定位问题快。比如运行结果异常时如果怀疑是场景生成的问题单独跑第二个脚本检查生成的场景曲线就行不用每次都在整个工程里翻。4.2 场景生成的Matlab实现细节场景生成这一块是整个项目的基石。我简化了核心逻辑把最关键的几行列出来采样的时候风电先对威布尔分布抽样得到风速再通过出力曲线转换成功率光伏直接对贝塔分布抽样负荷做正态分布抽样。每个时段的抽样互相独立把T个时段拼接起来就得到一个完整场景。场景削减的部分我用了同步回代法。中间有个关键细节计算场景距离前要对不同维度的变量做归一化因为风电功率数值大、热负荷数值可能相差更大如果不归一化距离计算会被数值大的维度主导削减出的场景可能丧失多样性。这个坑我一开始踩过后来加上归一化处理后场景质量明显改善。4.3 核心优化模型的Matlab代码实现接下来是重头戏怎么把数学模型变成可求解的Yalmip代码。我按代码段来讲解每一段都有对应的工作机理。首先是决策变量定义。两类变量第一类是第一阶段变量在随机优化里表示需要提前决定、不能随场景改变的量。本项目里包括% 第一阶段变量燃气轮机启停状态、各机组出力、储能充放电计划 z_chp binvar(1, T); % 燃气轮机启停机1表示开机 P_chp sdpvar(1, T); % 燃气轮机发电出力 P_gb sdpvar(1, T); % 燃气锅炉出力 soc_ess sdpvar(1, T); % 电储能荷电状态 soc_hs sdpvar(1, T); % 蓄热罐状态第二类变量属于第二阶段对每个场景单独定义表示运行动态调整量P_wt_s sdpvar(N, T); % 风电实际消纳出力N是场景个数 P_pv_s sdpvar(N, T); % 光伏实际消纳出力 P_grid_s sdpvar(N, T); % 与电网交换功率正为购电负为售电 L_e_cut_s sdpvar(N, T); % 电负荷削减量这里有个初学者容易混淆的概念第一阶段变量在所有场景下取值相同第二阶段变量随场景变化。打个比方第一阶段变量是你出发前就定好的路线第二阶段变量是路上根据实际天气做的调整。建模时这个区别必须对应到代码里否则整个随机优化的结构就瓦解了。其次是约束条件。核心的功率平衡约束用循环来构建Scenarios sdpvar(1, 1); % 用于构建场景集 Constraints []; for t 1:T for k 1:N % 电功率平衡 Constraints [Constraints, ... P_wt_s(k,t) P_pv_s(k,t) P_chp(t) ... P_dis_ess_s(k,t) P_grid_s(k,t) ... L_e(k,t) P_ch_ess_s(k,t) - L_e_cut_s(k,t)]; end end储能SOC的动态递推是另一个核心约束表达的是充电量和放电量与SOC变化的关系。这里要特别注意充放电效率的方向充电时存入的能量不是全额存进去的要乘以充电效率放电时能放出的能量也不是全额放出来电要先乘放电效率再输送到负荷侧。然后还有备用约束。这部分我把每个场景的备用需求表示为一个不等式约束要求系统在该场景下的可用出力之和大于等于负荷与备用之和。这么做能直接体现随机优化对可靠性的保障。最后把这些约束和场景概率加权的期望目标函数一起送入优化求解器Objective W1 W2; % W1为第一阶段成本W2为所有场景下的期望运行成本 ops sdpsettings(solver, cplex, verbose, 2); result optimize(Constraints, Objective, ops);4.4 参数设置与灵敏度分析我实测中发现场景削减后的场景数量对结果的影响不是线性的。场景数从5个增加到20个时目标函数值变化明显从20个增加到50个时变化幅度变得很小而求解时间几乎线性增长。也就是说存在一个“性价比拐点”。针对典型算例20个削减后的场景可以用较快的求解速度和较高的结果稳定性实现平衡。电价结构的设定同样关键。峰谷电价差越大储能系统越有“低储高放”的套利空间。假如峰谷价差从0.4元/kWh扩大到0.8元/kWh储能的日充放电循环次数会从不到1次增加到接近2次整个运行成本的变化也非常明显。这个结论可以和项目里的分时电价表对应上。需要提醒的是初始场景数是另一个容易忽视的参数。我建议初始抽样数不低于1000个否则削减后的场景集可能覆盖不到极端情况导致调度结果过于乐观。5. 结果分析与对比验证5.1 典型日调度结果解读我在标准算例下跑了一次完整调度典型日的预测电负荷峰值出现在晚上7点左右热负荷峰值出现在清晨6点到8点。优化结果显示燃气轮机在大部分时段维持较高出力水平因为它的发电成本低于高峰时段从电网购电的成本。电储能在凌晨谷电时段充满在上午和傍晚的高峰时段放电。蓄热罐的运行策略和电储能正好形成互补——燃气轮机在夜间产生大量余热蓄热罐先把热量存起来白天热负荷上升时再放热实现了热电解耦。特别要关注的是风电和光伏的消纳情况。在大多数场景下风电和光伏的出力都被全额消纳了说明备用约束和储能配置基本合理。但在少数极端场景中比如风速很低而负荷很高的时刻出现了少量切负荷惩罚成本被记入期望成本。真实调度中这种情况的发生概率不高经济上可以被接受。5.2 随机优化与确定性优化结果对比为了验证随机优化的价值我做了对照组把风电、光伏和负荷固定为预测值跑一个确定性优化模型再把这个确定性方案放到随机场景里回测计算它在各个场景下的期望成本。结果很有意思。确定性方案的“名义成本”比随机优化方案低5%到8%看着好像更省钱。但在随机场景回测中确定性方案因为缺少备用裕度经常触发高额的切负荷惩罚或高价购电实际的期望运行成本反而比随机优化方案高出10%到15%。这个对比鲜明地说明了“名义最优”和“真实最优”的差距。如果项目里只做确定性优化最终得到的调度方案在实际运行中可能不只是次优甚至会直接失效。5.3 不同场景数量的灵敏度分析为了确定合理的场景削减数量我做了一组实验分别把削减后场景数设为10、20、30、50记录目标函数值和求解耗时。场景数目标函数值元求解耗时秒结果差异相对50个场景10184203.22.10%20181357.80.52%301807213.50.17%501804128.6基准从表格结果来看场景数从10增加到20带来的精度提升明显求解时间也在可接受范围内20到30的边际收益开始递减到50个场景时结果趋于稳定但耗时翻倍。综合精度和效率20个削减场景是这个测试系统的甜点值。你的系统如果负荷曲线波动更剧烈可以把场景数适当调到30以获取更多稳健性。6. 常见问题与排查技巧实录6.1 Yalmip报错与求解器配置问题Q1Yalmip报“No suitable solver found”怎么办A这通常是系统没装Cplex或者路径没配置好。Cplex安装后在Matlab中执行addpath把对应路径添加上去然后运行yalmiptest进行验证。注意Cplex版本必须和Matlab版本兼容我之前在R2022b上装老版Cplex就出现过兼容警告。如果没有Cplex可以先换用intlinprog顶住。不过为了让求解大型场景更顺畅建议还是解决Cplex的授权和路径问题。Q2求解提示infeasible problem但我觉得模型没问题A大概率是约束矛盾导致的。建议逐段注释约束定位重点检查储能初始电量和最后时段电量约束是否合理、爬坡约束是否和出力上下限冲突、备用约束是否要求过高等。我调试过程中最常碰到的坑是SOC的初始设定和第一时段电平衡约束不匹配把初始SOC从0.2调成0.35后问题就迎刃而解了。6.2 求解时间长或内存溢出随机优化场景数一旦增大变量和约束的规模会瞬间膨胀。如果发现求解时间超过10分钟小规模系统参考先看场景数是不是设得过高。可以从50个场景削减到20个试试。如果模型里包含很多非线性项检查是否做完线性化。还要打开Cplex的输出信息直接看它卡在哪个阶段。有时候瓶颈不在模型大小而在于某个约束写得低效。6.3 结果看起来“不合理”怎么排查结果不合理分两种。第一种是储能基本不工作。原因通常是峰谷价差太小储能套利收益覆盖不了充放电损耗。你可以试着把峰谷价差加大再跑一次如果储能开始工作了说明模型逻辑正确只是经济驱动不足。第二种是燃气轮机一直满发。可能原因是热负荷需求很高“以热定电”模式下为了满足热需求CHP电出力被迫推高。此时需要检查热负荷数据是否过大或者蓄热罐容量是不是偏小、无法在低热负荷时段蓄热来削减CHP出力。还有一种情况是切负荷量明显偏大而且总在同一个时段。这时候去查该时段的备用约束是不是覆盖到了这个场景以及电储能是否在该时段已处于SOC下限无法提供额外支撑。6.4 实操心得汇总从我多次跑这个项目的经验看有几个细节值得专门拿出来说道。场景削减后一定要检查削减前后所有场景的平均期望值是否接近如果偏差超过5%说明削减算法实现有问题或者初始场景数太少生成的代表性场景已经失真。我早期写削减算法时犯过一个错误距离矩阵计算忘了除以样本数导致距离被夸大削减结果偏差很大排查了好一阵才意识到是归一化问题。目标函数的量级要提前预估。各个成本项最好控制在相近数量级如果切负荷惩罚系数比其他成本项大出一百倍以上有时候会导致求解器数值稳定性变差出现奇怪的迭代不收敛问题。我不是说惩罚系数不能大而是说可以先用较小的惩罚系数做一次试算确认模型逻辑正确后再逐步加大到目标值。MATLAB的随机数种子必须固定。记得在调试或者做对比实验前写一行rng(0)否则每次运行的场景都不同对比结果就没有意义了。尤其是做“随机优化vs确定性优化”对照实验时更要在同一组场景下回测不然出来的差异无法归因。最后关于项目扩展我建议对两阶段随机优化熟悉之后可以考虑加入鲁棒优化方法做对比。场景法需要假设概率分布已知如果你手里的历史数据不足分布假设就不太可靠。此时用盒式不确定性集合驱动的鲁棒优化会更有优势。两种方法各有侧重结合起来看问题会更全面。