简介面向电力系统调度与优化研究者的MATLAB源码包围绕考虑特性分布的储能电站接入电网场景实现日前-日内-实时多时间尺度源储荷协调调度并融合需求响应机制可用于教学实验与课题验证。压缩包内含12个m脚本文件约21KB由多个主程序与储能约束、机组组合、母线导纳矩阵求解等配套函数构成结构清晰便于按模块理解和二次开发。从中可掌握多时间尺度调度模型的构建方法、储能特性分布处理逻辑以及需求响应与日前日内实时计划的衔接思路适合电力系统专业学生、科研人员及算法工程师参考。已有315人学习代码可直接运行调试也可作为论文复现、课程设计或工程项目验证的起点。1. 多时间尺度源储荷协调调度一天为什么需要三层决策如果只做一次日前调度把24小时的光伏、风电预测当成精准数据那么在午间光伏大发时储能可能会提前充满等到晚高峰再用。可现实是光伏预测每过一小时就会偏差5%到20%负荷也在变日前给出的储能动作到了下午往往已经不再最优。另一个反直觉的点是调度收益的80%来自日前但95%的违规风险来自接近实时的那段时间。所以多时间尺度源储荷协调调度的核心不是在一天里加更多约束而是把决策拆成日前-日内-实时三层每一层只处理自己带宽内能确定的误差。这里我直接用Matlab默认的优化工具箱linprog/intlinprog搭一套最小可运行框架并把需求响应嵌进三个时间断面适合正在做园区微电网或综合能源系统调度的人用来快速验证边界、算例和论文数据。2. 日前-日内-实时三层模型的约束设计与滚动衔接2.1 三层时间尺度的划分逻辑和各自承担的任务常见做法是把调度周期分成三层日前调度Day-ahead以1小时为步长覆盖24小时决策储能充放电计划、购售电计划和可中断负荷的预安排日内滚动Intra-day以15分钟为步长窗口4到6小时每15分钟重新优化一次修正预测误差和日前计划的偏差实时控制Real-time以5分钟为步长不做大范围优化只对储能功率做局部修正并读取AMI采集的分钟级负荷。为什么这样分因为优化问题的规模和解的可靠性需要权衡。如果直接用5分钟步长优化7天决策变量会爆炸到一个可笑的规模而且约束矩阵条件数变差linprog在稀疏矩阵下也未必很快。三层递进之后日前提供全局可行解日内提供可执行计划实时只负责消除残差。我一般用下面这套参数作为默认值。调度层步长滚动周期更新频率主要决策变量需求响应深度日前1h24h每日一次储能24点计划、购售电、可中断预安排日前签约量日内15min4h每15min储能修正序列、外购电修正日内实际削减实时5min15min每5min储能功率微调仅紧急响应日内和实时之间的配合是很多工程项目的分水岭。日内可以看作带滚动优化实时用比例-积分反馈或一个简单的二次规划做校正。我遇到过直接把日内优化结果下发到EMS结果储能动作过于频繁后来在实时层加了一个不动作死区命令变化小于5kW就不下发。这就是多时间尺度存在的另一个意义给控制设备留出喘息空间。2.2 日前调度全天约束集与目标函数怎么写日前调度是三层模型的基础所有时段的耦合约束都在这层建立。最常见的模型是混合整数线性规划MILP但为了快速验证我建议先写成线性规划LP。目标函数一般取最小化全周期运行成本加上一个惩罚项min F 购电成本 可中断负荷补偿 储能退化惩罚 - 售电收益其中购电成本用分时电价向量乘购电变量售电收益是上网电价乘售电变量。储能退化惩罚我习惯用一个很小的系数乘以储能出力绝对值避免储能无谓频繁充放。在LP里用绝对值只要把储能净功率拆成正负两个变量或者用一对不等式引入辅助变量。这里我先用净功率正负可变的做法因为演示代码足够清晰。功率平衡约束是每个时段的硬约束也是所有协调调度的核心购电 光伏 风电 储能净放电 需求响应削减量 负荷 售电。这个等式必须逐时刻成立。储能SOC约束把全天串起来这也是时序耦合的主要来源。SOC递推方程写成离散形式SOC(t) SOC(t-1) - Pb(t) * dt / Cap其中Pb(t)为正表示放电负表示充电。这是源储荷协调调度里最容易写错的一环尤其是符号约定。我统一用“储能对母线出力为正”的参考方向这样负荷平衡里的正负关系一目了然。2.3 日内与实时模型预测控制式滚动更新日内滚动优化的模型和日前几乎一样但有三处变化。第一预测曲线换成15分钟分辨率并只保留未来4小时窗口第二储能SOC的初始值不再假设而是取当前实际量测值第三需求响应的削减量上限从日前签约值缩小到日内允许值因为用户在执行日当天能响应的能力是有限的。每一步滚动完成后只执行窗口内第一个15分钟的计划然后推进到下一个时间点重新优化这就是模型预测控制的标准套路。实时层我一般不再调用linprog而是用一个极小规模的二次规划或者甚至直接用查表法。原因是实时层只有15分钟窗口决策变量少但要求响应必须在几百毫秒内完成。如果硬要在实时层跑完整MILPEMS的通讯周期会被拖垮。有一个不错的做法是把日内的第一点储能功率作为基准实时层只计算一个修正项 ΔPb目标函数是 ΔPb^2 最小约束是 ΔPb 在死区内并且修正后的SOC不超过分钟级硬限。这样既保证了跟踪又不会引入新的组合变量。这里要特别强调滚动衔接中的“硬约束”和“软约束”划分。日前层只允许超限的惩罚项在目标函数里出现不能直接出现在等式约束里日内层则可以把重要越限量作为松弛变量加上罚函数。实时层严禁增加新的等式约束否则容易无解。我们做需求响应时也遵循同一个原则日前定契约日内调边界实时只微调。3. 用Matlab linprog搭建源储荷协调调度的最小代码骨架3.1 数据结构规划把设备参数、预测曲线和状态量分开存储在用Matlab写多时间尺度调度之前我习惯先把数据结构定下来不然到日内滚动阶段变量索引会越写越乱。常见做法是定义一个结构体数组opt包含以下几个键T时间尺度点数、dt步长小时数、pv和wind曲线、load曲线、价格向量、储能参数、DR参数。另设一个state结构体保存当前SOC和上一层下发的计划值。为什么不用全局变量因为滚动优化时每次要改变SOC初值和预测窗口用结构体传入函数更安全。% 数据结构定义示例不含真实曲线只展示骨架 opt struct(); opt.T 24; % 日前24点 opt.dt 1; % 1小时步长 opt.Pv zeros(1, opt.T); % 光伏预测单位kW opt.W zeros(1, opt.T); % 风电预测 opt.Pload zeros(1, opt.T); % 负荷预测 opt.cBuy 0.8 * ones(1, opt.T); % 分时购电价 元/kWh opt.cSell 0.4 * ones(1, opt.T); % 上网电价 opt.PbuyMax 800; opt.PsellMax 500; % 电网交互限值 opt.Cap 1000; opt.PbMax 200; % 储能容量和功率上限 opt.SOC0 0.5; opt.SOCmin 0.2; opt.SOCmax 0.9; opt.PdrMax 80; opt.cDr 1.2; % 需求响应上限和补偿单价这段代码的意图是把所有可调参数集中到一处。后续不管是做日前还是日内只要复制这个结构体再改变T、dt和曲线即可而不必重写约束矩阵。参数说明中cBuy和cSell是行向量因为要与变量长度对应PbMax是绝对值上限实际上下限设置为负的同一个值用于充电限值。注意这里我特意把SOC初始值和储能的物理参数分开存避免在滚动循环里被覆盖。3.2 约束矩阵组装能直接跑通24小时日前调度的linprog调用下面给出一段可以放在脚本里直接跑的日前调度代码。为了演示方便预测曲线我用随机数占位真正使用时应换成读表数据或预测接口。决策变量顺序是x [Pbuy(1..T); Psell(1..T); Pb(1..T); Pdr(1..T)]其中Pb为正放电、负充电。% 变量索引偏移 n opt.T * 4; offSell opt.T; offPb 2 * opt.T; offPdr 3 * opt.T; % 目标函数系数 f [opt.cBuy, -opt.cSell, zeros(1, opt.T), opt.cDr * ones(1, opt.T)]; % 功率平衡等式 Pbuy - Psell Pb Pdr Pload - Pv - W Aeq zeros(opt.T, n); beq zeros(opt.T, 1); for i 1:opt.T Aeq(i, i) 1; Aeq(i, offSell i) -1; Aeq(i, offPb i) 1; Aeq(i, offPdr i) 1; beq(i) opt.Pload(i) - opt.Pv(i) - opt.W(i); end % SOC上下限约束转化为线性不等式 A*x b % SOC_i SOC0 - dt/Cap * cumsum(Pb(1:i)) A zeros(2*opt.T, n); b zeros(2*opt.T, 1); for i 1:opt.T idxPb offPb (1:i); % 约束1: SOC_i SOCmin (dt/Cap)*sum(Pb) SOC0 - SOCmin A(i, idxPb) opt.dt / opt.Cap; b(i) opt.SOC0 - opt.SOCmin; % 约束2: SOC_i SOCmax -(dt/Cap)*sum(Pb) SOCmax - SOC0 A(opt.Ti, idxPb) -opt.dt / opt.Cap; b(opt.Ti) opt.SOCmax - opt.SOC0; end % 变量边界 lb [zeros(1, opt.T), zeros(1, opt.T), -opt.PbMax*ones(1,opt.T), zeros(1,opt.T)]; ub [opt.PbuyMax*ones(1,opt.T), opt.PsellMax*ones(1,opt.T), ... opt.PbMax*ones(1,opt.T), opt.PdrMax*ones(1,opt.T)]; % 求解 options optimoptions(linprog, Algorithm, dual-simplex, Display, off); [x, fval] linprog(f, A, b, Aeq, beq, lb, ub, options); % 结果拆解 Pbuy x(1:opt.T); Psell x(offSell1:offSellopt.T); Pb x(offPb1:offPbopt.T); Pdr x(offPdr1:offPdropt.T);这段代码的核心是先把时间和设备顺序想清楚再组装矩阵。linprog求解时要保证Aeq的行数是时间点数T变量顺序与之前定义一致SOC不等式使用了累计矩阵在24点时相当于对全天累计SOC范围约束避免能量越限。参数里我在optimoptions中用了dual-simplex算法它在目标函数非光滑时比默认的interior-point更稳定且在大矩阵下内存占用更少。如果你用的是R2023b后的版本这个选项仍然保留。3.3 结果落盘与边界判断求解完不能只看fval我一般会立即做三件事校验功率平衡残差、检查SOC是否越限、和上一层计划绘制对比曲线。功率平衡残差可以用Aeq*x-beq的无穷范数来检查SOC可以从Pb递推出来并画图。如果SOC出现锯齿状频繁到边界说明储能退化惩罚系数太小或者日前预测曲线有明显突变这时不是调整求解器而是要回看数据和惩罚系数。% 校验功率平衡 resid full(Aeq * x - beq); if norm(resid, inf) 1e-6 error(日前调度功率平衡残差超限); end % 递推SOC并校验 SOC zeros(opt.T,1); SOC(1) opt.SOC0 - Pb(1)*opt.dt/opt.Cap; for i 2:opt.T SOC(i) SOC(i-1) - Pb(i)*opt.dt/opt.Cap; end if any(SOC opt.SOCmin - 1e-6) || any(SOC opt.SOCmax 1e-6) warning(SOC越限请检查惩罚系数); end这段代码的价值在于把“能跑”变成“有信心”。实际工程里我遇到过约束矩阵写错但fval仍然收敛的情况只有残差校验能抓住符号颠倒。最后一章中我会再给出更高效的稀疏矩阵做法和实时校正技巧。4. 需求响应建模可平移负荷、可中断负荷与实时电价参数怎么设4.1 可平移负荷的时间窗约束写法可平移负荷时间型需求响应是指洗衣机、热水器这类可以在一个时间窗内任意时点运行的负荷但总电量必须保持固定。在源储荷协调调度中这对应一组整数或连续变量对每个可平移设备i引入T维变量loadShift(i,t)表示该设备在时段t是否供电。约束为sum_t loadShift(i,t) E_i并且loadShift(i,t)在允许窗口 [t_start, t_end] 之外必须为0。在Matlab中用linprog时如果连续可调直接加等式约束如果设备是开关型需要整数变量就用intlinprog。下面给出连续型可平移负荷在矩阵里的构造片段设备总数M配合前面变量的顺序向后扩展。% 假设每个设备i的可运行窗口为 [winStart(i), winEnd(i)] M 5; % 可平移设备数 winStart ones(M,1) * 6; % 最早运行时刻 winEnd ones(M,1) * 22; % 最晚运行时刻 E 3 * ones(M,1); % 总电量 kWh shiftIdx n 1; % 新变量起始列 n n M * opt.T; AeqAdd zeros(M, n); beqAdd zeros(M,1); for i 1:M colStart shiftIdx (i-1)*opt.T; AeqAdd(i, colStart (winStart(i):winEnd(i))) 1; beqAdd(i) E(i); end % 合并到原Aeq/beq末尾同时lb中对应变量为0ub为1代码中关键参数是窗口和电量E。注意这里的loadShift变量单位可以是kW功率在单时段内恒定T个时段累加后得到kWh。将可平移负荷并入总功率平衡等式时需要在平衡方程里减去这些负荷变量否则会造成供大于求。实际操作中我会把这些变量累加写成行向量再加到Aeq的原功率平衡矩阵中而不是像上面示例这样直接扩在末尾否则两个等式之间缺少耦合。也就是说功率平衡矩阵要同时包含储能、购售电和可平移负荷所有变量在同一个等式里。4.2 可中断负荷和价格型需求响应怎么接入三层模型可中断负荷建模比较简单变量Pdr(t)表示在t时刻被削减的负荷功率上限PdrMax(t)由日前合同确定补偿单价c_dr(t)进入目标函数。在日前层PdrMax是全天固定值或分时段值日内层需要把PdrMax缩小到滚动窗口内剩余可削减量实时层一般不再新增削减除非日前合同里包含紧急削减条款。价格型需求响应则不同它不增加变量而是通过弹性矩阵修改负荷预测曲线。最简单的方法是用自弹性系数实际负荷P_DR P0 .* (1 epsilon .* (price_ref - price_actual) / price_ref)。其中P0是原始预测epsilon为自弹性通常取-0.1到-0.3。这样做的好处是不增加决策变量直接把修改后的负荷作为已知参数传入日前/日内模型。缺点是没法精确描述跨时段替代因此我一般把它作为场景分析工具而不放进在线优化模型。下面给出可中断负荷与价格型需求响应在日前目标函数中的参数给法参数日前典型值日内典型值说明PdrMax最大负荷的10%最大负荷的5%剩余可削减量随执行进度衰减c_dr1.5倍购电均价2倍购电均价越接近实时补偿越高epsilon-0.2-0.1价格弹性绝对值随临近时间变小SOC下限0.20.3日内必须为实时留更多空间从表格可以看到越靠近实时需求响应能力越弱、价格越高这是符合动态定价逻辑的。工程上经常犯的错是三层都用同一个PdrMax导致日内优化把需求响应全部用完实时层遇到冲击时只能弃光或切负荷。因此在日前层计算完后我会用代码把Pdr序列保存到state节点日内层读取时自动扣减已执行部分。4.3 用灵敏度分析确定需求响应参数参数设置不能拍脑袋。我一般会固定其他条件单独扫描epsilon在-0.05到-0.4之间的目标值变化画一张二维曲线观察什么时候系统运行成本开始下降平缓。这一步用Matlab脚本做非常简单不需要重新写模型只要把改好的P0传入调度函数循环得到成本和SOC曲线。如果发现epsilon从-0.25变到-0.3时成本只降了0.2%说明该场景需求响应的边际效益已经很低再加大弹性只会让用户舒适度恶化。可中断负荷的补偿价格同样需要扫描。一个值得注意的经验是当c_dr低于购电峰值电价的一半时调度结果基本不会削减负荷当c_dr高于峰值电价的1.5倍时削减量占满上限。也就是说设备真正起作用的区间很窄。此时可以用二分法找到临界点然后设定c_dr临界价的1.2倍这样既能让DR参与又不会让市场成本失控。5. 求解效率与工程落地的三个实战技巧5.1 用sparse矩阵替代全矩阵避免大数据量卡死多时间尺度调度里日内滚动16小时每15分钟的约束矩阵已经很大如果还和日前一样用零矩阵组装Matlab很容易吃内存。我的做法是一开始就构造索引向量然后用sparse函数组装。例如在第二节代码中Aeq可以在循环中收集行列值和值最后一行生成稀疏矩阵。linprog内部对稀疏矩阵的处理效率远高于全矩阵特别是在终端条件数较大时。% 稀疏矩阵组装示例替代全零矩阵循环 rows []; cols []; vals []; for i 1:opt.T % 功率平衡行i的各列系数 rows [rows; i; i; i; i]; cols [cols; i; offSelli; offPbi; offPdri]; vals [vals; 1; -1; 1; 1]; end Aeq sparse(rows, cols, vals, opt.T, n);这里要注意rows/cols/vals三个数组用分号拼接避免在循环内动态增长太慢可以先预分配。对于4小时窗口的日内优化这样拼接的时间可以忽略不计。实时层的问题规模小不需要用sparse。5.2 intlinprog热启动与整数容差设置如果你决定把可平移负荷建模为整数变量就需要用intlinprog而不是linprog。intlinprog支持提供初始可行解x0这在滚动优化中特别有用因为当前窗口的最优解往往和前一个窗口的最优解非常接近。我一般从优化结果中取出前一窗口对应的变量平移映射到新变量索引中作为x0并设置IntegerTolerance,1e-4来加快求解。options optimoptions(intlinprog, IntegerTolerance, 1e-4, ... LPMaxIterations, 200, Display, off); x0 zeros(n,1); % 将上一轮窗口解的对应片段放到x0中 x0(oldIdx) x_prev(newIdx); [x, fval] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, x0, options);热启动的代价是必须自己维护变量索引映射否则x0错位会导致求解器放弃初始点。这个技巧在日内滚动中效果明显求解时间可以缩短一半以上。如果你对索引映射不熟宁可不传x0让求解器自己找也不要传递错误初值。5.3 实时层的最小二乘校正不用再跑一次优化最后一层我推荐用最小二乘校正而不是重新求解调度。设日前或日内下发的储能参考功率为Pb_ref当前量测误差为deltaP由负荷突变引起则储能实际功率命令为Pb_cmd Pb_ref K * deltaPK是限幅系数通常取0.1到0.3。如果用Matlab的lsqlin做带约束最小二乘约束是Pb_cmd不超过储能功率上下限同时一分钟级的SOC不越限。这个二次规划极小算一次在1毫秒以内非常适合嵌入式控制器。% 实时层二次规划min (Pb_cmd - Pb_ref)^2 % 使用quadprog求Pb_cmd约束在功率上下限内 H 1; f_r -Pb_ref; lb_r -opt.PbMax; ub_r opt.PbMax; Pb_cmd quadprog(H, f_r, [], [], [], [], lb_r, ub_r, Pb_ref);注意上面这个调用没有包括SOC约束实际工程中需要把SOC越限量转成一个小的功率区间再叠加到lb_r和ub_r上。我经常看到有人把实时层也做成MILP结果一个周期超过2秒直接被EMS超时踢掉。正确的做法是让实时层尽量简单只解决“误差修正”不重新做经济调度。本文还有配套的精品资源点击获取