t_{j}^{+}=t_{j}^{-}(1-z_{j}),\quad\forall j\in\mathcal{H}\setminus\\{0,n+1\\}$
(48) $\displaystyle w_{ij}\leq x_{ij},\quad\forall(i,j)\in\mathcal{A}$ (49)
$\displaystyle w_{ij}\leq y_{ij},\quad\forall(i,j)\in\mathcal{A}$ (50)
$\displaystyle w_{ij}\leq z_{i},\quad\forall(i,j)\in\mathcal{A}$ (51)
$\displaystyle w_{ij}\leq z_{j},\quad\forall(i,j)\in\mathcal{A}$ (52)
$\displaystyle x_{ij}+y_{ij}+z_{i}+z_{j}\leq
3+w_{ij},\quad\forall(i,j)\in\mathcal{A}$ (53) $\displaystyle t_{j}^{-}\leq
t_{i}^{+}(1-z_{i})+\tau_{ij}^{D}(1-w_{ij})+M_{2}(1-x_{ij}),\quad\forall(i,j)\in\mathcal{A}$
(54) $\displaystyle t_{j}^{-}\geq
t_{i}^{+}(1-z_{i})+\tau_{ij}^{D}(1-w_{ij})-M_{2}(1-x_{ij}),\quad\forall(i,j)\in\mathcal{A}$
(55) $\displaystyle l_{j}\leq
T_{\text{bl}}z_{j}+M_{2}\sum_{i:(i,j)\in\mathcal{A}}w_{ij},\quad\forall
j\in\mathcal{H}$ (56) $\displaystyle l_{j}\geq
t_{j}^{-}-M_{2}(1-z_{j}),\quad\forall j\in\mathcal{H}$ (57) $\displaystyle
l_{j}\geq\sum_{i:(i,j)\in\mathcal{A}}\Big{(}\tau_{ij}^{T}y_{ij}+\big{(}\max(\tau_{ij}^{T},T_{\text{s}})-\tau_{ij}^{T}\big{)}w_{ij}\Big{)}-M_{1}(1-z_{j}),\quad\forall
j\in\mathcal{H}$ (58) $\displaystyle z_{0}=1,t_{0}^{+}=0$ (59) $\displaystyle
x_{ij}\in\\{0,1\\},y_{ij}\in\\{0,1\\},w_{ij}\in\\{0,1\\},\quad\forall(i,j)\in\mathcal{A}$
(60) $\displaystyle
z_{i}\in\\{0,1\\},u_{i}\in[0,n+1],t_{i}^{-},t_{i}^{+}\in[0,T_{\text{bl}}],l_{i}\in\mathcal{R}_{+},\quad\forall
i\in\mathcal{H}$ (61)

The ZONVRP model is presented as follows:

$\displaystyle(\text{ZONVRP})\quad\min\quad$
$\displaystyle\sum_{i\in\mathcal{H}}l_{i}+T_{s}\sum_{i\in
H}z_{i}-T_{s}\sum_{(i,j)\in A}2w_{ij}$ s.t. constraints: (1)–(14), (48)–(58)
$\displaystyle z_{0}=1,t_{0}^{+}=0$ (62) $\displaystyle
x_{ij}\in\\{0,1\\},y_{ij}\in\\{0,1\\},w_{ij}\in\\{0,1\\},\quad\forall(i,j)\in\mathcal{A}$
(63) $\displaystyle
z_{i}\in\\{0,1\\},u_{i}\in[0,n+1],t_{i}^{-},t_{i}^{+}\in[0,T_{\text{bl}}],l_{i}\in\mathcal{R}_{+},\quad\forall
i\in\mathcal{H}$ (64)

To make the proof self-contained, we introduce the notations that are applied
in TDTL model and reproduce the TDTL model from González-Rodríguez et al.,
(2020).

Table A.1: Notations used in the TDTL MIP formulation.  Sets
---
$\mathcal{N}$ | Set of nodes of graph $G=(\mathcal{N},\mathcal{A})$.
$\mathcal{A}$ | Set of directed links in $G=(\mathcal{N},\mathcal{A})$.
$o/e$ | The origin/ending node of the mission, $o,e\in\mathcal{N}$.
$\delta^{+}(i)$ | Nodes that can be reached from $i$, $i\in\mathcal{N}$.
$\delta^{-}(i)$ | Nodes that can reach to node $i$, $i\in\mathcal{N}$.
Parameters
Q | Battery capacity expressed in time units.
$t_{ij}^{T}/$ | Truck travel time at link $(i,j)\in\mathcal{A}$.
$t_{ij}^{D}$ | Drone travel time at link $(i,j)\in\mathcal{A}$.
M | A big enough constant.
Variables
$u_{ij}$ | If $u_{ij}=1$, link $(i,j)$ is traversed by the truck.
$v_{ij}$ | If $v_{ij}=1$, link $(i,j)$ is traversed by the drone.
$s_{i}$ | The earliest departure time from node $i\in\mathcal{N}$.
$b_{i}^{-}$ | Drone battery level when if drone is just coming to the node $i\in\mathcal{N}$.
$b_{i}^{+}$ | Drone battery level when the drone is just departing from node $i\in\mathcal{N}$.

