最近复现了一个比较有代表性的项目考虑双重低碳需求响应的电力系统优化调度研究并把整套Matlab代码实现完整跑通了。所谓双重低碳需求响应通俗讲就是让负荷不再是一条“只等着被满足”的曲线而是变成能够在削峰填谷的同时也响应碳减排信号的灵活性资源。这个思路很实际因为传统调度只算煤耗、不算碳账时经济最优不一定低碳只做普通需求响应又可能因为削峰补谷的时机不对把负荷补到煤机出力很低的高碳时段。如果你正准备复现同类论文或者需要拿Matlab做低碳电力调度这篇文章可以当一份调试手册来用。我不会绕弯子从模型假设、目标函数、约束写法、代码结构到算例对比和最容易被坑的参数一次性摊开讲清楚。复现过程会涉及线性化处理和商业求解器配置我也会把实际踩过的坑说明白。1. 项目定位这个调度问题到底在优化什么1.1 双重低碳需求响应不是普通DR换了个名字很多文献里的“需求响应”主要是两类价格型需求响应即用户看到分时电价后调整用电时段激励型需求响应即用户与调度机构签协议在特定时段被调用削减负荷并获得补偿。这两类机制都能改善系统负荷曲线但它们自身并没有主动回答“这样调整是否低碳”的问题。“双重低碳需求响应”要把两个信号同时送到负荷侧一个是电价信号引导用户削峰填谷另一个是碳成本信号引导用户避开高碳发电时段。换句话说第一重需求响应服务系统经济性第二重需求响应服务系统低碳性。我在复现时没有把两个机制简单相加而是让它们在同一个调度模型里耦合可调负荷在峰谷和碳排两个维度上同时参与优化最终让系统在成本可控的前提下降低碳排放。一个很典型的例子晚高峰电价很高但此时风电出力可能很小机组组成里高煤耗、高排放机组占比高。普通需求响应会把晚高峰负荷削减并且允许用户在低谷时段补回一部分电量。问题在于低谷时段如果是高排放火电机组的最小出力区间这部分补回电量反而会推高碳排放。加入碳信号后模型会尽可能把负荷转移到新能源出力高、碳排放强度低的时段而不是只看电价低谷。这套逻辑落到优化代码里本质上是一个机组组合叠加经济调度问题只是多了负荷侧调节变量和碳交易成本项。决策变量一部分是火电机组的启停状态这部分是0-1整数变量另一部分是机组出力、负荷削减量、转移量、风光实际出力等连续变量。因此模型会是一个混合整数线性规划或者经线性化后的MILP。1.2 复现前要先确定哪些基本假设论文复现最怕的是模型边界不清楚。我这次复现采用的是典型日24时段集中调度发电机包含火电、风电、光伏负荷侧有常规负荷和可调节负荷忽略输电网网损和节点电压不讨论线路阻塞。这样可以把研究重点集中在“优化调度”本身的机理上代码也更适合在个人电脑上跑通。风电机光伏出力按预测曲线输入暂时不做随机场景。如果后面想扩展成鲁棒优化或随机优化可以在此基础上增加不确定性集合或场景集但基础版本不建议一上来就上随机规划否则遇到无解时很难判断是模型问题、数据问题还是算法问题。还有一个重要假设与碳排放配额有关我采用“免费配额 实际排放差额交易”的碳成本模型。火电机组根据发电量获得一定免费碳配额超出的排放需要购买配额盈余部分可以卖出。这个机制的好处是碳成本会进入目标函数系统会自动在提高煤电出力与购买碳排放权之间做经济比较而不是拍脑袋设定一个硬性排放上限。1.3 为什么用Matlab做这个复现更合适这个项目需要处理整数变量、多时段约束、数据可视化和参数修改Matlab加YALMIP工具箱是比较顺手的一种组合。YALMIP可以用很接近数学模型的语法写约束避免手工展开一大堆矩阵。对比直接用Matlab自带的intlinprogYALMIP代码可读性高很多。后续更换求解器比如从CPLEX换成Gurobi只需要改一行sdpsettings配置。对于需要反复调整约束和查论文公式的复现场景这个优势非常明显。2. 从头搭模型目标函数与约束逐项拆解2.1 目标函数里究竟放了多少项我在模型里把目标函数定义为系统总成本最小化。这里不追求最小化单个指标而是把发电成本、启停成本、需求响应成本、碳交易成本、弃风弃光惩罚都放进同一个目标函数。第一项是火电燃料成本。机组i在时段t的煤耗成本一般写成出力P的二次函数ai乘Pi平方加bi乘Pi再加ci。因为二次函数会让模型变成MIQP很多复现论文会采用分段线性化把它变成MILP。如果机组数不多二次规划也能用CPLEX直接求解但我要强调一旦后面增加碳价分段函数、DR分段补偿等逻辑线性模型会更稳定。第二项是启停成本。火电机组热态启动和冷态启动成本不同但模型里通常简化为每次启动产生一个固定启动费用停机费用也按固定值处理。由于启停变量是0-1变量这项成本本身是线性的。第三项是需求响应补偿成本。这里分两部分可削减负荷按削减电量乘以补偿单价计算可转移负荷按实际触发的转移量计算少量调节补偿。补偿单价取多少很关键后面在第5章会专门讲参数坑。第四项是碳交易成本。火电机组实际碳排放量与免费配额做差后乘以碳价。差额为正就是购买配额成本差额为负时就是出售配额的收入。这项直接让低碳要求变成经济信号。第五项是弃风弃光惩罚成本也叫新能源消纳惩罚。实际风电光伏上网电量如果低于预测值会产生一个很大的罚系数。这个罚系数不是物理成本而是代表消纳新能源的优先级避免求解器为了省钱而随意弃风弃光只在极特殊情况下才允许削减新能源。2.2 功率平衡与旋转备用的约束写法约束条件是优化模型的血肉。第一类是节点功率平衡忽略网损后就是系统总出力等于总负荷加需求响应后的负荷变化。我习惯写成 火电出力之和加风电出力加光伏出力等于原始负荷减去削减负荷再加上转移负荷净变化。风电和光伏出力上限分别取预测值且不会为负。 风光的预测误差在确定性模型里不处理只通过备用约束间接保留一定调节能力。第二类是机组出力约束。当机组停机时出力必须是零当机组启动后出力要落在该机组的最小技术出力和最大出力之间。这个逻辑需要引入0-1运行状态变量u用u乘以Pmin和Pmax去限制实际出力不能只写一个Pmin小于P小于Pmax。第三类是爬坡约束。机组在当前时段和下一时段的出力变化必须控制在最大上升速率和最大下降速率以内。如果相邻时段有启停机动作爬坡约束需要小心处理不然机组刚启动就要立刻抬起很多出力会超出实际物理能力。第四类是最小开机和最小停机时间约束。火电机组启动成本高不能频繁启停。这个问题一般用一组滚动约束来表达代码里循环每个机组、每个时段检查前面若干时段的运行状态之和。如果增加了求解时间过长的问题可以适当弱化最小停机时间因为它对总成本影响通常不如启停成本和燃料成本大。第五类是系统旋转备用约束。每个时段所有在线机组的最大可增出力要大于当前负荷加需求响应后的预测负荷与安全备用之和。新能源占比越高备用需求越大。2.3 需求响应模型怎么拆成可计算的线性约束需求响应建模我采用了“可削减负荷 可转移负荷”两类变量。可削减负荷表示峰时段用户在中短时间段内减少用电这一部分电量不补回可转移负荷表示某些可延迟生产的工业负荷或电动汽车充电负荷从高峰挪到低谷但全天总电量必须守恒。用YALMIP定义时可以简单写成一个长度24的行向量。可削减变量Pcut的每个元素都应该有下限0和上限PcutMaxPcutMax由负荷基数和用户参与比例决定。可转移负荷则拆成“转入”和“转出”两个非负变量任意一个时段负荷先转出再在另一个时段转入因此全时段转出总量要等于转入总量。需要特别注意的是需求响应后的负荷不能只盯“削峰”还要防止低谷时段填谷过度。现实里调峰资源不足时低谷填谷过度会导致机组下压到最小技术出力以外。所以我在模型里另外加了约束要求可转移负荷造成的净增负荷不能超过低谷时段系统调节裕度的某个百分比这个约束经常是模型从不可行变为可行的关键。对于低碳需求响应我不单独建立一个复杂的用户行为公式而是把碳交易成本放进目标函数后让模型自身去权衡高碳时段削减负荷的边际收益。当某时段电网需要高碳机组出力负荷稍微削减就能显著降低碳排放额度缺口优化器会自动增加该时段的可削减量反过来在清洁能源主导时段削减负荷不划算优化器不会刻意削减。这个相互作用在代码里不增加额外约束但需要观察结果时把各时段削减量和对应碳排放量画在同一张图里验证。2.4 碳交易成本与阶梯碳价的线性化思路碳交易成本的基本表达式是实际排放量减去免费配额再乘以碳价。如果只设固定碳价这是一个线性表达式非常简单。但不少论文会用阶梯碳价来模拟“超排越多、碳价越高”的惩罚逻辑。此时不能直接把碳价当作常数乘进去因为超过配额的部分可能需要分段计价。我在代码里没有使用二进制变量做分段而是用低成本区间优先的特征来简化。设超过配额的总缺口为x把它拆成三段非负变量x1、x2、x3每段分别有容量上限并让分段价格从低到高递增。因为求解器做最小化它会优先填满低价区间然后才用高价区间。这样分段碳价可以用纯连续变量完成不增加整数变量求解速度更快。这种处理方式也保留了碳价机制的可解释性排放缺口较小只购买低价配额排放缺口越大边际碳成本越高系统就会更主动地用低碳机组替代高碳机组。3. Matlab代码实现从框架到关键代码3.1 先搭好代码目录和数据结构复现论文代码不能把所有内容堆在main脚本里。我习惯把工程拆成几个文件 主脚本负责定义输入文件路径、调用模型、输出结果参数文件负责记录所有机组、负荷、风光、碳参数建模函数负责创建决策变量、约束和目标函数求解脚本负责调用求解器并提取结果画图脚本负责绘制负荷曲线、机组出力和成本分布。这样拆最大的好处是定位问题快。如果模型不可行可以先在参数文件里把需求响应与碳成本全部屏蔽先把最基础的经济调度跑通再逐步把模块加回来。如果所有逻辑堆在一个脚本里每次调试都要从头运行几十行初始化代码非常浪费时间。3.2 用YALMIP定义决策变量和约束决策变量定义的核心代码通常是这样 P sdpvar(nUnits, nPeriods); u binvar(nUnits, nPeriods); xStart binvar(nUnits, nPeriods); Pcut sdpvar(1, nPeriods); PtrUp sdpvar(1, nPeriods); PtrDown sdpvar(1, nPeriods);P是火电出力矩阵u是运行状态xStart是启动动作。Pcut、PtrUp和PtrDown是负荷侧调节变量。约束添加用约束数组C保存 Constraints []; for t 1:nPeriods Constraints [Constraints, sum(P(:,t)) Pwind(t) Ppv(t) PbaseLoad(t) - Pcut(t) PtransferNet(t)]; end Constraints [Constraints, 0 Pcut PcutMax]; Constraints [Constraints, sum(PtrUp) sum(PtrDown)];功率平衡约束里的PtransferNet是PtrUp减PtrDown这样某一个时段的净转移量可为正也可为负。全时段转入等于转出保证可转移负荷不会凭空增加总用电量。对于机组运行状态和出力上限的约束我通常会写 Pmin .* u P Pmax .* u; 这一段虽然简单但YALMIP里Pmin和Pmax需要做维度广播。为了避免维度错误我建议把这两个参数写成与P同维的满矩阵而不是单纯的行向量或列向量。3.3 碳成本代码的具体写法碳成本模块最容易被单位坑到。我在代码里统一用MW和h为单位电量就是MWh碳排放因子用tCO2/MWh。若某台机组排放因子是0.46表示它每发1MWh电产生0.46吨二氧化碳。免费配额我按系统设定的配额系数乘以每台机组实际发电量来计算这样配额与发电量成正比配额系数取0.35时相当于0.35tCO2/MWh以内的排放不需要买碳配额。 E_quota quotaCoeff * sum(P, 2); E_emit co2Intensity * sum(P, 2); E_balance sum(E_emit) - sum(E_quota);如果碳价线性可以直接写成 C_co2 carbonPrice * E_balance; 此时E_balance为负代表