This section discusses the RT logistics center as a typical industry case study. Firstly, we analyze the original operation mode and the exposed problems. Then we propose an improved logistics terminal distribution mode. Finally, we construct a mathematical model based on the improved mode to solve the terminal distribution problem.
Describe the case and analyze the problem
The automotive scrap parts logistics center is responsible for receiving scrap parts from automotive service stations around the country and delivering them to the suppliers. RT logistics center mainly has the following functions: collection, sorting, inventory, storage and custody, and distribution. In the original distribution model, the RT logistics center adopts a pointtopoint batch pickup model with suppliers, as shown in Fig. 3. That is, when the suppliers’ materials in the logistics center reach the set inventory level, the logistics center will send a pickup notice to the supplier, and then the supplier will arrange its vehicle to pick up the materials and send them back to the factory. Through field research, we found the following problems in the logistics center under the original model:

(1)
The pointtopoint bulk pickup model leads to high inventory levels, which results in expensive inventory costs. At the same time, excessive inventory takes up a large amount of storage space, making it difficult to adapt to future growth within the rapidly expanding automotive market.

(2)
The limited scale of the supplier’s selfpickup model leads to high transportation costs and low logistics efficiency, resulting in high logistics costs and high levels of carbon emissions for the entire distribution network.

(3)
Since multiple suppliers are involved, it is difficult to coordinate the vehicle models and pickup times of each supplier, this can easily lead to confusion in the management of the logistics center, interfering with the normal operation status of the outbound storage link.
Designing improved distribution mode
Under the original distribution model, a direct way to reduce inventory levels was to increase the frequency of supplier pickups. However, the rational decisionmakers of the suppliers, are not willing to bear more logistics and transport costs while their interests remain unchanged. Therefore, this paper decides to improve the logistics system from the perspective of changing the pointtopoint transport mode to reduce the logistics inventory level and the overall logistics costs of the whole supply chain. It will also help to reduce the carbon emissions in the logistics process and support the implementation of green sustainable development policies.
The Milkrun model is a pointtogroup efficient delivery method, which many scholars^{27,28} have applied in production research. Studies have shown that the adoption of this model can reduce transport and inventory costs by increasing the loading rate of transport vehicles, thus achieving a reduction in the total costs of logistics. In view of its characteristics of “multifrequency, small batch, and fixed time window”, it can be a better solution to the problems existing in the RT logistics center under the original model. Therefore, this paper decides to establish a kind of circular distribution network based on Milkrun with the logistics center as the leader, as shown in Fig. 4, and then combines it into the MTGVRPTW for indepth study.
Model assumptions and symbol description
Model assumptions
In the actual transportation and distribution process, the vehicle will be affected by a variety of uncontrollable factors, so this paper makes the following assumptions about the mathematical model: (1) Sufficient capacity to take on the distribution needs of suppliers. (2) The transport process is not affected by weather, traffic control, travel peaks, etc., and always maintains the set average speed at an even pace. (3) After each trip delivery departs from the logistics center, it serves the customers in turn according to the optimization results and returns to the logistics center upon completing the delivery task. (4) Each distribution trip can serve multiple suppliers, but each supplier set cannot be delivered by multiple distribution trips. (5) Distribution vehicles are subject to the double limitation of carrying capacity and loading space, which cannot exceed the constraints limitations, and the loading space is expressed in the form of the number of loading units that can be loaded. (6) The vehicle stays at each site for the same time.
Symbol description
The symbols used in the MTGVRPTW model constructed in this paper and their related descriptions are shown in Tables 1 and 2.

(1)
Sets

(2)
Parameters

