(顶刊复现)基于配电网韧性提升的应急移动电源预配置和动态调度(下)—MPS动态调度

参考文献：

[1] Lei S , Chen C , Zhou H ,et al.Routing and Scheduling of Mobile Power Sources for Distribution System Resilience Enhancement[J].IEEE Transactions on Smart Grid, 2019:5650-5662.DOI:10.1109/TSG.2018.2889347.

        这篇博客是上述SCI一区论文的部分复现，即采用求解MPS动态调度策略的部分。

1.摘要

        移动式电源（MPSs）包括电动汽车（EV）车队、车载移动能量存储系统（MESSs）和移动应急发电机（MEGs），具有极大潜力增强配电系统（DS）对极端天气事件的韧性。然而，它们的派遣并没有得到深入的研究。本文通过一个两阶段框架实施MPS的韧性配置和调度。在第一阶段，即事件发生前，MPS被预置在DS中，以实现快速的预恢复，从而增强对关键负载的电力供应的可持续性。DS网络也被积极重新配置为受影响较小或压力较小的状态。构建了一个两阶段的鲁棒优化模型，并通过列-约束生成算法来得出第一阶段的决策。在第二阶段，即事件发生后，MPS在DS中被动态调度，以配合传统的恢复工作，从而增强系统的恢复能力。制定了一个新颖的混合整数规划模型，用于优化MPS的动态调度，解决了MPS派遣和DS运行等不同时间尺度之间的耦合，以及道路和电力网络的耦合等问题。在IEEE 33节点和123节点测试系统上进行的案例研究证明了所提出的方法在提高DS韧性方面在MPS路由和调度方面的有效性。

2.原理介绍

        本文详细阐述了配电网韧性曲线。在图2中，R是系统韧性水平的度量。伴随极端天气事件，配电系统经历以下状态：韧性状态 t0∼te，事件进展 te∼tpe，事件后退化状态 tpe∼tr，恢复状态 tr∼tpr，修复后状态 tpr∼tir，基础设施恢复 tir∼tpir。需要注意的是，事件发生后的系统恢复是一个多阶段的过程。在典型和主要的修复和恢复阶段之前，可以立即进行快速的预恢复，例如连接应急备用发电机。通过预恢复，事件后的状态 tpe∼tr 可以进一步分为三个阶段：事件后的状态 tpe∼tp，预恢复状态 tp∼tpp（进行预恢复时），和预恢复后状态 tpp∼tr（预恢复后的状态）。

        本文提出了一个两阶段框架，用于实施移动式电源（MPSs）的韧性路由和调度。在第一阶段，即事件发生前 t0∼te，MPSs被主动预置到DS的一些选定节点，以增强系统的生存能力。在这里，我们通过 tpe∼tr 期间的系统韧性水平来评估生存能力。预置的MPSs能够快速进行系统的预恢复。事件发生后不久，在 tp 时刻，MPSs可以连接到DS，迅速将系统的韧性水平从 Rpe 提升到 Rpp。DS的网络重新配置也被协同优化，将系统转变为受事件影响较小、压力较小的状态[34]。通过这些措施，在事件后的状态，即系统的生存能力水平在 tpp∼tr 期间从 Rpe 提升到 Rpp，如图2所示。相比之下，如果没有MPSs等，系统的事件后韧性水平必须一直保持在 Rpe 直到 tr 时刻。

        在第二阶段，即事件发生后，MPSs被动态调度，以配合修复和基础设施恢复工作，以增强系统在 tr∼tpir 期间的恢复能力。这里的恢复类似于[9]提出的度量方法，用于评估系统的快速恢复程度。MPSs的交通运输和电力调度、DS的动态网络重新配置和功率调度等都被协同优化。因此，在恢复后的状态，系统的韧性水平从 Rpr 提升到 Rpr，并且系统的负载可以在比 tpir 更早的 tpir 时刻得到完全恢复。

