# High-resolution building and road
detection from Sentinel-2

Wojciech Sirko, Emmanuel Asiedu Brempong, Juliana T. C. Marcos, Abigail Annkah
Abel Korme, Mohammed Alewi Hassen, Krishna Sapkota, Tomer Shekel,
Abdoulaye Diack, Sella Nevo, Jason Hickey, John Quinn
Google Research  Address for correspondence<EMAIL_ADDRESS>
###### Abstract

Mapping buildings and roads automatically with remote sensing typically
requires high-resolution imagery, which is expensive to obtain and often
sparsely available. In this work we demonstrate how multiple 10 m resolution
Sentinel-2 images can be used to generate 50 cm resolution building and road
segmentation masks. This is done by training a ‘student’ model with access to
Sentinel-2 images to reproduce the predictions of a ‘teacher’ model which has
access to corresponding high-resolution imagery. While the predictions do not
have all the fine detail of the teacher model, we find that we are able to
retain much of the performance: for building segmentation we achieve 78.3%
mIoU, compared to the high-resolution teacher model accuracy of 85.3% mIoU. We
also describe a related method for counting individual buildings in a
Sentinel-2 patch which achieves $R^{2}=0.91$ against true counts. This work
opens up new possibilities for using freely available Sentinel-2 imagery for a
range of tasks that previously could only be done with high-resolution
satellite imagery.

##  1 Introduction

Buildings and roads are important to map for a range of practical
applications. Models that do this automatically using high-resolution (50 cm
or better) satellite imagery are increasingly effective. However, this type of
imagery is difficult to obtain: there is little control over revisit times,
the cost can be prohibitive for larger analyses, and there may be limited or
no historical imagery for a particular area. This rules out certain analyses,
such as spatially comprehensive surveys of buildings and roads, or systematic
study of changes over time, e.g. for studying urbanisation, economic or
environmental changes.

Figure 1: Example operation of our model, where multiple frames of low-
resolution Sentinel-2 imagery are used to make a single frame of high-
resolution predictions for a variety of output types. A high-resolution image
of the same scene is shown for comparison.

The Sentinel-2 earth observation missions collect imagery globally every 2-5
days and depending on the band at up to 10 m ground resolution. This freely
available data source is commonly used for the adjacent task of land cover
mapping, and recent work on super-resolution with Sentinel-2 imagery has shown
that a higher level of detail can be obtained than the native image resolution
would suggest. This is intuitively possible because of the small variation in
spatial position across successive Sentinel-2 image frames, caused by
atmospheric disturbances, and even across bands in the same image due to
sensor layout. A slightly different $10^{2}$ m2 area is being captured each
time, and hence a sequence of such images can be used to reconstruct finer
detail than exists in any one frame. We review related work in Section 2,
which includes a number of experiments on obtaining 2.5m super-resolution from
Sentinel-2.

50 cm grayscale |  Sentinel-2 grayscale |  Sentinel-2 grayscale super-resolution |  Building detection |  Road detection
---|---|---|---|---

Figure 2: Examples of building and road detection from Sentinel-2 imagery,
each covering an area of $192^{2}$ m2. The panels on the left show high-
resolution satellite imagery of the scene for comparison; although Sentinel-2
imagery has much lower level of detail in each frame, we are able to predict
fine-scale features of buildings and roads.

Figure 3: Estimation of the number of buildings in a tile, based on predicting
building centroids (left: high resolution image for comparison, centre: top of
Sentinel-2 stack input to the model; right: predicted centroid mask). This
method can obtain $R^{2}=0.91$ with respect to true counts even though
individual buildings cannot be discerned in the source imagery.

In our work, we attempt to extend the limit of fine detail that can be
recreated from a set of Sentinel-2 images, assessing the quality of building
and road presence predictions made at 50 cm resolution. To do this, we use a
teacher model with access to 50 cm resolution imagery to create training
labels for a large dataset of worldwide imagery. We then train a student model
to reconstruct these labels given only a stack of Sentinel-2 images from the
corresponding places and times (see Figure 1). We find that we are able to
retain much – though not all – of the accuracy of the high-resolution teacher
model: our Sentinel-2-based building segmentation has 78.3% mIoU, compared to
the high-resolution-based teacher performance of 85.3% mIoU. We found that
this accuracy level was equivalent to what could be achieved by a high-
resolution model using a single frame of 4 m resolution imagery (see Table 5),
though also noting through visual inspection that the segmentation quality
sometimes far exceeds that (see Figure 2 and Section 10).

Finally, we describe a method using our framework for counting the number of
individual buildings in a patch, based on predicting the locations of building
centroids. We can obtain $R^{2}=0.91$ compared to true building counts, again
capturing much of the performance of the high-resolution teacher model
($R^{2}=0.95$). The predicted counts are surprisingly accurate even when
buildings are very small relative to the raw Sentinel-2 resolution, or close
together. Some examples are shown in Figures 3 and 15.

This work extends the range of analysis tasks that can be carried out with
freely available Sentinel-2 data. As well as using it to improve the Open
Buildings dataset111https://sites.research.google/open-buildings, we plan to
make model resources available for social benefit purposes in the near future.

##  2 Related work

Several researchers have noted the potential of using freely available but
relatively low-resolution remote sensing imagery, such as Sentinel-2, to
obtain insights about buildings and other features with previously
unobtainable spatial and temporal scope. As well as experimental results,
practical data is already being produced from such systems: for example,
Sentinel-2 super-resolution from 10 m to 2.5 m has been used to create
buildings data for 35 cities across China [1].

Super-resolution, the task of reconstructing a high-resolution image from one
or more low-resolution images, has been widely applied to photographic images,
commonly with generative models such as GANs. There is existing work on
applying this type of model in remote sensing imagery [2], although there is a
risk of increasing resolution at the expense of introducing spurious details.
Certain characteristics of remote sensing imagery, and Sentinel-2 in
particular, can be exploited to aid hallucination-free super-resolution. Alias
and shift between sensor channels are undesirable for the purposes of
visualization, but provide a useful signal for super-resolution [3, 4, 2]. It
is also possible to exploit aspects of the physical design of the Sentinel-2
satellites, to obtain 5 m super-resolved training data for specific areas with
detector overlap [5].

Other results in the literature indicate that it is not necessary to
explicitly model these alias and shift effects to obtain hallucination-free
super-resolution, and that this can be learned directly from data with
standard semantic segmentation models. U-Net is popular in remote sensing
analysis and has been successfully applied in this setting. Simply up-sampling
medium resolution imagery and trying to extract building detections gives
increased effective resolution, which can then be used for building detection
[6]. Super-resolution followed by Mask-RCNN was used to detect buildings in
locations across Japan [7]. Another two stage model, SRBuildingSeg, carries
out super-resolution followed by building segmentation [8]. A similar two
stage setup based on U-Net was used to detect buildings using Sentinel-2 in
Spain [9].

Another principle useful for remote sensing super-resolution is to take
multiple images to generate a single high-resolution output. Physical features
such as buildings and roads change slowly in comparison to the Sentinel-2
revisit time of 5 days, so an image stack is likely (though not guaranteed) to
show the same scene. HighResNet [10] is an architecture for fusing a temporal
sequence of lower-resolution remote sensing images to predict a single,
higher-resolution image. This works with a convolutional architecture to fuse
the images together, and a loss function which can account for differences in
alignment between the low-resolution input and high-resolution training
labels. PIUNet (permutation invariance and uncertainty in multitemporal image
super-resolution) is an architecture for multiple–image super-resolution, used
to increase Proba-V resolution from 300 m to 100 m [11]. Most of this existing
work is based on fully convolutional architectures, however transformer-based
models have also been used [12]. The approach in our work is also to use
multiple images to predict high-resolution features, although we do not have
an intermediate super-resolution step and instead train end-to-end.

##  3 Training setup

Figure 4: Teacher-student setup in this work. The student model is trained to
reproduce the same outputs as a high-resolution model using 50 cm resolution
imagery, but using only a stack of Sentinel-2 images at 10 m resolution.

We propose a method for semantic segmentation from a stack of low-resolution
Sentinel-2 images at a much higher resolution than the input resolution. Our
setup consists of a teacher model and a student model, with inputs and outputs
as shown in Figure 4. The teacher model is trained on high-resolution (50 cm)
satellite imagery and outputs high-resolution semantic segmentation confidence
masks for buildings, building centroids, and roads. The student model tries to
mimic the output of the teacher model using only a stack of low-resolution
imagery obtained from Sentinel-2 of the same location.

Our method is an end-to-end super-resolution segmentation model, so that
instead of first carrying out super-resolution on the image and then running
semantic segmentation on the output, we make the model predict a high-
resolution semantic segmentation mask directly from low-resolution inputs.

The input to the model is a set of low-resolution Sentinel-2 frames
$\mathrm{LR}_{i,t},\in\mathbb{R}^{H\times W\times C}$, arranged in a stack of
time frames from $t=1,\ldots,32$ (where $H\times W$ are spatial dimensions and
$C$ is the number of input channels). The input channels include both imagery
and metadata, as described in Section 4. The output of the model is
$\mathrm{SR}_{i}\in\mathbb{R}^{\gamma H\times\gamma W\times C^{\prime}}$,
where $\gamma$ is the upscaling factor and $C^{\prime}$ is the number of
output channels. We denote the label by $\mathrm{HR}_{i}\in\mathbb{R}^{(\gamma
H+\Delta y_{max})\times(\gamma W+\Delta x_{max})\times C^{\prime}}$, where
$\Delta x_{max}$ and $\Delta y_{max}$ are margins corresponding to the maximum
allowed translation in $x$ and $y$ direction respectively (Section 6).

Each output channel can potentially represent a separate task. We train models
with up to four output channels: building semantic segmentation, roads
semantic segmentation, building centroids, and super-resolution grayscale
image, as shown in Figure 4. We use the building centroid predictions with one
extra step to calculate building counts. We predict a super-resolution
grayscale image as one of the output channels as this helps with image
registration during training and evaluation.