$\displaystyle(\text{TDTL})\quad\min\quad$ $\displaystyle s_{e}$ s.t.
$\displaystyle\sum_{j\in\delta^{-}(i)}u_{ji}\leq 1,\quad\forall
i\in\mathcal{N}\setminus\\{o,e\\}$ (65)
$\displaystyle\sum_{j\in\delta^{+}(i)}u_{ij}-\sum_{j\in\delta^{-}(i)}u_{ji}=0,\quad\forall
i\in\mathcal{N}\setminus\\{o,e\\}$ (66)
$\displaystyle\sum_{j\in\delta^{+}(o)}u_{oj}=1$ (67)
$\displaystyle\sum_{i\in\delta^{-}(e)}u_{ie}=1$ (68)
$\displaystyle\sum_{j\in\delta^{-}(i)}v_{ji}\leq 1,\quad\forall
i\in\mathcal{N}\setminus\\{o,e\\}$ (69)
$\displaystyle\sum_{j\in\delta^{+}(i)}v_{ij}-\sum_{j\in\delta^{-}(i)}v_{ji}=0,\quad\forall
i\in\mathcal{N}\setminus\\{o,e\\}$ (70)
$\displaystyle\sum_{j\in\delta^{+}(o)}v_{oj}=1$ (71)
$\displaystyle\sum_{i\in\delta^{-}(e)}v_{ie}=1$ (72)
$\displaystyle\sum_{i\in\delta^{-}(j)}u_{ij}+\sum_{i\in\delta^{-}(j)}v_{ij}\geq
1,\quad\forall j\in\mathcal{N}\setminus\\{o\\}$ (73) $\displaystyle s_{j}\geq
s_{i}+t_{ij}^{T}u_{ij}-M\Big{(}1-u_{ij}\Big{)},\quad(i,j)\in\mathcal{A}$ (74)
$\displaystyle s_{j}\geq
s_{i}+t_{ij}^{D}v_{ij}-M\Big{(}1-v_{ij}+u_{ij}\Big{)},\quad(i,j)\in\mathcal{A}$
(75) $\displaystyle s_{o}=0$ (76) $\displaystyle b_{j}^{-}\leq
Q+M\Big{(}2-v_{ij}-u_{ij}\Big{)},\quad(i,j)\in\mathcal{A}$ (77) $\displaystyle
b_{j}^{-}\geq Q-M\Big{(}2-v_{ij}-u_{ij}\Big{)},\quad(i,j)\in\mathcal{A}$ (78)
$\displaystyle b_{j}^{+}\leq
Q+M\Big{(}2-v_{ij}-u_{ij}\Big{)},\quad(i,j)\in\mathcal{A}$ (79) $\displaystyle
b_{j}^{+}\geq Q-M\Big{(}2-v_{ij}-u_{ij}\Big{)},\quad(i,j)\in\mathcal{A}$ (80)
$\displaystyle b_{j}^{-}\leq
b_{i}^{+}-t_{ij}^{D}+M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (81) $\displaystyle b_{j}^{-}\geq
b_{i}^{+}-t_{ij}^{D}-M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (82) $\displaystyle b_{j}^{+}\leq
b_{j}^{-}+M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (83) $\displaystyle b_{j}^{+}\geq
b_{j}^{-}-M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (84) $\displaystyle b_{j}^{-}\leq
b_{i}^{+}-t_{ij}^{D}+M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (85) $\displaystyle b_{j}^{-}\geq
b_{i}^{+}-t_{ij}^{D}-M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (86) $\displaystyle b_{j}^{+}\leq
Q+M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (87) $\displaystyle b_{j}^{+}\geq
Q-M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)},\quad(i,j)\in\mathcal{A}$ (88) $\displaystyle b_{o}^{+}=Q$
(89)

Let the polyhedron of the linear relaxation of models ZONVRP and TDTL be
defined by

$\displaystyle
P(\text{ZONVRP})=\Big{\\{}(l,t^{-},t^{+},z,u,w,x,y)\in\mathcal{R}^{5(n+1)}_{\geq
0}\times[0,1]^{3n^{2}}\big{|}\text{constraints:
(\ref{nmodel:drone_1})--(\ref{nmodel:v_3})}\Big{\\}}$ $\displaystyle
P(\text{TDTL})=\Big{\\{}(s,b^{-},b^{+},u,v)\in\mathcal{R}^{3(n+1)}_{\geq
0}\times[0,1]^{2n^{2}}\big{|}\ \text{constraints:
(\ref{g:1})--(\ref{g:25})}\Big{\\}}$

Given the linear transformation $\Phi$ in (90)–(94) below, we will show that
each constraint in $P(\text{TDTL})$ is implied by that in $P(\text{ZONVRP})$.
We use the notation TDTL.(k) to refer to the constraint (k) in the TDTL model.

$\displaystyle x_{ij}=v_{ij},\quad\forall(i,j)\in\mathcal{A}$ (90)
$\displaystyle y_{ij}=u_{ij},\quad\forall(i,j)\in\mathcal{A}$ (91)
$\displaystyle z_{j}=\sum_{i:(i,j)\in\mathcal{A}}u_{ij},\quad\forall
j\in\mathcal{N}$ (92) $\displaystyle t_{i}^{-}=Q-b_{i}^{-},\quad\forall
i\in\mathcal{N}$ (93) $\displaystyle t_{i}^{+}=Q-b_{i}^{+},\quad\forall
i\in\mathcal{N}$ (94)

TDTL.(65)

$\sum_{j\in\delta^{-}(i)}u_{ji}=\sum_{j:(j,i)\in\mathcal{A}}y_{ji}=z_{i}\leq
1$

TDTL.(66)–(68) are implied by constraint (42)–(44).
TDTL.(69)–(72) are implied by constraints (37)–(41) directly.
TDTL.(73)

$\sum_{i\in\delta^{-}(j)}u_{ij}+\sum_{i\in\delta^{-}(j)}v_{ij}=z_{j}+\sum_{i:(i,j)\in\mathcal{A}}x_{ij}=z_{j}+1\geq
1$

TDTL.(74) is in charge of setting the departure time at node $j$ given that
arc $(i,j)$ is traversed by the truck.

$\displaystyle s_{j}-s_{i}-t_{ij}^{T}u_{ij}+M(1-u_{ij})$
$\displaystyle=s_{j}-s_{i}-\tau_{ij}^{T}y_{ij}+M(1-y_{ij})$ $\displaystyle\geq
s_{j}-s_{i}-\tau^{T}_{ij},\quad\text{implied by\ }y_{ij}=1$
$\displaystyle=\max\Big{\\{}\tau_{ij}^{T},\max\\{\tau_{ij}^{T},T_{s}\\}\Big{\\}}-\tau^{T}_{ij},\text{where
$s_{j}-s_{i}$ is no less than the truck travel time from location i to j}$
$\displaystyle\geq 0$

