源荷随机特征下热电联供微网日前调度：场景建模与Matlab实现

1. 项目核心热电联供微网为什么绕不开“随机性”做热电联供微网优化的人大概率都遇到过这个问题模型建得漂漂亮亮约束写得整整齐齐结果一放到实际场景里光伏出力对不上、负荷预测偏了系统运行成本直接飙升甚至出现切负荷的极端情况。究其根源就是我们在建模时把“源”和“荷”都当成了确定的量来对待。所谓“考虑源荷随机特征”本质上是承认这样一个事实光伏和风电出力受天气影响负荷受用户行为影响它们都不是一条固定的曲线而是一个随机过程。今天中午的光伏出力可能是满发也可能是云层遮挡后的三成出力明天早晨的负荷可能是尖峰也可能平缓得像条直线。如果优化模型只按一条预测曲线来安排机组出力和储能充放策略一旦实际值偏离预测值整个调度方案就失去了最优性甚至失去可行性——这可不是小事热电联供系统里电和热深度耦合电跟不上可以买电热跟不上就得启动备用锅炉成本翻着跟头往上涨。这个课题的典型应用场景包括园区型综合能源系统、楼宇冷热电联供、农村或海岛独立微网等。研究的核心产物是一套“能扛住不确定性”的日前调度方案通常在Matlab中完成建模与仿真验证。对于电气工程、控制工程或能源动力方向的研究生以及做微网能量管理的工程师来说这是一个非常经典的切入点——既有理论深度随机优化建模又有工程落地性Matlab代码可复现还可以向鲁棒优化、分布鲁棒、多智能体博弈等方向自然延伸。这篇文章我就从建模思路、随机特征刻画、优化模型构建、Matlab代码实现到调试经验把整个项目完整地拆一遍。2. 源荷随机特征到底难在哪三个不确定性来源的建模2.1 风电与光伏出力随机性的数学刻画先说源侧。风电出力取决于风速而风速的统计特性在工程上常用两参数威布尔分布来描述其概率密度函数为f(v) (k/λ) · (v/λ)^(k-1) · exp(−(v/λ)^k)其中k是形状参数λ是尺度参数。知道风速分布之后还要经过风机功率曲线映射到电功率。这里有三个关键区间切入风速以下不出力额定风速到切出风速之间满发中间段近似线性爬升。很多初学者直接把风速均值代入功率曲线算出力这是典型的确定性思维——风速均值对应的不是“平均出力”因为功率曲线的非线性会把偏差放大。光伏出力相对好处理一些。光照强度可以用Beta分布描述光伏电池的输出功率与光照强度近似呈线性关系再乘上温度修正系数和转换效率即可。实际操作中我会建议直接基于历史出力数据做非参数估计或者用分时段的光照强度Beta分布拟合比全局一个分布更精细。这里有个非常容易踩的坑源侧随机性往往存在时空相关性。一个园区里的两台风机距离不远风速高度相关同一时刻的光照和温度也有耦合关系。如果只做独立抽样生成的风光场景会过于“分散”优化结果会偏保守。处理相关性常用Copula函数或Cholesky分解前者适合描述非线性相关性后者简单高效但只处理线性相关。对入门项目来说先用相关性系数矩阵加Cholesky分解就够了没必要一上来就上Copula。2.2 负荷预测误差才是主战场负荷侧的随机性本质上来自负荷预测误差。绝大多数研究不会直接对负荷历史数据建立复杂的随机模型而是假设负荷预测误差服从正态分布即P_load_real P_load_forecast ε ε ~ N(0, σ²)σ的取值通常取预测点负荷值的1%~5%。这个思路背后有一个工程合理性负荷预测算法已经很成熟误差特性相对稳定用正态分布近似在大多数场景下是够用的。但要注意一点不同时段的负荷预测误差方差是不一样的。早晚高峰时段用户行为变化大预测误差方差也大深谷时段负荷稳定误差方差小。如果你全程用一个固定σ生成出来的场景在峰时段显得过于“乐观”在谷时段又过于“悲观”最终调度方案会在不必要的地方留过多备用。更讲究一点的做法是把一天按峰平谷划分成3~4个时段每个时段用不同的σ再用基于历史预测残差的统计方法估计。我在实际做这个项目时就吃过固定σ的亏换了分时段变方差之后场景质量明显提升目标函数值也更有区分度。2.3 场景生成与削减从1000个到10个有了概率分布模型之后如何把这些随机性“塞进”优化模型里最主流的方法是场景法思路很直接从分布中大量抽样生成N个场景然后把每个场景都当作一个确定性的“世界”优化模型要同时满足所有场景下的约束、最小化期望成本这就是所谓的随机规划。具体操作用蒙特卡洛抽样就行每个场景包含24小时的光伏出力、风电出力、电负荷和热负荷时序。N取多少合适如果直接把1000个原始场景丢进优化模型模型规模会爆炸——每个场景都是一套约束副本1000个场景意味着模型里有几万个变量和约束。所以必须做场景削减把1000个场景浓缩成10~20个典型场景同时尽量保留原始场景集的概率分布特征。场景削减最常用的方法是同步回代消除法大致逻辑是计算所有场景两两之间的距离把距离最近的那个场景“合并”进另一个场景被合并场景的概率累加到保留场景上迭代直到场景数降到目标值。Matlab里自己写这个逻辑也就二三十行代码也可以用聚类的方法比如k-means或者k-medoids把场景聚成若干类取聚类中心作为典型场景。区别在于同步回代保留的是原始场景中的点物理意义更清晰k-means得到的是“虚拟场景”可能不符合物理特性。我个人的建议是如果是写论文用同步回代法理由是方法经典、有据可查如果是做工程快速验证用聚类法更快。3. 优化模型怎么建热电联供微网的核心约束与目标函数3.1 决策变量与目标函数先把模型框架说清楚。以最常见的园区热电联供微网为例设备包括燃气轮机或内燃机作为热电联产机组、余热锅炉、燃气锅炉、电储能、蓄热罐、光伏和风机外加与外部电网的交互连接。决策变量分两类。连续变量包括各机组出力、储能充放功率、购售电功率等整数变量包括机组启停状态、储能充放状态防止同时充放等。有整数变量问题就从线性规划变成了混合整数线性规划MILP这也是热电联供微网优化最标准的数学形态。目标函数通常是单目标最小化综合运行成本包含以下几个部分燃料成本燃气轮机和燃气锅炉消耗天然气的费用是成本大头运行维护成本各设备按出力线性计算维护费用购电成本扣除售电收益从电网购电按分时电价计费如果允许余电上网售电电价单独计碳排放成本按系统总排放量折算或者设置碳排放配额弃风弃光惩罚为了鼓励消纳可再生能源可以在目标函数里对弃风弃光加惩罚项。目标函数写成数学形式是min Σ (C_fuel C_om C_grid C_carbon C_penalty)所有成本的期望值都在场景概率加权下计算。也就是每个场景算一个成本最终目标值等于各场景成本乘以该场景概率之后的累加值。这样优化出来的调度方案不是针对某一个场景最优而是对所有场景整体期望成本最优——这正是随机规划区别于确定性优化的核心。3.2 电功率平衡与热功率平衡约束条件分三大块功率平衡、设备模型、储能模型。先看电功率平衡P_gas P_pv P_wind P_dch − P_ch P_buy − P_sell P_load其中P_gas是热电联产机组的电出力P_pv和P_wind是风光出力场景给定P_dch和P_ch是电储能放电和充电功率P_buy和P_sell是购电和售电功率。注意这个平衡约束在每个场景的每个时段都要成立。热功率平衡同样关键H_gas H_boiler H_tank_dch − H_tank_ch H_load其中H_gas是机组的热出力H_boiler是燃气锅炉出力H_tank_dch和H_tank_ch是蓄热罐的放热和蓄热功率。热负荷同样有场景数据。热电联产机组的特殊之处在于它同时产生电和热而电出力和热出力之间存在耦合关系。以最常见的抽凝式机组为例运行可行域是一个由电出力、热出力构成的凸多边形可用一组线性不等式来描述。简化模型中通常用“热电比”来近似H_gas R_chp · P_gasR_chp是热电比。但这只适用于背压机组。抽凝机组的灵活性更大热出力可以在一定范围内调节代价是电出力的上限随热出力变化。严谨的做法是引入可行域的顶点组合建模即把机组运行区域划分成若干个子区域每个子区域用一组线性约束描述。这块如果做不好后续怎么调都容易出问题。3.3 储能约束最容易写错的环节储能模型是初学者最容易出问题的地方。电储能的约束包括容量状态递推SOC(t1) SOC(t) η_ch · P_ch(t) − P_dch(t)/η_dch充放电功率上限0 ≤ P_ch(t) ≤ P_ch_max · u_ch(t)0 ≤ P_dch(t) ≤ P_dch_max · u_dch(t)充放状态互斥u_ch(t) u_dch(t) ≤ 1SOC上下限SOC_min ≤ SOC(t) ≤ SOC_max周期约束SOC(24) SOC(0)调度周期末回到初始电量前四条大家一般不会漏但第五条“周期约束”经常被忽略。如果不加这一条优化器会倾向在周期末把储能放空因为这样成本最小但这不符合实际运行要求——第二天还得继续用。加上周期约束之后储能实际上是在一个“日循环”内做能量套利结论更有工程意义。蓄热罐的结构类似只是把充电效率换成热损失系数。蓄热罐的损耗比电储能小很多但储热密度低、体积大这也是热电联供系统里“热电解耦”的重要手段——在电价低谷时段多产热储起来高峰时段不启动高成本热源。很多论文的重点就是研究怎么通过蓄热罐把原本捆绑在一起的电出力和热出力“解耦”提高系统运行灵活性。除上述约束外机组爬坡约束、出力上下限、购售电功率上限、与电网交互的容量限制等都是必要组成。这些内容全部加进去之后模型就具备了一个标准MILP问题的完整形态。到这里理论上已经把“源荷随机特征”通过场景法“翻译”成了确定性的多场景模型下一步就是把它变成Matlab代码。4. Matlab代码实现从数据结构到求解流程4.1 代码整体架构拿到一个这样的课题不要上来就撸代码。先想清楚数据结构。我的习惯是把整个工程拆成四个模块参数初始化模块、场景生成模块、优化模型构建模块、结果分析模块。文件命名用数字前缀控制执行顺序比如01_param_init.m基础参数、电价、负荷数据、设备参数02_scenario_gen.m生成风光负荷的不确定场景并削减03_build_model.m调用建模工具箱建立优化模型04_solve_and_plot.m求解并把结果画出来四步之间通过工作区的变量衔接。虽然可以用脚本一把梭但长期维护的角度看拆成函数更合理——因为场景削减方法和模型约束经常会调整拆开之后改起来不容易出错。4.2 场景生成的Matlab实现场景生成模块是体现“源荷随机特征”的重头戏。以最简单的独立抽样为例伪代码如下% 参数设置 N_scen 1000; % 初始场景数 T 24; % 调度时段数(小时) N_reduce 10; % 削减后场景数 % 风速抽样 - 风电出力 % 假设风速服从威布尔分布shapek, scalelambda wind_speed wblrnd(k, lambda, N_scen, T); % 风机功率曲线映射(分段函数) P_wind wind_curve(wind_speed, P_wind_rated, v_in, v_r, v_out); % 光照抽样 - 光伏出力 % 假设光照强度服从Beta分布用betarnd抽样 irradiance betarnd(a_t, b_t, N_scen, T); % 简化写法 P_pv eff_pv * S_area .* irradiance; % eff_pv含温度和损耗修正 % 负荷抽样 load_error normrnd(0, sigma_t, N_scen, T); P_load repmat(P_load_forecast, N_scen, 1) load_error; H_load repmat(H_load_forecast, N_scen, 1) ... normrnd(0, sigma_heat, N_scen, T);场景削减的同步回代算法实现也不复杂。核心逻辑是计算场景距离矩阵循环找最小距离删除概率最小的那个场景把它的概率加到距离最近的场景上。用Matlab写出来大约就是这样一个循环function [scen_reduced, prob_reduced] scenario_reduction(scen, prob, N_target) % scen: n_scen x T 的矩阵, prob: 各场景概率 n_scen size(scen, 1); while n_scen N_target % 计算场景间欧氏距离矩阵 D pdist2(scen, scen); D(D 0) inf; % 排除自身 % 找距离最小的场景对 [~, idx] min(D(:)); [i, j] ind2sub([n_scen, n_scen], idx); % 合并: 将概率小的场景概率加到另一个场景上 if prob(i) prob(j) prob(j) prob(i) prob(j); scen(i, :) []; prob(i) []; else prob(i) prob(i) prob(j); scen(j, :) []; prob(j) []; end n_scen n_scen - 1; end scen_reduced scen; prob_reduced prob / sum(prob); % 归一化 end这段代码虽然能跑但是每次循环都重新算距离矩阵复杂度是O(n²×削减次数)1000个场景削减到10个循环将近一千次耗时很高。工程优化手段是只算一次距离矩阵之后每次合并只更新受影响的行和列。但是好在场景削减是一次性的离线工作跑个十几秒甚至几十秒都能接受所以够用就行不必过度优化。投影到图形上场景削减这件事可以用一句很直白的话总结把一大片场景云的“形状”用十几个代表性点子表示出来让这些点子的加权平均效果尽量逼近整片云。配了削减前后的对比图审稿人看一眼就知道你做的是随机优化而不是拿着一条曲线硬算。4.3 模型构建与求解模型构建我用的是Matlab里通用的建模工具箱加外接求解器的路子。对MILP而言三个要点决策变量声明、约束写入、求解器参数设置。先定义变量。二进制变量直接声明连续变量要注意给出上下界上下界写得好不好直接影响求解速度% 声明变量(以场景s下的t时段为例) % 连续变量 P_chp sdpvar(T, N_scen, full); % 机组电出力 H_chp sdpvar(T, N_scen, full); % 机组热出力 P_buy sdpvar(T, N_scen, full); % 购电功率 P_sell sdpvar(T, N_scen, full); % 售电功率 SOC sdpvar(T1, N_scen, full); % 电储能SOC % 二进制变量 u_chp binvar(T, N_scen, full); % 机组启停 u_ch binvar(T, N_scen, full); % 储能充电状态 u_dch binvar(T, N_scen, full); % 储能放电状态然后写约束。这里提醒一个关键点随机规划的约束必须对每个场景分别加不能只加一次。很多新手在for循环外面写约束结果模型只覆盖了第一个场景求解结果出来怎么看都不对检查半天才发现是循环范围的问题。正确写法是对每个场景循环加约束Constraints []; for s 1:N_scen for t 1:T % 电功率平衡(每个场景每个时段) Constraints [Constraints, ... P_chp(t,s) P_pv(t,s) P_wind(t,s) ... P_dch(t,s) P_buy(t,s) - P_ch(t,s) - P_sell(t,s) P_load(t,s)]; % 其他约束... end % 储能SOC递推与周期约束 Constraints [Constraints, SOC(1,s) SOC_0]; Constraints [Constraints, SOC(T1,s) SOC_0]; Constraints [Constraints, ... SOC(t1,s) SOC(t,s) eta_ch * P_ch(t,s) - P_dch(t,s)/eta_dch]; end目标函数用场景概率加权Objective 0; for s 1:N_scen cost_s sum(price_gas * (fuel_chp(:,s) fuel_boiler(:,s)) ... % 燃料费 C_om_chp .* P_chp(:,s) ... % 维护费 price_buy(t) .* P_buy(:,s) - price_sell(t) .* P_sell(:,s)); % 购售电费 Objective Objective prob(s) * cost_s; end求解时直接调用内置的整数线性规划求解器也可以处理中小规模问题如果问题规模比较大比如场景数超过20、时段细化到15分钟内置求解器会明显变慢。这时建议用外部商用求解器配合求解器参数把MIP gap设置为0.5%或1%能大幅缩短求解时间。需要注意外部求解器需要额外安装并且要配置好Matlab调用路径这一步很多人卡住其实就是环境变量的路径问题。5. 典型参数与结果解读怎么判断模型跑得对不对5.1 典型算例参数参数设计是论文或项目报告里最容易被追问的部分。说几个我认可的典型配置调度周期24小时步长1小时共24个时段初始场景数1000个削减后保留10~20个分时电价峰时段10:00~15:0018:00~21:00高电价谷时段23:00~7:00低电价其余为平段峰谷比约3:1燃气轮机额定容量比如1MW电效率35%热电比1.2简化模型燃气锅炉效率90%电储能容量500kWh最大充放电功率200kW充放电效率0.95蓄热罐容量1000kWh热损系数0.01/h负荷曲线取某典型季节的负荷数据电负荷峰值约800kW热负荷峰值约600kW。注意负荷数据不能直接拿一条曲线就完事。上文说过负荷场景要叠加预测误差而且峰平谷时段的误差方差要区别对待。这些负荷基础数据最好用自己实测或者公开的数据集不要凭空捏造一条形状好看但没有代表性的曲线。5.2 确定性方案与随机方案的对比跑完模型之后最重要的验证工作是做对比分析。最基本的对照是确定性模型只取负荷预测值和风光预测值不生成场景与随机规划模型的对比。两者的目标函数形式一样差异只在数据输入——确定性模型相当于场景数N1且概率为1。对比结果通常呈现两个规律第一随机规划得到的期望成本略高于确定性规划。这是正常的因为随机规划要兼顾所有场景必然会在部分场景下牺牲目标值的“最优性”来换取可行性和稳健性。学者把这叫做“随机解的期望价值损失”。你要做的是把这个差额控制在可接受范围比如5%以内说明模型既稳健又不过度保守。第二更具说服力的是“事后仿真”验证把两个模型求出的调度方案放到一组未参与优化的真实场景比如100个测试场景中回代计算看哪个方案的平均运行成本更低、哪些方案出现了失负荷或热不足。常见的“意外结果”是确定性方案在测试场景中有若干个场景无法满足功率平衡——因为它没有为不确定性留余地而随机方案的失负荷概率显著较低。这就是随机优化价值的直接证据。这个对比思路大家在写报告或论文时一定要保留它是整个项目最核心的产出。5.3 敏感性分析与模型验证敏感性分析同样重要。常见的分析维度包括预测误差方差大小、场景削减后的保留数量、电价波动幅度、储能容量等。做敏感性分析有个长远的好处可以验证模型是不是“对参数太敏感”如果一个模型在σ从2%变到3%时调度结果就发生剧烈变化说明模型存在稳定性问题需要检查约束或者定价机制是否合理。场景数量的敏感性有讲究。保留5个场景时期望成本可能偏乐观且不稳定保留15个场景后结果逐渐稳定再增加到30个结果几乎不变但求解时间成倍增加。做一组“期望成本—求解时间—场景数”的折中曲线选拐点作为最终场景数这是比较专业的做法。6. 常见问题与排查实录6.1 模型不可行MILP最常见的坑新模型第一次求解报“infeasible”大概率是以下原因之一电/热平衡约束永远不可能同时满足比如总发电容量小于所有场景中最大负荷或者蓄热罐容量小于热负荷峰谷差需求。排查方法很简单先解一个不带任何整数约束的LP松弛版本看松弛问题是否可行。如果LP都不可行那就不是整数变量搞的鬼而是连续约束本身矛盾。充放互斥约束写得不对。比如u_ch和u_dch分别乘在充电和放电功率上但忘了加“不能同为1”的约束这不会导致不可行但会让模型出现“同时充放”这种物理上不可能的“能量对消”现象浪费效率项还让SOC曲线看起来很奇怪。周期约束和初始SOC冲突。比如要求SOC(24)等于初始值但储能容量小、负荷又大在给定的充放功率范围内根本不可能在24小时内回到初始值。这时候要么放宽周期约束为“SOC(24)在初始值的±10%范围内”要么增加储能容量。遇到不可行建议先做不可行性诊断——用求解器自带的不可行性分析功能定位冲突约束。Matlab里可以通过合理注释约束二分法缩小冲突范围。我实测下来最有效的调试方法是先把所有约束全部注释掉然后一批一批加回来每加一批就求解一次。找到让模型从“可行”变“不可行”的那批约束问题就直接暴露了。6.2 求解时间爆炸如何降低复杂度热电联供微网模型一旦场景数和设备数多起来求解时间会非常难看。几个提升求解效率的实用手段变量域收紧给每个连续变量设置紧凑的上下界比如SOC的初始范围而不是默认的-inf到inf。变量界越紧分支定界搜索的剪枝效率越高求解速度提升非常明显。减少整数变量储能充放互斥状态可以用充电功率和放电功率本身非负且乘积为零来近似但更稳妥的是用连续变量加一个很小的惩罚项替代整型变量或者用“多参数规划”思想简化。如果二进制变量实在消不掉考虑是否能合并一些时段——并不是所有设备都需要小时级启停决策。使用求解器热启动把上一个场景迭代的解作为下一个场景的初始解可以大幅减少搜索时间。调整MIP gap从默认的1e-4放宽到1e-2或1%牺牲千分之一的精度换来几个数量级的提速做敏感性分析的时候尤其香。6.3 场景削减导致极端场景丢失这是一个非常隐蔽但后果严重的问题。同步回代削减倾向保留“主流”场景而把概率小但出力极端的场景合并掉。后果是优化结果对未来可能出现的极端天气连续阴雨天、极端高温毫无准备。我在一次复现中遇到的情况是削减到10个场景时模型给出的方案在正常天气下很漂亮但一遇到连续阴雨天的测试场景蓄电池直接放空、系统切负荷。应对方法有三种削减后人工注入几个极端场景赋予一个小概率比如0.1%让模型对极端情况保持敏感性在约束里额外加一个“失负荷率上限”约束从模型层面保证不能过度牺牲可靠性改用分布鲁棒优化框架用一组统计矩信息约束所有可能的分布而不是依赖具体场景。这个方法理论更硬但实现复杂度也上一个台阶前期可以先从场景法起步。6.4 热负荷与电负荷的时间错配问题热负荷的峰谷时段和电负荷往往不一致北方冬季的供热负荷全天较高白天电负荷高但热负荷可能有午间低谷晚上电负荷低谷但热负荷仍然很高。很多第一次做热电联供的人会轻视热平衡约束结果热负荷峰值时段计算出的蓄热罐容量非常大或者燃气锅炉需要常年高负荷运行整体经济性很差。解决方向是引入“热电解耦”思路蓄热罐在其中扮演桥梁角色在电负荷高峰时段让机组优先发电、多产热存入蓄热罐在电负荷低谷时段用储热来满足供热甚至可以让机组降出力运行把宝贵的低谷期电能用在给储能充电上。蓄热罐容量参数、热电比和分时电价之间的匹配关系是这类系统实现经济性的核心杠杆。7. 项目扩展方向与个人体会跑通一个基本的随机规划版本只是第一步。这个项目本身有相当好的扩展性我列出几条我自己尝试过或想尝试的方向供参考不确定集建模升级把场景法替换为鲁棒优化盒式/椭球式不确定集或分布鲁棒优化重点对比保守性与计算复杂度多时间尺度协同日前随机调度 日内滚动修正的双层框架日内层利用实时预测信息修正日前计划更贴近工程实际引入需求响应资源把可平移负荷、可中断负荷作为灵活性资源纳入模型源荷两侧协同优化电-热综合储能容量规划在运行优化基础上叠加投资规划模型做两阶段优化回答“储电和储热配多大比例最划算”碳交易机制建模把碳配额和碳交易价格纳入目标函数分析不同碳价对机组运行方式的影响。最后说一点个人体会。做这类课题切忌把注意力全放在算法迭代和解算技巧上。算法是“术”问题建模才是“道”——尤其是源荷随机特征的刻画它直接决定了你优化出来的方案在实际天气和负荷波动下是否真的扛得住。我自己踩过的最大坑就是把大量的时间花在调求解器参数上反复调gap、调初始解最后发现真正的问题是光伏场景生成时没有考虑云层遮挡的时序相关性生成的光伏场景过于乐观导致调度方案里光伏消纳率虚高。所以给后续做这个项目的朋友一个建议先花足够的时间把随机特征建准、把场景生成好验证场景的统计特性与历史数据一致再进入优化求解环节。这一步投入的精力回报率是最高的。从建模到Matlab实现再到结果分析与调试整套流程走下来你对“什么叫不确定性环境下做决策”这件事会有一个远比课本更深刻的理解。这个项目不复杂但足够完整非常适合作为综合能源系统方向的第一篇完整仿真研究。