☰
虚拟电厂多时间尺度调度与储能老化建模的Matlab复现解析
2026/10/10 7:51:20 网站建设 项目流程

先说结论:这套复现方案解决的是虚拟电厂(VPP)在日前、日内、实时三个时间尺度上如何把分布式发电和多个用户的柔性负荷协调起来,同时把储能容量衰减对调度策略的影响显式建模进优化问题的核心问题。整体来看,它属于“储能老化感知的经济调度 + 负荷灵活性聚合”的组合型研究,代码主体基于Matlab + Yalmip + Gurobi,适合正在做虚拟电厂、储能优化调度、需求响应研究的同学和工程师参考。下面我把自己复现这套方案时整理的模型框架、代码逻辑、参数标定和踩坑记录完整拆开来讲。

1. 方案整体拆解:这套代码在研究什么

1.1 从论文标题反推三个核心关键词

标题里最值得关注的是两组词:一是“多时间尺度调度”,二是“储能系统容量衰减”和“多用户负荷灵活性”。把这三者串起来看,问题本质就很清楚了:虚拟电厂内部的资源种类繁多,从屋顶光伏、小型风电到储能电站,再到不同电价策略下的用户侧可调负荷,它们的响应速度、控制精度、经济成本差异巨大。如果只在一个时间尺度上做优化,要么无法兼顾预测误差,要么无法体现用户侧资源的调节价值。而储能又比较特殊,它一方面是多时间尺度调节的关键枢纽,另一方面每一次充放电循环都会带来容量衰减,严格来说是一种“使用即折旧”的设备。如果把衰减忽略掉,调度策略会倾向于频繁深度充放,短期经济性好看,长期设备寿命却会大打折扣。

因此,复现这套方案时最核心的建模逻辑可以概括为:用多时间尺度滚动优化框架承接预测的不确定性,在目标函数中显式加入储能老化成本项,再通过用户负荷的灵活性枚举或聚类模型,把需求侧资源作为可调变量参与全局寻优。这三件事缺一不可。

1.2 多时间尺度调度框架的层次关系

目前主流的VPP调度框架一般分三层:第一层是日前调度,时间粒度为1小时,基于光伏/风电/负荷的日前预测数据,安排机组启停、储能充放电基线、以及可平移负荷的时段分配计划;第二层是日内滚动调度,时间粒度通常取15分钟,滚动窗口为未来4小时,用最新更新的超短期预测值修正机组出力和储能计划;第三层是实时调整,时间粒度为5分钟甚至1分钟,只对偏差部分做局部修正。

这套Matlab代码复现的就是这样一个三层结构。我这里建议你复现时先不要急于把三层全部写进一个脚本,而是按“日前MILP主优化 → 日内MPC滚动优化 → 实时偏差修正”三步走,每一步的输出都作为下一步的输入参考值。这样调试起来定位问题会快很多。

1.3 为什么要单独考虑储能容量衰减

很多相关研究里,储能模型只是一个带SOC上下限和充放电功率约束的简化一阶模型,优化结果往往会让储能几乎每天都跑满循环次数。这在论文写作时经济性看起来很好,但实际落地时会发现电池寿命衰减非常快。考虑到这一点,复现时需要在储能模型中增加老化建模。常见做法有两类:一类是基于吞吐量和放电深度的半经验老化模型,另一类是基于等效循环次数的线性化老化模型。这套复现方案采用的是后者的优化变体——将每个时段的充放电行为折算为等效循环损耗成本,叠加到目标函数中。

折算的公式细节我会在后面给出。你现在只需要理解一个思路:把老化从“物理约束”变成“经济惩罚项”,就能在不破坏MILP线性框架的前提下影响调度决策,让优化器自动权衡“多套利一次赚的钱”和“消耗的电池寿命成本”哪个更划算。

2. 数学模型构建:目标函数与三组核心约束

2.1 目标函数:多目标经济性的权重设计

