这两年凡是做能源优化调度的几乎都绕不开“经济-碳协调”这四个字。我最初接触这类问题的时候想法很简单给成本函数里加一个碳排放惩罚项不就行了真动手在 IEEE33 节点配电网模型上做综合能源系统最优调度才发现事情远没那么轻松。不同设备的时间尺度、电热气的耦合关系、碳交易价格对调度策略的牵引任何一个环节处理不好结果就会偏离实际很远。这篇内容是我基于一个综合能源系统调度项目整理出来的完整思路重点解决三件事第一如何建立经济-碳双目标下的调度模型第二如何在 IEEE33 节点系统上用 Matlab 实现求解第三通过灵敏度分析判断碳价、负荷波动、新能源渗透率对调度结果的影响。整套方法是可复现的尤其适合做园区微网、区域综合能源系统、虚拟电厂、配电网低碳运行方向的学生和工程师参考。我尽量把所有模型细节、代码结构和踩过的坑都说清楚你可以把它当成一套带注释的实践笔记来用。1. 问题到底在算什么经济-碳协调的技术拆解1.1 为什么单纯“成本最低”不够用很多人一开始会困惑既然传统经济调度已经在做为什么还要单独提“经济-碳协调”因为这两者在物理层面是存在冲突的。举个例子向上级电网买电很便宜但电网电力的碳排放因子很高你买得越多系统整体碳排放就越高。天然气热电联产机组单位发电成本高一些排放却要低不少。光伏和风电属于零碳电源但出力的随机性又强不能让系统完全依赖它们。所以在不引入碳排放相关约束或成本项之前优化算法一定会偏向最便宜的电源组合结果往往是碳排放量居高不下。反过来说如果只追求碳排放最小化会让系统拼命用昂贵的清洁能源总运行成本又会高得离谱。“经济-碳协调”的核心不是在两个极端里二选一而是把碳排放量通过碳交易机制、碳税或者多目标权重的方式转成一种可量化的成本让优化算法自己去权衡。这样得到的调度方案既不是经济最优也不是碳最优而是综合代价最低的折中结果。我在项目里采用的是碳交易机制的表达方式。模型会先给系统一个免费碳排放配额实际排放超过配额就需要购买碳配额产生碳交易成本如果实际排放低于配额则可以把富余配额出售获得收益。这样碳排放就从抽象的环境指标变成了目标函数里的一个显式项。1.2 为什么选 IEEE33 节点作为验证网络IEEE33 节点配电网是电力系统领域用得最多的辐射状配网测试系统之一。它由 33 个节点、32 条支路组成馈线电压等级 12.66kV系统总有功峰值约 3.715MW、总无功峰值约 2.3Mvar。这个系统最大的优势在于公开、标准、可对比几乎任何一篇配电系统优化论文都可以用同一套数据做验证。你可能觉得 33 节点规模太小但我想强调一点做调度算法研究时网络规模不是越大越好关键要看模型是否能暴露核心问题。IEEE33 节点有足够的拓扑复杂度让你研究电压约束、线路容量约束和多接入点协调同时规模又可控用 Matlab 跑模型时迭代速度快做灵敏度分析时可以批量多次求解。如果你直接用 IEEE123 节点甚至更大规模的网络单次求解时间会明显拉长反而不利于做大范围的参数扫描实验。在我的配置中电力主干就是标准 IEEE33 节点辐射网再将综合能源部分作为设备接入包括分布式风电、光伏、微燃机热电联产机组、燃气锅炉和电池储能。这样做的好处是电力潮流约束落在公开的配电网拓扑上热力部分则用热功率平衡来简化既保留了多能流耦合特征又不至于把模型搞得过于庞杂。1.3 项目的完整技术路线整体路线可以拆成四条线建模、求解、分析、验证。建模构建以总运行成本与碳交易成本之和最小为目标的目标函数决策变量包括各机组出力、储能充放电功率、向上级电网购电功率、CHP 热电出力等。求解把模型写成 Matlab YALMIP 的形式用 CPLEX 或 Gurobi 这类商用求解器解混合整数线性规划问题因为模型中含有机组启停、储能充放状态等 0-1 变量。分析对碳交易价格、免费配额比例、负荷预测偏差、新能源渗透率等关键参数做灵敏度分析观察调度方案的变化趋势。验证在最优调度结果基础上对关键时段的节点电压和支路潮流做一次校验确认调度方案在物理上可行。后面所有章节基本就是沿着这四条线展开的。2. 数学模型是怎么一步步搭起来的2.1 目标函数把碳排放翻译成货币成本整个模型的目标函数是一个标量最小化问题我把它写成下面这种形式min C_total C_fuel C_grid C_om C_carbonC_fuel天然气燃料成本主要来自微燃机 CHP 和燃气锅炉。C_grid向上级电网购电的成本按分时电价计算。C_om各设备的运行维护成本通常和设备出力成线性关系。C_carbon碳交易成本。如果系统实际碳排放量超过了免费配额这项为正如果低于配额这项为负相当于卖碳收益。碳交易成本的具体表达式是C_carbon lambda_c * (E_total - E_quota)其中 lambda_c 是碳交易价格E_total 是系统实际碳排放总量E_quota 是免费配额。这里有一个关键细节值得展开E_total 主要包含两部分碳源。第一部分是向上级电网购电对应的间接碳排放用购电量和电网碳排放因子相乘得到第二部分是天然气燃烧产生的直接排放用天然气消耗量乘排放因子得到。要注意如果系统里有 P2G 这种把电转化为气的设备还需要在碳流计算中扣掉这部分转化过程吸收的碳量否则会出现重复计算。我在基础模型里没有引入 P2G只保留燃气设备和外购电两个碳源这样模型更干净灵敏度分析结论也更容易解读。免费配额 E_quota 的设定方式在很大程度上决定了调度结果。常见做法是按系统的总电负荷乘以一个基准排放强度来给配额比如“每 MWh 电负荷发放 0.55 吨二氧化碳配额”。这个基准值如果取得高相当于对系统很宽松碳排放约束几乎不起作用取得低则等于强制系统使用低碳设备。所以做灵敏度分析时配额比例和碳价是同等重要的两个关键参数。2.2 约束条件全梳理目标函数只是问题的“发动机”真正让模型变得复杂的是约束条件。一个完整的经济-碳协调调度模型至少要包含下面几组约束多能平衡约束电功率平衡表示每个时段所有电源出力加上储能放电、电网购电必须等于所有电负荷加上储能充电、电锅炉等用电设备消耗的功率。热功率平衡则是 CHP 余热回收和燃气锅炉产热之和等于热负荷需求。这类约束是等式约束也是模型里最基础的部分。机组约束微燃机 CHP 有功出力有上下限并且有爬坡速率限制。因为同一台 CHP 同时产电和产热运行时还存在一个热电比运行区间不能像独立发电机一样随意调节。燃气锅炉相对简单只需要限制出力范围。储能约束这个最容易被忽略也最容易出问题。储能电池要满足三组关系一是容量和能量之间的关系即 SOC 递推方程二是充放电功率上限三是每个时段只能充电或者只能放电的互斥约束。后者必须引入两个 0-1 变量如果忘了加互斥条件求解器很可能让电池在同一时刻既充电又放电产生无效解甚至无解。此外通常还要求调度周期结束时 SOC 回到初始值否则优化算法会把电池能量在最后一个时段全部放光得到一个不可持续运行的“短期最优解”。配电网潮流约束这是让模型从“能量平衡计算”升级为“调度决策”的关键。我在模型里采用 DistFlow 分支流模型来描述辐射状配电网潮流。DistFlow 的核心思想是对每条支路写出有功、无功和电压降落的递推关系。在辐射型网络中每条支路的末端节点功率等于该节点下游所有节点负荷与注入功率之和。完整 DistFlow 是非线性的因为式中包含支路电流平方项和电压平方项。如果做严格的凸松弛可以转成二阶锥规划求解这样能保证得到全局最优解。但我做第一版复现时用了线性化 DistFlow做法是假设系统各节点电压都在 1.0 p.u. 附近忽略网损项从而把潮流约束变成线性等式。这样处理后整个模型就是一个标准的 MILPCPLEX、Gurobi 求解速度非常快。求解完成后再用完整潮流计算去做校验确认简化后的结果在电气上仍然合理。2.3 灵敏度分析要回答三类实际问题如果说最优调度回答的是“给定条件下怎么做最划算”那灵敏度分析回答的就是“条件变了方案该怎么调整”。具体到这套系统灵敏度分析主要分三条线第一条线是碳交易价格变化。碳价提高后外购电力的间接碳成本上升系统会主动减少外购电、增加本地燃气机组出力和新能源消纳。这个方法可以用来判断一个园区或配电网对碳价信号的响应强度。第二条线是负荷波动包括电负荷和热负荷的预测偏差。调度方案对负荷偏差有多敏感直接关系到运行人员在日前计划基础上需不需要做日内调整。第三条线是新能源渗透率变化。把风电、光伏的装机容量按比例上调观察系统碳排放下降的边际效果。当渗透率增加到一定水平后由于弃风弃光矛盾加剧碳排放下降速度会趋缓这时候就能找到“最佳新能源配置比例”的参考依据。除此之外免费配额比例通常也会作为灵敏度的辅助分析对象。不过它的作用和碳价几乎是对称的配额比例降低、碳价提高效果都是促使系统低碳化。两个参数同时变化时甚至会呈现很强的等价性。3. Matlab 实现从数据准备到求解器调用3.1 工具链选型YALMIP 求解器的搭配逻辑整套代码我之前是用 Matlab 来写的。先说明一下我的环境配置Matlab R2022a 以上版本安装了 YALMIP 工具箱求解器用的是 Gurobi备用选项是 CPLEX 或者开源的 CBC。YALMIP 本身不是一个求解器而是一层建模语言可以帮你把约束和目标函数翻译成求解器能识别的标准格式。做这类问题我不建议直接用 Matlab 自带的 linprog 或 intlinprog 去逐行建模。虽然 intlinprog 勉强能处理小规模 MILP但一旦模型里出现几十个变量和几百条约束写代码的效率就变得很低。YALMIP 允许你像写数学表达式一样直接在代码里写约束比如 0 P_chp 0.6 这样的关系它内部会自动完成变量匹配和矩阵拼接极大减少出错概率。求解器层面Gurobi 和 CPLEX 是商用软件里对 MILP 支持最好、速度最快的两个。学生和科研用户通常可以申请免费 license。如果没有条件安装商业求解器也可以先用 CBC 或 GLPK 顶着用只是求解大规模模型时会明显变慢。3.2 IEEE33 节点数据怎么组织IEEE33 节点系统在 Matlab 里通常用一个支路表和一个节点负荷表来描述。支路表包含起始节点、终止节点、支路电阻和电抗还有支路容量上限节点负荷表则是每个节点的有功和无功负荷。我在代码里是把这些数据按有名值存放的电阻单位为欧姆功率单位为 kW这样后续生成潮流约束时直观一些。一个容易踩的坑是节点编号方式。IEEE33 节点标准的头文件里支路数据的第一行是从节点 1 到节点 2但很多 Matlab 版本里矩阵索引是从 1 开始的而某些第三方数据文件又习惯从 0 开始编号。只要有一个节点对不齐后面的拓扑关系全部错位调度结果看起来可能有数值实际物理上完全不成立。我第一次跑的时候就因为这个原因导致电压分布完全异常。我的建议是在做潮流校验之前先单独拿支路数据计算一遍网络拓扑连通性确认 33 个节点全部联通、根节点编号正确再进入优化模块。3.3 核心建模代码骨架我用 YALMIP 建模时代码的大致骨架如下。先定义每个时段的决策变量% 定义变量时间为1到24小时 P_chp sdpvar(1, 24, full); % CHP电出力单位kW Q_chp sdpvar(1, 24, full); % CHP热出力单位kW P_gb sdpvar(1, 24, full); % 燃气锅炉热出力单位kW P_buy sdpvar(1, 24, full); % 向上级电网购电功率 P_chg sdpvar(1, 24, full); % 储能充电功率 P_dis sdpvar(1, 24, full); % 储能放电功率 SOC sdpvar(1, 25, full); % 储能荷电状态多一个端点便于约束 u_chp binvar(1, 24); % CHP启停状态 u_chg binvar(1, 24); % 充电状态标志 u_dis binvar(1, 24); % 放电状态标志然后再写约束。储能是约束最多的一块我把关键约束摘出来说明Constraints []; for t 1:24 % 充放电互斥同一时刻只能处于一种状态 Constraints [Constraints, u_chg(t) u_dis(t) 1]; % 充放电功率上下限与状态变量挂钩 Constraints [Constraints, 0 P_chg(t) P_max * u_chg(t)]; Constraints [Constraints, 0 P_dis(t) P_max * u_dis(t)]; % SOC递推关系 Constraints [Constraints, SOC(t1) SOC(t) eta_chg*P_chg(t)*dt - P_dis(t)*dt/eta_dis]; % SOC上下限 Constraints [Constraints, SOC_min SOC(t1) SOC_max]; end % 调度周期结束时回到初始SOC防止把电池放空 Constraints [Constraints, SOC(25) SOC(1)];由于储能充放电效率的存在放电功率 P_dis 表示实际输出到母线的功率而不是电池内部消耗的能量所以递推公式里充电项要乘效率、放电项要除效率。这个细节如果写反了储能会像一个永远能产出多余能量的“永动机”而且优化器很可能会利用这个漏洞来压低成本。对于配电网的线性化 DistFlow 约束它们在 YALMIP 里是逐条支路、逐时段循环写出来的。以一条连接节点 i 和节点 j 的支路为例节点 j 的注入有功可以写成% 以支路末端节点为对象功率平衡 % 这里需要提前建立支路-末端节点映射表 branch_end for t 1:24 for k 1:32 j branch_end(k); % 节点注入功率 该节点上级支路流入功率 - 下级支路流出功率 % 在此基础上去掉线性化网损项即可 P_inj(j,t) P_from_branch(k,t) - P_to_branch(k,t); end end严格说完整潮流约束会比这段示意代码复杂很多因为你得根据拓扑先算出每个节点对应的父支路和子支路集合然后在循环里拼装。实际项目我建议把“配电网潮流模块”单独写成一个函数输入各节点净注入功率输出支路潮流和节点电压这样便于单独测试。3.4 求解配置与输出组织约束组装完成后目标函数也得按 24 个时段累加。YALMIP 里最方便的做法是把成本函数直接写成向量内积的形式。比如购电成本可以写C_grid price_buy(:) * P_buy(:);这里的 price_buy 是 24×1 的分时电价向量。天然气成本、运维成本、碳交易成本都可以写成类似形式。最后调用求解器ops sdpsettings(solver, gurobi, verbose, 2); result optimize(Constraints, Objective, ops);求解完成后直接用 value() 取出各变量的数值保存成结构体或者写入 Excel方便后续绘图和做灵敏度分析。我习惯把每个工况的结果存成一个 CSV 文件文件名就是该工况的关键参数比如 carbon_price_30_quota_0.55.csv。这样批量扫描时不容易搞混数据。4. 灵敏度分析怎么做才有说服力4.1 基准场景先跑通后面才有可比性做灵敏度分析的第一原则是先固定一个基准场景然后一次只改变一个参数。如果同时动碳价、负荷、光伏出力最后结果出现变化了你根本说不清是哪个因素引起的。我在项目里设置的基准场景是这样的仿真周期为冬季典型日 24 小时调度步长为 1 小时风电装机 1.2MW光伏装机 0.3MW微燃机 CHP 额定电功率 600kW燃气锅炉额定热功率 800kW储能容量 1MWh最大充放电功率 250kW初始 SOC 为 0.5基准碳交易价格设为 60 元/吨免费碳配额按总电负荷乘以 0.55 吨/MWh 发放。这个基准场景下先跑一次优化记录三个关键指标总运行成本、总碳排放量、各时段的购电功率曲线。后面的灵敏度分析都基于这套基准结果做对照。4.2 碳交易价格扫描的典型结果与解读我将碳交易价格从 20 元/吨逐步增加到 160 元/吨每 20 元取一个点一共 8 个工况。每个工况都重新完整求解一次 MILP 模型记录成本与排放量。在我的参数体系下结果呈现出一条比较清晰的规律碳价从 20 提高到 60 的阶段总碳排放量下降幅度最明显主要原因是我配置的系统中外购电占比大、本地燃气机组有足够的调节空间碳成本提高后成本对比关系很快翻转系统主动减少了高排放的外购电力。碳价继续从 80 向 160 爬升时碳排放量下降速度逐渐放缓原因是燃气机组的热电比约束限制了进一步降低排放的能力。这个结果指向一个很重要的结论碳价对系统低碳化转型的边际作用不是恒定的。在某个区间内提价效果明显超过拐点后再往上加价环境收益会越来越小但对系统运行成本的冲击会越来越大。这个拐点位置和系统设备配置直接相关不同系统要分别分析。4.3 负荷扰动与新能源渗透率的影响分析负荷扰动灵敏度我做了两种方式第一种是整体按比例增加 5% 和 10%第二种是只调整负荷峰值时段模拟对晚高峰负荷预测偏低的场景。结果显示整体负荷增加时系统主要依靠增加购电量来支撑碳排放基本随负荷等比例上升系统内部的机组组合没有太大变化。但如果是峰值时段负荷上调情况就不同了因为峰值时段的购电价格通常更高优化结果会更倾向提前在低谷时段给储能充电、高峰时段放电以及提高 CHP 在高峰时段的出力。这种扰动分析对运行人员更有参考价值因为它直接告诉你日前计划里哪个时段的预测精度最值得关注。新能源渗透率灵敏度的做法则相对简单把风电、光伏的装机容量按 0.5、1.0、1.5、2.0 倍基准值的比例放大然后重新优化每个工况。由于配电网存在支路容量和节点电压约束渗透率不是越大越好。电动率渗透率增加到一定水平后某些时段光伏大发而负荷很轻会出现节点电压越上限的风险。在我的模型里这个现象通过潮流约束传递到调度结果中表现为弃风弃光增加。所以这个灵敏度的价值在于找到最优渗透率区间而不是一味追求“百分之百新能源”。对于结果呈现方式我用的是“参数-指标”二维表格并配了曲线图。曲线图比表格更适合看拐点表格则负责记录精确数值。发布结果时建议把对应的参数文件和结果文件一并归档方便其他人在你的数据基础上继续扩展。碳价元/t总运行成本相对基准%总碳排放量相对基准%外购电量占比变化趋势2096.8103.5上升4099.1101.2微升60100.0100.0基准80101.997.1下降100103.694.8持续下降120105.293.4趋于平缓140106.592.7基本饱和160107.892.3基本饱和提示上面的百分比只是我演示配置下的一组典型输出用来展示分析方法和结果形态。你换一套负荷数据或者机组参数后具体数值会变但趋势形态通常不会有本质变化。5. 实际跑模型中遇到的坑和排查经验5.1 储能SOC递推方向反了成本低得离谱我在第一版模型里把储能效率项写反了充电时除效率、放电时乘效率结果求解出来的总成本比预期低了一大截。原因是储能每充放一次不仅没有损耗反而“凭空多出”一部分能量优化器自然会让储能反复循环充放来套利。这类问题非常隐蔽因为最终结果仍然是一组合理范围内的数据电压不越限、功率不超限只有总成本明显偏低这个异常信号。排查方式很简单单独把储能从模型里拿出来做一个理想化测试给定一个分时电价曲线检查储能是否只会“低谷充、高峰放”并且单次循环的放电量是否严格小于充电量乘以充放电效率。如果发现系统出现一个完整循环后能量反而增加那一定是递推公式里的效率方向有误。5.2 IEEE33 节点数据编号不一致导致拓扑错位这个问题前面提到过但值得再强调一次。IEEE33 节点的支路表在网络上有多个版本有的从 0 开始编号有的从 1 开始编号。我遇到过节点编号整体错一位的情况表现为优化结果里节点 18 这个位置的负荷始终无法被正常供电求解器只能通过增加其他节点购电来弥补最终成本明显偏高很多。为了避免这种问题我后来每次建模前都会做一个拓扑可视化。用 Matlab 的 plot 函数把 33 个节点按坐标画出来支路用线段连接一眼就能看出网络结构是不是标准的“单馈线辐射状”形状。坐标可以在代码里手写不需要精确地理解地理分布只需要能反映节点间的前后连接关系就行。5.3 二进制变量过多导致求解时间失控我最初给每个时段、每台机组都定义了启停变量同时储能又定义了两个独立的充放电 0-1 变量一个 24 小时模型里的二进制变量数量轻松破百。对于这种规模的 MILPGurobi 通常还能处理但如果模型继续扩大比如把场景缩减法或鲁棒优化引入变量数量会成倍增长求解时间就会急剧上升。在纯日前调度模型里有一个比我实际经验更简单的改善办法如果系统内燃机台数不多且不需要考虑频繁启停带来的额外寿命损耗可以不给机组定义启停变量只限制出力上下限并按最小技术出力处理模型就会退化成线性规划问题求解速度能快一到两个数量级。当然这样做的代价是丢失启停优化这一层自由度到底要不要启停变量取决于你是否需要做机组组合层面的决策。5.4 求解器返回 infeasible 时的系统排查顺序跑这类模型经常会遇到约束不可行。Gurobi 会提示模型无解但不会告诉你哪条约束出了问题。我总结了一套排查顺序能覆盖绝大多数情况第一先检查电量平衡是否闭合。把每个时段的电源出力求和、负荷求和看目标函数里有没有漏掉某一个时段或某一个设备的用能项。第二检查储能终值约束。SOC(25) 等于 SOC(1) 这种约束在极端配置下可能与 SOC 上下限冲突。如果储能容量太小、负荷峰谷差又很大系统可能找不到一个可行轨迹。第三检查机组出力是否超过该时段的最大可用功率。比如 CHP 机组的热电比约束如果热负荷很小而电出力又要求很高就可能突破运行区间。第四用 YALMIP 的 diagnose 命令查看冲突约束分析报告它会给出可能导致不可行的约束索引帮助你快速定位。个人建议不要一上来就去改动目标函数的权重因为无解的本质原因是约束矛盾而不是权重取的不合适。5.5 灵敏度扫描实验的批处理与归档技巧最后一个经验是关于如何组织灵敏度分析过程的。手动改参数、手动记录结果这种方式非常容易出错尤其是当你需要扫 8 个碳价、5 个渗透率、3 个负荷比例的时候工况数量会达到几十组。如果手工修改代码里的参数很可能哪一组漏改了自己都未必记得。我用的是循环嵌套批处理的结构外层循环遍历各灵敏度参数组合内层重新构造参数结构体然后调用同一个优化函数把结果写入以参数命名的 CSV 文件。同时我还会把求解器输出的求解时间、目标值、可行性状态一并保存到记录表里后续写报告时直接用脚本汇总。另一个容易忽略的点是随机初始化。优化问题本身是确定性的但如果风光的预测值本身带有随机性做灵敏度分析时需要固定随机种子否则重复执行同一组参数时结果会发生变化就很难判断两组结果之间的差异到底是参数变化引起的还是随机场景引起的。固定随机种子是这类实验里非常小但是非常重要的细节。这套模型后续能扩展的方向还有很多比如把单目标碳交易模型改成双目标并计算 Pareto 前沿或者把日前调度改成模型预测控制的滚动调度框架再或者引入多场景随机规划来刻画风电出力的不确定性。但无论往哪个方向走经济-碳协调的建模逻辑、灵敏度分析的实验思路、以及对潮流约束的处理方式都是最需要打扎实的地基。我实际跑下来最大的体会是这类项目的瓶颈往往不在求解器而在建模时对每一处约束物理意义的理解多花一点时间把约束写清楚、写对后面所有分析都会省心很多。