At a high level, our model employs an encoder and a decoder. The encoder
encodes each low-resolution image independently. The decoder takes a fused
representation of these encodings and applies successive upsampling to it to
output at target resolution. The model architecture is described in more
detail in Section 5.

##  4 Dataset

The dataset for training and evaluation consists of low-resolution image
stacks and high-resolution label pairs
$\left\\{\mathrm{LR}_{i,t=1:32},\mathrm{HR}_{i}\right\\}$. To generate labels
for the temporal stack we fetch Sentinel-2 imagery stacks at locations where
we also have high-resolution imagery available. We fetch the stack of
Sentinel-2 imagery such that the high-resolution image corresponds to the
middle of the low-resolution stack, i.e. between 16th and 17th time frame in a
32 frame stack (see Figure 5). Labels for the primary tasks of buildings and
roads semantic segmentation are generated as follows: labels for the training
split are per pixel building and roads presence confidence scores output by
the teacher model, whereas labels for the evaluation split are binary
segmentation masks obtained from human drawn polygons of buildings and roads
on high-resolution images (see Figure 6). Since we generate training labels
using a teacher model we can generate a large number of training examples,
only limited by the amount of high-resolution imagery available and
computational constraints.

Figure 5: Generation of training and evaluation data for the low-res student
model. The low-resolution stack is constructed so that the time of the high-
resolution image falls between the 16th and 17th frames. Figure 6: Comparison
of the model-derived teacher label (centre) and human-derived label (right)
for one evaluation example.

The training split of the dataset consists of about 10 million stacks of
Sentinel-2 imagery sampled globally. Examples are sampled randomly at uniform
first and then subsampled by building presence using the high-resolution
teacher model so that approximately 90% of the examples have buildings. The
geographical distribution of examples is shown in Figure 8. A consequence of
our sampling strategy is that the distribution of our dataset is roughly
correlated with population density. The validation split consists of 1165
Sentinel-2 imagery stacks, which are all within the continent of Africa and
chosen to include a range of urban and rural examples of different densities,
as well as settings of humanitarian significance such as refugee settlements.
Training examples close to validation examples (within some radius) are
discarded. For simplicity we do not perform any sampling based on roads. Often
roads tend to be present in the vicinity of buildings, but this is not always
true.

Sentinel-2 imagery is fetched from the Harmonized Sentinel-2 Top-Of-Atmosphere
(TOA) Level-1C collection [13] on Earth Engine. Each Sentinel-2 timeframe
consists of 13 bands, each upsampled to 4 meters per pixel regardless of their
native resolution. Additionally, we add metadata for each frame, as described
in Section 4.5 in more detail.

###  4.1 Teacher labels

The teacher model which is used to label the large worldwide training set is
the same as used to generate the Open Buildings dataset, trained along the
lines described in [14], with a few refinements for improved accuracy: these
include the collection of extra training data in areas where earlier models
had poor performance, ensembling of models trained on different partitions of
the training data (distilled back to a single model) and use of the HRNet [15]
architecture instead of U-Net. The geographical distribution of human-labelled
training data for buildings is shown in Figure 7, reflecting the current Open
Buildings coverage across Africa, South and Southeast Asia, and Latin America.
Training is done in a similar way for road detection as for building
detection.

Figure 7: Geographical distributions of high-resolution, human labelled data
on which the teacher model was trained, and human-labelled evaluation data on
which we report all evaluation metrics. Counts per level 5 S2 cell. Figure 8:
Geographical distribution of teacher-labelled training data used to train the
student model. Counts per level 5 S2 cell.

###  4.2 Registration

Neither low-resolution Sentinel-2 images nor the high-resolution image are
registered to any reference frame, and as such, both the input imagery stack
as well as the labels are potentially misaligned relative to each other. Given
the resolution of Sentinel-2 imagery, misalignments of a few pixels in image
space amount to tens of meters on the ground. We rely on the model being able
to implicitly align the low-resolution input frames and to facilitate such
implicit alignment, inspired by [10], we tried pairing each input frame with
the frame closest in time to the high-resolution image used to derive labels,
i.e. the 17th frame. The model output however can still be misaligned relative
to the labels, and in order to not penalise the model for that we manually
align labels to model output during loss computation, as described in Section
6. We also do similar alignment in evaluation.

###  4.3 Clouds

To filter cloudy images we discard Sentinel-2 images that have one or more
pixels in the opaque cloud mask set. An opaque cloud mask is derived from the
QA60 band in the Sentinel-2 TOA collection (10th bit). QA60 band however is
inaccurate and does not always capture all clouds (see example in Figure 9).
We retain Sentinel-2 images even if QA60 indicates they have cirrus clouds,
because often these images still looked useful. Moreover, cloudy high-
resolution images are especially problematic since they introduce noise into
the labels. Therefore both the training and evaluation datasets are filtered
to have cloud-free high-resolution imagery. During inference, we only use the
QA60 band to filter out timeframes with opaque clouds.

###  4.4 Orthorectification

Label transfer from high resolution to low resolution works well only if both
the satellites have the same viewing angle or both the imagery collections are
orthorectified. The high-resolution imagery collection we use is however not
orthorectified (see example in Figure 10) and as such, labels generated using
such imagery do not always align with Sentinel-2 imagery. Particularly
problematic are areas with tall buildings where building roofs (which is what
our teacher model is trained to detect) appear in different positions
depending on the satellite viewing angle. Orthorectification of satellite
imagery however requires an accurate elevation model – something that is not
available globally. In the absence of orthorectification, labels tends to be
more accurate for rural areas than for urban areas with taller buildings.

Figure 9: Inaccurate QA60 band cloud mask. Figure 10: Non-orthorectified
image: tall structure.

###  4.5 Sentinel-2 data processing

Due to the way imagery assets are tiled in the Harmonized Sentinel-2 Top-Of-
Atmosphere (TOA) Level-1C collection on Earth Engine, the imagery stack
fetched for a given location can potentially contain duplicate imagery [16].
We therefore perform deduplication of the imagery stack as follows: we group
images by datatake ids and for each datatake id we take the highest processing
baseline. Datatakes correspond to a swath of imagery taken by the Sentinel-2
satellite and often cover a very large surface area of the earth (up to 15,000
square kilometers). The processing baseline is effectively a set of
configurations used to post process raw imagery acquired from the satellites
and are updated regularly by the European Space Agency (ESA).

We feed in metadata associated with each frame to the model using a simple
approach where scalar value for each metadata item is broadcast to image
spatial dimensions and then appended to the image in the channel dimension. As
such, each scalar metadata value adds an additional channel to the input. The
following metadata are used for each frame: normalized time relative to the
17th time frame, mean incidence azimuth angle, mean incidence zenith angle,
mean solar azimuth angle, mean solar zenith angle, latitude, and longitude.
The time of image acquisition is normalised by dividing the time duration
relative to the 17th frame by ten years in seconds, and scale the other
features to the [0, 1] range. We do not perform any data augmentation other
than random cropping, which we found is sufficient to prevent overfitting due
to the very large dataset size.

To account for the fact that, in certain cloudy parts of the world it is not
always possible to acquire a stack of 32 Sentinel-2 images in a given time
span, during training we randomly truncate and pad the image stack on both
ends. More formally, for both halves of the stack, with some probability $p$
we generate a padding of length $l$ distributed uniformly in $[1,16]$. This
helps to make the model robust to missing frames at inference time.

##  5 Model

At a high level our student model consists of a shared encoder module that
encodes each LR frame separately, a simple mean based encoding fusion and a
decoder module that up-samples the fused encoding (Figure 13). Below we
describe encoder and decoder modules in detail.

###  5.1 Encoder

We employ the HRNet [15] architecture (see Figure 11) as an encoder without
making any significant modifications, except for the root block. The original
root block in HRNet downsamples the input image by a factor of 4, which we
find to be too aggressive for LR input images that are already quite low
resolution. To address this, we remove the 4x downsampling by reducing stride
from 2 to 1 on the two 3x3 convolutions in the first block of HRNet. As a
result, the encoder outputs features at the same spatial resolution as the
input. We pre-train the HRNet encoder on ImageNet. To adapt filter weights of
the first convolution trained on 3 channel RGB ImageNet input to more
channels, for each filter we take the mean of the weights across channels and
replicate the mean across the target channel dimensions (=13+7, corresponding
to all bands in a Sentinel-2 image and all metadata we pass in).

Figure 11: HRNet architecture, as it would be applied to a single frame of
imagery. The network consists of 4 different stages and each stage consists of
blocks that correspond to features at different resolutions. Each Sentinel-2
input has 13 image channels and 7 metadata channels.

###  5.2 Cross-time information fusion

In order to make the model learn not just spatial features but also the
temporal relationship between spatial features we experimented with cross-time
information fusion as features go through the various stages of HRNet encoder.
For each stage and block (corresponding to features at various depths and
resolutions), spatial features are fused across time using a depth-wise
convolution in a residual fashion (Figure 12).

For a given block and stage, features for pixel $(u,v)$ across all 32 time
steps are passed through a $1\times 32$ depth-wise convolution (with depth
multiplier $32$) to transform a $1\times 32\times c$ feature tensor to
$1\times 1\times 32c$. The tensor is then reshaped to $1\times 32\times c$
before being passed to $1\times 1$ point-wise convolution to obtain $1\times
32\times c$ fused features. We reshape the tensor before the point-wise
convolution to reduce the number of convolution parameters. Otherwise, this
setup is equivalent to a 1D depth-wise separable convolution across time.

Note that this fusion approach is much more efficient than 3D convolution or
ConvLSTM in terms of the number of parameters. This is because our approach
only fuses features across time for a single pixel, and using depth-wise
convolution (with reshaping before point-wise convolution) helps to reduce the
number of parameters even further.

The original features are then added to the fused features as a residual
connection. More details about pairing schemes and experimental results are
given in Section 7.3.

