## 3 Materials and Methods
### Solving Algorithm
#### 3.3.1 S-PSO Algorithm  
PSO is a population-based stochastic optimization algorithm proposed by [PERSON] and [PERSON] in 1995 [44]. The algorithm was first inspired by the social interaction behavior of a bird community and fish education. In order to provide the best solution, each so-called \"particle\" updates its flight speed and current position iteratively according to its own flight experience and the flight experience of other particles. The speed and position update methods of traditional particle swarms are defined in Formulas (18) and (19). \[v_{i,i+1}^{d}=v_{i,t}^{d}+c_{1}*rand_{1}*\left(pbest_{i,t}^{d}-x_{i,t}^{d} \right)+c_{2}*rand_{2}*\left(\left(gbest_{i,t}^{d}-x_{i,t}^{d}\right)-x_{i,t}^{ d}\right) \tag{18}\]  
\[x_{i,i+1}^{d}=x_{i,t}^{d}+v_{i,i+1}^{d} \tag{19}\]\[x_{k}^{h}=\left[x_{k}^{h_{1}},\cdots,x_{k}^{h_{m}}\right] \tag{20}\]  
The method used in this paper is the set-based particle swarm optimization (S-PSO), which is a discrete particle swarm optimization algorithm designed specifically for solving the MDVRPTW. \(x\) represents the position of the particle, and \(u\) represents the velocity is represented by \(u\). If the number of particles in the swarm \(P\) is \(N\), \(h\) represents each particle in the swarm,\(h\in[1,n]\). \(k\) represents the current number of iterations. The position of the particle in VRPTW is expressed as \(x_{k}^{h}\), and the position of each particle in the m-dimension is expressed as \(x_{k}^{h_{j}},\delta\in(1,m)\). Therefore, the position of the particles can be expressed as Equation (20). The position of each particle represents a directed Hamilton circle in the complete graph, and the Hamilton circle can be converted into a specific solution for a set of vehicle routes using the constraint-based decoder in this paper. As shown in Figure 4, starting from the outer arc of the warehouse, if an arc in the Hamilton circle satisfies all the constraints in the model, it will be retained; otherwise, the warehouse will be inserted between two nodes and the arc replaced by two new arcs. Through the definition of particle position and a constraint-based decoder, each particle represents a feasible solution to the model of the collaborative and persistent delivery problem of SR in the algorithm. The position of the particle in each dimension is composed of two arcs adjacent to node d\(\delta\), as show in Equation (21), represents the current dimension; \(nb_{1}\) is the previous customer point; and \(,nb_{2}\) is the next customer point. among them, \(nb_{1},nb_{2}\in\{1,2,\ldots,\delta-1,\delta+1,\ldots,m\}\), \(nb_{1}\
eq nb_{2}\). \[x_{k}^{h_{j}}=\left[\langle nb_{1},\delta\rangle,\langle\delta,nb_{2}\rangle\right] \tag{21}\]  
The velocity of a particle is the probability of a set of arc lengths being selected. \(u_{k}^{h_{j}},\delta\in[1,m]\) represents the speed of the particle in the m-dimensional space. The speed in the m-dimensional space is defined in Equation (22). \[u_{k}^{h}=\left[u_{k}^{h_{1}},\cdots,u_{k}^{h_{m}}\right] \tag{22}\]  
\[u_{k}^{h_{j}}=\left\{\left\langle c_{i},c_{j}\right\rangle/p\left(c_{i},c_{j} \right)\right|\left\langle c_{i},c_{j}\right\rangle\in A^{\delta}\right\} \tag{23}\]  
\(P(c_{i},c_{j})\) represents the probability of edge \(\left\langle\left(c_{i},c_{j}\right)\right\rangle\) being selected for particle position update, \(P(c_{i},c_{j})\in[0.1]\). If \(P(c_{i},c_{j})=0\), it means that side \(\left\langle c_{i},c_{j}\right\rangle\) will not be selected. When  
Figure 4: Decoder method based on model constraints. the particle position is updated, the velocity u is converted into an edge, and the edge selection is shown in Equation (24). \[Cut\big{(}u_{k}^{h}\big{)}=\Big{\{}\left\langle c_{i},c_{j}\right\rangle\Big{|} \big{(}c_{i},c_{j}\big{)}/p\big{(}c_{i},c_{j}\big{)}\in u_{k}^{h}\Big{\}} \tag{24}\]  
The update of the particle position \(x_{k+1}^{h}\) adopts a constructed method; \(\Omega\) is used to represent the vehicle load and service time window.