# On Size and Hardness Generalization in Unsupervised Learning for the
Travelling Salesman Problem

Yimemg Min  Carla P. Gomes

###### Abstract

We study the generalization capability of Unsupervised Learning in solving the
Travelling Salesman Problem (TSP). We use a Graph Neural Network (GNN) trained
with a surrogate loss function to generate an embedding for each node. We use
these embeddings to construct a heat map that indicates the likelihood of each
edge being part of the optimal route. We then apply local search to generate
our final predictions. Our investigation explores how different training
instance sizes, embedding dimensions, and distributions influence the outcomes
of Unsupervised Learning methods. Our results show that training with larger
instance sizes and increasing embedding dimensions can build a more effective
representation, enhancing the model’s ability to solve TSP. Furthermore, in
evaluating generalization across different distributions, we first determine
the hardness of various distributions and explore how different hardnesses
affect the final results. Our findings suggest that models trained on harder
instances exhibit better generalization capabilities, highlighting the
importance of selecting appropriate training instances in solving TSP using
Unsupervised Learning.

Travelling Salesman Problem, Graph Neural Networks, Combinatorial
Optimization, Unsupervised Learning

##  1 Introduction

The goal of machine learning for Combinatorial Optimization (CO) is to enhance
or surpass handcrafted heuristics. Recently, there has been an increasing
trend in applying Machine Learning (ML) to tackle CO problems (Bengio et al.,
2021). Different from manually crafted heuristics, machine learning approaches
harness the power of data to uncover patterns in CO problems.

The Euclidean Travelling Salesman Problem (TSP) is one of the most famous and
intensively studied CO problems. TSP asks the following question: _Given a
list of cities and the distances between each pair of cities, what is the
shortest possible route that visits each city exactly once and returns to the
origin city?_ A variety of methods have been developed to solve TSP, including
the Lin-Kernighan-Helsgaun (LKH) heuristics, which is known for their
effectiveness in approximating solutions (Helsgaun, 2000), and the Concorde
solver, which guarantees optimality of the solutions. The application of ML
for TSP has primarily focused on Supervised Learning (SL) and Reinforcement
Learning (RL). However, SL methods encounter the challenge of expensive
annotations, while RL methods struggle with sparse reward problems.

Recently, (Min et al., 2024) proposes a new approach named UTSP that employs
Unsupervised Learning (UL) to build a data-driven heuristics for the TSP. This
unsupervised method does not depend on any labelled dataset and generates a
heatmap in a non-autoregressive manner, offering a distinct alternative to
traditional SL and RL models.

While the UL heuristics offer a promising approach, the challenge of
generalizing across varying sizes and distributions remains significant. In
particular, the model presented in (Min et al., 2024) requires retraining to
adapt to new sizes, indicating that a model trained on one size cannot
effectively generalize to different sizes.

This paper explores the generalization capabilities of unsupervised heuristics
for the TSP. Our findings indicate that the UL model is able to generalize
across different problem sizes. Regarding the generalization behavior of
different distributions, based on the hardness results by (Gent & Walsh,
1996), we relate different distributions to distinct levels of hardnesses.
This allows us to investigate the impact of the training data’s hardness on
the model’s performance.

Our primary contributions are outlined as follows: We propose a novel approach
for enabling a TSP model, once trained, to generalize effectively across
different problem sizes. We show that training with larger problem sizes can
enhance model performance. Furthermore, we investigate the impact of various
embedding dimensions on TSP performance, finding that larger embedding
dimensions can build more effective representations to guide the search
process. Additionally, we explore how the model performs when trained on
datasets of varying distributions. Our findings indicate that models trained
on harder instances exhibit better performance, which underscores the
importance of training instances’ distribution within the framework of UL for
solving CO problems like the TSP.

While recent research papers explored using data-driven techniques for CO
problems, most have focused on SL or RL. Very few have examined the
generalization behaviours, particularly how training data (different
distributions of TSP instances) influences final model performance (Bi et al.,
2022). Our work addresses this gap, offering insights into the significance of
training data selection and its direct impact on the effectiveness of ML
models for CO tasks. This exploration contributes to understanding ML models
in CO and provides practical guidelines for improving model generalization and
performance in solving TSP.

##  2 Related works

###  2.1 RL for TSP

The goal of using RL for CO is to train an agent capable of either maximizing
or minimizing the expected sum of future rewards, known as the return. For a
given policy, the expected return from a current state is defined as the value
function. In the context of TSP, RL typically focuses on minimizing the length
of the predicted route (Ye et al., 2024; Zhou et al., 2023; Chen et al., 2024;
Ma et al., 2024). For example, (Kool et al., 2019) proposes a model based on
attention layers and trains the model using RL using a deterministic greedy
rollout. (Bello et al., 2016) trains a recurrent neural network to predict
permutations of city coordinates and optimizes the parameters with a policy
gradient method using the negative tour length as a reward signal.

However, as the size of the TSP increases, the rewards become increasingly
sparse, necessitating long exploration steps before the agent achieves a
positive return. So the RL setting is challenging as it only learns once the
agent, randomly or through more sophisticated strategies, finds a better
solution. Additionally, within RL, the learning process is hard to converge,
and the process may become trapped in local minima, as discussed in (Bengio et
al., 2021).

###  2.2 SL For TSP

In SL, the model is trained with a dataset including input coordinates
alongside their corresponding optimal TSP solutions. The objective is to
identify a function that predicts outputs for any given input coordinates,
aiming for these predictions to approximate the optimal solutions (Li et al.,
2024; Sun & Yang, 2024; Fu et al., 2021). For example, (Xin et al., 2021)
trains a Sparse Graph Network using SL to evaluate edge scores, which are then
integrated with the Lin-Kernighan-Helsgaun (LKH) algorithm to guide its search
process. (Fu et al., 2021) uses a GNN to learn from solved optimal solutions.
The model is trained on a small-scale instances, which could be used to build
larger heat maps.

However, In SL, the generation of optimal solutions for training is time-
consuming. Finding optimal or near-optimal solutions for large TSP instances
requires significant computational resources and sophisticated algorithms.

In other words, an ideal model should circumvent these issues, avoiding the
sparse reward problem in RL and not relying on labelled optimal solutuons in
SL. Addressing this, a recent approach by (Min et al., 2024) uses unsupervised
learning (UL) and trains a GNN using a surrogate loss. The model generates
heat maps through a non-autoregressive process, without relying on labelled
optimal solutions or requiring the agents to explore better solutions, thereby
circumventing the need for expensive annotation and mitigating the sparse
reward problem.

This paper is structured as follows: Section 3 introduces the background of UL
for TSP. Section 4 presents a method for generalizing across various problem
sizes. Section 5 investigates the generalization behavior w.r.t. different
embedding dimensions and training sizes. Finally, Section 6 explores the
generalization across different distributions through the lens of instance
hardness.

##  3 UL for TSP

Let’s revisit the definition of the TSP. Essentially, the TSP can be
reinterpreted as identifying the shortest Hamiltonian Cycle that encompasses
all the cities. In UL for TSP, the authors first reformulate the TSP into two
constraints: the shortest path constraint and the Hamiltonian Cycle
constraint. Subsequently, they construct a proxy for each of these constraints
(Min et al., 2024).

In UTSP, given $n$ cities and their coordinates
$(x_{i},y_{i})\in\mathbb{R}^{2}$, UTSP first uses GNN to generate a soft
indicator matrix $\mathbb{T}\in\mathbb{R}^{n\times n}$ and use $\mathbb{T}$ to
build the heat map $\mathcal{H}\in\mathbb{R}^{n\times n}$. Row $i$ of
$\mathcal{H}$ represents the probability distribution of directed edges
originating from city $i$, while column $j$ corresponds to the probability
distribution of directed edges terminating in city $j$. This heat map is
subsequently used to direct a local search algorithm. As mentioned, the goal
of UTSP is to construct a proxy for two constraints. For the shortest
constraint, the authors optimize the distance term:
$\langle\mathbf{D},\mathcal{H}\rangle=\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbf{D}_{i,j}\mathcal{H}_{i,j}$,
where $\langle\cdot,\cdot\rangle$ is the Frobenius inner product,
$\mathbf{D}\in\mathbb{R}^{n\times n}$ is the distance matrix and
$\mathbf{D}_{ij}$ is the distance between city $i$ and city $j$. To address
the Hamiltonian Cycle constraint, the authors introduce the
$\mathbb{T}\rightarrow\mathcal{H}$ transformation, which is designed to
implicitly encode this constraint.

###  3.1 Understanding $\mathbb{T}\rightarrow\mathcal{H}$ transformation

$\mathbb{T}\rightarrow\mathcal{H}$ transformation is defined as:

$\mathcal{H}=\mathbb{T}\mathbb{V}\mathbb{T}^{T},$ (1)

where

$\mathbb{V}=\begin{pmatrix}0&1&0&0&\cdots&0&0&0\\\ 0&0&1&0&\cdots&0&0&0\\\
0&0&0&1&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\ddots&\ddots&\vdots&\vdots&\vdots\\\
0&0&0&0&\ddots&1&0&0\\\ 0&0&0&0&\cdots&0&1&0\\\ 0&0&0&0&\cdots&0&0&1\\\
1&0&0&0&\cdots&0&0&0\end{pmatrix}$

is the Sylvester shift matrix, $\mathbb{V}\in\mathbb{R}^{n\times n}$. We can
interpret $\mathbb{V}$ as a Hamiltonian Cycle starts from city $1\rightarrow
2\rightarrow 3\rightarrow\cdots\rightarrow n\rightarrow 1$ and $\mathbb{T}$ as
an approximation of a general permutation matrix. Given that our starting heat
map is $\mathbb{V}$, which represents a Hamiltonian Cycle, and considering
that the Hamiltonian Cycle is invariant under permutation operations, the
Hamiltonian Cycle constraint is implicitly encoded in this framework, we refer
more details to (Min & Gomes, 2023).

We can also write $\mathbb{T}\rightarrow\mathcal{H}$ transformation as:

$\mathcal{H}=\sum_{t=1}^{n-1}p_{t}p^{T}_{t+1}+p_{n}p_{1}^{T},$ (2)

where $p_{t}\in\mathbb{R}^{n\times 1}$ is the $t_{th}$ column of $\mathbb{T}$,
$\mathbb{T}=[p_{1}|p_{2}|...|p_{n}]$. Equation 2 provides another way of
understanding the $\mathbb{T}\rightarrow\mathcal{H}$ transformation. The
elements in $\mathcal{H}$ are defined using two nearest columns in
$\mathbb{T}$. As shown in Figure 1, $p_{1}=[1,0,0,0,0]^{T}$ and
$p_{2}=[0,0,1,0,0]^{T}$. Since the non-zero element in $p_{1}$ is located at
the first position and the non-zero element in $p_{2}$ is at the third
position, it indicates a directed edge from node 1 to node 3 in the heat map
$\mathcal{H}$. This is depicted as the purple edge in Figure 1. Similarly, the
presence of a non-zero element at the second position in $p_{3}$ implies that
there is a directed edge from node 3 to node 2 in the heat map $\mathcal{H}$,
represented by the yellow edge.

Figure 1: Illustration of $\mathbb{T}$ and the corresponding $\mathcal{H}$.
$p_{1}$[1] = $p_{2}$[3] = $p_{3}$[2] = $p_{4}$[5] = $p_{5}$[4] = 1. This means
there is a corresponding Hamiltonian Cycle: $1\rightarrow 3\rightarrow
2\rightarrow 5\rightarrow 4\rightarrow 1$.

###  3.2 Training UTSP

In UTSP, the author train the model using the following loss $\mathcal{L}$ is:

$\displaystyle\lambda_{1}\underbrace{\sum_{i=1}^{n}(\sum_{j=1}^{n}\mathbb{T}_{i,j}-1)^{2}}_{\text{Row-
wise
constraint}}+\lambda_{2}\underbrace{\sum_{i=1}^{n}\mathcal{H}_{i,i}}_{\text{No
self-
loops}}+\underbrace{\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbf{D}_{i,j}\mathcal{H}_{i,j}}_{\text{Minimize
the distance }}.$ (3)

Here, the _Row-wise constraint_ encourages $\mathbb{T}$ to behave like a
doubly stochastic matrix, thus serving as a soft relaxation of a permutation
matrix (Min & Gomes, 2023). The _No self-loops_ term discourages self loops in
$\mathcal{H}$, where $\lambda_{2}$ is the distance of _self-loop_ , the
_Minimize the Distance_ term acts as a proxy for minimizing the distance of a
Hamiltonian Cycle.

Although UTSP offers a promising unsupervised way to learn the heat maps, a
notable limitation of the model is its lack of generalization. Specifically, a
model trained on TSP instances with $n$ cities cannot be applied to other
instances, such as instances with $n+1$ or $n-1$ cities. This limitation
arises due to $\mathbb{T}$ having a fixed dimension of $\mathbb{R}^{n\times
n}$. Consequently, the model’s architecture is inherently tied to the size of
the training instances, restricting its adaptability to TSP instances of
varying city counts.

##  4 Size Generalization

Recall the understanding of $\mathbb{T}\rightarrow\mathcal{H}$ transformation
in Equation 2. We can interpret that the GNN generates a $n$-dimensional
embedding for each city. In our generalized model, given TSP instances with
different sizes, for each node in these instances, the GNN outputs an
embedding of dimension $m$. Following this, a Softmax activation function is
applied to each column of the embedding matrix, resulting in the generation of
$\mathbb{T}\in\mathbb{R}^{n\times m}$.

We then build $\mathcal{H}$ using111It is important to observe that when
$m\neq n$, $\mathcal{H}$ is not doubly stochastic. We also tried either
replacing $\mathbb{T}$ with $\sqrt{\frac{n}{m}}\mathbb{T}$ or substituting
$\mathcal{H}$ with $\frac{n}{m}\mathcal{H}$, both of which yield similar
outcomes.:

$\mathcal{H}=\sum_{t=1}^{m-1}p_{t}p^{T}_{t+1}+p_{m}p_{1}^{T},$ (4)

where $p_{t}\in\mathbb{R}^{n\times 1}$ is the $t_{th}$ column of $\mathbb{T}$.
Equation 4 can be reformulated analogously to Equation 1 with
$\mathbb{V}\in\mathbb{R}^{m\times m}$.

Figure 2: Training history of different $m$ (embedding dimension) on TSP-2000.
Figure 3: Training history of different $n$ (instance size) with same
embedding dimension $m=1500$.

In practice, we train our model under the loss $\mathcal{L}$:

$\displaystyle\lambda_{1}\underbrace{\bigl{\\{}\sum_{j=1}^{n}(1-\sum_{i=1}^{n}\mathcal{H}_{i,j})^{2}+\sum_{i=1}^{n}(1-\sum_{j=1}^{n}\mathcal{H}_{i,j})^{2}\bigl{\\}}}_{\text{Row
and column-wise constraint}}$ (5)
$\displaystyle+\underbrace{\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbf{D}_{i,j}\mathcal{H}_{i,j}}_{\text{Minimize
the distance }}.$

By letting the GNN to output an $m$-dimensional embedding for each city, the
model achieves generalization across different instances. This means that,
through Equation 2, the heat map $\mathcal{H}$ will consistently match the
size of the input cities ($n\times n$ ).

##  5 Experiment

Table 1: Results of our generalizable model + Local Search w.r.t. existing
baselines, tested on 128 instances with $n$ = 200, 500 and 1000.
Method | Type | TSP200 | TSP500 | TSP1000
---|---|---|---|---
Length | Gap (%) | Time | Length | Gap (%) | Time | Length | Gap (%) | Time
Concorde | Solver | 10.7191 | 0.0000 | 3.44m | 16.5458 | 0.0000 | 37.66m | 23.1182 | 0.0000 | 6.65h
Gurobi | Solver | 10.7036 | -0.1446 | 40.49m | 16.5171 | -0.1733 | 45.63h | - | - | -
LKH3 | Heuristic | 10.7195 | 0.0040 | 2.01m | 16.5463 | 0.0029 | 11.41m | 23.1190 | 0.0036 | 38.09m
GAT (Deudon et al., 2018) | RL, S | 13.1746 | 22.9079 | 4.84m | 28.6291 | 73.0293 | 20.18m | 50.3018 | 117.5860 | 37.07m
GAT (Kool et al., 2019) | RL, BS | 11.3769 | 6.1364 | 5.77m | 19.5283 | 18.0257 | 21.99m | 29.9048 | 29.2359 | 1.64h
GCN (Joshi et al., 2019) | SL, G | 17.0141 | 58.7272 | 59.11s | 29.7173 | 79.6063 | 6.67m | 48.6151 | 110.2900 | 28.52m
Att-GCRN(Fu et al., 2021) |  SL+RL MCTS | 10.7358 | 0.1563 | 20.62s + | 16.7471 | 1.2169 | 31.17s + | 23.5153 | 1.7179 | 43.94s +
1.33m | 3.33m | 6.68m
UTSP (Min et al., 2024) | UL, Search | 10.7289 | 0.0918 | 4.83s + | 16.6846 | 0.8394 | 7.28s + | 23.3903 | 1.1770 | 0.23m+
1.11m | 1.54m | 3.51m
Our Model | UL, Search | 10.7251 | 0.0558 | 4.94s + | 16.6820 | 0.8229 | 5.66s + | 23.3867 | 1.1616 | 0.24m+
1.11m | 1.54m | 3.51m

Here, we explore the impact of the generalized model on different problem
sizes. Specifically, we study TSP with 200, 500, and 1000 cities, each size is
evaluated using 128 test instances.

Different from previously UTSP setting, our new methodology involves training
models on larger datasets and testing them on smaller ones. Specifically, we
train a model on a TSP-2000 dataset with $m=1500$ and test it on a TSP-1000
dataset; another model is trained on TSP-1000 with $m=800$ and tested on
TSP-500; and finally, a model trained on TSP-400 with $m=320$ is tested on
TSP-200. The TSP-2000, 1000, and 400 training datasets are created by randomly
distributing points on a 2D plane, subject to a uniform distribution. For
TSP-200 and TSP-400, we train the model for 300 epochs, while for TSP-1000, we
train the model for 600 epochs. Each of these datasets consists of 5,000
training instances.

We train our model on one NVIDIA A100 Graphics Processing Unit, using the same
Graph Neural Network (GNN) architecture as described in (Min et al., 2024).
The model is trained on TSP instances of sizes 400, 1000, and 2000, using a
configuration of two hidden layers, with each layer comprising 128 hidden
units. The hyperparameter $\lambda_{1}$, as specified in Equation 4, is set to
100. Our test instances are taken from (Fu et al., 2021). Here, the
performance gap is calculated using the $\frac{l-l_{opt}}{l_{opt}}$, where $l$
represents the TSP length generated by our model and $l_{opt}$ denotes the
optimal length. We run the search algorithm on Intel Xeon Gold 6326.

In our approach, consistent with the existing UTSP framework, we employ the
same search methodology. The process begins with the generation of the heat
map $\mathcal{H}$, from which we extract the top $M$ largest values in each
row. This extraction leads to the formation of a new heat map, denoted as
$\tilde{H}$. We compute $\mathcal{H}^{\prime}=\tilde{H}+\tilde{H}^{T}$ to
symmetrize this updated heat map. $\mathcal{H}^{\prime}$ is then used to guide
the search process. We further calculate the overlap between non-zero edges in
$\mathcal{H}^{\prime}$ and the optimal solutions, where a higher overlap ratio
indicates that $\mathcal{H}^{\prime}$ more effectively covers the optimal
solution. For more detailed information, we refer to (Min et al., 2024).

Our results are shown in Table 1, in the case of TSP-200, our model achieves a
gap of 0.0558 %, when tackling TSP-500, the model continues to demonstrate its
robustness, with a gap of 0.8229%. The performance in both TSP-200 and TSP-500
suggests that our model’s approach to guiding the local search is effective
across various scales of the TSP.