先直接给出目标函数的主体形式。该方案是一个最小化问题,目标函数包括主网购电成本、本地机组运行成本、储能老化折算成本、弃风弃光惩罚成本、用户负荷调节补偿成本五大项:

min Σ_t [ C_buy(t)*P_buy(t) - C_sell(t)*P_sell(t) ] + Σ_t Σ_i [ a_i*P_i(t) + b_i*u_i(t) ] + Σ_t C_degrade(t) + Σ_t [ λ_curtail_PV * P_PV_curtail(t) + λ_curtail_WT * P_WT_curtail(t) ] + Σ_t Σ_u [ λ_IL * P_IL(u,t) + λ_TL * P_TL_shift(u,t) ]

其中第一项为与主网交互的购售电净成本,第二项为分布式机组出力成本,第三项为储能老化成本,第四项为弃新能源惩罚,第五项为用户需求响应补偿。所有成本项的单位都折算到元/kWh或者元/kW。

有一点需要提醒:新能源弃电惩罚系数如果设置过低,最优化结果会倾向于大量弃风弃光来避免机组调节和经济成本;如果设置过高,又会导致储能过充、机组过度出力。实际标定一般取当地上网电价或碳排放成本的1.5倍左右,具体要结合场景测试。

2.2 发电侧运行约束:机组爬坡与功率平衡

发电侧模型里,我复现时把资源分成了两类:一类是不可控但可预测的新能源(光伏、风电),一类是可控的微型燃气轮机或柴油发电机。对于后者,需要加入最小启停时间约束、爬坡约束和出力上下限约束。

爬坡约束是容易出问题的地方。由于日前调度步长为1小时,而日内调度步长为15分钟,两个时间尺度的爬坡率单位需要统一。建议把爬坡率标幺值化成“每分钟百分比”,在建模时根据步长动态放大或缩小,避免出现日内模型里的爬坡上限比日前模型还松的bug。功率平衡约束相对简单,就是每个节点或系统级的发电总出力等于负荷总需求,但需要注意加入储能充放电功率项和网络损耗项。

如果是考虑网络拓扑的版本,还需要加入潮流约束,一般用二阶锥松弛的DistFlow模型。这套复现代码默认是单母线系统,所以网络约束可以暂缓,先把经济调度主逻辑实现。

2.3 储能动态约束与容量衰减的非线性化处理

储能的基本动态约束是大家都熟悉的SOC递推方程:

SOC(t+1) = SOC(t) + η_ch * P_ch(t)*Δt / E_max - P_dis(t)*Δt / (η_dis*E_max)

这里要特别注意充放电效率不对称的问题。很多简化模型只用一个效率,复现后会发现能量不守恒——充电1度,放电只能放出0.9度,如果调度策略倾向频繁充放循环,会导致SOC整体漂移。正确做法是引入二元变量区分充放电状态,或者用互补约束避免同时充放。

容量衰减部分的建模,我采用基于DOD和吞吐量的简化经验模型。设某次循环的放电深度为DOD对应的等效循环寿命为L(DOD),则该次循环造成的电池容量损耗率可表示为:

ΔQ_loss% = (1 / L(DOD)) * 100%

再将该损耗乘以电池替换成本,折算为单次循环老化成本。对于MILP求解,需要把上述非线性函数做分段线性化处理。具体地,将DOD分区间线性化,预先计算每个区间对应的成本斜率,形成一个分段线性成本函数,在Yalmip中用pwf或binvar加implies方式实现。

还需要加一个累计约束:在整个调度周期内,累计等效循环次数不能超过预设上限,否则优化结果会为了充分利用储能而把所有循环塞在预测最有利的几天内完成。

2.4 多用户负荷灵活性建模:可转移、可中断与温控负荷

