把“多时间尺度滚动优化”“双层调度”和“MATLAB代码”这几个词放在一起绝大多数第一反应是论文复现第二反应是头疼。我去年用这套思路做了一个多能互补微网的仿真项目前前后后改了三版代码从单层静态优化一路做到日前计划加日内滚动踩过不少坑也积累了不少可以直接抄的套路。这篇就把整个模型的实现逻辑、代码结构和调试经验完整拆开讲一遍主要适合正在做微网调度、综合能源系统优化相关课题的研究生也适合刚接触滚动优化想模仿落地的工程师。你看完以后至少能搞清楚两件事这套模型在算什么以及MATLAB代码里那些循环、变量、约束到底怎么拼起来。我尽量不堆公式但该有的约束一个都不会少。1. 先把这套模型拆开看到底调了什么为什么搞两层1.1 微网里的“多能源”和“多时间尺度”分别指什么多能源微网不是只有电而是电、热、气多种能源互相耦合的小型供能系统。典型配置是这样光伏板和风机负责清洁发电微型燃气轮机或热电联供机组同时出电和出热燃气锅炉专供热负荷电锅炉把电力转成热力蓄电池存电蓄热罐存热。这些设备通过一条低压配电网和一个区域热网连在一起对外还能跟上级电网买卖电力。时间尺度则对应调度决策的更新频率。传统调度喜欢做24小时日前计划步长取1小时提前一天把明天的机组出力和储能充放定下来。但真实运行里光照、负荷、电价全是波动的一天前定的计划到当天下午早就失真了。所以需要在当天向前滚动地重新优化比如每15分钟执行一次优化每次只看未来4小时执行完第1个步长后下个周期再刷新状态重新算。这种“边走边看、只执行眼前一段”的思想就是模型预测控制里的滚动优化拿到微网调度上完全成立。1.2 双层调度模型到底在“双”什么很多论文把调度分成日前层和日内层两层各有使命。日前层解决的是“明天整体怎么安排更经济”的问题它掌握的是全天的预测数据做全局寻优决定机组启停、储能日内的储放策略这些慢变量。但日前层不会随实际情况调整预测的误差会累积所以加了日内层。日内层跑在当天它把日前层给出的计划作为参考基准以更短的步长、更短的预测窗口实时修正各设备的出力。这里存在一个约束衔接关系日内层优化出的各时刻功率必须尽量贴近日前计划同时允许在一个偏差带内灵活调整。两层的关系可以简单理解成日前层定大方向日内层做精细化纠偏。代码上体现得很直白日前层把输出存档日内层读取这些存档值作为参考然后追加偏差惩罚项。1.3 滚动优化和一次性全局优化的差异代码上体现得很直白普通静态优化是给定T个时段的所有预测数据一次求得T个时段的完整计划表。滚动优化则把整个调度拆成很多轮每一轮只求解一个长度为H的窗口。窗口内部的约束结构和静态优化一模一样区别在于每轮都要重新读取当前状态量把第一轮的第一个控制量执行然后窗口整体向前平移一个步长。不要小看这个改动。循环一旦写起来原来一次性建模的约束矩阵全部要在循环里重建变量也要重新声明而且初始SOC、蓄热量、爬坡基准这些上一轮结束时的终值必须准确传递给下一轮。很多人的代码跑起来结果异常并不是数学模型错而是滚动循环的状态传递写错了。2. 建模是代码的地基MATLAB里怎么描述设备和约束2.1 时序数据怎么组织才不容易乱我习惯把所有外部输入数据先读进一个结构体或表格里。比如data.Pload是电负荷序列data.Hload是热负荷序列data.Ppv是光伏预测出力data.price是分时电价统一用interp1或直接按索引取值对齐到15分钟步长。如果你要对比不同时间尺度比如日前用96点、日内用15分钟最好建一个时间轴数组再用ismember定位索引避免靠猜索引数写死。数据的时间对齐是滚动优化的地基推荐把数据按15分钟一个点统一存储。日前层做96点优化日内层则使用每轮窗口内的16个或32个点。这样代码逻辑清晰后面改成5分钟步长也只是换数据密度的事不用动模型结构。2.2 设备模型写成约束矩阵时的核心套路设备建模不复杂核心就是把每个设备的输入输出关系变成等式或不等式约束。光伏和风电直接按预测数据给定蓄电池用一阶差分方程描述电量变化SOC(t1)SOC(t)-Pbat(t)×Δt/容量充电时Pbat为负放电时Pbat为正SOC约束在0.2到0.9之间充放功率各有上下限。热电联供机组要同时满足电出力和热出力耦合关系热出力等于电出力乘以热电比燃料成本由电出力反算燃料消耗量得出。燃气锅炉则是热功率乘以锅炉效率得到燃料量。蓄热罐和电池的结构类似只是时间常数更长、损耗系数更大。把这些公式翻译成MATLAB约束时最忌讳一行一个约束去写。正确做法是向量化决策变量定义成sdpvar(1,96)的行向量约束批量追加在同一个约束数组里。例如蓄电池SOC的递推关系可以直接写成soc(2:end) soc(1:end-1) - Pbat(1:end-1)*dt/Emax一条语句替代96条循环。2.3 目标函数怎么设计才贴合双层结构日前层的目标函数是全天运行成本最小包含购电费用、购气费用、弃光弃风惩罚、电池储能退化成本。这里的退化成本一般做一个很小的惩罚系数乘上充放电功率绝对值用于防止储能频繁充放。日内层的目标函数则在成本基础上增加一项日前计划偏差惩罚偏差越大罚的越多这样日内优化器会不由自主地向日前计划靠拢。偏差惩罚系数是个关键参数。设得太小日内层会彻底放飞自我日前计划形同虚设设得太大日内层几乎没有调节能力滚动优化失去意义。我实际操作下来这个系数设置为该时段购电电价的20%到50%比较合理具体还需结合预测误差水平调。2.4 约束的矩阵维度必须预先想清楚做双层代码时最常见的报错是维度不匹配。YALMIP里sdpvar(1,96)是一维变量和标量约束直接相加没问题但两个不同长度的向量约束相加就会报错。建议所有跟设备容量相关的约束都写成向量形式所有常值参数通过引用同一个ones系数矩阵进行扩展。另外设备爬坡约束写起来很绕因为上一时刻和当前时刻同时出现在一个约束里我一般写成-Ramp x(2:end) - x(1:end-1) Ramp这样整段约束一次生成从第二点开始循环避免越界。3. 代码落地的完整链路从日前计划到日内滚动3.1 日前调度层的关键代码结构先用YALMIP定义决策变量然后累加约束和目标最后调用求解器。核心结构大致如下%% 日前调度层 T 96; % 24h15min一个点 Pbat sdpvar(1, T); Pchp sdpvar(1, T); Hgb sdpvar(1, T); S sdpvar(1, T); % 电池SOC cons []; % 功率平衡约束 cons [cons, Ppv Pchp Pbat Pgrid_buy Pload Pgrid_sell Pelec_heat]; % 设备上下限 cons [cons, 0 Pchp Pchp_max]; cons [cons, -Pbat_dis_max Pbat Pbat_ch_max]; % SOC递推 cons [cons, S(2:end) S(1:end-1) - Pbat(1:end-1)*dt/Emax]; cons [cons, 0.2 S 0.9]; obj sum(grid_price .* Pgrid_buy * dt) ... sum(gas_price .* (Pchp/eta_chp Hgb/eta_gb) * dt) ... sum(penalty_curtail .* (Ppv_max - Ppv_used)); ops sdpsettings(solver, gurobi, verbose, 2); optimize(cons, obj, ops);日前层计算完成后把每个时段的电池SOC、CHP出力、锅炉出力、联络线功率都存下来作为日内层的参考。存储方式建议用MATLAB的struct或array字段名写清楚否则后面半个月再回来看代码你自己都不知道哪个变量是日前值。3.2 日内滚动优化层的循环逻辑日内层是典型的多循环结构。外层循环是滚动次数内层是在当前窗口内构建并求解一个规模较小的优化问题。代码如下%% 日内滚动优化 H 16; % 未来4小时15min间隔 N T_H; % 实际运行96个点 soc_current s0; % 当前SOC for k 1:N % 读取窗口内的负荷、光伏、电价预测 pload_win load_forecast(k:kH-1); ppv_win pv_forecast(k:kH-1); price_win price(k:kH-1); % 声明窗口决策变量 pbat_win sdpvar(1, H); pchp_win sdpvar(1, H); S_win sdpvar(1, H1); % 初值也占一个位置 % 约束与目标 cons []; cons [cons, S_win(1) soc_current]; cons [cons, S_win(2:end) S_win(1:end-1) - pbat_win*dt/Emax]; cons [cons, 0.2 S_win 0.9]; % 功率平衡约束类似日前层... obj sum(price_win .* pgrid_win * dt) ... sum(gas_price .* (pchp_win/eta_chp) * dt) ... bias_penalty * sum(abs(pchp_win - ref_pchp(k:kH-1))); optimize(cons, obj, ops); % 只执行第一个控制量 soc_current value(S_win(2)); % 更新状态 result_record(k) value(pbat_win(1)); % 记录实际执行值 end这个循环看起来简单但实践时有个很容易忽略的点S_win同时包含初值位置和后续状态位置初值来自外部而不是优化变量所以窗口内S的索引和pbat索引存在一个滑动的错位关系。如果直接把S_win(2:end)和S_win(1:end-1)做递推必须确保pbat_win长度是H而不是H1否则维数对不上。另外每轮循环开头都要重新声明sdpvar旧变量会残留在工作区因此在处理时最好用函数封装或及时清理不需要的变量。3.3 上下层之间如何传参考值和初始值两层衔接是代码中比较微妙的部分。日前层的参考值在传递给日内层时通常会遇到时间分辨率不一致的问题。比如日前层步长1小时日内层步长15分钟日前参考必须用repelem或插值扩展到96个点再按窗口截取。如果严格保持两层步长一致直接索引引用即可但很多论文为了展示不同时间尺度故意让两层步长不同这时就要谨慎处理坐标对齐。还有一种做法是日内层不直接引用日前层的绝对出力而是引用日前计划的松紧度。比如让日内层在该时刻的联络线功率允许在日前值的正负15%内浮动这个浮动带本身也作为约束写进模型。这种方法比单一的偏差惩罚更稳健因为偏差惩罚是软约束极端情况下优化器宁愿交罚款也不执行日前计划硬约束则能保证功率外特性不跑偏。3.4 求解器配置和版本适配的实话MATLAB下跑这类模型目前最顺手的组合是YALMIP加Gurobi。Gurobi的MIP求解能力很强适合含启停整数变量的双层模型。如果你的模型全部是连续线性规划直接用Gurobi或Cplex都行甚至MATLAB自带linprog也只牺牲一点速度。但一旦加了整数变量比如机组的启停状态linprog就无能为力了。我建议在sdpsettings里设置mipgap为1e-3或1e-4并且开启feasibilityTol到合适范围不然整数规划求解时间会非常久。另外新版MATLAB的optimoptions和旧版YALMIP之间偶尔会有求解器接口不兼容的问题。如果你用Gurobi 11以上最好把YALMIP升级到最新版否则经常出现“could not locate solver”的诡异报错。4. 调试现场我踩过的坑和排查思路4.1 优化器报无解最常出现在哪些环节无解是所有调度代码的噩梦。我在这个项目里遇到的无解绝大多数不是模型本身不可行而是约束内部冲突。典型场景有两个一个是SOC的初始值设得太低但当前时刻又要求高功率放电导致递推约束直接矛盾另一个是热电联供的电热耦合关系跟热负荷约束打架电负荷要求CHP多发电但热负荷已经饱和CHP又不能弃热于是找不到可行解。排查无解问题有一个笨而有效的办法将目标函数改成常数零逐条注释掉约束组看哪组约束注释后优化器变可行。一般是先去掉热力平衡约束再单独检查SOC递推。更快的办法是开启YALMIP的showprogress调试模式它会在求解前报告哪个约束维度异常。大多数无解问题最终都指向约束里的系数符号比如放电时Pbat符号定义反了检查一下充放功率的正负约定能省很多调试时间。4.2 滚动循环里SOC曲线断裂哪里的索引出了问题滚动优化的曲线断裂往往不是模型问题而是状态更新代码写错。我第一版代码里SOC更新直接写成soc_current value(S_win(1))忽略了窗口内已经执行完第一个动作后状态应当是第二个值。这样每轮都把SOC重置到了初始值运行结果自然是错误的。正确的逻辑是SOC按照递推方程更新到下一时刻的值即value(S_win(2))并把该值传给下一轮作初值。如果中间有多个能量存储设备每个都要分别维护各自的当前状态变量不能偷懒共用同一个变量名。还有一类索引问题是参考数组的错位。比如日前计划的参考值数组是1到96而日内循环从第k步开始取值时如果直接用ref_pchp(k:kH)会导致循环末尾越界。我的做法是事先将参考数组预留H个点或者用min(kH-1, N)做裁剪并在裁剪后单独处理不足段直接终止下一轮。这个细节卡了我一晚上后来用断点观察索引才发现。4.3 滚动窗口设多长计算时间怎么平衡窗口长度直接决定单轮优化规模。窗口长了滚动次数少了但单次求解变慢窗口短了求解快但每轮接入的新预测信息少优化视野太短容易跑出局部不合理方案。经过我反复测试对15分钟步长来说4到8小时的预测窗口是最实用的区间。再长的话求解时间翻倍优化效果提升不明显而且预测误差本身也大长窗口意义有限。如果你嫌求解速度慢最简单的加速手段是放宽MIPGAP。默认1e-4和1e-2之间的求解时间可能差十倍而调度问题一般不需要精确到万分之一所以我会调节到1e-3。另外我还会主动把一些对整体目标影响极小的整数变量松弛成连续变量比如某些启停频繁但对成本影响低于0.1%的小设备直接允许部分负荷运行速度提升立竿见影。4.4 高频问题的速查表表现常见原因排查方法解决方案Gurobi报“Infeasible”约束冲突或初值不匹配逐条注释约束定位检查SOC初值与功率上下限调整约束范围求解时间过长整数变量过多看求解日志里MIP gap下降速度放宽mipgap松弛非关键整数变量滚动结果与日前计划差异巨大偏差惩罚系数设太小单独把偏差项值打印出来提高偏差惩罚或改偏差带为硬约束SOC曲线在窗口边界跳变状态更新取错索引打印每轮SOC初值与终值修正为取S_win(2)作为下轮初值运行速度越来越慢循环内累计了历史变量查看工作区变量数量每轮用函数封装或clearvars指定变量光伏出力被过度弃用弃光惩罚系数低于购电成本检查目标项符号和数值提高弃光惩罚确保不弃光更经济5. 想基于这套代码做二次开发和论文建议重点关注这几处5.1 从确定性调度升级到不确定性处理的替换思路当前模型基于预测值做确定性优化论文里最常见的扩展方向是加入不确定性。一个简单替换是把光伏和负荷预测改成一个区间范围然后在约束中增加鲁棒对偶约束。另一种替换是场景法生成若干个预测场景把目标函数改成各场景期望成本的加权和但这会成倍增加变量数量和约束结构。如果你想快速验证鲁棒模型而不是从头写我的建议是保留现在两层结构只把日内层目标函数中的单一预测值替换为场景平均值并给约束增加一个备用容量项。这样改动量小代码主体不变但仿真效果上能体现鲁棒性。等接受后再慢慢改造成标准鲁棒优化格式。5.2 代码结构上的可扩展性调整很多人写完一版代码后换一个网络拓扑就要改动大量引索非常痛苦。我后来吸取教训把设备参数全部塞进一个equip结构体例如equip.chp.Pmax、equip.bat.Ecap。这样不同的算例只需要复制一个基础配置文件设备数目变化时改动也局限在数据初始化和约束生成函数内部。更进一步可以把日前层和日内层分别封装成函数输入是数据结构和参考计划输出是调度结果和运行日志这样后期做敏感性分析、画对比图都不需要改动核心代码。5.3 结果可视化时最值得画的几张图论文配图里最常见的组合是三类图调度结果堆叠图、SOC和蓄热变化曲线、日前与日内出力的对比图。堆叠图用area函数画依次叠加光伏、风电、CHP和电网购入电量能直观看出能量来源结构。SOC曲线最好用阶梯图或带标记的折线图便于说明滚动更新的过程。日前和日内的对比图则可以画成分段折线两条线之间的包络区域用fill函数填充读者一眼就能看出纠偏效果。不要小看这些图很多审稿人先看图再决定要不要细看正文。画图的配色和线宽我习惯统一成一套模板保存成plotStyle()函数每次画图前调用省心很多。6. 最后说几句个人体会做这类仿真项目最有价值的部分其实不是把模型跑出结果而是通过代码把背后那套“预测、决策、反馈、再决策”的逻辑亲手搭出来。双层调度和滚动优化听起来是高大上的名词落到MATLAB里就是几个循环加一堆约束的排列组合但每次调试都是一次对物理逻辑的重新审视。我特别想留给大家一条建议运行代码之前先把时序图画出来。把负荷曲线、光伏曲线、电价曲线拉在同一张图上用眼睛找出那些“光伏中午过剩、电价高峰、负荷低谷”的关键时段再对照调度结果检查。模型可以出错但物理直觉永远是好帮手。如果你也正在折腾多能源微网的调度代码而且卡在滚动循环或者双层衔接上欢迎对照着这份“复盘”逐行检查。特别是SOC的初值传递和日前计划的偏差惩罚这两处理顺了整套代码基本就通了。