# Iterative Training: Finding Binary Weight Deep Neural Networks with Layer
Binarization

Cheng-Chou Lan
Rakuten Mobile, Inc.
<EMAIL_ADDRESS>

###### Abstract

In low-latency or mobile applications, lower computation complexity, lower
memory footprint and better energy efficiency are desired. Many prior works
address this need by removing redundant parameters. Parameter quantization
replaces floating-point arithmetic with lower precision fixed-point
arithmetic, further reducing complexity.

Typical training of quantized weight neural networks starts from fully
quantized weights. Quantization creates random noise. As a way to compensate
for this noise, during training, we propose to quantize some weights while
keeping others in floating-point precision. A deep neural network has many
layers. To arrive at a fully quantized weight network, we start from one
quantized layer and then quantize more and more layers. We show that the order
of layer quantization affects accuracies. Order count is large for deep neural
networks. A sensitivity pre-training is proposed to guide the layer
quantization order.

Recent work in weight binarization replaces weight-input matrix multiplication
with additions. We apply the proposed iterative training to weight
binarization. Our experiments cover fully connected and convolutional networks
on MNIST, CIFAR-10 and ImageNet datasets. We show empirically that, starting
from partial binary weights instead of from fully binary ones, training
reaches fully binary weight networks with better accuracies for larger and
deeper networks. Layer binarization in the forward order results in better
accuracies. Guided layer binarization can further improve that. The
improvements come at a cost of longer training time.

##  1 Introduction

Recent works using deep convolutional networks have been successfully applied
to a large variety of computer vision tasks, such as image recognition (He et
al., 2016), object segmentation (He et al., 2017) and scene segmentation (Chen
et al., 2018). These networks are large. For example, ResNet-152 has 60.2
million parameters (Zagoruyko and Komodakis, 2016) and requires 11.3 billion
FLOPs (He et al., 2016). A large number of parameters results in a large
memory footprint. At 32-bit floating-point precision, 229.64 MB is needed to
store the ResNet-152 parameter values.

In low-latency or mobile applications, lower computation complexity, lower
memory footprint and better energy efficiency are desired. Many prior works
address this need of lower computation complexity. In a survey paper Cheng et
al. (2018), efficient computation of neural networks is organized into four
categories: network pruning, low-rank decomposition, teacher-student network
and network quantization.

Network pruning removes redundant parameters which are not sensitive to
performance. Low-rank decomposition uses matrix or tensor decomposition
methods to reduce number of parameters. In teacher-student network, knowledge
transfer is exploited to train a smaller student network using a bigger
teacher network. In these three categories, their common theme is a reduction
of number of parameters. During forward propagation, one of the most
computationally intensive operation in a neural network is the matrix
multiplication of parameters with input. With reduced parameters, FLOPs and
memory footprint reduce. With reduced FLOPs, energy efficiency improves.

In the fore-mentioned categories, network parameters typically use floating-
point precision. In the last category, network quantization, the parameters
and, for some works, all computations are quantized. For many low-latency or
mobile applications, we typically train offline and deploy pre-trained models.
Thus, the main goal is the efficiency in forward propagation. It is desirable
to compute backward propagation and parameter updates in floating-point
precision. The seminal work Courbariaux et al. (2015) matches our scope. They
quantize network weights to binary values, e.g. -1.0 and 1.0, while also
keeping weight values in floating-point precision for backward propagation.
During forward propagation, instead of matrix multiplication of weights with
input, the sign of these binary weights specify addition or subtraction of
inputs. Memory footprint is dramatically reduced to 1 bit per weight. Energy
efficiency improves because addition is more energy efficient than
multiplication Horowitz (2014).

Prior works in network quantization (Courbariaux et al., 2015; Li et al.,
2016; Hubara et al., 2016; Zhou et al., 2016; Wu et al., 2018) typically start
training from quantizing all weights in the network. Quantization creates
error which is the difference between the original value and its quantized
value. In other words, actual weight value, $w_{q}$, is

$w_{q}=w-w_{error}$ (1)

To reduce the impact of $w_{error}$, we hypothesize that if we quantize some
weights while leaving others in floating-point precision, the latter ones
would be able to compensate for the error introduced by quantization. To reach
a fully quantized network, we propose an iterative training, where we
gradually quantize more and more weights. This raises two questions. First,
how to choose the grouping of weights to quantize together at each iteration.
Second, how to choose the quantization order across groups. A feedforward,
deep neural network has many layers. One natural grouping choice is one group
per layer. For the quantization order of groups, we propose a sensitivity pre-
training to choose the order. A random order and other obvious orders are
chosen as comparison.

### Contributions.

  * •

We propose an iterative training regime that gradually finds a full binary
weight network starting from an initial partial binary weight network.

  * •

We demonstrate empirically that starting from a partial binary weight network
result in higher accuracy than starting from a full binary weight one.

  * •

We demonstrate empirically that the forward order is best, compared to other
obvious orders. In addition, sensitivity pre-training can further improve
that.

  * •

Code is available at https://github.com/rakutentech/iterative_training.

In the sections that follow, we describe the iterative training algorithm in
detail. Next, we present the iterative training of fully connected networks
using the MNIST dataset (Lecun et al., 1998) and of convolutional neural
networks using the CIFAR-10 (Krizhevsky, 2009) and ImageNet (Russakovsky et
al., 2015) datasets. Then we present the sensitivity pre-training for
convolution neural networks. Finally, we discuss related work and conclusion.

##  2 Iterative Training

Algorithm 1 Iterative Training

Input: Input data and label
Parameter: Number of iterations, $N$
Parameter: Number of layers, $L$
Parameter: Total number of epochs, $T$
Parameter: $BinarizationOrder$ array, length $L$
Output: A trained neural network with binary weights

1: $BinarizationState\leftarrow zeros(L)$.

2: for $j\leftarrow 1$ to $L$ do

3: $layer\leftarrow BinarizationOrder[j]$.

4: $BinarizationState[layer]\leftarrow 1$.

5: for $i\leftarrow 1$ to $N$ do

6: BinarizeWeights($BinarizationState$).

7: ForwardPropagation().

8: BackwardPropagation().

9: UpdateParameters().

10: end for

11: end for

12: $i\leftarrow L*N$

13: while $i<T$ do

14: BinarizeWeights($BinarizationState$).

15: ForwardPropagation().

16: BackwardPropagation().

17: UpdateParameters().

18: end while

A feedforward, deep neural network has many layers, say, $L$. We study
iterative training by quantizing more and more weights layer-by-layer.
Iterative training starts from one quantized layer while all other layers are
in floating-point precision. Each iteration trains for a fixed number of
epochs, $N$. Next, we quantize the next layer and trains for another $N$
epochs. Iterative training stops when there are no more layers to quantize. In
the case of ResNet architectures, same as the original paper, we reduce
learning rate by 10 twice and continue training. Algorithm 1 illustrates the
iterative training regime. As the experiments will show, this regime
consistently finds fully quantized network with better accuracies than
starting from an initial fully quantized network (our baseline).

For quantization scheme, we follow weight binarization in Courbariaux et al.
(2015), but, for simplicity, without "tricks": no weight clipping and no
learning rate scaling. In addition, we use softmax instead of square hinge
loss. The inner for-loop in Algorithm 1 is the same as the training regime in
Courbariaux et al. (2015), except that a state variable is introduced to
control whether a layer needs binarization or not. We use the PyTorch
framework Paszke et al. (2019). ImageNet results in the biggest GPU memory
needs and longest training time, which are about 10 GB and about one day to
train one model on a Nvidia V100, respectively.

Network | 300-100-10 | 784-784-10 | Vgg-5 | Vgg-9 | ResNet-20 | ResNet-21
---|---|---|---|---|---|---
Convolutional layers |  |  | 64, 64 | 64, 64 | 16, 3x[16, 16] | 64, 4x[64]
128, 128 | 3x[32, 32] | 5x[128], 5x[256]
256, 256 | 3x[64, 64] | 5x[512]
Fully connected layers | 300, 100, 10 | 784, 784, 10 | 256, 256, 10 | 256, 256, 10 | 10 | 1000
Dataset | MNIST | MNIST | CIFAR-10 | CIFAR-10 | CIFAR-10 | ImageNet 2012
Train / Validation / Test | 55K / 5K / 10K | 55K / 5K / 10K | 45K / 5K / 10K | 45K / 5K / 10K | 45K / 5K / 10K | 1.2M / 0 / 50K
Batch size | 100 | 100 | 100 | 100 | 128 | 256
Optimizer | Adam | Adam | Adam | Adam | SGD | SGD
Momentum 0.9 | Momentum 0.9
Weight decay 1e-4 | Weight decay 1e-4
Pre-training epochs, $K$ | 150 | 150 | 200 | 450 | 300 | 20
Epochs per layer | 150 | 150 | 150 | 150 | 50 | 2
Layers, $L$ | 3 | 3 | 5 | 9 | 20 | 21
Total epochs, $T$ | 450 | 450 | 750 | 1350 | 1200 | 67
Order count | 6 | 6 | 120 | 362,880 | 2e+18 | 5e+19
Weight count | 266,200 | 1,237,152 | 4,300,992 | 2,261,184 | 268,336 | About 11e6

Table 1: Summary of network architectures and their hyper-parameters.

As shown by order count in Table 1, there is a large number of layer
binarization order for a deep neural network. In this work, we experiment with
random and obvious orders, to show that starting from a partially quantized
weight network is better than starting from fully quantized one. In a later
section, we introduce the proposed sensitivity pre-training to select a layer
binarization order.

For obvious orders, we experiment with the forward order, i.e., quantizing
layer-by-layer from input layer towards output layer and the reverse order,
i.e., from output layer towards input layer. We then compare to training when:
(1) all weights are quantized from start (baseline) (2) all weights are in
floating-point precision and stay so. As the experiments will show, for bigger
and deeper networks, the forward order consistently finds fully quantized
network with better accuracies than other orders.

In the following subsections, we discuss experimental results for fully
connected and convolutional networks.

###  2.1 Iterative Training for Fully Connected Networks

We investigate iterative training of fully connected networks with the MNIST
dataset, which has 60,000 training and 10,000 test images. We use the last
5,000 images from the training set for validation and the remaining 55,000 as
training images for all MNIST experiments. We use no data augmentation. We use
batch normalization (Ioffe and Szegedy, 2015), no drop-out and weight
initialization as He et al. (2015). We use softmax as classifier.

For iterative training, we train for 150 epochs per layer. For each network
architecture, the total number of training epochs is number of layers
multiplied by 150 epochs. Because there are three layers for the chosen
networks, all MNIST experiments are trained for 450 epochs. For all cases, we
find the best learning rate from the best error on the validation set. For
layer-by-layer binarization cases, the best error is selected from epochs when
all layers are binarized. We then use each corresponding best learning rate
for the error on the test set. We vary the seed for 5 training sessions and
report the learning curves of the average test errors in the figures. Table 1
reports other hyper-parameters.

