风电、光伏大规模并网之后电网的净负荷曲线变得越来越陡。午间光伏大发时全网负荷被压出一个大坑傍晚光伏退坡后负荷又急剧拉升火电机组跟着上蹿下跳AGC指令根本来不及响应。传统做法是靠火电硬扛调峰但火电的爬坡速率和最小技术出力摆在那里深度调峰不仅煤耗飙升机组寿命也受影响。这种背景下把具备快速调节能力的水电、储能拉进调峰主战场让风光水火储五个电源按照各自的物理特性互补协调就成了调度侧最现实的一条路。这篇文章就是从Matlab代码落地的角度把一个“计及调峰主动性的多能互补协调优化调度”问题从建模到求解完整过一遍。适合刚接触电力系统优化调度的研究生、做新能源并网研究的工程师以及对YALMIP和Cplex求解框架感兴趣的同行参考。1. 课题核心拆解什么是“调峰主动性”驱动的互补调度1.1 先从多能系统的耦合关系说起做多能互补调度之前必须先摸清每个电源的“脾气”。风电是反调峰特性夜间出力大、白天反而小而且预测误差随超前时间快速增大基本不可控光伏相对规律一些白天大发、中午到峰值、夜晚零出力但碰上云层遮挡可能出现短时间内出力滑坡。水电调节能力强、启停快是系统里最灵活的常规电源但它受来水约束和日调节库容限制一天能用的水量就那么多不能无节制地调。火电出力可调范围其实不小问题是爬坡速率慢最小技术出力一般要到额定出力的40%到50%还得忍受频繁调节带来的煤耗飙升和设备损耗。储能响应速度最快毫秒级就能从满充切到满放可容量有限充放电循环还有寿命成本。把五种电源放在一起看互补的逻辑其实很清晰风光提供电量主体但波动性大需要有人给它们“打圆场”水电机动灵活适合做快速调节的先锋火电守住基底负荷保证系统不缺电储能作为最后一道缓冲负责短时峰谷搬移和应急备用。但这只是静态描述。真正难的是在时序上安排“谁先动、谁后动、动多少”这就进入了优化调度的核心问题。1.2 调峰主动性的建模动机传统调度模式下调峰是被动的——调度中心下达指令火电机组跟着净负荷曲线走。净负荷一陡火电就得拼命爬坡来不及就切部分风电光伏这就是我们常说的弃风弃光。被动调峰的问题很具体一方面火电爬坡速率跟不上净负荷变化率尤其是傍晚光伏退坡和晚高峰负荷叠加的那一两个小时另一方面深度调峰时火电为了稳燃得投油成本非线性上涨频繁变负荷还会缩短检修周期。“调峰主动性”这个概念就是针对这个痛点提出来的。它不是说火电不调峰了而是把调峰压力按照调节成本和经济性在多个资源之间合理分摊调度模型在目标函数里显式加入调峰代价和净负荷平滑度的惩罚项让模型自己去决定“谁来调峰、调多少、什么时候调”。水电调节成本低、速度快主动承担主要跟踪任务储能响应快但容量有限负责高峰短时缺口火电输出相对平稳只在净负荷趋势性变化时才调整出力。这样做的好处不仅是弃风弃光率降下来了整个系统的运行压力也小很多。1.3 典型调度架构与时间尺度实际工程中的调度是分层分时的日前计划做24小时的大框架、以1小时为分辨率日内滚动修正用15分钟级到5分钟级的窗口实时调整则靠AGC等自动控制手段。这篇Matlab代码实现的是日前调度层也是优化调度课题里最核心、最常被拿来练手的一块。日前调度的本质是在已知次日负荷预测、风电光伏预测出力上限、来水预测等约束的前提下安排各机组次日24个时段的出力计划和启停计划使系统运行成本最低、新能源消纳最大、调峰效果最好。这个决策窗口足够长适合用混合整数规划来建模求解也方便后续扩展日内滚动和实时环节。2. 优化模型构建目标函数与约束条件的完整数学表达2.1 目标函数的设计成本和调峰效果双目标模型的数学形式决定了代码怎么写所以这一步必须想清楚。我用的目标函数是“运行成本弃能惩罚调峰主动性代价”三项相加转成单目标最小化问题。火电运行成本用二次函数拟合比较常见[ C_{th}\sum_{t1}^{T}\sum_{i1}^{N_g}\left(a_i P_{th,i}(t)^2 b_i P_{th,i}(t) c_i S_{th,i}(t)\right) ]其中 (a_i,b_i,c_i) 是煤耗特性系数(S_{th,i}(t)) 是启停状态二进制变量。这个二次项直接丢给求解器会让模型变成MIQP虽然Cplex能解但速度偏慢。我在工程实现里一般手工做分段线性化5段就够用既保持MILP结构又省求解时间。弃风弃光惩罚项是最常见的“软约束”写法[ C_{cut}\sum_{t1}^{T}\left[\rho_w\left(P_{w}^{pre}(t)-P_w(t)\right)\rho_s\left(P_{s}^{pre}(t)-P_s(t)\right)\right] ]惩罚系数取得很大比如500元/MWh这样模型在可行的情况下不会主动弃电只有确保证功率平衡不可避免时才会牺牲新能源出力。调峰主动性代价是这篇文章的重点。我在目标函数里加了一项净负荷波动惩罚[ C_{peak}\sum_{t1}^{T}\kappa \left(D(t)-D_{avg}\right)^2 ]其中 (D(t)P_{load}(t)-P_w(t)-P_s(t)) 是净负荷(D_{avg}) 是全天净负荷的平均值。这一项看着简单实际作用很直接它迫使调度结果把净负荷曲线“擀平”逼迫水电、储能这些源头调节资源主动出动而不是让火电在那儿硬扛。(\kappa) 的取值决定调峰的权重我一般取系统平均电价的0.1到0.3倍太大容易过度牺牲经济性太小又起不到平滑作用需要试算几轮。2.2 约束条件别漏了时间耦合约束条件是优化模型里最容易“阴沟翻船”的地方。功率平衡是最基本的[ \sum_{i1}^{N_g}P_{th,i}(t)P_h(t)P_w(t)P_s(t)P_d(t)-P_c(t)P_{load}(t) ]火电约束除了出力上下限 (P_{th,i}^{min}S_{th,i}(t)\le P_{th,i}(t)\le P_{th,i}^{max}S_{th,i}(t)) 之外爬坡约束是关键它把相邻时段的状态耦合在一起[ -R_{down,i}\le P_{th,i}(t)-P_{th,i}(t-1)\le R_{up,i} ]另外还有最小启停时间约束用二进制变量配合逻辑矩阵实现这个在后面代码部分会提到。水电的约束相对复杂一点主要是水量平衡和库容限制。日调节水库可以写成[ V(t1)V(t)q_{in}(t)-q_{out}(t) ]出力与发电流量的关系简化为线性(P_h(t)k_h q_{out}(t))(V(t)) 和 (q_{out}(t)) 都要限制在上下限之间。储能方面SOC递推是典型的时序耦合约束[ SOC(t1)SOC(t)\frac{\eta_c P_c(t)\Delta t}{E_{max}}-\frac{P_d(t)\Delta t}{\eta_d E_{max}} ]还要确保不能同时充放电所以引入了辅助二进制变量 (b_{st}(t))用大M法约束(P_c(t)\le M b_{st}(t))(P_d(t)\le M(1-b_{st}(t)))。这里大M取值很有讲究取太大会造成数值病态取太小又会带来约束不可行经验上取储能最大功率的2到5倍比较稳妥。2.3 非线性项的线性化处理心得上面这套模型里最容易出问题的就是非线性项。火电成本二次函数我用分段线性近似把出力区间切成5段每段用线性插值水电出力特性有人会用“水头-出力”二维查表但工程上一天内的水头变化不大简化为线性系数 (k_h) 完全够用。还有绝对值项和min/max操作比如净负荷峰谷差的计算在YALMIP里可以用辅助变量加不等式约束来替代不要直接写abs()或min()否则求解器会抱怨非线性。这些细节在代码实现时都要提前想好不然建模半天最后卡在求解器报错上就尴尬了。3. Matlab代码实现:从数学公式到可运行的MILP模型3.1 代码框架与变量定义Matlab里做优化调度的标准组合是YALMIPCplex。YALMIP负责建模、生成标准形式Cplex负责求解混合整数规划。用这套框架的好处是代码可读性好换求解器也方便——Gurobi、Mosek都能无缝切换。我习惯把整个代码模块化数据准备、变量定义、约束构建、目标函数、求解、结果后处理六个模块各管一段调起来不头疼。变量定义用YALMIP的sdpvar和binvar两条路分开定义连续变量和二进制变量T 24; n_fire 2; % 火电机组台数 % 连续变量火电出力、水电出力、风电出力、光伏出力 P_fire sdpvar(T, n_fire, full); P_hydro sdpvar(T, 1, full); P_wind sdpvar(T, 1, full); P_pv sdpvar(T, 1, full); % 储能变量充电功率、放电功率、SOC状态 P_ch sdpvar(T, 1, full); P_dis sdpvar(T, 1, full); SOC sdpvar(T1, 1, full); % 二进制变量火电启停、储能充放电状态 S_fire binvar(T, n_fire, full); b_st binvar(T, 1, full);这里的SOC定义成T1是为了让最后一个时段的SOC能延续到次日方便做日间耦合。做工程和做学术在这点上有区别学术文章里很多会直接把SOC终值设成初值但对于连续日的调度SOC初始状态应该来自上一日的边界条件不能拍脑袋定。3.2 约束构建与目标函数装配约束用加号拼接在同一个Constraints变量上这是YALMIP的惯用写法。先写功率平衡和各机组基本约束Constraints []; % 功率平衡火电水电风电光伏放电-充电 负荷 for t 1:T Constraints [Constraints, ... sum(P_fire(t,:),2) P_hydro(t) P_wind(t) P_pv(t) P_dis(t) - P_ch(t) P_load(t)]; end % 火电出力上下限状态变量约束 for t 1:T for i 1:n_fire Constraints [Constraints, ... P_fire_min(i) * S_fire(t,i) P_fire(t,i) P_fire_max(i) * S_fire(t,i)]; end end % 火电爬坡约束 for t 2:T for i 1:n_fire Constraints [Constraints, ... -P_ramp_down(i) P_fire(t,i) - P_fire(t-1,i) P_ramp_up(i)]; end end % 储能SOC递推与充放限制 for t 1:T Constraints [Constraints, SOC(t1) SOC(t) eta_c*P_ch(t)/E_max ... - P_dis(t)/(eta_d*E_max)]; Constraints [Constraints, 0 P_ch(t) P_ch_max*b_st(t)]; Constraints [Constraints, 0 P_dis(t) P_dis_max*(1-b_st(t))]; Constraints [Constraints, 0 SOC(t) 1]; end % 水电出力与库容约束 for t 1:T Constraints [Constraints, 0 P_hydro(t) P_hydro_max]; Constraints [Constraints, V(t1) V(t) q_in(t) - P_hydro(t)/k_h]; Constraints [Constraints, V_min V(t1) V_max]; end注意爬坡约束在t1时段没有定义因为端点没有历史值。如果你的机组有初始出力可以把第一时段的爬坡也约束上用初始出力当基准否则默认当天起点状态是自由的。目标函数按照2.1节设计的三项成本装配Objective 0; for t 1:T for i 1:n_fire % 分段线性化的燃料成本函数 cost_pre[i] 是输入向量 Objective Objective cost_fire_interp{i}(P_fire(t,i)) ... start_cost(i) * max(0, S_fire(t,i) - S_fire(max(1,t-1),i)); end Objective Objective rho_w * (P_w_pre(t) - P_wind(t)) ... rho_s * (P_s_pre(t) - P_pv(t)); end % 调峰主动性惩罚项 net_load P_load - P_wind - P_pv; Objective Objective kappa * sum((net_load - mean(net_load)).^2);如果做了分段线性化直接把分段插值函数写进去即可如果直接传二次多项式YALMIP会生成MIQP模型Cplex也能解但大场景下求解时间会肉眼可见地拉长。在这个例子里我可以放心用二次形式因为变量规模小但如果扩展到几十台机组就要慎重了。3.3 求解器配置与参数调试YALMIP调用Cplex的解法非常直接ops sdpsettings(solver, cplex, verbose, 2, ... cplex.mip.tolerances.mipgap, 1e-4, ... cplex.timelimit, 300); sol optimize(Constraints, Objective, ops); if sol.problem 0 disp(求解成功); else disp(sol.info); endmipgap设置成1e-4是我常用的标准既保证精度又不至于死磕最优性。5分钟求解时间的上限对于24时段小规模算例足够遇到大规模算例也可以先设300秒看gap收敛情况再做调整。verbose级别开到2可以输出每次迭代的gap变化方便判断方程压在哪里。早期调试时可以开3甚至4但正常跑批量仿真我一般设2否则控制台输出太多反而淹没了关键信息。4. 典型场景仿真从数据准备到结果分析4.1 系统参数与预测数据为了讲清楚这台代码能干什么我搭了一个简化的测试系统两台火电机组、一座日调节水电站、一个风电场、一个光伏电站、一套储能系统。主参数如下设备参数数值火电1额定/最小出力300MW / 120MW火电2额定/最小出力200MW / 90MW水电最大出力100MW风电预测出力峰值约180MW光伏预测出力峰值约150MW储能功率/容量50MW / 100MWh负荷曲线模拟典型夏季傍晚高峰晚峰最大负荷出现在19时峰值约680MW风电预测的形态是夜间大、白天小光伏预测从6时开始爬升、12时到峰值、18时归零。水电日来水量总量固定为80000m³库容上下限限制在某个范围内第一天初始库容取中间值。运行系数方面火电的煤耗系数取常规的二次函数形式弃风弃光惩罚单价取500元/MWh调峰主动性惩罚系数取50相当于把净负荷峰谷差也折算为成本放进目标函数。4.2 调度计划与出力曲线解读求解完成之后我习惯先画三个图各机组出力堆叠图、净负荷和储能充放曲线、SOC变化曲线。堆叠图能直观看出火电出力是否平稳、水电和储能是否在净负荷陡升时及时补位。典型结果是风电和光伏全部消纳弃电量为零。这是因为弃电惩罚高模型宁可调整火电出力也要把新能源吃进去。水电在傍晚爬坡时段出力明显加大替代了火电的升负荷压力储能则在晚高峰前低价充电、高峰时放电实现了峰谷套利的同时也压低了净负荷峰值。火电出力曲线相比没加调峰主动性惩罚的版本要平滑得多说明“调峰主动性”确实在起作用——模型倾向于让调节能力强的资源先动火电的任务被主动分担了。这里有个细节值得注意如果完全不考虑调峰主动性只优化成本储能可能只在电价低谷充、高峰放虽然经济性最优但净负荷尖峰仍然很陡火电承受的压力没有减轻。加入调峰惩罚之后储能和水电会提前布局、提前响应净负荷曲线平滑了很多。这就是“主动性”三个字的含义不是等指令来了再动而是优化模型根据全局预判提前把调节资源布置在最需要的位置。4.3 多方案对比调峰主动性的量化评估为了量化效果我跑过三个方案对比方案A完全不考虑调峰主动性方案B加入净负荷平滑惩罚方案C再加入调节资源利用率约束。评估指标主要有三个弃风弃光率、净负荷峰谷差、火电出力的标准差。指标方案A方案B方案C弃风弃光率2.3%0%0%净负荷峰谷差410MW336MW322MW火电出力标准差62MW41MW38MW方案B相比方案A峰谷差下降了约18%火电出力标准差下降了约34%。这说明调峰主动性惩罚项的增设效果非常显著火电不再频繁大幅调整系统调峰压力被摊薄到了水电和储能身上。方案C进一步限制调节资源利用率不低于某个阈值峰谷差进一步下降到322MW但代价是水电和储能的调节深度更大寿命损耗相应上升。实际应用中到底选方案B还是C取决于系统里调节资源的规模和调度员的偏好这篇文章的核心代码框架两种方案都能支持。5. 工程调试经验与常见问题速查5.1 求解不收敛或求解过慢的排查路径做这类MILP模型最容易遇到的情况就是用Cplex解了半天不收敛或者干脆报infeasible。我踩过几次坑之后总结了一套排查顺序。第一看变量规模。如果T取24、机组数量超过20台二进制变量数量就会上千模型规模已经不小。此时求解慢是正常的优先考虑聚合策略把同参数的机组聚合成一座虚拟机组能大幅压缩变量数量。第二看数值尺度。功率如果用实际MW值比如几百上千成本可能会到几百万目标函数数值差异跨好几个数量级Cplex的数值稳定性就会变差。我一般统一用标幺值基准功率100MW成本折算成万元/小时收敛速度快很多。第三查约束一致性。功率平衡导致矛盾是最隐蔽的某一时段负荷预测值太低低于所有机组最小出力之和加水电最小出力模型必然不可行。这种情况解法不是调求解器而是允许本时段按一定比例弃电或者增加可中断负荷也就是把“硬约束”降级成“软约束”。此外大规模模型的调试建议先用小规模算例跑通比如T取4、机组取1台确认建模逻辑无误后再放大。这个习惯救过我很多次不然一次丢几百行代码进求解器报错信息根本定位不到问题在哪。5.2 几个绕不开的典型坑SOC递推、启动成本、大M取值储能SOC递推看着简单写错代码很容易。特别是SOC初始值设置的合理性如果第一天的SOC是0.5最后一天的SOC理应回到0.5附近否则连续调度时储能会“贪便宜”把电量放光。常用做法是加上终端约束SOC(T1)SOC(1)但这会牺牲一定的自由度实际运行时未必是最经济的选择。我的折中方案是不加硬约束而是在目标函数里加一个很小的终端SOC偏差惩罚既保证可持续性又不至于过度扰动最优解。启动成本的计算也不能简单写代码因为 (S_{fire}(t)-S_{fire}(t-1)) 虽然可以取值为1表示启动但如果在目标函数里直接乘上启动成本构建出来的表达式包含二进制变量之差逻辑上没问题数值上却容易造成病态。稳妥的写法是引入辅助变量 (SU(t))加约束 (SU(t)\ge S_{fire}(t)-S_{fire}(t-1)) 和 (SU(t)\ge 0)目标函数用 (SU(t)) 乘启动成本这样求解器有明确的凸结构可用。大M法的M值我在储能充放电互斥约束里踩过坑最初取500结果SOC递推变得不稳定缩到50之后模型正常收敛。M值只要比该约束里可能出现最大差值大一点就行没必要取得太大。5.3 代码移植、数据管理经验做调度的课题最后的瓶颈常常不在算法而在数据整理。我建议从一开始就用结构体把所有输入参数封装起来input.time linspace(1, 24, T); input.load P_load; input.wind_pre P_w_pre; input.pv_pre P_s_pre; input.fire struct(max, P_fire_max, min, P_fire_min, ... ramp_up, P_ramp_up, ramp_down, P_ramp_down); input.hydro struct(max, P_hydro_max, inflow, q_in, ... Vmin, V_min, Vmax, V_max, k, k_h); input.storage struct(power, P_st_max, energy, E_max, ... eta_c, eta_c, eta_d, eta_d);这样批量跑场景时只需要改input里的数据不改模型代码。结果输出我用table类型存一次性把时段、火电出力、水电出力、风电、光伏、储能SOC全部存成一个CSV方便后面导入Python或者直接Excel做分析。画图建议用stairs函数画离散调度曲线比plot更符合电力调度的阶梯式出力习惯。写在后面做了一段时间的调度优化一个很深的体会是再漂亮的模型到了工程现场都要跟机组的实际响应速度对齐。储能不是你想充就能瞬间充满的火电也不是你说降负荷立马就能降下来水电还有生态流量的底线。所以这套Matlab代码的价值不是给你一个万能答案而是给你一套可以持续加码的框架——水电多了就加梯级耦合储能比例大了就加电池寿命损耗模型分布式光伏多了就加配电网潮流约束。把“调峰主动性”这个思路真正落地需要你对每一种电源的约束边界、每一处成本的物理意义都有感觉这就是这类课题最磨人也是最有趣的地方。