2.1 MPS的动态调度

        在极端天气事件发生后，系统的损坏和停电情况被评估，然后恢复和基础设施恢复过程开始。在这个过程中，MPSs也被动态调度，以配合恢复工作，以减少停电范围和持续时间。MPS动态调度的目标是：

        在式（22）中，第一项是 t=1∼T恢复负载的加权和，第二项是由于MPS行驶而产生的运输成本，第三项是由于MESSs和EVs的电池寿命衰减成本。第二项主要出现的两个原因是：1）如果在某个时间 t < T 所有负载都得到了恢复，那么在 t∼T 的时间内就不再需要运输MPSs。在这种情况下，添加第二项可以限制 t∼T 时间内的MPS行驶；2）可能存在多种MPS调度策略可以实现相同的最优恢复效果。添加第二项是为了确保选择能够以最小的MPS运输成本实现最佳效果的调度策略，因为MPS的冗余行驶会导致额外的成本。至于电池寿命衰减问题，由于增强可靠性或韧性的价值较高且突发事件发生频率较低，许多相关的出版物并未考虑这个问题。具体来说，虽然使用车辆对电网频率进行调节时由于频繁充电和放电而引起电池的磨损成为一个问题，但当使用EV进行恢复时，避免了电池快速衰减的问题，因为停电次数有限。然而，在这里，我们包括这一项是为了避免不必要的电池充电和放电。最大化恢复的负载的目标通常是主导的。

        本小节的其余部分介绍了MPS动态调度问题的约束条件。首先，强制执行以下类似于约束(2)-(3)的约束：

        请注意，车辆路径规划本质上是一个NP困难的组合优化问题，本身就具有挑战性。对于涉及路径规划的问题，比如在[42]中，通常会引入多个额外的二进制变量，用于制定路径-流平衡、起止于车场的限制等，这是基于可能进一步简化为仅包含感兴趣节点和链接的有向或无向图的道路网络。在这项工作中，MPSs的路径规划以一种新颖而简洁的方式建模，而不引入新的二进制变量。也就是说，我们只需要约束（26）来确保MPSs在不同节点之间的运输满足必要的行驶时间。 其他涉及的约束，例如路径距离限制、路径-流平衡和请求满足约束，在这里是隐含地得到满足或反映出来，这是由于MPS路径规划与其他决策和约束之间的相互依赖关系。我会用一个简单的例子来解释约束（26）：如果需要MPS 1 在节点1和节点2之间行驶2个时间段，我们会强制执行以下约束：

        因此，如果 β1,1,t = 1（MPS 1 在时间 t 连接到节点 1），那么 β2,1,t+1 = β2,1,t+2 = 0（由于从节点 1 到节点 2 的必要行驶时间，它在接下来的两个时间段内不能连接到节点 2）。反之亦然。很明显，约束（26）也可以等价地转化为:

        这样，就减少了约束条件的数量。关于MPS电源调度的约束条件如下：

        方程（30）和（31）表示了时间内MESSs和EV车队的SoC变化。MESS的SoC由其充电和放电行为决定，而EV车队的SoC也由其行驶行为决定。约束（32）对MESSs和EV车队的SoC范围施加限制。约束（33）和（34）分别指定它们的充电功率和放电功率限制。约束（35）确保在每个时间段内，MESS或EV车队的充电和放电是互斥状态，如果它没有连接到配电系统，则既不能充电也不能放电。约束（36）通过其容量限制了MEG的实际功率输出，并强制要求未连接到配电系统的MEG的实际功率输出为零。约束（37）类似地限制了MPS的无功功率输出。需要注意的是，MEG的燃料可以在配电系统中进行预分配[32]。便携式或可拖动的燃料箱和优化调度的燃料卡车也可以给MEG加油[4]。因此，我们遵循[22]等文献，并在这项工作中假设MEG有充足的燃料。

        上述约束大致可分为两个相互依赖的组。约束（38）-（43）主要涉及网络重构，通过开关切换来实现。约束（38）规定了在时间 t 时，如果支路损坏且仍未修复，则支路应该处于打开状态。约束（39）表示没有开关的未损坏支路保持其初始状态。类似于MPS预定位问题中的约束（4）-（8），约束（40）-（43）确保了动态重构过程中配电系统的径向拓扑结构。具体来说，所有的 \( l_{i,t} \) 被设为 1，约束（40）表明关闭的支路数量应该小于节点数量减去在时间 t 由受损和未修复支路引起的岛屿数量。类似于MPS预定位问题中的约束（11）-（21），约束（44）-（55）主要涉及电力分配。具体来说，方程（46）和（47）通过分别将连接到候选节点的MPS的实际和无功功率输出相加，计算了候选节点的实际和无功功率输出。约束（50）强制执行逐渐增加恢复以达到100%恢复的每个负载。也就是说，它防止对已恢复的负载停电。上述约束大致可以分为两个相互依赖的组。约束（38）-（43）主要涉及网络重配置，通过开关与支路的连接状态变化来实现。约束（38）限制了如果在时间 t 时支路损坏且尚未修复，则支路应保持开路状态。约束（39）表示没有开关的未损坏支路保持其初始状态。约束（40）-（43），类似于MPS预定位问题中的约束（4）-（8），确保DS在动态重配置过程中的径向拓扑结构。具体而言，所有的li,t都设置为1，约束（40）规定了在时间 t 时关闭的支路数量应小于节点数减去因损坏而未修复支路引起的岛屿数量。约束（44）-（55），类似于MPS预定位问题中的约束（11）-（21），主要涉及功率分配。具体而言，方程（46）和（47）通过将连接到候选节点的MPS的实际和无功功率输出相加来计算候选节点的实际和无功功率输出。约束（50）强制每个负载都按非递减的方式恢复以达到100%的恢复。也就是说，它防止了恢复的负载再次掉电。

