风电占比一上来调度问题就变味了。以前做经济调度看好负荷预测排一遍机组启停计划火电把尖峰扛住基本就稳了现在风电出力跟闹情绪似的早高峰来得出力小、晚高峰风大得往外送等到了大风天预测功率和实际出力差个15%—20%都是家常便饭。这种背景下“风电随机性动态经济调度”成了电力系统优化里绕不开的研究方向。它要解决的问题用一句话说就是在风电预测带有明显不确定性的前提下怎样提前安排好火电机组的启停状态、出力基准和旋转备用让系统运行成本尽量低的同时又能扛住各种可能出现的“预测失准”场景。这篇文章就是围绕这个方向给出一套完整可跑的Matlab代码实现框架。整套代码覆盖了风电场景生成、场景削减、YALMIP下单目标随机动态经济调度的建模、求解和结果可视化不依赖任何复杂的商业平台有Matlab和合适的求解器就能跑。适合刚接触随机优化的研究生做算法复现或论文对比也适合电力系统运行人员快速了解随机调度模型到底在解什么题、边界在哪。下面我按建模思路、代码实现、实验设计、调试经验四个维度拆开讲。1. 动态经济调度到底在解什么题——从确定性到随机性1.1 确定性调度的痛点把预测值当成真值用传统确定性经济调度的套路很直接把系统负荷视为已知曲线把风电出力固定在预测值上然后让所有机组在满足功率平衡、出力上下限和备用要求的前提下找到一个总发电成本最低的出力方案。模型在数学上大多是二次成本函数的优化用Matlab调用求解器几秒就能算完这也是教材里最基础的形式。问题是风电介入后预测值这个假设本身就站不住脚。我实测过一个风电渗透率在20%左右的系统单个时段风电预测误差可以达到装机容量的两成以上一天下来最大偏差甚至能到四成。你如果按预测值做调度实际运行时风电出力一旦偏小缺口就得靠火电快速顶上但火电有爬坡速度限制不是你说加出力就能立刻加。更麻烦的是这种缺额在一天里往往会持续好几个小时等备用用尽最终的代价可能就是切负荷或紧急拉闸限电。确定性模型最大的问题是它把不确定性当成不存在所有结论都是“拿预测值对正确答案”然而真实电网运行时根本不存在一个唯一的“正确答案”。1.2 “动态”二字的真正含义跨时段耦合约束为什么是“动态”经济调度而不直接叫经济调度因为实际运行里火电机组不是独立的单时段决策机组从一个出力点调整到另一个出力点需要时间这就是爬坡约束机组开起来之后不能立刻关关了也不能马上再开这对应最小开停机时间约束每次启停还要额外付出一笔不小的启停费用。这些约束把24个小时的决策串成了一条连续链条任何单时段解出来的最优方案放在整个时间尺度上可能完全不可行。举个例子凌晨负荷低谷时可能最优解是关掉两台小机组省燃料费用但如果白天负荷猛增时又需要这两台机组立刻顶上而它们的最小停机时间还没到那就只能眼睁睁看着缺功率。所以动态调度模型里必须同时处理机组状态变量开机/停机0-1变量和连续出力变量整张问题变成了混合整数二次规划。这个“动态”三个字看着不起眼却在建模复杂度和求解难度上都带来了数量级的提升。1.3 引入随机性后的两难成本与风险的权衡把风电随机性引入之后问题就变成了一个两难权衡如果调度方案设计得特别保守把旋转备用设置得很高那么成本一定上去严重的时候一天下来多烧出十几万的油钱和煤钱如果完全按期望值来做不预留裕度看起来经济性最优但真实运行中只要预测稍微偏一点就可能触发代价极高的切负荷事件。系统运行要做的不是寻找一个绝对最优解而是在“预期运行成本”和“极端场景下的失负荷风险”之间找到一个平衡点。这里面有个工程上特别容易理解的关系我常跟同行打比方做确定性调度就像开车只看路况预测随机调度则是不仅要看预测还要看这条路上有多少种可能的大雨、滑坡、堵车情况然后决定要不要绕路或者多留一点车距。只不过在电力系统里这个“车距”就是旋转备用容量和机组组合方案的鲁棒性。后面的建模就是围绕这件事展开的。2. 风电随机性建模三种处理思路与选型判断2.1 场景法用多个场景逼近未来的不确定性场景法是目前学术界和工程应用里最常用的随机性建模手段。核心思路很直观风电预测结果其实是一个期望值真实的出力可以看作期望值加上一个随机误差我们把这个误差分布进行大量采样得到成百上千个“如果明天是这种天气、这种风速风电会这样出力”的样本场景然后把这些场景作为模型的输入条件让调度方案在所有场景下都尽量可行且成本期望值最小。实际实现时有三个环节要注意。第一个是误差分布的选择很多论文直接假设预测误差服从正态分布实测下来风电功率预测误差普遍有偏峰厚尾特性一小时级以上的误差分布尾部比正态要厚得多所以有条件的话用历史预测数据和实际出力数据去拟合一个混合分布或直接采用拉普拉斯分布去近似效果会好不少。第二个是场景生成时要注意时间相关性不能每一时段独立采样否则生成的场景抖动太剧烈专业术语叫“场景不连续”做出来的备用需求失真。第三个是场景数量生成1000个场景让模型能“看到”足够多的可能性但求解时1000个场景会让问题规模暴涨所以必须先做场景削减。2.2 机会约束在概率约束下做权衡场景法是用有限个离散场景去近似连续分布的随机性机会约束则是直接在约束层面做成“大概率满足”的形式。比如系统旋转备用约束不必在每个可能的风电出力场景下严格成立而是允许以一定概率比如95%的置信度成立把约束写成P{R(t) 需要的备用} 0.95。这样模型不会为了那极端的5%场景过度配置备用经济性更好。这个思路在处理备用约束时特别实用。做法是把机会约束线性化假设系统净负荷预测误差大致服从正态分布那么给定置信度α需求对应的备用容量就等于期望值加上z_{α}乘以标准差。比如置信度取95%z就是1.645模型里直接加上这个确定性备用即可。如果误差分布不是正态可以用分位数法做更精确的处理。机会约束的优点是计算负担小缺点是它本质是做了解析近似对分布假设比较敏感极端场景下的失效概率需要额外校验。2.3 鲁棒优化对抗最坏场景鲁棒优化的思路完全不一样它不看概率分布而是给定一个不确定集把风电出力限制在这个集合里变化然后要求调度方案在集合里最不利的那个场景下也要满足约束。工程上常用的是盒式不确定集比如风电出力在预测值正负20%范围内波动模型就为这个最坏情况准备备用。鲁棒模型的求解难度比场景法低对误差分布的假设依赖也少是工程实用中相当讨喜的选择。但鲁棒优化有个绕不开的天性保守。纯盒式鲁棒模型守着最极端场景的结果往往比场景法算出来的成本高不少在风电占比较高的时候尤其明显。字符上更精细的做法是引入不确定性预算budget of uncertainty限制同时出现最坏场景的时段数量用参数Γ控制保守程度这样模型在极端保守和过度乐观之间通过一个参数连续调节工程上非常灵活。2.4 工程选型判断该用哪一种做实际项目时怎么选我给大家一个简单的判断规则如果手上只有预测曲线和典型误差统计数据优先用场景法因为它直观、可解释性强审查起来也容易如果调度对安全要求极高不允许出现备用不足的情况且计算资源有限用鲁棒优化更稳妥如果计算时间有限、模型规模大、又不需要特别精细的分布建模就用机会约束处理备用约束其他部分保留确定性混合建模。三者的对比我整理成一个表格方便直接对照。建模方式数据需求计算复杂度保守程度适用场景场景法需要误差分布或历史场景数据高场景数×问题规模中等学术研究、方案比选、多场景分析机会约束需要误差均值和方差分位数低较低备用约束、运行备用配置鲁棒优化只需要不确定集参数低高可调Γ保安全场景、工程实用3. 风电随机性动态经济调度Matlab实现全流程3.1 代码整体模块设计写代码之前先把模块想清楚。我的习惯是把这个工程拆成五个主模块数据输入模块、场景生成模块、模型构建模块、求解模块、结果后处理模块。这样做有个好处你后面想换数据、换场景数、换求解器只需要改对应模块的接口不用把整个脚本推倒重来。模块主要文件核心内容数据输入case_data.m机组参数、负荷曲线、风电预测数据场景生成scenario_generate.m风电误差采样生成初始场景场景削减scenario_reduction.m场景聚合并输出各场景概率模型构建build_model.mYALMIP变量定义、约束、目标函数求解与后处理run_dispatch.m / plot_result.m求解、成本统计、画图3.2 场景生成与削减数据准备环节场景生成模块是整套随机模型的输入端。下面的代码实现了一个基于正态误差扰动的场景生成函数实际使用时误差标准差建议按预测出力的比例取比如取15%左右并根据历史误差数据校准。生成场景后注意把负数钳位到零因为风电机组出力物理上不可能为负。function Scen scenario_generate(wind_fcst, nScen, err_std) % wind_fcst: T x N 矩阵T时段数N风电场数 % nScen: 采样场景数 % err_std: T x N 误差标准差与风电场预测同量纲 [T, N] size(wind_fcst); Scen zeros(T, N, nScen); for k 1:nScen for t 1:T Scen(t, :, k) wind_fcst(t, :) err_std(t, :) .* randn(1, N); end Scen(:, :, k) max(Scen(:, :, k), 0); end end生成1000个初始场景之后直接拿到优化模型里跑问题规模会非常夸张。这里推荐先做场景削减常用方法是快速前向选择Fast Forward Selection每次从场景集合里找出“与其余场景概率距离之和最小”的场景把它删掉同时把它的概率并给距离最近的保留场景直到场景数降到预设值。Matlab自带统计工具箱的kmeans也能做简化版聚类削减虽然削减效果不如FFS精细但胜在快速。function [ScenRed, ProbRed] scenario_reduction(ScenRaw, nKeep) % 用kmeans聚类做简化场景削减 % ScenRaw: T*N x nScen每列是一个场景 [nDim, nScen] size(ScenRaw); [idx, ~] kmeans(ScenRaw, nKeep, Replicates, 5); ScenRed zeros(nDim, nKeep); ProbRed zeros(nKeep, 1); for k 1:nKeep members idx k; ScenRed(:, k) mean(ScenRaw(:, members), 2); ProbRed(k) sum(members) / nScen; end end这个函数返回每个保留场景的出力曲线和对应概率后面建模时直接用期望累计。实测中把1000个场景削减到20个计算精度损失通常在可接受范围内但求解时间能从几十分钟降到几分钟性价比非常高。3.3 YALMIP建模目标函数与约束的代码化模型构建是整个代码库的核心。系统目标函数包含四部分火电燃料成本、机组启停成本、备用成本、弃风惩罚。燃料成本是P的二次函数启用0-1变量表示机组开停机状态爬坡约束把相邻时段的出力关联起来旋转备用约束负责兜底。用YALMIP声明优化变量的好处是可以直接用sdpvar和binvar两种类型模型读起来跟数学表达式基本一一对应。T 24; nG 10; nS 20; % 时段、机组数、场景数 P sdpvar(nG, T, nS, full); % 机组出力三维机组x时段x场景 U binvar(nG, T, full); % 机组启停状态所有场景共用 SU sdpvar(nG, T, full); % 启停成本变量 Pw sdpvar(nW, T, full); % 实际风电出力弃风后可调度注意这里U的维度是nG×T所有场景共用一套启停计划这才是“动态经济调度”里机组组合决策的真实约束在实际运行中调度的启停计划只能定一套不能每个场景各定各的。目标函数里的场景期望成本用概率加权实现objective 0; cost_fuel 0; cost_start 0; cost_wind 0; cost_res 0; for s 1:nS for g 1:nG for t 1:T cost_fuel cost_fuel ProbRed(s) * (a(g)*P(g,t,s)^2 b(g)*P(g,t,s) c(g)*U(g,t)); cost_start cost_start ProbRed(s) * suc(g) * max(0, U(g,t)-U(g,t-1)); end end % 弃风惩罚 for w 1:nW for t 1:T cost_wind cost_wind ProbRed(s) * lambda_w * (wind_fcst(t,w) - Pw(w,t,s)); end end end objective cost_fuel cost_start cost_wind;功率平衡约束写起来也比较直接。注意场景s下的负荷、机组出力、风电出力要满足功率平衡同时风电出力受限于该场景下生成的预测场景最大值Constraints []; for t 1:T for s 1:nS Constraints [Constraints, sum(P(:,t,s)) sum(Pw(:,t,s)) Load(t)]; for w 1:nW Constraints [Constraints, 0 Pw(w,t,s) ScenRed(t,w,s)]; end end end爬坡约束要跨时段关联这里特别容易写错索引。爬坡率的含义是机组在前后一小时内的出力变化上限所以在同一个场景内t时段和t-1时段的出力差值要同时小于向上爬坡速率和向下爬坡速率for t 2:T for s 1:nS Constraints [Constraints, P(:,t,s) - P(:,t-1,s) R_up(:)]; Constraints [Constraints, P(:,t-1,s) - P(:,t,s) R_dn(:)]; end end这里我把R_up和R_dn定义成nG×1的列向量用冒号展开保证维度对齐。如果你把R_up定义成1×nG的行向量就会在这里踩维度不一致的报错坑。3.4 结果输出该画哪些图、算哪些指标模型运行完成后后处理模块要输出三类东西一是机组组合方案各时段哪些机组开着、哪些关着二是机组出力堆叠图三是关键指标统计。我实测下来最有用的可视化是“负荷-出力堆叠图”用Matlab的area函数把各机组出力按从下到上的顺序堆叠起来顶端叠上风电出力旁边叠加负荷曲线一眼就能看出系统供需是否匹配、备用余量大概多少。figure; hold on; area(1:T, P_base(1,:)); % 第一台机组出力 for g 2:nG area(1:T, P_base(g,:), FaceAlpha, 0.8); end plot(1:T, Load, k-, LineWidth, 2); legend([cellstr(num2str((1:nG))), 负荷]); xlabel(时段); ylabel(功率/MW); grid on;指标方面重点统计总运行成本、弃风率、备用满足率、以及各场景下的失负荷频次。备用满足率是随机模型特有的指标把所有削减后场景逐一回代到调度方案里检查旋转备用约束是否有场景不满足统计满足约束的场景概率权重之和。这个指标能直观反映模型方案的可靠性工程汇报时也特别有说服力。4. 实验设计与对比验证我们关注什么指标4.1 测试系统参数一套可复现的10机组系统为了让代码有对比基准我这里用一套经典的10机组系统作为测试算例24时段负荷取某个典型日曲线两座风电场总装机容量分别占系统负荷峰值的10%左右。机组参数表如下成本系数a的单位是$/MW^2hb是$/MWhc是$/h爬坡率单位是MW/h实际用的时候注意单位对齐。机组Pmax/MWPmin/MWabcR_upR_dn最小开机时间/h最小停机时间/hG14551500.0004816.19100016016088G24551500.0003117.2697013013088G3130200.0012016.60700606055G4130200.0011016.50680606055G580250.0020019.704504040334.2 对比实验确定性模型 vs 随机模型在同样的系统参数和负荷数据下我做了三组实验。第一组是确定性模型把风电固定在预测曲线上不加任何随机场景第二组是场景法随机模型把风电预测误差按15%标准差生成1000个场景并削减到20个第三组是机会约束模型只用95%置信度写备用约束。三类模型都跑出机组组合方案和总成本。结果非常有意思。确定性模型的总成本看起来最低但它实际上严重“透支了未来的风险预算”回代真实场景时备用不足的概率高达三成这意味着这个方案在实操中几乎不可用。随机模型的总成本比确定性模型高出4%—7%这在工程上是完全可以接受的代价换来的是备用满足率提升到95%以上。机会约束模型虽然成本也比确定性模型高但比场景法低不少它的备用满足率刚好卡在95%附近算是一种“性价比”策略。模型类型总运行成本/万元备用满足率弃风率求解时间/s确定性模型72.568%8.2%3场景法随机模型76.897.4%4.1%260机会约束模型74.995.1%5.6%154.3 关键指标解读成本构成与可靠性这组实验反映了一个核心的trade-off总成本不是衡量随机调度模型唯一的标准还必须看备用满足率。确定性模型“省”下来的几万块钱是拿安全裕度换来的。实际项目汇报中我一般用“单位可靠性成本”这个指标辅助判断计算公式是随机模型成本-确定性模型成本/随机模型备用满足率-确定性模型备用满足率含义是为了提高一个百分点的可靠性系统平均要多付出多少成本。这个比值在不同风电渗透率下有显著差异渗透率越高提高同等可靠性的成本越大这也是为什么高比例风电系统需要花更大的代价去配置储能或灵活性资源。5. 调试实录我在跑这个模型时踩过的几个坑5.1 模型无解先查可行域而不是查求解器YALMIP给“无解”的报错信息相当简略就一个“Problem is infeasible”不告诉你哪个约束导致的。我第一次跑这个模型时就在这卡了很久。实测下来无解大概率是下面几个原因之一。第一个原因是爬坡约束和机组最小启停时间约束同时过紧。比如某台机组最小停机时间是3小时但模型为了让凌晨出力满足爬坡限制强行让它在2小时后开机约束直接打架。这类问题用Gurobi或CPLEX求解时会输出冲突约束的Irreducible Inconsistent SubsystemIIS信息YALMIP在R2023a之后的版本也支持通过export导出模型然后读mps文件分析IIS这是排查无解最快的手段。第二个原因是风电场景里存在极端值。生成场景时如果误差标准差设置过大某几个场景的风电出力会接近零同时负荷又很高备用就完全不够用。处理方法不是删掉极端场景而是在模型里给失负荷设计一个高额惩罚变量允许极端场景切点负荷而不是让整体模型死锁。第三个原因是约束写错索引尤其是爬坡约束里对t1时段的处理。第一时段没有上一时段出力理论上应该用初始运行状态约束或者干脆不写第一时段的爬坡约束。很多人复现论文代码时把爬坡约束从t1写到T第一时段引用U(:,0)直接报错或者用循环变量的时候把维度写歪了这类问题不仔细看根本发现不了。5.2 求解时间爆炸场景数和整数变量是主因随机动态经济调度的求解时间基本上由两个因素决定0-1整数变量的数量、场景数量。拿10机组24时段来说状态变量U就有10×24240个0-1变量每个场景的目标函数和约束都要展开场景数从20涨到50求解时间不是线性增长而是近似指数增长。实测中20个场景要跑4分钟左右50个场景直接奔着半小时去了。工程中如果求解时间超预算优先做三件事场景削减压缩到10个以内、把最小启停约束用线性不等式展开而不是用大M法慢慢迭代、把求解器MIP gap从0.01%放宽到0.5%。这三招组合下来求解时间通常能压到一个数量级以内成本精度损失可以控制在1%以内。5.3 结果不合理的典型排查有一类问题特别有迷惑性模型有解、结果也跑出来了但一看就觉得不对劲。比如弃风率异常高明明模型里目标函数中有弃风惩罚结果还是大片弃风。这时要去查弃风惩罚系数我见过有人把lambda_w设成很小的一个数比如0.001元/MWh这个数值比火电的边际成本低好几个数量级模型当然选择多烧煤、不消纳风电因为弃风的“代价”微不足道。惩罚系数量级必须高于火电边际成本才能在目标函数里真正起到约束作用。另一种情况是机组出力一直顶在上限看起来“满发”很合理但系统同时又存在大量弃风说明备用约束根本没有被激活模型里备用约束可能被写成了软约束或者根本没加。排查时可以先单独跑一个确定性模型逐时段打印功率平衡方程左右两侧的数值看约束残差是否是零如果不是零那必然是有某处约束没写上。5.4 数值病态量纲和标幺化的坑数值问题在配电网和微电网算例里尤其突出。比如微电网里机组功率只有几百千瓦而成本系数按元/MWh量级代入目标函数里就会出现P的平方乘十几万、再乘以1e6的系数数值跨度能到1e12求解器一旦遇到这种病态数收敛速度会慢到让人怀疑人生甚至直接给一个错误的最优解。我处理这类问题的标准做法是先做标幺化把功率基准定成系统峰值负荷比如1000MVA所有功率量纲都除以这个基准对应成本系数也做等比例换算让变量和系数都落在0到几十这个量级求解稳定性立刻好转。还有一个容易被忽略的细节是目标函数里二次项系数(a)的数值如果系统容量特别大a可能小到1e-7这一档建议把所有成本的计算单位统一成“万元”避免结果里出现一大堆科学计数法。6. 后续扩展方向这套模型还能往哪改6.1 加入储能与需求响应储能是最自然的一个扩展。加储能本质上是引入一组新的决策变量储能充放电功率、储能SOC状态变量新增的约束是SOC动态方程、充放电功率限制、以及最终SOC回到初始值。把这组变量和约束写进随机模型的框架里储能会起到非常明显的削峰填谷和备用替代作用。风电渗透率越高的系统储能在随机模型中的价值越突出因为储能能直接对冲风电误差带来的功率波动。代码层面改动也不大在目标函数里加一个储能的损耗成本项约束里加上时间序列递推方程就能得到一个含储能系统的随机经济调度模型。6.2 目标函数中考虑碳成本如果你有双碳课题背景扩展方向就是把碳排放成本直接塞进目标函数。做法不复杂给每台火电机组的燃料成本后面加一项碳价×排放强度×出力大概率就能让高排放机组在组合方案里被排挤掉。这个扩展非常值得做因为很多情况下高碳机组和低碳机组仅仅差在排放强度参数上但加了碳价之后机组组合方案可能有质的改变。实测中碳价设在40—80元/tCO2系统总成本和碳排放量能取得一个相当合理的平衡。6.3 从两阶段到多阶段多时间尺度协调再往深处走可以把单层两阶段随机模型扩展成多时间尺度协调模型。比如日前调度做机组组合日内滚动修正做机组出力调整实时调度再做备用激活。这本质上是一个三层递进的问题各层之间的变量通过场景树或者约束传递联系起来。Matlab里实现多时间尺度模型比较实用的方式是外层用一个循环控制时间窗口滚动内层求解单时段或短时段的优化子问题而不是一次性建模求解不然模型规模完全没法控制。这套代码作为框架的好处就在这里你只需要把build_model函数里对应的约束和目标改一下就能往这些方向延伸。最后再分享一个小技巧做这一整套随机经济调度代码时别把重点全放在“场景生成”上场景削减这一环同样值得花时间调。很多新手一开始就是生成几百个场景直接丢进YALMIP结果模型大得跑不动又回头改场景数其实只要做一个好的削减算法10个精心挑选的场景就能让模型具备很高的代表性计算量却小了一个量级。这算是我在多次复现和改造这套模型之后印象最深的一点希望对正在啃这个题目的你有帮助。