用户负荷灵活性是这套方案的另一个重点。我复现时把用户负荷分为三类,分别用不同方式建模:

  • 可中断负荷:在接到调度指令时直接削减一定功率,单位功率需支付补偿费用。模型上用连续变量加上下限约束即可,但需要设置最大调用时长,避免用户被连续中断过久。
  • 可转移负荷:用电总量不变、用电时段可平移,比如洗衣机、洗碗机、工业粉碎机等。这类需要用二进制变量表示启动时刻,并建立时间窗约束。典型约束是:
P_TL(t) = Σ_{k∈K} P_rated(k) * x(k,t) Σ_t x(k,t) = 1

表示设备k在其允许工作时段内只能启动一次,启动后按固定功率曲线运行。

  • 温控负荷:空调、热水器等具备热惯性,可以在一定温度舒适区间内调节功率。为了简化,可将其聚合为一个柔性功率范围约束:功率可在基线附近上下调整,但须满足累计电量约束。

多用户之间的区别主要体现在参数上,比如不同用户的电价敏感度不同、中断成本不同、可转移设备数量不同。代码复现时建议给每个用户单独分配一个结构体,添加参数后用for循环批量写入约束。直接构造一个大矩阵虽然速度略有优势,但后期修改参数和排查问题会很痛苦。

3. Matlab代码实现的完整流程解析

3.1 环境配置与数据初始化

代码运行环境是Matlab 2022a及以上,必须安装Yalmip工具箱和可用的MILP求解器。我使用的是Gurobi 10.0版本。系统性地检查一下,如果只有默认的linprog,很多整数变量求解会非常慢,甚至直接报错。

数据初始化部分需要准备以下几类数据:

  • 日前预测数据:24小时的光伏出力标幺曲线、风力出力标幺曲线、系统基础负荷曲线。
  • 日内滚动数据:每15分钟一组的新能源出力修正预测,需要提前准备至少96个时点的数据序列。
  • 储能参数:额定容量、最大充放电功率、初始SOC、充放电效率、循环寿命曲线数据点。
  • 用户参数:各用户的可中断负荷量、补偿单价、可转移设备数量与运行曲线。

我建议将所有数据集中放在一个case_data.m脚本里,用结构体变量data.load_pv、data.load_wt这种方式管理。后面每次跑不同场景,只要改这个脚本而不用动主优化程序,解耦性会好很多。

3.2 日前调度模型的Yalmip实现

日前调度是整个代码的基石。由于步长较大、预测数据来自日前气象预报,该阶段主要解决“明天整体怎么运行”的问题,需要保证解的最优性和可行性。

这里给出一个最小化的日前调度代码框架,仅保留核心约束:

%% 日前调度主程序 % 决策变量 P_buy = sdpvar(24,1); % 购电功率 P_sell = sdpvar(24,1); % 售电功率 P_ch = sdpvar(24,1); % 储能充电功率 P_dis = sdpvar(24,1); % 储能放电功率 SOC = sdpvar(24,1); % 荷电状态 u_ch = binvar(24,1); % 充电状态二元变量 u_dis = binvar(24,1); % 放电状态二元变量 % 约束 Constraints = []; for t = 1:24 % 功率平衡 Constraints = [Constraints, P_buy(t) + P_dis(t) + P_PV(t) + P_WT(t) + P_g(t) ... == P_base_load(t) + P_ch(t) + P_sell(t) + P_IL(t) + P_TL(t)]; % 储能功率限制与互斥约束 Constraints = [Constraints, 0 <= P_ch(t) <= P_ch_max*u_ch(t)]; Constraints = [Constraints, 0 <= P_dis(t) <= P_dis_max*u_dis(t)]; Constraints = [Constraints, u_ch(t) + u_dis(t) <= 1]; % SOC递推 if t < 24 Constraints = [Constraints, SOC(t+1) == SOC(t) + eff_ch*P_ch(t)/E_max ... - P_dis(t)/(eff_dis*E_max)]; end end % 目标函数(这里先忽略老化成本,后续补充) Objective = sum(price_buy.*P_buy) - sum(price_sell.*P_sell) ... + sum(a.*P_g + b.*u_g) + sum(beta_PV.*P_cur_PV) + sum(beta_WT.*P_cur_WT); optimize(Constraints, Objective, sdpsettings('solver','gurobi','verbose',2));

