风电随机性动态经济调度：Matlab+Yalmip建模与场景削减实战

风电随机性的动态经济调度，这个题目我前后玩了有一阵子。做电力系统优化的人应该都有感触，传统的经济调度模型，大多基于确定性负荷预测，给一组固定的机组出力。但风电一旦接入，情况就完全不一样了——风速本身是个随机过程，出力跟着波动，如果在模型里不处理这层随机性，算出来的调度方案轻则保守、重则不可行，到了实际运行里可能连负荷都平衡不了。

这篇文章把我用Matlab搭建风电随机性动态经济调度模型的完整过程记录下来，包括风电随机性怎么建模、动态约束怎么处理、旋转备用怎么跟置信水平挂钩，以及完整的Yalmip代码实现和调试经验。不管你是正在做毕业设计，还是刚接触随机优化想找个能跑的框架，这篇都能给你一条清晰可执行的路径。

1. 模型问题定义与整体思路

1.1 风电随机性到底在随机什么

风电出力的随机性，根源在于风速的随机性。天然风速的分布通常用两参数Weibull分布描述，尺度参数和形状参数直接决定风速分布的形态。但风电出力和风速之间是非线性关系，可不是简单的线性映射，风速太小发不出电，风速太大为了保护机组反而要切机，中间段才近似线性。这种非线性会放大风速的随机性，导致风电出力呈现出一种偏态、厚尾的分布特征。

风电随机性对调度的影响，不是简单的"出力预测不准"这么轻描淡写。一个很直接的问题是：如果实际风电出力比预测值低很多，常规机组的爬坡能不能跟上？如果实际风电出力偏高，系统能不能及时压减火电出力？这些都是动态约束层面的考量，需要模型显式地把风电随机性纳入优化范畴。

很多人第一次做这个题目，容易把风电处理成确定性场景或仅用期望值，然后用一个固定的旋转备用量来兜底。这么做不是不行，但本质上没有回答一个核心问题：风电出力在给定置信水平下的波动区间到底是多少？以及这个波动区间对应的备用需求是多少？正因为如此，我在建模时选择了随机场景法，把风电出力用多场景模拟，同时对旋转备用约束加上置信水平约束，让模型自己去寻找一个兼顾经济性与可靠性的调度方案。

1.2 动态经济调度和经济调度的区别

传统的静态经济调度，只看一个时间断面的负荷需求和发电成本，做出该时段的机组出力安排。这种静态方案在火电时代还算够用，因为火电机组可以近似认为在较短时间内调整到位，各时段间的耦合不强。风电接入后情况不同了，风电机组的出力波动性强，而且具有明显的时序相关性，上一时段的风电出力会影响下一时段的机组爬坡空间。

动态经济调度是在静态经济调度的基础上，引入跨时段约束。通俗地讲，就是机组不仅要满足本时段的负荷平衡，还要满足相邻时段之间的爬坡速率约束、最小启停时间约束、燃料约束等。这种跨时段耦合让问题的维度一下子上去，也让优化结果更贴近真实运行情况。

两者的本质区别在于，动态经济调度的决策变量是时序关联的，前一时段的最优解不一定在下一时段仍然可行，需要从整个调度周期全局优化。所以在Matlab实现里，模型规模会比静态调度大很多，需要把决策变量向量化构造，减少循环，提升求解效率。

1.3 建模思路：为什么用场景法而不是解析法

把风电随机性纳入优化模型，主流路线大约有这几类：机会约束规划、鲁棒优化、随机场景规划。三者各有适用场景，但面试里最常遇到、写论文最容易讲清楚的还是随机场景规划，也是我最推荐上手的方式。

场景法的核心思路很简单：用一组离散的风电出力场景来近似风电出力的连续概率分布。比如通过蒙特卡洛采样生成1000个风速样本，映射成风电出力，再用场景削减技术保留10个有代表性的场景。每个场景有对应的概率，优化目标变成所有场景下期望成本最小化，约束条件要求每个场景下系统都满足运行约束。

