小区微网的综合能源调度很多做优化的同学第一步就卡在不确定性三个字上。风电和光伏出力是波动的负荷也不是一成不变如果拍脑袋把新能源出力当成固定值算出来的调度方案一落地就偏差严重的会直接导致切负荷或弃风弃光。我最近在做一个电气设备综合能源系统的协同优化项目核心就是用Matlab把新能源出力的不确定性考虑进去让电、气、热多种能源在设备层面协同运行。这套代码跑通之后不光解决了调度方案不贴合实际的问题还顺带把系统运行成本降了一截今天就把建模思路、代码实现和踩过的坑完整梳理一遍给正在做能源优化、写论文、搞毕设的朋友当个参考。这个内容对三类人最有用一是电气工程、能源系统方向的研究生想复现带随机优化性质的协同调度模型二是做园区或微网能源管理的工程师需要让调度算法兼顾安全性和经济性三是刚接触MatlabYalmip优化建模的初学者想搞清楚场景生成、机会约束、设备协同这些概念到底怎么落成代码。1. 项目思路拆解不确定性与协同优化到底在解决什么问题1.1 新能源出力不确定性为什么不能简单忽略很多初版模型喜欢把风电、光伏在某时刻的出力用一个固定值代替然后去求解一个确定性的优化问题。这么做的确简单但结果往往很脆弱。风电出力受风速影响一天内的波动可以超过额定功率的八成光伏更是白天有、晚上没有碰上多云天气出力能在几分钟内掉一半。用平均值或者预测值代表真实出力相当于默认预测零误差这在工程上是不成立的。不确定性处理不好最直接的后果就是优化结果不可信。比如某个时段系统按预测的风电出力安排燃气轮机少发结果实际风速偏低风电实际出力比预测值小燃气轮机又因为启动慢顶不上最终只能切除负荷或者从上级电网高价买电调度成本远超预期。反过来如果预测偏保守明明风很大却让燃气轮机满发又会造成不必要的浪费和碳排放。所以计及不确定性的本质是让调度方案在面对多种可能的出力场景时要么保守到所有场景都能扛住要么在某些概率约束下保证可行而不是只针对一个单点预测值优化。这也是这个项目为什么把不确定性建模放到第一优先级的原因。1.2 综合能源系统协同优化的协同体现在哪这个项目说的电气设备综合能源系统不是单一的电网络而是包含电、气、热三种能量流的小型综合能源系统。典型的结构是风电、光伏作为清洁电源燃气轮机或热电联产机组提供电和热电锅炉或热泵把多余的电力转化成热能储能电池平衡短期功率波动储热罐配合热力系统平移供热压力再加上电转气设备把富余电力转化成天然气存入储气罐实现跨时段、跨能源的耦合。协同的含义就很清楚了电能、天然气、热能不是各管各的而是通过设备耦合在一起整体调度。比如夜里风电出力大、电价低可以用电锅炉多产热存进储热罐白天热负荷高峰再释放如果天然气价格便宜就让燃气轮机多发电顺便供热减少从电网买电。这些决策是互相牵制的——多用电锅炉就会减少燃气轮机的发电需求多用电转气又会反过来影响锅炉的用气量。如果只做单能源优化比如只优化电力系统那热负荷怎么满足就管不着热力系统的灵活性完全浪费了。综合能源协同优化的价值就是把这些耦合关系用一个统一的数学规划问题表达出来让电、气、热三个系统联合运行整体运行成本最低、弃风弃光最少、系统越限风险最小。1.3 优化模型总体框架整个模型可以概括为以系统总运行成本最小为目标在设备运行约束、能源平衡约束、网络约束和不确定性约束的共同限制下决定各设备各时段的出力计划。具体到数学形式目标函数通常是从上级电网购电成本天然气购买成本设备运行维护成本弃风弃光惩罚成本切负荷惩罚成本这几个部分加起来最小化。约束包括每个时段的电功率平衡风电光伏燃气轮机储能放电电网购电 电负荷电锅炉耗电电转气耗电储能充电热功率平衡燃气轮机余热电锅炉产热储热罐放热 热负荷储热罐充热以及气网平衡天然气网购气电转气产气 燃气轮机用气热负荷用气等。为了计及新能源不确定性需要把风电、光伏出力从固定值改成随机变量再通过场景法或鲁棒区间等思路把随机优化转化为可解的确定性等价问题。我下面会分别讲两种主流做法再给出Matlab代码实现。2. 建模细节从不确定性表达到设备模型2.1 不确定性建模场景法 vs 鲁棒优化 vs 机会约束处理新能源出力不确定性目前主流有三条技术路线。场景法是直观也最容易上手的方法。思路是生成若干组风电、光伏出力的可能场景每组都是一个完整时段的出力曲线然后用一个含概率的期望值目标函数替代原来的单场景目标。约束条件要求每组场景下系统的运行约束都要满足最后一个多场景优化的解对各个场景都有一定的适应性。场景法能精细刻画随机分布也能兼容时序相关性缺点是场景数量多时计算量成倍上涨。鲁棒优化则是把不确定性描述成一个区间或集合要求调度方案对集合内最恶劣的场景都可行。鲁棒优化不需要概率分布只要知道出力的上下边界求解结果偏保守但工程上很可靠。适合那些安全要求极高、宁愿多花钱也不能出事的情况。机会约束采用了折中思路允许约束以一定的概率违反。比如要求系统功率平衡被打破的概率不超过5%把硬约束变成概率约束然后再用给定分布的分位数或采样法把概率约束转化为确定性约束。机会约束可以在保守性和经济性之间取一个平衡点参数设置更灵活但实际求解需要一点技巧常常要配合场景削减或者解析近似。我在代码里选择的是场景法加机会约束的混合思路主体调度用场景法生成典型场景同时把可容忍的失负荷概率用机会约束表达再用场景削减控制计算量。这是目前科研和工程里兼顾效果和可实现性的方案。2.2 关键设备模型与运行约束综合能源系统里的设备每个都是一个带约束的数学模型。建模必须准确否则优化结果根本没法用。燃气轮机CHP输入天然气输出电功率和热功率。模型可以用一个热电比电/热输出之比来简化约束包括出力上下限、爬坡约束、最小启停时间。代码里通常用线性化约束把热电比处理成可变区间。电锅炉电能转化为热能转换效率基本恒定比如0.9。主要约束是耗电功率上限以及功率变化爬坡限制。储能电池状态变量是荷电状态SOC充放电效率不对称充0.9、放0.95还要避免同时充放。需要添加SOC连续性约束以及每个时段的充放电功率上下限。储热罐与储气罐与储能电池类似有容量上下限和充放能速率约束同时必须周期性约束比如一个调度周期结束时储能量要回到初始值否则优化会透支储能。电转气P2G消耗电能和二氧化碳或水蒸气反应生成天然气简化模型就是电能转化为气能有一个转换效率系数。约束主要是产气功率上限。设备模型的每一个参数都不能乱设需要根据实际设备手册或参考典型值。代码中我统一用标幺值建模方便调参。2.3 目标函数与约束条件的Matlab实现Matlab下做优化建模首选组合是Yalmip工具箱加外部求解器Gurobi或Cplex。Yalmip是个建模层把目标函数和约束写成接近数学表达式的形式然后交给求解器去算能省掉大量手写矩阵转换的功夫。以目标函数为例核心代码如下%% 目标函数总运行成本最小 Objective 0; for t 1:T % 购电成本 Objective Objective price_e(t) * P_buy(t) * dt; % 购气成本 Objective Objective price_g(t) * V_gas(t) * dt; % 运维成本简化按出力的比例系数 Objective Objective k_wt * P_wt(t) k_pv * P_pv(t) ... k_chp * P_chp(t) k_gb * Q_gb(t); % 弃风弃光惩罚 Objective Objective pen_cur * (P_wt_avail(t) - P_wt(t)) ... pen_cur * (P_pv_avail(t) - P_pv(t)); end % 若考虑多场景则对场景取期望Objective sum(scen_prob * Objective_scen);约束条件用Yalmip写起来也很直观。比如电功率平衡约束Constraints []; for t 1:T Constraints [Constraints, ... P_wt(t) P_pv(t) P_chp(t) P_es_d(t) P_buy(t) ... P_load(t) P_eb(t) P_p2g(t) P_es_c(t)]; end这里P_wt、P_pv如果是场景变量就要在场景循环里逐个定义。我建议先用单场景把模型调通再扩展到多场景调试效率会高很多。3. 代码实现与求解过程从数据生成到结果输出3.1 数据准备与场景生成蒙特卡洛/拉丁超立方场景生成是整个随机优化里最吃数据处理功底的部分。我用的方法是先对风电、光伏的预测误差构造随机分布再用蒙特卡洛抽样生成大量场景最后用场景削减算法选出一组代表性场景并把每个场景的概率算出来。风电出力我通常用Weibull分布描述风速再通过风速-功率转换曲线得到出力光伏出力则用Beta分布拟合光照强度乘以光电转换效率得到功率。预测误差可以设为正态分布均值0标准差取预测值的10%-20%。生成场景的核心代码%% 生成风速与光照场景示意 rng(2024); % 固定随机种子保证可复现 N_scen_total 500; % 总抽样场景数 % 风速抽样Weibull分布 wind_speed wblrnd(k_shape, c_scale, N_scen_total, T); % 光伏光照强度抽样Beta分布归一化 irradiance betarnd(a_beta, b_beta, N_scen_total, T); % 通过功率曲线得到出力场景 P_wt_scen P_wt_rated .* (wind_speed cut_in) .* ... ((wind_speed.^3 ./ (rated_speed^3 - cut_in^3)) ... - (cut_in^3 ./ (rated_speed^3 - cut_in^3))); P_pv_scen eta_pv * S_pv .* irradiance .* solar_watch;直接用500个场景求解会非常慢所以要场景削减。我用的同步回代削减SBR方法原理是迭代地删除与其他场景距离最近的场景并把被删场景的概率累加到最近邻的场景上。削减后的场景能保留原始分布的主要特征。削减后的场景数一般取10-20个是我试下来精度和计算速度的平衡点。如果场景太少不确定性信息丢失严重太多求解时间指数上升。3.2 基于Yalmip的模型搭建核心代码有了场景数据多场景协同优化模型的Yalmip代码结构可以这样组织%% 定义决策变量按场景 for s 1:N_scen for t 1:T % 场景s下t时刻的决策变量 P_wt{s}(t) sdpvar(1,1); % 实际风电出力≤场景可用出力 P_pv{s}(t) sdpvar(1,1); P_chp{s}(t) sdpvar(1,1); Q_chp{s}(t) sdpvar(1,1); P_eb{s}(t) sdpvar(1,1); Q_eb{s}(t) sdpvar(1,1); P_es_c{s}(t) sdpvar(1,1); P_es_d{s}(t) sdpvar(1,1); SOC_es{s}(t1) sdpvar(1,1); P_p2g{s}(t) sdpvar(1,1); V_gas_buy{t} sdpvar(1,1); % 购气决策可不分场景第一阶段决策 end end注意一个优化建模的重要细节哪些决策变量需要在实时运行前就确定第一阶段决策比如购气合同量哪些可以等不确定性实现后再调整第二阶段决策比如储能出力。在我的模型里把储能、燃气轮机出力这类可以快速响应的变量设为场景相关把与上级电网的购电、购气这类需要提前安排的设为一阶段共享变量。这种两阶段决策结构更贴近实际的日前调度加日内调整过程代码写起来会多一层变量索引但结果明显更合理。约束的循环添加方式Constraints []; for s 1:N_scen for t 1:T % 场景约束 Constraints [Constraints, ... P_wt{s}(t) P_wt_scen(s,t), ... P_pv{s}(t) P_pv_scen(s,t), ... P_chp{s}(t) P_chp_max, ... Q_chp{s}(t) eta_hp * P_chp{s}(t), ... P_es_d{s}(t) - P_es_c{s}(t) P_es_rated, ... SOC_es{s}(t1) SOC_es{s}(t) eta_c*P_es_c{s}(t)*dt ... - (1/eta_d)*P_es_d{s}(t)*dt, ... % 电功率平衡 P_wt{s}(t) P_pv{s}(t) P_chp{s}(t) ... P_es_d{s}(t) P_buy(t) ... P_load{t} P_eb{s}(t) P_es_c{s}(t) P_p2g{s}(t)]; end end这里有一个容易踩的坑Yalmip的约束拼接会随着循环次数增加变慢如果模型规模大建议先初始化空的约束单元数组再用中括号合并。另外务必避免把同一组约束重复添加两遍否则约束矩阵的维度会翻倍求解慢且容易内存溢出。3.3 求解器选择与参数调优模型编好之后选择合适的求解器能让效果天差地别。我的模型是混合整数线性规划MILP因为要处理燃气轮机的启停状态、储能充放状态的0-1变量所以不能只用线性规划求解器。在Matlab环境中推荐的组合是Gurobi性能最强MILP求解速度业界标杆适合中大规模综合能源系统模型。Cplex老牌商业求解器稳定但近年更新相对慢。intlinprogMatlab自带MILP求解器小模型够用大模型会非常慢。SCIP开源求解器免费但速度不如前两者。调用方式%% 求解设置 ops sdpsettings(solver,gurobi, ... gurobi.TimeLimit,3600, ... gurobi.MIPGap,0.01, ... gurobi.Presolve,1, ... gurobi.NumericFocus,2); optimize(Constraints, Objective, ops);参数调优的核心是MIPGap最优间隙和TimeLimit求解时间上限。工程上不需要把MIPGap设到0因为综合能源模型本身数据精度有限设到1%或0.5%已经能保证结果可靠性还能大幅缩短计算时间。我试过将一个500个场景的模型强行求解到0.1%的MIPGap跑了四个小时还没收敛改成1%后十分钟出结果目标函数值只差不到0.8%工程上完全可接受。NumericFocus调成2是提醒Gurobi加强数值处理的苛刻程度能减少因为约束数量级差异造成数值问题。尤其是储能SOC约束里涉及到充放电效率和时间步长数值可能跨好几个数量级不加这个参数容易误报不可行。3.4 结果可视化与敏感性分析求解完成后至少要输出三类结果设备出力曲线、系统能量平衡情况、成本构成。我习惯用几个子图放在同一张figure里figure; subplot(3,1,1); stairs(t, P_wt_result, g); hold on; stairs(t, P_pv_result,b); stairs(t, P_chp_result,r); stairs(t, P_load,k--); legend(风电,光伏,燃气轮机,电负荷,Location,best); title(电力平衡); subplot(3,1,2); stairs(t, Q_eb_result,m); hold on; stairs(t, Q_chp_result,r); stairs(t, Q_heat_load,k--); legend(电锅炉供热,燃气轮机余热,热负荷); title(热力平衡); subplot(3,1,3); stairs(t, SOC_es_result); hold on; stairs(t, P_p2g_result); legend(储能SOC,电转气功率); title(储能与P2G);敏感性分析也是论文和工作汇报里很关键的一环。我做了三组敏感性测试一是新能源渗透率从20%提到60%看系统运行成本和弃风率的曲线二是储能容量从0.5MWh到2MWh看整体成本变化三是场景数量从5个到30个看求解时间和目标函数差异的变化趋势。这三类结果能直观说明模型的鲁棒性和协同优化的价值。4. 实操中的常见问题与避坑指南4.1 场景数怎么选计算量爆炸怎么办场景法是场景越多越准确这个直觉是不完全对的。当场景数超过一定阈值求得的期望目标函数变化很小但求解时间却成倍增加。我建议先用10个场景跑通流程然后做消融实验分别用10、15、20、30个场景求解对比目标函数值。如果15个场景与30个场景的结果差异不到1%选15个就够了。计算量爆炸的另一个常见原因是约束冗余。比如储能SOC约束如果每个场景、每个时刻都重复定义同一个初始SOC连续性条件会产生大量无效变量。好的做法是在场景循环外定义一个基础的储能初始状态场景内只定义增量变化。此外削减场景时要注意不要为了追求均匀覆盖而把所有极端场景都保留那样模型会被少数极端场景主导结果过于保守。4.2 求解器报错与收敛性调优我遇到最多的报错是Infeasible problem不可行。这通常不是模型本身无解而是约束写错了。排查步骤先检查平衡约束左右两边的符号是否反了。我早期把电平衡约束写成右边减去左边等于0因为符号问题排查了很久。检查储能SOC初始值设置。如果SOC初始值是0.5但后续约束要求连续运行24小时后SOC回到0.5而储能容量和充放电功率限制不能同时满足就可能无解。检查设备出力上限和功率平衡是否自相矛盾。比如热负荷很大但燃气轮机最大余热加电锅炉最大产热仍然不够模型就会不可行。这时需要允许一定比例的切负荷惩罚让模型在极端情况下可以牺牲一点负荷而不是崩溃。收敛性不好表现为Gap一直降不下去可能有几个原因。一是模型本身存在对称性比如两个相同的储能设备优化结果可能反复切换造成分支定界效率极低可以加一个简单的0-1变量排序约束来打破对称二是MIPGap和TimeLimit的配合问题可以通过设置求解器输出诊断信息用ops.gurobi.OutputFlag 1查看分支情况。4.3 不确定性参数分布设错导致优化结果离谱场景生成的分布参数直接影响优化结果的合理性。比如风电出力我一开始直接套用Weibull分布的形状参数k2、尺度参数c6生成的风速序列经常出现风电机组满发导致优化结果严重依赖风电但实际场景中很少连续满发。后来发现不同地区的风速分布差异很大需要先用历史数据做参数拟合。用Matlab的fitdist函数可以方便地估计Weibull分布的参数不要凭经验拍脑袋。光伏的Beta分布参数同样要基于实际辐照数据拟合。如果不做拟合可以退一步用历史数据分时段的均值与方差代替理论分布用经验分布抽样效果也不错。另一个坑是预测误差的方差设置。方差太小生成场景几乎重合等于没加不确定性方差太大场景离群严重优化结果过度保守。我建议误差标准差设为预测值的15%左右然后根据敏感性分析适当调整。4.4 模型验证与代码复现的经验做优化模型最怕的是模型写对了但磨出来的方案无人验证。我的验证方法是先让模型在确定性条件下即去除不确定性只用预测值求解把得到的结果与已知的实际运行数据对比看购电、购气量是否在合理范围。确定性场景验证通过后再叠加不确定性对比不同场景下的成本差异确认不确定性约束真的在发挥作用。复现别人的代码时如果发现结果对不上先检查单位。综合能源领域最容易出问题的是功率MW与能量MWh的转换时间步长dt如果是小时功率乘以dt才是能量但不少代码里直接不乘dt导致储能SOC变化量偏大。我习惯在代码注释里写清楚每个变量的单位并统一使用小时作为时间单位。另外数据文件的格式必须严格对齐。场景生成脚本产生的矩阵维度是场景数×时段数如果数据读取时行列没转置就会导致循环索引越界排查起来相当耗时间。我建议在处理完数据后立即验证size()和assert()条件杜绝这类隐藏bug。5. 后续扩展方向与个人心得这个模型跑通之后后续可以扩展的空间很大。比如加入需求响应机制让部分电负荷、热负荷可以根据电价和热价灵活平移或者引入碳交易成本把碳排放配额和碳价纳入目标函数这样协同优化就不只是经济驱动还能兼顾低碳。更复杂一点可以在模型里加入配电网潮流约束用DistFlow模型细化电压和线路容量限制让调度结果更贴合真实电网。我个人在实际操作中的一个体会是综合能源协同优化最大的难点不是算法有多深而是怎么把不同能源设备的物理耦合关系梳理清楚。我总是先从简单的电-热耦合开始搭建跑通基础场景再逐步加入气网和P2G每新增一种设备就重新做一次单场景验证。这样做的好处是每个耦合关系都能被单独检验一旦出现问题能快速定位到是哪一类设备约束出了问题而不是在整个复杂模型里大海捞针。最后再分享一个小技巧在Matlab中调试多场景模型时可以把场景参数和求解结果以结构体形式保存到本地mat文件方便多次对比调试。我习惯每次求解完都存一个带时间戳的快照如果之后改了模型导致结果变得不合理可以直接对比新旧结果找到是哪一处约束或参数变化引起了问题。这套调试习惯在工程项目里能省下大量重复排查时间。