储能参与调峰的配置与经济性分析Matlab实现全流程

论文里那句“系统总成本最小”看着简单，真落到Matlab代码里才发现到处都是坑。我去年完整复现过一个参与调峰的储能配置与经济性分析算例，从容量定容、时序调度到经济性评估，前前后后折腾了两周多，这里把整个过程中的模型设计、代码实现和踩坑记录整理出来，给准备做类似课题的朋友一条相对顺的路。

这篇文章适合三类人：一是正在复现储能配置类论文的研究生，二是做新能源消纳或电网调峰项目的工程师，三是想用Matlab快速验证“储能到底该配多大”的规划人员。下面内容不涉及具体某篇论文的照搬，而是把这类问题通用的建模思路、代码骨架和计算细节拆开讲清楚。

1. 储能参与调峰到底在优化什么：先把问题翻译成数学语言

1.1 调峰问题的本质与储能的角色

电网的净负荷曲线是指“用户负荷扣除新能源出力之后”的剩余负荷。当光伏装机比例上来以后，午间净负荷可能被压得很低，傍晚又迅速爬升，形成所谓的“鸭子曲线”。这时候常规火电机组要么得压出力到很低的调峰深度，要么就得频繁启停，无论哪种都会带来煤耗上升和设备损耗。

储能在调峰场景里做的事情很朴素：负荷低谷时段充电，把电能临时存起来；负荷高峰时段放电，把电能还回去。但“朴素”不等于“简单”，因为储能系统的额定功率和额定容量是两个不同的物理量——功率决定你一小时能充放多少，容量决定你总共能存多少。这两个量如果配得不匹配，要么功率足够但容量不够、高峰期放不出那么多电，要么容量够但功率太小、充放电速度跟不上系统调节需求。配置优化就是要把这两个量放在一起决策。

1.2 决策变量分两层：配置层和调度层

储能配置问题里，决策变量天然分成两类，这一点在建模时一定要分清楚：

• 配置层变量：储能的额定功率 P_ess（MW）、额定容量 E_ess（MWh），有些文章还会把接入节点也作为变量，但大多数调峰经济性分析都是单节点或典型系统模型，接入节点可以在灵敏度分析里再做。
• 调度层变量：每一个时段 t 的充电功率 P_ch(t)、放电功率 P_dis(t)、荷电状态 SOC(t)，以及充放电状态标志 u(t)。

两层变量不是独立的。调度层变量会在约束里反复用到 P_ess 和 E_ess，比如“充电功率不能超过额定功率”“SOC不能超过容量上限”，所以配置优化本质是一个把容量大小和运行策略同时求解的联合优化问题。多数EI论文用的是单层优化一次性求出结果，这也是Matlab代码实现里最主流的方式。

1.3 目标函数与约束条件怎么取舍

目标函数通常是系统综合成本最小，典型形式是：

min C_total = C_inv + C_om + C_fuel + C_penalty

其中：

• C_inv 是储能投资成本的年值化结果（用资金回收系数把一次性投资折算到每年）
• C_om 是年运维成本
• C_fuel 是常规机组发电燃料成本
• C_penalty 是弃电惩罚项或切负荷惩罚项

约束条件这部分是建模的关键，我列一下调峰储能模型里必须有的几类约束：

有一个取舍经验：很多新手会把SOC递推写成 SOC(t)=SOC(t-1)+[P_ch(t)·η_ch-P_dis(t)/η_dis]·Δt/E_ess，这没问题，但要注意单位统一。P_ch 单位是MW，Δt 单位是小时，乘出来是MWh，除以 E_ess（MWh）以后才是无量纲的SOC变化量。一旦单位不一致，结果里SOC会出现奇怪的跳变甚至超限，这个问题我后面专门讲。

目标函数里，C_penalty 弃电惩罚项值得单独强调。储能配置优化如果完全不考虑弃电成本，结果往往会偏保守——因为系统宁可让风机弃掉也不愿意花钱配储能。加入弃电惩罚后，优化器才会在“新建储能投资”和“弃电损失”之间做真正的权衡，这正是“参与调峰”这个场景的核心。我在复现时的做法是在约束里把 P_w(t) 拆成“实际消纳风电”和“弃风功率”两个变量，后者进入目标函数的惩罚项。

