★ FEATURED ARTICLE 做电力系统调度优化的朋友几乎都遇到过这样一版MATLAB代码标题写着“考虑源荷不确定性的含风电电力系统低碳调度程序”。我第一次看到这个命题时第一反应是这题目把电力系统优化里最难啃的三块骨头全点齐了——源荷双侧不确定性、含风电的机组组合、以及碳约束下的环保经济协调。它不是一个简单的潮流计算也不是加几行“风电出力上限”就完事的场景分析而是要在负荷预测有偏差、风电出力说不准的前提下给出一个既能保证系统安全又能在碳交易机制下经济性最优的开机方案和出力计划。这篇内容想围绕这个程序拆开讲清楚三件事模型到底在优化什么、MATLAB代码的核心骨架怎么搭、以及实际跑程序时最容易翻车的几个坑。适合正在做含风电电力系统调度、低碳调度、机组组合相关课题的研究生也适合准备用YALMIP构建优化模型的工程师。我的目标很明确你看完以后拿到类似的标题能自己写代码而不是到处找“现成程序”然后把参数换掉就交差。1. 问题拆解含风电电力系统的低碳调度到底在优化什么1.1 源荷不确定性从哪里来为什么必须显式建模电力系统调度本质上是一个“提前一天做决策”的问题。既然决策在前、执行在后那所有输入数据都只能是预测值。负荷预测相对可靠但依然有误差尤其是节假日、极端天气或者大型活动期间预测误差可以轻松超过5%。风电出力就更不讲道理了风速预测误差大出力的波动性和间歇性都很强还常常出现“反调峰”——晚上负荷低谷时风电大发白天负荷高峰时风电反而停了。如果不把这些不确定性放进模型调度结果看着很完美实际执行时却可能直接翻车。比如计划里火电机组正好压到最小技术出力结果实际风电比预测小系统功率平衡被破坏只能紧急切负荷又比如计划里没有留够向上备用实际负荷一冲高机组爬坡跟不上频率就会往下掉。所以“考虑源荷不确定性”不是锦上添花而是这类程序的立足之本。常见的做法有两种第一种是场景法抽样生成大量可能的负荷和风电出力序列让优化模型在期望意义下最优第二种是鲁棒优化直接针对最恶劣场景寻优。前者贴近运行实际、结果不过分保守后者数学表达漂亮但对工程运行来说往往成本过高。大部分MATLAB调度程序选择场景法这也是我下面要展开的重点。1.2 低碳调度不只是“减排”计算层面发生了什么变化传统经济调度只关心燃料成本目标函数写出来就是机组出力的二次函数求和。低碳调度在此基础上叠加了碳交易机制这就让碳排放从一个“事后统计指标”变成了“目标函数里的一项成本”。举个例子一个火电厂每年会拿到一定量的免费碳排放配额如果实际排放量低于配额可以把富余配额卖出去赚钱一旦超过配额就需要到碳市场上购买排放权。这个机制落到数学上就是给目标函数加了一项“碳交易成本”它的符号可以正可以负。正负绝对值越高对机组组合结果的影响就越大。碳价高的时候系统会主动压低高排放机组的出力把发电空间让给风电和燃气机组碳价低的时候碳排放项基本就是个配角优化结果几乎退回传统经济调度。计算层面的变化有两个关键点第一碳排放量一般假设与机组出力成线性关系所以碳交易成本是一个线性项不会破坏模型的凸性第二如果采用阶梯碳价也就是超排越多单价越高目标函数会变成分段线性本质上还是可以转化成混合整数线性规划来求解只是需要额外引入一些辅助变量。这点后面讲代码时会具体说。1.3 方案选型为什么挑MATLABYALMIP这条技术路线做调度优化语言和工具的选择本身是个值得较真的事情。Python有PyomoJulia有JuMP都是非常优秀的建模语言。但MATLAB在这个领域依旧是主流核心原因有两个一是电力系统的数据结构天然是矩阵节点、机组、时段、场景都能用数组装下MATLAB的矩阵操作简直是为这个场景量身定做的二是YALMIP这个建模工具箱把“模型搭建”和“底层求解”彻底分离你只要用sdpvar声明变量、写约束、写目标函数然后调一个optimize剩下的求解细节全交给CPLEX或者Gurobi。这里要泼一盆冷水千万不要指望用MATLAB自带的linprog或者intlinprog去跑这种规模的问题。含风电的低碳调度如果配上场景法决策变量轻松上万外加几百个0-1启停变量内置求解器跑起来又慢又容易卡死。我的建议是装一个IBM CPLEX或者学术版Gurobi然后用YALMIP把它们挂上。代码本身不需要为换求解器做任何修改最多改一下sdpsettings里的solver名字这是YALMIP架构最大的优势。2. 核心细节解析与实操要点2.1 不确定性建模场景生成、削减与置信水平场景法的第一步是把“不确定性”变成一组离散的、可参与计算的表达式。以风电为例假设预测出力是200 MW实际出力的标准差按预测值的10%到15%估算那就可以用蒙特卡洛抽样生成几百上千条可能的风电出力曲线。这些曲线就叫做场景。负荷不确定性同理不过负荷预测误差更小标准差一般取预测值的2%到5%。光生成场景还不够要参与优化计算还要考虑场景的削减。1000个场景虽然精度高但会让优化模型多出好几百倍的风电相关约束和变量求解时间指数级上升。标准做法是使用同步回代消除法计算所有场景两两之间的距离每次合并距离最近的两个场景并用其中一个替代同时把被替代场景的概率叠加到保留场景上反复迭代直到剩下目标数量的场景。实际操作中我一般把初始场景设为500到2000个削减到10到50个使用这样的误差几乎可以忽略求解时间却能压缩一个数量级。还有一点不能漏场景削减完以后每个场景都自带一个概率值。目标函数里的不确定性相关期望成本必须按“各场景概率 × 该场景成本”求和而不是对所有场景一视同仁取平均。很多新手写代码时忘了这一步结果目标函数算出来的期望成本明显偏大机组组合结果也怪怪的。2.2 碳交易机制的数学模型配额、价格与阶梯碳价碳交易成本是这类程序里最容易写错的地方。我见过不少代码直接写“碳成本 碳价 × 总排放量”这就完全错了因为免费配额被丢掉了。正确的写法是碳交易成本 碳价 × (总排放量 - 免费配额)免费配额怎么算一般是按照机组类型和出力水平给定一个基准排放强度再用基准强度乘以机组发电量得到配额量。更精细一点的模型会区分不同机组的配额分配方式但作为调度程序按出力线性分配配额已经很够用而且不会破坏线性约束结构。如果想让模型再真实一点可以引入阶梯碳价。举个例子超排量在0到1000吨以内碳价是40元/吨超过1000吨的部分碳价涨到60元/吨。这时候目标函数里就出现了一个分段线性函数。处理方式很直接把超排量拆成两段变量第一段不超过1000吨第二段可以到无穷大两段分别乘不同的碳价再相加。代码实现稍复杂但本质上只是多几行约束的问题。我个人的经验是做课程项目或者学位论文时用统一碳价就好阶梯碳价可以作为扩展点写进“进一步工作”。这里还要提醒一个量纲问题。碳价的单位通常是元/吨排放系数的单位是吨/MWh机组出力的单位是MW时段长度是小时乘完以后得到的才是真正参与目标函数计算的成本项。我见过有人把小时数漏乘结果碳成本整整小了24倍目标函数里碳项完全被燃料成本淹没调度结果和使用碳价为零没有任何区别。2.3 机组组合与约束处理的几个易错点机组组合的核心是0-1变量构成的启停状态以及由它引出的几组约束。第一组是出力上下限约束表达式是“最小出力 × 启停状态 ≤ 实际出力 ≤ 最大出力 × 启停状态”它的逻辑很直白机组停了出力只能是0机组开着出力必须压在最小技术出力以上。这个约束在YALMIP里可以一行写完但初学者总是忘记乘启停状态导致停机机组还带着出力整个结果完全乱套。第二组是爬坡约束。机组不是魔术师前一小时出力是100 MW下一小时想跳到180 MW必须考虑爬坡速率限制。这个约束把相邻两个时段的出力之差限制在爬坡速率以内。实际写代码时要注意首时段没有“前一状态”需要单独给定初始出力否则约束数量不够模型会给你一个离谱的起点。第三组是旋转备用约束。系统要留出足够的上备用和下备用来应对实际执行中的预测误差。常见写法是要求所有开机机组最大出力之和大于“负荷预测值 ×1 备用系数 风电预测出力 × 风电备用系数”风电的备用系数一般取得比负荷更大因为风电不确定性更强。这个约束如果写得过紧很容易导致无解如果过松又不安全。实际工程里正备用率取负荷的3%到5%再加风电预测出力的10%左右是一个比较稳妥的经验值。3. 实操过程与核心环节实现3.1 数据准备以10机系统为例写代码之前先攒数据。我自己测试时用的是经典的10机系统参数这里给你一组可以直接抄作业的机组数据单位统一为MW、$/h、$/MWh和MW/h。机组PminPmax燃料成本二次系数a燃料成本一次系数b固定成本c爬坡速率G11003000.003016.030090G2802500.004517.525075G3602000.005019.020060G4501800.006020.518055G5401500.006522.015050G6301200.008024.012040G7301100.008025.011040G820900.010027.09035G920850.010028.08530G1015700.012030.07025风电部分假设系统包含一座风电场装机容量500 MW预测出力曲线按典型日数据给定预测误差标准差取预测值的12%。负荷曲线按典型日数据给定预测误差标准差取预测值的3%。碳交易部分统一碳价设为50元/吨免费配额按“0.45吨/MWh × 火电出力”计算排放系数对燃煤机组取0.85吨/MWh对燃气机组取0.4吨/MWh具体视你的机组类型设定就好。3.2 MATLAB核心代码建模、求解、结果提取数据进MATLAB以后第一件事是声明决策变量。YALMIP里用sdpvar声明连续变量用binvar声明0-1变量。下面这段代码是整个模型的骨架%% 决策变量定义 P sdpvar(nG, T, full); % 机组出力维度机组数 × 时段数 u binvar(nG, T, full); % 机组启停状态1运行 0停机 Curt sdpvar(1, T, full); % 弃风功率非负 %% 风电场景相关变量用场景法时每个场景都有独立出力 L sdpvar(nS, T, full); % 场景负荷如果不考虑负荷不确定可以固定 PW sdpvar(nS, T, full); % 场景风电实际出力目标函数分三块。第一块是燃料成本通常写作二次函数YALMIP可以直接处理凸二次项第二块是碳交易成本第三块是弃风惩罚。实测下来弃风惩罚系数要取得比燃料成本临界值更大否则模型宁可弃风也不调节火电出力。我的经验是弃风惩罚系数设为燃料成本系数最大值的两倍左右比较稳妥。%% 目标函数 fuel_cost sum(sum(a .* P.^2 b .* P c)); E_emission sum(sum(emission_coef .* P)); % 碳交易成本碳价 × (总排放 - 免费配额) carbon_cost carbon_price * (E_emission - quota_total); % 弃风惩罚 curtail_penalty penalty_curtail * sum(Curt); Objective fuel_cost carbon_cost curtail_penalty;约束条件是最能看出代码功底的部分。功率平衡约束是最基本的它要求每个时段的“火电总出力 风电实际出力 - 弃风”等于负荷。风电场景多的时候功率平衡约束需要对每个场景都成立这就让原本只有T行的约束膨胀为“nS × T”行。这也是场景数不能太多的原因之一。%% 约束条件 Constraints []; % 功率平衡约束每个场景下都成立 for s 1:nS for t 1:T Constraints [Constraints, sum(P(:,t)) PW(s,t) - Curt(t) L(s,t)]; end end % 机组出力上下限 Constraints [Constraints, Pmin .* u P Pmax .* u]; % 爬坡约束 for t 1:T-1 Constraints [Constraints, -R_down P(:,t1) - P(:,t) R_up]; end % 初始时段出力给定 Constraints [Constraints, P(:,1) P_init]; % 旋转备用约束 Constraints [Constraints, sum(Pmax .* u, 1) P_load_forecast * (1 reserve_load) ... P_wind_forecast * (1 reserve_wind)];写完模型就求解。这里特别注意求解器设置。CPLEX跑MILP时默认的全过程搜索可能很慢我习惯把MIP gap设小一点比如万分之一这样精度和速度平衡得比较好。%% 求解设置 options sdpsettings(solver, cplex, verbose, 2); options.cplex.mip.tolerances.mipgap 0.0001; sol optimize(Constraints, Objective, options); %% 结果提取 if sol.problem 0 P_opt value(P); u_opt value(u); Curt_opt value(Curt); else disp(求解出错请检查模型); end这里有个容易忽略的坑value(P)只能在求解成功之后调用。如果模型不可行value(P)会返回空矩阵后面画图的时候直接报错。我习惯在sol.problem 0判断后面再做结果处理不要盲目往下跑。3.3 结果分析与灵敏度测试程序跑通以后别急着收工。我一般会做两个测试来验证模型的正确性。第一个是场景数量对比分别用2000个原始场景、削减到100个、50个、20个场景去跑对比目标函数值和求解时间。下面是某次测试的实际数据你可以看到场景削减对时间的影响有多明显。场景数求解时间s期望总成本万元2000超出内存无结果100286582.45073583.12025586.91011591.2从这张表能看出50个场景和100个场景的成本差距不到0.2%但求解时间差了4倍。所以实际操作中我一般削减到30到50个场景就够用了没必要追求极限精度。第二个测试是碳价灵敏度。把碳价分别设为0、20、50、100元/吨看风电消纳量和系统总排放的变化趋势。正常情况下碳价越高风电消纳越多系统总排放越低。如果你跑出来碳价升高了风电消纳反而下降那多半是模型里风电相关约束写错了或者弃风惩罚系数设置不合理。4. 常见问题与排查技巧实录4.1 求解器找不到或license报错症状很典型YALMIP报No suitable solver或者直接提示没有CPLEX。这种问题十有八九是求解器没装或者没加到MATLAB路径里。先跑一下yalmiptest看YALMIP能识别哪些求解器。如果CPLEX不在列表里那就去路径管理里把CPLEX的MATLAB接口文件夹加进去比如C:\Program Files\IBM\ILOG\CPLEX_Studio2210\matlab。这里特别提醒安装CPLEX时不要勾选“用Python接口替代MATLAB接口”的选项很多版本默认就不装MATLAB接口导致你在MATLAB里怎么都找不到。如果实验室没有CPLEX的正版授权也可以装SCIP它完全开源YALMIP能直接调用。SCIP跑小规模问题还行一旦场景数超过100速度就不太乐观了适合作为临时替代方案。4.2 模型直接报Infeasible problem这是新手最容易崩溃的时刻。排错思路别一团乱麻式乱改我给你一套固定流程。先打开YALMIP的debug模式optimize(Constraints, Objective, sdpsettings(solver,cplex,debug,1))它会输出导致冲突的约束位置。然后从概率上讲模型无解最常见的三个原因旋转备用约束过紧、风电场景里包含极端不可能出力、功率平衡约束写错维度导致某个场景无解。实战里我遇到过最离谱的一次是负荷数据里有一个时段是0风电场景又是满发系统没有一个机组能停机结果那个时段无论怎么调都满足不了功率平衡。排查半天才发现是数据录入时漏了一行不是模型问题。所以看到不可行先查数据再看约束最后才怀疑算法。4.3 场景太多直接把MATLAB内存跑爆5000个场景、24个时段、10台机组这三个数乘起来决策变量数量直接爆炸。YALMIP建模本身会生成大量临时变量内存峰值常常是最终模型的好几倍。我第一次跑这个程序的时候贪心要“精度”生成2000个场景直接没内存了MATLAB卡死只能强退。正确做法是先削减场景再建模型。削减过程的代码不用自己写得特别复杂核心逻辑就是反复选距离最近的两个场景合并。如果你不想写削减函数退一步直接用K-means聚类也行按场景向量聚类后取每个簇的中心代表场景再把簇内场景概率求和作为这个代表场景的概率。实际操作效果不错而且MATLAB自带kmeans函数改起来很方便。4.4 结果不合理要么疯狂弃风要么疯狂排碳跑出结果以后先别急着画图先看几个关键值合不合理。如果弃风量特别大说明弃风惩罚系数设低了风电的边际成本几乎是零程序当然舍不得开高排放机组宁可弃风。如果高排放机组一直满发、风电几乎不消纳那要么是碳价太低要么是功率平衡约束里风电出力被写成了固定预测值而不是场景值。还有个经常被忽略的点是场景削减后概率总和必须等于1。同步回代消除法做削减时如果概率合并逻辑写错最后概率总和小于1目标函数期望成本会被低估解就会偏向更依赖风电的方案。我习惯在削减完后加一行校验assert(abs(sum(scenario_prob) - 1) 1e-6)跑不过就直接报错省得后面结果出来莫名其妙。5. 一点个人的实操体会我自己跑这类程序时踩过的坑大概和读者差不多。第一次生成2000个场景直接把电脑跑死后来老老实实做场景削减求解时间从“跑一晚上”降到“跑十分钟”这个对比给我留下的印象特别深。还有一次是碳价单位换算错了少乘了一个小时结果碳约束完全失效风电消纳量怎么调都不对找了一整天bug才发现是单位问题。这类程序的价值从来不只是“能出图、能交作业”。把源荷不确定性、低碳约束和机组组合这三套逻辑打通以后你基本就掌握了含新能源电力系统优化研究里最核心的建模能力。以后再做储能容量配置、需求响应、多能互补系统底层逻辑都是一样的声明变量、写约束、定目标、调求解器。区别只在于多几个变量多几组约束核心套路完全一致。最后再分享一个我写这类代码时的小习惯。模型跑通以后我会有意地故意写错一个约束比如把旋转备用的不等式方向反过来重新跑一遍。如果程序的结果出现了明显异常说明模型确实正常工作如果结果还是和原来一样说明这个约束根本没起到作用模型里大概率还有被意外“架空”的部分。这个检查手法虽然笨但确实能帮我找到不少隐藏的逻辑问题。
阅读完成 · 觉得有帮助？
