*Article* **A Multi-Period Vehicle Routing Problem for Emergency Perishable Materials under Uncertain Demand Based on an Improved Whale Optimization Algorithm**

**Xiaodong Li 1,†, Yang Xu 1,\*,†, Kin Keung Lai 2,†, Hao Ji 1, Yaning Xu <sup>1</sup> and Jia Li <sup>3</sup>**


**Abstract:** The distribution of emergency perishable materials after a disaster, such as an earthquake, is an essential part of emergency resource dispatching. However, the traditional single-period distribution model can hardly solve this problem because of incomplete demand information for emergency perishable materials in affected sites. Therefore, for such problems we firstly construct a multi-period vehicle path distribution optimization model with the dual objectives of minimizing the cost penalty of distribution delay and the total corruption during delivery, and minimizing the total amount of demand that is not met, by applying the interval boundary and most likely value weighting method to make uncertain demand clear. Then, we formulate the differential evolutionary whale optimization algorithm (DE-WOA) combing the differential evolutionary algorithm with the whale algorithm to solve the constructed model, which is an up-and-coming algorithm for solving this type of problem. Finally, to validate the feasibility and practicality of the proposed model and the novel algorithm, a comparison between the proposed model and the standard whale optimization algorithm is performed on a numerical instance. The result indicates the proposed model converges faster and the overall optimization effect is improved by 23%, which further verifies that the improved whale optimization algorithm has better performance.

**Keywords:** emergency material distribution; multi-period; uncertain demand; perishable materials; whale optimization algorithm; differential evolution algorithm

**MSC:** 90B06

#### **1. Introduction**

Large-scale sudden natural or man-made disasters occur frequently around the world every year, posing serious threats and impacts on society, human production, and life [1]. How to respond quickly effectively to these unpredictable emergencies has attracted much attention from governments and management at all levels, and has also placed high requirements on them from all aspects [2].

A scientific distribution and reasonable delivery of emergency relief materials, a key aspect of emergency relief work, can reduce the damage to property and casualties caused by disasters, improve the efficiency of rescue work and release the psychological pressure on the victims [3–5]. Due to the suddenness of disasters and the urgency of rescues, the demand for perishable emergency supplies for affected locations is often vague [6]. In addition, the longer the transport time, the more serious the spoilage phenomenon. Currently, this problem can be solved by a single-period delivery model. However, in this case,

**Citation:** Li, X.; Xu, Y.; Lai, K.K.; Ji, H.; Xu, Y.; Li, J. A Multi-Period Vehicle Routing Problem for Emergency Perishable Materials under Uncertain Demand Based on an Improved Whale Optimization Algorithm. *Mathematics* **2022**, *10*, 3124. https://doi.org/10.3390/ math10173124

Academic Editor: Yaping Ren

Received: 26 July 2022 Accepted: 27 August 2022 Published: 31 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

plenty of restrictions influence the solution's precision. For example, the actual demand for one disaster site is much greater than the maximum loading capacity of a vehicle, and the number of vehicles is limited. Thus, the single-period delivery cannot satisfactorily solve this problem. To more efficiently solve this problem, we consider a multi-period distribution model. Given the situation of sufficient supplies in the distribution center, one can use it to make decisions on the distribution vehicle's path and the distribution quantity in each period to minimize the cost penalty of distribution delay and total corruption during delivery, and minimize the total amount of demand that is not met, which is worthy of studying to improve relief work's efficiency and reduce losses in disaster areas.

The remainder of this paper is organized as follows. Section 2 performs literature reviews relevant to this study. Section 3 constructs the optimization model with the biobjective and multi-period vehicle path distribution, and proposes the improved differential whale optimization algorithm, which is a novel algorithm for solving the vehicle path problem with multi-objective optimization. Section 4 presents a practical example to verify the validity of the proposed model and algorithm. A comparison of the solution results of the algorithm before and after the improvement reveals that the improved differential evolutionary whale optimization algorithm optimizes better regarding the two objectives of minimizing the distribution delay penalty and corruption cost, and minimizing the unsatisfied degree of demand. Finally, conclusions and possible future research are given in Section 5.

#### **2. Literature Review**

The vehicle routing problem (VRP), as a classical problem in the field of operations research and combinatorial optimization, has been widely studied and played a significant role in transportation, logistics production and emergency rescue since its introduction in 1959 by Dantzig and Ramser [7]. A large number of experts and scholars have conducted in-depth research and analysis on it so far. Many variants of the VRP problem have been derived, and related theories and models have become relatively mature, among which, the multi-period vehicle path problem (PVRP) is one. Traditional vehicle paths and their derivatives are mostly deterministic vehicle path optimization problems, where the relevant variables are known in advance. However, in practice, uncertain information abounds whether in production transport or emergency relief, including demands, road conditions, casualties and so on. It can be divided into fuzzy information, random information and dynamic information concerning the properties of uncertain information. Therefore, the analysis and research of uncertain vehicle path optimization problems have become the focuses of experts and scholars. With the increasing development of intelligent optimization algorithms in recent years, a good research foundation has been laid for solving such problems.

The problem of optimizing the routes of emergency material distribution vehicles is a typical VRP problem in which the distribution center provides materials to some demand points with different quantities of materials, and vehicles are assigned to appropriate routes which form closed loops such that departure and final return are both the distribution center, so that the demand points' needs are met. Such goals of minimum transport costs, shortest driving distance and time spent under certain constraints should also be accounted for. Given the condition of the demand of distribution being known, the shortest driving distance is used as a goal to indicate the shortest resource allocation time, and a suitable distribution path is selected for a vehicle to satisfy the distribution demand of each affected location. Zhou et al. constructed a heterogeneous vehicle path optimization model for the vehicle path problem in which the pre-emergency transporters' vehicles are insufficient; the maximum system satisfaction and the minimum total time and the total cost were considered as the goals [8]. Li Zhuo et al., focusing on different interests of demand points and transporters, developed a multi-objective hybrid vehicle path optimization model with a soft time window, and a non-dominated sorting ant colony algorithm was proposed to

solve this model. An arithmetic case analysis indicated the effectiveness of the modified algorithm [9].

For the multi-site, open-emergency material distribution vehicle path problem considering secondary disasters, with the objective of shortest transportation time, Tan Jie et al. established two types of mathematical planning models that, under the one-sided fuzzy soft time window and fuzzy demand constraints, consider the risk of random failure at supply points, and designed an improved variable neighborhood search algorithm to solve the problem [10].

To solve the site-path problem of post-disaster emergency relief, different objectives and models were developed by scholars. With the objectives of maximizing rescue efficiency and minimizing the total cost, Gao Xinyu et al. developed a multi-stage site-path optimization model under the constraint of demand uncertainty and proposed an improved fast non-dominated genetic algorithm [11]. To maximize the matching degree of emergency demand at each dispatch point in the current stage, minimize the variance of the average matching degree of emergency demand at the previous k stages of dispatch and minimize the total travel time of the dispatch path, Liu Yang et al. constructed a multi-stage distribution and dispatching model for emergency relief supplies based on the historical travel time functions of road sections to portray the dynamics of the traffic on a road network [12]. In addition, an integrated optimization algorithm and coding adjustment strategy were made for the solution of multi-stage distribution and dispatching of disaster relief supplies. With the objective of minimizing the maximum distribution time, Zhou Yufeng et al. formulated an emergency facility siting-allocation model applicable to the initial post-earthquake relief phase by considering the phase characteristics, facility disruption scenarios, multi-species uncertain demand, facility capacity limitations, etc. The defuzzification of uncertain demand was processed through the expectation value formula of interval boundary, and on this basis, the result could be obtained by the proposed hybrid integer coding genetic algorithm [13].

The period vehicle routing problem (PVRP) was firstly proposed by Bekrami and Bodin in 1974 [14], arguing that different customers have different access frequencies for the recycling of industrial waste in New York City. Christofides et al. (1984) initially constructed a mathematical model of PVRP [15]. After nearly forty-five years of development, PVRP has been further extended in practical applications, such as the period vehicle routing problem with time window (PVRPTW), multidepot and periodic vehicle routing problems (MAPVRP) and the dynamic multi-period vehicle routing problem (DPVRPD) [16–19], and other existing studies mostly used heuristic algorithms to solve the extended PVRP model.

Wang et al. (2019) put forward a multi-stage model for distributing emergency supplies to multiple affected locations with the objectives of minimizing losses caused by shortages, total fixed costs of transportation and distribution costs. They designed and constructed a nonlinear utility function to reflect the negative utility caused by a lack of funding, and experiment results proved the feasibility of this model [20]. With the objectives of minimizing total delay time and total system loss for distribution of emergency supplies, Wang Yanyan et al. developed a dynamic distribution optimization decision model that uses fuzzy information conditions with multiple demand points, multiple distribution centers, multiple supplies, multiple periods and multiple objectives. After analyzing the clarification methods of the interval objective function, interval fuzzy constraints and triangular fuzzy constraints, they designed a two-dimensional Euclidean distance-based objective empowerment fuzzy algorithm to solve the model [21]. With the dual objectives of minimizing the risk of sending unsatisfying amounts of resources and minimizing the risk of resources not reaching disaster areas, Zhou et al. considered of the inherent nature of the multi-period dynamic emergency resource scheduling (ERS) problem to establish a multi-objective optimization model for the multi-period dynamic emergency resource scheduling (ERS) problem, and a decomposition-based multi-objective evolutionary algorithm (MOEA/D) was made to achieve great performance [22].