2. Matlab代码实现：从参数表到可复现的项目骨架

2.1 项目文件和目录组织

复现这种优化类问题，我强烈建议一开始就把代码按模块拆开，不要全部堆在一个脚本里。因为储能配置问题通常要反复调参、换场景、做灵敏度分析，一个两三百行的单一脚本到后期会非常痛苦。

我的目录结构是这样：

storag_config/ ├── 1_main.m % 主程序，串联整个流程 ├── 2_load_data.m % 数据读入与预处理 ├── 3_build_model.m % 构建优化模型 ├── 4_solve_model.m % 调用求解器 ├── 5_econ_analysis.m % 经济性指标计算 ├── 6_plot_results.m % 结果可视化 ├── data/ │ ├── load_profile.xlsx % 负荷时序数据 │ ├── wind_profile.xlsx % 风电时序数据 │ └── price_profile.xlsx % 分时电价数据 └── result/ └── (输出结果与图)

每个文件只干一件事，主程序里用函数调用串联。这样最大的好处是：当你把24时段算例换成8760时段全年算例时，只需要改数据文件和模型里的时间维度，不需要重写整个程序。

2.2 负荷和风电时序数据的读入与预处理

数据读入这块我用的是readtable，它比xlsread更稳，特别是在读取含时间戳和多列数据的Excel时。

%% 读取负荷、风电、分时电价数据 load_data = readtable('data/load_profile.xlsx', 'VariableNamingRule', 'preserve'); wind_data = readtable('data/wind_profile.xlsx', 'VariableNamingRule', 'preserve'); price_data = readtable('data/price_profile.xlsx', 'VariableNamingRule', 'preserve'); P_load = load_data.Load_MW; % 负荷序列，MW P_wind = wind_data.Wind_MW; % 风电序列，MW price = price_data.Price; % 分时电价，元/MWh T = length(P_load); % 时段数 dt = 1; % 步长，小时

注意VariableNamingRule这个参数，新版本Matlab如果不设它，遇到中文变量名或特殊字符会直接报错或者改名，这个坑我在R2023a以前遇到过好几次。另外建议数据文件里统一用英文列名，中文列名虽然能读，但后面索引容易出现编码问题。

时序数据的预处理通常包含三件事：

1. 缺失值处理：用前一天同时段数据填补，或者线性插值。
2. 数据缩放：如果原始数据是标幺值或百分比，要乘以基准值还原成MW。
3. 时间对齐：负荷、风电、电价三组数据必须严格逐时段对齐，最保险的方法是在Excel里就排成同样的行顺序。

我习惯在读完数据后加一个自检：

assert(length(P_load) == length(P_wind) && length(P_load) == length(price), ... '负荷、风电、电价数据长度不一致，请检查数据文件');

这一段小代码看着不起眼，但能避免后面模型求解时因为维度不匹配出现的各种莫名其妙错误。

2.3 用YALMIP把储能配置模型搬进Matlab

YALMIP是Matlab里做优化建模的利器。它的核心思路是让你用接近数学公式的方式定义变量、写约束，然后统一交给底层的求解器去算。对于储能配置这种混合整数线性规划（MILP）问题，YALMIP写起来比直接用intlinprog的矩阵形式简洁得多。

先定义决策变量：

%% 决策变量定义 % 配置层：储能额定功率和额定容量 P_ess = sdpvar(1, 1); % 额定功率 MW E_ess = sdpvar(1, 1); % 额定容量 MWh % 调度层：时序变量 P_ch = sdpvar(T, 1); % 充电功率 P_dis = sdpvar(T, 1); % 放电功率 u = binvar(T, 1); % 充放电状态，1充电/0放电 SOC = sdpvar(T+1, 1); % SOC序列，多一个初始时刻 P_g = sdpvar(T, 1); % 常规机组出力 P_w_use = sdpvar(T, 1); % 实际消纳风电功率 P_curt = sdpvar(T, 1); % 弃风功率