When the model is applied to the largest tested instance size, TSP-1000, it
achieves a gap of 1.1616%. which is the minimum one among all the methods.
More importantly, it underscores the model’s generalization to scale and
maintain a level of efficiency in large-scale TSP instances. Our results
across all three instance sizes illustrate that the model trained using
Equation 5 is able to generalize across instances of different sizes and
effectively enhances the search process.

Figure 4: Expansion

Figure 5: Explosion

Figure 6: Implosion

Figure 7: Uniform

###  5.1 Impact of Varying $m$ on Training Performance

As mentioned in Equation 4, $m$ represents the embedding dimension of each
node. In this study, we investigate the effect of the embedding dimension $m$
on the model’s performance. Specifically, we train models on TSP-2000
instances with varying embedding dimensions: $m=500,1000,$ and $1500$. We then
evaluate these models on TSP-1000 instances to assess their performance.

Table 2: Overlap ratios and the search results on 128 TSP-1000 instances in
(Fu et al., 2021) using different embedding dimension $m$. We select top 5
elements from each row in the heat maps.
$m$ | Overlap Ratio(%) | Performance Gap(%)
---|---|---
500 | 82.70 |  2.0746 $\pm$ 0.5457
1000 | 93.75 |  1.4832 $\pm$ 0.2305
1500 | 94.93 |  1.4145 $\pm$ 0.2005

The training curves for different embedding dimensions are shown in Figure 2.
We calculate the overlap ratios and search performance using models with
different embedding dimensions, the results are shown in Table 2, 3. Our
findings indicate that an increase in the embedding dimension contributes to
higher overlap ratios and enhanced search performance. For instance, the
overlap ratio improves from 82.70% to 94.93% when the embedding dimension $m$
is increased from 500 to 1500, based on the heat maps with top 5 elements from
each row. Correspondingly, the search performance also improves, with the gap
decreasing from 2.0746% to 1.4145%. This highlights the significance of
embedding dimension in increasing model efficacy. A larger embedding dimension
can better identify optimal or near-optimal solutions and narrow the gap.

Table 3: Overlap ratios and the search results on 128 TSP-1000 instances in
(Fu et al., 2021) using different embedding dimension $m$. We select top 20
elements from each row in the heat maps.
$m$ | Overlap Ratio(%) | Performance Gap(%)
---|---|---
500 | 99.99 |  1.1995 $\pm$ 0.1849
1000 | 100.00 |  1.1608 $\pm$ 0.1844
1500 | 100.00 |  1.1616 $\pm$ 0.1743

Specifically, it is noteworthy that when selecting the top 20 elements from
each row, both $m=1000$ and $m=1500$ achieve a 100.00% overlap ratio, whereas
$m=500$ does not cover all the optimal solutions, resulting in a larger gap.
Furthermore, we observe that $m=1000$ exhibits marginally better performance
compared to $m=1500$. This suggests that beyond a certain threshold,
increasing the embedding dimension yields diminishing returns in terms of
covering optimal solutions. It also implies that there might be an optimal
range for the embedding dimension, indicating a need for careful consideration
in the choice of $m$ to optimize model performance.

###  5.2 Impact of Varying $n$ on Training Performance

Our model can generalize across different sizes, meaning that training on one
size can effectively translate to performance on another, previously unseen
size. Here we investigate how varying the training size impacts the model’s
performance. We train the model using TSP-400, TSP-1000, and TSP-2000
instances, all with the same embedding dimension $m=1500$. The training
results are illustrated in Figure 3.

We then test how different training instances’ sizes can affect the overlap
ratio and the performance. The results are shown in Table 4, 5. We note that
training with larger instances enhances search performance under both top 5
and top 20 conditions. Specifically, when selecting the top 5 elements from
each row, the performance gap improves from 3.0762% to 1.4145%. Similarly,
when choosing the top 20 elements from each row, the gap shows a marked
improvement, decreasing from 1.1885% to 1.1616%.

Table 4: Overlap ratios and the search results on 128 TSP-1000 instances
instances in (Fu et al., 2021) with $m=1500$ using training instances with
different sizes. We select top 5 elements from each row in the heat maps. The
first column denotes different training sizes.
$n$ | Overlap Ratio(%) | Performance Gap(%)
---|---|---
400 | 68.83 |  3.0762 $\pm$ 1.3141
1000 | 93.48 |  1.5563 $\pm$ 0.2345
2000 | 94.93 |  1.4145 $\pm$ 0.2005
Table 5: Overlap ratios and the search results on 128 TSP-1000 instances in
(Fu et al., 2021) with $m=1500$ using training instances with different sizes.
We select top 20 elements from each row in the heat maps. The first column
denotes different training sizes.
$n$ | Overlap Ratio(%) | Performance Gap(%)
---|---|---
400 | 99.96 |  1.1885 $\pm$ 0.1927
1000 | 100.00 |  1.1763 $\pm$ 0.1743
2000 | 100.00 |  1.1616 $\pm$ 0.1743

Our results highlight the importance of selecting larger training instance
sizes to enhance model performance and efficiency.

##  6 Hardness Generalization

Previous studies suggest that UL can generalize across different sizes, guide
the search and reduce the search space, Here, we delve into how UL’s
capability to reduce the search space is influenced by different
distributions. Specifically, we explore the relationship between different
distributions and the efficiency of using UL for solving the TSP.

However, building a connection between various distributions and the efficacy
of UL in reducing the search space presents significant challenges. To address
this, we first focus on correlating different distributions with their
hardness levels.

#### Phase transition

A phase transition refers to a change in the solvability of NP-hard problems.
When some parameters of the problem is varied, for example, the density of
constraints in a Boolean satisfiability problem (SAT) problem (Mitchell et
al., 1992), the problem undergoes a transition from being almost solvable to
unsolvable. To be specific, The phase transition in SAT refers to a sharp
change in the solvability of these problems, depending on the ratio of the
number of clauses to the number of variables in the formula. When the ratio is
low (few clauses relative to variables), most instances of the problem are
easy to solve. This is because there are fewer constraints, making it more
likely to find a satisfying assignment. Conversely, when this ratio is high
(many clauses relative to variables), the problem becomes over-constrained,
and most instances are also easy to solve because they are almost certainly
unsatisfiable. The most interesting part occurs at a certain critical ratio,
typically around 4.3 for 3-SAT problems. At this ratio, the problems undergo a
phase transition and become extremely hard to solve. In other words, the
problems are most difficult around the phase transition point (Monasson et
al., 1999).

Phase transitions provides a powerful framework to study the properties of NP-
hard problems. However, the exact nature and location of these transitions can
be difficult to determine and may depend intricately on the structure of the
different problems. For TSP, (Gent & Walsh, 1996) suggest using the parameter
$\tau=l_{opt}/\sqrt{nA}$, where $A$ denotes the area covered by the TSP
instance, $l_{opt}$ represents the length of the optimal solution, and $n$ is
the number of cities. This approach is based on the observation that there is
a rapid transition in solvability around a fixed value of the parameter,
specifically at approximately $T_{c}=0.78$.

Figure 8: TSP phase transition and the $\tau$ values for different
distributions.

Here we study four different distributions and see how it can effect the
search space reduction, an illustration of these four distribution is shown in
Figure 7 $\sim$ 7. As mentioned earlier, around the phase transition point,
the problems often exhibits the greatest computational complexity (Hard).
Figure 8 illustrates the scheme of phase transition in the TSP. The $x$-axis
is the $\tau$ value, while the $y$-axis corresponds to the level of hardness.
The point at which $\tau$ equals the critical threshold $T_{c}=0.78$ marks the
peak of difficulty, exhibiting the highest hardness, we refer more details to
(Gent & Walsh, 1996).

Furthermore, we present the $\tau$ values for four different distributions,
where each $\tau$ is computed as an average from 100 instances, each with a
size of 200, 500 and 1000, detailed in Table 6 and Figure 8.

Table 6: $\tau=l_{opt}/\sqrt{nA}$ of Expansion, Explosion, Implosion, and
Uniform for different sizes.
Size | Expansion | Explosion | Implosion | Uniform
---|---|---|---|---
1000 | 0.4838 | 0.5629 | 0.7237 | 0.7460
500 | 0.5114 | 0.5905 | 0.7338 | 0.7515
200 | 0.5796 | 0.6337 | 0.7539 | 0.7745

As shown in Figure 8, the Uniform distribution is closest to the phase
transition point $T_{c}$. This indicates a highest level of hardness.
Consequently, in terms of transitioning from hard to easy, the order is
observed as follows: Uniform $\approx$ Implosion $>$ Explosion $>$ Expansion.
Following upon this concept, we examine how these distributions influence the
capacity of UL to efficiently reduce the search space and guide the search.

Figure 9: The training curves for TSP-400 with $m=320$ across four different
distributions are shown; the model is then tested on 128 TSP-200 instances.
Figure 10: The training curves for TSP-1000 with $m=800$ across four different
distributions are shown; the model is then tested on 128 TSP-500 instances.
Table 7: Overlap ratios and the search results on 128 TSP-200 instances in (Fu
et al., 2021) using different distributions. We select top 5 elements from
each row in the heat maps.
Dataset | Overlap Ratio(%) | Performance Gap(%)
---|---|---
Uniform | 95.64 |  0.0883$\pm$ 0.0885
Implosion | 95.50 |  0.0876$\pm$ 0.0920
Explosion | 95.29 |  0.0979$\pm$ 0.0907
Expansion | 94.00 |  0.1131$\pm$ 0.0973
Table 8: Overlap ratios and the search results on 128 TSP-500 instances in (Fu
et al., 2021) using different distributions. We select top 5 elements from
each row in the heat maps.
Dataset | Overlap Ratio(%) | Performance Gap(%)
---|---|---
Uniform | 95.47 |  0.9311$\pm$ 0.1638
Implosion | 95.40 |  0.9394$\pm$ 0.1732
Explosion | 94.99 |  0.9410$\pm$ 0.1764
Expansion | 94.03 |  1.0137$\pm$ 0.1800
Table 9: Overlap ratios and the search results on 128 TSP-1000 instances in
(Fu et al., 2021) using different distributions. We select top 5 elements from
each row in the heat maps.
Dataset | Overlap Ratio(%) | Performance Gap(%)
---|---|---
Uniform | 94.93 |  1.4145 $\pm$ 0.2005
Implosion | 94.71 |  1.4060 $\pm$ 0.2078
Explosion | 93.86 |  1.5274 $\pm$ 0.2632
Expansion | 93.38 |  1.5777 $\pm$ 0.2735

