基于MILP的风储深度调峰模型：Matlab实现与优化调度全流程

做电力系统优化调度的项目，这几年碰得最多的一类需求，就是把风电、储能和火电深度调峰放到同一个模型里算。风电出力一高，火电又不能随手就停，火电机组最低技术出力卡在那里，电网低谷时段很容易出现弃风。把储能加进去，等于多了一块能快速响应的缓冲，但储能容量终归有限，真正要兜底的还是得靠火电往下压。用Matlab实现风储深度调峰模型，难点往往不在“会不会调用优化函数”，而在于怎么把机组启停、深调区间、爬坡、储能SOC这些约束组织成一个能稳定求解的优化问题。这篇文章从建模思路开始，把混合整数线性规划（MILP）建模、Matlab代码实现、算例结果和调试经验完整地捋一遍，适合正在写调度类论文、做Matlab仿真，或者刚从普通机组组合转向深度调峰建模的读者。

1. 深度调峰模型该建什么：先把目标与约束理清楚

1.1 风电高占比下，调峰为什么“深”不起来

先明确一个基础概念：什么叫深度调峰。常规火电机组有最小技术出力，一般设计在额定容量的40%到50%左右，再往下锅炉燃烧稳定性就会出问题。深度调峰就是通过设备改造、低负荷稳燃等技术手段，把最小出力压到30%、甚至20%这个量级。代价是煤耗上升、设备损耗增加，所以调度里通常会把深度调峰当作一种特殊运行状态来对待，而不是简单地把下限往下一改。

这个动作放到优化模型里，立刻带来一个数学问题：机组不再是一个简单的0-1二值状态（停机或运行），运行状态还得分“常规区间”和“深度调峰区间”，而且两个区间之间有一个不连续的断层。比如一台600MW机组，额定区间可能是250~600MW，深调区间是140~250MW，那么140到250之间被切成两段，150MW这个出力点理论上合法，但要清楚它落在哪个区间。这个非凸特性，正是深度调峰模型比普通机组组合麻烦的地方。

风电出力高的时段通常集中在凌晨和午间。这个时候负荷偏低，风电又顶着预测值出力，火电全部压到常规下限仍然供大于求，只能让一部分机组继续下探，也就是深度调峰，同时储能充电把多余电量转移。如果这些手段都用完了还是压不下来，最终选择只能是弃风。所以风储深度调峰模型的本质，是在“火电深调、储能充电、弃风”三个动作之间做经济性权衡，目标是让总运行成本最低，而不是机械地追求某单一指标最小。

1.2 目标函数：从煤耗成本到弃风惩罚

我建议把目标函数设为综合运行成本最小，而不是简单地把弃风率最小或者煤耗最小作为唯一目标。原因很直接：如果只追求弃风率最小，模型会把火电压到极限甚至频繁启停，结果是风光消纳率很好看，实际运行成本却高得离谱；如果只追求煤耗最小，模型会倾向少开机，旋转备用约束又容易出问题。更符合工程习惯的目标函数应该是燃煤成本、机组启动停运成本、弃风惩罚成本和储能运行成本一起最小，权重通过价格系数来控制。

写成数学表达式就是这样：

[ \min \sum_{t=1}^{T}\sum_{g=1}^{G}\left[ C_g(P_{g,t})u_{g,t} + S_g^{up}y_{g,t} + S_g^{down}z_{g,t}\right] + \sum_{t=1}^{T}c_w P_t^{curt} + \sum_{t=1}^{T}c_{ess}(P_t^{ch}+P_t^{dis}) ]

其中 (C_g(P)) 是煤耗成本函数，通常用二次曲线 (aP^2+bP+c) 近似；(u) 是机组运行状态，(y) 和 (z) 分别是启动和停机标记变量；(P_t^{curt}) 是弃风功率，(P_t^{ch})、(P_t^{dis}) 是储能充放电功率。储能运行成本系数 (c_{ess}) 我一般设得很小，主要作用不是真实计费，而是防止储能一天到晚来回充放电，把多余损耗压掉。