这里 SOC 定义成 T+1 个变量，是因为递推公式里需要初始时刻和最终时刻两个端点的SOC。写的时候我就直接使用这个约定，不要省那一个变量，否则后面写周期约束时索引不好对。

再写约束：

%% 约束条件 Constraints = []; % 功率平衡约束 for t = 1:T Constraints = [Constraints, P_g(t) + P_w_use(t) + P_dis(t) == P_load(t) + P_ch(t)]; Constraints = [Constraints, P_w_use(t) + P_curt(t) == P_wind(t)]; end % 储能运行约束 for t = 1:T % 充放电功率不能超过额定功率，且不能同时充放 Constraints = [Constraints, 0 <= P_ch(t) <= P_ess * u(t)]; Constraints = [Constraints, 0 <= P_dis(t) <= P_ess * (1 - u(t))]; % SOC递推 Constraints = [Constraints, SOC(t+1) == SOC(t) + (eta_ch * P_ch(t) - P_dis(t) / eta_dis) * dt / E_ess]; end % SOC上下限与边界周期约束 Constraints = [Constraints, SOC_min <= SOC <= SOC_max]; Constraints = [Constraints, SOC(1) == SOC_init]; Constraints = [Constraints, SOC(T+1) == SOC_init]; % 周期约束 % 机组出力上下限 Constraints = [Constraints, P_g_min <= P_g <= P_g_max]; % 储能规模下限 Constraints = [Constraints, P_ess >= 0, E_ess >= 0];

目标函数：

%% 目标函数 C_inv = (c_p * P_ess + c_e * E_ess) * CRF; % 投资年值 C_om = om_rate * (c_p * P_ess + c_e * E_ess); % 年运维成本 C_fuel = sum(fuel_price * P_g * dt); % 燃料成本 C_curt = curt_price * sum(P_curt * dt); % 弃风惩罚 Objective = C_inv + C_om + C_fuel + C_curt;

求解：

%% 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'showprogress', 1); optimize(Constraints, Objective, ops); % 提取结果 P_ess_opt = value(P_ess); E_ess_opt = value(E_ess); SOC_opt = value(SOC); P_ch_opt = value(P_ch); P_dis_opt = value(P_dis);

YALMIP的好处是约束写法几乎和数学公式一一对应，代码不容易把约束写错。注意sdpvar(T+1,1)在约束循环里使用SOC(t+1)时不会越界，因为 SOC 有 T+1 个元素，t 从 1 到 T，刚好用到 SOC(2)到SOC(T+1)。

2.4 intlinprog与商用求解器在调峰场景下的取舍

很多朋友可能一开始不想装Gurobi，想直接用Matlab自带的intlinprog。我的经验是：24小时或48小时的典型日算例，intlinprog完全够用，求解时间基本在几秒到几十秒；一旦把时间尺度拉长到全年8760小时，变量数量会大幅增加，intlinprog的求解时间会指数级增长，有时候几个小时都跑不完。

我个人在复现阶段会装一个Gurobi，理由有三点：

• MILP求解速度比intlinprog快一到两个数量级，8760小时的算例差距非常明显。
• 求解数值稳定性更好，不容易出现“明明有可行解但求解器报告不可行”的误判。
• YALMIP对Gurobi支持完善，改动量很小。

如果确实不方便装商用求解器，也有一个折中方案：先用24小时典型日做配置决策和灵敏度分析，再把最优容量固定下来，用intlinprog对全年做运行模拟校核。这样既能得到较可信的容量结果，又不会因为求解器性能卡住整个项目。

求解器安装后要在YALMIP里确认路径有效，可以用yalmiptest检测。有一个常见的坑是：安装完Gurobi后直接运行optimize，YALMIP提示找不到求解器，原因是Gurobi的路径没有加进Matlab路径里。解决办法是在startup.m里添加addpath('C:\gurobi1000\win64\matlab')，或者在Gurobi安装目录下运行gurobi_setup脚本。

3. 经济性分析：储能不是省了电费就万事大吉

