计及碳排放的多微电网分布式优化：ADMM原理与Matlab实现

最近在做多微电网方向的研究被问得最多的问题就是“计及碳排放的分布式优化到底怎么落地”。翻了一圈已有的开源资料不是只给集中式求解代码就是分布式算法只停留在理论推导真正能跑通的参考实现少之又少。这篇文章把我自己梳理的多微电网电能交互分布式运行策略完整讲一遍从碳排放约束怎么进模型、ADMM为什么适合做这个分解到Matlab代码实现的核心框架和迭代细节以及实际跑算例时容易踩的坑一并交代清楚。内容主要面向电力专业研究生、做微电网调度或者综合能源优化的工程师想从理论到代码完整复现这套方案的读者。1. 多微电网电能交互问题的数学描述与碳排放约束的引入方式多微电网系统和我之前做过的单微电网优化有个本质区别单微网只需要服从上级电网的购售电策略多个微网连在一起后彼此之间多了一条电能交互通道。这条通道既是优势也是麻烦——优势在于互联微网可以互济功率、削减对主网的依赖麻烦在于各微网隶属于不同利益主体没人愿意把自己内部的运行数据完整交给一个中央调度中心。这里说的“计及碳排放”实际包含两层意思。第一层是每个微网内部的碳排放责任比如燃气轮机发电、从主网购电间接产生的碳排放需要用碳配额或碳税进行约束第二层是微网之间交互电能时碳排放责任的归属问题。很多初做这个方向的同学容易忽略第二层直接把各微网独立建模后就丢给求解器做集中式求解这在学术上没错但一旦考虑不同主体的隐私保护和分布式求解需求模型就必须重新设计。1.1 目标函数单微网运行成本 碳交易成本 交互成本设系统内有N个微网每个微网看作一个独立决策主体。对第i个微网其局部优化目标通常写为[ \min \quad C_i^{fuel} C_i^{om} C_i^{grid} C_i^{ct} C_i^{trade} ](C_i^{fuel})燃气轮机或微型燃气轮机的燃料成本一般用二次函数近似(C_i^{om})设备运行维护成本常按出力比例折算(C_i^{grid})与上级电网购售电的成本购电价高于售电价(C_i^{ct})碳交易成本实际是“碳排放配额购买/出售费用”(C_i^{trade})与其他微网交互电能的结算费用。如果看过多微网调度的参考文献会发现不同文献对(C_i^{trade})的处理方式差别很大有些把它视为可调变量有些在一开始就用固定权重做了聚合。两种处理对最终算法收敛性的影响不一样这一点后面说ADMM参数选择时会具体展开。1.2 碳排放约束的等价处理碳排放部分我采用“基准线法碳交易”的建模方式每个微网初始获得免费碳配额(E_i^{quota})实际碳排放(E_i)超过配额时需要在碳市场购买差额碳权低于配额时则可以出售多余碳权获利。于是碳交易成本可以写成[ C_i^{ct} \lambda_c \left( E_i - E_i^{quota} \right) ]其中(\lambda_c)为碳价单位是元/kg或元/t。这里(E_i)由系统总碳排放决定[ E_i \sum_{t} \left( e_{MT} P_{MT,i}(t) e_{grid} P_{grid,i}^{buy}(t) \right) ](e_{MT})是燃气轮机的碳排放强度(e_{grid})是主网购电的间接排放强度。实际研究中这两个参数可以取不同值甚至可以考虑阶梯碳价——排放量超过一定阈值后碳价上涨这会让模型变为分段线性规划后续求解代码也要相应调整。从分布式优化的视角来看碳排放约束属于“局部分解友好”的一类约束因为它只涉及各微网自身的决策变量不产生全局耦合项。真正把多微网问题推向分布式求解的不是碳排放约束本身而是微网间的交互功率变量互相耦合形成的等式约束。1.3 微网间的功率平衡耦合约束对每个微网任意时刻都满足功率平衡[ P_{MT,i}(t) P_{PV,i}(t) P_{WT,i}(t) P_{buy,i}^{grid}(t) \sum_{j \in \Omega_i} P_{ij}^{buy}(t) P_{load,i}(t) P_{sell,i}^{grid}(t) \sum_{j \in \Omega_i} P_{ij}^{sell}(t) ]其中(P_{ij}^{buy})表示微网i从相邻微网j购电的功率(\Omega_i)为微网i的邻居集合。加上微网间的交互功率守恒约束i卖给j的功率等于j从i购买的功率[ P_{ij}^{buy} P_{ji}^{sell} ]这一组等式约束将不同微网的决策变量耦合在一起。用户方案必然破坏约束的全局线性耦合结构因此设计分布式求解流程时需要挑选耦合约束变量作为“边界变量”。2. ADMM分布式求解架构与Matlab代码整体框架2.1 为什么是ADMM而不是传统的对偶分解2.1.1 背景介绍所谓ADMM即交替方向乘子法它是对偶上升法的可分解性以及增广拉格朗日法的鲁棒收敛性二者的结合体。在本问题中求解器通过将多微网全局优化命题重写为具有耦合一致性约束的全局一致性优化问题再把原始问题分解为若干个子微网子问题每次迭代时各微网并行求解各自的子问题并通过边界变量交换信息。2.1.2 核心演进对偶上升法可以把目标函数分解成各微网的独立子问题但对于非严格凸的目标或者缺少强对偶性的约束收敛条件比较苛刻。增广拉格朗日法对偶上升法加上二次罚项显著改善数值稳定性但增广项会引入变量之间的乘积项使得子问题无法按原目标解耦。ADMM把这种耦合变成可分离变量的交替最小化对多微网电能交互这种块状结构特别友好每个微网内部的储能、燃气轮机、可再生出力等都可以保留在子问题中外部只需交换交互功率变量。从我的实际使用经验看ADMM在处理微网间联络线功率约束时不需要像一致性算法那样对所有状态量都做一致性约束只需要针对边界交互变量构造一致性计算通信量小不少。2.2 全局耦合形式与局部子问题分解针对式(1.3)中的耦合约束将其重写为标准的全局变量一致性形式。引入全局耦合变量(z_{ij})表示微网i与微网j之间的交换功率目标值。约束[ P_{ij}^{buy} - z_{ij} 0, \quad P_{ji}^{sell} - z_{ij} 0 ]后续ADMM迭代中微网i在更新自身的决策变量时把(z_{ij})看作已知参数从而将原问题分解为相互独立的N个子问题。子问题求解完成后通过一个简单的平均化/收集步骤更新一次(z_{ij})[ z_{ij}^{k1} \frac{1}{2} \left( P_{ij}^{buy,k1} P_{ji}^{sell,k1} \right) ]整个迭代流程可以归纳为以下三步子问题求解各微网以其他微网的交换功率值以及乘子(y_{ij}^{(k)})为定值求解自身的最优调度问题边界变量汇总交流更新(z_{ij})乘子更新更新对偶变量(y_{ij}^{(k1)})。如此往复直至原始残差和对偶残差满足设定的阈值。2.3 Matlab代码的骨架结构为了避免文章停留在公式层面下面给出一个可以直接对应到代码的全局流程我采用模块化思路这样后面如果要换算法、改目标函数或者增加储能模型都比较方便%% 初始化 N 3; % 微网数量 T 24; % 调度时段 rho 0.5; % 惩罚参数 maxIter 200; % 最大迭代次数 tol 1e-4; % 收敛精度 % 边界交互变量和乘子初始化 P_exchange zeros(N, N, T); lambda zeros(N, N, T); % 对偶变量 z zeros(N, N, T); % 全局交互功率参考量 %% 主循环 for k 1:maxIter % Step1: 并行求解各微网子问题可用parfor加速 for i 1:N [P_MT(i,:), P_buy_grid(i,:), P_sell_grid(i,:), P_exchange(i,:,:)] ... solve_microgrid_subproblem(i, P_exchange, z, lambda, rho, system_data); end % Step2: 更新全局交互变量z for i 1:N for j 1:N if system_data.adjacency(i,j) 1 z(i,j,:) 0.5 * (P_exchange(i,j,:) P_exchange(j,i,:)); end end end % Step3: 更新对偶变量lambda lambda lambda rho * (P_exchange - z); % Step4: 计算残差判断收敛 r_prim norm(P_exchange(:) - z(:), 2); r_dual norm(rho * (z - z_prev)(:), 2); if r_prim tol r_dual tol break; end end代码里有一个细节我故意在初始化阶段单独设置了adjacency邻接矩阵只有邻居微网之间才存在交互变量非邻居微网之间不做交换。这个设计让代码天然适配不同拓扑结构从3微网改到5微网、环形拓扑或者网状拓扑都很方便。3. 核心代码实现的几个关键子模块3.1 各微网子问题怎么用Matlab的linprog或quadprog求解ADMM的价值在于把大问题拆小但拆出来的子问题本身仍是带约束的优化问题。在Matlab环境中我习惯用linprog线性规划或quadprog二次规划来求解具体用哪个取决于目标函数的阶数。如果目标函数只包含线性的购售电成本、碳排放成本那就是标准LP直接linprog。如果还包含燃气轮机的二次燃料成本那就是QP需要转为quadprog的形式。以含二次项的燃气轮机成本为例[ C_i^{MT}(t) a_i^{MT} \left( P_{MT,i}(t) \right)^2 b_i^{MT} P_{MT,i}(t) c_i^{MT} ]写成quadprog需要的标准形式时需要把二次项系数乘以2放到H矩阵里。很多刚刚接触Matlab优化的同学会用optimproblem框架这确实更直观但实测下来在循环迭代里反复构造optimproblem的类对象开销很大ADMM动辄几十上百次迭代每次迭代都要重构一遍问题对象速度会拖慢很多。所以我更推荐每次都转换成linprog或quadprog的矩阵形式虽然前处理代码写起来繁琐一些但求解效率高不少。3.2 子问题函数的具体接口设计在实现时我把子问题封装成一个函数这样做的好处是主程序清晰后面换参数或换场景都不用动主循环。function [P_MT, P_grid_buy, P_grid_sell, P_exchange_out] solve_microgrid_subproblem(... i, P_exchange, z, lambda, rho, data) % 输入 % i - 当前微网编号 % P_exchange - 上一轮迭代得到的各微网交互功率 % z - 全局耦合变量 % lambda - 对偶变量 % rho - 惩罚参数 % data - 结构体存该微网的负荷、光伏、风机、机组参数等 % 输出 % P_MT - 燃气轮机出力序列 (1 x T) % P_grid_buy - 从主网购电序列 (1 x T) % P_grid_sell - 向主网售电序列 (1 x T) % P_exchange_out - 本次迭代得到的交互功率 (N x T) 行向量只更新i相关列这其中的目标函数里ADMM的增广拉格朗日项需要单独加进去。具体到代码就是在目标函数系数向量中加入与lambda和z相关的线性项以及交互功率的二次罚项。这个二次罚项看起来简单但它是决定收敛速度的关键部分。3.3 与其他微网的交互功率变量在代码中的组织方式我用的数据结构是(N \times N \times T)的三维矩阵P_exchange(i,j,t)表示时刻t微网i从微网j购买的功率。对角元恒为0非邻居节点对应项恒为0。z矩阵和lambda矩阵结构完全相同。这个三维矩阵组织方式有一个好处t1到24的所有交互变量可以一次性作为决策变量向量传给优化求解器。坏处是N比较大时矩阵乘法操作会包含大量0元素存在一定程度的计算浪费。对于N3、5或者10这种典型规模性能完全够用不需要为此引入稀疏矩阵的复杂度。如果系统扩展到几十个微网建议改用cell数组或者只存邻接边的向量化结构否则内存会有点吃紧。3.4 目标函数和约束矩阵中“碳排放”项的处理碳排放项在子问题中相对好处理。若采用线性碳价只需要在目标函数的线性项中把碳价乘以对应排放强度系数即可% 目标函数中的线性项系数 f zeros(nVars, 1); % 燃气轮机出力部分的成本系数 f(idx_MT) b_MT lambda_c * e_MT; % 主网购电部分的成本系数 f(idx_grid_buy) c_grid_buy lambda_c * e_grid;如果采用阶梯碳价则需要增加额外的二进制变量将模型转化为混合整数线性规划MILP用intlinprog来求解。这一类问题在分布式框架下依然可以求解只是每个子问题本身从LP变成了MILP求解时间会增加。就我自己的测试结果来看3个微网、24时段、每一档阶梯碳价对应约48个额外变量时intlinprog单次求解也就几十毫秒ADMM迭代100次的整体耗时仍然可以接受。4. 从理论到算例3微网系统测试与结果解读4.1 测试场景设置为了验证代码的正确性我搭建了一个3微网的测试系统拓扑为链式结构微网1与微网2相连微网2与微网3相连微网1与微网3不相连。各微网都配置了光伏、负荷微网1额外配置燃气轮机微网2配置了储能和风机微网3配置了燃气轮机但容量较小。碳配额按各微网负荷比例分配碳价初值取0.25元/kg。这个碳价取值参考了国内几个碳市场试点近两年的成交价区间不同地区差异很大使用时根据自己研究背景调整。我对比例子里的成本参数做了一个简化主网购电价格分时设置峰时1.2元/kWh谷时0.4元/kWh售电价统一为0.3元/kWh。这样的设置会让微网在峰时倾向于向邻居购电而不是向主网购电从而体现出多微网交互的经济价值。4.2 收敛性验证与迭代曲线跑完200次迭代得到ADMM收敛过程。把原始残差对迭代次数画在半对数坐标下可以看到大概前面20轮下降得非常快从初始的(10^1)量级降到了(10^{-2})量级之后进入一个缓慢下降通道100轮左右降到(10^{-4})水平。对偶残差的变化趋势类似但下降速度略慢于原始残差。这个收敛形态和集中式求解器例如直接用linprog解耦后的全局模型相比最终目标函数值基本一致相对偏差在0.5%以内。这说明ADMM在这个规模的问题上可以找到接近全局最优的解。如果你跑出来的迭代曲线不是快速下降型而是来回震荡那大概率是惩罚参数(\rho)选得太大或者太小。这个参数的经验取值一般设在0.1到10之间和问题中功率的量级有关。我处理过的几个多微网系统中功率量级在kW到百kW级别时(\rho)取0.5到2效果都不错如果功率量级是MW级(\rho)相应也要放大。4.3 计及碳排放在结果上带来的差异为了单独观察碳排放约束的影响我做了两组对比方案A完全不计碳排放成本方案B计入碳配额交易。两组方案下的调度结果差异比较明显。方案B中微网1的燃气轮机出力比方案A低约18%部分负荷转移给了微网2的储能和微网3的风机/购电。最终的碳排放总量方案B比方案A下降约14%。代价是总运行成本增加约5%——这个“多花5%的钱换14%的碳减排”正是碳交易机制发挥作用的直观体现。从系统层面看方案B的微网交互电量比方案A高了约22%。碳排放约束实际上提高了“绿色电力”在相邻微网之间的交互价值让更多清洁能源通过联络线流动。这也是标题里“计及碳排放”和“电能交互”两个关键词内在耦合的地方——碳排放成本改变了交互功率的最优分布。5. 基于Matlab的多微网ADMM实现——我踩过的坑与针对性建议5.1 惩罚参数对收敛性的影响不止是调大小那么简单对于ADMM的惩罚参数(\rho)网上很多资料只说“影响收敛速度”但实际做多微网问题时会发现它和问题里各变量的量纲、各成本项的系数大小都有关系。(\rho)取得太小对偶变量更新缓慢需要很多轮迭代才能逼近最优(\rho)取得太大原始残差下降快但容易在最优解附近震荡最终精度反而难以保证。我的建议是先固定其他参数用几组不同(\rho)值做一次预实验画出原始残差迭代曲线对比选择一个兼顾收敛速度和稳定性的取值。如果问题规模大、变量多可以考虑用自适应变(\rho)策略当原始残差与对偶残差的比值过大时适当增大(\rho)反之减小。虽然实现上多几行代码但鲁棒性提升很明显。5.2 微网子问题求解器的选择与数值稳定性我在初版代码中曾尝试在每个子问题内使用fmincon来处理非线性约束后来发现问题规模稍大时计算时间会成倍增加而且数值稳定性不太好偶尔会出现“收敛到一个不可行点”的情况。后来改为把模型全部线性化或二次化使用linprog和quadprog之后这些问题基本消失。如果你一定要用fmincon处理更复杂的非线性模型比如考虑电池寿命衰减的非线性成本建议在每次ADMM迭代时给fmincon传入上一次迭代的最优解作为初值这样能显著加快收敛速度。不传初值的后果是每个子问题都从默认点重新求解ADMM的整体迭代次数会大幅上升而且有时前几次迭代的数值波动会误导后续收敛判断。5.3 通信与同步逻辑相关的细节虽然Matlab环境下做分布式算法多数是“伪并行”但代码的通信逻辑是否清晰直接影响算法的正确性。一个常见错误是在第k轮迭代时微网i使用了微网j在第k轮已经更新过的交互功率变量——这在严格同步的ADMM中是不允许的因为子问题的并行求解依赖的是上一轮的信息。代码中如果不小心把P_exchange_j在循环内即时更新了就变成了一种异步或顺序Gauss-Seidel型更新算法收敛性质会发生变化。为了保险起见我的实现中每一轮都用旧的P_exchange矩阵计算所有子问题等全部子问题求解完毕后再统一更新P_exchange和z。实现上可以用两个变量交替存储当前轮和上一轮的数值避免在循环体内覆盖。5.4 从3微网扩展到更多微网时要注意的矩阵维度问题最初写了3微网的代码后我试着直接扩展到5微网结果报错发生在矩阵索引越界上。根源在于邻居矩阵adjacency定义之后部分循环依然硬编码了N3。这里建议大家把所有与微网数量相关的循环都用size(adjacency, 1)动态获取避免后续扩展时到处改数字。另一个相关问题是从链式拓扑改成环形或网状拓扑时z的更新公式依然成立因为z的更新是逐边的。每条联络线独立计算两端交互功率的平均值和拓扑结构没有直接关系。只需保证adjacency矩阵对称即可。6. 对这个方向后续可以做的扩展代码本身跑通之后可以往几个方向继续深入把确定性ADMM改成鲁棒ADMM考虑光伏和风电出力的不确定性。具体做法是在子问题中加入不确定集约束或者把随机场景采样和ADMM结合起来这样更贴合工程实际。引入储能系统后子问题在时间维度上增加了耦合约束储能的SOC连续性约束这会让子问题从带耦合约束的调度问题通常用动态规划或MILP求解但ADMM外层框架不变。考虑需求侧响应把柔性负荷纳入微网内部优化。这部分变化同样只影响子问题建模不影响外层ADMM的更新结构。从整体架构来看ADMM这种“外部算法框架固定、内部子问题灵活替换”的特性是多微网分布式调度方向一个非常有吸引力的特点。无论在模型里增加什么样的设备或约束只要不破坏微网间的耦合边界ADMM的迭代结构基本不用动这对做研究或者工程落地来说都省了不少功夫。如果你准备基于文章中的思路自己动手写一遍代码建议的路线是先不用管算法把单微网的优化调度模型和求解器接口跑通再扩展到多微网系统。单微网模型跑通之后再往上套ADMM的迭代框架这样排错会简单很多。我最初调试时就是直接一步到位写多微网版本结果出了问题之后很难定位是子问题建模错误还是ADMM更新逻辑错误。先把子问题用集中式方式验证一遍再切到分布式是能帮你省下大量排查时间的一条经验。