这个方法比机会约束规划更容易在Matlab中实现，而且不依赖对特定概率分布的假设，对风电预测误差甚至实际风速分布的拟合都适用。鲁棒优化虽然更保守，但对不确定集合的参数设置比较敏感，对新手来说调参成本较高。场景法最大的优势是直观，生成场景、削减场景、跑优化、看结果，每一步都能可视化，方便理解和调试。

2. 风电随机性建模与场景生成

2.1 风速概率模型与实际风电出力计算

风速模型我选用两参数Weibull分布，这是业内最常见的风速概率描述方式。概率密度函数是：

[ f(v) = \frac{k}{c}\left(\frac{v}{c}\right)^{k-1} \exp\left[-\left(\frac{v}{c}\right)^k\right] ]

其中 (k) 是形状参数，(c) 是尺度参数。风电场的实测风速拟合通常得到 (k) 在1.8到2.5之间，(c) 在6到10之间。但如果手头没有实测数据，也可以设定一个合理值然后生成符合这个分布的样本。

风速转化为风电出力，用的是风电机组的标准功率曲线。工程上常用三段式近似：

[ P_w(v) = \begin{cases} 0, & v < v_{ci} \text{ 或 } v > v_{co} \ \frac{v^3 - v_{ci}^3}{v_r^3 - v_{ci}^3}P_r, & v_{ci} \le v < v_r \ P_r, & v_r \le v \le v_{co} \end{cases} ]

这里的 (v_{ci})、(v_r)、(v_{co}) 分别是切入风速、额定风速和切出风速。常见参数是切入3m/s、额定12m/s、切出25m/s。公式并不复杂，但要注意在程序实现时对分段条件做向量化处理，别写成for循环逐点判断，否则速度会很难看。

2.2 场景生成：蒙特卡洛采样实操

场景生成的过程，说白了就是用随机数模拟风速的不确定性，然后映射成风电出力场景。这里有个关键点，风电出力在不同时段之间不是独立的，相邻时段的风速具有相关性，如果不处理这个相关性直接采样，生成的场景会过度震荡，最终优化的结果也会偏乐观或者偏悲观。

处理时段相关性的一个实用做法，是先对风速序列建立一阶自回归模型，也就是AR(1)模型：

[ v_t = \mu + \rho(v_{t-1} - \mu) + \varepsilon_t ]

(\rho) 根据实际风速序列的自相关系数来定，通常在0.8到0.9之间。这样采样时，先采第一时段的风速，后面每个时段的风速都基于前一时刻的值生成，自然带上时序相关性。在Matlab里生成风电出力场景序列的代码可以这样写：

% 风电随机场景生成 - 基于AR(1)风速模型 % 参数设置 k = 2.2; % Weibull形状参数 c = 8.5; % Weibull尺度参数 rho = 0.85; % 风速自相关系数 H = 24; % 调度时段数 N_scenarios = 500; % 初始场景数 % 风电机组参数 v_ci = 3; v_r = 12; v_co = 25; Pr = 300; % 额定功率MW % 生成风速场景矩阵 (N_scenarios x H) v = zeros(N_scenarios, H); for s = 1:N_scenarios % 第一时段直接按Weibull采样 v(s,1) = wblrnd(c, k); % 后续时段基于AR(1)递推 for t = 2:H eta = wblrnd(c, k) - c*gamma(1+1/k); % 均值归零的随机项 v(s,t) = c * (1-rho) + rho * v(s,t-1) + eta; if v(s,t) < 0, v(s,t) = 0; end end end % 风速转出力 P_wind = zeros(N_scenarios, H); for s = 1:N_scenarios for t = 1:H vt = v(s,t); if vt < v_ci || vt > v_co P_wind(s,t) = 0; elseif vt >= v_r P_wind(s,t) = Pr; else P_wind(s,t) = Pr * (vt^3 - v_ci^3) / (v_r^3 - v_ci^3); end end end

上面的代码胜在逻辑清晰、没绕弯，新手对照公式就能看懂。它把风速生成和出力映射分成两个阶段，调试起来也很方便，哪个环节出问题一目了然。

