简介：面向电力系统调度与优化研究者的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在稀疏矩阵下也未必很快。三层递进之后，日前提供全局可行解，日内提供可执行计划，实时只负责消除残差。我一般用下面这套参数作为默认值。
| 调度层 | 步长 | 滚动周期 | 更新频率 | 主要决策变量 | 需求响应深度 |
|---|---|---|---|---|---|
| 日前 | 1h | 24h | 每日一次 | 储能24点计划、购售电、可中断预安排 | 日前签约量 |
| 日内 | 15min | 4h | 每15min | 储能修正序列、外购电修正 | 日内实际削减 |
| 实时 | 5min | 15min | 每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分钟窗口，决策变量少，但要求响应必须在几百毫秒内完成。如果硬要在实时层跑完整MILP，EMS的通讯周期会被拖垮。有一个不错的做法是把日内的第一点储能功率作为基准，实时层只计算一个修正项 Δ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.T+i, idxPb) = -opt.dt / opt.Cap; b(opt.T+i) = 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(offSell+1:offSell+opt.T); Pb = x(offPb+1:offPb+opt.T); Pdr = x(offPdr+1:offPdr+opt.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中对应变量为0，ub为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_dr | 1.5倍购电均价 | 2倍购电均价 | 越接近实时，补偿越高 |
| epsilon | -0.2 | -0.1 | 价格弹性绝对值随临近时间变小 |
| SOC下限 | 0.2 | 0.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; offSell+i; offPb+i; offPdr+i]; 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 * deltaP，K是限幅系数，通常取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超时踢掉。正确的做法是让实时层尽量简单，只解决“误差修正”，不重新做经济调度。
本文还有配套的精品资源，点击获取