因为目标函数里带了二次煤耗成本，很多人第一反应是用 fmincon。但我更建议先把 (C_g(P)) 分段线性化，然后转成 MILP 用 intlinprog 求解。理由有三点：第一，机组组合本身含0-1整数变量，二次目标会让问题变成MIQP，intlinprog不直接支持，用遗传算法工具箱处理这种离散优化又不如分支定界稳定；第二，MILP有全局收敛性质，不像非线性求解器那样容易陷入局部最优；第三，求解速度快，论文里常见的十几台机组、24到96时段算例，MILP求解器非常成熟。分段线性化的具体做法是在 (P_{min}) 到 (P_{max}) 之间取若干个分段点，比如取五个点，用折线去拟合二次成本曲线。

这里要特别提醒一个细节：很多教材写分段线性成本时直接用分段函数表达，但在MILP里要避免“多个分段函数按区间激活”造成的非凸问题。对于凸的成本曲线，可以写成叠加形式，把总出力表达为最小出力加上各段出力增量，并要求这些增量按顺序累加，这样不需要引入额外的0-1变量。如果忽略这一点，求解结果会出现成本曲线跳变，或者干脆得不出合理解。

1.3 约束条件：机组、风电与储能的核心限制

先看电力平衡。这个约束没什么悬念，每个时段火电总出力、风电实际出力和储能放电功率之和等于负荷加上充电功率：

[ \sum_g P_{g,t} + P_{w,t}^{use} + P_t^{dis} = L_t + P_t^{ch} ]

风电实际出力 (P_{w,t}^{use}) 和弃风量之间满足 (P_{w,t}^{use} = P_t^{forecast} - P_t^{curt})。

机组运行区间是最容易写错的地方。假设常规最小出力是 (P_{min})，深调能压到 (P_{min}^{deep})，额定出力是 (P_{max})。机组状态分三档：停机（0）、深调区间（(P_{min}^{deep}\sim P_{min})）、常规区间（(P_{min}\sim P_{max})）。为了用MILP表达这个非凸可行域，我习惯给每台机组引入两个0-1变量：(x_g^{deep}) 表示是否运行在深调区间，(x_g^{norm}) 表示是否运行在常规区间，再通过 (x_g^{deep}+x_g^{norm} \le 1) 把状态隔离开。出力约束写成：

[ P_{min}^{deep}\cdot x_g^{deep} + P_{min}\cdot x_g^{norm} \le P_{g,t} \le P_{min}\cdot x_g^{deep} + P_{max}\cdot x_g^{norm} ]

仔细看这个式子：如果机组在深调区间，(x^{deep}=1, x^{norm}=0)，出力就被限制在深调下限和常规下限之间；如果机组在常规区间，就限制在常规下限和额定出力之间；如果两个状态变量都为0，则出力被强制为0，也就是停机。这组约束是线性的，非常干净。如果论文里还要分30%、40%、50%多个深调档位，思路完全一样，多准备几个0-1变量即可。

爬坡约束在MILP里比较麻烦，因为上一时段机组可能是停机状态，出力跳变和启停过程是耦合的。常见写法是引入启动标记 (y_{g,t}) 和停机标记 (z_{g,t})：

[ P_{g,t}-P_{g,t-1} \le R^{up}\cdot u_{g,t-1} + P_{max}\cdot y_{g,t} ]

[ P_{g,t-1}-P_{g,t} \le R^{down}\cdot u_{g,t} + P_{max}\cdot z_{g,t} ]

这样允许机组在启动和停机时段出现较大的出力跳变。很多初学者漏掉右侧的第二项，结果机组启动那一小时爬坡约束怎么都满足不了，模型报不可行。

最小启停时间约束同样不能省。它要求机组启动后至少维持运行 (m_{up}) 个小时，停机后至少维持停运 (m_{down}) 个小时：

[ \sum_{i=t}^{t+m_{up}-1} u_{g,i} \ge m_{up}\cdot y_{g,t} ]

[ \sum_{i=t}^{t+m_{down}-1} (1-u_{g,i}) \ge m_{down}\cdot z_{g,t} ]

在Matlab里构建这个约束时要小心下标越界，(t+m_{up}-1) 可能超过总时段数T，需要通过截断处理。

旋转备用约束也不能漏，否则模型会把所有机组压到最低而不考虑系统可靠性。备用需求可以写成每个时段在线机组可调容量上限之和减去当前出力，大于等于系统备用需求。储能部分用SOC动态方程串联：

