多微电网的电能交互研究这几年热度一直不低,但很多人一开始上手就是做集中式调度,等到算例规模变大、参与主体变多之后才发现问题:不同微网属于不同利益主体,谁都不想把内部运行数据全部交给一个中心;集中式模型传到第三方求解器里也容易遇到计算瓶颈。我自己做这个方向时,把算例从集中式改成交替方向乘子法(ADMM)分布式求解,同时把碳排放交易机制加了进去,前后调了两周参数才稳定收敛。这篇文章就完整梳理一下整个建模思路和Matlab实现过程,给正在做多微网电能交互方向的同学一个可以直接参考的框架。
这套内容解决的是这样一类问题:多个微电网通过公共联络线相连,彼此之间可以买卖电能,每个微网内部有分布式电源、储能和负荷,同时还需要考虑碳排放配额和碳交易成本。目标是在满足各微网功率平衡、储能约束、联络线传输限制等条件的前提下,让整个多微网系统的运行成本最低。ADMM的价值在于把全局优化问题拆成多个子问题,各微网独立求解,只交换交互功率信息,既能保护各主体隐私,又能通过迭代达到与集中式接近的最优解。适合电力系统方向的研究生、做微电网能量管理算法的工程师,以及刚入门分布式优化、想快速把ADMM落到实际算例的读者。
1. 多微网模型与碳排放建模
1.1 为什么需要电能交互
先看一个实际场景。假设园区A里装了大量光伏,中午发电充裕,自己的负荷根本用不完;旁边的园区B是商业楼宇,白天负荷高但屋顶光伏装机很小。如果两个微网完全独立运行,园区A只能把多余的光伏弃掉,园区B则要以较高的价格从大电网买电。一买一弃,从整个系统角度看就是双重浪费。
这时候如果有一条联络线把两个微网连起来,园区A把多余的绿色电力卖给园区B,两边都能受益:A获得了卖电收益,B减少了外购电成本,系统整体的新能源消纳率也上去了。微网之间的交互功率就是在这种背景下引入的。多个微网之间的交互还需要考虑线路容量、交互电价、调度主体之间如何结算等一系列问题,这就构成了典型的优化调度问题。
实际算例中我一般设置3个微网来验证算法。MG1配光伏、储能和柴油机,MG2配风电和储能,MG3配光伏、燃气轮机和较重的基础负荷。这样设计的好处是各微网的特征差异足够大,交互功率的曲线会比较有看头,也方便验证不同电源组合下的调度策略。调度周期取24小时,步长1小时,一共24个时段。
1.2 碳排放如何进入目标函数
碳排放在多微网模型里不是一个简单的罚函数,通常要结合碳交易机制来做。常见的做法是采用配额制:每个微网根据历史排放强度或基准线法获得一定数量的免费碳排放配额,实际排放量超过配额的部分需要到碳市场购买,如果实际排放低于配额,富余部分可以出售获利。
这样碳排放成本就变成:
C_carbon_i = λ_c × (E_actual_i - E_quota_i)
λ_c是碳交易价格,E_actual_i是微网i的实际碳排放量,E_quota_i是分配给微网i的免费配额。当 (E_actual_i - E_quota_i) 为正时是支出成本,为负时是收益,目标函数里直接相加就可以。注意这里不需要分段判断,一个线性表达式就同时覆盖了买入和卖出两种情况,模型实现起来很方便。
实际碳排放量需要分来源计算。柴油发电机和燃气轮机的碳排放按燃料消耗量乘以排放系数折算;从大电网购入的电量也隐含着排放,因为电网侧的电力来自火电、水电、新能源等多种能源的混合,工程上通常用区域电网的平均排放因子估算。向大电网售电的部分,对应的碳排放归属需要在建模时小心处理,一般有两种思路:一种是不计入微网自身排放,因为这部分电力的碳排放已经在电网侧核算过;另一种是只计算净购电量的排放。我在代码里采用的是按净交互功率计算间接排放的方式,具体实现时要区分正负方向,否则容易算错。
1.3 整体优化数学模型
综合起来,多微网系统的目标函数可以写为:
min Σ_i [ C_fuel_i + C_om_i + C_购售电_i + C_carbon_i ]
其中C_fuel_i是微网i内柴油机或燃气轮机的燃料成本,通常表示为输出功率的二次函数;C_om_i是运行维护成本;C_购售电_i是微网与大电网之间的购电成本和售电收益,同时也包括微网间的交易结算成本;C_carbon_i是上面提到的碳排放成本。
约束条件分为几个层次。第一层是各微网内部的功率平衡约束,即各电源出力、储能充放电、交互功率之和等于负荷需求。第二层是储能约束,包括充放电功率上下限、SOC状态转移关系、SOC上下限,以及调度周期首末SOC相等。第三层是交互功率约束,包括微网与公共联络线的交换功率上下限,以及各微网与大电网之间允许的最大购售电功率。第四层是分布式电源出力上下限和爬坡约束,燃气轮机或柴油机通常有爬坡率限制。
这个模型直接交给集中式求解器没问题,但一旦涉及多个微网、多个利益主体,集中式方案的弊端就很明显了。下一节详细说ADMM怎么把这个模型拆开求解。
2. ADMM分布式求解思路
2.1 为什么不用集中式求解
集中式优化在算例小、主体单一的时候确实最简单:把所有变量集中起来,用Yalmip配合Gurobi或Cplex一次性求解就结束了。但放到多微网场景下,集中式至少有三个问题绕不开。
第一个问题是隐私保护。各微网属于不同运营主体,内部负荷数据、储能配置、电源成本参数都被看作是商业信息。集中式优化要求每个微网把完整信息上传给调度中心,这在工程实际中很难有说服力。
第二个问题是计算效率。当微网数量增多、每个微网内部节点变多之后,整个优化问题的维度会快速膨胀,而且不同微网的约束条件耦合在一起,求解器的计算时间增长很快。虽然平时的3微网算例还能接受,但扩展到十几个微网、几百个节点时就吃力了。
第三个问题是鲁棒性。集中式架构要求调度中心必须可靠在线,一旦中心节点出故障,整个系统的优化调度就瘫痪了。分布式架构里单个微网故障不影响其他微网的本地优化,整体韧性更好。
ADMM在处理这类模型时有一个非常合适的特性:它天然支持把带耦合约束的优化问题拆成多个子问题,各子问题并行求解,然后通过交换少量的边界变量实现全局协调。实现难度比完全分布式的次梯度类算法低,收敛效果比很多原始对偶类算法稳定,是工程实践里比较稳妥的选择。
2.2 ADMM的分解原理与迭代流程
ADMM解决的是如下形式的优化问题:
min f(x) + g(z) s.t. Ax + Bz = c
在多微网问题里,关键耦合变量是微网间的交互功率。定义变量x_ij表示微网i向微网j输送的功率,那么必然有x_ij = -x_ji,或者用一种更便于ADMM处理的写法:引入一个全局参考变量z_ij,要求x_ij = z_ij且x_ji = -z_ij。
这样可以构造增广拉格朗日函数。各微网在第k+1轮迭代时,需要求解自己的本地子问题:
x_i^{k+1} = argmin [ f_i(x_i) + Σ_j (ρ/2) × (x_ij - z_ij^k + u_ij^k)^2 ]
其中u_ij是对偶变量,ρ是惩罚因子。求解完各子问题后,各微网把交互功率值上传到协调层,协调层更新全局参考变量:
z_ij^{k+1} = (x_ij^{k+1} - x_ji^{k+1}) / 2
然后更新对偶变量:
u_ij^{k+1} = u_ij^k + x_ij^{k+1} - z_ij^{k+1}
这个过程持续迭代,直到满足收敛条件。可以看到,每个微网的子问题里只包含自身变量和交互功率变量,不涉及其他微网的内部变量,隐私保护和并行计算的要求都满足了。
2.3 收敛判据与参数经验
ADMM的收敛判据通常使用原始残差和对偶残差。原始残差反映的是交互功率与全局参考变量之间的偏差,公式如下:
r_k = sqrt(Σ_ij (x_ij^k - z_ij^k)^2)
对偶残差反映的是全局参考变量的变化幅度:
s_k = ρ × sqrt(Σ_ij (z_ij^k - z_ij^{k-1})^2)
当 r_k < ε_pri 且 s_k < ε_dual 时判定收敛。我用的收敛阈值为1e-4,这个精度足够得到与集中式结果接近的解。
关于惩罚因子ρ的调整,这是整个ADMM调参里最影响成败的部分。ρ太小会导致对偶变量更新缓慢,迭代次数显著增多;ρ太大则会使目标函数里的增广项过强,子问题求解时过度牺牲原目标。实际操作中可以这样把握:先运行一次小规模算例,观察原始残差和对偶残差的下降速度,如果原始残差下降慢而对偶残差下降快,就应该调大ρ;反过来则应调小ρ。对于本文规模的多微网算例,ρ通常在0.1到10这个区间内取值,但要注意目标函数中成本量级的归一化,各成本项的量级差别很大时,ρ也需要做相应缩放。
3. Matlab代码实现全流程
3.1 代码框架与测试算例设计
我用的环境是Matlab R2022b加Yalmip工具箱,求解器是Gurobi。Gurobi对二次约束规划问题的求解速度很快,ADMM子问题每轮迭代都能在很短时间内完成。
代码整体分为以下文件:
| 模块 | 功能 |
|---|---|
| main.m | 主程序,负责参数初始化、循环调用子问题求解、更新对偶变量、判断收敛 |
| subproblem_mg1.m | 微网1的本地优化子问题 |
| subproblem_mg2.m | 微网2的本地优化子问题 |
| subproblem_mg3.m | 微网3的本地优化子问题 |
| update_duals.m | 协调层更新全局参考变量和对偶变量 |
| plot_results.m | 绘制交互功率、SOC、各微网出力和碳排放结果 |
测试算例参数方面,微网互联的线路功率上限设为300kW,微网与大电网的交换功率上限设为500kW。柴油机和燃气轮机的出力上限分别为200kW和500kW。储能容量设为600kWh,最大充放电功率150kW,充放电效率0.95,SOC运行范围0.1到0.9,初始SOC设为0.5。碳配额按照各微网历史负荷的一定比例确定,碳交易价格参考试点碳市场,取80元/吨。
负荷和新能源出力数据我用的是一组典型的夏季日曲线,光伏出力在12点到15点达到峰值,风电出力在凌晨和晚间较大,负荷呈早晚高峰特征。这些数据没有用真实某地数据,而是根据典型曲线人工构造的,但不影响算法验证。
3.2 核心代码实现
先看主循环代码。整个ADMM迭代就在这个循环里完成:
% 初始化 P_ex = {zeros(1, 24), zeros(1, 24), zeros(1, 24)}; z = zeros(3, 24); % 全局参考变量 u = zeros(3, 24); % 对偶变量 rho = 1; tol = 1e-4; max_iter = 200; for k = 1:max_iter % 各微网并行求解本地子问题 P_ex{1} = subproblem_mg1(z, u, rho, param); P_ex{2} = subproblem_mg2(z, u, rho, param); P_ex{3} = subproblem_mg3(z, u, rho, param); % 协调层更新全局参考变量 z_new = zeros(3, 24); z_new(1, :) = (P_ex{1}(1, :) - P_ex{2}(2, :)) / 2; z_new(2, :) = (P_ex{2}(2, :) - P_ex{1}(1, :) + P_ex{2}(1, :) - P_ex{3}(3, :)) / 2; z_new(3, :) = (P_ex{3}(3, :) - P_ex{2}(1, :)) / 2; % 更新对偶变量 u = u + rho * (cell_of_Pex_to_matrix(P_ex) - z_new); % 计算残差 r_pri = norm(cell_of_Pex_to_matrix(P_ex) - z_new, 'fro'); s_dual = rho * norm(z_new - z, 'fro'); z = z_new; if r_pri < tol && s_dual < tol break; end end这段代码里的z_new更新方式值得注意。对于两个微网之间的交互,全局参考变量等于两个微网各自求解出的交互功率的平均值,这是ADMM里最常用也最稳定的更新策略。三个微网互联时,每个微网上报的交互功率需要按互联拓扑关系整理成统一矩阵,代码里用cell数组存储各微网的P_ex,再通过一个转换函数得到统一的矩阵形式。这个映射关系在写代码时最容易出错,建议先画一个简单的拓扑图,把每个微网的交互功率对应到矩阵的哪个位置标清楚再写。
再看子问题的Yalmip建模。以微网1为例:
function P_ex = subproblem_mg1(z, u, rho, param) % 提取参数 P_load = param.load1; P_pv = param.pv1; dt = 1; % 决策变量 P_dg = sdpvar(1, 24); % 柴油机出力 P_ch = sdpvar(1, 24); % 储能充电功率 P_dis = sdpvar(1, 24); % 储能放电功率 SOC = sdpvar(1, 25); % 荷电状态 P_ex_sell = sdpvar(1, 24); % 向外输出功率 P_ex_buy = sdpvar(1, 24); % 从外部买入功率 % 约束 Constraints = []; Constraints = [Constraints, P_dg >= 0, P_dg <= 200]; Constraints = [Constraints, P_ch >= 0, P_ch <= 150]; Constraints = [Constraints, P_dis >= 0, P_dis <= 150]; Constraints = [Constraints, SOC >= 0.1, SOC <= 0.9]; Constraints = [Constraints, SOC(1) == 0.5, SOC(25) == 0.5]; Constraints = [Constraints, SOC(2:25) == SOC(1:24) + 0.95*P_ch*dt - P_dis/0.95*dt]; Constraints = [Constraints, P_ex_sell >= 0, P_ex_sell <= 300]; Constraints = [Constraints, P_ex_buy >= 0, P_ex_buy <= 300]; Constraints = [Constraints, P_dg + P_pv + P_dis - P_ch + P_ex_sell - P_ex_buy == P_load]; P_ex_tot = P_ex_sell - P_ex_buy; % 净交互功率,正为输出 % 目标函数:燃料成本 + 运维成本 + 碳排放成本 + ADMM增广项 Objective = sum(0.3*P_dg + 0.002*P_dg.^2) ... + sum(0.02*(P_ch + P_dis)) ... + 80 * (0.872*sum(P_dg) + 0.581*sum(P_ex_buy) - param.quota1) ... + rho/2 * sum((P_ex_tot - z(1,:) + u(1,:)).^2); ops = sdpsettings('solver', 'gurobi', 'verbose', 0); optimize(Constraints, Objective, ops); P_ex = value(P_ex_tot); end这个子问题里混入了碳排放成本,注意间接排放部分用的是P_ex_buy,也就是从外部购入的功率,而不是净交互功率。这是一个容易算错的地方:如果微网同时存在购入和售出,用净交互功率计算碳排放会低估或高估排放量。但在实际建模中,同一时段既买入又卖出通常不是最优选择,所以很多简化模型直接用净交互功率的正负来判断购售方向,结果差异不大。不过如果想做得严谨,还是建议把P_ex_sell和P_ex_buy分别建模。
3.3 结果展示与后处理
迭代完成后,可以画出各微网的交互功率、SOC曲线和柴油机出力。一个明显能看到的现象是:MG1在中午光伏大发时段向外输出功率,MG2在凌晨风电大发时段向外输出功率,MG3因为负荷重且本地电源有限,大部分时段需要从其他微网购电。这说明交互功率的方向和大小直接反映了各微网的发电特性与负荷特性的互补程度。
与集中式结果对比时,我通常计算两个指标:一是目标函数值的相对偏差,一般在0.5%以内;二是各微网交互功率的最大偏差。如果偏差超过2%,重点检查收敛阈值是否太宽松,或者ρ是否需要调整。实测下来,ADMM迭代大约需要80到120次收敛,设置200次上限足够。
4. 调试中的常见问题与避坑指南
4.1 不收敛与解精度问题
我在调试过程中遇到最多的问题是迭代振荡,典型表现是残差曲线上下波动不平滑,或者交互功率在两个值之间反复横跳。原因有时不在算法本身,而在于子问题里二次项的系数与目标函数其他项量级不匹配。
举一个具体例子。如果燃料成本按元/千瓦时计算,数值在0.3左右,而碳排放成本按吨二氧化碳折算后可能在几十甚至上百的量级,两者差距很大。此时ADMM增广项里的ρ/2乘上功率平方,数值可能远大于实际成本项,导致子问题只关注增广项而忽略原目标。解决办法有两种:一是把成本量纲统一,比如把碳成本折算到元/千瓦时;二是把ρ调到与目标函数量级匹配的水平。我习惯在生产成本单位上统一用元,功率单位用千瓦,时间步长用小时,这样各成本项自然在同一个量级上。
另一个常见问题是储能SOC的循环约束处理不当。如果强制要求SOC(25) == SOC(1),在子问题里这个约束会让储能系统提前预留充电量,可能导致最后几个时段的调度策略偏离集中式解。实际中我的做法是把这个约束保留,但把初始SOC设置为0.5,并允许结束SOC与初始SOC有微小偏差,用不等式约束SOC(25) >= SOC(1) - 0.01 来处理。这样既保证了储能运行的周期性,又给优化留了足够灵活性。
4.2 碳排放参数选择技巧
碳排放模型里最敏感的参数是碳价λ_c和配额E_quota。碳价定得越高,系统就越倾向于减少化石能源出力、增加微网间的电能交互和新能源消纳。在写论文做灵敏度分析时,我通常会把碳价从0逐步提高到150元/吨,观察系统总碳排放量的下降趋势和交互功率的变化幅度。
碳配额的分配方式对结果影响也很大。如果配额设置过松,所有微网都有富余配额,碳排放成本变成收益项,可能出现微网为了卖配额而故意增加基础发电的情况,这就不太合理了。建议配额按基准线法设定,即以历史平均排放强度的80%到90%作为分配依据,这样既给了减排压力,又不至于让约束过于苛刻。
关于外购电的排放因子,不同区域电网的数值差异较大,可以从国家或地区发布的最新电网排放因子数据中取一个当前值。需要注意的是这个数值是逐年更新的,论文写作时应标注数据来源和年份。类似地,柴油和天然气的排放系数也应以权威数据库为准,算例里用近似值可以,但要注明。
4.3 与集中式解的一致性验证
写分布式算法之前,我非常建议先做一步:把同一个多微网模型用集中式方法完整求解一遍,把目标函数值和各时段交互功率存下来。这个集中式解就是后续校验分布式算法性能的基准。
具体操作是:在Yalmip里把三个微网的变量和约束全部合并到一个模型里求解,约束里直接用P_ex_i + P_ex_j == 0的形式表示微网间的交互平衡。这样解出来的结果在数学上是全局最优的。然后运行ADMM代码,每轮迭代后计算当前目标值与集中式最优解的偏差。正常情况下列表会呈现快速下降—缓慢逼近—稳定收敛的趋势。
如果你发现ADMM最终解与集中式解的偏差始终降不下去,不要急着调ρ。先检查全局参考变量的更新公式是否正确,特别是多个微网互联时,每个交互方向的参考变量是否都做了平均处理。一个很小的符号错误可能导致偏差永远存在。我在最初写代码时就犯过这个错:P_ex(2,1)和P_ex(1,2)的方向反了,导致两个微网之间的交互功率始终无法达到一致。
4.4 一个容易踩坑的实现细节
最后特别提一个关于Matlab数组操作的细节。在Yalmip中定义24维变量时,price向量和变量向量相乘要小心维度。例如,碳价80元/吨乘以排放总量时,如果排放量是按24小时累加的总吨数,计算结果是一个标量;而如果写成80乘以每小时排放向量,结果就是一个1×24的向量。两种写法在Yalmip里都能触发求解器建模,但含义完全不同,会导致成本计算错误。
用debug工具检查约束和目标函数的维度是简单有效的方法。如果在optimize调用时求解器报“Model is unbounded”或“Non-convex QP detected”,绝大多数情况不是算法问题,而是变量或约束定义时某个维度出错导致模型语义变了。在怀疑ADMM本身之前,先确认模型是正确的,这一条经验能帮你省下大量排查时间。
其实做这套分布式多微网调度,最大的收获是把算法变成可运行代码的过程中,对模型里每个公式的理解会深很多。纸上推ADMM迭代公式很顺畅,但真到写代码时,变量的方向、约束的耦合、参数的量级都会逼着你重新审视问题。如果你也在调类似的算例,建议先跑通两微网,再扩展到三微网或多微网,每次只增加一个复杂度维度。这样即使出了问题,也知道从哪个环节开始排查。后续如果感兴趣,还可以在这个框架上加入需求响应、多时间尺度协调、通信拓扑变化等扩展场景,ADMM的分解结构对这类扩展天然友好。