Figure 1:  Left: Test errors for 300-100-10 network. Right: Test errors for 784-784-10 network.  Case | 300-100-10 | 784-784-10 | Improvement
---|---|---|---
Binary | 0.023 | 0.021 | 0.002
Float | 0.015 | 0.014 | 0.001
Reverse | 0.024 | 0.018 | 0.006
Forward | 0.026 | 0.016 | 0.010
Table 2: Error improvement from using 784-784-10 over 300-100-10.

For network architectures, we study the 300-100-10 network (Lecun et al.,
1998) and a bigger variant, 784-784-10. Figure 1 shows test errors for the
300-100-10 and 784-784-10 networks. The float case is training where all
weights are in floating-point precision and stay so. The binary case
(baseline) is training where all weights are binarized from the start. The
forward case is training where layer binarization is in the forward order, the
reverse in the reverse order. The solid lines are the mean across multiple
runs and the matching shaded color is one standard deviation.

For the smaller network, 300-100-10, the binary case reaches a lower error
than forward and reverse orders. Next best is the reverse order then the
forward one. This shows that order of layer binarization matters for accuracy.
On the contrary, for the bigger network, 784-784-10, the forward and reverse
cases does better than the binary one. Binarization operation is not
differentiable. According to Equation 1, it injects a random error signal into
the network. During iterative training, some of the weights are in floating-
point precision. We hypothesize that they are compensating for the random
error. At the same time, we think bigger networks are more robust due to more
parameters.

The error improvement of upgrading to a bigger network is given in Table 2.
The forward and reverse orders have significantly higher improvement than
float and binary, showing that iterative training is beneficial. In addition,
the forward order has a higher improvement than reverse. We observe the same
pattern for subsequent network architectures. Namely, for bigger and deeper
networks, starting from partial binary weight network, instead of full binary
weight network, iterative training with forward weight quantization order
finds full binary weight network with higher accuracies.

###  2.2 Iterative Training for Convolutional Networks

We investigate iterative training of convolutional networks with the CIFAR-10
dataset, which has 50,000 training and 10,000 test images. We randomly choose
5,000 images from the training set as the validation set and the remaining
45,000 as training images for all CIFAR-10 experiments. We use the same data
augmentation as He et al. (2016): 4 pixels are padded on each side, and a
32x32 crop is randomly sampled from the padded image or its horizontal flip.
We use batch normalization, no drop-out and weight initialization as He et al.
(2015). We use softmax as classifier.

We experiment with VGG (Simonyan and Zisserman, 2015) and ResNet (He et al.,
2016) architectures. For iterative training of VGG architectures, we train for
150 epochs per layer. For iterative training of ResNet-20 architecture, we
train for 50 epochs per layer. Same as the original paper, we reduce learning
rate by a factor 10 twice, once at 1000 epochs and a second time at 1100
epochs. Then stop training at 1200 epochs. Using same methodology as MNIST
experiments, for all cases, we use the validation set to tune the learning
rate and test set to report errors. Table 1 reports other hyper-parameters.

Figure 2:  Left: Test errors for Vgg-5 network. Right: Test errors for Vgg-9
network.

For VGG architecture, we study a shallower, Vgg-5, and a deeper network,
Vgg-9. As their names suggest, Vgg-5 has 5 layers and Vgg-9, 9. Figure 2 shows
test errors for Vgg-5 and Vgg-9 networks. The float case is training where all
weights are in floating-point precision and stay so. The binary case
(baseline) is training where all weights are binarized from the start. The
forward case is training where layer binarization is in the forward order, the
reverse in the reverse order and the random case, a randomly selected order.

For both network architectures, the binary case has the highest error and the
float case the lowest error. In the same pattern as the larger MNIST network,
starting from partial binary weight networks, iterative training finds full
binary weight networks that have lower error than the binary cases. For Vgg-5,
a shallower network, the ascending error ranking is reverse, forward then
random. For Vgg-9, a deeper network, the ranking is forward, random then
reverse. This shows again that layer binarization order matters.

Case | Vgg-5 | Vgg-9 | Improvement
---|---|---|---
Binary | 0.30 | 0.28 | 0.02
Float | 0.16 | 0.08 | 0.08
Reverse | 0.22 | 0.1126 | 0.1074
Forward | 0.22 | 0.1025 | 0.1175
Table 3: Error improvement from using Vgg-9 over Vgg-5.

The error improvement of upgrading to Vgg-9 from Vgg-5 is summarized in Table
3. There is a small improvement for the binary case. The float case has a
significantly higher improvement than binary. Next higher is the reverse case.
Finally, the forward case has the highest improvement. In the same pattern as
in the MNIST experiments, favoring iterative training and the forward order.

As shown in Table 1, although Vgg-9 has a smaller number of weight parameters
than Vgg-5, it has more layers. Iterative training continues to be beneficial.
We hypothesize that this is due to a more gradual rate of total binarization.
For Vgg-9, as each layer is binarized, relatively more weights stay in
floating-point precision to compensate for the random noise injected by the
binarization operation.

Figure 3:  Left: Test errors for ResNet-20 network. Right: Zoom to final
epochs. Learning rate is reduced by 10x at 1000 and again at 1100 epochs.

For an even deeper network, we study ResNet-20 from He et al. (2016), which
has 20 layers, as its name suggests. Figure 3 shows test errors for the
ResNet-20 network. The binary case has the highest error and the float case
has the lowest error. In the same pattern as other network architectures,
iterative training finds full binary weight networks that have lower error
than the binary case. In increasing error order is forward, random and
reverse. Again, showing that the order of binarization matter and the forward
order has advantage. In the next section, we propose a sensitivity pre-
training to select a binarization order.

##  3 Sensitivity Pre-training

In prior sections we demonstrated empirically that starting from a partial
binary weight network results in higher accuracy than starting from a fully
binary weight one for larger and deeper networks. In this section, we describe
the proposed sensitivity pre-training to choose a the binarization order.

For shallower neural networks like the 3-layer fully connected networks for
the MNIST dataset, exhaustive search for the best binarization order is
possible. For deeper neural networks such as Vgg-5, Vgg-9 and Resnet-20, it is
impractical to do so, as shown by order count in Table 1. However, we can
obtain a measure of error sensitivity to layer quantization. Then let the
sensitivity be a guide for binarization ordering.

Sensitivity is computed as follows. We train $L$ models, where in each model
only the weights of the L-th layer is binarized while others are in floating-
point precision. We train for $K$ epochs and, as before, use validation set to
tune the learning rate to get the best validation error for each sensitivity
model. $K$ for Vgg-5 is 200 and for Vgg-9 is 450. $K$ for ResNet-20 is 300.
For ResNet, same as the original paper, we reduce learning rate by 10 twice,
one at epoch 200 and again at epoch 250. Then we rank these $L$ best
validation errors in ascending order. This becomes the ascending layer
binarization order for iterative training.

Figure 4:  Left: Test errors for Vgg-5 network. Right: Test errors for Vgg-9
network. Ascending and descending orders are chosen by sensitivity pre-
training.

During iterative training using ascending order, the layer that had the lowest
error will be binarized first, while the layer that had the highest error
last, meaning the latter stays in floating-point precision the longest during
training. As shown in Figure 4 for Vgg-5 and Vgg-9, the ascending order
results in a fully binary weight network with the lowest error, beating the
forward ones. Also shown is the descending order, which is the reverse of the
ascending one. For both networks, the descending order results in error higher
than ascending, showing again that binarization order matters. In the case of
Vgg-5, the random order is worst while descending one follows closely behind.
In the case of Vgg-9, the descending one is the worst of all. In short, the
lower the error for one order, the higher its reverse order would be.

Figure 5:  Left: Test errors for Resnet-20 network. Right: Zoom to final
epochs. Ascending and descending orders are chosen by sensitivity pre-
training.

For ResNet-20, Figure 5 shows the test errors with ascending and descending
orders. Unlike for Vgg-5 and Vgg-9, the forward order reaches accuracy better
than both ascending and descending orders. The proposed sensitivity pre-
training considers binarization of layers independently. We hypothesize that
there may be interactions between multiple layers.

Figure 6:  Left: Test errors for ResNet-21 network. Right: Zoom to final
epochs. Ascending order are chosen by sensitivity pre-training.

For ImageNet, we experiment with ResNet-18 He et al. (2016). Since it has 21
layers, we will refer to it as ResNet-21. The optimizer is SGD with momentum
0.9 and weight decay 1e-4. For sensitivity pre-training, $K$ is 20 epochs. For
each layer, we sweep 3 learning rates and use the last-epoch errors of the
test set to choose the ascending order. In the full training, we choose 2
epochs per layer. The starting learning rate, 0.01, comes from the best
learning rate in sensitivity pre-training. Same as the orginal paper, we
reduce learning rate by 10 twice, after 42 epochs and again after 57 epochs.
We stop training after 67 epochs. The floating-point training is just one run,
while all other binarization training are from 5 random-seeded runs. Figure 6
shows the test errors with forward and ascending orders. The ascending order
has a lower mean error than forward. Both of which are better than binary.
Again, binarization order matters and ascending order is better than the
forward one.

###  3.1 Exhaustive Search

Figure 7:  Left: Test errors for 300-100-10 network. Right: Test errors for
784-784-10 network. 132 means binarization order is layer 1, layer 3 then
layer 2. Forward order is 123.

For shallower neural networks like the 3-layer fully connected network for the
MNIST dataset, exhaustive search for the best binarization order is possible.
Figure 7 shows result for all combinations of layer binarization order for
300-100-10 and 784-784-10 networks. For the former, a smaller network, the
ascending order turns out to be same as reverse. Errors for all combinations
are very close. The best order is not the ascending one, but 132 and 312, both
of which are better than binary by a small margin. 132 means binarization
order is layer 1, layer 3 then layer 2. Thus, also for 300-100-10, starting
from partial weight binarization is better than from full weight binarization.

For the bigger network, 784-784-10, the ascending order is better than forward
and reverse ones. The descending order is worst of all others. This is
consistent with the results from convolutional networks. Here, the ascending
one shares with another in best accuracy.

In summary, we proposed using sensitivity pre-training as a guide for layer
binarization order. For 784-784-10, Vgg-5, Vgg-9 and ResNet-21, we have shown
empircally that better accuracies are achieved. This improvement comes at a
cost of pre-training additional $L$ models.

##  4 Related Work

Our work introduces an iterative layer-by-layer quantization training regime.
Although we demonstrated the results using weight binarization, the regime is
independent of quantization schemes. We think other schemes, e.g., Li et al.
(2016) (where weights are ternary: -1.0, 0 and 1.0), may yield similar trends.