(3)
Decision variables
$$x_{ijk}^{m} = \left\{ {\begin{array}{*{20}l} {1,} \hfill & {{\text{A}}\;{\mathrm{direct}}\;{\mathrm{path}}\;{\mathrm{exists}}\;{\mathrm{between}}\;{\mathrm{arc}}\;{(}i{, }j{)}\;{\mathrm{and}}\;{\mathrm{belongs}}\;{\mathrm{to}}\;{\mathrm{the}}\;{\mathrm{trip}}\;k\;{\mathrm{of}}\;{\mathrm{vehicle}}\;m} \hfill \\ {0,} \hfill & {{\mathrm{otherwise}}} \hfill \\ \end{array} } \right.$$
$$y_{km} = \left\{ {\begin{array}{*{20}l} {1,} \hfill & {{\mathrm{Distribution}}\;{\mathrm{task}}\;{\mathrm{for}}\;{\mathrm{trip}}\;k\;{\mathrm{is}}\;{\mathrm{carried}}\;{\mathrm{by}}\;{\mathrm{vehicle}}\;m} \hfill \\ {0,} \hfill & {{\mathrm{otherwise}}} \hfill \\ \end{array} } \right.$$
$$r_{m} = \left\{ {\begin{array}{*{20}l} {1,} \hfill & {{\mathrm{The}}\;{\mathrm{vehicle}}\;m\;{\mathrm{of}}\;{\mathrm{the}}\;{\mathrm{set}}\;{\mathrm{of}}\;{\mathrm{vehicles}}\;{\mathrm{is}}\;{\mathrm{in}}\;{\mathrm{service}}} \hfill \\ {0,} \hfill & {{\mathrm{otherwise}}} \hfill \\ \end{array} } \right.$$
$$ds_{i} = \left\{ \begin{gathered} 0,\quad ET_{i} \le t_{a\_i} < LT_{i} \, \hfill \\ 1,\quad {\mathrm{otherwise}} \hfill \\ \end{gathered} \right.$$
Mathematical model
Objective functions
In the context of cleaner production and sustainable development strategies, when building the MTGVRPTW model for automotive scrap parts distribution, it is necessary to consider the impacts of various factors, including: (1) Reducing ecological impacts. (2) Achieving cost reduction and efficiency in logistics. (3) Ensuring the timely delivery of distribution services. On this basis, this study proposes three different optimization objectives.

(1)
Minimizing carbon dioxide emissions
The ecological impact of the vehicle distribution process is usually caused by the fact that driving a vehicle consumes fuel and produces carbon dioxide. Therefore, the first optimization objective in the model is set to minimize carbon dioxide emissions during the delivery process. Zhou et al. summarized most of the estimation methods on carbon emissions from automobiles^{29}, considering the difficulty of obtaining data, this paper decided to use the fuel consumption rate measure to calculate. The specific formula is shown in Eq. (1).
$$Min\;{\mathrm{Z}}_{1} = \sum\limits_{i} {\sum\limits_{j} {\sum\limits_{k} {\sum\limits_{m} {e_{0} x_{ijk}^{m} D_{ij} \left( {\rho_{ok} + \frac{{\left( {\rho_{k}^{*} – \rho_{ok} } \right)}}{{W_{V} }}W_{ijk}^{m} } \right)} } } }$$
(1)

(2)
Minimizing overall logistics costs
In the process of logistics distribution, a variety of transportation resources need to be deployed, which will generate several costs, including internal preparation, vehicle rental, driver labor, and distance transportation. The goal of the actual operation of the enterprise is to reduce the overall logistics costs and improve the efficiency of resource utilization. Therefore, this paper sets the lowest overall logistics costs as the second optimization objective. The specific formula is shown in Eq. (2).
$$Min\;{\mathrm{Z}}_{2} = C_{preparation} + C_{vehicle} + C_{driver} + C_{transport}$$
(2)

(i)
Preparation costs for departure
Before the departure of each trip within the enterprise involved in the shelves, moving storage, loading and other logistics arrangements, including a large number of logistics equipment, material and human resources, the unified deployment. The costs of this item is shown in Eq. (3).
$$C_{preparation} = C_{p} \sum\limits_{{j \in {\varvec{N}}^{\prime}}} {x_{0jk}^{m} }$$
(3)

(ii)
Vehicle rental costs
Distribution vehicles required for distribution are obtained by leasing with thirdparty companies, which generates vehicle leasing costs. The final decision on the number of vehicles to be rented is related to the result of combining multiple trips. The costs of this item is shown in Eq. (4).
$$C_{vehicle} = C_{v} \sum\limits_{m} {r_{m} }$$
(4)

(iii)
Driver labor costs
Each vehicle is provided with a driver, who is employed on a temporary basis. If the actual delivery time of the vehicle m is less than the legal working hours of half a working day (4 h), a halfday contract is concluded with the driver of the vehicle; otherwise, a fullday contract is concluded. The costs of this item is shown in Eq. (5).
$$C_{driver} = \sum\limits_{m} {CH_{m} }$$
(5)
where CH_{m} represents the labor costs of equipping the vehicle m with a driver and is related to the travel time of each vehicle. It is shown in Eq. (6).
$$CH_{m} = \left\{ {\begin{array}{*{20}l} {100,\quad T_{m} \le 240\,{\mathrm{(minute)}}\;{\mathrm{\& }}\;r_{m} = 1 \, } \\ 200,\quad T_{m} > 240\,{\mathrm{(minute)}}\;{\mathrm{\& }}\;r_{m} = 1 \, \\ {0},\quad {\mathrm{otherwise }} \\ \end{array} } \right.$$
(6)
where the formula for the travel time of each vehicle is as described in Eq. (7).
$$T_{m} = \sum\limits_{k} {y_{km} } t_{DC} + \sum\limits_{k} {y_{km} } Q_{k} t_{CS} + \sum\limits_{i} {\sum\limits_{j} {\sum\limits_{k} {\left( {x_{ijk}^{m} D_{ij} } \right)/v} } } \,$$
(7)

(iv)
Vehicle transportation costs
Transportation costs will be generated when transporting goods, which include fuel costs, etc. And it is directly proportional to the distance traveled. The costs of this item is shown in Eq. (8).
$$C_{transport} = \sum\limits_{i} {\sum\limits_{j} {\sum\limits_{k} {\sum\limits_{m} {c_{1} x_{ijk}^{m} D_{ij} } } } }$$
(8)

(3)
Minimizing delayed delivery rates
If the delivery vehicle arrives at the destination earlier or later than the time window required by suppliers, it may disrupt the supplier’s normal work situation. Therefore, we expect the delayed delivery rate to be minimized to avoid disrupting the supplier’s work schedule due to unfavorable delivery. The specific formula is shown in Eq. (9).
$$Min\;Z_{3} = \frac{{\sum\limits_{i} {ds_{i} } }}{n}$$
(9)
subject to
$$\sum\limits_{m} {\sum\limits_{k} {\sum\limits_{{j \in {{\varvec{\Delta}}}^{ + } (i)}} {x_{ijk}^{m} } } = 1\quad \forall i \in {\varvec{N}}}$$
(10)
$$\sum\limits_{j} {x_{0jk}^{m} } = \sum\limits_{i} {x_{i0k}^{m} } \quad \forall k \in {\varvec{K}}{;}\quad \forall m \in {\varvec{M}}$$
(11)
$$ \sum\limits_{i \in {\Delta_{ – }} (j)} {x_{ijk}^{m} } – \sum\limits_{i \in {\Delta^{ + }} (j)} {x_{jik}^{m} } = 0\quad \forall j \in {\varvec{N}};\quad \forall k \in {\varvec{K}};\quad \forall m \in {\varvec{M}}$$
(12)
$$ \sum\limits_{i} {\sum\limits_{j} {q_{i} x_{jik}^{m} } \le Q_{V} } \quad \forall k \in {\varvec{K}};\quad \forall m \in {\varvec{M}}$$
(13)
$$\sum\limits_{i} {\sum\limits_{j} {W_{ijk}^{m} x_{jik}^{m} } \le } \, W_{V} \quad \forall k \in {\varvec{K}};\quad \forall m \in {\varvec{M}}$$
(14)
$$\sum\limits_{k} {\sum\limits_{i} {\sum\limits_{j} {D_{ij} x_{jik}^{m} } \le D} \quad \forall m \in {\varvec{M}}}$$
(15)
$$\sum\limits_{i} {\sum\limits_{j} {\sum\limits_{k} {x_{ijk}^{m} } } \le I \times r_{m} \quad \forall m \in {\varvec{M}}}$$
(16)
$$\sum\limits_{j} {x_{0jk}^{m} } \ge \sum\limits_{j} {x_{0jk}^{m + 1} } \quad \forall k \in {\varvec{K}}{;}\quad \forall m \in {\varvec{M}}/\left\{ {M^{*} } \right\}$$
(17)
$$x_{ijk}^{m} \in \{ 0,1\}$$
(18)
$$r_{m} \in \{ 0,1\}$$
(19)
$$ds_{i} \in \{ 0,1\}$$
(20)
where constraint Eq. (10) restricts the allocation of each site to only one trip. Constraint Eq. (11) limits the number of times a single vehicle can enter and exit the logistics center under multitrip distribution to the same number of times. Constraint Eq. (12) restricts the number of times a vehicle can drive in and out of the same stop to remain equal. Constraint Eq. (13) restricts that no distribution task on each trip exceeds the volume limit. Constraint Eq. (14) restricts each distribution task on each trip to not exceeding the load limits. Constraint Eq. (15) limits the number of miles traveled by each carrier vehicle to no more than the vehicle’s range limit. Constraint Eq. (16) restricts only vehicles that are in service to distribution duties. Constraint Eq. (17) restricts the order in which vehicles are put into service from m to m + 1. The constraint Eqs. (18–20) defines the decision variable as 0 or 1.
Multiobjective processing
There are three objectives of different properties in the optimization model. Considering the complexity of the solution and the decision preference, this study uses the weighted sum method to integrate these objectives into a single objective function^{30}. Setting the weights of the three objectives as convex combinations, i.e.,\(\lambda_{1} ,\lambda_{2} ,\lambda_{3} \ge 0\) and \(\lambda_{1} + \lambda_{2} + \lambda_{3} = 1\). Due to the large number of members in a distribution network and the fact that distribution costs are shared by all members, decisionmaking must take into account the opinions of all members. Groupanalytic hierarchy process (GAHP) is a comprehensive evaluation method developed based on hierarchical analysis, which can effectively integrate the knowledge and experience of multiple experts for group decisionmaking^{31}. Therefore, this paper will apply GAHP to determine the weights of the three objective functions, drawing on the reviews of a team of experts consisting of representatives from logistics centers and suppliers.

(1)
Calculate the weight vector of each expert’s judgment matrix
In the formulation of the weight standard, a total of 5 experts participate in the group decision, and the judgment matrix constructed by the review opinion of each expert are \(A_{1} ,A_{2} ,A_{3} ,A_{4} ,A_{5}\). The expression form of each judgment matrix is \(A_{l} = (a_{ijl} ); \, i,j = 1,2,3; \, l = 1,2, \ldots ,5\), where \(a_{ijl}\) denotes the relative importance of factor i over factor j as perceived by the expert l. Separately solve their weight vectors \(\tilde{w}_{il} \, = \,(\,\tilde{w}_{1l} ,\tilde{w}_{2l} \,,\tilde{w}_{3l} \,)^{T}\), where \(\tilde{w}_{il}\) represents the judgment weight value of the expert l for the objective function i. To ensure the consistency of the weight calculation results, a consistency check is performed, aiming for \(CR_{l} = CI_{l} /RI_{l} < 0.1\).

(2)
Calculate the group’s composite weight vector
In this paper, considering the fairness and balance, under the condition that \(\sum {\lambda_{l} = 1,\quad (\lambda_{l} > 0,\quad l = 1,2, \ldots ,5)}\), the weight of each review expert’s opinion is set to be \(\lambda_{l} = 1/5 = 0.2\). Performing the weighted arithmetic mean calculation on the respective components of each weight vector, as shown in Eq. (21).
$$\tilde{w}_{i} = \sum\limits_{l} {\lambda_{l} \cdot \tilde{w}_{il} }$$
(21)
Following normalization of \(\tilde{w}_{i}\) as Eq. (22), the weight vectors for the three objectives can be derived as \(w_{i} = (0.164,0.539,0.297)^{T}\).
$$w_{i} = \frac{{\tilde{w}_{i} }}{{\sum\limits_{i} {\tilde{w}_{i} } }}$$
(22)
To address the dimensional differences among \(Z_{1}\) and \(Z_{2}\), the method of min–max normalization is employed to transform the multiobjective problem into a scalar optimization problem. Where \(\underline{{Z_{1} }} ,\underline{{Z_{2} }} ,\underline{{Z_{3} }}\) are the lower bounds of \(Z_{1} ,Z_{2} ,Z_{3}\), and \(\overline{{Z_{1} }} ,\overline{{Z_{2} }} ,\overline{{Z_{3} }}\) are the upper bounds of \(Z_{1} ,Z_{2} ,Z_{3}\). The lower and upper bounds are calculated by the box constraint. The transformed single objective optimization model is shown in Eq. (23).
$$\begin{aligned} Min\;Z & = w_{1} \left( {\frac{{Z_{1} – \underline{{Z_{1} }} }}{{\overline{{Z_{1} }} – \underline{{Z_{1} }} }}} \right) + w_{2} \left( {\frac{{Z_{2} – \underline{{Z_{2} }} }}{{\overline{{Z_{2} }} – \underline{{Z_{2} }} }}} \right) + w_{3} \left( {\frac{{Z_{3} – \underline{{Z_{3} }} }}{{\overline{{Z_{3} }} – \underline{{Z_{3} }} }}} \right) \\ & = 0.164\left( {\frac{{Z_{1} – \underline{{Z_{1} }} }}{{\overline{{Z_{1} }} – \underline{{Z_{1} }} }}} \right) + 0.539\left( {\frac{{Z_{2} – \underline{{Z_{2} }} }}{{\overline{{Z_{2} }} – \underline{{Z_{2} }} }}} \right) + 0.297\left( {\frac{{Z_{3} – \underline{{Z_{3} }} }}{{\overline{{Z_{3} }} – \underline{{Z_{3} }} }}} \right) \\ \end{aligned}$$
(23)