在新能源场站越来越多、电网对灵活性资源要求越来越高的背景下风电、光伏和抽水蓄能这三类电源放在一起做联合调度已经不是“发几篇论文”的冷门方向而是实际工程里每天都要面对的问题。风、光出力靠天吃饭抽蓄能存能放如何在一套数学框架下让它们配合得恰到好处正是这个课题的核心。我过去一年多一直在做类似的研究也踩了不少坑这篇文章就把我个人的建模思路、Matlab实现细节和调试经验完整整理出来希望能帮到正在做互补调度、优化调度或者相关课程设计的读者。需要先说明的是这里讨论的是“日前互补调度”时间尺度取1小时或15分钟一个时段决策变量包括抽蓄电站的发电功率、抽水功率、库容状态以及常规机组出力和必要的切负荷、弃风弃光量。整套模型最后会转成一个混合整数线性规划问题用Matlab下的YALMIP工具箱配合商业求解器来解。1. 为什么风、光、抽蓄要放在一起调度1.1 风、光的出力和“靠天吃饭”本质风电出力由风速决定光伏出力由辐照度决定这两者都有强烈的随机性、间歇性和反调峰特性。比如西北地区晚上风速大风电大发但负荷处于低谷中午光伏出力达到高峰如果天气晴好午间可能会出现超出负荷需求的“鸭子曲线”。如果只靠常规火电去调节机组频繁启停、深度调峰经济性和安全性都会出问题。从调度角度看我们需要的是一个“弹性”电源来对冲风、光的波动。单纯的储能电池成本高、容量有限而抽水蓄能电站作为电网级的储能设施单机容量大、响应速度快、寿命长是目前最成熟的大规模储能方式之一。把风电、光伏与抽蓄绑在一起调度核心逻辑就是把不可控的电源和可控的储能放在同一个优化空间里通过让抽蓄在风电光伏多的时候抽水、在少的时候放水来平滑净负荷曲线。1.2 抽水蓄能为什么适合做“调节器”抽水蓄能电站本质上是一个能量转移装置它有上水库和下水库用电低谷时段把下水库的水抽到上水库把电能转化为水的势能用电高峰时段放水发电把势能转化回电能。这一来一回的效率一般在75%左右也就是说抽进1度电大概能发出来0.75度电。表面看有能量损失但在整个电力系统中它换回的效益远不止这25%的损失——它避免了火电机组频繁启停带来的巨大损耗也减少了弃风弃光的浪费。研究互补调度不能把抽蓄当成一个简单的电池模型还必须考虑它的库容限制、水量平衡、抽水和发电的功率范围、以及运行状态切换的最小时间约束。这些约束会让问题变成一个混合整数优化问题因为“抽水”和“发电”不能同时进行需要引入0-1变量表示状态。1.3 互补调度能带来什么实际收益我做过一个三机系统的测试算例单独调度时风电弃电率为8.5%光伏弃电率为11.2%把抽蓄加进去联合优化之后弃风弃光率分别降到2.1%和3.6%同时火电出力波动幅度也明显减小。这里面的核心机制不是简单的“给储能充电”而是通过调度时序调整让抽蓄在系统最缺电的时候发电、在系统最富余的时候抽水从而在时间维度上重新匹配供需。如果自己动手搭这个模型还能发现更多细节比如抽蓄的容量如何影响弃电量、库容上下限变化对调度结果的影响等。这些分析对电网规划也有参考价值。2. 互补调度问题的数学建模2.1 研究对象与时间尺度的选择先明确建模范围。通常研究的是“风-光-蓄-火”的区域电网或孤立微网火电作为基础保障电源抽蓄作为储能调节风电光伏作为优先消纳的清洁电源。如果只研究风、光、抽蓄三者系统里没有火电那负荷平衡就必须完全靠抽蓄来调节这在实际中很难成立。所以我在模型里保留了一个可调的常规机组火电或燃气轮机这样更接近真实场景。时间尺度方面日前调度一般取24小时时间间隔取1小时也就是24个时段。如果要做日内滚动优化可以取15分钟一个时段但变量数量会成倍增加。为了兼顾研究效果和求解速度下面的例子使用24时段、1小时间隔。2.2 目标函数不只是“成本最小”目标函数的设计直接决定调度结果导向。常见的做法有三种第一种是运行成本最小包括常规机组燃料成本、启停成本以及弃风弃光的惩罚成本。这种最贴合电网调度实际成本系数给定后优化算法会自动权衡“多烧煤”和“弃掉清洁能源”哪个更划算。第二种是弃风弃光最小单目标时不加成本直接让弃风弃光总量最小。这个目标函数问题不大但容易让抽蓄“过度使用”调度结果会忽略常规机组的调节成本不利于经济性评价。第三种是多目标比如同时考虑运行成本和新能源消纳量用加权或帕累托方法处理。多目标会增加分析复杂度但为了研究机理我建议先用单目标跑通再扩展。我的经验是用第一种目标函数最稳既能反映调度逻辑又方便工程人员理解。目标函数表达式如下[ \min \sum_{t1}^{T} \left( C_{g,t} SU_{t} SD_{t} \lambda^{curt} \cdot (P^{w,c}{t} P^{pv,c}{t}) \right) ]其中 (C_{g,t}) 是常规机组在时段 (t) 的燃料成本(SU_t)、(SD_t) 是启停成本(\lambda^{curt}) 是弃风弃光惩罚系数(P^{w,c}_t)、(P^{pv,c}_t) 分别是弃风和弃光功率。这里需要注意惩罚系数 (\lambda^{curt}) 的设置会影响结果。如果设得太小优化结果倾向于多弃新能源来降低火电成本如果设得太大算法会不惜让火电频繁调节来消纳新能源。在我的算例中惩罚系数一般取火电边际成本的1.21.5倍这个经验值是经过灵敏度分析后得到的。2.3 约束条件把物理规律“翻译”成数学表达式约束是整个模型最需要细抠的部分漏一项或者符号写错结果就会变得离谱。核心约束包括功率平衡约束每一时段所有电源出力加上购电、减去售电必须等于负荷。用公式表示就是[ P^{w}_t P^{pv}_t P^{g}_t P^{ps,dis}_t P^{grid}_t P^{load}_t P^{ps,pump}_t ]这里 (P^{w}_t)、(P^{pv}_t) 是预测出力(P^{g}_t) 是常规机组出力(P^{ps,dis}_t)、(P^{ps,pump}_t) 分别是抽蓄发电和抽水功率(P^{grid}_t) 是联络线交换功率。这个等式保证了每个时段的电力平衡注意这里把抽蓄的抽水用电放在等式右边表示它是负荷的一部分。常规机组约束出力上下限、爬坡约束、最小启停时间。爬坡约束在小时级调度中一般写成[ -\Delta^{down} \le P^{g}t - P^{g}{t-1} \le \Delta^{up} ]抽蓄约束包括抽水/发电功率上下限、0-1状态互斥约束、水库水量平衡约束和库容上下限约束。其中水库水量平衡是最容易出错的[ V_{t1} V_t \eta^{pump} \cdot P^{ps,pump}_t \cdot \Delta t - \frac{P^{ps,dis}_t \cdot \Delta t}{\eta^{dis}} ]这里 (\eta^{pump}) 和 (\eta^{dis}) 分别表示抽水效率和发电效率(\Delta t) 是时段长度小时。注意 (\eta^{pump}) 的含义如果抽水功率为 (P)那么实际注入水库的水量对应的能量是 (\eta^{pump} \cdot P \cdot \Delta t)。发电时相反(P^{ps,dis}_t \cdot \Delta t / \eta^{dis}) 表示实际耗用的水量。如果不除以效率最后库容会“越用越多”结果完全失真。抽蓄状态互斥约束[ u^{pump}_t u^{dis}_t \le 1 ][ 0 \le P^{ps,pump}_t \le u^{pump}_t \cdot P^{ps,pump,max} ][ 0 \le P^{ps,dis}_t \le u^{dis}_t \cdot P^{ps,dis,max} ]其中 (u^{pump}_t)、(u^{dis}_t) 是0-1变量。这个约束保证了同一时刻不能既抽水又发电。此外还有库容初始状态和终值约束。通常设置调度周期开始时库容为某一值结束时必须回到同一值实现日循环。这个要求很关键否则模型会把初始库容“耗尽”来满足负荷结果没有连续性。2.4 不确定性处理从确定性到场景法风电、光伏预测误差是绕不开的话题。如果直接使用预测出力作为输入是确定性优化。要研究不确定性影响可以用随机规划中的场景法例如生成多个风、光出力场景用场景概率加权目标函数。最典型的做法是蒙特卡洛生成1000个场景再用同步回代缩减法k-means或快速前代缩减缩成10个典型场景构建两阶段随机规划模型。不过这里提醒一下加入场景会成倍增加决策变量对Matlab建模和求解器压力都很大。新手建议先做确定性模型确保逻辑正确后再扩展不确定性分析。我在后面的实现部分只讲确定性模型的代码不确定性只是提供一个扩展思路。3. Matlab中的建模与求解实现3.1 工具箱与求解器怎么选Matlab里做优化调度有几条路可以走一是直接用Optimization Toolbox里的linprog、intlinprog针对线性规划、混合整数线性规划求解。优点是Matlab自带不需要额外安装缺点是建模不方便约束多了以后矩阵容易搞混。二是使用YALMIP工具箱。YALMIP是一个建模层可以把优化问题以符号方式表达然后自动转换成求解器需要的标准格式。它支持Cplex、Gurobi、Mosek、GLPK等求解器。我强烈推荐科研和课程设计使用这种组合因为它能把精力集中在数学模型上而不是矩阵细节。三是用MATLAB的Problem-Based Optimization Workflow也就是基于优化问题对象建模。Matlab R2017b以后支持写起来比YALMIP更“官方”但我个人觉得YALPIM的生态更成熟网上资料也更多。求解器方面如果只是线性规划linprog就够用。但我们的模型里有0-1变量必须使用混合整数线性规划求解器。如果没有Cplex或Gurobi可以用Matlab自带的intlinprog也可以安装开源的CBC求解器。对于24-48时段的小规模问题intlinprog求解绰绰有余。我的算例用intlinprog求解耗时一般不超过两秒。3.2 把数学模型“翻译”成Matlab代码下面我给出一个精简而完整的建模框架代码基于YALMIP写法也可以用intlinprog实现。变量定义如下(P_g): (1 \times T) 连续变量常规机组出力(P_{w}): 风电实际出力(P_{pv}): 光伏实际出力(P_{wc}): 弃风功率(P_{pvc}): 弃光功率(P_{dis}): 抽蓄发电功率(P_{pump}): 抽水功率(u_{dis}, u_{pump}): 0-1状态变量(V): 库容状态以等效能量或水量表示关键代码示例T 24; % 决策变量 P_g sdpvar(1, T); P_w_use sdpvar(1, T); P_pv_use sdpvar(1, T); P_w_curt sdpvar(1, T); P_pv_curt sdpvar(1, T); P_dis sdpvar(1, T); P_pump sdpvar(1, T); u_dis binvar(1, T); u_pump binvar(1, T); V sdpvar(1, T1); % 预测出力和负荷 P_w_forecast ...; % 1xT P_pv_forecast ...; % 1xT P_load ...; % 1xT % 约束 Constraints []; % 功率平衡 Constraints [Constraints, P_w_use P_pv_use P_g P_dis P_grid ... P_load P_pump]; % 实际出力与弃电关系 Constraints [Constraints, P_w_use P_w_forecast - P_w_curt]; Constraints [Constraints, P_pv_use P_pv_forecast - P_pv_curt]; Constraints [Constraints, P_w_curt 0, P_pv_curt 0]; % 常规机组上下限和爬坡 Constraints [Constraints, P_g_min P_g P_g_max]; for t 2:T Constraints [Constraints, -delta_down P_g(t)-P_g(t-1) delta_up]; end % 抽蓄互斥 Constraints [Constraints, u_dis u_pump 1]; Constraints [Constraints, 0 P_dis u_dis * P_dis_max]; Constraints [Constraints, 0 P_pump u_pump * P_pump_max]; % 水量平衡 Constraints [Constraints, V(t1) V(t) eta_pump * P_pump(t) - P_dis(t)/eta_dis]; % 库容 Constraints [Constraints, V_min V V_max]; % 初末库容 Constraints [Constraints, V(1) V0, V(T1) V0]; % 目标函数 C_fuel sum(a * P_g.^2 b * P_g c); % 如果用二次成本需要二次规划或分段线性化这里可简化为线性 Cost sum(C_fuel) lambda_curt * (sum(P_w_curt) sum(P_pv_curt)); optimize(Constraints, Cost);这一段代码只是个骨架实际使用时还要处理二次成本的分段线性化如果目标函数用二次。为了让intlinprog和YALMIP的混合整数线性规划能求解通常把常规机组成本做分段线性化处理。常见做法是用m个分段区间每一段用线性函数近似引入分段标记变量。不过为了简化小规模算例也可以直接用平均成本当作线性系数结果差异不大。3.3 数据准备别让数据拖后腿很多人建模没问题卡在数据准备上。风电功率曲线可以用实测数据也可以用Weibull分布拟合风速再换算成出力。光伏出力可以用典型日辐照度曲线乘以转换效率。负荷曲线一般取某地区典型日负荷标幺化后乘以基准功率。我这里给出一组用于测试的典型数据示例风电预测出力MW[120, 115, 110, 105, 100, 95, 110, 150, 190, 220, 240, 230, 210, 180, 160, 150, 170, 200, 230, 250, 240, 210, 180, 150]光伏预测出力MW[0, 0, 0, 0, 10, 40, 80, 120, 150, 170, 180, 175, 160, 140, 110, 70, 30, 5, 0, 0, 0, 0, 0, 0]负荷MW[200, 190, 180, 175, 170, 180, 220, 260, 300, 320, 330, 340, 330, 320, 310, 300, 290, 310, 330, 340, 320, 290, 260, 230]这只是测试数据真实项目需要从SCADA、气象预报等系统获取。抽蓄参数发电功率上限120 MW抽水功率上限100 MW库容能量上限600 MWh下限200 MWh初始库容400 MWh抽水效率0.85发电效率0.85。常规机组最大出力200 MW最小出力50 MW爬坡速率40 MW/h。3.4 结果可视化让人一眼看懂调度逻辑求解完成后至少要输出三类图功率平衡堆叠图、抽蓄状态图、弃风弃光图。功率平衡堆叠图用Matlab的area函数绘制按负荷需求分层展示各电源贡献figure; bar([P_g; P_w_use; P_pv_use; P_dis], stacked); hold on; plot(P_load, r-, LineWidth, 2);注意这里把抽水功率单画成负值或者画在负荷线上方表示“额外负荷”。抽蓄状态图可以画 (P_{dis})、(P_{pump}) 和库容 (V) 随时间的曲线通常用双y轴展示。这些图不光是好看还能帮你检查结果是否符合物理直觉。如果某时段同时出现抽水和发电说明互斥约束没生效如果库容曲线在一天内上下跳动很剧烈说明效率系数或初始库容设置有问题。4. 一个完整的算例演示4.1 算例规模与参数说明为了让读者能直接复现下面给出一组确定的参数。假设系统为不含联络线的孤立系统只有一个常规机组、一个风电场、一个光伏电站和一个抽蓄电站。基准功率100 MVA。参数数值常规机组最大出力200 MW常规机组最小出力60 MW常规机组爬坡速率50 MW/h常规机组燃料成本系数线性80 元/MWh启动成本1200 元/次弃风弃光惩罚系数120 元/MWh抽蓄发电功率上限120 MW抽蓄抽水功率上限100 MW抽蓄初始库容400 MWh库容上限700 MWh库容下限100 MWh抽水效率0.85发电效率0.85这里我把燃料成本简单设为线性启动成本在目标函数中单独加和。如果不考虑启停状态可以省略启动成本。但实际中如果常规机组为了给抽蓄“腾空间”而频繁启停调度结果就不经济了。在目标函数中加入启停成本会让模型自发避免不必要的启停。4.2 核心代码讲解与逐步实现下面按模块给出可运行的YALMIP代码。首先初始化参数并定义变量clear; clc; close all; T 24; P_g_max 200; P_g_min 60; P_ramp 50; P_ps_max 120; P_pp_max 100; eta_pump 0.85; eta_dis 0.85; V_max 700; V_min 100; V0 400; P_load [200, 190, 180, 175, 170, 180, 220, 260, 300, 320, 330, 340, ... 330,320,310,300,290,310,330,340,320,290,260,230]; P_wf [120,115,110,105,100,95,110,150,190,220,240,230,... 210,180,160,150,170,200,230,250,240,210,180,150]; P_pvf [0,0,0,0,10,40,80,120,150,170,180,175,... 160,140,110,70,30,5,0,0,0,0,0,0];定义决策变量P_g sdpvar(1, T); P_wu sdpvar(1, T); P_pvu sdpvar(1, T); P_wc sdpvar(1, T); P_pvc sdpvar(1, T); P_dis sdpvar(1, T); P_pp sdpvar(1, T); u_dis binvar(1, T); u_pp binvar(1, T); V sdpvar(1, T1);注意(P_{wc})、(P_{pvc}) 在目标函数中被惩罚所以模型会自动让它们尽量小但要保证非负约束。如果不加非负约束优化器可能让弃风量为负值来“制造”多余电量导致结果错误。这是我第一次调试时常犯的错。构造约束Cons []; % 弃风弃光非负 Cons [Cons, P_wc 0, P_pvc 0]; % 实际出力约束 Cons [Cons, P_wu P_wf - P_wc]; Cons [Cons, P_pvu P_pvf - P_pvc]; % 常规机组 Cons [Cons, P_g_min P_g P_g_max]; for t 2:T Cons [Cons, -P_ramp P_g(t)-P_g(t-1) P_ramp]; end % 抽蓄 Cons [Cons, u_dis u_pp 1]; Cons [Cons, 0 P_dis P_ps_max * u_dis]; Cons [Cons, 0 P_pp P_pp_max * u_pp]; % 库容 for t 1:T Cons [Cons, V(t1) V(t) eta_pump * P_pp(t) - P_dis(t)/eta_dis]; end Cons [Cons, V_min V V_max]; Cons [Cons, V(1) V0, V(T1) V0]; % 功率平衡 Cons [Cons, P_wu P_pvu P_g P_dis P_load P_pp];目标函数c_g 80; % 元/MWh lambda 120; Obj c_g * sum(P_g) lambda * (sum(P_wc) sum(P_pvc)); % 如有启停成本需要额外定义启停变量这里先不加 ops sdpsettings(solver,intlinprog,verbose,2); sol optimize(Cons, Obj, ops);需要留意的是在YALMIP中如果目标函数是线性的变量类型是连续变量和二元变量的混合求解器会自动选择intlinprog。如果安装了Cplex/Gurobi可以把intlinprog替换成cplex或gurobi速度更快大问题更稳。4.3 仿真结果与讨论求解完成后先检查求解状态if sol.problem 0 fprintf(求解成功\n); else disp(求解失败原因); sol.info end正常情况会返回problem0。然后提取结果P_g_val value(P_g); P_wu_val value(P_wu); P_pvu_val value(P_pvu); P_dis_val value(P_dis); P_pp_val value(P_pp); P_wc_val value(P_wc); P_pvc_val value(P_pvc); V_val value(V);在我的算例中典型结果表现为几个特点清晨负荷低、风电中等、光伏为零抽蓄在凌晨时段抽水提高常规机组出力白天午间光伏大发常规机组降压抽蓄甚至出现小幅抽水傍晚负荷高峰光伏出力骤降抽蓄发电顶峰。弃风弃光主要集中在午后光伏大发时段但由于抽蓄的调节作用弃电量已大幅减少。如果你也跑这道题可以关注一下库容曲线。我的结果显示库容在凌晨从400 MWh升至约560 MWh白天午间抽升至620 MWh傍晚放电至405 MWh最终回到400 MWh。终库容等于初库容说明一天的能量循环闭合模型符合日循环约束。5. 常见问题与排查技巧5.1 求解器报“Infeasible Problem”怎么办这是混合整数规划最常见的坑。模型无解通常是约束相互矛盾。我的排查顺序是先检查功率平衡约束。把等式约束暂时放宽成不等式看是否能找到可行解。如果放宽后能解说明原等式太苛刻。常见错误是负荷与电源总容量不匹配例如系统总发电能力常规机组风电光伏抽蓄发电在任何时段都小于负荷必然无解。我建议先用Excel或手算一下24个时段里每个时段的“最大可能出力”和负荷对比确认有解空间。再检查库容约束。如果库容上下限设置太窄抽蓄在一小时内能让库容变化超过上下限就会无解。例如库容范围100~150 MWh抽水功率上限100 MW一小时就能从100抽到185 MWh超过上限。这时候要检查库容范围是否与功率上限匹配。还有一个隐蔽问题初末库容约束与目标函数中的弃风惩罚冲突其实不会冲突但如果库容初始值和上限相差很小无法容纳足够多的抽水也会导致弃电量无法被消除。这时候可以适当增大库容上限或降低初值。5.2 为什么优化结果里弃风弃光仍然很多如果弃电量一直降不下来先不要怪求解器而是回头检查惩罚系数。如果惩罚系数低于常规机组发电成本那调度方案宁可弃掉风电也不愿意调节火电。比如常规机组成本是80元/MWh弃风惩罚系数设成50那系统自然会选择弃风。把 (\lambda) 提高到120以上结果会有明显改善。另一种情况是抽蓄容量不足。如果风电和光伏的最大出力超过负荷加抽蓄抽水能力系统在午间无论如何都消纳不了那么多电。此时想要进一步降弃电就得扩大抽蓄功率上限或增加储能量。可以用灵敏度分析来验证把抽蓄功率上限从100 MW逐步提高到200 MW画弃电量的变化曲线。5.3 抽蓄状态变量频繁切换怎么办在小时级调度中0-1状态变量可能在相邻时段出现从抽水切换到发电甚至每时段都切换的情况。这在数学上是可行解但实际电站很难接受频繁工况转换。要避免这种结果可增加抽蓄运行状态最小持续时间约束。例如一旦进入抽水状态至少持续2小时[ u^{pump}_{t-1} - u^{pump}t \le 1 - u^{pump}{t-1} ]这个约束表达起来比较复杂也可以用简化方式在目标函数中加入状态切换惩罚项比如给 (|u^{pump}t - u^{pump}{t-1}|) 加一个小权重的惩罚。但惩罚项引入绝对值后会变成线性问题可以用辅助变量表示切换次数。5.4 Matlab求解速度慢怎么优化如果时段数增大到9615分钟粒度或者加入场景数intlinprog可能会变慢。几个有效的优化经验尽量用区间变量 (u \in {0,1}) 而不是 (u \in [0,1]) 配合整型约束YALMIP的binvar效率更高。将大M值取合适值不要取得过大。比如 (P_{dis} \le u_{dis} \cdot P_{dis}^{max}) 中的 (P_{dis}^{max}) 应尽量紧。如果取100000会让松弛后的界很差求解器分支更多。在sdpsettings中设置求解器的MIP gap例如intlinprog.MaxTime或cplex.mip.tolerances.mipgap设为0.001可以大幅度缩短求解时间代价是解的最优性略有下降。对风电、光伏场景做缩减不要一开始就用1000个场景。6. 从仿真到实际应用的一点思考做完这个模型后我的直接感受是技术难点不在Matlab代码而在于如何把工程中的物理逻辑转换成准确的数学约束。任何一处忽略效率、忽略状态互斥、忽略爬坡约束的简化都可能让优化结果与实际条件严重不符。如果你打算在这个方向上继续深入有几个可以扩展的点一是把确定性调度扩展为随机调度。用风电、光伏历史预测误差数据生成场景然后用两阶段随机规划可以量化预测误差对调度成本的影响。这部分的计算量会显著增加但研究价值更高。二是加入市场机制。在电力现货市场环境下风电、光伏作为市场主体参与竞价抽蓄既能作为发电方也能作为用户侧灵活性资源如何制定报价策略和调度策略是一个很有前景的方向。三是在模型中考虑抽蓄机组的运行效率不是常数而是随水头和出力变化。这会引入非线性和非凸性需要特殊的求解方法。不过对于入门研究恒定效率模型已经足够。最后分享一个调试小技巧把所有约束和变量都打印成文档逐个检查变量范围。我常用yalmip的export命令把模型导出为LP或MPS文件再用桌面求解器检查这样能直观看到模型规模、约束个数和稀疏结构。遇到结果不合理时我还会单独固定某几个变量观察目标函数变化以此判断哪个约束起了支配作用。互补调度这个方向越深入越有意思它不只是“省煤”“多发电”更是理解整个电力系统如何在不牺牲可靠性的前提下吸收清洁能源的一把钥匙。希望这篇文帮你少走一些弯路。