1. 从确定性到不确定性为什么源荷两侧都不能再用点预测做电力系统调度的朋友应该都有过这种憋屈时刻用点预测值算好了一版调度计划结果实际运行中风速突然起来风电出力比预测高出一截或者晚高峰负荷比预测值低了上百兆瓦最后只能靠AGC自动发电控制和备用硬扛。单靠点预测做日前调度预测误差越大结果偏离实际就越多这是所有调度模型都绕不开的痛点。风电出力预测误差的来源太多了。数值天气预报本身有误差地形影响没有完全建模风电场的尾流效应、风机切入切出风速附近的不连续出力都会让预测值和实际值偏离。短时尺度预测误差通常可以达到额定出力的10%到20%在极端天气过程前后甚至可能更高。如果系统里风电渗透率再高一点这么大体量的不确定性已经不是实时调整一下能解决的了。负荷侧同样不省心。电力负荷受天气、经济活动和用户行为影响虽然大趋势可预测但存在明显的日内波动和随机成分。特别是近年来分布式光伏、电动汽车充电桩等新元素进入负荷侧使得纯负荷曲线本身也变得更加难以琢磨。负荷短期预测误差在1%到3%之间听起来不大但放到区域电网就可能是几十上百兆瓦的量级跟风电误差叠加在一起调度裕度被大幅压缩。更麻烦的是风电出力和负荷误差有时是同向叠加的——高估风电的同时低估负荷或者反过来。这使得净负荷负荷减去风电出力的不确定性往往比两个变量各自的不确定性简单相加还要复杂。考虑源荷两侧不确定性的意义就在于调度计划需要有足够的灵活容量来同时应对两类偏差而不是只盯住某一边。只用确定性模型做日前调度通常只出一条基准曲线真到了运行时刻风电不出力或者负荷偏大系统就得靠旋转备用和快速启停机组兜底。备用不足就面临切负荷风险备用过多则经济性下降。这个问题在低碳调度里更加棘手因为碳排放约束和碳交易成本会改变机组出力的优先次序常规备用机组的碳排放特性也会影响整体碳指标。1.1 风电预测误差不只是小扰动风电出力预测误差的工程意义很大程度上取决于风电渗透率。渗透率低时一点预测误差靠系统惯性和快速调节机组就能消化渗透率高了风电预测误差可能直接触发备用容量不足甚至影响频率安全。用一个直观的数字说明一个100MW风电场预测误差按15%算就是15MW。如果系统里常规机组的调节速率平均是每分钟3%额定出力一组额定容量500MW的机组也要好几分钟才能完全弥补这个缺口。如果是日内滚动调度这个恢复时间窗口内净负荷一旦再快速变化风险就很明显。更麻烦的是风电出力的概率分布往往是非对称、多峰的。风速在切入风速附近时出力对风速变化非常敏感风速稍微波动出力就大幅变化风速超过额定风速后出力又趋于饱和。用一个对称分布比如简单的高斯误差去描述风电出力预测误差在某些时段会有系统性偏差。常见的处理思路有两个一是用风速的威布尔分布结合风电功率曲线生成出力场景二是在风电功率预测值基础上叠加一个符合Beta分布的误差项。Beta分布定义在[0,1]区间恰好对应风电出力的归一化范围灵活性很好实际拟合效果也不错。1.2 负荷预测偏差同样不可忽略负荷预测在城市和区域电网的短期误差通常在1%到3%绝对量不容小觑。负荷曲线往往有明显的早晚双峰结构高峰时段的预测偏差对备用和调度的影响最大因为高峰时段的调节空间本来就小。负荷不确定性的描述相对规范一般直接用正态分布。负荷预测中心值用短期负荷预测的结果标准差取预测点历史误差的标准差也可以分时段设置不同的标准差因为高峰时段的负荷误差绝对量明显更大。源荷两侧同时考虑时有一个比较关键的问题相关性。风电出力和负荷在全时段上相关性通常不强但在同一地理区域的极端天气事件中比如寒潮导致风电静风和负荷攀升同时出现两者会出现明显的尾部相关。严格的研究应该用联合分布甚至Copula函数建模但工程上为了可操作性普遍假设两者独立各自生成场景再做组合约束。我的代码也沿用这个假设在后面的灵敏度分析里可以观察这个假设对结果的影响。2. 低碳调度的核心机制碳交易成本怎么进目标函数低碳调度和普通经济调度的本质区别不在于约束方程多复杂而在于目标函数里多了一项碳相关的成本并且因此改变了机组出力的优先级排行。先理解碳交易机制的基本逻辑。在碳排放权交易体系下发电企业会被分配一定数量的碳排放配额。配额可以覆盖全部或部分排放量实际排放超出配额的部分需要从市场上购买碳排放权实际排放低于配额的部分则可以出售获利相当于给低碳机组一笔额外收入。2.1 碳配额怎么算目前常用的配额分配方法有历史强度法按历史排放强度分配和基准线法按行业先进水平分配。在电力系统调度模型中简化处理时一般用基准线法配额 装机容量 × 年等效利用小时数 × 基准排放强度但调度模型是日前或小时级问题所以需要把配额按时间尺度折算。常见做法是给每个机组设定一个单位电量碳排放配额比如0.6到0.8 tCO2/MWh实际单位排放强度高于这个值就得购买碳配额低于这个值就能出售。这个值在不同国家和地区的碳市场中差异很大具体研究时需要按目标市场的政策参数来设。2.2 碳交易成本的表达式假设机组i在t时刻出力为P(i,t)碳排放强度为E(i)tCO2/MWh碳配额强度为R(i)碳价为λ那么机组i的碳排放量为C(i,t) E(i) × P(i,t)配额为Q(i,t) R(i) × P(i,t)碳交易成本为T(i,t) λ × [C(i,t) - Q(i,t)] λ × [E(i) - R(i)] × P(i,t)从式子里能看出一个非常重要的现象如果某台机组的碳排放强度高于配额强度碳交易成本是正的出力越多需要买的碳排放权越多反过来如果排放强度低于配额强度碳交易成本为负值相当于给低碳机组一种绿色补贴出力越多越有利。这就改变了火电机组之间的经济竞争关系——一台低效但燃料便宜的老机组在碳价提高之后可能反而不如高效但燃料稍贵的机组有调度优先级。2.3 目标函数完整形式低碳调度的目标函数可以写为min F Σ(t) Σ(i) [ F_fuel(i,t) ST(i,t) ] Σ(t) Σ(i) [ λ × (E(i) - R(i)) × P(i,t) ]其中F_fuel(i,t)是燃料成本通常表达为二次函数 F a bP cP²ST(i,t)是启停成本由0-1变量控制。实际代码里燃料成本的常数项a可以省略它对最优解的出力分布没有影响只影响目标函数绝对值对场景对比和灵敏度分析也没影响。约束条件与常规经济调度类似包括系统功率平衡约束所有机组出力加风电出力等于负荷机组出力上下限约束机组爬坡速率约束最小启停时间约束旋转备用约束碳排放总量约束可选政策面设置系统碳预算时加入碳总量约束如果加上相当于给整个系统设定了一个碳排放红线模型会变得更严格也更容易看出碳价政策替代硬约束的力度差异。我在代码里把这条做成可选项注释里说明了开关位置。3. 不确定性建模为什么我选了场景法以及场景怎么来3.1 三种不确定性处理方法的工程化对比不确定性建模的常见手段有场景法随机规划、鲁棒优化和机会约束规划三种各有适用场景方法数学模型输出计划优点缺点场景法随机规划多场景期望最小化确定性出力计划概念直观经济指标明确场景数量影响计算量鲁棒优化不确定集合内最坏情况保守出力计划安全性强集合参数易调结果偏保守经济性差机会约束规划概率约束允许小概率越限出力计划置信水平风险与经济平衡求解复杂分布假设敏感对于需要复现和做教学演示的Matlab代码来说场景法是最合适的选择。它把随机性问题转化成多个确定性问题的加权组合YALMIP配合CPLEX或者Gurobi可以直接求解不需要自己写Benders分解或列约束生成这类复杂算法后续做灵敏度分析和结果可视化也很方便。3.2 源荷场景生成的具体做法我用的方法分三步。第一步生成风电出力场景。设风电预测出力为w_f(t)预测误差用Beta分布描述每个场景的风电出力为w_s(t) w_f(t) Δw(t)Δw按Beta分布采样Beta分布的形状参数根据历史误差统计拟合。注意Beta分布的支撑区间是[0,1]所以需要结合风电装机容量做量纲换算和上下界修正采样后要强制截断避免场景里出现负出力或超过装机容量的不合理值。第二步生成负荷场景。负荷预测中心值为L_f(t)误差服从正态分布N(0, σ²_L(t))σ_L(t)取对应时段历史预测误差的标准差。峰谷时段的σ_L可以做差异化处理高峰时段放大、低谷时段缩小更贴合实际情况。第三步场景削减。我首先生成1000到2000个初始场景然后用K-means聚类把场景数降到10到20个。K-means聚类质心要按概率加权每个簇里包含的原场景数量作为权重这样削减后的期望成本和原始场景集期望成本偏差会小很多。有些同行用同步回代消除法backward reduction实测效果也不错但K-means实现起来更简洁Matlab内置了kmeans函数不需要额外写概率距离计算的循环代码量少很多跑起来也快。3.3 场景削减后的一个隐蔽坑K-means聚类完成后别忘记对削减后的场景做平滑化处理。因为聚类质心可能让风电出力在相邻时段之间出现跳变这在功率平衡约束里没有问题但在爬坡约束里可能导致本来不该触发的爬坡被触发了这种失真。我处理的办法是对聚类后的场景序列做时间维度的中值滤波或者约束每个时段出力变化不超过预测场景的最大变化率。这个小细节直接关系到最后结果里爬坡约束是否真实有效。4. Matlab代码实现从数学模型到可运行代码4.1 代码结构与数据准备我的代码文件组织如下main.m 主程序串联场景生成、建模、求解、出图 case30.m 系统数据机组参数、风电场参数、负荷曲线 scenario_gen.m 源荷场景生成与削减 build_model.m YALMIP建模与求解 plot_result.m 结果绘图与分析指标计算case30.m里存放的主要数据包括机组数、爬坡率、出力上下限、燃料成本系数、碳排放强度、碳配额强度、风电装机容量、负荷时序数据。数据来源用的是IEEE 30节点系统的标准算例风电场接在某个节点上替换掉一部分常规机组容量。设备参数表大致如下数值可按研究需要调整机组出力上限(MW)爬坡率(MW/h)燃料成本系数b碳排放强度(tCO2/MWh)碳配额强度G12005020.50.920.70G21504018.20.850.70G31003024.11.050.70风电80-000风电不产生碳排放燃料成本为零但引入不确定性后需要额外配置备用容量这个隐形代价会体现在约束和成本计算中。4.2 场景生成核心代码场景生成部分的核心逻辑大致如下function [scenLoad, scenWind, prob] scenario_gen(load_f, wind_f, windCap, numClu) % load_f, wind_f: 预测序列长度为T % windCap: 风电场装机容量 % numClu: 聚类后的场景数 T length(load_f); nInit 1500; % 初始场景数 alpha_w 2.0; beta_w 2.5; % Beta分布形状参数按历史误差拟合 sigma_l 0.02 * load_f; % 负荷误差标准差按比例设定 % 生成初始场景 initLoad zeros(nInit, T); initWind zeros(nInit, T); for k 1:nInit delta_w windCap * (betarnd(alpha_w, beta_w, 1, T) - 0.5) * 1.8; initWind(k,:) wind_f delta_w; initWind(k, initWind(k,:) 0) 0; initWind(k, initWind(k,:) windCap) windCap; initLoad(k,:) load_f normrnd(0, sigma_l); initLoad(k, initLoad(k,:) 0) 0; end % K-means聚类削减 X [initLoad initWind]; [idx, C] kmeans(X, numClu, Replicates, 10); prob accumarray(idx, 1) / nInit; % 各场景概率 scenLoad C(:, 1:T); scenWind C(:, T1:2*T); % 时间维平滑避免场景内相邻时段跳变过大 for s 1:numClu scenWind(s,:) medfilt1(scenWind(s,:), 3); scenLoad(s,:) medfilt1(scenLoad(s,:), 3); end end需要注意Beta分布的形状参数直接决定了误差分布是偏左还是偏右。如果场址的历史数据显示高估概率更大就把α调小β调大反之则调大α调小β。这个参数不是随便拍脑袋的最好用实际风电场的出力历史数据做一次最大似然估计。4.3 YALMIP建模变量定义与约束实现YALMIP建模的关键在于把场景索引编进变量里。机组出力P是三维变量机组数×时段数×场景数而0-1启停变量u是二维的因为机组的开停机决策在所有场景下保持一致——这是场景法的标准处理方式叫非预期性约束。核心代码示意%% 决策变量定义 P sdpvar(nG, T, numScen, full); % 机组出力 u binvar(nG, T, full); % 开停机状态场景间一致 startup binvar(nG, T, full); % 启动标志 shutdown binvar(nG, T, full); % 停机标志 %% 约束集合 Constraints []; % 功率平衡每个场景下全部机组出力风电负荷 for s 1:numScen for t 1:T Constraints [Constraints, sum(P(:,t,s), 1) scenWind(s,t) scenLoad(s,t)]; end end % 出力上下限与启停耦合 for s 1:numScen for t 1:T Constraints [Constraints, P(:,t,s) u(:,t) .* Pmax]; Constraints [Constraints, P(:,t,s) u(:,t) .* Pmin]; end end % 爬坡约束 for s 1:numScen for t 2:T Constraints [Constraints, P(:,t,s) - P(:,t-1,s) ramp_up]; Constraints [Constraints, P(:,t,s) - P(:,t-1,s) -ramp_down]; end end % 旋转备用约束常规机组的可调容量应对风电和负荷场景偏差 for t 1:T Constraints [Constraints, sum(P(:,t,1), 1) sum(scenWind(:,t)) ... scenLoad(:,t) * (1 reserve_rate)]; end %% 目标函数场景期望成本 Objective 0; for s 1:numScen for t 1:T % 燃料成本 Objective Objective prob(s) * sum(fuel_b .* P(:,t,s) fuel_c .* P(:,t,s).^2); % 碳交易成本 Objective Objective prob(s) * sum(carbon_price .* (E_int - R_int) .* P(:,t,s)); % 启停成本 Objective Objective prob(s) * sum(startup_cost .* startup(:,t)); end end %% 求解 ops sdpsettings(solver, cplex, verbose, 2, debug, 0); result optimize(Constraints, Objective, ops);这段代码有几个细节值得展开。一是启停变量u所有场景共用这保证了日前调度给出的机组开停计划是确定性的不会随场景变化而变卦。如果让u也按场景独立实际上就变成看人下菜碟在实际运行中做不到。二是爬坡约束需要特别注意时段边界。t1时段没有历史出力我给的约束从t2开始避免边界效应。部分文献会在初始时段假设机组以初始状态连续出力那是另一种处理方式。三是目标函数中的二次项fuel_c * P²在CPLEX里求解混合整数二次规划MIQP是可以的但如果数据里的fuel_c很小数值上可能出现求解器报告数值问题的现象最好是先做归一化处理或者把成本归一到百万美元级别减少量级差距。四是场景概率的权重必须加对。每个场景的目标函数贡献要乘以prob(s)否则所有场景等权叠加聚类削减就白做了。4.4 求解器配置与常见坑YALMIP后端求解器我推荐CPLEX或Gurobi学术界用得最多。如果机器上没装也可以用MATLAB自带的intlinprog求解但大规模场景下速度会明显变慢。实测下来20个场景、24时段、6台机组的混合整数二次规划问题CPLEX大概几秒到十几秒能出结果intlinprog可能要几分钟。版本兼容性是一个容易踩的坑。Matlab版本和CPLEX版本必须匹配YALMIP官方文档有兼容性表格安装时建议先查一下。我自己的环境是Matlab 2023b配CPLEX 12.10跑常规70维左右的变量没有问题。如果遇到无法识别cplex求解器的报错多半是路径没添加正确用addpath(genpath(C:\...\cplex\matlab))重新加载一遍即可。再提一个可能有用的技巧如果初始求解时发现某些二进制变量无法收敛先检查约束里有没有等式两端同时乘以二进制变量这种非线性写法YALMIP对这类问题会自动引入大M法如果大M值设得过大或过小都会导致松弛问题病态表现为求解器输出infeasible。5. 算例验证改装后的IEEE 30节点系统能说明什么5.1 基础数据与参数设置我用的算例是在IEEE 30节点系统基础上改装的。常规机组取3台火电总装机450MW风电场装机80MW负荷峰值230MW负荷曲线取典型冬季峰谷形状。碳价格λ设置为25元/tCO2碳配额强度统一取0.70 tCO2/MWh旋转备用率取5%。场景数量取15个初始场景1500个K-means削减后各场景概率分布相对均匀。对比方案是确定性调度不考虑任何不确定性用预测值直接求解和考虑源荷两侧不确定性的场景法调度。5.2 确定性调度和不确定性调度的对比先说总成本。确定性调度的期望成本比场景法调度低约4%到6%这个差异看起来不大但注意两点一是确定性方案没有预留足够的爬坡和备用容量实际运行中遭遇预测偏差时的二次调频成本和外购高价电力没有计入门内二是如果计算确定性方案在15个真实场景下的模拟运行成本它会比它声称的成本高出一截而场景法方案在真实场景下的模拟成本与其优化目标基本一致。换句话说确定性方案省下来的是纸面成本不是真实成本。再看碳排放。场景法调度下的碳排放总量比确定性调度高约1%到2%原因在于为了应对不确定性效率稍低但爬坡能力更强的机组会有更多的出力和启停机会。这个时候如果碳价进一步上升系统会趋向于让高效机组多带基荷、低效机组只做调节碳排放总量反而可能下降——这就是碳价对调度结构的影响在结果上的最终体现。5.3 碳价与风电比例的敏感性分析把碳价从20元/tCO2逐步提高到60元/tCO2观察系统总碳排放和总成本的变化可以发现一个有趣的非线性现象碳价在20到40元区间时碳排放下降较慢因为此时低效机组依然有一部分的基荷份额当碳价超过40元后碳排放下降速度明显加快说明临界碳价促使低效机组的出力占比大幅压缩。这个临界碳价的位置取决于高效和低效机组之间的燃料成本差和排放强度差。再调整风电渗透率从10%、20%到30%场景法调度的优势会越来越明显。渗透率30%时确定性方案面对的不确定性已大到不加处理就很难保证安全约束的满足程度而场景法方案虽然总成本上升但所有场景下均满足功率平衡和爬坡约束。这正好印证了开篇的观点高比例新能源接入之后不确定性不再是小修小补的问题而是调度建模的前提条件。6. 个人实操中的几个体会代码跑通只是第一步真正花时间的往往在数据清洗和结果调校上。我把自己反复踩过的一些教训列在这里供后来者参考。第一场景数量和求解时间的平衡点。一开始我以为场景越多越好结果30个场景的模型求解时间比15个场景翻了三倍不止而优化结果的期望成本差异很小。实测下来15到20个场景对IEEE 30节点这样的算例已经足够再往上的边际收益微乎其微。如果系统规模更大可以先用少量场景快速验证模型逻辑等逻辑无误后再加大场景数跑最终版。第二约束的可行域要留一点呼吸空间。有些文献里功率平衡约束写成严格等式实际求解时一旦场景生成存在数值误差很容易出现infeasible。我的做法是在功率平衡约束中引入一个很小的松弛量比如0.01MW既不影响结果的合理性又避免了很多无谓的不可行问题。这个松弛在论文里可以注明在实际工程里本来就存在这样量级的计量误差。第三K-means聚类的初始化对结果有不可忽略的影响。我设置了Replicates, 10让Matlab重复跑10次选最优的聚类结果这样能有效避免陷入局部最优。如果时间紧张至少要设3次否则个别场景会明显偏离原始分布特征。第四结果可视化不要只画机组出力曲线要把风电场景带、负荷场景带和系统备用容量放在同一张图上看。场景带直观展示不确定性的量级调度结果落在场景带内才说明调度方案对不确定性免疫这一步对判断模型的可靠性非常重要。最后如果后续想把模型扩展得更精细可以考虑加入储能系统、需求响应资源或者碳捕集装置。存储的充放电策略天然适合用不确定性场景框架来优化——正是这套场景法框架的扩展场景储能配置问题的建模逻辑和本文完全一致只是决策变量再多一个储能SOC的状态转移约束而已。 计及源荷不确定性的低碳调度场景法建模与MATLAB实现
1. 从确定性到不确定性为什么源荷两侧都不能再用点预测做电力系统调度的朋友应该都有过这种憋屈时刻用点预测值算好了一版调度计划结果实际运行中风速突然起来风电出力比预测高出一截或者晚高峰负荷比预测值低了上百兆瓦最后只能靠AGC自动发电控制和备用硬扛。单靠点预测做日前调度预测误差越大结果偏离实际就越多这是所有调度模型都绕不开的痛点。风电出力预测误差的来源太多了。数值天气预报本身有误差地形影响没有完全建模风电场的尾流效应、风机切入切出风速附近的不连续出力都会让预测值和实际值偏离。短时尺度预测误差通常可以达到额定出力的10%到20%在极端天气过程前后甚至可能更高。如果系统里风电渗透率再高一点这么大体量的不确定性已经不是实时调整一下能解决的了。负荷侧同样不省心。电力负荷受天气、经济活动和用户行为影响虽然大趋势可预测但存在明显的日内波动和随机成分。特别是近年来分布式光伏、电动汽车充电桩等新元素进入负荷侧使得纯负荷曲线本身也变得更加难以琢磨。负荷短期预测误差在1%到3%之间听起来不大但放到区域电网就可能是几十上百兆瓦的量级跟风电误差叠加在一起调度裕度被大幅压缩。更麻烦的是风电出力和负荷误差有时是同向叠加的——高估风电的同时低估负荷或者反过来。这使得净负荷负荷减去风电出力的不确定性往往比两个变量各自的不确定性简单相加还要复杂。考虑源荷两侧不确定性的意义就在于调度计划需要有足够的灵活容量来同时应对两类偏差而不是只盯住某一边。只用确定性模型做日前调度通常只出一条基准曲线真到了运行时刻风电不出力或者负荷偏大系统就得靠旋转备用和快速启停机组兜底。备用不足就面临切负荷风险备用过多则经济性下降。这个问题在低碳调度里更加棘手因为碳排放约束和碳交易成本会改变机组出力的优先次序常规备用机组的碳排放特性也会影响整体碳指标。1.1 风电预测误差不只是小扰动风电出力预测误差的工程意义很大程度上取决于风电渗透率。渗透率低时一点预测误差靠系统惯性和快速调节机组就能消化渗透率高了风电预测误差可能直接触发备用容量不足甚至影响频率安全。用一个直观的数字说明一个100MW风电场预测误差按15%算就是15MW。如果系统里常规机组的调节速率平均是每分钟3%额定出力一组额定容量500MW的机组也要好几分钟才能完全弥补这个缺口。如果是日内滚动调度这个恢复时间窗口内净负荷一旦再快速变化风险就很明显。更麻烦的是风电出力的概率分布往往是非对称、多峰的。风速在切入风速附近时出力对风速变化非常敏感风速稍微波动出力就大幅变化风速超过额定风速后出力又趋于饱和。用一个对称分布比如简单的高斯误差去描述风电出力预测误差在某些时段会有系统性偏差。常见的处理思路有两个一是用风速的威布尔分布结合风电功率曲线生成出力场景二是在风电功率预测值基础上叠加一个符合Beta分布的误差项。Beta分布定义在[0,1]区间恰好对应风电出力的归一化范围灵活性很好实际拟合效果也不错。1.2 负荷预测偏差同样不可忽略负荷预测在城市和区域电网的短期误差通常在1%到3%绝对量不容小觑。负荷曲线往往有明显的早晚双峰结构高峰时段的预测偏差对备用和调度的影响最大因为高峰时段的调节空间本来就小。负荷不确定性的描述相对规范一般直接用正态分布。负荷预测中心值用短期负荷预测的结果标准差取预测点历史误差的标准差也可以分时段设置不同的标准差因为高峰时段的负荷误差绝对量明显更大。源荷两侧同时考虑时有一个比较关键的问题相关性。风电出力和负荷在全时段上相关性通常不强但在同一地理区域的极端天气事件中比如寒潮导致风电静风和负荷攀升同时出现两者会出现明显的尾部相关。严格的研究应该用联合分布甚至Copula函数建模但工程上为了可操作性普遍假设两者独立各自生成场景再做组合约束。我的代码也沿用这个假设在后面的灵敏度分析里可以观察这个假设对结果的影响。2. 低碳调度的核心机制碳交易成本怎么进目标函数低碳调度和普通经济调度的本质区别不在于约束方程多复杂而在于目标函数里多了一项碳相关的成本并且因此改变了机组出力的优先级排行。先理解碳交易机制的基本逻辑。在碳排放权交易体系下发电企业会被分配一定数量的碳排放配额。配额可以覆盖全部或部分排放量实际排放超出配额的部分需要从市场上购买碳排放权实际排放低于配额的部分则可以出售获利相当于给低碳机组一笔额外收入。2.1 碳配额怎么算目前常用的配额分配方法有历史强度法按历史排放强度分配和基准线法按行业先进水平分配。在电力系统调度模型中简化处理时一般用基准线法配额 装机容量 × 年等效利用小时数 × 基准排放强度但调度模型是日前或小时级问题所以需要把配额按时间尺度折算。常见做法是给每个机组设定一个单位电量碳排放配额比如0.6到0.8 tCO2/MWh实际单位排放强度高于这个值就得购买碳配额低于这个值就能出售。这个值在不同国家和地区的碳市场中差异很大具体研究时需要按目标市场的政策参数来设。2.2 碳交易成本的表达式假设机组i在t时刻出力为P(i,t)碳排放强度为E(i)tCO2/MWh碳配额强度为R(i)碳价为λ那么机组i的碳排放量为C(i,t) E(i) × P(i,t)配额为Q(i,t) R(i) × P(i,t)碳交易成本为T(i,t) λ × [C(i,t) - Q(i,t)] λ × [E(i) - R(i)] × P(i,t)从式子里能看出一个非常重要的现象如果某台机组的碳排放强度高于配额强度碳交易成本是正的出力越多需要买的碳排放权越多反过来如果排放强度低于配额强度碳交易成本为负值相当于给低碳机组一种绿色补贴出力越多越有利。这就改变了火电机组之间的经济竞争关系——一台低效但燃料便宜的老机组在碳价提高之后可能反而不如高效但燃料稍贵的机组有调度优先级。2.3 目标函数完整形式低碳调度的目标函数可以写为min F Σ(t) Σ(i) [ F_fuel(i,t) ST(i,t) ] Σ(t) Σ(i) [ λ × (E(i) - R(i)) × P(i,t) ]其中F_fuel(i,t)是燃料成本通常表达为二次函数 F a bP cP²ST(i,t)是启停成本由0-1变量控制。实际代码里燃料成本的常数项a可以省略它对最优解的出力分布没有影响只影响目标函数绝对值对场景对比和灵敏度分析也没影响。约束条件与常规经济调度类似包括系统功率平衡约束所有机组出力加风电出力等于负荷机组出力上下限约束机组爬坡速率约束最小启停时间约束旋转备用约束碳排放总量约束可选政策面设置系统碳预算时加入碳总量约束如果加上相当于给整个系统设定了一个碳排放红线模型会变得更严格也更容易看出碳价政策替代硬约束的力度差异。我在代码里把这条做成可选项注释里说明了开关位置。3. 不确定性建模为什么我选了场景法以及场景怎么来3.1 三种不确定性处理方法的工程化对比不确定性建模的常见手段有场景法随机规划、鲁棒优化和机会约束规划三种各有适用场景方法数学模型输出计划优点缺点场景法随机规划多场景期望最小化确定性出力计划概念直观经济指标明确场景数量影响计算量鲁棒优化不确定集合内最坏情况保守出力计划安全性强集合参数易调结果偏保守经济性差机会约束规划概率约束允许小概率越限出力计划置信水平风险与经济平衡求解复杂分布假设敏感对于需要复现和做教学演示的Matlab代码来说场景法是最合适的选择。它把随机性问题转化成多个确定性问题的加权组合YALMIP配合CPLEX或者Gurobi可以直接求解不需要自己写Benders分解或列约束生成这类复杂算法后续做灵敏度分析和结果可视化也很方便。3.2 源荷场景生成的具体做法我用的方法分三步。第一步生成风电出力场景。设风电预测出力为w_f(t)预测误差用Beta分布描述每个场景的风电出力为w_s(t) w_f(t) Δw(t)Δw按Beta分布采样Beta分布的形状参数根据历史误差统计拟合。注意Beta分布的支撑区间是[0,1]所以需要结合风电装机容量做量纲换算和上下界修正采样后要强制截断避免场景里出现负出力或超过装机容量的不合理值。第二步生成负荷场景。负荷预测中心值为L_f(t)误差服从正态分布N(0, σ²_L(t))σ_L(t)取对应时段历史预测误差的标准差。峰谷时段的σ_L可以做差异化处理高峰时段放大、低谷时段缩小更贴合实际情况。第三步场景削减。我首先生成1000到2000个初始场景然后用K-means聚类把场景数降到10到20个。K-means聚类质心要按概率加权每个簇里包含的原场景数量作为权重这样削减后的期望成本和原始场景集期望成本偏差会小很多。有些同行用同步回代消除法backward reduction实测效果也不错但K-means实现起来更简洁Matlab内置了kmeans函数不需要额外写概率距离计算的循环代码量少很多跑起来也快。3.3 场景削减后的一个隐蔽坑K-means聚类完成后别忘记对削减后的场景做平滑化处理。因为聚类质心可能让风电出力在相邻时段之间出现跳变这在功率平衡约束里没有问题但在爬坡约束里可能导致本来不该触发的爬坡被触发了这种失真。我处理的办法是对聚类后的场景序列做时间维度的中值滤波或者约束每个时段出力变化不超过预测场景的最大变化率。这个小细节直接关系到最后结果里爬坡约束是否真实有效。4. Matlab代码实现从数学模型到可运行代码4.1 代码结构与数据准备我的代码文件组织如下main.m 主程序串联场景生成、建模、求解、出图 case30.m 系统数据机组参数、风电场参数、负荷曲线 scenario_gen.m 源荷场景生成与削减 build_model.m YALMIP建模与求解 plot_result.m 结果绘图与分析指标计算case30.m里存放的主要数据包括机组数、爬坡率、出力上下限、燃料成本系数、碳排放强度、碳配额强度、风电装机容量、负荷时序数据。数据来源用的是IEEE 30节点系统的标准算例风电场接在某个节点上替换掉一部分常规机组容量。设备参数表大致如下数值可按研究需要调整机组出力上限(MW)爬坡率(MW/h)燃料成本系数b碳排放强度(tCO2/MWh)碳配额强度G12005020.50.920.70G21504018.20.850.70G31003024.11.050.70风电80-000风电不产生碳排放燃料成本为零但引入不确定性后需要额外配置备用容量这个隐形代价会体现在约束和成本计算中。4.2 场景生成核心代码场景生成部分的核心逻辑大致如下function [scenLoad, scenWind, prob] scenario_gen(load_f, wind_f, windCap, numClu) % load_f, wind_f: 预测序列长度为T % windCap: 风电场装机容量 % numClu: 聚类后的场景数 T length(load_f); nInit 1500; % 初始场景数 alpha_w 2.0; beta_w 2.5; % Beta分布形状参数按历史误差拟合 sigma_l 0.02 * load_f; % 负荷误差标准差按比例设定 % 生成初始场景 initLoad zeros(nInit, T); initWind zeros(nInit, T); for k 1:nInit delta_w windCap * (betarnd(alpha_w, beta_w, 1, T) - 0.5) * 1.8; initWind(k,:) wind_f delta_w; initWind(k, initWind(k,:) 0) 0; initWind(k, initWind(k,:) windCap) windCap; initLoad(k,:) load_f normrnd(0, sigma_l); initLoad(k, initLoad(k,:) 0) 0; end % K-means聚类削减 X [initLoad initWind]; [idx, C] kmeans(X, numClu, Replicates, 10); prob accumarray(idx, 1) / nInit; % 各场景概率 scenLoad C(:, 1:T); scenWind C(:, T1:2*T); % 时间维平滑避免场景内相邻时段跳变过大 for s 1:numClu scenWind(s,:) medfilt1(scenWind(s,:), 3); scenLoad(s,:) medfilt1(scenLoad(s,:), 3); end end需要注意Beta分布的形状参数直接决定了误差分布是偏左还是偏右。如果场址的历史数据显示高估概率更大就把α调小β调大反之则调大α调小β。这个参数不是随便拍脑袋的最好用实际风电场的出力历史数据做一次最大似然估计。4.3 YALMIP建模变量定义与约束实现YALMIP建模的关键在于把场景索引编进变量里。机组出力P是三维变量机组数×时段数×场景数而0-1启停变量u是二维的因为机组的开停机决策在所有场景下保持一致——这是场景法的标准处理方式叫非预期性约束。核心代码示意%% 决策变量定义 P sdpvar(nG, T, numScen, full); % 机组出力 u binvar(nG, T, full); % 开停机状态场景间一致 startup binvar(nG, T, full); % 启动标志 shutdown binvar(nG, T, full); % 停机标志 %% 约束集合 Constraints []; % 功率平衡每个场景下全部机组出力风电负荷 for s 1:numScen for t 1:T Constraints [Constraints, sum(P(:,t,s), 1) scenWind(s,t) scenLoad(s,t)]; end end % 出力上下限与启停耦合 for s 1:numScen for t 1:T Constraints [Constraints, P(:,t,s) u(:,t) .* Pmax]; Constraints [Constraints, P(:,t,s) u(:,t) .* Pmin]; end end % 爬坡约束 for s 1:numScen for t 2:T Constraints [Constraints, P(:,t,s) - P(:,t-1,s) ramp_up]; Constraints [Constraints, P(:,t,s) - P(:,t-1,s) -ramp_down]; end end % 旋转备用约束常规机组的可调容量应对风电和负荷场景偏差 for t 1:T Constraints [Constraints, sum(P(:,t,1), 1) sum(scenWind(:,t)) ... scenLoad(:,t) * (1 reserve_rate)]; end %% 目标函数场景期望成本 Objective 0; for s 1:numScen for t 1:T % 燃料成本 Objective Objective prob(s) * sum(fuel_b .* P(:,t,s) fuel_c .* P(:,t,s).^2); % 碳交易成本 Objective Objective prob(s) * sum(carbon_price .* (E_int - R_int) .* P(:,t,s)); % 启停成本 Objective Objective prob(s) * sum(startup_cost .* startup(:,t)); end end %% 求解 ops sdpsettings(solver, cplex, verbose, 2, debug, 0); result optimize(Constraints, Objective, ops);这段代码有几个细节值得展开。一是启停变量u所有场景共用这保证了日前调度给出的机组开停计划是确定性的不会随场景变化而变卦。如果让u也按场景独立实际上就变成看人下菜碟在实际运行中做不到。二是爬坡约束需要特别注意时段边界。t1时段没有历史出力我给的约束从t2开始避免边界效应。部分文献会在初始时段假设机组以初始状态连续出力那是另一种处理方式。三是目标函数中的二次项fuel_c * P²在CPLEX里求解混合整数二次规划MIQP是可以的但如果数据里的fuel_c很小数值上可能出现求解器报告数值问题的现象最好是先做归一化处理或者把成本归一到百万美元级别减少量级差距。四是场景概率的权重必须加对。每个场景的目标函数贡献要乘以prob(s)否则所有场景等权叠加聚类削减就白做了。4.4 求解器配置与常见坑YALMIP后端求解器我推荐CPLEX或Gurobi学术界用得最多。如果机器上没装也可以用MATLAB自带的intlinprog求解但大规模场景下速度会明显变慢。实测下来20个场景、24时段、6台机组的混合整数二次规划问题CPLEX大概几秒到十几秒能出结果intlinprog可能要几分钟。版本兼容性是一个容易踩的坑。Matlab版本和CPLEX版本必须匹配YALMIP官方文档有兼容性表格安装时建议先查一下。我自己的环境是Matlab 2023b配CPLEX 12.10跑常规70维左右的变量没有问题。如果遇到无法识别cplex求解器的报错多半是路径没添加正确用addpath(genpath(C:\...\cplex\matlab))重新加载一遍即可。再提一个可能有用的技巧如果初始求解时发现某些二进制变量无法收敛先检查约束里有没有等式两端同时乘以二进制变量这种非线性写法YALMIP对这类问题会自动引入大M法如果大M值设得过大或过小都会导致松弛问题病态表现为求解器输出infeasible。5. 算例验证改装后的IEEE 30节点系统能说明什么5.1 基础数据与参数设置我用的算例是在IEEE 30节点系统基础上改装的。常规机组取3台火电总装机450MW风电场装机80MW负荷峰值230MW负荷曲线取典型冬季峰谷形状。碳价格λ设置为25元/tCO2碳配额强度统一取0.70 tCO2/MWh旋转备用率取5%。场景数量取15个初始场景1500个K-means削减后各场景概率分布相对均匀。对比方案是确定性调度不考虑任何不确定性用预测值直接求解和考虑源荷两侧不确定性的场景法调度。5.2 确定性调度和不确定性调度的对比先说总成本。确定性调度的期望成本比场景法调度低约4%到6%这个差异看起来不大但注意两点一是确定性方案没有预留足够的爬坡和备用容量实际运行中遭遇预测偏差时的二次调频成本和外购高价电力没有计入门内二是如果计算确定性方案在15个真实场景下的模拟运行成本它会比它声称的成本高出一截而场景法方案在真实场景下的模拟成本与其优化目标基本一致。换句话说确定性方案省下来的是纸面成本不是真实成本。再看碳排放。场景法调度下的碳排放总量比确定性调度高约1%到2%原因在于为了应对不确定性效率稍低但爬坡能力更强的机组会有更多的出力和启停机会。这个时候如果碳价进一步上升系统会趋向于让高效机组多带基荷、低效机组只做调节碳排放总量反而可能下降——这就是碳价对调度结构的影响在结果上的最终体现。5.3 碳价与风电比例的敏感性分析把碳价从20元/tCO2逐步提高到60元/tCO2观察系统总碳排放和总成本的变化可以发现一个有趣的非线性现象碳价在20到40元区间时碳排放下降较慢因为此时低效机组依然有一部分的基荷份额当碳价超过40元后碳排放下降速度明显加快说明临界碳价促使低效机组的出力占比大幅压缩。这个临界碳价的位置取决于高效和低效机组之间的燃料成本差和排放强度差。再调整风电渗透率从10%、20%到30%场景法调度的优势会越来越明显。渗透率30%时确定性方案面对的不确定性已大到不加处理就很难保证安全约束的满足程度而场景法方案虽然总成本上升但所有场景下均满足功率平衡和爬坡约束。这正好印证了开篇的观点高比例新能源接入之后不确定性不再是小修小补的问题而是调度建模的前提条件。6. 个人实操中的几个体会代码跑通只是第一步真正花时间的往往在数据清洗和结果调校上。我把自己反复踩过的一些教训列在这里供后来者参考。第一场景数量和求解时间的平衡点。一开始我以为场景越多越好结果30个场景的模型求解时间比15个场景翻了三倍不止而优化结果的期望成本差异很小。实测下来15到20个场景对IEEE 30节点这样的算例已经足够再往上的边际收益微乎其微。如果系统规模更大可以先用少量场景快速验证模型逻辑等逻辑无误后再加大场景数跑最终版。第二约束的可行域要留一点呼吸空间。有些文献里功率平衡约束写成严格等式实际求解时一旦场景生成存在数值误差很容易出现infeasible。我的做法是在功率平衡约束中引入一个很小的松弛量比如0.01MW既不影响结果的合理性又避免了很多无谓的不可行问题。这个松弛在论文里可以注明在实际工程里本来就存在这样量级的计量误差。第三K-means聚类的初始化对结果有不可忽略的影响。我设置了Replicates, 10让Matlab重复跑10次选最优的聚类结果这样能有效避免陷入局部最优。如果时间紧张至少要设3次否则个别场景会明显偏离原始分布特征。第四结果可视化不要只画机组出力曲线要把风电场景带、负荷场景带和系统备用容量放在同一张图上看。场景带直观展示不确定性的量级调度结果落在场景带内才说明调度方案对不确定性免疫这一步对判断模型的可靠性非常重要。最后如果后续想把模型扩展得更精细可以考虑加入储能系统、需求响应资源或者碳捕集装置。存储的充放电策略天然适合用不确定性场景框架来优化——正是这套场景法框架的扩展场景储能配置问题的建模逻辑和本文完全一致只是决策变量再多一个储能SOC的状态转移约束而已。