2.3 场景削减：同步回代消除法

初始场景数越大，随机性刻画越精确，但模型求解规模也越大。风电随机性动态经济调度场景数量翻倍，决策变量和约束条件也跟着接近倍增，求解时间可能从几十秒激增到几十分钟。为了平衡精度和计算效率，需要对场景削减。

场景削减的经典算法是同步回代消除法。思路是：迭代删除对场景集合概率分布影响最小的场景，同时把这个场景的概率累加到离它最近的场景上。这个算法的Matlab实现不复杂，核心代码如下：

% 场景削减 - 同步回代消除法 function [reduced_scen, reduced_prob] = scenario_reduction(scenarios, probs, keep_num) % scenarios: 每个场景是一个行向量，矩阵维度为 场景数 x 时段数 % probs: 每个场景的概率，列向量 current_scen = scenarios; current_prob = probs; n_scen = size(current_scen, 1); while n_scen > keep_num % 计算场景两两之间的距离（欧氏距离） D = zeros(n_scen, n_scen); for i = 1:n_scen for j = 1:n_scen if i ~= j D(i,j) = norm(current_scen(i,:) - current_scen(j,:)); end end end % 找到要删除的场景，使得概率加权距离最小 weighted_dist = zeros(n_scen, 1); for i = 1:n_scen % 找到距离场景i最近的其他场景 [min_dist, ~] = min(D(i,:), [], 2); weighted_dist(i) = current_prob(i) * min_dist; end [~, idx_del] = min(weighted_dist); % 找到与被删场景最近的场景，把概率转移过去 [~, idx_near] = min(D(idx_del,:), [], 2); current_prob(idx_near) = current_prob(idx_near) + current_prob(idx_del); % 删除该场景 current_scen(idx_del, :) = []; current_prob(idx_del) = []; n_scen = n_scen - 1; end reduced_scen = current_scen; reduced_prob = current_prob; end

这个算法的for循环嵌套虽然效率一般，但对场景数据量（几百个）来说完全够用。如果追求更快的速度，可以在距离计算部分使用pdist2函数优化。削减后的场景数一般取10到20个比较合适，太少丢随机性，太多模型求解慢。

3. 动态经济调度数学模型

3.1 目标函数与约束条件

动态经济调度模型的目标函数是调度周期内所有时段的总运行成本最小化，包括火电机组的燃料成本和风电的运维成本。风电的运维成本很低，通常只有火电成本的零头，但还是要计入模型，否则可能出现风电占比虚高的调度方案。目标函数数学形式如下：

[ \min \sum_{t=1}^{T} \sum_{i=1}^{N_g} \left[ a_i P_{i,t}^2 + b_i P_{i,t} + c_i \right] + \sum_{t=1}^{T} c_w P_{w,t} ]

(a_i, b_i, c_i) 是火电机组煤耗特性系数，(c_w) 是风电单位运维成本。火电成本用二次函数描述，在Yalmip里用二次目标函数求解效率很不错。

约束条件包含几个方面。首先是系统功率平衡约束：

[ \sum_{i=1}^{N_g} P_{i,t} + P_{w,t}^{scen} = D_t, \quad \forall t ]

这里 (P_{w,t}^{scen}) 是每个场景下风电的预测出力。然后是机组出力上下限约束，以及动态经济调度关键的爬坡约束：

[ -R_{i}^{down} \le P_{i,t} - P_{i,t-1} \le R_{i}^{up} ]

爬坡约束是逐时段耦合的，直接导致模型求解难度上升，也是动态模型和静态模型最本质的区别所在。

3.2 旋转备用约束的随机处理

风电随机性带来的最大挑战在旋转备用上。传统确定性模型通常设置备用容量不小于最大负荷的某个比例或最大单机容量，这个做法对高比例风电系统来说过于粗糙。

我采用的方案是把旋转备用约束和置信水平结合。考虑风电预测误差的概率分布，备用约束写成：

