不确定环境下多周期动力电池回收网络优化

Optimization of Multi-Cycle Power Battery Recycling Network under Uncertain Environment

摘要: “双碳”背景下，动力电池的产销量随着新能源汽车的发展呈快速增长的趋势。由于动力电池退役时仍然有能量剩余，所以梯次利用是针对废旧电池回收再利用的有效措施。而对考虑梯次利用的动力电池逆向物流回收网络进行规划，是提高动力电池回收效率，缓解资源和环境压力的有效方法。本文考虑动力电池回收量与梯次利用市场需求量的不确定性，建立总成本最小的多周期动力电池回收网络的鲁棒优化模型，并构造了测试算例对模型进行求解，并对鲁棒控制参数进行灵敏度分析。

Abstract: Under the background of “dual carbon”, the amount of power battery is showing a rapid growth trend with the development of new energy vehicles. Since there is still energy left when power battery is retired, cascade utilization is an effective measure for the recycling and reuse of waste batteries. Planning a reverse logistics recycling network for power battery that takes cascade utilization into consideration is an effective way to improve the efficiency of power battery recycling and alleviate resource and environmental pressures. This paper considers the uncertainty of the amount of power battery recycling and the market demand for cascade utilization, establishes a robust optimization model for a multi-period power battery recycling network with the minimum total cost, constructs a test example to solve the model and performs sensitivity analysis on the robust control parameters.

1. 引言

在“双碳”背景下，新能源汽车因其具备节能减排的特性而逐渐成为我国汽车行业的主流趋势。根据国际能源署的数据统计，我国自2012年起，电动汽车销量与保有量都呈快速增长的趋势，截止2022年底我国电动汽车销量已达到590万辆，预计2025年销量将达到1460万辆，年销量整体保持50%左右的增速。动力电池作为电动汽车的核心部件，其产销量也呈快速增长的趋势。然而当动力电池的容量衰减至不足80%时，由于无法满足汽车的动力要求将面临退役，而退役动力电池中含有可利用的剩余电量和有价金属，若回收处理不当，会造成环境与资源问题。动力电池的再生利用与梯次利用可以有效地缓解动力电池回收带来的资源浪费与环境污染问题，同时将其剩余电量降级利用在其他领域。所以研究考虑梯次利用的动力电池逆向物流网络优化问题具有实际意义，同时由于动力电池回收起步较晚，回收量与梯次利用需求量存在不确定性，并且动力电池退役已经形成规模化退役的常态，其回收工作是长期动态的决策，所以考虑不确定性与多周期性是有必要的。基于此在不确定环境下研究多周期动力电池回收网络规划问题具有重要意义。

2. 文献综述

目前针对动力电池回收主要聚焦在回收模式选择[1] [2]与定价策略[3]-[5]等领域，针对动力电池回收网络优化的研究较少[6]-[8]，并且其中较少考虑多周期与不确定性。其中Guan等[9]基于第三方回收模式，构建由消费者、第三方回收中心与处理中心组成的逆向物流网络，以此建立考虑居民与环境影响的双目标模型，采用遗传算法求解并对社会负效应进行灵敏度分析。焦建玲等[10]构建了以动力电池与整车制造企业和第三方回收拆解企业为回收主体的Stackelberg博弈模型，求解了最优决策与最优利润并分析再生材料收益对回收模式选择以及闭环供应链系统的碳排放的影响。杨康康等[11]同时考虑了再生利用与梯次利用，研究三种政府干预机制下的回收激励政策，对比无政府干预和政府补贴与奖惩机制下，回收网络各主体的最优决策与废旧电池回收率。杨玉香[12]等为废旧动力电池回收设计了一个考虑经济、环境和社会平衡发展的可持续逆向物流网络，为此建立了动力电池双渠道回收模式的多目标模型，并采用改进混合协同进化遗传算法求解。刘娟娟[13]等考虑动力电池回收量与处理水平的不确定性，采用三角模糊数表示处理成本、运输成本、单位价格、回收量与再循环利用率，采用模糊机会约束建模并对以上参数灵敏度分析。Hao和Sun等[14]以电动汽车电池回收为研究背景，以最小化安全和环境风险，最大化社会责任和经济效益为目标，建立了一个逆向物流网络选址模型。孙强[15]等将逆向物流网络中不同周期的回收数量、再生产成本与市场需求作为不确定参数，建立了多周期多产品的两阶段逆向物流鲁棒优化模型。

综上所述，大多数针对动力电池回收网络优化的研究没有考虑决策的动态性与不确定性，往往建立的是单周期与确定性模型，并且对于不确定的处理较少采用鲁棒优化的方法。基于此，本文构建了不确定环境下多周期动力电池回收网络的鲁棒优化模型，并采用Budget鲁棒优化与对偶理论将模型转化为易于求解的鲁棒对等模型。

3. 问题描述与数学模型