3.1 费用项到底有哪些，别漏算

储能经济性分析最怕两件事：一是收益算过头，二是成本算漏项。先看成本侧，完整清单是这样的：

• 初始投资：PCS（功率变换系统）按功率计价，单位元/kW；电池本体按容量计价，单位元/kWh。此外还有土建、接入、并网等一次性费用，通常按初始投资的百分比估算。
• 年运维成本：固定运维按额定容量算（元/kWh/年），可变运维按年充电电量算（元/MWh）。
• 电池置换成本：锂电池寿命通常在8000到10000次循环，或8到12年。项目周期如果超过电池寿命，要计入置换成本。这个很多人会漏掉。
• 折现率：把未来各年的现金流折算到当前。电力项目折现率一般取6%到8%。

年值化的数学处理是资金回收系数（CRF）：

CRF = r·(1+r)^Ny / [(1+r)^Ny - 1]

其中 r 是折现率，Ny 是项目年限。用CRF把一次性投资折算成等额年值，才能和每年的运行收益公平对比。

3.2 收益项怎么识别，别重复计算

储能的收益来源在调峰场景里主要有四块：

1. 峰谷套利收益：低谷充电、高峰放电，赚取峰谷电价差。这是最直接的收益。
2. 容量补偿或调峰辅助服务收益：储能参与调峰后，为系统提供了调节容量，部分地区或市场会对这种辅助服务付费。
3. 减少弃电带来的收益：如果系统本来要弃风弃光，储能消纳这部分电量后，相当于挽回了这部分电量的价值。
4. 延缓新建调峰机组的投资：储能容量替代了一部分新建机组的容量需求，这部分“避免投资”也可以算作收益。

在计算时要注意：峰谷套利和容量补偿不能同时全额计算，因为储能同一时段放电，既赚了电量价差，又提供了容量服务，这两块如果都按全额算，等于同一笔放电被算了两次。论文里常见的处理方式是：套利收益按实际峰谷价差计算，容量收益按容量市场或辅助服务市场规则单独核算，两者属于不同结算规则，但要在论文里说明口径。

3.3 NPV、IRR和回收期的Matlab实现

经济性评估我单独写成一个函数，输入是储能配置结果和价格参数，输出是净现值、内部收益率和回收期。

function [NPV, IRR, payback] = econ_eval(P_ess, E_ess, P_dis_opt, P_ch_opt, price, params) % 参数 r = params.r; % 折现率 Ny = params.Ny; % 项目年限 c_p = params.c_p; % 功率成本 元/kW c_e = params.c_e; % 容量成本 元/kWh om_fixed = params.om_fixed; % 固定运维 元/kWh/年 om_var = params.om_var; % 可变运维 元/MWh % 初始投资 C_init = c_p * P_ess * 1000 + c_e * E_ess * 1000; % 单位换算成元 % 年运行收益（峰谷套利按全年模拟结果折算） annual_energy = sum(P_dis_opt * params.dt) * 365; % 假设典型日代表全年，MWh/年 annual_revenue = annual_energy * params.avg_price_diff; % 平均峰谷价差 % 年运维成本 C_om_annual = om_fixed * E_ess * 1000 + om_var * sum(P_ch_opt * params.dt) * 365; % 年净现金流 CF_annual = annual_revenue - C_om_annual; % 现金流序列 CF = zeros(1, Ny); CF(1) = -C_init; for k = 2:Ny CF(k) = CF_annual; end % NPV NPV = sum(CF ./ ((1+r).^(0:Ny-1))); % IRR：令NPV=0求折现率 IRR = fzero(@(x) sum(CF ./ ((1+x).^(0:Ny-1))), 0.1); % 回收期：累计现金流首次转正的年份 cum_CF = cumsum(CF); payback = find(cum_CF >= 0, 1, 'first'); if isempty(payback) payback = Inf; end end

用典型日代表全年计算收益是一种简化做法，适合初步估算。如果要做全年精确计算，最好还是跑8760小时的运行模拟，再把全年电量代入经济性模型。另外IRR的求解用fzero时初始值设为0.1比较稳妥，IRR在合理范围时收敛很快，但如果项目亏得太厉害，fzero可能找不到零点，要加错误处理。