Figure 12: cross-time fusion of features. For each stage and block
(corresponding to features at various depths and resolutions) spatial features
for pixel $(u,v)$ across all 32 time steps are fused using depth-wise
convolution in a residual fashion.

###  5.3 Decoder

The decoder module in our model consists of a number of decoder blocks that
each upsample the input by a factor of 2. We use the same decoder block as in
our previous technical report [14] which consists of x2 (batch normalization,
ReLU, convolution) followed by another (batch normalization, ReLU), fusion
with residual connection to the input and finally an upconvolution, as
illustrated in Figure 14. We use transposed convolution or deconvolution for
up convolution. The fused encoding is passed through 3 decoder blocks with
widths 360, 180 and 90 respectively. Finally the upsampled features are passed
through a (convolution, batch normalization, ReLU) block before being passed
through a final convolution that outputs the desired number of output
channels. There is one output channel for each task: building segmentation,
road segmentation, centroid detection and image super-resolution. Since the
input Sentinel-2 frames are the same for each of these tasks, a single model
can be trained in a multi task fashion to perform these tasks simultaneously.
Apart from reducing cost and time, the different tasks could possibly serve as
regularization for the other tasks to produce better results.

Figure 13: Overall model architecture. Each of the input images is encoded
separately using a shared encoder. Encodings are then averaged together (which
collapses the time dimension) and passed on to the decoder that upsamples the
fused encoding by a factor of 8. Figure 14: Residual Decoder block with width
$n$ that upsamples the input by a factor of 2.

###  5.4 Building counts via centroid prediction

In some practical situations, people need building counts in an area rather
than a segmentation mask of building presence. We investigated various ways of
trying to estimate the count of buildings in an image tile from the Sentinel-2
image stack. The most effective method we found is based on first predicting
the positions of building centroids, inspired by [17] and other work on object
counting e.g. for crowds. For the centroid detection task, the labels are
Gaussian ‘splats’, as shown in Figure 15, where each splat is the same size
regardless of the size of the building. This means that we expect a roughly
constant response for each building, and we can derive an estimate count by
simply summing the model output over the spatial dimensions and then dividing
by a constant scaling factor.

We therefore derive the count without requiring any post-processing such as
peak detection or non-maxima suppression. For a given tile $i$, the building
count, $\hat{n}_{i}$ can be estimated as follows:

$\hat{n}_{i}=\dfrac{1}{K}\sum_{u,v}C_{u,v}^{i}\ ,$ (1)

where $C_{u,v}^{i}$ is model output at pixel $u,v$ (for the channel
corresponding to centroid labels) and $K$ is the scaling factor. $K$ is
derived from the training split by least squares regression.

The centroid detection task is incorporated into the overall model as one of
the output channels. Building instances from the teacher model, above a
certain score threshold, are used to generate centroid labels. See Figure 15
for example labels and model output. We detail some of the results obtained
with this approach in Section 7.7.

Figure 15: From the left to the right high-resolution image, the first images
of the corresponding Sentinel-2 stack, the instance masks derived from teacher
model, corresponding centroid labels, and model output.

##  6 Loss function

We employ per pixel Kullback-Leibler Divergence (KLD) loss between the label
and the prediction, defined as follows:

$KLD(y_{i},\hat{y}_{i})=\left(y_{i}\log\left(\frac{y_{i}}{\hat{y}_{i}}\right)+\epsilon\right)^{\gamma}\
,$ (2)

where $y_{i}$ and $\hat{y_{i}}$ are the label and prediction for pixel $i$
respectively and $\gamma$ is the focal term. Both $y_{i}$ and $\hat{y_{i}}$,
originally $\in[0,1]$, are clipped to $[\epsilon,1-\epsilon]$, where
$\epsilon=1e-7$, to avoid division by zero errors. The focal term $\gamma$
varies the importance given to hard (misclassified) examples. Following a
hyper-parameter sweep we set $\gamma$ to 0.25. We also add $\epsilon$ to KLD
before exponentiation by $\gamma$ to guard against an undefined gradient.

Since the label and the model output could be misaligned, and in order to not
penalize the model due to misalignment, we do an exhaustive neighborhood
search based registration between the label and the prediction. For this we
assume a translation model and find the $(\Delta x,\Delta y)$ translation that
minimizes the mean squared error (MSE) between the label and prediction. The
label is then shifted by $(\Delta x,\Delta y)$ to register it with the model
output, before the loss in (2) is computed. This exhaustive neighborhood
search is however computationally expensive. A possible alternative would be
to also learn registration as proposed in [10], where a separate registration
module takes in (prediction, label) pair as input and predicts a $k\times k$
kernel which when convolved with the model output aligns it with the label.

We crop the label such that we keep a margin around it whose size is set to
the maximum translation $(\Delta x_{max},\Delta y_{max})$ allowed in a given
direction (Figure 16). Once the label is shifted relative to the model output
we crop the label to the output dimension so that both the label and model
output have the same size.

Figure 16: Alignment at loss: Label with margin whose size is equal to the
maximum allowed translation along a given image dimension.

To make alignment more robust, we also make the model predict a super-resolved
grayscale image. The label for that task is the 50 cm image converted to
grayscale. We use similar registration logic in both training and evaluation.

##  7 Experiments

In our experiments, we report the mean Intersection over Union (mIoU) for a
binary segmentation setup, where the foreground is either buildings or roads.
Instead of using a fixed threshold (such as 0.5) to convert the per-pixel
confidence values output by the model into a segmentation mask, we use a
threshold in the range $[0,1]$ and a dilation kernel size that maximizes the
mIoU. We then report the maximum mIoU obtained. For the auxiliary task of tile
level building count, we report mean absolute error (MAE) and the coefficient
of determination, R2.

The numbers between different tables are not necessarily comparable due to
slight differences in configurations. Our evaluations focus on the building
detection task and we do not report metrics for the road detection and image
super-resolution tasks.

###  7.1 Training Details

Our models are trained with the Adam optimizer [18] using a constant learning
rate of $3\times 10^{-4}$ for 500,000 steps with batch size 256. During
training images are cropped randomly to size of 512$\times$512\. For
evaluation, the random crop is replaced by a center crop but in the end
metrics are calculated on a further center crop of 384$\times$384\. All the
models are initialized with the checkpoints of a model trained on ImageNet
[19].

###  7.2 Handling of multiple timeframes

We explored two options for handling multiple Sentinel-2 timeframes in the
model. In the first one we concatenated all timeframes in the channel
dimension. In the second approach we passed each timeframe through a shared
encoder (see Section 5). The second approach seems to require cross-time
information fusion (see Section 5.2) and/or a pairing scheme for timeframes
(see Section 7.3) to exceed the performance of the first approach (see Table
1). Both approaches used slightly modified off-the-shelf HRNet [15] and U-Net
[20]. See Section 5.1 for a description of how HRNet was modified, in case of
U-Net the modifications were equivalent. Unless stated otherwise, in all
remaining experiments we use the second approach with cross-time information
fusion and no pairing scheme.

Table 1: Comparison of two approaches to handling multiple timeframes: concatenating the timeframes in the channel dimension and separate encoding using a shared encoder. The second approach requires additional modifications to exceed the performance of the first approach. Model | mIoU (Buildings)
---|---
U-Net concatenate timeframes in channel dimension | 72.7
HRNet concatenate timeframes in channel dimension | 73.8
HRNet separate timeframe encoding | 70.7
HRNet separate timeframe encoding + cross-time fusion | 76.7
HRNet separate timeframe encoding + Pairing (17th frame) | 76.9
HRNet separate timeframe encoding + Cross-time fusion + Pairing (17th frame) | 76.4

###  7.3 Pairing schemes

Inspired by HighResNet [10], where each low-resolution frame is paired with
the per-pixel median of the stack, we experimented with various pairing
schemes. We find that pairing each timeframe with the 17th timeframe (that is
closest in time to the teacher label) provides a significant boost in
performance over model without pairing. In comparison, model trained with no
pairing but with cross-time communication (Section 5.2) performs slightly
worse (see Table 1). Surprisingly enough, model trained with both pairing
(17th frame) and cross-time communication performs worse than that with just
pairing. We also experimented with pairing with either the median (as is done
in HighResNet), mean or second darkest timeframe and find none of these
pairing schemes to be better than pairing with the 17th frame.

###  7.4 Resolution sensitivity analysis

To understand the influence of input, output and target resolutions on model
performance we carry out the following resolution sensitivity experiments:

#### Input resolution

We altered the resolution of the input while keeping the model output and
label at 50 cm resolutions. The effective resolution of the Sentinel-2 input
is 10 m or lower depending on the specific band, but we find that upsampling
it with a simple image resize has non-trivial impact on model performance. For
example, training on 4 meter inputs provides a sizable improvement over
training on 8 meter inputs (Table 2). However, increasing the input resolution
to 2 meters does not provide similar benefit. The models in Table 2 are
trained on eight Sentinel-2 timeframes because of the computational
constraints associated with training at 2 meter input resolution. For training
on the stack of 32 images, we only explored training on 8 and 4 meter inputs.
The performance improvement associated with training on 4 meter inputs is also
observed in this case.

Table 2: Sensitivity analysis of the Sentinel-2 input resolution (8 timeframes). Input resolution | mIoU (Buildings)
---|---
8 meters | 74.3
4 meters | 75.6
2 meters | 75.8

#### Output resolution

For output resolution sensitivity experiments we alter the effective
resolution of the output while keeping the model input at 4 m and the label at
50 cm resolutions.

As described in Section 5, to produce super-resolved outputs we add residual
upsampling blocks akin to the decoder blocks used in U-Net to progressively
upsample the input. We investigate the usefulness of these blocks by
progressively replacing them with a bi-linear resize operation. Although the
output of the model is at the same resolution as the label, the effective
resolution is actually much lower since the resize operation does not add any
new details. Table 3 shows that using the upsampling block to increase the
effective resolution leads to an increase in the performance of the model. At
4 meters, no upsampling blocks are used, resulting in lower performance
compared to other models that make use of them.