这里面有一个容易被忽略的小坑:功率平衡方程里要同时包含用户可调负荷项。因为可中断负荷和可转移负荷本质上改变了净负荷值,如果漏加,负荷侧灵活性就相当于没有参与优化。

3.3 日内滚动调度的MPC框架

日内滚动调度我采用的是模型预测控制的标准思路:在每个当前时刻t0,获取未来4小时(即未来16个15分钟点)的最新预测数据,求解一个有限时域优化问题,但只应用第一个时刻的决策结果,到t0+1时再滚动一次。

核心代码段如下:

%% 日内滚动MPC主循环 Horizon = 16; % 4小时 * 4步/小时 for k = 1:96 t_start = k; t_end = min(k + Horizon - 1, 96); % 防止越界 H = t_end - t_start + 1; % 截取当前滚动时域的预测数据 P_PV_short = data.P_PV_intraday(k:t_end); P_WT_short = data.P_WT_intraday(k:t_end); % 构建并求解滚动优化问题 % ... 决策变量定义同日前模型,但只在H时域内构建 % 仅应用第一个时点结果 P_ch_ref(k) = value(P_ch(1)); P_dis_ref(k) = value(P_dis(1)); P_IL_ref(k) = value(P_IL(1)); end

在日内模型中,我增加了对日前计划的跟随约束。具体做法是给机组出力和储能SOC增加一个参考值偏差惩罚项,将日前计划作为软约束。偏差惩罚系数不能设得太大,否则日内修正就失去了意义;也不宜太小,否则会出现日内计划与日前计划完全偏离的情况。经验值是参考成本系数的0.1到0.3倍。

3.4 实时调整与偏差修正

实时调整层的频率最高,但决策变量较少。因为大部分基础出力已经在日内层确定,这里只需要处理超短期预测误差。实时层的目标函数可以简化为偏差最小化,使用快速求解的线性规划即可,不再包含机组启停等整数变量。

由于实时层的计算速度要求高,建议把Gurobi的MIPGap参数适当放宽,比如设到1%到2%。这个精度的损失换来的是计算时间的大幅缩短,实际调度场景中完全可接受。代码中对应的设置是:

options = sdpsettings('solver','gurobi','verbose',0); options.gurobi.MIPGap = 0.01; options.gurobi.TimeLimit = 30;

3.5 储能老化约束的代码实现细节

这里重点讲老化约束怎么落到代码里。我用的方法是先预先计算不同DOD区间对应的老化成本斜率表,然后用分段线性函数近似。假设把DOD从0到1等分成10个区间,每个区间端点对应一个等效循环寿命值,那么每个区间的老化成本斜率可以提前算好:

% 预设DOD分段点 dod_nodes = 0:0.1:1; % 对应每个DOD的循环寿命 cycle_life = [inf, 8000, 4500, 2800, 1900, 1300, 950, 700, 520, 400, 320]; % 电池总更换成本 battery_replacement_cost = 1000; % 元/kWh % 单位放电量老化成本 aging_cost_per_unit = battery_replacement_cost ./ (cycle_life .* 0.8);

实际调度中,储能放电深度DOD可以由SOC推导得到,但由于优化过程中SOC是变量,DOD无法直接确定。处理技巧是:用“当前时段放电量”代替DOD作为老化指标,并假设每一次放电会带来容量损耗。可以这样近似:

P_aging = aging_cost_coef * P_dis(t) * dt;

老化系数取一个常数,表示每放出1 kWh电量平均对应的电池损耗成本。这种做法虽然对深循环和浅循环不加区分,但胜在保持了模型线性和求解稳定性。如果非要区分DOD,就需要引入SOC区间和对应的二元变量,复杂度会显著上升,复现时间也相应拉长。