The vehicle routing problem for perishable goods (VRPFPG) is one of the vehicle routing problems (VRP) [23]. Large quantities of perishable goods around the world are transported from suppliers to consignees on a daily basis. Perishable goods, such as food and pharmaceuticals, require special handling during transportation due to their limited lifespans, and they must be transported as fast as possible before they deteriorate. Besides transport time constraints, the high frequency of transport can generate high transportation costs, which makes the optimization of perishable materials particularly vital. With the multiple objectives of minimizing operational costs, spoilage costs and carbon emissions, and maximizing customer satisfaction, Zulvia et al. paid attention to time windows, different travel times during peak and off-peak hours and working hours to develop a green VRP model and design a multi-objective gradient evolution (MOGE) algorithm whose results showed great performance [24].

To solve the perishable with uncertain demand material distribution vehicle path problem, researchers constructed various models with different objectives. With minimum total cost, maximum product freshness and minimum carbon emission as objectives, Qian Zhang et al. established a multi-objective optimization model for distribution path planning, and designed the main objective method and fruit fly algorithm based on robust optimization to deal with the uncertainty problem [25]. With the objective of minimizing the operating cost and emission cost, Babagolzadeh et al. constructed a two-stage stochastic planning model to determine the optimal replenishment strategy and transportation plan in the presence of carbon tax controls and uncertain demand, and an improved result was obtained by the proposed mathematical algorithm with respect to iterative local search (ILS) algorithm and mixed integer programming [26]. With the objectives of minimizing costs, minimizing environmental impacts and maximizing customer service levels, Talouki et al. formulated a dynamic green vehicle path model for perishable material under green transportation conditions in view of time window implementation, and then designed an algorithm based on a new augmented-constrained heuristic for solutions [27]. With the goal of profit maximization, Wu et al. developed a variable fractional inequality distribution path optimization model considering the uncertainty of perishable food demand for high speed rail and designed an augmented Lagrangian with the Euler algorithm based on the pairwise algorithm [28].

For the problem of uncertainty in demand and return of perishable goods with different periods, with the objective of minimizing the total cost of the system, Guo Jiangyan et al. constructed a multi-period closed-loop logistics network for perishable goods and figured out a mixed-integer linear programming (MILP) model solving by a proposed genetic algorithm [29].

For the problems of high-frequency distribution, uncertain demand and return of fresh goods due to perishability, with the objective of minimizing the total cost of the system, Yang et al. constructed the corresponding fuzzy mixed-integer linear programming (FMILP) model for the system and designed genetic algorithm (GA) and particle swarm optimization (PSO) algorithms to solve it [30].

The vehicle path problem is considered as an NP-hard problem, so it may be timeconsuming and ineffective to use ordinary mathematical methods, such as exact algorithms, to deal with it. Most scholars nowadays use intelligent optimization algorithms for solving such problems. The whale optimization algorithm (WOA) is a biomimetic metaheuristic algorithm developed by Mirjalili et al. in 2016 to simulate the feeding mode of whales [31]. In recent years, it has been successfully applied to some large-scale optimization problems with the advantages of few artificial parameters and simple operation, such as resource scheduling problems [32], workflow planning for construction sites, site selection and path planning [33] and neural network training [34]. However, because the traditional WOA has the disadvantages of slow late convergence and easily falling into a local optima, some scholars have combined other algorithms with it to improve its performance in operation speed. Rohit Salgotra et al. addressed the problems of poor search performance and easily falling into a local optimum of the WOA algorithm. Three different improved

versions, including WOA-adversarial-based learning, exponentially reduced parameters and worst-particle elimination and reinitialization methods, have been proposed to improve its exploratory capabilities. These properties have been exploited to improve the exploration capabilities of WOA by maintaining the diversity among search agents [35]. Shang Mang et al. proposed a WOA-based vehicle path optimization method for the distribution logistics of the VRP problem; modified the WOA algorithm using random inertia weights and a non-uniform variation strategy; and verified the effectiveness of the improved algorithm by testing functions. The verification results showed that the improved whale optimization algorithm can efficiently optimize the distribution path for vehicles and reduce the distribution cost of logistics [36].

As a novel algorithm, the WOA algorithm has attracted extensive attention from scholars in various fields since a basis has been built for the research, development and improvement of the algorithm, and application studies have been conducted regarding engineering, scheduling, optimization and site selection. Additionally, there is a richness in algorithm improvement. However, there are fewer applications in vehicle path research, so further development and utilization are needed.

A great deal of research has been carried out in the existing literature on the optimization of vehicle paths for the distribution of emergency and perishable materials. In addition to large demands for emergency supplies, such as communication equipment, quilts and tents, in the early stage of post-disaster relief, there also would be large demands for life-saving and living emergency supplies, such as medicines and foodstuffs. As for the perishable characteristics of these emergency supplies, along with the likelihood of severe damage to some roads, there is often uncertainty about the needs of the affected sites, making it difficult for these emergency perishable supplies to be delivered quickly and meet demand requirements at once. Therefore, in order to improve the optimization efficiency, this paper combines the differential evolutionary algorithm with the whale optimization algorithm to solve the vehicle path problem for the distribution of emergency perishable materials with dual objectives, which is rarely studied at present. Finally, further verification of the effectiveness of the improved whale optimization algorithm at solving the realistic vehicle path problem through examples shows convincing performance.

#### **3. Multi-Period Optimization Model and DE-WOA Algorithm**

The distribution vehicle path problem for emergency perishable materials has special characteristics compared to the same problem for general emergency materials, which negatively affect the solving process. The standard whale optimization algorithm greatly improves the operation efficiency of the algorithm because of the relatively simpler process and searching mechanism. Thus, it is suitable for solving the problem of emergency perishable material distribution optimization. This sub-section, while taking the uncertain demand situation into account, analyses the situation of adequate supplies in distribution centers and constructs a multi-period vehicle path distribution optimization model with the dual objectives of minimizing the cost penalty of distribution delay and total corruption during delivery, and minimizing the total amount of demand that is not met. The improved differential evolutionary whale algorithm is designed to solve the model by combining the features of the differential evolutionary algorithm with the standard whale optimization algorithm with strong global search capability.

#### *3.1. Description of the Problem*

Given a simple discrete undirected road traffic network *G* = (*V*, *E*), where *V* = (*v*0, *v*1, *v*2, ... , *vn*) is the set of nodes and *v*<sup>0</sup> denotes the distribution center for emergency perishable relief supplies, *v*1, *v*2, ... , *vn* denotes the affected point and *E* = - *e*(*vi*, *vj*)|*vi*, *vj* ∈ *V* is the set of edges. Assume that the supplies are sufficient. The demand for emergency perishable supplies at the affected point *vi* is represented by interval boundary *Di*(*i* = 1, 2, ... , *n*). The distribution center *v*<sup>0</sup> possesses a sufficient emergency perishable supply, and the total quantity is *S*. The spoilage rate of emergency perishable supplies during the distribution

process will linearly change along with the distribution time, and it is *θtj*. The demand for emergency perishables at each site can be hardly met at once due to uncertain information on demand and limited vehicle capacity.

*v*0: The distribution center now has *k* vehicles available with a maximum capacity of *R* for each; *dij* represents the distance between any two points; *vij*, depending on the road conditions, represents the actual speed of the vehicle on edge (*i*, *j*) during transportation, and *vij* represents the average speed of the vehicle so that the actual time for the vehicle to reach the disaster site *j* is *tj* = *dij*/*vij* and the ideal time is *tj* = *dij*/*vij*. *cj* is defined as the delay penalty, relying on the demand and the degree of damage at the disaster site. The distribution service will deliver emergency perishable materials to each disaster site and back to the distribution center until all the needs of the disaster sites are met. The most probable value weighting method is used to identify the uncertain demand, and the distribution route and the amount of each demand point in each period are decided with the dual objectives of minimizing the cost penalty of distribution delay and the total corruption during delivery, and minimizing the total amount of demand that is not met.

The model was constructed based on the following assumptions.


#### *3.2. Model Building and Notation Definition*

The symbols and parameters used in this model are defined in Table 1. Decision variables are identified in Table 2.

**Table 1.** The symbols and parameters used in the model.


**Table 2.** Decision variables.


Taking into account all the objectives and constraints, the model is developed.

$$\min \sum\_{i \in \mathcal{Z}} \sum\_{p \in P} (t\_j - t\_0) c\_{1j} \mathbf{x}\_{ijk}^p + c\_{2j} d\_{ipk} \theta t\_j, \text{ while } t\_j \le t\_0, \ t\_j - t\_0 = 0 \tag{1}$$

$$\min \left\{ \sum\_{p \in P} 1 - d\_{ipk} (1 - \widetilde{\theta} t\_{\widetilde{\jmath}}) / \widetilde{D}\_{ip\prime} \,\,\,\mathrm{i} \in Z, \,\,\, k \in K \right\} \tag{2}$$

$$\text{s.t. } \sum\_{i \in \mathcal{Z}} d\_{ipk} y\_{ik}^p \le R, \ k \in \mathcal{K}, \ p \in P; \tag{3}$$

$$\sum\_{i \in Z} \sum\_{k \in K} d\_{ipk} \le S\_\prime \ p \in P; \tag{4}$$

$$0 \le \sum\_{k \in K} d\_{ipk} \theta t\_j \le \vec{D}\_{ip}, \; i \in Z, \; p \in P; \tag{5}$$

$$\sum\_{i \in S} \sum\_{j \in S} x\_{ijk}^p \le |S| - 1, \ k \in \mathbb{K}, \ p \in P; \tag{6}$$

