做热电联供微电网优化的人,应该都碰过这么一个问题:模型里明明写得清清楚楚,一到实际工程里就变了味。光伏出力今天飘了,风说停就停,热负荷也跟着天气和人活动习惯一起抖。要是把这些随机性全丢到确定性的调度模型里,算出来的结果往往只有两个走向——要么成本预测虚低,要么指令发下去根本执行不了。我这次要拆解的项目,标题已经点了题:“考虑源荷随机特征的热电联供微网优化研究(Matlab代码实现)”。核心思路也不是搭一个静态的典型日算完就交差,而是把光伏、风电出力以及电、热负荷的不确定性全数纳入调度决策,用场景化的随机规划框架去求解,最后落到一套可以直接跑的 Matlab 实现上。
这个内容适合谁?一种是刚接触微电网调度的研究生,想弄明白随机优化到底怎么建模、怎么从随机变量走到可求解的数学规划;另一种是已经在做综合能源系统、虚拟电厂或者园区供能设计的工程师,手里有实际负荷曲线和光伏出力数据,想给自己的项目加一层“抗不确定性”的调度逻辑。下面我就按我从建模到写代码、再从调参到复盘的全过程,把这套东西一层层剥开讲。
1. 先搞清楚源荷随机特征到底“随机”在哪
1.1 热电联供微网的基本构成
热电联供微网和普通微网最大的区别在于,它同时承担电力和热力的供应,常见的设备组合包括热电联产机组(CHP)、燃气锅炉、电锅炉、蓄电池、储热罐,再加上光伏和风电这类可再生能源,以及电网交互关口。CHP机组是这张网的核心,它烧气同时产电和余热,热电比这个参数把电出力和热出力捆在一起,所以电和热两套平衡不能各自独立做。
传统调度模型里会把CHP简单当成一个“出力可调、效率固定”的电源,但实际上它的电出力和热出力有强耦合区间,爬坡约束、最小出力约束、启停约束都比普通火电机组复杂。再加上热负荷本身有惯性、储热罐能缓冲,热电联供微网的调度问题天然是一个多能互补、时序耦合的优化问题。
1.2 电源侧的随机性来自哪里
电源侧随机性主要来自两部分:光伏和风电。
光伏出力本质上跟辐照度强相关,而辐照度在一天内服从明显的波动规律。工程上常用 Beta 分布描述晴空指数的概率特性,或者直接用历史出力数据做统计,拟合出预测误差分布。做日前调度的时候,你手里只有光伏功率的预测曲线,实际出力会在预测值附近波动,这个波动幅度在不同季节、不同天气下差异很大。所以不能简单设定一个固定误差带,更好的做法是把预测误差建模成正态分布,再结合天气类型修正标准差。
风电的随机性更直接,风速服从 Weibull 分布,风机出力又是风速的分段非线性函数。很多简化模型直接把风电出力当作一个带误差的随机变量,然后用场景集去覆盖可能的出力区间。我做了几个算例之后的感觉是:风电的随机特性比光伏更难“驯服”,因为它可能在一个调度时段内大幅摆动,对系统备用容量和爬坡能力的要求更高。
1.3 负荷侧的随机性同样不能轻视
电负荷和热负荷都有明显的时序性,但随机特征不太一样。电负荷波动快,受用户行为影响大,尤其是工商业混合的园区,尖峰时段的不确定性非常明显。热负荷则受天气温度影响更大,热惯性使得它在短时间内的波动相对平滑,但一旦遇到温度骤降,热负荷预测偏差会突然放大。
如果忽视负荷随机性,模型很容易出现“平衡约束过紧”的问题。比如确定性模型把某一时段的电负荷定死为 1000kW,实际波动到 1080kW,系统如果没有留备用,就只能切负荷或者从电网高价购电。把负荷随机特征纳入优化,本质上是给调度决策留出“安全裕度”,但裕度留多大又会影响经济性,这就是随机优化的核心矛盾。
2. 从随机变量到可计算的场景集合
2.1 用场景法还是鲁棒优化
描述随机特征大体有三条路:随机规划场景法、鲁棒优化区间法、分布鲁棒优化。这个项目采用的是场景法,原因很朴素:场景法能给出具体的概率分布信息,能算出期望成本,而且它的结果可以直观地翻译成“在某种出力情况下我该怎么调整机组出力”。鲁棒优化虽然更保守、能保证极端情况下不越限,但它把所有不确定性都压进一个区间盒子里,结果往往过于保守,经济性差很多。
场景法的思路是:把随机变量的概率分布离散成有限个代表性子场景,每个场景附带一个发生概率,然后把所有场景同时代入优化模型求解。这样,目标函数就变成了各个场景下成本的加权期望值,约束里的功率平衡量也会带上场景下标。
2.2 场景生成怎么做
场景生成这一步直接决定后面优化的质量。常见的做法有蒙特卡洛抽样和拉丁超立方抽样。蒙特卡洛简单直接,但需要的样本量很大,场景一多优化规模就爆炸。拉丁超立方抽样则先把每个随机变量的累积概率分布分成等间隔区间,再从每个区间里取样本,能用少得多的场景覆盖相近的概率空间,这个项目里我用的是后者。
% 生成光伏出力场景(简化示意) n_scenarios = 200; t_horizon = 24; % pv_forecast 是日前预测曲线,sigma_pv 是预测误差标准差序列 pv_scenarios = zeros(t_horizon, n_scenarios); for t = 1:t_horizon error_t = normrnd(0, sigma_pv(t), 1, n_scenarios); pv_scenarios(t, :) = max(0, pv_forecast(t) + error_t); end这段代码有个容易踩的坑:直接用相加的方式可能让场景出现负值和超过装机容量的数值,所以必须做截断,把出力压回 [0, 装机容量] 区间。另外,各时段的误差之间不应该完全不相关,更精细的做法是加入时间相关性,比如用一阶自回归模型让相邻时段的误差互相牵动,这样场景才更贴近真实波动曲线。
2.3 场景削减:为了能求解
200个场景已经能让非线性模型卡死了,何况热电联供微网模型本身带整数变量,场景一多就是指数级爆炸。所以必须做场景削减。工程上最常用的是基于概率距离的快速前向选择法和同步回代消除法,思路是不断剔除那些与其他场景“距离最近”的场景,把被剔除场景的概率累加到保留场景上。
场景削减的度也要掌握好。削减到太少,比如只剩 5 个,随机特征就被磨平了,优化结果跟确定性模型差别不大;削减到太多,比如 50 个以上,求解时间又会让人崩溃。就24小时调度问题来说,我实践下来留 20 到 30 个场景是比较舒服的区间,既能保留随机分布的形状,求解时间也还能控制在分钟级别。
3. 优化模型的框架:目标、约束和两阶段结构
3.1 目标函数:先把要算的钱列清楚
热电联供微网优化的目标函数一般是综合成本最小化,主要有这么几块:
| 成本项 | 说明 | 备注 |
|---|---|---|
| 购电费用 | 从上级电网买电,按分时电价计算 | 峰谷电价差大时,模型会优先用储能套利 |
| 购气费用 | CHP和燃气锅炉的燃料成本 | 热负荷权重大的园区,这是大头 |
| 运维成本 | 各设备单位出力的维护费用 | 光伏风电运维成本低,CHP高一些 |
| 启停成本 | CHP启停的热应力损耗折算 | 带整数变量,建模别漏 |
| 惩罚成本 | 失负荷惩罚和弃风弃光惩罚 | 给模型一个“软退路”,避免无解 |
如果项目还要求考虑环保效益,可以把碳排放成本也加进去,也就是在气耗成本后叠加二氧化碳排放因子。这样模型在比较“多用电锅炉产热”和“多用燃气锅炉产热”时,会自动权衡电价、气价和碳排放价格之间的相对大小。
3.2 两阶段随机规划的结构
把随机场景引入之后,调度决策分成两个层次:第一阶段决策是日前制定的、不依赖具体场景的“现在就要定”的变量,比如机组启停状态、是否与电网交换电力这些需要提前宣告的动作;第二阶段决策是可以在场景发生后调整的变量,比如各机组的具体出力、储能充放电功率、切负荷量。
为什么要这么分?因为实际运行里,调度员不可能等到光伏出力精确知道了再决定CHP启停,机组开机是提前定的,但开机后发多少电可以实时调整。两阶段结构正好对应这种“事前决策+事后调整”的逻辑,这也是它比确定性模型更贴近工程实际的原因。Matlab 里用 YALMIP 建模时,第一阶段变量和第二阶段变量都可以直接用 sdpvar 声明,差异只体现在约束条件和目标函数的下标组织上。
3.3 关键约束与建模中的坑
功率平衡约束要分电和热两条线分别写。电功率平衡式是:电网购入功率加光伏出力加热电联供机组电出力加蓄电池放电,等于电负荷加蓄电池充电加电锅炉耗电。热功率平衡式是:CHP余热加燃气锅炉产热加储热罐放热,等于热负荷加储热罐蓄热。
CHP建模是整篇代码里最容易出错的地方。它的可运行区域不是一个简单矩形,而是由电出力上下限、热出力上下限、热电比上下限共同围成的多边形。常见做法是用一组线性不等式去逼近这个运行区域,比如:
% CHP运行区域约束(示意,数值需按实际机组参数标定) H_chp_min <= H_chp <= H_chp_max; P_chp_min <= P_chp <= P_chp_max; H_chp <= k_upper * P_chp; % 热电比上限 H_chp >= k_lower * P_chp; % 热电比下限实际运行中还有个容易忽略的点:CHP刚启动时热出力爬坡跟不上电出力爬坡,这会破坏瞬时热平衡。如果模型里没加供热爬坡约束,求解器会给出一个理论上可行、实际上不可能执行的出力轨迹。所以建议在CHP热出力上也加一段爬坡约束,数值可以比电出力爬坡限制更严格。
储能约束也有讲究。蓄电池的荷电状态递推式是 SOC(t+1) = SOC(t) + Pch·ηch·Δt − Pdis/ηdis·Δt,但这里必须有充放电互斥约束,否则求解器会在同一时段既充电又放电,白白损耗能量,算出来的成本还显得挺“划算”。互斥约束一般用两个二元变量和一个大M量来处理。大M的取值不能拍脑袋,取小了会错误地砍掉可行域,取大了又会影响求解器数值稳定性,我习惯取设备最大功率的3到5倍。
燃气锅炉相对简单,只有出力上下限和爬坡约束。储热罐跟蓄电池一样有放热和蓄热两个方向,但它的自损耗系数通常比电池小,而且热负荷的短时波动用储热罐平抑的效果非常好。
3.4 旋转备用的随机化处理
考虑源荷随机特征后,光靠功率平衡约束还不够,系统还得留旋转备用来应对“预测场景之外”的突发偏差。项目里常见做法是加一条备用约束:常规机组可调出力加上储能可放电容量,要大于当前时段所有场景下新能源出力的最大预测偏差再加上负荷的最大预测偏差。
这条约束的实质,是把“随机特征”从目标函数里的概率期望,扩展到约束条件里的概率覆盖。两阶段优化与备用约束配合,模型才能既考虑经济期望,又保证安全底线。我在测试中发现,备用约束设太紧会让求解时间暴涨,建议先把备用系数设成 5% 到 10% 跑通流程,再逐步收紧观察成本变化趋势。
4. Matlab 代码实现的完整流程
4.1 整体程序结构
这个项目落到代码层面,我习惯按五个模块组织程序:数据输入模块、场景生成模块、优化建模模块、求解配置模块、结果后处理模块。数据输入模块把设备参数、分时电价、气价、负荷预测曲线、光伏风电预测曲线单独放一个结构体里,不要散落在脚本各处,否则后面参数标定能找死人。
场景生成模块按上一节的方法生成并削减场景。优化建模模块是核心,用 YALMIP 声明变量和约束。求解配置模块里设置求解器参数,比如精度、时间上限、MIP gap。结果后处理模块负责把每个时段的电功率平衡、热功率平衡、机组出力、储能SOC画出来,顺便输出各场景下的成本分布直方图。
4.2 用 YALMIP 搭模型的基本骨架
YALMIP 是个好东西,能把 Matlab 的优化建模效率提高一个量级。先用 sdpvar 声明决策变量,再把约束和目标函数写成 YALMIP 表达式,最后调用 optimize 求解。骨架大概长这样:
%% 声明变量(以24时段、20个场景为例) T = 24; S = 20; P_chp = sdpvar(T, S, 'full'); % CHP电出力 H_chp = sdpvar(T, S, 'full'); % CHP热出力 P_gb = sdpvar(T, S, 'full'); % 燃气锅炉热出力 P_buy = sdpvar(T, S, 'full'); % 从电网购电 SOC = sdpvar(T+1, S, 'full'); % 蓄电池SOC u_chp = binvar(T, 1, 'full'); % CHP启停状态,第一阶段变量 %% 约束集合 Constraints = []; % ... 在此处加入功率平衡、CHP区域、储能递推、备用约束等 %% 目标函数 Objective = ...; %% 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, ... 'savesolveroutput', 1); result = optimize(Constraints, Objective, ops);这里有个细节:第一阶段变量 u_chp 只有时间维度 T,没有场景维度 S;第二阶段变量 P_chp、H_chp 则带 T×S 两维。这样写才能体现“启停决策在知道随机结果之前就定了下来”的非预期性约束。
4.3 求解器选择和参数经验
YALMIP 底层可以接多个求解器,常见的就两类:一类是内点法为主的连续优化求解器,另一类是分支定界法为主的混合整数求解器。这个项目因为带启停变量和充放电互斥变量,属于混合整数线性规划(MILP),所以必须选支持整数规划的求解器。学术圈常用 Gurobi 或 CPLEX,它们的求解速度确实比开源求解器快很多,尤其场景数上来之后,差距是数量级的。
我们自己测试过,几十个场景的 MILP 模型,用默认参数可能要跑十几分钟,但把 MIP gap 容忍度从默认值放宽到 1%,求解时间能压缩到两三分钟,而且成本差距基本在千分之一以内。调度问题又不需要站点收敛到小数点后八位,所以实战中我都是直接设'mipgap', 0.01。
如果手头没有商业求解器授权,Matlab 自带的 intlinprog 也可以凑合,但它的符号运算没有 YALMIP 方便,需要把约束手动整理成矩阵形式,代码量大不说,还容易写错系数。所以我还是推荐 YALMIP 建模、商业求解器求解的组合。
4.4 结果的输出和分析
求解完成后,第一件事不是写论文,而是做平衡校验。把每个时段的电功率平衡和热功率平衡都算一遍,看残差是否为零。如果残差有数值噪声,大概率不是模型问题,而是求解精度设置太低,调一下求解器的可行容差就行。
校验通过后,输出关键图:电功率平衡堆叠图、热功率平衡堆叠图、CHP机组电热出力轨迹、蓄电池和储热罐的 SOC 曲线、各场景总成本直方图。这些图放在一起,评审一眼就能看出模型是否合理:如果热功率平衡堆叠图里燃气锅炉出力整天满负荷,而储热罐几乎不用,那说明经济参数设置有问题,电锅炉或者CHP的产热成本应该再调一调。
5. 调试经验与常见问题排查
5.1 模型无解时先查什么
模型无解是这个项目里最高频的问题,没有之一。我自己的排查顺序是:先看功率平衡约束的符号有没有写反,再看所有设备出力上限是否比下限大,然后看储能递推的初值 SOC(1) 是否设了,最后检查有没有漏掉“储热罐在同一时段不能既蓄热又放热”这个约束。
还有一个非常隐蔽的坑:CHP的热电比上下界约束写得过紧,加上爬坡约束后,在热负荷低谷时模型找不到解。这是因为CHP为了满足电负荷必须最低出电,但对应的热出力下界已经超过了热负荷需求。解决办法有两种,一是调低CHP的最小电出力,二是加一个散热器(或电锅炉)作为热负荷的“安全阀”,允许多余的热量被消纳掉。
5.2 求解时间太长怎么压
场景数上来之后,整数变量的总数会线性增长,求解难度却是超线性增长。压缩求解时间的路径优先级我排一下:先削减场景数,再看备用约束能不能放松,再检查储能互斥约束是不是用了太多大M,最后才是调求解器的 MIP gap。如果这些都不够,可以考虑把部分决策变量从单时段自由度改成分段常数,比如把热负荷分成峰、平、谷三段分别设平抑目标,变量维数会大幅下降。
另外一个常被忽略的技巧是给二元变量设初始解。用确定性模型的启停结果作为热启动值传给随机优化,YALMIP 里可以用assign和optimize(...,'usex0',1)实现。这一步操作简单,但求解速度往往能快 30% 以上,尤其在场景数多、机组数多的模型里,效果立竿见影。
5.3 结果不合理时如何定位
如果优化结果出现“某个时段电网购电量极高但CHP停机”的现象,多半不是模型错,而是电价和气价的比值设置不合理,导致经济上确实买电比烧气划算。这时候要回到数据校核,而不是改模型。如果结果出现“储能整天充满放净好多次”的锯齿状曲线,那就是没有给储能加循环次数限制或者寿命损耗成本,模型在疯狂套利,实际工程里电池会很快报废。处理办法是在目标函数里加一项储能充放电的累积能量惩罚,数值虽然小,但足以让模型停止“高频低买高卖”的投机行为。
我还有一个习惯:每次调参之后,把核心输出指标存到一个 Excel 或者 CSV 里,包括总成本、购电购气比例、CHP利用小时数、弃风弃光量。跑上几十轮对比,哪个参数影响大、哪个参数不敏感,一目了然。这个方法帮我在项目验收阶段省了很多解释成本。
6. 从期望成本到风险管控的扩展
6.1 期望成本之外的尾部风险
纯期望值目标有一个天然短板:它可能挑中“平均来看很好,但极端场景下会巨亏”的方案。比如某个调度策略在95%的场景里成本很低,但在5%的极端场景里切负荷惩罚高得离谱。期望值一平均,看起来挺美,但实际运行碰到那次极端情况就麻烦了。
这个项目如果想进一步,可以在目标函数里加 CVaR(条件风险价值)项,也就是对最差的那部分场景求条件期望,把它也纳入优化。这样决策者可以通过调 CVaR 权重,在“平均经济最优”和“极端风险可控”之间滑杆,这也正是现实调度里运营方最关心的取舍。如果不加这项,模型在随机特征很强、误差分布尾部分散的场景下,给出的调度策略可能偏冒险。
6.2 滚动时域和日内修正
日前模型算完,还有一步值得做:把调度周期缩短成滚动时域。具体做法是每次只执行前几个时段的指令,到下一个时段用最新的光伏、负荷实测数据重新生成场景、重新求解。这本质上是用“新信息”不断收缩随机性的范围,比单纯把日前优化结果一口气执行到底要稳得多。工程落地的时候,很多项目就是靠这种“日前计划+日内滚动修正”的框架把随机优化真正用起来的。
6.3 往多微网协同扩展的思考
单微网做完随机优化之后,自然想到的问题是:多微网之间能不能通过共享备用容量来对冲各自的随机波动?不同微网的光伏出力曲线和负荷曲线在时间上不完全同相位,A微网的光伏低谷时段,B微网的热负荷高峰还没来,如果有一条联络线,两个微网就可以互补余缺,各自需要的备用容量都小一些。这也是这类研究往后延伸最有价值的方向之一。
我个人在实际操作中的体会是:随机优化最难的从来不是求解器,而是对不确定性的建模分寸。场景生成太粗,结果不可信;场景生成太细,模型算不动。这个度需要反复试。建议你先拿一个小系统、少场景数把整个流程跑通,再逐步加压,不要一上来就奔着几百个场景去。另外,所有随机参数的标准差和分布类型,尽量从实际数据里去拟合,别拍脑袋设值,否则模型再精巧,输入数据失真,输出也只能是自欺欺人的结果。这套代码框架跑顺之后,换设备参数、换场景生成方式、换风险度量指标都比较容易,算是性价比很高的一个项目起点。