Hu et al. (2018) transforms weight binarization as a hashing problem. Same as
ours, their iterative method also operates layer-by-layer, from input layer
towards output layer. However, they start from a pre-trained network and,
after weight quantization without fine-tuning, fine-tune the biases. Ours
starts from an untrained network and gradually trains a full binary weight
network, which we believe allows the network to adapt to the random noise
created by the quantization operation. In addition, their final weights are
not pure binary, but power-of-2 multiples. When constrained to pure binary,
they report non-convergence. Our iterative training does not require pure
binary weights. For future work, we can binarize using power-of-2 multiples.

Zhou et al. (2017) iterates over both pruning and quantization techniques.
First, weights are partitioned into two groups. Then, weights in first group
are quantized to power-of-2 multiples or zero. Next, weights in the second
groups are fine-tuned, while the first group receives no parameter updates. In
the next iteration, some of weights in second group is assigned to the first
group. The process is repeated until all weights are members of the first
group. In this partitioning scheme, the first group contains weights from all
layers. It is possible to merge both methods because their partitioning is
orthogonal to ours. Once weights join the first group, their values stay
unchanged for the rest of the fine-tuning. Because our binarization is based
on (Courbariaux et al., 2015), floating-point weights prior to quantization
are saved for parameter updates. Thus, during iterative training of later
layers, weights of prior layer are allowed to adapt and flip signs. However,
the disadvantage is more memory are required during training.

In low-rank decompositions and teacher-student network, weights are still in
floating-point precision. For low-rank decomposition, the implementation
requires decomposition operation, which is computationally expensive, and
factorization requires extensive model retraining to achieve convergence when
compared to the original model (Cheng et al., 2018). Similarly, due to
iterative nature of our proposed training regime, training time is also
lengthened.

##  5 Conclusions and Further Work

In this work, we proposed a simple iterative training that gradually trains a
partial binary weight network to a full binary weight network layer-by-layer.
We showed empirically that this regime results in higher accuracy than
starting training from a fully binarized weight network. The order of layer
binarization matters. We show empircally that, for larger and deeper neural
networks, the forward order achieves better accuracies than other binarization
orders. We proposed a sensitivity pre-training for selection of binarization
order. For 784-784-10, Vgg-5, Vgg-9 and ResNet-21, this guided order achieve
better accuracies than the forward order.

Iterative training has a cost, which is lengthened training. This trade-off
may be acceptable in many applications where pre-trained models are deployed,
because efficiency in only forward propagation is needed. A binary weight
neural network dramatically reduces computation complexity, memory footprint
and, thus, increases energy efficiency. For future work, we would like to
understand analytically why layer quantization works and the optimal
quantization order.

## References

  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep Residual Learning for Image Recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778. IEEE, 2016.
  * He et al. [2017] Kaiming He, Georgia Gkioxari, Piotr Dollar, and Ross Girshick.  Mask R-CNN.  In _2017 IEEE International Conference on Computer Vision (ICCV)_ , pages 2980–2988. IEEE, oct 2017.  ISBN 978-1-5386-1032-9.  doi:10.1109/ICCV.2017.322.
  * Chen et al. [2018] Liang-chieh Chen, Yukun Zhu, George Papandreou, Florian Schroff, and Hartwig Adam.  Encoder-Decoder with Atrous Separable Convolution for Semantic Image Segmentation.  In _European Conference on Computer Vision_ , pages 833–851, feb 2018\.
  * Zagoruyko and Komodakis [2016] Sergey Zagoruyko and Nikos Komodakis.  Wide residual networks.  In Edwin R. Hancock Richard C. Wilson and William A. P. Smith, editors, _Proceedings of the British Machine Vision Conference (BMVC)_ , pages 87.1–87.12. BMVA Press, September 2016.  ISBN 1-901725-59-6.  doi:10.5244/C.30.87.
  * Cheng et al. [2018] Yu Cheng, Duo Wang, Pan Zhou, and Tao Zhang.  A Survey of Model Compression and Acceleration for Deep Neural Networks.  _IEEE Signal Processing Magazine_ , 35(1):126–136, 2018.
  * Courbariaux et al. [2015] Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David.  BinaryConnect: Training Deep Neural Networks with binary weights during propagations.  In _Neural Information Processing Systems_ , nov 2015.
  * Horowitz [2014] Mark Horowitz.  1.1 Computing’s energy problem (and what we can do about it).  In _2014 IEEE International Solid-State Circuits Conference Digest of Technical Papers (ISSCC)_ , pages 10–14. IEEE, feb 2014.  ISBN 978-1-4799-0920-9.  doi:10.1109/ISSCC.2014.6757323.
  * Li et al. [2016] Fengfu Li, Bo Zhang, and Bin Liu.  Ternary Weight Networks.  In _NIPS 2016, 1st International Workshop on Efficient Methods for Deep Neural Networks_ , 2016.
  * Hubara et al. [2016] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Binarized neural networks.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 29. Curran Associates, Inc., 2016.
  * Zhou et al. [2016] Shuchang Zhou, Yuxin Wu, Zekun Ni, Xinyu Zhou, He Wen, and Yuheng Zou.  DoReFa-Net: Training Low Bitwidth Convolutional Neural Networks with Low Bitwidth Gradients, Jun 2016.
  * Wu et al. [2018] Shuang Wu, Guoqi Li, Feng Chen, and Luping Shi.  Training and Inference with Integers in Deep Neural Networks.  In _International Conference on Learning Representations_ , pages 1–14, feb 2018.
  * Lecun et al. [1998] Yann Lecun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  volume 86, pages 2278–2324, 1998.
  * Krizhevsky [2009] Alex Krizhevsky.  Learning Multiple Layers of Features from Tiny Images.  Technical report, 2009.
  * Russakovsky et al. [2015] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei.  ImageNet Large Scale Visual Recognition Challenge.  In _International Journal of Computer Vision (IJCV)_ , volume 115, pages 211–252, 2015.  doi:10.1007/s11263-015-0816-y.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift.  In _International Conference on Machine Learning_ , volume 37, pages 448–456, 2015.
  * He et al. [2015] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification.  In _2015 IEEE International Conference on Computer Vision (ICCV)_ , pages 1026–1034. IEEE, 2015.
  * Simonyan and Zisserman [2015] Karen Simonyan and Andrew Zisserman.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  In _International Conference on Learning Representations_ , 2015.
  * Hu et al. [2018] Qinghao Hu, Peisong Wang, and Jian Cheng.  From Hashing to CNNs: Training Binary Weight Networks via Hashing.  In _Thirty-Second AAAI Conference on Artificial Intelligence (AAAI-18)_ , pages 3247–3254, 2018.
  * Zhou et al. [2017] Aojun Zhou, Anbang Yao, Yiwen Guo, Lin Xu, and Yurong Chen.  Incremental Network Quantization: Towards Lossless CNNs with Low-Precision Weights.  In _International Conference on Learning Representations_ , pages 1–14, feb 2017.


# Closed-Loop Model Identification and MPC-based Navigation of Quadcopters: A
Case Study of Parrot Bebop 2

Mohsen Amiri and Mehdi Hosseinzadeh  School of Mechanical and Materials
Engineering, Washington State University, Pullman, WA 99164, USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The growing potential of quadcopters in various domains, such as aerial
photography, search and rescue, and infrastructure inspection, underscores the
need for real-time control under strict safety and operational constraints.
This challenge is compounded by the inherent nonlinear dynamics of quadcopters
and the on-board computational limitations they face. This paper aims at
addressing these challenges. First, this paper presents a comprehensive
procedure for deriving a linear yet efficient model to describe the dynamics
of quadrotors, thereby reducing complexity without compromising efficiency.
Then, this paper develops a steady-state-aware Model Predictive Control (MPC)
to effectively navigate quadcopters, while guaranteeing constraint
satisfaction at all times. The main advantage of the steady-state-aware MPC is
its low computational complexity, which makes it an appropriate choice for
systems with limited computing capacity, like quadcopters. This paper
considers the Parrot Bebop 2 as the running example, and experimentally
validates and evaluates the proposed algorithms.

###### keywords:

Closed-Loop Model Identification, Quadcopter, Model Predictive Control, Parrot
Bebop 2.

##  1 Introduction

Motivation: Quadcopters have demonstrated proficiency in hovering and have
found diverse applications across fields such as aerial photography,
surveillance, precision agriculture, and package delivery (Gandhi and Ghosal,
2018; Hermand et al., 2018; Duggal et al., 2016; Kumar et al., 2018).
Nevertheless, their nonlinear dynamics pose significant challenges for
modeling and control, and specifically for addressing safety and operational
constraints (Shauqee et al., 2021; Bouabdallah et al., 2004; Idrissi et al.,
2022). In particular, the inherent nonlinearity of quadcopter dynamics,
arising from intricate aerodynamic interactions and dynamic coupling between
motion axes, complicates the accurate identification of models and the
development of safe and reliable control strategies.

A further significant challenge in navigating quadcopters, in practical, in
terms of precise position tracking, is that they are often designed to keep
weight and energy consumption as low as possible (Ichnowski et al., 2019;
Hosseinzadeh et al., 2022, 2023a); thus, the computing power of quadcopters is
naturally limited, while multiple computationally expensive diagnostics and
control functions must be executed at any time instant. This underlines the
importance of identifying a simple, yet accurate, model to describe the
dynamics of quadcopters, and of developing a control scheme that is capable of
guaranteeing constraint satisfaction despite limitations of computing power.

