说实话第一次在搜索框里敲下“风电随机性动态经济调度 Matlab”的时候我以为这又是个套娃模型把风电功率预测曲线往确定性经济调度里一插加两个约束跑个24时段的优化完事。真正动手以后我发现“随机性”这三个字带来的复杂度跳跃远比想象中大——同一个算例场景怎么生成、备用怎么留、机会约束怎么转化直接决定结果是“合理调度”还是“一个花哨的玩具”。这篇文章是我从问题建模到 Matlab 代码落地的一次完整复盘。面向的读者有三类正在写电力系统课程设计或论文复现的研究生、刚接触随机优化的工程师以及那些下载了别人代码但跑不通、一换数据就崩的朋友。我会先把模型讲透再给可直接修改的 Yalmip 写法最后分享我复现过程中踩过的坑。标题里的“风电随机性动态经济调度模型”听起来很短但背后要处理的问题其实足够写一本小册子。1. 风电随机性为什么让传统经济调度“失灵”1.1 传统确定性调度的“一条曲线思维”传统动态经济调度DED的逻辑非常直白给定未来24小时的负荷曲线和风电短期预测曲线再给定机组参数求一个满足所有约束的火电出力计划使总煤耗成本最低。数学上就是min Σᵢ Σₜ (aᵢ·Pᵢ,ₜ² bᵢ·Pᵢ,ₜ cᵢ)约束包括功率平衡、机组出力上下限、爬坡约束和备用约束。这套模型用了好几十年并且在风电渗透率不高的时候它确实够用。比如系统里风电装机只占总负荷的5%即使风电预测误差达到20%对系统功率平衡的影响也只有总负荷的1%左右系统自动发电控制AGC、一次调频完全能兜住。风电渗透率超过20%、甚至30%以后情况就变了。一条预测曲线背后的误差分布不再是可以忽略的小扰动而是可能在某个时段造成几百兆瓦功率缺口的“大象”。确定性调度只针对一条曲线求解它在真实风电场景下经常给出两种结果要么备用不足导致切负荷风险要么因为对风电过于乐观火电爬坡跟不上等风电骤降时机组加出力加不上去。1.2 随机性来源从风速预报误差到功率预测误差风电随机性严格来说来自风速。风速预报误差受气象条件、地形、预测时效影响再通过风机功率曲线 P(v) 映射成功率误差。工程落地时不方便做“风速到功率”的两层抽样因为功率曲线往往是一张查表数据不同风机类型差异也大。更常用的做法是直接对风电功率预测误差建模。最常见的假设是功率预测误差近似服从正态分布 N(0, σ²)σ 与预测提前时间、风电场规模正相关。也有的模型用 Beta 分布因为 Beta 分布定义在 [0,1] 区间适合描述功率在 0 到装机容量之间变化的有界性。实际算例里把预测误差取为预测值的10%~20%或者装机容量的某个比例基本都在合理范围。在调度模型中我们更关心的是“风电实际出力落在什么范围”。有三种描述方式区间描述Wₜ ∈ [Ŵₜ − εₜ, Ŵₜ εₜ]这是鲁棒优化的输入场景描述抽样 N 个可能的 Wₜ 向量每个带概率这是随机规划/场景法的输入分位数描述用预测误差的 α 分位数确定备用需求这是工程上最通用的输入。三种方式各有代价区间最保守场景最灵活但计算量大分位数最简单但难以刻画多时段相关性。初学阶段建议从分位数方案入手跑通后再扩展到场景模型。1.3 随机调度真正要回答的问题确定性模型回答的是“给定风电曲线谁该发多少电。”随机模型回答的是“面对一组可能的风电曲线现在就要决定机组出力计划使得在所有或绝大多数可能场景下系统都能安全运行且期望成本最低。”注意“现在就要决定”这个表述。实际调度中机组出力计划是提前一天定的风电实现在第二天才能知道。所以随机调度模型必须处理一个决策时序问题哪些变量现在定日前计划哪些变量明天根据风电实测再调实时调整。把这个时序理清楚后面的建模才不会乱。我自己的判断标准是如果一个随机调度模型里所有变量都不分场景、只有一个风电期望值参与约束那它本质上是确定性等价模型不是随机模型如果所有变量都分场景那它假设调度员在事前就知道风电场景属于“后验优化”同样失真。真正合理的随机模型一定在这两者之间。2. 数学模型的落地选择机会约束、备用分位数还是两阶段场景2.1 目标函数三项成本怎么配权重风电随机动态经济调度的目标函数不光是煤耗成本。因为引入了随机性必须考虑风电消纳和系统可靠性。常见目标函数是min Σₜ Σᵢ (aᵢ·Pᵢ,ₜ² bᵢ·Pᵢ,ₜ cᵢ) λᵥᵥ·Σₜ(Ŵₜ − Wₜ) λ꜀·Σₜ(Dₜ − ΣᵢPᵢ,ₜ − Wₜ)⁺后两项分别是弃风惩罚和切负荷惩罚。这里有个经验切负荷惩罚系数一定要远大于弃风惩罚系数通常要相差一个数量级以上。否则优化器会“聪明”地选择切一点负荷来满足煤耗最优结果就是把系统频率风险转移到用户侧这在工程上是绝对不允许的。弃风惩罚系数怎么定一个常用做法是取煤耗成本曲线边际成本的1.5~2倍也可以直接设为一个常数比如500元/MWh。只要它显著高于边际煤耗模型就会尽量消纳风电。我自己常用 500 这个量级因为它在大多数30机算例中能保证风电优先消纳又不会极端到让模型为了消纳最后一度电而触发机组的启停振荡。2.2 约束全集别漏掉爬坡与备用以下是一套24时段动态经济调度模型的完整约束清单每个约束我都标注了在 Matlab 里组装时需要注意的维度。功率平衡Σᵢ Pᵢ,ₜ Wₜ Dₜ。这里的 Wₜ 是火电实际消纳的风电功率不是风电预测值同时满足 0 ≤ Wₜ ≤ Ŵₜ。出力上下限Pᵢᵐⁱⁿ ≤ Pᵢ,ₜ ≤ Pᵢᵐᵃˣ。爬坡约束−Rᵢᵈᵒʷⁿ ≤ Pᵢ,ₜ − Pᵢ,ₜ₋₁ ≤ Rᵢᵘᵖ。特别注意 t1 时段要与初始出力 Pᵢ,₀ 比较别把初始出力漏掉。旋转备用Σᵢ min(Pᵢᵐᵃˣ − Pᵢ,ₜ, Rᵢᵘᵖ) ≥ Rₜⁿᵉᵉᵈ。Rₜⁿᵉᵉᵈ 是考虑风电预测误差后的备用需求。弃风约束0 ≤ Wₜ ≤ Ŵₜ当 Wₜ Ŵₜ 时目标函数里就会产生弃风惩罚。很多人第一次搭模型会漏掉备用里“爬坡受限”的部分。风电功率突然下降300MW如果某台机组虽然还有200MW的出力上升空间但爬坡速率只有20MW/时段那它在15分钟内只能提供80MW而不是200MW。Σ min(Pᵢᵐᵃˣ − Pᵢ,ₜ, Rᵢᵘᵖ) 这个写法就是为了防止备用容量被高估。虽然表达式麻烦一点但它才是物理上真实可用的旋转备用。2.3 随机性嵌入方式三种主流模型的比较这是整篇最值得花时间的部分。同一个问题有三种主流处理方式确定性等价CEE直接把 Wₜ 替换为预测值 Ŵₜ约束照常。优点是简单到极点缺点是忽略分布结果偏乐观。风电渗透率低时可用渗透率高时千万别用。机会约束 / 备用分位数要求在置信水平 α 下系统能承受风电波动即 Pr(ΣᵢPᵢ,ₜ Wₜ ≥ Dₜ Rₜ) ≥ α。工程上常转化为更朴素的形式把旋转备用需求设为预测误差的 α 分位数例如 Rₜ 1.645σₜ对应95%置信度。这是我后面主代码用的方案简单、稳健、物理意义清晰非常适合工程复现。场景两阶段随机规划生成 N 个风电场景第一阶段决定基准出力 Pᵢ,ₜ第二阶段根据场景 s 决定调整量 Δᵢ,ₜˢ期望总成本最小。表达式大致是min ΣₜΣᵢcostᵢ(Pᵢ,ₜ) Σₛ probₛ·ΣₜΣᵢcostᵢ(Pᵢ,ₜ Δᵢ,ₜˢ)这个模型信息含量最高但变量维度会扩大 N 倍。实际算例中 N 取 20~50 已经能覆盖大部分风险取 500 会导致求解器内存和时间双双爆炸。处理方式随机性描述计算量保守程度Matlab落地难度确定性等价期望值小低偏乐观最低备用分位数α分位数小中低场景两阶段N个场景大中高较高3. 场景生成与削减把“随机性”变成一个有限集合3.1 正态假设下的蒙特卡洛抽样如果选场景两阶段模型第一步是生成 N 个风电出力场景。做法通常是以预测曲线每时段的 Ŵₜ 作为均值假设误差 εₜ 服从 N(0, σₜ²)σₜ 取预测值的10%~20%或装机容量的某比例蒙特卡洛抽样得到场景 s 的 Wₜˢ Ŵₜ εₜˢ对越界值截断Wₜˢ 不能小于 0也不能超过风电场装机容量。这里有个容易翻车的细节如果对每个时段独立抽样生成的场景会非常“毛糙”相邻时段的风电出力跳变巨大与风电场实际出力的连续变化特性不符。更好的做法是引入时间相关性常用 AR(1) 模型εₜ ρ·εₜ₋₁ ηₜ ηₜ ~ N(0, σ·√(1−ρ²))ρ 取 0.8~0.95。ρ 越大场景越平滑也更符合实际风速的持续性特征。我在算例里常用 ρ0.9各时段残差标准差 σ 设为预测值的15%效果比较稳。3.2 场景削减从上千个到二十个的两种主流算法蒙特卡洛抽1000个场景直接把两阶段模型变量扩大1000倍求解器多半吃不消。场景削减的目标是用少量场景逼近原始场景集的概率分布特征。主流有两种K-means/K-medoids 削减。把每个场景看成24维向量聚成 K 类用每类中心代表该类概率为该类场景数除以总数。优点实现简单Matlab 里一行 kmeans 就能跑缺点均值聚类容易抹掉尾部分布极端场景可能被聚掉而极端场景恰恰是随机调度最关心的。快速前代削减Forward Selection。这是随机规划文献最常用的方法。核心思路是迭代选场景使得被选场景集合与剩余场景之间的某种距离度量Kantorovich 距离或 Wasserstein 距离减少最多选够 K 个后再重新分配概率。效果上比 K-means 更保尾代价是计算量约 O(N²K)原始场景多时会慢。我实际测试的经验是原始场景数在 200 以内时快速前代削减效果明显好于 K-means原始场景上千时先用 K-means 粗削减到 100再用快速前代削减到 20速度和效果兼顾。这部分代码建议自己写网上很多现成脚本是小规模示例直接套大规模场景矩阵很容易溢出。3.3 削减完必须做的质量检验别削减完就闷头跑优化。先用三个指标检查场景集质量均值对齐度削减后场景集的均值曲线与原始场景均值曲线的偏差应小于风电装机容量的2%。偏差大说明削减过程把主体趋势弄歪了。覆盖率统计削减后场景在各时段的最大值/最小值与原始场景集对比。如果削减后区间明显收窄说明尾部风险被丢掉运行结果会偏乐观。概率集中度削减后的场景概率分布应该比较均匀。如果某个场景概率超过 0.3说明削减集太集中建议增加场景数或换削减方法。这个环节容易被忽略但恰恰是决定随机调度可靠性的关键。模型再漂亮输入场景集是歪的结果就是精致的错误。我见过不少论文复现代码场景数取10个就敢上两阶段模型跑出来的成本比确定性等价还低这就是场景质量崩掉的典型信号。4. Matlab 实现从数据结构到 Yalmip 求解4.1 输入数据结构设计我强烈建议把全部参数放进 struct方便换算例、换数据。核心字段可以这样组织% gen 每行代表一台机组列约定如下 % 列1: a煤耗二次系数, 列2: b煤耗一次系数, 列3: c煤耗常数 % 列4: Pmin, 列5: Pmax, 列6: 上行爬坡速率 Rup % 列7: 下行爬坡速率 Rdown, 列8: 初始出力 P0 caseData.gen [...]; % ng x 8 caseData.load [...]; % 24 x 1负荷曲线 caseData.wind [...]; % 24 x 1风电预测均值 caseData.windSigma [...]; % 24 x 1预测误差标准差 caseData.scen [...]; % 24 x Ns风电场景矩阵每列一个场景 caseData.prob [...]; % 1 x Ns场景概率用 struct 的好处是函数接口干净写约束组装函数时不会出现几十个形参。后面不管调 K-means 削减、做敏感性分析还是换 IEEE 算例只要替换 caseData 字段就行。4.2 主模型代码备用分位数版本先给最实用、最容易跑通的版本——备用分位数方案。这个模型不需要场景集只需要预测曲线和误差标准差计算量小适合作为第一版跑通和对照基线。%% 参数初始化 T 24; ng size(gen, 1); P sdpvar(ng, T, full); % 各机组各时段出力 W sdpvar(1, T, full); % 各时段实际消纳风电 % 旋转备用需求取风电预测误差的95%分位数 负荷预测误差贡献 Rneed 1.645 * caseData.windSigma 0.02 * caseData.load; Constraints []; % 功率平衡 Constraints [Constraints, sum(P, 1) W caseData.load]; % 风电消纳上下限 Constraints [Constraints, 0 W caseData.wind]; % 出力上下限 Constraints [Constraints, repmat(gen(:,4), 1, T) P repmat(gen(:,5), 1, T)]; % 爬坡约束含第1时段与初始出力比较 Constraints [Constraints, -gen(:,7) P(:,2:T) - P(:,1:T-1) gen(:,6)]; Constraints [Constraints, -gen(:,7) P(:,1) - gen(:,8) gen(:,6)]; % 旋转备用考虑机组爬坡上限的实际可用备用 availReserve min(repmat(gen(:,5), 1, T) - P, repmat(gen(:,6), 1, T)); Constraints [Constraints, sum(availReserve, 1) Rneed]; % 目标函数煤耗 弃风惩罚 obj sum(sum(repmat(gen(:,1), 1, T).*P.^2 ... repmat(gen(:,2), 1, T).*P ... repmat(gen(:,3), 1, T))) ... 500 * sum(caseData.wind - W); % 求解 ops sdpsettings(solver, cplex, verbose, 2); optimize(Constraints, obj, ops);注意几个细节min函数对矩阵逐元素取最小值得到 ng×T 矩阵再用sum(...,1)按列求和得到 1×T 的可用备用向量。这个写法依赖 Matlab 2016b 及以上版本的隐式扩展如果还在用老版本把-gen(:,7)改成-repmat(gen(:,7), 1, T)之类的显式扩展。4.3 场景两阶段模型的关键写法如果要用场景模型核心是区分基准出力变量和场景调整量。示意代码如下P0 sdpvar(ng, T, full); % 第一阶段基准出力不依赖场景 delta sdpvar(ng, T, Ns, full); % 第二阶段各场景下的调整量 W sdpvar(1, T, Ns, full); % 各场景下实际消纳风电 Constraints []; for s 1:Ns % 场景 s 的功率平衡 Constraints [Constraints, sum(P0, 1) sum(delta(:,:,s), 1) ... W(:,:,s) caseData.load]; % 场景 s 的风电消纳上限 Constraints [Constraints, 0 W(:,:,s) caseData.scen(:,s)]; % 场景 s 下的机组限幅 Constraints [Constraints, repmat(gen(:,4), 1, T) P0 delta(:,:,s) ... repmat(gen(:,5), 1, T)]; % 场景 s 下的爬坡 Constraints [Constraints, -gen(:,7) (P0(:,2:T) delta(:,2:T,s)) ... - (P0(:,1:T-1) delta(:,1:T-1,s)) gen(:,6)]; end % 目标第一阶段煤耗 各场景概率加权调整成本 obj sum(sum(repmat(gen(:,1), 1, T).*P0.^2 ... repmat(gen(:,2), 1, T).*P0 ... repmat(gen(:,3), 1, T))); for s 1:Ns obj obj caseData.prob(s) * sum(sum(repmat(gen(:,1), 1, T).*(P0 delta(:,:,s)).^2 ... repmat(gen(:,2), 1, T).*(P0 delta(:,:,s)) ... repmat(gen(:,3), 1, T))); end ops sdpsettings(solver, cplex, verbose, 2); optimize(Constraints, obj, ops);这个循环写法比三维隐式扩展更易读而且不容易把维度弄错。代价是多写几行但对调试非常友好。注意第二阶段调整量不一定非要用二次煤耗工程上也常用线性惩罚替代比如正向调整价格高一点、负向调整价格低一点这样目标函数变成线性求解压力更小。4.4 求解器选型实测Matlab 自带的linprog/quadprog也能解但有两个明显短板。第一quadprog内部求解器对稀疏性和预处理的利用不如商业求解器规模一大就慢第二随机两阶段模型变量多、约束多自带求解器经常需要手动调容差否则会出现数值警告。我实测过同样一个 1500 个约束、5000 个变量的算例quadprog要 30~60 秒CPLEX 只要 2~5 秒Gurobi 和 CPLEX 接近。如果是课程作业quadprog完全够用如果要跑论文仿真、做敏感性分析强烈建议装 CPLEX 或 Gurobi然后用 Yalmip 统一建模。一个非常实际的坑很多小白卡在“Yalmip 显示 No suitable solver”其实是没装外部求解器或者装了但没添加路径。CPLEX 有学术版许可Gurobi 有学术版和免费限制版本。实在都没有可以用 SCIPYalmip 也支持。先跑通模型再优化求解器这个顺序别反了。5. 30机24时段算例从搭建到结果解读5.1 算例参数怎么设才能不翻车很多复现失败是算例参数不自洽。最典型的是机组出力上限加起来小于负荷峰值或者爬坡速率与备用时间尺度不匹配。我自己调试 30 机算例的参数原则总装机容量 负荷峰值 × 1.6~1.8给足备用空间单机 Pmin 不能太高否则夜间负荷低谷时被迫高煤耗风电装机取总负荷的 25%~30%这样随机性影响足够明显每台机组的爬坡速率按“15 分钟内能爬满备用”来反推即 Rup ≥ (Pmax − 当前出力)/4。没有标准算例时用 IEEE 30 机数据配合一条合理的负荷曲线就行。关键是把模型先跑通再慢慢调参数。不要在第一次运行时就追求精确结果先把约束全拼对数值不出警告再去做敏感性分析。5.2 随机性对成本与出力的影响固定场景数和备用水平运行后重点看三类结果总煤耗、弃风量、最大备用缺额。一个很典型的结论随着备用需求从 0 增加到预测误差的 95% 分位数总煤耗成本大约增加 3%~8%但切负荷风险显著下降。这个“成本换可靠性”的权衡曲线是论文里最常用的一张图。有些代码把备用系数调得过高比如取 3σ成本增加超过 15%这时候就要提醒自己过度保守的调度和过度乐观的调度一样糟糕。场景两阶段模型下的对比更有意思。固定场景数 20 个期望成本通常比确定性等价高 2%~5%但最坏场景下的功率平衡违反量几乎为 0。换句话说随机模型多花一点钱买的是极端场景下的安全性。风电渗透率越高这个价差越明显。5.3 结果图怎么读我习惯输出三张图机组出力堆叠图每时段各机组出力叠加上方留出的空白是风电消纳空间。如果堆叠顶部和负荷曲线之间的空隙很大说明火电出力压得不够弃风可能性高。风电消纳柱状图预测风电和实际消纳风电对比差值就是弃风。如果弃风集中在负荷低谷时段说明是系统调峰能力不足而不是风电本身的问题。备用裕度曲线每时段可用备用减去备用需求出现负值的位置就是要重点关注的风险时段。读图时有一个容易误判的点备用裕度为负不代表系统一定出事只是说明在该时段的极端风电偏差下AGC 可能来不及恢复频率。但备用裕度大面积为负说明模型约束没建对或者备用需求设得太高需要回去检查参数。6. 复现随机调度代码时我踩过的坑6.1 场景数太少优化器开始“自欺欺人”我第一次只生成 10 个场景就跑两阶段模型结果目标成本比确定性等价还低。检查后发现10 个场景的风电均值偏差了 8%而且没有覆盖任何极端高/低风电场景。相当于让优化器在一个被“平均化”的世界里做计划自然成本低。后来把场景加到 50 个再削减到 20 个结果才合理。经验是场景数小于 15 时先做场景质量检验削减后如果某个场景概率超过 0.2视为警告。随机优化最怕的不是随机性本身而是把随机性平均到看不见。6.2 爬坡约束第一时段的老问题用矩阵化方式写爬坡约束时最容易犯的错是漏掉 t1 与初始出力的约束。我一开始直接用P(:,2:end) - P(:,1:end-1)结果第一时段完全不受爬坡限制优化器把第一时段出力直接拉到最高后面所有爬坡约束都被“带偏”。加一行Constraints [Constraints, -gen(:,7) P(:,1) - gen(:,8) gen(:,6)];就解决了。这个坑非常隐蔽因为模型能正常求解只是结果从第一时段开始就是错的。6.3 备用约束被爬坡上限卡死备用约束写成 Σ(Pmax − P) ≥ R 时系统看似备用充足真到风电骤降时某些机组爬坡速率根本供不上。改成 Σ min(Pmax − P, Rup) ≥ R 后备用容量缩水不少模型“被迫”增加开机或提高低负荷机组的出力这才接近物理真实状态。这个坑的教训是模型里的“可用性”一定要考虑时间尺度15 分钟能调出来的功率才是真备用。6.4 二次项数值尺度问题煤耗二次项系数 aᵢ 可能小到 0.0001一次项系数 bᵢ 到几十目标函数总量级 10⁶。如果各系数跨了太多数量级求解器对偶算法容易出现数值警告甚至报“伪不可行”。解决办法有两个一是做目标归一化把整个目标除以 1000 再优化二是设sdpsettings(cplex.optimization.simplex.display,off)这类容差配置但更直接的办法还是保持系数同数量级别让小系数项被容差吃掉。6.5 独立抽样丢了时间相关性前面提过独立对每时段抽样会让场景毛糙备用需求虚高。我换用 AR(1) 模型ρ0.9 之后场景平滑度明显改善极端场景也保留得不错。另外我习惯把随机种子固定比如rng(2024)这样复现实验时结果可重复。做对比实验时固定同一份场景集才能说清楚“是调度模型改变了结果而不是随机数改了结果”。我因为这个问题吃过亏现在所有随机调度实验都强制固定随机种子和场景文件。到这里一套完整的风电随机性动态经济调度模型就算真正落地了。我最后采用的配置是原始场景 50 个、快速前代削减到 20 个、备用分位数取 95%、弃风惩罚 500 元/MWh。这个组合在成本和可靠性之间取得了不错的平衡。你也可以直接套这个配置但换算例时一定先跑一遍场景质量检验。如果你的模型里风电渗透率更高或者有储能参与可以把储能变量作为一个额外决策块接入同一套 Yalmip 框架原理完全一致。祝跑通。
