我读研刚接触电力系统调度那会儿最头疼的其实不是最优潮流和机组组合模型本身而是“源荷不确定性”这五个字。负荷明明是个时变曲线风电光伏更是一阵风一阵云可教科书里的调度模型却总把它们当成固定值去优化。后来自己动手用Matlab做实战才意识到真正难的不是公式推导而是怎么把不确定性装进代码里让调度结果在真实场景里经得起检验。这篇文章就围绕“电力系统调度之源荷不确定性Matlab实战”展开用一个简化火电加风电系统作为载体把负荷预测误差、风电出力随机性、场景生成与削减、随机经济调度建模、求解与结果分析这整条链路完整拆开讲。适合正在做调度优化课题的研究生、刚入门电力系统方向的技术人员以及所有想用Matlab把论文模型落到可运行代码里的人。1. 源荷不确定性到底指什么Matlab里该怎么描述1.1 负荷不确定性不是一条曲线而是一个分布族传统调度里负荷常被当成确定性的预测曲线比如每小时一个数值所有机组按这个数值去安排出力。但实际运行中负荷预测总存在误差尤其是午后和晚上高峰时段误差更大。这个误差可以看成随机变量实践中最常用的是正态分布[ D(t)D_0(t)\varepsilon_D(t),\quad \varepsilon_D(t)\sim N(0,\sigma_D^2(t)) ]其中 (D_0(t)) 是t时段的负荷预测值(\sigma_D(t)) 是预测标准差单位MW。负荷预测精度随时间尺度不同而不同日前尺度标准差一般在2%到5%左右日内滚动会低一些。但如果每个时段都独立抽样会出现一个很反常识的问题相邻时段的负荷误差完全不相关生出来的场景曲线会像高频噪声一样上下乱跳这在现实中几乎不可能发生。所以做源荷不确定性建模时要考虑误差的时序相关性最简单的做法是AR(1)模型sigma_d 0.03 * D0; % 预测标准差取负荷的3% phi 0.8; % 自相关系数 eps_d zeros(NS, NT); for s 1:NS e_prev 0; for t 1:NT e phi * e_prev sigma_d(t) * randn(1); eps_d(s, t) e; e_prev e; end end D_scn D0 eps_d; % 负荷场景矩阵NS行NT列这里NS是场景数NT是时段数D0是1行NT列的基础负荷。AR(1)模型里的phi决定了相邻时段误差的连续程度phi0.8表示误差会缓慢回转不会突变。实际业务中还有更复杂的指数平滑模型但对研究来说AR(1)已经能抓住主要矛盾。1.2 风电出力不确定性风速分布与功率曲线的叠加风电不确定性通常从风速开始建模。风速常用两参数Weibull分布描述概率密度函数为[ f(v)\frac{k}{c}\left(\frac{v}{c}\right)^{k-1}\exp\left(-\left(\frac{v}{c}\right)^k\right) ]其中c是尺度参数k是形状参数。典型的风电场风速场景c可以取6到10k取2左右。Matlab里一句话就能生成风速样本v wblrnd(c_v, k_v, NS, NT);不过风速本身也有时间相关性不能简单独立抽。和负荷误差一样可以对风速时序做AR(1)处理或者用风速的历史样本叠加上随机扰动。生成风速场景后再通过风机功率曲线转成出力function Pw windPowerCurve(v, v_in, v_r, v_out, P_r) Pw zeros(size(v)); idx1 v v_in v v_r; idx2 v v_r v v_out; Pw(idx1) P_r .* (v(idx1) - v_in) / (v_r - v_in); Pw(idx2) P_r; % 低于切入风速或高于切出风速时出力为0 endv_in是切入风速通常3m/s左右v_r是额定风速12到15m/sv_out是切出风速一般25m/sP_r是风机额定功率。注意风速到功率是非线性映射除非你直接采样风电功率否则不能绕开风速分布。1.3 为什么不能把期望值直接当实际值用很多初学者会先算风电出力期望、负荷期望然后把这两个期望值丢进确定性模型里求最优解。这个做法不是完全不能用但结果往往偏“乐观”。原因很简单期望值只是分布的中心不代表真实场景的极端情况。风电大发的深夜场景和风电很小的白天高峰场景对机组出力和备用的需求完全不同。如果把所有不确定性都压成一个点优化器会觉得“不需要那么多备用”因为期望场景下风电刚好够用。等到实际系统遇到风电出力大幅低于预测时才发现爬坡不够、备用不足只能切负荷。所以源荷不确定性建模的意义就是让调度决策在多种可能场景下都能站得住脚而不是只在期望场景下最优。2. 场景生成与削减把连续概率分布变成有限的离散场景2.1 蒙特卡洛采样与拉丁超立方采样的取舍最简单粗暴的方法是蒙特卡洛采样直接用randn和wblrnd生成大量场景比如NS1000。但纯蒙特卡洛有一个问题样本分布会有随机团簇和稀疏区需要有足够多样本才能稳定逼近原始分布。拉丁超立方采样是改进方案先把每个变量的累计分布均匀分N层再从每层里抽取一个样本这样边缘分布覆盖更均匀在同等样本量下误差更小。Matlab可以用lhsdesign生成均匀拉丁超立方数据再通过逆变换映射到目标分布。对负荷误差这种正态分布做法是U lhsdesign(NS, NT); eps_d norminv(U, 0, 1) .* sigma_d;但要提醒的是变量之间独立抽样同样会丢掉时序相关性。处理思路是先抽样独立扰动再通过排序法或Cholesky分解引入相关性。我的经验是如果只是做课程设计和论文仿真AR(1)之后再套拉丁超立方已经够用如果追求工业级精度才需要上Copula或Iman-Conover方法。2.2 同步回代削减从上千场景压缩到十来个典型场景场景太多会让优化问题规模爆炸。一个24时段的场景每个场景都加一套约束如果机组数量和约束条数再多一点Matlab连着求解器很容易卡死。实战中很少直接优化上千个场景而是先做场景削减。我比较常用的是同步回代削减思路选择并保留少数场景让削减前后的概率分布距离尽可能小。核心步骤是计算所有场景两两之间的欧氏距离或加权距离距离定义可以用时空向量差。每次找距离最近的一对场景删掉其中一个把它的概率叠加到另一个场景上。重复直到场景数达到目标值K。Matlab里可以写一个函数输入完整场景矩阵和对应概率输出削减后的场景索引和权重function [idx_keep, prob_keep] ScenarioReduction(scn, prob, K) NS size(scn, 1); dist zeros(NS, NS); for i 1:NS for j i1:NS d sqrt(sum((scn(i,:) - scn(j,:)).^2)); dist(i,j) d; dist(j,i) d; end end keep 1:NS; while length(keep) K min_d inf; rm_i 0; rm_j 0; for i 1:length(keep) for j i1:length(keep) if dist(keep(i), keep(j)) min_d min_d dist(keep(i), keep(j)); rm_i keep(i); rm_j keep(j); end end end prob(rm_j) prob(rm_j) prob(rm_i); keep(keep rm_i) []; end idx_keep keep; prob_keep prob(keep) / sum(prob(keep)); end这个实现属于教学版本O(N^2)复杂度和三重循环在NS上千时会很慢。数据量大时建议用K-means聚类或者调用现成的Fast Forward Selection工具包。但理解同步回代的“合并概率”思想很重要它保证删除场景不是简单丢弃而是把概率转移到相近场景上。我自己的经验是从1000个场景削减到10个典型风电场景的包络和均值都保留得不错再把削减前和削减后分别丢进调度模型求期望成本误差通常能控制在2%以内。如果误差偏大多半是K取得太少或者距离定义时把负荷和风电放在一起时量纲不统一记得先归一化。2.3 场景质量校验别等跑完优化才发现场景有问题很多人做完场景削减直接进优化结果发现调度结果出现一些很不合理的“尖峰”出力。其实问题往往出在场景质量上而不是优化模型上。建议在生成和削减之后先做三件事第一画场景带。把原始场景和削减后场景叠加在一张图上看曲线包络是否连续、是否存在明显突变。风电场景如果出现相邻时段功率跳变几十MW说明风速相关参数可能没调好。第二看极端场景。找风电出力最低和负荷最高的场景单独画出来确认它们不是数值异常而是真实可信的边界情况。调度模型的备用和爬坡决策很大程度上就是被这些极端场景逼出来的。第三多次更换随机种子重复生成检验场景均值方差。运行几次后观察期望成本是否稳定如果波动很大说明场景数不够或削减比例过高。3. 调度优化模型考虑源荷不确定性的随机经济调度3.1 模型设计两阶段随机规划的基本框架源荷不确定性下的调度问题最自然的数学框架是两阶段随机规划。第一阶段在今天做决策确定各机组基点出力、向上备用和向下备用第二阶段在实际运行中根据风电和负荷场景的发生情况在备用范围内调整机组出力必要时切负荷或弃风。第一阶段变量不依赖场景第二阶段变量依赖场景。这个区分非常关键。如果所有变量都依赖场景那就等于对所有场景各算一套独立调度失去了“提前安排备用”的意义。反过来如果所有变量都不依赖场景又过于保守永远无法应对实际波动。本篇文章使用的简化系统包括3台火电机组、1个风电场、1个集中负荷点忽略电网网络潮流约束只保留功率平衡和机组运行约束。这适合入门跑通之后再往模型里加直流潮流或者安全约束。3.2 目标函数与约束条件怎么写目标函数由三部分组成燃料成本、备用成本、期望失负荷和弃风惩罚。燃料成本按基点出力计算备用成本按预留的备用容量计算失负荷和弃风则对所有场景加权期望后加入目标[ \min \sum_{t}\sum_{i}\left(a_i p_{0,i,t}^2 b_i p_{0,i,t} c_i\right) \sum_{t}\sum_{i}\left(c_{ru,i}r_{u,i,t}c_{rd,i}r_{d,i,t}\right) \sum_{s}\pi_s\sum_{t}\left(VOLL\cdot L_{s,t}^{shed}WC\cdot P_{s,t}^{curt}\right) ]约束中功率平衡是每个场景每个时段都要满足的。对于场景s、时段t有[ \sum_i p_{i,s,t}P_{w,s,t}-P_{s,t}^{curt}L_{s,t}^{shed}D_{s,t} ]这个约束把风电出力、弃风、失负荷、系统负荷全部放在一个等式里。实际出力需要在第一阶段预留的备用范围内调节[ -r_{d,i,t}\le p_{i,s,t}-p_{0,i,t}\le r_{u,i,t} ]此外还有机组出力上下限约束、备用容量上下限约束、爬坡约束。爬坡约束必须针对实时场景出力加这样才能体现“遇到风电骤降时机组能不能在1小时内追上来”。实践中很多论文只在基准点上加爬坡约束最后算出来的场景调整量会严重违反物理规律。3.3 YALMIP建模核心代码解析Matlab里做优化建模我一般直接上YALMIP。它能把变量、约束、目标写得接近数学表达式不用自己去拉稀疏矩阵。完整脚本太长这里给关键部分组成。先初始化变量p0 sdpvar(ng, NT, full); % 基点出力 ru sdpvar(ng, NT, full); % 向上备用 rd sdpvar(ng, NT, full); % 向下备用 p_s sdpvar(ng, NS, NT, full); % 场景s下机组实际出力 lshed sdpvar(NS, NT, full); % 失负荷 pcur sdpvar(NS, NT, full); % 弃风然后搭约束关键部分是场景循环里的功率平衡和备用约束Constraints []; for s 1:NS for t 1:NT Constraints [Constraints, sum(p_s(:,s,t)) Pw(s,t) - pcur(s,t) lshed(s,t) D(s,t)]; for i 1:ng Constraints [Constraints, Pmin(i) p_s(i,s,t) Pmax(i)]; Constraints [Constraints, -rd(i,t) p_s(i,s,t) - p0(i,t) ru(i,t)]; Constraints [Constraints, 0 ru(i,t) Ru_max(i), 0 rd(i,t) Rd_max(i)]; end if t 1 for i 1:ng Constraints [Constraints, -Ramp(i) p_s(i,s,t) - p_s(i,s,t-1) Ramp(i)]; end end end end Constraints [Constraints, lshed 0, pcur 0];目标函数这样写Objective sum(sum(repmat(a,1,NT).*p0.^2 repmat(b,1,NT).*p0 repmat(c,1,NT))); Objective Objective sum(sum(c_ru .* ru c_rd .* rd)); for s 1:NS Objective Objective prob(s) * (VOLL * sum(lshed(s,:)) WC * sum(pcur(s,:))); end求解时用ops sdpsettings(solver, gurobi, verbose, 1); sol optimize(Constraints, Objective, ops);a、b、c是燃料成本二次项、一次项和常数项的列向量c_ru、c_rd是备用成本矩阵。没有Gurobi时YALMIP会自动尝试本机的quadprog但模型规模稍大后还是推荐装一个商用求解器速度差距非常明显。3.4 为什么备用成本不能省如果目标里只留燃料成本不考虑备用成本优化器会倾向于把机组出力压到很贴近上限这样一旦负荷增加或风电变小就没有上调空间。实际调度中备用的价值要体现在目标函数里这样才能在“多留备用多花钱”和“不够备用切负荷”之间做合理权衡。失负荷惩罚系数VOLLValue of Lost Load设得越高系统越愿意多留备用。这个参数带有很强的政策属性工程上一般取几千到上万美元每MWh我自己在仿真里常用5000代表极端事件下停电代价远高于发电成本。4. 完整案例流程从基础数据到调度结果一次跑通4.1 案例系统参数与数据准备为了复现方便我给出一个3机1风电场的小系统。3台火电参数如下表机组Pmin(MW)Pmax(MW)a($/MW²h)b($/MWh)c($/h)爬坡(MW/h)G11003000.005830060G2502000.0101020040G3501500.0151215030风电场额定装机200MW切入风速v_in3m/s额定风速v_r12m/s切出风速v_out25m/s。负荷预测曲线设为典型白天两峰一谷8点爬升12点小高峰19点晚高峰凌晨低谷300MW左右最高660MW。风电场出力场景采用NS1000生成再削减到K10个场景。4.2 主流程代码框架整个脚本按四步组织rng(20240101); % 固定随机种子保证可复现 % Step1 生成场景 [D_scn, Pw_scn] GenScenarios(NS, NT, D0, windPara); % Step2 场景削减 [D_red, Pw_red, prob_red] ScenarioReduction(D_scn, Pw_scn, K); % Step3 建立调度模型并求解 [p0_opt, ru_opt, rd_opt, detail] SolveStochasticDispatch(D_red, Pw_red, prob_red, unitPara); % Step4 结果可视化 PlotDispatchResult(p0_opt, ru_opt, rd_opt, D_red, Pw_red);GenScenarios负责生成带相关性的负荷误差与风速时序场景ScenarioReduction返回削减后的代表场景和概率SolveStochasticDispatch内部用YALMIP建模并求解。真正工程化项目里建议把这些步骤拆成函数而不是放在一个脚本里方便后期换数据、换模型、换求解器。4.3 典型结果长什么样以这套参数跑出来的结果几个特征非常典型。第一凌晨风电大发时段G3这类小机组会减少出力甚至接近下限给风电腾空间但向下的备用不能留太少否则风电突然变大时没法消纳。第二傍晚负荷爬升又快又猛G1作为爬坡能力最强的机组会被安排带一部分不是最优经济点的出力换取的代价是备用空间充足。第三削减后的场景里总会有个别场景风电出力很低这时lshed会有非零值但期望失负荷量被VOLL限制住了不会出现大面积切负荷。我习惯把结果画成三张图第一张是负荷和风电场景带第二张是各机组基点出力与备用堆叠面积图第三张是失负荷和弃风期望时序柱状图。这三张图基本能满足论文和汇报的展示需求。如果进一步做敏感性分析可以把VOLL从1000改到10000观察备用总量和失负荷期望的变化曲线这就是稍加扩展就能出成果的方向。5. 拿Matlab做这类问题最容易踩的坑和自查清单5.1 场景数和随机种子的“玄学”有一段时间我跑同一套代码两次结果差距大得离谱后来发现是随机数状态没固定。生成场景前写上rng(整数)是基本操作但更严谨的做法是跑多组种子看目标函数期望值是否稳定。一组种子只能说明“在你抽到的这些场景下结果是这样”不能说明系统整体风险。把种子从1跑到20统计目标函数均值和方差比单独跑一次1000场景更有参考意义。此外场景削减后的K值不是越大越好。K太小会漏掉极端情况K太大会让优化模型变量爆炸。以我的经验24时段3机系统K取10到20就够如果系统节点数很多K通常取5到10配合备用约束已经能反映主要不确定性。5.2 约束拼接速度慢到怀疑人生YALMIP的约束Constraints[Constraints, newConstraint]在循环里写起来很自然但在场景数和时段数变大后这种逐条拼接会让Matlab慢得令人抓狂。原因是每次拼接都要复制整个约束对象。优化方法是先一次性构建所有约束放进cell数组最后再合并consCell cell(NS*NT, 1); idx 0; for s 1:NS for t 1:NT idx idx 1; consCell{idx} [...]; end end Constraints [consCell{:}];另一个经验是把所有的p0、ru、rd当成矩阵变量处理避免循环内单个元素的约束。能用矩阵运算是Matlab的天然优势代码可读性和求解器构建效率都会高很多。5.3 丢失时序相关性的场景一堆废数据很多初学者只对每个时段独立采样完全不管时序生成的风电场景白天黑夜一个样或者负荷场景抖动特别剧烈。把这种场景投入调度模型爬坡约束会把机组折腾得极其保守备用容量无脑拉高目标成本虚高。这个问题肉眼能看出来画一张场景带图如果曲线像雪花一样乱颤说明相关性建模有问题。修正方法就是前文提到的AR(1)或指数平滑。对风电来说如果风速序列已经有历史数据更推荐用历史样本加小扰动的办法生成场景保留原始相关性。5.4 结果合理性检查不要只盯着目标函数值目标函数值降下来了不代表调度可行。模型里约束写漏了一个可能直接导致机组爬坡不够或功率不平衡。在出结果后我至少会做三层检查。第一层逐时段检查功率平衡。把最优解代回平衡等式看残差是不是在1e-6以下。这个检查最简单也最能发现维度错误。第二层检查爬坡量。把每个机组相邻时段出力相减看有没有超出爬坡约束的隐性越界。第三层检查备用利用率。如果某机组向上备用全程为0说明这个备用容量根本没参与决策大概率是约束没激活需要回头看约束关系。如果做了这三层检查仍然觉得不踏实还可以把随机调度结果和确定性调度做对比。确定性调度模型把风电设成期望值、负荷设成预测值求出来的目标函数值通常更低但把所有场景套上去做滚动校验时会发现失负荷期望很高。这个对比正好解释了为什么必须做不确定性优化。最后再补两句我自己绕过的弯路做这个案例时我在场景削减上耽误的时间比建模还多。最初图省事直接用K-means聚类结果聚类中心是“平均形态”极端的风电低谷场景被完全吸收掉了调度结果总在个别时段失负荷。后来换成同步回代并保留每个极端场景结果立刻合理很多。还有一件事值得说。YALMIP配合Gurobi求解时如果目标函数是二次的默认会走QP或MIQP路径模型规模变大后内存占用会涨得很快。我的习惯是先跑一个线性成本版本把整体流程调通确认场景和约束都没问题再切到二次成本。这样出错时你能很快判断是模型问题还是求解器问题。源荷不确定性调度是一个越做越深的主题这篇内容覆盖的只是从不确定性建模到随机调度的完整闭环后面你可以继续把网络潮流约束、机组组合0/1变量、储能和需求响应逐步加进去每一步扩展都会让模型离真实系统更近。