前阵子一个师弟做毕业设计拿了一个热电联供微网优化模型跑了很久在仿真里结果一直很漂亮可一换数据就崩最后来问我我一看就明白了——他那套模型里,源和荷全用的是确定值。这个方向其实很有意思也很折磨人热电联供微网优化本身就是一个多变量、多约束、强耦合的问题再把源荷随机特征扔进来模型瞬间从“麻烦”变成“复杂系统”。但恰恰是这一步才是从“能发论文”走向“能落地”的分水岭。今天我就以自己复现和改造相关模型的经验把“考虑源荷随机特征的热电联供微网优化”这件事从头到尾拆一遍包括建模思路、不确定性怎么处理、MATLAB里怎么实现以及我在实际调试中踩过的坑。需要说明的是下面这部分内容基于该领域目前的主流研究框架和通用工程实践结合我个人对这类模型的复现经验来写代码和参数只代表一种可用的方案不一定是最优雅的但足够你跑通并且在这个基础上往上加东西。1. 为什么说源荷不确定性是热电联供微网调度的“隐形天花板”1.1 确定性调度在真实场景中的尴尬先讲一个很直观的对比。确定性调度做起来很爽给定一条负荷曲线给定一个风电光伏出力序列优化器把各台机组的出力、购电量、储热罐充放热一次全算出来目标函数最小约束全满足画出来的图也是一条条平滑漂亮的曲线连审稿人看着都舒服。但你想过没有这条“最优调度曲线”存在的根基是——你给优化器的负荷预测和新能源出力预测是百分百准确的。现实里可能吗不可能。风电的随机波动、光伏的云层遮挡、用户负荷在一天里的突发变化这些不确定性会直接导致一个结果按照确定性优化算出来的调度计划执行实际运行中会出现功率不平衡、热负荷供应不足、频率电压越限甚至被迫弃风弃光、切负荷。我和不少同行交流过一个共识确定性优化在纯理论推演、案例分析、教科书示范里有价值但只要你手里的微网真实接了风机、光伏、电锅炉、储热罐确定性结果只能当参考不能当执行方案。1.2 热电联供系统的双层不确定性链条热电联供微网比纯电力微网更麻烦的地方在于——它有两条能量流电和热。这两条能量流不是独立的CHP机组发多少电往往就带出多少热取决于机组类型。你调整供电策略供暖平衡就跟着动你调蓄热罐电平衡又受影响。这种强耦合关系下不确定性会在两条链上同时传导电源侧不确定性风机出力随风速波动光伏出力因云层变化呈间歇性这两个变量直接冲击电功率平衡。负荷侧不确定性电负荷随用户行为随机波动热负荷受天气、建筑围护结构、供热方式影响时间和空间上都有明显的随机特征。更麻烦的是负荷侧的电和热本身还有相关性。比如冬天的傍晚大家下班回家电负荷上升热负荷也上升夏天的中午空调制冷负荷飙升热负荷反而处在低谷。这种相关性如果处理不好会导致你生成的不确定性场景失真优化结果出现系统性偏差。1.3 三类主流处理方法随机优化、鲁棒优化、区间优化要处理源荷随机特征现在学术界和工程界常用三套框架各有各的适用场景。随机优化Stochastic Optimization对随机变量的概率分布进行采样生成大量场景用期望值做目标做两阶段或多阶段决策。优点是贴近实际能利用概率信息缺点是计算量大且需要一个“可信”的概率分布。鲁棒优化Robust Optimization不依赖具体概率分布只给不确定性变量一个“集合”区间、盒式、椭球式等优化结果在所有集合内的最坏情况下依然可行。优点是稳健、对概率分布不敏感缺点是偏保守最坏情况不常发生经济性会有损失。区间优化Interval Optimization是鲁棒优化的一种简化形态直接用区间上下界描述不确定量思路简单、求解快但信息利用率低结果往往更粗糙。我做复现时最常用的是“随机场景 两阶段”的方式因为它在经济性和鲁棒性之间比较好平衡而且MATLAB里用YALMIP配合商业求解器处理起来非常顺手。具体怎么搭往下看。2. 热电联供微网的建模底层目标函数与电热耦合约束2.1 从“电热解耦”到“电热耦合”的变化很多刚开始接触这个方向的人容易陷入一个思维定式把热电联供系统拆成“电网子问题 热网子问题”两个独立模块分别优化后再合并。这个问题我劝你从一开始就别这么做。热电联供的核心特征是“以热定电”或“以电定热”具体取决于机组类型和运行策略。抽凝式机组的热电比在一定范围内可调背压式机组则是“发多少电就排多少热”强制耦合。你把电热拆开等于把最关键的耦合约束扔掉了算出来的“最优”方案在物理上根本可能不成立。正确的做法是用一组耦合变量把两条能量流绑在一起在同一个优化模型里同时决策。说白了就是在一个大的优化问题里同时写电平衡、热平衡、CHP出力关系、储热罐状态、购电交互等约束让求解器自己去协调。2.2 热电联产机组的可行域抽凝式与背压式的不同CHP机组的建模是整个热电联合优化的基础也是新手最容易出错的环节。背压式机组汽轮机排汽全部进入热网加热器放热发电功率和热功率之间存在一个近似固定的比例关系[ P_{chp} c·Q_{chp} ]其中c对某一台确定机组是一个常数如0.4~0.7之间浮动。这类机组模型简单但因为没有调节自由度在运行中灵活性很差。我在实际项目里一般不单独用背压式除非场景固定、热负荷稳定。抽凝式机组可以从汽轮机中间级抽取部分蒸汽用于供热其余蒸汽继续发电。它的电出力和热出力是一个二维可行域通常用多边形的顶点来刻画[ \begin{aligned} P_{chp}^{min} \leq P_{chp} \leq P_{chp}^{max} \ 0 \leq Q_{chp} \leq Q_{chp}^{max} \ P_{chp} \geq \alpha_{chp} Q_{chp} \beta_{chp} \end{aligned} ]实际建模时根据机组的不同工况可行域可能是一个五边形、六边形甚至更复杂的多边形。我建议你先把机组的技术资料里的运行工况点拿到然后用顶点枚举法把可行域描述出来再转成线性约束。一个要注意的细节是很多论文为了简化把不等式的斜率/截距写成一个固定的常系数。但真实机组的可行域往往是非凸的直接线性化会引入误差。对于硕士阶段的项目线性化精度通常够了如果是博士课题或工程应用还是建议用分段线性化甚至混合整数建模把非凸区域尽量逼近。2.3 热网动态约束与热负荷平衡条件热网和电网不一样电是“瞬时平衡”热是有“惯性”的。热水在管道里流动有延迟建筑本身有热存储效应所以热负荷平衡不一定需要每一时刻严格等于供热端出力。但很多初版模型会直接把热平衡写成静态等式[ Q_{chp} Q_b Q_{hs,out} - Q_{hs,in} Q_{load} ]这个写法在较长调度时间段比如24小时步长1小时里可以接受但如果你的步长缩小到15分钟甚至更短就必须考虑管道延迟和建筑热惯性。否则算出来的热出力曲线会剧烈跳动和实际系统完全对不上。我在实际项目里处理这个问题通常分两步走在规划级优化里用静态热平衡步长1小时先求整体的调度策略。在执行级日内滚动里引入热网延迟系数和储热罐SOC做精细化校核。这样做的好处是既不会让规划模型因热网动态约束维数爆炸而不可解又能在执行层面把热网惯性利用起来。3. 源荷随机特征的数学抽象概率、场景与不确定性集合3.1 风光出力随机特征从时序波动到概率分布先说风电。风速本身服从威布尔分布风电出力与风速之间又有非线性关系切入风速、额定风速、切出风速三段所以直接对出力做概率统计更实用。你可以用历史运行数据通过核密度估计拟合出风电场出力的概率密度函数然后生成大量随机出力曲线。光伏相对简单一些主要受光照强度影响。同一时刻的光照强度可以用Beta分布近似描述光伏出力近似等于光照强度乘以一个转换效率。但注意光伏还有一个特点日内的时序相关性很强上午到中午上升、下午下降你不可能把每个时刻都当作独立随机变量来处理。我的做法是先对24小时的风光出力分别建一个“基准场景”和“误差场景”。基准场景取预测值误差场景从历史预测误差的概率分布中抽样然后把二者叠加。这样既保留了时序特性又引入了随机性。3.2 负荷预测误差的统计特征与相关性电负荷和热负荷也都不是确定值。对负荷建模我通常用“预测值 误差项”的结构[ L_{load} L_{load}^{forecast} \varepsilon ]误差项ε的分布一般假设为零均值正态分布标准差与负荷量级和历史预测精度相关。如果你手里有微网的历史负荷数据和对应预测数据可以直接算出误差序列再统计出标准差。这里有一个我踩过的坑负荷误差不是每个时刻独立的。比如晚间高峰出现预测偏差往往连着好几个小时都偏高或偏低因为造成偏差的天气原因通常会持续一段时间。所以单纯用独立正态分布生成场景场景会过分“毛躁”和真实情况的持续性偏差不符。处理办法是引入时间相关性可以用一阶自回归模型AR(1)来生成误差序列或者先对误差做主成分分析保留主要模态再加随机扰动。后一种在数据量不够时更稳。3.3 场景生成与场景缩减从一万条曲线到几条代表曲线有了随机变量的分布下一步就是生成场景。蒙特卡洛抽样是首选比如对风、光、电负荷、热负荷四个随机变量每个变量抽1000次样组合起来就产生1000个“四维场景”。每个场景就是一组24小时的风出力、光出力、电负荷、热负荷曲线附带一个概率值等概率的话就是1/1000。但1000个场景直接代入优化模型别说是非线性哪怕是线性规划求解器也会被拖垮。所以必须做场景缩减把1000个场景聚合成少量的代表性场景比如5~10个然后每个场景分配一个权重等于被削减到该场景的原始场景概率之和。常用的场景缩减方法有快速前向选择法同步回代消除法为主每次找一对“最相似”的场景把其中一个的权重合并到另一个重复直到场景数量达标。K-means聚类型方法把每个场景看作高维空间的一个点聚类中心就是代表性场景。MATLAB里这两个方法都有现成实现也可以自己写。我在实际项目中更倾向用同步回代消除法因为它在保持概率分布形状方面效果好并且实现起来也不复杂。4. 考虑随机特征的两阶段优化框架从“调度计划”到“再调整”4.1 两阶段决策变量怎么划分两阶段随机优化是处理源荷不确定性的一个非常自然的框架。它的核心思想是把决策分成“现在必须做、还没看到随机变量实现值”的决策和“看到随机变量实现值后、可以调整”的决策。在热电联供微网里通常这样划分第一阶段决策Here-and-Now机组的启停状态0/1变量与上级电网的购电/售电的中长期计划储热罐的充放热策略基准值这些决策必须提前确定且在调度周期内不能轻易改变。第二阶段决策Wait-and-See各机组在每个场景下的实际出力调整量切负荷量/弃风弃光量的实时修正储能设备的实时充放电调整这些决策是在每个具体场景即不确定性实现后内做出的为的是保证系统在该场景下可行同时尽可能经济。模型写成紧凑形式就是[ \min_{x} \left( c^T x \sum_{s1}^{S} \pi_s Q(x, \xi_s) \right) ]其中(Q(x, \xi_s))是场景(s)下的第二阶段最优成本(\pi_s)是场景概率x代表第一阶段决策。4.2 不确定性集合在阶段二怎么发挥作用如果你用的是随机场景方法第二阶段就是在每个场景(\xi_s)下分别求解一个确定性优化问题然后取加权期望。如果你决定用鲁棒优化框架那就不需要场景了取而代之的是一个不确定性集合U。优化问题变成[ \min_{x} \left( c^T x \max_{\xi \in U} Q(x, \xi) \right) ]内层的“max”找最坏情况外层“min”优化第一阶段决策使最坏情况下的成本最小或可行性最好。最常用的不确定性集合是“盒式 预算约束”的组合也就是每个随机变量在预测值(\bar{\xi})附近可波动但整个调度周期内波动总量受一个预算控制。这样做的目的是避免模型无限保守——允许每一天、每一时刻都取最坏值那优化结果会保守到没法用。预算参数Γ读作Gamma是一个关键参数它实质上是“你愿意为鲁棒性牺牲多少经济性”。Γ越大越保守总成本越高Γ越小越接近确定性模型。实际调参时我一般从Γ0开始逐步增大观察总成本和“违背约束风险”的曲线拐点选择一个经济性损失小于5%——10%的最大Γ值。4.3 结合经典研究思路的参考框架王锐等人在含可再生能源热电联供型微网方面的研究思路大致可以概括为“不确定性建模 鲁棒/随机优化 电热耦合调度”这个框架在业界被广泛借鉴。我自己复现这类模型时通常会再加一层“运行风险评估”也就是在优化完成后统计各个不确定性场景下的越限情况把“约束被违背的概率”作为输出指标。这样做有一个很大的好处审稿人或者领导不会只盯着总成本数字他们会问“你这个方案如果遇到极端天气扛不扛得住”。你把风险评估结果拿出来比任何解释都更有说服力。5. MATLABYALMIP实现从搭建模型到求解的完整流程5.1 为什么用YALMIP而不是手写求解器MATLAB里做优化建模我强烈建议用YALMIP。原因有三第一YALMIP的语法非常接近数学表达变量定义、约束写入、目标函数设置几乎和论文公式一一对应调试效率极高。第二YALMIP后端可以无缝切换Gurobi、CPLEX、Mosek等商业求解器你可以先拿Gurobi跑再拿CPLEX验证结果一致性不用改模型代码。第三YALMIP对混合整数线性规划MILP、混合整数二阶锥规划MISOCP等复杂模型的支持很成熟而热电联供微网优化恰恰经常要处理整数变量和非线性项。如果你非要用MATLAB自带的linprog、intlinprog也不是不行但处理大规模场景、多个整数变量时求解速度会慢很多。一个24小时、10个场景、20台设备的热电联供优化问题用intlinprog可能要好几分钟Gurobi一般几秒上下。5.2 求解器选型与参数设置这里单独说下求解器的选择。模型是纯线性规划LP或者只有连续变量的二次规划QP用Gurobi/CPLEX都行差距不大。但如果包含0/1启停变量就必须上MILP求解器并且需要注意设置时间限制防止极端情况卡死。设置MIP Gap阈值比如1e-3或0.1%避免求解器为了追求无穷精度一直跑下去。开启多线程充分利用电脑CPU。在YALMIP里设置这些参数很简单ops sdpsettings(solver, gurobi, gurobi.TimeLimit, 300, ... gurobi.MIPGap, 0.001, gurobi.Threads, 8); optimize(constraints, objective, ops);如果你没有商业求解器授权也可以先用MATLAB自带的intlinprog做验证模型小的时候完全够用。5.3 一个最小可运行案例的核心代码与结果解读下面我写一个极简的、但结构完整的热电联供微网“场景法两阶段”优化示例方便你对整体流程有个直观感受。这个示例还谈不上“完整工程”但把最关键的点都涵盖了。假设调度周期24小时步长1h系统含1台抽凝式CHP机组1台燃气锅炉1台储能电池可选风电场出力含随机误差电负荷预测值误差与上级电网可交互购电%% 基础参数 H 24; % 时段数 S 5; % 场景数 c_gas 0.35; % 天然气价格元/kWh热值 c_buy 0.6; % 购电价元/kWh c_sell 0.4; % 售电价元/kWh Pi ones(S, 1) / S; % 场景概率均分 % CHP机组参数 Pchp_max 200; Pchp_min 40; Qchp_max 180; % 效率系数简化的热电比线性关系 eta_chp_e 0.35; eta_chp_h 0.45; %% 场景数据风电出力以每小时为单位的预测值误差 wind_forecast 50 30 * sin((1:H)/24 * pi); % 预测曲线 % 生成场景每个场景是 1xH 的风电曲线在预测值上叠加噪声 wind_scen zeros(S, H); for s 1:S wind_scen(s, :) wind_forecast 10 * randn(1, H); wind_scen(s, :) max(wind_scen(s, :), 0); % 风电不为负 end %% 负荷场景 load_forecast 300 80 * sin((1:H)/24 * pi 0.5); load_scen zeros(S, H); for s 1:S load_scen(s, :) load_forecast 20 * randn(1, H); load_scen(s, :) max(load_scen(s, :), 0); end %% 热负荷假设比电负荷低一些形状类似 heat_forecast 200 50 * sin((1:H)/24 * pi 0.8); heat_load repmat(heat_forecast, S, 1); %% 建立优化模型 yalmip(clear) X sdpvar(S, H); % CHP发电量场景s时刻h Y sdpvar(S, H); % CHP供热量场景s时刻h B sdpvar(S, H); % 锅炉供热量 G sdpvar(S, H); % 从电网购电量 R sdpvar(S, H); % 弃风量 % 目标函数期望的总成本 obj 0; for s 1:S cost_chp_gas c_gas * ( X(s,:) / eta_chp_e Y(s,:) / eta_chp_h ) / 2; % 简化燃气成本 cost_buy c_buy * G(s,:); obj obj Pi(s) * sum(cost_chp_gas cost_buy); end % 约束 C []; for s 1:S for h 1:H % CHP可行域约束简化 C [C, Pchp_min X(s,h) Pchp_max]; C [C, 0 Y(s,h) Qchp_max]; % 电功率平衡CHP 风电 购电 负荷 弃风 C [C, X(s,h) wind_scen(s,h) G(s,h) load_scen(s,h) R(s,h)]; % 热功率平衡CHP 锅炉 热负荷 C [C, Y(s,h) B(s,h) heat_load(s,h)]; C [C, B(s,h) 0]; end end ops sdpsettings(solver, gurobi, verbose, 0); optimize(C, obj, ops);这里CHP的燃气成本我做了高度简化真实项目里要按气耗量曲线来写通常是关于电出力的二次函数用分段线性近似处理。跑完后你可以画一张图把场景1~5下的CHP电出力、购电量画出来看看figure; plot(1:H, value(X(1,:)), r-, LineWidth, 1.5); hold on; plot(1:H, value(X(2,:)), b--, LineWidth, 1.5); plot(1:H, value(X(3,:)), g-., LineWidth, 1.5); legend(场景1, 场景2, 场景3); xlabel(时刻/h); ylabel(CHP电出力/kW);正常情况下不同场景下的CHP出力曲线应该有一些差别这正体现了“考虑随机特征”的作用——每个场景对应一种可能的实现调度方案用期望总成本来平衡不会只针对某一条曲线做到最优。6. 实际项目里的坑与调优经验6.1 不确定性集合边界定多大才不假保守前面提到鲁棒优化的不确定性集合这里具体展开一下边界大小的选择。我见过很多刚接触鲁棒优化的人直接把不确定量的上下限设成年最大偏差/年最小偏差结果优化出来一堆机组全开、储能全满的“堡垒式调度”成本高得离谱实际根本用不上。正确的做法是用“一定置信水平的预测区间”作为波动范围而不是历史极值。比如对风电预测误差统计出误差的90%分位数把不确定性集合边界设为±1.28σ反正态分布的90%分位数对应约1.28倍标准差这样既覆盖了绝大多数情况又不会因为极端小概率事件让模型过度保守。预算参数Γ的调节可以参考我前面提到的“成本—风险曲线”法逐点扫描。如果你懒得扫描有一个经验值Γ取调度时段数H的1/4到1/3。比如24小时调度Γ取6~8一般的工程问题效果就不错。6.2 热网模型的时间步长与动态约束处理这个坑我前面提过但值得再强调一次。很多人为了精细化把热网管道延迟精确到分钟级建模然后和小时的调度模型耦合结果模型规模爆炸求解时间从几秒变成几十分钟而且整数变量一多经常半天出不来最优解。我的实践建议是优化层面步长取1小时热平衡用静态等式 储热罐SOC修正。校核层面针对优化结果用15分钟步长做热网动态仿真检查是否有温度越限、流量越限问题。如果确实需要在优化里加动态热网约束尽量采用“管道传输延迟 热损失系数”的简化和线性化模型并且减少整数变量否则模型稳定性很难保证。6.3 求解时间爆炸与预归约当场景数增加到50个以上设备数超过20台MILP模型的规模会迅速膨胀求解时间可能从秒级跳到分钟级甚至小时级。几个有效的手段减少整数变量很多启停变量在部分时段可以通过逻辑关系提前固定比如热负荷高峰期CHP必开这类变量可以直接赋固定值减少分支定界的搜索空间。场景预处理如果两个场景之间相似度过高做一次场景缩减再进优化比直接丢100个场景进去要快得多。给求解器一个热启动初值先用确定性模型算一个解作为两阶段模型的整数变量初值传入能大幅加快MILP收敛。松弛检验先把整数变量松弛成连续变量跑一遍LP看看目标和约束是否合理。如果LP都解不出来那就不是求解器的问题而是模型本身有bug。6.4 不收敛和奇异解的排查清单最后给一个排查清单都是我实际调试时遇到过的问题约束写成了双向不等式且两边方向反了检查每个等式的左右单位是否一致。储能SOC出现非物理值忘了加SOC范围约束或者在循环中SOC递推写反了正负号。场景权重没归一化导致目标函数量级失真优化器给出的解偏向某个权重异常大的场景。热平衡约束缺少松弛变量当负荷和供热量在某个场景下无论如何都配不平尤其是极端场景模型会直接报不可行。此时需要检查场景生成有没有产生明显不合理的极限值必要时加少量松弛变量并惩罚。求解器设置过于激进MIPGap设得过于接近0求解器为了证明最优性会一直跑实际解早已稳定。工程上设0.1%或0.5%就足够了。我自己复现这类“考虑源荷随机特征的热电联供微网优化”模型最大的体会是模型不是越复杂越好关键是把不确定性的特征描述准确、把电热耦合关系写对、把求解环节调稳这三点做到位你的模型就比大多数“套模板”的实现要靠谱得多。尤其是场景缩减和不确定性预算这两个环节看似小实际上决定了整个模型的经济性和可解释性值得多花点时间琢磨。