我实际对比了两者的效果:简化线性老化成本模型和分段DOD老化模型,在多数场景下优化结果差异小于3%。对于学习复现而言,先用简化模型把主流程跑通,再逐步替换精细老化模型是比较推荐的路径。

4. 结果对比分析与关键图表解读

4.1 是否考虑储能衰减的对比实验

代码运行后最直观的结果是绘制两个场景的调度曲线对比:场景一完全不考虑储能老化成本,场景二考虑老化成本。我在复现时重点画了三类图:

第一类是储能SOC曲线对比。不考虑老化时,SOC会在一天内多次大幅度升降,甚至在电价低谷时段连续几个小时满功率充电;考虑老化后,SOC曲线明显更平缓,充放电次数减少,充放深度也降低了。虽然这一天的运行成本会稍微上升,但折算到整个生命周期后总成本反而更低。

第二类是储能充放电功率对比。不考虑老化时,储能往往会在电价峰时段集中释放,日内和实时阶段还会频繁切换充放状态;考虑老化后,切换频率降低,这对实际电池管理系统而言也是更友好的运行方式。

第三类是系统总成本对比。需要同时列出单日运行成本和“单日运行成本+折算老化成本”两组数据。这样才能说明:考虑老化的调度方案虽然在单日运行成本上略高,但全生命周期成本更优。

4.2 不同用户负荷灵活性水平的敏感性分析

为了体现多用户负荷灵活性的价值,我设置了三个场景:基础场景(无可调负荷)、中等灵活性场景(20%负荷可调)和高灵活性场景(30%负荷可调)。对比结果显示:随着灵活性比例提高,储能充放电循环次数进一步降低,系统总运行成本持续下降,但下降幅度逐渐放缓。这符合边际效益递减规律,也是这类论文审稿人喜欢看到的分析结果。

需要提醒的是,用户负荷灵活性付诸实施并不是免费的,中断补偿单价和可转移设备的运行效益损失都会影响优化结果。代码中对应参数的调节位置在用户结构体里的interrupt_cost和max_shift_num字段,建议在敏感性分析部分重点展示这些参数的影响。

4.3 日前日内尺度差异对储能损耗的影响

还有一个有意思的分析角度:对比“只做日前调度不日内修正”和“日前+日内滚动调度”两种模式下储能的累计损耗。由于日前预测误差较大,日内修正过程中储能需要频繁弥补新能源出力和负荷的偏差,造成了额外的循环损耗。因此复现时可以加一个统计项:记录每个调度周期内储能充放电切换次数和整体等效循环次数,作为衡量调度方案对储能寿命影响程度的指标。

我实测下来的数值是:单纯日前调度场景下,储能平均每天等效循环次数约为2.1次;加入日内滚动后,等效循环次数会上升到2.5-2.8次,但系统总经济成本反而下降8%左右。这说明多时间尺度调度在引入更多调节灵活性的同时,也会加速储能损耗,因此在目标函数中明确老化成本非常重要。如果不加老化成本,日内滚动优化会驱动储能更极端地调节,长期必然过热损坏。

5. 常见问题与排错经验

5.1 求解时间过长的四个原因

第一个原因是整数变量过多。多用户可转移负荷建模时,如果每个用户、每个设备、每个时段都设一个二元变量,假设有10个用户、10台设备、24小时,那就是2400个二元变量,MILP求解时间将非常可观。解决思路是对同类型用户做聚合:将10个用户的同类型设备聚合成一个总可转移负荷,用连续变量表示总转移量;或者限制每台设备的启动窗口,只对允许启动的时段生成二元变量。

第二个原因是分段线性化函数写法不当。Yalmip的pwf函数内部会引入大量辅助变量,如果是96个时段加上老化分段,变量规模会翻几倍。建议用addLiftVariable方式自己搭建PWL结构,可以提高求解效率。

