In this paper we deal with an extension of the crossing postman problem to design Hamiltonian routes that have to visit different shapes of dimensional elements (neighborhoods or polygonal chains) rather than edges. This problem models routes of drones that must visit a number of geographical elements to deliver some good or service and then move directly to the next target element using straight line displacements. We present two families of mathematical programming formulations. The first one is time-dependent and captures a number of actual characteristics of real applications at the price of using three indexes variables. The second one are not referenced to the stages of the route. We compare them on a testbed of randomly generated instances with different shapes of elements: second order cone representable (SOC) and polyhedral neighborhoods and polygonal chains. The computational result reported in this paper show that our models are useful and can solve to optimality medium size instances of sizes similar to other combinatorial problems with neighborhoods. To address larger instances we also present a heuristic algorithm that runs in two phases: clustering and VNS. This algorithm performs very well in quality of solutions provided and can be used to initialize the exact methods with promising initial solutions.

§ INTRODUCTION

Drones, or UAVs (unmanned aerial vehicles), provide new opportunities for improving logistics in a variety of settings. We refer them as drones since the main characteristic that we wish to exploit is their aerial displacement using straight lines.
Recent technological improvements as battery life, better communication devices and reduction in manufacturing costs have increased the use of drones in  logistics. Drones have been used in many different fields as disaster management in remote regions ([11]), parcel delivery [12], communication coverage ([1]), traffic monitoring, infrastructure inspection, coastal surveying and many other applications. The reader is referred to the review [18] for further references.

The availability of this new technology has brought  new business opportunities and at the same time has  opened a lot of new challenges in the Operations Research field to propose solutions to new emerging problems in the areas logistics and routing. As drones play a growing role in business operations, questions of planning and optimization increase in practical and
academic importance. Actual characteristics of drone's displacement make most previous routing models in literature not readily applicable. Unlike standard ground vehicles, that must follow paths,  one of the specific characteristic of drones is their ability  to use direct connections by straight lines between targeted destinations because they can fly across areas.

In 1962, Meigu Guan introduced the undirected Chinese Postman Problem (CPP) whose aim is to determine a least-cost closed route that traverses all edges of the graph. Orloff [17] extended the CPP to travel through a subset of required edges that is known as the Rural Postman Problem (RPP). Based on this idea, Garfinkel [7] relaxed the RPP to the case in which it is permitted to leave the edges of the network and cross from one edge to another at points other than the original vertices. These ARPs are studied in depth in [5]. On the other hand, this problem takes the structure of the well-known Traveling Salesman Problem that is studied in [8] using convex sets. Yuan [22] presented in his work a hybrid framework in which metaheuristics and classical TSP solvers are combined strategically to produce high quality solutions for TSPN with arbitrary neighborhoods.

The aim of this paper is motivated by the design of drones' routes that must connect a number of targets with given shapes,  that we will call from now on elements,  that are spread on an area. These elements are dimensional. This means that they can not be represented, in general, by points. Moreover, they have some requirements of service beyond its simple visit: either one has to visit a percentage of its total length (assuming that their dimension is one) or one has to cross at least some distance over them (if they are two dimensional). Some previous attempt that can be applied to this problem is to approximate  all the shapes by polygons with sufficiently large number of breakpoints (see [4]). Nevertheless, we would like to exploit some new features of MINLP to address an alternative approach. Obviously, we have to impose some limit to the shapes of the considered elements to achieve tractable models. As a first building block, we restrict ourselves to two main types of elements (see Figure <ref>): convex bodies and piecewise linear chains (including segments). This limitation will be extended to more general shapes at the end of the paper.

An example of convex sets and polygonal chains considered in the problem

We assume a structure of costs trying to capture actual elements in these situations. The costs of moving between elements may change over time: it may be cheaper to go from A to B at time 1 than at time 2. The cost of moving on the elements may be cheaper or more expensive than moving freely over the ground: controlling the drone over polygonal chains to do some inspection may be more expensive than flying directly between targets. On the other hand, one can obtain some discount for flying over some large area (parks, lakes, natural reserves...) because the drone can do a secondary job, as reporting information, while doing its primary delivery job.

The goal is to find a Hamiltonian route that connects all the elements satisfies the visits' requirements and is of minimum total cost.

The contribution of this paper is the combination on the same model of different elements that have never been put together before: design of Hamiltonian routes without underlying graph structure, required targets defined on elements (like in the RPP), free entry and exit points over the elements, use of dimensional elements (as in the TSP with neighborhoods). Combining these features altogether gives rise to a challenging new problem that is analyzed for the first time in this paper.

The paper is structured in 8 sections. The first section is the introduction. In the second section we describe the problem and set the notation followed in the rest of the paper. Section <ref> is devoted to present different valid formulations of the problem.
In Section <ref> we present a heuristic algorithm for solving XPPN. This heuristic has two phases clustering and VNS. The results show that it provides good quality solution in very limited computation time. Section <ref> deals with some strengthening of our formulations: pre-processing variables and deriving valid inequalities to be added to the formulations. Next, in Section <ref> we present a decomposition algorithm 'a la' Benders that can be also applied to solve the problem. We derive all the details of that decomposition and show preliminary computational results.

An extensive computational experience is reported in Section <ref>. There, we compare the different formulations in terms of final gaps and computing time. The paper ends with a section devoted to conclusions and extensions, where we list some interesting open lines of research connected with the problems addressed in this paper.

§ DESCRIPTION OF THE PROBLEM

Let $G=(V, E)$ be a connected undirected graph, whose vertices are embedded in $\mathbb R^2$. (The reader may note that extensions to the three dimensional space are possible at the price of increasing the models' complexity.) Associated with each vertex $v\in V$, we attach an element that can belong to two different types: either a convex set or a polygonal chain. In the former case, let $\mathcal C_v\subset\mathbb R^2$ denote the convex set associated to $v$ that must contain $v$ in its interior. In the latter, let  $\mathcal P_v\subset\mathbb R^2$ denote the polygonal chain attached to $v$ that we assume to be parameterized by its breakpoints $A_v^1,\ldots,A_v^{n_v+1}$, where $n_v$ is the number of line segments of the polygonal chain. We denote

$$V_{\mathcal C}=\{v\in V : v \text{ is associated with a convex set} \}, $$

$$V_{\mathcal P}=\{v\in V : v \text{ is associated with a polygonal chain} \}. $$

Feasible solutions to the Crossing Postman Problem with Neighborhoods (XPPN) problem consists of
a set of pairs of points, $X=\{(x_v^1, x_v^2) : v\in V\}$, that denotes the access and exit points to the component associated with vertices $v\in V$, together with a Hamiltonian tour $\mathcal T$ on the graph $G' = (X, E')$, with edge set $E'=E_{\text{out}} \cup E_{\text{in}}$, where:

$$E_\text{out}= \{(x_v^1, x_w^2) : (v, w)\in E\},\quad E_\text{in} = \{(x_v^1, x_v^2) : v \in V\}.$$

Edges in the set $E_\text{out}$ are links between different elements whereas those in $E_\text{in}$ are those that define the part of the tour that is traveled within the convex neighborhoods or the polygonal chains.
Edges lengths are given by the Euclidean distance, $\|\cdot\|_2$, between their endpoints.

The overall cost of a feasible solution $(X, \mathcal T)$ is then given by the overall sum of outside edges plus the discounted sum of the inner edges:

$$d(X,\mathcal T) = \sum_{e_{vw}=(x_v^1, x_w^2)\in\mathcal T} d(x_v^1, x_w^2) + \sum_{e_v=(x_v^1, x_v^2)\in\mathcal T} f_v\,d(x_v^1, x_v^2),$$

where $f_v$ is a discount factor for traveling within the neighborhoods.

Throughout this paper we employ the following notation:

* $\mathcal{T}_G$ as the set of incidence vectors associated with Hamiltonian tours
on $G$, i.e., $\mathcal{T}_G=\{z\in\mathbb{R}_+^{|E|}:z\text{ is a Hamiltonian tour on } G\}$.

* $\mathcal{X}=\left(\prod_{v\in V_{\mathcal{C}}}\mathcal C_v\right)\cup \left(\prod_{v\in V_{\mathcal{P}}}\mathcal P_v\right)$, where $\mathcal C_v$ (resp. $\mathcal P_v$) is the neighborhood (resp. polygonal chain) associated to vertex $v$, which contains the possible sets of vertices for the Hamiltonian tours of XPPN.

The goal of XPPN is to find a feasible solution $(X, \mathcal T)$ of minimal total cost. Then, it can be expressed as:

∑_e∈E_out d_e z_e + ∑_e∈E_in f_e d_e
z∈𝒯_G,  x∈𝒳

Some observations can be stated from the formulation above:

* Fixing $x\in\mathcal X$ results in a Rural Postman Problem that has been analyzed in the literature and the state of current research provides good heuristics to solve it efficiently.

* Because of the expression of its objective function, XPPN is not separable. It is suitable to represent this problem as a MINLP.

In this paper, we focus on the case where the sets $\mathcal C_v$ are second order cone (SOC) representable, that is, the sets can be expressed by using second-order cone constraints as follows:

\begin{equation}\label{C-C}\tag{$\mathcal{C}$-C}
\|B_ix + b_i\|\leq c_i^Tx + d_i,\quad i=1,\ldots,N,
\end{equation}

where $x\in\mathbb R^2$ is the decision variable and $B_i, b_i, c_i$ and $d_i$ are parameters of the constraint. Note that these constraints can model linear constraints (for $B_i, b_i\equiv 0$), ellipsoids and hyperbolic constraints (see [13] for more details).

These type of elements could be extended further to unions of SOC representable sets. This type of  neighborhood is obtained introducing binary variables, whose meaning is similar to those in disjunctive programming. Thus, we can determine in which set of the union happen the access or the departure points in those unions of sets.

Let $\{\mathcal C_v^1, \ldots, \mathcal C_v^{N_v}\}$ be the second order cone representable sets that makes the neighborhood associated to the vertex $v$ and let $\mathcal U_v = \displaystyle \bigcup_{j=1}^{N_v} \mathcal C_v^j$ denote the union of these sets. Consider the binary variable  $\chi_v^{ij}$ that assumes the value of one if $x_v^i$ is located in the set $\mathcal C_v^j$ and zero otherwise. Thus, for each $v\in V$, we can model that $x_v^i \in \mathcal U_v$ by using the following inequalities for each $i=1,2$:

\begin{equation}\label{U-C}\tag{$\mathcal U$-C}
x_v^i\in \mathcal U_v \Longleftrightarrow
\left\{
\begin{array}{cclr}
\|A_v^j x_v^i + b_v^j\|& \leq & (c_v^j)^T x_v^i + M_v^j(1-\chi_v^{ij}), & j=1,\ldots,N_v, \\
\sum_{j = 1}^{N_v} \chi_v^{ij} & =    & 1.
\end{array}
\right.
\end{equation}

The reader may observe that one can replace <ref> by <ref> in all our formulations without altering their validity. Therefore, our model can deal easily with these more general forms of neighborhoods.

On the other hand, the second type of elements are the piecewise linear constraints. Let $n_v$ be the number of line segments of the polygonal chain $v$. Since we need to refer to interior points of the segment, these continuum of points  is parametrized by the two endpoints of the segment: $x\in [A_v^j,A_v^{j+1}]$ if and only if $\exists \gamma\in [0,1]$ such that $x=\gamma A_v^j+(1-\gamma) A_v^{j+1}$. In order to deal with them, we  introduce the following variables for each vertex $v\in V_\mathcal{P}$ and $i=1,2$:

* $u_v$: Binary variable that determines the traveling direction in the polygonal chain $v$.

* $\gamma_v^{ij}$: Continuous variable in $[0,1]$ that represents the parameter value of the $x_v^i$ variable in the line segment $j$ of the polygonal chain $v$, $j=1,\ldots,n_v$.

* $\mu_v^{ij}$: Binary variable that is one when $x_v^i$ is located in the line segment $j$ of the polygonal chain $v$, and zero otherwise, for $j=1,\ldots,n_v$.

Using these variables, we can model the polygonal chain introducing the following inequalities for each $i = 1, 2$:

\begin{equation}\label{P-C}\tag{$\mathcal P$-C}
\small
x_v^i\in \mathcal P_v \Longleftrightarrow
\left\{
\begin{array}{cclr}
\lambda_v^i - j                    & \geq & \gamma_v^{ij} - (n_v+1)(1-\mu_v^{ij}),                            & j=2,\ldots,n_v+1 \\
\lambda_v^i - j                    & \leq & \gamma_v^{ij} + (n_v+1)(1-\mu_v^{ij}),                            & j=2,\ldots,n_v+1 \\
\gamma_v^{i1}                      & \leq & \mu_v^{i1} & \\
\gamma_v^{ij}                      & \leq & \mu_v^{ij-1} + \mu_v^{ij}                                         & j=2,\ldots,n_v\\
\gamma_v^{in_v}                    & \leq & \mu_v^{in_v} \\
\sum_{j = 1}^{n_v} \mu_v^{ij}      & =    & 1 \\
\sum_{j = 1}^{n_v+1} \gamma_v^{ij} & =    & 1 \\
x_v^i                              & = & \sum_{j=1}^{n_v+1}\gamma_v^{ij}A_v^j \\
\end{array}
\right.
\end{equation}

Observe that the first and second inequalities determine the lower and upper limits of each segment of $\mathcal{P}_v$. If $\mu_v^{ij}=0$ the inequalities are always fulfilled and there is no distinguished (entry or exit) point in the $j$-th segment. The third and fourth inequalities link $\mu$ and $\lambda$ variables. They state that the variable $\gamma^{ij}$ that gives the representation of a point $x_i$ on the line segment $j$ is active (non-null)  only if this line segment is chosen (to enter of exit), i.e., $\mu^{ij}=1$. The fifth equation sets that only one line segment is chosen for entering or leaving each polygonal chain.  Finally, the sixth equation and seventh inequality sets the representation of $x^i$ as a convex combination of the extreme points of the adequate line segment.

In addition, we assume that the tour must traverse at least some given percentage $\alpha_v$  of each polygonal chain  total length. Denoting by $\lambda_i$ the parameter value of $x_v^i$ in the parametrization of the polygonal chain $\mathcal{P}_v$ and $\lambda^{\text{min}}_v$ and $\lambda^{\text{max}}_v$ the parameter values $\lambda$ of the access (entry) and exit points to  $\mathcal{P}_v$, respectively, we can model that condition by the following absolute value constraint:

\begin{equation}\label{alpha-C}\tag{$\alpha$-C}
|\lambda_v^1-\lambda_v^2|\geq \alpha_v \Longleftrightarrow
\left\{
\begin{array}{ccl}
\lambda_v^1 - \lambda_v^2                       & =    & \lambda^{\text{max}}_v - \lambda^{\text{min}}_v \\
\lambda^{\text{max}}_v + \lambda^{\text{min}}_v & \geq & \alpha_v n_v                                    \\
\lambda^{\text{max}}_v                          & \leq & n_v(1-u_v)                                      \\
\lambda^{\text{min}}_v                          & \leq & n_v u_v,                                        \\
\end{array}
\right.
\end{equation}

The above modeling assumptions are sufficient to model a wide range of actual situations that appear in the routing problems with drones that we want to model. Obviously, they could be more general at the price of not being easily implemented with off-the-shelf solvers.

§.§ Some interesting particular cases

Three very interesting well-known models appear as particular cases of the problems that can be modeled within our framework. If the element associated with each vertex is the empty set the problem reduces to the standard traveling salesman problem. If the element associated with each vertex $v\in V$ is a segment $C_v= [x_v^1,x_v^2]$, then XPPN becomes the classical Rural Postman Problem in which the edges $(x_v^1, x_v^2)$ are required, in the complete graph induced by these vertices with edge lengths given by the Euclidean norm distance.
On the other hand, if the considered neighborhoods are big enough so that $\cap_{v\in V} \mathcal C_v\neq\emptyset$, then the problem reduces to finding a degenerate one-vertex tour and the solution to the XPPN is that vertex with cost 0. Finally, if all neighborhoods are line segments and no condition is imposed on required percentage of visit, we obtain the Crossing Postman Problem which is studied in [7].

Figure <ref> shows an example of the solution obtained for a case in which the elements are circles, triangles and we also have two polygonal chains to visit in our required route.

An example with 9 elements: 7 convex sets and 2 polygonal chains

The discussion above allows us to state the complexity of the XPPN.

The decision version of the problem XPPN, given a length $L$ to decide whether the graph has a XPPN tour of length at most $L$,  is NP-complete.

The proof follows using a reduction from TSP that as shown above is a particular case of this problem.

§ MIXED INTEGER NON LINEAR PROGRAMMING FORMULATIONS

In this section we present alternative MINLP formulations for the XPPN that
will be compared computationally in later sections. First, we start with a time dependent formulation that allows us to include a number of actual requirements in the modeling phase such as time dependent travel distances, time windows or time dependent discount factors. Then, we give another formulation that does not make reference to stages in the routes and that simplifies the model at the price of losing some of the above mentioned characteristics.

§.§ A Time Dependent Formulation

One way to model the routes in our problem is to make variables dependent on the index of the stage when an element is visited in the sequence of visited elements. Thus, this formulation requires binary variables depending on the index order when they are chosen. Since variables depend on time also parameters in the problem as discount factors for traveling the neighborhoods  $(b_v^t)$ and distances $d^t$ (as proxy for travel times $\beta^t$) can be dependent on the instant when they are traveled.

To model the problem, we introduce a binary variable $y_v^t$ to select that the element associated with vertex $v$ is visited at stage $t$. In addition, we define the following  variables:

* $y_v^t$: Binary variable whose value is one when $v$ is visited at the $t$ position in the route sequence and zero otherwise.

* $z_{vw}^t$:	Binary variable that is one when $v$ and $w$ are visited consecutively, assuming that $v$ is visited at the instant $t$ and zero otherwise.

* $z_{vw}^t = y_v^t\,y_{w}^{t+1}$, $v\neq w$.

* $d_{vw}^t$: Continuous variable that represents the distance between  pairs of chosen points $v,w$ from different components at the instant $t$.

* $d_v^t$: Continuous variable that represents the distance between two consecutive points within the same component associated with $v\in V$ at the instant $t$.

* $\lambda^1_v,\,\lambda^2_v$: Continuous variables determining the position of  $x_v^1$ and $x_v^2$, respectively,  in the polygonal chain $\mathcal{P}_v$.

Using these variables, the first formulation follows:
\beginmini!|s|
∑_t=1^|V|∑_v≠w d_vw^tz_vw^t+∑_t=1^|V|∑_v∈Vf_v^t d_v^t
d_vw^t≥β_uv^t x_v^2 - x_w^1,    ∀v≠w
d_v^t≥β_v^t x_v^1 - x_v^2,    ∀v∈V
∑_v∈Vy_v^t= 1,     ∀t
∑_t=1^|V|y_v^t= 1,     ∀v∈V
y_v^t + y_w^t+1-1≤z_vw^t,     ∀v≠w, t=1,…,|𝒞|-1
(<ref>), (<ref>), (<ref>)

\endmini!

The first addend of the objective function (<ref>) includes the traveling distance among different elements while the second one accounts for the distances between the entry and exit points of each component taking into account the discount factor for travelling within this component at the instant $t$. Constraints (<ref>) and (<ref>) state, respectively, that in each instant the route visits one element and each component is traversed once and only once. Constraint (<ref>) is obtained by linearizing $z_{vw}^t$ and ensures that if we travel from $v$ to $w$ assuming that we are in $v$ at the instant $t$, then we visit $v$ in $t$ and $w$ in $t+1$. Constraint (<ref>) refers to the domain of the entry and exit points of each element in the problem, as well as the minimal required percentage of the polygonal chain length that must be traversed. They were defined in Section <ref>.

Despite the versatility of this formulation for capturing actual characteristics of drone routes, its drawback comes from the three index dimension of its variables which makes it difficult to handle medium size instances. In the next section, we shall simplify this formulation making it independent of time at the price of losing some of its real-world characteristics. Because of that, the following
step consists on reducing the dimension of z variables and simplifying useless variables.

§.§ Non-Time Dependent Formulations

The simplification mentioned above can be performed based on the rationale of ensuring connectivity through different sets of inequalities. In particular, we compare Miller-Tucker-Zemlin (MTZ) inequalities and  subtour elimination constraints (SEC).
All formulations use the following sets of decision variables:

* Binary variables $z_e\in\{0,1\},\,e\in E$, to represent the edges of the tours.

* Continuous variables $d_e\geq 0,\,e=\{v, w\}\in E_\text{out}\subseteq E'$, to
represent the distance $d(x^1_v, x^2_w)$ between the pairs of selected points
of different elements (neighborhoods) and $d_v\geq 0,\,v\in V$, to represent the distance
$d(x^1_v, x^2_v)$ between the pairs of points of the same element.

$$\mathcal D_e=\{d\in\mathbb R_+^{|E_\text{out}|}:d_e\geq d(x_v^1, x_w^2), \forall e=(v, w)\in E_\text{out},\,x\in\mathcal X\},$$

$$\mathcal D_v=\{d\in\mathbb R_+^{|E_\text{in}|}:d_v\geq d(x_v^1, x_v^2), \forall v\in V,\,x\in\mathcal X\},$$

denote the domains for the feasibility of the $d$ variables. Then, a generic
bilinear formulation for XPPN is

∑_e∈E_out d_e z_e + ∑_e∈E_in f_v d_vPdz
z∈𝒯_G,d_e∈𝒟_e, d_v∈𝒟_v

(<ref>), (<ref>), (<ref>)

The reader should observe that, as already mentioned, the above formulation is bilinear since the first term of the objective function contains products of variables of the form $d_ez_e$, for $e\in E_\text{out}$.

Next, we use McCormick's envelopes [14] for the linearization of those bilinear terms of the objective function. We define additional variables $p_e\geq 0,\,e\in E_\text{out}$ that stand for that product.

Replacing the products by the new variables and introducing a new set of constraints enforcing the correct representation, we obtain the following formulation:
P = ∑_e∈E_out p_e + ∑_v∈V f_v d_vRL-XPPN
p_e≥d_e - M_e(1 - z_e)  ∀e∈E_out
(<ref>), (<ref>), (<ref>)
z∈𝒯_G,  d_e∈𝒟_e, d_v∈𝒟_v  p_e≥0,  ∀e∈E_out.

Here $M_e$ denotes an upper bound of the distance between the sets that are joined by $e$.

In addition, we describe the sets $\mathcal D_e$ and $\mathcal D_v$ using the constraints

\begin{align}
\|x_v^1 - x_w^2\|_2 & \leq d_e,\quad\forall e=\{v, w\}\in E_\text{out}, \label{D1}\tag*{(D$_1$)} \\
\|x_v^1 - x_v^2\|_2 & \leq d_v,\quad\forall v\in V, \label{D2}\tag*{(D$_2$)}                     \\
x                   & \in  \mathcal X, \tag*{(D$_3$)}
\end{align}

which set the distance values and impose that $x$ belongs to its suitable neighborhood.

Furthermore, this formulation can be reinforced by adding the valid inequalities:
$p_e\geq m_e z_e,\forall e\in E_\text{out}$ and $d_v\leq M_v,\forall v\in V$, where $m_e$ and $M_v$ are bounds that are adjusted in Section <ref>. The first family of valid inequalities sets lower bounds on the values for $p_e$ whereas the second ones sets upper bounds on the distances traveled within neighborhoods.

The above discussion leads us to strengthen a generic formulation for XPPN. This formulation will be particularized once the connectivity condition of the solutions is specifically introduced in the model.

P = ∑_e∈E_out p_e + ∑_v∈V f_v d_v
p_e≥d_e - M_e(1 - z_e)  ∀e∈E_outLIN-Mc
p_e≥m_e z_e  ∀e∈E_out
d_v≤M_v  ∀v∈V
(<ref>), (<ref>), (<ref>)

The two formulations that we present below differ from one another in the family of constraints used to enforce connectivity. One of them is by the family of subtour elimination constraints (SEC) [6]. The other one relies on a compact formulation based on the well-known Miller-Tucker-Zemlin (MTZ) constraints [15].

§.§.§ A valid formulation for XPPN based on SECs.

The family of SEC is well-known in combinatorial optimization. It enforces connectivity by imposing that the number of edges among any subset of vertices can not exceed its cardinality minus one. Augmenting these constraints into the generic formulation presented above we obtain the following valid formulation for XPPN:

P = ∑_e∈E_out p_e + ∑_v∈V f_v d_vSEC-XPPN
(<ref>), (<ref>), (<ref>)
∑_w∈V∖{v} z_vw=1,  ∀v ∈VC$_1$
∑_w∈V∖{v} z_wv=1,  ∀v ∈V
∑_e=(v, w):v,w∈S z_e≤|S| - 1,  ∀S⊂VSEC
(<ref>), (<ref>), (<ref>)

Assignment Constraints (<ref>) and (<ref>) ensure that any feasible solution found enters and exits each component of the problem exactly once. Constraint (<ref>) prevents the existence of subtours. This constraint forces that in any subset of nodes $S$ included in $V$ there can not be more edges between nodes in $S$ than the number of nodes that contains minus one, thus avoiding the existence of cycles.

Since there is an exponential number of SEC constraints, when we implement this formulation we need to perform a row generation procedure including constraints whenever they are required by a separation oracle. To find SEC inequalities, as usual, we search for disconnected components in the current solution. Among them, we choose the shortest subtour found in the solution to be added as a lazy constraint to the model.

If the considered distance between components is symmetric, we obtain the symmetric formulation based on SECs, denoted by (sSEC-XPPN). In this formulation, we can halve the number of binary variables and replace constraints (<ref>) and (<ref>) in (<ref>) by the following connectivity restrictions:

\begin{equation*}
\sum_{w\in V\setminus\{v\}}z_{wv} = 2,\qquad \forall v\in V.
\end{equation*}

§.§.§ XPPN formulation based on the Miller Tucker Zemlin.

This section addresses an alternative formulation that results replacing SEC inequalities by the so called Miller-Tucker-Zemlin constraints [15]. In this formulation, we introduce the integer variable $s_v$ to generate an alternative formulation that eliminates the subtours and the exponential number of inequalities of (<ref>).
P = ∑_e∈E_out p_e + ∑_v∈V f_v d_vMTZ-XPPN
(<ref>), (<ref>), (<ref>)
∑_w∈V∖{v} z_vw=1,  ∀v ∈VC$_1$
∑_w∈V∖{v} z_wv=1,  ∀v ∈VC$_2$
|V|z_vw + s_v - s_w≤|V| - 1,  ∀e=(v,w)∈E_outMTZ$_1$
s_1= 1MTZ$_2$
2≤s_v≤|V|,   ∀v∈VMTZ$_3$
s_v-s_w + |V|z_wv≤|V|-1,  ∀e=(v,w)∈E_out, w>1MTZ$_4$
s_v-s_w + (|V|-2)z_wv≤|V|-1,  ∀e=(v,w)∈E_out, v>1MTZ$_5$
(<ref>), (<ref>), (<ref>)

Again constraints (<ref>) and (<ref>) require that in each feasible solution only one edge departs from node $v$ and only one edge enters at node $v$ for any $v\in V$, respectively. The constraints (<ref>)-(<ref>) enforce connectivity, i.e., that there is only a single tour covering all vertices. The constraints (<ref>) and (<ref>) define the intermediate conditions for the  tour that may improve performance of subtour elimination constraints (see [20] for more details).

For the sake of completeness, we include as a remark a well-known explanation of these last constraints [15].

Constraints (<ref>)-(<ref>) model the elimination of subtours.

To show that, we must see that every feasible solution contains only one closed sequence of cities and that for every single tour covering all cities, there are values for the dummy variables $s_{v}$ that satisfy the constraints.

To prove that every feasible solution contains only one closed sequence of cities, it suffices to show that every subtour in a feasible solution passes through city 1. If we sum all the inequalities corresponding to $z_{vw}=1$ for any subtour of $k$ steps, we obtain:

$$|V|k\leq (|V|-1)k, $$

which is not possible.

It now must be shown that for every single tour covering all cities, there are values for the dummy variables $u_{i}$ that satisfy the constraints. If we define $s_{v}=t$ if city $v$ is visited in step $t$. Then

$$|V|z_{vw} + |V| - 2 \leq |V| z_{vw} + s_v - s_w \leq n-1,$$

forces $z_{vw}=0$. For $z_{vw}=1$, we have:

$$s_{v}-s_{w}+|V|z_{vw}= t-(t+1)+|V|=|V|-1,$$

satisfying the constraint.

Now we state a result related to the relationship between the SEC and MTZ polytopes of our formulations of the XPPN, that is, the feasible regions of the respective LP relaxations of these models.

The SEC polytope is contained in the MTZ polytope for the XPPN.

Notice that the only difference between these polytopes is the constraint that ensures the elimination of subtours. Therefore, it is enough to see that the (<ref>) constraints are stronger than those given in (<ref>)-(<ref>), that is proved in [21].

§ A HEURISTIC ALGORITHM FOR XPPN

In this section we present a heuristic algorithm for solving XPPN. This algorithm has two different applications. On the one hand, it provides good quality feasible solutions for XPPN that become a promising alternative to exact methods whenever the size of the problems is large. On the other hand, it also helps in solving exactly XPPN by feeding the exact formulations with a good initial solution which in turns speeds up the branch and bound search. The considered algorithm is composed by two phases: the Clustering Phase and the Variable Neighborhood Search (VNS) Phase. The so called clustering phase determines some points in each dimensional element (polygonal chain or neighborhood) and then the VNS phase finds a heuristic tour on the complete graph spanned by the previously obtained points.

The Clustering Phase

The first phase of the heuristic algorithm is based on solving a relatively easy single facility location problem: the Weber or median Problem. The solution of this problem looks for a prototype point (a representative) of the dimensional elements in the problem (neighborhoods and polygonal chains) that minimizes the sum of the Euclidean distance from all those elements, i.e., we want to find a point $C$ and some points $x_v$ in the element associated with the vertex $v\in V$ such that:

∑_v∈Vx_v - C
(<ref>), (<ref>), (<ref>)

The idea of this approach is to find some points that are likely to be close to the true chosen  points in each element in the final optimal route. Figure <ref> shows an example that combines three neighborhoods and seven polygonal chains. Red points represent the points of each set and the green point is the proposed 1-median obtained after solving the corresponding Weber problem described above.

Illustration of the first phase of the heuristic algorithm

The Variable Neighborhood Search Phase

Once the points of each set have been chosen,  the idea is to find the minimal cost route that joins these points. To obtain this route, we have used the well-known Variable Neighborhood Search metaheuristic developed in [16]. The code has been taken from [19].

Using the example depicted in Figure <ref>, we generate a tour considering this VNS approach with a maximum number of 25 attempts, a neighborhood of size 5 and 10 iterations. The final result is shown in Figure <ref>.

Application of the VNS phase to the example of Figure <ref>

Finally, in order to build a feasible solution for XPPN we take into account the position of the points (represented by $x_v^1$ and $x_v^2$) and the order in which they have to be visited in the  tour obtained by VNS phase of our heuristic (represented by $z_{vw}$). Once the solution is built, it can also be taken as an initial solution for any of the exact formulations presented above.

§ STRENGTHENING THE FORMULATION OF XPPN

§.§ Pre-processing

In this section we explore the geometry of the neighborhoods that appear in the problem to fix a priori some variables and to increase the efficiency of the model.

Firstly, we consider two special cases that relate the position of the entry and exit points of each neighborhood with the coefficient $f_v$ of the objective function.

If the problem verifies that $f _v= 0$ for all $v\in V _{\mathcal C}$, then the entry and exit points $x_v^1$ and $x_v^2$ selected in each neighborhood are the same that the ones obtained by minimizing the distance between the neighborhoods. (No sé cómo escribirlo, lo hablamos a ver cómo se podría poner)

If the problem verifies that $f_v\geq 1$ for all $v\in V _{\mathcal C}$, then, in the optimum, $x_v^1=x_v^2$.

Let $x_u^2$ and $x_w^1$, the exit (resp. entry) point of the neighborhood $u\in V_{\mathcal C}$ (resp. $w\in V_{\mathcal C}$) and suppose that in the optimum $x_v^1 \neq x_v^2$, i.e., $d(x_v^1, x_v^2)>0$. Let $x'_v=\frac{x_v^1+x_{v}^2}{2}$ the middle point of these points. If we compare the length of the new built path $p'=[x_u^2,x_v', x_w^1]$ with the path $p=[x_u^2, x_v^1, x_v^2, x_w^1]$, we obtain by using the triangular inequality that

\begin{align*}
length(p') &= d(x_u^2, x_v') + d(x_v', x_w^1) \\
&\leq d(x_u^2, x_v^1) + f_vd(x_v^1, x_v')+f_vd(x_v',x_v^2)+d(x_v^2, x_w^1)\\
&= d(x_u^2, x_v^1) + f_vd(x_v^1, x_v^2)+d(x_v^2, x_w^1)\\
\end{align*}

This fact contradicts that $p'$ is the optimum.

From now on, we assume in the rest of this section that $f_v\geq 1$ for all $v\in V$. The following outcome restricts the domain where the selected point are located.

Any point selected in an optimal solution of the XPPN must be place in the boundary of the neighborhoods.

By induction in the number of neighborhoods.

* $n=2:$ it is known that the minimum distance between two convex sets occurs in the boundary of these sets.

* $n\Rightarrow n + 1:$ let assume that the problem has $n+1$ sets and there exists a neighborhood $N_k$ whose selected point in the optimum is in $int(N_k)$. Let $T$ be the triangle formed by the point $A$ of the previous neighborhood, the point $B$ in the neighborhood $k$ and the point $C$ of the next neighborhood in the optimum edge sequence.  Let also $\overline{AC}$ the line segment that joins $A$ and $C$. The idea of the proof is to construct a point in the border of $N_k$ whose distance to $A$ and $C$ is lower. We have two cases:

* If $A, B, C$ are aligned, we consider the point $B'=\overline{AC}\cap\partial N_k$. Since $B'$ is also aligned with $A$ and $C$, $\overline{AB'} + \overline{B'C}=\overline{AB} + \overline{BC}$.

* If $A, B, C$ are not aligned, we consider the point $B'=r^{\bot}\cap \partial N_k$, where $r^\bot$ is the perpendicular line to the line segment $\overline{AC}$ and $\partial N_k$ is the boundary of $N_k$. Note that this intersection produces two points in $\partial N_k$, we take the closest one to $\overline{AC}$. If we call $T'$ the triangle generated by $A$, $B'$ and $C$, then the height of $T'$ to $\overline{AC}$ is lower than the one of $T$, hence by the Pythagorean Theorem, $\overline{AB'} < \overline{AB}$ and $\overline{B'C} < \overline{BC}$, which is a contradiction.

The special case in which all the neighborhoods are circles, allows us to limit more the location of the points based on the construction given in the Proposition <ref>.

Any point selected in an optimal solution of the XPPN when all the neighborhoods are circles is placed in some arc of one of the circumferences inside of the convex hull generated by the center of the circles.

By induction in the number of circles.

* $n=2:$ the points are located in the line segment that joins the center of the circles.

* $n\Rightarrow n + 1:$ let assume that the problem has $n+1$ circles and there exists a circle $N_k$ whose selected point in the optimum is not in the convex hull of the center of the neighborhoods. Let $T$ be the triangle formed by the point $A$ of the previous neighborhood, the point $B$ in the neighborhood $k$ and the point $C$ of the next neighborhood in the optimum edge sequence.  Let also $\overline{AC}$  be the line segment that joins $A$ and $C$. The idea of the proof is to construct a point in the border of the convex hull $\mathcal C$ whose distance to $A$ and $C$ is lower. We have two cases depending on the location of $B$ in the neighborhood:

* If $A, B, C$ are aligned, we consider the point $B'=\overline{AC}\cap\partial N_k$. Since $B'$ is also aligned with $A$ and $C$, $\overline{AB'} + \overline{B'C}=\overline{AB} + \overline{BC}$.

* If $A, B, C$ are not aligned, we split the circumference on two arcs $A_R$ and $A_B$. These arcs are built by taking the perpendicular line to the edge of the convex hull:

* If $B\in A_R$, we take $B'$ the projection to the convex hull and it produces a triangle $T'$ with lower height to $\overline{AC}$. Then, we can use the Proposition <ref> to construct a point in the boundary of $N_k$ that lies in the convex hull. (See Figure )

* If $B\in A_B$, we construct $B'$ the diametrically opposite point of $B$ in $N_k$ and it also produces a smaller height that contradicts the assumption that $B$ gives the shortest tour.

Now we give a result to eliminate some neighborhoods and simplify the problem without modifying the objective value.

Given two neighborhoods $A$ and $B$, if $B\supset A$, then $B$ can be deleted in the problem.

Starting from the optimal solution of problem without $B$, we are going to build an optimal solution including $B$ that is essentially the same. Let $z^*$ the optimal tour by deleting the neighborhood $B$ in the problem. By connectivity, there exist two neighborhoods $A_{-1}$ and $A_{+1}$ that are connected with $A$, i.e., such that $z^*_{A_{-1}A}=z^* _{AA_{+1}}= 1$. In addition, let $x^*_A$ be the point chosen to visit the neighborhood $A$. If we include $B\supset A$ in the problem and we fix $x_B=x^*_A$ and $z_{AB}=z_{BA_{+1}}=1$. This solution is also a Hamiltonian path whose objective value is the same because $d(A, B) =0$ and $d(B, A_{+1})=d^*(A, A_{+1})$.

§.§ Valid inequalities

The different models that we have proposed include in one way or another big-M constants. We have defined different big-M constants along this work. In order to strengthen the formulations we provide good upper bounds for those constants. In this section we present some results that adjust them for each kind of set considered in our models.

The first bound we need to adjust is $M_e$ that denotes an upper bound of the distance between the sets joined by the edge $e$. We have three cases that depend of the shape of the sets $A$ and $B$:

* If $A$ and $B$ are both ellipsoids, we cannot easily compute the maximum distance between $A$ and $B$, but we can generate an upper bound of this distance by taking diametrically opposite points of minimum radius circles that contain each entire ellipsoid. When both ellipsoids are circles, this bound coincides with the maximum distance.

* If $A$ is an ellipsoid and $B$ is a polygon or a polygonal chain, we can set this bound by the maximum of the distances of each vertex of $B$ to the center of $A$ plus the radius of the minimum circle that contains the ellipsoid $A$.

* If $A$ and $B$ are both polygons or polygonal chains, this bound can be computed exactly by taking the maximum of the distances between vertices of $A$ and $B$.

The second bound to be adjusted is $m_e$. It denotes a lower bound of the distance of the sets joined by the edge $e$. In this case, we can compute this distance exactly by solving a convex program the minimizes the distance of the sets $A$ and $B$.

In the Figures <ref>, <ref> and <ref> we show the selected maximal (red) and minimal (blue) bounds depending on the shape of the sets:

Upper and lower bound when both sets are ellipsoids

Upper and lower bound when a set is a polygon and the other is an ellipsoid

Upper and lower bound when both sets are polygons

In addition, the third bound represent the maximal distance between two points within a given neighborhood.  We can compute this upper bound according to the shape of this set (see Figure <ref>):

* If the set is an ellipsoid, we can take diametrically opposite points of the minimum radius circle that contains this ellipsoid.

* If the set is a polygon, we can compute the maximum of the distances between each pair of vertices.

* If the set is a polygonal chain, this bound equals the length of the polygonal.

Upper bound on the maximal  distance within a set

§ A DECOMPOSITION ALGORITHM

In this section we present an alternative row generation approach to solve the XPPN based on a Benders decomposition of the problem. The general method is based on the following observation: If we fix $z\in\mathcal T_G$ in the generic formulation of XPPN, we obtain a continuous SOC problem, which is well-known to be convex. On the other hand, the
objective function that we are considering is bilinear. Hence, we can use a Benders-like decomposition approach (see [2]) to generate
an iterative algorithm that solves this problem.

For a given $\bar{z} \in \mathcal T_G$, the “optimal” vertices and distances of its associated XPPN can be computed by solving the following subproblem:

d(z̅)=∑_e∈E_out d_e z̅_e + ∑_v∈V f_v d_vPd$\bar z$
d_e∈𝒟_e, d_v∈𝒟_v.

Note that the number of $d$ variables in (<ref>) is $2n$, because only distances with nonzero $\bar z_e$ variables need to be calculated. Thus, Benders decomposition is
a good approach for the  XPPN model (see [3]). The explicit form of the Benders cuts is the following:

\begin{equation}
P \geq d(\bar z) + \sum_{e:\bar z_e = 1}M_e(z_e - 1) + \sum_{e:\bar z_e = 0}m_e z_e
\end{equation}

where $P = \sum_{e\in E_\text{out}} p_e + \sum_{v\in V} f_v d_v$ with $p_e\geq 0$ and $M_e$ and $m_e$ are the upper and lower bounds estimated in the above section.

Therefore, the relaxed master problem at the $K$-th iteration of the  row-generation solution algorithm can be stated as:

\begin{align}
P^* = \min \qquad &  P \nonumber\\
& \hspace*{-1cm}  P \geq  d(\bar z^k) + \sum_{e:\bar z_e^k = 1}M_e(z_e^k - 1) + \sum_{e:\bar z_e^k = 0}m_e z_e^k,\; k=1, \ldots, K,\label{xmatheu000}\\
& P = \sum_{e\in E_\text{out}} p_e + \sum_{v\in V} f_v d_v \\
& z \in \mathcal{T}_G.\nonumber
\end{align}

Adding  the above cuts sequentally gives rise to the solution scheme described in Algorithm <ref>:

Let $z^0\in \mathcal{T}_G$ be an initial  solution and $\varepsilon$ a given threshold value.

Set $LB=0$, $UB=+\infty$, $\bar z=z^0$.

* Solve (<ref>) for $\overline z$ to get $d(\overline{z})$.

* Add the cut $ P \geq  d(\bar z) + \sum_{e:\bar z_e = 1}M_e(z_e - 1) + \sum_{e:\bar z_e = 0}m_e z_e$ to the current master problem.

* Obtain the optimal value $\bar{P}$ to the current master problem, and its associated solution $\bar z$.

* Update $LB=\max\{LB, \bar{P}\}$ and $UB=\min\{UB, \sum_{e\in E} d(\overline{z})_e\overline{z}_e +\sum_{v\in V} f_v d_v\}$

Decomposition Algorithm for solving XPPN.

The stopping criterion is that the gap between the upper and lower bound does not exceed the fixed threshold value $\varepsilon$.

Theorem 2.4 in [9] states the finite convergence of the decomposition approach under the following assumptions: convexity and finiteness of the feasible domains, closeness of the “linking” constraints between the sets, and convexity of the objective functions. In our case, the finiteness of the number of underlying Hamiltonian tour of $\mathcal T_G$, the convexity of (<ref>) for any $\overline z \in \mathcal T_G$,  and the linear separability of the problem allows to apply the above result, which assures that Algorithm <ref> terminates in a finite number of steps (for any given $\varepsilon\geq 0$).

To avoid the enumeration of all Hamiltonian tours of $\mathcal T_G$, we have started with a non-empty set of cuts which  give a suitable initial representation of the lower envelope of $P$.

Given that the master problem exhibits a combinatorial nature, we have embedded the cut generation mechanism within a branch-and-cut scheme.

The computational results obtained by our implementation are included in  table <ref>. This table compares the results of the Final Gap, cpu time and number of cuts added applying the  decomposition algorithm versus those obtained with formulation MTZ. From these results we conclude that the decomposition algorithm performs worse than formulation MTZ even for small size problems. To reinforce our conclusion, we have also included a performance profile of number of solved instances versus time for formulations sSEC, MTZ and the decomposition algorithm (see Figure <ref>). The reader can observe that the number of solved instances within the time limit is approximately one half comparing Benders decomposition with MTZ and sSEC. These results lead us not to include this algorithm in the final computational experience for larger problem sizes presented in the next section.

Size     Radii     Mode     Final Gap (Benders)     Time (Benders)     #Cuts     Final Gap (MTZ)     Time (MTZ)
10     1     1     0.0     15.72     19.0     0.0     1.93
10     1     2     0.0     23.64     55.4     0.0     0.75
10     1     3     0.0     13.22     25.8     0.0     0.72
10     1     4     0.0     33.96     29.8     0.0     1.52
10     2     1     76.1     6430.08     1209.0     0.0     38.83
10     2     2     56.14     4777.58     1009.6     0.0     14.14
10     2     3     0.0     1766.06     380.6     0.0     2.23
10     2     4     10.57     5993.57     804.6     0.0     2.52
10     3     1     96.21     7208.56     1481.2     0.0     487.94
10     3     2     92.16     7203.99     1352.2     0.0     35.81
10     3     3     9.29     5832.51     520.4     0.0     13.28
10     3     4     84.41     7214.86     921.8     0.0     133.81
10     4     1     98.79     7205.35     2283.0     19.28     3513.38
10     4     2     95.53     7207.51     1343.4     0.0     238.98
10     4     3     19.55     7220.14     499.0     0.0     20.25
10     4     4     82.69     7211.46     789.8     0.0     1142.17

Computational comparison between MTZ formulation and Benders algorithm for problems with up to 10 neighbors

Performance profile: Time vs #Solved

§ COMPUTATIONAL EXPERIMENTS

In this section we have performed a series of experiments to compare the non-time dependent formulations presented in Section <ref>. Based on the work of Blanco, Fernández and Puerto (see [3] for more details) we have generated five instances of size $m\in[5, 10]$ and we report the average results. We have considered three different types of neighbors to be visited:

* Circles of radii $r$.

* Regular polygons of radii $r$ with a random number of vertices $|V|\in[3, 10]$.

* Polygonal chains parameterized by its breakpoints that are distanced in $r$ and some random percentage $\alpha\in [0, 1]$.

In addition, the centers or breakpoints of these elements have been generated uniformly in the square [0, 100]. On the one hand, we have studied four different scenarios to generate the radii to define the elements:

* Small size Neighborhoods $(r=1)$: Radii randomly generated in $[0, 5]$.

* Small-Medium Neighborhoods $(r=2)$: Radii randomly generated in $[5, 10]$.

* Medium-Large size Neighborhoods $(r=3)$: Radii randomly generated in $[10, 15]$.

* Large size Neighborhoods $(r=4)$: Radii randomly generated in $[15, 20]$.

On the other hand, we have also considered four modes depending on the nature of the neighborhoods:

* Mode 1: All neighborhoods are circles.

* Mode 2: All neighborhoods are regular polygons.

* Mode 3: All neighborhoods are polygonal chains.

* Mode 4: Mixture of the three previously  considered neighborhoods.

All the formulations were coded in Python 3.7, and solved using Gurobi 8.1.0. [10] in a Intel(R) Xeon(R) E-2146G CPU @ 3.50 GHz and 64GB of RAM. A time limit of 2 hours was set in all the experiments.

Our preliminary test is devoted to decide whether initializing the solution process with an initial solution helps is solving the problem or not. In this regard, we have performed a comparison between formulation MTZ with and without the initial solution provided by the VNS heuristic. The results are summarized in Table <ref>. This table reports average results for instances of sizes 5,10 and 15 neighbors with all combinations of radii and modes. It compares the final gap and running times for the formulation with initial solution (Final Gap (Init), Opt. Time (Init)) and without (Final Gap (NoInit), Opt. Time (NoInit)). The results are also depicted in the boxplox diagrams in Figure <ref>. Both, table and figure, clearly show that loading an initial solution helps in reducing the gap and the cpu time: all the instances up to 10 neighbors are solved to optimality with and without initial solution but for 15 neighbors the final gap in the first case is always better than in the latter (blue boxes are always below orange ones). Based on this results, in the following, we have always solved the instances loading an initial solution.

Size     Radii     Mode     Final Gap (Init)     Final Gap (NoInit)     Opt. Time (Init)     Opt. Time (NoInit)
5     1     1     0.0     0.0     0.61     0.45
5     1     2     0.0     0.0     0.12     0.08
5     1     3     0.0     0.0     0.21     0.13
5     1     4     0.0     0.0     0.25     0.27
5     2     1     0.0     0.01     0.27     0.4
5     2     2     0.0     0.0     0.13     0.11
5     2     3     0.0     0.0     0.17     0.14
5     2     4     0.0     0.0     0.17     0.19
5     3     1     0.0     0.0     0.44     0.42
5     3     2     0.0     0.01     0.16     0.12
5     3     3     0.0     0.0     0.17     0.16
5     3     4     0.0     0.0     0.3     0.35
5     4     1     0.0     0.01     0.34     0.4
5     4     2     0.0     0.0     0.14     0.34
5     4     3     0.0     0.0     0.16     0.16
5     4     4     0.0     0.0     0.28     0.3
10     1     1     0.0     0.0     1.93     4.53
10     1     2     0.0     0.0     0.75     0.84
10     1     3     0.0     0.0     0.72     1.77
10     1     4     0.0     0.0     1.52     2.95
10     2     1     0.0     0.0     38.83     61.53
10     2     2     0.0     0.0     14.14     44.93
10     2     3     0.0     0.0     2.23     4.65
10     2     4     0.0     0.0     2.52     7.5
10     3     1     0.0     1.09     487.94     1049.86
10     3     2     0.0     0.0     35.81     153.37
10     3     3     0.0     0.0     13.28     29.43
10     3     4     0.0     0.0     133.81     510.58
10     4     1     19.28     10.0     3513.38     4134.31
10     4     2     0.0     0.0     238.98     1253.21
10     4     3     0.0     0.0     20.25     82.39
10     4     4     0.0     11.75     1142.17     3490.88
15     1     1     0.0     0.0     18.58     196.71
15     1     2     0.0     0.0     3.38     23.37
15     1     3     0.0     0.0     314.57     44.56
15     1     4     0.0     0.0     10.94     90.1
15     2     1     6.69     23.17     3135.29     7200.72
15     2     2     0.0     14.06     2460.78     7200.49
15     2     3     0.0     0.0     14.58     49.51
15     2     4     0.0     5.88     1052.16     4660.88
15     3     1     46.33     59.74     5760.56     7200.28
15     3     2     20.79     31.2     5760.84     7200.8
15     3     3     0.0     0.0     322.77     599.25
15     3     4     14.07     19.17     5865.82     6896.78
15     4     1     100.0     100.0     7200.47     7200.98
15     4     2     20.2     36.9     4421.25     7200.51
15     4     3     0.19     0.72     2195.2     3566.82
15     4     4     21.6     27.71     7200.42     7200.5

Computational comparison between MTZ formulation with and without initial solution

Comparison of the final gap between MTZ formulation with and without initial solution after 7200 seconds

The remaining information of our computational experiments can be found in tables <ref>, <ref> and <ref>. The first one reports our results for formulation SEC, the second one for sSEC (symmetric version of SEC) and the third one for MTZ. Information in all the three tables is organized in the same way. Each row shows averages of five instances for different combinations of factors (Size, Radii and Mode) each one with four different levels. Our tables have 9 columns. The first three (Size, Radii and Mode) describe the parameters of the problem. Then, we report the final gap (% Final Gap), time required by the exact method (Exact Time), time required by the heuristic (Heur. Time) and the % improvement of the gap with respect to the initial solution (% Improved Gap).

Size     Radii     Mode     Final Gap     Exact Time     Heur. Time     % Improved Gap
5     1     1     0.0     0.11     2.29     0.83
5     1     2     0.0     0.06     2.2     0.57
5     1     3     0.0     0.14     3.82     0.37
5     1     4     0.0     0.12     2.8     0.92
5     2     1     0.0     0.09     2.35     2.51
5     2     2     0.0     0.06     2.37     2.25
5     2     3     0.0     0.15     3.51     1.37
5     2     4     0.0     0.09     2.68     2.53
5     3     1     0.0     0.11     2.26     3.79
5     3     2     0.0     0.08     2.34     3.85
5     3     3     0.0     0.16     3.72     2.05
5     3     4     0.0     0.15     2.82     2.29
5     4     1     0.0     0.13     2.31     4.42
5     4     2     0.0     0.26     2.18     5.18
5     4     3     0.0     0.16     3.48     3.69
5     4     4     0.0     0.14     2.89     8.05
10     1     1     0.0     0.95     4.31     3.25
10     1     2     0.0     0.47     4.37     2.56
10     1     3     0.0     1.72     7.72     1.28
10     1     4     0.0     4.81     5.6     1.53
10     2     1     0.0     21.74     4.73     6.48
10     2     2     0.0     8.45     4.36     6.37
10     2     3     3.23     2949.61     9.25     2.21
10     2     4     1.38     2188.88     6.14     3.22
10     3     1     0.0     522.23     5.08     10.0
10     3     2     0.0     84.64     4.9     9.26
10     3     3     3.47     4324.06     10.54     1.69
10     3     4     0.0     404.51     5.29     7.48
10     4     1     5.32     2484.47     5.17     7.76
10     4     2     0.0     539.03     4.85     9.21
10     4     3     3.57     3656.07     10.11     5.47
10     4     4     16.69     6548.93     6.18     7.86
15     1     1     0.0     14.12     5.63     2.74
15     1     2     0.0     4.63     5.58     4.59
15     1     3     12.12     7200.55     12.9     0.6
15     1     4     0.46     1597.94     7.44     4.2
15     2     1     29.42     7200.43     5.7     11.07
15     2     2     17.89     6178.28     5.6     8.74
15     2     3     13.91     7200.95     12.55     0.1
15     2     4     23.44     7200.6     7.25     3.11
15     3     1     70.59     7200.52     5.77     12.59
15     3     2     35.67     7200.65     5.78     12.85
15     3     3     9.7     5828.43     12.44     2.16
15     3     4     45.94     7200.79     9.95     0.49
15     4     1     0.0     10.68     5.65     11.74
15     4     2     43.12     7200.7     5.67     7.02
15     4     3     7.6     5789.95     12.67     3.58
15     4     4     41.1     7200.6     8.48     6.57
20     1     1     2.58     3189.58     6.9     2.72
20     1     2     1.78     2894.34     6.47     4.59
20     1     3     10.17     5797.93     13.57     1.78
20     1     4     11.01     6434.25     11.34     1.47
20     2     1     63.7     7200.95     6.41     10.96
20     2     2     39.3     7201.34     6.77     10.83
20     2     3     11.82     6050.34     14.24     4.17
20     2     4     37.99     7200.84     10.26     2.74
20     3     1     95.47     7200.71     6.7     15.29
20     3     2     55.89     7201.07     6.62     16.12
20     3     3     17.75     7201.0     14.43     2.0
20     3     4     45.88     7200.79     11.44     11.6
20     4     1     0.0     855.97     6.72     10.99
20     4     2     60.12     7201.0     6.43     1.85
20     4     3     10.06     7201.0     14.34     4.05

Asymmetric SEC results with initial solution

Size     Radii     Mode     Final Gap     Exact Time     Heur. Time     % Improved Gap
5     1     1     0.0     0.08     2.38     0.0
5     1     2     0.0     0.04     2.29     0.0
5     1     3     0.0     0.08     3.74     0.51
5     1     4     0.0     0.05     2.88     0.09
5     2     1     0.0     0.06     2.35     0.01
5     2     2     0.0     0.06     2.38     0.0
5     2     3     0.0     0.08     3.59     1.5
5     2     4     0.0     0.06     2.72     1.26
5     3     1     0.0     0.05     2.28     0.01
5     3     2     0.0     0.05     2.36     2.47
5     3     3     0.0     0.09     3.67     0.81
5     3     4     0.0     0.07     2.83     3.36
5     4     1     0.0     0.06     2.32     0.01
5     4     2     0.0     0.05     2.31     0.13
5     4     3     0.0     0.07     3.62     5.95
5     4     4     0.0     0.08     2.95     1.9
10     1     1     0.0     0.29     4.52     0.01
10     1     2     0.0     0.14     4.66     0.01
10     1     3     0.0     0.2     7.8     0.0
10     1     4     0.0     0.23     5.51     0.04
10     2     1     0.0     3.72     5.41     0.05
10     2     2     0.0     1.88     5.69     2.28
10     2     3     0.0     1.34     9.34     0.0
10     2     4     0.0     1.36     9.68     1.17
10     3     1     0.0     45.54     5.27     0.15
10     3     2     0.0     9.37     4.9     0.53
10     3     3     0.0     490.19     10.53     0.67
10     3     4     0.0     7.78     5.33     4.51
10     4     1     0.0     520.74     5.22     0.26
10     4     2     0.0     36.72     4.99     8.1
10     4     3     0.54     1474.78     11.41     0.0
10     4     4     0.0     1461.93     6.84     2.91
15     1     1     0.0     2.1     6.83     0.0
15     1     2     0.0     0.86     6.46     1.11
15     1     3     3.14     2881.8     13.26     0.0
15     1     4     0.0     1.17     7.92     0.0
15     2     1     12.93     5840.0     7.19     0.0
15     2     2     8.9     4560.43     7.0     0.75
15     2     3     11.18     5761.15     15.14     0.11
15     2     4     8.3     4642.84     8.36     0.0
15     3     1     64.34     7200.42     7.39     0.0
15     3     2     28.89     7200.44     7.45     0.0
15     3     3     5.59     5765.79     16.7     0.42
15     3     4     24.3     7200.44     10.87     2.94
15     4     1     0.0     271.8     7.55     0.0
15     4     2     35.34     7200.52     7.41     3.21
15     4     3     12.84     7200.49     248.17     0.96
15     4     4     35.59     7200.53     10.82     0.8
20     1     1     0.0     175.65     11.47     0.81
20     1     2     0.95     1632.15     10.98     0.03
20     1     3     0.0     466.74     18.07     1.28
20     1     4     8.59     2887.48     16.38     1.36
20     2     1     43.77     7200.49     11.95     0.0
20     2     2     26.72     7200.65     11.51     0.0
20     2     3     4.65     4354.1     27.26     0.0
20     2     4     26.11     7200.5     16.42     0.37
20     3     1     81.51     7200.5     12.52     0.0
20     3     2     47.27     7200.95     12.13     0.0
20     3     3     16.84     7200.81     37.43     0.26
20     3     4     44.18     7200.59     19.15     0.0
20     4     1     0.0     483.78     13.51     0.0
20     4     2     55.27     7200.73     12.43     0.0
20     4     3     15.84     7200.84     3191.77     0.28
20     4     4     40.08     7200.79     20.47     0.0

Symmetric SEC results with initial solution

Size     Radii     Mode     Final Gap     Exact Time     Heur. Time     % Improved Gap
5     1     1     0.0     0.61     1.47     0.02
5     1     2     0.0     0.12     1.3     0.0
5     1     3     0.0     0.21     2.01     0.08
5     1     4     0.0     0.25     1.64     0.06
5     2     1     0.0     0.27     1.33     2.38
5     2     2     0.0     0.13     1.25     0.01
5     2     3     0.0     0.17     1.9     0.64
5     2     4     0.0     0.17     1.49     2.22
5     3     1     0.0     0.44     1.28     1.01
5     3     2     0.0     0.16     1.28     4.02
5     3     3     0.0     0.17     1.95     0.73
5     3     4     0.0     0.3     1.68     0.44
5     4     1     0.0     0.34     1.28     15.36
5     4     2     0.0     0.14     1.26     8.52
5     4     3     0.0     0.16     1.98     1.55
5     4     4     0.0     0.28     1.7     4.39
10     1     1     0.0     1.93     2.63     0.09
10     1     2     0.0     0.75     2.3     0.01
10     1     3     0.0     0.72     4.49     0.3
10     1     4     0.0     1.52     5.11     0.05
10     2     1     0.0     38.83     2.33     0.01
10     2     2     0.0     14.14     2.11     2.29
10     2     3     0.0     2.23     15.6     0.97
10     2     4     0.0     2.52     3.28     0.77
10     3     1     0.0     487.94     2.44     0.16
10     3     2     0.0     35.81     2.22     1.8
10     3     3     0.0     13.28     14.76     2.94
10     3     4     0.0     133.81     3.1     1.16
10     4     1     19.28     3513.38     2.39     0.47
10     4     2     0.0     238.98     2.28     13.46
10     4     3     0.0     20.25     21.99     5.45
10     4     4     0.0     1142.17     3.57     5.52
15     1     1     0.0     18.58     5.98     0.15
15     1     2     0.0     3.38     2.98     0.35
15     1     3     0.0     314.57     9.09     0.2
15     1     4     0.0     10.94     8.58     2.26
15     2     1     6.69     3135.29     3.7     0.67
15     2     2     0.0     2460.78     3.8     4.81
15     2     3     0.0     14.58     87.78     3.68
15     2     4     0.0     1052.16     8.27     2.28
15     3     1     46.33     5760.56     4.12     4.04
15     3     2     20.79     5760.84     4.51     3.57
15     3     3     0.0     322.77     420.37     4.26
15     3     4     14.07     5865.82     7.74     2.37
15     4     1     100.0     7200.47     5.23     3.13
15     4     2     20.2     4421.25     4.35     9.72
15     4     3     0.19     2195.2     237.9     6.28
15     4     4     21.6     7200.42     6.55     11.68
20     1     1     0.0     743.32     13.95     2.78
20     1     2     0.0     110.91     62.75     9.0
20     1     3     1.6     2896.62     17.67     0.13
20     1     4     1.16     3112.33     20.05     1.19
20     2     1     37.26     7200.45     90.1     9.42
20     2     2     19.43     7200.68     5.23     4.37
20     2     3     0.0     1051.22     254.06     5.93
20     2     4     17.15     7200.48     19.33     2.06
20     3     1     78.16     7200.35     6.05     4.07
20     3     2     34.44     5763.72     5.63     6.17
20     3     3     0.73     4530.27     299.19     5.04
20     3     4     30.71     7200.52     22.32     7.34
20     4     1     100.0     7200.63     8.71     6.0
20     4     2     41.32     7200.67     35.68     25.64
20     4     3     1.83     7200.52     307.92     7.45
20     4     4     29.84     7200.52     33.81     9.7

MTZ results with initial solution

To have a clearer view of the results we also present some comparative boxplots obtained from the tables above. First of all, we report the final gap after two hours of running time. We have gathered all the information in Figure <ref>. It is organized in four rows corresponding with the different modes: row $i$ shows results for Mode $i$, $i=1,\ldots,4$. Within each row, there are four columns one per radius size. Then, each graph within this $4 \times 4$ grid contains comparative diagrams for the four different problem sizes considered, namely $n=5,10,15,20$ neighbors. Finally, for each problem size we compare the results obtained for our three different formulations Miller-Tucker-Zemlin (MTZ), Subtour elimination (SEC) and the symmetric version of SEC (sSEC). For instance, looking at the second row, third column (Mode 2, Radius 3) one can see that for $n=5$ and $10$, which correspond to the first two boxes the gap of the three formulations is zero in all the instances (actually the boxes are collapsed to lines). However, for $n=15$ and $20$ MTZ seems to outperforms SEC and sSEC, and moreover, sSEC is also better than SEC since the green boxes lie below the orange ones. As a general comment, one can observe that for all combinations of factors MTZ (the blue boxes) outperforms SEC (orange) and sSEC (green) and also sSEC reports smaller gaps than SEC, with the only exception of Mode 3 where SEC seems to work better than sSEC.

Finally, we have included in Figure <ref> a performance profile graph of number of instances solved versus time. This figure shows that SEC formulation is the weakest since it solves less number of instances in the same cpu time. The comparison between MTZ and sSEC is not that clear although in the long run MTZ seems to outperforms sSEC since the former solves more instances than the latter.

Final gap after 7200 seconds

Performance profile: Time vs #Solved

We also compare next, the behavior of SEC and sSEC in number of cuts required by these two formulations to solve the corresponding problems. As before, we have organized the information in a $4 \times 4$ grid of boxplox graphs. The reader can easily observe that sSEC always requires less number of cuts (blue boxes  corresponding to SEC are always above orange ones corresponding to sSEC) showing that this formulation is more accurate than SEC: reports smaller gaps (see Figure <ref>) and needs less number of cuts.

Number of SEC added in the execution time

§ CONCLUDING REMARKS

This paper has analyzed a novel version of the crossing postman problem with neighbors. We have shown that the problem can be cast within the framework of the family of mixed integer second order cone programming and several exact formulations are presented and computationally tested on an extensive testbed of instances. Additionally, we have presented a heuristic algorithm providing good quality solutions and that can be considered for large scale problems and also as a procedure to obtain initial solutions to be loaded into exact solvers with the exact formulations. Computational results show that the problem is very hard and already for problems with 20 neighbors exact approaches fail to find optimal solution within two hours of cpu time.

This research opens up  several research lines and extensions of the basic problem that can be included in the model. Among them we mention finding better formulations or decomposition schemes that help in solving exactly larger instance sizes; and alternative heuristic algorithms that allow tackling large scale problems. Other extensions of the proposed models considered in this paper are the consideration of unions of second order cone representable sets that can fit more general neighborhoods,  barriers that represent some buildings that the tour can not cross or including conditions that  control the displacement on the border of nonlinear neighborhoods like circles. Some of these topics will be the subject of a follow up paper.

§ ACKNOWLEDGEMENTS

This research has been partially supported by Spanish Ministry of Education and Science/FEDER grant number  MTM2016-74983-C02-(01-02), and projects FEDER-US-1256951, CEI-3-FQM331 and  NetmeetData: Ayudas Fundación BBVA a equipos de investigación científica 2019.

[1]
Amorosi, L., Chiaraviglio, L., D'Andreagiovanni, F., and
Blefari-Melazzi, N.
Energy-efficient mission planning of uavs for 5g coverage in rural
In 2018 IEEE International Conference on Environmental
Engineering (EE) (2018), pp. 1–9.

[2]
Benders, J. F.
Partitioning procedures for solving mixed-variables programming
Numerische Mathematik 4 (1962), 238–252.

[3]
Blanco, V., Fernández, E., and Puerto, J.
Minimum spanning trees with neighborhoods: Mathematical programming
formulations and solution methods.
European Journal of Operational Research 262, 3 (nov 2017),

[4]
Campbell, J. F., Corberán, Á., Plana, I., and Sanchis, J. M.
Drone arc routing problems.
Networks 72, 4 (oct 2018), 543–559.

[5]
Corberán, Á., and Laporte, G., Eds.
Arc Routing: Problems, Methods, and Applications.
Society for Industrial and Applied Mathematics, Philadelphia (USA),
feb 2015.

[6]
Edmonds, J.
Combinatorial Optimization — Eureka, You Shrink!
Springer, 2003, ch. Submodular Functions, Matroids, and Certain
Polyhedra, pp. 11–26.

[7]
Garfinkel, R., and Webb, J.
On crossings, the crossing postman problem, and the rural postman
Networks: An International Journal 34, 3 (1999), 173–180.

[8]
Gentilini, I., Margot, F., and Shimada, K.
The travelling salesman problem with neighbourhoods: Minlp solution.
Optimization Methods and Software 28, 2 (2013), 364–378.

[9]
Geoffrion, A. M.
Generalized Benders decomposition.
Journal of Optimization Theory and Applications 10, 4 (1972),

[10]
Gurobi Optimization, L.
Gurobi Optimizer Reference Manual, Version 8.1.0, 2019.

[11]
Knight, R.
Drones deliver healthcare, 2016.

[12]
Lavars, N.
Amazon to begin testing new delivery drones in the US.
New Atlas 13 (2015).

[13]
Lobo, M. S., Vandenberghe, L., Boyd, S., and Lebret, H.
Applications of second-order cone programming.
Linear Algebra and its Applications 284, 1 (1998), 193 – 228.
International Linear Algebra Society (ILAS) Symposium on Fast
Algorithms for Control, Signals and Image Processing.

[14]
McCormick, G. P.
Computability of global solutions to factorable nonconvex programs:
Part i — convex underestimating problems.
Mathematical Programming 10 (1976), 147–175.

[15]
Miller, C. E., Tucker, A. W., and Zemlin, R. A.
Integer programming formulation of traveling salesman problems.
J. ACM 7, 4 (Oct. 1960), 326–329.

[16]
Mladenović, N., and Hansen, P.
Variable neighborhood search.
Computers & Operations Research 24, 11 (nov 1997),

[17]
Orloff, C. S.
A fundamental problem in vehicle routing.
Networks 4, 1 (1974), 35–64.

[18]
Otto, A., Agatz, N., Campbell, J., Golden, B., and Pesch, E.
Optimization approaches for civil applications of unmanned aerial
vehicles (uavs) or aerial drones: A survey.
Networks 72, 4 (2018), 411–458.

[19]
Pereira, V.
Metaheuristic-Local_Search-Variable_Neighborhood_Search, 2018.

[20]
Sawik, T.
A note on the Miller-Tucker-Zemlin model for the asymmetric
traveling salesman problem.
Bulletin of the Polish Academy of Sciences Technical Sciences
64, 3 (2016), 517–520.

[21]
Velednitsky, M.
Short combinatorial proof that the DFJ polytope is contained in the
MTZ polytope for the asymmetric traveling salesman problem.
Operations Research Letters 45, 4 (jul 2017), 323–324.

[22]
Yuan, B., and Zhang, T.
Towards Solving TSPN with Arbitrary Neighborhoods: A Hybrid
In ACALCI (2017), pp. 204–215.

[1]
Amorosi, L., Chiaraviglio, L., D'Andreagiovanni, F., and
Blefari-Melazzi, N.
Energy-efficient mission planning of uavs for 5g coverage in rural
In 2018 IEEE International Conference on Environmental
Engineering (EE) (2018), pp. 1–9.

[2]
Benders, J. F.
Partitioning procedures for solving mixed-variables programming
Numerische Mathematik 4 (1962), 238–252.

[3]
Blanco, V., Fernández, E., and Puerto, J.
Minimum spanning trees with neighborhoods: Mathematical programming
formulations and solution methods.
European Journal of Operational Research 262, 3 (nov 2017),

[4]
Campbell, J. F., Corberán, Á., Plana, I., and Sanchis, J. M.
Drone arc routing problems.
Networks 72, 4 (oct 2018), 543–559.

[5]
Corberán, Á., and Laporte, G., Eds.
Arc Routing: Problems, Methods, and Applications.
Society for Industrial and Applied Mathematics, Philadelphia (USA),
feb 2015.

[6]
Edmonds, J.
Combinatorial Optimization — Eureka, You Shrink!
Springer, 2003, ch. Submodular Functions, Matroids, and Certain
Polyhedra, pp. 11–26.

[7]
Garfinkel, R., and Webb, J.
On crossings, the crossing postman problem, and the rural postman
Networks: An International Journal 34, 3 (1999), 173–180.

[8]
Gentilini, I., Margot, F., and Shimada, K.
The travelling salesman problem with neighbourhoods: Minlp solution.
Optimization Methods and Software 28, 2 (2013), 364–378.

[9]
Geoffrion, A. M.
Generalized Benders decomposition.
Journal of Optimization Theory and Applications 10, 4 (1972),

[10]
Gurobi Optimization, L.
Gurobi Optimizer Reference Manual, Version 8.1.0, 2019.

[11]
Knight, R.
Drones deliver healthcare, 2016.

[12]
Lavars, N.
Amazon to begin testing new delivery drones in the US.
New Atlas 13 (2015).

[13]
Lobo, M. S., Vandenberghe, L., Boyd, S., and Lebret, H.
Applications of second-order cone programming.
Linear Algebra and its Applications 284, 1 (1998), 193 – 228.
International Linear Algebra Society (ILAS) Symposium on Fast
Algorithms for Control, Signals and Image Processing.

[14]
McCormick, G. P.
Computability of global solutions to factorable nonconvex programs:
Part i — convex underestimating problems.
Mathematical Programming 10 (1976), 147–175.

[15]
Miller, C. E., Tucker, A. W., and Zemlin, R. A.
Integer programming formulation of traveling salesman problems.
J. ACM 7, 4 (Oct. 1960), 326–329.

[16]
Mladenović, N., and Hansen, P.
Variable neighborhood search.
Computers & Operations Research 24, 11 (nov 1997),

[17]
Orloff, C. S.
A fundamental problem in vehicle routing.
Networks 4, 1 (1974), 35–64.

[18]
Otto, A., Agatz, N., Campbell, J., Golden, B., and Pesch, E.
Optimization approaches for civil applications of unmanned aerial
vehicles (uavs) or aerial drones: A survey.
Networks 72, 4 (2018), 411–458.

[19]
Pereira, V.
Metaheuristic-Local_Search-Variable_Neighborhood_Search, 2018.

[20]
Sawik, T.
A note on the Miller-Tucker-Zemlin model for the asymmetric
traveling salesman problem.
Bulletin of the Polish Academy of Sciences Technical Sciences
64, 3 (2016), 517–520.

[21]
Velednitsky, M.
Short combinatorial proof that the DFJ polytope is contained in the
MTZ polytope for the asymmetric traveling salesman problem.
Operations Research Letters 45, 4 (jul 2017), 323–324.

[22]
Yuan, B., and Zhang, T.
Towards Solving TSPN with Arbitrary Neighborhoods: A Hybrid
In ACALCI (2017), pp. 204–215.


# Robust and Efficient Single-Pixel Image Classification with Nonlinear Optics

Santosh Kumar  Department of Physics, Stevens Institute of Technology,
Hoboken, NJ, 07030, USA  Center for Quantum Science and Engineering, Stevens
Institute of Technology, Hoboken, NJ, 07030, USA  These authors contributed
equally  Ting Bu  Department of Physics, Stevens Institute of Technology,
Hoboken, NJ, 07030, USA  Center for Quantum Science and Engineering, Stevens
Institute of Technology, Hoboken, NJ, 07030, USA  These authors contributed
equally  He Zhang  Department of Physics, Stevens Institute of Technology,
Hoboken, NJ, 07030, USA  Center for Quantum Science and Engineering, Stevens
Institute of Technology, Hoboken, NJ, 07030, USA  Irwin Huang  Department of
Physics, Stevens Institute of Technology, Hoboken, NJ, 07030, USA  Yuping
Huang  Department of Physics, Stevens Institute of Technology, Hoboken, NJ,
07030, USA  Center for Quantum Science and Engineering, Stevens Institute of
Technology, Hoboken, NJ, 07030, USA  Corresponding author<EMAIL_ADDRESS>
###### Abstract

We present a hybrid image classifier by mode-selective image upconversion,
single pixel photodetection, and deep learning, aiming at fast processing a
large number of pixels. It utilizes partial Fourier transform to extract the
signature features of images in both the original and Fourier domains, thereby
significantly increasing the classification accuracy and robustness. Tested on
the MNIST handwritten digit images, it boosts the accuracy from 81.25% to
99.23%, and achieves an 83% accuracy for highly contaminated images whose
signal-to-noise ratio is only -17 dB. Our approach could prove useful for fast
lidar data processing, high resolution image recognition, occluded target
identification, atmosphere monitoring, and so on.

††journal: ol

Machine learning techniques based on deep neural networks (DNNs) have scored
considerable success in image classification [1, 2], speech recognition [3],
image generation [4], image reconstruction [5], and so on. They use a feed-
forward multi-layer neural architecture to achieve high performance for
complex operations [6, 7]. Despite notable success for their digital
implementations, it remains a challenge to adopt them directly for embedded
systems due to the high memory-demand, limited scalability, and excessive
energy budget [8]. On the other hand, optical systems naturally offer extreme
parallelism and multiplexing yet with little energy consumption and fast data-
processing [9, 10, 11]. Thus far, optical Fourier transformation, diffraction,
interference, and filtering have been exploited for optical pattern
recognition [12, 13] and restoration [14], phase retrieval [15, 16], and
optical information processing [17]. Recently, optically assisted image
processing has become a vivid pursuit in the field of artificial intelligence
[18], computer vision [19], robotics [20, 21], medical science [22], and so
on.

While most of these studies are based on linear optics, nonlinear optics is
poised to enable even richer and more complex operations and lift our
processing capability to yet another level. To this end, optical neural
network with nonlinear activation functions was implemented for object
identification and classification of ordered and disordered phases of Ising
models [23]. Spontaneous parametric down-conversion was utilized for quantum-
correlated pattern recognition with spatially structured photons [24].
Artificial neural networks based on nonlinear optics have shown advantages in
reconstructing the amplitude and phase profiles of ultrashort pulses [25, 26]
and neuromorphic computing [27, 28]. A super-Ising emulator was demonstrated
with unprecedented four-body interaction using spatial light modulation and
second-harmonic generation [29].

Recently, we have proposed and experimentally demonstrated a nonlinear-optics
approach to pattern recognition with single-pixel imaging and a deep neural
network [30]. It employs mode-selective image up-conversion to project a raw
image with up to mega pixels onto a set of coherent spatial modes, whereby its
signature features are extracted optically in a nonlinear manner. Our
experimental results promise distinct applications in online classification of
large-size images, fast lidar data analyses [31], complex pattern recognition,
and so on.

In this work, we continue to study the aforementioned hybrid machine learning
technique, focusing on further improving the efficiency and effectiveness in
the feature extraction. In our previous system [30], the mode-selective
conversion is applied to the fully Fourier-transformed images, which samples
the signature features in the Fourier domain, but overlooks some important
identifying features in the original image domain. Here, we study a new
approach where the upconversion is applied to the images undergoing only
partial Fourier transform. While challenging to calculate partial Fourier
transform numerically [32, 33], it is easily realized in the current system by
displacing the nonlinear crystal from the focal point of the Fourier lens for
the images. This simple adjustment allows both the original and Fourier
features of the images to be sampled at the same time, leading to high
classification accuracy using only a limited number of pump modes for the
upconversion. Tested on the MNIST (Modified National Institute of Standards
and Technology) handwritten digit images, a 40 mm crystal displacement
significantly increases the accuracy from 81.25% to 99.23%. The same technique
enables robust pattern recognition when the images are mixed with strong
noise, where a high accuracy $\sim 83\%$ is achieved even when the signal-to-
noise (SNR) is as low as -17 dB. All of our experimental results agree well
with the simulated results, confirming this effect. Finally, we numerically
compare two choices of the pump mode bases for the upconversion that can
affect the feature extraction, showing better performance with the Laguerre
Gaussian (LG) mode bases than the Hermite-Gaussian (HG) modes, likely due to a
lower cross-talk for the former [34]. Our results highlight a viable approach
to image recognition through feature extraction in both the original and
Fourier domains by mode-selective image upconversion, for high accuracy, high
efficiency, and exceptional robustness.

Figure 1: Experimental setup for the feature extraction via frequency up-
conversion. The PPLN crystal is moved away from the focal point to the focal
lens by a distance $\Delta$, as shown in the lower right inset. SLM: Spatial
Light Modulator, BS: Beamsplitter, WDM: Wavelength-division multiplexing,
PPLN: Magnesium-doped Periodic Poled Lithium Niobate crystal; PM: power meter.

The experimental setup for the present nonlinear optical pattern-recognition
scheme is outlined in Fig. 1. The phase patterns of the Gaussian beams in
signal and pump arms are modulated by two separate spatial light modulators
(SLM). In the signal arm, SLM1 is used to upload phase patterns of input
images while the SLM2 is applied to map phase patterns of spatial LG (or HG)
modes on pump beam. The phase value is wrapped in the interval between 0 and
$2\pi$ to express on the SLMs [35, 36]. The phase modulated two beams are then
merged at a beam splitter (BS). A Fourier lens ($f$=200 mm) is used to focus
the coupled beam in a temperature-stabilized periodic poled lithium niobate
(PPLN) crystal for sum-frequency (SF) generation process. The power readings
of the SF generation are extracted features from images, which are inputs to a
deep neural network for digits recognition. Further details of our setup can
be found in our previous work [30]. The difference between the current setup
and the previous one is the position of the PPLN. As illustrated in the lower
right inset in Fig. 1, $\Delta$ is the of the displacement of PPLN crystal
towards the Fourier lens.

In simulation, for each pair of the signal beam $E_{s}$ and pump beam $E_{p}$,
their corresponding SF field $E_{f}$ is solved under the slowly-varying-
envelope approximation as:

$\displaystyle
2ik_{s}\partial_{z}E_{s}+(\partial_{x}^{2}+\partial_{y}^{2})E_{s}=-2\frac{\omega_{s}^{2}}{c^{2}}\chi^{(2)}E_{p}^{*}E_{f}e^{i\triangle
kz},$ (1) $\displaystyle
2ik_{p}\partial_{z}E_{p}+(\partial_{x}^{2}+\partial_{y}^{2})E_{p}=-2\frac{\omega_{p}^{2}}{c^{2}}\chi^{(2)}E_{s}^{*}E_{f}e^{i\triangle
kz},$ (2) $\displaystyle
2ik_{f}\partial_{z}E_{f}+(\partial_{x}^{2}+\partial_{y}^{2})E_{f}=-2\frac{\omega_{f}^{2}}{c^{2}}\chi^{(2)}E_{p}E_{s}e^{-i\triangle
kz},$ (3)

where $\omega_{s}$, $\omega_{p}$, and $\omega_{f}$ are the beam waists of
signal, pump and SF light, respectively. $k_{s}={n_{s}\omega_{s}}/{c}$ ,
$k_{p}={n_{p}\omega_{p}}/{c}$ and $k_{f}={n_{f}\omega_{f}}/{c}$ are their wave
numbers. $\Delta k=k_{s}+k_{p}-k_{f}-2\pi/\Lambda$ gives their phase mismatch,
$\Lambda=19.36$ $\mu$m is the poling period of the PPLN. $z$ is the distance
in propagation direction with $z=0$ being the focal point. $\chi^{(2)}$ is the
second-order nonlinear susceptibility. We can use the standard split-step
Fourier and adaptive step size methods [37] to numerically solve Eqs.(1)-(3)
with an initial location $z=-\frac{L}{2}-\Delta$, where $L=1$ $cm$ is the
length of the PPLN.

Figure 2: Normalized confusion matrices for testing the performance of the
hand-written digits with three different crystal positions. The top panels
(a-c) are the simulated results using 40 LG pump modes with $\Delta=0$, 27 mm
and 40 mm, respectively. The bottom panels (d-f) are the corresponding
experimental results.

In this paper, one of the most popular MNIST handwritten digit database is
selected as our benchmark, which is a 10-class database including handwritten
images representing digits from "0" to "9" [38]. A subset of the original
database consisting of the first 200 handwritten images of each digit is
applied as our database. These images in our database are shuffled and
separated into training (1600 images) and testing (400 images) sets. The
resolution of these 28$\times$28 images are increased to 400$\times$400 pixels
to match with the input signal beam size and SLM pitch. The SF power readings
of each handwritten image are normalized between 0 and 1, after collecting all
results for our training and testing data sets. This normalization is
different from our previous work in [30], where the upconversion takes place
in the Fourier plane and data are normalized within each LG mode. While the
original and Fourier planes offer equivalent descriptions [39], the amount of
information that can be extracted using a small number of pump modes is quite
different. In [30], the upconversion was performed in the Fourier domain,
which resulted in similar SF power distributions for the first few pump modes.
The signature features are contained in the higher order modes where the SF
power is relatively low. To distinguish those digits therefore requires
“amplifying” the SF power difference in the higher order modes, which was
realized through SF power normalization by mode. While high performance is
permissible, that method renders the system prone to noise and fluctuations.
Also, after adding new training or testing data, the whole data set needs to
be normalized again and new training process is warranted for the best
performance.

Figure 3: (a) An example of random noise mask with $\textrm{SNR}=-10$ dB. The
zoomed windows are representing the noise clusters with cluster size S equal
to the length of 10 SLM pixels (i.e. one noise pixel size, $\textrm{S}=10$ SLM
pixels or 104 $\mu$m). (b) Experimental and simulated accuracy vs SNR. Blue
and green dotted curves are the experimental and simulated accuracy for
$\Delta$ = 27 mm, respectively, and red and orange dashed curves are for
$\Delta$ = 40 mm.

In contrast, when the upconversion is applied to partially Fourier transformed
images, the signature features of the images are extracted in both the
original and Fourier domains. They are readily manifested in SF power for the
first few LG modes of high readings, resulting in distinct power distribution
over the 40 LG modes. This eliminates the need to amplify the SF power for
higher-order modes, and allows normalization for each image. It makes the
system more robust and efficient, while avoiding the need of repeated
normalization within each mode over the entire data set when new data are
added.

The structure of the deep neural network is same as the one in our previous
work in [30], which is a combination of one convolutional layer and five fully
connected layers. In the training process, the adaptive moment estimation
(ADAM) gradient descent algorithm is selected to minimize the loss between
predictions and their ground truth. The activation function softmax is applied
to normalize the output values between 0 and 1 so that they can represent the
probabilities of each class. Then the predicted class of each input image will
be chosen as the one with the largest probability. Finally, the classification
accuracy can be calculated as the fraction of correctly classified images.

To optimally extract the spatial and Fourier features at the same time, we
displace the PPLN crystal closer to the Fourier lens with its center away from
the focal point by $\Delta$, as shown in Fig.1. While this leads to efficient
and robust classification, a large $\Delta$ de-focuses the light and
suppresses the SF generation. In light of this trade-off, we choose
$\Delta=0$, 27 mm and 40 mm in this study. A larger $\Delta$ may be
accommodated by using a higher pump power to compensate for the SF efficiency
loss. Also, here we follow [30] to use a set of 40 LG modes, with $l\in[-2,2]$
and $p\in[0,7]$. Figure 2 lists the normalized confusion matrices using the 40
LG pump modes for features extraction. Figure 2(a-c) are the simulated results
and (d-f) are the experimental results with $\Delta=0,$ 27 mm and 40 mm,
respectively. When the crystal is centered at the Fourier plane with $\Delta$
= 0, the setup is the same with that in [30]. In this case, the simulated
accuracy is only 77.25% and experimental accuracy is 81.25%. This discrepancy
is likely attributed to the measurement errors in the upconverted SF power
readings for high-order pump projection modes. When the crystal is displaced
by $\Delta=27$ mm, both the simulated and experimental accuracy reached above
97% as shown in Fig. 2(b) and (e), respectively. This shows the advantage of
feature extraction in both spatial and Fourier domains. Figure 2(c) and (f)
show that the accuracy can be further increased over 99% with $\Delta=40$ mm.

To test the robustness of our technique against noise, we next mix the images
with uniformly distributed random noise. The noise level is quantified by the
signal-to-noise ratio (SNR), defined as:

$\textrm{SNR}=10\leavevmode\nobreak\
\textrm{log}_{10}(\sigma_{s}^{2}/\sigma_{n}^{2}),$ (4)

where $\sigma_{s}^{2}$ and $\sigma_{n}^{2}$ are the variance of per-pixel
phase values in the digit images and the added noise masks, respectively.
Figure 3(a) shows an example of one uniformly distributed noise mask with SNR
= -10 dB and cluster size $\textrm{S}=104$ $\mu$m, where the cluster size is
the number of pixels in one noise cluster (which is 10 SLM pixels in this
case). The white noise masks are added to the original handwritten digit
images and are uploaded on SLM1 as the inputs. The results in Fig. 3(b) show
that, with a larger SNR, the accuracy is higher and the difference is smaller
between the simulated and experimental results. Similar to those in Fig. 2 for
clear images, for the noise-contaminated images, both results reveal better
performance with $\Delta$ = 40 mm than $\Delta$ = 27 mm. The advantage reaches
12.18% when the SNR is -15 dB. This again validates our approach for efficient
and robust information extraction.

Figure 4: The accuracy for different sets of LG modes with different azimuthal
($l$) and radial ($p$) indexes, with $\textrm{SNR}=-15$ dB and $\textrm{S}=52$
$\mu$m.

In addition to SNR, the cluster size S is also an import factor to the
classification accuracy. When the noise feature size is comparable with the
image’s characteristic features, it is difficult to distinguish them in both
original and Fourier domains. The simulated and experimental results with
SNR=-15 dB and $\Delta$ = 40 mm in Fig. 3 (b) are 59% and 72.68%,
respectively. However, when we reduce the cluster size by half
($\textrm{S}=52$ $\mu$m), the accuracy increases to 91.25% and 93.4% in
simulation and experiment, respectively.

Thus far, 40 LG modes {LG${}^{p}_{l}$} are chosen for the pump, with
$l\in[-2,2]$ and $p\in[0,7]$, to achieve the accuracy for clear images while
aiming at high efficiency. In practice, one may further optimize the mode
bases for a particular application. To illustrate this opportunity, Fig. 4
plots the simulated accuracy using different sets of LG mode basis with
$\textrm{SNR}=-15$ dB and $\textrm{S}=52$ $\mu$m. It shows that with more
modes, the accuracy can be significantly higher when more than 40 LG modes are
used. With 176 LG modes ($l\in[-5,5]$ and $p\in[0,15]$), the accuracy improved
to 96% from 91%. On the other hand, when the number of LG modes is smaller
than the 40 modes, the accuracy could drop quickly from 87% with $l\in[-2,2]$
and $p\in[0,4]$ to 48% with $l=0$ and $p\in[0,4]$).

Figure 5: Simulated normalized confusion matrices using 49 HG pump modes with
(a) $\Delta$ = 27 mm and (b) $\Delta$ = 40 mm.

Lastly, we compare the choices of LG and HG pump modes to drive the
upconversion for image classification. HG modes are another popular class of
orthogonal mode basis which, unlike LG modes of rotational symmetric
intensity, exhibit rectangular symmetry in intensity. Figure 5 shows the
simulated results with the 49 lowest order HG modes, {HGmn} with $m$ and
$n\in[0,7]$. The classification accuracy of the same data base is 88.5% with
$\Delta=27$ mm and 96.75% with $\Delta=40$ $mm$, which are all inferior to the
simulated results with 40 LG modes in Fig. 2(b) and (c). Therefore, LG modes
are more efficient and suitable for the feature extractions of the handwritten
digits than HG modes. Of course, in practice the optimum mode choice shall
depend on the spatial characteristics of the target images.

In summary, we have demonstrated a new approach to efficient and robust image
classification using a hybrid system integrating nonlinear optics and a deep
neural network. It combines mode selective frequency up-conversion and partial
Fourier transformation to sample the signature features in the original and
Fourier-transformed images through a single nonlinear optical stage. The
partial Fourier transform is conveniently implemented by moving the up-
conversion crystal away from the focal point. It significantly improves the
classification accuracy for both clear and noise-contaminated images. Our
results tap on an inherent advantage of coherent optics for complex machine
learning tasks: the ease and flexibility of linear algebra operations for high
volume data. The robustness of this approach against strong noise invites its
future applications in challenging tasks such as identifying occluded targets,
wide-field surveillance, and remote sensing under low visibility.

Disclosures. The authors declare no conflicts of interest.

## References

  * [1] O. A. Rosso, R. Ospina, and A. C. Frery, PloS one 11, e0166868 (2016).
  * [2] A. T. Sahlol, D. Yousri, A. A. Ewees, M. A. A. Al-qaness, R. Damasevicius, and M. A. Elaziz, Scientific Reports 10, 15364 (2020).
  * [3] L. Deng, G. Hinton, and B. Kingsbury, “New types of deep neural network learning for speech recognition and related applications: An overview,” in _2013 IEEE international conference on acoustics, speech and signal processing,_ (IEEE, 2013), pp. 8599–8603.
  * [4] A. Krizhevsky, I. Sutskever, and G. E. Hinton, Communications of the ACM 60, 84 (2017).
  * [5] G. Wang, J. C. Ye, and B. De Man, Nature Machine Intelligence 2, 737 (2020).
  * [6] J. Schmidhuber, Neural Networks 61, 85 (2015).
  * [7] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ (MIT Press, 2016). http://www.deeplearningbook.org.
  * [8] A. A. Moreno, J. Olivito, J. Resano, and H. Mecha, IEEE Transactions on Very Large Scale Integration (VLSI) Systems 28, 1993 (2020).
  * [9] I. Chakraborty, G. Saha, and K. Roy, Phys. Rev. Applied 11, 014063 (2019).
  * [10] R. Hamerly, L. Bernstein, A. Sludds, M. Soljačić, and D. Englund, Phys. Rev. X 9, 021032 (2019).
  * [11] J. Feldmann, N. Youngblood, M. Karpov, H. Gehring, X. Li, M. Stappers, M. Le Gallo, X. Fu, A. Lukashchuk, A. S. Raja, J. Liu, C. D. Wright, A. Sebastian, T. J. Kippenberg, W. H. P. Pernice, and H. Bhaskaran, Nature 589, 52 (2021).
  * [12] S. Jiao, J. Feng, Y. Gao, T. Lei, Z. Xie, and X. Yuan, Opt. Lett. 44, 5186 (2019).
  * [13] Z. Zhang, X. Li, S. Zheng, M. Yao, G. Zheng, and J. Zhong, Opt. Express 28, 13269 (2020).
  * [14] D. Brunner, A. Marandi, W. Bogaerts, and A. Ozcan, Nanophotonics 9, 4053 (2020).
  * [15] A. Goy, K. Arthur, S. Li, and G. Barbastathis, Phys. Rev. Lett. 121, 243902 (2018).
  * [16] J. Qian, S. Feng, Y. Li, T. Tao, J. Han, Q. Chen, and C. Zuo, Opt. Lett. 45, 1842 (2020).
  * [17] G. R. Steinbrecher, J. P. Olson, D. Englund, and J. Carolan, npj Quantum Information 5, 1 (2019).
  * [18] G. Wetzstein, A. Ozcan, S. Gigan, S. Fan, D. Englund, M. Soljačić, C. Denz, D. A. B. Miller, and D. Psaltis, Nature 588, 39 (2020).
  * [19] P. Latorre-Carmona, V. J. Traver, J. S. Sánchez, and E. Tajahuerce, Image and Vision Computing 86, 28 (2019).
  * [20] Y. J. jae Shin, W. Shin, T. Taniguchi, K. Watanabe, P. Kim, and S.-H. Bae, 2D Materials (2020).
  * [21] P. A. Xu, A. K. Mishra, H. Bai, C. A. Aubin, L. Zullo, and R. F. Shepherd, Science Robotics 4, eaaw6304 (2019).
  * [22] D. P. Sullivan, C. F. Winsnes, L. Åkesson, M. Hjelmare, M. Wiking, R. Schutten, L. Campbell, H. Leifsson, S. Rhodes, A. Nordgren, K. Smith, B. Revaz, B. Finnbogason, A. Szantner, and E. Lundberg, Nature Biotechnology 36, 820 (2018).
  * [23] Y. Zuo, B. Li, Y. Zhao, Y. Jiang, Y.-C. Chen, P. Chen, G.-B. Jo, J. Liu, and S. Du, Optica 6, 1132 (2019).
  * [24] X. Qiu, D. Zhang, W. Zhang, and L. Chen, Phys. Rev. Lett. 122, 123901 (2019).
  * [25] R. Ziv, A. Dikopoltsev, T. Zahavy, I. Rubinstein, P. Sidorenko, O. Cohen, and M. Segev, Optics Express 28, 7528 (2020).
  * [26] T. Zahavy, A. Dikopoltsev, D. Moss, G. I. Haham, O. Cohen, S. Mannor, and M. Segev, Optica 5, 666 (2018).
  * [27] G. Marcucci, D. Pierangeli, and C. Conti, Phys. Rev. Lett. 125, 093901 (2020).
  * [28] D. Ballarini, A. Gianfrate, R. Panico, A. Opala, S. Ghosh, L. Dominici, V. Ardizzone, M. De Giorgi, G. Lerario, G. Gigli _et al._ , Nano Letters 20, 3506 (2020).
  * [29] S. Kumar, H. Zhang, and Y.-P. Huang, Communications Physics 3, 108 (2020).
  * [30] T. Bu, S. Kumar, H. Zhang, I. Huang, and Y.-P. Huang, Optics Letters 45, 6771 (2020).
  * [31] P. Rehain, Y. M. Sua, S. Zhu, I. Dickson, B. Muthuswamy, J. Ramanathan, A. Shahverdi, and Y.-P. Huang, Nature Communications 11, 921 (2020).
  * [32] J. C. Bowman and Z. Ghoggali, Journal of Scientific Computing 76, 1578 (2018).
  * [33] Y. chan Park, J.-G. Jang, and U. Kang, “Fast partial fourier transform,” (2021).
  * [34] S. Restuccia, D. Giovannini, G. Gibson, and M. Padgett, Opt. Express 24, 27127 (2016).
  * [35] S. Kumar, H. Zhang, S. Maruca, and Y.-P. Huang, Opt. Lett. 44, 98 (2019).
  * [36] H. Zhang, S. Kumar, and Y.-P. Huang, Scientific Reports 9, 17481 (2019).
  * [37] G. Agrawal, _Nonlinear Fiber Optics_ (Elsevier Science, 2007).
  * [38] Y. LeCun and C. Cortes, “Mnist handwritten digit database,” (2010).
  * [39] J. W. Goodman, _Introduction to Fourier optics_ (McGraw-hill, 2008).

ref