4. 算例验证与结果呈现：怎么证明配置结果是可信的

4.1 典型算例场景设定

算例参数不要随便拍脑袋，尽量参考同领域论文或典型的系统测试数据。我这里给一套常用的参数组合，你可以作为起点：

数据来源方面，负荷曲线可以用IEEE RTS系统的典型日数据，风电曲线可以用某个风电场的历史出力数据归一化后放大。关键是论文里要交代清楚，审稿人会盯这一块。

4.2 优化结果解读：SOC曲线和充放电行为是否合理

求解完成后先别急着看经济性指标，先检查优化结果是否符合物理直觉。我复现时第一个检查项是SOC曲线：

figure; plot(1:T, SOC_opt(1:T), 'LineWidth', 1.5); xlabel('时段/h'); ylabel('SOC'); title('储能SOC变化曲线'); grid on;

合理的SOC曲线应该满足三个特征：

1. 在谷时充电段SOC上升，在峰时放电段SOC下降，和电价曲线的联动关系清晰。
2. SOC始终在上下限之间，不会出现振荡或突变。
3. 起始和结束的SOC相等，满足周期约束。

如果SOC曲线出现“充满又放空又充满”的反复跳变，很可能是容量配置过小或者峰谷价差不规则，需要仔细核实。另外一个合理性的检查是把充放电功率曲线和负荷曲线画在一起，看看储能是不是真的在削峰填谷，而不是在随机充放。

4.3 敏感性分析：储能造价、峰谷价差的影响

配置优化不能只给一个静态结果，要做容量配置对关键参数的灵敏度分析。最常做的两个维度的扫描是储能单位造价和峰谷价差。

%% 容量成本扫描 c_e_range = [600:200:2000]; % 元/kWh E_opt_record = zeros(size(c_e_range)); for i = 1:length(c_e_range) c_e = c_e_range(i); % 重新运行优化模型 ... E_opt_record(i) = value(E_ess); end figure; plot(c_e_range, E_opt_record, 'o-', 'LineWidth', 1.5); xlabel('电池容量成本/(元/kWh)'); ylabel('最优储能容量/MWh'); grid on;

这类曲线是最有力的论文素材：横轴是成本参数，纵轴是最优配置，展示出“随着储能成本下降，最优配置容量如何增长”的规律。通常结果是一条近似线性的增长曲线，拐点出现在储能成本下降到接近峰谷价差收益的临界值时。我做这个分析时发现了一个在实际工程中很重要的规律：储能配置对峰谷价差的敏感度远高于对储能成本的敏感度，价差从400元/MWh涨到600元/MWh，最优容量能翻一倍多，而降成本带来的边际增长要平缓得多。这说明调峰储能项目的经济性瓶颈很大程度在收益价格机制上，单纯降设备成本解决不了根本问题。

5. 复现途中踩过的坑：从求解失败到论文出图

5.1 模型无解、结果NaN与不合理SOC的排查链路

复现这类优化问题，绕不开的坎就是求解器报“无解”或给出的结果全是NaN。我总结了一套排查顺序：

第一步，检查不可行约束。先把SOC上下限放宽到0和1，把机组出力范围放宽，看看模型能不能求解。如果放宽后能解，说明是某个约束条件之间的冲突。我遇到最多的是周期约束和SOC初值冲突——SOC(1)设为某个值的同时又要求SOC(T+1)也等于这个值，但充电量根本不足以让SOC在周期内回到初值，解决办法是先不加周期约束，求解一次看终值SOC在哪，再据此调整初值。

第二步，检查数值尺度。储能额定容量动辄几十上百MWh，但如果某种电价或成本参数用的是元/MWh、数百的数值，目标函数各项的量级可能相差10^6倍以上，求解器数值稳定性会变差。解决办法是对大数值参数做缩放，比如把成本单位改成万元，或者对目标函数做归一化。