3.1. 问题描述

本文建立了如图1所示的废旧动力电池逆向物流回收网络，该回收物流主要由新能源汽车客户、收运点、分类检测中心、梯次利用中心、回收拆解中心与梯次利用客户组成。

该回收网络的物流环节如下：新能源汽车客户点会产生不满足汽车动力需求的退役动力电池；收运点即为新能源汽车零售商，负责一定范围内的废旧动力电池回收，并将回收到的废旧电池集中运输到分类检测中心；分类检测中心一般为电池制造商，负责对电池剩余容量进行检测分类，满足梯次利用要求的运输至梯次利用中心进行进一步的重组加工，不满足要求的则送往回收拆解中心；回收拆解中心一般由电池回收或电池制造企业负责，回收其中的金属材料进行循环再利用；梯次利用中心主要是专注于废旧动力电池回收的企业，将废旧电池重组为梯次产品根据梯次利用客户的需求进行分配与运输。

Figure 1. Power battery reverse logistics recycling network

图1. 动力电池逆向物流回收网络

该回收网络面临的决策问题是：在每周期产生的退役电池数量与梯次利用市场的需求量不确定的情况下，对收运点和分类检测中心的位置与数量进行决策，以及决策网络中各设施之间每周期的运输量，使得回收网络总成本最小。

根据动力电池回收网络的特点，为简化模型，作出如下假设：

(1) 网络节点之间的单位运输成本与距离成正比；

(2) 不考虑动力电池类型与动力电池的重组损耗；

(3) 不考虑选址的多周期决策，只考虑每个周期流量分配的变化；

(4) 为统一梯次利用的单位回收收益，本文的梯次利用市场仅面向某一单一领域。

3.2. 符号说明

本文所使用的数学符号如表1，其中参数中带有“~”的表示不确定参数。

Table 1. Explanation of symbols

表1. 符号说明

类别	符号	定义
集合	$I$	新能源汽车客户
	$J,K$	设施选址备选点

续表

	$L,N$	回收拆解与梯次利用中心
	P	梯次利用客户
	Z	回收网络设施集合
	T	运营周期
	$U_q,U_p$	参数${\tilde{q}}_i^t,{\tilde{q}}_p^t$ 的不确定集合
模型参数	$f_j,f_k$	开放设施的成本
	$i,j,i,j\in Z$	动力电池回收网络节点
	${\theta }_{ij},d_{ij},i,j\in Z$	回收网络节点间的距离与单位运输成本
	$C_i,i\in Z\backslash \left(I\cup P\right)$	设施的最大容量
	$h_i,,i\in K\cup N$	设施的单位库存成本
	$b_i,i\in N$	梯次利用的单位处理成本
	${\mu }_p,p\in P$	梯次市场的单位缺货成本
	M	一个足够大的常数
	$\alpha$	梯次利用率
	${\tilde{q}}_i^t,i\in I$	废旧动力电池回收量
	$r_m$	收运点的最大服务半径
	$P_r^1,P_r^2$	单位梯次利用与回收拆解收入
	${\tilde{q}}_p^t,p\in P$	梯次利用市场的需求
决策变量	$X_j$	如果收运点j选址，则为1，否则为0
	$X_k$	如果分类检测中心k选址，则为1，否则为0
	$Y_{ij}^t$	t周期将新能源客户i分配给收运点j，则为1，否则为0
	$Q_{ij}^t,i,j\in Z/I$	t周期网络节点间的废旧动力电池运输量
	$Q_{k1}^t$	t周期由分类检测中心k分拣可以梯次利用的电池数量
	$Q_{k2}^t$	t周期由分类检测中心k分拣只能回收拆解的电池数量
	$I_i^t,i\in K\cup N$	t周期梯次利用中心与分类检测中心的梯次产品库存量
	$S_p^t$	t周期梯次利用市场的缺货量
	$W_{ij}^t,i,j\in N,P$	t周期提供给梯次市场的动力电池量

3.3. 数学模型

不确定环境下的多周期动力电池回收网络的鲁棒优化模型如下：

