基于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_{t1}^{T}\sum_{g1}^{G}\left[ C_g(P_{g,t})u_{g,t} S_g^{up}y_{g,t} S_g^{down}z_{g,t}\right] \sum_{t1}^{T}c_w P_t^{curt} \sum_{t1}^{T}c_{ess}(P_t^{ch}P_t^{dis}) ]其中 (C_g(P)) 是煤耗成本函数通常用二次曲线 (aP^2bPc) 近似(u) 是机组运行状态(y) 和 (z) 分别是启动和停机标记变量(P_t^{curt}) 是弃风功率(P_t^{ch})、(P_t^{dis}) 是储能充放电功率。储能运行成本系数 (c_{ess}) 我一般设得很小主要作用不是真实计费而是防止储能一天到晚来回充放电把多余损耗压掉。因为目标函数里带了二次煤耗成本很多人第一反应是用 fmincon。但我更建议先把 (C_g(P)) 分段线性化然后转成 MILP 用 intlinprog 求解。理由有三点第一机组组合本身含0-1整数变量二次目标会让问题变成MIQPintlinprog不直接支持用遗传算法工具箱处理这种离散优化又不如分支定界稳定第二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_{it}^{tm_{up}-1} u_{g,i} \ge m_{up}\cdot y_{g,t} ][ \sum_{it}^{tm_{down}-1} (1-u_{g,i}) \ge m_{down}\cdot z_{g,t} ]在Matlab里构建这个约束时要小心下标越界(tm_{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个调度时段做演示参数G1G2G3额定容量 PmaxMW600400250常规最小出力 PminMW250180100深调最小出力 Pmin_deepMW1409050爬坡速率MW/h12010080煤耗二次系数 a0.00110.00150.0020煤耗一次系数 b14.316.819.5煤耗常数 c604530启动成本元1200900600最小运行时间h432最小停机时间h432储能参数我设置成额定功率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}xb_{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*T1:2*G*T, T, G); v.xd reshape(2*G*T1:3*G*T, T, G); v.xn reshape(3*G*T1:4*G*T, T, G); v.y reshape(4*G*T1:5*G*T, T, G); v.z reshape(5*G*T1:6*G*T, T, G); base 6*G*T; v.pw base (1:T); v.pc base (T1:2*T); v.pd base (2*T1:3*T); v.s base (3*T1: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只放二进制变量索引即可连续变量不要混进去整型变量数量越少求解越快。约束矩阵我习惯逐步拼接。第一步是状态互斥约束 (xdxn\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最稳如果只是想快速看趋势可以先把机组最后的启停方案通过启发式定下来再用二次规划算出力分配。最终调用intlinprogintcon []; 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储能。三个场景共用同一套负荷曲线和风电预测数据只改相应的参数开关。我跑完典型数据后结果趋势大概是这样场景总运行成本万元弃风电量MWh火电深调电量MWh最大深调量MWA 常规无储能321.66800—B 深度调峰无储能306.82151180220C 深调储能294.560760180具体数值会因负荷曲线和风电序列不同而变化但规律基本一致。从场景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}) 在等式左边作为电源项连乘效率以后更容易写错。然后检查状态变量的组合约束。如果 (xdxn\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变量的逻辑关系都会直接影响求解质量。把建模基本功打扎实后面无论换场景、换参数、换更复杂的约束都会从容很多。