TDTL.(75) regulates the departure time at node $j$ given that arc $(i,j)$ is
traversed by the drone.

$\displaystyle s_{j}-s_{i}-t_{ij}^{D}v_{ij}+M(1-v_{ij}+u_{ij})$
$\displaystyle=s_{j}-s_{i}-\tau_{ij}^{D}x_{ij}+M(1-x_{ij}+y_{ij})$
$\displaystyle\geq s_{j}-s_{i}-\tau_{ij}^{D},\quad\text{implied by\
}x_{ij}=1,y_{ij}=0$ $\displaystyle=\tau_{ij}^{D}-\tau_{ij}^{D}$
$\displaystyle\geq 0$

TDTL.(76) is implied by constraint (62).

TDTL.(77)–(80) describe the drone’s battery level at the time it just arrives
at or departs from a node $j$.

$\displaystyle
b_{j}^{-}-Q-M(2-v_{ij}-u_{ij})=-t_{j}^{-}-M(2-x_{ij}-y_{ij})\leq 0$
$\displaystyle
b_{j}^{-}-Q+M(2-v_{ij}-u_{ij})=-t_{j}^{-}+M(2-x_{ij}-y_{ij})\geq
0,\quad\text{implied by constraints (\ref{nmodel:w_1})--
(\ref{nmodel:time_4})}$ $\displaystyle
b_{j}^{+}-Q-M(2-v_{ij}-u_{ij})=-t_{j}^{+}-M(2-x_{ij}-y_{ij})\leq 0$
$\displaystyle
b_{j}^{+}-Q+M(2-v_{ij}-u_{ij})=-t_{j}^{+}+M(2-x_{ij}-y_{ij})\geq
0,\quad\text{implied by\ }\text{constraints (\ref{nmodel:time_5})--
(\ref{nmodel:time_4})}$

TDTL.(81), $\forall(i,j)\in\mathcal{A}$

$\displaystyle
b_{j}^{-}-b_{i}^{+}+\tau_{ij}^{D}-M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)}$
$\displaystyle=t_{i}^{+}-t_{j}^{-}+\tau_{ij}^{D}-M(1-x_{ij}+z_{j})$
$\displaystyle\leq t_{i}^{+}-t_{j}^{-}+\tau_{ij}^{D},\quad\text{implied by\
}x_{ij}=1,z_{j}=0$ $\displaystyle=-\tau_{ij}^{D}+\tau_{ij}^{D}$
$\displaystyle\leq 0$

TDTL.(82), $\forall(i,j)\in\mathcal{A}$

$\displaystyle
b_{j}^{-}-b_{i}^{+}+\tau_{ij}^{D}+M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)}$
$\displaystyle=t_{i}^{+}-t_{j}^{-}+\tau_{ij}^{D}+M(1-x_{ij}+z_{j})$
$\displaystyle\geq t_{i}^{+}-t_{j}^{-}+\tau_{ij}^{D},\quad\text{implied by\
}x_{ij}=1,z_{j}=0$ $\displaystyle=-\tau_{ij}^{D}+\tau_{ij}^{D}$
$\displaystyle\geq 0$

TDTL.(83), $\forall(i,j)\in\mathcal{A}$

$\displaystyle b_{j}^{+}-b_{j}^{-}-M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)}$ $\displaystyle=t_{j}^{-}-t_{j}^{+}-M(1-x_{ij}+z_{j})$
$\displaystyle\leq t_{j}^{-}-t_{j}^{-}(1-z_{j}),\quad\text{implied by\
}x_{ij}=1,z_{j}=0,\text{constraint (\ref{nmodel:time_5})}$ $\displaystyle\leq
0$

TDTL.(84), $\forall(i,j)\in\mathcal{A}$

$\displaystyle b_{j}^{+}-b_{j}^{-}+M\Big{(}1-v_{ij}+u_{ij}+\sum_{k\neq
i}u_{kj}\Big{)}$ $\displaystyle=t_{j}^{-}-t_{j}^{+}+M(1-x_{ij}+z_{j})$
$\displaystyle\geq t_{j}^{-}-t_{j}^{-}(1-z_{j}),\quad\text{implied by\
}x_{ij}=1,z_{j}=0,\text{constraint (\ref{nmodel:time_5})}$ $\displaystyle\geq
0$

TDTL.(85)

$\displaystyle
b_{j}^{-}-b_{i}^{+}+\tau_{ij}^{D}-M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)}$
$\displaystyle=-t_{j}^{-}+t_{i}^{+}+t_{ij}^{D}-M\Big{(}1-x_{ij}+y_{ij}+1-\sum_{k\neq
i}y_{kj}\Big{)}$ $\displaystyle\leq-
t_{j}^{-}+t_{i}^{+}+t_{ij}^{D},\quad\text{implied by\
}x_{ij}=1,y_{ij}=0,z_{j}=1$ $\displaystyle=-t_{ij}^{D}+t_{ij}^{D}$
$\displaystyle\leq 0$

TDTL.(86), $\forall(i,j)\in\mathcal{A}$

$\displaystyle
b_{j}^{-}-b_{i}^{+}+\tau_{ij}^{D}+M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)}$
$\displaystyle=-t_{j}^{-}+t_{i}^{+}+t_{ij}^{D}+M\Big{(}1-x_{ij}+y_{ij}+1-\sum_{k\neq
i}y_{kj}\Big{)}$ $\displaystyle\geq-
t_{j}^{-}+t_{i}^{+}+t_{ij}^{D},\quad\text{implied by\
}x_{ij}=1,y_{ij}=0,z_{j}=1$ $\displaystyle=-t_{ij}^{D}+t_{ij}^{D}$
$\displaystyle\geq 0$

TDTL.(87), $\forall(i,j)\in\mathcal{A}$