$\begin{array}{l}\min \underset{U_q,U_p}{\max }{\displaystyle \sum }_{j\in J}{f}_j{X}_j+{\displaystyle \sum }_{k\in K}{f}_k{X}_k+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{i\in I}{\displaystyle \sum }_{j\in J}{\theta }_{ij}{d}_{ij}{\tilde{q}}_i^t{Y}_{ij}^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{j\in J}{\displaystyle \sum }_{k\in K}{\theta }_{jk}{d}_{jk}{Q}_{jk}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{\theta }_{kl}{d}_{kl}{Q}_{kl}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{\theta }_{kn}{d}_{kn}{Q}_{kn}^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{\theta }_n{d}_n{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{b}_n{Q}_{kn}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{h}_k{I}_k^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{h}_n{I}_n^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{p\in P}{\mu }_p{S}_p^t-\left({\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{P}_r^1{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{P}_r^2{Q}_{kl}^t\right)\end{array}$ (1)

${\displaystyle \sum }_{j\in J}{Y}_{ij}^t=1;\forall i\in I,t\in T$ (2)

$Y_{ij}^t\le X_j;\forall i\in I,\forall j\in J,t\in T$ (3)

${\displaystyle \sum }_{i\in I}{\tilde{q}}_i^t{Y}_{ij}^t\ge {\displaystyle \sum }_{k\in K}{Q}_{jk}^t;\forall j\in J,t\in T$ (4)

$Q_{jk}^t\le M{X}_k;\forall j\in J,\forall k\in K,t\in T$ (5)

${\displaystyle \sum }_{j\in J}{Q}_{jk}^t=Q_{k1}^t+Q_{k2}^t;\forall k\in K,t\in T$ (6)

$Q_{k1}^t=\alpha {\displaystyle \sum }_{j\in J}{Q}_{jk}^t;\forall k\in K,t\in T$ (7)

$I_k^t=Q_{k1}^t+I_k^{t-1}-{\displaystyle \sum }_{n\in N}{Q}_{kn}^t;\forall k\in K,t\in T$ (8)

$Q_{k2}^t={\displaystyle \sum }_{l\in L}{Q}_{kl}^t;\forall k\in K,t\in T$ (9)

$I_n^t={\displaystyle \sum }_{k\in K}{Q}_{kn}^t+I_n^{t-1}-{\displaystyle \sum }_{p\in P}{W}_{np}^t;\forall n\in N,t\in T$ (10)

${\displaystyle \sum }_{i\in I}{\tilde{q}}_i^t{Y}_{ij}^t\le C_j{X}_j;\forall j\in J,t\in T$ (11)

$I_k^t\le X_k{C}_k;\forall k\in K,t\in T$ (12)

${\displaystyle \sum }_{k\in K}{Q}_{kl}^t\le C_l;\forall n\in N,t\in T$ (13)

$I_n^t\le C_n;\forall n\in N,t\in T$ (14)

${\displaystyle \sum }_{n\in N}{W}_{np}^t+S_p^t\ge {\tilde{q}}_p^t;\forall t\in T,p\in P$ (15)

$d_{ij}{Y}_{ij}^t\le r_m;\forall i\in I,\forall j\in J,t\in T$ (16)

$X_j,X_k,Y_{ij}^t\in \left\{0,1\right\},\forall i\in I,\forall j\in J,\forall k\in K,t\in T$ (17)

$Q_{jk}^t,Q_{k1}^t,Q_{k2}^t,Q_{kl}^t,Q_{kn}^t\ge 0;\forall j\in J,\forall k\in K,\forall m\in M,\forall n\in N,\forall t\in T$ (18)

$I_k^t,I_n^t,S_p^t,W_{np}^t\ge 0;\forall k\in K,\forall n\in N,\forall t\in T$ (19)

其中，式(1)表示由固定成本、运输成本、处理成本、库存成本、缺货成本以及回收收益构成的网络总成本最小。约束(2)表示新能源汽车客户只能由一个收运点负责回收。约束(3)与(5)是逻辑约束，表示只有该候选点被选址了才有对应的服务关系。约束(4)是收运点的流量平衡约束。约束(6)与(7)表示在分类检测中心将废旧电池按照是否可以梯次利用来分类，满足容量要求的送往梯次利用中心进行重组利用，不满足要求则送往回收拆解中心。约束(8)与(10)是分类检测中心与梯次利用中心的库存平衡约束。约束(9)表示当前周期不满足梯次利用条件的废旧电池全部进行拆解回收。约束(11)、(12)、(13)与(14)表示设施的容量约束。约束(15)表示每个周期的梯次利用客户产生的市场需求被满足，未满足部分按照缺货来处理。约束(16)表示收运点的服务能力约束。约束(17)、(18)与(19)是决策变量的取值约束。

4. 鲁棒优化处理

4.1. 目标函数的不确定性处理

对于目标函数，通过引入变量$\tau$ 来将不确定性放入约束条件中处理：

$\begin{array}{l}\min \underset{U_q,U_p}{\max }{\displaystyle \sum }_{j\in J}{f}_j{X}_j+{\displaystyle \sum }_{k\in K}{f}_k{X}_k+\tau \\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{j\in J}{\displaystyle \sum }_{k\in K}{\theta }_{jk}{d}_{jk}{Q}_{jk}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{\theta }_{kl}{d}_{kl}{Q}_{kl}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{\theta }_{kn}{d}_{kn}{Q}_{kn}^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{\theta }_n{d}_n{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{b}_n{Q}_{kn}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{h}_k{I}_k^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{h}_n{I}_n^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{p\in P}{\mu }_p{S}_p^t-\left({\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{P}_r^1{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{P}_r^2{Q}_{kl}^t\right)\end{array}$

$\tau \ge {\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{i\in I}{\displaystyle \sum }_{j\in J}{\theta }_{ij}{d}_{ij}{\tilde{q}}_i^t{Y}_{ij}^t$

所以基于上述处理，上一节的鲁棒模型可以表示为如下形式：

$\begin{array}{l}\min \underset{U_q,U_p}{\max }{\displaystyle \sum }_{j\in J}{f}_j{X}_j+{\displaystyle \sum }_{k\in K}{f}_k{X}_k+\tau \\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{j\in J}{\displaystyle \sum }_{k\in K}{\theta }_{jk}{d}_{jk}{Q}_{jk}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{\theta }_{kl}{d}_{kl}{Q}_{kl}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{\theta }_{kn}{d}_{kn}{Q}_{kn}^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{\theta }_n{d}_n{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{b}_n{Q}_{kn}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{h}_k{I}_k^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{h}_n{I}_n^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{p\in P}{\mu }_p{S}_p^t-\left({\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{P}_r^1{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{P}_r^2{Q}_{kl}^t\right)\end{array}$ (20)

$\tau \ge {\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{i\in I}{\displaystyle \sum }_{j\in J}{\theta }_{ij}{d}_{ij}{\tilde{q}}_i^t{Y}_{ij}^t$ (21)

约束(2)~(19)。

4.2. 动力电池回收量的不确定性处理

对于模型中的废旧动力电池回收量不确定，本节将采用前述的Budget鲁棒优化来对不确定参数${\tilde{q}}_i^t$ 进行处理，以得到鲁棒对等模型。

以上鲁棒优化模型中，约束(4)、(11)与(21)含有不确定参数${\tilde{q}}_i^t$ ，首先对约束(4)进行鲁棒优化处理。定义${\tilde{q}}_i^t$ 的取值范围为$\left[{\overline{q}}_i^t-{\widehat{q}}_i^t,{\overline{q}}_i^t+{\widehat{q}}_i^t\right]$ ，其中${\overline{q}}_i^t$ 与${\widehat{q}}_i^t$ 分别表示为退役动力电池回收量的名义值与最大偏离值，且${\widehat{q}}_i^t=\mu {\overline{q}}_i^t$ ，$\mu \in \left[0,1\right]$ ，为扰动系数。引入参数$\left|Q^t\right|$ 来控制鲁棒优化模型的保守程度，其中 $Q^t=\left\{i|{\widehat{q}}_i^t>0\right\}$ 表示回收量参数中发生不确定的个数。通过参数$\left|Q^t\right|$ 可以调节鲁棒模型的保守性，当$\left|Q^t\right|=0$ 时即当前周期模型不存在不确定参数，模型退化为确定性模型；当$\left|Q^t\right|>0$ 且不断增加时，模型的保守性提高，但相应的最优解由于考虑了所有不确定性会变得很差。

引入budget控制参数${\text{Γ}}_1^t$ ，其中${\text{Γ}}_1^t\in \left[0,\left|Q^t\right|\right]$ ，${\text{Γ}}_1^t$ 并不一定是整数，可以取范围内的任意实数。令${\text{Γ}}_1^t$ 的整数部分即$\lfloor {\text{Γ}}_1^t\rfloor$ (不超过${\text{Γ}}_1^t$ 的最大整数)个回收量参数取最差情况，即这部分的参数的变动值为${\widehat{q}}_i^t$ ，那么余下的回收量的变动值为$\left({\text{Γ}}_1^t-\lfloor {\text{Γ}}_1^t\rfloor \right){\widehat{q}}_i^t$，所以在这种情况下，最多会出现$\lfloor {\text{Γ}}_1^t\rfloor +1$ 个不确定回收量。

根据上述分析，回收量的最大偏离值${\widehat{q}}_i^t$ 的鲁棒控制参数为${\text{Γ}}_1^t$ ，则约的鲁棒表达式可转化为：

${\displaystyle \sum }_{i\in I}{\overline{q}}_i^t{Y}_{ij}^t-\underset{\left\{S_t\cup \left\{r\right\}|S_t\subseteq Q^t,\left|S_t\right|=\lfloor {\text{Γ}}_1^t\rfloor ,r\in Q^t\backslash S_t\right\}}{\max }\left\{{\displaystyle \sum }_{i\in S_t}{\widehat{q}}_i^t{Y}_{ij}^t+\left({\text{Γ}}_1^t-\lfloor {\text{Γ}}_1^t\rfloor \right){\widehat{q}}_r^t{Y}_{rj}^t\right\}\ge {\displaystyle \sum }_{k\in K}{Q}_{jk}^t.$ (22)

引入保护函数$\gamma \left(Y,{\text{Γ}}_1^t\right)$ ，其中${\text{Γ}}_1^t$ 是回收量最大偏移值的鲁棒控制参数。

$\gamma \left(Y,{\text{Γ}}_1^t\right)=\underset{\left\{S_t\cup \left\{r\right\}|S_t\subseteq Q^t,\left|S_t\right|=\lfloor {\text{Γ}}_1^t\rfloor ,r\in Q^t\backslash S_t\right\}}{\max }\left\{{\displaystyle \sum }_{i\in S_t}{\widehat{q}}_i^t{Y}_{ij}^t+\left({\text{Γ}}_1^t-\lfloor {\text{Γ}}_1^t\rfloor \right){\widehat{q}}_r^t{Y}_{rj}^t\right\}$ (23)

按照式(23)，将式(22)表示为：

${\displaystyle \sum }_{i\in I}{\overline{q}}_i^t{Y}_{ij}^t-\gamma \left(Y,{\text{Γ}}_1^t\right)\ge {\displaystyle \sum }_{k\in K}{Q}_{jk}^t;\forall j\in J,t\in T$ (24)

根据Budget鲁棒优化方法，通过对偶理论将式(22)含有max的表达式等价转换为以下模型，以下模型中${\rho }_i^t$ 与${\lambda }_{ij}^t$ 为对偶变量。

${\displaystyle \sum }_{i\in I}{\overline{q}}_i^t{Y}_{ij}^t-{\rho }_j^t{\text{Γ}}_1^t-{\displaystyle \sum }_{i\in Q^t}{\lambda }_{ij}^t\ge {\displaystyle \sum }_{k\in K}{Q}_{jk}^t;\forall j\in J,t\in T$ (25)

${\rho }_j^t+{\lambda }_{ij}^t\ge {\widehat{q}}_i^t{Y}_{ij}^t,i\in Q^t,\forall j\in J,t\in T$ (26)

${\rho }_j^t\ge 0,{\lambda }_{ij}^t\ge 0,i\in Q^t,\forall j\in J,t\in T$ (27)

下面证明模型(25)~(27)等价于约束(4)。首先引入实数决策变量${\psi }_{it}\in \left[0,1\right],i\in Q^t$ ，该变量表示是否从$Q^t$ 中选取参数${\widehat{q}}_i^t$ ，即该变量可以用于表示不确定参数的波动情况，所以(23)的保护函数可以等价为：

$\begin{array}{l}\gamma \left(Y,{\Gamma }_1^t\right)=\max {\displaystyle \sum }_{i\in I}{\widehat{q}}_i^t{\psi }_{it}{Y}_{ij}^t\\ {\displaystyle \sum }_{i\in Q^t}{\psi }_{it}\le {\Gamma }_1^t\\ 0\le {\psi }_{it}\le 1,\forall i\in Q^t\end{array}$ (28)

以上模型可知引入的决策变量${\psi }_{it}$ 表示了不确定参数的变化。对于$Q^t$ 内的任意i，当${\psi }_{it}=0$ 时，鲁棒控制参数${\text{Γ}}_1^t=0$ ，即所有参数未出现不确定，此时模型退化为确定性模型；当${\psi }_{it}=1$ 时，此时${\text{Γ}}_1^t=\left|Q^t\right|$ 时，即所有的不确定参数都发生了最大偏移，此时的模型最优解为最坏情况下的，保守程度过高。

按照对偶理论，将模型(28)转化为：

$\begin{array}{l}\min {\rho }_j^t{\Gamma }_1^t+{\displaystyle \sum }_{i\in Q^t}{\lambda }_{ij}^t\\ {\rho }_j^t+{\lambda }_{ij}^t\ge {\widehat{q}}_i^t{Y}_{ij}^t,i\in Q^t\\ {\rho }_j^t\ge 0,{\lambda }_{ij}^t\ge 0,i\in Q^t\end{array}$ (29)

模型中${\rho }_i^t$ 与${\lambda }_{ij}^t$ 为对偶变量，根据强对偶理论，模型(28)与(29)的最优目标值是相等的。所以将模型(29)的目标代入到(24)中的保护函数，等价转换后即得到模型(25)，即模型(25)与(24)等价，而(24)是由(4)推导得到的，所以三者都是等价的。

同理证明可得，式(11)等价于以下模型(30)~(32)：

${\displaystyle \sum }_{i\in I}{\overline{q}}_i^t{Y}_{ij}^t+{\xi }_j^t{\text{Γ}}_1^t+{\displaystyle \sum }_{i\in Q^t}{\kappa }_{ij}^t\le C_j{X}_j;\forall j\in J,t\in T$ (30)

${\xi }_j^t+{\kappa }_{ij}^t\ge {\widehat{q}}_i^t{Y}_{ij}^t,\forall i\in Q^t,\forall j\in J,t\in T$ (31)

${\kappa }_{ij}^t\ge 0,{\xi }_j^t\ge 0,\forall i\in Q^t,\forall j\in J,t\in T$ (32)

其中，${\xi }_j^t$ 与${\kappa }_{ij}^t$ 为对应保护函数max表达式的对偶变量。

类似地，可以证明式(21)等价于以下模型(33)~(35)，其中${\phi }^t$ 与${\vartheta }_i^t$ 为对应的对偶变量。

$\tau \ge {\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{i\in I}{\displaystyle \sum }_{j\in J}{\theta }_{ij}{d}_{ij}{\overline{q}}_i^t{Y}_{ij}^t+{\displaystyle \sum }_{t\in T}{\phi }^t{\text{Γ}}_1^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{i\in Q^t}{\vartheta }_i^t$ (33)

${\phi }^t+{\vartheta }_i^t\ge {\widehat{q}}_i^t{\displaystyle \sum }_{j\in J}{\theta }_{ij}{d}_{ij}{Y}_{ij}^t,\forall i\in Q^t$ (34)

${\vartheta }_i^t\ge 0,{\phi }^t\ge 0,\forall i\in Q^t,\forall t\in T$ (35)

4.3. 梯次市场需求量的不确定性处理

与上一小节类似，以上鲁棒优化模型中，约束(15)含有不确定参数${\tilde{q}}_p^t$ 。定义${\tilde{q}}_p^t$ 的不确定集为 $U_p=\left\{{\tilde{q}}_p^t\in \left[{\overline{q}}_i^t-{\widehat{q}}_i^t,{\overline{q}}_i^t+{\widehat{q}}_i^t\right],\forall p\in P\right\}$ ，其中${\overline{q}}_p^t$ 与${\widehat{q}}_p^t$ 分别表示为梯次利用市场的需求量的名义值与最大偏离值，且${\widehat{q}}_p^t=\theta {\overline{q}}_p^t$ ，$\theta \in \left[0,1\right]$ ，为扰动系数。引入参数$\left|D^t\right|$ 来控棒优化模型的保守程度，其中$D^t=\left\{i|{\widehat{q}}_p^t>0\right\}$ 表示市场需求量参数中发生不确定的个数。通过参数$\left|D^t\right|$ 可以调节鲁棒模型的保守性，当$\left|D^t\right|=0$ 时即当前周期模型不存在不确定参数，此时为确定性模型；当$\left|D^t\right|>0$ 且不断增加时，模型的保守性提高。类似地，引入budget控制参数${\text{Γ}}_2^t$ ，其中${\text{Γ}}_2^t\in \left[0,\left|D^t\right|\right]$ ，${\text{Γ}}_2^t$ 并不一定是整数，可以取范围内的任意实数。与之前的处理类似，令${\text{Γ}}_2^t$ 的整数部分即${\text{Γ}}_1^t$ 个参数取最差情况，即变动值为${\widehat{q}}_p^t$ ，那么余下的回收量的变动值为$\left({\text{Γ}}_2^t-\lfloor {\text{Γ}}_2^t\rfloor \right){\widehat{q}}_p^t$。因为最大偏离值${\widehat{q}}_p^t$ 的鲁棒控制参约${\Gamma }_2^t$ ，约束(15)的鲁棒表达式可转化为：

${\displaystyle \sum }_{n\in N}{W}_{np}^t+S_p^t\ge {\overline{q}}_p^t+\frac{{\Gamma }_2^t}{\left|D^t\right|}{\widehat{q}}_i^t;\forall t\in T,p\in D^t$ (36)

4.4. 鲁棒对等模型

经过上述budget鲁棒优化的推导，得到考虑客户废旧电池回收量与梯次市场需求量不确定的逆向物流回收网络的鲁棒对等模型：

$\begin{array}{l}\min {\displaystyle \sum }_{j\in J}{f}_j{X}_j+{\displaystyle \sum }_{k\in K}{f}_k{X}_k+\tau \\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{j\in J}{\displaystyle \sum }_{k\in K}{\theta }_{jk}{d}_{jk}{Q}_{jk}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{\theta }_{kl}{d}_{kl}{Q}_{kl}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{\theta }_{kn}{d}_{kn}{Q}_{kn}^t\\ +{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{\theta }_n{d}_n{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N}{b}_n{Q}_{kn}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{h}_k{I}_k^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{h}_n{I}_n^t\\ \begin{array}{c}+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{p\in P}{\mu }_p{S}_p^t-\left({\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{n\in N}{\displaystyle \sum }_{p\in P}{P}_r^1{W}_{np}^t+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{l\in L}{P}_r^2{Q}_{kl}^t\right)\end{array}\end{array}$ (37)

s.t. 式(2)~(3), (5)~(10), (12)~(14), (16)~(19), (25)~(27), (30)~(36)

5. 算例分析

5.1. 算例设置

由于引入不确定参数，模型求解难度增加，所以本节构造包含24个新能源汽车客户点、12个备选收运点、8个备选分类检测中心、3个梯次利用中心、3个回收拆解中心与15个梯次利用客户点在内的逆向物流回收网络算例。该算例研究3个周期下的选址与流量分配决策，假定模型的两个不确定参数的名义值已知，且在区间内按照均匀分布随机生成，其余相关参数设置如表2所示。

Table 2. Model parameter setting

表2. 模型参数设置

参数	符号	单位	数值
固定成本	$f_j$	元	$U\left[5,6.5\right]\times {10}^5$
	$f_k$		$U\left[3,8\right]\times {10}^5$
库存成本	$h_k$	元/kg	$U\left[5,10\right]$
	$h_n$		$U\left[5,10\right]$

续表

最大容量	$C_j$	组	800
	$C_k$		600
	$C_n$		1000
	$C_l$		1000
处理成本	$b_n$	元/kg	25
惩罚成本	${\mu }_p$	元/kg	100
网络收益	$P_r^1$	元/kg	77
	$P_r^2$		18

5.2. 决策分析

通过上述Budget鲁棒优化方法处理，去除了目标函数与约束条件的不确定性，推导得到的鲁棒对等模型为易于求解的混合整数线性规划模型。针对鲁棒对等模型，使用Python调用Gurobi求解器对模型进行求解。

令鲁棒控制参数${\text{Γ}}_1^t$ 与${\text{Γ}}_2^t$ 为0，即所有不确定参数都为名义值，此时得到确定性模型。如图2所示，表示当鲁棒控制参数都为0时得到的动力电池回收网络选址结果，收运点的选址决策为1，3，4，6，8，10，12，分类检测中心的选址决策为1，4，8，该逆向物流回收网络的总成本为31717959.44。

Figure 2. Site selection result of multi-cycle power battery recycling network

图2. 多周期动力电池回收网络选址结果

设置模型不确定参数的鲁棒控制参数与扰动系数，求解不确定情况下的模型。

首先令梯次市场需求量是确定的，即${\text{Γ}}_2^t=0$ ，求解扰动系数$\mu =0.5$ 下的模型，确定性与不确定性模型选址结果与网络成本如表3所示。

Table 3. Robust optimization decision results with $\mu =0.5$

表3. $\mu =0.5$ 的鲁棒优化决策结果

${\text{Γ}}_1^t$	收运点选址结果	分类检测中心选址结果	回收网络总成本
0	1, 3, 4, 6, 8, 10, 12	1, 4, 8	31717959.44
8	1, 3, 4, 6, 8, 10, 11, 12	1, 3, 4	60436502.98
12	1, 3, 4, 6, 8, 10, 11, 12	3, 4	62924225.10
24	1, 2, 3, 4, 6, 8, 10, 11, 12	3, 4	66943709.81
18	1, 3, 4, 6, 8, 10, 11, 12	3, 4	66074897.87
24	1, 2, 3, 4, 6, 8, 10, 11, 12	3, 4	66943709.81

令动力电池回收量是确定的，即${\text{Γ}}_1^t=0$ ，求解扰动系数$\theta =0.5$ 下的模型，确定性与不确定性模型选址结果与网络成本如表4所示。

Table 4. Robust optimization decision results with $\theta =0.5$

表4. $\theta =0.5$ 的鲁棒优化决策结果

${\text{Γ}}_2^t$	收运点选址结果	分类检测中心选址结果	回收网络总成本
0	1, 3, 4, 6, 8, 10, 12	1, 4, 8	31717959.44
5	1, 3, 4, 6, 8, 10, 12	1, 4, 8	33722767.24
10	1, 3, 4, 6, 8, 10, 11, 12	1, 4, 8	35754078.85
15	1, 3, 4, 6, 8, 10, 11, 12	1, 3, 4, 8	38936122.15

当模型引入了不确定性时，网络总成本相较于确定性模型显著增加，且收运点与分类检测中心的选址结果也发生了变化，这是因为为了应对动力电池回收量与梯次利用市场需求量的不确定性，决策者倾向于做出保守的选择，而保守的结果也带来了成本的增加，这是为了应对模型不确定性所付出的代价。

当${\text{Γ}}_1^t$ 增大时，意味着每周期产生的动力电池回收量的不确定性增加，为应对增加的回收量不确定性，决策者需要调整收运点与分类检测中心的选址规划，相应地设施间的流量分配也会发生变化，因此带来了回收网络的总成本不断增加。当${\text{Γ}}_1^t=8$ 时，增加了收运点11，分类检测中心的选址从1，4，8变为1，3，4。有8个新能源汽车客户点的回收量发生波动，决策时需要增加收运点以满足不确定回收量的回收需求。而随着收运点的选址决策的改变，分类检测中心也做出了选址调整，关闭了原先的8，增设了4，由于确定性模型中分类检测中心的容量尚有剩余，所以总的数量并未发生变化，而因为原先的分类检测中心8无法将引入不确定带来的回收量全部处理掉，所以在收运点11附近新增了分类检测中心4。

当${\text{Γ}}_2^t$ 增大时，意味着每周期梯次利用市场的需求量发生变化的数量增多，为应对可能增加的梯次市场的需求，决策者需要增加设施数量，调整选址决策以提高对废旧动力电池的回收能力，同样地决策趋于保守，回收网络总成本也将随着${\text{Γ}}_2^t$ 的变化而不断增大。当${\text{Γ}}_2^t\le 5$ 时，鲁棒优化模型的决策结果为：收运点选址1，3，4，6，8，10，12，分类检测中心选址1，4，8；当${\text{Γ}}_2^t\ge 10$ 时收运点选址发生变化为1，3，4，6，8，10，11，12，并且当${\text{Γ}}_2^t=15$ 时分类检测中心选址也发生了变化为1，3，4，8。当${\text{Γ}}_2^t\le 5$ 时，相比于确定性模型，选址决策上并没有发生变化，即当前的选址方案可以应对变化的梯次市场的需求波动，但是变化的需求仍然会带来总成本的增加，因为设施间的流量会随着鲁棒控制参数的增大而增加。随着鲁棒控制参数${\text{Γ}}_2^t$ 的继续增加，收运点的选址决策首先发生改变，为了波动的梯次市场需求，需要增设更多的收运点以提升废旧电池的回收能力，接着分类检测中心的选址数量也增加了，增加8是为了提高废旧电池的处理能力，增加梯次利用电池的数量以满足可能变化的梯次利用客户需求。

5.3. 灵敏度分析

以上决策分析部分，是在固定参数扰动比例的情况下，探究鲁棒水平对于选址决策的影响，为进一步研究不确定性对于鲁棒优化模型的影响程度，下面对不确定参数的扰动系数与鲁棒控制参数进行灵敏度分析。

令${\text{Γ}}_2^t=0$ ，$\theta =0$ ，模型最优值的鲁棒性主要由鲁棒控制参数${\text{Γ}}_1^t$ 与扰动系数$\mu$ 决定，取 ${\text{Γ}}_1^t=0,3,6,12,15,18,24$ ，$\mu =0.05,0.1,0.2$ ，分别求解模型得到最优目标值，如图3所示。同理，令${\text{Γ}}_1^t=0$ ，$\mu =0$ ，取${\text{Γ}}_2^t=0,3,6,9,12,15$ ，$\theta =0.05,0.1,0.2$ ，分别求解模型得到最优目标值，如图4所示。

Figure 3. Total network cost under different robust levels of recycling

图3. 不同回收量的鲁棒水平下的网络总成本

Figure 4. Total network cost under different robustness levels of demand

图4. 不同需求量的鲁棒水平下的网络总成本

分析图3与图4可知，随着鲁棒控制参数的增大，不论是${\text{Γ}}_1^t$ 抑或${\text{Γ}}_2^t$ ，网络总成本都呈现增长趋势，并且成本增长速度都在逐渐放缓。这是由于随着参数不确定性的扩大，模型趋于达到全部参数发生波动的情况，此时模型的最优目标值的增大空间会越来越小。所以决策者可以权衡鲁棒水平与自身风险偏好，选择合适的鲁棒水平，以获得满意的最优目标值。其次在相同的鲁棒控制参数下，扰动系数越大，目标函数也即网络总成本也越大。这是由于波扰动系数越大，不确定性增加的越大，可能会带来选址的增加与库存和流量的增加，所以最终带来了网络总成本的增加。最后相同情况下，动力电池回收量的不确定性所带来的总成本的增加会比梯次市场需求的不确定性要明显。这是由于模型中目标函数不仅仅有成本项还有回收收益项，虽然为应对不确定风险，两者都会改变选址策略，相关的成本都会有所增加，但是梯次利用市场需求量不确定性还会带来一部分额外的需求，相对于回收量的不确定性，这会导致梯次利用的收益增加，会抵消一部分的回收网络成本的增加。

6. 结论

在“双碳”背景下，针对动力电池回收网络规划问题，本文考虑了新能源汽车客户动力电池回收量与梯次利用市场需求的不确定性，建立了多周期动力电池逆向物流网络的鲁棒优化模型。针对模型中的两个不确定参数，采用鲁棒优化方法进行处理。按照Budget鲁棒优化理论建立参数的不确定集合，并将含有不确定参数的目标函数与约束条件进行等价转化，通过对偶理论证明鲁棒对等模型与原模型的等价关系。然后构造随机算例，通过Python调用Gurobi求解模型，并对鲁棒控制参数与扰动系数进行灵敏度分析，分析不确定性对于决策结果的影响，以验证鲁棒优化方法在处理不确定优化问题上的有效性。本文在求解规模与效率上还有待完善，未来可以进一步设计高效的算法求解大规模算例下的鲁棒模型。

NOTES

^*通讯作者。