Prior Work: Several previous studies have attempted to model quadcopter
dynamics, considering aerodynamic effects, rotor dynamics, and environmental
disturbances; see e.g., (Wang et al., 2016; Khodja et al., 2017; Zhao and Go,
2014; Saif and Eminoglu, 2022; Tang et al., 2017). Identification methods
based on artificial intelligence have been considered in (Kantue and Pedro,
2018; Pairan et al., 2020; Wang et al., 2023). The above-mentioned efforts
provide a nonlinear model to describe the dynamical behavior of quadcopters.
One inevitable issue with the nonlinear models is that controllers designed
based on them tend to be complex and computationally intensive (Wang et al.,
2017), challenging their practical implementation, in particular, when the
available computing power is limited. To overcome these limitations,
specifically for navigation purposes, a shift towards simpler nonlinear models
(e.g., (Pinto et al., 2020)) and linear models (e.g., (Hermand et al., 2018;
Santos et al., 2017a, b; Santana et al., 2014; Santos et al., 2019) has been
made in recent years. A linear model offers a more tractable framework for
modeling and control, facilitating the development of precise navigation
algorithms while addressing safety constraints. To the best of our knowledge,
prior work on identifying a linear model to describe the dynamical behavior of
quadcopters lacks a detailed procedure and/or fails to effectively account for
output constraints; in particular, most of the existing articles perform an
open-loop identification that is not appropriate for indoor identification.

Once a dynamical model is identified, one possible approach (Huang and Tamayo,
2000; Huang et al., 2003) to evaluate its reliability is to investigate the
performance of a Model Predictive Control (MPC) scheme developed based on that
model. Ensuring constraint satisfaction is crucial for safe and efficient
operation of quadcopters (Dubay and Pan, 2018), and MPC has shown promise in
this domain by enabling online decision-making based on predictive models
(Camacho and Bordons, 2007). However, the computational challenges associated
with MPC implementation pose significant barriers. Several approaches have
been proposed to mitigate this challenge, including explicit MPC (Alessio and
Bemporad, 2009; Kvasnica et al., 2019) and self- and event-triggered MPC (Sun
et al., 2019; Wang et al., 2020; Feller and Ebenbauer, 2018). Anytime MPC
(Feller and Ebenbauer, 2018) ensures stability with minimal iterations,
although it cannot guarantee constraint satisfaction at all times. Converting
MPC problem into the evolution of a continuous-time system has been introduced
in (Nicotra et al., 2018; Hosseinzadeh et al., 2023b), while its discrete-time
implementation is still challenging.

Shortening the prediction horizon, as the most intuitive approach to reduce
the computational cost of MPC, has been considered in (Sawma et al., 2018;
Pannek and Worthmann, 2011). To address feasibility issues with small
prediction horizons, (Ferramosca et al., 2009; Limon et al., 2008; Ferramosca
et al., 2008) introduced modifications to the conventional MPC; note that to
ensure output tracking, such modifications need prior knowledge about the
steady state and steady input corresponding to any given reference, which is
computationally demanding (Yousfi and Tournier, 1991; Muske and Rawlings,
1993; Shead et al., 2010) and practically unrealistic, in particular, when the
reference is unknown a priori.

Despite the above-mentioned efforts, achieving computational efficiency, while
addressing performance objectives and enforcing constraint satisfaction
remains a key research focus. Recently, the authors have proposed a steady-
state-aware MPC framework in (Amiri and Hosseinzadeh, 2024) which is
appropriate for systems with limited computing power. The main feature of the
steady-state-aware MPC is that, given a desired steady state and input, it
guarantees that the state and input of the system converge to the desired
steady values, while ensuring output tracking and constraint satisfaction at
all times, without requiring any prior knowledge about the reference.
Simulation studies demonstrated the capability of the steady-state-aware MPC
in guaranteeing output tracking and steady-state convergence, as well as
enforcing constraint satisfaction at all times, all under the limitations on
available computing power. However, experimental validations are still
required to establish the steady-state-aware MPC as a standard approach for
safe and reliable control of safety-critical systems with limited
computational resources.

Goal: In a nutshell, the main goal of this paper is to systematically address
the intricacies of quadcopter dynamics and computational limitations. First,
we advocate for a physically informed modeling approach, and identify a simple
linear yet efficient model for a quadcopter to enable efficient onboard
computation and safe control. Then, we develop a steady-state-aware MPC to
safely navigate the quadcopter, while managing constraints effectively despite
limited available computational resources. We consider a Parrot Bebop 2 (see
Fig. 1) as a running example throughout the paper, and experimentally evaluate
and validate the proposed methodology.

Organization: The rest of this paper is organized as follows. The presentation
of the model structure followed by identification procedure is discussed in
Section 2. Section 3 develops a steady-state-aware MPC and experimentally
assesses its effectiveness on Parrot Bebop2 drone. Finally, Section 4
concludes this paper.

Notation: We denote the set of real numbers by $\mathbb{R}$, the set of
positive real numbers by $\mathbb{R}_{>0}$, and the set of non-negative real
numbers by $\mathbb{R}_{\geq 0}$. Similarly, we use $\mathbb{Z}$ to denote the
set of integer numbers, $\mathbb{Z}_{>0}$ to denote the set of positive
integer numbers, and $\mathbb{Z}_{\geq 0}$ to denote the set of non-negative
integer numbers. We denote the transpose of matrix $A$ by $A^{\top}$. Also,
$A\succ 0$ ($A\succeq 0$) indicates that $A$ is positive definite (positive
semi-definite). Given $x\in\mathbb{R}^{n}$ and $Q\succeq 0$, we have
$\left\|x\right\|_{Q}^{2}=x^{\top}Qx$. We use $t$ to denote continuous time
and $k$ to denote sampling instants. The notation
$\text{diag}\\{a_{1},\cdots,a_{n}\\}$ indicates a $n\times n$ diagonal matrix
with elements $a_{1},\cdots,a_{n}$ on the main diagonal. Finally, $\mathbf{0}$
denotes the zero matrix with appropriate dimension.

##  2 Closed-Loop Model Identification

Figure 1: Parrot Bebop 2 and considered reference and body frames; note that
to make the measurements more relevant, the considered body frame is different
than the one in the “Parrot Drone Support from MATLAB” package (MATLAB, 2019).

###  2.1 Model Structure

Suppose that the yaw angle of the quadcopter is very small; see Subsection 2.2
for details on how one can regulate the yaw angle to zero. Taking inspiration
from (Pinto et al., 2020; Santos et al., 2019), we propose the following
structure to describe the position dynamics of a quadcopter:

$\displaystyle\ddot{p}_{x}+\alpha_{x}\dot{p}_{x}=$
$\displaystyle\beta_{x}u_{x},$ (1a)
$\displaystyle\ddot{p}_{y}+\alpha_{y}\dot{p}_{y}=$
$\displaystyle\beta_{y}u_{y},$ (1b)
$\displaystyle\ddot{p}_{z}+\alpha_{z}\dot{p}_{z}=$
$\displaystyle\beta_{z}u_{z},$ (1c)

where $p_{x},p_{y},p_{z}\in\mathbb{R}$ are X, Y, and Z positions of the
quadcopter in the global Cartesian coordinate system,
$\alpha_{x},\alpha_{y},\alpha_{z},\beta_{x},\beta_{y},\beta_{z}\in\mathbb{R}$
are system parameters, and $u_{x},u_{y},u_{z}\in\mathbb{R}$ are control inputs
on X, Y, and Z directions.

We use the “Parrot Drone Support from MATLAB” (MATLAB, 2019) to send control
commands to the Parrot Bebop 2 drone. In this framework, $u_{x}$ is the pitch
angle (in [rad]), $u_{y}$ is the roll angle (in [rad]), and $u_{z}$ is the
vertical velocity (in [m/s]); see Fig. 1.

The dynamical model (1) can be expressed via the following state-space
equations:

$\displaystyle\dot{x}=$ $\displaystyle\begin{bmatrix}0&1&0&0&0&0\\\
0&-\alpha_{x}&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&-\alpha_{y}&0&0\\\
0&0&0&0&0&1\\\ 0&0&0&0&0&-\alpha_{z}\end{bmatrix}x+\begin{bmatrix}0&0&0\\\
\beta_{x}&0&0\\\ 0&0&0\\\ 0&\beta_{y}&0\\\ 0&0&0\\\
0&0&\beta_{z}\end{bmatrix}u,$ (2a) $\displaystyle y=$
$\displaystyle\begin{bmatrix}1&0&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&0&1&0\end{bmatrix}x,$ (2b)

where
$x=\begin{bmatrix}p_{x}&\dot{p}_{x}&p_{y}&\dot{p}_{y}&p_{z}&\dot{p}_{z}\end{bmatrix}^{\top}$
and $u=\begin{bmatrix}u_{x}&u_{y}&u_{z}\end{bmatrix}^{\top}$.

Note that model (2) is consistent with practical observations. First, dynamics
along X, Y, and Z directions are decoupled, as one can control the position of
the quadcopter along X, Y, and Z directions separately. Second, when
$\dot{p}_{x}\neq 0$ ($\dot{p}_{y}\neq 0$ and $\dot{p}_{z}\neq 0$) and
$u_{x}=0$ ($u_{y}=0$, $u_{z}=0$), the quadcopter continues moving in x
direction (y and z directions) while its velocity along that direction
decreases to zero; from this observation, we anticipate that $\alpha_{x}$,
$\alpha_{y}$, and $\alpha_{z}$ are positive (i.e.,
$\alpha_{x},\alpha_{y},\alpha_{z}\in\mathbb{R}_{>0}$). Third, when pitch and
roll angles are small, the quadcopter behaves as a second-order system.

###  2.2 Identification Procedure

According to model (1), it is possible to obtain the system parameters (i.e.,
$\alpha_{i},\beta_{i},~{}i\in\\{x,y,z\\}$) by applying a control input signal
and measuring its response; this method is called the open-loop identification
and is followed widely in prior work, e.g., (Pinto et al., 2020; Santana et
al., 2014; Santos et al., 2019). Although the quadcopters are stable, open-
loop identification is not appropriate for indoor identification, as the
spatial constraints and limits cannot be enforced. A different approach, which
is called the closed-loop identification (Forssell and Ljung, 1999, 2000), is
to identify the system parameters when a controller is utilized to control the
position and angles of the quadcopter. In addition to enforcing output
constraints, closed-loop identification allows us (den Hof, 1998) to limit the
operating points of the quadcopter to a region in which the quadcopter
presents a linear behavior (i.e., keeping the pitch and roll angles small).

Note that since the measurement noise in motion capture systems used in indoor
applications is very small, the control input due to the measurement noise
will be small as well. Thus, one can use the closed-loop identification
approach without being worried about the correlation between the measurement
noise and the control input, which is discussed in (Forssell and Ljung, 1999).

In what follows, we provide a systematic procedure for closed-loop
identification of quadcopters.

  * •

Network Establishment: The experimental setup should consist of a motion
capture system that measures the position of the quadcopter, a computing unit
that analyzes the camera data and computes the control commands, and a
communication channel that sends the control commands to the quadcopter.

Our experimental setup is shown in Fig. 2. We use OptiTrack system with ten
Prime${}^{\text{x}}$ 13 cameras operating at a frequency of 120 Hz; these
cameras provide a 3D accuracy of $\pm 0.02$ millimeters which is acceptable
for identification purposes. The computing unit is a $13^{\text{th}}$ Gen
$\text{Intel}^{\text{\textregistered}}$ $\text{Core}^{\text{\texttrademark}}$
i9-13900K processor with 64GB RAM, on which the software Motive is installed
to analyze and interpret the camera data. We use the “Parrot Drone Support
from MATLAB” package (MATLAB, 2019), and send the control commands to the
Parrot Bebop 2 via WiFi and by using the command move($\cdot$). It should be
mentioned that the communication between Motive and MATLAB is done through
User Datagram Protocol (UDP) communication using the NatNet service.

Figure 2: Overview of the experimental setup utilized to evaluate and validate
the proposed methods.

  * •

Nominal Controller: One can use a Proportional-Derivative (PD) control law to
control the position of the quadcopter and its yaw angle, as follows:

$\displaystyle u_{x}=$ $\displaystyle-
k_{p}^{x}\left(p_{x}-p_{x}^{d}\right)-k_{d}^{x}\left(\dot{p}_{x}-\dot{p}_{x}^{d}\right),$
(3a) $\displaystyle u_{y}=$ $\displaystyle-
k_{p}^{y}\left(p_{y}-p_{y}^{d}\right)-k_{d}^{y}\left(\dot{p}_{y}-\dot{p}_{y}^{d}\right),$
(3b) $\displaystyle u_{z}=$ $\displaystyle-
k_{p}^{z}\left(p_{z}-p_{z}^{d}\right)-k_{d}^{z}\left(\dot{p}_{z}-\dot{p}_{z}^{d}\right),$
(3c) $\displaystyle u_{\psi}=$
$\displaystyle-k^{\psi}_{p}\left(\psi-\psi^{d}\right)-k^{\psi}_{d}\left(\dot{\psi}-\dot{\psi}^{d}\right),$
(3d)

where $\psi$ is the yaw angle, $p_{x}^{d},p_{y}^{d},p_{z}^{d}\in\mathbb{R}$
are the desired X, Y, and Z trajectories, $\psi^{d}=0$ is the desired yaw
angle, $k^{x}_{p},k_{p}^{y},k^{z}_{p},k^{\psi}_{p}\in\mathbb{R}_{>0}$ are
proportional gains, and
$k^{x}_{d},k_{d}^{y},k^{z}_{d},k^{\psi}_{d}\in\mathbb{R}_{>0}$ are derivative
gains.

Table 1: Parameters of the PD controllers used to control the Parrot Bebop 2 drone. Parameter | Identified Value
---|---
$k_{p}^{x}$ | 0.050
$k_{p}^{y}$ | 0.050
$k_{p}^{z}$ | 1.700
$k_{p}^{\psi}$ | 1.000
$k_{d}^{x}$ | 0.065
$k_{d}^{y}$ | 0.065
$k_{d}^{z}$ | 0.200
$k_{d}^{\psi}$ | 0.800

In our experiments, we select the parameters of the PD controllers as in Table
1. Note that this selection ensures that the yaw angle $\psi$ remains around
zero while Parrot Bebop 2 moves toward the desired position
$(P_{x}^{d},P_{y}^{d},P_{z}^{d})$. Also, note that since the OptiTrack system
does not provide the velocities along X, Y, and Z directions, to implement the
PD controllers (3), we use the Newton’s difference quotient to compute the
velocities along all directions.

  * •

Trajectory Determination For Closed-Loop Identification: The most intuitive
way to identify a system is to apply a combination of sinusoidal inputs at
different frequencies and observe the output of the system. It is well-known
that if the input to a linear system is sinusoidal, the output will also be
sinusoidal, and at exactly the same frequency as the input. This property,
which is called sinusoidal fidelity, provides insights on how to determine the
desired output trajectory for closed-loop identification and on how to
identify the operating region in which the system presents a linear behavior.

For quadcopters, one can follow the following steps for each direction
separately to determine the desired output trajectory that is appropriate for
closed-loop identification: i) consider a sinusoidal trajectory at small
frequencies and apply the PD controllers given in (3); ii) compute Fourier
transform of the control input and of the position of the quadcopter; iii)
gradually increase the frequency of the sinusoidal trajectory until the
corresponding component disappears in the Fourier transform of the
quadcopter’s position (this step would determine the bandwidth of the system);
and iv) gradually increase the magnitude of the sinusoidal trajectory (while
satisfying the spatial limits) until frequency components other than that of
the considered trajectory appear in the Fourier transform of the control input
(this step would determine the saturation levels for the control inputs so as
to ensure a linear behavior).

