拿到“考虑特性分布的储能电站接入的电网多时间尺度源储荷协调调度策略(Matlab代码实现)”这个题目时,我的第一反应是:项目名虽然绕口,但信息量很足。拆开来看就是三件事——储能电站内部和外部的“特性分布”、日前/日内/实时三个“时间尺度”、源电源储能负荷的“协调调度”,最后全部落到Matlab上。储能电站不是单台电池,电站内部有成百上千个电池簇,SOC、SOH、内阻各不相同;源和荷也不是常数,风电光伏和负荷预测都有误差分布。把这些分布特性考虑进去,再放到多时间尺度上做优化,用Matlab跑通,就是整个项目的核心。这篇文章我会把模型设计、Matlab实现、算例结果和踩过的坑都摊开讲,适合刚接触电力系统优化调度的研究生,也适合想用Matlab搭调度模型的工程师。
1. 项目拆解:特性分布到底在“分布”什么
1.1 储能电站内部的特性分布
先说“特性分布”这个词。很多人一看到储能电站,习惯把它当成一个大电池处理,给一个总功率、总容量、一个SOC就完事。但在实际工程里,一个100MW/200MWh储能电站通常由几十个2.5MW/5MWh的储能单元组成,每个单元又由好多电池簇并联。这些单元的出厂批次、运行温度、循环次数不同,导致SOC可用区间、SOH衰减速度、最大充放电功率都不一样。
举个例子,我处理过的一组脱敏数据里,同一电站40台PCS对应储能单元,SOH最低的只有90%,最高的还有95%。如果把电站当成一个等效SOC为50%的大电池去做调度,下达一个100MW的放电指令,等到真正执行时,SOH偏低的那部分单元会先触到SOC下限,实际能放出的功率远小于计划值。这就是特性分布不建模会带来的直接后果。
所以项目里的第一层“分布”,是指储能电站内部单元之间的SOC、SOH、可调容量分布。调度策略不能只给电站一个总功率指令,还要考虑如何在单元之间分配功率,避免个别单元过充过放。这也是后面“实时分配”环节要解决的问题。
1.2 电源与负荷的随机特性分布
第二层“分布”来自电源和负荷侧。风电出力、光伏出力、负荷大小在调度时刻并不是已知常数,而是带预测误差的随机量。实际建模时,我会把风电功率预测误差近似成正态分布或贝塔分布,负荷误差近似成正态分布。误差的标准差随着预测时间尺度缩短而减小,这也是为什么要用多时间尺度框架。
源储荷协调调度的核心,就是把这三个带不确定性的对象放在同一个优化模型里决策。常规机组要留有功功率备用,储能要参与调峰和调频,负荷侧可以切一部分可中断负荷作为备用支撑。如果不考虑这些随机分布,只按预测值做单点确定性优化,运行方案往往会在实际执行时出现功率缺额或弃风弃光。
1.3 多时间尺度为什么是必选项
预测时间越长,误差越大。日前24小时的风电预测误差可能达到15%~20%,但未来1小时的预测误差能降到5%以下。如果只用日前计划去应对实时扰动,要么备用留太多,经济性差;要么备用不足,安全风险高。多时间尺度调度的思路是层层修正:
- 日前层:以小时为步长,制定机组启停和储能充放电的基准计划,决策周期24小时;
- 日内层:以15分钟到1小时为步长,滚动修正机组出力和储能计划,把最新的超短期预测装进模型;
- 实时层:以分钟甚至秒为周期,只做偏差调整和储能单元功率再分配,把前两层留下的误差吃掉。
这就像出门前看一周天气预报决定带不带厚衣服,出门前再看小时级天气预报决定要不要加伞,走到楼下发现下雨了再掏伞。每一层都只处理当前尺度能处理的不确定性,不把远期误差硬扛下来。
2. 调度策略整体设计:三层框架与数学模型
2.1 三层时间尺度分工
项目采用经典的“日前计划—日内滚动—实时调整”三层框架。各层的时间尺度、决策内容、求解频率各不相同,具体分工如下表所示。
| 层 | 时间粒度 | 决策变量 | 目标侧重点 | 求解频率 |
|---|---|---|---|---|
| 日前 | 1小时,共24点 | 机组启停、机组出力、储能充放电计划、弃风弃光量 | 全周期运行成本最小,启动次数合理 | 每天一次 |
| 日内 | 15分钟或1小时,滚动窗口4~8小时 | 机组出力修正量、储能充放电修正量、负荷削减量 | 跟踪日前计划并平抑预测偏差 | 每15分钟或每1小时一次 |
| 实时 | 1~5分钟 | 储能内部各单元功率分配、AGC机组调节量 | 消除短时偏差,考虑SOC/SOH差异 | 每个控制周期一次 |
这种分工的现实意义是:日前优化可以慢一点、精细一点,求解时间长一点也能接受;日内滚动必须快,所以模型要适当简化,比如把机组启停状态固定为日前得到的值;实时层则更轻量,甚至可以退化成规则分配。
2.2 日前计划的数学模型
日前计划是整个策略的“压舱石”。目标函数我通常写成运行成本最小,包含四部分:常规机组燃料成本、机组启停成本、储能寿命损耗折算成本、弃风弃光惩罚成本。如果负荷侧有可中断负荷,还要加负荷削减成本。
用公式表示就是:
[ \min ; \sum_{t=1}^{T}\sum_{g \in G} \left( C_{g}^{P} P_{g,t} + C_{g}^{U} u_{g,t} + C_{g}^{SD} v_{g,t} \right) + \sum_{t=1}^{T}\sum_{e \in E} C_{e}^{deg} \left( P_{e,t}^{ch} + P_{e,t}^{dis} \right) + \sum_{t=1}^{T} C^{W} \left( P_{w,t}^{fore} - P_{w,t}^{use} \right) ]
其中 (u_{g,t}) 是机组启动状态二进制变量,(v_{g,t}) 是停机状态二进制变量,(P_{g,t}) 是机组出力,(P_{e,t}^{ch})、(P_{e,t}^{dis}) 是储能充放电功率,(P_{w,t}^{fore})、(P_{w,t}^{use}) 是风电预测和实际消纳功率。
约束条件里最基础的是功率平衡约束:
[ \sum_{g \in G} P_{g,t} + P_{w,t}^{use} + \sum_{e \in E} P_{e,t}^{dis} = P_{t}^{load} + \sum_{e \in E} P_{e,t}^{ch} ]
然后是常规机组上下限约束、爬坡约束、最小开停机时间约束,储能SOC动态约束:
[ SOC_{e,t+1} = SOC_{e,t} + \eta_{e}^{ch} P_{e,t}^{ch} \Delta t - \frac{P_{e,t}^{dis} \Delta t}{\eta_{e}^{dis}} ]
以及SOC上下限约束、储能充放电功率上下限约束、充放电互斥约束。对于网络约束,如果项目只做经济调度,可以暂时忽略网架;如果要考虑潮流,我一般用直流潮流模型,把线路有功功率限制写成线性约束。
2.3 日内滚动修正模型
日内滚动是基于模型预测控制思想。上一层的日前计划不是“圣旨”,而是“参考线”。日内层每隔一段时间,利用最新超短期预测数据重新优化未来几个小时的计划,但只执行当前时段的决策。
日内层决策变量用修正量表示,比如机组出力修正量 (\Delta P_{g,t})、储能功率修正量 (\Delta P_{e,t})。目标函数包含两部分:一是对日前计划的偏离惩罚,防止日内把日前计划改得面目全非;二是对预测偏差的平衡,最小化弃风和切负荷量。
目标函数可以写为:
[ \min ; \sum_{t \in H} \left( \lambda_{g} |\Delta P_{g,t}| + \lambda_{e} |\Delta P_{e,t}| + C^{W}\Delta P_{w,t}^{curt} + C^{LS}\Delta P_{t}^{LS} \right) ]
这里面向量 (H) 是滚动优化窗口。日内层通常把机组启停变量固定为日前解,把模型变成LP或QP,求解速度很快。这一步是整个策略能否落地的关键,因为日内求解频率高,模型太重就跑不动。
2.4 实时调整层的设计
实时层不做全局优化,重点是“怎么分”。储能电站收到实时有功指令后,需要按照各单元SOC可用容量、SOH健康状态、功率响应能力等实际情况分配功率。最常用的分配原则是等比例可调容量分配:
[ P_{e,i}^{set} = P_{e}^{order} \cdot \frac{SOC_{i} - SOC_{i}^{min}}{\sum_{j} \left( SOC_{j} - SOC_{j}^{min} \right)} ]
这样SOC偏高的单元在放电时多承担一些功率,充电时少承担一些,避免个别单元率先触界。如果实时层还需要调节常规机组,我习惯用最简单的PID或AGC分配逻辑,再叠加储能的快速功率补偿。
3. Matlab实现:从数据到可跑程序
3.1 工程代码模块划分
整个Matlab工程我按功能拆成六个模块:主程序、数据生成、日前优化、日内滚动、实时分配、结果绘图。每个模块独立成文件,方便调试和替换数据。
| 文件 | 功能 | 关键函数/脚本 |
|---|---|---|
| main.m | 总入口,按顺序调用各模块 | 无 |
| setup_case.m | 读入系统参数、负荷和新能源数据 | loadData, initialize |
| build_dayahead_model.m | 构建日前优化模型并求解 | buildConstraints, solveDayAhead |
| run_rolling_mpc.m | 日内滚动优化主循环 | updateForecast, solveMPC |
| dispatch_allocation.m | 实时储能功率分配 | allocatePower |
| plot_results.m | 画SOC、机组出力、弃风率结果 | 无 |
这种结构的优点是:数据模块独立,换一个算例系统时只需要改setup_case.m;优化模型和求解逻辑分开,想换求解器只需改一处;绘图单独放,方便最后调图。
3.2 数据准备和参数设置
我习惯把所有原始参数集中到一个xlsx文件里,由setup_case.m统一读取。参数至少包括:常规机组容量、爬坡率、启停成本、燃料成本系数;储能单元个数、额定功率、额定容量、初始SOC、SOC上下限、充放电效率;负荷曲线、风电预测曲线、光伏预测曲线;以及各层时间尺度和滚动窗口长度。
注意单位统一。我踩过最大的坑就是功率用MW、能量用MWh、时间用小时时很容易混,特别是储能SOC更新公式里 (\Delta t) 的单位。我在代码里统一用MWh作为能量单位,所有功率乘以时间后得到能量,这样SOC约束不会差数量级。
3.3 日前优化模型核心代码
如果安装了YALMIP和Gurobi/CPLEX,日前模型写起来会非常清爽。下面是我常用的日前优化核心代码框架。
% 变量定义 P_g = sdpvar(n_gen, T, 'full'); % 常规机组出力 u_g = binvar(n_gen, T); % 机组启动状态 P_ch = sdpvar(n_ess, T, 'full'); % 储能充电功率 P_dis = sdpvar(n_ess, T, 'full'); % 储能放电功率 SOC = sdpvar(n_ess, T+1, 'full'); % SOC状态 P_w_use = sdpvar(n_wind, T, 'full'); % 风电消纳功率 % 目标函数 Objective = sum(sum(C_g * P_g)) + sum(sum(StartCost .* u_g)) + ... sum(sum(deg_cost * (P_ch + P_dis))) + ... sum(sum(curtail_penalty * (P_w_fore - P_w_use))); % 约束 Constraints = []; % 功率平衡约束 for t = 1:T Constraints = [Constraints, sum(P_g(:,t)) + sum(P_w_use(:,t)) + sum(P_dis(:,t)) ... == load(t) + sum(P_ch(:,t))]; end % 机组出力上下限和爬坡约束 for g = 1:n_gen for t = 1:T Constraints = [Constraints, P_g_min(g)*u_g(g,t) <= P_g(g,t) <= P_g_max(g)*u_g(g,t)]; if t > 1 Constraints = [Constraints, P_g(g,t) - P_g(g,t-1) <= ramp_up(g)]; Constraints = [Constraints, P_g(g,t-1) - P_g(g,t) <= ramp_down(g)]; end end end % 储能SOC及功率约束 for e = 1:n_ess Constraints = [Constraints, SOC(e,1) == SOC_init(e)]; for t = 1:T Constraints = [Constraints, SOC(e,t+1) == SOC(e,t) + eta_ch(e)*P_ch(e,t) - P_dis(e,t)/eta_dis(e)]; Constraints = [Constraints, SOC_min(e) <= SOC(e,t+1) <= SOC_max(e)]; Constraints = [Constraints, 0 <= P_ch(e,t) <= P_ch_max(e)]; Constraints = [Constraints, 0 <= P_dis(e,t) <= P_dis_max(e)]; Constraints = [Constraints, P_ch(e,t) + P_dis(e,t) <= P_ess_max(e)]; % 互斥约束松弛 end end % 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(Constraints, Objective, ops);这里充放电互斥我用的是 ( P_{ch} + P_{dis} \leq P_{max} ) 来代替二进制互斥变量。因为一个储能单元在同一时刻不可能一边充一边放,但这个松弛约束在目标函数有成本惩罚时,大部分情况下解出的结果不会同时充放电。如果非得严格互斥,再引入二进制变量:
z = binvar(n_ess, T); Constraints = [Constraints, P_ch <= P_ch_max * z]; Constraints = [Constraints, P_dis <= P_dis_max * (1 - z)];注意引入二进制变量会把LP问题变成MILP,求解时间会上去。如果你的储能单元数量多、调度周期长,先用松弛版本跑通,再对照二进制版本的差别。
3.4 日内滚动优化主循环
日内滚动层用for循环按时间窗口推进。核心思路是每步读取最新的预测数据,构建一个子优化问题,求解后只取第一个时段的修正量执行,然后窗口平移。
for k = 1:num_steps_intraday % 获取最新预测 wind_forecast_now = get_latest_forecast(k); load_forecast_now = get_latest_load_forecast(k); % 构建日内优化模型 [Constraints, Objective] = build_intraday_model(P_g_dayahead, wind_forecast_now, ...); % 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 0); optimize(Constraints, Objective, ops); % 取当前时段结果执行 delta_P_g(:,k) = value(delta_P_g_var(:,1)); delta_P_ess(:,k) = value(delta_P_ess_var(:,1)); end日内模型里最重要的一点是加“偏向日前计划”的惩罚项,比如 ( \lambda_g |\Delta P_g| )。如果不加这个项,优化器可能会重新规划出一个完全不同的出力曲线,导致机组频繁调整,实际工程里很难接受。
3.5 实时分配逻辑实现
实时分配逻辑我推荐写成一个独立函数,输入是当前SOC、SOH、储能单元上下限和总功率指令,输出是每个单元的功率指令。最经典的按可用容量等比例分配,代码只有几行:
function P_set = allocate_power(SOC_now, SOH_now, SOC_min, SOC_max, P_order) avail_dis = (SOC_now - SOC_min) .* SOH_now; avail_ch = (SOC_max - SOC_now) .* SOH_now; if P_order >= 0 % 放电需求 weights = avail_dis / sum(avail_dis); P_set = P_order * weights; else % 充电需求 weights = avail_ch / sum(avail_ch); P_set = P_order * weights; end end这种分配方式没有用优化算法,胜在快、容易理解。如果要做得更精细,可以加一个带SOC均衡目标的实时优化分配,用quadprog求解,但实际工程里,规则分配在大多数时候已经够用。
4. 针对“特性分布”的关键处理技术
4.1 储能SOC一致性建模
前面说了储能电站内部单元SOC分布不均衡,如果模型里把这些单元都当成独立变量,目标函数里最好能体现对“SOC一致性”的偏好。最直接的做法是在目标函数里加一个SOC方差惩罚项:
[ J_{soc} = \lambda_{soc} \sum_{t=1}^{T} \sum_{e=1}^{E} \left( SOC_{e,t} - \overline{SOC}_{t} \right)^2 ]
其中 (\overline{SOC}_{t}) 是所有储能单元在t时刻的平均SOC。这个二次项会让优化器倾向于让各单元SOC尽量接近,避免出现某些单元满载、某些单元空载的情况。
不过二次项会让MILP变成MIQP,求解变慢。工程上可以用分段线性近似替代方差项:把SOC偏差绝对值罚掉,用线性项逼近方差效果。如果用的是二次规划求解器,直接保留平方项也能接受。我跑过60个储能单元、24个时段的算例,引入二次项后求解时间从1.8秒涨到3.2秒,还在可接受范围内。
4.2 随机场景生成与缩减
考虑电源和负荷不确定性的常用做法是场景法。先用Monte Carlo抽样生成大量预测误差场景,再用场景缩减技术挑出少数代表性场景。Matlab自带kmeans,做场景聚类非常方便。
具体步骤:
- 生成N个误差场景,每个场景包含24个时段的误差序列;
- 以场景向量作为特征,用
kmeans聚成K类; - 每类取聚类中心作为典型场景,权重为该类场景数占总数的比例;
- 把典型场景代入两阶段随机优化模型,目标函数变为各场景目标值的加权和。
K值很关键。太小了覆盖不了尾部风险,太大了计算量爆炸。根据我的经验,K取5~10个场景就能覆盖95%以上的波动特性,再往上增加场景数量,精度提升很小,计算时间却线性增长。我一般先跑K=5,再试K=10对比稳定性。
4.3 储能寿命损耗折算
储能参与调度会影响寿命,如果目标函数里不考虑寿命,优化器会让储能频繁深度充放,虽然运行成本低了,但实际换电池成本很高。常用的办法是把寿命损耗折算成经济成本。
工程上常用DoD(放电深度)和循环次数关系曲线:( N_{cyc} = a \cdot DoD^{-b} )。比如某磷酸铁锂电池,100%DoD时循环3000次,50%DoD时循环6000次,拟合得到a和b。每次等效循环的损耗成本是:
[ C_{deg} = \frac{C_{replace}}{N_{cyc}(DoD) \cdot E_{rated}} ]
调度模型里用SOC变化绝对值近似等效循环深度:
[ DoD_{e,t} \approx |SOC_{e,t} - SOC_{e,t-1}| ]
然后在目标函数加 ( C_{deg} \cdot DoD ) 惩罚。这个线性近似不算特别精确,但胜在可导、能嵌进优化模型,而且能抑制储能的“过度活跃”问题。
5. 算例分析与调优心得
5.1 算例设置与结果对比
我搭了一个简化算例做验证:常规机组4台,总装机600MW;风电装机200MW;储能电站由20个5MW/10MWh单元组成,总容量100MW/200MWh;负荷峰值为500MW。日前时间尺度为24点,日内滚动窗口4小时,步长15分钟。数据采用某地区典型日负荷和风电曲线。
跑完三组策略对比:
| 指标 | 不含特性分布 | 考虑SOC分布+随机场景 | 考虑SOC分布+随机场景+寿命折损 |
|---|---|---|---|
| 运行成本(万元/日) | 158.2 | 160.4 | 162.1 |
| 弃风率(%) | 4.8 | 3.2 | 3.0 |
| 储能SOC方差 | 0.036 | 0.012 | 0.009 |
| 等效循环次数 | 1.12 | 1.03 | 0.91 |
从表格能明显看出,考虑特性分布后运行成本小幅上升,但弃风率明显下降,储能SOC更健康,循环次数也更低。这说明特性分布模型本质上是拿一点经济成本换运行安全性和设备寿命。
5.2 关键参数敏感性分析
调节惩罚系数是调优的主要手段。弃风惩罚系数 (C^{W}) 设太高,优化器会为了消纳风电让常规机组压到最低出力,甚至牺牲备用;设太低,又会出现大量弃风。我常用做法是先按上网电价加补贴定一个基准值,再在50%~200%之间做敏感性扫描,找到成本和弃风率的拐点。
储能寿命折损系数 (C_{deg}) 同理。系数太小时,储能会被当作“免费调节资源”频繁充放;系数太大时,储能干脆趴着不动,发挥不了调节作用。我在实际调参时,会把寿命折损系数设为储能更换成本除以全生命周期等效循环次数的量级,这样调度结果比较接近真实运行。
5.3 计算效率问题与提速
日前模型如果含机组启停二进制变量,用YALMIP+Gurobi求解4台机组、24时段的MILP,通常1~3秒能出结果。当储能单元数量增加到60个,且每个单元都单独建模时,变量规模翻倍,求解时间可能涨到30秒。这时候需要做聚合简化:
- 将同型号、同SOH水平的储能单元合并成一组,用“聚合SOC”建模;
- 只在实时分配层保留单元级模型;
- 使用滚动窗口缩短日内优化时域,不优化全天。
我一般保留最多4~5组储能聚合体,既能反映SOC分布差异,又不至于让MILP变得太重。实时分配层的单元级计算量很小,完全不影响整体性能。
6. 常见问题与实战避坑
6.1 求解器报“infeasible”怎么办
如果YALMIP返回Infeasible problem,大概率是约束之间产生了矛盾。我最常见的诱因是:储能初始SOC设置不合理,导致接下来几个时段无论怎么充放都违反SOC上下限;还有功率平衡约束中负荷和新能源数据不匹配,出现了长时间缺电。
排查方法是从简单到复杂:先关掉机组启停变量,固定机组出力,只跑储能调度,看是否可行;再逐步加入爬坡约束、SOC约束。哪一步开始不可行,问题就在哪一步。另外检查一下SOC更新公式里的充电效率,有些资料里效率是乘在放电侧,写反了也会导致SOC越界。
6.2 数值尺度病态问题
优化模型里如果变量范围差异过大,比如机组出力从0到600MW,储能SOC从0.1到0.9,而惩罚成本有小数,求解器很容易出现数值问题。解决办法是统一单位基准:要么全系统用标幺值,要么把功率、能量、成本都转成同一个数量级。
我自己的习惯是:功率用MW,能量用MWh,成本用万元。SOC限制写成“0.1~0.9”这种小数,而机组出力写成“100~600”这种整数,条件数确实比较差。后来我把SOC上下限改成以10MWh为基准的可用电量值,比如“20~180MWh”,数值尺度就顺了。
6.3 充放电互斥约束导致求解变慢
严格互斥约束需要二进制变量。储能单元多的时候,二进制变量会成倍增加,MILP求解时间急剧上升。这里可以分两种情况处理:
- 如果储能参与的是谷充峰放这种长时间尺度调度,直接用 (P_{ch}+P_{dis} \leq P_{max}) 松弛就足够,解出来几乎不会同时充放;
- 如果储能还参与调频,实时功率波动大,那就得用真正的互斥约束,但可以在实时层用MPC或规则分配,不一定在日前层加二进制变量。
6.4 日内滚动模型“飘”得太远
日内滚动只优化未来4小时,求解器会为了当前窗口利益把机组出力大幅调整,导致与日前计划脱节。解决方法是给修正量加上二次惩罚,并且每隔一段时间强制刷新日前计划。更工程化的做法是,日内层只允许储能和可中断负荷动作,常规机组保持日前基准,除非出现严重预测偏差才允许调整机组出力。
实际操作中,我发现惩罚系数 (\lambda_g) 取燃料成本系数的1.5倍左右比较合适,太小则机组修正频繁,太大则日内层形同虚设。
6.5 Matlab版本与求解器接口
YALMIP支持Gurobi、CPLEX、Mosek、OSQP等多种求解器。如果电脑没装商业求解器,可以用Matlab自带的intlinprog配合YALMIP求解MILP,也可以直接手写linprog做LP模型。我建议先装OSQP或Gurobi,求解速度快一个数量级。但要注意,YALMIP和求解器的版本兼容性偶尔会出问题,常见报错是“solver not found”或“license error”,这种时候先把LS接口和求解器路径检查一遍。
7. 一些实操体会
这个项目最考验人的不是写代码,而是想清楚“特性分布”和“多时间尺度”怎么在模型里互相咬合。我最初的版本只考虑了储能SOC单点约束,结果实时仿真时频繁出现某几个单元过放的告警;后来把SOC分布方差和单元级分配加进去,结果才真正能落地。调试过程中我发现,SOC一致性指标比总运行成本更能反映方案的鲁棒性,所以现在做类似项目,我会把SOC方差、SOH均衡度放在和成本同等的地位去汇报,而不是只报一个总费用。
如果你也想复现这个策略,建议从单储能单元+日前调度开始,跑通以后再扩展到多单元和日内滚动。一开始别急着加随机场景和寿命折损,等基础模型稳定了,再一层层往上面叠复杂度。最后一个小技巧:所有参数统一放Excel表里,代码里不写死任何数字,这样换一组算例、调一组惩罚系数会轻松很多。