$$\sum\_{j \in \mathbb{Z}} \mathbf{x}\_{\text{ojk}}^p = \sum\_{j \in \mathbb{Z}} \mathbf{x}\_{\text{jok}}^p \le \mathbf{1}, \; i \in \mathbb{Z}, \; k \in \mathbb{K}, p \in \mathcal{P}; \tag{7}$$

$$\sum\_{j \in Z, i \neq j} \mathbf{x}\_{ijk}^p \le \mathbf{1}, \; i \in Z, \; k \in K, p \in P; \tag{8}$$

$$\mathbf{x}\_{ijk}^p \le y\_{ik'}^p \tag{9}$$

$$0 \le \frac{d\_{ipk}\theta t\_j}{d\_{ipk}} \le \delta;\tag{10}$$

$$\frac{d\_{ipk}}{R\_k} \ge \omega;\tag{11}$$

$$x\_{ijk}^p = 0 \text{ or } 1, \ y\_{ik}^p = 0 \text{ or } 1, \ (i, j) \in Z, \ p \in P; \tag{12}$$

$$d\_{ipk} \ge 0;\tag{13}$$

Equations (1) and (2), respectively, represent the dual objectives that minimize the cost penalty of distribution delay and the total corruption during delivery and minimize the total amount of demand that is not met. Equation (3) guarantees the load amount of each vehicle does not exceed the maximum capacity per vehicle. Equation (4) ensures that the total amount of distribution in each period is less than the available amount in the distribution center. Equation (5) is aimed at restricting the amount of emergency perishable supplies delivered to the disaster site in each period that does not exceed its ideal demand. Equation (6) indicates that the sub-loop in the distribution process is broken. Equation (7) guarantees each vehicle starts and ends transportation at the distribution center. Equation (8) presents a vehicle does not pass through the same path twice or more in any period. Equation (9) ensures that the vehicle serves a disaster site before passing through. Equation (10) indicates the degree of spoilage of emergency perishable materials during distribution should be less than a given rate. Equation (11) represents that each vehicle's utilization rate for each period should be more than a given rate. Equation (12) is related to the integer variable constraint. Equation(13) represents the non-negative constraint.

Where *<sup>D</sup>ip* = *<sup>D</sup>ip* + |*Di*(*p*−1) − <sup>∑</sup> *dik*(*p*−1)| when *<sup>p</sup>* ≥ 2.

#### *3.3. Clarity of Ambiguous Needs*

In this paper, uncertain demand is expressed as interval boundary:

$$D\_{ip} = [q\_{1i}, q\_{2i}, q\_{3i}]\_\prime \ q\_{1i} \stackrel{\ast}{\leq} q\_{2i} \stackrel{\ast}{\leq} q\_{3i} \tag{14}$$

The affiliation function is:

$$\mu\_{\bar{D}\_i}(\mathbf{x}) = \begin{cases} 0 & \mathbf{x} \le q\_{1i}, \mathbf{x} \ge q\_{3i} \\ (\mathbf{x} - q\_{1i}) / (q\_{2i} - q\_{1i}) & q\_{1i} < \mathbf{x} < q\_{2i} \\ (q\_{3i} - \mathbf{x}) / (q\_{3i} - q\_{2i}) & q\_{2i} < \mathbf{x} < q\_{3i} \end{cases} \tag{15}$$

where *q*1*i*, *q*2*i* and *q*3*i* represent the left boundary, the point with affiliation 1 (most likely value) and the right boundary of the interval boundary, respectively. The interval boundary is constant with the weights given by experts or decision-makers. *Dip* = [*q*1*i*, *q*2*i*, *q*3*i*] is expressed by the Equation (16)

$$
\tilde{D}\_{ip} = \omega\_1 q\_{1i} + \omega\_2 q\_{2i} + \omega\_3 q\_{3i} \,\text{.}\tag{16}
$$

*ω*<sup>1</sup> is the weight of the lower boundary, *ω*<sup>2</sup> is the weight of the most probable value and *ω*<sup>3</sup> is the weight of the upper boundary.

Such methods that determine weights by experience and knowledge of experts or decision makers are relatively subjective. The results thus are influenced by strong human elements. Some more objective methods to identify fuzzy weights were developed, such as the same weight method and hierarchical analysis. The most common is the most likely method. The most likely value of the interval boundary is given the highest weight, as it is most accurate. The value of boundary is less accurate; thus, they are assigned smaller weights.

To indicate differences between the three estimates *q*1,*<sup>q</sup>* and *q*3, the weights of them are consequently determined by *ω*<sup>1</sup> = *ω*<sup>3</sup> = 1/6 and *ω*<sup>2</sup> = 4/6. Therefore, Equation (15) can be converted into Equation (16).

$$
\hat{D}\_{ip} = \frac{q\_{1i} + 4q\_{2i} + q\_{3i}}{6} \tag{17}
$$

After replacing Equation (5) with Equation (17), the updated constraint is shown as Equation (18):

$$0 \le \sum\_{k \in K} d\_{ijk} \le \frac{1}{6} q\_{1i} + \frac{4}{6} q\_{2i} + \frac{1}{6} q\_{3i}, \ i \in Z, \ p \in P \tag{18}$$

$$\min \left\{ \sum\_{p \in P} 1 - d\_{ikp} (1 - \widetilde{\theta} t\_j) / \frac{1}{6} q\_{1i} + \frac{4}{6} q\_{2i} + \frac{1}{6} q\_{3i}, \ i \in Z, \ k \in K \right\} \tag{19}$$

#### *3.4. Handling of Dual Targets*

The conventional method aims to convert a muti-objective problem into a singleobjective optimization problem by linear weighting. However, because of the non-uniformity of the objective magnitude, the solution of the original problem and that of the converted problem are not in simple one-to-one correspondence. The weights of each objective may largely affect the accuracy of solutions. This paper takes advantage of the idea of the constraint method (Haimes et al. 1971), combining it with an improved differential evolutionary whale algorithm.

In this case, two single-objective optimization problems are solved by converting one of the dual objectives into the other's constraints based on the importance of the objectives in each period in turn and solving them separately to obtain the Pareto solution set of the model.

1. Construct a single objective optimization problem with objective A and objective B, respectively, and find the value domain (upper and lower bounds) of the two objective functions. Objective A.

$$\begin{aligned} \min & \sum\_{j \in Z} \sum\_{p \in P} (t\_j - t\_0) c\_{1j} x\_{ijk}^p + c\_{2j} d\_{ij \text{pk}} \theta t\_{j\text{\textquotedbl{}}} \text{ while } t\_j \le t\_0, \ t\_j - t\_0 = 0\\ \text{s.t.} & \operatorname{constant}(1) - \operatorname{constant}(13) \end{aligned} \tag{20}$$

Objective B.

$$\min \left\{ \sum\_{p \in P} 1 - d\_{ipk} (1 - \tilde{\theta} t\_j) / \tilde{D}\_{ip\prime} \, \, \, \mathbf{i} \in \mathbf{Z} \, \, k \in \mathbf{K} \right\} \tag{21}$$
  $\text{s.t.} \, \text{ constraint}(1) - \text{constant}(13)$ 

2. Step 1 finds the minimum value of objective A as *m*, and then adds *ZA* ≤ *a* as a constraint to get the result of objective B. Construct a single-objective optimization problem for objective B. If the problem has a feasible solution, find the optimal solution for objective B as *Z*∗ *<sup>B</sup>*, and go to Step 3; if there is no feasible solution, go to Step 4.

3. Then, add *ZB* ≤ *Z*<sup>∗</sup> *<sup>B</sup>* as a constraint to objective A to construct a single objective optimization problem for objective A. Similarly, if there is a feasible solution to the problem, the optimal solution to objective A is found at *Z*∗ *<sup>A</sup>*, at which point the solution obtained in the above step is counted in the Pareto solution set; if there is no feasible solution, then go to step 4.

4. Make *a* = *a* + . is a fixed step; go to step 2 to continue solving.

5. Stop when *a* is greater than the maximum value of target A.

*3.5. The Basic Process of DE-WOA*

The improved differential evolutionary whale algorithm is computed as Figure 1.

**Figure 1.** Flowchart of the DE-WOA algorithm.

Step 1: Uncertain demand clarification. The demand parameters in the model are the interval boundary. The most probable weight method is used to convert the interval boundary into definite values and replace them in the model.

Step 2: Initialize parameters. Assign values to parameters, such as population size *pop*, the maximum number of iterations *M*, the logarithmic spiral shape constant *b*, the scaling factor *F* and the crossover probability *CR*.

Step 3: Calculate the individual fitness function at *F* - and the population average fitness function at *F* - *average*. Based on the obtained fitness function values, record the location of the global optimal solution in the initial population *xbest*, where the global optimal value is *F* - *best*.

Step 4: When *F* - ≤ *F* - *average*, iteratively update the solution and calculate the values of parameters such as *a*, *A*, *p*, *C* and *l*; otherwise adapt it for global exploitation using *Xi*(*t* + 1) = *Xbest*(*t*) − *A* ∗ *D*, *A* = 2*a* ∗ *r* − *a*, *C* = 2 ∗ *r* to expand the population diversity.

Step 5: When *P* < 0.5 and |*A*| < 1, the whale position is updated using *D* = |*C* ∗ *Xbest*(*t*) − *xi*(*t*)|; when *P* ≥ 0.5 and |*A*| < 1, the whale position is updated using *D*- = |*Xbest*(*t*) − *xi*(*t*)|; when *P* < 0.5 and |*A*| ≥ 1, the whale position is updated using *D* = |*C* ∗ *Xrand*(*t*) − *xi*(*t*)|.

Step 6: Update the global optimal solution *xbest* and the global optimal value *F* - *best*.

Step 7: Stop the iteration if the algorithm stopping condition is met; otherwise, repeat step 4–step 7.

#### **4. Analysis of Numerical Examples and Computational Results**

In this section, we report the results of numerical experiments that were applied to verify the feasibility and effectiveness of the constructed model and proposed algorithm. All experiments were tested on a PC equipped with an Intel(R) Core(TM) i7-9750H CPU @ 2.60 GHz 2.59 and 8 GB of RAM. The model programming was solved by Python 3.8.1.

#### *4.1. Parameter Setting*

There are one distribution center and ten disaster locations labeled in order from 0 to 10. Related information, including coordinate values, is shown in Table 3. The network topology between the distribution center and the affected points is shown in Figure 2.


**Table 3.** Coordinates of the distribution center and locations of affected points.

**Figure 2.** Network topology in the affected area.

Due to the lack of information on data from the affected areas, the demand for emergency perishable goods and the speed of vehicle movements at each affected location need to be estimated based on published information, such as local casualties and the probability of secondary disasters. The specific demand parameters *q*<sup>1</sup> (pessimistic value), *q*<sup>2</sup> (most likely value) and *q*<sup>3</sup> (optimistic value) are shown in Table 4.

**Table 4.** Demand parameters.


The distribution center has three small trucks of the same type. In order to obtain a more accurate distribution time, the actual distance between any two points is measured according to the latitude and longitude of the map, and the maximum speed of the vehicles traveling on each road is estimated according to the road damage. The transport network parameters (*a*, *b*) and vehicle parameters are shown in Tables 5 and 6.

**Table 5.** Transport network parameters.




In Table 5, *a* denotes the distance between two points (km), *b* denotes the actual speed of the vehicle traveling between this path *vk* (km/h) and "-" denotes that this section is impassable, resulting in the delivery time parameters shown in Table 7.

**Table 7.** Distribution time parameters.


#### *4.2. Results*

After several trials, the parameters of the improved whale algorithm (DE-WOA) based on the difference algorithm were set as shown in Table 8.

**Table 8.** DE-WOA parameter settings.


After five trials, a Pareto frontier solution set for the problem was obtained and is shown in Figure 3. The horizontal coordinates and vertical coordinates, respectively, represent the value of objective A (the distribution delay penalty and corruption cost) and the value of objective B (total amount of demand that is not met). Each point represents a distribution solution that satisfies the Pareto optimum. The decision maker can choose a relative compromise by weighing the relationship between multiple objectives according to the situation in practice.

**Figure 3.** Pareto frontier solution set.

The relationship between the transport volume and the optimal solution at each affected point is obtained, and the optimal path is output as shown in Tables 9 and 10.

**Table 9.** Relationship between transport volumes at affected sites and optimal path.


**Table 10.** Distribution vehicle paths and objective function values.


#### *4.3. Algorithm Comparison*

Two algorithms, the standard whale algorithm (WOA) and the improved differential evolutionary whale algorithm (DE-WOA), were used to solve the algorithms, resulting in the set of Pareto front solutions under both algorithms shown in Figure 4. It can be seen from Figure 4 that the Pareto ranks of the solutions of the improved differential evolutionary whale algorithm are lower than the Pareto ranks of the solutions of the standard whale algorithm, thereby indicating that the improved differential evolutionary whale algorithm

can effectively improve the local search capability and increase the diversity of solutions in the population.

**Figure 4.** Set of Pareto front solutions under both algorithms.

After 100 runs of both algorithms, the optimal objective values were obtained as shown in Table 8, and the convergence of the algorithms under the two objectives was obtained as shown in Figures 5 and 6.

1. Comparison of target results.


2. Analysis of convergence effects.

**Figure 5.** Convergence diagram of distribution delay penalty and corruption costs.

**Figure 6.** Convergence diagram of total amount of demand that is not met.

By comparison, it can be seen that the improved differential evolutionary whale algorithm outperformed the standard whale algorithm in terms of solution results and converged faster when solving. This shows that the improved differential evolutionary whale algorithm outperforms the standard whale algorithm in solving the dual objective model of this paper, thereby also verifying the validity of the model and algorithm.

#### **5. Discussion and Conclusions**

For the optimization problem of the multi-cycle distribution of emergency perishable materials under a uncertain demand, this paper draws the following conclusions.


The main purpose of this paper was to provide a set of scientific distribution scheme for emergency rescue, through the reasonable distribution of emergency perishable materials and reasonable arrangement of vehicles, so as to effectively reduce the damage caused by an earthquake, reduce casualties, improve the rescue work efficiency, etc. The research model and algorithm proposed in this paper can be applied not only to the disaster scenario, but also to the logistics distribution in urban and rural areas in practical daily life, which can effectively improve the operation efficiency among supply chains.

The model we proposed in this paper also has a few shortcomings. This case operates under the assumption that the distribution center has sufficient supplies and only distributes a single variety of perishable materials, though in practice the distribution center is often short of supplies and the demand for emergency perishable materials at the disaster site is often multi-species. In future work, we will take this shortcoming into account and consider how to combine and distribute multiple species of emergency perishable materials and improve the model to take more factors into account and build a more realistic emergency material distribution model. At the same time, as the complexity of the model increases, more efficient algorithms should be designed to correspondingly solve the model.

**Author Contributions:** Conceptualization, X.L.; Data curation, X.L.; Formal analysis, X.L.; Funding acquisition, Y.X. (Yang Xu); Investigation, X.L.; Methodology, X.L.; Project administration, Y.X. (Yang Xu), K.K.L., H.J. and Y.X. (Yaning Xu); Resources, Y.X. (Yang Xu) and K.K.L.; Software, X.L. and H.J.; Supervision, Y.X. (Yang Xu), K.K.L., H.J. and Y.X. (Yaning Xu); Validation, Y.X. (Yang Xu) and K.K.L.; Visualization, K.K.L. and Jia Li; Writing—original draft, X.L.; Writing—review & editing, X.L. and Y.X. (Yang Xu). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by General Research Project on Major Theoretical and Practical Issues in Philosophy and Social Sciences of Shaanxi Province of funder grant number 2022ND0185. The APC was funded by Yang Xu.

**Data Availability Statement:** Not applicable

**Conflicts of Interest:** The authors declare no conflict of interest

#### **Abbreviations**


#### **References**


## *Article* **Multi-AGV Flexible Manufacturing Cell Scheduling Considering Charging**

**Jianxun Li 1, Wenjie Cheng 1, Kin Keung Lai 2,\* and Bhagwat Ram <sup>3</sup>**


**\*** Correspondence: mskklai@outlook.com

**Abstract:** Because of their flexibility, controllability and convenience, Automated Guided Vehicles (AGV) have gradually gained popularity in intelligent manufacturing because to their adaptability, controllability, and simplicity. We examine the relationship between AGV scheduling tasks, charging thresholds, and power consumption, in order to address the issue of how AGV charging affects the scheduling of flexible manufacturing units with multiple AGVs. Aiming to promote AGVs load balance and reduce AGV charging times while meeting customer demands, we establish a scheduling model with the objective of minimizing the maximum completion time based on process sequence limitations, processing time restrictions, and workpiece transportation constraints. In accordance with the model's characteristics, we code the machine, workpiece, and AGV independently, solve the model using a genetic algorithm, adjust the crossover mutation operator, and incorporate an elite retention strategy to the population initialization process to improve genetic diversity. Calculation examples are used to examine the marginal utility of the number of AGVs and electricity and validate the efficiency and viability of the scheduling model. The results show that the AVGs are effectively scheduled to complete transportation tasks and reduce the charging wait time. The multi-AGV flexible manufacturing cell scheduling can also help decision makers to seek AGVs load balance by simulation, reduce the charging times, and decrease the final completion time of manufacturing unit. In addition, AGV utilization can be maximized when the fleet size of AGV is 20%-40% of the number of workpieces.

**Keywords:** AGV scheduling; flexible manufacturing cell; AGV charging; genetic algorithm

**MSC:** 68T20

#### **1. Introduction**

An AGV is a transport vehicle that can navigate autonomously and carry out spontaneous or controlled transportation along a prescribed route [1]. It assists intelligent factories to realize workpiece production tasks under the condition of unmanned driving. In order to efficiently use AGVs to participate in intelligent operations, it is necessary to integrate information such as processes, equipment, and workpieces into workshop scheduling. This kind of scheduling is actually the Flexible Job-Shop Scheduling Problem (FJSP) that has evolved from the traditional Job-Shop Scheduling (JSP). There is an exceptionally rich domain of research that covers many aspects of AGV scheduling and proposes numerous approaches. The majority of interest in AGV scheduling stems from its potential impact in practice by increasing efficiency and lowering costs. However, due to its complexity and numerous features, it remains a challenging problem [2]. The literature review and discussion in this paper are limited to the number of available AGVs, the power threshold, and the charging time in AGV and machine task scheduling issues with single or multiple objectives AGVs. The scheduling optimization problem of AGVs is the primary focus of

**Citation:** Li, J.; Cheng, W.; Lai, K.K.; Ram, B. Multi-AGV Flexible Manufacturing Cell Scheduling Considering Charging. *Mathematics* **2022**, *10*, 3417. https://doi.org/ 10.3390/math10193417

Academic Editor: Yaping Ren

Received: 22 August 2022 Accepted: 17 September 2022 Published: 20 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

related research. Fu et al. [3] elaborated the analysis process, scheduling rules and optimization methods for the AGV scheduling optimization problem. Bilge and Ulusoy [4] proposed a pseudo-polynomial-time algorithm to obtain optimal machine and vehicle schedules. Abdelmaguid et al. [5] studied a hybrid genetic algorithm approach to schedule the machines and used a heuristic technique to obtain a vehicle assignment. The authors addressed the operations scheduling and provided the vehicle assigning heuristic with a starting time for each operation, as well as its predecessor operation in the job sequence. Zhang et al. [6] proposed a mixed integer linear programming model for the joint production and transportation scheduling in flexible manufacturing system (FMS) without considering the transportation tasks associated with each job from the lower and upper area to the machine. Fontes and Homayouni [7] established a new mixed integer linear programming model to solve the joint production transportation scheduling problem, which integrated the machine scheduling problem and the AGV scheduling problem, and used two sets of chain decision strategies (machine and AGV). Using the hybrid taboo bat algorithm, Yonglai et al. [8] carried out in-depth analysis on AGV material distribution scheduling and provided an effective solution for workshop material distribution scheduling by multiple AGVs under certain constraints. In addition, Liu et al. [9] also adopted the improved pollen algorithm to study the co-integration AGV manufacturing unit scheduling problem. Regarding the average delay in flexible manufacturing units with AGVs, Heger [10,11] studied the dynamic priority assignment rules of AGVs using sorting and routing rules. Xu and Guo [12] discussed the multi-objective and multi-dimensional green scheduling method of FMS with AGV by means of the evolutionary algorithm of segmented coding, and Zhang et al. [13] studied the AGV allocation problem for mixed-flow assembly lines considering the load capacity of AGVs. Umar et al. [14] considered the path conflict problem of AGV during transportation, and solved the problem with an improved genetic algorithm, but the algorithm has a long running time and low efficiency. Mousavi et al. [15] compared the genetic algorithm with the particle swarm algorithm that combined the multi-objective AGV scheduling in order to solve the proposed model in the flexible manufacturing system, which was different from the conventional intelligent algorithm Nouri et al. [16] proposed; a hybrid meta-heuristic algorithm based on the clustering multi-agent model, which simultaneously scheduled the machines and AGVs, treated AGVs as special machines for coding. Zhang et al. [17] proposed a two-stage solution approach and a particle sworn optimization method to solve an energy-efficient path planning model of a single-load AGV. The authors analyzed the AGV energy consumption characteristics by motion state and vehicle structure, where it was found that energy consumption was an independent optimization objective in AGV path planning.

The above research has effectively advanced the exploration of AGVs' participation in intelligent job scheduling. However, there are relatively few studies on the AGV charging problem in the workshop scheduling problem in the existing literature, and only a few experts and scholars have discussed the related issues. A typical example was given by Dehnavi et al. [18], who studied the optimization problem of AGV charging station location in manufacturing units through GAMES simulation. Zhengfeng and Yangyang [19] studied the job-shop workshop scheduling problem under the consideration of multiple charging AGVs. Wang et al. [20] studied an AGV scheduling method that integrated AGV energy consumption and workshop complexity. However, this type of research adopted the AGV to perform a task and then returned immediately, and did not perform the task continuously. This was not in line with the reality that the AGV could continue to perform transportation multiple times after being charged that directly affected the determination of the number and scale of AGVs in the workshop, their work efficiency, and its evaluation. Fazlollahtabar and Saidi-Mehrabad [21] examined literature related to different methodologies to optimize AGV systems for the problems of scheduling and routing at manufacturing, distribution, trans-shipment, and transportation systems. The authors categorized the methodologies into simulation studies, metaheuristic techniques, artificial intelligent, exact, and heuristics mathematical methods. Further, De Ryck et al. [22] overviewed a number of AGV-related

control algorithms and techniques that were employed in not only the early stages of AGVs, but also the algorithms and techniques used in the most recent AGV-systems, as well as the algorithms and techniques with high potential.

Extant literature shows that experts and scholars have fully realized the significant role played by intelligent robots in the manufacturing industry, and have carried out research on the scheduling problem of flexible manufacturing cells with AGVs. The model realizes the solution to the problem and optimizes the scheduling optimization objective. However, the charging ability and continuous working ability of AGVs are rarely addressed in the model in literature. In particular, most of the existing literature assumes that the charging process is completed immediately and there is no waiting. Moreover, the selection of AGVs to transport workpieces is based purely on the customer's demand for completion time without regard to AGV priority, resulting in load imbalance in AGVs. Therefore, the calculation result is quite different from the actual scheduling AGV result, and sometimes a part of the workpieces cannot be completed according to customer demands. The load imbalance even caused some AGVs not to undertake any transportation tasks, but some AGVs are always in the working state. Then, the core question of research is: how to schedule AGVs to balance the distribution of transportation tasks and minimize the maximum completion time, while meeting customer demands, and considering the charging process and waiting time of AGVs?

Aiming to seek the marginal utility and load balance of AGV with charging, this paper studies the allocation and handling process scheduling problems of AGVs on the basis of FJSP. The AGV comprehensively considers factors such as workpiece scheduling, transportation task scheduling, scheduling priority, and charging constraints, to build models, and we solve these through genetic algorithms and complete use of AGV capabilities to obtain optimal production efficiency. To achieve this objective, we describe the problem and its formulation in Section 2. Further, algorithm development and solutions are given in Section 3. Section 4 presents numerical experiments followed by conclusion and future remarks in Section 5.

#### **2. Mathematical Modelling**

#### *2.1. Problem Description*

Given the number of machines, workpieces, and customer demands, the multi-AGV flexible manufacturing cell scheduling model explores the optimal distribution of multiple AGVs with charging, so as to complete all the processing tasks and minimize the final completion time of the manufacturing unit. In the model, AGVs must transport workpiece to machine according to the corresponding handling procedure. The total time consumption includes not only AGV transportation time and machine processing time, but also AGV charging time and waiting time. The problem of charging multiple AGVs in FJSP can be described as: there are *n* workpieces to be machined in a machining cell, expressed as *J* = {*J*1, *J*2, *J*3,..., *Jn*}, each workpiece *Ji*(*Ji* ∈ *J*) has *p*(*p* ≥ 1) different machining operations. Each operation *Oij* can be machined by an optional set of machines, the processing time of it is different for different operations *Oi*(*Oi* ∈ *O*) of each workpiece. The processing time of the same process varies with the change of the machine, and the processing sequence of the workpiece has been specified in the process route file in advance. Workpiece *Ji* contains *p* + 1 handling procedures *Tij*, and the distance between each processing machine is measured by the AGV travel time. In the course of processing, *k* AGV handling robots *R* = {*R*1, *R*2, ···, *Rk*} take out the workpiece blanks from the loading station and transport the workpiece *J* to *m* different machines *M* = {*M*1, *M*2, ···, *Mm*} for processing until the last process of the workpiece is completed, and the AGV transports the processed workpiece to the unloading station. The AGV charging method is to stay and charge at the charging station. The charging time cannot be ignored; it affects the final scheduling decision. Therefore, the charging time is taken into account in the AGV's handling time. The issues that need to be considered are the selection of the processing machine where the process is located, the sequence of each process in the machine, and the

AGV allocation problem that considers charging. In order to build a multi-AGV flexible manufacturing cell scheduling model with charging, this paper introduces the symbol definitions in Table 1, and assumes the following conditions:

(1) The machine follows the principle of first-come, first-served. Each machine can process only one workpiece at a time, and each AGV can transport only workpiece for one operation at a time; (2) Each machine and workpiece have the same priority except for user requirements, and the processing and AGV transport process are uninterrupted; (3) There is a workpiece buffer area next to each machine. Ignoring the workpiece congestion, the buffer capacity is unlimited; (4) Ignore the time the AGV takes to load and unload workpieces in the warehouse or machine buffer; (5) Initially, the AGV position and the initial workpiece position are located in the loading station, and the workpieces that have completed all processes are placed in the unloading station; (6) If two adjacent processes of the same workpiece are in the same machine, the AGV handling resources will not be occupied; (7) The AGV handling speed remains unchanged during transportation, is not affected by the load, and does not consider path conflicts; (8) The difference between power consumption of the AGV with or without a load is not considered, the unit power consumption is fixed; and (9) After the AGV completes the handling of a process, if the power is less than the threshold, it returns to the charging pile for charging. When the power is greater than the threshold and meets the AGV to perform the handling task, it waits for the next handling task. The following table provides symbols that satisfy assumptions and model requirements.

**Table 1.** Symbol definitions.


There are multiple AGVs in the manufacturing unit, and they stop working when the power stored in their batteries is exhausted, they are then required to be charged again. Therefore, it is necessary to define the available time parameters of the AGVs to ensure that when the AGVs are dispatched for handling, the AGVs are in the available state, not the charging state. The driving path of the AGV can be divided into three sections: the no-load time to the task location; the waiting time when waiting for the task profile to complete the previous process; and the time for the loaded profile to go to the destination machine. The waiting time is sometimes 0. Because the calculation of AGV power consumption is based on the driving power consumption, it is necessary to define the three states of AGV driving. The sum of the above three parts is the total time *Rr*−*ij* required to complete the task. In the case when waiting time is 0, there is *Rr*−*ij* = *Rr*−*total*.

Due to the flexibility of selecting machines for processing in FJSP, the decision variable *PMl ij* needs to be set. It means assigning machine *Ml* to complete the machining work of *Oij*. Similarly, the handling task needs to assign one AGV from multiple AGVs to perform the task, so it is also necessary to define decision variables to specify the AGV to perform the work. The specific decision variable symbols and meanings are as follows:

> *∂<sup>M</sup> ij*−*ij*- = 1,*Oij* is machined on the same machine before *Oij*- 0,*else ∂R ij*−*ij*- = 1, *Tij* is being transported by the same AGV before *Tij*- 0,*else PMl ij* = 1, Operation *Oij* is processed on machine *Ml* 0,*else PRr ij* = 1, Conveying process *Tij* is transported by AGV(*Rr*) 0,*else*

#### *2.2. Model Formulations*

The multi-AGV flexible manufacturing cell scheduling considering charging problem is actually a typical Flexible Job-shop Scheduling Problem. Because of the continuity of the system, it is not suitable to be solved by Mixed Integer Linear Programming or Multistage Stochastic Programming. Considering the complexity of Deep Q-learning in the calculation of partial derivatives, we employed Multivariate Nonlinear Programming, so as to comprehensively describe system constraints and process status. The scheduling objective is to minimize the maximum time cost of production. The completion time of the last process of the workpiece is *tipe*, then the objective function that can be expressed as:

$$f = \min\_{1 \le i \le n} (\max(C\_i))\tag{1}$$

The flexible manufacturing cell scheduling constraints considering multi-AGV charging are:

$$\mathbb{C}\_{i} \ge t\_{i\text{ps}} + T\_{l-i\text{j}} + T\_{r-l-l\text{ll}},\\\mathbb{R}\_{r-i(j+1)\text{s}} \ge t\_{i\text{js}} + T\_{l-i\text{j}},\\t\_{i(j+1)\text{s}} \ge R\_{r-i\text{js}} + R\_{r-i\text{j}}\tag{2}$$

$$
\partial\_{ij - i'j'}^{M\_l} + \partial\_{i'j' - ij}^{M\_l} = 1,\\
\sum\_{M\_l \in (M\_1, M\_2, \dots, M\_m)} P\_{ij}^{M\_l} = 1 \tag{3}
$$

$$t\_{i\bar{j}s} + T\_{l-\bar{i}\bar{j}} \le t\_{i'\bar{j}'s} + (1 - \mathfrak{d}^{M\mathfrak{l}}\_{i\bar{j}-\bar{i}'\bar{j}'}) \times Q,\\ P^{M\mathfrak{l}}\_{i\bar{j}} \times P^{M\mathfrak{l}}\_{i'\bar{j}'} \times (t\_{i'\bar{j}'s} + T\_{l-\bar{i}'\bar{j}'}) \le t\_{i\bar{j}s} + \mathfrak{d}^{M\mathfrak{l}}\_{i\bar{j}-\bar{i}'\bar{j}'} \times Q \tag{4}$$

$$
\partial\_{i\underline{j}-\dot{i}'\underline{j}'}^{R\_r} + \partial\_{\dot{i}'\dot{j}'-\dot{i}\underline{j}}^{R\_r} = 1 \tag{5}
$$

$$R\_{r-i\bar{j}s} + R\_{r-i\bar{j}} \le R\_{r-\bar{i}'\bar{j}'s} + (1 - \partial\_{\bar{i}\bar{j}-\bar{i}'\bar{j}'}^{R\_{r-\bar{i}}}) \times Q\_{r} P\_{i\bar{j}}^{R\_{r}} \times P\_{i'\bar{j}'}^{R\_{r}} \times (R\_{r-\bar{i}'\bar{j}'s} + R\_{r-\bar{i}'\bar{j}'}) \le t\_{\bar{i}\bar{s}} + \partial\_{\bar{i}\bar{j}-\bar{i}'\bar{j}'}^{R\_{r}} \times Q \tag{6}$$

$$\sum\_{\substack{R\_r \in \left(R\_1, R\_{2^\*}, \dots, R\_k\right)}} P\_{ij}^{R\_r} = 1 \tag{7}$$

$$R\_{r-ijs} = \max\left\{ R\_{r-\text{ablle}} + T\_{r-\text{pos}-M\_{i(j-1)}}, t\_{i(j-1)c} \right\} \tag{8}$$

$$R\_{r-i\text{jc}} = R\_{r-i\text{js}} + R\_{r-i\text{j}} \tag{9}$$

$$R\_{r-i\circ j} = T\_{\text{agv}} \times \left( R\_{r-\text{ablc}} + T\_{r-\text{pos}-M\_{i(j-1)}} - t\_{i(j-1)\circ} \right) + T\_{r-O\_{ij}-\Delta} + T\_{r-M\_{i(j-1)}} - M\_{ij}$$

$$T\_{\text{agv}} = \begin{cases} -1\_r R\_{r-\text{ablc}} + T\_{r-\text{pos}-M\_{i(j-1)}} < t\_{i(j-1)\circ} \\ 0, \text{else} \end{cases} \tag{10}$$

$$R\_{r-total} = T\_{r-pos-O\_{ij}-\Lambda} + T\_{r-O\_{ij}-\Lambda} + T\_{r-\Lambda-M\_{ij}} \tag{11}$$

$$R\_{r-\text{actual}} = E\_{\text{max}} - (R\_{r-\text{total}} \times E\_{\text{average}}), \\ R\_{r-\text{ct}} = (E\_{\text{max}} - R\_{r-\text{actual}})/E\_v \tag{12}$$

$$\Theta\_{\rm I} = \Theta \times \frac{1}{\sum\_{\text{L}^{\prime} \text{th} \leq t} \mathcal{U}\_{\rm rh} + 1} \times \frac{E\_{\rm rt}}{E\_{\rm max} - E\_{\rm min}} \tag{13}$$

$$
\Omega\_r = \Omega \times T\_{r-O\_{ij}-\Lambda} \times \tau\_{j\epsilon} \times \frac{E\_{\text{max}} - E\_{\text{min}}}{E\_{rt}} \tag{14}
$$

$$\forall t\_{ijc} \stackrel{\ast}{\preceq} \ x\_{jc}, \forall t\_{ijs} \stackrel{\ast}{\preceq} R\_{r-ijc}, \forall \downarrow \\ I\_{rt} \stackrel{\ast}{\preceq} \downarrow \mathcal{U}\_{\max} \tag{15}$$

$$
\forall E\_{rt} \leq E\_{\text{max}}, \forall E\_{rt} \geq E\_{\text{min}} \tag{16}
$$

$$[\left(\forall R\_{r-i\bar{j}s} (r \in R\_r)\right) \cap \left(\neg \exists R\_{r'-i\bar{j}s} (r' \in R\_r)\right)] \cap [\neg \exists r \neq r' (T\_{r-O\_{i\bar{j}}-\Lambda} = T\_{r'-O\_{i\bar{j}}-\Lambda})] \tag{17}$$