$\displaystyle b_{j}^{+}-Q-M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)}$
$\displaystyle=-t_{j}^{+}-M\Big{(}1-x_{ij}+y_{ij}+1-\sum_{k\neq
i}y_{kj}\Big{)}$ $\displaystyle\leq-t_{j}^{-}(1-z_{j}),\quad\text{implied by\
}x_{ij}=1,y_{ij}=0,z_{j}=1$ $\displaystyle\leq 0$

TDTL.(88), $\forall(i,j)\in\mathcal{A}$

$\displaystyle b_{j}^{+}-Q+M\Big{(}1-v_{ij}+u_{ij}+1-\sum_{k\neq
i}u_{kj}\Big{)}$
$\displaystyle=-t_{j}^{+}-M\Big{(}1-x_{ij}+y_{ij}+1-\sum_{k\neq
i}y_{kj}\Big{)}$ $\displaystyle\geq-t_{j}^{-}(1-z_{j}),\quad\text{implied by\
}x_{ij}=1,y_{ij}=0,z_{j}=1$ $\displaystyle\geq 0$

TDTL.(89) is implied by constraint (62).

We have finally shown that $\Phi(P(\text{ZONVRP}))\subseteq P(\text{TDTL})$.
To further establish that the feasible region of $\Phi(P(\text{ZONVRP}))$ is
strictly contained in $P(\text{TDTL})$, it is sufficient to give a solution
that is valid in TDTL model but infeasible in ZONVRP model.

Figure A.1: Partial solutions that are avoided in ZONVRP.

In Figure A.1 (a), we depict a nested unit in which the truck visits location
$k_{n}$ before drives to rendezvous location $j$. This nested unit could be
part of any TDTL solution. However, such a solution is infeasible to the
ZONVRP model since constraints (1)–(2) require that each location has to be
visited by the drone exactly once. In this case shown in Figure A.1 (b), the
truck travels for more than $T_{\text{bl}}$ time units to the designated
destination $j$; and here the drone has depleted its battery and has been
forced to land on the ground. But the ZONVRP excludes this situation via the
constraint (12). $\Box$

##  Appendix B Proof for proposition 1

The following analysis applies to general cases where the total number of
locations $n$ is more than two.

Consider the left hand side of constraints (DL-1) with the coefficient of term
$x_{j,n+1}$ being replaced by $a_{j,n+1}$, we compute the largest possible
value for $a_{j,n+1}$ such that the inequality remains valid.

$1+(1-x_{0j})+a_{j,n+1}x_{j,n+1}\leq u_{j}$ (95)

case 1. $x_{0j}=1$. Since location $j$ is right after location $0$, we have
$u_{j}=1$. Meanwhile, location $j$ should not be connected to the designated
location $n+1$ when $n>2$ (i.e., $x_{j,n+1}=0$), otherwise, the drone misses
out locations before returning back to the depot. In this case, constraint
(95) becomes $0\leq 1$ which is valid.

Case 2. $x_{0j}=0$. Since the location $j$ is not the one right after $0$, the
order $u_{j}$ of location $j$ must satisfy $2\leq u_{j}\leq n+1$. Furthermore,
if $x_{j,n+1}=0$, we have $2\leq u_{j}$ which is valid. However, if
$x_{j,n+1}=1$, we have $u_{j}=n$ and the constraint (95) becomes
$a_{j,n+1}\leq u_{j}-2=n-2$. Therefore, we set $a_{j,n+1}=n-2$.

Next, we consider the right hand side of constraint (DL-1) with the
coefficient of term $x_{0j}$ being replaced by $b_{0j}$. Similarly, we compute
the largest possible value of $b_{0j}$.

$u_{j}\leq(n+1)-b_{0j}x_{0j}-(1-x_{j,n+1})$ (96)

Case 1.  $x_{j,n+1}=1$. Since location $j$ is the last location visited before
the drone coming back to depot $n+1$, we know the rank of location $j$ is $n$
(i.e., $u_{j}=n$). Moreover, location $j$ must not be right after depot $0$
Thus, the inequalities (96) become $n\leq n+1$ which are satisfied at all
times.

Case 2.  $x_{j,n+1}=0$. If $x_{0j}=0$, since location $j$ could be anywhere
not right after $0$ and right before $n+1$, the rank of location $j$ must
satisfies $2\leq u_{j}\leq n$ inequalities (96) become $u_{j}\leq n$ which is
satisfied at all times. However, if $x_{0j}=1$, we have $u_{j}=1$ and
$b_{0j}\leq n-1$. Therefore, we choose $b_{0j}=n-1$.

Consider the constraint (DL-2) with the coefficient of term $x_{ji}$ being
replaced by $c_{ji}$, if $c_{ji}=0$, the constraint is the same with
constraint (5)

$u_{i}-u_{j}+(n+1)x_{ij}+c_{ji}x_{ji}\leq n$ (97)

Case 1. $x_{ji}=0$. The constraint is the same with constraint (5).

Case 2. $x_{ji}=1$. This means that $x_{ij}=0$ and $u_{i}\geq u_{j}+1$ in the
drone order. Thus, constraint (DL-2) becomes $u_{i}-u_{j}+c_{ji}\leq n$. We
have $c_{ji}\leq n-1$.

##  Appendix C Proof for proposition 2

We would like to show that constraints (SD-1) – (SD-5) are valid restatement
of constraints (4) and (5). To achieve the goal, we follow the RTL technique
proposed in Sherali and Driscoll, (2002) which consists of a reformulation and
linearization steps.

We restate the MTZ subtour elimination constraints as follows:

$\displaystyle u_{j}x_{ij}=(u_{i}+1)x_{ij},\quad\forall(i,j)\in A,i\neq 0$
(98) $\displaystyle u_{j}x_{0j}=x_{0j},\quad\forall j\in H\setminus\\{0\\}$
(99) $\displaystyle u_{j}x_{j,n+1}=nx_{j,n+1},\quad\forall j\in
H\setminus\\{n+1\\}$ (100) $\displaystyle 1\leq u_{j}\leq n+1,\quad\forall
j\in H\setminus\\{0\\}$ (101)