[ SOC_t = SOC_{t-1} + \left(\eta_{ch}P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}}\right)\cdot \frac{\Delta t}{E} ]

再加上充放电功率上限、SOC上下限和初末值等式。对于“不能同时充电和放电”这个约束，工程上不一定非要用0-1变量限制，因为同时充放电会带来能量损耗和多余成本，模型在理性情况下不会主动这么干，我一般只给很小的 (c_{ess}) 系数，让这个逻辑自然成立。

2. Matlab代码实现：优化问题怎么落到矩阵上

2.1 数据组织与参数初始化

我习惯把数据集中到一个脚本文件里，用结构体组织，不要散落着一堆全局变量。原因很简单：后面做参数敏感性分析时，改一个容量参数不用翻三个脚本。下面给出一套三台火电机组的典型参数，用24个调度时段做演示：

储能参数我设置成额定功率50MW、容量200MWh、充放电效率0.95、SOC范围0.1到0.9，初始SOC为0.5，同时要求调度周期结束时回到0.5。负荷曲线按典型的峰谷形态设置，峰谷比大概1.7；风电预测数据取一个日出力高、凌晨也偏高的典型日序列，模拟弃风高风险场景。

在Matlab里用结构体保存这些数据，代码清晰很多：

this.unit.Pmax = [600; 400; 250]; this.unit.Pmin = [250; 180; 100]; this.unit.Pmin_deep = [140; 90; 50]; this.unit.ramp = [120; 100; 80]; this.unit.a = [0.0011; 0.0015; 0.0020]; this.unit.b = [14.3; 16.8; 19.5]; this.unit.c = [60; 45; 30]; this.unit.S_up = [1200; 900; 600]; this.unit.m_up = [4; 3; 2]; this.unit.m_down = [4; 3; 2]; this.ess.Pmax_ch = 50; this.ess.Pmax_dis = 50; this.ess.E = 200; this.ess.eta_ch = 0.95; this.ess.eta_dis = 0.95; this.ess.soc_min = 0.1; this.ess.soc_max = 0.9;

2.2 决策变量编组与约束矩阵搭建

intlinprog的标准输入形式是目标向量f、整数变量索引intcon、不等式约束 (A_{ineq}x\le b_{ineq})、等式约束 (A_{eq}x=b_{eq})，以及变量上下界。变量怎么排，直接决定后面所有约束好不好写。我习惯按“时段×机组”的矩阵顺序展开，然后用一个结构体记录每个子变量的起始索引，这样后期检查错误方便很多。

G = numel(this.unit.Pmax); T = 24; % 变量顺序：p, u, xd, xn, y, z, pw, pc, pd, soc v.p = reshape(1:G*T, T, G); v.u = reshape(G*T+1:2*G*T, T, G); v.xd = reshape(2*G*T+1:3*G*T, T, G); v.xn = reshape(3*G*T+1:4*G*T, T, G); v.y = reshape(4*G*T+1:5*G*T, T, G); v.z = reshape(5*G*T+1:6*G*T, T, G); base = 6*G*T; v.pw = base + (1:T); v.pc = base + (T+1:2*T); v.pd = base + (2*T+1:3*T); v.s = base + (3*T+1:4*T); nvars = base + 4*T; lb = zeros(nvars, 1); ub = 800 * ones(nvars, 1);

第1台机组变量实际上不需要把 (u) 直接设为二进制？需要，u/xd/xn/y/z是0-1变量，P、pw、soc等都是连续变量。intcon只放二进制变量索引即可，连续变量不要混进去，整型变量数量越少求解越快。

约束矩阵我习惯逐步拼接。第一步是状态互斥约束 (xd+xn\le 1)：

Aineq = []; bineq = []; temp = zeros(G*T, nvars); b_temp = ones(G*T, 1); row = 0; for t = 1:T for g = 1:G row = row + 1; temp(row, v.xd(t,g)) = 1; temp(row, v.xn(t,g)) = 1; end end Aineq = [Aineq; temp]; bineq = [bineq; b_temp];

第二步是机组出力区间约束。这一步把每台机组每个时段的上下界写进约束矩阵，上下界用 (xd)、(xn) 的线性表达式表示：