$$i \in \{1, 2, \ldots, n\}, j \in \{1, 2, \ldots, p\}, r \in \{1, 2, \ldots, k\} \tag{18}$$

In the above inequalities and equations, constraint (2) represents the process sequence constraint of the same workpiece, which means that the maximum completion time of the workpiece is jointly constrained by the start processing time of the last process, the processing time, and the AGV handling time, and the AGV start handling time must be greater than the completion time of the previous process. Equation (3) represents the constraints of processing machine resources, requiring that one machine can process only one workpiece at a time, and the workpiece can be allocated to only one machine for processing. Here, *Q* is a positive integer that represents the elasticity of the system. Equation (4) indicates that the two machining processes on the same machine satisfy the sequence relationship. Equations (5)–(12) represent the AGV resource constraints, where Equation (5) indicates that an AGV can handle only one workpiece at a time, and Equation (6) indicates that the two handling processes completed by the same AGV satisfy the sequence relationship. Equation (7) means that a workpiece can be transported by only one AGV at a time. Equation (8) indicates that the start time of the AGV is constrained by the available time of the AGV, the time from the AGV to the processing equipment of the previous process, and the completion time of the previous process of the workpiece. Formula (9) indicates that the completion time of the AGV handling task is constrained by the time when the AGV starts to transport and the time from the current position to the processing machine. Equation (10) indicates that the actual transportation time is constrained by the available time of the AGV, the current position of the AGV, the start position of the transportation, and the completion time of the previous process. Equation (11) represents the constraint on the total transportation time of the AGV. Equation (12) represents the constraint on the actual power of the AGV and the charging time of the AGV. Equation (13) represents the priority of selecting the AGV, which is composed of two parts, one is the load size of the completed task, and the other is the current power situation. To achieve load balancing, the AGV with sufficient power and small historical load is preferentially selected to undertake transportation task. Θ is a constant for equivalent conversion. Equation (14) represents the priority of the AGV loading workpiece for *Oij*. The AGV with earlier completion time and