[ \sum_{i=1}^{N_g} R_{i,t}^{up} \ge D_t \times \epsilon + \text{VaR}{\alpha}(P{w,t}) ]

其中 (\text{VaR}{\alpha}(P{w,t})) 是风电出力的在险价值，表示在置信水平 (\alpha) 下的最大预测误差。实际计算中，可以从场景集合中统计不同置信度下的风电出力分位数。用Matlab的prctile函数就能算出，比如置信水平95%对应的风电出力下分位数：

% 计算风电出力在置信水平下的分位数 alpha = 0.95; P_wind_low = prctile(P_wind, (1-alpha)*100, 1); % 备用需求 = 负荷比例备用 + 风电不确定性备用 reserve_req = 0.05 * D + (mean(P_wind,1) - P_wind_low);

这段代码的逻辑是，在95%置信水平下，风电出力最低可能低到 (P_{wind_low})，所以需要额外备用来弥补期望出力与该低分位数之间的差值。这种做法比拍脑袋定备用比例要科学得多，而且思路在论文里也容易解释得通。

3.3 模型整体架构

模型整体架构分三层。最底层是输入数据层，包括火电机组参数表、负荷曲线、风速分布参数、风电场参数。中间层是场景处理层，负责生成和削减风电场景。最高层是优化调度层，把场景数据转化为优化模型的参数，调用求解器求解。

这种分层架构的好处是各层可以独立修改。比如想换一种风速模型，只需修改场景处理层的代码，优化调度层的代码不用动。想增加机组数，只需修改数据层的参数表，其他部分自动适应。

4. Matlab代码实现细节

4.1 数据准备与参数设置

我用一个6机测试系统做算例，这是电力系统经济调度领域很经典的标准测试系统。6台机组参数覆盖不同类型，有燃煤、燃气，机组容量从50MW到200MW不等，成本系数也差异明显，能真实反映不同机组间的调度博弈。机组参数表如下：

负荷曲线采用典型日负荷数据，峰值负荷650MW，谷值负荷480MW，用一条平滑的24时段曲线模拟。

4.2 核心求解代码（Yalmip建模）

建模用Yalmip工具箱，求解器选用Cplex或Gurobi都可以。Yalmip的优势在于建模语言贴近数学表达式，写出来的代码可读性强，而且Cplex这类商业求解器和开源的求解器之间的切换，只需要改一行代码就行。核心建模代码：

% 动态经济调度 - Yalmip建模求解 % 决策变量 P = sdpvar(6, H, 'full'); % 火电机组出力 R_up = sdpvar(6, H, 'full'); % 上旋转备用 R_down = sdpvar(6, H, 'full'); % 下旋转备用 % 目标函数：期望成本最小化 objective = 0; for t = 1:H for i = 1:6 objective = objective + a(i)*P(i,t)^2 + b(i)*P(i,t) + c(i); end end % 加入风电成本（取场景期望） objective = objective + c_w * mean(P_wind_reduced, 1) * ones(H,1); % 约束条件 constraints = []; for t = 1:H % 功率平衡约束（对每个场景） for s = 1:N_keep constraints = [constraints, sum(P(:,t)) + P_wind_reduced(s,t) == D(t)]; end % 机组出力上下限 constraints = [constraints, P(:,t) >= P_min, P(:,t) <= P_max]; % 备用容量约束 constraints = [constraints, R_up(:,t) >= 0, R_up(:,t) <= P_max - P(:,t)]; constraints = [constraints, R_down(:,t) >= 0, R_down(:,t) <= P(:,t) - P_min]; constraints = [constraints, sum(R_up(:,t)) >= reserve_req(t)]; end % 爬坡约束（跨时段耦合） for t = 2:H constraints = [constraints, P(:,t) - P(:,t-1) <= R_up_rate]; constraints = [constraints, P(:,t-1) - P(:,t) <= R_down_rate]; end % 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); optimize(constraints, objective, ops);