Figure 11: The training curves on TSP-2000 with $m=1500$ across four different
distributions are shown; the model is then tested on 128 TSP-1000 instances.

We first train the models using 4 different distributions with the same
parameters in Section 4. We calculate the overlap ratio of these models for
TSP-200, 500, and 1000. The training results are shown in Figure 9, 10 and 11.
We observe that models trained with harder instances consistently exhibit a
lower loss. Specifically, the loss curves for models trained using the Uniform
distribution consistently show the lowest loss, while those trained with
Expansion and Explosion distributions demonstrate higher losses. This suggests
that the hardness level of training instances plays a significant role in the
effectiveness of the model training, directly impacting the loss metrics. It
is important to note that throughout our training process, all other
hyperparameter settings remained constant. Therefore, the observed variations
in loss can be attributed solely to the differences in training distributions.

Table 10: Overlap ratios and the search results on 128 TSP-1000 instances in
(Fu et al., 2021) using different distributions. We select top 20 elements
from each row in the heat maps.
Dataset | Overlap Ratio(%) | Performance Gap(%)
---|---|---
Uniform | 100.00 |  1.1616$\pm$ 0.1743
Implosion | 100.00 |  1.1844$\pm$ 0.1572
Explosion | 100.00 |  1.1937$\pm$ 0.1764
Expansion | 100.00 |  1.1797$\pm$ 0.2025

We then evaluate how different distributions can affect the search results. We
pick the top 5 element each row and build the heat maps. The overlap ratio and
the search results are shown in Table 7, 8 and 9. When training on easier
distributions such as Explosion and Expansion, we observe low overlap ratios
and larger performance gaps. This indicates that models trained on simpler
distributions may struggle to generalize effectively to more challenging
instances of the problem. The lower overlap ratios suggest that the solutions
generated by these models are less aligned with the optimal solutions, and the
larger performance gaps highlight a significant disparity in effectiveness
when these models are applied to the test TSP instances. Training on harder
distributions, such as Uniform, yields higher overlap ratios and improved
search performance. This indicates that models trained on harder distributions
can build a better representation of the search space, which enables the
search to perform more effectively. It is also observed that the _plateaus_
during training are more pronounced when training on harder instances,
suggesting that the optimization landscape becomes more complex when the
hardness level increases.

We evaluate the model’s performance on TSP-1000 instances by utilizing the top
20 elements from each row for each distribution, as detailed in Table 10. We
observe that by selecting the top 20 elements, $\mathcal{H}^{\prime}$ is able
to cover 100.00% of the optimal solutions. Overall, the performance gaps
across the distributions are similar, with training on uniform distributions
continuing to exhibit the lowest performance gap.

##  7 Conclusion

This work introduces a new methodology that allows a trained, unsupervised TSP
model to generalize across different problem sizes. Our results demonstrate
that training on larger problem instances can improve performance compared to
training with smaller instances. Additionally, we delve into the influence of
embedding dimensions on TSP results, showing that larger embedding dimensions
are important in constructing more effective representations that guide the
search process more efficiently. Moreover, we investigate the model’s
performance using training datasets with different levels of hardnesses. We
show that training on harder instances can improve model performance,
emphasizing the importance of selecting training instances with appropriate
difficulty levels. We train our models on different TSP distributions to
understand their impact on the effectiveness of UL models. Our study indicates
a clear relationship between the inherent hardness of distribution and the
model’s capacity to generalize and effectively solve TSP instances. To our
knowledge, this is the first study to systematically investigate and
demonstrate this connection.

Our results highlight the relationship between the characteristics of training
instances (size and hardness), embedding dimensions, and model performance in
UL, particularly when addressing CO problems such as the TSP. We anticipate
that these findings — emphasizing the benefits of training on larger, harder
instances with increased embedding dimensions — can inspire further research
in the application of Unsupervised Learning to Combinatorial Optimization
tasks.

##  8 Acknowledgement

This project is partially supported by the Eric and Wendy Schmidt AI in
Science Postdoctoral Fellowship, a Schmidt Futures program; the National
Science Foundation (NSF) and the National Institute of Food and Agriculture
(NIFA); the Air Force Office of Scientific Research (AFOSR); the Department of
Energy; and the Toyota Research Institute (TRI).

## References

  * Bello et al. (2016) Bello, I., Pham, H., Le, Q. V., Norouzi, M., and Bengio, S.  Neural combinatorial optimization with reinforcement learning.  _arXiv preprint arXiv:1611.09940_ , 2016.
  * Bengio et al. (2021) Bengio, Y., Lodi, A., and Prouvost, A.  Machine learning for combinatorial optimization: a methodological tour d’horizon.  _European Journal of Operational Research_ , 290(2):405–421, 2021.
  * Bi et al. (2022) Bi, J., Ma, Y., Wang, J., Cao, Z., Chen, J., Sun, Y., and Chee, Y. M.  Learning generalizable models for vehicle routing problems via knowledge distillation.  _Advances in Neural Information Processing Systems_ , 35:31226–31238, 2022.
  * Chen et al. (2024) Chen, J., Wang, J., Zhang, Z., Cao, Z., Ye, T., and Chen, S.  Efficient meta neural heuristic for multi-objective combinatorial optimization.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Deudon et al. (2018) Deudon, M., Cournut, P., Lacoste, A., Adulyasak, Y., and Rousseau, L.-M.  Learning heuristics for the tsp by policy gradient.  In _International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research_ , pp. 170–181. Springer, 2018.
  * Fu et al. (2021) Fu, Z.-H., Qiu, K.-B., and Zha, H.  Generalize a small pre-trained model to arbitrarily large tsp instances.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(8):7474–7482, 2021.
  * Gent & Walsh (1996) Gent, I. P. and Walsh, T.  The tsp phase transition.  _Artificial Intelligence_ , 88(1-2):349–358, 1996\.
  * Helsgaun (2000) Helsgaun, K.  An effective implementation of the lin–kernighan traveling salesman heuristic.  _European journal of operational research_ , 126(1):106–130, 2000.
  * Joshi et al. (2019) Joshi, C. K., Laurent, T., and Bresson, X.  An efficient graph convolutional network technique for the travelling salesman problem.  _arXiv preprint arXiv:1906.01227_ , 2019.
  * Kool et al. (2019) Kool, W., van Hoof, H., and Welling, M.  Attention, learn to solve routing problems!  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=ByxBFsRqYm.
  * Li et al. (2024) Li, Y., Guo, J., Wang, R., and Yan, J.  From distribution learning in training to gradient search in testing for combinatorial optimization.  volume 36, 2024.
  * Ma et al. (2024) Ma, Y., Cao, Z., and Chee, Y. M.  Learning to search feasible and infeasible regions of routing problems with flexible neural k-opt.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Min & Gomes (2023) Min, Y. and Gomes, C.  Unsupervised learning permutations for tsp using gumbel-sinkhorn operator.  In _NeurIPS 2023 Workshop Optimal Transport and Machine Learning_ , 2023.
  * Min et al. (2024) Min, Y., Bai, Y., and Gomes, C. P.  Unsupervised learning for solving the travelling salesman problem.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Mitchell et al. (1992) Mitchell, D., Selman, B., Levesque, H., et al.  Hard and easy distributions of sat problems.  In _Aaai_ , volume 92, pp. 459–465, 1992.
  * Monasson et al. (1999) Monasson, R., Zecchina, R., Kirkpatrick, S., Selman, B., and Troyansky, L.  Determining computational complexity from characteristic ‘phase transitions’.  _Nature_ , 400(6740):133–137, 1999.
  * Sun & Yang (2024) Sun, Z. and Yang, Y.  Difusco: Graph-based diffusion solvers for combinatorial optimization.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Xin et al. (2021) Xin, L., Song, W., Cao, Z., and Zhang, J.  Neurolkh: Combining deep learning model with lin-kernighan-helsgaun heuristic for solving the traveling salesman problem.  _Advances in Neural Information Processing Systems_ , 34:7472–7483, 2021.
  * Ye et al. (2024) Ye, H., Wang, J., Cao, Z., Liang, H., and Li, Y.  Deepaco: Neural-enhanced ant systems for combinatorial optimization.  volume 36, 2024.
  * Zhou et al. (2023) Zhou, J., Wu, Y., Song, W., Cao, Z., and Zhang, J.  Towards omni-generalizable neural methods for vehicle routing problems.  In _International Conference on Machine Learning_ , pp. 42769–42789. PMLR, 2023.


# UAV-Enabled Mobile Edge Computing: Offloading Optimization and Trajectory
Design