temp = zeros(2*G*T, nvars); b_temp = zeros(2*G*T, 1); row = 0; for t = 1:T for g = 1:G row = row + 1; % 下限 P >= Pmin_deep*xd + Pmin*xn temp(row, v.p(t,g)) = 1; temp(row, v.xd(t,g)) = -this.unit.Pmin_deep(g); temp(row, v.xn(t,g)) = -this.unit.Pmin(g); row = row + 1; % 上限 P <= Pmin*xd + Pmax*xn temp(row, v.p(t,g)) = 1; temp(row, v.xd(t,g)) = -this.unit.Pmin(g); temp(row, v.xn(t,g)) = -this.unit.Pmax(g); b_temp(row) = 0; end end Aineq = [Aineq; temp]; bineq = [bineq; b_temp];

注意我写的全是“小于等于”形式，所以下限约束写成 (-P + Pmin_deep\cdot xd + Pmin\cdot xn \le 0)，这样符号不容易乱。

第三步把启动停机标记的绑定关系写进去：启动标记必须大于等于 (u_t - u_{t-1})，停机标记必须大于等于 (u_{t-1} - u_t)：

temp = zeros(3*G*T, nvars); b_temp = zeros(3*G*T, 1); row = 0; for t = 1:T for g = 1:G % y >= u_t - u_{t-1} row = row + 1; temp(row, v.y(t,g)) = 1; temp(row, v.u(t,g)) = -1; if t > 1 temp(row, v.u(t-1,g)) = 1; end % z >= u_{t-1} - u_t row = row + 1; temp(row, v.z(t,g)) = 1; temp(row, v.u(t,g)) = 1; if t > 1 temp(row, v.u(t-1,g)) = -1; end % u_t <= 1 row = row + 1; temp(row, v.u(t,g)) = 1; b_temp(row) = 1; end end Aineq = [Aineq; temp]; bineq = [bineq; b_temp];

最小启停时间约束看在Matlab里要循环 (t) 和 (g)，把求和区间逐项填进矩阵。这部分代码不复杂，只是容易越界，写的时候统一判断 (t + m - 1 \le T)，超过的时段取到T就截断。

等式约束有两类，第一类是功率平衡，第二类是储能SOC递推。功率平衡每个时段一行：

Aeq = zeros(T, nvars); beq = zeros(T, 1); for t = 1:T for g = 1:G Aeq(t, v.p(t,g)) = 1; end Aeq(t, v.pw(t)) = 1; Aeq(t, v.pd(t)) = 1; Aeq(t, v.pc(t)) = -1; beq(t) = load_curve(t) - wind_forecast(t); end

我这里直接把风电预测移到等式右边，再用 (pw) 表示实际消纳的风功率。如果想显式保留弃风变量，可以写成 (pw + pcurt = forecast)，两种方式等价，看个人习惯。

SOC递推方程写成等式形式：

Aeq_soc = zeros(T, nvars); beq_soc = zeros(T, 1); for t = 1:T Aeq_soc(t, v.s(t)) = 1; if t > 1 Aeq_soc(t, v.s(t-1)) = -1; else beq_soc(t) = 0.5; % SOC(0)=0.5 end Aeq_soc(t, v.pc(t)) = -this.ess.eta_ch * dt / this.ess.E; Aeq_soc(t, v.pd(t)) = dt / (this.ess.eta_dis * this.ess.E); end Aeq = [Aeq; Aeq_soc]; beq = [beq; beq_soc];

2.3 目标函数向量化和求解调用

目标函数里涉及分段线性的煤耗成本，我在前面说过用叠加方式处理。为了简化演示，这里给出一个常用的做法：把每台机组出力相对最小出力的增量部分拆成几个连续段，每个分段用连续变量表示，分段系数是斜率差。实现时需要额外引入分段长度变量，代码量多一些，但逻辑不复杂。

另一种更省事的方式是直接用二次成本函数配合 fmincon 做非线性规划，前提是你把机组状态变量固定以后再做连续优化，否则混合整数非线性问题用内置工具很难收拾。我自己的经验是：论文复现用MILP最稳，如果只是想快速看趋势，可以先把机组最后的启停方案通过启发式定下来，再用二次规划算出力分配。

最终调用intlinprog：