arrival time, and smaller power, has a higher priority for loading workpiece and a lower wait time. Here, Ω is a constant for equivalent conversion. Equation (15) indicates the time constraint and load constraint of any processing operation required by the customer. Equation (16) indicates the power constraint on which the AGV can continue to work. Equation (17) indicates the constraint to avoid deadlocks. The AGV must exclusively serve the operation *Oij* in order. At the same time, only one AGV can serve one operation, and only one AGV can pick up workpiece for an operation from loading station. Equation (18) are the space sizes of all entities in the system.

#### **3. Solution Algorithms**

The scheduling problem of flexible manufacturing units with AGVs considering charging can be divided into three sub-problems, namely, assigning processes to machines, sequencing processes assigned to each machine, and AGV allocation for transportation tasks. Ant Colony Algorithm, Gray Wolf Algorithm, and Deep Learning Network, is applied to obtain optimal strategy, but all of these received local solutions. However, Genetic Algorithm is more inclined to global search, which can quickly retrieve all feasible solutions and prevent the convergence of local optimal solutions too fast. Thus, aiming at the above three problems, this paper introduces the elite retention strategy and uses the genetic algorithm to solve it. First, all feasible solutions are divided into three segments of chromosomes: process, machine, and the AGV, then the population is obtained after initialization, and then decoding is performed after performing crossover and mutation. Finally, the fitness value is calculated to complete the selection until it reaches the threshold of the number of iterations. The core algorithm part is as follows.