Constraint (98) – (101) excludes any subtour in the drone route. To see that,
if $x_{ij}=1$, location $j$ increases rank by 1 as compared to that of the
location $i$. Specifically, if $x_{0j}=1$, location $j$ must have rank 1. If
$x_{j,n+1}=1$, location $j$ has rank $n$. Note that Constraint (98) – (100)
contain quadratic terms, we therefore perform linearization process by
introducing new variables.

Reformulation: We reformulate the constraints (98) – (101) by generating
additional implied constraints:

$\displaystyle(\text{r1}):u_{i}\big{(}\sum_{j:(i,j)\in
A}x_{ij}-1\big{)}=0,\quad\forall i\in H\setminus\\{n+1\\}$
$\displaystyle(\text{r2}):u_{j}\big{(}\sum_{i:(i,j)\in
A}x_{ij}-1\big{)}=0,\quad\forall j\in H\setminus\\{0\\}$
$\displaystyle(\text{r3}):(u_{i}-1)x_{ij}\geq 0,\quad\forall(i,j)\in A,i\neq
0$ $\displaystyle(\text{r4}):(n+1-u_{i})x_{ij}\geq 0,\quad\forall(i,j)\in
A,i\neq 0$ $\displaystyle(\text{r5}):(u_{i}-1)(1-x_{ij}-x_{ji})\geq
0,\quad\forall(i,j)\in A,i\neq 0$
$\displaystyle(\text{r6}):(n+1-u_{i})(1-x_{ij}-x_{ji})\geq
0,\quad\forall(i,j)\in A,i\neq 0$
$\displaystyle(\text{r7}):(u_{i}-2)(1-x_{0i}-x_{i,n+1})\geq 0,\quad\forall
i\in H\setminus\\{0,n+1\\}$
$\displaystyle(\text{r8}):(n-1-u_{i})(1-x_{0i}-x_{i,n+1})\geq 0,\quad\forall
i\in H\setminus\\{0,n+1\\}$

Linearization:

Let $y_{ij}=u_{i}x_{ij},\forall(i,j)\in A,i\neq 0$ and
$z_{ij}=u_{j}x_{ij},\forall(i,j)\in A,i\neq 0$. Constraints (98)–(100) become:

$\displaystyle\begin{cases}&z_{ij}=y_{ij}+x_{ij},\quad\forall(i,j)\in A,i\neq
0\\\ &u_{j}x_{0j}=x_{0j},\quad\forall j\in H\setminus\\{0\\}\\\
&u_{j}x_{j,n+1}=nx_{j,n+1},\quad\forall j\in H\setminus\\{n+1\\}\end{cases}$
(102)

We then linearize constraints (r1)–(r8) one by one by leveraging the equations
in (102).

r1:

$\forall i\in H\setminus\\{n+1\\}$:

$\displaystyle\qquad u_{i}\big{(}\sum_{j:(i,j)\in A}x_{ij}-1\big{)}=0$
$\displaystyle\Longleftrightarrow
u_{i}x_{i1}+u_{i}x_{i2}+\cdots+u_{i}x_{i,n+1}-u_{i}=0$
$\displaystyle\Longleftrightarrow\sum_{j:(i,j)\in A,j\neq
n+1}y_{ij}+nx_{i,n+1}-u_{i}=0$ $\displaystyle\implies\text{constraint\
(\ref{sd_1})}$

r2:

$\forall j\in H\setminus\\{0\\}$:

$\displaystyle\qquad u_{j}\big{(}\sum_{i:(i,j)\in A}x_{ij}-1\big{)}=0$
$\displaystyle\Longleftrightarrow\sum_{i:(i,j)\in A}(x_{ij}+y_{ij})-u_{j}=0$
$\displaystyle\Longleftrightarrow\sum_{i:(i,j)\in A}y_{ij}+1-u_{j}=0$
$\displaystyle\implies\text{constraint\ (\ref{sd_2})}$

r3:

$\forall(i,j)\in A,i\neq 0$:

$\displaystyle\qquad(u_{i}-1)x_{ij}\geq 0$ $\displaystyle\Longleftrightarrow
x_{ij}\leq y_{ij}$ $\displaystyle\implies\text{constraint\
(\ref{sd_3})}\text{\ lower bound}$

r4:

$\forall(i,j)\in A,i\neq 0$:

$\displaystyle\qquad(n+1-u_{i})x_{ij}\geq 0$
$\displaystyle\Longleftrightarrow(n+1)x_{ij}-y_{ij}\geq 0$
$\displaystyle\Longleftrightarrow y_{ij}\leq(n+1)x_{ij}$
$\displaystyle\implies\text{constraint\ (\ref{sd_3})}\text{\ upper bound}$

r5:

$\forall(i,j)\in A,i\neq 0$:

$\displaystyle\qquad(u_{i}-1)(1-x_{ij}-x_{ji})\geq 0$
$\displaystyle\Longleftrightarrow u_{i}-y_{ij}-z_{ij}-1+x_{ij}+x_{ji}\geq 0$
$\displaystyle\Longleftrightarrow y_{ij}+y_{ji}\leq u_{i}-(1-x_{ij})$
$\displaystyle\implies\text{constraint\ (\ref{sd_4})}\text{\ upper bound}$

r6:

$\forall(i,j)\in A,i\neq 0$:

$\displaystyle\qquad(n+1-u_{i})(1-x_{ij}-x_{ji})\geq 0$
$\displaystyle\Longleftrightarrow(n+1)-(n+1)x_{ij}-(n+1)x_{ji}-u_{i}+y_{ij}+x_{ji}+y_{ji}\geq
0$ $\displaystyle y_{ij}+y_{ji}\geq u_{i}+(n+1)(x_{ij}-1)+nx_{ji}$
$\displaystyle\implies\text{constraint\ (\ref{sd_4})}\text{\ lower bound}$

r7:

$\forall i\in H\setminus\\{0,n+1\\}$:

