做风电随机性动态经济调度这个课题一开始我是被随机性三个字折腾得够呛。单看动态经济调度无非是多时段滚动优化机组出力把煤耗曲线、爬坡约束、功率平衡一股脑塞进求解器。但一旦把风电扯进来问题性质就变了——风电出力不是调度员能决定的值它是一个随机变量而所有确定性模型里功率平衡必须严格相等的写法在风电面前基本站不住脚。这篇文章就是把我自己在Matlab里从零搭这套模型的过程完整复盘一遍。从风速的Weibull分布怎么变成场景集到目标函数和约束怎么设计才不至于让求解器报infeasible再到YalmipGurobi的代码骨架和调参坑全部摊开讲。适合两类人一是课程作业或毕业设计要做含风电经济调度的学生二是刚入行做电力系统优化、想快速把随机调度跑通的工程师。1. 为什么风电一多传统静态经济调度就不够用了1.1 静态经济调度的问题每个时段都在各扫门前雪传统的经济调度Economic DispatchED在单一时间断面内做优化给定当前时段的负荷、机组状态、风电预测值求解各机组出力使得总煤耗成本最低同时满足出力上下限和功率平衡。这类模型每个时段独立求解时段之间唯一的联系是机组当前出力状态但算法本身并不感知上一时刻机组能出多少、下一时刻能不能跟上。这么做在火电主导、负荷平稳的年代问题不大。但风电接入后场景立刻变得麻烦午间风电大发机组要快速压出力傍晚风停了机组又要迅速顶上。如果每个时段独立优化不考虑爬坡能力模型给出的方案往往是理论上成本最低、实际上无法执行的——因为相邻时段间的出力跳变量超出了机组爬坡速率。注意静态ED最容易出现的问题是时段间出力剪刀差。比如第t时段风大最优解是把某台大机组压到200MW第t1时段风小最优解又要它立刻拉到280MW。单看任一时段约束都是满足的但两台机组80MW/小时的爬坡能力根本做不到这个跨度。这也是动态二字的关键意义所在把所有时段放进同一个优化问题里让爬坡约束真正成为连接相邻时段的桥梁。1.2 动态经济调度到底动在哪里动态经济调度Dynamic Economic DispatchDED把整个调度周期典型是24小时作为一个整体来优化。决策变量是所有机组在每个时段的出力约束里除了常规的功率平衡、上下限还必须显式加入相邻时段出力的爬坡约束[ -RD_i \le P_{i,t} - P_{i,t-1} \le RU_i ]其中 (RU_i)、(RD_i) 分别是机组i的向上、向下爬坡速率。这行约束一加时段之间的计划就不再是孤岛而是被物理爬坡能力捆在一起。从这个意义上说动态不是指时间推进的算法而是指优化模型本身包含了时间耦合的约束。如果进一步考虑机组启停模型就升级为动态经济调度机组组合DEDUC需要引入二进制启停变量、最小开停机时间等约束。本文核心聚焦在出力优化层但代码结构里我会说明怎么扩展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比较常见。我自己调试时发现如果手里没有实测风速数据直接用 (k2.3, c8.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_{t1}^{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或CPLEXMIQP问题也能直接解如果只能用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表示infeasible2表示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。机组参数是常规的煤耗系数风电场额定容量100MW24小时负荷曲线有明显早晚高峰。结果很有代表性。确定性模型的总成本略低但调度方案非常脆夜间风电大发时段方案把两台大机压到接近下限一旦实际风电比预测低20%系统就必须紧急上调出力——而爬坡约束根本不允许瞬间做到。随机模型的总成本高出约1.5%~2.5%但机组出力曲线更平滑每个时段预留的向上备用明显更多即使在最差场景下也没有出现切负荷或严重弃风。这个成本差可以理解为购买鲁棒性的价格。做汇报或者写论文时把这张对比表放出来很有说服力指标确定性DED随机DED20场景总运行成本万元356.2362.8弃风率9.8%3.1%最差场景切负荷风险有无求解时间秒2.118.65.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代码的骨架搭好之后这些扩展都只是在这个框架上增加变量和约束的问题不算难。