Table 3: Sensitivity analysis of the Sentinel-2 output resolution demonstrating usefulness of upsampling blocks (32 timeframes). Effective output resolution | mIoU (Buildings)
---|---
4 meters | 76.1
2 meters | 76.6
1 meter | 76.8
50 centimeters | 77.3

#### Label resolution

For label resolution sensitivity experiment we alter the effective resolution
of the label while keeping the model input at 4 m and output at 50 cm
resolutions. Unsurprisingly, we see that the model performs better when
trained on higher resolution labels, as shown in Table 4.

Table 4: Sensitivity analysis of label resolution (32 timeframes). Effective label resolution | mIoU (Buildings)
---|---
4 meters | 76.2
2 meters | 75.6
50 centimeters | 76.8

#### Comparison to single-frame detection with varying image resolution

We wished to understand the performance of our model in comparison to a high-
resolution building detection model operating on a single frame of imagery.
That is, if only one image was available, what resolution would it need to be
in order to get the same accuracy of detection as we obtain with 32 Sentinel-2
frames.

To do this, we trained a model similar to our teacher model on downsampled
high-resolution images. The images are downsampled to a target lower
resolution and resized back to the original image dimensions. Thus the
dimensions of the images are the same but the information content is reduced.
We find that at 4 meter resolution, this model’s performance matches the
performance of the Sentinel-2 super-resolution model (Table 5). The 50 cm
model has comparable metrics to the model that was used to generate teacher
labels for training the Sentinel-2 models.

Table 5: Resolution sensitivity for a model trained on a single high-resolution image. For comparison, the last row represents our best Sentinel-2 super-resolution segmentation model. It has comparable mIoU to the 4 meter model. Effective image resolution (cm) | mIoU (Buildings)
---|---
1000 | 67.4
500 | 75.8
400 | 78.1
300 | 80.5
200 | 82.9
100 | 84.8
50 | 85.3
Best Sentinel-2 model | 78.3

###  7.5 Training set size

To quantify the impact of dataset size on model performance we trained our
model on subsets of training data of different sizes. Table 6 shows monotonic
improvement in performance as training data size increases. The model trained
on the full set of 10 million images outperforms the model trained on 10% of
the training set by 2 percentage points.

Table 6: Building segmentation performance as a function of fraction of training data size used (32 timeframes). Fraction of training data | mIoU (Buildings)
---|---
1% | 69.1
5% | 73.4
10% | 74.7
100% | 76.6

###  7.6 Number and position of Sentinel-2 timeframes

Each training example contains 32 Sentinel-2 timeframes with the corresponding
high-resolution image located temporally somewhere between the 16th and 17th
timeframes. The timeframes are sorted by time. This means that the model gets
to see 16 frames in the past and 16 in the future around the high-resolution
image time. We carry out a series of experiments to determine the relationship
between the number of timeframes and the performance of the model.

We find that increasing the number of timeframes leads to a monotonic increase
in the performance of the model. In Table 7, the model trained on all 32
timeframes outperforms the model trained on a single timeframe by 5 percentage
points. As a corollary to this, we train a model on a single timeframe
duplicated 32 timeframes. This model achieves a mean intersection over union
of 71.7% compared to 76.7% of the model trained on all timeframes. This shows
that the additional timeframes provide useful information to the model. The
conditions at the time a snapshot is taken is different for the 32 timeframes,
thus there must be features that are captured in some frames but absent in
others. Using multiple frames allows the model to draw from the information
available in each time frame to make a single prediction. The use of multiple
frames can also be understood as producing a ‘dither’ effect. Since noise
should be randomly distributed across the image, shifts between consecutive
frames allow the model to isolate the signal from the noise to produce a
super-resolved output.

Table 7: Model performance as a function of the number of timeframes in the Sentinel-2 stack. Number of timeframes | mIoU (Buildings)
---|---
1 | 71.7
2 | 73.2
4 | 74.5
8 | 75.8
16 | 76.1
32 | 76.7

We also explore how sensitive the model is to having access to future or past
timeframes by training the model on only past or future frames. In Table 8,
the model that only gets to see 16 past timeframes is better than a model that
sees only 16 future timeframes, and roughly comparable to one that sees 8 past
and 8 future timeframes.

###  7.7 Building counts via centroid prediction

On the auxiliary task of tile level count prediction via centroid prediction
as outlined in Section 5.4 we observe a very high correlation between
predicted building count and the true count (see Figure 17). We were able to
achieve a coefficient of determination ($R^{2}$) of 0.912 and MAE of 5.67 on
human labels. For comparison, model trained to do the same but using 50 cm
imagery was able to achieve $R^{2}$ of 0.955 and MAE of 4.42 on human labels.

Figure 17: Regression plots (in linear and logarithmic scale) and metrics for
the buildings count from just a stack of low-resolution images. Counts are
made on tiles of size $192^{2}$ m2 on the ground.

##  8 Discussion

In this work, we present an end-to-end super-resolution segmentation framework
to segment buildings and roads from Sentinel-2 imagery at a much higher
effective resolution than the input imagery. To this end we demonstrate label
transfer from a high-resolution satellite to a low-resolution satellite that
both image the same surface of the earth. We show that such label transfer not
only allows generation of more accurate and fine detailed labels but also
enables automatic label generation through a teacher model trained to perform
the same task but using imagery from a high-resolution satellite. This
significantly reduces the amount of high-resolution imagery needed for certain
analysis tasks, such as large scale building mapping.

Unfortunately, our approach has some limitations. First, it still assumes that
one has access to a good amount of high-resolution images to first train a
teacher high-resolution model and then also to generate a large dataset to
train a student model as described in Section 4. However, this is a one-time
cost. Additionally, with the advancement in deep learning, models can now be
trained in a much more data efficient manner [21]. Furthermore, there are many
publicly available datasets and models that can be used to bootstrap different
tasks [22]. Second, it assumes that for a given location one can assemble a
deep stack (i.e. 32) of cloud-free Sentinel-2 images. This actually can be
quite problematic if the stack has to be centered timewise around a fairly old
or recent date, and especially if the location is cloudy; in humid locations
such as Equatorial Guinea we noticed that for some very cloudy locations a 32
timeframe stack can span more than 2 years.

The results we present in this technical report are preliminary and we believe
that there is still room for improvement. For example, we did not evaluate the
model on change detection, and do not have metrics on how many Sentinel-2
timeframes with change does the model need to see to detect these changes
consistently. Additionally, recent advances in super-resolution using
generative AI models such as GANs [2] and diffusion models [23] have shown
potential. These approaches have however been limited to single image super
resolution of natural images and have not been extensively used for remote
sensing tasks where multiple low-resolution images of the same location are
available. As such, a promising area of future research could be exploration
of use of some of these approaches for multi-frame super-resolution
segmentation tasks in remote sensing.

###  8.1 Social and ethical considerations

We believe that timely and accurate building information, particularly in
areas which have few mapping resources already, are critically important for
disaster response, service delivery planning and many other beneficial
applications. However, there are potential issues with improvements in remote
sensing analysis, both in terms of unintended consequences and from malicious
use. Where such a model is used as a source about information on human
population centres, for example during emergency response in a poorly mapped
and inaccessible area, any false negatives could lead to settlements being
neglected, and false positives lead to resources being wasted. Particular
risks for the kind of Sentinel-2 based analysis we describe include
settlements consisting of very small buildings made of natural materials, and
buildings in deserts, both of which are challenging.

Table 8: Effect of the temporal relationship between the high-resolution label and the Sentinel-2 stack, with 16 frames either at the end of the stack (after the time of the high-resolution label), the start (before the label), or split across both. Method | mIoU (Buildings)
---|---
Only future | 75.8
Only past | 76.6
Equal past & future | 76.3

##  9 Acknowledgements

We thank Sergii Kashubin, Maxim Neumann and Daniel Keysers for feedback which
helped to improve this paper.

