1. 项目概述与整体设计思路
1.1 核心需求解析
前阵子帮一个做综合能源方向的同学看代码,发现他卡在“源荷随机性”这个概念上很久。他手里有一份热电厂出力的确定性调度代码,跑起来没问题,但只要一涉及风光预测误差、负荷波动,整个模型就不知道该怎么改。这个项目标题“考虑源荷随机特征的热电联供微网优化研究(Matlab代码实现)”其实点得很明白:不是单纯做个微网调度,而是要把不确定性“请进”优化模型里,用随机优化的思路去处理源侧和荷侧的波动。
先交代一下背景。热电联供(Combined Heat and Power,CHP)微网是目前园区级能源系统的主流形态,核心设备是燃气轮机或内燃机,它在发电的同时回收余热供给热负荷,能效比单纯的电热分产高出不少。微网里还会配风电、光伏这类可再生能源,以及电储能、蓄热罐、燃气锅炉等辅助设备。整个系统的运行目标是在满足电负荷和热负荷的前提下,让总运行成本最低。
但问题来了——风电和光伏出力是随机波动的,负荷也有预测误差,如果优化调度时把这些都当成固定值,那调度方案在实际执行时可能严重偏离预期,极端情况下会导致切负荷或者弃风弃光。所以就有了“考虑源荷随机特征”这个研究方向。
这个项目适合谁看?两类人。一类是正在做微电网/综合能源系统方向毕业论文的硕士生,另一类是刚接触随机优化、想知道怎么用Matlab落地实现的工程师。我会把随机建模的原理、场景生成与削减、优化模型的数学表达、以及Yalmip+Cplex的完整实现逻辑全部讲透,最后附上我在调试过程中踩过的坑。
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 Reduction,SBR)。核心思想很简单:每次找出一对距离最近的场景,删掉其中一个,同时把被删场景的概率累加到保留场景上,让总概率分布尽量不变。
具体算法流程是这样:
- 计算所有场景对之间的概率距离,常用的是欧氏距离乘以场景概率的加权值
- 找到距离最小的一对场景(假设场景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.2 | 0.5 MW/h |
| 燃气锅炉 | 5 MW | 热效率85% | 0.8 MW/h |
| 电储能 | 2 MWh | 充放电效率95%,SOC范围0.1-0.9 | 0.5 MW |
| 蓄热罐 | 4 MWh | 储放热效率90%,容量范围0.15-0.95 | 0.6 MW |
| 风电 | 2 MW | - | - |
| 光伏 | 1 MW | - | - |
3.2 目标函数:最低运行成本到底在优化什么
优化目标是最小化系统在一个调度周期(通常是24小时)内的总运行成本。总成本的组成项比较多,但每一项都有明确的工程意义,我列个清单:
第一项是购电成本(或售电收益)。微网与上级电网存在功率交换,用分时电价结算。峰时段电价高,谷时段电价低,这直接影响储能的充放电策略。第二项是燃料成本。燃气轮机和燃气锅炉都烧天然气,成本按单位热值价格乘以消耗量计算,而消耗量与设备的出力存在明确的效率换算关系。第三项是运维成本。我按设备出力的一定比例折算,比如燃气轮机的运维费率取0.05元/kWh,风光的运维费率低一些,取0.01元/kWh。第四项是弃风弃光惩罚。当系统无法消纳全部可再生能源时,会以惩罚系数的形式计入成本,这样优化器会尽量避免弃风弃光。最后一项是切负荷惩罚。如果系统功率不足导致必须削减电负荷或热负荷,会产生很高的惩罚成本,我把系数设置成购电价的几十倍,确保优化器只在极端情况下才舍得切负荷。
把这几项加总写成一个目标函数,就得到了我们的优化目标。特别提醒一下:在第二阶段的期望成本计算里,每个场景都会算一次从上级电网购电、弃风弃光、切负荷的成本,再按场景概率加权求和。这种“第一阶段决策 + 第二阶段期望调整”的建模方式正是随机规划中典型的两阶段随机优化的结构。
3.3 约束条件的体系与细节
约束条件是优化模型的主体,我拆成几类来说。
第一类是功率平衡约束。对于每个时段和每个场景,系统的电功率必须满足平衡关系:风电出力加光伏出力加CHP发电加电储能放电加购电,等于电负荷加电储能充电加售电。热功率同理:CHP产热加燃气锅炉产热加蓄热罐放热,等于热负荷加蓄热罐充热。平衡约束是所有调度问题的骨架,表达式本身不复杂,但容易在编写代码时搞混正负号的方向。
第二类是设备出力上下限约束。每台设备都有最小技术出力和最大技术出力限制。特别要说的是CHP机组还有最小开机电出力要求,低于这个值时机组无法稳定运行,反映出燃气轮机不能频繁启停、不能在极低负荷下运行的工程实际。
第三类是爬坡约束。燃气轮机和燃气锅炉在相邻时段的出力变化不能超过规定速率。这一点在随机优化里尤其重要,因为场景中可能出现相邻时段风电出力大幅跳变的情况,如果没有爬坡约束,调度方案在物理上根本执行不了。
第四类是储能系统约束。电储能的荷电状态(State of Charge,SOC)在每个时段都有动态递推关系,并且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个场景) |
|---|---|---|---|
| 10 | 18420 | 3.2 | 2.10% |
| 20 | 18135 | 7.8 | 0.52% |
| 30 | 18072 | 13.5 | 0.17% |
| 50 | 18041 | 28.6 | 基准 |
从表格结果来看,场景数从10增加到20带来的精度提升明显,求解时间也在可接受范围内;20到30的边际收益开始递减;到50个场景时结果趋于稳定但耗时翻倍。综合精度和效率,20个削减场景是这个测试系统的甜点值。你的系统如果负荷曲线波动更剧烈,可以把场景数适当调到30以获取更多稳健性。
6. 常见问题与排查技巧实录
6.1 Yalmip报错与求解器配置问题
Q1:Yalmip报“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确定性优化”对照实验时,更要在同一组场景下回测,不然出来的差异无法归因。
最后,关于项目扩展,我建议对两阶段随机优化熟悉之后,可以考虑加入鲁棒优化方法做对比。场景法需要假设概率分布已知,如果你手里的历史数据不足,分布假设就不太可靠。此时用盒式不确定性集合驱动的鲁棒优化会更有优势。两种方法各有侧重,结合起来看问题会更全面。