2.2求解算法

        方程（46）-（47）包含双线性项。这些非线性项可以被线性化。例如，我们可以用下面的约束来替换方程（46）中的 βim,t×cpm,t 项：

        在进行了这样的线性化之后，通过算法1解决MPS预定位问题和MPS动态调度问题时，涉及的MP和SP都变成了混合整数二阶锥规划（MISOCP）模型。如果约束（18）和（52）也线性化了，那么遇到的优化问题就变成了混合整数线性规划（MILP）模型。MISOCP和MILP问题都可以通过许多现成的求解器，如Gurobi，进行高效求解。

3.编程思路

3.1参数和变量定义

表1 相关参数

表2 决策变量

3.2编程思路

        根据对文献内容的解读，可以设计下面的编程思路：

步骤1：输入所需数据

        大部分所需数据文中都已给出，部分没有提供的数据可以自己假设一下。需要注意的是，原文献中并没有给出参数(t时刻的虚拟电源节点集合)和It（在t时刻形成的孤岛数）的取值，但这两个参数并不能随意设定，需要根据不同时刻配电网的故障情况进行设定。为进一步说明，首先根据原文献中表3所提供的故障修复方案，画出初始时刻配电网的故障情况：

        由上图可知，初始时刻共有8条支路发生故障，形成了9个孤立节点，其中孤立节点1形成一个孤岛，孤立节点25，29形成一组孤岛，孤立节点30-32形成一组孤岛，孤立节点17，18，32形成一组孤岛，其余正常节点形成一组孤岛，因此故障发生的初始时刻配电网中共有五个孤岛，可以从每个孤岛中任意选择一个节点作为虚拟电源节点，例如可选择节点1，2，25，30，33作为虚拟电源节点。对于其余时刻的故障情况分析思路完全一致，此处不再赘述。

步骤2：定义第二阶段MPS动态调度的决策变量

        这一步比较简单，按照表2初始化决策变量即可，同时每个决策变量的维度以及类型(sdpvar还是binvar)不要出错。需要注意的是，EV车队、MESS和MEG的数学模型和可连接节点都是不同的，为了方便编程，可以用不同的变量进行定义。

步骤3：写第二阶段MPS动态调度的目标函数和约束条件

        这一步需要按照给定的数据和定义的变量，分别写出第二阶段MPS动态调度问题的目标函数和约束条件。需要重点关注的是约束条件46-47，这几个约束均为非线性约束，包含0-1变量和连续变量相乘的非线性项，可以通过引入一个中间变量进行线性化。Gurobi求解器中也可以直接对这种形式的非线性变量进行处理。

        原文献中的约束条件比较复杂，写约束条件时需要比较细心，最好每新增一组约束，就检查一下优化问题是否可行，确定无误后再列下下一组约束，避免一次性对大量约束进行检查。

步骤4：使用yalmip工具箱对优化问题进行求解

        文中构建的MPS动态调度优化问题是混合整数二阶锥规划形式，可采用Yalmip工具箱进行编程，并调用gurobi求解器进行求解。

4.完整Matlab代码

%% 配电网MPS动态调度优化问题求解

%% 清除内存空间
clc
clear
close all
warning off