intcon = []; for t = 1:T for g = 1:G intcon = [intcon, v.u(t,g), v.xd(t,g), v.xn(t,g), v.y(t,g), v.z(t,g)]; end end options = optimoptions('intlinprog', ... 'MaxTime', 600, ... 'Display', 'final', ... 'CutGeneration', 'basic'); [xopt, fval, exitflag] = intlinprog(f, intcon, Aineq, bineq, Aeq, beq, lb, ub, options);

f向量需要逐一填：机组煤耗线性段斜率、启动停机成本、弃风惩罚、储能成本。intlinprog返回后，把xopt按之前的结构体还原成果矩阵：

P = reshape(xopt(v.p), T, G); U = reshape(xopt(v.u), T, G); XD = reshape(xopt(v.xd), T, G); XN = reshape(xopt(v.xn), T, G); SOC = xopt(v.s);

3. 仿真算例与结果分析：模型到底算出了什么

3.1 算例场景设计

把模型跑起来之前，我先设计了三个场景做对比。场景A：火电只能按常规最小出力运行，不开深调，也不加储能，作为基准。场景B：火电允许进入深调区间，储能不参与。场景C：允许深调，同时投入200MWh储能。三个场景共用同一套负荷曲线和风电预测数据，只改相应的参数开关。

我跑完典型数据后，结果趋势大概是这样：

具体数值会因负荷曲线和风电序列不同而变化，但规律基本一致。从场景A到场景B，弃风量大幅下降、总成本下降，说明火电有了深调选项以后，低谷时段大部分弃风压力被机组自身消化了。再从场景B到场景C，弃风进一步被压到60MWh，同时深调电量反而减少，说明储能替代了一部分深度调峰功能，火电不需要一直压在极限低位，系统运行变得更加灵活。

3.2 结果出来以后先别急着分析

拿到结果以后我通常会做两件“笨”事情。第一，把每个时段的电力平衡残差打印出来，也就是用求解出来的变量反算 ( \sum P_g + P_w + P_dis - L - P_ch )，如果某个时段残差很大，说明约束矩阵的符号或者索引有问题，这个检查能拦住一大批低级错误。第二，检查火电出力轨迹有没有真的避开不可运行区间。比如G1的常规最小出力是250MW，深调下限是140MW，那么G1在深调状态下出力应该落在140到250之间，如果看到某时段G1出力跑到120MW，肯定是状态变量组合写错了。

绘图是最直观的检查方式。我习惯把负荷、火电总出力、风电实际出力、储能SOC画在同一张图上，用堆积面积图显示火电和风电的分担比例。如果看到SOC曲线在某个时段出现不合理的锯齿状突变，多半是充放电效率写反了，或者SOC递推方程里dt和容量E的单位没有统一。

3.3 储能容量变化对调峰效果的影响

对储能容量做敏感性分析是这类模型里最有实用价值的部分。我把储能能量从100MWh逐步增加到400MWh，步长50MWh，每跑一次记录弃风率和总成本。典型结果是：弃风率下降的斜率前段很陡、后段明显平缓。原因不复杂，低谷时段能被转移的负荷总量有限，储能又受50MW功率上限约束，能量再大也只有那几个小时能充电，边际效益自然递减。

更值得关注的是储能和深度调峰之间的替代关系。在小容量储能下，系统仍然需要火电深度调峰来压低谷；当储能容量增大以后，深调电量明显下降，火电平均出力抬高，整体煤耗成本也在降低。但从“减少1MWh弃风所需增加的储能投资”角度看，容量超过某个阈值以后，再加大电池容量就不划算了。这个阈值正是这套模型能帮你算出来的东西——把火电深调改造和储能投资放在同一个成本函数里做比较，而不是拍脑袋定容量。

类似地，把G1的深调下限从140MW往下降到110MW，也能看到弃风下降、成本下降的边际递减效应。深调越深，煤耗二次项上升得越快，设备损耗越大，最终谁更经济，取决于深调改造费用、补偿电价和储能成本。这套模型的价值就是把两边放到同一个天平上比。

4. 常见问题、调试经验与避坑指南

4.1 模型不可行：先按这个顺序排查

intlinprog返回exitflag非1，提示找不到可行解，这是最常见的失败模式。我的排查顺序很固定。

