我做这个课题的时候最直观的感受就是风电、光伏、储能调度这块资料不少但大多数要么只讲理论公式要么直接甩一个封装好的simulink模型很难看到一套从数学建模到Matlab代码再到结果分析的完整链路。这次正好借着一个包含电池储能和废弃矿井小型抽水蓄能的互补调度项目把整条路线用Matlab从头到尾走了一遍。这篇文章就是我个人的完整记录包含各单元建模、优化调度模型、核心代码、调试经验以及一些常规文档里不会写的心得体会。这个课题的价值在于风电和光伏出力天然波动负荷也在不断变化单一储能手段总有一些力不从心的地方。电池响应快但容量贵抽水蓄能容量大但选址苛刻。把电池和抽蓄组合起来、并且把抽蓄选址放宽到废弃矿井改造是一个既有学术价值又贴近工程实际的思路。如果你正在做新能源微电网、源网荷储一体化、储能优化调度相关的研究这篇文章可以当一份带代码的参考。不需要很强的优化理论功底我会尽量把每个关键选择背后的原因讲清楚。1. 项目背景与整体设计思路1.1 为什么风电光伏必须配储能风电、光伏的出力特征是间歇性强、随机波动大而且峰谷时段和负荷曲线经常错位。风大的时候可能正是半夜负荷低谷光伏大发的中午负荷可能还没上来。如果这些电直接并网系统需要准备大量的旋转备用去平衡波动否则电网频率和电压质量都会受影响。配上储能之后逻辑就顺了多余的清洁电能先存起来等风小、光弱、负荷高峰的时候再放出来。这就是移峰填谷加平滑波动两个功能的叠加。但储能本身也有分类思维电池储能响应速度可以从毫秒级到分钟级适合做短时功率平抑、调频和快速爬坡支撑但容量做到很大之后成本呈线性上涨而且循环寿命有限天天深度充放几年就得换。抽水蓄能则恰好反过来单机容量大、可长时间运行、寿命以几十年计但响应速度慢、启动流程复杂而且常规抽蓄需要合适的地形建设上下两个水库选址极其受限。所以电池抽蓄的组合本质上是用两类时间尺度完全不同的储能去配合电池负责高频次、短时、快速响应抽蓄负责大容量、长周期、跨时段的能量搬移。这种搭配在风光基地、源网荷储一体化等场景里越来越常见也是我这个课题最开始立意的出发点。1.2 废弃矿井小型抽水蓄能一个被低估的方案常规抽蓄选址难但很多资源枯竭型矿区恰恰有大量废弃巷道和竖井。把地面水池作为上水库、把井下巷道空间作为下水库就能改造成一个小型抽水蓄能电站。相比新建常规抽蓄这种方案土建成本低、不占用新的土地资源还能给矿区带来产业转型方向。德国、美国都有相关示范工程国内这几年也在做前期研究和试点。不过在调度建模里废弃矿井抽蓄比常规抽蓄多几个现实约束容量小。矿井巷道容积决定了库容通常只能做成小型抽蓄单机功率可能也就几兆瓦到几十兆瓦。水头可能不稳定。水位变化会带来净水头和发电效率的变化简化建模可以用平均水头加固定效率近似精细研究要考虑H-Q特性曲线。地下空间安全约束。运行水位上下限要严格限制防止巷道结构失稳或者涌水倒灌。响应慢。机组启动、工况转换都要时间调度模型里通常用爬坡约束和最小运行/停运时间约束来表达。我对矿井抽蓄的处理方式是先把上下水库抽象成两个水箱上库对应地面蓄水池下库对应井下巷道空间两者之间的水量守恒用线性方程描述水头影响先简化成固定效率后期如果有实际机组参数再升级成变效率模型。这样既保证调度模型可解又不丢库容有限、不能同时抽发这些核心物理约束。1.3 整体技术路线整个研究的技术路线分四步输入数据准备风电、光伏出力和负荷的时序曲线可以用历史数据、典型日场景也可以通过随机模拟生成。单元建模把风电、光伏、电池、抽蓄分别写成状态量到出力的数学模型。互补调度优化以系统经济运行和新能源消纳为双重目标构造优化模型用Matlab调用求解器求解。结果评估统计弃风弃光率、购电成本、储能利用率用图表对比不同配置方案。这四步里最容易出问题的是第2步和第3步的衔接。模型写太薄丢掉关键约束结果没有物理意义模型写太厚约束矩阵庞大求解器跑不动。后面我会重点讲这个平衡怎么把握。2. 各单元数学模型搭建2.1 风电出力模型风电调度模型通常分两步先产生风速序列再用风机功率曲线映射到出力。风速序列我用威布尔分布随机模拟加时间相关性修正。威布尔分布的两个参数形状参数k和尺度参数c可以用历史风速数据通过极大似然估计拟合出来。比如某个内陆风电场k约等于2.0、c约等于7.5 m/s这是比较典型的参数。但这里有个细节容易被忽略如果直接用独立随机序列作为风速相邻时段风速跳变会很剧烈映射出来的风电出力是高频抖动的白噪声调度优化出来的结果会显得非常假。我加了简单的一阶自回归滤波让风速具备时间连续性rng(42); T 24; % 时段数单位1小时 k 2.0; c 7.5; % 威布尔分布参数 v_iid wblrnd(c, k, [T, 1]); % 独立随机风速 rho 0.85; % 一阶自回归系数 v_ar zeros(T, 1); v_ar(1) v_iid(1); for t 2:T v_ar(t) rho * v_ar(t-1) (1 - rho) * v_iid(t); end v_ar max(v_ar, 0); % 风机功率曲线映射切入风速、额定风速、切出风速 v_in 3; v_r 12; v_out 25; Pr 1.5; % 单位MW Pw zeros(T, 1); idx_linear (v_ar v_in) (v_ar v_r); idx_rate (v_ar v_r) (v_ar v_out); Pw(idx_linear) Pr * (v_ar(idx_linear) - v_in) / (v_r - v_in); Pw(idx_rate) Pr;实际项目里如果有风电场的实测出力序列可以直接读数据不一定要走风速映射。但代码里保留一个数据模式/模拟模式的开关能让没有实测数据的人也能复现整个流程这个设计很值得保留。2.2 光伏出力模型光伏出力本质上取决于水平面总辐照度和组件温度。工程上常用的简化表达式是P_pv(t) G(t) * A * η_pv * [1 - β * (T_c - T_ref)]其中G是辐照度A是组件总面积η_pv是标准测试条件下的效率β是温度系数单晶硅组件一般在0.004左右/°CT_c是组件温度T_ref是参考温度25°C。辐照度模拟我一般用两种办法一种是直接取典型日辐照曲线形状类似钟形再加上云层遮挡的随机扰动另一种用Beta分布采样因为Beta分布可以较好地刻画辐照度的偏态特性。调度研究里常见做法是在晴空模型基础上叠加天气类型随机系数比如多云天折减系数在0.3到0.7之间随机。一个必须注意的点光伏出力序列和负荷、风电的时序相关性不能忽略。比如同一个地区白天光照强的时候可能恰恰是负荷高峰的下午时段夜晚光伏出力为0。如果完全独立随机生成三类曲线调度结果虽然能跑通但经济性指标和真实情况差了十万八千里。正确做法是先确定基准负荷曲线再让风电和光伏围绕各自典型日形状做随机扰动。H 24; t_h (0:H-1); G_clear 800 * max(sin(pi * (t_h - 6) / 12), 0).^2; % 6点到18点的晴空辐照度 cloud 0.5 0.5 * rand(H, 1); % 云量系数0-1 G G_clear .* (0.7 0.3 * cloud); A 5000; % 组件总面积单位m2 eta_pv 0.18; beta 0.004; Tc 30; Tref 25; Ppv G .* A * eta_pv .* (1 - beta * (Tc - Tref)) / 1e3; % 单位kW转MW注意这里的单位换算。调度模型里所有功率量纲最好统一成MW否则后面构造约束矩阵时很容易出现十倍、千倍的数值差直接导致求解器精度问题。2.3 电池储能系统模型电池储能的核心状态量是SOC本质上是一阶差分方程SOC(t) SOC(t-1) η_ch * P_ch(t) * Δt / E_bat - P_dis(t) * Δt / (η_dis * E_bat)其中η_ch、η_dis分别是充电、放电效率E_bat是电池额定容量Δt是调度步长。充电功率和放电功率物理上互斥处理互斥关系有两条路加二元变量u_ch u_dis ≤ 1问题变成混合整数规划求解严谨但计算量增大。不加二元变量靠目标函数自然排斥只要充电成本和放电成本都为正同一时段同时充放电会使功率平衡等式两端同时增加不可能让目标函数更小。实际项目中不加互斥约束跑出来通常也合理但存在风险如果某些价格信号极端比如充电补贴高到充电成本为负模型就真可能出现同时充放电的假象。所以我的建议是研究性质的项目老老实实加二元变量用MILP求解只是快速估算或教学演示可以省略。电池还要考虑SOC上下限、充放电功率上限。另外每完成一次充放电循环电池容量会有衰退这在日前的调度模型里常简化成固定充放电循环成本加在目标函数里用来体现减少电池深度充放的意图。E_bat 20; % 额定容量 MWh SOC_min 0.2; SOC_max 0.9; P_ch_max 5; P_dis_max 5; % MW eta_ch 0.95; eta_dis 0.92; C_bat 12; % 充放电循环成本 元/MWh2.4 废弃矿井小型抽水蓄能模型抽蓄模型和电池在数学形式上很像都是储能罐思想但物理约束差异很大能量状态变量是上库水量或等效电量而不是SOC。抽水和发电两个方向效率不同启动过程有时间迟滞。库容受上下水库水位限制不是简单一个SOC区间。水头变化影响出力范围简化模型用固定水头加固定效率精细模型要引入水头-功率耦合。我用的简化模型是E_res(t) E_res(t-1) η_pump * P_pump(t) * Δt - P_gen(t) * Δt / η_gen其中E_res是抽蓄电站的等效电量P_pump是抽水功率P_gen是发电功率η_pump、η_gen分别是抽水和发电效率。等效电量的思路是把上库水量乘上重力势能换算成MWh这样就能把抽蓄和电池放进同一个储能母线框架里比较。在容量不大、水头变化不剧烈的废弃矿井小型抽蓄场景下这个简化是合理的。和电池类似抽水和发电互斥、等效电量有上下限、功率有上限。但抽蓄还多两个约束最小运行/停运时间约束因为机组启停慢日抽发循环次数限制因为频繁启停损伤机组。如果研究重点是互补调度逻辑可以先放下这层用软约束在目标函数里加启停惩罚项代替。E_res_max 60; % 等效电量上限 MWh E_res_min 6; % 等效电量下限对应死水位 eta_pump 0.75; eta_gen 0.82; P_pump_max 4; P_gen_max 5; % MW3. 优化调度模型目标函数与约束体系3.1 目标函数经济性和消纳率如何平衡互补调度的核心目标通常不是单一的成本最低也不是单一的消纳最大而是一个组合目标。我这里选的是最常见的加权和形式min F Σ C_buy(t) * P_grid(t) * Δt Σ C_bat * (P_ch(t) P_dis(t)) * Δt Σ λ_w * Curt_w(t) * Δt λ_pv * Curt_pv(t) * Δt Σ λ_pump * (P_pump(t) P_gen(t)) * Δt第一项是购电成本第二项是电池循环成本第三项是弃风弃光惩罚第四项是抽蓄运行成本。这个目标里真正需要刻意设计的是弃风弃光惩罚系数λ_w、λ_pv。设得太大模型会不计成本地消纳新能源哪怕购电价格很低也选择不买电设得太小模型宁可弃风弃光也不愿意让储能设备动作消纳率惨不忍睹。一般标幺化后λ取购电电价的1.5到3倍比较合适能体现优先消纳新能源的政策导向。3.2 约束条件功率平衡是灵魂约束分四类功率平衡约束这是所有调度模型的核心P_w(t) - Curt_w(t) P_pv(t) - Curt_pv(t) P_dis(t) P_gen(t) P_grid(t) L(t) P_ch(t) P_pump(t)左边是供电源风电实际出力、光伏实际出力、电池放电、抽蓄发电、电网购电右边是需求侧负荷、电池充电、抽蓄抽水。所有量纲统一为MW。储能单元自身约束电池SOC差分方程和上下限、抽蓄等效电量差分方程和上下限、充放功率上下限、互斥关系。外购电约束P_grid(t)有最大允许值模拟联络线功率极限。爬坡约束抽蓄相邻时段功率变化率有限电池一般不用爬坡约束因为电池爬坡能力远快于调度步长。这里插一句容易踩的坑功率平衡方程里每一类电源的实际出力是决策变量不是预测值。比如风电预测是10MW但弃风了2MW那么平衡方程左边只能用实际出力8MW。很多初写调度代码的人直接把预测值当常数塞进平衡方程结果弃风变量怎么设都对功率平衡没影响这就是模型逻辑出问题的典型症状。3.3 求解算法为什么我选MILP由于引入了电池充放互斥、抽蓄抽发互斥这些二元变量这个问题天然是混合整数线性规划。Matlab自带的intlinprog可以直接求解中小规模问题。我的算例是24时段、单元数量有限的单场景决策变量200个左右intlinprog几秒到十几秒就能出结果完全够用。如果是多场景比如8760小时全年时序或上千个随机场景或者要考虑非线性水头效率曲线那就得用YALMIP作为建模层后端接Gurobi或Cplex求解效率会高很多。YALMIP的好处是符号化建模代码可读性和可维护性明显优于裸写intlinprog的大矩阵。我在项目里先用YALMIP搭第一版验证模型逻辑正确后再用intlinprog重写一版加速两版结果对比一致后以intlinprog版本作为交付代码。4. Matlab代码实现从建模到求解的完整流程4.1 模块化代码结构整个项目的代码结构如下main.m % 主程序调度求解总入口 data_prepare.m % 数据生成与读取 unit_model.m % 各单元参数初始化 build_schedule.m % 构建优化模型并求解 plot_results.m % 画图与结果分析代码分层有一个很现实的好处换数据、换参数、换求解器时不需要动核心逻辑。我那会儿为了对比有抽蓄/无抽蓄电池与抽蓄配比不同等好几组方案如果所有代码写在一个脚本里每次都要复制粘贴大段代码还容易改错变量。模块化之后只是改data_prepare里的参数主程序一行不用动。4.2 核心代码段YALMIP建模版本用YALMIP建模最直观代码和数学表达式几乎一一对应% 变量定义 P_grid sdpvar(T, 1); P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); u_ch binvar(T, 1); u_dis binvar(T, 1); SOC sdpvar(T, 1); P_pump sdpvar(T, 1); P_gen sdpvar(T, 1); u_pump binvar(T, 1); u_gen binvar(T, 1); E_res sdpvar(T, 1); Curt_w sdpvar(T, 1); Curt_pv sdpvar(T, 1); % 约束集合 C []; % 功率平衡 C [C, Pw - Curt_w Ppv - Curt_pv P_dis P_gen P_grid ... L P_ch P_pump]; % 电池SOC递推 C [C, SOC(1) SOC_init (eta_ch * P_ch(1) - P_dis(1)/eta_dis) * dt / E_bat]; C [C, SOC(2:end) SOC(1:end-1) ... (eta_ch * P_ch(2:end) - P_dis(2:end)/eta_dis) * dt / E_bat]; % 电池容量与功率限制 C [C, SOC_min SOC SOC_max]; C [C, 0 P_ch P_ch_max * u_ch]; C [C, 0 P_dis P_dis_max * u_dis]; C [C, u_ch u_dis 1]; % 抽蓄等效电量递推 C [C, E_res(1) E_res_init (eta_pump * P_pump(1) - P_gen(1)/eta_gen) * dt]; C [C, E_res(2:end) E_res(1:end-1) ... (eta_pump * P_pump(2:end) - P_gen(2:end)/eta_gen) * dt]; % 抽蓄限制 C [C, E_res_min E_res E_res_max]; C [C, 0 P_pump P_pump_max * u_pump]; C [C, 0 P_gen P_gen_max * u_gen]; C [C, u_pump u_gen 1]; % 联络线功率限制 C [C, 0 P_grid P_grid_max]; % 目标函数 Objective sum(C_buy .* P_grid) * dt ... C_bat * sum(P_ch P_dis) * dt ... lambda_w * sum(Curt_w) * dt ... lambda_pv * sum(Curt_pv) * dt ... lambda_pump * sum(P_pump P_gen) * dt; % 求解 options sdpsettings(solver, gurobi, verbose, 1); diagnostics optimize(C, Objective, options);这段代码基本能直接跑前提是装了YALMIP和Gurobi。如果没有外部求解器把solver改成sedumi或者直接用intlinprog版本。几个细节SOC递推式右边用了SOC(1:end-1)向量这是YALMIP支持的向量化写法比for循环快得多。sdpvar的向量维度必须与T一致否则维度报错。如果YALMIP版本较旧对向量等式支持不够好可能得老老实实写for循环但绝大多数现代版本没问题。4.3 无YALMIP的intlinprog实现思路有些读者暂时装不了YALMIP或者公司电脑对安装工具有权限限制那就用Matlab自带的intlinprog。思路是把全部决策变量堆成一个大列向量x把约束写成Ax ≤ b和Aeqx beq的矩阵形式。好处是不依赖额外工具箱缺点是写矩阵非常痛苦特别是SOC递推这种跨时段耦合约束需要构造稀疏矩阵一不小心维度就错了。我的建议是先用YALMIP把模型跑通把结果保存下来作为基准再用intlinprog复现一遍对比两者的目标函数值和决策变量曲线是否一致。这样即使一直用YALMIP交付也能在论文里写为保证可复现性使用Matlab intlinprog与YALMIP交叉验证。4.4 结果可视化调度结果出来之后至少画三张图功率平衡堆叠图横轴时间纵轴功率风电、光伏、电池放电、抽蓄发电、购电从上往下堆叠负荷曲线用黑色实线叠加能直观看出各类电源在每个时段的出力贡献。储能状态曲线电池SOC和抽蓄等效电量的时序变化放同一张图的双Y轴能看出两类储能的时间尺度配合。弃风弃光面积图单独画弃风弃光功率直观体现调度策略在消纳率上的效果。figure; bar(t_h, [Pw_actual, Ppv_actual, P_dis, P_gen, P_grid], stacked); hold on; plot(t_h, L, k-, LineWidth, 2); xlabel(时间/h); ylabel(功率/MW); legend(风电实际出力, 光伏实际出力, 电池放电, 抽蓄发电, 购电, 负荷);这套可视化代码看着简单但我实际研究里吃过亏堆叠图如果不把负荷画成明显的线而是也做成柱状图整个图挤在一起根本看不清。所以负荷曲线单独用plot叠加并且图例顺序和堆叠顺序保持一致别让读者猜颜色。5. 实测中的常见问题与排查记录5.1 求解器报Infeasible怎么办这是做调度优化的人见面第一问。我的排查顺序看输出日志里哪条约束被标红YALMIP的diagnostics里有约束编号信息可以定位矛盾点。检查功率平衡等式是不是负荷曲线和电源出力曲线完全不匹配比如深夜光伏为0、风电也极小、储能全放完、购电又设了上限。检查SOC初值和末值约束如果SOC初值设了0.5但E_bat过小负荷高峰又特别高储能很快放空之后怎么都平衡不了。检查整数变量和连续变量的耦合比如P_ch≤P_ch_max*u_ch这条约束如果P_ch_max设成了infMILP会有数值病态问题求解器直接给不可行或NaN。打过几轮之后我习惯把约束用硬约束/软约束分类功率平衡这种物理绝对定律必须硬SOC上下限这种工程运行限制可以加一个小松弛量比如把SOC_min从0.2放宽到0.19很多时候能快速定位问题是不是出在边界上。5.2 向量化递推的维度错误用sdpvar构造时序变量时最容易出现维度报错多半原因是变量是行向量还是列向量不一致。我习惯全局统一用列向量T×1。负荷、预测、价格数据读进来后马上转一下data data(:)。这行代码能省掉70%的维度报错。另外SOC(2:end) SOC(1:end-1) ...这种写法右边加法的每一项都要确认方向一致不然MATLAB自动广播扩张出来的维度会吓人一跳。5.3 抽蓄模型里的电量和水量混用问题抽蓄建模时如果一会儿用等效电量MWh一会儿用水量m³单位转换很容易混乱。我一开始用水量结果效率和库容的单位老对不上数值量纲乱七八糟。后来统一改成等效电量即把上水库蓄水的重力势能等效成MWh所有约束和状态递推都在电能维度上完成。这样电池和抽蓄都变成同一个储能罐模型主程序处理起来非常统一。代价是如果研究侧重点在水力特性比如水头变化对效率的影响等效电量模型就比较粗糙那只能再回去用水量模型这是项目要求权衡的结果。5.4 求解时间突然变长有一次我把调度时段从24小时扩到168小时intlinprog跑了20多分钟没结束。原因是整数变量数量翻了7倍分支定界树的规模指数增长。解决办法先松弛二元变量改成0-1连续变量跑一遍LP看看目标函数下界再限制求解器的节点数或时间上限最后实在不行把模型拆成分日滚动调度每24小时一个子问题带状态传递也能得到近似最优。研究初期强烈建议先用24小时场景验证逻辑确定没问题再扩时段。6. 结果分析与个人体会6.1 算例结果能说明什么我用一组典型日数据跑了一遍风电装机30MW、光伏20MW、电池20MWh/5MW、废弃矿井抽蓄等效库容60MWh/5MW负荷峰值35MW联络线购电上限25MW。结果显示加入储能后系统购电成本下降约12%弃风弃光率从没有储能时的19.6%降到4.8%。其中抽蓄承担了大约70%的跨时段能量搬移电池主要做平滑和快速调节。这两类储能的分工从SOC和E_res曲线上一目了然电池的SOC一天内频繁起伏说明它一直在追功率波动抽蓄的等效电量则是夜间抽水蓄能、白天放水发电的大日循环节奏。这验证了我最初电池跑短时、抽蓄跑长时的设计思路。6.2 对废弃矿井抽蓄的认知变化做这个课题之前我也觉得抽蓄就是大型水电工程的代名词动辄百万千瓦级别。但把废弃矿井小型抽蓄放进调度模型后看法变了它真正适合的场景是分布式源荷储一体化比如矿区自备电网、工业园区的绿色供能系统。这种场景里负荷不大、源侧也不大大型抽蓄根本吃不下一个5MW的小型矿井抽蓄反而能匹配得刚刚好。而且矿井巷道本身是现成的地下空间改造成本远低于开挖新洞室。当然也必须承认这类项目目前的工程案例还不多技术经济性、地质安全性都还需要更多示范项目验证。但从调度模型层面看它和常规抽蓄在数学形式上几乎无缝衔接意味着未来如果某个矿区的实际数据出来可以直接套用这套调度框架不需要推翻重来。6.3 最后分享一个后来一直沿用的习惯所有调度模型的代码我都尽量写成参数-数据-模型-求解-可视化五段式参数和模型分离。以前我把参数写在脚本开头每次改参数都要全局搜索哪些地方引用了它后面改成单独的参数结构体比如params.bat.E 20所有函数通过params传参逻辑清晰很多。这个习惯让我后来换数据、做敏感性分析时省了几天的返工时间。做研究尤其是多方案对比这个参数集中管理的习惯建议早一点养成。