$\displaystyle\qquad(u_{i}-2)(1-x_{0i}-x_{i,n+1})\geq 0$
$\displaystyle\Longleftrightarrow u_{i}\geq
2+x_{0i}+nx_{i,n+1}-2x_{0i}-2x_{i,n+1}$ $\displaystyle\Longleftrightarrow
u_{i}\geq 1+(1-x_{0i})+(n-2)x_{i,n+1}$ $\displaystyle\implies\text{constraint\
(\ref{sd_5})}\text{\ lower bound}$

r8:

$\forall i\in H\setminus\\{0,n+1\\}$:

$\displaystyle\qquad(n-1-u_{i})(1-x_{0i}-x_{i,n+1})\geq 0$
$\displaystyle\Longleftrightarrow
n-1-(n-1)x_{0i}-(n-1)x_{i,n+1}-u_{i}+x_{0i}+nx_{i,n+1}$
$\displaystyle\Longleftrightarrow u_{i}\leq n-(n-2)x_{0i}-(1-x_{i,n+1})$
$\displaystyle\implies\text{constraint (\ref{sd_5})}\text{\ upper bound}$

##  Appendix D Proof of Theorem 3

Given a Nested-VRP described in graph $\mathcal{G}=(\mathcal{H},\mathcal{A})$,
the optimal solution to the problem consists of the drone route
$\mathcal{X}=\\{(i,j)\ |\ (i,j)\in\mathcal{A},\ \text{and}\ x_{ij}=1\\}$ and
the truck route $\mathcal{Y}=\\{(i,j)\ |\ (i,j)\in\mathcal{A},\ \text{and}\
y_{ij}=1\\}$. The optimal solution can also be described as a set of non-
overlapping nested units $U$. Mathematically, a nested unit $u\in U$ can be
viewed as a sub-graph of $\mathcal{G}$. The sub graph contains a set of
locations $V(u)$ to be observed by the drone and the corresponding drone path
$E(u)$. Most importantly, one can see that $\mathcal{X}\equiv\\{(i,j)\ |\
(i,j)\in E(u),u\in U\\}$ and $\mathcal{H}\equiv\\{k\ |\ k\in V(u),u\in U\\}$.

In a nested unit $u$, when the drone is about to meet with the truck, the
drone arrives at the rendezvous location either earlier than the truck and
idles for $\Delta_{u}$ time units or later than the truck without idling. We
define an indicator function $\mathbbm{1}_{u}$ for each nested unit $u$ that
forms the optimal Nested-VRP solution. In a nested unit $u$, if the drone
arrives earlier than the truck at rendezvous, $\mathbbm{1}_{u}=1$, otherwise
$\mathbbm{1}_{u}=0$. After two vehicles meet up successfully, the drone
relinquishes all remaining battery life before it obtains a new battery. This
portion of unused battery life is denoted as $\delta_{u}$.

A battery is either used for drone surveillance (e.g., routing between
locations, surveying locations, and possibly waiting for the truck) or wasted
after the two vehicles meet up. To derive the lower bound of the mission
makespan of a Nested-VRP solution, we first investigate the battery
consumption associated with a Nested-VRP solution. Denote the IBR $l_{u}$ of a
nested unit as in Equation (103). We summarize the battery consumption
breakdowns in Equation (104).

$\displaystyle l_{u}$ $\displaystyle=$ $\displaystyle\sum_{(i,j)\in
E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in V(u)}o_{k}+\mathbbm{1}_{u}\Delta_{u}$
(103) $\displaystyle T_{\text{bl}}$ $\displaystyle=$ $\displaystyle
l_{u}+\delta_{u}$ (104)

Due to the conservation of energy, the amount of energy extracted from all
battery replacements scheduled en route should be able to balance off the
total battery consumption needed for the mission. Therefore, we can derive the
necessary number of battery swaps $N_{\text{s}}$ as follows:

$\displaystyle N_{\text{s}}T_{\text{bl}}$ $\displaystyle=\sum_{u\in
U}(l_{u}+\delta_{u})=\sum_{u\in U}\Big{(}\sum_{(i,j)\in
E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in
V(u)}o_{k}+\mathbbm{1}_{u}\Delta_{u}+\delta_{u}\Big{)}$
$N_{\text{s}}=\frac{1}{T_{\text{bl}}}\sum_{u\in U}\Big{(}\sum_{(i,j)\in
E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in
V(u)}o_{k}+\mathbbm{1}_{u}\Delta_{u}+\delta_{u}\Big{)}$

Since the mission makespan is the sum of the IBRs $l_{u},\forall u\in U$ and
battery swap service times:

makespan $\displaystyle=\sum_{u\in U}l_{u}+N_{\text{s}}T_{\text{s}}\quad$
$\displaystyle=\sum_{u\in U}\Big{(}\sum_{(i,j)\in
E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in
V(u)}o_{k}+\mathbbm{1}_{u}\Delta_{u}\Big{)}+\frac{T_{\text{s}}}{T_{\text{bl}}}\sum_{u\in
U}\Big{(}\sum_{(i,j)\in E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in
V(u)}o_{k}+\mathbbm{1}_{u}\Delta_{u}+\delta_{u}\Big{)}$

We relax constraints imposed on how battery swap could happen and allow the
drone to charge itself at any time when its battery has depleted. Therefore,
the mission makespan of a relaxed version of Nested-VRP instance is given by:

makespan $\displaystyle=\sum_{u\in U}\Big{(}\sum_{(i,j)\in
E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in
V(u)}o_{k}\Big{)}+\frac{T_{\text{s}}}{T_{\text{bl}}}\sum_{u\in
U}\Big{(}\sum_{(i,j)\in E(u)}\tau_{ij}^{D}x_{ij}+\sum_{k\in V(u)}o_{k}\Big{)}$
$\displaystyle=\sum_{(i,j)\in\mathcal{X}}\tau_{ij}^{D}x_{ij}+\sum_{k\in\mathcal{H}}o_{k}+\frac{T_{\text{s}}}{T_{\text{bl}}}\Big{(}\sum_{(i,j)\in\mathcal{X}}\tau_{ij}^{D}x_{ij}+\sum_{k\in\mathcal{H}}o_{k}\Big{)}$
(105)
$\displaystyle\geq\sum_{(i,j)\in\mathcal{S}}\tau_{ij}^{D}+\sum_{k\in\mathcal{H}}o_{k}+\left\lfloor\frac{1}{T_{\text{bl}}}\Big{(}\sum_{(i,j)\in\mathcal{S}}\tau_{ij}^{D}+\sum_{k\in\mathcal{H}}o_{k}\Big{)}\right\rfloor
T_{\text{s}},$ (106)

