搞电力系统调度的朋友应该都清楚风电并网之后最让人不踏实的就是它那股“说变就变”的出力。天气一变化风机出力从800MW掉到200MW可能只需要十几分钟负荷那边也不是个省心的主儿早晚高峰的预测误差经常能到好几千千瓦。传统的确定性调度把预测值当成真实值用定出来的出力计划常常在实时运行中严重失配轻则弃风重则切负荷。这篇博文要聊的就是怎么把电源侧风电出力的随机性和负荷侧负荷预测误差两方面的不确定性一起塞进低碳调度模型里并给出完整的Matlab实现思路。代码我已经跑通了文章会按“问题建模 → 不确定性场景生成 → Yalmip建模求解 → 结果分析 → 常见坑位”的顺序来写。适合正在做电力系统优化调度课题的研究生、刚接触不确定性建模的工程师还有想把论文模型落到代码上的同学。整个过程不抖理论包袱尽量用能直接抄作业的细节讲透。1. 项目概述与问题建模1.1 为什么必须同时考虑源荷两侧不确定性很多初学者一开始只把风电出力当随机量处理负荷取一个确定预测值然后心安理得地开始写代码。我最早也这么干过结果算完一看备用约束怎么算都觉得别扭。后来想明白了系统调度关心的核心其实是“净负荷”也就是负荷减去风电出力。你光把风电视为随机量负荷却当成确定值相当于只考虑了净负荷波动的一半。风电出力随机性主要来源于风速的间歇性。对同一个风电场来说短期预测误差的标准差通常能到预测出力的10%20%极端天气下更大。负荷侧则是另一回事它受气象、生产生活节律、突发事件影响预测误差普遍在1%3%之间绝对值不小。一个600MW的风电场预测出力300MW实际可能只有几十MW甚至接近零配上一座有2000MW负荷的城市负荷预测偏差2%就是40MW。两边叠加净负荷的波动范围可能超过600MW。如果调度模型里完全没考虑这个波动范围旋转备用和机组组合安排得再漂亮实时一跑大概率会出问题。同时考虑双侧不确定性的意义在于它让备用配置和机组组合结果更贴近真实运行需求。你可以把净负荷当成一个随机变量来处理也可以把风电和负荷分别建模再通过功率平衡方程耦合。实际工程里前一种操作更简洁后一种物理意义更清楚。我写代码时更喜欢分别建模因为这样便于做源荷相关性的灵敏度分析也方便后续扩展储能的充电策略。1.2 低碳调度模型的目标函数与约束体系“低碳调度”不是单纯把碳排放量最小化而是在保证经济性和安全性的前提下把碳排放成本或碳排放上限纳入决策。国内现在做相关课题基本都会结合碳交易机制系统根据配额免费获得一部分碳排放额度实际排放低于配额时可以在市场卖出超出配额时就得掏钱买。这个机制放进调度模型里非常自然碳价就是碳排放的边际成本。目标函数我习惯写成四项之和燃煤成本通常用二次函数表示a·P² b·P c机组启停成本用0-1变量配合启动变量来算碳交易成本等于碳价乘以“实际碳排放量减去配额”弃风惩罚成本等于弃风量乘以惩罚单价。为什么要加弃风惩罚因为模型优化时可能为了躲开机组最低出力限制或者减少启停次数而故意压低风电出力导致大量弃风。加上惩罚项之后系统会在“多烧煤”和“弃风电”之间做个折中算出来的结果更合理。惩罚单价一般取得比煤电边际成本高一些我常用的是50100元/MWh。约束方面功率平衡是硬约束Σ P_fire(i,t) P_wind(t) P_load(t)这个等式在每个场景、每个时段都要满足。另一个关键约束是旋转备用容量我一般同时配正备用和负备用。正备用要求所有在线机组的剩余可调容量之和大于负荷预测偏差的向上分位点负备用则要求能往下调节的容量足够覆盖风电突然增加的场景。在场景法里这些备用约束可以写成对所有场景都成立的形式也可以放宽为置信水平约束。2. 不确定性建模方法选型与场景生成2.1 为什么选场景法而不是鲁棒优化或机会约束处理不确定性主要有三条路随机规划场景法、鲁棒优化、机会约束规划。鲁棒优化的思路是不给概率分布只给一个不确定集让方案对集合内所有可能都可行。优点是不需要分布假设缺点是结果通常非常保守成本高出一大截而且不确定集的边界怎么定本身就是个难题。机会约束规划允许约束在一定程度上失稳比如“功率平衡以95%概率满足”。它直接利用概率分布来压缩保守性但求解时如果碰到非正态分布或非线性约束转换过程非常麻烦。Matlab里做小规模还可以规模一大就卡住。场景法本质是蒙特卡洛思想采样一批可能的风电功率和负荷功率把随机优化问题等价成一个大规模确定性问题来解。它最大的优点就是用起来直观代码实现结构清晰商业求解器能直接处理。论文里面对审稿人也有一定说服力因为你提供了具体的场景集和削减算法。我日常做课题基本首选场景法。2.2 风电出力场景生成拉丁超立方采样风电出力场景的核心是构造“预测出力 随机误差”的样本。误差的分布有很多选择工程里最常用的是正态分布但要注意两个问题一是误差可能在0和装机容量附近截尾否则会采样出负的风电功率或超过装机的值二是风电出力在0附近的概率密度实际上并不表现为正态更准确可以用Beta分布拟合。我个人的常用做法是先用正态分布采样再对越界值做截断处理实操上完全够用。标准差取预测出力的12%18%。如果风电场数量多于一个还要考虑空间相关性我一般用Cholesky分解协方差矩阵来生成相关样本。Matlab里有现成的lhsnorm函数可以生成拉丁超立方多维正态随机数没有统计工具箱的话自己写也不难。% 生成单风电场出力场景含截断 P_wf 300; % 预测出力(MW) P_cap 600; % 装机容量(MW) sigma 0.15 * P_wf; % 误差标准差 n_scene 500; % 采样场景数 e normrnd(0, sigma, n_scene, 1); % 正态误差 P_wscen P_wf e; P_wscen max(P_wscen, 0); P_wscen min(P_wscen, P_cap);拉丁超立方采样比简单蒙特卡洛好在哪里它把每个维度分成等概率的区间在每个区间里只采一个点保证样本在整个分布空间里铺得均匀。我做过对比场景数只有100时LHS得到的优化成本方差就显著小于普通随机采样做灵敏分析时结论更稳定。2.3 负荷场景生成与源荷相关性负荷场景更简单误差一般取预测负荷的1%3%同样用正态分布生成。需要注意的是一天24小时负荷预测误差的标准差并不恒定——早晚高峰预测难误差大凌晨平稳时段误差小。建模时可以按时段设置不同的标准差我试过恒定标准差和分时段两种分时段对结果的影响在备用配置上体现得比较明显。关于源荷相关性严格来说风电出力与负荷之间是有耦合的比如大风降温天气风电大发而负荷飙升但这两者相关性往往弱于风电场的空间相关性。对于绝大多数课程项目和期刊论文独立采样已经能说明问题。如果你立志要做精细化建模可以采用Copula函数构造联合分布但代码量会明显上去求解环节也得多一层转换。我建议先把独立采样的版本跑通再考虑相关性扩展。2.4 场景削减同步回代消除法实现500个场景直接塞进优化模型大多数个人电脑跑不动。场景削减是必须的。同步回代消除法Simultaneous Backward Reduction是最经典的做法思路很直白先给每个场景等概率然后每轮删除一个概率最小的场景把它的概率累加到距离它最近的场景上直到场景数降到目标值。% 同步回代消除法核心逻辑示意 function [s_reduced, p_reduced] sc_reduction(s_in, w_in, K_target) s s_in; w w_in; while length(w) K_target n length(w); % 计算场景两两距离按时间序列欧氏距离 D zeros(n, n); for i 1:n for j i1:n D(i,j) norm(s(i,:) - s(j,:)); D(j,i) D(i,j); end end % 找每个场景的最近邻距离 D(D0) inf; [d_min, idx_nb] min(D, [], 2); % 最小删除代价场景(1-w)*d_min cost (1 - w) .* d_min; [~, del_idx] min(cost); % 概率叠加 w(idx_nb(del_idx)) w(idx_nb(del_idx)) w(del_idx); % 删除场景 s(del_idx,:) []; w(del_idx) []; end s_reduced s; p_reduced w; end削减的要点是每次删除代价最小的场景代价等于该场景概率乘以到最近邻的距离。这个算法的Matlab实现网上有很多版本自己写的时候注意距离矩阵的对称性和对角线置零否则很容易把“自己跟自己最近”选进去。我用下来把500个场景削减到20个优化结果和直接用500个场景相差不到5%但求解时间能缩小一个数量级。3. Matlab代码实现核心逻辑3.1 总体框架与主函数结构完整的代码框架大概长这样数据准备 → 场景生成 → 场景削减 → Yalmip建模 → 求解 → 结果整理与画图。我习惯把每个环节封装成独立函数方便改参数重跑。数据部分包括机组参数、24小时负荷预测曲线、风电预测曲线、碳价、备用系数。%% 主函数骨架 clear; close all; clc; addpath(genpath(yalmip)); % Yalmip工具箱 addpath(genpath(gurobi)); % 求解器路径 % 1. 基础数据 [gen, load_pred, wind_pred] init_data(); % 2. 生成场景并削减 [load_scen, wind_scen, prob] gen_scenarios(load_pred, wind_pred, 500, 20); % 3. 建立优化模型 [x, u, z, objective] build_model(gen, load_scen, wind_scen, prob); % 4. 求解 ops sdpsettings(solver, gurobi, mipgap, 0.01, verbose, 0); result optimize(x, objective, ops); % 5. 结果整理与画图 plot_results(value(u), value(x), load_scen, wind_scen);主函数逻辑很直白。我习惯把参数设置全部放在init_data里这样后面做灵敏度分析时只需要改一个文件不用动主体代码。3.2 Yalmip建模的变量与约束写法Yalmip是Matlab环境下最顺手的建模语言。关键点是决策变量的维度安排。假设有3台火电机组、24个时段、20个场景出力变量就是三维数组P sdpvar(3, 24, 20, full)。启停变量是和场景无关的第一阶段变量所以是二维0-1变量u binvar(3, 24)。功率平衡约束写起来要小心索引顺序Constraints []; for t 1:24 for k 1:20 Constraints [Constraints, ... sum(P(:,t,k), 1) wind_use(t,k) load_scen(t,k)]; end end这里wind_use是决策变量表示每个场景下实际使用的风电功率它的上限是场景削减后的风电场景值。所以还要写Constraints [Constraints, 0 wind_use(t,k) wind_scen(t,k)];机组出力上下限约束要跟启停变量耦合写成Constraints [Constraints, ... repmat(u(:,t),1,20) .* Pmin P(:,t,:) repmat(u(:,t),1,20) .* Pmax];3.3 目标函数中碳成本与弃风惩罚的写法目标函数我用sum叠加。燃煤成本是二次函数Yalmip可以直接写二次项Cplex和Gurobi都支持MIQP求解。fuel_cost 0; for i 1:3 for t 1:24 for k 1:20 fuel_cost fuel_cost prob(k) * (a(i)*P(i,t,k)^2 b(i)*P(i,t,k) c(i)); end end end碳成本需要先算每个时段的实际碳排放。碳排放强度可以简化为每台机组的固定排放系数t/MWh乘以出力得到总排放量再乘以碳价并减去配额。E_emit sum(sum(sum(prob .* carbon_intensity .* P))); % 期望碳排放 carbon_cost carbon_price * (E_emit - E_quota);弃风惩罚先算期望弃风量curtail sum(sum(sum(prob .* (wind_scen - wind_use)))); curtail_cost 80 * curtail; % 弃风惩罚单价80元/MWh目标函数最后写成objective fuel_cost startup_cost carbon_cost curtail_cost;。我提醒一句碳成本、燃煤成本和弃风惩罚的量纲必须统一否则碳价取太大时模型会“研究性”地疯狂弃风看起来很科学实际上是目标函数失衡。3.4 求解器配置与调试建议求解器首选Gurobi或Cplex学术版都是免费授权的。装了之后要确保Yalmip找得到常用做法是addpath(genpath(C:\gurobi1100\matlab)); savepath;用yalmiptest快速检测求解器是否加载成功。如果没有商业求解器Matlab自带的intlinprog也能跑但场景数一旦超过10个、时段到24小时速度会非常酸爽。我建议至少装一个GurobiMILP的求解速度能差几十倍。调试阶段有个技巧先把场景数设成1误差标准差设成0这样等价于确定性模型跑通之后再加随机场景。我踩过很多次坑大部分模型错误其实是索引写错确定性版本能帮你快速定位问题在建模逻辑还是随机因素。4. 实操过程与运行结果分析4.1 算例系统与基础参数我用一个简化系统跑通整个流程3台火电机组、1个风电场、1个聚合负荷。机组参数如下表碳价基准值取80元/吨风电装机600MW24小时预测负荷峰值1200MW。机组出力下限(MW)出力上限(MW)煤耗系数a(元/MW²h)煤耗系数b(元/MWh)碳排放强度(t/MWh)爬坡速率(MW/h)G11003500.00022200.5580G2803000.00032000.7270G3501500.00041800.9145注意机组容量越大碳排放强度越低这是我们给“大机组”的低碳偏好。煤耗系数a和b的量纲一定要统一我最初写代码时就把a的单位搞错过导致成本结果偏了几个量级。4.2 场景生成与求解规模设置采样500个原始场景削减到20个等概率场景。24个时段每个场景下3台机组出力加上20个场景下的风电使用量决策变量数量在几百个0-1变量是3×2472个。这种规模Gurobi通常几十秒到几分钟就能收敛。如果用的是老版本Matlab和Gurobi我建议把mipgap设为0.5%或1%不要追求默认精度速度快很多。另外Yalmip求解时如果遇到“连续变量过多导致内存爆炸”优先检查是不是约束里出现了norm函数它会引入很多辅助变量能不用就不用。4.3 结果可视化与对比维度跑完以后必画的图有三个。第一张是24小时机组出力堆叠图把G1、G2、G3和风电出力叠起来与负荷曲线对照能直观看到风电多发时段火电怎么压出力。第二张是各场景下的旋转备用容量曲线看是不是所有场景下备用都满足要求。第三张是碳排放柱状图把“无碳约束”“有碳成本”“碳价提高”三种情况对比这是论文里最常用的结果图。对比维度上我建议做三组确定性调度 vs 随机调度看风电利用率和系统总成本差异不同场景数比如5、20、50下的目标函数值看场景削减带来的近似误差不同碳价80、160、320元/吨下的碳排放总量验证碳价对低碳调度的引导作用。4.4 参数灵敏度分析与扩展方向参数灵敏度分析是课题和论文的加分项。最常见的是扫描碳价你会发现碳价从0开始逐步上升时碳排放总量并不是线性下降而是存在一个阶梯效应——碳价高到一定值后最优方案会切换机组组合方式。这个拐点很值得写在论文里。还有一个可以做的扩展是把网络约束加进来。系统级调度不关心线路潮流但如果题目要求考虑网络安全就得引入直流潮流方程用节点注入功率求解支路潮流再加线路传输容量约束。这个扩展在Yalmip里实现不算太难模型规模会大很多但结论也更可信。储能和碳捕集也是后续很自然的扩展方向模型框架基本不用改只需要在功率平衡里加几个变量和约束。5. 常见问题与排查技巧实录5.1 模型不可行的必查项模型不可行是最容易遇到的问题。我通常按这个顺序排查功率平衡约束的左右格式是否一致是不是把预测值用成了场景值旋转备用约束在风电出力极端高和极端低时是否都能满足启停变量和出力上下限的耦合条件是否写反比如Pmin ≤ P还是Pmin ≤ u*P负荷场景和风电场景是否出现负数截断处理是否生效。还有一个百试百灵的技巧在约束组里临时加一个虚拟松弛变量比如把功率平衡写成sum(P)wind_use load_scen slack然后看求解出来的slack分布在哪些时段和场景。slack不为0的位置就是问题所在。5.2 求解时间爆炸怎么办场景数太多、0-1变量交互复杂都会导致求解慢。减场景数是最直接的方案从20个减到10个可能只要一半时间结果差别不大。其次可以调大mipgap到1%2%在做灵敏度分析时这个误差完全可接受。如果还是慢检查目标函数里有没有引入大量非凸二次项。燃煤成本二次函数如果写成a*P^2对求解器是友好的但如果你手滑写成a*P*P会把矩阵维度搞错求解器可能被迫用非线性算法速度慢得离谱。我用Yalmip遇到最夸张的一次是连续变量只有几百个但求解器跑了三个小时没结果最后发现就是某一行把sum写丢导致目标函数里多了一整张矩阵乘法。5.3 结果“弃风为零”或“成本为负”异常弃风量正好是零不一定是错但要检查是不是因为弃风惩罚单价太低导致模型“懒得弃风”。比较合理的现象是在风电大发时段若机组已经压到最小技术出力弃风量会非常小但有正数。如果是确定性模型下弃风量严格为零且所有机组都有调峰空间那也说得通但随机模型下因为每个场景风电波动几乎不可能所有场景都零弃风。成本为负十有八九是配额的碳交易机制算反了。注意碳交易成本表达式如果实际排放小于配额目标函数是减项总成本可以变低但不应该出现“总成本远小于煤耗成本”的情况。仔细检查配额的单位如果配额设成总排放的上限而非差值表达式会完全乱掉。5.4 版本兼容与系统环境问题Matlab版本和Yalmip、Gurobi的兼容性是个常年坑。Matlab R2023a之后有些内建函数行为变化旧版Yalmip可能出现“找不到sdpvar构造函数”的问题。我的建议是统一用新版本Yalmip并安装对应求解器的最新Matlab接口。Gurobi的License服务器偶尔会被杀毒软件误杀启动报License Manager Error -8时先关掉安全软件白名单再重启License服务。另外提醒一下用笔记本电脑跑大算例时一定要把系统休眠关掉。之前我跑一个2400变量的大模型中途电脑休眠导致求解状态丢失白白浪费了三个小时。现在我的习惯是写代码前先执行system(powercfg -change -standby-timeout-ac 0)或者直接在电源设置里选了“从不睡眠”。这套代码和模型框架我已经在多个算例上重复用过。做课题最怕的不是模型复杂而是问题定位困难。把场景生成、建模、求解三个环节拆开各自测试确认独立模块没问题再组装几乎能避开九成以上的调试痛苦。