## References

  * [1] Lin Feng, Penglei Xu, Hong Tang, Zeping Liu, and Peng Hou.  National-scale mapping of building footprints using feature super-resolution semantic segmentation of sentinel-2 images.  GIScience & Remote Sensing, 60(1):2196154, 2023.
  * [2] F Pineda, V Ayma, and C Beltran.  A generative adversarial network approach for super-resolution of sentinel-2 satellite images.  The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 43:9–14, 2020.
  * [3] Ngoc Long Nguyen, Jérémy Anger, Lara Raad, Bruno Galerne, and Gabriele Facciolo.  On the role of alias and band-shift for sentinel-2 super-resolution.  arXiv preprint arXiv:2302.11494, 2023.
  * [4] Mario Beaulieu, Samuel Foucher, Dan Haberman, and Colin Stewart.  Deep image-to-image transfer applied to resolution enhancement of sentinel-2 images.  In IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium, pages 2611–2614. IEEE, 2018.
  * [5] Ngoc Long Nguyen, Jérémy Anger, Axel Davy, Pablo Arias, and Gabriele Facciolo.  L1bsr: Exploiting detector overlap for self-supervised single-image super-resolution of sentinel-2 l1b imagery.  arXiv preprint arXiv:2304.06871, 2023.
  * [6] Jonathan Prexl and Michael Schmitt.  The potential of sentinel-2 data for global building footprint mapping with high temporal resolution.  In 2023 Joint Urban Remote Sensing Event (JURSE), pages 1–4, 2023\.
  * [7] Shenglong Chen, Yoshiki Ogawa, Chenbo Zhao, and Yoshihide Sekimoto.  Large-scale individual building extraction from open-source satellite imagery via super-resolution-based instance segmentation approach.  ISPRS Journal of Photogrammetry and Remote Sensing, 195:129–152, 2023.
  * [8] Lixian Zhang, Runmin Dong, Shuai Yuan, Weijia Li, Juepeng Zheng, and Haohuan Fu.  Making low-resolution satellite images reborn: a deep learning approach for super-resolution building extraction.  Remote Sensing, 13(15):2872, 2021.
  * [9] C. Ayala, C. Aranda, and M. Galar.  Pushing the limits of sentinel-2 for building footprint extraction.  In IGARSS 2022 - 2022 IEEE International Geoscience and Remote Sensing Symposium, pages 322–325, 2022.
  * [10] Michel Deudon, Alfredo Kalaitzis, Israel Goytom, Md Rifat Arefin, Zhichao Lin, Kris Sankaran, Vincent Michalski, Samira E Kahou, Julien Cornebise, and Yoshua Bengio.  Highres-net: Recursive fusion for multi-frame super-resolution of satellite imagery.  arXiv preprint arXiv:2002.06460, 2020.
  * [11] Diego Valsesia and Enrico Magli.  Permutation invariance and uncertainty in multitemporal image super-resolution.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–12, 2021\.
  * [12] Tai An, Xin Zhang, Chunlei Huo, Bin Xue, Lingfeng Wang, and Chunhong Pan.  Tr-misr: Multiimage super-resolution based on feature fusion with transformers.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 15:1373–1388, 2022.
  * [13] Harmonized Sentinel-2 Top-Of-Atmosphere (TOA) Level-1C Earth Engine collection.  https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_HARMONIZED.
  * [14] Wojciech Sirko, Sergii Kashubin, Marvin Ritter, Abigail Annkah, Yasser Salah Eddine Bouchareb, Yann N. Dauphin, Daniel Keysers, Maxim Neumann, Moustapha Cissé, and John Quinn.  Continental-scale building detection from high resolution satellite imagery.  CoRR, abs/2107.12283, 2021.
  * [15] Jingdong Wang, Ke Sun, Tianheng Cheng, Borui Jiang, Chaorui Deng, Yang Zhao, Dong Liu, Yadong Mu, Mingkui Tan, Xinggang Wang, Wenyu Liu, and Bin Xiao.  Deep high-resolution representation learning for visual recognition, 2020\.
  * [16] Jan Jackson.  Deduplicate Datatake.  https://forum.step.esa.int/t/duplicate-observations-in-the-same-archive/1793/4, 2016\.  [Online; accessed 21-August-2023].
  * [17] Xingyi Zhou, Dequan Wang, and Philipp Krähenbühl.  Objects as points.  arXiv preprint arXiv:1904.07850, 2019.
  * [18] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  International Conference on Learning Representations, 2015.
  * [19] J. Deng, Wei Dong, R. Socher, Li-Jia Li, K. Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 248–255, 2009.
  * [20] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  MICCAI, 2015.
  * [21] Yezhen Cong, Samar Khanna, Chenlin Meng, Patrick Liu, Erik Rozi, Yutong He, Marshall Burke, David Lobell, and Stefano Ermon.  Satmae: Pre-training transformers for temporal and multi-spectral satellite imagery.  Advances in Neural Information Processing Systems, 35:197–211, 2022\.
  * [22] SpaceNet Challenge GitHub repository.  https://github.com/SpaceNetChallenge.
  * [23] Chitwan Saharia, Jonathan Ho, William Chan, Tim Salimans, David J Fleet, and Mohammad Norouzi.  Image super-resolution via iterative refinement.  arXiv:2104.07636, 2021.

##  10 Appendix

50 cm grayscale |  Sentinel-2 grayscale |  Sentinel-2 grayscale super-resolution |  Building detection |  Road detection |  Teacher building detection |  Teacher road detection
---|---|---|---|---|---|---

Figure 18: Urban examples of building and road detection output, each covering
an area of $192^{2}$ m2.

50 cm grayscale |  Sentinel-2 grayscale |  Sentinel-2 grayscale super-resolution |  Building detection |  Road detection |  Teacher building detection |  Teacher road detection
---|---|---|---|---|---|---

Figure 19: Rural examples of building and road detection output, each covering
an area of $192^{2}$ m2.

(a) |
---|---
(b) |
(c) |
(d) |
(e) |
(f) |
(g) |

|  50 cm grayscale |  Sentinel-2 grayscale |  Sentinel-2 grayscale super-resolution |  Building detection |  Road detection |  Teacher building detection |  Teacher road detection
---|---|---|---|---|---|---|---

Figure 20: Interesting examples of building and road detection output, each
covering an area of $192^{2}$ m2: (a) some regions can develop very rapidly,
the 50 cm image is 3 weeks older than the next available cloud-free Sentinel-2
image but new structures have already appeared during that timespan; (b) the
model struggles with separating many small building structures from each
other; (c) the model can be sensitive to roof color, in this example lighter
colored roofs have higher confidence; (d) model struggles on small buildings
interspersed with trees; (e) buildings in desert environment can be
challenging; (f) tall structures can have blurry predictions due to 50 cm
imagery often being at high off-nadir angle; (g) the road prediction can
sometimes be better than the teacher road prediction, in this case the model
is less confident the beach is a road.

2016-11-17 |  |  |
---|---|---|---
2019-02-07 |  |  |
2023-01-04 |  |  |

|  50 cm RGB |  Building detection |  Road detection
---|---|---|---

Figure 21: Urban growth in Bulbula in Oromiya Region in Ethiopia over time.
Leftmost column is high-resolution imagery for reference.

2016-01-31 |  |  |
---|---|---|---
2016-07-15 |  |  |
2017-02-28 |  |  |

|  50 cm RGB |  Building detection |  Road detection
---|---|---|---

Figure 22: The construction of Pagirinya refugee settlement in Eastern
Adjumani District, Uganda. Leftmost column is high-resolution imagery for
reference.


# Concurrent Learning Based Tracking Control of Nonlinear Systems using
Gaussian Process

