做综合能源系统调度最让人头疼的不是模型本身有多复杂而是“预测永远不准”。光伏和风电出力跟着天气走电负荷、热负荷跟着人的作息走你拿早上的预测数据做完全天优化下午就可能发现原计划已经和实际脱节。所以行业内基本形成了共识先做日前计划日内再滚动修正两阶段各管一段。这个项目就是用Matlab和Yalmip把完整的“日前-日内两阶段综合能源系统调度”跑通包含模型搭建、代码实现、仿真分析和常见坑位整理。这套方案适合两类人一是做园区能源管理、微电网调度、综合能源优化控制的研究生和工程师二是刚接触Yalmip想知道新能源消纳、储能调度、热电联产约束怎么落到代码里的人。文章里所有参数、目标函数、约束和代码都可以直接照着自己改不用从零开始啃文献。1. 项目背景与整体思路为什么非要做两阶段调度先聊清楚为什么不能只做一次调度。综合能源系统里光伏、风电这类可再生能源的预测误差随预测时域缩短而快速减小。日前预测提前24小时光伏出力的预测误差可能到20%以上但到了当天超短期预测误差可以压到5%以内。如果只靠日前计划一旦出现云层遮光、负荷突变系统就得临时切负荷或者弃电经济性和可靠性都受影响。如果完全放弃日前计划只做实时优化又容易陷入“只见树木不见森林”——不同时段之间缺少经济协调储能总是没有跨时段套利的意识整体成本反而更高。1.1 单阶段调度为什么不够用单阶段调度通常指只基于某个固定预测场景做一次全局优化。这种方法的优点是简单缺点是缺乏修正能力。我最早试过只做日前优化把24小时计划下发之后第二天实际光伏出力比预测低了30%结果中午时段只能大量购电而当天峰时电价恰好最贵一个午间时段就吃掉了一整天的成本结余。有人可能会说那把预测做准一点行不行问题在于中长期气象预报的天然不确定性就摆在那里投入再多算力也消不掉。更合理的做法不是拼命提升预测精度而是设计一个能容忍误差、能在运行中逐步修正的调度框架。这也是两阶段调度的核心动机。1.2 两阶段分工预测精度与经济性的权衡两阶段调度把决策拆成两个时间尺度。日前阶段以小时为间隔基于未来24小时的日前预测数据以系统运行成本最小为目标确定机组出力计划、储能充放策略、购售电计划。这个阶段解决的是“经济上最优的方向在哪”。日内阶段则是滚动修正通常以15分钟或5分钟为间隔把预测时域缩短到未来几小时基于更新的超短期预测和实时量测数据以“尽量跟踪日前计划、同时最小化调整成本”为目标重新求解各设备的出力增量。这个阶段解决的是“实际运行怎么应对不确定因素”。两阶段的好处在于经济性和鲁棒性不再互斥。日前阶段充分挖掘了峰谷价差、储能转移、热电联供的好处日内阶段又把预测偏差逐步消化掉而且因为修正周期很短偏差被限制在一个小范围内调整代价远小于推倒重来。1.3 为什么是MatlabYalmip调度模型本质上是优化问题而且是带约束的线性/混合整数规划。Matlab的优势在于数据处理和可视化成熟做负荷曲线、光伏曲线、储能SOC的曲线绘制非常顺手。Yalmip的作用则是把优化建模的复杂度降下来你不需要手动把约束写到求解器的标准矩阵形式里只要用数学表达式直接声明变量、写约束、定义目标就能交给Gurobi、CPLEX、Mosek这些商业求解器去解。实际工程里Gurobi和CPLEX对线性规划和混合整数规划的求解性能远好于Matlab自带的linprog/intlinprog尤其是MILP规模一上来差距非常明显。Yalmip就是在Matlab和求解器之间做翻译层既保留了Matlab的灵活又拿到了商业求解器的性能。这也是这个项目选择这套技术栈的根本原因。2. 调度模型目标函数、约束与热电耦合处理模型是整套调度的核心。我在项目里建立的综合能源系统包含光伏、风电、燃气轮机热电联产、燃气锅炉、电储能、热储能以及与电网的购售电交互。负荷侧考虑电负荷和热负荷系统通过电母线和热母线分别实现能量平衡。2.1 系统拓扑与能量流系统按母线分成两部分。电母线连接光伏、风电、燃气轮机、电储能、电网和电负荷热母线连接燃气轮机余热、燃气锅炉、热储能和热负荷。燃气轮机是热电耦合点发电的同时产生余热通过换热装置供给热母线。这种结构在商业园区、学校园区、工业厂区都非常典型。建模时最关键的一点是别把电和热割裂开燃气轮机的电出力和热出力存在耦合约束不能分别单独决定。后面我会详细讲热电耦合怎么处理。2.2 日前调度目标函数日前阶段的目标函数是系统总运行成本最小包含四项燃气轮机燃料成本、燃气锅炉燃料成本、电网购电成本、弃风弃光惩罚成本。售电收入在目标函数中表现为负成本也就是卖出电能抵减购电支出。目标函数写成$C_{total} \sum_{t1}^{24} \left[ C_{gas}\cdot (F_{gt,t} F_{gb,t}) price_{buy,t}\cdot P_{buy,t} - price_{sell,t}\cdot P_{sell,t} \lambda_{curtail}\cdot (P_{pv}^{fore} P_{wt}^{fore} - P_{pv,t} - P_{wt,t}) \right]$其中燃料耗量根据出力和效率折算。燃气轮机综合效率取35%余热回收效率按80%估算燃气锅炉效率取90%。天然气单价3.0元/m³低位热值9.7 kWh/m³。购电价格采用分时电价峰谷价差明显储能就可以通过谷充峰放来套利。注意弃风弃光惩罚项的作用是让优化倾向多消纳新能源但惩罚系数不能设太高或太低否则会扭曲经济信号。我这里取0.5元/kWh略低于购电均价、高于上网电价逻辑是优先消纳本地新能源其次才是购电。2.3 日内滚动调度目标函数日内阶段目标函数和日前不完全一样。理论上如果完全沿用日前目标滚动求解时可能出现“每次都在重新优化全局”导致出力计划在时段之间剧烈跳变。我在日内阶段采用“偏差最小化为主、运行成本微调为辅”的目标结构$C_{intra} \sum_{j1}^{N_h} \left[ \lambda_{dev}\cdot (\Delta P_{gt,j}^2 \Delta P_{ch,j}^2) \text{燃料成本增量} \right]$其中 $\Delta P$ 是各设备出力相对日前计划的偏差量。二次项是为了让偏差惩罚趋近于小偏差避免出现大范围调整。实际求解时如果是MILP二次项可以用线性化代替或者直接用一次项配合分段权重。核心思想是日前计划是“经济锚点”日内只做必要修正。比如光伏突然少发30kW日内就优先让电储能多放电、燃气轮机轻微爬坡而不是重新安排全天的储放策略。2.4 核心约束条件约束按性质分四组。功率平衡约束是硬约束电和热分别平衡设备容量约束限制各出力变量的上下限爬坡约束限制燃气轮机单位时间内的出力变化幅度储能约束包括SOC递推关系、容量范围和充放功率限制。电功率平衡$P_{gt,t} P_{pv,t} P_{wt,t} P_{buy,t} P_{dis,t} - P_{ch,t} - P_{sell,t} P_{load,t}$热功率平衡$H_{gt,t} H_{gb,t} H_{dis,t} - H_{ch,t} H_{load,t}$储能SOC递推$SOC_{e,t} SOC_{e,t-1} \eta_{ch}\cdot P_{ch,t}\cdot \Delta t - P_{dis,t}\cdot \Delta t / \eta_{dis}$这里要特别注意效率放在充电还是放电侧。我习惯把充电效率乘在充电项把放电效率放在分母上这样SOC模型更符合实际但代码里一定保证单位统一否则会出现SOC虚增。2.5 热电耦合与非线性的处理思路燃气轮机的热电耦合是综合能源系统调度里最需要注意的建模点。工程上常用“热电比”近似描述$H_{gt,t} r_{th}\cdot P_{gt,t}$。我取热电比1.2意味着发100kW电的同时输出120kW热。这个关系让电出力和热出力绑定在一起。如果热电比固定模型是线性的非常好解。但实际微型燃气轮机存在最小技术出力限制在低负荷率时热电比会变化。更精细的做法是用可行域多边形约束描述把若干个离散工况点线性组合成凸包这样电和热的可行组合区域更真实。我项目里为了先跑通框架采用了固定热电比方案同时在代码里保留了可行域多边形的注释方便后续扩展。燃料成本函数严格来说也不是线性的燃气轮机在部分负荷率下效率下降燃料成本是出力的二次函数。如果直接用二次函数模型变成MIQPGurobi也能解但求解时间明显增加。工业界更常用分段线性化把出力区间切成几段每段用线性函数逼近模型转成MILP求解快且精度足够。3. MatlabYalmip建模实现全流程模型想清楚之后代码实现反而是相对机械的工作。但Yalmip建模有几个细节如果不注意后面排查会非常痛苦最典型的就是变量维度不匹配和约束拼装错误。这一节把我踩过的坑按流程梳理一遍。3.1 环境准备Yalmip与求解器安装配置Yalmip本质是一个Matlab工具箱下载后把文件夹放到Matlab路径里即可。注意要用addpath把Yalmip目录和求解器目录都加进去否则调用optimize时提示找不到求解器。addpath(D:\Tools\yalmip); addpath(D:\Tools\gurobi\win64);求解器方面Gurobi和CPLEX对学生都有免费学术许可学校邮箱申请即可。如果只是小规模测试Matlab自带的linprog也能跑但MILP求解能力较弱不建议在正式项目里依赖它。安装Gurobi后要在命令行执行一次授权激活然后在Matlab里用gurobi_setup完成对接。一个容易忽略的点Yalmip在调用求解器前会先做预处理如果模型里存在无界变量或约束矩阵稀疏度过高求解器可能直接报“QCP not supported”之类的错误这时不要急着怀疑求解器坏了先回头检查建模。3.2 变量定义与维数管理Yalmip用sdpvar定义决策变量。我的习惯是按照“设备类别时间维度”批量定义。以24时段日前调度为例变量定义为P_gt sdpvar(24, 1, full); % 燃气轮机电出力 H_gt sdpvar(24, 1, full); % 燃气轮机热出力 P_buy sdpvar(24, 1, full); % 购电 P_sell sdpvar(24, 1, full); % 售电 P_ch sdpvar(24, 1, full); % 电储能充电 P_dis sdpvar(24, 1, full); % 电储能放电 SOC_e sdpvar(24, 1, full); % 电储能SOCsdpvar的第三个参数默认是full表示普通实数变量。如果定义矩阵变量时省略这个参数Yalmip默认也会按全矩阵处理但在二次规划里如果漏写symmetric会导致约束不对称所以养成显式写参数的习惯比较好。维数管理的核心原则是所有 24×1 的变量在拼装约束时左右两边必须同样是24×1否则Matlab会做隐式扩展或直接报错。我最常犯的错误是电价向量写成了1×24和24×1的购电变量逐项相乘前没有转置结果目标函数里出现维度错误。3.3 约束和目标函数的代码化约束拼装用方括号一直累加即可。以电功率平衡、爬坡约束和储能初值为例Constraints []; % 电功率平衡每时段都满足 Constraints [Constraints, ... P_gt P_pv P_wt P_buy - P_sell P_dis - P_ch P_load]; % 燃气轮机爬坡约束 for t 2:24 Constraints [Constraints, ... -ramp_rate P_gt(t) - P_gt(t-1) ramp_rate]; end % 储能初值和SOC递推 Constraints [Constraints, SOC_e(1) 0.5]; for t 2:24 Constraints [Constraints, ... SOC_e(t) SOC_e(t-1) eta_ch * P_ch(t) * dt - P_dis(t) * dt / eta_dis]; end目标函数我习惯先算成本矩阵再做向量点乘和求和。这样逻辑清晰也方便后续改价格参数C_fuel c_gas / LHV .* ( P_gt / eta_gt_elec H_gt / eta_gt_heat H_gb / eta_gb ); C_buy price_buy .* P_buy; C_sell price_sell .* P_sell; C_curtail lambda_curtail * ( (P_pv_fore P_wt_fore) - (P_pv P_wt) ); Objective sum( C_fuel C_buy - C_sell C_curtail );注意这里C_fuel是24×1向量sum对向量求和得到总成本。Gas提单算燃料成本时如果c_gas/LHV*(P_gt/eta_gt_elec)算出来的单位是元验证一下量纲就能避免很多低级错误。3.4 求解配置与结果后处理求解器设置使用sdpsettings。Gurobi下常用配置ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; ops.gurobi.TimeLimit 300; res optimize(Constraints, Objective, ops);verbose2会让求解器输出迭代日志方便判断是否找到最优解。如果res.problem返回0表示求解成功返回1表示求解器内报错返回15表示模型不可行。拿到结果后用value函数提取变量数值P_gt_opt value(P_gt); SOC_e_opt value(SOC_e);一个非常实用的细节提取结果后先做简单校验比如把电功率平衡约束左侧的数值代入计算残差。如果残差在1e-6量级以上多半是某些变量的单位或者索引没对齐这时候直接画曲线对比会非常容易暴露问题。3.5 日前计划与日内滚动的数据衔接两阶段调度在代码层面最关键的衔接是把日前结果作为日内优化的参考基准。我的做法是把日前决策结果保存为结构体日内滚动时逐时段读取对应值DA_plan.P_gt value(P_gt); DA_plan.SOC_e value(SOC_e); % 日内滚动每15分钟决策一次预测时域4小时共16个点 Horizon 16; for k 0:95 idx k 1 : k Horizon; ref_P_gt DA_plan.P_gt(ceil((k1)/4)) * ones(Horizon, 1); % 用超短期预测替换日前预测 P_pv_intra P_pv_ultra(idx); ... optimize(Con_intra, Obj_intra, ops_intra); % 取第一个时段的控制指令下发 u_next(:, k1) value(P_gt(1)); end这个循环结构是MPC模型预测控制的典型写法。每次只实施当前时段的指令下个时段重新求解。注意在滚动循环开始前把日前计划按照日内的时间分辨率插值否则15分钟分辨率的变量和1小时的参考值对不齐。4. 仿真算例日前-日内两阶段调度结果分析光有模型和代码还不够需要每个项目都有可复现的算例。我这边设定了一个典型冬季日场景工业与商业混合负荷参数都有依据方便你替换成自己的数据后对比。4.1 算例场景与参数表系统规模和设备参数如下设备参数数值单位光伏装机容量300kW风电装机容量150kW燃气轮机额定电出力200kW燃气轮机热电比1.2-燃气轮机发电效率0.35-燃气轮机爬坡速率60kW/h燃气锅炉额定热出力300kW电储能容量200kWh电储能最大充/放功率100kW电储能充/放效率0.95-热储能容量300kWh热储能最大充/放热功率150kW购电价格采用典型工商业分时电价。峰时段为10:00-15:00、18:00-22:00电价1.0元/kWh平时段为07:00-10:00、15:00-18:00电价0.6元/kWh谷时段为22:00-次日07:00电价0.35元/kWh。售电价格统一取0.35元/kWh。电负荷峰值约310kW出现在晚间热负荷峰值约350kW出现在清晨和傍晚。光伏预测曲线在中午达到300kW但日内实际出力在13:00-15:00出现明显下降模拟短时云层遮挡造成的预测偏差。4.2 日前调度计划解读求解后日前计划显示了非常清晰的运行规律。夜间电价低燃气轮机维持在较低出力水平电储能低谷充电到早晨SOC接近满值。上午光伏出力上升燃气轮机降出力电储能开始放电满足负荷爬坡。中午光伏大发电储能充电多余电量少量上网。傍晚光伏退出热负荷上升燃气轮机进入全天最高出力状态电储能放电支撑晚高峰。深夜电价再次进入谷值储能重新充电。这个结果的主要特点就是“峰谷联动”储能充放与电价变化严格对应燃气轮机的出力峰值也基本落在电价峰值和热负荷高峰重叠的时段。这说明目标函数和约束的设计是合理的经济信号传导到了设备出力上。4.3 日内修正一次典型预测偏差的消纳过程我在算例中人为构造了一个预测偏差场景日前预测13:00-15:00光伏出力分别为285kW、290kW、295kW而日内实际出力只有180kW、160kW、150kW相当于出现约40%的短时偏差。如果只按日前计划执行系统就会在这三个时段缺电约300kWh只能被迫高价购电或者切负荷。日内滚动调度在13:00时刻检测到光伏低出力后立即重新求解。结果显示电储能提前增加放电功率从原计划的30kW提升到70kW燃气轮机在爬坡速率允许的范围内每小时增加40kW出力同时削减电储能充电计划。三路响应叠加最终电功率平衡被精确满足且没有触发切负荷。这组结果说明日内修正机制是有效的。关键不在于每一路调整了多少而在于调整方向与日前经济方向保持一致储能优先放、气机适度爬没有出现破坏性的出力震荡。4.4 经济性结果对比为了量化两阶段调度的价值我做了三组对比调度方式运行成本元弃光弃风量kWh缺电风险纯日前计划开环执行14320120高日前日内两阶段1386515低理想预测全知调度135100无两阶段相比纯日前开环运行成本节省约3.2%弃光弃风量下降87.5%而且完全没有出现缺电时段。相比理想预测两阶段只高出约2.6%这个差距主要来自日内调整的燃料和储能损耗属于为鲁棒性付出的合理代价。从工程角度看这3%左右的成本节省非常可观。园区级系统一年运行成本几百万节省3%就是十几万级别的收益这也是综合能源调度项目能商业化落地的核心逻辑。5. 常见问题与排查经验实录这部分是我最想分享的。模型和算法原理看起来不难但实际跑起来总会遇到各种莫名其妙的报错、不可行、结果异常。下面的经验都是我反复踩过的坑。5.1 模型报不可行怎么快速定位Yalmip在不可行时会提示“infeasible problem”但你根本不知道是哪个约束出了问题。我推荐一套三步定位法。第一步把二进制变量全部固定为初始值或者把MILP松弛成LP看是不是整数变量导致的不可行。很多时候是储能充放互斥逻辑写错了导致某个时段同时要求充电和放电。第二步用check函数检查约束残差找出残差最大的约束。Yalmip里check(Constraints)会返回每个约束的残差接近0的说明约束被满足远离0的很可能就是问题来源。第三步逐条删除约束做二分查找没跑过的人觉得麻烦实际上在24小时模型里五分钟就能定位。最常见的原因有三个储能初值和边界条件冲突、SOC变量在滚动优化中被错误重置、某个时段负荷值明显高于所有设备最大出力之和。第三种情况最隐蔽我会在输入数据前先做一个简单的“最大供给能力曲线”检查。5.2 求解器与MILP规模控制两阶段调度里如果加入了机组启停状态或者储能充放互斥状态模型就是MILP求解时间会显著上升。Gurobi默认的MIPGap设置为1e-4时会比较慢实际工程中可以把MIPGap放宽到5e-3或者1e-3运行成本误差在0.1%以内但求解时间可能缩短一个数量级。另外控制变量数量的核心在于减少冗余的二进制变量。比如储能充放互斥用状态变量引入后模型规模可能成倍增长。我后来发现很多情况下储能SOC的递推约束已经能够自然抑制同时充放因为同时充放意味着无效循环目标函数会主动排斥。直接在SOC递推里加一个微小的充放互斥惩罚项比引入状态变量更省事。5.3 非线性约束线性化与收敛技巧我之前在燃料成本里直接用了二次函数Gurobi配Yalmip求解MIQP小规模24小时模型要跑好几分钟而且结果一波动画曲线就会出现毛刺。后来切成段线性化8段近似后误差低于1%求解时间降到十几秒。段线性化有两种实现方式一种是用二进制变量选择段模型变成MILP另一种是用特殊有序集Gurobi支持SOS2约束Yalmip中可以直接用sos2声明求解效率更高。我实测下来SOS2方式代码更简洁MIPGap收敛也更快。二次项还有一个替代方案就是改用绝对值惩罚并通过逻辑变量分段。对于日内偏差最小化的目标函数用一次项配合不同权重往往比二次项更稳定不容易让求解器陷入数值振荡。5.4 Matlab/Yalmip环境运行典型报错Yalmip本身对环境依赖较小但Matlab版本的兼容性问题确实存在。不同年份的Matlab对Matlab内置函数结果有细微差异比如sdpvar与optimize的交互在新版本里偶尔会提示“Unable to create a simplex”之类的内部错误这时候优先检查求解器是否是最新版再检查Yalmip是否需要更新。还有一个非常常见的坑Matlab许可证验证失败导致优化求解器崩溃。刚装完Matlab或者更换电脑后如果报错信息里带有“license manager error”字样一般是许可证文件与当前机器绑定信息不匹配。解决办法是重新激活许可证用学校或单位提供的正版授权确保Matlab和Gurobi的授权都能正常访问再回来跑模型。这个问题和模型本身无关但容易误导你以为是建模错误先把环境确认好再排查。5.5 结果异常排查锤有一次我跑出来的电储能SOC曲线在凌晨时段不断下降但对应时段电储能明明没有放电。最后发现是因为我用的电负荷数据是“系统净负荷”已经扣除了光伏出力但模型里我又把光伏单独放在平衡方程中等于重复计入了一次新能源电量。这是典型的“数据口径没对齐”造成的异常。排查结果异常时我的习惯是先画能量平衡曲线把“供能总和-负荷总和”画成一条折线。正常情况这条折线应该始终在0附近波动一旦出现明显偏移就能快速定位是某个时段的数据口径或约束出了问题再深入看设备级变量。6. 这套框架还能往哪些方向扩两阶段调度框架跑通之后扩展空间非常大。我在实际项目中已经探索过几个方向顺手整理出来供你做技术选型时参考。6.1 从园区级到虚拟电厂聚合在园区级框架中把所有可调设备建模为单一优化问题本质上是一种集中式调度。当系统延伸到虚拟电厂时多个分布式资源分属不同主体就需要引入分布式优化常见的是交替方向乘子法或分布式梯度投影。Yalmip建模部分可以继续复用但求解循环要自己写分布式算法中的本地子问题仍然适合用Yalmip表达。6.2 引入随机规划/鲁棒优化两阶段调度的下一步是处理更极端的不确定性。随机规划需要给光伏和负荷构建场景集每个场景权重不同模型规模会增大很多。鲁棒优化则要构建不确定集把约束写成对最坏情况的鲁棒形式。Yalmip本身提供了鲁棒优化的高级接口可以直接对带不确定参数的约束声明uncertain内部自动重构。这条路我还在实践但初体验是建模代码改动不大求解性能是主要瓶颈。6.3 与数据驱动方法结合最近不少人在研究把深度学习预测和优化调度串成一个大流程。思路是用神经网络先给出光伏和负荷预测的置信区间再把这些区间信息转化为调度模型的不确定集参数。Matlab的深度学习工具箱与Yalmip在同一环境下协同代码衔接非常自然。我之前试过用LSTM预测光伏短期波动再把预测误差分位数输入日内调度模型整体效果比单纯滚动优化又提升了一截。我在实际项目中最大的体会是调度系统的好坏取决于三件事模型是否说清楚了物理约束代码是否把数据和变量口径对整齐求解配置是否匹配了实际规模。这三点做到位方案离落地就不远了。后续你自己改场景时建议先把算例参数按真实设备数据校准一遍再逐步增加设备类型和不确性建模每一步都留出前后对比的曲线这样既能验证模型正确性也能让工程决策更有底气。