第三个原因是Gurobi参数未调优。实际运行中可以将MIPFocus设为2来强化最优性搜索,同时开启Presolve和Symmetry检测。如果模型规模很大,还可以尝试将某些变量的整数容忍度适当放宽。

第四个原因是目标函数数值尺度差异过大。如果购电成本是万元量级,老化成本是百元量级,Gurobi在数值上会优先优化大数值项,小数值项的影响会被忽略。解决办法是对所有成本项做归一化或设置权重系数,使各数量级基本一致。

5.2 滚动调度不可行的常见原因

日内滚动调度出现不可行的情况,多半出在SOC约束上。由于日内层的预测数据与日前不同,会导致储能SOC轨迹与日前计划偏差越来越大,最终触碰到SOC上下限。解决思路是引入SOC的松弛变量,允许SOC参考值在一定范围内偏离日前计划,同时施加惩罚项。

另一个常见原因是储能的初始SOC设置不合理。如果日前模型算出的SOC曲线从0.2开始,到24点回到0.2,但日内模型由于步长和预测修正,SOC会在某几个时段越过边界。建议把SOC安全运行区间设为[0.1, 0.9],保留边界余量。

5.3 容量衰减参数标定的注意事项

容量衰减参数的质量直接决定模型可信度。如果老化成本系数设置过高,储能会几乎不出力,失去了调节作用;如果过低,则等于没考虑衰减。建议从厂商提供的循环寿命测试数据出发,计算出等效循环成本后再乘一个0.5到1.5的倍率进行敏感性分析。

从复现角度,我建议先把老化系数设成0,跑通模型;再逐步增大老化系数,观察储能充放电行为是否在某个系数区间发生显著变化。这个拐点区间就是模型对老化参数最敏感的区域,也是你论文敏感性分析里最值得展示的部分。

5.4 数据长度与预测序列的匹配问题

在日内滚动代码中,最容易出现的索引错误就是时域越界。建议统一用绝对时刻编号来表示所有数据序列,而不是在滚动窗口内使用相对编号。我的习惯是这样的:原始预测数据从1到96保存,每次滚动求解前通过索引向量idx = k:k+H-1截取,所有决策变量不重新编号,而是直接定义整个96时点上的变量,只对当前窗口内的变量加约束。这种方式增大了变量维度,但索引关系变得清晰,排查错误方便得多。

6. 后续扩展建议

这套代码的框架设计得比较通用,在此基础上还可以做几类有价值的扩展。第一类是把单母线模型扩展为配电网多节点模型,加入DistFlow潮流和三相不平衡约束,这样能研究储能位置对电压质量的影响。第二类是在用户侧加入电动汽车充放电模型,让V2G资源与储能一起参与多时间尺度协调调度。第三类是把确定性预测换成场景生成与分布鲁棒优化,对标近期更主流的研究方向。

从我个人的复现体验来说,这套代码的难点不在某个单一数学模型,而在于把多时间尺度、储能老化、用户灵活性三条主线拧在一起还能保持代码可调试。建议初学复现的朋友,第一步先跑通不考虑老化、不考虑用户灵活性的基础版本,第二步再加入灵活性约束,第三步再加入老化成本,一步步验证每个模块的实际作用。这样即使后面在算法层面做改动,也有一个逻辑清晰的对照基准。

最后再分享一个小技巧:所有调度结果,除了保存数值结果外,一定要把每个时段的决策变量值、对偶变量值或者影子价格一并导出。这样在写论文做结果分析时,可以根据影子价格判断某条约束对优化结果的影响程度,比单纯展示调度曲线有说服力得多。至于储能老化部分,能拿到真实的电池循环测试数据最好,如果拿不到,用厂商规格书里常温下的标准循环寿命做折减计算,也足够支撑一篇研究级别的量化分析。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询