最近在做一个风光蓄互补调度的小项目把风电、光伏和抽水蓄能电站放进同一个优化框架里用Matlab跑日前调度。这套模型我前后改了五版从最开始的风光单独消纳到后来把抽蓄的启停、抽发状态全部用整数变量建模才算是把弃风弃光和削峰填谷的效果同时抓了出来。今天把完整的调度思路、数学模型、Matlab实现和调试过程中踩过的坑整理一遍给做电力系统优化、微电网调度或者新能源消纳研究的同学参考。这套代码的思路并不复杂但里面涉及的目标函数设置、约束矩阵搭建、求解器参数选择很多都是教科书不会细讲的东西。1. 互补调度的整体设计思路1.1 为什么是风电、光伏加抽水蓄能这个组合风电和光伏在时间尺度上的出力特性其实挺互补的。风电往往在夜间和凌晨出力偏大光伏则集中在白天中午前后两者叠加以后净负荷曲线的波动幅度依然很大尤其是早晚负荷高峰时段新能源出力往往帮不上忙甚至会出现光伏中午大发、负荷却处于午谷需要大量调峰容量的情况。这时候如果有一个能够快速充电和放电的调节电源就能把低谷时段的富余新能源电量存起来等到高峰时段再放出来这就是抽水蓄能在这个模型里扮演的角色。抽水蓄能相比电化学储能最大的优势是容量大、单位容量成本低、寿命长适合日内甚至周尺度的能量搬移。响应速度虽然不如电池快但用来平抑小时级的净负荷波动完全够用。而且抽蓄电站本身具有同步电机特性还能提供惯量和无功支撑这些附加价值在纯经济调度里不容易体现但在实际电网调度中是实打实的收益。所以做风电、光伏与抽水蓄能互补调度不是简单的“加一个储能”而是要把抽蓄的抽水、发电、停机三种状态和库容的日循环过程一起纳入优化。1.2 调度研究的核心命题是什么这个课题的落脚点是回答一个很实际的问题在已知风电、光伏预测出力和负荷预测曲线的前提下抽水蓄能电站该怎么抽、什么时候发、一天之内要抽多少水才能让新能源消纳最多、系统峰谷差最小、运行也最经济。时间尺度上常规做法是给定一个调度周期典型的是24小时时间分辨率取15分钟或1小时。15分钟对应96个时段1小时对应24个时段。分辨率越细越能反映风光出力的波动细节但优化问题的变量数和约束数会成倍增加。我做项目时默认用96时段既保留了足够的时间细节又不至于让求解器跑不动。调度过程本身是开环的即用预测数据进行日前决策。实际工程中还需要日内滚动修正但作为研究模型先把确定性日前调度做扎实再往后考虑不确定性是比较稳妥的技术路线。很多人一上来就想把鲁棒优化、随机优化全部塞进去结果模型复杂到根本找不到头绪反而失去了一次把基础逻辑吃透的机会。1.3 建模方案选型混合整数线性规划为什么是首选互补调度的核心难点在于抽蓄电站的离散运行状态同一台机组不能同时抽水和发电甚至同一时段整个电站往往只允许处于一种工况。这种“互斥逻辑”天然需要0-1整数变量来表达。再加上功率平衡、库容动态、出力上下限都是线性关系整个问题就自然落到了混合整数线性规划MILP这个框架里。MILP的好处有三个。第一有成熟的商业求解器比如Gurobi、Cplex开源的有CBCMatlab自带的intlinprog也能处理中小规模问题。第二线性模型目标函数和约束都容易调试出问题的时候可以逐条检查。第三求解结果具有全局最优性不会像启发式算法那样陷入局部最优对做研究发表结论来说更有说服力。有些同学会问用动态规划或者遗传算法行不行动态规划处理单台抽蓄机组可以但一旦扩展成多机组、多水库状态维数就会爆炸。遗传算法等启发式算法在论文里经常看到但缺乏最优性保证参数调起来也麻烦。我的建议是除非问题规模大到MILP实在解不动否则优先上MILPMatlab里直接调用求解器效率和可靠性都高得多。2. 数学模型与约束条件2.1 目标函数设计消纳优先还是经济优先调度目标决定了优化方向。如果系统里有火电通常目标会写成火电煤耗成本最小化或购电成本最小化新能源和抽蓄则作为可调资源进入约束。但在纯风光蓄系统里没有火电做比较基准比较自然的做法是最小化弃风弃光量同时给抽蓄的运行成本加一个很小的惩罚系数防止求解器无意义地频繁启停抽蓄机组。目标函数可以这样写[ \min \sum_{t1}^{T} \left[ \alpha_w \cdot P_{w,curt}(t) \alpha_v \cdot P_{v,curt}(t) \beta \cdot (P_{pump}(t) P_{gen}(t)) \right] ]其中 (P_{w,curt}(t)) 和 (P_{v,curt}(t)) 分别表示t时段风电和光伏的弃电功率(\alpha_w)、(\alpha_v) 是弃电惩罚系数(\beta) 是抽蓄运行成本惩罚系数。弃风弃光惩罚系数要远大于 (\beta)这样才能保证优化器优先消纳新能源而不是图省事直接把预测出力砍掉。还有一种做法是直接最大化新能源实际并网电量即目标函数写成负的 ( (P_{w,use} P_{v,use}) ) 之和。两种写法在数学上等价但用弃电量作为变量更直观后续统计弃风率、弃光率时也方便。2.2 功率平衡与新能源出力约束不管目标函数怎么定系统约束都必须严格满足。首先是每个时段的功率平衡[ P_{load}(t) P_{w,use}(t) P_{v,use}(t) P_{gen}(t) - P_{pump}(t) ]其中 (P_{load}(t)) 是系统净负荷风电、光伏的实际并网功率加上抽蓄发电功率减去抽水消耗的功率必须时刻等于负荷。这个约束用等式形式写入优化模型是所有变量联动的主线索。新能源侧实际并网功率不能超过预测功率[ 0 \le P_{w,use}(t) \le P_{w,fore}(t) ] [ 0 \le P_{v,use}(t) \le P_{v,fore}(t) ]弃电量定义为预测值减去实际值[ P_{w,curt}(t) P_{w,fore}(t) - P_{w,use}(t) ] [ P_{v,curt}(t) P_{v,fore}(t) - P_{v,use}(t) ]只要预测出力上限给得合理弃电量会自动落在非负区间。2.3 抽水蓄能电站运行约束最关键的建模环节抽蓄的约束是整个模型的核心也是最容易写错的地方。首先运行状态互斥同一时段不能既抽水又发电[ u_{pump}(t) u_{gen}(t) \le 1 ]其中 (u_{pump}(t)) 和 (u_{gen}(t)) 是0-1变量。然后功率上下限与状态变量绑定[ 0 \le P_{pump}(t) \le P_{pump}^{max} \cdot u_{pump}(t) ] [ 0 \le P_{gen}(t) \le P_{gen}^{max} \cdot u_{gen}(t) ]如果抽蓄电站由多台机组组成可以按机组聚合或单机建模。单机建模更精细但整数变量数量会成倍增加求解压力大。作为研究模型通常可以先按单机等效处理后续再细化。库容动态是另一个重要约束。如果用水量体积表示库容还需要换算成功率的积分关系容易写出非线性项。这里建议直接把水库储存量归一化为“可发电能量” (E(t))单位取MWh这样公式更简洁[ E(t1) E(t) \eta_p \cdot P_{pump}(t) \cdot \Delta t - \frac{P_{gen}(t) \cdot \Delta t}{\eta_g} ](\eta_p) 是抽水效率(\eta_g) 是发电效率(\Delta t) 是单个时段的小时数。这个方程的含义是抽水时把电能转化为水的势能存入水库发电时再把势能转化为电能释放出来转换过程有损耗。考虑到抽蓄电站综合效率一般在0.75到0.8左右抽水和发电效率可以分别取0.9和0.85乘积约0.765。库容上下限约束[ E_{min} \le E(t) \le E_{max} ](E_{min}) 不能取0因为水库存在死库容不能真正放空(E_{max}) 则由水库有效库容决定。日循环约束是研究日内调度时常用的一招[ E(T1) \ge E(1) ]或者更严格地取等号要求调度周期结束时库容不小于初始库容从而保证第二天还能继续调度。如果允许跨日调节可以把末库容做成决策变量但要加入惩罚或给定初末值。2.4 决策变量汇总为了方便编写代码通常把所有决策变量整理成一个向量。完整问题的决策变量包括变量类型含义(P_{w,use}(t))连续风电实际并网功率(P_{v,use}(t))连续光伏实际并网功率(P_{w,curt}(t))连续风电弃电功率(P_{v,curt}(t))连续光伏弃电功率(P_{pump}(t))连续抽水功率(P_{gen}(t))连续发电功率(E(t))连续水库储能水平(u_{pump}(t))0-1抽水状态(u_{gen}(t))0-1发电状态如果系统还有联络线功率、火电出力就需要增加对应变量。变量的索引方式决定了约束矩阵的构造方式建议先画一张变量布局图再写矩阵能省掉很多调bug的时间。3. Matlab代码实现与求解细节3.1 代码结构脚本模块化别全塞一个文件里Matlab实现求解MILP有两种常见方式一种是直接用内置的intlinprog用矩阵和向量描述目标与约束另一种是安装Yalmip工具箱用符号化的方式搭建优化模型再把模型交给Gurobi、Cplex或intlinprog求解。我强烈建议研究阶段用Yalmip因为模型可读性高改约束方便但如果要做成自动化程序或嵌入其他系统intlinprog的更底层写法反而少一层依赖。代码整体可以分成几个模块loadData.m读取风电、光伏、负荷预测数据处理缺失值生成时段序列。buildModel.m定义决策变量、目标函数和约束条件。solveModel.m调用求解器设置求解参数输出结果。plotResults.m绘制功率平衡、库容变化、弃电情况等曲线。runMain.m主程序按顺序执行上述流程。数据读取这一步整理成结构体或表格比较方便。比如用readtable读取Excel文件后再通过interp1把数据插值成需要的分辨率。预测数据如果是从别处仿真生成的可以提前绘制出来检查曲线是否合理避免后期求解结果被脏数据带偏。3.2 基于Yalmip的快速建模代码直观半小时跑通模型下面给出一段简化但完整可运行的Yalmip建模代码时间分辨率取15分钟T96目标是最大化新能源消纳并附加抽蓄运行惩罚。T 96; dt 0.25; % 小时数15分钟 P_load load_profile; % 1x96 负荷预测MW P_w_fore wind_forecast; % 1x96 风电预测 P_v_fore pv_forecast; % 1x96 光伏预测 P_p_max 400; % 最大抽水功率MW P_g_max 400; % 最大发电功率MW E_min 200; % 最小储能MWh E_max 2000; % 最大储能MWh E_init 1000; % 初始储能MWh eta_p 0.9; eta_g 0.85; P_w sdpvar(1, T); % 风电实际出力 P_v sdpvar(1, T); % 光伏实际出力 P_p sdpvar(1, T); % 抽水功率 P_g sdpvar(1, T); % 发电功率 E sdpvar(1, T1); % 储能水平 u_p binvar(1, T); % 抽水状态 u_g binvar(1, T); % 发电状态 Constraints []; for t 1:T Constraints [Constraints, 0 P_w(t) P_w_fore(t)]; Constraints [Constraints, 0 P_v(t) P_v_fore(t)]; Constraints [Constraints, P_load(t) P_w(t) P_v(t) P_g(t) - P_p(t)]; Constraints [Constraints, 0 P_p(t) P_p_max * u_p(t)]; Constraints [Constraints, 0 P_g(t) P_g_max * u_g(t)]; Constraints [Constraints, u_p(t) u_g(t) 1]; Constraints [Constraints, E_min E(t) E_max]; Constraints [Constraints, E(t1) E(t) eta_p * P_p(t) * dt - P_g(t) * dt / eta_g]; end Constraints [Constraints, E(1) E_init]; Constraints [Constraints, E(T1) E_init]; alpha 10; % 弃电惩罚系数 beta 0.01; % 抽蓄运行惩罚系数 Objective alpha * sum(P_w_fore - P_w) alpha * sum(P_v_fore - P_v) ... beta * sum(P_p P_g); ops sdpsettings(solver, gurobi, verbose, 2); diagnostic optimize(Constraints, Objective, ops);这里有几个细节值得注意。约束写成循环虽然直观但T很大时循环效率会下降。Yalmip会自动把循环拆成向量运算所以性能尚可。如果坚持用intlinprog就需要把所有约束手动拼成大矩阵代码长很多但运行速度会更快。目标函数中alpha和beta的量级关系很重要。弃电惩罚必须远大于抽蓄运行成本否则求解器会在弃电和抽蓄耗电之间摇摆不定。我一般设alpha10beta0.01这样每弃1MWh相当于抽蓄运行1000MWh的代价优化器自然优先消纳新能源。3.3 如果不用Yalmipintlinprog矩阵化怎么写不使用Yalmip时所有约束都要转化成 (A x \le b) 和 (A_{eq} x b_{eq}) 的形式。决策变量向量x需要按固定顺序拼接比如% 变量分段索引P_w(1:T), P_v(1:T), P_p(1:T), P_g(1:T), E(1:T1), u_p(1:T), u_g(1:T) n_w T; n_v T; n_pump T; n_gen T; n_e T1; n_up T; n_ug T; idx_w 1:T; idx_v idx_w(end) (1:T); idx_p idx_v(end) (1:T); idx_g idx_p(end) (1:T); idx_e idx_g(end) (1:T1); idx_up idx_e(end) (1:T); idx_ug idx_up(end) (1:T); n idx_ug(end);然后分别初始化f、A、b、Aeq、beq、lb、ub。比如功率平衡约束在每一行对应的Aeq中把P_w、P_v、P_g位置的系数置为1P_p位置置为-1右侧等于P_load(t)。库容动态约束则把E(t1)位置置1E(t)位置置-1P_p(t)位置减去eta_p*dtP_g(t)位置加上dt/eta_g右侧为0。状态耦合约束通过A不等式实现比如P_p(t) - P_p_max*u_p(t) 0。最后调用intcon [idx_up, idx_ug]; % 整数变量索引 options optimoptions(intlinprog, Display, final, AbsoluteGapTolerance, 1e-4, RelativeGapTolerance, 1e-4); [x, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, options);intlinprog把整数变量索引单独作为第二个输入参数其余连续变量默认可以取实数。需要特别注意变量数量越多矩阵越容易写错建议每加一个类型的约束就用小规模随机数据验证一次Aeq的维度是否匹配。3.4 结果可视化画图不只是给论文用的调度结果画图不仅是展示更是排查模型问题的重要手段。我习惯至少画三张图第一张是功率平衡堆叠图可以看出风电、光伏、抽蓄发电和负荷之间的时间匹配第二张是库容变化曲线检查储能是否触顶触底第三张是抽蓄状态图用stairs或分段线段来表示抽水、发电、停机状态。堆叠图可以直接用Matlab的area命令figure; area(t, [P_w_use, P_v_use, P_gen], LineWidth, 0.5); hold on; plot(t, P_load, k-, LineWidth, 1.5); legend(风电, 光伏, 抽蓄发电, 负荷); xlabel(时刻 (h)); ylabel(功率 (MW));库容图用plot和ylim就可以。导出图片时建议用exportgraphics(gcf, result.png, Resolution, 300)比直接print更清晰。如果需要在论文里用矢量图导出成eps或pdf效果更好。4. 仿真结果与效果评估4.1 算例参数设定先定好边界再跑数据为了检验模型我设置了一个中等规模算例风电装机500MW光伏装机300MW抽蓄电站最大抽水/发电功率400MW水库有效库容对应2000MWh可发电能量初始库容1000MWh抽水效率0.9发电效率0.85。负荷曲线取典型工商业日负荷峰值约800MW谷值约450MW。参数数值说明风电装机500 MW预测出力上限光伏装机300 MW预测出力上限抽蓄功率400 MW抽水和发电最大功率库容上限2000 MWh等效可发电能量库容下限200 MWh死库容等效初始库容1000 MWh调度起点抽水效率0.9电能转势能发电效率0.85势能转电能时间分辨率15 min96时段/天数据生成时可以用历史出力曲线的缩放也可以模拟一组带随机波动的出力序列。模拟的好处是能控制风光的互补特性比如让风电夜间出力大、光伏中午出力大。使用真实数据时需要注意数据的时间戳是否对齐负荷、风电、光伏必须处于同一时区、同一分辨率。4.2 典型日调度结果抽蓄到底在什么时段动作跑完模型后最直观的结果是抽蓄的“低抽高发”模式。在一组典型日数据下中午12点到14点光伏大发、负荷相对平稳抽蓄开始抽水把多余的光伏电量存起来晚上18点到20点负荷快速攀升、光伏出力已经降为零抽蓄转入发电弥补晚高峰的电力缺口。风电出力大但负荷较低的凌晨时段抽蓄同样会抽水储备能量。从数据看加入抽蓄后弃风弃光电量从无抽蓄时的280MWh下降到约45MWh弃电率从8%左右降到1.3%左右。系统净负荷的峰谷差也明显减小说明抽蓄确实起到了能量搬移的作用。曲线图上最明显的特征就是负荷线不再大起大落抽蓄发电时段恰好顶在负荷爬坡最陡的位置。4.3 互补效果评估指标不要只看弃电率弃风弃光率是核心指标但不是全部。调度效果通常还需要从削峰填谷、净负荷波动、抽蓄利用率几个维度综合评估。削峰率可以这样计算[ R_{peak} \frac{(P_{load,max} - P_{load,min}) - (P_{net,max} - P_{net,min})}{P_{load,max} - P_{load,min}} \times 100% ]其中 (P_{net}) 是扣除了风光实际出力后的净负荷即 (P_{load} - P_{w,use} - P_{v,use})。抽蓄的目标是把净负荷峰谷差压下来。净负荷标准差则能反映波动程度标准差的下降比例比峰谷差更稳健。某次仿真中典型日原始负荷峰谷差约350MW经抽蓄调度后净负荷峰谷差约195MW削峰率约44%。抽蓄日发电量约1800MWh对应等效利用约4.5小时还算是比较健康的运行方式。如果抽蓄频繁启停或发电量过高反而会造成水库频繁穿越上下限需要在目标函数里加大对状态切换的惩罚。4.4 敏感性分析与扩展思路模型跑通后最值得做的事情是敏感性分析。固定其他参数分别改变抽蓄功率容量、库容上限、综合效率观察弃电率如何变化。我试过把库容从1000MWh增加到3000MWh弃电率下降非常明显但继续增加到5000MWh后收益变缓原因是弃电主要发生在中午光伏高峰和夜间风电高峰储能容量一旦覆盖这两个时段的多余电量再大的库容就没有额外作用了。同理提升抽水效率或发电效率可以直接降低储能损耗但效率改善在仿真中体现为更多电量被有效利用。还有一个容易忽略的维度是初始库容。初始库容过高抽蓄白天没有足够库容接纳光伏可能导致中午弃光初始库容过低晚高峰又放不出电来。所以日循环约束中初末库容的设置直接影响调度结果的合理性和可信度。5. 常见问题与调试技巧5.1 模型无解或报警不可行先查约束再查数据用intlinprog或Yalmip时经常会遇到“infeasible problem”或者exitflag极小的结果。遇到这种问题我的排查顺序是先检查功率平衡等式再检查库容动态和初始库容最后检查变量边界和数据单位。最常见的原因是库容动态约束和末库容约束相互矛盾。比如初始库容设得过高白天又要大量抽水导致库容顶到上限求解器发现无论如何都无法满足所有约束就会判成不可行。这时可以先把末库容约束放宽为 E_min等模型能跑通后再逐步加严。还有一种隐蔽问题是功率平衡中抽水功率前面少写了负号。抽蓄抽水时是耗能元件功率平衡右侧必须减去P_pump。有些同学把它当成发电功率加进去模型照样可能给出一个“看起来合理”的结果实际上抽蓄永远在发电库容却不会下降。最后检查结果时才发现库容曲线是条直线。5.2 求解速度慢与收敛性调优MILP是NP-hard问题变量一多求解时间可能从几秒变成几小时。如果模型规模不大但速度很慢首先要检查是否生成了大量冗余整数变量。比如把每台机组每个时段各设一个抽水和一个发电状态100台机组96时段就有19200个整数变量求解压力非常大。研究阶段可以用聚合模型代替单机模型显著降低问题规模。另外求解器参数对收敛速度影响很大。Gurobi可以设置MIPGapYalmip中通过sdpsettings(gurobi.MIPGap, 0.01)来指定1%的MIP间隙。实际调度并不需要严格全局最优1%以内的可行优解完全够用而且速度快很多。Matlab的intlinprog同样可以设置AbsoluteGapTolerance和RelativeGapTolerance。如果求解器长时间卡住不返回可能是初始解质量太差。可以先用启发式方法生成一个可行调度比如人工指定抽蓄在全天固定时段抽水、固定时段发电然后作为MIP start传入求解器。Gurobi对这类初始解很友好能明显加快分支定界过程。5.3 常见报错速查表调试过程中我整理了一张常见报错对照表很多问题都是数据结构或参数设置引起的。报错信息可能原因解决方法Index exceeds the number of array elements变量索引错位尤其是矩阵化建模时索引计算错误打印idx分段索引检查各段是否超出总变量数YALMIP: no solver found未安装求解器或路径未配置安装Gurobi/Cplex并添加路径或设置solver为intlinprogInfeasible problem数据矛盾或约束过强增加松弛变量先定位冲突约束或调整初始/末库容Solver timed out求解时间超限缩小模型规模设置MIPGap提供MIP startNaN in constraints预测数据含有NaN或Inf用isfinite检查并fillmissing处理Undefined function optimizeYalmip未正确安装运行yalmiptest检查安装还有一点要注意单位和量纲的统一。功率用MW时间用小时能量就用MWh。如果时间分辨率是15分钟dt是0.25千万别写成1。很多结果异常比如库容变化曲线严重突跳就是因为时间单位换算错了。收尾一点实际体会把这套模型从无到有跑通之后我最深的感受是数学模型写出来是一回事能在Matlab里稳定求解是另一回事。前几版代码我沉迷于把约束写得很完整结果模型规模太大求解时间动辄半个小时根本没法做敏感性分析。后来我把目标函数里的惩罚系数重新调整把机组聚合、时间分辨率改成96点速度一下子提了上来。所以如果你也在做类似的调度研究我的建议是先从一个典型日、单机组模型开始用简单的预测数据跑通流程再逐步加入机组细节、滚动修正和不确定性。数据质量永远是第一位的如果输入的风光负荷曲线本身问题百出后面所有的优化结果都会失真。希望这篇记录能帮你少走几步弯路后面如果有更深入的滚动调度或鲁棒优化实践我再继续更新。