先检查变量上下界。看是不是某个二进制变量的ub设成了0，或者风电、储能的lb、ub把可行域裁掉了一大块。接着检查等式约束：功率平衡方向有没有写反，储能SOC初末值是否一致。最常见的问题是充电功率符号混乱，(P_{ch}) 应该在等式右边作为负荷项，(P_{dis}) 在等式左边作为电源项，连乘效率以后更容易写错。

然后检查状态变量的组合约束。如果 (xd+xn\le 1) 被误写成等于1，那么机组要么深调要么常规，却永远不能停机，系统低谷时段备用容量会凭空少一大块，很可能直接不可行。

调试期的小技巧是把约束逐批去掉，逐步缩小范围。先不启用最小启停时间和爬坡，只保留电力平衡和出力区间，看模型能不能出解。能出解就把约束一组一组加回来，加哪组挂了，问题就锁定在哪组。SOC的末值等式在调试期可以先放宽为区间 ([0.4, 0.6])，跑通以后再收紧。

4.2 求解太慢：问题多半出在变量设计

MILP的求解速度和0-1变量数量强相关。第一要务是检查intcon里有没有混进连续变量，比如把出力P、SOC也设置成整数，那规模会直接爆炸，24时段三台机组可能就要跑几分钟。第二，机组状态变量其实可以用更紧凑的方式表达，例如只引入“深调”这一个0-1变量，常规区间用另一个二进制变量的补逻辑来表示，变量数量能省下一些。

如果模型规模确实很大，比如上百台机组、几百个时段，建议给intlinprog设置合理的停止条件，限制MaxTime或者MaxNodes。实际工程中并不需要每次都求到全局最优，一个gap在0.5%以内的次优解已经完全够用。另外，intlinprog支持通过start point传入初始可行解，用上一轮的调度结果作为热启动，能显著减少分支定界的搜索时间。

代码层面也不要忽略向量化。我在2.2节演示的循环版本对初学者很友好，但同样的约束用稀疏矩阵和矩阵运算一次构建，求解前的预处理时间能差一个数量级。第一版用循环保证正确性，性能瓶颈出现了再向量化，这是写优化模型的通用策略。

4.3 结果合理但“不自然”：从指标反推模型缺陷

如果模型顺利求出解，但结果看起来不符合实际，问题往往出在目标函数权重或者遗漏约束上。

机组频繁启停，一小时内开机、一小时后又停机，大概率是最小启停时间约束没写，或者启动成本给得太低。给 (S^{up}) 一个足够大的值，让模型把“多开一台机组连续运行”和“频繁启停”之间的成本权衡关系表达清楚。

储能一整天几乎不动，先问弃风惩罚系数 (c_w) 是不是给得太低。如果弃风惩罚低于储能运行成本，模型当然会直接选择弃风，储能的边际价值体现不出来。反之，如果储能充满以后还在充电，说明SOC上限约束或者充电功率约束没生效。

风功率曲线出现明显的弃风但火电还维持在高出力，说明风电出力上限约束写错了。检查 (P_w^{use}) 是否被ub限制，以及弃风变量的非负性。还有一种是风功率预测值本身就设置得过低，模型不知道有多余的风可以消纳，这一点在做数据时要单独核验。

4.4 我的几个实操小习惯

最后分享几条自己摸索出来的经验。

第一，模型参数用有名值（MW、MWh、元）而不是标幺值，虽然标幺值在电力系统里很常见，但优化模型的量纲太乱时，约束矩阵的病态程度会上升。用有名值写清楚，intlinprog的默认数值处理反而更稳。

第二，风电出力和负荷曲线先做平滑处理。如果原始数据里有个别跳变点，会导致某些时段负荷突变、爬坡约束和备用约束同时吃紧，模型被迫启停机组，结果看起来会很折腾。

第三，建议先用2小时的小规模算例验证约束矩阵。小规模时所有中间变量都能直接打印出来，哪里索引错了马上能定位。我一般会构造一个已知最优解的小场景，比如只有一台机组和一台储能，手动算一遍再和求解器对照，确认一致后才跑完整24小时算例。

这套风储深度调峰模型的工程实现，核心其实不在“深度学习”或者“智能化算法”，而在把复杂的物理约束和运行逻辑准确翻译成MILP矩阵。每一个索引、每一个符号、每一条0-1变量的逻辑关系，都会直接影响求解质量。把建模基本功打扎实，后面无论换场景、换参数、换更复杂的约束，都会从容很多。