Vedant Bhandari and Erkan Kayacan *This work was supported by UQ Early Career
Researcher grant (UQECR2057325).The authors are with the School of Mechanical
and Mining Engineering, The University of Queensland, 4072 Brisbane,
Australia. {v.bhandari, e.kayacan<EMAIL_ADDRESS>
###### Abstract

This paper demonstrates the applicability of the combination of concurrent
learning as a tool for parameter estimation and non-parametric Gaussian
Process for online disturbance learning. A control law is developed by using
both techniques sequentially in the context of feedback linearization. The
concurrent learning algorithm estimates the system parameters of structured
uncertainty without requiring persistent excitation, which are used in the
design of the feedback linearization law. Then, a non-parametric Gaussian
Process learns unstructured uncertainty. The closed-loop system stability for
the nth-order system is proven using the Lyapunov stability theorem. The
simulation results show that the tracking error is minimized (i) when true
values of model parameters have not been provided, (ii) in the presence of
disturbances introduced once the parameters have converged to their true
values and (iii) when system parameters have not converged to their true
values in the presence of disturbances.

##  I INTRODUCTION

The performance of feedback linearization depends on the accuracy of the
system model, and its parameter estimates in the linearization process and all
generalized external forces that change in different environments [1, 2, 3].
It is often difficult to exactly know the model parameters for any practical
system, with uncertainty in measurement always present [4, 5]. Uncertainty in
model parameters leads to a decreasing reference tracking performance [6].
Furthermore, changing operation environments introduces varying disturbances,
which cannot always be accounted for in the controller design [7].

These challenges mentioned above require feedback linearization to be combined
with more robust design techniques to account for model uncertainties and
adapt to various environments with varying disturbances. The proposed law by
[8] combines feedback linearization (FBL) control with Gaussian Process (GP)
to account for model uncertainty. Using online learning GP, the controller
learns the model mismatch between the actual model and the estimated model
used for the feedback linearization, reducing the model uncertainty and
leading to better reference tracking. However, no attempt is made to estimate
the system model’s true parameters and the effect of any external disturbances
is directly captured by the GP as a model mismatch. This is not ideal as the
model mismatch and disturbance cannot be distinguished. A similar proposal by
[9, 10] also uses GP-based tracking control, where the Gaussian process is
used offline to learn the model mismatch and any modelled dynamics together.
This presents the same issues as mentioned previously.

Figure 1: Schematic diagram for the developed control algorithm. Along with
the use of the traditional feedback linearization technique (FBL and state-
feedback controller blocks), a concurrent learning algorithm is used to
estimate the model parameters to be used in the feedback linearization
process, and an online non-parametric Gaussian Process and Robustness term are
used to account for unknown disturbances.

Concurrent learning approach was used and demonstrated significantly with
model reference adaptive control methods as initially proposed in [11].
Recently, it has also been used as an online parameter estimator [12]. It can
be used in any other nonlinear control methods. This paper will demonstrate
its applicability in the scheme of the FBL method – this is unique and has not
been shown in the literature.

The developed control law in this paper provides a two-stage controller
design, as shown in Fig.1. In the first stage, the system model parameters are
estimated to provide optimal tracking control with feedback linearization
using concurrent learning in a known and controlled environment. The learned
model parameters are then used with a GP-based online learning approach to
operate in environments with unknown disturbances in the second stage. This
approach’s key benefits are that the system model parameters can be accurately
estimated, which is critical in a model-based controller design process. The
unstructured uncertainties can then be captured by a non-parametric GP; thus,
the control system can also work in unstructured environments.

The main contributions of this study are as follow:

  * •

An approach that uses concurrent learning to estimate system model parameters
in the framework of FBL in a known environment accurately.

  * •

A control system that learns the disturbances and model mismatch in the second
stage, where a non-parametric GP is used to compensate for the disturbances in
unknown, unstructured environments.

  * •

The stability of the nth-order closed-loop system is proven using the Lyapunov
stability theorem.

  * •

Tracking error convergences to zero with and without exact parameter
estimation convergence by using GP.

This paper is organized as follows: The problem statement is described in
Section II. Concurrent learning and Gaussian process are given in Section III.
In Section IV, the details of the developed control algorithm is explained,
the system stability for an n-th order system is proven. Simulations results
have verified the developed control approach leading to minimized tracking
error in Section V. Finally, the study is concluded in Section VI.

##  II PROBLEM STATEMENT

Single-input single-output (SISO) nonlinear system dynamics that can be
linearly parameterized are given by:

$\mathbf{\dot{x}}=\mathbf{Ax}+\mathbf{b}\big{(}u+\mathbf{w}^{*T}{\mathbf{\bm{\phi}(x)}}+d\big{)}$
(1)

where $\mathbf{A}\in\mathbb{R}^{n\times n}$ exploits the integrator chain
structure, $\mathbf{b}=[0,0,...,1]^{T}\in\mathbb{R}^{n}$, $u\in\mathbb{R}$ is
the control input, $\mathbf{x}\in\mathbb{R}^{n}$ is the state vector, the
system dynamics are given by $\mathbf{w}^{*T}\mathbf{\bm{\phi}(x)}$, and $d$
is an unknown disturbance. The system model regressor vectors
$\bm{\phi}(\mathbf{x})\in\mathbb{R}^{m}$ are known and
$\mathbf{w}^{*}\in\mathbb{R}^{m}$ is the unknown ideal weight vector.
Expanding the system equations (1), the state-derivative is given by the last
equation:

$\dot{x}_{n}={\mathbf{w}^{*T}\bm{\phi}(\mathbf{x})}+u+d$ (2)

where $\dot{x}_{n}$ is the state-derivative and $n$ is the order of the
system.

In addition to being linearly parameterizable, the system considered in this
paper must have no internal dynamics, and the state-derivative in the system
model must be measurable or estimated. Moreover, it is assumed there exists a
training environment where the dynamics of the system can be accurately
modelled in the regressor vector $\bm{\phi}(\mathbf{x})$ before the system is
operated in an environment where disturbance $d$ may exist. This does not
necessarily include the weightings of the system dynamics, $\mathbf{w}^{*}$,
but only the dynamics representing the system.

Feedback linearization algebraically transforms a non-linear system (fully or
partially) into a linear one, allowing the use of linear control techniques
[13, 14]. This is achieved by exact state transformations and feedback control
law. Feedback linearization requires a mathematical model of the system and
considers a fixed model estimate traditionally [15]. The larger the error in
the model estimate, the greater the tracking error. Let the estimated weight
vector be given by vector $\mathbf{w}$. To linearize a nonlinear system using
feedback linearization, the control input is given by:

$u=-\mathbf{w}^{T}{\bm{\phi}}(\mathbf{x})+v$ (3)

where $\mathbf{w}$ is the estimate of the model weights and $v$ is a selected
linear controller. This results in the following system:

$\dot{x}_{n}=-\tilde{\mathbf{w}}\bm{\phi}(\mathbf{x})+d+v$ (4)

where $\mathbf{\tilde{w}=w-w^{*}}$ is the modelling error. Note that if
$\mathbf{w\rightarrow w^{*}}$, then $\mathbf{\tilde{w}}=0$ and also if $d=0$,
then a linear system is obtained as

$\dot{x}_{n}=v$ (5)

Considering the closed-loop dynamics presented by (5), the linear controller
$v$ is chosen such that the closed-loop system is stable. Let the tracking
error be defined as :

$\mathbf{e=x_{ref}-x}$ (6)

where $\mathbf{e}$ is the vector of the error states and $\mathbf{x_{ref}}$ is
the vector of the state references where

$\mathbf{e}=\begin{bmatrix}e_{1}&e_{2}&\cdots&e_{n}\end{bmatrix}^{T},\enspace
e_{i}=x_{i_{\mathbf{ref}}}-x_{i}$ (7)

The linear controller, i.e., $v$, is designed such that:

$v=\dot{x}_{n_{\mathbf{ref}}}+k_{1}e_{1}+k_{2}e_{2}+\cdots+k_{n}e_{n}$ (8)

Substituting the linear controller $v$ (8) in (5) gives:

$0=\dot{e}_{n}+k_{n}e_{n}+\cdots+k_{2}e_{2}+k_{1}e_{1}$ (9)

The values of the polynomial coefficients $k_{1},k_{2},\cdots,k_{n}\geq 0$ are
such that the system poles are on the left half of the plane, leading to
exponentially stable dynamics [6]. This demonstrates that the control input
given by (3) results in stable closed-loop dynamics when the system is
perfectly modelled and in the absence of any disturbances [16].

However, in the common case that $\mathbf{w\neq w^{*}}$, and in the presence
of a disturbance, (9) becomes:

$\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d=\dot{e}_{n}+k_{n}e_{n}+\cdots+k_{2}e_{2}+k_{1}e_{1}$
(10)

Eq. (10) shows that the tracking error cannot converge to zero due to the
presence of modelling error and disturbance under the traditional feedback
linearization control law proposed in (3). As the control law is deemed
insufficient in this case, it needs to be redesigned to account for system
modelling errors and disturbances.

The goal of the developed control law in this paper is to design the control
input $u$ such that the tracking error is minimized and converges to zero in
finite-time in the presence of modelling error and disturbance. This goal is
to be achieved using concurrent learning to accurately estimate the system
parameters to minimize the modelling error and Gaussian Process to learn the
disturbance and cancel its effects by compensating in the control input.

##  III CONCURRENT LEARNING AND GAUSSIAN PROCESS

###  III-A Concurrent Learning

Traditionally, an inaccurate model parameter estimate in the feedback
linearization control method leads to an error in the linearization process,
and it is not easy to obtain the model parameter estimate accurately. It would
be ideal to update the model parameter weights using state feedback until the
model parameters converge to their true values. Concurrent learning allows for
model parameter estimation without requiring persistent excitation instead of
providing a user-selected model parameter estimate. It uses the current error
and a set of recorded data to converge parameter estimates to their true
values, as visible in the second half of (11).

The concurrent learning method is initially introduced in the context of model
reference adaptive control. The same concept is integrated with the feedback
linearization framework in this paper. The following weight update law is
used:

$\dot{\mathbf{w}}=-\Gamma_{W}\bm{\phi}(\mathbf{x})\mathbf{e}^{T}\mathbf{Pb}-\sum_{j=1}^{p}\Gamma_{W}\bm{\phi}(\mathbf{x})\epsilon_{j}$
(11)

with the conditions and parameters detailed in [11]. The algorithm to store
data points is given in [12]. The Lyapunov equation is written as follows:

$\mathbf{(A-bk)}^{T}\mathbf{P}+\mathbf{P(A-bk)}+\mathbf{Q}=0$ (12)

where $\mathbf{(A-bk)}$ is introduced in Section IV-C,
$\mathbf{Q}\in\mathbb{R}^{n\times n}$ is any positive definite matrix and
$\mathbf{P}\in\mathbb{R}^{n\times n}$ is the positive definite solution to the
Lyapunov equation [17]. It was shown the weight update law to be exponentially
stable within the model reference adaptive control context [11]. The error
dynamics of the proposed closed-loop system in the context of feedback
linearization is formulated in Section IV-C and is shown to be in the form as
in [11], meaning the stability results apply directly to the proposed closed-
loop system as well. The stability results for the parameter convergence are
not repeated here and can be found in [11, 18].

###  III-B Gaussian Process

Gaussian Processes are used in the control law to compensate for disturbances
when systems exposed to an unknown environment [19, 10]. The following part of
this subsection outlines the key concepts used in the controller design in
this paper.

Gaussian Process (GP) Regression allows for a nonlinear mapping from a set of
input values $X=${$x_{1},x_{2},\cdots,x_{n}$} to output values
$Y=${$y_{1},\cdots,y_{n}$} [20]. The GP is provided with a set of training
data $D=${$X,Y$}, from which a nonlinear map is generated, so an unseen input
$X$ returns a predicted output $Y$ [21]. This makes GP a powerful prediction
tool and useful addition to the control law discussed in Section IV-B when
dealing with disturbances. The idea of Bayesian inference is used to perform
regression on a set of training data, in which the current hypothesis is
updated based on new information. If the greater number of training data is
provided around a point to be estimated, the better GP prediction is obtained.

Two priors are required to define the Gaussian distribution along with the
training data: a mean function $\mu$, which describes the expected value of
the distribution, generally set to zero, and a covariance matrix or kernel,
$\sum$, which provides a measure of similarity between two points and their
values, describing the shape of the distribution. These two priors can
significantly impact the nonlinear mapping generated by the GP and must be
selected carefully. Each kernel has various hyperparameters that are also
user-selected and significantly impact the overall input-output mapping.
However, these parameters can be optimized by solving the maximum log-
likelihood problem [21].

Mathematically, consider the joint Gaussian distribution of the training
outputs, $Y$ corresponding to inputs $X$, and the test outputs, $Y_{*}$
corresponding to test inputs $X_{*}$ [21]:

$\begin{bmatrix}Y\\\
Y_{*}\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}K(X,X)&K(X,X_{*})\\\
K(X_{*},X)&K(X_{*},X_{*})\end{bmatrix}\right)$ (13)

In the controller, the test inputs, $X_{*}$, are the time-dependent system
states. For $n$ training points and $n_{*}$ test points, $K(X,X_{*})$
represents the $n\times n_{*}$ matrix of covariances evaluated at all pairs of
training and test points. The same applies to the other entries of the
covariance function. The squared exponential covariance function is used in
this paper and is given by [22]:

$k_{SE}(X,X_{*})=\sigma^{2}_{f}\exp{\left(-\frac{(X-X_{*})^{2}}{2l^{2}}\right)}$
(14)

Conditioning the joint prior distribution on the training data, $D$, the mean
and variance at query points $X_{*}$ is given by [21]:

$\mu(X_{*})=K(X_{*},X)K(X,X)^{-1}Y$ (15)

$\sigma^{2}(X_{*})=K(X_{*},X_{*})-K(X_{*},X)K(X,X)^{-1}K(X,X_{*})$

(16)