%% 系统参数
Copyright;                                  % 配电网参数，脆弱支路和非脆弱支路，负荷曲线
SB = mpc.baseMVA;                           % 基准功率，MVA
VB = mpc.bus(1,10);                         % 基准电压,kV
Yb = mpc.bus(:,1);                          % 节点集合
Nb = length(Yb);                            % 节点数目
Ys = mpc.gen(:,1);                          % 变电站集合
Ns = length(Ys);                            % 变电站数目
Ym = [15,21,29;5,8,33];                     % 第m个MPS可以连接的节点集合
Nm = [3;3];                                 % 第m个MPS可以连接的节点数目
NL = length(mpc.branch(:,1));               % 支路数目
YL = 1:NL;                                  % 支路集合
wi = 5*ones(Nb,1);                          % 所有节点权重均设为5，可调整
Pl_it = 1;                                  % 虚拟负荷
pload0 = mpc.bus(:,3)/SB;                   % 节点 i负荷有功功率需求
qload0 = mpc.bus(:,4)/SB;                   % 节点 i负荷无功功率需求
Vmax = 1.06;                                % 节点电压最大值
Vmin = 0.94;                                % 节点电压最小值
Capi = 1;                                   % 每个节点最大允许接入的MPS
rij = mpc.branch(:,3);                      % 支路电阻
xij = mpc.branch(:,4);                      % 支路电抗
sij_max = 0.5;                              % 支路功率最大值
Pm_max = [300,500,800]/(SB*1e3);            % 第m个MPS的输出有功功率上限
Qm_max = [500,776,600]/(SB*1e3);            % 第m个MPS的输出无功功率上限
G1 = 5;                                     % 脆弱支路的最大故障数
G2 = 3;                                     % 非脆弱支路的最大故障数
K1 = Nb - Ns;                               % 足够大的正数K1
K2 = Nb;                                    % 足够大的正数K2
YT = 0.5:0.5:24;                            % 控制时段集合
NT = 48;                                    % 控制时段数
dt = 0.5;                                   % 一个控制时段的持续时间
tr_ij = ...                                 % MPS从节点i到节点j的通行时间。
    {[0,4,2;4,0,2;2,2,0],[0,2,2;2,0,4;2,4,0]};
SOC0 = 0.5*[300;776]/(SB*1e3);              % 初始荷电状态
nc = 0.95;                                  % 充电效率
nd = 0.95;                                  % 放电效率
SOC_max = 0.9*[300;776]/(SB*1e3);           % 荷电状态上限
SOC_min = 0.2*[300;776]/(SB*1e3);           % 荷电状态下限
tp = 0.25/(SB*1e3);                         % EV车队能耗率
YL_RCS = 1:37;                              % 拥有自动开关的支路集合
dm = 0.002;                                 % MPS的运输成本系数
Zm = 0.000015;                              % MPS电池老化速度
CBm = 0.01;                                 % MPS单位充放电价格
Lij0 = ones(NL,1);                          % 初始支路开断状态
Lij0([2,10,13,15,24]) = 0;


branch_to_node = zeros(Nb,NL);              % 节点的上游支路
branch_from_node = zeros(Nb,NL);            % 节点的下游支路
for k = 1:NL
    branch_to_node(mpc.branch(k,2), k) = 1;
    branch_from_node(mpc.branch(k,1), k) = 1;
end

Ysrc = cell(1,NT);                          % t时刻的虚拟电源节点集合
L_off = cell(1,NT);                         % t时刻仍未修复的支路
It = zeros(1,NT);                           % 在t时刻形成的孤岛数
for t = 1:NT
    if t >= 1 && t <= 6
        L_off{t} = [1,19,33,16,32,29,28,24];
        Ysrc{t} = [1,2,25,30,33];
        It(t) = 5;
    elseif t >= 7 && t <= 12
        L_off{t} = [19,33,16,32,29,28,24];
        Ysrc{t} = [1,25,30,33];
        It(t) = 4;
    elseif t >= 13 && t <= 18
        L_off{t} = [33,16,32,29,28,24];
        Ysrc{t} = [1,25,30,33];
        It(t) = 4;
    elseif t >= 19 && t <= 26
        L_off{t} = [16,32,29,28,24];
        Ysrc{t} = [1,25,30,33];
        It(t) = 4;
    elseif t >= 27 && t <= 32
        L_off{t} = [32,29,28,24];
        Ysrc{t} = [1,25,30];
        It(t) = 3;
    elseif t >= 33 && t <= 40
        L_off{t} = [29,28,24];
        Ysrc{t} = [1,25];
        It(t) = 2;
    elseif t >= 41 && t <= 44
        L_off{t} = [28,24];
        Ysrc{t} = 1;
        It(t) = 1;
    elseif t >= 45 && t <= 48
        L_off{t} = [24];
        Ysrc{t} = 1;
        It(t) = 1;
    end
