## 2 Problem Description and Formulations
### Problem Description  
In this problem, we define a directed graph \(G=(V,A)\), where \(V=\{0,\ldots,n+1\}\) is a vertex set, and \(A=\{(i,j)\in V,i\
eq j\}\) is an arc set. The vertex 0 and \(n+1\) represent the same central depot. For each route, the starting vertex is 0 and the ending vertex is \(n+1\). The vertex set \(V^{\prime}=\{1,\ldots,n\}\) represents demand locations. For convenience, we further define \(V_{0}=V^{\prime}\cup 0\) and \(V_{n+1}=V^{\prime}\cup n+1\). The arc \((i,j)\) represents the shortest route between vertices \(i\) and \(j\).  
The planning horizon in each cycle is \(T\). At the beginning of the planning cycle, check the initial inventory level \(I_{0}\). Then, a set of land vehicles and a set of aerial crafts depart from the central depot. They must return to the central depot before the end of the planning cycle. The unit transportation cost of a land vehicle is \(c^{l}\), and the unit transportation cost of an aerial craft is \(c^{a}\), where \(c^{a}>c^{l}\). Similarly, for each kind of vehicle, there are a capacity \(Q^{l}\), \(Q^{a}\) and a fixed cost \(f^{l}\), \(f^{a}\).  
We assume that the replenishment lead time of each demand location under each replenishment mode only considers the time consumption of transportation. We use \(l^{a}\) and \(l^{l}\) to represent the aerial and land replenishment lead time of each demand location, respectively. \(q^{a}\) and \(q^{l}\) represent the aerial and land replenishment quantity. Then we further define the lead time and replenishment quantity as 0 if the replenishment mode does not visit the demand location. For a demand location, the land replenishment mode may arrive earlier than the aerial mode due to the delivery order. Therefore, let \(l^{f}=min\left\{l^{a},l^{l}\right\}\) and \(l^{s}=max\left\{l^{a},l^{l}\right\}\) be defined as the lead time of the faster arrival replenishment method and  
Figure 1: A brief sketch for the problem (two planning cycles).  
the slower one when the demand location is replenished by both two modes. Similarly, \(q^{f}\) and \(q^{s}\) are faster and slower replenishment quantities, respectively. For a demand location, the inventory cost before the first replenishment arrival is not dependent on the replenishment quantity in the planning cycle, only depends on the initial inventory level, as shown in Figure 1. The initial inventory level is determined by the replenishment quantities of the last planning cycle. So, the objective should minimize the routing cost and the inventory cost from the beginning of this planning cycle to the fast replenishment arrival time in the next planning cycle. However, it is impossible to know the lead time information of the next planning cycle when planning the current cycle. Thus, we use the average fast lead time \(l_{p-1}^{qS}\) in the previous planning cycle to estimate it. This strategy is the planning horizon crossing.  
Due to incomplete information, the distribution of demand is unknown. We assume demands \(\hat{d}_{it}\) are independent, and their values are in the positive symmetric interval \([\hat{d}_{it}-\hat{d}_{it},\hat{d}_{it}+\hat{d}_{it}]\), where \(\hat{d}_{it}\) is the estimated value of the demand and \(\hat{d}_{it}\) the maximum deviation of the demand. At the time point \(t\), shelters deduct the demand between \(t-1\) to \(t\) from the inventory held. If the available inventory does not meet the demand, it will have a backlog of demand to be met in the future, and a shortage cost will occur; otherwise, the on-holding inventory will be stored in the demand location to meet future demand, and lead to a holding cost. Each demand location has a unit holding cost \(h\) and a unit shortage cost \(s\), where \(s>h\).