#### *3.1. Coding and Initialization*

#### 3.1.1. Encoding Operation

The genetic algorithm involves three segments of chromosome coding: process coding, machine coding, and AGV coding. The process code describes the sequence of the production process, and its gene length is *<sup>L</sup>*<sup>1</sup> <sup>=</sup> *<sup>n</sup>* ∑ *i*=1 (*Pi* + 1). The machine code indicates the processing machine where the process of each workpiece is located, and the length of the gene string is *<sup>L</sup>*<sup>2</sup> <sup>=</sup> *<sup>n</sup>* ∑ *i*=1 *Pi*. The AGV code corresponds to the task code, and the length of the gene string is the length of the AGV task code string. For example, in the code shown in Figure 1, process codes are sorted from left to right, the first number 3 represents the first process of workpiece No. 3, and the fourth number 2 represents the workpiece No. 2. The fifth number 3 of the machine code indicates that the second process of the workpiece No. 2 is processed by machine No. 3, and the tenth number 4 indicates that the unloading task is virtualized by the virtual equipment, that is, the unloading station. The fourth number of the AGV code is 1. The fourth task representing the process code, that is, the second process of workpiece No. 2, is transported by the AGV numbered 1. If the handling process is assigned to the corresponding AGV code, a sequence with process handling will be generated: (*O*31, *M*3, *T*1−03),(*O*11, *M*1, *T*2−01),(*O*21, *M*2, *T*2−02), The corresponding handling time series is [1, 1, 3].