end

%% 设决策变量
bimt_EV = binvar(Nm(2),NT);                 % 0-1变量，如果EV预配置到节点i，取值为1，否则取值为0。
bimt_MESS = binvar(Nm(1),NT);               % 0-1变量，如果MESS预配置到节点i，取值为1，否则取值为0。
bimt_MEG = binvar(Nm(1),NT);                % 0-1变量，如果MEG预配置到节点i，取值为1，否则取值为0。
ymt_EV = binvar(1,NT);                      % 0-1变量，如果EV在t时刻正在行驶，则取值为1，否则取值为0。
ymt_MESS = binvar(1,NT);                    % 0-1变量，如果MESS在t时刻正在行驶，则取值为1，否则取值为0。
ymt_MEG = binvar(1,NT);                     % 0-1变量，如果MEG在t时刻正在行驶，则取值为1，否则取值为0。
cmt_EV = binvar(1,NT);                      % 0-1变量，如果EV在t时刻正在充电，则取值为1，否则取值为0。
cmt_MESS = binvar(1,NT);                    % 0-1变量，如果MESS在t时刻正在充电，则取值为1，否则取值为0。
dmt_EV = binvar(1,NT);                      % 0-1变量，如果EV在t时刻正在放电，则取值为1，否则取值为0。
dmt_MESS = binvar(1,NT);                    % 0-1变量，如果MESS在t时刻正在放电，则取值为1，否则取值为0。
cpmt_EV = sdpvar(1,NT);                     % EV在t时刻充电功率。
cpmt_MESS = sdpvar(1,NT);                   % MESS在t时刻充电功率。
dpmt_EV = sdpvar(1,NT);                     % EV在t时刻放电功率。
dpmt_MESS = sdpvar(1,NT);                   % MESS在t时刻放电功率。
SOC_mt_EV = sdpvar(1,NT);                   % EV在t时刻的荷电状态
SOC_mt_MESS = sdpvar(1,NT);                 % MESS在t时刻的荷电状态
gp_mt = sdpvar(3,NT);                       % MPS在t时刻的有功输出
gq_mt = sdpvar(3,NT);                       % MPS在t时刻的无功输出
Lijt = binvar(NL,NT);                       % 0-1变量，取值为1时表示支路闭合，取值为0时表示支路断开
git = cell(1,NT);                           % 虚拟电源的出力
for t = 1:NT
    git{t} = sdpvar(It(t),1);
end                      
fijt = sdpvar(NL,NT);                       % 支路ij的虚拟功率
pload = sdpvar(Nb,NT);                      % 节点i的实际有功负荷
qload = sdpvar(Nb,NT);                      % 节点i的实际无功负荷
Pi = sdpvar(Nb,NT);                         % 节点i注入有功功率
Qi = sdpvar(Nb,NT);                         % 节点i注入无功功率
pf_ijt = sdpvar(NL,NT);                     % 支路ij的有功功率
qf_ijt = sdpvar(NL,NT);                     % 支路ij的无功功率
vit = sdpvar(Nb,NT);                        % 节点i的电压平方

%% 约束条件
省略。。。。。。。
%% 目标函数

objective = sum(wi.*sum(pload, 2)) - dm*sum(ymt_EV + ymt_MESS + ymt_MEG) - Zm*CBm*sum(cpmt_EV + dpmt_EV + cpmt_MESS + dpmt_MESS);



%% 设求解器
% gurobi求解器
ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1, 'debug', 1);
ops.gurobi.TimeLimit = 600;                 % 运行时间限制为10min
ops.gurobi.MIPGap = 0.01;                   % 收敛精度限制为0.01

sol = optimize(Constraints, -objective, ops);


%% 分析错误标志
if sol.problem == 0
    disp('求解成功');
else
    disp('运行出错');
    yalmiperror(sol.problem)
end

%% 展示运行结果
show_result;

        以上仅为主函数的部分代码，完整代码获取链接如下：

(顶刊复现)基于配电网韧性提升的应急移动电源预配置和动态调度(下)-MPS动态调度资源-CSDN文库

        全文完整复现代码获取链接如下：

(顶刊复现)基于配电网韧性提升的应急移动电源预配置和动态调度matlab代码资源-CSDN文库

5.代码运行结果

        原文中数据提供不全，所以代码复现结果和原文献相比会有偏差，但原理完全一样。