源荷随机特征下热电联供微网优化调度与Matlab实现

项目标题里最核心的四个词“源荷随机特征”“热电联供微网”“优化”“Matlab实现”分开看都不算冷门但组合在一起就是典型的综合能源系统方向的研究课题了。这两年做微电网优化的朋友应该都能感觉到纯电微网的套路已经非常成熟真正容易出问题也容易出成果的恰恰是这种带热负荷、带热电联供机组、还要把风电光伏和负荷随机性一起塞进优化模型的场景。这篇文章我就围绕这个题目把从模型拆解到Matlab代码落地的整个思路捋一遍也会把我自己调试过程中踩过的坑一并写出来给正在做相关课题或者准备复现类似工作的同学一个参考。这个研究到底解决什么问题一句话来说在风电、光伏出力以及电/热负荷都存在不确定性的情况下怎么安排热电联供机组、燃气锅炉、储能等设备的出力计划让系统在满足供电和供热需求的同时运行成本最低、排放尽量小而且还能扛得住源荷波动。适合正在做微电网优化调度、综合能源系统规划或运行优化的同学看也适合刚接手MatlabYalmip这类工具、想快速跑通一个完整算例的入门者。1. 项目解读热电联供微网优化到底在做什么1.1 源荷随机特征这句话怎么拆解“源”指的是风电机组和光伏阵列这种出力受天气影响的分布式电源出力曲线不是人为可控的来一阵云、起一阵风出力就变了。“荷”指的是用户侧的电负荷和热负荷居民用电有早晚高峰工厂用电跟着生产班次走冬天的热负荷还跟室外温度强相关。这些源和荷都有一个共同特点用确定性曲线来描述它们是不充分的实际运行中它们一直在随机波动。做热电联供微网优化研究时如果无视这种随机性把光伏出力设成一个固定的预测曲线然后对着这个曲线去做优化最后得到的结果通常是非常乐观的。遇到实际运行时风光出力比预测低或者负荷比预测高原先的调度方案就会出现功率缺额供热温度也可能掉出允许区间。所以课题名里专门强调“考虑源荷随机特征”本质上就是要让优化模型在面对不确定性时依然有足够的鲁棒性和经济性。随机特征具体怎么刻画一般有两层意思。第一层是概率分布比如风速近似服从Weibull分布太阳辐照度可以用Beta分布描述负荷预测误差通常假设为正态分布。第二层是时间相关性源荷的预测误差在相邻时段之间不是相互独立的而是有自相关的特性。这两层信息最终都可以转化成优化模型中的随机场景用多个场景来代表可能发生的源荷情况。1.2 为什么需要联合考虑源和电负荷、热负荷的随机性很多初接触这个课题的朋友会有一个疑问既然热负荷相对平稳供热需求的变化也比电负荷缓慢那能不能只考虑电负荷和风光出力的随机性热负荷直接用确定值我在实际做课题时发现这种简化在纯供热系统里问题不大但在热电联供微网里会带来明显偏差。原因在于热电联供机组存在“以热定电”的运行约束也就是机组的发电功率和供热功率是耦合在一起的。热负荷一旦出现波动机组的电出力上下限就会跟着变电功率平衡也会受影响。举个直观的例子一个严寒天的早晨热负荷预测值假设是10MW但实际因为气温骤降变成了11.5MW。机组为了多供1.5MW热量其电出力就会被抬高或者压低具体要看机组运行区间和运行模式。如果优化模型里没有把这种热负荷的随机波动考虑进去系统中其他电源的计划就可能完全偏离实际需求储能充放电计划也会错位。因此至少要把电、热两类负荷以及风光出力都纳入随机场景的构建才算真正做到了“源荷随机特征”。联合考虑还有一个实际收益就是可以利用电热负荷之间的时序互补性。比如夜间风电出力大、电负荷低这时候可以通过电锅炉或者热泵把多余电能转化为热能存储起来白天再放热。但如果电负荷和热负荷的随机性没有被描述在同一个框架里这种协调优化的机会就很难被模型捕捉到做出来的调度方案自然不够经济。2. 系统建模与随机性处理2.1 热电联供微网的典型结构与设备模型要搭建优化模型首先要把微网里有哪些设备、怎么连接、能量怎么流动用数学语言描述出来。一个典型的热电联供微网包括以下模块热电联产机组CHP这是系统的核心供能设备可以同时输出电功率和热功率。它的技术特点是热电出力比在一定范围内可调通常用多边形的可行运行域来描述。简化建模时也可以用固定热电比但精度稍差。燃气锅炉作为补充热源在CHP供热不足时启动特点是响应快、热效率高但消耗燃气。风电机组和光伏阵列可再生能源电源出力具有随机性通常作为负的负荷或者可控程度有限的电源。预测值取自气象数据随机性在后面用场景描述。电储能电池和蓄热罐实现电、热两种能量在时间维度上的平移。电池有充放电功率限制和容量约束蓄热罐有储/放热功率限制和热容量上限。电负荷和热负荷用户需求分为预测值和实际值。预测值用于场景生成实际值用于仿真验证。设备模型的核心是输入-输出关系。以CHP机组为例常见建模方式是[ C_{chp,t} a \cdot P_{chp,t} b \cdot H_{chp,t} c ]其中(C_{chp,t})是燃料成本(P_{chp,t})是电功率(H_{chp,t})是热功率(a)、(b)、(c)是效率系数。这个二次函数在求解时往往线性化成分段线性函数。电池模型需要考虑充放电效率蓄热罐模型需要考虑热损失系数。把这些数学模型写清楚优化模型的基础才算打牢。2.2 风光出力与负荷的随机特征建模方法源荷随机特征的建模本质上是把不确定性从“不可描述”变成“可以用场景描述”。我在实际项目中用过的比较可行的方法有这么几种第一种是参数化概率分布法。风速用两参数Weibull分布拟合太阳辐照度用Beta分布拟合负荷预测误差假设为均值为0的正态分布。这种方法的优点是简单缺点是无法表达复杂的时空相关性。第二种是场景生成法。基于历史数据通过蒙特卡洛抽样或者拉丁超立方抽样生成大量包含时序相关性的源荷场景。注意这里不能直接按独立分布抽样否则生成的负荷曲线会非常毛糙和真实负荷曲线的平滑特性不符。一般做法是先对预测误差的时间序列用ARMA模型或者ARIMA模型建模再从这个模型里抽样得到误差场景叠加到预测曲线上。第三种是数据驱动的场景生成比如用生成对抗网络去拟合源荷联合分布。这种方法适合数据量比较大的场合我在论文复现里用得不多因为代码复杂度和训练成本都比较高对课题组没有专门算力的情况不太友好。实际工程中最稳的组合是“预测曲线 ARMA误差抽样 拉丁超立方采样”。既能保证时序相关性又能用较少样本覆盖大概率空间后续的场景削减也方便。2.3 场景生成与削减让随机问题变成可计算问题有了随机场景之后优化问题就变成多场景下的随机优化问题。理论上场景数越大不确定性描述得越精准但计算量也越大。一个24小时调度问题如果场景数取1000每个场景有24个时段决策变量就会膨胀到几十万个尤其加入二进制变量之后求解难度会急剧上升。所以实践中几乎都会做场景削减。比较常用的削减方法是快速前向选择法或者同步回代消除法。核心思想是在保留原始场景集合概率特征的前提下删除那些与其余场景距离较近、代表性弱的场景把1000个场景削减到50个或者20个。每个保留场景赋予一个概率权重所有保留场景的权重之和等于1。用Matlab实现场景削减时我偏向直接用现成的函数。如果不想依赖额外工具箱可以自己写一个简单的贪心削减每次找一对欧式距离最小的场景删除其中一个把被删场景的概率加到保留场景上。虽然不如同步回代消除法那么精确但对大多数课题研究来说已经足够了。场景削减后的一个细节是削减后的集合所对应的期望成本会和原始场景集合的期望成本有偏差。做结果分析时要留意这个偏差的规模一般控制在5%以内就可以接受。如果偏差太大说明初始场景数不够或者削减比例过高需要回调参数。3. 优化模型与求解策略3.1 目标函数经济性、环保性与供能可靠性的取舍热电联供微网优化的目标函数通常不是单一的。最经典的是运行成本最小成本项包括燃气成本CHP机组和燃气锅炉消耗的燃料费用。购电成本微网与配电网交换功率的费用。如果允许向配电网卖电还要加入售电收益项。设备启停成本CHP机组和锅炉的启动成本这个虽然单次不大但对机组频繁启停有一定抑制作用。环保性作为第二个目标时可以引入碳排放成本将CO2排放量按碳价折算到目标函数中。也可以用多目标优化的形式同时优化经济性和排放量但多目标问题的求解复杂度比单目标高不少。研究性课题里把排放折算成经济成本放入单目标函数处理起来简洁得多。供能可靠性作为第三个目标时最常见的处理方式是引入失负荷惩罚。在随机优化框架下如果一个场景下实际发生的负荷超过了可用出力那么缺额部分就会被切掉切负荷量乘以惩罚系数进入目标函数。这个系数如果设得太低优化结果会保持很小的备用经济性很好但风险大设得太高备用过大成本上升。我自己的经验是惩罚系数至少取单位购电成本的10到20倍才能有效抑制切负荷的发生。3.2 约束条件功率平衡、爬坡、储能、热负荷平衡约束条件是整个模型里最容易出bug的地方。我列一下必须覆盖的几类并说明它们各自在模型中起什么作用。电功率平衡约束要求每个时段内所有电源出力加上储能放电、加上从配电网购入的电功率等于电负荷加储能充电加向配电网售出的功率。如果模型允许切负荷右侧还需要加上切负荷变量。这条约束如果写得不对求解出的结果要么明显违反物理规律要么干脆无解。热功率平衡约束类似但要注意的是CHP机组的高温出水经过换热站后供给热负荷蓄热罐可以吸收或释放热量。蓄热罐的蓄热状态量还必须满足时序递推方程并受容量的上下限约束。机组爬坡约束描述的是CHP机组在相邻时段之间电出力的变化速率限制。如果不加这条约束优化模型可能会给出一个锯齿状出力曲线现实中机组完全跟不上。同类约束也适用于燃气锅炉的出力变化。储能约束则包括充放电功率上限、容量上下限以及同一时段不能同时充放电的逻辑约束这个逻辑约束通常用二进制变量实现。还有一个容易被忽略的约束是CHP运行区间约束。实际机组的电出力和热出力不是任意组合都可以的存在一个可行域需要用一组线性不等式去逼近。不少初学的同学只在模型里写了电出力和热出力的上下限结果得到的调度点落在了机组的不可运行区域仿真时才发现不对。3.3 求解思路混合整数线性规划与场景法的结合优化模型经过线性化处理后因为里面含有二进制变量所以是一个混合整数线性规划MILP问题。场景法把随机性问题转化为一个大规模确定性MILP问题然后用成熟的商业求解器求解。在Matlab环境里用Yalmip作为建模工具底层调Cplex、Gurobi或者开源的Cbc都可以。我个人最常用的是YalmipGurobi的组合建模代码简洁求解效率也高。Gurobi对MILP问题的剪枝效果在同类求解器里数一数二对场景规模较大的算例能明显感觉到时间优势。没有Gurobi许可证的话Cplex也可以Cbc免费但求解速度要慢不少适合模型比较小的验证场景。求解思路的细节是把所有场景下的决策变量统一放进一个大的优化问题里场景之间的变量通过“非预期性约束”耦合。最简单的耦合方式是让第一个时段的决策变量在所有场景下保持一致因为实际运行时在第一个时段还没有获取未来信息调度决策不能因为场景不同而变化。如果忽略这一点就会得到“事前知道未来”的非预期性解结果过于乐观。这一点在做随机优化研究时非常关键。4. Matlab代码实现要点4.1 代码整体结构与数据准备拿到这个课题之后我建议先把代码的结构想清楚不要一上来就写优化模型。一个可维护的Matlab代码通常分成四个部分数据输入、场景生成、模型构建、求解与结果输出。数据输入部分需要准备的基础数据包括24小时的电负荷预测曲线和热负荷预测曲线风速和光照强度的历史序列或预测序列CHP机组、燃气锅炉、电池、蓄热罐的技术参数天然气价格、购电电价等经济参数。这些参数用结构体或者表格存起来就行我用得比较习惯的是把所有参数放在一个结构体参数表params里后续模型里引用起来非常方便。场景生成部分把源荷不确定性转化成场景集合。这里要留意的是用来做随机优化的风电、光伏和负荷场景需要来自同一个随机源不能各编各的。比如用ARMA模型生成误差序列时得把风、光、电负荷、热负荷的误差序列一起模拟保留它们之间的相关性。一种简便做法是假设各随机变量之间的相关系数矩阵然后用经Cholesky分解处理后的随机数序列来抽样。模型构建部分Yalmip的语法我个人觉得是能提高效率的。每个变量用sdpvar定义二进制变量额外指定类型约束用循环填入一个cell数组目标函数直接写表达式最后调用optimize。代码可读性比直接用Cplex的底层API要高很多。求解与结果输出部分除了打印出目标函数值和各设备出力还应把电平衡、热平衡的残差输出来方便定位约束写错的地方。仿真结果如果画成图建议把预测场景和实际场景画在同一条时间轴上能直观看出随机优化的价值。4.2 核心函数与求解器调用Yalmip/Gurobi举例下面给一个Yalmip建模的核心代码片段覆盖CHP和储能的基本约束。这个是我调试过很多次的模板写给大家做个参考。% T为时段数Ns为场景数x可能包含三维变量时段、场景、设备 P_chp sdpvar(T, Ns, full); % CHP电出力 H_chp sdpvar(T, Ns, full); % CHP热出力 P_bat sdpvar(T, Ns, full); % 电池出力正为放电负为充电 S_bat sdpvar(T1, Ns, full); % 电池荷电状态 u_chp binvar(T, Ns); % CHP启停状态 Constraints []; % CHP运行区间约束简化 Constraints [Constraints, P_chp 5000*u_chp, P_chp 0]; Constraints [Constraints, H_chp 3000*u_chp, H_chp 0]; % 电池充放电约束 Constraints [Constraints, S_bat(:,11) S_bat(:,1) ...]; % 实际写的时候按“上一时刻SOC 充电功率*效率 - 放电功率/效率”递推 Constraints [Constraints, -200 P_bat 200]; Constraints [Constraints, 0 S_bat 1000]; % 电功率平衡 Constraints [Constraints, P_chp P_wind_scene - P_load_scene P_bat P_grid 0]; % 目标函数 Objective sum(sum(C_gas * ... )) sum(sum(C_grid .* P_grid));这段代码虽然简化核心意思很清楚所有设备变量都按“时段 x 场景”两个维度展开约束逐条添加。注意binvar声明二进制变量的位置如果场景数很大二进制变量会非常多这时候尽量在Yalmip里用维度矩阵方式声明比用循环一件一件定义要快得多。调用Gurobi之前不要忘了在代码开头设置求解器参数options sdpsettings(solver,gurobi,gurobi.MIPGap,0.001,verbose,2); optimize(Constraints, Objective, options);MIPGap设为0.001表示允许相对最优间隙在千分之一既保证精度又不至于卡死。如果是调试阶段可以放宽到0.01会快很多。4.3 参数标定与代码调试心得参数标定这件事想认真做起来比写优化模型还费时间。设备参数一般在手册上都能查个大概但实际使用时要根据自己系统的规模调整。比如CHP的燃料系数a、b、c如果直接把某区域电网的平均供电煤耗套进去成本和设备实际消耗的对应关系就会失真。我的习惯是先做一次确定性优化的基准运行把结果和历史运行数据对照一下如果各设备出力比例和实际相差很大就反推是哪个效率参数偏了。代码调试的第一个重灾区是量纲一致性。Matlab本身不自动检查单位最容易出现的是电功率用了kW但热功率还停留在MW导致热平衡约束完全失效。我建议在数据准备阶段就统一单位比如全部用kW和kWth天然气消耗量用kW或者m³单独转换确认。第二个重灾区内置变量名冲突比如循环变量i和表达式里面的虚数单位i冲突这个在复数运算多的脚本里特别容易出问题。第三个是Yalmip的x、y等函数容易和变量名撞车所以自定义变量名尽量用带语义的字符串。调试的时候建议先跑一个单场景的确定性优化确认代码逻辑没有报错、结果趋势合理再扩展成多场景随机优化。这样定位问题会快很多不然几十个场景一起跑一旦出问题很难判断是模型问题还是场景数据问题。5. 仿真案例与结果分析5.1 典型日场景设置为了说明随机优化的价值我做了一个典型日仿真系统包括1台CHP机组、1台燃气锅炉、一组电储能和一组蓄热罐。风电和光伏各配置了一定容量并给它们设置了比较强的波动场景。电负荷预测曲线设成了传统的双峰曲线早高峰在8点到11点晚高峰在18点到22点。热负荷则是夜间和清晨偏高白天略低这是北方供暖期典型特征。风电场景设成了夜间大发、白天减少的模式光伏则相反集中在中午出力。这个时间是典型案例晚上风电多但电负荷相对低正好可以用来蓄热或者给电池充电白天光伏大发但晚高峰时光伏出力已经衰减需要靠储能和机组顶上。场景生成阶段我用ARMA(1,1)模型生成了误差序列初始场景数取500削减后保留20个。这样模型规模不至于太大又能比较充分地代表不确定性。为了对比还设置了一个“确定性优化”对照组即以预测曲线为唯一场景进行优化。5.2 对比实验确定性优化 vs 随机优化两种方式下优化得到的机组调度计划有显著差异。确定性优化给出的方案在一些时段非常激进比如它会在夜间风电预测出力较大时尽量压低CHP电出力把多余的电用来给电池充电或蓄热罐蓄热。这个计划在预测曲线成立的条件下成本很低但一旦实际风电出力达不到预测值系统的电功率就会出现缺口。多场景模拟回测时确定性方案出现切负荷和切热负荷的时段数明显高于随机优化方案。随机优化因为同时考虑了20个场景会在调度时主动留出一部分备用空间。例如在风电出力预测较大但波动也大的时段它不会把所有风电都视为可靠出力而是把储能的充电功率稍微压低或者让CHP保持在一个更稳妥的电出力水平。代价是随机优化期望成本比确定性优化高出4%到6%但最恶劣场景下的失负荷量大幅降低而且供热可靠性指标明显更好。这个“成本升高但风险下降”的权衡正是随机优化研究的核心价值所在。做结果分析时不能只看成本绝对值还要把可靠性指标和风险指标放到一起看。我习惯用条件风险价值CVaR来衡量风险把最差5%场景的成本平均值作为风险指标这样优化时就可以在目标中加入CVaR项得到更稳健的调度方案。当然这个属于扩展内容基础版本只做期望成本最小化也完全可以。5.3 关键结果解读与能效指标仿真结果输出后需要用几个关键指标来衡量系统的运行水平。第一个是CHP平均热电比用来判断机组工作点是否落在较优区间。第二个是电储能和蓄热罐的日循环次数如果一天内频繁满充满放说明储能被过度调用实际运行寿命会受影响。第三个是综合能效比即系统输出电热之和与输入燃料热值之比。热电联供系统的优越性就体现在这个综合能效比上一般可以做到0.8以上。我在结果里还专门检查了热平衡约束的残差。随机优化后由于蓄热罐在其中起了缓冲作用大多数时段热平衡都是严格满足的。但在极端场景下出现过蓄热罐达到容量上限导致多余热量无法吸收的时段这时候系统不得不通过降低CHP热出力来配合热负荷若仍不足就必须由燃气锅炉顶上。这种边界情况如果不做仿真回测单看优化结果是察觉不到的。做结果分析时另一个值得关注的细节是随机优化中不同场景下同一时段的CHP出力会有一定散布范围这个散布范围其实是系统应对不确定性的“策略空间”。可以把它画成概率带图和确定性优化的单点曲线作对比视觉冲击力很强也便于在论文或汇报中说明随机优化的必要性。6. 常见问题与排查技巧实录6.1 场景数量怎么选才合理场景数量的选择没有绝对标准但有一条经验区间如果采用蒙特卡洛抽样初始场景数取200到1000之间削减到10到30个效果通常都不错。太少比如只留5个场景随机性表达不够优化结果会和确定性方案差不多太多比如留100个计算时间成倍增长而且贡献有限。我自己试过初始场景500、削减到50和削减到20两种情况期望成本差异不超过2%但求解时间后者快了一倍多。所以如果只是做课程设计或者论文中的一个算例我建议削减到15到25个即可。如果是更复杂的热电联供系统需要加入机组组合变量建议场景数控制在20以内否则MILP的整数变量会撑爆内存。6.2 求解慢、内存占用高怎么办场景法最麻烦的问题就是求解时间太长。慢的时候一个24小时、20场景、含二进制变量的MILP问题可能也要跑几分钟。如果模型规模大跑几个小时也正常。我常用的加速手段有三个。第一尽量用连续变量近似二进制变量这需要仔细分析哪些逻辑约束是必须的。比如电池充放电互斥约束可以先用一个很小的惩罚项代替让模型允许极小概率下同时充放电求解速度会有提升但要在结果检查时确认不会有实际问题。第二在Yalmip里定义变量时用维度矩阵不要用cell矩阵后者的索引开销很大。第三给Gurobi设置MIPGap比如0.005而不是默认的0这样求解器可以早早返回一个可接受的精化解能省下大量时间。一个容易被忽略的点是场景削减后的场景权重总和是否严格等于1。如果削减代码里有数值误差导致权重和不等于1目标函数中的成本项就会整体失真求解器会出现一些奇怪的解。我建议在调用optimize之前加一行断言检查权重和与1的偏差。6.3 热负荷与电负荷的时间尺度错配问题热电联供微网里最隐蔽的问题是电负荷和热负荷的波动特性完全不同。电负荷响应速度是秒级到分钟级热负荷因为有建筑热惯性的缓冲变化相对平缓而且供热管网的传输延迟也让热负荷的“实际感受”滞后于热源出力。如果优化模型用同样的时间分辨率去处理电和热比如都用一小时一个时段热负荷延迟效应就会被完全忽略。这样得到的蓄热罐放热策略可能在真实系统中过早地把热量放完了等热负荷真正升高的时候蓄热罐里已经没有余量。解决这个问题有两种常见手段一是把热负荷加上一段滞后模型比如把当前时段的热负荷与上一时段的平均供热输出绑定二是在热平衡约束中引入一个热惯性系数模拟建筑允许的温度波动范围。这些措施会增加模型复杂度但能明显提高调度结果在实际运行中的可行性。我还有一个经验如果条件允许热负荷模型的时段可以设成电负荷的两倍即热力部分的调度周期拉长到2小时这样能自然缓解热惯性问题。当然这要看具体系统的实际情况不是所有场景都适用。除了上面这个热惯性问题还有一类典型坑来自蓄热罐的离散状态。有些同学把蓄热罐的热损失系数设成常数但实际运行时蓄热罐顶部和底部的温度分层不一样热损失不是一个固定比例。做研究时用常数近似可以接受但如果做工程应用建议蓄热罐模型至少分两层温度节点否则长期运行的累计误差会非常可观。最后再分享一点我的实际体会一个研究课题的价值不完全体现在算法有多高深更多在于能不能把问题背后的物理特性看清楚。源荷随机特征这件事看起来只是数学建模但真正理解了热电联供机组“以热定电”的运行边界理解了蓄热罐在时间维度上的缓冲作用再回头看那些随机场景、多场景优化就都有了落地的抓手。如果想把这个课题向下延伸比较值得做的是把模型扩展到多微网协同或者把运行优化和容量配置放到同一个框架里做两阶段优化。前者需要面对更复杂的荷侧不确定性后者则要处理设备投资与运行成本的耦合。无论哪条路Matlab加Yalmip加成熟求解器的这套技术栈都还能继续用模型改起来也比较顺手。最后给一个小建议拿到任何优化模型的时候务必先花半小时检查约束的物理意义再去盯求解器的设置很多看起来是“求解失败”的问题本质上都是约束里藏着一两个量纲或者索引错误。