实际运行时要特别注意一个细节：功率平衡约束应该只对基准场景或期望场景严格等值成立，而对每个单独场景严格等值，会导致不同场景间的机组出力互相打架，模型从数学上看甚至可能是不可行的。我的做法是把功率平衡约束的等号改为期望场景下的平衡约束，而各场景的差异通过备用约束来吸收，这样模型就有解了。

4.3 结果可视化与分析

求解完成后，把机组出力计划画出来，最容易看出调度方案的合理性。用这么一段代码做可视化：

% 结果可视化 figure; bar(P', 'stacked'); hold on; plot(D, 'r-', 'LineWidth', 2); plot(mean(P_wind_reduced,1), 'g-', 'LineWidth', 2); legend('G1','G2','G3','G4','G5','G6','负荷','风电出力'); xlabel('时段/h'); ylabel('功率/MW'); title('动态经济调度结果'); % 各机组出力曲线 figure; for i = 1:6 subplot(3,2,i); plot(P(i,:), 'LineWidth', 1.5); title(['机组 G', num2str(i)]); xlabel('时段/h'); ylabel('出力/MW'); ylim([0, max(P_max(i)*1.2, 1)]); end

从可视化结果能直观看出：高峰时段大机组（G1、G2）出力增加，低谷时段小机组压荷，风大的时段火电出力相应减小。如果模型结果出现某台机组出力曲线抖得厉害，多半是爬坡约束或备用约束设置不合理，需要回头检查参数的合理性。

5. 算例分析与结果讨论

5.1 不同置信水平下的调度结果对比

旋转备用约束的置信水平从85%调到99%，对调度结果的影响非常显著。我跑了一组对比实验，结果如下：

置信水平从85%提高到99%，总成本增加了约5.9%，但最坏情况下的备用容量从128.5MW增加到了286.4MW，增幅超过一倍。这说明，追求更高的供电可靠性，代价不只是备用的线性增加，还会因为备用抬高了机组的运行区间，导致煤耗成本非线性上涨。从实例里基本能得出一个结论：95%置信水平是一个比较划算的折中点，再往上走边际成本就明显偏高了。

5.2 风电渗透率对调度成本的影响

把风电场容量从100MW逐步提高到500MW，观察系统总成本的变化。这组实验用来回答一个关键问题：风电多了，系统是不是真的省钱了？

风电容量增加后，总成本确实在下降，但降幅越来越小。原因有两方面：一是风电出力小时段，火电机组已经压到最低技术出力，没法再让出更多空间；二是备用需求跟着风电容量增加而增加，这部分成本抵消了一部分燃料成本的节省。当风电容量到500MW时，弃风率接近10%，说明系统接纳风电的能力已经接近极限。

5.3 模型收敛性与求解性能

在Matlab里跑这个模型，我观察到的求解时间分布大致是：场景生成和削减耗时不到1秒，Yalmip建模耗时1秒左右，Cplex求解耗时3到5秒。整体来说，对一个24时段、6机组、10个场景的模型，总耗时在10秒内，完全能满足学习和研究用途。

如果要把规模扩大到几十台机组、数百个节点，建议在建模时使用稀疏矩阵构造约束，或者先对模型做预处理，看看能不能合并冗余约束。另外一个实用技巧是，为爬坡约束设置一个好的初始解，可以显著缩短求解时间。

6. 常见问题与调试经验

6.1 求解失败：不可行问题排查思路

遇到infeasible problem，基本上就是约束条件互相矛盾。排查路径按顺序来：先单独检查功率平衡约束，把风电场景全取期望值，看能不能求解；再逐步加入爬坡约束和备用约束，哪一层加了以后模型不可行，问题就出在哪一层。

我踩过一个很典型的坑：备用约束里，我把旋转备用上限和机组出力上限绑在一起，同时爬坡约束要求机组在相邻时段间快速升降出力。当备用容量需求设置得比机组爬坡能力还大时，模型直接不可行。解决办法是把机组备用分解为上备用和下备用，同时加入爬坡备用约束，保证备用的调节速度跟得上系统需求。

6.2 风电场景数目与计算复杂度的平衡

场景削减数量是另一个影响模型成败的关键因素。削减到5个场景，模型求解很快，但结果偏差可能达到8%以上。削减到50个场景，结果精度提升有限，求解时间却从几秒涨到几十秒。我测下来10到15个场景是比较合理的配置，结果精度和求解效率兼得。

如果你发现在10个场景下，不同场景的调度结果差异很大，这说明场景代表性不够，可以从削减后的场景集合里检查一下场景间的距离，看看是否存在过于相似的冗余场景，或者削减时把异常场景强行删除了，导致覆盖不到极端情况。

6.3 Matlab代码调试的独家技巧

调试动态经济调度模型，我自己的习惯是先跑一个简化版本：取消风电随机性，把风电出力固定为期望值，模型退化成一个确定性动态经济调度。这一步跑通了，再逐步加入场景和随机约束，问题定位会容易很多。

再分享一个细节：Yalmip里二次成本函数（(a_i P_{i,t}^2)）容易让Cplex选择求解QP问题，如果机组数量大，求解时间可能明显上升。可以尝试用分段线性化来近似二次成本函数，精度损失很小但求解速度提升明显。很多商用电力系统软件在内部也是这么处理的。

6.4 参数敏感性分析建议

做参数敏感性分析时，重点关注三个参数：Weibull分布的形状参数、风速自相关系数、备用置信水平。这三个参数分别影响风电出力的波动幅度、时序关联性和备用需求水平，是模型里最有"调控手感"的旋钮。

形状参数从2.0调到2.6，风电出力的标准差会下降约15%，相应的备用需求也会下降，总成本可能变化3%左右。自相关系数从0.8调到0.9，进入模型的风速序列波动方式会不同，风电场景的峰谷位置变化，对机组爬坡压力有明显影响。如果能把这些参数敏感性分析的图表放进论文，审稿人通常会觉得研究做得很扎实。

7. 扩展方向与进阶思路

7.1 考虑储能系统的联合调度

风电随机性调度模型已经跑通以后，下一步可以往模型里加入储能系统。储能可以削峰填谷，把风电多发时段的多余电量储存起来，在风电出力不足时释放，能有效缓解系统对旋转备用的依赖。模型需要新增的决策变量包括储能的充放电功率、荷电状态（SOC），以及充放电效率约束。储能系统的SOC跨时段递推关系是动态经济调度模型很好的扩展维度。

7.2 多风电场相关性建模

多个风电场接入同一系统，如果地理位置较近，风电出力之间存在正相关性。把多个风电场的随机性独立建模会低估系统的整体风险，高估风电的消纳能力。处理手段是引入Copula函数描述多个风电场出力之间的相关性结构，或者直接用历史数据构造联合场景。对Matlab来说，Copula建模有现成工具箱，模型改造也不算复杂，学到这一层基本可以把这个课题写成一篇有一定深度的期刊论文。

7.3 从离线优化到滚动调度

动态经济调度在实操里通常不是一次算完，而是用模型预测控制的思想做滚动优化，也就是每次只执行当前时段的结果，下一时段重新求解。这种滚动模式对风电随机性的动态变化更加鲁棒。代码实现上，只需要在外层加一个for循环，每次更新风预测信息和系统状态，重新构建约束和优化目标就行。这个方向的工程感更强，放简历上是一个不错的加分项。

我在这个模型的调试过程中，最深的体会是：风电随机性的动态经济调度，难点不在于求解器能不能解，也不在于某个约束公式怎么写精准，而在于你是否真正理解了"随机性"在你的模型里到底扮演什么角色。如果只是把随机性当成一个扰动项，随便加一个正态分布就算处理了，那出来的结果往往经不起推敲。把场景生成、削减、备用约束、置信水平这条链路上的每个环节都做到有据可依，再回头审视调度结果，你会发现它对系统的描述能力和决策支持价值，是确定性模型无法比拟的。做这类研究，多花时间在数据分析和参数标定上，比单纯调调代码参数，收益要高得多。