高比例可再生能源下虚拟电厂多时间尺度调度与储能衰减建模的Matlab复现指南

前阵子我们课题组打算复现一篇电源领域顶刊上的储能调度论文题目正好就是“高比例可再生能源并网下虚拟电厂多时间尺度调度及衰减建模”。那篇摘要写得很漂亮又是灵活性、又是储能成本、又是衰减但真要自己动手在Matlab里跑通才知道坑有多深。我前后折腾了差不多两周把多时间尺度调度框架、储能容量衰减的成本化处理、求解器对接、结果可视化的逻辑全部捋顺了。这篇博文就把整个复现过程掰开揉碎讲一遍包含数学建模、代码关键片段、参数调优和调试心得给想复现类似工作的同学一个可以直接上手的参考。先说清楚这东西到底是干嘛的。所谓高比例可再生能源并网说白了就是风、光在系统里占了很大一块它们出力的随机性和波动性会让电网措手不及。过去靠火电顶上去就行现在调峰压力太大必须把储能、柔性负荷、可中断负荷这些“灵活性资源”打通。虚拟电厂干的事就是把这些分散资源聚合成一个整体参与调度。而储能恰恰是其中最灵活、但也是最花银子的环节。频繁充放电电池寿命掉得快换电池又是一大笔钱。所以这篇文章的科学问题就是在保证系统不缺电、不弃风弃光的硬约束下找到一个储能充放电策略让“运行成本储能老化成本弃电惩罚”的总账单最小。这种问题用数学语言讲就是一个带时序约束的优化问题常见的解法是混合整数规划。我复现下来的整体感受是表面上难点是调度算法真正考验人的其实是储能衰减建模和怎么把这一坨非线性关系塞进优化模型。很多论文对衰减只是一笔带过实际跑的时候你会发现在Matlab里要么表达式写错了要么求解器不给面子报错要么模型规模太大直接算不动。下面按我的复现顺序把每一步的设计逻辑、代码写法、调试过程都写清楚。1. 项目整体设计与核心思路拆解1.1 为什么高频次储能充放会加速退役先建立一个直觉。手机电池用一两年就不行了是因为每天充放电化学活性物质一直在损耗。对应到电网储能一天可能充放电十几个循环再加上大功率快充快放储能寿命衰减会比想象中快得多。传统调度论文通常把储能当成一个“能量中转池”只约束荷电状态在0.1到0.9之间却完全不考虑每次充放电循环对寿命的影响。结果算出来的方案确实让系统运行成本很低但实际执行一两年后电池就报废了这笔账根本划不来。这也是那篇论文的立论基础。它把容量衰减和调度决策耦合起来在目标函数里显式加入“每充放一笔电电池就损失几厘钱寿命成本”的项。一旦有这个成本优化器就会自动选择合适的放电深度和充放电频率不会为了省电费去拼命折磨电池。复现时要想清楚我们不是要去预测电池的化学状态而是要把工程上的寿命经验公式变成优化器看得懂的数学表达式。1.2 为什么要分三个时间尺度调度风光预测误差有个特点远期预测偏差大近期预测偏差小。如果只看未来24小时一次性把计划做死遇到突发风速变化就只能眼睁睁看系统失衡。所以论文采用日前、日内、实时三层结构。日前计划算一遍确定机组启停、储能初始SOC曲线、购电计划日内每隔固定时段滚动更新修正风光功率的短期预测实时层只调节小范围内的偏差比如储能快速响应和切负荷。这个结构和实际电网调度层级是一致的。我复现的时候三层之间要传递状态变量。尤其是SOC日前结束时的SOC就是日内开场时的初始SOC不能各算各的。如果这里接错了后文所有结果都会对不上。1.3 衰减建模为什么是这次复现的核心衰减建模直接决定模型复杂度。最简单的做法是用“吞吐量法”认为总寿命折损只与累计充放电电量有关每千瓦时一定钱。工程上更讲究的是考虑放电深度浅充浅放对电池友好深充深放则加速老化。论文里用的是“等效循环寿命法”把不同放电深度换算成等效全循环次数从而估算每一次动作造成的寿命损耗。但问题来了放电深度是一个带绝对值的非线性函数直接丢进优化问题会非常难解。我的做法是把它线性化把充放电状态拆成二进制变量SOC与端到端的差额用绝对值约束成线性再用分段线性函数逼近寿命曲线。这样保持了论文核心意图又保证模型能跑。后面会详细展开。2. 关键模型与数学公式解析2.1 目标函数里的三笔账怎么平衡总目标函数大概是最小化系统总运行成本 火电机组燃料费用 从上级电网购电成本 储能充电和放电导致的寿命损耗成本 弃风弃光惩罚。这几项的量纲完全不同机组成本是每兆瓦时几十块弃电惩罚是每兆瓦时几百块储能衰减成本更是可能只有每兆瓦时几块钱。如果直接相加衰减那项会被完全淹没等于白建模。所以复现时必须做归一化我采用的方式是先把各项除以各自的基准值变成无量纲相对值再乘以权重系数。权重系数怎么设参考论文里的表。没有表就自己调先设储能权重为零跑通调度再逐步加上去看SOC轨迹是否明显变平缓。一个容易忽略的点购电成本和售电价格在日内不同时段有峰谷差。目标函数里必须区分买电和卖电两个单价如果只用单一电价虚拟电厂的经济效益会被严重高估。那篇论文标配是日内和实时电价曲线我复现时用了模拟的峰谷电价并加噪后续替换真实市场数据也方便。2.2 储能衰减模型的数学表达与线性化我采用的等效循环寿命法表达式为当单次循环放电深度为DoD时对应等效全循环次数为DoD乘以一个系数而实际寿命折损等于有效电量除以电池额定全循环次数下可处理的总电量。简化成优化目标时每一项衰减成本近似等于储能充电或放电功率乘以一个随SOC变化的系数。这个系数随放电深度增大而变大体现了“深充深放更伤电池”。处理方式是把SOC分成若干固定区间每个区间给一个近似的寿命损耗系数。充放电过程中SOC如果跨越多个区间就按区间累加损耗。我在Matlab里用的大M法来实现这个分段逻辑。下面这段是核心写法% SOC分段线性化示例 % s_break [0.1 0.3 0.5 0.7 0.9]; % lambda [0.8 1.0 1.2 1.5]; % 每段斜率影响系数 % 用二进制变量选择激活段 seg_active binvar(T, n_seg, full); soc_seg sdpvar(T, n_seg, full); % 保证每个时段只激活一个段 for t 1:T constraints constraints [sum(seg_active(t,:)) 1]; constraints constraints [soc_seg(t,:) min_seg * seg_active(t,:)]; constraints constraints [soc_seg(t,:) max_seg * seg_active(t,:)]; end % 损耗成本为各段 SOC 对应系数的加权和 degradation_cost sum(sum(lambda .* soc_seg)) * cost_per_unit;这样衰减成本和SOC的关系就是线性的优化器可以正常处理。我最开始偷懒直接用SOC的二次函数来近似衰减成本结果求解器跑一个小时都不收敛改成分段线性后几分钟就出结果了。2.3 约束条件里最容易出错的三个地方第一功率平衡约束永远是等式加上网损修正。但也别忘了虚拟电厂内部还有光伏和风电的出力上限输出功率不能超过预测值否则就是凭空造电。第二储能SOC更新约束里存在充放电同时为1的可能必须加一个二进制互斥约束。第三如果目标函数里有弃风弃光惩罚那弃掉的量不能小于零也不能大于总预测功率这两端都要封死否则优化器会给一个负的弃电值来平账。互斥约束的常见写法是充电标识加放电标识小于等于1再配两个非负连续变量表示充/放电功率。这样充放电功率不会同时为正。我见不少初学者在这里跳坑说储能莫名其妙边充边放就是没加互斥约束。3. Matlab实现流程与代码关键片段3.1 仿真的基础数据准备先设置基础参数调度周期24小时日内滚动周期4小时实时间隔15分钟。更细的调度间隔会显著增加变量规模。在复现初期我建议先用24个时段跑通逻辑再加密到96时段验证精度。风光负荷曲线我用了一段模拟数据风速用威布尔概率密度生成光照用正弦函数叠加云层扰动负荷按工作日形状生成并人为叠加预测误差。要注意的是确定性场景下优化结果只是一个快照论文通常会用蒙特卡洛跑出多组结果取均值。我复现时生成了一百组场景每组场景都求解一次。3.2 不确定性场景生成与压缩风光预测误差用正态分布近似均值是预测值方差按预测时尺度平方增长。有了误差之后还需要对场景进行缩减不然一百个场景丢进优化器模型变量数量爆炸。我用的场景缩减方法是概率距离聚合简单说就是把距离相近的场景合并成一个权重叠加。你也可以直接调用统计工具箱里的pdist和层次聚类效果足够。场景生成之后重点是把每个场景对应成优化模型里的一组参数。不同场景下系统必须满足功率平衡约束但储能SOC可以在场景间独立演化还是需要全局统一取决于论文设定。通常要求“非预期性”也就是调度决策不能看到未来只根据当前和过去的信息。复现时我用滚动时域方式处理每个时段只优化当前到未来四小时决策变量中只有当前时段的动作会被执行。3.3 求解器的选择与代码组织模型规模比较大场景一多就变成大规模MILP问题。我用Matlab的Yalmip建模底层接的求解器是我本科毕设时的老朋友了现在也还在用。如果实验室没有商业求解器可以用开源的Cbc或Gurobi试用版处理规模稍小的模型也够用。求解速度受几个因素影响二进制变量数量、大M系数、场景压缩率。我处理方式是发电机启停变量设为二进制储能充放电变量也设为二进制其余连续变量用sdpvar定义。% 建立模型 yalmip clear T 24; % 决策变量 u_gen binvar(T, 1, full); % 火电机组启停 p_gen sdpvar(T, 1, full); % 火电出力 p_ch sdpvar(T, 1, full); % 储能充电功率 p_dis sdpvar(T, 1, full); % 储能放电功率 soc sdpvar(T1, 1, full); % 荷电状态 % 目标函数 cost_gen sum(alpha .* p_gen beta .* u_gen); cost_grid sum(price_buy .* p_buy - price_sell .* p_sell); cost_degrad sum(lambda_seg .* soc(2:end)); objective cost_gen cost_grid cost_degrad; % 约束堆叠 constraints []; for t 1:T constraints [constraints, soc(t1) soc(t) eff_ch*p_ch(t) - (1/eff_dis)*p_dis(t)]; constraints [constraints, u_ch(t) u_dis(t) 1]; 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, p_gen(t) p_dis(t) p_wind(t) p_load(t) p_ch(t) p_buy(t)]; end % 求解 ops sdpsettings(solver, gurobi, verbose, 1, debug, 1); sol optimize(constraints, objective, ops);代码里最需要注重的是SOC状态转移方程。我特意在soc定义上多留了一个索引soc(t1)表示时段末的状态这样初始SOC可以单独设置不会和目标函数里的index错位。另外一个很关键的点是充电和放电标识必须和功率对应起来否则充电功率大于0但标识为0约束就直接崩了。3.4 结果可视化与指标计算优化跑完之后不能只看目标函数值。我习惯把SOC曲线、充放电功率、机组出力堆在一张图里再单独画一张成本和衰减成本随调度时长变化的面积图。这样可以直观看出加了衰减成本后SOC曲线是不是真的更平缓充放电是不是更少出现满充满放的尖峰。指标计算上论文一般看三样东西弃风弃光率、系统备用缺失概率、储能的循环寿命估计。循环寿命估计要结合实际SOC曲线一个个循环去数或者用雨流计数法。我自己写了个粗略的辅助函数统计SOC越过阈值的次数给个近似值。这样做出来的表格和论文里的曲线对得上才有说服力。4. 复现中的常见问题与排查技巧4.1 求解器报错“infeasible”怎么排查我最开始也踩了这个坑。模型一加储能互斥约束立刻显示不可行。排查思路分三步第一步把衰减成本项去掉只保留功率平衡约束确认模型有没有解第二步加入储能约束但只保留SOC转移看到底是哪个约束把可行域边界逼到空集。第三步检查SOC初始值如果初始SOC定死是0而充电功率被限制为最小出力为正那当然无解。一个非常实用的调试技巧把目标函数改成极简单的最小化常数然后让求解器输出“修复不可行集”Yalmip的debug模式可以直接指向冲突条目。用这个方式我不夸张地说对排查问题能节省一半时间。4.2 衰减模型参数不知道怎么标定论文里不同品牌电池、不同温度、不同放电倍率下的循环寿命参数全不一样。你先确定仿真对象是哪类电池锂离子还是铅酸然后去厂商数据手册或公开论文里找真实曲线。标定时要把曲线转成“放电深度-等效循环次数”表格再用线性差值填入我们模型的分段表达式。没有数据的话可以按泄漏数基准整定比如100%放电深度对应4000次循环50%放电深度对应8000次循环浅充浅放按比例延寿。我实操中碰到的另一个情况是衰减成本单位设置不当导致衰减成本在总成本里占比极小。最后我用对数坐标观察成本变化把衰减系数调高了100倍才让SOC曲线变得平滑。4.3 代码运行速度和内存占用过高怎么办处理96时段加几百个场景时Yalmip会瞬间生成几万个变量求解器内存很容易爆。我的解决办法是分层求解先跑日前主问题只包含机组组合和储能SOC计划再把日内滚动作为一个子问题子问题的决策逐时段带入。这种方法损失一点全局最优性但速度快很多。严格复现顶刊结果时可以采用Benders分解或把场景并行求解。Matlab的parfor在这里很好用每个场景单独调一个solve最后汇总。还有一个小tip将矩阵运算矢量化少用for循环。我算过一笔账把约束构建从一堆For循环改成矩阵形式求解初始化时间从3分钟降到20秒推荐所有变量下标常量都提前算好向量再一次性构建sdpvar表达式。4.4 结果图对不上原文需要如何调整论文里给的图一定是某个特定典型场景下的结果而不是所有场景的平均值。如果复现出来的功率曲线跟原文差异大先看看时间尺度是否一致、负荷是否归一化、数据是否经过光滑处理。很多时候数据根本不公开你需要根据图形反推曲线特征。比如负荷曲线峰值出现时段、风电的反调峰特性这些都是可以人为设置的。对比SOC曲线时如果形状类似但时间轴偏移很可能是时基设置有问题把t0对齐到调度首时段的凌晨很多论文会把横轴改成调度序号要注意区分。唯一准则就是确保坐标轴含义和原论文一致。5. 复现结果分析框架与扩展方向5.1 结果分析框架从三条曲线看系统平衡结果分析我建议直接从下面表格里的三项指标入手。指标定义说明灵活供给率实际灵活出力占系统最大灵活需求的比例越高越不容易甩负荷储能容量衰减率调度周期末SOH与初始健康度的差值越低电池越耐用综合经济指标总运行成本与衰减成本之和越低策略综合越经济复现时先跑一组“基准场景”再跑“加衰减建模的场景”对比表格中的数值。加了衰减建模后灵活供给率可能会有小幅下降但衰减率大幅下滑综合经济指标更好。这种“以略微降低短时灵活性换取更优全生命周期成本”的结果正好呼应标题里“平衡灵活性与储能成本”的主题。配上SOC曲线对比图能明显看出不加衰减成本时SOC上下波动剧烈经常冲顶打底加了衰减成本后SOC曲线平滑很多充放电次数减少但总电量调度序贯性更好。把这两条曲线放在同一张图里听我一句劝审稿人爱看这个。5.2 可以怎么做扩展和落地复现不是终点能把模型用到实际工程里才算掌握。我刚跑通的时候做的第一个扩展是把固定储能替换成“锂电液流电池”的混合储能。液流电池循环寿命长、成本高锂电池能量密度高、成本适中。这个扩展只需要把同一套模型里的储能参数替换成两种电池分别加上容量和功率约束目标函数里再增加一个投资成本和运维成本项就能对比不同技术路线的经济性。另一个可做的落地方向是接入现货市场的实时电价。虚拟电厂的外特性是参与市场交易而不是单纯自我平衡。把价格信号从固定参数变成随外部场景变化的输入模型就能回答什么时候买电、什么时候卖电、储能在什么时段套利最划算。我拿模拟的真实市场电价曲线试过调度结果确实出现了“低充高放”的经典套利形状和电网真实运行经验一致。还可以把电动汽车有序充电并进来。电动汽车既是负荷又是一个移动储能电池组但充放电决策受到用户出行需求的强约束。加入这类约束后还需要考虑车主最低电量需求以及集中式调度的可接受性。这一条扩展到实际工程项目中很有说服力模型脑子里的“灵活性来源”从一个固定储能扩展成移动分布式储能。6. 写在最后复现论文这件事方法和运气各占一半。方法在于系统性拆解模型把目标函数、约束、参数、求解步骤一层层剥离运气在于不同求解器、不同版本、不同机器对同一模型的执行表现完全不同。我个人经验是先把退化场景跑通比如只有两小时、三个节点、一台机组验证逻辑无误再逐步加节点数和时间尺度。这个策略帮我避开了绝大多数“模型正确但算不出来”的尴尬。最后再分享一个小技巧每次修改模型都要把改动前后目标函数和约束数打印出来记录到一个文本里。我复现期间最喜欢用Yalmip自带报告来跟踪每个约束的数量变化。出现问题时对照前后差异比重新一行行盯代码高效得多。只要把这个套路走完整不管是这篇虚拟电厂顶刊还是其他储能优化论文复现难度都能降一个数量级。希望大家都能顺利跑通自己的仿真。