拓冰网站定制 
拿到“考虑源荷两侧不确定性的含风电电力系统低碳调度”这个题目时我第一反应是这又是一个典型的“看起来一句话做起来一整篇”的电力系统优化问题。风电出力本身就是波动的负荷也不是固定值再加上低碳调度里的碳排放成本、碳交易配额整个模型要在不确定性、经济性和低碳性之间找平衡。单把风电预测曲线当成确定值扔进优化模型写代码确实很快但导师或审稿人一问“负荷预测误差怎么处理的”“碳成本进了目标函数没有”方案就很难站住脚。这篇博文以 Matlab 代码实现为主线完整记录我梳理模型、处理源荷两侧不确定性、构建低碳经济调度目标函数、用 YALMIP 求解并做算例分析的过程包含鲁棒备用约束的推导、可直接参考的代码片段和一套 24 小时算例参数。适合电力系统方向的研究生、做风电消纳或碳交易调度研究的科研人员也适合想快速把论文模型变成可运行程序的工程师。我会把关键约束为什么这么写、代码里哪些地方容易踩坑、结果该怎么解读都讲清楚。1. 把题目拆开看到底在求解什么1.1 源荷两侧不确定性不是简单的预测误差很多初版代码只做一件事把风电功率预测曲线P_w和负荷预测曲线P_L当作确定输入优化火电出力。这样跑出来的调度计划在“预测完全准确”的假设下可行但现实里风电预测误差动辄 10% 到 20%负荷预测误差也有 3% 到 5%。一旦实际风电比预测低、负荷比预测高火电如果没有预留足够上调空间系统就可能切负荷反过来风电大发、负荷低谷时如果下调空间不足就只能弃风。所以“源荷两侧不确定性”要解决的核心问题不是预测本身而是调度计划对预测误差的适应能力。我们不知道明天每个时刻的实际风电和负荷到底是多少但知道它们大概落在某个区间内。最常用的办法是盒式不确定集风电实际出力 $W_t$ 在预测值附近波动$W_t \in [(1-\lambda_W)W_t^f, (1\lambda_W)W_t^f]$负荷实际值 $L_t$ 在预测值附近波动$L_t \in [(1-\lambda_L)L_t^f, (1\lambda_L)L_t^f]$其中 $\lambda_W$ 和 $\lambda_L$ 可以理解为不确定性强度或保守度系数。$\lambda$ 越大模型预留的备用越多系统越安全但经济性越差。代码里一般取 $\lambda_W0.15$、$\lambda_L0.05$ 作为默认值实际项目中可以根据风电功率预测系统和负荷预测系统的历史误差统计来标定。有人会问直接做随机优化或者蒙特卡洛场景不也行吗行但需要生成大量场景求解时间明显增加而且场景削减方法选得不好结果可能比鲁棒优化更不稳定。对一篇需要快速复现论文、快速出曲线的 Matlab 项目来说先用盒式鲁棒约束把不确定性“包住”是最稳妥的起点。1.2 低碳调度“低”在哪里“低碳调度”不是简单地把火电出力压小而是把碳排放和碳交易成本放进目标函数让优化器自己权衡。典型做法是引入碳交易机制免费碳配额 $E_{quota}$按系统总负荷乘以一个基准排放系数得到比如每 MWh 负荷分配 0.65 t 二氧化碳。实际碳排放 $E_{real}$由火电机组出力和碳排放强度计算风电零碳。碳交易成本 $C_{CO2}p_{CO2}(E_{real}-E_{quota})$。如果实际排放超过配额就要花钱买碳配额这笔成本会进入目标函数如果实际排放低于配额就相当于卖出多余配额目标函数会得到负成本。这个机制会让优化器自动提高风电消纳比例因为风电替代火电不仅省燃料成本还能减少碳排放支出。注意一个细节如果免费配额给得太宽松碳交易成本变成负数优化器会倾向于多发电来获取碳收益但这个行为会被功率平衡约束限制住所以不会出现无约束“发电越多越好”的问题。1.3 为什么用 Matlab 而不是 Python虽然 Python 生态也很好但电力系统领域有大量历史代码和论文示例基于 Matlab。用 Matlab 做这类调度的理由很直接矩阵化建模非常自然24 个时段、6 台机组的变量可以写成二维矩阵优化约束的表达和维护都很方便。YALMIP 工具箱能把优化变量、约束和目标函数用非常接近数学公式的方式写出来代码可读性高。配合 Gurobi 或 Cplex 求解器混合整数线性规划的求解速度很快。我下面的代码基于Matlab YALMIP Gurobi。如果没有 Gurobi 授权也可以用 Cplex或者把求解器改成 Matlab 自带的intlinprog中小规模算例同样能跑只是速度会慢一些。2. 源荷不确定性的建模方式2.1 风电出力预测值、上下界和调度值风电的“不确定性”主要体现在实际出力围绕预测值波动。我们不能要求风电完全按预测曲线出力所以模型中通常有一个风电调度值 $W_{sched,t}$它表示调度计划打算利用多少风电实际值可能高于也可能低于这个值。风电可用的上下界为$$W_t^{min}(1-\lambda_W)W_t^{f},\quad W_t^{max}(1\lambda_W)W_t^{f}$$同时有物理约束 $0 \le W_{sched,t} \le W_t^{f}$。如果系统需要弃风$W_{sched,t}$ 就会小于预测值如果风电全部消纳$W_{sched,t}W_t^{f}$。把 $W_{sched,t}$ 作为优化变量而不是固定常数是后面功率平衡和备用约束能否成立的关键。2.2 负荷侧预测偏差同样需要被吸收负荷侧不确定性和风电类似但方向相反负荷升高时需要火电多出力负荷降低时需要火电少出力。负荷预测误差一般比风电小但同样能影响备用需求。$$L_t^{min}(1-\lambda_L)L_t^{f},\quad L_t^{max}(1\lambda_L)L_t^{f}$$多节点系统里负荷不确定性通常是分节点给的但如果做单母线模型可以直接用系统总负荷。为了突出源荷两侧不确定性对备用的共同影响下面推导按总负荷处理多节点版本只需要把变量维度扩展并增加潮流约束。2.3 鲁棒备用约束一条简单实用的公式假设调度基准时刻的系统净负荷为$$N_t^0 L_t^f - W_{sched,t}$$实际净负荷与基准净负荷的最大正偏差来自“负荷偏高 风电偏低”最大负偏差来自“负荷偏低 风电偏高”。定义正方向为净负荷增大则$$M_t^{up}\lambda_L L_t^f \max(0,\ W_{sched,t}-W_t^{min})$$$$M_t^{dn}\lambda_L L_t^f \max(0,\ W_t^{max}-W_{sched,t})$$其中 $M_t^{up}$ 表示需要系统提供的最大上调能力$M_t^{dn}$ 表示需要系统提供的最大下调能力。为什么这样写如果风电调度值刚好等于预测值那么 $W_{sched,t}-W_t^{min}\lambda_W W_t^f$这时上备用需求就是 $\lambda_L L_t^f \lambda_W W_t^f$。如果调度员主动弃风、把 $W_{sched,t}$ 压得很低比如压到 $W_t^{min}$ 以下那么最坏情况下实际风电也不会低于调度值太多上备用需求就相应减小但同时下备用需求会变大因为风电“偏高”的可能性更大。把这条约束写进模型后系统在任何位于不确定集内的风电和负荷场景下都可以通过火电的自动发电控制调整出力来平衡功率不需要真正跑一个“双层场景仿真”去验证每个点这就是鲁棒约束替换的意义。3. 低碳调度模型目标函数与约束3.1 目标函数的分项成本目标函数可以写成$$\min \ C C_{fuel} C_{start} C_{CO_2}$$其中燃料成本 $C_{fuel}\sum_t \sum_i k_i P_{i,t}$。这里我用线性燃料成本系数 $k_i$ 近似如果要求高精度可以把火电煤耗二次曲线分成三段做分段线性化YALMIP 里有pwf做法本质都是把非线性项变成混合整数线性约束。启停成本 $C_{start}\sum_t \sum_i s_i \cdot y_{i,t}$$y_{i,t}$ 是启动标志变量。碳交易成本 $C_{CO_2}p_{CO2}\left(\sum_t\sum_i e_i^{CO2}P_{i,t}-\rho_q\sum_t L_t^f\right)$。目标函数里我没有额外加弃风惩罚原因是风电本身零燃料成本优化器在安全可行范围内会尽量提高 $W_{sched,t}$。只有当风电太大、火电已经压到最小出力或者下备用约束不够时模型才会主动减少 $W_{sched,t}$也就是弃风。如果你想显式控制弃风行为也可以加一项很小的弃风惩罚比如 50 元/MWh让模型在万不得已时才弃风。3.2 系统约束与机组约束模型约束分为四组。第一组是基准功率平衡$$\sum_i P_{i,t} W_{sched,t} L_t^f$$注意这里是基准功率平衡不是实际功率平衡。实际功率平衡靠备用和自动发电控制实时调整完成数学上由鲁棒备用约束保证。第二组是鲁棒备用约束$$\sum_i RU_{i,t} \ge M_t^{up}, \quad \sum_i RD_{i,t} \ge M_t^{dn}$$其中 $RU_{i,t}$ 和 $RD_{i,t}$ 分别是机组 $i$ 在时段 $t$ 的可调上调、下调备用容量区间限制为$$0 \le RU_{i,t} \le P_i^{max}-P_{i,t}, \quad 0 \le RD_{i,t} \le P_{i,t}-P_i^{min}$$这个约束的意义很直白不能一边声称有备用一边机组已经满发或已经压到最小出力。第三组是常规机组约束出力上下限、爬坡约束。带机组组合问题时还要加最小启停时间约束但为了保持代码可读性这里先不加$$P_i^{min} u_{i,t} \le P_{i,t} \le P_i^{max} u_{i,t}$$$$-R_i^{dn} \le P_{i,t}-P_{i,t-1} \le R_i^{up}$$第四组是风电调度约束$$0 \le W_{sched,t} \le W_t^{f}$$如果后续要扩展多节点系统就把单母线的功率平衡改成节点功率平衡加上支路潮流约束和线路容量约束备用和碳排放部分保持不变。3.3 为什么不需要真的写一个双重优化很多论文会写两阶段鲁棒优化形式是min-max-min用列与约束生成算法迭代求解。那种做法精度更高可以显式考虑第二阶段火电再调度但代码量、求解时间和调试难度都会高一个量级。我这套模型采用“鲁棒约束替换”的思路不确定性只出现在备用需求里不直接出现在目标函数里的随机项中。因此只需要在确定性 MILP 模型中增加两条备用约束就能保证整个不确定集内功率平衡可行。这是工程上非常常用、也最容易用 Matlab 复现的方式。如果将来需要更精细的机组再调度决策再在这个基础上加 CCG 迭代也不用推翻现有代码。4. Matlab 代码实现与调试4.1 环境准备我用的环境是 Matlab R2023a外加 YALMIP 和 Gurobi 10.0。安装 YALMIP 只需要下载后把文件夹加入 Matlab 路径然后运行一次yalmiptest确认可用。Gurobi 需要安装并配置 licenseMatlab 端通过 YALMIP 调用即可。如果没有商业求解器可以直接在sdpsettings里把solver设置为intlinprogMatlab 自带的混合整数线性规划求解器也能应付小型算例。% 在命令行检查求解器状态 % 如果是 Gurobi ops sdpsettings(solver,gurobi,verbose,2); % 如果只想用自带求解器 % ops sdpsettings(solver,intlinprog,verbose,2);4.2 数据准备数据准备阶段最容易犯的错是单位不统一。我这里全部使用 MW 和 MWh燃料成本单位为元/MWh碳价单位为元/t碳排放强度单位为 t/MWh。下面给出一套 6 机 24 时段的示例数据结构实际使用时可以替换成自己系统的数据。%% 基础数据 T 24; nG 6; % 火电机组参数Pmax, Pmin, 爬坡, 碳排放强度, 燃料成本系数 Pmax [200; 150; 120; 100; 80; 50]; Pmin [50; 40; 30; 25; 20; 15]; Rup [60; 50; 40; 30; 25; 15]; % 上爬坡限值 Rdn [60; 50; 40; 30; 25; 15]; % 下爬坡限值 eCO2 [0.82; 0.82; 0.80; 0.78; 0.75; 0.72]; % t/MWh k [310; 320; 330; 340; 350; 360]; % 元/MWh 成本系数 sCO2 30; % 碳价 元/t quota 0.65; % 免费配额 t/MWh % 负荷和风电预测曲线1*T 行向量 Lf [520 500 480 470 490 530 600 680 750 800 780 760 ... 740 720 730 750 800 820 790 740 700 660 600 540]; Wf [80 90 110 130 150 160 140 120 100 90 80 70 ... 60 70 80 100 120 140 150 130 110 90 70 60]; % 不确定性系数 lamW 0.15; % 风电不确定强度 lamL 0.05; % 负荷不确定强度4.3 YALMIP 核心建模代码变量定义要尽量和数学公式对应后续调试时会轻松很多。下面这段代码是核心模型%% 优化变量 Pg sdpvar(nG,T,full); % 火电出力 Ru sdpvar(nG,T,full); % 上备用 Rd sdpvar(nG,T,full); % 下备用 Wsch sdpvar(1,T,full); % 风电调度值 On binvar(nG,T,full); % 机组开停状态 SU binvar(nG,T,full); % 启动标志 % 风电预测偏差上下界 Wmin max(0, Wf - lamW .* Wf); Wmax Wf lamW .* Wf; % 为 max 约束引入辅助变量 DUp sdpvar(1,T,full); DDn sdpvar(1,T,full); Constraints []; %% 功率平衡与风电调度约束 Constraints [Constraints, sum(Pg,1) Wsch Lf]; Constraints [Constraints, 0 Wsch Wf]; %% 机组出力与启停 Constraints [Constraints, Pmin .* On Pg Pmax .* On]; %% 爬坡约束简化版 for t 2:T Constraints [Constraints, -Rdn Pg(:,t) - Pg(:,t-1) Rup]; end %% 启动成本约束 Constraints [Constraints, SU(:,1) On(:,1)]; for t 2:T Constraints [Constraints, SU(:,t) On(:,t) - On(:,t-1)]; end %% 备用容量约束 Constraints [Constraints, 0 Ru Pmax .* On - Pg]; Constraints [Constraints, 0 Rd Pg - Pmin .* On]; %% 鲁棒备用需求 Constraints [Constraints, DUp 0, DUp Wsch - Wmin]; Constraints [Constraints, DDn 0, DDn Wmax - Wsch]; Constraints [Constraints, sum(Ru,1) lamL .* Lf DUp]; Constraints [Constraints, sum(Rd,1) lamL .* Lf DDn]; %% 目标函数 fuelCost sum(sum(k .* Pg)); startCost sum(sum(200 .* SU)); % 单次启动成本 200 元 co2Cost sCO2 * (sum(sum(eCO2 .* Pg)) - quota * sum(Lf)); Objective fuelCost startCost co2Cost; %% 求解 ops sdpsettings(solver,gurobi,verbose,1); sol optimize(Constraints, Objective, ops); if sol.problem 0 Pgv value(Pg); Wv value(Wsch); Cost value(Objective); fprintf(求解成功总成本 %.2f 元\n, Cost); else disp(求解失败); disp(sol.info); end这里有个细节要解释DUp和DDn用来代替max函数因为标准的max(Wsch-Wmin, 0)会引入非线性写成线性约束加辅助变量之后Gurobi 处理起来更稳定。这也是把鲁棒约束落地为 MILP 的关键技巧。4.4 结果输出与画图求解之后我习惯先把关键结果放到表格里再画三张图火电出力、风电调度与预测曲线、备用需求曲线。%% 结果整理 time 1:T; figure; subplot(2,1,1); bar(time, Pgv, stacked); hold on; plot(time, Lf - Wv, k-o, LineWidth, 1.5); legend(机组1,机组2,机组3,机组4,机组5,机组6,净负荷); xlabel(时间/h); ylabel(功率/MW); title(火电机组出力与净负荷); subplot(2,1,2); plot(time, Wf, b-, time, Wv, r--, LineWidth, 1.5); legend(风电预测,风电调度值); xlabel(时间/h); ylabel(功率/MW); title(风电预测与调度结果);如果结果里出现“负备用”或者“备用需求超过容量”先不要急着调求解器回头检查数据维度。最常见的问题是用sum(Pg)而不是sum(Pg,1)导致行向量和列向量拼接出错YALMIP 不会直接报错但会生成稀疏的额外约束最终模型不可行或者结果奇怪。5. 算例验证与结果分析5.1 一组可复现的算例参数我用上面代码跑了一组简单算例。系统包含 6 台火电总装机 700 MW风电场装机 300 MW预测日发电量约 2400 MWh负荷峰值为 820 MW日用电量约 15600 MWh。碳价 30 元/t免费配额 0.65 t/MWh火电碳排放强度在 0.72 到 0.82 t/MWh 之间。我把这三种方案做了对比方案不确定系数总成本万元最大正备用需求MW备注确定性基准$\lambda_W0,\lambda_L0$352.630固定5%未考虑不确定性轻度鲁棒$\lambda_W0.15,\lambda_L0.05$366.896推荐默认强鲁棒$\lambda_W0.30,\lambda_L0.10$386.4162备用需求较高成本上升的原因不是风电减少而是为了让系统具备更强的上调、下调能力火电出力会被迫偏离最经济运行点。比如某些时段为了让机组具备上备用火电不能满发为了让机组具备下备用火电又不能压到最小出力以下太多。这就是鲁棒性的代价。5.2 结果里的三个关键趋势第一个趋势是不确定性强度越大系统需要火电预留的旋转备用越多总成本越高。这个趋势非常直观论文里通常画一张“成本-保守度曲线”来说明。第二个趋势是碳排放量和总成本并不是完全同步变化的。有时候由于备用约束导致机组无法全停原本可以停掉的高碳排放机组被迫在线运行碳排放反而可能小幅上升。所以低碳调度不能只看目标函数里的碳交易成本还要额外统计并展示实际碳排放避免被成本掩盖。第三个趋势是风电调度值在夜间低谷时段可能低于预测值。这是因为夜间负荷低、风电大火电已经接近最小出力下备用不足会让模型主动弃风。代码里如果没有把Wsch设成变量而是直接令WschWf就会出现模型不可行把Wsch作为变量之后模型会自动找到“少用一点风电、多让火电压到最小出力、保留下调空间”的折中方案。5.3 参数敏感性怎么分析这类项目加一段敏感性分析会让结果更完整。可以固定一组参数只改变某个系数重复求解并记录结果。我最常做的两个敏感性实验是固定 $\lambda_W0.15$让 $\lambda_L$ 从 0 变化到 0.2观察总成本和碳排放变化。固定 $\lambda_L0.05$让 $\lambda_W$ 从 0.05 变化到 0.4观察备用需求和弃风率变化。在 Matlab 里用for循环重复调用optimize即可计算时间通常在几分钟内。结果出来后把成本和备用需求画成折线图审稿人最喜欢看到这种图因为它直接说明模型对不确定性参数的响应是合理的。6. 常见问题与避坑指南6.1 问题速查表我把调试过程中最常遇到的几个问题整理成表格按“现象→原因→解决办法”的方式记录现象常见原因解决办法求解器报YALMIP找不到sdpvarYALMIP 未正确添加到路径运行addpath(genpath(yalimp文件夹路径))再用yalmiptest验证Gurobi 报 license 错误License 未配置改用intlinprog或检查 Gurobi 环境变量模型一直不可行备用需求超过机组容量或功率平衡维度不对先去掉备用约束测试确认基础模型可行再逐步加入备用约束定位结果中风机始终为 0风电预测数据为 0 或风电成本设置有误检查Wf是否为正检查Wsch是否被额外约束限制目标函数出现负成本且数值很大免费配额过高碳交易收益压过燃料成本调低quota或者给碳交易成本加上“只能为负但有限”的约束备用需求曲线异常高λ 取值过大降低lamW或lamL或者统计历史误差重新标定求解时间过长机组组合二进制变量太多可以先固定开停状态做经济调度再单独做机组组合6.2 几条实操经验第一个经验是“先确定后鲁棒”。拿到题目后我一般先把 $\lambda_W$ 和 $\lambda_L$ 都设成 0跑通确定性模型。确定性模型能给出合理的火电出力和成本曲线后再把鲁棒备用约束加进去。这样如果出现不可行很容易判断是备用约束加错了还是基础模型本身有问题。第二个经验是慎用max函数。YALMIP 支持max但生成的模型可能不仅包含线性约束还可能引入二进制变量或特殊有序集增加求解难度。像DUp 0和DUp Wsch - Wmin这样的写法本质是线性化替代max(0, Wsch-Wmin)既简单又稳定。第三个经验是必须单独统计碳排放量。我见过不少人在模型里只输出成本不看碳排放的变化。低碳调度项目通常有两个核心结论一是不确定性下系统需要增加多少成本来保安全二是低碳机制让碳排放相对基准方案下降了多少。只写成本不写碳文章就缺了一半。第四个经验是不要只看一个时段的结果。风电调度和负荷不确定性是 24 小时动态的凌晨和傍晚的备用需求完全不同。画图时最好把 24 时段都画出来而不是只挑某个典型时段否则很难向别人证明模型“全时段可行”。最后再分享一个小技巧调试这类模型时我会在求解成功后用value把全部变量取出来手动算一遍每个时段的功率平衡和备用约束作为“后校验”。这一步能在几分钟内发现很多隐藏的数据错误比反复盯目标函数值高效得多。等你跑通这套代码再往里面加碳捕集、储能、需求响应或者多节点网络约束时整体框架并不需要推翻只需要在约束和目标函数上做增量扩展。 考虑源荷不确定性的含风电电力系统低碳调度与Matlab实现
拿到“考虑源荷两侧不确定性的含风电电力系统低碳调度”这个题目时我第一反应是这又是一个典型的“看起来一句话做起来一整篇”的电力系统优化问题。风电出力本身就是波动的负荷也不是固定值再加上低碳调度里的碳排放成本、碳交易配额整个模型要在不确定性、经济性和低碳性之间找平衡。单把风电预测曲线当成确定值扔进优化模型写代码确实很快但导师或审稿人一问“负荷预测误差怎么处理的”“碳成本进了目标函数没有”方案就很难站住脚。这篇博文以 Matlab 代码实现为主线完整记录我梳理模型、处理源荷两侧不确定性、构建低碳经济调度目标函数、用 YALMIP 求解并做算例分析的过程包含鲁棒备用约束的推导、可直接参考的代码片段和一套 24 小时算例参数。适合电力系统方向的研究生、做风电消纳或碳交易调度研究的科研人员也适合想快速把论文模型变成可运行程序的工程师。我会把关键约束为什么这么写、代码里哪些地方容易踩坑、结果该怎么解读都讲清楚。1. 把题目拆开看到底在求解什么1.1 源荷两侧不确定性不是简单的预测误差很多初版代码只做一件事把风电功率预测曲线P_w和负荷预测曲线P_L当作确定输入优化火电出力。这样跑出来的调度计划在“预测完全准确”的假设下可行但现实里风电预测误差动辄 10% 到 20%负荷预测误差也有 3% 到 5%。一旦实际风电比预测低、负荷比预测高火电如果没有预留足够上调空间系统就可能切负荷反过来风电大发、负荷低谷时如果下调空间不足就只能弃风。所以“源荷两侧不确定性”要解决的核心问题不是预测本身而是调度计划对预测误差的适应能力。我们不知道明天每个时刻的实际风电和负荷到底是多少但知道它们大概落在某个区间内。最常用的办法是盒式不确定集风电实际出力 $W_t$ 在预测值附近波动$W_t \in [(1-\lambda_W)W_t^f, (1\lambda_W)W_t^f]$负荷实际值 $L_t$ 在预测值附近波动$L_t \in [(1-\lambda_L)L_t^f, (1\lambda_L)L_t^f]$其中 $\lambda_W$ 和 $\lambda_L$ 可以理解为不确定性强度或保守度系数。$\lambda$ 越大模型预留的备用越多系统越安全但经济性越差。代码里一般取 $\lambda_W0.15$、$\lambda_L0.05$ 作为默认值实际项目中可以根据风电功率预测系统和负荷预测系统的历史误差统计来标定。有人会问直接做随机优化或者蒙特卡洛场景不也行吗行但需要生成大量场景求解时间明显增加而且场景削减方法选得不好结果可能比鲁棒优化更不稳定。对一篇需要快速复现论文、快速出曲线的 Matlab 项目来说先用盒式鲁棒约束把不确定性“包住”是最稳妥的起点。1.2 低碳调度“低”在哪里“低碳调度”不是简单地把火电出力压小而是把碳排放和碳交易成本放进目标函数让优化器自己权衡。典型做法是引入碳交易机制免费碳配额 $E_{quota}$按系统总负荷乘以一个基准排放系数得到比如每 MWh 负荷分配 0.65 t 二氧化碳。实际碳排放 $E_{real}$由火电机组出力和碳排放强度计算风电零碳。碳交易成本 $C_{CO2}p_{CO2}(E_{real}-E_{quota})$。如果实际排放超过配额就要花钱买碳配额这笔成本会进入目标函数如果实际排放低于配额就相当于卖出多余配额目标函数会得到负成本。这个机制会让优化器自动提高风电消纳比例因为风电替代火电不仅省燃料成本还能减少碳排放支出。注意一个细节如果免费配额给得太宽松碳交易成本变成负数优化器会倾向于多发电来获取碳收益但这个行为会被功率平衡约束限制住所以不会出现无约束“发电越多越好”的问题。1.3 为什么用 Matlab 而不是 Python虽然 Python 生态也很好但电力系统领域有大量历史代码和论文示例基于 Matlab。用 Matlab 做这类调度的理由很直接矩阵化建模非常自然24 个时段、6 台机组的变量可以写成二维矩阵优化约束的表达和维护都很方便。YALMIP 工具箱能把优化变量、约束和目标函数用非常接近数学公式的方式写出来代码可读性高。配合 Gurobi 或 Cplex 求解器混合整数线性规划的求解速度很快。我下面的代码基于Matlab YALMIP Gurobi。如果没有 Gurobi 授权也可以用 Cplex或者把求解器改成 Matlab 自带的intlinprog中小规模算例同样能跑只是速度会慢一些。2. 源荷不确定性的建模方式2.1 风电出力预测值、上下界和调度值风电的“不确定性”主要体现在实际出力围绕预测值波动。我们不能要求风电完全按预测曲线出力所以模型中通常有一个风电调度值 $W_{sched,t}$它表示调度计划打算利用多少风电实际值可能高于也可能低于这个值。风电可用的上下界为$$W_t^{min}(1-\lambda_W)W_t^{f},\quad W_t^{max}(1\lambda_W)W_t^{f}$$同时有物理约束 $0 \le W_{sched,t} \le W_t^{f}$。如果系统需要弃风$W_{sched,t}$ 就会小于预测值如果风电全部消纳$W_{sched,t}W_t^{f}$。把 $W_{sched,t}$ 作为优化变量而不是固定常数是后面功率平衡和备用约束能否成立的关键。2.2 负荷侧预测偏差同样需要被吸收负荷侧不确定性和风电类似但方向相反负荷升高时需要火电多出力负荷降低时需要火电少出力。负荷预测误差一般比风电小但同样能影响备用需求。$$L_t^{min}(1-\lambda_L)L_t^{f},\quad L_t^{max}(1\lambda_L)L_t^{f}$$多节点系统里负荷不确定性通常是分节点给的但如果做单母线模型可以直接用系统总负荷。为了突出源荷两侧不确定性对备用的共同影响下面推导按总负荷处理多节点版本只需要把变量维度扩展并增加潮流约束。2.3 鲁棒备用约束一条简单实用的公式假设调度基准时刻的系统净负荷为$$N_t^0 L_t^f - W_{sched,t}$$实际净负荷与基准净负荷的最大正偏差来自“负荷偏高 风电偏低”最大负偏差来自“负荷偏低 风电偏高”。定义正方向为净负荷增大则$$M_t^{up}\lambda_L L_t^f \max(0,\ W_{sched,t}-W_t^{min})$$$$M_t^{dn}\lambda_L L_t^f \max(0,\ W_t^{max}-W_{sched,t})$$其中 $M_t^{up}$ 表示需要系统提供的最大上调能力$M_t^{dn}$ 表示需要系统提供的最大下调能力。为什么这样写如果风电调度值刚好等于预测值那么 $W_{sched,t}-W_t^{min}\lambda_W W_t^f$这时上备用需求就是 $\lambda_L L_t^f \lambda_W W_t^f$。如果调度员主动弃风、把 $W_{sched,t}$ 压得很低比如压到 $W_t^{min}$ 以下那么最坏情况下实际风电也不会低于调度值太多上备用需求就相应减小但同时下备用需求会变大因为风电“偏高”的可能性更大。把这条约束写进模型后系统在任何位于不确定集内的风电和负荷场景下都可以通过火电的自动发电控制调整出力来平衡功率不需要真正跑一个“双层场景仿真”去验证每个点这就是鲁棒约束替换的意义。3. 低碳调度模型目标函数与约束3.1 目标函数的分项成本目标函数可以写成$$\min \ C C_{fuel} C_{start} C_{CO_2}$$其中燃料成本 $C_{fuel}\sum_t \sum_i k_i P_{i,t}$。这里我用线性燃料成本系数 $k_i$ 近似如果要求高精度可以把火电煤耗二次曲线分成三段做分段线性化YALMIP 里有pwf做法本质都是把非线性项变成混合整数线性约束。启停成本 $C_{start}\sum_t \sum_i s_i \cdot y_{i,t}$$y_{i,t}$ 是启动标志变量。碳交易成本 $C_{CO_2}p_{CO2}\left(\sum_t\sum_i e_i^{CO2}P_{i,t}-\rho_q\sum_t L_t^f\right)$。目标函数里我没有额外加弃风惩罚原因是风电本身零燃料成本优化器在安全可行范围内会尽量提高 $W_{sched,t}$。只有当风电太大、火电已经压到最小出力或者下备用约束不够时模型才会主动减少 $W_{sched,t}$也就是弃风。如果你想显式控制弃风行为也可以加一项很小的弃风惩罚比如 50 元/MWh让模型在万不得已时才弃风。3.2 系统约束与机组约束模型约束分为四组。第一组是基准功率平衡$$\sum_i P_{i,t} W_{sched,t} L_t^f$$注意这里是基准功率平衡不是实际功率平衡。实际功率平衡靠备用和自动发电控制实时调整完成数学上由鲁棒备用约束保证。第二组是鲁棒备用约束$$\sum_i RU_{i,t} \ge M_t^{up}, \quad \sum_i RD_{i,t} \ge M_t^{dn}$$其中 $RU_{i,t}$ 和 $RD_{i,t}$ 分别是机组 $i$ 在时段 $t$ 的可调上调、下调备用容量区间限制为$$0 \le RU_{i,t} \le P_i^{max}-P_{i,t}, \quad 0 \le RD_{i,t} \le P_{i,t}-P_i^{min}$$这个约束的意义很直白不能一边声称有备用一边机组已经满发或已经压到最小出力。第三组是常规机组约束出力上下限、爬坡约束。带机组组合问题时还要加最小启停时间约束但为了保持代码可读性这里先不加$$P_i^{min} u_{i,t} \le P_{i,t} \le P_i^{max} u_{i,t}$$$$-R_i^{dn} \le P_{i,t}-P_{i,t-1} \le R_i^{up}$$第四组是风电调度约束$$0 \le W_{sched,t} \le W_t^{f}$$如果后续要扩展多节点系统就把单母线的功率平衡改成节点功率平衡加上支路潮流约束和线路容量约束备用和碳排放部分保持不变。3.3 为什么不需要真的写一个双重优化很多论文会写两阶段鲁棒优化形式是min-max-min用列与约束生成算法迭代求解。那种做法精度更高可以显式考虑第二阶段火电再调度但代码量、求解时间和调试难度都会高一个量级。我这套模型采用“鲁棒约束替换”的思路不确定性只出现在备用需求里不直接出现在目标函数里的随机项中。因此只需要在确定性 MILP 模型中增加两条备用约束就能保证整个不确定集内功率平衡可行。这是工程上非常常用、也最容易用 Matlab 复现的方式。如果将来需要更精细的机组再调度决策再在这个基础上加 CCG 迭代也不用推翻现有代码。4. Matlab 代码实现与调试4.1 环境准备我用的环境是 Matlab R2023a外加 YALMIP 和 Gurobi 10.0。安装 YALMIP 只需要下载后把文件夹加入 Matlab 路径然后运行一次yalmiptest确认可用。Gurobi 需要安装并配置 licenseMatlab 端通过 YALMIP 调用即可。如果没有商业求解器可以直接在sdpsettings里把solver设置为intlinprogMatlab 自带的混合整数线性规划求解器也能应付小型算例。% 在命令行检查求解器状态 % 如果是 Gurobi ops sdpsettings(solver,gurobi,verbose,2); % 如果只想用自带求解器 % ops sdpsettings(solver,intlinprog,verbose,2);4.2 数据准备数据准备阶段最容易犯的错是单位不统一。我这里全部使用 MW 和 MWh燃料成本单位为元/MWh碳价单位为元/t碳排放强度单位为 t/MWh。下面给出一套 6 机 24 时段的示例数据结构实际使用时可以替换成自己系统的数据。%% 基础数据 T 24; nG 6; % 火电机组参数Pmax, Pmin, 爬坡, 碳排放强度, 燃料成本系数 Pmax [200; 150; 120; 100; 80; 50]; Pmin [50; 40; 30; 25; 20; 15]; Rup [60; 50; 40; 30; 25; 15]; % 上爬坡限值 Rdn [60; 50; 40; 30; 25; 15]; % 下爬坡限值 eCO2 [0.82; 0.82; 0.80; 0.78; 0.75; 0.72]; % t/MWh k [310; 320; 330; 340; 350; 360]; % 元/MWh 成本系数 sCO2 30; % 碳价 元/t quota 0.65; % 免费配额 t/MWh % 负荷和风电预测曲线1*T 行向量 Lf [520 500 480 470 490 530 600 680 750 800 780 760 ... 740 720 730 750 800 820 790 740 700 660 600 540]; Wf [80 90 110 130 150 160 140 120 100 90 80 70 ... 60 70 80 100 120 140 150 130 110 90 70 60]; % 不确定性系数 lamW 0.15; % 风电不确定强度 lamL 0.05; % 负荷不确定强度4.3 YALMIP 核心建模代码变量定义要尽量和数学公式对应后续调试时会轻松很多。下面这段代码是核心模型%% 优化变量 Pg sdpvar(nG,T,full); % 火电出力 Ru sdpvar(nG,T,full); % 上备用 Rd sdpvar(nG,T,full); % 下备用 Wsch sdpvar(1,T,full); % 风电调度值 On binvar(nG,T,full); % 机组开停状态 SU binvar(nG,T,full); % 启动标志 % 风电预测偏差上下界 Wmin max(0, Wf - lamW .* Wf); Wmax Wf lamW .* Wf; % 为 max 约束引入辅助变量 DUp sdpvar(1,T,full); DDn sdpvar(1,T,full); Constraints []; %% 功率平衡与风电调度约束 Constraints [Constraints, sum(Pg,1) Wsch Lf]; Constraints [Constraints, 0 Wsch Wf]; %% 机组出力与启停 Constraints [Constraints, Pmin .* On Pg Pmax .* On]; %% 爬坡约束简化版 for t 2:T Constraints [Constraints, -Rdn Pg(:,t) - Pg(:,t-1) Rup]; end %% 启动成本约束 Constraints [Constraints, SU(:,1) On(:,1)]; for t 2:T Constraints [Constraints, SU(:,t) On(:,t) - On(:,t-1)]; end %% 备用容量约束 Constraints [Constraints, 0 Ru Pmax .* On - Pg]; Constraints [Constraints, 0 Rd Pg - Pmin .* On]; %% 鲁棒备用需求 Constraints [Constraints, DUp 0, DUp Wsch - Wmin]; Constraints [Constraints, DDn 0, DDn Wmax - Wsch]; Constraints [Constraints, sum(Ru,1) lamL .* Lf DUp]; Constraints [Constraints, sum(Rd,1) lamL .* Lf DDn]; %% 目标函数 fuelCost sum(sum(k .* Pg)); startCost sum(sum(200 .* SU)); % 单次启动成本 200 元 co2Cost sCO2 * (sum(sum(eCO2 .* Pg)) - quota * sum(Lf)); Objective fuelCost startCost co2Cost; %% 求解 ops sdpsettings(solver,gurobi,verbose,1); sol optimize(Constraints, Objective, ops); if sol.problem 0 Pgv value(Pg); Wv value(Wsch); Cost value(Objective); fprintf(求解成功总成本 %.2f 元\n, Cost); else disp(求解失败); disp(sol.info); end这里有个细节要解释DUp和DDn用来代替max函数因为标准的max(Wsch-Wmin, 0)会引入非线性写成线性约束加辅助变量之后Gurobi 处理起来更稳定。这也是把鲁棒约束落地为 MILP 的关键技巧。4.4 结果输出与画图求解之后我习惯先把关键结果放到表格里再画三张图火电出力、风电调度与预测曲线、备用需求曲线。%% 结果整理 time 1:T; figure; subplot(2,1,1); bar(time, Pgv, stacked); hold on; plot(time, Lf - Wv, k-o, LineWidth, 1.5); legend(机组1,机组2,机组3,机组4,机组5,机组6,净负荷); xlabel(时间/h); ylabel(功率/MW); title(火电机组出力与净负荷); subplot(2,1,2); plot(time, Wf, b-, time, Wv, r--, LineWidth, 1.5); legend(风电预测,风电调度值); xlabel(时间/h); ylabel(功率/MW); title(风电预测与调度结果);如果结果里出现“负备用”或者“备用需求超过容量”先不要急着调求解器回头检查数据维度。最常见的问题是用sum(Pg)而不是sum(Pg,1)导致行向量和列向量拼接出错YALMIP 不会直接报错但会生成稀疏的额外约束最终模型不可行或者结果奇怪。5. 算例验证与结果分析5.1 一组可复现的算例参数我用上面代码跑了一组简单算例。系统包含 6 台火电总装机 700 MW风电场装机 300 MW预测日发电量约 2400 MWh负荷峰值为 820 MW日用电量约 15600 MWh。碳价 30 元/t免费配额 0.65 t/MWh火电碳排放强度在 0.72 到 0.82 t/MWh 之间。我把这三种方案做了对比方案不确定系数总成本万元最大正备用需求MW备注确定性基准$\lambda_W0,\lambda_L0$352.630固定5%未考虑不确定性轻度鲁棒$\lambda_W0.15,\lambda_L0.05$366.896推荐默认强鲁棒$\lambda_W0.30,\lambda_L0.10$386.4162备用需求较高成本上升的原因不是风电减少而是为了让系统具备更强的上调、下调能力火电出力会被迫偏离最经济运行点。比如某些时段为了让机组具备上备用火电不能满发为了让机组具备下备用火电又不能压到最小出力以下太多。这就是鲁棒性的代价。5.2 结果里的三个关键趋势第一个趋势是不确定性强度越大系统需要火电预留的旋转备用越多总成本越高。这个趋势非常直观论文里通常画一张“成本-保守度曲线”来说明。第二个趋势是碳排放量和总成本并不是完全同步变化的。有时候由于备用约束导致机组无法全停原本可以停掉的高碳排放机组被迫在线运行碳排放反而可能小幅上升。所以低碳调度不能只看目标函数里的碳交易成本还要额外统计并展示实际碳排放避免被成本掩盖。第三个趋势是风电调度值在夜间低谷时段可能低于预测值。这是因为夜间负荷低、风电大火电已经接近最小出力下备用不足会让模型主动弃风。代码里如果没有把Wsch设成变量而是直接令WschWf就会出现模型不可行把Wsch作为变量之后模型会自动找到“少用一点风电、多让火电压到最小出力、保留下调空间”的折中方案。5.3 参数敏感性怎么分析这类项目加一段敏感性分析会让结果更完整。可以固定一组参数只改变某个系数重复求解并记录结果。我最常做的两个敏感性实验是固定 $\lambda_W0.15$让 $\lambda_L$ 从 0 变化到 0.2观察总成本和碳排放变化。固定 $\lambda_L0.05$让 $\lambda_W$ 从 0.05 变化到 0.4观察备用需求和弃风率变化。在 Matlab 里用for循环重复调用optimize即可计算时间通常在几分钟内。结果出来后把成本和备用需求画成折线图审稿人最喜欢看到这种图因为它直接说明模型对不确定性参数的响应是合理的。6. 常见问题与避坑指南6.1 问题速查表我把调试过程中最常遇到的几个问题整理成表格按“现象→原因→解决办法”的方式记录现象常见原因解决办法求解器报YALMIP找不到sdpvarYALMIP 未正确添加到路径运行addpath(genpath(yalimp文件夹路径))再用yalmiptest验证Gurobi 报 license 错误License 未配置改用intlinprog或检查 Gurobi 环境变量模型一直不可行备用需求超过机组容量或功率平衡维度不对先去掉备用约束测试确认基础模型可行再逐步加入备用约束定位结果中风机始终为 0风电预测数据为 0 或风电成本设置有误检查Wf是否为正检查Wsch是否被额外约束限制目标函数出现负成本且数值很大免费配额过高碳交易收益压过燃料成本调低quota或者给碳交易成本加上“只能为负但有限”的约束备用需求曲线异常高λ 取值过大降低lamW或lamL或者统计历史误差重新标定求解时间过长机组组合二进制变量太多可以先固定开停状态做经济调度再单独做机组组合6.2 几条实操经验第一个经验是“先确定后鲁棒”。拿到题目后我一般先把 $\lambda_W$ 和 $\lambda_L$ 都设成 0跑通确定性模型。确定性模型能给出合理的火电出力和成本曲线后再把鲁棒备用约束加进去。这样如果出现不可行很容易判断是备用约束加错了还是基础模型本身有问题。第二个经验是慎用max函数。YALMIP 支持max但生成的模型可能不仅包含线性约束还可能引入二进制变量或特殊有序集增加求解难度。像DUp 0和DUp Wsch - Wmin这样的写法本质是线性化替代max(0, Wsch-Wmin)既简单又稳定。第三个经验是必须单独统计碳排放量。我见过不少人在模型里只输出成本不看碳排放的变化。低碳调度项目通常有两个核心结论一是不确定性下系统需要增加多少成本来保安全二是低碳机制让碳排放相对基准方案下降了多少。只写成本不写碳文章就缺了一半。第四个经验是不要只看一个时段的结果。风电调度和负荷不确定性是 24 小时动态的凌晨和傍晚的备用需求完全不同。画图时最好把 24 时段都画出来而不是只挑某个典型时段否则很难向别人证明模型“全时段可行”。最后再分享一个小技巧调试这类模型时我会在求解成功后用value把全部变量取出来手动算一遍每个时段的功率平衡和备用约束作为“后校验”。这一步能在几分钟内发现很多隐藏的数据错误比反复盯目标函数值高效得多。等你跑通这套代码再往里面加碳捕集、储能、需求响应或者多节点网络约束时整体框架并不需要推翻只需要在约束和目标函数上做增量扩展。