Following the above-mentioned procedure, we select the following desired
sinusoidal trajectory to identify the dynamics of the Parrot Bebop 2:

$\displaystyle\text{Desired Trajectory}=\frac{1}{10}\Big{(}3\sin(0.2\pi t)$
$\displaystyle+0.6\sin(0.4\pi t)+0.1\sin(0.7\pi t)+0.1\sin(\pi t)\Big{)},$ (4)

which has frequency components at 0.1, 0.2, 0.35, and 0.5 Hz. We consider the
desired trajectory given in (• ‣ 2.2) separately for each direction, while the
desired trajectory for the other two directions is set to zero. Fig. 3
presents the Fourier transform of control input and of the position of the
Parrot Bebop 2 along each direction. As seen in this figure, the selected
desired trajectory ensures that the Parrot Bebop 2 behaves as a linear system.
According to the Fourier analysis, we set the saturation levels for the
control inputs $u_{x}$ and $u_{y}$ to 0.06 [rad], and for the control input
$u_{z}$ to 0.6 [m/s] (i.e., maximum pitch and roll angles are $\pm$4 [deg],
and maximum vertical velocity is $\pm$0.6 [m/s]). Note that although these
numbers are smaller than the manufacturer’s recommendation, they ensure a
linear behavior for the Parrot Bebop 2 drone. Furthermore, we conclude from
the Fourier analysis that -3B-bandwidth of the Parrot Bebop 2 in all
directions is lower than 0.5 Hz.

Figure 3: Fourier transform of positions and corresponding control inputs in
X, Y, and Z direction.

  * •

Data Collection: According to model (1), to identify the system parameters
$\alpha_{i},\beta_{i},~{}i\in\\{x,y,z\\}$, one needs to observe velocity and
acceleration in each direction for a given control input. Since motion capture
systems usually provide only position data, one would need to use numerical
derivation to compute velocity and acceleration data. Note that although
measurement noise in motion capture systems is small, numerical derivation can
amplify the impact of the noise on velocities and accelerations, and
consequently can compromise the accuracy of the identified parameters. Hence,
it becomes imperative to properly filter data.

In our experiments, we use the PD control law given in (3) and the desired
identification trajectory given in (• ‣ 2.2) to collect training data. Fig. 4
shows the obtained data in all directions, where the Newton’s difference
quotient is used to obtain the velocity and acceleration data. Note that in
the identification process, we use the unfiltered position and velocity data.
Whereas, to obtain the acceleration data, we first use a low-pass filter with
a cutoff frequency of 5 Hz to filter the velocity data, then employ the
Newton’s difference quotient to generate acceleration data, and finally use a
low-pass filter with a cutoff frequency of 10 Hz to smoothen the acceleration
data.

Figure 4: Obtained data in X, Y, and Z directions and control input along each
direction

.

  * •

Parameter Identification: Once position, velocity, and acceleration data is
collected, model (1) suggests that one can use the linear regression
techniques (see, e.g., (Kleinbaum et al., 2013)) to identify the system
parameters $\alpha_{i},\beta_{i},~{}i\in\\{x,y,z\\}$.

Table 2: Identified parameters for the Parrot Bebop 2 drone. Parameter | Identified Value
---|---
$\alpha_{x}$ | 0.0527
$\alpha_{y}$ | 0.0187
$\alpha_{z}$ | 1.7873
$\beta_{x}$ | -5.4779
$\beta_{y}$ | -7.0608
$\beta_{z}$ | -1.7382

We use the least square method to identify the parameters of the Parrot Bebop
2 drone. The obtained values are reported in Table 2. Note that the obtained
parameters satisfy our practical observations. For instance, let $p_{x}(0)=0$,
$u_{x}(t)=0$, and $\dot{p}_{x}(0)=\sigma$ for some $\sigma\in\mathbb{R}_{>0}$;
that is, the Parrot Bebop 2 has an initial velocity along X direction, while
the roll angle is zero. According to model (1), $p_{x}(t)$ is given by:

$\displaystyle
p_{x}(t)=\frac{\sigma}{\alpha_{x}}\left(1-e^{-\alpha_{x}t}\right),~{}t\geq 0,$
(5)

which implies that the Parrot Bebop 2 moves toward the point
$\frac{\sigma}{\alpha_{x}}$ in X direction, while its velocity
$\dot{p}_{x}(t)$ converges to zero (note that $\alpha_{x}\in\mathbb{R}_{>0}$);
this behavior is consistent with our practical observations. Also, setting
$p_{x}(0)=\dot{p}_{x}(0)=0$ and $u_{x}(t)=\sigma$ for some
$\sigma\in\mathbb{R}_{>0}$ yields:

$\displaystyle
p_{x}(t)=\frac{\sigma}{\alpha_{x}}t+\frac{\sigma}{\alpha_{x}^{2}}\left(e^{-\alpha_{x}t}-1\right),~{}t\geq
0,$ (6)

which indicates that when the roll angle is kept at $\sigma$ radian, the
Parrot Bebop 2 will stay in motion (with a constant steady-state velocity) in
X direction; this behavior is also consistent with our practical observations.

The obtained parameters satisfy our frequency-domain observations as well.
Fig. 5 compares the magnitude of the obtained transfer function at frequencies
0.1, 0.2, 0.35, and 0.5 Hz with the one computed from the Fourier transform of
practical data shown in Fig. 3. As seen in Fig. 5, the identified model
approximately fits the data in all X, Y, and Z directions.

##  3 Experimental Validation: MPC-Based Navigation

Suggested by prior work, e.g., (Huang and Tamayo, 2000; Huang et al., 2003),
this section aims at experimentally validating the identified model by
investigating the performance of a MPC control scheme developed based upon the
identified model. For this purpose, this section develops a steady-state-aware
MPC (Amiri and Hosseinzadeh, 2024) for system (2) with parameters given in
Subsection 2.2. Our motivation to use the steady-state-aware MPC is that it is
capable of addressing output tracking, steady-state convergence requirements,
and constraint satisfaction, despite the limitations of computing power.

Figure 5: Bode plot of the identified model for each direction (blue lines)
and the magnitudes computed based on Fourier transforms (red asterisks).

###  3.1 System Discretization

The first step to design a steady-state-aware MPC is to discretize system (2)
with an appropriate sampling period, to obtain the following discrete-time
linear system:

$\displaystyle x(k+1)=$ $\displaystyle A_{d}x(k)+B_{d}u(k),$ (7a)
$\displaystyle y(k)=$ $\displaystyle C_{d}x(k),$ (7b)

where the interval $[k,k+1)$ is equal to the considered sampling period.

For the Parrot Bebop 2 drone, we use a sampling period of 0.2 seconds. Note
that the bandwidth of the system is 0.5 Hz (see Subsection 2.2), and thus, the
considered sampling frequency satisfies the Nyquist criterion (Ogata, 1995).
The obtained $A_{d}$, $B_{d}$, and $C_{d}$ matrices are:

$\displaystyle A_{d}=\begin{bmatrix}1&0.19895&0&0&0&0\\\ 0&0.98952&0&0&0&0\\\
0&0&1.000&0.19963&0&0\\\ 0&0&0&0.99627&0&0\\\ 0&0&0&0&1.000&0.16816\\\
0&0&0&0&0&0.69946\end{bmatrix},$ (8a) $\displaystyle
B_{d}=\begin{bmatrix}-0.10917348&0&0\\\ -1.08982035&0&0\\\ 0&-0.141040918&0\\\
0&-1.409531141&0\\\ 0&0&-0.030967224\\\ 0&0&-0.292295416\end{bmatrix},$ (8b)
$\displaystyle C_{d}=\begin{bmatrix}1&0&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&0&1&0\end{bmatrix}.$ (8c)

