**Notation of Variables:**

*Yijk*: loaded quantity of vehicle *k* from pickup trip *i* to *j*;

*Zijk*: unloaded quantity of vehicle *k* from delivery trip *i* to *j*;

*tcijk*: the transportation cost of vehicle *k* from customer *i* to *j*;

*etijk*: time for vehicle *k* to move from *i* to *j*;

*δik*: service time required by vehicle *k* to load/unload the quantity demanded at *i*;

*m*: number of vehicles;

*n*: number of customers;

*ck*: fixed cost of vehicle *k*;

*qk*: maximum capacity for each vehicle *k*;

*T*: planning horizon;

*P*: set of unit demand from each pickup stop;

*D*: set of unit demand from each delivery stop;


**Objective Function:**

$$\text{Minimize } Z = \sum\_{i=0}^{n} \sum\_{j=0}^{n} \sum\_{k=1}^{m} t c\_{ij} X\_{ijk} + \sum\_{k=1}^{m} \sum\_{j=1}^{n} c\_k X\_{0jk} \tag{1}$$

**Subject to:**

$$\mathop{\rm s\sum}\_{i=0}^{n}\sum\_{k=1}^{m}X\_{ijk} = 1,\text{ for } j=1,\text{ }2,\text{ }\dots,n;\tag{2}$$

$$\sum\_{j=0}^{n} \sum\_{k=1}^{m} X\_{ijk} = 1, \text{ for } i = 1, \ 2, \ \dots, n; \tag{3}$$

$$\sum\_{i=1}^{n} X\_{ihk} = \sum\_{j=1}^{n} X\_{hjk}, \text{ for } k = 1, 2, \dots, m; \ h = 1, 2, \dots, \dots, n; \tag{4}$$

$$\sum\_{j=1}^{n} X\_{0jk} \le \mathbf{1}, \text{ for } k = 1, \mathbf{2}, \dots, m; \tag{5}$$

$$\sum\_{i=1}^{n} X\_{i0k} \le 1,\text{ for } k = 1,\text{ }2,\dots,m;\tag{6}$$