**Figure 1.** Individual coding gene string.

#### 3.1.2. Initialization Operation

During initialization, depending upon the requirement of minimizing the maximum completion time, the machine with the shortest processing time is selected with a probability of 0.2, and the AGV code is initialized according to the principle of uniform distribution to balance the usage efficiency of the AGV and avoid idleness. After obtaining the fitness according to the objective function, the roulette method is used to select individuals from the population of size *N*. The probability that individual *i* can be retained for the next generation is related to the fitness value *fi*, and the probability is *pi* = *fi*/∑*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *fi*. In addition, considering that after a series of genetic manipulations in the parent generation, individuals with high fitness have a certain probability of being eliminated, which may lead to a decrease in the average fitness of the group. To this end, we use an elite retention strategy to ensure the top 10% of the parent population do not participate in crossover and mutation. When there are multiple optimal solutions at the same time, the overall handling time of the AGV is used as the second evaluation criterion to calculate the total

running time *k* ∑ *i*=1 *Rr*−*total* of the AGV. The individuals are sorted in sequence according to formula: *<sup>F</sup>* <sup>=</sup> *<sup>F</sup>* <sup>+</sup> 0.0001 <sup>×</sup> <sup>∑</sup>*<sup>k</sup> <sup>i</sup>*=<sup>1</sup> *Rr*−*total*, and the individual with the smallest AGV total running time is taken as the optimal solution. After participating in the crossover mutation, a new species population of size P excepting the elite is generated, and then 10% of the worst solutions are replaced with the optimal solution stored by the parent population to maintain the next generation of the parent population and ensure excellent genetic continuity and diversity.

#### *3.2. Crossover and Mutation Operations*

#### 3.2.1. Crossover Operation

The interleaving operation of the process code uses the POX method. Select all processes for a certain workpiece from the parent P1, keep the positions of all processes of the workpiece unchanged on the chromosome, and keep them separate from the child individual C1. Randomly generate the parent P2, and encode the parent individual P2. (Except for the process selected in the previous step). Fill in the blank positions of C1 in turn. Similarly, the corresponding parts of P1 and P2 are combined in order to obtain the offspring individual C2. As shown in Figure 2a, all processes of workpiece 3 in the parent P1 are randomly selected and inherited to the offspring C1, and are kept in the position in the gene string remains unchanged. Then, all the processes of the parent P2, except for workpiece 3, are sequentially filled in the remaining gene string vacancies of the offspring C1 to form a complete offspring C1. Similarly, the remaining parts of P1 and P2 can be merged. The offspring C2 is obtained. The crossover selection of machine gene strings is operated using the Multi-point Preservative Crossover (MPX) method, which is used for machine crossover for chromosome process assignment, and the process order is reserved for the offspring, as shown in Figure 2b. For the parent gene string P1 And P2, the randomly generated 0\_1 sequence, exchange the genes located at positions 2, 5, 6, 8, and 12, of the chromosome, and cross to generate offspring C1 and C2, that is, 1 corresponds to the same position of the parent P1 and P2 and the rest of the positions remain unchanged. The crossover of the AGV gene string is consistent with the MPX crossover method of the machine gene string. As shown in Figure 2c, the genes at positions 3, 6, 8, 10, and 12 are exchanged.

**Figure 2.** Chromosome crossover operation and mutation operation. (**a**) POX crossover operation of process code; (**b**) MPX interleaving operation of machine code; (**c**) MPX crossover operation of AGV code; (**d**) Mutation operation of process code; (**e**) Mutation operation of machine code.

#### 3.2.2. Mutation Operation

For process code, the extended insertion mutation method is adopted, that is, a gene at a given position is randomly selected in the chromosome encoded by the process, and it is randomly inserted into another position in the chromosome under the constraint that the sequence of processes of the same workpiece is fixed. The sequence of the gene after the insertion is set back one space. As shown in Figure 2d, insert the second process of workpiece No. 3 at chromosome position 6 of parent P1 into the second process of workpiece No. 2 at chromosome position 4, the second process of workpiece 2 originally at position 4. The process needs to be moved one bit backwards. As shown in Figure 2e, the mutation operation of the machine code randomly selects one of the two machine codes, and selects the machine with the shortest processing time in the machine group. The gene at position 7 is replaced, that is, the processing machine of the first process of workpiece 3 is replaced, and the original processing machine No. 3 and processing time 5 are mutated into processing machine No. 2 and processing time 3. The AGV coding mutation operation is the same as the machine coding. For example, the AGV No. 1 at position 6 in the parent P1 chromosome is randomly selected according to the probability, and replaced by AGV No. 2, and the feasible solution progeny C1 for handling the AGV is obtained.

#### *3.3. Fit Function and Decoding*

The fitness function solved by the genetic algorithm for scheduling of multi-AGV flexible manufacturing units considering charging is *f* = *max*( 1≤*i*≤*n Ci*), the maximum completion time. The smaller the *f* is, the better is the individual. The decoding operation of chromosomes is actually the process of converting chromosomes into feasible solutions for

scheduling, which can obtain the process profile, the machine used in the process and its processing time, the AGV trolley transported, and its handling time, etc.

**Step 1:** Select a process chromosome from the process and machine code in turn to judge the task type, convert it into the corresponding process *Oij*; the machine *Ml* corresponding to the process; the corresponding processing time *Tl*−*ij* and the completion time *ti*(*j*−1)*<sup>e</sup>* of the previous process *Oi*(*j*−1) in the machine.

**Step 2:** Obtain the transportation task sequence *Tij*, and the transportation start point (L or machine), end point (U or machine) and transportation time.

**Step 3:** Read the AGV chromosome code in sequence, obtain the AGV status information, and assign the handling tasks to the corresponding AGVs in sequence.

**Step 4:** Calculate the time *Tr*−*pos*−*Mi*(*j*−1) from the AGV to the starting point of the handling task and the time *Rr*−*ijs* when the transportation starts, update time *Rr*−*ij* it takes for the AGV to handle the *j* process of Profile *i*, and calculate the total transportation time *Rr*−*total*, current actual power *Rr*−*actual*, and available time *Rr*−*able* for the AGV to perform each handling task.

**Step 5:** Based on the AGV status information, see whether the actual power of the AGV is less than the threshold; if it is less than the threshold, it will be charged at the loading station L, and calculate the charging time *Rr*−*ct*, and then update the available time *Rr*−*able* of the AGV.

#### **4. Analysis of Examples**

In order to verify the performance of the model and the algorithm, genetic algorithm and MATLAB programming are used to solve the FJSP 10×8 problem, FJSP 15×8 problem and FJSP 25 × 8 problem. We conduct hyper-parameter fine-tuning by GeatPy. Based on Values Exchanged Recombination and Mutation for Binary Chromosomes, the automated machine learning method is applied to obtain the optimal crossover and mutation operation probabilities. With the constructed field vector, objective vector, and fitness vector, we use the Elite Copy Selection method for fitness selection, and use the constraint violation value matrix to store the degree of individual violation constraints. After 1628 times evolution by Evolution Tracker, we acquire the parameters, as shown in Table 2. In addition, to reduce the problem scale, the maximum battery power is set to 100, the power threshold to 10, the number of charging piles to 1, the power consumption per unit time to 2, and the AGV charging rate to 6.

**Parameters Values** maximum battery power 100 power threshold 10 power consumption 2 AGV charging rate 6 population size 100 iteration threshold 10,000 mutation probability 0.208 crossover probability 0.846

**Table 2.** The parameters of the algorithm model.