###  3.2 Control Structure

Given $r=[r_{1}~{}r_{2}~{}r_{3}]^{\top}\in\mathbb{R}^{3}$ as the reference
position of the quadcopter, since $(A_{d},B_{d})$ is controllable111Since
$(A_{d},B_{d})$ is controllable, for any given $r$, there exists (Muske and
Rawlings, 1993) a steady state $\mathbf{x}_{f}$ and steady input
$\mathbf{u}_{f}$ such that
$\mathbf{x}_{f}=A_{d}\mathbf{x}_{f}+B_{d}\mathbf{u}_{f}$ and
$r=C_{d}\mathbf{x}_{f}$., it is always possible (Amiri and Hosseinzadeh, 2024;
Limon et al., 2008) to characterize the reference signal $r$, corresponding
steady state $\mathbf{x}_{f}$, and steady input $\mathbf{u}_{f}$ as
$r=L\theta$, $\mathbf{x}_{f}=M\theta$ and $\mathbf{u}_{f}=W\theta$, where
$L\in\mathbb{R}^{3\times 3}$, $M\in\mathbb{R}^{6\times 3}$,
$W\in\mathbb{R}^{3\times 3}$, and $\theta\in\mathbb{R}^{3}$ is the
characterizing vector.

Given the prediction horizon $N\in\mathbb{Z}_{>0}$, the steady-state-aware MPC
computes the optimal characterizing vector $\theta^{\ast}(k)$ and the optimal
control sequence
$\mathbf{u}^{\ast}(k):=\left[\begin{matrix}\left(u^{\ast}(0|k\right)^{\top}&\cdots&\left(u^{\ast}(N-1|k)\right)^{\top}\end{matrix}\right]^{\top}\in\mathbb{R}^{3N}$
at any time instant $k$ by solving the following optimization problem:

$\displaystyle\theta^{\ast}(k),\mathbf{u}^{\ast}(k)=\arg\,\min_{\mathbf{u},\theta}\Big{(}\sum_{s=0}^{N-1}\left\|\hat{x}(s|k)-M\theta\right\|_{Q_{x}}^{2}$
$\displaystyle+\sum_{s=0}^{N-1}\left\|u(s|k)-W\theta\right\|_{Q_{u}}^{2}+\left\|\hat{x}(N|k)-M\theta\right\|_{Q_{N}}^{2}$
$\displaystyle+\left\|r-L\theta\right\|_{Q_{r}}^{2}+\left\|M\theta-
x_{des}\right\|_{Q_{fx}}^{2}+\left\|W\theta-
u_{des}\right\|_{Q_{fu}}^{2}\Big{)},$ (9a) subject to the following
constraints222Note that since the purpose of this section is to validate the
identified model, we do not consider any constraints on the states, like
spatial constraints. It is noteworthy that the steady-state-aware MPC is able
to handle both state and input constraints.:
$\displaystyle\hat{x}(s+1|k)=A_{d}\hat{x}(s|k)+B_{d}u(s|k),~{}\hat{x}(0|k)=x(k),$
(9b) $\displaystyle u(s|k)\in\mathcal{U},~{}s\in\\{0,\cdots,N-1\\},$ (9c)
$\displaystyle(\hat{x}(N|k),\theta)\in\Omega,$ (9d)

where $Q_{x}=Q_{x}^{\top}\succeq 0$ ($Q_{x}\in\mathbb{R}^{6\times 6}$),
$Q_{u}=Q_{u}^{\top}\succ 0$ ($Q_{u}\in\mathbb{R}^{3\times 3}$),
$Q_{fx}=Q_{fx}^{\top}\succeq 0$ ($Q_{fx}\in\mathbb{R}^{6\times 6}$),
$Q_{fu}=Q_{fu}^{\top}\succeq 0$ ($Q_{fu}\in\mathbb{R}^{3\times 3}$),
$x_{des}\in\mathbb{R}^{6}$ and $u_{des}\in\mathbb{R}^{3}$ are desired steady
state and steady input determined by the designer, $Q_{N}\succeq 0$
($Q_{N}\in\mathbb{R}^{6\times 6}$), $Q_{r}=Q_{r}^{\top}\succ 0$
($Q_{r}\in\mathbb{R}^{3\times 3}$), $\mathcal{U}$ is the constraint set on the
control inputs (i.e., $\mathcal{U}=\\{u||u_{x}|\leq 0.06,|u_{y}|\leq
0.06,|u_{z}|\leq 0.6\\}$), and $\Omega$ is the terminal constraint set (see
subsection 3.4).

For quadcopters with dynamical model described in (2), the desired steady
state and steady input are
$x_{des}=[r_{1}~{}0~{}r_{2}~{}0~{}r_{3}~{}0]^{\top}$ and
$u_{des}=[0~{}0~{}0]^{\top}$, respectively. Also, matrices $M$, $L$, and $W$
are computed as:

$\displaystyle M=\begin{bmatrix}1&0&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&0&0&1&0\end{bmatrix}^{\top},~{}L=\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{bmatrix},~{}W=\begin{bmatrix}0&0&0\\\ 0&0&0\\\ 0&0&0\end{bmatrix}.$
(10)

Note that since $W$ is the zero matrix and $u_{des}$ is set to zero, for
quadcopters modeled as in (2), the last term (i.e.,
$\left\|\mathbf{u}_{f}-u_{des}\right\|_{Q_{fu}}^{2}$) in the cost function
(9a) is zero.

###  3.3 Determination of $Q_{N}$

Drawing inspiration from conventional MPC such as those discussed in (Nicotra
et al., 2018; Hosseinzadeh et al., 2023b), it is advantageous to set $Q_{N}$
in (9a) to the solution to the algebraic Riccati equation
$Q_{N}=A_{d}^{\top}Q_{N}A_{d}-(A_{d}^{\top}Q_{N}B_{d})(Q_{u}+B_{d}^{\top}Q_{N}B_{d})^{-1}(B_{d}^{\top}Q_{N}A_{d})+Q_{x}$.

###  3.4 Determining the Terminal Constraint Set $\Omega$

Let the terminal control law be $\kappa(x(k),\theta)=K(x(k)-M\theta)$, where
the gain $K$ is
$K=-(Q_{u}+B_{d}^{\top}Q_{N}B_{d})^{-1}(B_{d}^{\top}Q_{N}A_{d})$, where
$Q_{N}$ is as in Subsection 3.3. It is easy to show that $A_{d}+B_{d}K$ is
Schur.

Following the arguments similar to (Amiri and Hosseinzadeh, 2024), the
terminal constraint set $\Omega$ can be computed as follows:

$\displaystyle\Omega=\\{(x,\theta)|\hat{u}(\omega|x,\theta)\in\mathcal{U},\omega=0,\cdots,\omega^{\ast}\\},$
(11)

where $\omega^{\ast}\in\mathbb{Z}_{>0}$ can be obtained by solving a sequence
of mathematical programming problems detailed in (Gilbert and Tan, 1991), and

$\displaystyle\hat{u}(\omega$
$\displaystyle|x,\theta)=K(A_{d}+B_{d}K)^{\omega}\hat{x}(N|k)$
$\displaystyle-K\sum_{j=1}^{\omega}(A_{d}+B_{d}K)^{j-1}B_{d}KM\theta-
B_{d}KM\theta.$ (12)

Note that the terminal constraint set $\Omega$ given in (11) is convex and
positively invariant, and can be easily computed at any time instant $k$.

###  3.5 Theoretical Properties

This subsection proves the properties of the steady-state-aware MPC designed
to navigate the quadcopter. First, we show that the proposed steady-state-
aware MPC is recursively feasible.

######  Theorem 3.1 (Recursive Feasibility)

Consider system (7), which is subject to the above-mentioned constraints on
the control inputs. Suppose that (9) is feasible at $k=0$. Then, it remains
feasible for all $k>0$.

Suppose that MPC (9) is feasible at $k$, where the optimal characterizing
vector and the optimal control sequence are denoted by $\theta^{\ast}(k)$ and
$\mathbf{u}^{\ast}(k)$, respectively. Also,
$(\hat{x}(N|k),\theta^{\ast}(k))\in\Omega$, where $\hat{x}(N|k)$ is the
terminal state given the optimal control sequence $\mathbf{u}^{\ast}(k)$.

Since the terminal constraint set $\Omega$ is positively invariant, it is
concluded that
$\kappa\big{(}\hat{x}(N|k),\theta^{\ast}(k)\big{)}=K(A_{d}+B_{d}K)\hat{x}(N|k)-KB_{d}KM\theta^{\ast}(k)-B_{d}KM\theta^{\ast}(k)\in\mathcal{U}$,
and $\left(\hat{x}(N+1|k),\theta^{\ast}(k)\right)\in\Omega$. Thus, the
characterizing vector $\theta^{\ast}(k)$ and the control sequence
$\big{[}\left(u^{\ast}(1|k)\right)^{\top}~{}\cdots$
$\left(u^{\ast}(N-1|k)\right)^{\top}~{}\left(\kappa\big{(}\hat{x}(N|k),\theta^{\ast}(k)\big{)}\right)^{\top}\big{]}^{\top}$
construct a feasible solution for the proposed steady-state-aware MPC at time
instant $k+1$. Therefore, feasibility at time instant $k$ implies feasibility
at time instant $k+1$. This implies that the steady-state-aware MPC given in
(9) is recursively feasible, which completes the proof. ∎

Next, we show that the proposed steady-state-aware MPC ensures closed-loop
stability.

######  Theorem 3.2 (Closed-Loop Stability)

Suppose that the steady-state-aware MPC given in (9) is used to control system
(7). Then, $x(k)\rightarrow x_{des}$, $u(k)\rightarrow u_{dex}$, and
$y(k)\rightarrow r$ as $k\rightarrow\infty$, where
$x_{des}=[r_{1}~{}0~{}r_{2}~{}0~{}r_{3}~{}0]^{\top}$ and
$u_{des}=[0~{}0~{}0]^{\top}$.

Let $J\big{(}\theta,\mathbf{u}|x(k)\big{)}$ and
$J\big{(}\theta,\mathbf{u}|x(k+1)\big{)}$ be the cost functions of the steady-
state-aware MPC given in (9) at time instants $k$ and $k+1$, respectively.
Also, let $\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)\big{)}$ and
$\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)\big{)}$ be the optimal
solutions at time instants $k$ and $k+1$, respectively. According to the
optimality of the solution
$\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)\big{)}\in\mathbb{R}^{3}\times\mathbb{R}^{3}$
at time instant $k+1$, we have:

$\displaystyle
J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\leq
J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}.$ (13)