$$X\_{ijk} + Z\_{ijk} \le \sum\_{i=1}^{n} Q\_i \sum\_{j=1}^{n} X\_{ijk'} \text{ for } k = 1, \ 2, \dots, \dots, m;\tag{7}$$

$$Y\_{jik} - Y\_{ijk} = \begin{cases} P\_{j\prime} & \text{if } j \in P\_{\prime} i = 1, \, 2, \, \dots, \dots, n, \\ 0, & \text{if } j \in D\_{\prime} i = 1, \, 2, \, \dots, \dots, n, \\ -\sum\_{i=1}^{n} P\_{i\prime} & \text{if } j \in 0, i = 1, \, 2, \dots, \dots, n; \end{cases} \tag{8}$$

$$Z\_{ijk} - Z\_{jik} = \begin{cases} 0, & \text{if } j \in P\_\prime i = 1, \, 2, \, \dots, \dots, n, \\ d\_{i\prime} & \text{if } j \in D\_\prime i = 1, \, 2, \, \dots, \dots, n, \\ \sum\_{i=1}^n d\_{i\prime} & \text{if } j \in 0, i = 1, \, 2, \dots, \dots, n; \end{cases} \tag{9}$$

$$\sum\_{i=0}^{n} \sum\_{j=0}^{n} \delta\_{ik} X\_{ijk} + \sum\_{i=0}^{n} \sum\_{j=0}^{n} \epsilon t\_{ijk} X\_{ijk} \le T, \text{ for } k = 1, 2, \dots, m;\tag{10}$$

$$DT\_{jk} = \left(\varepsilon t\_{\text{ij}} + DT\_{\text{ik}} + \delta\_{\text{j}}\right) X\_{\text{ijk}}, \text{ for } k = 1, 2, \dots, \dots, m;\tag{11}$$

$$AT\_k = (\varepsilon t\_{i0} + DT\_{i\bar{k}}) X\_{i0\bar{k}\prime} \text{ for } k = 1, \ 2, \dots, \dots, m;\tag{12}$$

$$AT\_d = AT\_{b\prime} \text{ for } a \neq b. \tag{13}$$

Equation (1) is the objective function whose goal is to minimize both transportation costs and fixed costs. The constraints of each customer being serviced by only one vehicle are described in Equations (2) and (3), while Equation (4) represents that each vehicle arriving at the node must also leave from the same node. Equations (5) and (6) state the constraints of each vehicle being only permitted to start from and return to the crossdocking depot and being used to serve at most one route, respectively. The loaded and uploaded demand from pickup and delivery processes cannot exceed the vehicle quantity limit detailed by Equation (7), while Equations (8) and (9) each represent the quantity limit for pickup and delivery en route. The total distance visited and time traveled cannot exceed the planning horizon, which is limited by Equation (10). The departure and arrival times are functioned by Equations (11) and (12), respectively. The vehicles' simultaneous arrival at the CD hub is stated by Equation (13).

#### **4. Proposed sAIS Algorithm**

#### *4.1. Clustering Method*

In order to minimize total cost, each dispatched vehicle must serve the maximum number of customers within its capacity limit while traveling the shortest distance possible. In this research, we simulate this practical strategy to solve the problem of generating vehicle route sequences with one CD operation in the supply chain. Moreover, for faster convergence, we use a two-phase algorithm to solve VRPCD.

In the initial route-generating phase, we use the sweep method as the first phase to provide the initial solution for the second phase of sAIS. Usually, the sweep method is applied to a polar coordinate system, and the center of the coordinate system serves as the depot of the VRP. In this case, for each traveling distance of route *Oi* (*i* = 1, 2, . . . , *n*, where *n* denotes the number of customers), the *i*th customer's *Ci* served as the first customer in the cluster. Next, search for the closest customer from *Ci* + 1 to *Ci* by increasing the angle, and add it to the same cluster without exceeding vehicle capacity. When a customer cannot be added to the current cluster due to a vehicle capacity limit, this customer becomes the first customer in the next cluster. After all customers are clustered, we calculate the objective value of this route. After *O*1, *O*2,..., *On*, were calculated, we chose the minimum distance of *n* routes as the initial solution of the AIS algorithm. Figure 4 illustrates an example of the clustering process for 12 nodes.

**Figure 4.** An illustration of clustering by the sweep method.

The total cost is accumulated by the fixed cost of vehicles used (operational cost) and transportation costs as vehicles travel from one customer to another. In some cases, the aggregated service time consumed by each customer is also considered. Overall, the shorter the distance traveled by vehicles, the lower the total cost incurred. Figure 5 exemplifies the hierarchies of the proposed two-phase algorithm for 10 customers with homogeneous demand served by 2 vehicles as well as the chromosome encoding format in this study.

**Figure 5.** An illustration of 10 nodes with homogeneous demand served by 2 vehicles.

#### *4.2. AIS Procedure*

As mentioned previously, the AIS was successfully implemented to solve combinatorial optimization and abnormal detection problems. For data representation, we used chromosomes to denote routing sequences. All suppliers and retailers were assigned a unique number. We made the variation field of each cell, and its dimensions (length of chromosome) equivalent to the suppliers' (for pickup routes) or customers' (for delivery routes) numbers to be served, as shown in Figure 5. Therefore, the cells' chromosomes represent the nodes of the positions to be visited. In the meantime, the suppliers' or customers' sequences also mean the routing of orders. First, somatic hypermutation will be selected to make cells evolve, which can search for any feasible cells in the solution field. Next, we calculate the affinities between the antigen-antibodies and antibodies-antibodies according to the objective function. After iterations, the relationship between antigens and antibodies has a greater affinity. Moreover, elimination is a mechanism that will delete worse cells with each iteration. With affinity maturation and elimination, the cells will mutate toward the better-quality field of solution. In this case, each iteration individually will preserve the best cell, which means a higher affinity value.

In the proposed sAIS optimization algorithm, the principles of each antibody mutation have four parts. The basic components of the sAIS applied in the VRPCD are shown as follows:


The proposed sAIS optimization approach has amechanismwith theinside operation processes:

1. Somatic Recombination

Somatic recombination is one of the gene rearrangement processes that involves cutting out small regions of DNA and then putting the remaining pieces of DNA back together in an error-prone way in the adaptive immune system. For every iteration of our proposed algorithm, we randomly create a new set of antibodies by using the somatic recombination process in order to search for optimized solutions;

2. Somatic Hypermutation

There are five types of immune mutations that we use in our proposed sAIS algorithm. Each type has a different function. If the cost of objection function of antibodies in IgM is larger than the cost of the old ones, we randomly selected one of four kinds of immune mutations to operate in the next step of mutation in the same iteration. In addition to the traditional four types of immune mutations, we proposed a new immune mutation, called IgG2, to further improve the performance of the sAIS optimizations:


**Figure 6.** The operation of IgM mode.

**Figure 7.** The operation of IgA mode.

**Figure 8.** The operation of IgG2 mode.

#### 3. Affinity Maturation

Affinity means the total cost of the vehicles routed from the objection function of the VRPCD model. Some mutant antibodies with lower costs may be generated after several iterations of somatic hypermutation. Then, from iteration to iteration, the cost of producing antibodies will steadily decrease. This phenomenon is called affinity maturation;

## 4. Elimination Process

The antibodies will be stored in a limited list per iteration. Except for the best antibody, the antibodies that were produced and mutated per iteration were deleted from previous iterations. As a result, antibodies will search for alternative solution spaces for the VRPCD;

5. Escape Criterion

We designed an escape criterion for searching the antibody in the global optimization to prevent the search from entering the same solution space and remaining in the local optimum. Figure 9 and Algorithm 1 show the proposed procedure of the sAIS algorithm.


**Figure 9.** The processes of the proposed sAIS flowchart.

#### **5. Computational Experiments**

First, a preliminary investigation of the optimal parameter settings for the proposed algorithm is presented here. Next, 60 VRPPD benchmark problems were retrieved from the VRP Web (http://www.bernabe.dorronsoro.es/vrp/ accessed on 3 February 2023) as the criteria to evaluate the performance of the proposed heuristics. The comparison set is based on the GA algorithm, which is functioned by the customized software Palisade Evolver (https://www.palisade.com/evolver/) industrial edition 5.5.1, the Genetic Algorithm Solver plug-in for Microsoft Excel 2016.
