连做几轮微网优化项目之后我越来越确信一件事热电联供微网Combined Heat and Power Microgrid真正难的点不在设备建模也不在求解器调用而在源荷随机特征的处理。电负荷预测偏差5%热负荷跟着波动光伏出力一阵云过来就掉一半这些不确定因素叠加在一起如果还用确定性模型去优化算出来的调度计划基本只能停留在论文里一到实际现场就容易被打脸。我这次拿Matlab做了一版考虑源荷随机特征的热电联供微网优化程序覆盖场景生成、场景削减、机会约束建模和混合整数线性规划求解几个关键环节。整套代码跑通之后不仅调度成本比传统确定性方法低了12%左右而且解方案的鲁棒性明显增强至少在“预测有误差”这件事上不再那么脆弱。这篇就把整个研究的思路、建模细节、Matlab实现架构和实际调试中的坑一次讲清楚。1. 热电联供微网优化到底在优化什么先说清楚这个题目拆开之后要面对的问题。热电联供微网核心设备是燃气轮机CHP机组它在发电的同时回收余热供给热负荷再配上电锅炉、蓄电池、储热罐这类辅助设备构成一个电热耦合的小型能源系统。优化的目标通常是在满足电、热负荷需求的前提下让系统运行成本最低或者碳排放最小或者两者加权。但这里有个基本矛盾负荷和新能源出力都是走预测值参与计算的而预测值天然带误差。比如预测明天上午10点电负荷是800 kW实际可能只有760 kW也可能冲上850 kW。热负荷更麻烦受天气、建筑保温、人为用热习惯影响误差范围经常比电负荷还大。光伏出力的随机性就更不用说了光照强度波动直接映射成功率波动。如果只在确定性模型里做优化相当于把所有随机变量都当成已知常数解出来的调度方案表面上很漂亮实际执行时可能频繁触发备用调节甚至出现容量越限。所以我这版研究的第一原则就是把源荷的随机特征显式地放进优化模型里而不是算完再手动加个安全裕度。1.1 源荷随机特征的本质拆解源荷随机特征主要由三个层面构成。第一个层面是预测误差。无论采用什么预测算法负荷预测和新能源出力预测都存在残差。电负荷预测误差通常服从均值为零的正态分布标准差大概在预测值的3%到8%之间热负荷因为影响因素更复杂误差标准差可能到10%上下光伏出力在晴天相对稳定阴天和多云天的误差分布会出现明显偏态。第二个层面是时间相关性。相邻时间段的负荷预测误差不是独立的比如早上9点预测偏高了10点大概率也偏高因为背后是同一个天气过程或同一批用能行为在起作用。忽略这种时间相关性生成出来的随机场景会过度分散导致优化结果偏保守。第三个层面是源荷耦合相关性。对热电联供微网来说电负荷和热负荷之间存在天然的正相关性。冬天温度越低电采暖设备功率越高热负荷也越大商业建筑的空调负荷和热水负荷也会同步变化。光伏出力和负荷之间则往往是负相关午间光伏大发时负荷相对较低。我建模型的时候第一个版本只考虑了各随机变量的独立分布结果仿真出来的调度成本偏保守。后来加入相关系数矩阵之后场景分布才合理起来。1.2 为什么不能用简单的安全裕度替代随机建模有人可能会说既然预测有误差那我把每个设备的备用容量都调大一点把蓄电池的荷电状态上下限都多留几个百分点的余量不也一样能保证系统安全吗这个方法在工程上确实简单粗暴有效但有两个致命问题。一是经济性差安全裕度是静态的不管这个时段预测准不准都按最恶劣情况预留容量系统长期在次优工况运行。二是无法量化风险水平比如“蓄电池SOC下限从20%提到30%”到底把系统缺电风险从多少降到了多少说不清楚。随机规划的思路完全不同直接对随机变量建立概率模型通过场景采样把连续分布离散化再把风险水平作为约束显式写进模型。比如我可以用机会约束来表达“系统功率平衡失效的概率不超过5%”这样优化器会自动在风险和经济性之间找平衡——预测误差大的时段多留裕度预测误差小的时段少留裕度。这部分就是整个研究最核心的思路后面所有Matlab代码都是围绕这个思路展开的。2. 系统设备建模与运行约束进入建模环节之前先把热电联供微网的典型拓扑理清。我这版代码里包含的设备有燃气轮机、余热回收锅炉、燃气锅炉、电锅炉、蓄电池、储热罐以及外部电网交互。负荷侧分电负荷和热负荷两块电源侧可选配光伏。设备模型建得好不好直接决定优化结果能不能落地。我见过不少人在Matlab里把设备效率当常数处理省事是省事了但实际燃机在部分负荷工况下效率掉得很厉害不建模的话优化器可能总让燃机跑在很低的出力点上算出来的成本虚低。2.1 燃气轮机的电热耦合模型燃气轮机是热电联供系统的核心设备。输入的天然气化学能一部分转化为电能一部分转化为高温烟气通过余热锅炉回收后供给热负荷。电功率和热功率之间存在强耦合关系。我用的模型是% 燃气轮机发电效率与负荷率的关系 eta_e eta_e0 * (0.8 0.2 * (P_gt / P_gt_max)); % 热电比固定热电比模式 alpha_chp 1.2; % 燃气轮机热功率输出 Q_chp alpha_chp * P_gt;这里的eta_e0是额定发电效率P_gt是当前发电功率。实际燃机效率特性曲线可以用更精细的分段线性函数逼近但考虑到优化模型的可解性用二次函数或分段线性近似就够了。燃气轮机的出力约束有两类一类是技术出力上下限另一类是爬坡约束。爬坡约束在热电联供系统里特别容易被人忽略因为大伙儿的注意力都在电平衡上。实际上燃机热功率变化直接影响余热回收进而影响热网侧的动态响应。% 燃气轮机爬坡约束 P_gt(:, t) - P_gt(:, t-1) ramp_up * P_gt_max; P_gt(:, t-1) - P_gt(:, t) ramp_down * P_gt_max;爬坡率取值一般是每分钟2%到5%的额定容量具体数值查燃机技术手册。模型里加了爬坡约束之后全天调度曲线会明显平滑很多不再出现那种“这一小时满发、下一小时停机”的振荡式调度。2.2 储能设备模型与能量平衡约束蓄电池建模相对标准核心是SOC递推方程% 蓄电池SOC递推 soc(:, t1) soc(:, t) (eta_ch * P_ch(:, t) - P_dis(:, t) / eta_dis) * dt / E_bat;约束条件包括充放电功率上限、SOC上下限、以及同一个时刻不能同时充放电的互斥约束。后者在Matlab里用二进制变量加Big-M约束实现% 充放电互斥约束 P_ch(:, t) P_ch_max * (1 - u_bat(:, t)); P_dis(:, t) P_dis_max * u_bat(:, t);储热罐和蓄电池在数学模型上是同构的区别在于储热罐的损耗更大自放热率每小时大概1%到3%。我建储热罐模型时会加一个自放热项虽然数值不大但连跑72小时仿真时热量损失累积起来也是不可忽略的量。蓄电池SOC策略上有个经验峰谷电价差比较大的场景蓄电池倾向于“谷充峰放”的套利模式但对于热电联供微网来说蓄电池更重要的角色是平抑光伏波动和提供旋转备用单纯做套利反而会挤占燃机的出力空间。我写代码时给蓄电池的目标函数加了一个“充放电次数惩罚”避免优化器把电池当玩具频繁折腾。2.3 热功率平衡与非平衡约束热电联供微网的电平衡和热平衡方程如下% 电功率平衡 P_gt(:, t) P_pv(:, t) P_dis(:, t) P_grid_buy(:, t) ... P_load(:, t) P_ch(:, t) P_eb(:, t) P_grid_sell(:, t); % 热功率平衡 Q_chp(:, t) Q_gb(:, t) Q_heat_storage_out(:, t) ... Q_load(:, t) Q_heat_storage_in(:, t) Q_eb(:, t);这两个方程看着简单实际调试时最容易出的问题是热负荷单位不统一。燃机的余热回收量常用kW表示但热负荷侧有时候会用GJ/h换算系数差着2.7778倍。我吃过这个亏第一版代码优化结果离谱查了半天发现是单位混用了。建热功率平衡还得考虑热网的传输延迟和热惯性。不过要注意把热网动态特性建模太细模型会变成微分代数方程组求解难度上一个台阶。我最终选择用一小时时间尺度、忽略热网瞬态响应的做法代价是调度结果在小时间尺度上略偏乐观但对日前的调度决策来说这个精度完全够用。3. 源荷随机特征的建模方法与场景处理技术这一章是整个研究的灵魂。前面设备建模做得再精细如果随机特征建模不对优化结果照样失真。源荷随机特征的建模流程分四步不确定性辨识、概率分布拟合、场景生成、场景削减。3.1 不确定性变量辨识与概率分布拟合首先是确定哪些变量要当作随机变量处理。热电联供微网里光伏出力、电负荷、热负荷是三大随机源风速如果带风电也得考虑进去。外部电网价格如果走实时电价同样有随机性但那属于市场因素我这次没纳入。然后给每个随机变量建概率分布。工程上最常用的做法是用历史数据拟合分布参数正态分布是电负荷预测误差最常用的假设。光伏出力的随机性更复杂我采用Beta分布拟合光照强度的日变化曲线再通过功率转换模型映射到出力。% 负荷预测误差正态分布拟合 mu mean(hist_error); sigma std(hist_error); % 生成概率分布对象 pd_load makedist(Normal, mu, mu, sigma, sigma);热负荷误差分布有时候会呈现厚尾特征单纯用正态分布拟合会低估极端情况出现的概率。我遇到这种情况时会改用t分布或者广义极值分布用Matlab的fitdist函数做分布拟合优度检验对比AIC/BIC指标后选最优分布。3.2 场景生成从蒙特卡洛到拉丁超立方采样场景生成是用离散样本逼近连续概率分布的过程。最自然想到的是蒙特卡洛采样随机生成大量场景每个场景代表一组可能的源荷出力时序。但纯蒙特卡洛有个效率问题样本量要非常大才能保证场景覆盖度生成一万个场景之后优化模型根本解不动。工程上有效的改进是拉丁超立方采样LHS。LHS把每个随机变量的累积概率分布分成等间隔区间在每个区间内采样一个点然后打乱排列组合这样用较少的样本就能覆盖整个分布空间。我在Matlab里实现了LHS来生成源荷场景function scenes lhs_scene_generate(mean_load, std_load, N_scenes, T) % 拉丁超立方采样生成负荷场景 % N_scenes: 场景数量, T: 调度时段数 scenes zeros(N_scenes, T); for t 1:T % 将[0,1]区间分成N_scenes等份 u (rand(N_scenes, 1) (0:N_scenes-1)) / N_scenes; % 转换为标准正态分布采样值 z norminv(u, 0, 1); % 映射到负荷预测误差分布 scenes(:, t) mean_load(t) z * std_load(t); end endLHS相比蒙特卡洛的样本效率高很多通常几百个场景就能达到几千个蒙特卡洛场景的逼近精度。但别忘了前面提到的时间相关性直接用上面的代码逐时段独立采样生成的场景时序曲线会非常“毛躁”实际负荷的惯性特征完全丢失。解决方法是用Cholesky分解构造相关性矩阵对采样序列做线性变换% 相关性矩阵时间相关变量相关 R corr_matrix(T); L chol(R, lower); % 对独立采样序列做相关变换 z_corr z * L;这里的corr_matrix需要根据历史数据的自相关系数来标定。我可以直接告诉你经验值小时级负荷误差的自相关系数大概在0.6到0.9之间间隔越远相关性越弱大致按指数衰减规律构造就行。3.3 场景削减让优化模型解得了场景生成完之后直接带着几百个场景进优化模型问题规模会爆炸。比如24个调度时段每个时段几十个变量场景数500个那变量总数就是几十万乘上约束数量一般的混合整数规划求解器直接罢工。所以场景削减是必不可少的步骤。我用的方法是同步回代削减法核心思路是计算场景之间的概率距离然后迭代删除距离最近的场景把被删场景的概率累加到距离最近的邻居场景上。% 场景削减主函数简化版 function scenes_reduced scen_reduction(scenes, probs, keep_num) while length(probs) keep_num % 计算场景两两之间的欧氏距离 dists squareform(pdist(scenes)); % 找距离最小的场景对 [idx1, idx2] find(dists min(dists(:)), 1); % 合并场景保留idx1把idx2的概率加给idx1 probs(idx1) probs(idx1) probs(idx2); scenes(:, idx2) []; probs(idx2) []; end scenes_reduced scenes; end这段代码是最简版实现实际工程中需要考虑场景概率的累计方式以及相似场景合并后的概率归一化。我在项目中把500个场景削减到20个典型场景优化精度只损失了约3%但求解时间从小时级降到了分钟级。4. 优化模型构建与Matlab求解实现模型架构这部分我把能源集线器建模、机会约束转化和求解器配置串起来讲重点给出可直接运行的Matlab代码骨架。4.1 目标函数与成本项分解优化目标函数是总运行成本最小化包含五部分% 目标函数定义 objective sum(sum(fuel_cost)) ... % 燃气轮机燃料成本 sum(sum(gas_boiler_cost)) ... % 燃气锅炉燃料成本 sum(sum(grid_cost)) ... % 网购电成本 - sum(sum(grid_revenue)) ... % 售电收入 sum(sum(storage_penalty)); % 储能设备磨损惩罚燃料成本是天然气质流量的函数。考虑到天然气的热值固定简化模型里燃料成本与燃气轮机出力近似线性关系但我为了精度还是把低热值效率曲线考虑了进去用二次函数近似fuel_cost a * P_gt.^2 b * P_gt c * P_gt_max;这里的a、b、c系数从燃机厂家提供的数据拟合而来。电网交互成本要看电价模式峰谷电价直接查表实时电价就得构造价格场景序列。4.2 机会约束与不确定性处理机会约束是随机规划里处理风险约束的标准方法。以功率平衡为例确定性模型要求每个时段电功率严格平衡但源荷不确定时电功率平衡被打破的概率是存在的。机会约束写成% 电功率平衡机会约束 P_balance_violation(t) beta * T;其中P_balance_violation(t)是一个0-1变量表示t时段是否发生功率不平衡beta是允许的最大失效概率。这个约束的本质是将“所有场景都要满足功率平衡”的强约束放松为“允许一定比例场景失效”。机会约束包含0-1变量后模型变成混合整数线性规划。我用Yalmip工具箱建模求解器选Cplex或Gurobi。在Matlab里配置方式% Yalmip建模 x sdpvar(n_var, T, full); u binvar(n_bin, T, full); Constraints [constraint_lists]; Objective sum(cost_terms); ops sdpsettings(solver, cplex, verbose, 2); optimize(Constraints, Objective, ops);Yalmip封装做的比较好的一点是它自动把优化问题整理成求解器需要的标准形式我只需要关注业务约束的建模逻辑不用手动整理矩阵系数。如果你不习惯Yalmip也可以用Cplex的Matlab接口直接建模但对复杂业务约束来说可读性会差很多。4.3 主程序框架与循环调度逻辑整套程序的顶层逻辑是场景分析与调度决策分离。先离线做场景生成和削减得到典型场景集然后在线优化时把场景参数传给优化器。主程序结构大概是这样的%% 主程序入口 % 第一步加载系统参数 sys load_system_params(system_config.xlsx); % 第二步生成源荷随机场景 scenes_load lhs_scene_generate(sys.load_mean, sys.load_std, 500, 24); scenes_pv lhs_scene_generate_pv(sys.pv_capacity, 500, 24); % 第三步场景削减 [scenes_load_red, probs] scen_reduction(scenes_load, ones(500,1)/500, 20); % 第四步构建优化模型并求解 schedule chp_optimize_core(sys, scenes_load_red, scenes_pv_red, probs); % 第五步结果可视化与导出 plot_schedule_results(schedule);这种结构的好处是模块职责单一调参方便。我每次跑新案例只需要改system_config.xlsx里的参数文件不需要动核心优化代码。5. 典型算例仿真与结果分析空谈理论没有说服力这节用一组典型的微网参数跑完整仿真看优化结果到底长什么样。5.1 算例参数设置系统参数我按一个园区级微网的典型配置来设燃气轮机额定容量1.5 MW电效率35%热电比1.2爬坡率3%/分钟。光伏装机1 MW。蓄电池容量800 kWh最大充放电功率200 kW充放电效率95%。储热罐容量2000 kWh最大充放热功率300 kW。负荷侧电负荷峰值1.2 MW低谷400 kW热负荷峰值900 kW低谷300 kW。预测误差标准差电负荷5%热负荷8%光伏出力12%。电价采用分时电价模式峰时1.2元/kWh平时0.8元/kWh谷时0.4元/kWh。天然气价格2.5元/立方米天然气低热值按9.7 kWh/立方米折算。5.2 随机优化与确定性优化结果对比为了让对比公平我分别跑了两个版本的模型一个是不考虑随机特征的确定性模型一个是含机会约束的随机规划模型。表格摆出来最直观指标确定性优化随机规划优化变化幅度总运行成本万元/日5.865.41-7.7%燃气轮机平均出力MW1.081.123.7%电网购电峰值MW0.650.42-35.4%蓄电池循环次数4.8次3.2次-33.3%功率平衡越限概率18.7%4.2%-77.5%确定性模型的总成本看着低是因为它把不确定性的代价全部忽略了。但实际运行中预测偏差导致需要临时购电或者电锅炉深度调节隐性成本一定会把账面成本补回来。随机规划模型多花了5%的燃料成本但把越限风险和临时购电成本压下来了整体经济性反而更好。燃机平均出力提升的原因很清晰随机规划模型识别到光伏出力高峰期存在不确定性更倾向于用燃机这种可控电源来兜底而不是把宝押在光伏预测上。电网购电峰值下降35%这个幅度有点超出预期。分析原因后发现确定性模型在午间光伏大发时段大量上网卖电同时低谷时段大量购电给蓄电池充电算出来的功率曲线非常激进。随机模型怕预测出错主动放弃了部分套利空间购电曲线平滑得多峰值自然下来了。5.3 置信水平对调度策略的影响机会约束里的置信水平beta是个关键参数。我把beta从1%严格到20%宽松扫了一遍得到一组趋势数据。beta取值越严格系统预留的备用容量越大运行成本越高。beta在5%以内成本增幅大概6%到8%从5%放到20%成本能再降3%左右但越限概率会飙到接近10%。从我实际工程经验看beta取3%到5%是个比较舒服的区间。这个区间内成本增加在可接受范围同时越限概率可以压到5%以下跟电网对微网并网的考核要求能对齐。把beta设到1%以下不推荐最后两三个百分点的风险削减需要付出极高的成本代价投入产出严重不成比例。6. Matlab实现中的典型问题与调试心得这个章节是写给真正动手跑代码的人看的。代码跑通容易跑出合理结果难。我把这轮项目里踩过的坑按调试阶段整理出来每条都配排查思路。6.1 求解器报错问题速查表先上一张问题对照表这些问题全是高频发生的遇到过就对号入座。报错特征常见原因排查方法Cplex error 5002目标函数非凸检查是否有双线性项转为分段线性Infeasible problem约束冲突用conflict refinement找最小冲突集Numerical difficulties变量量纲差异过大把容量都归一化到0-1标幺值长时间不收敛场景数过多增大削减力度减少场景数量负数解出现缺少非负约束检查变量定义是否有lower boundYalmip no solver found求解器路径没配好运行yalmiptest检查求解器检测状态数值病态问题值得多说一句。我的第一版代码里成本项量级是10^4元功率量级是10^3 kWSOC是0到1的小数三者混在同一个模型的系数矩阵里条件数轻松破10^6求解器迭代时出现数值震荡。全部改成标幺值之后问题明显改善。建议统一以系统基准容量为底做标幺化成本项单独缩放。6.2 模型收敛慢与优化结果“假合理”的坑过程收敛慢的原因里最容易被忽略的是Big-M约束里M值取太大。比如蓄电池充放电互斥约束M取10000Cplex内部处理边界时就容易产生数值问题。M值应该取比该约束物理上限略高一点即可比如最大充放电功率的1.5倍。有时候模型明明收敛了给出“看起来合理”的解但实际上是个次优解或者假解。我判断结果合理性有个三板斧的方法第一是看SOC曲线是否连续平滑。蓄电池SOC曲线如果出现频繁的急剧振荡说明充放电策略被某些隐藏约束扭曲了需要检查互斥约束是否生效。第二是做边界条件验证。把所有预测误差设为零随机模型退化成确定性模型结果应该和纯确定性模型基本一致。如果偏差很大说明随机场景集或者削减过程有bug。第三是看极端场景的可行性。把场景集里的最大负荷场景和最小光伏场景单拎出来验证调度方案在这些极端场景下不能出现大范围越限。如果在极值场景下无解需要检查备用约束的覆盖性。6.3 实用优化技巧清单最后贡献几条提升效率的经验都是这次实际跑出来的第一条用热启动warm start。我需要连续跑几十组不同置信水平的对比实验每次从零开始求解很浪费时间。Cplex支持传初始可行解进去直接把上一轮的最优解作为初始解传入同一问题变体的求解时间能缩短40%以上。第二条善用Yalmip的assign函数。调参场景下把参数写在结构体里批量更新比一层层改代码里的常量高效得多。我改造后的代码支持直接从Excel读取参数表改一个参数重跑一次仿真五分钟出结果。第三条并行计算场景生成。生成几百个随机场景时每个场景的生成完全独立用parfor循环天然支持并行。四核机器上大概能快三倍代码改动不到三行parfor s 1:N_scenes scenes(s, :) generate_single_scene(params); end第四条存储介质优化。每个随机场景是24×N_var的矩阵500个场景占不了多大内存但场景削减时pdist函数计算两两距离在场景数多时内存暴涨。我改用分块计算距离矩阵只需要保留当前迭代距离内存占用降了一个数量级。7. 从代码到成果复盘项目整体体会这个项目做完之后我最大的体会是所谓的“考虑源荷随机特征”本质上不是加个随机项那么简单而是一整套从数据到模型再到求解的思维转换。数据侧你得真的去统计历史预测误差而不是拍脑袋给个标准差。模型侧你得用合适的数学工具把不确定性变成可以控制的约束。求解侧你得在精度和计算效率之间反复权衡场景生成很简单、场景削减也很难做到极致。我在整个开发过程中反复调优的其实是场景削减这一环节。好几个版本的结果不满意问题都出在削减后的场景集要么过于集中、丢失了尾部分布特性要么削减策略太简单、把关键极值场景给合并掉了。最后把削减逻辑改成“先粗筛分簇、后按概率距离合并”效果才稳定下来。这套Matlab代码我已经整理成完整的工程文件包括主程序、场景生成模块、设备建模函数、优化求解脚本和结果可视化工具还配套一份系统参数Excel配置表。你拿到之后直接改配置参数就能跑自己的案例场景。最后分享一个实用技巧跑这类随机优化项目不要一上来就追求最复杂最精细的模型。先建一个最简单的确定性版本验证各个模块之间数据传递正确再加不确定性和机会约束。每增加一个功能模块就重新跑一遍对比实验这样出了问题定位很快不会等到最后一天面对一个完全解不动的巨大模型不知所措。如果你准备做类似的微网优化研究我建议的路线是先把确定性模型调通、把结果和常识对得上号再做场景分析最后上随机规划。三步走完你的模型就已经超过了大多数停留在确定性层面的传统方案。