By subtracting $J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}$
from both sides of (3.5), one can get the following inequality:

$\displaystyle
J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}$
$\displaystyle~{}~{}~{}\leq
J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}.$
(14)

According to (9a), and since $u^{\ast}(s+1|k)=u^{\ast}(s|k+1)$ and
$\hat{x}(s+1|k)=\hat{x}(s|k+1),~{}s=0,\cdots,N-2$, we obtain that:

$\displaystyle
J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}$
$\displaystyle\leq\left\|\hat{x}(N|k+1)-M\theta^{\ast}(k)\right\|_{Q_{N}}^{2}-\left\|\hat{x}(N|k)-M\theta^{\ast}(k)\right\|_{Q_{N}}^{2}$
$\displaystyle+\left\|\hat{x}(N-1|k+1)-M\theta^{\ast}(k)\right\|_{Q_{x}}^{2}+\left\|u(N-1|k+1)\right\|_{Q_{u}}^{2}$
$\displaystyle-\left\|\hat{x}(0|k)-M\theta^{\ast}(k)\right\|_{Q_{x}}^{2}-\left\|u^{\ast}(0|k)\right\|_{Q_{u}}^{2}.$
(15)

As shown in (Mayne et al., 2000; Amiri and Hosseinzadeh, 2024), when the
matrix $Q_{N}$ is determined as in Subsection 3.3, we have
$\big{\|}\hat{x}(N|k+1)-M\theta^{\ast}(k)\big{\|}_{Q_{N}}^{2}-\big{\|}\hat{x}(N|k)-M\theta^{\ast}(k)\big{\|}_{Q_{N}}^{2}+\big{\|}\hat{x}(N-1|k+1)-M\theta^{\ast}(k)\big{\|}_{Q_{x}}^{2}+\big{\|}u(N-1|k+1)\big{\|}_{Q_{u}}^{2}\leq
0$. Thus, it follows from (3.5) that:

$\displaystyle
J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}$
$\displaystyle\leq-\big{\|}\hat{x}(0|k)-M\theta^{\ast}(k)\big{\|}_{Q_{x}}^{2}-\big{\|}u^{\ast}(0|k)\big{\|}_{Q_{u}}^{2}\leq
0.$ (16)

According to (3.5), it is obvious that imposing
$J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}\equiv
0$ implies that $u(k)=u_{des}=\mathbf{0}$. In the following, we show that
imposing
$J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}\equiv
0$ implies that $x(k)=x_{des}$ as well.

On the one hand, as shown in (Amiri and Hosseinzadeh, 2024, Appendix II),
imposing
$J\big{(}\theta^{\ast}(k+1),\mathbf{u}^{\ast}(k+1)|x(k+1)\big{)}-J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}=0$
for $k\geq k^{\dagger}$ (for some $k^{\dagger}\in\mathbb{Z}_{\geq 0}$) implies
that the optimal cost function
$J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}$ takes the
following form for $k\geq k^{\dagger}$:

$\displaystyle J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}=$
$\displaystyle\left\|r-L\theta^{\ast}(k)\right\|_{Q_{r}}^{2}$
$\displaystyle+\left\|M\theta^{\ast}(k)-x_{des}\right\|_{Q_{fx}}^{2}.$ (17)

On the other hand, from optimality of the solution
$\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)\big{)}$ for $k\geq k^{\dagger}$,
we have:

$\displaystyle J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}\leq
J\big{(}r,\mathbf{u}^{\ast}(k)|x(k)\big{)},~{}k\geq k^{\dagger},$ (18)

From (3.5) and (18), we obtain that:

$\displaystyle J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}\leq$
$\displaystyle\left\|r-Lr\right\|_{Q_{r}}^{2}+\left\|Mr-
x_{des}\right\|_{Q_{fx}}^{2}.$ (19)

which according to matrices $M$ and $L$ given in (10), we have:

$\displaystyle J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}\leq$
$\displaystyle 0.$ (20)

Since $J\big{(}\theta^{\ast}(k),\mathbf{u}^{\ast}(k)|x(k)\big{)}$ cannot be
negative, it follows from (20) that
$\left\|r-L\theta^{\ast}(k)\right\|_{Q_{r}}^{2}=\left\|M\theta^{\ast}(k)-x_{des}\right\|_{Q_{fx}}^{2}=0$
for $k\geq k^{\dagger}$, or equivalently $L\theta^{\ast}(k)=r$ and
$M\theta^{\ast}(k)=x_{des}$ for $k\geq k^{\dagger}$. This means that the
system is at steady-state for $k\geq k^{\dagger}$, and the steady output and
state are $y(k)=L\theta^{\ast}(k)=r$ and $x(k)=M\theta^{\ast}(k)=x_{des}$,
respectively.

According to the above-mentioned discussion, we conclude that $x(k)\rightarrow
x_{des}$, $u(k)\rightarrow u_{des}$, and $y(k)\rightarrow r$ as
$k\rightarrow\infty$, which completes the proof. ∎

###  3.6 Experimental Setting

To conduct our experimental analysis, we set $N=10$,
$Q_{x}=\text{diag}\\{5,5,5,5,5\\}$, $Q_{u}=\text{diag}\\{35,20,1\\}$,
$Q_{r}=\text{diag}\\{500,500,500\\}$, and
$Q_{fx}=\text{diag}\\{1,1,1,1,1,1\\}$. We implement the steady-state-aware MPC
with one sample delay (Roy et al., 2021; Hosseinzadeh et al., 2022); that is,
the control signal computed based on the measurements at sampling instant $k$
is applied to the Parrot Bebop 2 at sampling instant $k+1$. We use the YALMIP
toolbox (Lofberg, 2004) to perform the computations of the developed steady-
state-aware MPC.

###  3.7 Experimental Results

Constant Piecewise Reference Tracking: We assume that desired positions
$p_{x}^{d}(t)$, $p_{u}^{d}(t)$, and $p_{z}^{d}(t)$ are piecewise constant.
Experimental results are presented in Fig. 6. As seen in this figure, the
developed steady-state-aware MPC effectively steers the Parrot Bebop 2 to the
desired location even in the presence of abrupt changes in the reference
signal. The obtained results can serve to validate the identified parameters
model for Parrot Bebop 2 which is reported in Subsection 2.2.

Figure 6: Set points traveled by the quadrotor Bebop 2 in the navigation
experiment.

Trajectory Tracking: Let $p_{z}^{d}(t)=1.5$, and $p_{x}^{d}(t)$ and
$P_{y}^{d}(t)$ be such that they construct the lemniscate of Bernoulli curve.
Fig. 7 shows the trajectory traveled by the Parrot Bebop 2 in 3D. As one can
see, the task of trajectory tracking with the developed steady-state-aware MPC
is accomplished, implying that the identified model for the Parrot Bebop 2 is
reliable.

Discussion: Experimental results shown in Fig. 6 and 7 affirm that although
the model (2) is a simplified approximation of the dynamics of the Parrot
Bebop 2, it can be safely used to design a model-based control scheme for it.
Also, the computing time of the developed steady-state-aware MPC is
0.0518$\pm$0.0055 seconds, which indicates that the steady-state-aware MPC can
be employed to address control objectives while ensuring constraint
satisfaction even on low-end hardware.

Figure 7: Path followed by the quadrotor when tracking the lemniscate of
Bernoulli trajectory.

##  4 Conclusion

This paper proposed to systematically address the complexities of quadcopters’
dynamics and computational limitations. First, a physically informed modeling
approach was advocated, leading to the identification of a simple yet
efficient linear model for enabling onboard computation efficiency and
ensuring safe navigation. Second, a steady-state-aware MPC framework was
developed to effectively navigate quadcopters, while guaranteeing constraint
satisfaction at all times within the confines of limited computational
resources. Using a Parrot Bebop 2 as a running example, proposed methods and
algorithms were evaluated and validated experimentally.