In the context of this paper, a GP is trained with $n$ data pairs, with the
system states, x as inputs, and the control input mismatch $(u-\hat{u})$, as
the output, representing the torque compensation for the disturbance. This can
be performed using two methods; online or offline learning displayed in Fig.
2.

Figure 2: Gaussian Process. Offline vs Online Gaussian Process learning.

In offline learning, the GP is trained across a set of generic inputs, and the
pre-recorded GP is used to predict the compensation for the disturbance in
simulation. This differs from online learning in the sense that the GP is not
updated during the simulation. In online learning, the GP model is
continuously updated throughout tracking the reference signal. This allows the
GP to be trained to meet the experimental conditions and update according to
any model mismatch during the simulation. Using an online GP requires greater
processing power as the prediction hypothesis is updated throughout the
simulation compared to only obtaining a prediction from the GP.

##  IV METHODOLOGY

###  IV-A Reference Model

In matrix form, the reference system can be represented as:

$\mathbf{\dot{x}_{ref}}=\mathbf{Ax_{ref}}+\mathbf{b}u$ (17)

where $\mathbf{A}\in\mathbb{R}^{n\times n}$ exploits the integrator chain
structure,
$\mathbf{x_{ref}}=[x_{1_{\mathbf{ref}}},x_{2_{\mathbf{ref}}},\cdots,x_{{n-1}_{\mathbf{ref}}}]$
is the reference output,
$\mathbf{b}\in\mathbb{R}^{n},\mathbf{b}=[0,0,...,1]^{T}$ as before and $u$ is
the reference input equal to the reference output,
$u=\dot{x}_{n_{\mathbf{ref}}}$. The reference signal is assumed to be
continuously differentiable.

###  IV-B Control Law

Let the control law be defined as:

$u=u_{FBL}+u_{SFB}+u_{ref}-u_{GP}-u_{rob}$ (18)

where $u_{FBL}$ is the feedback linearization law to linearize the system,
$u_{SFB}$ is the chosen linear controller, $u_{ref}$ is the reference signal
required to cancel the input from the reference model, $u_{GP}$ is the
Gaussian process-based estimate of the disturbance, and $u_{rob}$ is
robustness term with its formulation detailed in Section IV-D to ensure
closed-loop system stability. They are formulated as follows:

$\displaystyle u_{FBL}$ $\displaystyle=$
$\displaystyle-\mathbf{w}^{T}\bm{\phi}(\mathbf{x})$ (19) $\displaystyle
u_{SFB}$ $\displaystyle=$
$\displaystyle\mathbf{ke},\enspace\mathbf{k}=[k_{1},k_{2},\cdots,k_{n}]$ (20)
$\displaystyle u_{ref}$ $\displaystyle=$
$\displaystyle\dot{x}_{n_{\mathbf{ref}}}$ (21) $\displaystyle u_{GP}$
$\displaystyle=$
$\displaystyle\mu(\mathbf{x^{*}})=u(\mathbf{x^{*}})-\hat{u}(\mathbf{x^{*}})$
(22) $\displaystyle u_{rob}$ $\displaystyle=$
$\displaystyle-m\frac{\mathbf{b}^{T}\mathbf{Pe}}{\lVert\mathbf{b}^{T}\mathbf{Pe}\rVert}$
(23)

$u_{GP}$ is determined at a query state $\mathbf{x^{*}}$. The estimate is
given by learning the difference in the actual control input, i.e., $u$, and
the control input estimated from the known system model, i.e., $\hat{u}$. If
the model parameters converge to their true value during the learning stage,
only the disturbance is learnt by the Gaussian Process. The plant model is
given by:

$u=\dot{x}_{n}-\mathbf{w}^{*T}\bm{\phi}(\mathbf{x})-d$ (24)

The estimated model is given by:

$\hat{u}=\dot{x}_{n}-\mathbf{w}^{T}\bm{\phi}(\mathbf{x})$ (25)

Using the same state feedback, $\mathbf{x}$, $u-\hat{u}=-d$ if
$\mathbf{w\rightarrow w^{*}}$. If the parameters do not converge, the GP
learns the disturbance and the model convergence error.

###  IV-C Tracking Error

The error dynamics are given by:

$\dot{\mathbf{e}}=\mathbf{\dot{x}_{ref}}-\dot{\mathbf{x}}$ (26)

Substituting the reference and system model in (26), the error dynamics can be
written considering (17) as:

$\dot{\mathbf{e}}=\mathbf{Ax_{ref}}+\mathbf{b}\dot{x}_{n_{ref}}-\Big{(}\mathbf{A}\mathbf{x}+\mathbf{b}\big{(}u+\mathbf{w}^{*T}\bm{\phi}(\mathbf{x})\big{)}\Big{)}$
(27)

Substituting in for the control law (18) and letting
$\mathbf{\tilde{w}=w-w^{*}}$ be the model parameter estimation error, the
error dynamics can be simplified to:

$\dot{\mathbf{e}}=(\mathbf{A-bk})\mathbf{e}+\mathbf{b}\big{(}\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d+\mu(\mathbf{x^{*}})+u_{rob}\big{)}$
(28)

###  IV-D System Stability

The closed-loop system is shown to be stable using the Lyapunov stability
theorem in the possibility of concurrent learning not learning the system
parameters correctly as well. An approach similar to [8] is used below to show
the closed-loop system is stable. The common quadratic Lyapunov function is
proposed below:

$V=\mathbf{e}^{T}\mathbf{Pe}$ (29)

where $\mathbf{e}$ is the tracking error vector and $\mathbf{P}$ is the
positive definite matrix that solves the algebraic Ricatti equation:

$\mathbf{(A-bk)}^{T}\mathbf{P+P(A-bk)+\tilde{S}}=0$ (30)

where $\mathbf{\tilde{S}=Q+k^{T}Rk>0}$ with $\mathbf{Q,R}>0$. The derivative
of the Lyapunov function is given by:

$\dot{V}=\mathbf{e}^{T}\mathbf{P\dot{e}}$ (31)

The time-derivative of Lyapunov function is equivalently using the algebraic
Ricatti equation and by inserting (28) into the equation above:

$\dot{V}=\mathbf{-e}^{T}\mathbf{\tilde{S}e}+2\mathbf{(b}^{T}\mathbf{Pe)}^{T}\big{(}\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d+\mu(\mathbf{x^{*}})+u_{rob}\big{)}$
(32)

The robustness term $u_{rob}$ is now designed such that the system can be
shown to be stable. Using $u_{rob}$ defined in (23) and substituting into
(32):

$\displaystyle\dot{V}=$
$\displaystyle\mathbf{-e}^{T}\mathbf{\tilde{S}e}+2\mathbf{(b}^{T}\mathbf{Pe)}^{T}$
$\displaystyle\Big{(}\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d+\mu(\mathbf{x^{*}})-m\frac{\mathbf{b}^{T}\mathbf{Pe}}{\lVert\mathbf{b}^{T}\mathbf{Pe}\rVert}\Big{)}$
(33)

Using the Lyapunov theorem, $\dot{V}$ needs to be shown to be strictly
negative for the system to be closed-loop stable. The first term
$\mathbf{-e^{T}\tilde{S}e<0}$ is strictly negative. The second term also needs
to be shown to be strictly negative to show system stability. Let the second
term be represented by $V_{2}$. Using the Cauchy-Schwarz inequality, $V_{2}$
can be rewritten as:

$V_{2}\leq
2\lVert\mathbf{b}^{T}\mathbf{Pe}\rVert\big{(}-m+\lvert\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d+\mu(\mathbf{x^{*}})\rvert\big{)}$
(34)

$V_{2}$ is negative if and only if:

$m>\lvert\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d+\mu(\mathbf{x^{*}})\rvert$
(35)

Therefore, as $V_{1}$ is strictly negative, and $V_{2}$ is negative in the
presence of a disturbance or no disturbance, $\dot{V}=V_{1}+V_{2}<0$ and the
closed-loop system is stable.

###### Remark 1

The term above
$\lvert\mathbf{\tilde{w}}\bm{\phi}(\mathbf{x})-d+\mu(\mathbf{x^{*}})\rvert$ is
equal to $\lvert\dot{x}_{n_{measured}}-\dot{x}_{n}\rvert$.

##  V Simulation

The simulation results for the proposed control law applied to the system
presented in [11] are presented in this section. The following cases are
simulated for the system and referred to below:

  1. Case )

FBL without any learning, with model mismatch $\mathbf{\tilde{w}}\neq 0$ in
the absence of disturbance, $d=0$.

  2. Case )

FBL with only concurrent learning for parameter estimation without model
mismatch $\mathbf{\tilde{w}}=0$ in the absence of disturbance, $d=0$.

  3. Case )

FBL with only concurrent learning for parameter estimation without model
mismatch $\mathbf{\tilde{w}}=0$ in the presence of disturbance, $d\neq 0$.

  4. Case )

FBL with concurrent learning for parameter estimation with model mismatch
$\mathbf{\tilde{w}}\neq 0$ and Gaussian process-based online learning for
disturbance in the presence of disturbance, $d\neq 0$.

  5. Case )

FBL with concurrent learning for parameter estimation without model mismatch
$\mathbf{\tilde{w}}=0$ and Gaussian process-based online learning for
disturbance in the presence of disturbance, $d\neq 0$.

All simulations use gains $k_{1}=k_{2}=20$, $\rho=0.01$ and a reference
tracking of $x_{ref}=A\sin{\omega t}$ with $A=0.5$ and $\omega=1\enspace
rad/s$. All initial conditions are zero, i.e., $\textbf{x}(0)=\textbf{0}$. A
learning rate of $\Gamma_{W}=3$ is used for the concurrent parameter
estimation. For the GP, the squared exponential kernel is used as each new
training point reduces the posterior variance [23]. The hyperparameters are
optimized using the maximum log likelihood function and all GP processing is
handled by using the fitrgp function in MATLAB.