Recall that $\mathcal{S}$ is the collection of arcs that are in the TSP route.
Since the union of the drone paths in all nested units will produce a
Hamiltonian cycle containing all of the locations and whose total length is no
shorter than the TSP route, Equation (105) will be lower bounded by specifying
the drone route as the TSP route as well as rounding down the total number of
battery swaps to an integer. The result is that, the objective function value
of the original Nested-VRP model is lower bounded by term (106). $\Box$

##  Appendix E Computational results on small data set

Table E.2: Results from solving instances from the uniform pattern in the
small data set.

uniform | MIQCP |  | MILP+SD |  | LB |  | NS Heuristic
---|---|---|---|---|---|---|---
$N$ | $C_{\text{MIQCP}}(s)$ | $\gamma_{\text{MIQCP}}(\%)$ | $T_{\text{MIQCP}}(s)$ |  | $C_{\text{MILP+SD}}(s)$ | $\gamma_{\text{MILP+SD}}(\%)$ | $T_{\text{MILP+SD}}(s)$ |  | $\gamma_{\text{lb}}(\%)$ |  | $N^{\star}$ | $T_{\text{NS}}(s)$
$\alpha=1$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1615.5 | 0.0 | 0.4 |  | 1615.5 | 0.0 | 0.4 |  | -0.8 |  | 10/10 | 0.7
6 | 1795.1 | 0.0 | 2.2 |  | 1795.1 | 0.0 | 2.7 |  | -0.8 |  | 10/10 | 1.6
7 | 2313.9 | 0.0 | 26.4 |  | 2313.9 | 0.0 | 16.3 |  | -1.2 |  | 10/10 | 2.0
8 | 2484.2 | 0.0 | 246.4 |  | 2484.2 | 0.0 | 59.8 |  | -1.5 |  | 10/10 | 3.5
9 | 2901.8 | 36.6 | 901.7 |  | 2900.7 | 1.6 | 446 |  | -1.4 |  | 10/10 | 4.1
10 | 3004.9 | 54.8 | 901.8 |  | 3002.5 | 54.4 | 900.4 |  | -1.6 |  | 10/10 | 5.9
$\alpha=2$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1431.9 | 0.0 | 1.0 |  | 1431.9 | 0.0 | 0.5 |  | -3.7 |  | 10/10 | 0.7
6 | 1577.6 | 0.0 | 4.5 |  | 1577.6 | 0.0 | 2.1 |  | -2.2 |  | 10/10 | 1.4
7 | 1915.1 | 0.0 | 43.8 |  | 1915.1 | 0.0 | 11.1 |  | -2.9 |  | 10/10 | 2.5
8 | 2305.1 | 0.0 | 374.8 |  | 2305.1 | 0.0 | 72.2 |  | -2.8 |  | 10/10 | 3.5
9 | 2596.9 | 36.8 | 902.2 |  | 2591.8 | 5.1 | 475.4 |  | -3.0 |  | 10/10 | 4.9
10 | 2364.1 | 60.2 | 901.9 |  | 2361.1 | 52.8 | 900.4 |  | -2.9 |  | 10/10 | 7.6
$\alpha=3$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1474.8 | 0.0 | 1.1 |  | 1474.8 | 0.0 | 0.4 |  | -4.2 |  | 10/10 | 0.9
6 | 1414.8 | 0.0 | 4.6 |  | 1414.8 | 0.0 | 0.8 |  | -2.4 |  | 10/10 | 1.4
7 | 1880 | 0.0 | 30.5 |  | 1880 | 0.0 | 5.5 |  | -4.1 |  | 10/10 | 2.3
8 | 2227.1 | 0.0 | 231.1 |  | 2227.1 | 0.0 | 62.8 |  | -1.9 |  | 10/10 | 3.1
9 | 2562.5 | 40.4 | 902.0 |  | 2562.4 | 3.5 | 500.9 |  | -3.2 |  | 9/10 | 4.7
10 | 2449 | 62.6 | 901.8 |  | 2429.1 | 47.0 | 900.3 |  | -3.0 |  | 9/10 | 6.5
Summary | 120/180 of instances reach optimality |  | 145/180 instances reach optimality |  |  |  |  |

Table E.3: Results from solving instances from the single-center pattern in
the small data set.