第三步，检查二进制变量与连续变量的耦合。充放电互斥约束 P_ch(t)≤P_ess·u(t) 里，如果 P_ess 是决策变量，这其实是一个非线性约束（变量乘变量）。YALMIP会自动处理成线性形式吗？实际上这块要小心。如果 P_ess 是变量而 u(t) 也是变量，P_ess·u(t) 就是双线性项，不是线性约束，求解器会报错或进入非凸优化。常规做法是引入辅助变量替换，或者把 P_ess 的参数化：外层循环枚举候选功率等级，内层只优化调度变量，这是最稳妥的处理方式。我复现时发现不少新手在这里栽跟头，看论文公式以为可以直接写进YALMIP，实际上需要MILP化处理。

再补充一种排查经验：检查结果时如果发现SOC曲线出现轻微越界但求解器没报错，很可能是容差设置问题。求解器默认的MILP容差在1e-4级别，SOC约束的越界量可能被容忍。这时可以收紧sdpsettings('gurobi.MIPGap', 1e-4)或对应求解器的MIPGap参数。默认的MIPGap是1e-4还是别的值，不同版本不一样，最好显式设置。

5.2 8760小时全年算例的内存与计算时间问题

全年8760小时的算例，决策变量数量大约是：P_ch、P_dis、u、SOC、P_g、P_w_use、P_curt，各自8760个，再加上P_ess、E_ess，总共接近6万个变量，其中二进制变量8760个。这个规模对内存和求解器都有压力。

我的建议是分两步走。第一步，先把全年数据按季节性分成几个典型日场景（比如冬季、夏季、过渡季各取一个典型日），在典型日上做配置优化决策。第二步，把第一步得到的最优 P_ess 和 E_ess 固定下来，对全年8760小时只做运行模拟，校验全年收益和SOC是否始终在允许范围内。这种“配置优化用典型日，运行校核用全年”的做法在工程上足够可靠，计算复杂度也低很多。如果非要一次性做全年优化，建议至少把二进制变量 u(t) 通过“充放电时段预先给定”的方式变成线性规划，但这需要对运行策略有较强的先验判断。

5.3 论文级图表导出的细节处理

论文里需要导出的图主要是SOC曲线、削峰填谷效果对比图、灵敏度分析图、经济性指标对比表。Matlab画图默认的分辨率和字体在投期刊时往往不够，我用exportgraphics导出高清图：

exportgraphics(gcf, 'result/SOC_curve.eps', 'ContentType', 'vector');

EPS格式对应的是热词里“matlab 2025导出eps”这种搜索，说明很多人在卡这一步。新版Matlab官方推荐的是exportgraphics，比老的print -depsc更稳，并且支持矢量输出。导出PNG时设置分辨率：

exportgraphics(gcf, 'result/result.png', 'Resolution', 600);

字号方面，正文图里的文字大小建议10到12磅，坐标轴刻度标签8到9磅。中文字体问题：Matlab默认的Helvetica不支持中文，中文会变成方块，解决办法是画图时先设置：

set(groot, 'defaultAxesFontName', 'Times New Roman'); set(groot, 'defaultTextFontName', 'Times New Roman');

如果图里需要出现中文，就用'SimHei'或'Microsoft YaHei'，但投稿英文期刊时图里最好不要放中文，全部用英文标注更省事。线宽统一2磅，颜色用色盲友好配色（比如蓝、橙、绿三色），这样打印黑白版时也有区分度。

结尾想说的话

整套流程跑下来我最大的体会是：储能配置优化这个题目，数学建模只占三成，数据与代码工程占七成。跟着论文公式写YALMIP约束是最容易的部分，真正的功夫在数据处理、单位统一、求解器配置、结果合理性校验这些不起眼的细节上。建议你先拿24小时典型日把整条链路跑通，确认SOC曲线和经济性指标都合理了，再逐步扩展到全年算例，这样排查问题会顺手很多。最后再分享一个我自己一直在用的小技巧：遇到任何奇怪的求解结果，先把约束逐条注释掉、逐条加回来，找到影响结果的那条约束，比盯着整个模型猜要高效得多。