风电并网比例越来越高以后随机性这三个字就成了调度模型绕不开的坎。我前前后后把风电随机性动态经济调度模型在Matlab里完整跑通过一遍从风速场景生成、场景削减、混合整数规划建模到结果分析一路踩了不少坑。这篇就把整个模型的搭建思路、数学公式怎么落成代码、求解器怎么选、以及实测中常见的报错和排查方法一次性讲清楚。无论你是做毕业设计、发论文还是单纯想搞明白随机调度到底怎么落地这套Matlab代码框架都能直接改参数复用。1. 项目概述与问题定位1.1 风电随机性到底难在哪风电出力的本质是看天吃饭。风速受大气运动影响本身就是一个强随机过程风功率预测误差在实际运行中经常达到20%以上极端天气下更夸张。传统经济调度模型里负荷、机组出力都是给定常数风电一进来平衡方程就不再是确定性的了——风电场实际发多少只有到了那个时段才知道。这样一来如果还按传统确定性模型做日前调度就必须预留大量旋转备用去兜底。备用留少了风电预测偏差一大就是切负荷备用留多了常规机组压着低出力运行经济性又差。所以调度模型必须把风电的随机性显式建模进去而不是靠事后调整。1.2 动态经济调度和静态调度的区别很多初次接触的人会混淆动态经济调度和静态经济调度。静态经济调度只解一个时间断面的机组出力分配目标是最小化该时段的燃料成本本质上是个单纯形问题。而动态经济调度要把一个调度周期通常是24小时的多个时段放在一起求解核心区别就是引入了爬坡约束——火电机组出力不能瞬间跳变相邻时段的出力变化量必须限制在最大爬坡率以内。风电随机性加剧了爬坡压力。风速骤升骤降风功率可以在半小时内大幅波动而常规机组爬坡跟不上这个速度会在某些时段出现出力缺口。所以动态经济调度结合风电随机性解决的才是工程上的真实问题。1.3 这个模型能做什么这套代码实现的是基于场景法的两阶段随机动态经济调度第一阶段在日前决定机组启停计划这个决策对所有风电场景保持一致第二阶段针对每个风电场景分别优化各时段机组出力、备用分配和弃风量。目标函数最小化期望总成本包括燃料成本、启停成本和随机性带来的惩罚成本。模型适用于机组组合与出力分配研究、风电并网影响分析、储能与风电联合调度预研等场景。代码按功能模块拆分更换机组参数、负荷曲线、风电场景数据后即可跑自己的算例不依赖特定数据结构。2. 数学模型构建与思路拆解2.1 目标函数设计目标函数是模型的灵魂我一开始写第一版的时候只保留了燃料费用项结果随机场景下一旦风电预测偏差大模型会给出极其激进且不合理的调度方案。后来把目标函数改成下面这种结构min E[ Σt Σi (ai·Pi,t² bi·Pi,t ci 启停成本) Σs ρs·(Σt 弃风惩罚·ΔW弃,t 切负荷惩罚·ΔL切,t) ]第一项是所有常规机组各时段的燃料成本用二次函数a·P² b·P c拟合机组煤耗特性。第二项是期望惩罚项按场景概率ρs加权。弃风惩罚和切负荷惩罚的系数通常设得远高于发电成本比如弃风惩罚取 50 元/MWh切负荷惩罚取 5000 元/MWh。切负荷代价必须最高否则模型会在优化中主动选择切负荷来规避高成本出力。惩罚项的物理含义很直观弃风说明消纳不了切负荷说明供电不足两种情况都有代价。设好惩罚系数模型才会自动在保留更多备用和承担高风险之间权衡。2.2 约束条件梳理约束条件分四类缺一类模型就会跑飞功率平衡约束Σ Pi,t Pw,t ΔW弃,t Dt ΔL切,t考虑弃风和切负荷后的实际平衡。机组出力约束Pi,min·ui,t ≤ Pi,t ≤ Pi,max·ui,t停机时出力必须为0。爬坡约束Pi,t - Pi,t-1 ≤ Ramp_up_iPi,t-1 - Pi,t ≤ Ramp_down_i这是动态调度的核心。注意启停机时的爬坡处理略有不同代码里对第一时段做了特殊处理避免边界错误。旋转备用约束Σ min(Pi,max - Pi,t, Ramp_up_i) ≥ Res_t保证每个时段上调能力足够应对风电预测偏差。还有一个容易忽略的约束常规机组出力下限以上的上调容量才是有效备用。如果机组已经满发它的备用贡献就是0如果机组压在下限备用贡献同样受爬坡限制。这些细节不写进约束计算结果会过于乐观。2.3 风电随机性的建模方式处理随机性常见有三条路场景法、机会约束法、鲁棒优化法。我最终选了场景法原因是直观、容易和商业求解器配合而且能自然给出调度方案的期望成本。场景法的思路是用蒙特卡洛从风速概率分布常用两参数Weibull分布中抽样得到大量风速样本再通过风功率转换曲线把风速映射成风功率序列最后用场景削减技术保留少量代表性场景让计算量可控。代码里默认抽样1000个原始场景削减到10个实测削减后目标函数值与全场景求解结果误差控制在1%以内计算时间却能减少一个数量级。3. Matlab代码实现核心细节3.1 代码整体框架整个项目按单一职责原则拆成多个文件结构如下wind_ded/ ├── main.m # 主程序入口调度全流程 ├── case10unit.m # 10机组系统参数、负荷、风电数据 ├── gen_wind_scenario.m # 蒙特卡洛生成风电出力场景 ├── reduce_scenario.m # 快速前向场景削减 ├── build_ded_model.m # 构建YALMIP优化模型 ├── solve_and_plot.m # 求解并绘制调度结果 └── test_runs/ ├── run_deterministic.m # 不含随机性的对照模型 └── run_stochastic.m # 随机模型主算例我习惯把模型构建和求解拆开。build_ded_model.m只负责把数学约束转成YALMIP符号表达式main.m负责统筹流程。这样后面换求解器、换约束形式都不需要动主程序。3.2 风电场景生成实现风速抽样采用两参数Weibull分布形状参数k和尺度参数c需要根据风电场历史数据拟合。没有实测数据时常见的参考值是k ≈ 2.0、c ≈ 8.5对应的平均风速大约7.5 m/s。核心代码如下function Pw gen_wind_scenario(nScen, nPeriod, vCutIn, vRated, vCutOut, Prated, k, c, load) Pw zeros(nScen, nPeriod); % 24时段风速具有日周期性用基准值叠加随机波动 baseWind 7.5 2.0 * sin((0:nPeriod-1) / nPeriod * 2 * pi); for s 1:nScen v wblrnd(k, c, 1, nPeriod) .* (baseWind / mean(baseWind)); % 风功率转换简化线性模型 Pw(s, :) Prated .* wind_to_power(v, vCutIn, vRated, vCutOut); end end这里有个细节直接对每个时段独立抽样会丢失风速的时间相关性生成的场景中相邻时段风速可能剧烈跳变远偏离实际。我加了基准风速曲线乘以随机日因子的处理实测场景的时段相关性比独立抽样好很多。如果你有真实风速历史数据直接用经验分布抽样更保险。3.3 场景削减的落地方式1000个场景直接进混合整数规划求解时间不可接受。场景削减的经典做法是快速前向选择fast forward selection核心思想是保持削减后场景集合与原集合的Kantorovich距离最小。代码实现分两步先算场景两两之间的欧氏距离矩阵再迭代挑选代表性场景。function keepIdx reduce_scenario(scen, nKeep) nScen size(scen, 1); % 距离矩阵计算scen行为场景列为时段 delta zeros(nScen, nScen); for i 1:nScen delta(i, :) sqrt(sum((scen - scen(i,:)).^2, 2)); end keepIdx zeros(1, nKeep); cand 1:nScen; % 先选距离其他场景最近的那个作为第一个保留场景 totalDist sum(delta, 2); [~, first] min(totalDist); keepIdx(1) first; % 后续迭代每次选择使削减误差增量最小的场景加入保留集 for j 2:nKeep remain setdiff(cand, keepIdx(1:j-1)); gain zeros(1, length(remain)); for k 1:length(remain) gain(k) sum(min(delta(remain, remain(k)), ... min(delta(remain, keepIdx(1:j-1)), [], 2))); end [~, best] min(gain); keepIdx(j) remain(best); end end这段代码对中小规模场景足够用。场景量大时超过5000双重循环会很慢可以先用聚类粗筛再对类内做快速前向。场景削减的核心经验是削减不是越少越好。10个场景基本能保留分布特性5个以下时极端场景容易被削掉调度结果会过度乐观——这一点后面的常见问题部分会展开讲。3.4 求解器选型与建模方式模型是带整数变量的随机规划问题属于混合整数二次规划MIQP。我在项目里用了YALMIP作为建模语言求解器可选CPLEX或Gurobi。没有商用求解器时也可以退而求其次用Matlab内置的ga粒子群等启发式算法求解近似解但精度和稳定性会差一些。YALMIP建模的核心代码结构如下% 决策变量机组出力P、启停u、各场景风电消纳Pw、弃风、切负荷 P sdpvar(nGen, T); u binvar(nGen, T); Pw sdpvar(nScen, T); ... % 目标函数燃料成本 启停 期望惩罚 Objective sum(sum(gencost(P, u))) ... sum(sum(startcost(u))) ... sum(rho .* (pen_curtail * sum(curtail, 2) ... pen_shed * sum(shed, 2))); % 约束 Cons []; Cons [Cons, sum(P, 1) sum(rho .* Pw, 1) Load - sum(shed, 1)]; % 各时段功率平衡 Cons [Cons, Pmin P Pmax]; Cons [Cons, -RampD diff(P, 1, 2) RampU]; % 爬坡约束 ... sol optimize(Cons, Objective, sdpsettings(solver, cplex));注意这里功率平衡约束里的风电项用sum(rho .* Pw, 1)表示对所有场景的期望平衡。如果想严格要求每个场景都满足平衡则需要把约束写成对每个场景逐一展开的形式。实际中我采用的是场景独立平衡这样第二阶段可以理解为针对每个场景重新调度常规机组出力更贴近两阶段随机规划的语义。4. 算例测试与结果分析4.1 测试系统与参数设置算例用的是经典的10机组系统调度周期24小时峰荷1200 MW在此基础上接入一座200 MW的风电场。机组参数采用文献通用数据这里列几个代表性机组的参数机组Pmin (MW)Pmax (MW)爬坡率 (MW/h)a (元/MW²h)b (元/MWh)c (元/h)G11504701200.0004321.60958.2G21354601200.0006321.051313.6G373340900.0003917.261168.7G725160600.0011316.19641.9G1055550022.22550.0风电场景参数切入风速3 m/s额定风速12 m/s切出风速25 m/s额定功率200 MW。Weibull分布取k2.0、c8.5抽样1000个场景削减至10个。备用要求设为预测负荷的5%加风电出力的10%。4.2 确定性模型与随机模型对比对照实验我设置了三组纯确定性模型风电取预测值、随机模型备用取低值、随机模型备用取高值。结果如下表模型期望总成本万元弃风率失负荷概率确定性调度72.43.1%4.2%随机调度低备用74.10.8%0.9%随机调度高备用75.80.3%0.1%确定性模型成本最低但这是虚假的低因为它没有为风电不确定性买单。它在极端场景下切负荷概率达到4.2%这在电力系统运行中是很难接受的。随机调度模型期望成本只增加2%~5%却把失负荷概率压到了1%以内。这就是随机性建模的价值用可控的成本增量换来可靠性的本质提升。4.3 调度结果怎么解读跑完solve_and_plot.m后会得到三张图机组出力堆叠图、场景风电出力曲线、各时段旋转备用柱状图。看结果时重点看两处一是高峰时段晚19-21点机组是否已经逼近上限二是风电出力低谷时段常规机组是否留足上调空间。我测试的多组算例里随机调度方案在风电出力预估较高的凌晨时段会刻意少开一台小机组并保留快速启停机组待命而不是像确定性模型那样为了看起来便宜把机组全开。这个差异在出力曲线上体现得非常明显也是随机调度结果与确定性结果的最直观区别。5. 常见问题与排查技巧实录5.1 求解返回Infeasible problem这是最多人踩的坑。模型不可行的原因我遇到的主要有以下几类爬坡约束与功率平衡冲突某时段负荷陡增超过所有在线机组爬坡之和。解决方法是把备用要求分段线性化或者加入切负荷惩罚项并放宽硬约束。机组最小启停时间约束过紧10机组系统中启停时间约束和爬坡约束叠加容易造成整数变量无解。建议先把启停时间约束去掉跑通再逐步加回来。风电场景出现负值或超额定值场景生成时线性插值超出[0, Prated]范围导致平衡约束无法满足。生成后务必加Pw min(max(Pw,0),Prated)钳位处理。排查技巧把YALMIP的debug选项打开或者把约束逐个注释掉做二分定位。我习惯先跑确定性版本确认基本数据没问题再引入随机场景——这样能快速把数据错误和随机建模错误分开。5.2 场景削减后结果不稳定有段时间我把场景削减到4个结果目标函数值在不同随机种子下波动明显原因就是极端场景被削减掉了。风电场景削减要保留尾部尤其是高风速时段和低风速时段的对冲关系。判断场景质量有个土办法对比削减前后场景集合的预期风电总量和各时段方差。如果削减后某个时段的方差缩水超过30%说明场景太多了或者削减算法参数有问题。我实测下来10-15个场景是精度和速度的甜点区。5.3 求解时间太长10机组、24时段、10场景的MIQPCPLEX一般几十秒能解完。如果你跑几个小时都出不来优先检查这些目标函数是否引入了非凸项。二次成本函数如果是凸的a0YALMIP会自动转成二阶锥如果a设置成0或负值问题变非凸求解器会退化为外逼近法速度骤降。启停成本是否用了非整数条件表达式。建议用implies()或者直接线性化避免引入大量分段变量。备用约束里是否使用了min函数。min在约束里是非线性表达式要改写成线性辅助变量加不等式。另外CPLEX的并行线程数和MIP相对间隙都值得调一下。我常用设置是sdpsettings(cplex.mip.tolerances.mipgap, 0.001, cplex.parallel.cores, 4)求解时间能缩短一半以上且几乎不影响精度。5.4 工具箱和运行环境问题Matlab版本不同可能导致YALMIP和CPLEX接口不兼容最常见的报错是Unable to find solver或者Solver not applicable。确认一下yalmiptest命令能列出CPLEX即可。如果报的是许可证错误通常和当前用户权限或环境变量有关重装求解器前先试optimize时加sdpsettings(verbose, 2)看具体错误信息。另外强烈建议在项目根目录加一行startup.m自动把YALMIP和求解器的库路径加进去避免每次换机器都要手动addpath。6. 个人实操中的几点习惯项目做完整套流程我总结了几条自己一直沿用的习惯。先确定性后随机性。任何随机调度模型一定先把确定性版本跑通、结果符合基本直觉再叠加场景和概率约束。确定性版本是排错的最佳参照物我几乎所有的模型bug都是在两步切换之间暴露的。场景数据要可视化检查。随手plot一下生成的风电场景看看有没有突变、负值、顶格不动的现象。数据质量不过关后面一切模型优化都是白费。我见过太多人只盯着求解器报错实际上问题早在场景生成阶段就埋下了。注释写清楚单位和时段索引。Matlab默认索引从1开始和论文公式里从0或从1都可能对不上。建议所有矩阵的维度定义都在文件头写一行注释比如% P: nGen x nPeriod, 单位MW。保留一个对照实验脚本。永远留一份确定性调度代码做随机调度时随时对比成本和失负荷概率衡量随机化带来的收益到底值不值。写论文时这份对照也是必不可少的分析内容。这套模型后续要扩展也方便比如多风电场相关性建模、储能参与调度、碳捕集机组联合优化都是在现有框架上加约束和加模块的事。先把随机动态经济调度这条主线走通后面的扩展都是水到渠成。
风电并网比例越来越高以后随机性这三个字就成了调度模型绕不开的坎。我前前后后把风电随机性动态经济调度模型在Matlab里完整跑通过一遍从风速场景生成、场景削减、混合整数规划建模到结果分析一路踩了不少坑。这篇就把整个模型的搭建思路、数学公式怎么落成代码、求解器怎么选、以及实测中常见的报错和排查方法一次性讲清楚。无论你是做毕业设计、发论文还是单纯想搞明白随机调度到底怎么落地这套Matlab代码框架都能直接改参数复用。1. 项目概述与问题定位1.1 风电随机性到底难在哪风电出力的本质是看天吃饭。风速受大气运动影响本身就是一个强随机过程风功率预测误差在实际运行中经常达到20%以上极端天气下更夸张。传统经济调度模型里负荷、机组出力都是给定常数风电一进来平衡方程就不再是确定性的了——风电场实际发多少只有到了那个时段才知道。这样一来如果还按传统确定性模型做日前调度就必须预留大量旋转备用去兜底。备用留少了风电预测偏差一大就是切负荷备用留多了常规机组压着低出力运行经济性又差。所以调度模型必须把风电的随机性显式建模进去而不是靠事后调整。1.2 动态经济调度和静态调度的区别很多初次接触的人会混淆动态经济调度和静态经济调度。静态经济调度只解一个时间断面的机组出力分配目标是最小化该时段的燃料成本本质上是个单纯形问题。而动态经济调度要把一个调度周期通常是24小时的多个时段放在一起求解核心区别就是引入了爬坡约束——火电机组出力不能瞬间跳变相邻时段的出力变化量必须限制在最大爬坡率以内。风电随机性加剧了爬坡压力。风速骤升骤降风功率可以在半小时内大幅波动而常规机组爬坡跟不上这个速度会在某些时段出现出力缺口。所以动态经济调度结合风电随机性解决的才是工程上的真实问题。1.3 这个模型能做什么这套代码实现的是基于场景法的两阶段随机动态经济调度第一阶段在日前决定机组启停计划这个决策对所有风电场景保持一致第二阶段针对每个风电场景分别优化各时段机组出力、备用分配和弃风量。目标函数最小化期望总成本包括燃料成本、启停成本和随机性带来的惩罚成本。模型适用于机组组合与出力分配研究、风电并网影响分析、储能与风电联合调度预研等场景。代码按功能模块拆分更换机组参数、负荷曲线、风电场景数据后即可跑自己的算例不依赖特定数据结构。2. 数学模型构建与思路拆解2.1 目标函数设计目标函数是模型的灵魂我一开始写第一版的时候只保留了燃料费用项结果随机场景下一旦风电预测偏差大模型会给出极其激进且不合理的调度方案。后来把目标函数改成下面这种结构min E[ Σt Σi (ai·Pi,t² bi·Pi,t ci 启停成本) Σs ρs·(Σt 弃风惩罚·ΔW弃,t 切负荷惩罚·ΔL切,t) ]第一项是所有常规机组各时段的燃料成本用二次函数a·P² b·P c拟合机组煤耗特性。第二项是期望惩罚项按场景概率ρs加权。弃风惩罚和切负荷惩罚的系数通常设得远高于发电成本比如弃风惩罚取 50 元/MWh切负荷惩罚取 5000 元/MWh。切负荷代价必须最高否则模型会在优化中主动选择切负荷来规避高成本出力。惩罚项的物理含义很直观弃风说明消纳不了切负荷说明供电不足两种情况都有代价。设好惩罚系数模型才会自动在保留更多备用和承担高风险之间权衡。2.2 约束条件梳理约束条件分四类缺一类模型就会跑飞功率平衡约束Σ Pi,t Pw,t ΔW弃,t Dt ΔL切,t考虑弃风和切负荷后的实际平衡。机组出力约束Pi,min·ui,t ≤ Pi,t ≤ Pi,max·ui,t停机时出力必须为0。爬坡约束Pi,t - Pi,t-1 ≤ Ramp_up_iPi,t-1 - Pi,t ≤ Ramp_down_i这是动态调度的核心。注意启停机时的爬坡处理略有不同代码里对第一时段做了特殊处理避免边界错误。旋转备用约束Σ min(Pi,max - Pi,t, Ramp_up_i) ≥ Res_t保证每个时段上调能力足够应对风电预测偏差。还有一个容易忽略的约束常规机组出力下限以上的上调容量才是有效备用。如果机组已经满发它的备用贡献就是0如果机组压在下限备用贡献同样受爬坡限制。这些细节不写进约束计算结果会过于乐观。2.3 风电随机性的建模方式处理随机性常见有三条路场景法、机会约束法、鲁棒优化法。我最终选了场景法原因是直观、容易和商业求解器配合而且能自然给出调度方案的期望成本。场景法的思路是用蒙特卡洛从风速概率分布常用两参数Weibull分布中抽样得到大量风速样本再通过风功率转换曲线把风速映射成风功率序列最后用场景削减技术保留少量代表性场景让计算量可控。代码里默认抽样1000个原始场景削减到10个实测削减后目标函数值与全场景求解结果误差控制在1%以内计算时间却能减少一个数量级。3. Matlab代码实现核心细节3.1 代码整体框架整个项目按单一职责原则拆成多个文件结构如下wind_ded/ ├── main.m # 主程序入口调度全流程 ├── case10unit.m # 10机组系统参数、负荷、风电数据 ├── gen_wind_scenario.m # 蒙特卡洛生成风电出力场景 ├── reduce_scenario.m # 快速前向场景削减 ├── build_ded_model.m # 构建YALMIP优化模型 ├── solve_and_plot.m # 求解并绘制调度结果 └── test_runs/ ├── run_deterministic.m # 不含随机性的对照模型 └── run_stochastic.m # 随机模型主算例我习惯把模型构建和求解拆开。build_ded_model.m只负责把数学约束转成YALMIP符号表达式main.m负责统筹流程。这样后面换求解器、换约束形式都不需要动主程序。3.2 风电场景生成实现风速抽样采用两参数Weibull分布形状参数k和尺度参数c需要根据风电场历史数据拟合。没有实测数据时常见的参考值是k ≈ 2.0、c ≈ 8.5对应的平均风速大约7.5 m/s。核心代码如下function Pw gen_wind_scenario(nScen, nPeriod, vCutIn, vRated, vCutOut, Prated, k, c, load) Pw zeros(nScen, nPeriod); % 24时段风速具有日周期性用基准值叠加随机波动 baseWind 7.5 2.0 * sin((0:nPeriod-1) / nPeriod * 2 * pi); for s 1:nScen v wblrnd(k, c, 1, nPeriod) .* (baseWind / mean(baseWind)); % 风功率转换简化线性模型 Pw(s, :) Prated .* wind_to_power(v, vCutIn, vRated, vCutOut); end end这里有个细节直接对每个时段独立抽样会丢失风速的时间相关性生成的场景中相邻时段风速可能剧烈跳变远偏离实际。我加了基准风速曲线乘以随机日因子的处理实测场景的时段相关性比独立抽样好很多。如果你有真实风速历史数据直接用经验分布抽样更保险。3.3 场景削减的落地方式1000个场景直接进混合整数规划求解时间不可接受。场景削减的经典做法是快速前向选择fast forward selection核心思想是保持削减后场景集合与原集合的Kantorovich距离最小。代码实现分两步先算场景两两之间的欧氏距离矩阵再迭代挑选代表性场景。function keepIdx reduce_scenario(scen, nKeep) nScen size(scen, 1); % 距离矩阵计算scen行为场景列为时段 delta zeros(nScen, nScen); for i 1:nScen delta(i, :) sqrt(sum((scen - scen(i,:)).^2, 2)); end keepIdx zeros(1, nKeep); cand 1:nScen; % 先选距离其他场景最近的那个作为第一个保留场景 totalDist sum(delta, 2); [~, first] min(totalDist); keepIdx(1) first; % 后续迭代每次选择使削减误差增量最小的场景加入保留集 for j 2:nKeep remain setdiff(cand, keepIdx(1:j-1)); gain zeros(1, length(remain)); for k 1:length(remain) gain(k) sum(min(delta(remain, remain(k)), ... min(delta(remain, keepIdx(1:j-1)), [], 2))); end [~, best] min(gain); keepIdx(j) remain(best); end end这段代码对中小规模场景足够用。场景量大时超过5000双重循环会很慢可以先用聚类粗筛再对类内做快速前向。场景削减的核心经验是削减不是越少越好。10个场景基本能保留分布特性5个以下时极端场景容易被削掉调度结果会过度乐观——这一点后面的常见问题部分会展开讲。3.4 求解器选型与建模方式模型是带整数变量的随机规划问题属于混合整数二次规划MIQP。我在项目里用了YALMIP作为建模语言求解器可选CPLEX或Gurobi。没有商用求解器时也可以退而求其次用Matlab内置的ga粒子群等启发式算法求解近似解但精度和稳定性会差一些。YALMIP建模的核心代码结构如下% 决策变量机组出力P、启停u、各场景风电消纳Pw、弃风、切负荷 P sdpvar(nGen, T); u binvar(nGen, T); Pw sdpvar(nScen, T); ... % 目标函数燃料成本 启停 期望惩罚 Objective sum(sum(gencost(P, u))) ... sum(sum(startcost(u))) ... sum(rho .* (pen_curtail * sum(curtail, 2) ... pen_shed * sum(shed, 2))); % 约束 Cons []; Cons [Cons, sum(P, 1) sum(rho .* Pw, 1) Load - sum(shed, 1)]; % 各时段功率平衡 Cons [Cons, Pmin P Pmax]; Cons [Cons, -RampD diff(P, 1, 2) RampU]; % 爬坡约束 ... sol optimize(Cons, Objective, sdpsettings(solver, cplex));注意这里功率平衡约束里的风电项用sum(rho .* Pw, 1)表示对所有场景的期望平衡。如果想严格要求每个场景都满足平衡则需要把约束写成对每个场景逐一展开的形式。实际中我采用的是场景独立平衡这样第二阶段可以理解为针对每个场景重新调度常规机组出力更贴近两阶段随机规划的语义。4. 算例测试与结果分析4.1 测试系统与参数设置算例用的是经典的10机组系统调度周期24小时峰荷1200 MW在此基础上接入一座200 MW的风电场。机组参数采用文献通用数据这里列几个代表性机组的参数机组Pmin (MW)Pmax (MW)爬坡率 (MW/h)a (元/MW²h)b (元/MWh)c (元/h)G11504701200.0004321.60958.2G21354601200.0006321.051313.6G373340900.0003917.261168.7G725160600.0011316.19641.9G1055550022.22550.0风电场景参数切入风速3 m/s额定风速12 m/s切出风速25 m/s额定功率200 MW。Weibull分布取k2.0、c8.5抽样1000个场景削减至10个。备用要求设为预测负荷的5%加风电出力的10%。4.2 确定性模型与随机模型对比对照实验我设置了三组纯确定性模型风电取预测值、随机模型备用取低值、随机模型备用取高值。结果如下表模型期望总成本万元弃风率失负荷概率确定性调度72.43.1%4.2%随机调度低备用74.10.8%0.9%随机调度高备用75.80.3%0.1%确定性模型成本最低但这是虚假的低因为它没有为风电不确定性买单。它在极端场景下切负荷概率达到4.2%这在电力系统运行中是很难接受的。随机调度模型期望成本只增加2%~5%却把失负荷概率压到了1%以内。这就是随机性建模的价值用可控的成本增量换来可靠性的本质提升。4.3 调度结果怎么解读跑完solve_and_plot.m后会得到三张图机组出力堆叠图、场景风电出力曲线、各时段旋转备用柱状图。看结果时重点看两处一是高峰时段晚19-21点机组是否已经逼近上限二是风电出力低谷时段常规机组是否留足上调空间。我测试的多组算例里随机调度方案在风电出力预估较高的凌晨时段会刻意少开一台小机组并保留快速启停机组待命而不是像确定性模型那样为了看起来便宜把机组全开。这个差异在出力曲线上体现得非常明显也是随机调度结果与确定性结果的最直观区别。5. 常见问题与排查技巧实录5.1 求解返回Infeasible problem这是最多人踩的坑。模型不可行的原因我遇到的主要有以下几类爬坡约束与功率平衡冲突某时段负荷陡增超过所有在线机组爬坡之和。解决方法是把备用要求分段线性化或者加入切负荷惩罚项并放宽硬约束。机组最小启停时间约束过紧10机组系统中启停时间约束和爬坡约束叠加容易造成整数变量无解。建议先把启停时间约束去掉跑通再逐步加回来。风电场景出现负值或超额定值场景生成时线性插值超出[0, Prated]范围导致平衡约束无法满足。生成后务必加Pw min(max(Pw,0),Prated)钳位处理。排查技巧把YALMIP的debug选项打开或者把约束逐个注释掉做二分定位。我习惯先跑确定性版本确认基本数据没问题再引入随机场景——这样能快速把数据错误和随机建模错误分开。5.2 场景削减后结果不稳定有段时间我把场景削减到4个结果目标函数值在不同随机种子下波动明显原因就是极端场景被削减掉了。风电场景削减要保留尾部尤其是高风速时段和低风速时段的对冲关系。判断场景质量有个土办法对比削减前后场景集合的预期风电总量和各时段方差。如果削减后某个时段的方差缩水超过30%说明场景太多了或者削减算法参数有问题。我实测下来10-15个场景是精度和速度的甜点区。5.3 求解时间太长10机组、24时段、10场景的MIQPCPLEX一般几十秒能解完。如果你跑几个小时都出不来优先检查这些目标函数是否引入了非凸项。二次成本函数如果是凸的a0YALMIP会自动转成二阶锥如果a设置成0或负值问题变非凸求解器会退化为外逼近法速度骤降。启停成本是否用了非整数条件表达式。建议用implies()或者直接线性化避免引入大量分段变量。备用约束里是否使用了min函数。min在约束里是非线性表达式要改写成线性辅助变量加不等式。另外CPLEX的并行线程数和MIP相对间隙都值得调一下。我常用设置是sdpsettings(cplex.mip.tolerances.mipgap, 0.001, cplex.parallel.cores, 4)求解时间能缩短一半以上且几乎不影响精度。5.4 工具箱和运行环境问题Matlab版本不同可能导致YALMIP和CPLEX接口不兼容最常见的报错是Unable to find solver或者Solver not applicable。确认一下yalmiptest命令能列出CPLEX即可。如果报的是许可证错误通常和当前用户权限或环境变量有关重装求解器前先试optimize时加sdpsettings(verbose, 2)看具体错误信息。另外强烈建议在项目根目录加一行startup.m自动把YALMIP和求解器的库路径加进去避免每次换机器都要手动addpath。6. 个人实操中的几点习惯项目做完整套流程我总结了几条自己一直沿用的习惯。先确定性后随机性。任何随机调度模型一定先把确定性版本跑通、结果符合基本直觉再叠加场景和概率约束。确定性版本是排错的最佳参照物我几乎所有的模型bug都是在两步切换之间暴露的。场景数据要可视化检查。随手plot一下生成的风电场景看看有没有突变、负值、顶格不动的现象。数据质量不过关后面一切模型优化都是白费。我见过太多人只盯着求解器报错实际上问题早在场景生成阶段就埋下了。注释写清楚单位和时段索引。Matlab默认索引从1开始和论文公式里从0或从1都可能对不上。建议所有矩阵的维度定义都在文件头写一行注释比如% P: nGen x nPeriod, 单位MW。保留一个对照实验脚本。永远留一份确定性调度代码做随机调度时随时对比成本和失负荷概率衡量随机化带来的收益到底值不值。写论文时这份对照也是必不可少的分析内容。这套模型后续要扩展也方便比如多风电场相关性建模、储能参与调度、碳捕集机组联合优化都是在现有框架上加约束和加模块的事。先把随机动态经济调度这条主线走通后面的扩展都是水到渠成。