Fuhui Zhou§‡, Yongpeng Wu∗, Haijian Sun§, Zheng Chu†
§Utah State University, USA, ‡Nanchang University, Nanchang, China,
∗Shanghai Jiao Tong University, China, †University of Surrey, Guildford, U.K
Email: _{<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>The research was supported by the U.S. National
Science Foundation grant EARS-1547312, the National Natural Science Foundation
of China (61701214, 61701301, 61661028, 61631015, and 61561034), the Young
Natural Science Foundation of Jiangxi Province (20171BAB212002), the China
Postdoctoral Science Foundation (2017M610400) and the Postdoctoral schedule
fund of Jiangxi Province(2017RC17).

###### Abstract

With the emergence of diverse mobile applications (such as augmented reality),
the quality of experience of mobile users is greatly limited by their
computation capacity and finite battery lifetime. Mobile edge computing (MEC)
and wireless power transfer are promising to address this issue. However,
these two techniques are susceptible to propagation delay and loss. Motivated
by the chance of short-distance line-of-sight achieved by leveraging unmanned
aerial vehicle (UAV) communications, an UAV-enabled wireless powered MEC
system is studied. A power minimization problem is formulated subject to the
constraints on the number of the computation bits and energy harvesting
causality. The problem is non-convex and challenging to tackle. An alternative
optimization algorithm is proposed based on sequential convex optimization.
Simulation results show that our proposed design is superior to other
benchmark schemes and the proposed algorithm is efficient in terms of the
convergence.

###### Index Terms:

Mobile edge computing, resource allocation, unmanned aerial vehicle
communications, trajectory optimization, wireless power transfer.

##  I Introduction

WITH the development of Internet of Things (IoT), the emerging diverse mobile
applications (augmented reality, face recognition, mobile online gaming, etc.)
enable mobile users to enjoy a high quality of experience [1]. However, these
applications are latency-sensitive and need a high computation capability. Due
to the limited battery and low computation capability, it is challenging for
mobile devices to execute these applications [2], [3]. Fortunately, mobile
edge computing (MEC) has been recognized as a promising technique to tackle
this challenge [4]. It provides the edge network with cloud computing service.
Mobile users can offload their computation tasks into the edge network. Unlike
mobile cloud computing (MCC), network edge devices in MEC, such as access
points, can perform cloud-like computing and are deployed in close proximity
to users. MEC has received increasing attention in both academia and industry
since it has advantages of saving energy for users, providing low latency
services and achieving security for mobile applications [5]-[9].

On the other hand, wireless power transfer (WPT) techniques are promising for
prolonging the operational time of energy-limited mobile devices [10], [11].
Particularly, radio frequency (RF) signals are used as energy sources for
energy harvesting (EH). Compared to the conventional EH techniques, such as
solar charging, WPT techniques can provide a controllable and stable power.
They are important for energy-limited mobile devices, which are required to
execute local low-computation tasks and offload computation-extensive tasks
[6]-[9]. However, the harvested power level is greatly influenced by the
severe propagation loss. Recently, an unmanned aerial vehicle (UAV)-enabled
WPT architecture has been proposed to improve the energy transfer efficiency
[12]. It utilizes an UAV as an energy transmitter for powering the ground
mobile users. It was shown that the harvested power level can be significantly
improved due to the higher chance of short-distance line-of-sight (LoS) energy
transmit links [13].

Motivated by the UAV-enabled WPT architecture, a UAV-enabled wireless powered
MEC system is studied in this paper. In the system, the UAV transmits energy
to multiple ground users and the ground users exploit the harvested energy for
local computing and computation tasks offloading. _To the authors’ best
knowledge, this is the first work that establishes an UAV-enabled wireless
powered MEC system and studies the joint optimization of computation
offloading and trajectory design_. The related works are summarized as
follows.

In [5], the revenue of the wireless cellular networks with MEC was maximized
by jointly optimizing the computation offloading decision and resource
allocation. The authors in [6]-[9] extended the resource allocation problems
into wireless powered MEC systems. Specifically, in [6], an energy-efficient
resource allocation strategy was proposed by jointly optimizing the number of
the local computation bits and the offloading computation bits under the
causal energy harvesting constraint. The authors in [7] proposed an efficient
reinforcement learning-based resource management scheme in an MEC system with
energy harvesting. It was shown that the computation performance can be
significantly improved by using the proposed algorithm. In [8] and [9], a
joint optimization framework was proposed in wireless powered MEC systems with
different operational paradigms, namely, binary and partial offloading,
respectively. In the binary offloading paradigm, the computation task is
completely executed in mobile devices or offloaded into network edge devices
for computing. In the partial offloading paradigm, the computation task can be
divided into two parts, one for local computing and the other for offloading.
The computation rate was maximized in [8] and the transmit power was minimized
in [9] by jointly optimizing the computing frequency and the transmit power.

In [6]-[9], the access point or the energy transmitter equipped with an MEC
server is deployed at the fixed location. It results in a low energy transfer
efficiency due to the severe propagation loss [10], [11]. In order to tackle
this issue, the authors designed UAV-enabled wireless powered systems and
jointly optimized the resource allocation and trajectory of the UAV [12],
[13]. However, MEC was not considered in [12] and [13]. Recently, an UAV-
enabled MEC system was designed in [14] and the transmit power of users was
minimized by jointly optimizing the number of the local computation bits, the
offloading computation bits and the downloading bits.

Different from [12]-[14], an UAV-enabled wireless powered MEC system is
studied in this paper. A power minimization problem is formulated by jointly
optimizing the number of the offloading computation bits, the local
computation frequencies of users and the UAV, and the trajectory of the UAV.
It is challenging to solve the formulated non-convex problem due to the
existing couple among the optimized variables. An alternative optimization
algorithm is proposed to solve it by using the sequential convex optimization
(SCA) techniques. Simulation results show that our proposed resource
allocation scheme outperforms other benchmark schemes.

The remainder of this paper is organized as follows. The system model is
presented in Section II. Section III presents the energy minimization problem.
Section IV presents simulation results. The paper concludes with Section V.

##  II System Model

Figure 1: The system model.

An UAV-enabled wireless powered MEC system is considered in Fig. 1, where an
UAV equipped with an MEC server transmits energy to $K$ users and provides MEC
services for these users. In this paper, the partial offloading paradigm is
applied. Similar to [8] and [9], users can simultaneously perform energy
harvesting, local computing and computation offloading. Without loss of
generality, a three-dimensional (3D) Euclidean coordinate is adopted. Each
user is fixed at the ground. The location of the $k$th ground user is denoted
by $\mathbf{q}_{k}$, where $\mathbf{q}_{k}=[x_{k},y_{k}]$, $k\in{\cal K}$ and
${\cal K}=\left\\{1,2,\cdots,K\right\\}$. Boldface lower case letters
represent vectors, and $x_{k}$ and $y_{k}$ are the horizontal plane coordinate
of the $k$th ground user. It is assumed that the positions of users are known
to the UAV for designing trajectory. A finite time horizon with during $T$ is
considered. During the finite time, the UAV flies at a fixed altitude ($H>0$).
A block fading channel model is applied. During the finite time, the channel
is unchanged.

For ease of exposition, the finite time $T$ is discretized into $N$ equal-time
slots, denoted by $n=1,2,\cdots,N$. At the $n$th slot, it is assumed that the
horizontal plane coordinate of the UAV is
$\mathbf{q}_{u}\left[n\right]=[x_{u}[n],y_{u}[n]]$. Similar to [13]-[15], it
is assumed that the wireless channel between the UAV and each user is
dominated by the LoS channel. Thus, the channel power gain between the UAV and
the $k$th user is denoted by ${h_{k}}\left[n\right]$, given as

$\displaystyle\
{h_{k}}\left[n\right]={\beta_{0}}d_{k,n}^{-2}=\frac{{{\beta_{0}}}}{{{H^{2}}+{{\left\|{{\mathbf{q}_{u}}\left[n\right]-{\mathbf{q}_{k}}}\right\|}^{2}}}},k\in{\cal
K},n\in{\cal N},$ (1)

where $\beta_{0}$ is the channel power gain at a reference distance $d_{0}=1$
m; $d_{k,n}$ is the horizontal plane distance between the UAV and the $k$th
user at the $n$th slot, $n\in{\cal N}$, ${\cal
N}=\left\\{1,2,\cdots,N\right\\}$, and $\left\|\cdot\right\|$ denotes its
Euclidean norm. In order to reach meaningful insights into the design of an
UAV-enabled wireless powered MEC system, similar to [6]-[10], a linear EH
model is applied. Thus, the harvested energy at the $k$th user during $n$ time
slots denoted by ${E_{k}}\left[n\right]$, is given as

$\displaystyle\
{E_{k}}\left[n\right]=\sum\limits_{i=1}^{n}{\frac{{T\eta{h_{k}}\left[i\right]{P_{u}}}}{N}},$
(2)

where $\eta$ denotes the energy conservation efficiency, $0<\eta\leq 1$ and
${P_{u}}$ is the transmit power of the UAV. In this paper, the UAV employs a
constant power transmission [12]-[15]. During the $n$th slot, all users
perform energy harvesting, local computing and computation offloading. In
order for all users to offload their bits to the UAV for computation, a time
division multiple access (TDMA) protocol is applied. The time interval $T/N$
is divided into $K$ time slots with duration $\lambda=T/(NK)$ and $K$ users
offload their computation bits to the UAV one by one. Similar to [6]-[9], the
received energy and the energy for transmitting the computed results of the
UAV are ignored.

Let ${{l_{k}}\left[n\right]}$ and ${{f_{k}}\left[n\right]}$ denote the number
of the offloading bits and the central processing unit (CPU) frequency
(cycles/s) of the $k$th user at the $n$th slot, respectively. Thus, the
transmit power of the $k$th user for offloading ${{l_{k}}\left[n\right]}$
computation bits denoted by $P_{k}[n]$, is given as

$\displaystyle\
{P_{k}}\left[n\right]=\frac{{\Gamma{\sigma^{2}}\left({{2^{\frac{{{l_{k}}\left[n\right]}}{{B\lambda}}}}-1}\right)}}{{{h_{k}}\left[n\right]}},$
(3)

where $B$ is the communication bandwidth and $\sigma^{2}$ denotes the noise
power at the user. In $\left(3\right)$, $\Gamma$ is a constant related to the
gap from the channel capacity owning to a practical coding and modulation
scheme. It is assumed that $\Gamma=1$ in this paper for simplicity. Let
${{f_{u}}\left[n\right]}$ denote the CPU frequency of the UAV at the $n$th
slot. According to [6]-[9], the energy consumed for the local computation at
the $k$th user and that for the offloading computation at the UAV in the $n$th
slot are denoted by ${E_{k,l}}\left[n\right]$ and ${E_{u,o}}\left[n\right]$,
respectively given as

$\displaystyle{E_{k,l}}\left[n\right]={\gamma_{c}}K\lambda{\left[{{f_{k}}\left[n\right]}\right]^{3}},$
(4a)
$\displaystyle{E_{u,o}}\left[n\right]={\gamma_{c}}K\lambda{\left[{{f_{u}}\left[n\right]}\right]^{3}},$
(4b)

where ${\gamma_{c}}$ is the effective switched capacitance of the CPU. Similar
to the works in [14] and [16], the propulsion energy consumption model at the
UAV due to the flying in the $n$th slot denoted by $E_{s}\left[n\right]$, is
given as

$\displaystyle{E_{s}}\left[n\right]={\kappa{{\left\|{{v_{u}}\left[n\right]}\right\|}^{2}}},$
(5a)
$\displaystyle{{v_{u}}\left[n\right]}=\frac{{\left\|{{\mathbf{q}_{u}}\left[{n+1}\right]-{\mathbf{q}_{u}}\left[n\right]}\right\|}}{{K\lambda}},$
(5b)

where $\kappa=0.5WT/N$ and $W$ is the mass of the UAV. Note that the
propulsion energy consumption model employed in this paper only depends on the
velocity. In future work, we will exploit a more general model that considers
both the velocity and acceleration of the UAV.

##  III Energy Minimization Design

###  III-A The Energy Minimization Problem Formulation

In the UAV-enabled wireless powered MEC system, in order to minimize the
energy consumed at the UAV while guaranteeing the computation bits of all
users, the number of the offloading computation bits and the CPU frequency of
the users, the CPU frequency and the trajectory of the UAV are jointly
optimized. The energy minimization problem can be formulated as
$\text{P}_{{1}}$, given as

$\displaystyle\text{P}_{{1}}:{\mathop{\min}\limits_{\Xi}}\
{\sum\limits_{n=1}^{N}{{\kappa{{\left\|{{v_{u}}\left[n\right]}\right\|}^{2}}}}+{{T{P_{u}}}}+\sum\limits_{n=2}^{N}{{\gamma_{c}}K\lambda{{\left[{{f_{u}}\left[n\right]}\right]}^{3}}}},$
(6a) $\displaystyle\text{s.t.}\ \ C1:\sum\limits_{n=1}^{N}{\frac{{\lambda
K{f_{k}}\left[n\right]}}{M}}+\sum\limits_{n=1}^{N-1}{{l_{k}}\left[n\right]}={R_{k}},k\in{\cal
K},$ (6b) $\displaystyle
C2:\sum\limits_{i=1}^{n}{{E_{k,l}}\left[i\right]}+\sum\limits_{i=1}^{n}{\lambda{P_{k}}\left[i\right]}\leq\sum\limits_{i=1}^{n}{{E_{k}}\left[i\right]},k\in{\cal
K},n\in{\cal N},$ (6c) $\displaystyle
C3:\sum\limits_{j=2}^{n}{\frac{{{f_{u}}\left[j\right]K\lambda}}{M}}\leq\sum\limits_{k=1}^{K}{\sum\limits_{i=1}^{n-1}{{l_{k}}\left[i\right]}},k\in{\cal
K},n\in{\cal N}_{-N},$ (6d) $\displaystyle
C4:\sum\limits_{j=2}^{N}{\frac{{{f_{u}}\left[j\right]K\lambda}}{M}}=\sum\limits_{k=1}^{K}{\sum\limits_{i=1}^{N-1}{{l_{k}}\left[i\right]}},k\in{\cal
K},$ (6e) $\displaystyle
C5:{l_{k}}\left[N\right]=0,{f_{u}}\left[1\right]=0,k\in{\cal K},$ (6f)
$\displaystyle
C6:{\left\|{\mathbf{q}_{u}\left[{n+1}\right]-\mathbf{q}_{u}\left[n\right]}\right\|}\leq{V_{\max}}K\lambda,n\in{\cal
N},$ (6g) $\displaystyle
C7:\mathbf{q}_{u}\left[1\right]={\mathbf{q}_{0}},\mathbf{q}_{u}\left[{N+1}\right]={\mathbf{q}_{F}},$
(6h) $\displaystyle C8:f_{u}\left[n\right]\geq 0,f_{k}\left[{n}\right]\geq
0,k\in{\cal K},n\in{\cal N},$ (6i)

where $\Xi$ denotes the variable set consisting of
${f_{u}}\left[n\right],\mathbf{q}_{u}\left[n\right],{l_{k}}\left[n\right],{f_{k}}\left[n\right]$;
${R_{k}}$ denotes the total number of the computation bits of the $k$th user;
$V_{\max}$ is the maximum flying speed of the UAV; ${\mathbf{q}_{0}}$ and
${\mathbf{q}_{F}}$ are the destined initial and final locations of the UAV,
respectively; $M$ denotes the number of CPU cycles required for computing one
bit at the user and the UAV. ${\cal N}_{-N}$ denotes the set ${\cal N}$ other
than $N$. The constraint $C1$ is the total computation bits required at the
$k$th user; $C2$ is the energy causal constraint that the energy consumed for
the local computation and offloading computation bits cannot be higher than
the harvesting energy; $C3$ represents that the number of the computation bits
at the UAV in the $n$th slots cannot be higher than the total number of the
offloading computation bits of all users before the $n-1$th slot. Note that
the UAV starts to compute the offloading bits at the $n$th slots only when all
users finish offloading the computation bits of the $n-1$th slot; $C4$ denotes
that all the offloading computation bits of users should be computed; $C5$
represents that the UAV does not execute the computation task in the first
slot and all users do not offload their computation tasks in the last slot;
$C6$ is the flying speed constraint and $C7$ is the initial and final
locations constraint related to the UAV.

It is challenging to solve the non-convex problem $\text{P}_{{1}}$ due to the
presence of the couple among the optimization variables
$\mathbf{q}_{u}\left[n\right]$, ${l_{k}}\left[n\right]$ and
${f_{k}}\left[n\right]$. An alternative algorithm is proposed to solve
$\text{P}_{{1}}$ in the following subsection.

###  III-B Computation Offloading And CPU Frequency Optimization

It is seen from $\text{P}_{{1}}$ that $\text{P}_{{1}}$ is convex for a given
trajectory $\mathbf{q}_{u}\left[n\right]$. Thus, for a given
$\mathbf{q}_{u}\left[n\right]$, $\text{P}_{{1}}$ can be transformed as
$\text{P}_{{2}}$, given as

$\displaystyle\text{P}_{{2}}:\
{\mathop{\min}\limits_{{\scriptstyle{f_{u}}\left[n\right],\scriptstyle{l_{k}}\left[n\right],{f_{k}}\left[n\right]\hfill}}}\
{\sum\limits_{n=2}^{N}{{\gamma_{c}}K\lambda{{\left[{{f_{u}}\left[n\right]}\right]}^{3}}}},$
(7a) $\displaystyle\text{s.t.}\ \ C1-C5\ \text{and}\ C8.$ (7b)

Since $\text{P}_{{2}}$ is convex, it can be solved by using the Lagrange
duality method [17], [18]. By solving $\text{P}_{{2}}$, Theorem 1 can be
obtained as follows.

###### Theorem 1

For a given trajectory $\mathbf{q}_{u}\left[n\right]$, the optimal offloading
computation bits and the CPU frequency of the users, and the CPU frequency of
the UAV denoted by ${l_{k}}^{opt}\left[n\right],{f_{u}}^{opt}\left[n\right]$
and ${f_{k}}^{opt}\left[n\right]$, can be respectively given as

$\displaystyle{l_{k}}^{opt}\left[n\right]{\rm{=}}{B\lambda{\rm{lo}}{{\rm{g}}_{2}}\left\\{{\frac{{B{h_{k}}\left[n\right]\left[{\sum\limits_{j=n{\rm{+}}1}^{N{\rm{-}}1}{{\theta_{j}}}+{\mu_{k}}-{\theta_{N}}}\right]}}{{\sum\limits_{j=n}^{N}{{\nu_{k,j}}}\Gamma{\sigma^{2}}\ln
2}}}\right\\}},$ (8a)
$\displaystyle{f_{u}}^{opt}\left[n\right]=\left\\{\begin{array}[]{l}\begin{array}[]{*{20}{c}}{0,}&{n=1}\end{array}\\\
\begin{array}[]{*{20}{c}}{\sqrt{\frac{{{\theta_{N}}-\sum\limits_{j=n}^{N-1}{{\theta_{j}}}}}{{{3\gamma_{c}}M}},}}&{n=2,\cdots,N-1}\end{array}\\\
\begin{array}[]{*{20}{c}}{\sqrt{\frac{{{\theta_{N}}}}{{{3\gamma_{c}}M}},}}&{n=N}\end{array}\end{array}\right.$
(8h)
$\displaystyle{f_{k}}^{opt}\left[n\right]=\sqrt{\frac{{{\mu_{k}}}}{{3{\gamma_{c}}MK\sum\limits_{j=n}^{N}{{\nu_{k,j}}}}}},k\in{\cal
K},n\in{\cal N}$ (8i)

where $\mu_{k}\geq 0$, $\nu_{k,n}\geq 0$ and $\theta_{n}\geq 0$ are the dual
variables associated with the constraints $C1$, $C2$, $C3$ and $C4$,
respectively.

###### Proof:

See Appendix A. ∎

###### Remark 1

Theorem 1 indicates that the CPU frequency of the UAV increases with the time
slots since $\theta_{N}>0$ and $\theta_{n}\geq 0$ when $n=2,3,\cdots,N-1$. It
means that the total number of the offloading computation bits increases with
the time slots. Thus, in order to decrease the total energy consumed at the
UAV, users need to allocate a high energy for local computation so that the
number of the offloading computation can be decreased. It is also seen that
the number of the offloading computation bits is increased when the channel
condition between the UAV and users is improved. This indicates that the
number of the offloading computation bits of users increases with the decrease
of the distance between the user and the UAV. Finally, the dual variables can
be obtained by using the subgradient algorithm [19].

###  III-C Trajectory Optimization

For any given number of the offloading computation bits, the CPU frequencies
of users and the UAV, the trajectory optimization problem can be formulated as
$\text{P}_{{3}}$, given as

$\displaystyle\text{P}_{{3}}:\
{\mathop{\min}\limits_{\mathbf{q}_{u}\left[n\right]}}\
{\sum\limits_{n=1}^{N}{\kappa{{\left\|{{v_{u}}\left[n\right]}\right\|}^{2}}}}$
(9a) $\displaystyle\text{s.t.}\ \ C2,C6\ \text{and}\ C7.$ (9b)

Due to the constraint $C2$, $\text{P}_{{3}}$ is non-convex. In order to tackle
$C2$, the SCA technique is exploited. It can guarantee that the obtained
solutions satisfy the Karush-Kuhn-Tucker (KKT) conditions of $\text{P}_{{3}}$.
By using the SCA technique, Theorem 2 is given as follows.

###### Theorem 2

For any local trajectory $\mathbf{q}_{u}^{j}\left[n\right],n\in{\cal N}$ at
the $j$th iteration, one has

$\displaystyle\sum\limits_{i=1}^{n}{\frac{{K\lambda\eta{P_{u}}{\beta_{0}}}}{{{H^{2}}+{{\left\|{{\mathbf{q}_{u}}\left[i\right]-{\mathbf{q}_{k}}}\right\|}^{2}}}}}\geq
K\lambda\eta{P_{u}}{\beta_{0}}\overline{{h_{k}}}\left[n\right],$ (10a)
$\displaystyle\overline{{h_{k}}}\left[n\right]=\sum\limits_{i=1}^{n}\left\\{{\frac{{{H^{2}}+2{{\left\|{\mathbf{q}_{u}^{j}\left[i\right]-{\mathbf{q}_{k}}}\right\|}^{2}}-{{\left\|{{\mathbf{q}_{u}}\left[i\right]-{\mathbf{q}_{k}}}\right\|}^{2}}}}{{{{\left({{H^{2}}+{{\left\|{\mathbf{q}_{u}^{j}\left[i\right]-{\mathbf{q}_{k}}}\right\|}^{2}}}\right)}^{2}}}}}\right\\},$
(10b)

where the equality holds when
$\mathbf{q}_{u}\left[n\right]=\mathbf{q}_{u}^{j}\left[n\right]$.

###### Proof:

Let $f\left(z\right)=\frac{a}{{b+z}}$, where $a$ and $b$ are positive
constants, and $z\geq 0$. Since $f\left(z\right)$ is convex with respect to
$z$, the following inequality can be obtain:

$\displaystyle\
\frac{a}{{b+z}}\geq\frac{a}{{b+{z_{0}}}}-\frac{a}{{{{\left({b+{z_{0}}}\right)}^{2}}}}\left({z-{z_{0}}}\right),$
(11)

where $z_{0}$ is a given local point. By using eq. $\left(11\right)$, Theorem
2 is obtained. ∎

Using Theorem 2, $\text{P}_{{3}}$ can be solved by iteratively solving the
approximate problem $\text{P}_{{4}}$, given as

$\displaystyle\text{P}_{{4}}:\
{\mathop{\min}\limits_{\mathbf{q}_{u}\left[n\right]}}\
{\sum\limits_{n=1}^{N}{\kappa{{\left\|{{v_{u}}\left[n\right]}\right\|}^{2}}}},$
(12a) $\displaystyle\text{s.t.}\ \ \ C6\ \text{and}\ C7,$ (12b)
$\displaystyle\sum\limits_{i=1}^{n}{{E_{k,l}}\left[i\right]}+\sum\limits_{i=1}^{n}{\lambda{P_{k}}\left[i\right]}\leq
K\lambda\eta{P_{u}}{\beta_{0}}\overline{{h_{k}}}\left[n\right],k\in{\cal
K},n\in{\cal N}.$ (12c)

It is seen that $\text{P}_{{4}}$ is convex and can be readily solved by using
the software CVX [10]. Based on solving $\text{P}_{{2}}$ and $\text{P}_{{4}}$,
an alternative optimization algorithm denoted by Algorithm 1 is given to solve
$\text{P}_{{1}}$. The details for Algorithm 1 can be found in Table 1. In
Table 1, $E_{u}^{i}$ denotes the value of the objective function of
$\text{P}_{{1}}$.

$\displaystyle\ L\left({{\Xi_{1}}}\right)=$
$\displaystyle\sum\limits_{n=2}^{N}{{\gamma_{c}}K\lambda{{\left[{{f_{u}}\left[n\right]}\right]}^{3}}}+\sum\limits_{k=1}^{K}{{\mu_{k}}}\left[{{R_{k}}-\sum\limits_{n=1}^{N}{\frac{{\lambda{f_{k}}\left[n\right]}}{M}}-\sum\limits_{n=1}^{N-1}{{l_{k}}\left[n\right]}}\right]$
$\displaystyle+\sum\limits_{k=1}^{K}{\sum\limits_{n=1}^{N}{{\nu_{k,n}}\left\\{{\sum\limits_{i=1}^{n}{{E_{k,l}}}\left[i\right]+\sum\limits_{i=1}^{n}{\lambda{P_{k}}\left[i\right]}-\sum\limits_{i=1}^{n}{{E_{k}}\left[i\right]}}\right\\}}}+\sum\limits_{n=1}^{N-1}{{\theta_{n}}\left\\{{\sum\limits_{j=2}^{n}{\frac{{{f_{u}}\left[j\right]K\lambda}}{M}}-\sum\limits_{k=1}^{K}{\sum\limits_{i=1}^{n-1}{{l_{k}}\left[i\right]}}}\right\\}}$
$\displaystyle+{\theta_{N}}\left\\{{\sum\limits_{k=1}^{K}{\sum\limits_{i=1}^{N-1}{{l_{k}}\left[i\right]}-\sum\limits_{j=2}^{N}{\frac{{{f_{u}}\left[j\right]K\lambda}}{M}}}}\right\\}+\sum\limits_{k=1}^{K}{{\rho_{k}}{l_{k}}\left[N\right]}+\vartheta{f_{u}}\left[1\right]$
(13)
TABLE I: The alternative optimization algorithm Algorithm 1: The alternative
optimization algorithm for $\text{P}_{{1}}$
---
1: Setting:
$R_{k}$, $k\in{\cal K}$, $P_{u}$, $T$, $N$, ${V_{\max}}$, $\mathbf{q}_{0}$,
$\mathbf{q}_{F}$, and the tolerance errors $\xi$, $\xi_{1}$;
2: Initialization:
The iterative number $i=1$ and $\mathbf{q}_{u}^{i}\left[n\right]$;
3: Repeat 1:
calculate ${l_{k}}^{opt,i}\left[n\right],{f_{u}}^{opt,i}\left[n\right]$ and
${f_{k}}^{opt,i}\left[n\right]$
using eq. $\left(8\right)$ for given $\mathbf{q}_{u}^{i}\left[n\right]$;
update $\mu_{k}$, $\nu_{k,n}$ and $\theta_{n}$ using the subgradient
algorithm;
initialize the iterative number $j=1$;
Repeat 2:
solve $\text{P}_{\textbf{4}}$ by using CVX for the given
${l_{k}}^{opt,i}\left[n\right]$,
${f_{u}}^{opt,i}\left[n\right]$ and ${f_{k}}^{opt,i}\left[n\right]$;
update $j=j+1$, and $\mathbf{q}_{u}^{j}\left[n\right]$;
if
$\sum\limits_{n=1}^{N}{\left\|{\mathbf{q}_{u}^{j}\left[n\right]-\mathbf{q}_{u}^{j}\left[n\right]}\right\|}\leq\xi$
$\mathbf{q}_{u}^{i}\left[n\right]=\mathbf{q}_{u}^{j}\left[n\right]$ ;
break;
end
end Repeat 2
update the iterative number $i=i+1$;
if $\left|E_{u}^{i}-E_{u}^{i-1}\right|\leq\xi_{1}$
break;
end
end Repeat 1
4: Obtain solutions:
${l_{k}}^{opt}\left[n\right]$, ${f_{u}}^{opt}\left[n\right]$ and
${f_{k}}^{opt}\left[n\right]$ and $\mathbf{q}_{u}^{opt}\left[n\right]$.

##  IV Simulation Results

In this section, simulation results are presented to compare the performance
obtained by using our proposed design with that achieved by using two
benchmark schemes, denoted by Scheme 1 and Scheme 2, respectively. In Scheme
1, the UAV flies straight with a constant speed from the initial position to
the final position. In Scheme 2, the UAV flies along the trajectory that is a
semi-circle with its diameter being
${\left\|{\mathbf{q}_{F}-{\mathbf{q}_{0}}}\right\|}$. The converge performance
of the proposed algorithm is also evaluated by simulation results. The
simulation settings are based on the works in [9] and [14]. The positions and
the total number of the computation bits of users are set as:
$\mathbf{q}_{1}=[0,0]$, $\mathbf{q}_{2}=[0,10]$, $\mathbf{q}_{3}=[10,10]$,
$\mathbf{q}_{4}=[10,0]$, $R_{1}=2$ Mbits, $R_{2}=4$ Mbits, $R_{3}=6$ Mbits,
and $R_{4}=3$ Mbits, respectively. The detail settings are given in Table II.

TABLE II: Simulation Parameters Parameters | Notation | Typical Values
---|---|---
Numbers of Users | $K$ | $4$
The height of the UAV | $H$ |  $10$ m
The time length of the UAV flying | $T$ |  $2$ sec
Numbers of CPU cycles | $M$ |  $10^{3}$ cycles/bit
Energy conversation efficiency | $\eta$ | $0.8$
Communication bandwidth | $B$ |  $40$ MHz
The receiver noise power | $\sigma^{2}$ |  $10^{-9}$ W
The number of time slots | $N$ | $50$
The mass of the UAV | $W$ |  $9.65$ kg
The effective switched capacitance | $\gamma_{c}$ | $10^{-28}$
The channel power gain | $\beta_{0}$ |  $-50$ dB
The tolerance error | $\xi,\xi_{1}$ | $10^{-4}$
The initial position of the UAV | $\mathbf{q}_{0}$ | $[0,0]$
The final position of the UAV | $\mathbf{q}_{F}$ | $[10,0]$
The maximum speed of the UAV | $V_{\max}$ |  $10$ m/s
The transmit power of the UAV | $P_{u}$ |  $100$ dBm
Figure 2: The trajectories of the UAV under different schemes with $T=2$
seconds.

Fig. 2 shows the trajectories of the UAV under different schemes. The time
length of the UAV flying is set as $T=2$ seconds. The trajectories of the UAV
under Scheme 1 and Scheme 2 are also presented. As shown in Fig. 2, under our
proposed optimal trajectory, the UAV firstly flies smoothly and tends to User
2 and User 3, and then the UAV flies smoothly with a higher speed to the final
position. The reason is that the UAV needs to provide more energy to User 2
and User 3, which has a larger number of computation bits to be offloaded.
Moreover, in order to control the total number of the computation bits of all
users offloaded to the UAV, the UAV flies with a higher speed in the end of
the flying time so that the harvested energy of users used for offloading the
computation bits can be compromised.

Figure 3: The total energy consumed at the UAV versus the time length of the
UAV flying under different schemes.

Fig. 3 shows the total energy consumed at the UAV versus the time length of
the UAV flying under different schemes. It is seen that the energy consumed at
the UAV by using our proposed scheme is the smallest among those by using the
benchmark schemes. This demonstrates that our proposed scheme that jointly
optimizes the number of the offloading computation bits, the CPU frequency of
users and the UAV, and the trajectory of the UAV can is more efficient in
terms of the energy minimization of the UAV. It is also seen that the total
energy consumed at the UAV decreases with the increase of the time length of
the UAV flying, irrespective of the adopted scheme. It can be explained by the
fact that the total energy consumed at the UAV is dominated by the flying
speed and the CPU frequency of the UAV, and the flying speed and the CPU
frequency can be decreased when the flying time is increased.

Figure 4: The total energy consumed at the UAV versus the number of iterations
$T=2,2.2$ or $2.4$ seconds.

Fig. 4 is presented to verify the efficiency of our proposed alternative
algorithm. It can be seen that only several number of iterations are required
for Algorithm 1 to converge.

##  V Conclusion

An UAV-enabled wireless powered MEC system was studied where the UAV provides
multiple ground users with computation offloading and sustainable operation
opportunities. The number of the offloading computation bits and the CPU
frequency of users, the CPU frequency and the trajectory of the UAV were
jointly optimized in order to minimize the energy consumed at the UAV. An
alternative algorithm was proposed based on the SCA techniques. Simulation
results show that our proposed design outperforms other benchmark schemes and
that the proposed algorithm only requires several number of iterations to
converge.

##  Appendix A Proof of Theorem 1

The Lagrangian of $\text{P}_{{2}}$ related to the proof is given by eq.
$\left(13\right)$ at the top of the previous page, where ${\rho_{k}}\geq 0$
and $\vartheta\geq 0$ are the dual variables related to the constraint $C5$;
$\Xi_{1}$ is the set consisting of all optimization and dual variables. Thus,
the derivations of the Lagrangian of $\text{P}_{{2}}$ with respect to
${l_{k}}\left[n\right]$ and ${f_{k}}\left[n\right]$, can be respectively given
as

$\displaystyle\frac{{\partial
L\left(\Xi\right)}}{{\partial{l_{k}}\left[n\right]}}$
$\displaystyle=\frac{{\Gamma{\sigma^{2}}\left({{2^{\frac{{{l_{k}}\left[n\right]}}{{B\lambda}}}}}\right)\ln
2}}{{B{h_{k}}\left[n\right]}}\sum\limits_{j=n}^{N}{{\nu_{k,j}}}-\left[{\sum\limits_{j=n{\rm{+}}1}^{N-1}{{\theta_{j}}}+{\mu_{k}}-{\theta_{N}}}\right],$
(14a) $\displaystyle\frac{{\partial
L\left(\Xi\right)}}{{\partial{f_{k}}\left[n\right]}}=-\frac{{{\mu_{k}}\lambda}}{M}+3{\gamma_{c}}K\lambda\sum\limits_{j=n}^{N}{{\nu_{k,j}}{{\left[{{f_{k}}\left[n\right]}\right]}^{2}}}.$
(14b)

Let their derivations be zero. Thus, eq. $\left(8\rm{a}\right)$ and eq.
$\left(8\rm{c}\right)$ are obtained. Let the derivation of the Lagrangian of
$\text{P}_{{2}}$ with respect to ${f_{u}}\left[n\right]$ be zero. One has

$\displaystyle
3{\gamma_{c}}{\left[{{f_{u}}\left[n\right]}\right]^{2}}{\rm{+}}\frac{{\left({\sum\limits_{j=n}^{N-1}{{\theta_{j}}}-{\theta_{N}}}\right)}}{M}{\rm{=}}0,n=2,\cdots,N-1,$
(15a) $\displaystyle
3{\gamma_{c}}{\left[{{f_{u}}\left[N\right]}\right]^{2}}{\rm{-}}\frac{{{\theta_{N}}}}{M}{\rm{=}}0,n=N.$
(15b)

From eq. $\left(15\right)$, eq. $\left(8\rm{b}\right)$ is obtained. The proof
for Theorem 1 is completed.

## References

  * [1] L. Wei, R. Q. Hu, Y. Qian, G. Wu, “Enabling device-to-device communications underlaying cellular networks: challenges and research aspects,”  _IEEE Commun. Mag._ , vol.52, no.6, pp.90-96, June 2014.
  * [2] F. Zhou, Y. Wu, R. Q. Hu, Y. Wang, and K. K. Wong, “Energy-efficient NOMA enabled heterogeneous cloud radio access networks,”  _IEEE Network_ , to be published, 2017.
  * [3] R. Q. Hu and Y. Qian, “An energy efficient and spectrum efficient wireless heterogeneous network framework for 5G systems,”  _IEEE Commun. Mag._ , vol.52, no.5, pp.94-101, May 2014.
  * [4] L. Wei, R. Q. Hu, Y. Qian, G. Wu, “Key elements to enable millimeter wave communications for 5G wireless systems,”  _IEEE Wireless Commun. Mag._ , vol.21, no.6, pp. 136-143, Dec. 2014.
  * [5] C. Wang, _et al._ , “Computation offloading and resource allocation in wireless cellular networks with moible edge computing,”  _IEEE Trans. Wireless Commun._ , vol. 16, no. 8, pp. 4924-4938, Aug. 2017.
  * [6] C. You, K. Huang, and H. Chae, “Energy efficient mobile cloud computing powered by wireless energy transfer,”  _IEEE J. Sel. Areas Commun._ , vol. 34, no. 5, pp. 1757-1771, May, 2016.
  * [7] J. Xu, L. Chen, and S. Ren, “Online learning for offloading and autoscaling in energy harvesting mobile edge computing,”  _IEEE Trans. Cogn. Commun. Netw._ , vol. 3, no. 3, pp. 361-373, March, 2017.
  * [8] S. Bi and Y. Zhang, “Computation rate maximization for wireless powered mobile egde computing with binary computation offloading,” submitted to _IEEE Trans. Wireless Commun._ , https://arxiv.org/abs/1708.08810.
  * [9] F. Wang, J. Xu, X. Wang, and S. Cui, “Joint offloading and computing optimization in wireless powered mobile-edge computing systems,” submitted to _IEEE Trans. Wireless Commun._ , https://arxiv.org/abs/1702.00606.
  * [10] F. Zhou, _et al._ , “Robust AN-aided beamforming and power splitting design for secure MISO cognitive radio with SWIPT,”  _IEEE Trans. Wireless Commun._ , vol. 16, no. 4, pp. 2450-2464, April 2017.
  * [11] X. Lu, P. Wang, D. Niyato, D. I. Kim, and Z. Han, “Wireless networks with RF energy harvesting: A contemporary survey,”  _IEEE Commun. Surveys Tuts._ , vol. 17, pp. 757-789, Second Quarter, 2015.
  * [12] J. Xu, Y. Zeng, and R. Zhang, “UAV-enabled wireless power transfer: Trajectory design and energy optimization,” submitted to _IEEE Trans. Wireless Commun._ , https://arxiv.org/abs/1709.07590.
  * [13] J. Xu, Y. Zeng, and R. Zhang, “UAV-enabled wireless power transfer: Trajectory design and energy region charaterization,” in _Proc. IEEE Global Commun. Conf._ Singapore, 2017,
  * [14] S. Jeong, O. Simeone, and J. Kang, “Mobile edge computing via a UAV-mounted cloudlet: Optimization of bit allocation and path planning,”  _IEEE Trans. Vehicular Technol._ , to be published, 2017.
  * [15] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,”  _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 3747-3760, June 2017.
  * [16] N. Xue, “Design and optimization of lithiumion batteries for electricvehicle applications,” 2014, Doctoral dissertation, University of Michigan.
  * [17] S. P. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge, U.K.: Cambridge Univ. Press, 2004.
  * [18] Q. Li, Y. Xu, R. Q. Hu, and Y. Qian, “Optimal fractional frequency reuse and power control in the heterogeneous wireless networks,”  _IEEE Trans. Wireless Commun._ , vol.12, no.6, pp. 2658- 2668, June 2013.
  * [19] F. Zhou, N. C. Beaulieu, Z. Li, J. Si, and P. Qi, “Energy-efficient optimal power allocation for fading cognitive radio channels: Ergodic capacity, outage capacity and minimum-rate capacity,”  _IEEE Trans. Wireless Commun._ , vol. 15, no. 4, pp. 2741-2755, Apr. 2016.