1. 为什么源荷两侧不确定性必须放在一个模型里如果你正在做电力系统优化调度尤其是含风电的低碳调度那么“源荷两侧不确定性”这几个字一定会出现在开题报告或者项目需求里。这个题目看起来不大但真正动手用Matlab实现一遍后你会发现90%的工作量在建模和场景处理上剩下那10%才是写求解器约束。这篇分享我就把整个实现过程拆开讲包括模型怎么建、场景怎么生成和削减、YALMIP里怎么写约束以及实际调参时踩过的一些坑。先说这个题目的核心矛盾。传统调度里机组出力可以按确定性的负荷预测来安排顶多留一个固定旋转备用率。但一旦加入大规模风电源头出力变得不确定负荷预测也不可能完全准两个不确定量叠在一起原来的“确定”调度就变得非常脆弱。可能风功率预测值明明够用实际风小了系统就要拉闸限电也可能负荷预测偏低实际用电猛涨机组爬坡跟不上频率直接跌落。所以我一直建议做含风电调度项目的人不要只在风电侧做文章负荷侧的不确定性必须一起进模型否则算出来的结果到实际运行里很容易“翻车”。1.1 风电出力的不确定性到底难在哪风电出力本质上是一个强随机过程取决于风速。研究里最常用Weibull分布去拟合风速的统计特性再用风机功率曲线把风速映射成出力。这个映射不是线性的切入风速、额定风速、切出风速三个断点会让场景分布产生明显的截断效应。再加上风电还有反调峰特性白天负荷高峰风速可能很低深夜负荷低谷反而风大。你要是只用确定性预测曲线做调度大概率会得到“低谷时段大量弃风、高峰时段备用不足”的结果实际并网运行根本不敢这么执行。还有一个大家容易忽略的问题不同时间断面之间的风电出力具有很强的时序相关性。比如凌晨西北风系统过境很可能连续十几个小时风电出力都偏高而不是每个小时独立波动。如果场景生成时忽略了这种时间相关性调度结果的备用配置会被严重低估。我早期用独立抽样生成风速场景每个时段单独抽分布结果模型算出来的成本很低一加上时序检验就崩了后来改成时序相关的抽样才正常。1.2 负荷侧的波动不只是“加个误差”负荷预测误差一般用正态分布描述均值为预测值标准差取预测值的2%到5%。但负荷侧不确定性不只是“数值上偏一点”这么简单。它同样存在日内的连续变化和突发突变比如气温骤降、大型工业用户临时启停、节假日负荷特性偏移等。更关键的是风电出力误差和负荷预测误差之间存在相关性如果风速预测偏低往往意味着天气系统判断错误温度变化也会带动负荷预测偏差。所以做源荷联合场景的时候不能简单地把两个独立随机变量拼在一起最好用协方差矩阵描述它们之间的相关性再通过Cholesky分解生成联合场景。这一步看起来只是数学处理但对调度结果的保守程度影响很大。用一个生活类比出门前看天气预报说“下午有雨”于是带伞但又看到商场大促可能很多人排队这是两件独立的事。如果你把“下雨”和“商场人少”误当成独立事件可能带伞却人挤人或者不带伞却被淋湿。源荷两侧不确定性也是同一个逻辑必须当成一个联合随机过程来考虑而不是各算各的。1.3 低碳调度不只是“少烧煤”低碳调度的本质是在传统经济调度的基础上把碳排放变成一个可量化的成本项或硬约束。现在国内研究里最常用的是碳交易机制给每台火电机组分配一个免费的碳排放配额实际排放量超出配额要花钱买低于配额可以把差额卖掉。这相当于给“低碳”标了个价格调度目标从单纯省煤耗变成了“煤耗成本碳交易成本弃风惩罚”综合最小化。如果只把碳排放当成一个固定上限约束问题会简单很多但实际效果并不好。因为碳上限设得太紧系统可能无解或被迫大量上高成本机组设得太松又起不到减排作用。碳交易机制的好处是它给了系统一个弹性在负荷顶峰允许用高排放但低边际成本的机组顶上多买碳配额在负荷低谷让高排放机组少发省下来的配额卖出去赚钱。我在模型里把碳交易成本直接放进目标函数同时保留一个系统总碳排放的软约束场景化处理这样既体现低碳目标又不会因为刚性约束导致无解。2. 低碳调度模型怎么搭目标、约束与不确定性方法这个部分我说一下模型的主干。这里不贴完整推导先把思路捋清楚因为完整程序里最复杂的其实就是两件事一是目标函数里怎么把成本、碳交易、弃风切负荷惩罚揉在一起二是每个时段的功率平衡和备用约束怎么在不确定性场景下成立。2.1 目标函数怎么设计假设系统里有常规火电机组、风电场和负荷调度周期为24小时。目标函数我通常写成三个部分之和。第一是运行成本包含煤耗成本和启停成本。煤耗成本一般用二次函数表示但在混合整数规划里二次项会让求解变慢实际程序里我会用分段线性近似。第二是碳交易成本。采用基准线法时碳交易成本等于碳价乘以实际碳排放量减去免费配额。这里的实际碳排放量用机组出力乘排放强度累加得到配额可以按机组容量或历史排放水平设定。第三是弃风惩罚和失负荷惩罚。风电消纳不是无限度的场景化调度里某些极端场景可能需要弃风或者切除少量负荷否则功率平衡可能无解。给这些操作设一个较高的惩罚系数让优化器只在极端场景下才允许使用。目标函数的表达式不复杂但要正确处理“期望成本”。因为场景法里每个场景有不同风电出力、不同负荷优化目标应该是对所有场景概率加权后的总期望成本而不是拿某一个确定场景算。这样算出来的调度方案才是面对未来不确定性时的平均最优而不是某个“预测值”下的最优。2.2 约束条件清单约束条件是整个模型里最容易漏项也最容易导致无解的地方。我按类别列一下功率平衡约束每个时段所有机组出力加上风电消纳量要等于负荷需求。在场景法中这个约束需要针对每个场景单独成立。机组出力上下限常规机组的出力不能越过技术最小出力和最大出力这个不难但要注意启动状态变量对出力的限制。爬坡约束机组在相邻时段的出力变化不能超过爬坡速率风电波动越大这里越容易成为瓶颈。旋转备用约束为了保证可靠性系统需要预留一定上备用和下备用。上备用用来应对风电突然减少或负荷突增下备用用来应对风电突增或负荷突减。场景法中备用通常取场景偏差的一定比例或固定值。碳排放约束除了目标函数里的碳交易成本我会再加一个系统总碳排放量的上限约束。这个上限不要设成硬性死值可以设置成基准排放的百分比否则极端场景下容易无解。弃风和切负荷变量约束弃风量和切负荷量必须非负且不能超过该时段的可用风电或负荷值。这一步给求解器留了“应急出口”。每条约束在Matlab里都对应一组矩阵或者YALMIP的表达式真正写起来并不难。难的是场景数量和约束数量的增长关系每增加一个场景功率平衡和备用约束就扩展一份。如果不做场景削减500个场景直接进MIP内存和求解时间都会爆炸。2.3 不确定性建模场景法与鲁棒优化怎么选含不确定性的调度建模主流方法大致分三类。第一类是随机规划核心思路是生成大量场景用期望值做目标代表方法是两阶段随机规划。第二类是鲁棒优化不知道概率分布时用区间或盒式不确定集追求最坏情况下的安全但结果通常偏保守。第三类是分布鲁棒优化介于两者之间用模糊集描述分布不确定理论漂亮但实现复杂度高。如果只是做Matlab代码实现我个人推荐先用场景法。原因很实在场景法概念直观代码容易调试而且能够充分利用风电出力和负荷预测误差的历史统计信息。你生成500个初始场景削减到10个左右求解结果不仅有调度计划还能看到不同场景下的风电消纳情况和切负荷风险这对写分析报告非常有帮助。鲁棒优化虽然不用生成场景但是不确定集合边界怎么定、对偶约束怎么转化都要花不少功夫而且最终解的保守性很难向导师或甲方解释清楚。场景法里还有一个细节两阶段决策结构。第一阶段决策是机组启停和基本出力这些必须在看到实际风电和负荷之前确定第二阶段是弃风、切负荷和机组调整量这些可以根据场景实现后确定。在代码中第一阶段变量用binvar和sdpvar声明第二阶段变量只在对应场景约束里使用。目标函数里的期望成本就是对所有场景的第二阶段惩罚项做概率加权。这部分写代码时很容易犯错我后面会详细说。3. Matlab实现从场景生成到求解出图这一部分是实操量最大的地方我把自己的程序框架和核心代码思路都放出来。由于完整代码有上千行这里只保留最关键的逻辑片段。你完全可以按照这个框架自己搭一套数据换成自己系统的就行。3.1 程序架构与数据准备我习惯把程序分成四个模块参数设置、场景生成、模型求解、结果输出。参数设置模块包括机组参数、负荷预测曲线、风速/负荷不确定性参数、碳交易价格和配额等。机组参数建议放在表格里读入后方便修改。为了测试方便我用了6台火电机组加1个风电场的简化系统24小时调度周期初始场景数500削减后保留10个场景。机组参数表大致长这个样子机组最大出力/MW最小出力/MW煤耗系数a/($/MW²h)煤耗系数b/($/MWh)排放强度/(tCO2/MWh)爬坡速率/(MW/h)G1200500.0012300.7040G2150400.0018280.6535G3100250.0022320.7530G4100250.0020290.7230G580200.0025350.8025G680200.0024330.7825风电场额定容量300MW切入风速3m/s额定风速12m/s切出风速25m/s。负荷预测数据我直接取了一个典型日曲线这里不逐小时列出。碳价设为40元/吨免费配额按火电装机容量的一定比例分配实际排放超过配额的部分按碳价计入成本。3.2 用Matlab生成源荷联合场景的三个关键步骤第一步是生成初始场景。风速用Weibull分布抽样负荷误差用正态分布抽样。需要注意的是风速是时序相关的所以我会用马尔可夫链或者简单的一阶自回归模型来生成风速序列。自回归系数根据历史风速数据拟合通常设成0.85左右表示当前时段风速对下一时段影响较强。代码如下% 关键变量 N 500; % 初始场景数 T 24; % 调度时段 v_shape 2.0; % Weibull形状参数 v_scale 6.0; % Weibull尺度参数 rho 0.85; % 风速自回归系数 wind_scn zeros(N, T); load_base load_forecast; % 24x1负荷预测列向量 load_err_std 0.03; % 负荷预测标准差比例 load_scn zeros(N, T); % 生成时序相关风速场景 for i 1:N v zeros(1, T); for t 1:T innovation wblrnd(v_shape, v_scale); % Weibull噪声 if t 1 v(t) innovation; else v(t) rho * v(t-1) sqrt(1 - rho^2) * innovation; end end wind_scn(i, :) v; load_scn(i, :) load_base .* (1 normrnd(0, load_err_std, 1, T)); end第二步是把风速转成风电出力。这一步对应风机功率曲线可以用分段函数处理。切入风速以下或切出风速以上出力为0切入风速到额定风速之间近似线性上升额定风速到切出风速之间保持额定出力% 风机功率曲线 v_ci 3; v_r 12; v_co 25; P_r 300; wind_power zeros(N, T); idx_linear (wind_scn v_ci) (wind_scn v_r); wind_power(idx_linear) P_r .* (wind_scn(idx_linear) - v_ci) / (v_r - v_ci); idx_rated (wind_scn v_r) (wind_scn v_co); wind_power(idx_rated) P_r;第三步是考虑源荷相关性。严格的实现是两个随机变量通过协方差矩阵联合抽样但工程上可以先独立抽样再用Cholesky分解对已生成的风电出力和负荷误差做线性变换。要注意变换之后风电功率范围可能越界所以最后需要再做一次截断或重新归一化。这一步我一般放在场景削减之前避免把相关性破坏掉。3.3 场景削减同步回代消除法初始场景往往有几百个如果全塞进MIP变量数量和约束数量会剧增。我实测过500个场景、24时段、6台机组的目标函数展开之后YALMIP可能要卡上几分钟甚至报内存错误。场景削减的标准做法很多我比较推荐同步回代消除法它实现简单且效果稳定。核心思路是每次删掉一个场景并把被删场景的概率累加到离它最近的场景上使剩余场景集合与原始场景分布之间的概率距离增量最小。判断“最近”的距离可以用欧氏距离也可以把风速场景和负荷场景合并成一个高维向量再计算。代码如下scn [wind_power, load_scn]; % 合并源荷维度 probs ones(N, 1) / N; K_keep 10; % 削减后保留场景数 D pdist2(scn, scn, euclidean); D(1:N1:end) Inf; while size(scn, 1) K_keep min_val inf; for i 1:size(scn, 1) [v_j, idx_j] min(D(i, :)); inc probs(i) * v_j; if inc min_val min_val inc; del_i i; del_j idx_j; end end probs(del_j) probs(del_j) probs(del_i); probs(del_i) []; scn(del_i, :) []; D pdist2(scn, scn, euclidean); D(1:size(scn,1)1:end) Inf; end这个削减算法在实际运行时有个小坑pdist2每次重算距离矩阵如果场景数很大循环次数很多会非常慢。我处理的办法是先用K-means之类的粗聚类把场景缩到50个以下再用同步回代精确削减到10个。这样精度损失很小但速度能快一两个数量级。3.4 在YALMIP中建立优化模型并求解求解部分我用YALMIP做建模层求解器用Gurobi或者Cplex都行。如果没装商业求解器也能用SCS或OSQP跑部分模型但混合整数问题最好还是用商业求解器单纯形分支定界的性能差别很大。模型变量分成三组机组启停变量u二值变量24×6。机组出力变量P连续变量24×6。弃风变量和切负荷变量每个场景下24×1的连续变量。注意第二阶段变量是按场景设置的比如弃风变量就是24×6×场景数的大矩阵。YALMIP里我喜欢用细胞数组存这样约束循环写起来更清晰。核心建模代码如下% 变量定义 u binvar(nG, T); % 启停状态 P sdpvar(nG, T); % 机组出力 for s 1:K_keep curt{s} sdpvar(1, T); % 弃风 shed{s} sdpvar(1, T); % 切负荷 end % 目标函数 fuel_cost sum(sum(coe_a .* repmat(P.^2, [1 1]) coe_b .* P)); % 实际代码中二次项替换为分段线性 carbon_cost carbon_price * ... sum(sum(emission_factor .* P)) - free_quota; obj fuel_cost carbon_cost ... sum(probs .* (sum(curt{s}, 2) * penal_curt sum(shed{s}, 2) * penal_shed));约束添加的时候功率平衡约束写成每个场景下成立Constraints []; for s 1:K_keep Constraints [Constraints, ... sum(P, 1) wind_power_sce(s, :) - curt{s} load_sce(s, :) - shed{s}]; Constraints [Constraints, ... 0 curt{s} wind_power_sce(s, :)]; Constraints [Constraints, ... 0 shed{s} load_sce(s, :)]; end比较关键的一点是备用约束不要写成“所有场景都要满足”的确定性约束否则场景法就变成鲁棒优化了结果会非常保守。更合理的做法是让备用需求与负荷水平绑定并考虑风电出力预测误差的标准差reserve_up 0.05 * load_forecast 0.08 * (pred_wind_power - min_wind_power); Constraints [Constraints, ... sum(P_max .* u, 1) load_forecast reserve_up];然后调用求解器ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.01, gurobi.TimeLimit, 300); optimize(Constraints, obj, ops);MIPGap设成0.01就够了没必要追求严格最优。很多场景下相对误差从1%降到0.1%要花几倍时间但对工程结果没有本质区别。我第一次跑的时候把MIPGap设成0.0001结果一个小时都没跑完后来改成0.01两分钟内出解。3.5 结果分析与可视化求解之后我习惯画三张图。第一张是机组出力堆叠图横轴是24小时纵轴是各机组出力风电部分也堆叠上去。第二张是不同场景下的风电消纳对比图把削减后10个场景的消纳区间画出来能直观看到预测误差对风电消纳的影响。第三张是碳排放成本随碳价变化的敏感性曲线改变碳价从20到100元/吨看系统总成本和碳排放量的变化趋势。画图用Matlab自带bar和plot就行。要注意把YALMIP里sdpvar求解后的结果用value()取出来直接画会出错。另外堆叠图需要把风电放在最上面或最下面不然负值会导致图很难看。图里的坐标轴、图例都写好论文直接用这个图基本不用再加工。4. 跑代码时容易踩的坑和排查手册做这类调度模型最大的敌人不是数学推导而是模型写好了之后求解器给你一个“无解”或者“结果明显不合理”的输出。我把自己调试过程中遇到的典型问题整理成一张速查表希望能省掉你半天排查时间。现象常见原因排查方法提示无解或Infeasible功率平衡约束或备用约束过紧场景里有极端风电或负荷先固定场景逐个检查每个场景的约束是否可行把备用系数调小或增加惩罚变量求解时间过长场景数太多二次项没线性化MIPGap设得过严削减场景数二次成本改分段线性放松MIPGap出力结果出现负值机组变量定义中缺少下界约束检查P_min约束是否加了P P_min .* u弃风变量一直为0弃风惩罚设置太高或太低检查目标函数里弃风惩罚项是否真的被累加了边界条件是否正确场景削减后分布失真只用欧氏距离没考虑场景概率用概率距离或1-范数加权距离削减后重新归一化概率YALMIP报错“No solver available”没配置求解器路径或求解器未安装yalmiptest查看可用求解器确认Gurobi/Cplex的Matlab接口已配置碳交易成本为负碳配额设置高于实际排放碳价乘以配额确实可能为负这属于正常但如果长期为负说明配额过松调整配额系数4.1 场景数选多少合适场景数是个绕不开的权衡。初始场景太少分布覆盖不足削减后保留场景太少极端风电情况可能被削掉。我的经验值是初始场景300到500个削减后保留10到15个。如果你用了K-means粗聚类可以保留20个再精削减。这样既能保留典型极端情况又能让MIP在2分钟内收敛。如果场景里包含强相关性不要用随机抽样尽量用LHS或Sobol序列生成初始场景分布覆盖率会好很多。4.2 风电功率序列越界和负值问题风速转功率时最容易出现的就是功率曲线函数写得不够严谨导致风速低于切入风速时功率非零或者风速超过额定风速时功率超过额定容量。我建议转功率之后立即加一行截断wind_power max(0, min(P_r, wind_power));另外如果做了场景相关性变换变换后的功率值可能不在正常范围内也必须在变换后重新截断。但要注意截断会破坏一部分相关性所以协方差矩阵设计时要预留一点裕度不要让两个变量强相关到0.95以上否则截断之后相关性和预期差很远。4.3 碳交易参数的敏感性碳价和配额会影响优化结果的方向。碳价太低系统发现买碳配额比调整出力结构更便宜就不会主动减排风电消纳率可能不高。碳价太高系统会过度弃火电导致煤耗成本上升但碳排放确实下降。我建议在一开始先用一个中等碳价跑通模型然后做敏感性分析观察碳排放总量和总成本的帕累托前沿。如果你的场景里碳配额是按历史排放设置的记得每个场景的实际排放不同而配额通常是一个固定值不要不小心把配额写成场景相关变量。4.4 备用系数怎么取值备用需求可以用负荷比例加风电出力波动区间来表示。我常用的是上备用5%负荷8%风电装机容量×当前预测出力比例下备用5%负荷。这个系数没有统一标准要根据你系统的实际情况调。原则上备用系数越大模型越稳健但调度成本越高系数太低可能在极端场景下出现切负荷。可以先跑几组敏感性实验找到成本和风险的拐点。4.5 二次煤耗成本的处理很多教科书上煤耗成本用二次函数但直接放进MIP会带来凸二次项。Gurobi/Cplex能直接处理二次目标但和线性约束混在一起时求解速度会慢。实际做工程时我通常把机组出力区间分成三段用分段线性函数逼近二次成本曲线。分段点取在技术出力范围的1/3和2/3附近每段的斜率先算好然后用一组连续变量表示各段出力加上总和约束。这样目标函数变成线性求解速度明显加快精度损失不到1%。如果你的场景数很多这一步非常值得做。说了这么多最后聊一点个人体会。我一开始也迷信复杂的不确定性建模方法觉得鲁棒优化、分布鲁棒才是“高水平”。但真到了写代码、调参、出结果的时候发现把场景法做得扎实把源荷相关性处理好已经把95%的工程问题解决了。这个题目后续其实还有很多可以扩展的方向比如加入需求响应、碳捕集电厂、储能或者把调度尺度从日前延伸到日内滚动都是在这个框架上做加法。如果你也在做类似的东西建议先跑通这个基本版再一步一步加复杂度。模型不是越复杂越好能稳定出结果、能解释清楚才是眼下最要紧的。