single-center | MIQCP |  | MILP+SD |  | LB |  | NS Heuristic
---|---|---|---|---|---|---|---
$N$ | $C_{\text{MIQCP}}(s)$ | $\gamma_{\text{MIQCP}}(\%)$ | $T_{\text{MIQCP}}(s)$ |  | $C_{\text{MILP+SD}}(s)$ | $\gamma_{\text{MILP+SD}}(\%)$ | $T_{\text{MILP+SD}}(s)$ |  | $\gamma_{\text{lb}}(\%)$ |  | $N^{\star}$ | $T_{\text{NS}}(s)$
$\alpha=1$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1546 | 0.0 | 0.4 |  | 1546 | 0.0 | 0.3 |  | -1.2 |  | 10/10 | 0.7
6 | 2156.7 | 0.0 | 2.0 |  | 2156.7 | 0.0 | 2.3 |  | -1.6 |  | 10/10 | 1.2
7 | 2233.6 | 0.0 | 19.4 |  | 2233.6 | 0.0 | 19.0 |  | -1.3 |  | 10/10 | 1.8
8 | 2562.4 | 0.0 | 194.4 |  | 2562.4 | 0.0 | 178.5 |  | -1.3 |  | 10/10 | 2.9
9 | 2964.2 | 34.9 | 901.7 |  | 2960.3 | 15.0 | 778.5 |  | -1.6 |  | 10/10 | 4.0
10 | 3471.4 | 51.0 | 901.8 |  | 3461.2 | 32.3 | 900.5 |  | -2.3 |  | 9/10 | 5.8
$\alpha=2$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1340.6 | 0.0 | 0.9 |  | 1340.6 | 0.0 | 0.5 |  | -3.1 |  | 10/10 | 1.2
6 | 1741.3 | 0.0 | 5.0 |  | 1741.3 | 0.0 | 1.2 |  | -4.9 |  | 10/10 | 1.6
7 | 1917.1 | 0.0 | 32.5 |  | 1917.1 | 0.0 | 9.1 |  | -2.3 |  | 10/10 | 2.3
8 | 2241.6 | 0.0 | 363.3 |  | 2241.6 | 0.0 | 67.3 |  | -4.0 |  | 9/10 | 3.2
9 | 2348.9 | 37.7 | 902.4 |  | 2338.8 | 6.4 | 301.9 |  | -3.3 |  | 8/10 | 3.8
10 | 3213.4 | 57.5 | 901.6 |  | 3191.1 | 1.8 | 406.2 |  | -3.3 |  | 7/10 | 6.0
$\alpha=3$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1299.6 | 0.0 | 1.1 |  | 1299.6 | 0.0 | 0.4 |  | -2.4 |  | 10/10 | 1.3
6 | 1484.9 | 0.0 | 4.5 |  | 1484.9 | 0.0 | 0.7 |  | -3.8 |  | 10/10 | 1.3
7 | 1874.1 | 0.0 | 25.8 |  | 1874.1 | 0.0 | 3.0 |  | -3.0 |  | 10/10 | 2.2
8 | 2418 | 0.0 | 268.8 |  | 2418 | 0.0 | 27.0 |  | -2.8 |  | 9/10 | 2.2
9 | 2475.5 | 37.5 | 902.0 |  | 2468.8 | 1.8 | 104.7 |  | -2.8 |  | 10/10 | 5.2
10 | 2713.9 | 61.8 | 901.8 |  | 2673.5 | 0.9 | 245.8 |  | -1.6 |  | 4/10 | 6.6
Summary | 120/180 of instances reach optimality |  | 156/180 instances reach optimality |  |  |  |  |

Table E.4: Results from solving instances from the double-center pattern in
the small data set.

double-center | MIQCP |  | MILP+SD |  | LB |  | NS Heuristic
---|---|---|---|---|---|---|---
$N$ | $C_{\text{MIQCP}}(s)$ | $\gamma_{\text{MIQCP}}(\%)$ | $T_{\text{MIQCP}}(s)$ |  | $C_{\text{MILP+SD}}(s)$ | $\gamma_{\text{MILP+SD}}(\%)$ | $T_{\text{MILP+SD}}(s)$ |  | $\gamma_{\text{lb}}(\%)$ |  | $N^{\star}$ | $T_{\text{NS}}(s)$
$\alpha=1$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 2645.8 | 0.0 | 0.4 |  | 2645.8 | 0.0 | 0.2 |  | -1.0 |  | 10/10 | 0.8
6 | 2889.1 | 0.0 | 1.2 |  | 2889.1 | 0.0 | 0.9 |  | -0.6 |  | 10/10 | 1.1
7 | 3254 | 0.0 | 9.8 |  | 3254 | 0.0 | 6.2 |  | -0.8 |  | 10/10 | 1.9
8 | 3687.1 | 0.0 | 78.0 |  | 3687.1 | 0.0 | 41.3 |  | -1.2 |  | 9/10 | 2.6
9 | 3568.4 | 15.5 | 746.6 |  | 3567.3 | 5.0 | 394.9 |  | -0.9 |  | 10/10 | 4.7
10 | 4179.2 | 31.4 | 901.6 |  | 4176.4 | 7.6 | 900.5 |  | -1.4 |  | 10/10 | 5.0
$\alpha=2$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1855.7 | 0.0 | 0.5 |  | 1855.7 | 0.0 | 0.2 |  | -3.9 |  | 10/10 | 1.1
6 | 2429.1 | 0.0 | 3.5 |  | 2429.1 | 0.0 | 0.3 |  | -4.8 |  | 10/10 | 1.4
7 | 2354.9 | 0.0 | 23.4 |  | 2354.9 | 0.0 | 0.9 |  | -3.2 |  | 10/10 | 2.0
8 | 3027.8 | 0.0 | 372.5 |  | 3027.8 | 0.0 | 6.7 |  | -4.5 |  | 10/10 | 3.0
9 | 3296.2 | 31.2 | 901.8 |  | 3289.1 | 0.0 | 73.6 |  | -4.5 |  | 10/10 | 4.5
10 | 3362.6 | 50.5 | 901.6 |  | 3353.2 | 0.0 | 417.9 |  | -4.6 |  | 8/10 | 6.4
$\alpha=3$ |  |  |  |  |  |  |  |  |  |  |  |
5 | 1612.6 | 0.0 | 0.4 |  | 1612.6 | 0.0 | 0.1 |  | -2.6 |  | 9/10 | 1.0
6 | 1711.9 | 0.0 | 1.8 |  | 1711.9 | 0.0 | 0.3 |  | -2.7 |  | 10/10 | 1.4
7 | 2350.4 | 0.0 | 20.8 |  | 2350.4 | 0.0 | 0.6 |  | -2.2 |  | 9/10 | 2.2
8 | 2429.4 | 0.0 | 142.7 |  | 2429.4 | 0.0 | 2.6 |  | -2.3 |  | 10/10 | 3.3
9 | 2806.4 | 35.7 | 901.8 |  | 2794.9 | 0.0 | 23.3 |  | -3.0 |  | 9/10 | 4.8
10 | 3273.3 | 53.8 | 902.3 |  | 3251.7 | 0.0 | 93.9 |  | -2.6 |  | 6/10 | 6.3
Summary | 123/180 of instances reach optimality |  | 173/180 instances reach optimality |  |  |  |  |