## References

  * Alessio and Bemporad (2009) Alessio, A. and Bemporad, A. (2009).  A survey on explicit model predictive control.  In L. Magni, D.M. Raimondo, and F. Allgower (eds.), _Nonlinear Model Predictive Control: Towards New Challenging Applications_ , 345–369. Springer Berlin, Heidelberg.
  * Amiri and Hosseinzadeh (2024) Amiri, M. and Hosseinzadeh, M. (2024).  Steady-state-aware model predictive control for tracking in systems with limited computing capacity.  _IEEE Control Syst. Lett._ DOI:10.1109/LCSYS.2024.3370266.
  * Bouabdallah et al. (2004) Bouabdallah, S., Noth, A., and Siegwart, R. (2004).  PID vs LQ control techniques applied to an indoor micro quadrotor.  In _Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems_ , volume 3, 2451–2456. Sendai, Japan.
  * Camacho and Bordons (2007) Camacho, E.F. and Bordons, C. (2007).  _Model predictive controllers_.  Springer.
  * den Hof (1998) den Hof, P.V. (1998).  Closed-loop issues in system identification.  _Annu. Rev. Control._ , 22, 173–186.
  * Dubay and Pan (2018) Dubay, S. and Pan, Y.J. (2018).  Distributed mpc based collision avoidance approach for consensus of multiple quadcopters.  In _Proc. IEEE 14th Int. Conf. Control and Automation_ , 155–160.
  * Duggal et al. (2016) Duggal, V., Sukhwani, M., Bipin, K., Reddy, G.S., and Krishna, K.M. (2016).  Plantation monitoring and yield estimation using autonomous quadcopter for precision agriculture.  In _Proc. IEEE Int. Conf. Robotics and Automation_ , 5121–5127. Stockholm, Sweden.
  * Feller and Ebenbauer (2018) Feller, C. and Ebenbauer, C. (2018).  Sparsity-exploiting anytime algorithms for model predictive control: A relaxed barrier approach.  _IEEE Trans. Control Syst. Technol._ , 28(2), 425–435.
  * Ferramosca et al. (2008) Ferramosca, A., Limon, D., Alvarado, I., Alamo, T., and Camacho, E.F. (2008).  MPC for tracking with optimal closed-loop performance.  In _Proc. 47th IEEE Conf. Decision and Control_. Cancun, Mexico.
  * Ferramosca et al. (2009) Ferramosca, A., Limon, D., Alvarado, I., Alamo, T., and Camacho, E.F. (2009).  MPC for tracking with optimal closed-loop performance.  _Automatica_ , 45(8), 1975–1978.
  * Forssell and Ljung (1999) Forssell, U. and Ljung, L. (1999).  Closed-loop identification revisited.  _Automatica_ , 35(7), 1215–1241.
  * Forssell and Ljung (2000) Forssell, U. and Ljung, L. (2000).  A projection method for closed-loop identification.  _IEEE Trans. Autom. Control_ , 45(11), 2101–2106.
  * Gandhi and Ghosal (2018) Gandhi, D.A. and Ghosal, M. (2018).  Novel low cost quadcopter for surveillance application.  In _Proc. Int. Conf. Inventive Research in Computing Applications_ , 412–414. Coimbatore, India.
  * Gilbert and Tan (1991) Gilbert, E.G. and Tan, K.T. (1991).  Linear systems with state and control constraints: the theory and application of maximal output admissible sets.  _IEEE Trans. Autom. Control_ , 36(9), 1008–1020.
  * Hermand et al. (2018) Hermand, E., Nguyen, T.W., Hosseinzadeh, M., and Garone, E. (2018).  Constrained control of UAVs in geofencing applications.  In _Proc. 26th Mediterranean Conf. Control and Automation_ , 217–222. Zadar, Croatia.
  * Hosseinzadeh et al. (2022) Hosseinzadeh, M., Sinopoli, B., Kolmanovsky, I., and Baruah, S. (2022).  ROTEC: Robust to early termination command governor for systems with limited computing capacity.  _Syst. Control Lett._ , 161.  Art. no. 105142.
  * Hosseinzadeh et al. (2023a) Hosseinzadeh, M., Sinopoli, B., Kolmanovsky, I., and Baruah, S. (2023a).  Robust to early termination model predictive control.  _IEEE Trans. Autom. Control_.  DOI: 10.1109/TAC.2023.3308817.
  * Hosseinzadeh et al. (2023b) Hosseinzadeh, M., Sinopoli, B., Kolmanovsky, I., and Baruah, S. (2023b).  Robust to early termination model predictive control.  _IEEE Trans. Autom. Control_.  DOI: 10.1109/TAC.2023.3308817.
  * Huang et al. (2003) Huang, B., Malhotra, A., and Tamayo, E.C. (2003).  Model predictive control relevant identification and validation.  _Chem. Eng. Sci._ , 58(11), 2389–2401.
  * Huang and Tamayo (2000) Huang, B. and Tamayo, E.C. (2000).  Model validation for industrial model predictive control systems.  _Chem. Eng. Sci._ , 55(12), 2315–2327.
  * Ichnowski et al. (2019) Ichnowski, J., Prins, J., and Alterovitz, R. (2019).  The economic case for cloud-based computation for robot motion planning.  In _Proc. 18th Int. Symp. Robotics Research_ , 59–65. Puerto Varas, Chile.
  * Idrissi et al. (2022) Idrissi, M., Salami, M., and Annaz, F. (2022).  A review of quadrotor unmanned aerial vehicles: applications, architectural design and control algorithms.  _J. Intell. Robot. Syst._ , 104(2), 22.
  * Kantue and Pedro (2018) Kantue, P. and Pedro, J.O. (2018).  Nonlinear identification of an unmanned quadcopter rotor dynamics using RBF neural networks.  In _Proc. 22nd Int. Conf. System Theory, Control and Computing_ , 292–298. IEEE.
  * Khodja et al. (2017) Khodja, M.A., Tadjine, M., Boucherit, M.S., and Benzaoui, M. (2017).  Experimental dynamics identification and control of a quadcopter.  In _Proc. 6th Int. Conf. Systems and Control_ , 498–502. Batna, Algeria.
  * Kleinbaum et al. (2013) Kleinbaum, D.G., Kupper, L.L., Nizam, A., and Rosenberg, E.S. (2013).  _Applied Regression Analysis and Other Multivariable Methods_.  Cengage Learning, 5 edition.
  * Kumar et al. (2018) Kumar, K.V.V.M.S., Sohail, M., Mahesh, P., and Nelakuditi, U.R. (2018).  Crowd monitoring and payload delivery drone using quadcopter based UAV system.  In _Proc. Int. Conf. Smart Systems and Inventive Technology_. Tirunelveli, India.
  * Kvasnica et al. (2019) Kvasnica, M., Bakaráč, P., and Klaučo, M. (2019).  Complexity reduction in explicit MPC: A reachability approach.  _Syst. Control Lett._ , 124, 19–26.
  * Limon et al. (2008) Limon, D., Alvarado, I., Alamo, T., and Camacho, E.F. (2008).  MPC for tracking piecewise constant references for constrained linear systems.  _Automatica_ , 44(9), 2382–2387.
  * Lofberg (2004) Lofberg, J. (2004).  YALMIP: a toolbox for modeling and optimization in MATLAB.  In _Proc. IEEE Int. Conf. Robotics and Automation_ , 284–289. Taipei, Taiwan.
  * MATLAB (2019) MATLAB (2019).  Parrot drone support from matlab.  https://www.mathworks.com/hardware-support/parrot-drone-matlab.html.  [Accessed ].
  * Mayne et al. (2000) Mayne, D.Q., Rawlings, J.B., Rao, C.V., and Scokaert, P.O. (2000).  Constrained model predictive control: Stability and optimality.  _Automatica_ , 36(6), 789–814.
  * Muske and Rawlings (1993) Muske, K.R. and Rawlings, J.B. (1993).  Model predictive control with linear models.  _AIChE J._ , 39(2), 262–287.
  * Nicotra et al. (2018) Nicotra, M.M., Liao-McPherson, D., and Kolmanovsky, I.V. (2018).  Embedding constrained model predictive control in a continuous-time dynamic feedback.  _IEEE Trans. Autom. Control_ , 64(5), 1932–1946.
  * Ogata (1995) Ogata, K. (1995).  _Discrete-Time Control Systems_.  Pearson, 2 edition.
  * Pairan et al. (2020) Pairan, M.F., Shamsudin, S.S., and Zulkafli, M.F. (2020).  Neural network based system identification for quadcopter dynamic modelling: A review.  _J. Adv. Mech. Eng. Appl._ , 1(2), 20–33.
  * Pannek and Worthmann (2011) Pannek, J. and Worthmann, K. (2011).  Reducing the prediction horizon in NMPC: An algorithm based approach.  In _Proc. IEEE Int. Conf. Robotics and Automation_ , 7969–7974. Taipei, Taiwan.
  * Pinto et al. (2020) Pinto, A.O., Marciano, H.N., Bacheti, V.P., Moreira, M.S.M., Brandão, A.S., and Sarcinelli-Filho, M. (2020).  High-level modeling and control of the bebop 2 micro aerial vehicle.  In _Proc. Int. Conf. Unmanned Aircraft Systems_ , 939–947. Athens, Greece.
  * Roy et al. (2021) Roy, D., Hobbs, C., Anderson, J.H., Caccamo, M., and Chakraborty, S. (2021).  Timing debugging for cyber—physical systems.  In _Proc. 2021 Design, Automation and Test in Europe Conference and Exhibition_ , 1893–1898. Grenoble, France.
  * Saif and Eminoglu (2022) Saif, E. and Eminoglu, I. (2022).  Modelling of quad-rotor dynamics and hardware-in-the-loop simulation.  _J. Eng._ , 2022(10), 937–950.
  * Santana et al. (2014) Santana, L.V., Brandao, A.S., Sarcinelli-Filho, M., and Carelli, R. (2014).  A trajectory tracking and 3D positioning controller for the AR.Drone quadrotor.  In _Proc. Int. Conf. Unmanned Aircraft Systems_ , 756–767. Orlando, FL, USA.
  * Santos et al. (2019) Santos, M.C.P., Rosales, C.D., Sarapura, J.A., Sarcinelli-Filho, M., and Carelli, R. (2019).  An adaptive dynamic controller for quadrotor to perform trajectory tracking tasks.  _J. Intell. Robot. Syst._ , 93, 5–16.
  * Santos et al. (2017a) Santos, M.C.P., Rosales, C.D., Sarcinelli-Filho, M., and Carelli, R. (2017a).  A novel null-space-based UAV trajectory tracking controller with collision avoidance.  _IEEE/ASME Trans. Mechatron._ , 22(6), 2543–2553.
  * Santos et al. (2017b) Santos, M.C., Santana, L.V., Brandao, A.S., Sarcinelli-Filho, M., and Carelli, R. (2017b).  Indoor low-cost localization system for controlling aerial robots.  _Control Eng. Pract._ , 61, 93–111.
  * Sawma et al. (2018) Sawma, J., Khatounian, F., Monmasson, E., Ghosn, R., and Idkhajine, L. (2018).  The effect of prediction horizons in MPC for first order linear systems.  In _Proc. IEEE Int. Conf. Industrial Technology_ , 316–321. Lyon, France.
  * Shauqee et al. (2021) Shauqee, M.N., Rajendran, P., and Suhadis, N.M. (2021).  Quadrotor controller design techniques and applications review.  _INCAS Bulletin_ , 13(3), 179–194.
  * Shead et al. (2010) Shead, L.R.E., Muske, K.R., and Rossiter, J.A. (2010).  Conditions for which linear MPC converges to the correct target.  _J. Process Control_ , 20(10), 1243–1251.
  * Sun et al. (2019) Sun, Z., Dai, L., Liu, K., Dimarogonas, D.V., and Xia, Y. (2019).  Robust self-triggered MPC with adaptive prediction horizon for perturbed nonlinear systems.  _IEEE Trans. Autom. Control_ , 64(11), 4780–4787.
  * Tang et al. (2017) Tang, Y.R., Xiao, X., and Li, Y. (2017).  Nonlinear dynamic modeling and hybrid control design with dynamic compensator for a small-scale uav quadrotor.  _Meas._ , 109, 51–64.
  * Wang et al. (2020) Wang, B., Huang, J., Wen, C., Rodriguez, J., Garcia, C., Gooi, H.B., and Zeng, Z. (2020).  Event-triggered model predictive control for power converters.  _IEEE Trans. Ind. Electron._ , 68(1), 715–720.
  * Wang et al. (2023) Wang, M., Zhou, J., Duan, X., Zhao, D., Cai, P., Zhai, J., Liu, X., and Ren, C. (2023).  Nonlinear system identification for quadrotors with neural ordinary differential equations.  In _Proc. IEEE Int. Conf. Unmanned Systems_ , 317–322.
  * Wang et al. (2016) Wang, P., Man, Z., Cao, Z., Zheng, J., and Zhao, Y. (2016).  Dynamics modelling and linear control of quadcopter.  In _Proc. Int. Conf. Advanced Mechatronic Systems_ , 498–503. Melbourne, VIC, Australia.
  * Wang et al. (2017) Wang, Y., Ramirez-Jaime, A., Xu, F., and Puig, V. (2017).  Nonlinear model predictive control with constraint satisfactions for a quadcopter.  _J. Phys. Conf. Ser._ , 783(1), 012025.
  * Yousfi and Tournier (1991) Yousfi, C. and Tournier, R. (1991).  Steady state optimization inside model predictive control.  In _Proc. American Control Conf._ , 1866–1870. Boston, MA, USA.
  * Zhao and Go (2014) Zhao, W. and Go, T.H. (2014).  Quadcopter formation flight control combining MPC and robust feedback linearization.  _J. Frank. Inst._ , 351(3), 1335–1355.