## 3 Model Development
### Mathematical Model  
Maximize  
\[\sum_{r=1}^{R}z^{r}\quad\text{ \ where }z^{r}\leq s_{f}^{kr}\quad\forall r, \forall f,\forall k \tag{1}\]  
Subject to  
\[s_{f}^{kr}\leq\beta_{1,f}^{k,r}-\beta_{2,f}^{k,r}.t^{r}+\left(1-\sum_{n=1}^{N }\sum_{m=1}^{M}\sum_{d=1}^{D}y_{n,m}^{k,r}.a_{f,d}^{r}\right).B\quad\forall r,\forall f,\forall k \tag{2}\]  
Equation (2) calculates the customer satisfaction according to the elapsed time to deliver when an order \(r\) consisting of a type of goods \(f\) is processed to a means of delivery \(k\). In this constraint, the value of customer satisfaction is greater than 1 by multiplying a very large number \(B\) when an order \(r\) consisting of type of goods \(f\) is not processed with a means of delivery \(k\). Therefore, as shown in Equation (1), the smallest value among the customer satisfaction values of order \(r\) becomes the actual customer satisfaction of order \(r\). \[\sum_{k=1}^{K}\sum_{n=1}^{N}\sum_{m=1}^{M}y_{n,m}^{k,r}=1\quad\forall r \tag{3}\]  
\[\sum_{r=1}^{R}y_{n,m}^{k,r}\leq 1\quad\forall k,\forall n,\forall m \tag{4}\]  
\[\sum_{r=1}^{R}y_{n,m}^{k,r}\geq\sum_{r=1}^{R}y_{n,m+1}^{k,r^{\prime}}\quad \forall k,\forall n,\ m=1,\ldots,M-1 \tag{5}\]  
\[\sum_{m=1}^{M}\sum_{r=1}^{R}y_{n,m}^{k,r}\geq\sum_{m^{\prime}=1}^{M}\sum_{r=1 }^{R}y_{n+1,m^{\prime}}^{k,r^{\prime}}\quad\forall k,\ n=1,\ldots,N-1 \tag{6}\]  
Equation (3) indicates that each order \(r\) can be assigned to only one delivery sequence \(m\) of a specific round \(n\) of a specific means of delivery \(k\), while Equation (4) also represents that only one order can be assigned to a specific delivery sequence \(m\) in a specific round \(n\) of means of delivery \(k\). In addition, as shown in Equation (5), in the \(n\)th round of means of delivery \(k\), orders can be assigned to the \((m+1)\)th delivery sequence only when orders are assigned at the \(m\)th delivery sequence. In the similar manner, Equation (6) is restricted that the \((n+1)\)th round can be created only when the \(n\)th round is made by order allocation. \[\sum_{m=1}^{M}\sum_{r=1}^{R}g^{r}.y_{n,m}^{k,r}\leq G^{k}\quad\forall k,\ \forall n \tag{7}\]\[\sum\limits_{m=1}^{M}\sum\limits_{r=1}^{R}v^{r}\cdot y_{n,m}^{k,r}\leq V^{k} \quad\forall k,\ \forall n \tag{8}\]  
Equations (7) and (8) represent the weight and volume constraints that means of delivery \(k\) can carry at once in specific round \(n\), respectively. In this study, it is assumed that UGVs and UAVs can visit at the depot multiple times to newly load the goods within period of time. Therefore, the capacity constraints are presented according to the weight and volume. \[0\leq t^{r}\leq T\quad\forall r \tag{9}\]  
\[\sum\limits_{d=1}^{D}\sum\limits_{f=1}^{F}y_{1,1}^{k,r}a_{f,d}^{r}\cdot t_{0,d }^{k}\leq t^{r}\quad\forall k,\ \forall r \tag{10}\]  
The elapsed time taken to process the order \(r\) should be greater than 0 and less than the maximum period of time \(T\), as shown in Equation (9). In addition, Equation (10) indicates that if order \(r\) is assigned to the first delivery sequence of the first round at the means of delivery \(k\), then the delivery arrival time of order \(r\) is equal or greater to the elapsed time for moving from depot to destination node \(d\). \[ \begin{array}{l}t^{\prime}+t_{s}^{k}+\sum\limits_{d=1}^{D}\sum\limits_{d=1} ^{D}\sum\limits_{f=1}^{F}\sum\limits_{f=1}^{F}y_{n,m-1}^{k,r^{\prime}}a_{f,d} ^{r}\cdot t_{d,d}^{k}\cdot a_{f,d}^{r}\\ -B\cdot\left(2-y_{n,m}^{k,r}-y_{n,m-1}^{k,r^{\prime}}\right)\leq t^{r}\quad \forall k,\ \forall r^{\prime},\forall r,\forall n,m=2,\ldots,\end{array}  \tag{11}\]  
\[ \begin{array}{l}t^{\prime}+t_{s}^{k}+\sum\limits_{d=1}^{D}\sum\limits_{d=1} ^{D}\sum\limits_{f=1}^{F}\sum\limits_{f=1}^{F}y_{n-1,m}^{k,r^{\prime}}a_{f,d} ^{r^{\prime}}\cdot\left(t_{d,0}^{k}+t_{s}^{k}+t_{0,d}^{k}\right)\cdot a_{f,d} ^{r}\\ -B\cdot\left(2-y_{n,1}^{k,r}-y_{n-1,m}^{k,r^{\prime}}\right)\leq t^{r}\quad \forall k,\ \forall r^{\prime},\forall r,\forall m,n=2,\ldots,N\end{array}