#### *4.1. Analysis of Results*

According to the conversion time of 1:5, that is, one unit of time represents 5 min, we input the resource list and routing list according to the customer order. The final completion time of the manufacturing unit is optimally 227.5 min, which is about 3.79 h. The specific experimental results are shown in Figure 3. The black line in the figure is the driving path of AGV No. 1, and the red line is the driving path of AGV No. 2. The maximum completion time of the 10 × 8 scale case is a minimum of 45.5. The final completion time is measured by AGV removing all the workpieces. The production planning department can determine an appropriate time for different batches of production orders and draw actual schedules according to the scheduling plan. In the 10 × 8 calculation example, the two AGVs did not reach the power threshold, so the charging state wherein the AGV cannot be used does not appear in the Gantt chart. However, after the AGV completes the previous task and the handling time is greater than 90 units of time, it will fall below the power threshold and shall need to stop work for charging. In the case of 15 × 8 scale, the two AGVs reached the power threshold once and in the case of 20 × 8 scale, the two AGVs reached the power threshold twice. In both cases, the charging wait time is 0 and all transportation tasks are completed. However, in the 25 × 8 scale case, the charging wait time is 6 with all tasks completed. AGV1 reached the power threshold three times, and AGV2 reached the power threshold twice in this case. It means that the AGVs need to be recharged five times in total, which takes 438.06 units of time. Here, the current powers of AGVs when they go to charge are 11.25, 12.08, 12.67, 12.90, and 13.04, respectively. Although the current powers of the AGVs are higher than the minimum power *E*min, they have to be charged because they cannot complete a single transportation task. Obviously, the larger the number of transportation tasks for workpieces, the more AGV charging times, and the longer the waiting time. From the perspective of AGV load, AGV1 and AGV2 each completed 17 transportation tasks in the case 10 × 8 scale, that is, the load balance ratio reaches the optimal state of value of 0. The load balance ratio is defined as (Ψmax − Ψmin)/Ψ*avg*. Here, Ψmax, Ψmin and Ψ*avg* are the maximum, minimum, and average number of tasks completed by AGVs. As for the cases of 15 × 8 scale, 20 × 8 scale, and 25 × 8 scale, the load balance ratios are all less than 0.100, which are 0.074, 0.057, and 0.071, respectively. This shows that our model can effectively avoid AGV resource waste while meeting customer demands.

**Figure 3.** FJSP Gantt chart with 2 AGVs.

In general, the final completion time of manufacturing unit increases with the size of the problem. As shown in Table 3, when the number of AGVs is sufficient, the completion time grows in proportion to the problem scale. On average, for every 10 additional customer demands, the completion time increases by no more than 40 according to calculation examples. In contrast, if the number of AGVs is too small, the completion time will increase extremely with problem scale, and even some demands cannot be accomplished on time, as shown in the data marked with "\*" in Table 3. Only if the AGVs exceeds 4 can all the problems of different scales from 20 × 8 to 100 × 8 be solved. Obviously, increasing the number of AGVs will improve the ability to solve problems. As the amount of AGVs increases exponentially, the completion time first decreases rapidly and then becomes stable. For example, in the case of 80 × 8 scale problem, when the number of AGVs is doubled from 1 to 32 successively, the completion time is improved by 128.57, 85.46, 50.07, 24.45, and 21.49, respectively. This shows that blindly increasing AGVs cannot significantly improve the scheduling effect, but will cause part of the AGVs to be idle. It can be found that when the number of AGVs exceeds 1/5 of the scale of workpieces, the improvement space of completion time is very small. In this case, each additional AGV reduces the completion

time by a maximum of 3.17 unites. In order to obtain a more reasonable AGV fleet size, the marginal utility of the AGV is discussed in Section 4.2.


**Table 3.** Final completion time under different scales of problem and number of AGVs.

"\*" means some demands cannot be accomplished on time.

#### *4.2. AGV Marginal Utility*

The number of AGVs is closely related to transportation allocation, and it also has a greater impact on optimization goals. Usually, the AGV is subject to transportation constraints and a dedicated AGV is equipped for various workpieces, but it may increase the cost of the enterprise. In order to deeply study the marginal utility of AGVs, this paper selects 10 × 8 examples, takes different numbers of AGVs, and calculates the data 10 times. The results are shown in Table 4. When its number is about 20% of the number of workpieces, the optimization effect is prominent and appears as an inflection point. On the contrary, the optimization effect gradually decreases. The reason for the inflection point is that when the AGV input is less than 20% of the number of workpieces, the turnover between the AGV and workpiece is complicated, and the AGV needs more. It takes at least 15 units of time to charge during the next shutdown. It can be seen that AGV transportation and charging greatly affects production efficiency. When the number of AGVs is greater than 80% of the number of workpieces, although a better solution is achieved, the optimal performance remains unchanged and the economic cost increases. Enterprises can assign AGVs according to the economic situation of the workshop and planning of the completion time. As shown in Tables 3 and 4, when the number of AGVs is greater than 2/5 of the number of workpieces, the marginal utility by increasing the AGV is less than 2.5. It means that when the AGV is 40% of the workpieces, the optimal schedule is almost reached and it is uneconomical to increase the number of AGVs again. Thus, according to the results of different scales and considering the fixed cost of AGVs, it is recommended that the number of AGVs should account for 20–40% of the scale of the workpieces, and the marginal utility is the most reasonable.



#### *4.3. AGV Load and Charging*

In order to verify the optimization of our model on the charging process and load balance, we conducted experiments in the case of 100 × 8 scale. The rapid increase in the scale of the problem makes it difficult to complete all transportation tasks when AGVs are insufficient. As shown in Table 5, while the number of AGVs is 1 and 2, respectively, 4 customer demands and 3 customer demands are uncompleted on time. Until the AGV exceeds 4, all the workpiece transportation tasks can be effectively completed. Thereafter, as

the number of AGVs increases, the total charging times decreases almost exponentially and the waiting time drops sharply. When all transportation tasks are completed, the average wait time remains between 4 and 6. This indicates that the AVGs are effectively scheduled to accomplish all tasks, while maximizing the reduction of charging times and charging wait time. Especially when the AGVs exceeds 8, the charging times tend to be flat, and the utility brought by increasing AGV is not obvious. From the perspective of the AGV load, the load balance ratio is always less than 0.25, which means that the difference between the maximum and minimum load is less than 1/4 of the average load. This indicates that the loads of different AGVs are similar and tend to be balanced. In addition, the load balance ratio shows a bell-shaped trend of increasing first and then decreasing. The inflection point occurs when AGVs is 8, where the minimum load of AGVs changes by 20.7% and the maximum load of AGVs changes by only 17.3%. It means that when the charging times of AGVs do not decrease dramatically, the loads of AGVs will gradually become balanced based on scheduling priority. Furthermore, the rational load balance of AGVs is achieved when the charging times remains stable (the number of AGVs is greater than 12). That is, when considering charging, multi-AGV flexible manufacturing cell scheduling can assist decision makers to select a reasonable number of AGVs through simulation to meet the needs of all customers, reduce the charging times, and promote load balance.

**Table 5.** AGV load and charging.


#### **5. Conclusions**

AGVs are powered by batteries, they move along the planned path, and have automatic guidance equipment such as magnetic strips, rails, or lasers, which can assist the workshop to complete a series of transportation tasks, providing a strong support for intelligent manufacturing enterprises to reduce the time consumption of the production process. Considering the influence of the available number of AGVs, the power threshold, and the charging time in the workshop scheduling, we adopted minimization of the maximum completion time as the optimization objective, and have constructed a multi-AGV flexible manufacturing unit scheduling model considering charging. The genetic algorithm was used to solve the problem, and the optimized process arrangement, the AGV driving path, and the scheduling plan after the AGV reached the power threshold were obtained, which further confirmed the impact of AGV transportation and charging on production efficiency. Compared with the models given in the literature [18,19], our model more comprehensively considers the charging time and charging waiting. The AGV scheduling priority is also designed to drive AGV load balance. The AGV with earlier completion time, arrival time, and smaller power, has a higher priority for loading workpiece and a lower wait time. To achieve load balancing, the AGV with sufficient power and a small historical load is preferentially selected to undertake the transportation task. After analyzing the marginal utility of the AGV under different problem scales, the results shows that the AGVs are effectively scheduled to complete transportation tasks, while reducing the charging times and charging wait time. It is clear that AGV utilization can be maximized when the number of AGV scales is 20–40% of the number of workpieces. Furthermore, the scheduling model of multi-AGV flexible manufacturing cell when considering charging can help decision makers minimize the maximum completion time by simulation, and seek load balance, while meeting customer demands. The limitations of the model research in this paper are as follows: (1) considering the complexity of research caused by diverse workpiece weights, we ignored the time-consuming differences generated by AGV loads in order to simplify model building; and (2) the objective function focuses on the final completion time cost of manufacturing unit without considering the fixed cost of the AGV. Therefore, future research can try to establish a multi-objective nonlinear programming model including time cost and fixed cost, and discuss the influence of different loads on AGV power consumption.

**Author Contributions:** Conceptualization, J.L. and W.C.; methodology, W.C.; software, W.C.; validation, J.L. and W.C.; formal analysis, J.L. and W.C.; investigation, W.C.; resources, W.C.; data curation, W.C.; writing—original draft preparation, W.C.; writing—review and editing, J.L., K.K.L. and B.R.; visualization, W.C.; supervision, K.K.L. and B.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The program code and data that support the plots discussed within this paper are available from the corresponding author upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

