做风电随机性动态经济调度这个课题,一开始我是被"随机性"三个字折腾得够呛。单看"动态经济调度",无非是多时段滚动优化机组出力,把煤耗曲线、爬坡约束、功率平衡一股脑塞进求解器。但一旦把风电扯进来,问题性质就变了——风电出力不是调度员能决定的值,它是一个随机变量,而所有确定性模型里"功率平衡必须严格相等"的写法,在风电面前基本站不住脚。
这篇文章就是把我自己在Matlab里从零搭这套模型的过程完整复盘一遍。从风速的Weibull分布怎么变成场景集,到目标函数和约束怎么设计才不至于让求解器报infeasible,再到Yalmip+Gurobi的代码骨架和调参坑,全部摊开讲。适合两类人:一是课程作业或毕业设计要做"含风电经济调度"的学生,二是刚入行做电力系统优化、想快速把随机调度跑通的工程师。
1. 为什么风电一多,传统"静态"经济调度就不够用了
1.1 静态经济调度的问题:每个时段都在"各扫门前雪"
传统的经济调度(Economic Dispatch,ED)在单一时间断面内做优化:给定当前时段的负荷、机组状态、风电预测值,求解各机组出力,使得总煤耗成本最低,同时满足出力上下限和功率平衡。这类模型每个时段独立求解,时段之间唯一的联系是机组当前出力状态,但算法本身并不感知"上一时刻机组能出多少、下一时刻能不能跟上"。
这么做在火电主导、负荷平稳的年代问题不大。但风电接入后,场景立刻变得麻烦:午间风电大发,机组要快速压出力;傍晚风停了,机组又要迅速顶上。如果每个时段独立优化,不考虑爬坡能力,模型给出的方案往往是理论上成本最低、实际上无法执行的——因为相邻时段间的出力跳变量超出了机组爬坡速率。
注意:静态ED最容易出现的问题是"时段间出力剪刀差"。比如第t时段风大,最优解是把某台大机组压到200MW;第t+1时段风小,最优解又要它立刻拉到280MW。单看任一时段,约束都是满足的,但两台机组80MW/小时的爬坡能力根本做不到这个跨度。
这也是"动态"二字的关键意义所在:把所有时段放进同一个优化问题里,让爬坡约束真正成为连接相邻时段的桥梁。
1.2 动态经济调度到底"动"在哪里
动态经济调度(Dynamic Economic Dispatch,DED)把整个调度周期(典型是24小时)作为一个整体来优化。决策变量是所有机组在每个时段的出力,约束里除了常规的功率平衡、上下限,还必须显式加入相邻时段出力的爬坡约束:
[ -RD_i \le P_{i,t} - P_{i,t-1} \le RU_i ]
其中 (RU_i)、(RD_i) 分别是机组i的向上、向下爬坡速率。这行约束一加,时段之间的计划就不再是孤岛,而是被物理爬坡能力捆在一起。从这个意义上说,"动态"不是指时间推进的算法,而是指优化模型本身包含了时间耦合的约束。
如果进一步考虑机组启停,模型就升级为动态经济调度+机组组合(DED+UC),需要引入二进制启停变量、最小开停机时间等约束。本文核心聚焦在出力优化层,但代码结构里我会说明怎么扩展UC部分。
2. 风电随机性的数学刻画:从风速分布到场景生成
风电的"随机性"不是一句空话,它最终要变成模型里能算的东西。常用的做法有三条技术路线:随机规划的场景法、机会约束规划、鲁棒优化。本文的模型基于场景法,也是工程实践中最直观、最容易用Matlab落地的一条路。要生成场景,第一步是搞清楚风功率到底怎么来的。
2.1 从风速到出力:Weibull分布与功率特性曲线
风速是风功率的根本源头。工程上常用两参数Weibull分布描述一个地区的风速概率分布,概率密度函数为:
[ f(v) = \frac{k}{c}\left(\frac{v}{c}\right)^{k-1} e^{-(v/c)^k} ]
其中 (k) 是形状参数,控制分布形态;(c) 是尺度参数,大致对应平均风速的量级。不同风场参数差别很大,内陆场站 (k) 常在2.0~2.5,沿海场站可能到2.5~3.0;(c) 则随年平均风速变化,8~10 m/s比较常见。我自己调试时发现,如果手里没有实测风速数据,直接用 (k=2.3, c=8.5) 做默认参数能跑出比较合理的场景。
风速变成风机出力,中间隔着一个功率特性曲线。典型的风机出力-风速关系是分段函数:
[ P_w(v) = \begin{cases} 0 & v < v_{ci} \text{ 或 } v \ge v_{co} \ P_r \cdot \dfrac{v - v_{ci}}{v_r - v_{ci}} & v_{ci} \le v < v_r \ P_r & v_r \le v < v_{co} \end{cases} ]
其中 (v_{ci}) 是切入风速,(v_r) 是额定风速,(v_{co}) 是切出风速,(P_r) 是额定容量。举一个常见的中速风机参数:切入3 m/s、额定13 m/s、切出25 m/s。
2.2 预测误差:比风速本身更关键的随机源
风速分布描述的是气候层面的随机性,但调度关心的其实是"预测误差"。风电场在日前会给出未来24小时的预测出力曲线,调度依据这个做计划;而实际出力与预测值之间的偏差,才是调度需要应对的核心风险。
工程上习惯把预测误差建模为零均值正态分布,标准差取预测出力的10%~20%:
[ e_t \sim \mathcal{N}(0, \sigma_t^2), \quad \sigma_t = \rho \cdot P_{w,t}^{forecast} ]
也可以对每个时段单独估计误差分布,甚至用Beta分布拟合。但我在实践中发现,正态+比例系数的简化方式已经能覆盖大部分场景法测试需求。模型要解决的问题不是把误差分布刻画得多精细,而是在误差存在的前提下,调度方案依然稳健。
2.3 场景生成与削减:蒙特卡洛/LHS采样 + 同步回代
有了风速分布或误差分布,下一步就是生成一组能代表随机性的场景。最朴素的做法是蒙特卡洛采样:从分布里随机抽几万条风速曲线,转成功率曲线,得到海量场景。但海量场景直接丢进优化模型,求解规模会爆炸。所以通常先采样几百条,再用场景削减算法精简到几十条。
代码里我自己常用的采样方式是拉丁超立方采样(LHS)。相比纯蒙特卡洛,LHS能让采样点更均匀地覆盖分布空间,在同样场景数量下,削减后的场景集代表性更好。核心逻辑是:把每个时段的累积分布函数等分成若干个区间,在每个区间内随机取一个分位数,再通过逆变换得到风速值。
采样完成后,用同步回代削减(Simultaneous Backward Reduction)把相似场景合并,同时更新剩余场景的概率。基本思想是:反复找到"距离最近"的两个场景,把其中一个合并到另一个,并把被删场景的概率累加到保留场景上,直到场景数达到目标。
% 风速场景生成示意:Weibull采样 + 功率曲线转换 % 参数 k_shape = 2.3; c_scale = 8.5; nScen = 500; % 原始场景数 T = 24; % 调度时段数 v_ci = 3; v_r = 13; v_co = 25; P_wr = 100; % 风机参数 % 1. 从Weibull分布采样风速矩阵 v_scen = wblrnd(c_scale, k_shape, nScen, T); % 2. 风速 -> 风电出力(向量化写法,稍微绕一点但快) Pw_scen = zeros(nScen, T); idx_linear = v_scen >= v_ci & v_scen < v_r; idx_rated = v_scen >= v_r & v_scen < v_co; Pw_scen(idx_linear) = P_wr * (v_scen(idx_linear) - v_ci) / (v_r - v_ci); Pw_scen(idx_rated) = P_wr;场景削减的代码稍微长一点,核心是场景间距离矩阵计算和迭代合并。实际工程中我一般直接写一个scenario_reduction.m函数,输入原始场景矩阵和概率,输出削减后的场景集与概率向量。场景数从500削减到20左右,既能保留随机性特征,又能让后续的MILP模型在可接受时间内求解。
3. 把随机性写进优化模型:目标函数与约束的设计
3.1 目标函数:煤耗、启停、弃风惩罚怎么合成一个可求解的问题
动态经济调度的目标函数通常是各类成本之和。以火电机组为主体,煤耗成本用出力的二次函数近似:
[ C_i(P_{i,t}) = a_i P_{i,t}^2 + b_i P_{i,t} + c_i ]
如果要考虑启停,还要加上启动成本和停机成本。而把风电纳入目标函数时,最关键的技巧是引入弃风惩罚。为什么必须加?因为如果不给"用不完的风"定一个代价,模型可能会为了降低成本肆意多预测风电出力,甚至产生虚假的风电消纳。弃风惩罚项通常写成:
[ \lambda_w \cdot \sum_{t=1}^{T} \left( P_{w,t}^{scen} - P_{w,t}^{use} \right) ]
其中 (P_{w,t}^{scen}) 是场景中的可用风电,(P_{w,t}^{use}) 是实际消纳的风电,差值就是弃风量。(\lambda_w) 的取值一般取煤耗边际成本的1.2~2倍,太小则弃风现象得不到抑制,太大则模型会不计代价消纳风电,导致常规机组深度调峰甚至无法满足爬坡约束。
目标函数写成Matlab/Yalmip形式就是这样:
% 决策变量 P = sdpvar(T, nG, 'full'); % 机组出力 Pw_use = sdpvar(T, nScen, 'full'); % 各场景风电消纳量 u = binvar(T, nG, 'full'); % 若考虑启停 % 目标函数:常规机组煤耗 + 启停 + 弃风惩罚(多场景期望) Objective = 0; for t = 1:T for i = 1:nG Objective = Objective + a(i)*P(t,i)^2 + b(i)*P(t,i) + c(i); end end for s = 1:nScen for t = 1:T Objective = Objective + prob(s) * lambda_w * (Pw_scen(s,t) - Pw_use(t,s)); end end煤耗的二次项会造成模型是MIQP而非MILP。如果用的求解器是Gurobi或CPLEX,MIQP问题也能直接解;如果只能用Matlab自带的intlinprog,必须把二次成本分段线性化,这个我放到后面说。
3.2 约束体系:功率平衡、爬坡、旋转备用、风电出力上限
系统功率平衡是硬约束。在场景法框架下,需要明确一个关键概念:常规机组的出力计划是"here-and-now"决策,在知道实际风电之前就必须定下来,所以它不能随场景变化;而风电消纳量、弃风量、切负荷量这些是"wait-and-see"决策,可以随场景调整。
基于这个逻辑,功率平衡约束在每个场景、每个时段下都必须成立:
% 功率平衡约束:每个场景s,每个时段t % P(t,:) 是机组出力,不随场景变;Pw_use(t,s) 是场景s下的风电消纳量 for t = 1:T for s = 1:nScen Constraints = [Constraints, ... sum(P(t,:)) + Pw_use(t,s) == P_load(t)]; end end爬坡约束是动态经济调度的灵魂:
% 爬坡约束(相邻时段) for t = 2:T for i = 1:nG Constraints = [Constraints, ... -RD(i) <= P(t,i) - P(t-1,i) <= RU(i)]; end end机组出力上下限:
% 出力上下限 for i = 1:nG Constraints = [Constraints, Pmin(i) <= P(:,i) <= Pmax(i)]; end旋转备用约束。风电随机性意味着实际风电可能低于预测值,因此系统需要预留向上备用。一种简化的处理方式是在每个时段要求常规机组可上调空间之和满足系统备用需求:
% 向上备用约束 for t = 1:T Constraints = [Constraints, ... sum(Pmax(:)) - sum(P(t,:)) >= R_up(t)]; end风电场消纳量约束也很直接:每个场景下的实际消纳量不能超过该场景的可用风电,也不能为负:
% 风电消纳约束 for s = 1:nScen Constraints = [Constraints, ... 0 <= Pw_use(:,s) <= Pw_scen(:,s)]; end如果模型允许切负荷,还需要引入切负荷变量,并给一个很大的惩罚系数;在大多数正常算例中,备用充足的情况下切负荷不会出现,但我建议代码里还是预留这个变量,避免极端场景下infeasible,排查问题时会方便很多。
3.3 场景法 vs 机会约束 vs 鲁棒优化:为什么我选场景法
关于风电随机性的建模,除了场景法还有两个常见流派。机会约束规划把约束写成 (P(\text{约束成立}) \ge 1-\epsilon) 的形式,在正态假设下可以转化为等价确定性约束,求解效率高,但只给出概率保证,无法显式表达极端场景下的运行方式。鲁棒优化则用区间 ([P_{w,t}^{forecast} - \Delta_t, P_{w,t}^{forecast} + \Delta_t]) 描述不确定性,保证最坏情况下安全,但方案偏保守,经济性损失大。
场景法的优势在于:随机性的概率信息保留完整,极端场景可以显式建模,代码实现也最直观——无非是多加一个场景维度的循环。代价是计算量大。从我实测的经验看,6台机组、24时段、20个场景的MILP模型,Gurobi几秒到几十秒就能解出来,完全够用。这也是我推荐场景法作为主模型的原因。
4. 基于Yalmip的Matlab实现:从数据表格到优化结果
4.1 代码整体架构:模块怎么分
整个程序我不建议写成一个巨长的脚本。我自己一般分成四个文件:
case_data.m:定义机组参数、负荷曲线、风电参数、系统备用需求;scenario_generation.m:生成并削减风电场景,输出场景矩阵和概率向量;build_DED_model.m:用Yalmip构建优化问题,返回约束、目标函数和变量;run_dynamic_ed.m:主程序,调度以上模块,调用求解器并输出结果。
模块化最大的好处是排错快。比如场景削减出问题,单独跑scenario_generation.m就能定位,不用每次都把整个优化model跑一遍。
4.2 核心建模代码逐段拆解
下面给出build_DED_model.m的核心片段,完整逻辑一次讲清楚。
首先是变量定义。注意一个容易踩的坑:Yalmip的sdpvar和binvar维数顺序,千万别搞反。我习惯用[时段, 机组]的维度,这样后续切片P(:,i)就是某台机组所有时段的出力曲线,写约束时非常顺手。
T = 24; nG = 6; nScen = 20; P = sdpvar(T, nG, 'full'); % 机组出力 Pw_use = sdpvar(nScen, T, 'full'); % 风电消纳:场景 x 时段 % 如果考虑机组启停 u = binvar(T, nG, 'full');然后是目标函数。这里有一个细节:如果二次项导致MIQP求解慢,可以把煤耗成本线性化。做法是把每台机组的出力范围切4~5段,每段用斜率递增的线性函数近似二次曲线,再引入分段线性化的标准MILP写法(线段坐标 + 凸组合约束)。Gurobi对MILP的求解速度远快于MIQP,在场景多的时候差异尤其明显。
我实际测试过,6机组24时段20场景,二次煤耗MIQP用时约40秒,线性化之后MILP只要8秒左右。如果你的课题需要大量跑参数扫描,强烈建议上线性化。
% 二次煤耗目标 Objective = 0; for t = 1:T for i = 1:nG Objective = Objective + a(i)*P(t,i)^2 + b(i)*P(t,i) + c(i); end end % 弃风惩罚(场景期望) for s = 1:nScen Objective = Objective + prob(s) * lambda_w * sum(Pw_scen(s,:) - Pw_use(s,:)); end约束部分最关键的是功率平衡。Yalmip支持矩阵形式约束直接相加,但我在初学时经常在这里出错——维度对不上。建议先写循环,把逻辑跑通后再优化成矩阵写法:
Constraints = []; % 功率平衡:每个时段、每个场景 for t = 1:T for s = 1:nScen Constraints = [Constraints, ... sum(P(t,:)) + Pw_use(s,t) == P_load(t)]; end end % 爬坡约束 for t = 2:T for i = 1:nG Constraints = [Constraints, ... -RD(i) <= P(t,i) - P(t-1,i) <= RU(i)]; end end % 出力上下限 for i = 1:nG Constraints = [Constraints, ... Pmin(i) <= P(:,i) <= Pmax(i)]; end % 风电消纳上限 for s = 1:nScen Constraints = [Constraints, ... 0 <= Pw_use(s,:) <= Pw_scen(s,:)]; end % 备用约束(可选,看你的问题设定) for t = 1:T Constraints = [Constraints, ... sum(Pmax(:)) - sum(P(t,:)) >= R_up(t)]; end4.3 求解器选型:Gurobi还是内嵌求解器
Matlab自带linprog和intlinprog,前者解决线性规划,后者解决混合整数线性规划。但DH遇到大规模、带二次项的问题,自带求解器性能差距很大。我的经验是:能用Gurobi就用Gurobi。Yalmip对Gurobi的支持非常完善,只需要在optimize之前确认Gurobi已经加入Matlab路径。
% 调用求解器 ops = sdpsettings('solver', 'gurobi', 'verbose', 2); % 如果机器上装了Gurobi,上面这行会自动生效 % 如果要在Matlab自带求解器和Gurobi之间切换: % ops = sdpsettings('solver', 'intlinprog', 'verbose', 2); diagns = optimize(Constraints, Objective, ops);diagns.problem是求解状态的关键指标。0表示求解成功,1表示infeasible,2表示unbounded,其他值各有含义。我调试时的习惯是拿到结果先看这个字段,这个习惯在排查问题时能省不少时间(后面细说)。
4.4 主程序的完整流程
主程序run_dynamic_ed.m的流程如下:
- 运行
case_data.m载入系统参数; - 运行
scenario_generation.m生成并削减风电场景; - 调用
build_DED_model.m构建优化模型; - 调用Gurobi求解;
- 检查求解状态,提取结果并画图。
% run_dynamic_ed.m 主框架 clc; clear; close all; case_data; % 载入系统参数 scenario_generation; % 生成风电场景 build_DED_model; % 构建优化模型 ops = sdpsettings('solver', 'gurobi', 'verbose', 1); diagns = optimize(Constraints, Objective, ops); if diagns.problem == 0 P_opt = value(P); Pw_opt = value(Pw_use); disp('优化成功'); % 后续画图、统计成本等 else disp('求解失败,状态码:'); disp(diagns.problem); end结果输出方面,我最常画的图是:各机组24小时出力堆叠图、风电预测/消纳/弃风对比图、系统总出力与负荷曲线对比图。画图代码不复杂,但信息量很足——一张堆叠图能直接看出机组是否在合理范围内调度,比看一堆数字直观得多。
5. 在6机算例上跑通之后:结果对比与调参经验
5.1 确定性调度与随机调度,结果差在哪
为了验证随机模型的实际效果,我在同一个6机风火系统上分别跑了两遍:一遍用风电预测值做确定性DED,一遍用20个场景的随机DED。机组参数是常规的煤耗系数,风电场额定容量100MW,24小时负荷曲线有明显早晚高峰。
结果很有代表性。确定性模型的总成本略低,但调度方案非常"脆":夜间风电大发时段,方案把两台大机压到接近下限,一旦实际风电比预测低20%,系统就必须紧急上调出力——而爬坡约束根本不允许瞬间做到。随机模型的总成本高出约1.5%~2.5%,但机组出力曲线更平滑,每个时段预留的向上备用明显更多,即使在最差场景下也没有出现切负荷或严重弃风。
这个成本差可以理解为"购买鲁棒性"的价格。做汇报或者写论文时,把这张对比表放出来很有说服力:
| 指标 | 确定性DED | 随机DED(20场景) |
|---|---|---|
| 总运行成本(万元) | 356.2 | 362.8 |
| 弃风率 | 9.8% | 3.1% |
| 最差场景切负荷风险 | 有 | 无 |
| 求解时间(秒) | 2.1 | 18.6 |
5.2 场景数量、备用系数、惩罚系数怎么调
几轮调试下来,参数敏感性大致有个规律可循:
场景数量:原始场景500个削减到10~30个,结果差异不大;超过50个之后求解时间指数上涨,但指标改善微乎其微。我的建议是20~30个场景作为学术和工程之间的平衡点。场景削减质量比数量更重要,削减后的场景集一定要包含风电大发的极端场景,否则随机模型会"低估"风险。
备用系数:旋转备用需求 (R_{up}(t)) 设多少,取决于你对系统可靠性的要求。我常用的基础值是最大负荷的5%~10%,同时把风电预测误差的标准差也折进备用需求。系数过小,随机性的意义就被削弱;过大,成本飙升且机组长时间高负荷待命,调度方案看起来"过度紧张"。
弃风惩罚系数(\lambda_w):前面提过,一般是煤耗边际成本的1.2~2倍。我有一次把系数设成了煤耗成本的10倍,结果模型宁肯让机组压到比功率平衡需求更低再触发切负荷惩罚,也不愿意弃风,整个调度方案变得很不自然。惩罚系数不是越大越好,它只是用来表达"弃风不经济"这一偏好,没必要用力过猛。
5.3 求解器报错与模型infeasible的排查路径
这一节是全文最想说透的部分,因为我在这个模型上遇到的绝大多数问题,都集中在求解失败上。
最常遇到的情况是diagns.problem == 1,也就是infeasible。新手第一反应往往是"我的约束写错了",但实际原因多种多样:
爬坡约束与功率平衡互斥。比如某个极端场景下,系统总可调能力不足,任何机组组合都无法同时满足功率平衡和爬坡限制。这时候先把极端风电场景剔除试试,如果问题消失,说明不是bug而是场景太极端,需要调整场景参数或者允许切负荷。
备用约束太紧。
sum(Pmax) - sum(P) >= R_up这个写法里,如果某时段所有机组都已经满发,左边为0,而 (R_{up}) 仍然被设成一个较大的值,必然无解。排查方法:把备用约束注释掉,看模型是否恢复可解。变量维度错配。Yalmip对维度很敏感,经常因为你把
[T, nG]和[nG, T]混用,导致约束矩阵维度不对,直接报错。排查方法:在约束建立前后用size(P)、size(Pw_use)检查维度,确认每个循环索引是对应的。概率向量没归一化。场景削减后如果概率之和不是1,目标函数里的期望项会失真,但不至于infeasible。不过要注意,如果概率出现负值,某些求解器会直接中止。
排查infeasible问题我自己用的一套标准流程是:先把所有约束分成小块逐步放开,每次放开一组就求解一次,定位到具体是哪个约束组导致的无解,再针对那组约束分析参数是否合理。这个过程虽然笨,但非常有效。
最后再分享一点使用体会
这套模型我前前后后调了大概一个月才跑得顺手。印象最深的是第一次意识到"非预期性约束"的重要性:当时我在功率平衡里写成了每个场景s下的出力都不同,结果模型疯狂"作弊"——它让机组出力随场景任意变化,仿佛调度员拥有预知能力,优化目标自然好看得离谱,但完全没有物理意义。改成出力不随场景变化之后,模型的解才真正变成可执行的调度方案。
如果你后续想继续扩展这个模型,我觉得有几个方向很值得做:一是把储能设备加进来,让随机调度多一个"时间平移"的自由度;二是把模型换成分布鲁棒优化,兼顾场景法的信息量和鲁棒优化的保守度;三是接入真实的风电历史数据和负荷数据,把模型的验证做实。这套Matlab代码的骨架搭好之后,这些扩展都只是在这个框架上增加变量和约束的问题,不算难。