The online GP implementation has a ten-second learning period where the model
is initially learnt, and then the model is updated with every previous 100
data points throughout the simulation. Since the disturbance is unknown to a
control system in practice, it is not feasible to use offline GP learning
unless the same simulation was run under the same circumstances, the model was
trained and then deployed again. Using an online GP allows the learning and
training to happen simultaneously and is used below for all simulations.

###  V-A Systems

The system used for simulations is the same as that presented in [11], and is
given by:

$\ddot{\theta}=u+\sin{\theta}-\lvert\dot{\theta}\rvert\theta+0.5\exp{(\theta\dot{\theta})}$
(36)

Rewriting in the form of (2):

$\ddot{\theta}=\mathbf{w}^{*T}\bm{\phi}(\bm{\theta})+u+d$ (37)

where the ideal weights, $\mathbf{w}^{*}=[1,-1,0.5]^{T}$, the regressor vector
values are
$\bm{\phi}(\bm{\theta})=[\sin{\theta},\lvert\dot{\theta}\rvert\theta,\exp{(\theta\dot{\theta})}]^{T}$,
and $u$ is the control input. For Case a), the weight estimate is given by
$\mathbf{w}=[0.5,-1.3,0.75]^{T}$. The added external disturbance is given by:

$d=\begin{cases}0&0<t<10\\\ \cos{\theta}+\dot{\theta}&10\leq t\leq
30\end{cases}$ (38)

The simulation has three stages. In the first ten seconds, there is no
disturbance and concurrent learning is used to estimate the model parameters.
From 10 to 20 seconds, the disturbance is introduced and data is collected to
train the GP. After 20 seconds, the GP begins compensating for the disturbance
and the GP is continuously updated.

###  V-B Simulation Results

(a) Average tracking error (%).

(b) Reference tracking Case e).

(c) Reference tracking Case d).

(d) Parameter estimation: Weights.

(e) Gaussian estimate of the disturbance.

(f) Output of the robustness term.

Figure 3: Simulation results: (a) The average tracking error for each case
outlined at the beginning of Section V. Case d) and Case e) show the minimized
tracking error when using the proposed control law. (b) Reference tracking
when concurrent learning converges to the exact model parameters
($\mathbf{\tilde{w}}=0$) ( Case e) ). (c) Reference tracking when concurrent
learning does not converge to the true model parameters
($\mathbf{\tilde{w}}\neq 0$) ( Case d) ). (d) Exact parameter convergence
using concurrent learning in Case b) and Case e). (e) Gaussian estimate of the
disturbance when the model is exactly parameterized and when a model error is
present ( Case d) and Case e) ). (f) The robustness value when the model is
exactly parameterized and when a model error is present ( Case d) and Case e)
).

The average tracking error for the cases outlined at the beginning of Section
V is shown in Figure 3(a). The proposed control law leads to the reduced
tracking error in both scenarios, Case d) and Case e), even in the presence of
large continuous disturbances. With the Gaussian Process learning and
predicting the system disturbance, the average tracking errors for Case b) and
Case e) are approximately equal. Both cases would have the same tracking error
given that the GP exactly learns the disturbance. The effect of using a fixed
model estimate versus a concurrent learning model is visible in Case a) and
Case b), respectively, in which using concurrent learning drastically reduces
the error. If a case was to exist where the model parameters were not
estimated accurately using concurrent learning, the GP would learn the model
error as part of the disturbance and provide appropriate compensation, as is
shown by Case d).

The three stages for Case e) and Case (d) are shown in Figure 3(b) and 3(c),
respectively. In the first stage, the model parameters are being learnt in a
controlled environment for the first ten seconds only – this provides
sufficient time for the parameters to converge to their parameter estimates.
In the second stage, a disturbance is introduced, the GP also begins learning
the disturbance as a torque mismatch for the 10 to 20 second period. During
this time, the GP provides no compensation and only learns. After 20 seconds,
the GP begins compensating for the effects of the disturbance and the tracking
error is significantly reduced. The ability of the GP to compensate for large
disturbances is apparent by comparing Case c) and Case e). In Case e), the
Gaussian process captures the disturbance only as the model parameters are
exactly estimated, whereas in Case d), the Gaussian process accounts for the
disturbance and the system model error, as shown in Fig. 3(e). The exact
parameter convergence of concurrent learning in Case b) and Case e) is shown
in Fig. 3(d). The robustness value for Case d) and Case e) is also shown in
Fig. 3(f). As the GP compensates for the control input mismatch, the
robustness term has minimal input.

##  VI Conclusion

This paper has developed the use of concurrent learning for exact model
parameterization in the context of feedback linearization, combined with the
use of Gaussian Processes for learning the disturbance for uncertain nonlinear
systems. A robustness term was developed and added to the control law to
ensure closed-loop system stability, which has been proven by using the
Lyapunov stability theorem. Simulation results have verified the effectiveness
of the control law in achieving minimized tracking error in the presence of
model mismatch and disturbances, and ensured more accurate tracking
performance.

## References

  * [1] T. Westenbroek, D. Fridovich-Keil, E. Mazumdar, S. Arora, V. Prabhu, S. S. Sastry, and C. J. Tomlin, “Feedback linearization for unknown systems via reinforcement learning,” 2020.
  * [2] E. Kayacan, “Closed-loop error learning control for uncertain nonlinear systems with experimental validation on a mobile robot,” _IEEE/ASME Trans. Mechatronics_ , vol. 24, no. 5, pp. 2397–2405, 2019.
  * [3] J. Umlauft and S. Hirche, “Feedback linearization based on gaussian processes with event-triggered online learning,” _IEEE Trans. Autom. Control_ , vol. 65, no. 10, pp. 4154–4169, 2020.
  * [4] L. Oliveira, A. Bento, V. J. Leite, and F. Gomide, “Evolving granular feedback linearization: Design, analysis, and applications,” _Applied Soft Computing_ , vol. 86, p. 105927, 2020.
  * [5] E. Kayacan, “Sliding mode learning control of uncertain nonlinear systems with lyapunov stability analysis,” _Transactions of the Institute of Measurement and Control_ , vol. 41, no. 6, pp. 1750–1760, 2019.
  * [6] J.-J. E. Slotine and W. Li, _Applied nonlinear control_. Prentice-Hall, 1991.
  * [7] J. Chai, E. Medagoda, and E. Kayacan, “Adaptive and efficient model predictive control for booster reentry,” _Journal of Guidance, Control, and Dynamics_ , vol. 43, no. 12, pp. 2372–2382, 2020.
  * [8] M. Greeff and A. P. Schoellig, “Exploiting differential flatness for robust learning-based tracking control using gaussian processes,” _IEEE Control Systems Letters_ , vol. 5, no. 4, pp. 1121–1126, 2021.
  * [9] T. Beckers, J. Umlauft, D. Kulic, and S. Hirche, “Stable gaussian process based tracking control of lagrangian systems,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017, pp. 5180–5185.
  * [10] T. Beckers, D. Kulić, and S. Hirche, “Stable gaussian process based tracking control of euler–lagrange systems,” _Automatica_ , vol. 103, pp. 390–397, 2019.
  * [11] G. Chowdhary and E. Johnson, “Concurrent learning for convergence in adaptive control without persistency of excitation,” in _49th IEEE Conference on Decision and Control (CDC)_ , 2010, pp. 3674–3679.
  * [12] E. Kayacan, S. Park, C. Ratti, and D. Rus, “Online system identification algorithm without persistent excitation for robotic systems: Application to reconfigurable autonomous vessels,” in _IEEE/RSJ Int. Conf. Intell Robots Syst._ , 2019, pp. 1840–1847.
  * [13] L. Martins, C. Cardeira, and P. Oliveira, “Feedback linearization with zero dynamics stabilization for quadrotor control,” _Journal of Intelligent & Robotic Systems_, vol. 101, no. 1, p. 7, 2020.
  * [14] J. Moreno–Valenzuela, J. Montoya–Cháirez, and V. Santibáñez, “Robust trajectory tracking control of an underactuated control moment gyroscope via neural network–based feedback linearization,” _Neurocomputing_ , vol. 403, pp. 314–324, 2020.
  * [15] G. Li, A. Luo, Z. He, F. j. Ma, Y. Chen, W. Wu, Z. Zhu, and J. M. Guerrero, “A dc hybrid active power filter and its nonlinear unified controller using feedback linearization,” _IEEE Trans. Ind. Electron._ , pp. 1–1, 2020.
  * [16] M. Mehndiratta, E. Kayacan, M. Reyhanoglu, and E. Kayacan, “Robust tracking control of aerial robots via a simple learning strategy-based feedback linearization,” _IEEE Access_ , vol. 8, pp. 1653–1669, 2020.
  * [17] D. Boley and B. Datta, “Numerical methods for linear control systems,” 1994.
  * [18] R. Kamalapurkar, B. Reish, G. Chowdhary, and W. E. Dixon, “Concurrent learning for parameter estimation using dynamic state-derivative estimators,” _IEEE Trans. Autom. Control_ , vol. 62, no. 7, pp. 3594–3601, 2017.
  * [19] E. Bradford, L. Imsland, D. Zhang, and E. A. del Rio Chanona, “Stochastic data-driven model predictive control using gaussian processes,” _Computers & Chemical Engineering_, vol. 139, p. 106844, 2020.
  * [20] Y. Cao, J. Huang, H. Ru, W. Chen, and C. H. Xiong, “A visual servo-based predictive control with echo state gaussian process for soft bending actuator,” _IEEE/ASME Trans. Mechatronics_ , vol. 26, no. 1, pp. 574–585, 2021.
  * [21] C. E. Rasmussen and C. K. I. Williams, _Gaussian processes for machine learning_ , ser. Adaptive computation and machine learning. MIT Press, 2006.
  * [22] D. Duvenaud, “Automatic model construction with gaussian processes,” Ph.D. dissertation, 2014.
  * [23] I. Steinwart and A. Christmann, _Support vector machines_ , 1st ed., ser. Information science and statistics. Springer, 2008.