11institutetext: Visual Geometry Group, Department of Engineering Science,
University of Oxford 22institutetext: Plexalis Ltd

# A Convolutional Approach to Vertebrae Detection and Labelling in Whole Spine
MRI

Rhydian Windsor 11 Amir Jamaludin 11 Timor Kadir 22 Andrew Zisserman 11

###### Abstract

We propose a novel convolutional method for the detection and identification
of vertebrae in whole spine MRIs. This involves using a learnt vector field to
group detected vertebrae corners together into individual vertebral bodies and
convolutional image-to-image translation followed by beam search to label
vertebral levels in a self-consistent manner. The method can be applied
without modification to lumbar, cervical and thoracic-only scans across a
range of different MR sequences. The resulting system achieves 98.1% detection
rate and 96.5% identification rate on a challenging clinical dataset of whole
spine scans and matches or exceeds the performance of previous systems of
detecting and labelling vertebrae in lumbar-only scans. Finally, we
demonstrate the clinical applicability of this method, using it for automated
scoliosis detection in both lumbar and whole spine MR scans.

###### Keywords:

Vertebral bodies Whole spine MRI Scoliosis.

##  1 Introduction

The objective of this paper is automated vertebrae detection and
identification of vertebral levels. This is an important task for several
reasons. Firstly, automated diagnosis of many spinal diseases such as disc
degeneration [9, 10] or spinal stenosis [13] relies on accurate localisation
of vertebral structures or, in the case of pathological scoliosis, lordosis
and kyphosis [8], analysing the geometry of the spinal column. Secondly,
vertebral bodies can be used to infer other spinal structures of interest such
as the spinal cord or ribs. Finally, vertebrae can act as points to allow
registration between different scans.

There are several issues that make this task challenging. One of the most
obvious is that vertebrae are highly repetitive and hence distinguishing
between different levels can be hard. Labelling by simply counting down from
the C2 vertebra is problematic as it assumes that all vertebrae have been
detected, C2 is visible, and that every patient has the same number of
vertebrae which is not always true. Furthermore, for clinical use, labelling
must be robust to: variations in spinal anatomy (such as collapsed vertebrae,
hemivertebrae and fused vertebra); vertebrae numbers – around 11.3% of the
population have one more or one less mobile vertebra [18]); different imaging
parameters including MR weighting (e.g. T1, T2, STIR, TIRM and FLAIR); fields
of view (e.g. lumbar, whole spine scans); scan resolution and also
number/thickness of slices in the scan.

This paper proposes a new approach to this challenge, in particular in the
case of 3-D sagittal whole spine clinical MRIs which are important for
diagnosing several diseases such as ankylosing spondylitis and multiple
myeloma. We make the following contributions to the tasks of vertebral body
detection and vertebral level labelling in clinical MRIs: We propose a new
convolutional method of detection based on localising the corners and
centroids of vertebrae and then grouping them together (Section 2); We
reformulate the labelling task as a convolution enhancement followed by a
language modelling inspired sequential correction, removing the need for a
recurrent network and showing robustness to variations in vertebra numbers
(Section 3); We show that the resulting system is robust to a variety of
fields of view and pathologies by evaluating on a large clinical dataset of MR
lumbar scans, and also, for the first time, on whole spine scans. We achieve
state-of-the-art performance at vertebra identification in both (Section 4);
We demonstrate a clinical application of this system by using it to
automatically detect cases of scoliosis (Section 5).

Related work: There have been several approaches to automated detection and
labelling of vertebrae in MRIs although most focus on fixed fields of view,
e.g. lumbar or cervical scans only [3, 9, 12, 13]. Previous vertebrae
labelling methods tend to rely on either heuristic-based graphical models [3,
11] or assuming the bottom vertebra is S1 and counting up [13]. Zhao et al.
[21] perform labelling in MRIs with arbitrary fields of view, but only in the
lower spine (from S1-T12 to L4-T10). Windsor and Jamaludin [19] report a
method of detecting, though not labelling, vertebrae in full spine scans
iteratively but also require the location of the S1 vertebra for
initialization. Cai et al. [2] also perform detection and labelling by a
hierarchical deformation model and even report success in a single full spine
scan but do not evaluate the performance quantitatively on a whole spine
dataset. Furthermore such models are slow to apply and make strong assumptions
on spinal geometry such as a fixed number of vertebrae and lack of major
pathology (fused/collapsed vertebrae). Greater progress has been made in
whole-body CT images where the task is more straightforward; CT imaging
protocols are highly standardised due to image intensities being consistent,
representing X-Ray absorptions of voxel locations. Also, CT scans tend to have
more 3-D information, with higher resolution in the sagittal plane than that
typical of clinical MR. Approaches using graphical models and recurrent neural
networks for labelling have been reported to achieve 70-89% identification
rate in mixed field of view CT scans [4, 5, 20].

##  2 Detecting Vertebrae

Detection of the vertebrae proceeds in three stages, as illustrated in Figure
1. First, corners and centroids of each vertebra are predicted in each
sagittal slice. Second, each detected corner is assigned to a centroid by
predicting a vector field for each corner type (e.g. top left, bottom right)
which points to the corresponding centroid. The detected centroid and the four
corners which point to it define the bounding quadrilateral for that vertebra
in that slice. Third, these quadrilaterals are grouped across slices to define
detected volumes for each vertebra.

Inference is performed here by a U-Net architecture (shown in detail in the
appendix) which ingests an MRI (one slice) and outputs 13 channels: 5 channels
are the landmark detection channels and each is used to detect the centroid
and 4 corners respectively of all vertebral bodies appearing in that image;
the remaining 8 channels are the landmark grouping channels corresponding to
the $x$ and $y$ coordinates of the 4 vector fields used in grouping.
In more detail, the locations of corners and centroids are identified as modes
of the heatmaps in the landmark detection channels. For each corner detected,
there exists a corresponding grouping vector field from the landmark grouping
channels; 4 corners equals 4 vectors fields. A group is then formed from the 4
closest corners pointing to a single centroid; a group here forms a
quadrilateral. This is done for each detected centroid. If two or more
centroids are assigned the same corner, the centroid to which the corner’s
vector points closest remains, and the other centroid is discarded to stop
double detection of a single vertebra. If there is no detected corner within a
fixed range of a centroid, it is also discarded, making the system robust to
spurious centroid detections. This is performed in each sagittal slice.
Finally, the vertebral bodies detected are grouped across slices by measuring
the IoU between quadrilaterals and assigning them to the same vertebra if they
have an overlap of greater than 0.5.

Discussion. Our approach of detecting vertebrae as a series of points and then
grouping them differs from other methods which have used region proposal
networks [21] pixel/voxel-wise segmentation [23], deformable-part models [2]
or simply detected vertebrae centroids [3, 20]. However, this approach is
receiving increasing attention in the computer vision literature, with Zhou et
al. [22] showing state-of-the-art results by detecting objects as a series of
keypoints. The advantage of this method is its high speed (1 second inference
on GPU-enabled hardware) combined with more accurate bounding regions than
using standard bounding boxes.

Figure 1: The pipeline used to detect vertebrae in whole spine scans. The
output detection channels are thresholded and each resulting connected volume
becomes a corner or a centroid. Each corner has corresponding vector in the
grouping channels at the point-of-detection and ‘points’ with a magnitude and
direction according to that vector. Each centroid is assigned the 4 corners
that point closest to it.

Training: To train the network, we use separate loss functions for the
detection and grouping channels. The ground truth annotations are coordinates
of the corners of the vertebral bodies (centroids are computed from these).
Target detection channel outputs are constructed by overlaying a Gaussian
kernel on each annotated ground truth landmark in the detection channels with
variance proportional to the square root of the area of the detection. These
channels are trained using the weighted L1 loss;
$\mathcal{L}_{detect}(Y,\hat{Y})=\sum_{k=1}^{5}\sum_{i,j}\alpha_{ijk}|y_{ijk}-\hat{y}_{ijk}|$
where $Y$ is the response map output by the network and $\hat{Y}$ is the
target response map. $y_{ijk}$ and $\hat{y}_{ijk}$ are the value of response
and ground truth map respectively at image coordinate $(i,j)$ in detection
channel $k$ and $a_{ijk}$ is a weighing factor given by
$\alpha_{ijk}=\frac{N_{k}}{N_{k}+P_{k}}\mbox{if }\hat{y}_{ijk}\geq T$ or
$\frac{P_{k}}{N_{k}+P_{k}}\mbox{if }\hat{y}_{ijk}<T$, where $P_{k}$ and
$N_{k}$ are the number of pixels respectively above and below threshold $T$ in
channel $k$. This weighing factor balances the loss from false positive
responses in the heatmap and false negatives, speeding up training. In the
experiments in this paper, $T=0.01$ was used. Landmark grouping channels are
trained by using the L2 loss;
$\mathcal{L}_{group}=\sum_{l=1}^{4}\sum_{b}\sum_{(i,j)\in\mathcal{N}_{bl}}||\mathbf{v}^{l}_{ij}-\mathbf{r}^{k}_{ij}||_{2}^{2}.$
Here $l$ indexes each corner type/vector field (e.g. top left, bottom right),
$b$ indexes labelled ground truth vertebral body, $\mathcal{N}_{bl}$ is a
neighbourhood of pixels surrounding the $l^{th}$ corner of vertebral body $b$.
$\mathbf{v}^{l}_{ij}$ is the value of the vector field corresponding to corner
$l$ at the pixel location $(i,j)$ and $\mathbf{r}^{b}_{ij}$ is the
displacement from the centroid of vertebral body $b$ to $(i,j)$. As suggested
by [15], heavy augmentation is used during training. Scans are padded,
rotated, zoomed and flipped in the coronal plane. Non-square scans are split
into squares overlapping by 40% and resized to $224\times 224$ to ensure
constant sized input to the network.

##  3 Labelling Vertebrae

Once the vertebrae are detected, the next task is to label each vertebrae with
its level (e.g. S1, L5, L4 etc.). There are two types of information to
consider when labelling a vertebra; the appearance of the vertebrae – its
intensity pattern, shape and size – and its context – its position in relation
to other vertebrae that have been detected in the scan. As such, we train two
networks; an appearance network to infer the level of vertebra from its
appearance alone, and a context network which takes as input the predictions
of the appearance network along with the spatial configuration of the
detections in the scan to improve the predictions. The final stage is to
search for a consistent labelling of the vertebra using a sequential ‘language
model’ that builds in ordering constraints. The labelling pipeline is outlined
in Figure 2. It should be noted that both networks are fully convolutional.
This differs from the approaches outlined in [20] and [12] which use recurrent
neural networks and graphical models respectively.

Figure 2: Overview the labelling pipeline. Volumes surrounding each of the
vertebrae detected are extracted and input to an appearance network. This is
then used to construct an input image to the context network which gives the
probability of a given label at a given height in the image. The resulting
probability map is shown in the far right of the image, followed by a beam
search to generate the final level sequence. The architectures of both
networks are shown in the appendix.

The labelling pipeline proceeds as follows: A volume around each detected
vertebra is extracted and given as input to the appearance model which
predicts a softmax probability vector over the labels. Input volumes are
created by first fitting a bounding box around the detection, then expanding
the box by 100% to include nearby anatomical features and resampling its size
to $224\times 224\times 16$, where 16 is the number of sagittal slices. For
the appearance network, we use a simple 3D CNN (details of this and the
context network are given the appendix). Then, the predictions of the
appearance networks are used to make a probability-height map where the
probability vector at the height of each vertebra is equal to the output
vector of the appearance network for that vertebra. In detail: the
probability-height map, $P\in\mathbb{R}^{H\times N}$ where $H$ is the height
of the image and $N=24$ is the number of vertebrae level classes (C2 to S1).
If vertebra detection $v$ spans from heights $h_{v_{1}}$ to $h_{v_{2}}$ in the
scan then $P_{h^{\prime}n}=p^{a}_{vn}\forall\;h_{v_{1}}\leq h^{\prime}\leq
h_{v_{2}}$ and 0 otherwise where $p^{a}_{vn}$ is the probability $v$ has level
$n$, given by the appearance network. Next, the context network refines this
probability-height map, using detections around each vertebra to update its
predicted class. Finally, the output map is decoded into a logical sequence of
vertebra using a language-model inspired beam search which imposes global
constraints such as no repetitions. This is described in detail in Section
3.1.

Training: The appearance network is trained using a cross-entropy loss
function. For re-calibration [7], a softmax operation with temperature $T=10$
is applied to the logits layer. $P$ is given to the context net which produces
$P^{\prime}$, a refined version of the same map as shown on the far right of
Figure 2. The cross-entropy between the output probability map at the centroid
height of each vertebra and its ground truth label is used as the loss
function. When training the context network, augmentation is applied at
training time by randomly removing between 0 to 4 of the highest and lowest
detections from the input image and further removing each remaining vertebra’s
probability map from the input image with probability 0.2. This ensures the
labelling pipeline is robust to missed detections and not reliant on distinct
looking vertebrae at the top and bottom of the spine.

###  3.1 Enforcing Monotonicity and Dealing with Numerical Variations

To predict the labels for each vertebrae detected, $P^{\prime}$ shown in
Figure 2 must be decoded into a sequence of level predictions. A naive method
to do this would be to take the argmax of the probability map at each vertebra
centroid height and use this as a label. However, this allows implausible
sequences such as repeated vertebral levels to be predicted. If the detection
$v$ at height $h_{v}$ has label index $l$, then we know $v^{\prime}$ at
$h_{v^{\prime}}>h_{v}$ detection should have label $l^{\prime}>l$ with the
greatest probability of being $l+1$(e.g. L4 should be followed by L3). The
challenge of imposing such constraints is analogous to that faced in automatic
speech recognition where CTC training is followed by language modelling to
predict a valid character or word sequence [16, 6]. We take inspiration from
this and use a beam search to obtain a valid labelling: beginning at the
highest detection and searching down, sequences are generated and scored by
selecting the highest probability labels for each vertebra. The $k$ most
likely sequences of levels are stored in memory at each step, where $k$ is the
beam width. Sequences with repetitions of levels are given probability 0 and
those with skipped levels are penalised by multiplying the sequence
probability by a penalty score. The method can also incorporate numerical
variations, with $\pm 1$ lumbar vertebrae given a small sequence probability
penalty reflecting the low incidence of this.

##  4 Detection and Labelling Results

We evaluate the system at the task of vertebrae detection and labelling in
whole spine and lumbar scans. Following [11, 20], we define correct detections
to be when the ground truth vertebra centroid is contained entirely within a
single bounding quadrilateral. For detection, we report precision, recall, and
the localisation error (LE), defined to be the mean distance of ground truth
centroids from the closest detection quadrilateral centroid. For labelling, we
report identification rate (IDR), the fraction of vertebrae detected and
labelled correctly.

###  4.1 Datasets

Three datasets are used in this work: OWS, Genodisc and Zukić. OWS is a
dataset of 710 sagittal whole spine scans across 196 patients from the Picture
Archiving and Communication System (PACS) of an orthopaedic centre. The
dataset exhibits a wide range of pathologies such as hemivertebra, fused
vertebrae, numerical variations of vertebrae and scoliosis. Scans are taken
from different scanners with a range of MR parameters (T1, T2, FLAIR, TIRM and
STIR). The dataset is split into training, validation and testing sets with an
60/20/20% split at the patient level. Corners of vertebrae from S1 to C2 are
annotated and vertebrae levels were marked by a radiologist in one scan for
each patient, with S1 being the first vertebra attached to the pelvis. In the
case of 25 vertebrae between S1 and C2 instead of the normal 24, an extra
lumbar vertebrae is labelled (L6). Networks were trained on the OWS training
set, using an Adam optimizer with $\beta_{1}=0.9$, $\beta_{2}=0.999$ and a
learning rate of 0.001. The Genodisc and Zukić datasets are used only for
testing. Genodisc’s test set has 421 clinical lumbar MRIs used by Lootus et
al. [12]. Zukić [23] is a small dataset of 17 mostly lumbar sagittal MRIs
available on the online SpineWeb platform.

###  4.2 Results

The results of detection and labelling on all datasets are shown in Table 1
with comparisons to other methods reported on the same datasets where
available. We also compare our convolutional labelling pipeline results to a
baseline recurrent approach for vertebra labelling, training a bidirectional
LSTM on the appearance features extracted from each detected volume. Example
predicted detection and labelling sequences across a range of pathologies are
given in Figure 3. The LSTM baseline used is detailed in the appendix.

Dataset | No. Scans | No. Vert | Method | Prec. (%) | Rec. (%) | IDR(%) | IDR$\pm$ 1(%) | LE (mm)
---|---|---|---|---|---|---|---|---
OWS (Whole Spine) |  |  | Windsor† [19] | 99.4 | 99.4 | - | - | 1.0 $\pm$ 0.9
37 | 888 | Label Baseline | - | - | 86.9 | 93.4 | -
|  | Ours | 99.0 | 98.1 | 96.5 | 97.3 | 2.4 $\pm$ 1.3
Genodisc (Lumbar) |  |  | Lootus [11] | - | - | 86.9 | - | 3.5 $\pm$ 3.3
421 | 2947 | Label Baseline | - | - | 90.1 | 97.4 | -
|  | Ours | 99.7 | 99.7 | 98.4 | 99.7 | 1.6 $\pm$ 1.1
Zukić (Lumbar) |  |  | Zukić [23] | 98.7 | 92.9 | - | - | 1.6 $\pm$ 0.8
17 | 154 | Label Baseline | - | - | 87.0 | 94.3 | -
|  | Ours | 99.3 | 98.7 | 90.9 | 98.7 | 2.0 $\pm$ 1.5
Table 1: Performance of the pipeline on the three datasets. Our approach is
compared with other methods using the same datasets and also a LSTM labelling
baseline. Results are reported on a per-vertebra level. Higher is better for
detection precision (Prec.), detection recall (Rec.), and level identification
rate (IDR). Lower is better for localisation error (LE). We also report the
percentage of vertebrae within one level of their ground truth value (IDR$\pm
1$). Lootus [11] is tested on a subset of 291 scans from the Genodisc
dataset.Note, Windsor† [19] requires manual initialization by providing the
location of the S1 vertebra, so is not directly comparable Figure 3: Example
detection and labelling of vertebrae for a range of whole spine scans. a) and
b) are examples of typical spines. The other images shows examples of the
system dealing with pathologies; c) shows an example of a spine with fused C3
and C4 vertebrae; d) is an example of a spine with a collapsed vertebra; e)
shows a spine with an extra vertebra between S1 and C2 and f) shows an
extremely kyphotic spine with a hemivertebra at the bottom of the thoracic
spine.

Discussion: For whole spine detection on the OWS dataset, the proposed method
achieves a high precision and recall of 99.0% and 98.1% respectively. It
achieves a level identification rate of 96.5%, significantly exceeding the
LSTM baseline. The few errors the system made for labelling are generally due
to S2 being detected as S1, meaning all labels are out by one. In practice,
this is a mistake radiologists often make as it can be difficult to tell which
is the first sacral vertebra without looking at axial scans to see which bone
is joint to the pelvis. This also explains why [19] achieves slightly higher
precision and recall at detection than ours; it is a semi-automated algorithm
given the location of S1 at initialization and thus bypasses this difficult
problem of S1 recognition. On the Genodisc lumbar spine dataset, our method
again significantly outperforms the baseline, and it is also outperforms the
prior method of Lootus et al. [11] (98.4% compared to 86.9%). Importantly, OWS
and Genodisc are scans of patients with a wide range of pathologies imaged
using typical clinical MR protocols; hence strong performance here gives
evidence of the clinical usefulness of this approach. We show example results
for pathological spines in the appendix. In the appendix we also report
results by other groups for vertebrae detection and labelling in MRIs.
However, these are for different datasets to which we could not get access,
with different scanning protocols, fields of view (FoV) and patient sets, and
thus cannot be compared to directly.

##  5 Automated Scoliosis Detection

Finally, as an illustration of a potential application of the proposed
approach, we explore the ability of the system to classify cases of scoliosis
from sagittal MR scans. In clinical practice this is usually determined by
measuring the Cobb angle in coronal views of X-ray scans [1], however
measuring scoliosis in the supine position has also been shown to be possible
[17, 8] and MRI can be useful for understanding disease etiology and symptoms
[14]. This is a more difficult task in sagittal scans as it requires sensitive
detection of the sides of vertebra with a clear decision boundary as to when a
vertebra is present or not present in a slice which can be difficult in cases
of partial visibility. In the entirety of the Genodisc dataset, scoliosis was
reported by a radiologist in 198 of 3542 scans, across 2009 patients. By
measuring statistics of a quintic polynomial fit through vertebra centroids
and using them as predictive features we develop classifiers for this label.
Specifically, we measure the maximum curvature of the polynomial, and the
maximum deviation of the curve from a straight vertical line fit through the
vertebrae, assuming that a low curvature vertebral column with little
deviation from the centreline corresponds to a non-scoliotic spine. While
ground truth scoliosis labels were not available for the whole spine scans, we
give qualitative results, comparing the features of curves fit though
vertebrae of a scoliotic and non-scoliotic scan from Zukić. Results of these
experiments are shown in Figure 4.

Figure 4: Results of scoliosis classification in sagittal scans using the
proposed vertebrae detection system: (a) ROC curve for simple scoliosis
classifiers based on statistics of polynomial curves fitted through vertebrae
centroids in Genodisc scans; (b) A qualitative comparison of curves fit
through detected vertebrae in full spine scans. Curves are overlaid on coronal
slices synthesised from the sagittal slices.

Discussion: The results for automated scoliosis detection from sagittal scans
are promising. Using simple classifiers AUCs of 0.636-0.689 are achieved in a
highly class-imbalanced problem. Of the features measured, distance from the
vertical centreline of the vertebrae performed best. The system is also shown
to capture scoliotic curves in full spine scans too. These experiments
illustrate that the proposed method produces a strong geometric representation
of the spine, which can be used in further downstream tasks.

##  6 Conclusion

We introduce a novel method for vertebrae detection and labelling in whole
spine sagittal MRIs. It shows state-of-the-art results for vertebra
identification in lumbar scans, with little performance drop in the whole
spine case, and is robust to a range of spinal defects, numerical variations
of vertebrae and different scanning protocols. We also demonstrate a potential
diagnostic application: automated detection of scoliosis from sagittal MRIs.
Future work will include integrating automated detection of different spinal
pathologies and implementing the system for CT scans.

Acknowledgements. The authors would like to thank Dr. Sarim Ather for useful
discussions on spinal anatomy and clinical approaches to diagnosing disease,
as well as assistance labelling the data. Rhydian Windsor is supported by
Cancer Research UK as part of the EPSRC CDT in Autonomous Intelligent Machines
and Systems (EP/L015897/1). Amir Jamaludin is supported by EPSRC Programme
Grant Seebibyte (EP/M013774/1). The Genodisc data was obtained during the EC
FP7 project GENODISC (HEALTH-F2-2008-201626).

## References

  * [1] Aebi, M.: The adult scoliosis. European Spine Journal 14(10), 925–948 (2005)
  * [2] Cai, Y., Osman, S., Sharma, M., Landis, M., Li, S.: Multi-Modality Vertebra Recognition in Arbitrary Views Using 3D Deformable Hierarchical Model. IEEE Transactions on Medical Imaging 34(8), 1676–1693 (2015)
  * [3] Forsberg, D., Sjöblom, E., Sunshine, J.L.: Detection and Labeling of Vertebrae in MR Images Using Deep Learning with Clinical Annotations as Training Data. Journal of Digital Imaging 30(4), 406–412 (2017)
  * [4] Glocker, B., Feulner, J., Criminisi, A., Haynor, D.R., Konukoglu, E.: Automatic Localization and Identification of Vertebrae in Arbitrary Field-of-View CT Scans. In: Medical Image Computing and Computer-Assisted Intervention (2012)
  * [5] Glocker, B., Zikic, D., Konukoglu, E., Haynor, D.R., Criminisi, A.: Vertebrae Localization in Pathological Spine CT via Dense Classification from Sparse Annotations. In: Medical Image Computing and Computer-Assisted Intervention. pp. 262–270 (2013)
  * [6] Graves, A., Fernández, S., Gomez, F., Schmidhuber, J.: Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks. In: International Conference on Machine learning. ICML ’06 (2006)
  * [7] Guo, C., Pleiss, G., Sun, Y., Weinberger, K.Q.: On Calibration of Modern Neural Networks. In: International Conference on Machine Learning (2017)
  * [8] Jamaludin, A., Kadir, T., Clark, E., Zisserman, A.: Predicting spine geometry and scoliosis from DXA scans. In: MICCAI Workshop: Computational Methods and Clinical Applications in Musculoskeletal Imaging (2019)
  * [9] Jamaludin, A., Kadir, T., Zisserman, A.: SpineNet: Automated classification and evidence visualization in spinal MRIs. Medical Image Analysis 41, 63–73 (2017)
  * [10] Jamaludin, A., Lootus, M., Kadir, T., Zisserman, A., Urban, J., Battié, M.C., Fairbank, J., McCall, I.: Automation of reading of radiological features from magnetic resonance images (mris) of the lumbar spine without human intervention is comparable with an expert radiologist. European Spine Journal (2017)
  * [11] Lootus, M., Kadir, T., Zisserman, A.: Vertebrae detection and labelling in lumbar MR images. In: MICCAI Workshop: Computational Methods and Clinical Applications for Spine Imaging (2013)
  * [12] Lootus, M., Kadir, T., Zisserman, A.: Radiological grading of spinal MRI. In: MICCAI Workshop: Computational Methods and Clinical Applications for Spine Imaging (2014)
  * [13] Lu, J.T., Pedemonte, S., Bizzo, B., Doyle, S., Andriole, K.P., Michalski, M.H., Gonzalez, R.G., Pomerantz, S.R.: Deep spine: Automated lumbar vertebral segmentation, disc-level designation, and spinal stenosis grading using deep learning. In: Machine Learning for Healthcare Conference (2018)
  * [14] Ozturk, C., Karadereler, S., Ornek, I., Enercan, M., Ganiyusufoglu, K., Hamzaoglu, A.: The role of routine magnetic resonance imaging in the preoperative evaluation of adolescent idiopathic scoliosis. International Orthopaedics 34(4), 543–546 (Apr 2010)
  * [15] Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional Networks for Biomedical Image Segmentation. In: Medical Image Computing and Computer-Assisted Intervention (2015)
  * [16] Scheidl, H., Fiel, S., Sablatnig, R.: Word Beam Search: A Connectionist Temporal Classification Decoding Algorithm. In: 2018 16th International Conference on Frontiers in Handwriting Recognition (ICFHR) (2018)
  * [17] Taylor, H.J., Harding, I., Hutchinson, J., Nelson, I., Blom, A., Tobias, J.H., Clark, E.M.: Identifying scoliosis in population-based cohorts: development and validation of a novel method based on total-body dual-energy x-ray absorptiometric scans. Calcified Tissue International 92(6), 539–547 (Jun 2013)
  * [18] Tins, B.J., Balain, B.: Incidence of numerical variants and transitional lumbosacral vertebrae on whole-spine MRI. Insights into Imaging 7(2), 199–203 (2016)
  * [19] Windsor, R., Jamaludin, A.: The ladder algorithm: Finding repetitive structures in medical images by induction. In: IEEE International Symposium on Biomedical Imaging (2020)
  * [20] Yang, D., Xiong, T., Xu, D., Huang, Q., Liu, D., Zhou, S.K., Xu, Z., Park, J., Chen, M., Tran, T.D., Chin, S.P., Metaxas, D., Comaniciu, D.: Automatic Vertebra Labeling in Large-Scale 3D CT Using Deep Image-to-Image Network with Message Passing and Sparsity Regularization. In: Information Processing in Medical Imaging (2017)
  * [21] Zhao, S., Wu, X., Chen, B., Li, S.: Automatic Vertebrae Recognition from Arbitrary Spine MRI Images by a Hierarchical Self-calibration Detection Framework. In: Medical Image Computing and Computer Assisted Intervention (2019)
  * [22] Zhou, X., Wang, D., Krähenbühl, P.: Objects as points. In: arXiv preprint arXiv:1904.07850 (2019)
  * [23] Zukić, D., Vlasák, A., Egger, J., Hořínek, D., Nimsky, C., Kolb, A.: Robust Detection and Segmentation for Diagnosis of Vertebral Diseases Using Routine MR Images. Computer Graphics Forum 33(6), 190–204 (2014)


# Data-driven Self-triggered Control via Trajectory Prediction

Wenjie Liu, Jian Sun, , Gang Wang, ,
Francesco Bullo, , and Jie Chen The work was supported in part by the National
Key R&D Program of China under Grant 2021YFB1714800, in part by the National
Natural Science Foundation of China under Grants 62088101, 62173034, 61925303,
and in part by the Chongqing Natural Science Foundation under Grant
2021ZX4100027. (_Corresponding author_ : Gang Wang.) W. Liu, J. Sun, and G.
Wang are with the State Key Lab of Intelligent Control and Decision of Complex
Systems and the School of Automation, Beijing Institute of Technology, Beijing
100081, China, and also with the Beijing Institute of Technology Chongqing
Innovation Center, Chongqing 401120, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>gangwang@bit.edu.cn). F. Bullo is with the Mechanical
Engineering Department and the Center of Control, Dynamical Systems and
Computation, UC Santa Barbara, CA 93106-5070, USA (e-mail: bullo@ucsb.edu). J.
Chen is with the Department of Control Science and Engineering, Tongji
University, Shanghai 201804, China, and also with the State Key Lab of
Intelligent Control and Decision of Complex Systems, Beijing Institute of
Technology, Beijing 100081, China (e-mail: chenjie@bit.edu.cn).

###### Abstract

Self-triggered control, a well-documented technique for reducing the
communication overhead while ensuring desired system performance, is gaining
increasing popularity. However, existing methods for self-triggered control
require explicit system models that are assumed perfectly known _a priori_. An
end-to-end control paradigm known as data-driven control learns control laws
directly from data, and offers a competing alternative to the routine system
identification-then-control method. In this context, the present paper puts
forth data-driven self-triggered control schemes for unknown linear systems
using data collected offline. Specifically, for output feedback control
systems, a data-driven model predictive control (MPC) scheme is proposed,
which computes a sequence of control inputs while generating a predicted
system trajectory. A data-driven self-triggering law is designed using the
predicted trajectory, to determine the next triggering time once a new
measurement becomes available. For state feedback control systems, instead of
capitalizing on MPC to predict the trajectory, a data-fitting problem using
the pre-collected input-state data is solved, whose solution is employed to
construct the self-triggering mechanism. Both feasibility and stability are
established for the proposed self-triggered controllers, which are validated
using numerical examples.

###### Index Terms:

Data-driven control, data-driven MPC, self-triggered control, predicted
control.

##  I Introduction

Thanks to recent advances in data acquisition and computing technologies,
data-driven control has attracted considerable attention in the past years.
Designing control laws directly from data without resorting to any system
identification step, offers an appealing alternative to the traditional model-
based control paradigm [1, 2]. This is because in real-world applications, it
is always difficult or even impossible to acquire an accurate system model [3,
4, 5]. Indeed, a number of publications are devoted to studying data-driven
control. These were mainly inspired by the celebrated Fundamental Lemma
proposed in [6], which lays the theoretical foundation for data-driven
control. Several control problems have been addressed under the new framework,
including stabilization and optimization in [7, 8, 9, 10, 11], linear
quadratic regulation in [12], robust control in [13], quantized control in
[14], model predictive control (MPC) in [15, 16, 17, 18, 19], and control of
complex networks in [20].

Yet, the aforementioned works employ periodic transmission protocols, which
may be resource-inefficient for real-world systems in terms of processor
usage, communication bandwidth, and energy. In cyber-physical networked
systems [21, 22] for instance, whose communication network is shared by many
devices, the communication bandwidth is always restricted for each device. To
tackle this issue, a resource-efficient scheduling approach for data
transmissions, known as event-based control, has been widely studied in the
context of model-based control.

Generally, there are two event-based approaches that have been proven
effective, namely, event-triggered control and self-triggered control [23]. In
the former, an event e.g., a data transmission, is triggered only after some
designed triggering condition is met. This condition should be tested at each
state or output, thus requiring continuous monitoring of the system. While a
level of robustness against uncertainties and unmodeled dynamics can be
expected, having the sensor continuously operating results in waste of
resources. Related work can be found in [24, 25, 26, 27, 28]. Self-triggered
control, on the other hand, determines the next sampling time and transmission
once a sampled measurement is received, which does not need to continuously
sample the outputs; see, e.g., [29, 30, 31, 32]. Notably, the sensors in self-
triggered control can be completely shut off between sampling times, which
saves energy and prolongs service life of the sensor. This feature appears
promising in the model-based case and motivates the generalization of self-
triggered control from model-based to data-driven settings. The key idea of
traditional self-triggered control is to predict the future trajectory using
the system model. It remains unclear how a trajectory can be obtained in the
absence of a model.

The goal of this paper is to design data-driven self-triggered controllers for
unknown linear systems using only some data acquired off-line. The challenge
here lies in how to obtain a predicted trajectory using only data. To this
end, we begin by considering unknown output feedback systems. In this setup, a
new data-driven MPC scheme is proposed, which generates a sequence of optimal
control inputs as well as the associated system trajectory. Leveraging the
predicted trajectory, a data-driven self-triggering mechanism is designed so
that the next triggering time can be dictated without using a system model. We
further extend this method to state feedback control systems, where the state
is sampled to construct the control input based on a state feedback control
law. Again, a norm minimization problem is formulated to predict the system
trajectory using the pre-collected input-state data, whose solution enables
design of the data-driven self-triggering mechanism.

In succinct form, the contribution of this work is threefold, summarized as
follows.

  * c1)

To predict the future system trajectory based on input-output data, a data-
driven MPC scheme accounting for noisy online outputs is developed for output
feedback control systems.

  * c2)

Leveraging the solution of the data-driven MPC, a data-driven self-triggering
mechanism is designed which determines the next transmission time once a new
measurement is received, along with its rigorous and comprehensive stability
analysis.

  * c3)

For state feedback control systems, a norm maximization problem is suggested
to predict the system trajectory, whose solution is then employed to construct
the data-driven self-triggering law.

_Notation:_ Denote the set of real numbers (natural numbers) by $\mathbb{R}$
($\mathbb{N}$), and define $\mathbb{N}_{0}:=\mathbb{N}\cup\\{0\\}$. For a
matrix $M$, if it has full column rank (full row rank), its left pseudo-
inverse (right pseudo-inverse) is denoted by $M^{\dagger}$ ($M^{\ddagger}$).
Given a vector $x\in\mathbb{R}^{n_{x}}$, $\|x\|$ is its Euclidean norm,
$\|x\|_{\infty}$ is its infinity norm, and for a positive definite matrix
$P=P^{T}\succ 0$, define the weighted norm $\|x\|_{P}=\sqrt{x^{\top}Px}$. Let
further $\|M\|$ be the spectral norm of matrix $M$. Let
$\underline{\lambda}_{P}$ ($\overline{\lambda}_{P}$) represent the minimum
(maximum) eigenvalue of matrix $P$. Let $[t_{1},t_{2}]$ denote the time
interval from $t_{1}$ to $t_{2}$ with discrete times. A function
$\alpha:[0,\infty)\rightarrow[0,\infty)$ is said to be of class $\mathcal{K}$
if it is continuous, strictly increasing, and $\alpha(0)=0$. A function
$\alpha:[0,\infty)\rightarrow[0,\infty)$ is said to be of class
$\mathcal{K}_{\infty}$ if it is of class $\mathcal{K}$ and also unbounded. A
function $\beta:[0,\infty)\times[0,\infty)\rightarrow[0,\infty)$ is said to be
of class $\mathcal{KL}$ if $\beta(\cdot,t)$ is of class $\mathcal{K}$ for each
fixed $t\geq 0$ and $\beta(s,t)$ decreases to $0$ as $t\rightarrow\infty$ for
each fixed $s\geq 0$.

The Hankel matrix associated with the sequence $\\{x_{t}\\}_{t=0}^{N-1}$, is
denoted by

$H_{L}(x):=\left[\begin{matrix}x_{0}&x_{1}&\ldots&x_{N-L}\\\
x_{1}&x_{2}&\ldots&x_{N-L+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{L-1}&x_{L}&\ldots&x_{N-1}\end{matrix}\right].$ (1)

In this paper, we consider experiments of length $N$, so the dependence of
$H_{L}(x)$ on $N$ will be omitted. A stacked window of a sequence, say
$\\{x_{t}\\}_{t=t_{1}}^{t_{2}}$ is given by

$x_{[t_{1},t_{2}]}=\left[\begin{matrix}x_{t_{1}}\\\ \vdots\\\
x_{t_{2}}\end{matrix}\right].$ (2)

##  II Preliminaries and Problem Formulation

###  II-A Networked control systems

Consider the following discrete-time linear system

$\displaystyle x_{t+1}$ $\displaystyle=Ax_{t}+Bu_{t},\quad t\in\mathbb{N}_{0}$
(3a) $\displaystyle y_{t}$ $\displaystyle=Cx_{t}+Du_{t}$ (3b)

where $x_{t}\\!\in\\!\mathbb{R}^{n_{x}}$, $u_{t}\\!\in\\!\mathbb{R}^{n_{u}}$,
and $y_{t}\\!\in\\!\mathbb{R}^{n_{y}}$ are the state, control input, and
output, respectively. In this paper, we make the following standing
assumptions on the system (3).

######  Assumption 1 (_Controllability and observability_).

The pair $(A,B)$ is controllable, and $(C,A)$ is observable.

######  Assumption 2 (_Unknown system model_).

The system matrices $(A,B,C,D)$ in (3) are unknown. Instead, some pre-
collected input-output data, i.e., $\\{u^{p}_{t},y^{p}_{t}\\}_{t=0}^{N-1}$ are
available.

Regarding the two assumptions, a remark comes ready.

######  Remark 1.

Assumption 1 is standard for stability analysis of linear systems. As a matter
of fact, the data-driven self-triggering controllers along with associated
theoretical results developed in this paper can be generalized to linear
systems which are stablizable and detectable. In real-world setups, it is
almost always impossible to have perfect knowledge of the system matrices
$(A,B,C,D)$ due to modeling and/or instrumentation inaccuracies. On the other
hand, acquiring input-output data pairs
$\\{u^{p}_{t},y^{p}_{t}\\}_{t=0}^{N-1}$ by exciting the system with some
control inputs is often practical and doable.

We further recall the definition of observability index [33].

######  Definition 1 (_Observability index_).

The observability index of linear systems as in (3) is defined to be the
smallest integer $\eta\in\\{1,\cdots,n_{x}\\}$ such that the observability
matrix $\Theta$ has full rank $n_{x}$, i.e.,

${\rm rank}(\Theta)\stackrel{{\scriptstyle\triangle}}{{=}}{\rm
rank}\big{(}[C^{\top}\,(CA)^{\top}\,\cdots\,(CA^{\eta-1})^{\top}]^{\top}\big{)}=n_{x}.$
(4)

Figure 1: Pictorial description of system (3) with data-driven self-triggered
MPC.

In networked control systems, the plant is connected with a sensor which takes
the output measurement at sampling times. Sampled measurements are transmitted
to a remote controller to construct control inputs. Frequent sampling and
transmission inevitably results in waste of energy and network resources. To
save the network bandwidth as well as prolong the life-cycle of the sensor, we
incorporate a self-triggering module running at the controller side, to
dictate when a new output measurement should be sampled and transmitted.

Specifically, we consider that the output is transmitted over a network of
limited bandwidth to the remote controller, while the controller-to-plant
channel is assumed perfect either over wired lines or networks with sufficient
bandwidth. Furthermore, to model the network effect, the output is corrupted
by bounded additive noise $\|n_{t}\|\leq\bar{n}$ for all $t\in\mathbb{N}_{0}$
when transmitted over the network, i.e.,

$\zeta_{t}=y_{t}+n_{t}.$ (5)

See Fig. 1 for a pictorial depiction of the system structure.

Let $t_{\ell}\in\mathbb{N}_{0}$ stand for the time at which the $\ell$-th
event (i.e., sampling and transmission) is triggered by the self-triggering
module. Assume the sensor is equipped with a buffer of size $\eta$, which
records the most recent $\eta$ historical measurements. Once an event is
triggered, the data packet consisting of the $\eta$ buffered measurements is
sent to the controller. The self-triggering module computes the next
triggering time by evaluating a triggering function on a predicted trajectory
along with the received data packet. For example, at time $t_{\ell}$, the
$\ell$-th transmission has just been made, and our self-triggering module
receives the data packet
$y_{[t_{\ell}-\eta,t_{\ell}-1]}:=\\{y_{t_{\ell}-\eta},y_{t_{\ell}-\eta+1},\cdots\\!,y_{t_{\ell}-1}\\}$,
and computes the next triggering time $t_{\ell+1}$ which is fed back to the
sensor. During the interval $[t_{\ell-1},t_{\ell})$, it is only required that
the sensor samples output measurements for the sub-interval
$[t_{\ell}-\eta,t_{\ell})$, which is often much smaller than
$[t_{\ell-1},t_{\ell})$, therefore saving considerable energy as well as
expanding life expectancy of the sensor. Of course, if the inter-triggering
time $t_{\ell}-t_{\ell-1}$ is smaller than $\eta$, i.e.,
$t_{\ell}-t_{\ell-1}\leq\eta$, one only needs to send the packet of the new
$t_{\ell}-t_{\ell-1}$ measurements as the old $\eta-(t_{\ell}-t_{\ell-1})$
ones have already been sent to the controller in the last packet.

###  II-B Fundamental Lemma

In this section, we briefly review the data-driven system representation based
on the so-called Fundamental Lemma [6], which is key to derive our data-driven
self-triggering predictive controllers. Before presenting the lemma, the
definition of persistence of excitation is introduced first.

######  Definition 2 (_Persistence of excitation_).

A sequence $\\{u_{t}\in\mathbb{R}^{n_{u}}\\}_{t=0}^{N-1}$ is persistently
exciting of order $L$ if ${\rm{rank}}(H_{L}(u))=n_{u}L$.

Based on Def. 2, it has been shown in [6] that any trajectory of system (3)
can be expressed as a linear combination of pre-collected input-output data
$\\{u^{p}_{t},y^{p}_{t}\\}_{t=0}^{N-1}$, provided that the input sequence
$\\{u^{p}_{t}\\}_{t=0}^{N-1}$ is persistently exciting. This result, formally
summarized below, is also known as the Fundamental Lemma [6].

######  Lemma 1 (_Fundamental Lemma_).

Consider a persistently exciting input sequence $\\{u^{p}_{t}\\}_{t=0}^{N-1}$
of order $L+n_{x}$. A trajectory $\\{\bar{u}_{t},\bar{y}_{t}\\}_{t=0}^{L-1}$
is generated by system (3) if and only if there exists a vector
$g\in\mathbb{R}^{N-L+1}$ such that the following holds

$\left[\begin{matrix}H_{L}\\!\left(u^{p}\right)\\\
H_{L}\\!\left(y^{p}\right)\end{matrix}\right]g=\left[\begin{matrix}\bar{u}_{[0,L-1]}\\\
\bar{y}_{[0,L-1]}\end{matrix}\right].$ (6)

It can be concluded from Lemma 1 that to generate a predicted trajectory of
length $L$, the pre-collected data should be persistently exciting of at least
order $L$. To validate this requirement, the next assumption is imposed.

######  Assumption 3 (_Pre-collected data_).

Let $\\{u^{p}_{t}\\}_{t=0}^{N-1}$ be any sequence persistently exciting of
order $L+n_{x}+\eta$. The output sequence $\\{y^{p}_{t}\\}_{i=0}^{N-1}$ is
generated from system (3) offline with any initial condition $x^{p}_{0}$, and
input sequence $\\{u^{p}_{t}\\}_{t=0}^{N-1}$.

For notational brevity, we use $u^{p}$ and $y^{p}$ to represent sequences
$\\{u^{p}_{t}\\}_{t=0}^{N-1}$ and $\\{y^{p}_{t}\\}_{t=0}^{N-1}$ in the
following, respectively. Leveraging the observability index $\eta$ in Def. 1,
we construct an extended state for $t\geq\eta$ as follows

$\xi_{t}:=\left[\begin{matrix}u_{[t-\eta,t-1]}\\\
y_{[t-\eta,t-1]}\end{matrix}\right]\in\mathbb{R}^{n_{\xi}}$ (7)

with $n_{\xi}:=(n_{u}+n_{y})\eta$. Then, similar to [34], system (3) can be
converted into the ensuing linear system

$\displaystyle\xi_{t+1}$ $\displaystyle=\tilde{A}\xi_{t}+\tilde{B}u_{t},\quad
t\in\mathbb{N}_{\eta}$ (8a) $\displaystyle y_{t}$
$\displaystyle=\tilde{C}\xi_{t}+\tilde{D}u_{t}$ (8b)

for suitable matrices $(\tilde{A},\tilde{B},\tilde{C},\tilde{D})$ depending on
$(A,B,C,D)$. For subsequent analysis, let us denote the equilibrium point of
the new system (8) by

$\xi^{e}:=\left[\begin{matrix}u_{[0,\eta-1]}\\\
y_{[0,\eta-1]}\end{matrix}\right]$ (9)

in which $u_{i}=u^{e}$ and $y_{i}=y^{e}$ for all $i=0,1,\cdots,\eta-1$. As a
matter of fact, such an augmented system (state) has been commonly employed in
studies of data-driven MPC to simplify the stability analysis; see e.g., [15,
35].

##  III Data-Driven self-triggering MPC

This section advocates a data-driven approach to designing a controller and a
self-triggering mechanism for unknown linear systems as in (3). Commonly, a
self-triggering mechanism determines the next transmission by comparing the
current output with a predicted future output, and triggers a transmission if
they differ considerably. However, when the system matrices are unknown, three
grand challenges are there: c1) how to obtain a predicted system trajectory
for the self-triggering mechanism using only pre-collected data? c2) how to
design stablizing data-driven control inputs? and, c3) how to perform
performance analysis for the data-driven self-triggered closed-loop control
system? We first address c1) and c2) by putting forward a data-driven MPC
scheme, which can compute a window of control inputs and its corresponding
predicted outputs once a new output measurement is transmitted and received.

###  III-A Data-driven MPC scheme

Although Lemma 1 characterizes system trajectories in some sense, the
correspondence between the input-output trajectory and the coefficient vector
is not unique in general. In addition, the Fundamental Lemma does not account
for any noise, so it cannot be directly employed for control and trajectory
prediction. Building on existing works [15, 19], we propose a new data-driven
MPC that, at each triggered time $t_{\ell}$, computes a window of optimal
control inputs and associated outputs, using the received noisy data
$\zeta_{[t_{\ell}-\eta,t_{\ell}-1]}$. The control inputs are subsequently
employed one by one at times $t\in[t_{\ell},t_{\ell+1})$. Let vectors
$\bar{u}(t_{\ell}):=[\bar{u}^{\top}_{-\eta}(t_{\ell})\,\cdots\,\bar{u}^{\top}_{L-1}(t_{\ell})]^{\top}$
and
$\bar{y}(t_{\ell}):=[\bar{y}^{\top}_{-\eta}(t_{\ell})\,\cdots\,\bar{y}^{\top}_{L-1}(t_{\ell})]^{\top}$
denote the predicted input and output from time $t_{\ell}-\eta$ to
$t_{\ell}+L-1$, and
$\bar{\xi}(t_{\ell})=[\bar{\xi}_{0}^{\top}(t_{\ell})\cdots\bar{\xi}_{L}^{\top}(t_{\ell})]^{\top}$
the corresponding extended state with

$\displaystyle\bar{\xi}_{i}(t_{\ell})=\left[\begin{matrix}\bar{u}_{[i-\eta,i-1]}(t_{\ell})\\\
\bar{y}_{[i-\eta,i-1]}(t_{\ell})\end{matrix}\right]$

where $i\in\\{0,1,\cdots\\!,L\\}$. Mathematically, the following optimization
problem is solved at each self-triggered time $t_{\ell}$

$\displaystyle J^{*}_{L}(u_{t_{\ell}},$ $\displaystyle y_{t_{\ell}}):=$
$\displaystyle\underset{\\{g(t_{\ell}),h(t_{\ell})\atop\bar{y}(t_{\ell}),\bar{u}(t_{\ell})\\}}{\min}$
$\displaystyle\sum_{i=0}^{L-1}\|\bar{u}_{i}(t_{\ell})-u^{e}\|^{2}_{R}+\|\bar{y}_{i}(t_{\ell})-y^{e}\|^{2}_{Q}+\frac{\lambda_{h}}{\bar{n}}\|h(t_{\ell})\|^{2}$
$\displaystyle+\lambda_{g}\bar{n}\|g(t_{\ell})\|^{2}+\|\bar{\xi}_{L}(t_{\ell})-\xi^{e}\|^{2}_{P}$
(10a) $\displaystyle{\rm s.t.}\quad$
$\displaystyle\left[\begin{matrix}\bar{u}(t_{\ell})\\\
\bar{y}(t_{\ell})+h(t_{\ell})\end{matrix}\right]=\left[\begin{matrix}H_{L+\eta}(u^{p})\\\
H_{L+\eta}(y^{p})\end{matrix}\right]g(t_{\ell})$ (10b)
$\displaystyle\left[\begin{matrix}\bar{u}_{[-\eta,-1]}(t_{\ell})\\\
\bar{y}_{[-\eta,-1]}(t_{\ell})\end{matrix}\right]=\left[\begin{matrix}u_{[t_{\ell}-\eta,t_{\ell}-1]}\\\
\zeta_{[t_{\ell}-\eta,t_{\ell}-1]}\end{matrix}\right]$ (10c)
$\displaystyle\bar{\xi}_{L}(t_{\ell})\in\Xi_{\epsilon}$ (10d)
$\displaystyle\bar{u}_{i}\in\mathbb{U},\quad i\in\\{1,2,\cdots\\!,L-1\\}.$
(10e)

The data-driven MPC formulation (10) generalizes that of [35] by accommodating
noisy outputs. Preselected weighting matrices $R\succ 0$ and $Q\succ 0$
penalize the distances from the predicted input and output to their
equilibrium points. Constraint (10b) reflects the data-driven system
representation in Lemma 1, where $g(t_{\ell})\in\mathbb{R}^{N-L-\eta+1}$
stands for the coefficient vector at time $t_{\ell}$, and the vector
$h(t_{\ell})=[h^{\top}_{-\eta}(t_{\ell})~{}\cdots~{}h^{\top}_{L-1}(t_{\ell})]^{\top}$
is added to compensate for the unknown network-induced noise $n_{t}$ in
predicting future outputs. Penalties are also imposed on the use of
$g(t_{\ell})$ and $h(t_{\ell})$ in the objective function, to guarantee system
stability and improve robustness against noise, with weights $\lambda_{h}>0$
and $\lambda_{g}>0$ balancing between minimizing different costs. Feasibility
sets $\mathbb{U}$ and $\Xi_{\epsilon}$ are prescribed and convex. Finally,
(10d) is a terminal constraint, which together with the last summand having
$P\succ 0$ in (10a) ensures that the last predicted extended state
$\bar{\xi}_{L}(t_{\ell})$ (at the end of the window) stays in $\Xi_{\epsilon}$
(i.e., a small neighborhood of the desired equilibrium point). Overall,
problem (10) is convex and can be efficiently solved using off-the-shelf
solvers.

The terminal constraint (10d) and the terminal cost
$\|\bar{\xi}_{L}(t_{\ell})-\xi^{e}\|^{2}_{P}$ are also critical ingredients in
the model-based MPC, see e.g., [36]. To proceed, we make the following
assumption on the terminal ingredients $\Xi_{\epsilon}$ and the corresponding
matrix $P$.

######  Assumption 4.

There exist matrices $P=P^{\top}\succ 0$, $K\in\mathbb{R}^{n_{u}\times
n_{\xi}}$, and a set
$\Xi_{\epsilon}:\\{\xi\in\mathbb{R}^{n_{\xi}}|\|\xi-\xi^{e}\|_{P}\leq\epsilon\\}\subseteq\mathbb{U}^{\eta}\times\mathbb{R}^{\eta}$
such that for all $\xi\in\Xi_{\epsilon}$, $u=u^{e}+K(\xi-\xi^{e})$, and
$y=(\tilde{C}+\tilde{D}K)\xi$, the following statements hold true

  1. 1.

$u\in\mathbb{U}$, $\tilde{A}\xi+\tilde{B}u\in\Xi_{\epsilon}$ and,

  2. 2.

the following inequality holds

$\|(\tilde{A}+\tilde{B}K)\xi\|_{P}^{2}\leq\|\xi\|_{P}^{2}-\|K\xi\|_{R}^{2}-\|y\|_{Q}^{2}.$
(11)

A feasible data-based method for selecting the matrix $P$ and the set
$\Xi_{\epsilon}$ has been discussed in [35, Proposition 10]. It will be shown
in the next subsection that these terminal ingredients have influence on the
inter-triggering time, and should be handled with care. According to (10c), a
requirement is imposed on the prediction window as follows.

######  Assumption 5 (_Prediction horizon_).

The predict horizon of the data-driven MPC satisfies $L\geq\eta+1$.

Without loss of generality, we consider for simplicity the equilibrium
$\xi^{e}=[(u^{e})^{\top}~{}(y^{e})^{\top}]^{\top}=0$ in the remainder of the
paper. Let
$(\bar{y}^{*}(t_{\ell}),\bar{u}^{*}(t_{\ell}),g^{*}(t_{\ell}),h^{*}(t_{\ell}))$
denote any optimal solution of problem (10) at $t_{\ell}$, thereby yielding
the predicted extended state $\bar{\xi}^{*}(t_{\ell})$. Let $u_{t_{\ell}}$,
$y_{t_{\ell}}$ and $\xi_{t_{\ell}}$ represent the actual input, output and
extended state, respectively. The inter-triggering time between two
consecutive self-triggered times is defined as
$\tau_{\ell}:=t_{\ell+1}-t_{\ell}$ for $\ell\in\mathbb{N}_{0}$, with
$t_{0}=0$. With these definitions, the following lemma indicates that the
error between the predicted extended state $\bar{\xi}^{*}(t_{\ell})$ and the
actual extended state $\xi_{t}$ is bounded, and it can be rigorously
quantified using the optimizer of (10).

######  Lemma 2.

Let Assumptions 1—5 hold. Consider system (3) whose control input is generated
by solving problem (10) at each self-triggered time $t_{\ell}$. For every
$\ell\in\mathbb{N}_{0}$, the error between $\xi_{t_{\ell+1}}$ and
$\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})$ satisfies

$\displaystyle\|\xi_{e}(t_{\ell+1})\|$
$\displaystyle:=\|\xi_{t_{\ell+1}}-\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|$
(12a)
$\displaystyle\leq\big{(}\sqrt{\eta}\bar{n}+\|h^{*}_{[-\eta,-1]}(t_{\ell})\|\big{)}\sqrt{\sum_{i=\tau_{\ell}}^{\tau_{\ell}+\eta-1}\rho^{i}}$
(12b)
$\displaystyle~{}~{}~{}+\|h^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\|$
(12c)

where $\rho^{i}:=\|CA^{i+\eta}\Phi^{\dagger}\|$ for
$i=\tau_{\ell},\tau_{\ell}+1,\cdots,\tau_{\ell}+\eta-1$.

###### Proof.

It follows from (10d) that

$\displaystyle\xi_{e}(t_{\ell}+\tau_{\ell})$
$\displaystyle=\left[\begin{matrix}\bar{u}^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\\\
y_{[t_{\ell}+\tau_{\ell}-\eta,t_{\ell}+\tau_{\ell}-1]}\end{matrix}\right]\\!-\\!\left[\begin{matrix}\bar{u}^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\\\
\bar{y}^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\end{matrix}\right]$
$\displaystyle=\left[\begin{matrix}0\\\
y_{[t_{\ell}+\tau_{\ell}-\eta,t_{\ell}+\tau_{\ell}-1]}-\bar{y}^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\end{matrix}\right].$

Therefore, an upper bound on the extended state error
$\|\xi_{e}(t_{\ell}+\tau_{\ell})\|$ can be obtained by calculating
$y_{e}(t_{\ell}+\tau_{\ell}):=y_{t_{\ell}+\tau_{\ell}}-\bar{y}^{*}_{\tau_{\ell}}(t_{\ell})$.
It follows from (10b) that

$\displaystyle\bar{y}^{*}_{\tau_{\ell}}(t_{\ell})=I_{\tau_{\ell}}H_{\eta+L}(y^{p})g^{*}(t_{\ell})-h^{*}_{\tau_{\ell}}(t_{\ell})$
(13)

where $I_{\tau_{\ell}}$ denote the corresponding $\tau_{\ell}$-th block of
matrix $H_{\eta+L}(y^{p})g^{*}(t_{\ell})$ for the optimizer
$\bar{y}^{*}_{\tau_{\ell}}(t_{\ell})$. Let
$\hat{y}_{\tau_{\ell}}(t_{\ell}):=I_{\tau_{\ell}}H_{\eta+L}(y^{p})g^{*}(t_{\ell})$,
which is a trajectory of system (3) according to Lemma 1. Based on (3), one
can deduce that

$\displaystyle\hat{y}_{\tau_{\ell}}(t_{\ell})\\!-\\!y_{t_{\ell}+\tau_{\ell}}\\!=\\!CA^{\tau_{\ell}+\eta}\Theta^{{\dagger}}(\hat{y}_{[-\eta,-1]}(t_{\ell})\\!-\\!y_{[t_{\ell}-\eta,t_{\ell}-1]}).$

It follows from (13) and (10c) that

$\displaystyle\hat{y}_{[-\eta,-1]}(t_{\ell})-y_{[t_{\ell}-\eta,t_{\ell}-1]}$
$\displaystyle~{}~{}=y_{[t_{\ell}-\eta,t_{\ell}-1]}+n_{[t_{\ell}-\eta,t_{\ell}-1]}+h^{*}_{[-\eta,-1]}(t_{\ell})-y_{[t_{\ell}-\eta,t_{\ell}-1]}$
$\displaystyle~{}~{}=n_{[t_{\ell}-\eta,t_{\ell}-1]}+h^{*}_{[-\eta,-1]}(t_{\ell}).$

Therefore, the error between the actual output and the predicted output obeys

$\displaystyle
y_{e}(t_{\ell}\\!+\\!\tau_{\ell})\\!=\\!CA^{\tau_{\ell}+\eta}\Theta^{\dagger}\big{(}n_{[t_{\ell}-\eta,t_{\ell}-1]}\\!+\\!h^{*}_{[-\eta,-1]}(t_{\ell})\big{)}\\!-\\!h^{*}_{\tau_{\ell}}(t_{\ell}).$
(14)

We have that

$\displaystyle\|\xi_{e}(t_{\ell}+\tau_{\ell})\|$
$\displaystyle=\Big{\|}\Big{[}(CA^{\tau_{\ell}+\eta})^{\top}~{}\cdots~{}(CA^{\tau_{\ell}+2\eta-1})^{\top}\Big{]}^{\top}$
$\displaystyle\quad\times\big{[}n_{[t_{\ell}-\eta,t_{\ell}-1]}+h^{*}_{[-\eta,-1]}(t_{\ell})\big{]}$
$\displaystyle\quad-\big{[}(h^{*}_{\tau_{\ell}}(t_{\ell}))^{\top}~{}\cdots~{}(h^{*}_{\tau+\eta-1}(t_{\ell}))^{\top}\big{]}^{\top}\Big{\|}$
$\displaystyle\leq\big{(}\sqrt{\eta}\bar{n}+\|h^{*}_{[-\eta,-1]}(t_{\ell})\|\big{)}\sqrt{\sum_{i=\tau_{\ell}}^{\tau_{\ell}+\eta-1}\rho^{i}}$
$\displaystyle\quad+\|h^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\|$

which completes the proof. ∎

According to Lemma 2, an upper bound on the error between the predicted
extended state and the actual extended state is characterized in terms of
known or computable parameters $h^{*}(t_{\ell})$, $\bar{n}$, and unknown
parameters $\rho^{i}$ for $i=\tau_{\ell},\cdots,\tau_{\ell}+\eta-1$. Before
running the system online, upper bounds of parameters $\rho^{i}$ can be
obtained offline by solving the optimization problem proposed in [18, Section
V. B] for $i=1,\cdots,L-1$

$\displaystyle J_{i}^{\prime}:=$
$\displaystyle\max_{y,g^{\prime}}~{}\|y_{i}\|_{\infty}$ (15a)
$\displaystyle\,~{}~{}{\rm s.t.}~{}\|y_{[0,\eta-1]}\|_{\infty}\leq 1$ (15b)
$\displaystyle\,\quad\quad\left[\begin{matrix}H_{i+1}(u^{p})\\\
H_{i+1}(y^{p})\end{matrix}\right]g^{\prime}=\left[\begin{matrix}0\\\
y_{[0,i]}\end{matrix}\right]$ (15c)

where $H_{i+1}(u^{p})$ denotes the first $(i+1)n_{u}$ rows of matrix
$H_{L+\eta}(u^{p})$, and likewise for $H_{L+\eta}(y^{p})$; vector
$g\prime\in\mathbb{R}^{N-L-\eta+1}$, and $y_{[0,i]}$ consists of the $i+1$
predicted outputs. It has been shown in [18] that
$J_{i}^{\prime}\geq\|CA^{i+\eta}\Phi^{\dagger}\|$, and hence, an upper bound
for $\rho^{i}$ is obtained, i.e., $J_{i}^{\prime}\geq\rho^{i}$. Therefore, one
gets an upper bound on $\|\xi_{e}(t_{\ell+1})\|$.

###  III-B Self-triggering mechanism

Observe that the prediction horizon constrains at most $L-1$ future outputs at
time $t_{\ell}$ can be obtained by solving (10). Therefore, the inter-
triggering time between any two consecutive self-triggered times obeys
$\tau_{\ell}\leq L-1$ for all $\ell\in\mathbb{N}_{0}$. Leveraging the upper
bound in Lemma 2, our self-triggering mechanism is given in the following
lemma.

######  Lemma 3.

Let Assumptions 1—5 hold. Consider system (3) with control inputs generated by
the data-driven MPC (10) at each self-triggered time $t_{\ell}$ for all
$\ell\in\mathbb{N}_{0}$. Assume further that (10) is feasible at $t_{0}$. For
appropriate $\lambda_{g}>0$, $\lambda_{h}>0$, $\epsilon\geq 0$, and $P\succ 0$
satisfying Assumption 4, there exists a constant $\bar{\sigma}\in(0,1)$ such
that for all $\sigma\in(0,\bar{\sigma})$, problem (10) is feasible at all
$t_{\ell}\in\mathbb{N}_{0}$, if i) the constant $r$, matrices $P_{r}$ and
$K_{r}$ are chosen obeying Assumption 4 and

$\displaystyle\frac{\bar{\lambda}_{P_{r}}}{\underline{\lambda}_{P}}\Big{(}1-\frac{\underline{\lambda}_{K_{r}^{\top}RK_{r}}}{\bar{\lambda}_{P_{r}}}\Big{)}^{L}r^{2}\leq\epsilon^{2}\leq
r^{2}$ (16)

and ii) the inter-triggering time satisfies

$\displaystyle\tau_{\ell}=\min\big{\\{}\hat{\tau}_{\ell},\,\check{\tau}_{\ell},\,L-1\big{\\}}$
(17)

with

$\displaystyle\hat{\tau}_{\ell}$
$\displaystyle:=\sup\bigg{\\{}\tau_{\ell}\Big{|}(\sqrt{\eta}\bar{n}+\|h^{*}_{[-\eta,-1]}(t_{\ell})\|)\sqrt{\sum_{i=\tau_{\ell}}^{\tau_{\ell}+\eta-1}\rho^{i}}$
$\displaystyle~{}\quad\quad\quad+\|h^{*}_{[\tau_{\ell}-\eta,\tau_{\ell}-1]}(t_{\ell})\|\leq\frac{r}{\bar{\lambda}_{P_{r}}}-\|\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|\bigg{\\}}$
(18)

and

$\displaystyle\check{\tau}_{\ell}$
$\displaystyle:=\sup\bigg{\\{}\tau_{\ell}\Big{|}2\Big{[}\bar{\lambda}_{Q}+2\bar{\lambda}_{P_{r}}+\lambda_{g}\bar{n}\|H^{\ddagger}_{u\xi}\|^{2}\Big{(}1+\frac{\bar{\lambda}_{P_{r}}}{\underline{\lambda}_{R}}\Big{)}\Big{]}$
$\displaystyle\times\Big{[}\big{(}\eta\bar{n}^{2}+\|h^{*}_{[-\eta,-1](t_{\ell})}\|^{2}\big{)}\sum_{i=0}^{\tau_{\ell}-1}(\rho^{i+\eta})^{2}+\sum_{i=0}^{\tau_{\ell}-1}\|h^{*}_{i}(t_{\ell})\|^{2}\Big{]}$
$\displaystyle+\lambda_{h}\eta\bar{n}^{2}+\epsilon^{2}+2\Big{[}\bar{\lambda}_{P_{r}}+\lambda_{g}\bar{n}\|H_{u\xi}^{\ddagger}\|^{2}\Big{(}1+\frac{\bar{\lambda}_{P_{r}}}{\underline{\lambda}_{R}}\Big{)}\Big{]}$
$\displaystyle\times\|\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|^{2}\leq\sigma\sum_{i=0}^{\tau_{\ell}-1}\|\bar{\xi}^{*}_{i}(t_{\ell})\|^{2}\bigg{\\}}$
(19)

where matrix $H_{u\xi}$ is defined by

$\displaystyle H_{u\xi}=\left[\begin{matrix}H_{L+\eta}(u^{p})\\\
H_{1}(\xi^{p}_{N-L-\eta+1})\end{matrix}\right].$ (20)

###### Proof.

First, for some $\ell\in\mathbb{N}_{0}$, we assume that the problem is
feasible at $t_{\ell}$ and construct a candidate solution at $t_{\ell+1}$.
Then, leveraging a carefully designed Lyapunov function, recursive feasibility
of problem (10) is proven by showing the decrease of the Lyapunov function.

Suppose that (10) is feasible at $t_{\ell}$ for some $\ell\in\mathbb{N}_{0}$.
Denote a candidate solution of (10) at $t_{\ell+1}$ by $\bar{y}(t_{\ell+1})$,
$\bar{u}(t_{\ell+1})$, $g(t_{\ell+1})$, and $h(t_{\ell+1})$. Using (12), one
gets from (17)–(18) that

$\displaystyle\|\xi_{t_{\ell+1}}\|_{P_{r}}\leq\|\xi_{e}(t_{\ell+1})\|_{P_{r}}+\|\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|_{P_{r}}\leq
r.$

Thus, for terminal region
$\Xi_{r}:\\{\xi\in\mathbb{R}^{n_{\xi}}|\|\xi-\xi^{e}\|_{P_{r}}\leq r\\}$ there
exist matrices $P_{r}$ and $K_{r}$ such that conditions in Assumption 4 hold.
Let $\bar{u}_{i}(t_{\ell+1})=K\bar{\xi}_{i}(t_{\ell+1})$ and
$\bar{y}_{i}(t_{\ell+1})=y_{i+t_{\ell+1}}$ for $i\in\\{0,\cdots\\!,L-1\\}$.
Capitalizing on (10b), the slack variable $h(t_{\ell+1})$ can be chosen as
$h_{[-\eta,-1]}(t_{\ell+1})=n_{[t_{\ell+1}-\eta,t_{\ell+1}-1]}$ and
$h_{i}(t_{\ell+1})=0$ for $i\in\\{0,1,\cdots,L-1\\}$. Set

$\displaystyle{g}(t_{\ell+1})=H^{\ddagger}_{u\xi}\left[\begin{matrix}u_{[0,L-1]}(t_{\ell+1})\\\
\xi(t_{\ell+1})\end{matrix}\right]$ (21)

with $H_{u\xi}$ in (20). It follows from (11) that

$\displaystyle\|\bar{\xi}_{i+1}(t_{\ell+1})\|^{2}_{P_{r}}$
$\displaystyle\leq\|\bar{\xi}_{i}(t_{\ell+1})\|^{2}_{P_{r}}-\|\bar{\xi}_{i}(t_{\ell+1})\|_{K_{r}^{\top}RK_{r}}^{2}$
$\displaystyle\leq\big{[}1-\underline{\lambda}_{K_{r}^{\top}RK_{r}}/\bar{\lambda}_{P_{r}}\big{]}\|\bar{\xi}_{i}(t_{\ell+1})\|_{P_{r}}^{2}$

for all $i\in\\{0,\cdots,L-1\\}$. Recursively, we arrive at

$\displaystyle\|\bar{\xi}_{L}(t_{\ell+1})\|_{P_{r}}^{2}\leq\big{[}1-\underline{\lambda}_{K_{r}^{\top}RK_{r}}/\bar{\lambda}_{P_{r}}\big{]}^{L}\|\bar{\xi}_{0}(t_{\ell+1})\|_{P_{r}}^{2}.$
(22)

Combining (16) with $\|\bar{\xi}_{0}(t_{\ell+1})\|_{P_{r}}^{2}\leq r^{2}$,
inequality (22) implies that
$\|\bar{\xi}_{L}(t_{\ell+1})\|_{P}^{2}\leq\epsilon^{2}$. Therefore, if problem
(10) is feasible at $t_{\ell}$, it is feasible at $t_{\ell+1}$.

Next, we construct a Lyapunov function to show the recursive feasibility.
Since $(C,A)$ is observable, $(\tilde{C},\tilde{A})$ is detectable. According
to [37], there exists an input-output-to-state stability Lyapunov function
$W(\xi)=\xi^{\top}P_{\xi}\xi$ such that

$W(\tilde{A}\xi+\tilde{B}u)-W(\xi)\leq-\frac{1}{2}\|\xi\|^{2}+c_{1}\|u\|^{2}+c_{2}\|y\|^{2}$
(23)

for suitable $c_{1},c_{2}>0$, $P_{\xi}\succ 0$, and all
$u\in\mathbb{R}^{n_{u}}$, $\xi\in\mathbb{R}^{n_{\xi}}$,
$y=\tilde{C}\xi+\tilde{D}\xi$. Consider the Lyapunov function
$V(\xi_{t})=J_{L}^{*}(\xi_{t})+\gamma W(\xi_{t})$ having
$\gamma=\min\\{\underline{\lambda}_{Q},\underline{\lambda}_{R}\\}/{\max\\{c_{1},2c_{2}\\}}$.
Recalling the candidate solution
$(\bar{u}(t_{\ell+1}),\bar{y}(t_{\ell+1}),g(t_{\ell+1}),h(t_{\ell+1}))$ for
problem (10) in Lemma 3. Denote the corresponding cost function by
$\bar{J}_{L}(\xi_{t_{\ell}})$. The optimizer of (10) at $t_{\ell}$ is denoted
by $(\bar{u}^{*}(t_{\ell}),\
\bar{y}^{*}(t_{\ell}),g^{*}(t_{\ell}),h^{*}(t_{\ell}))$. The difference of the
Lyapunov function at two consecutive self-triggered times satisfies

$\displaystyle V(\xi_{t_{\ell+1}})-V(\xi_{t_{\ell}})$ (24) $\displaystyle=\;$
$\displaystyle\bar{J}_{L}(\xi_{t_{\ell+1}})+\gamma
W(\xi_{t_{\ell+1}})-J^{*}_{L}(\xi_{t_{\ell}})-\gamma W(\xi_{t_{\ell}})$ (25)
$\displaystyle\leq\;$
$\displaystyle\underbrace{\big{(}\bar{J}_{L}(\xi_{t_{\ell+1}})-J^{*}_{L}(\xi_{t_{\ell}})\big{)}}_{\stackrel{{\scriptstyle\triangle}}{{=}}\mathcal{T}_{1}}+\underbrace{\gamma\big{(}W(\xi_{t_{\ell+1}})-W(\xi_{t_{\ell}})\big{)}}_{\stackrel{{\scriptstyle\triangle}}{{=}}\mathcal{T}_{2}}.$
(26)

Next, we derive upper bounds on the two terms.

_Part a: Upper-bounding $\mathcal{T}_{1}$._ According to (10a), it can be
obtained that

$\displaystyle\mathcal{T}_{1}$
$\displaystyle\leq\sum_{i=0}^{L-1}\big{(}\|\bar{u}_{i}(t_{\ell+1})\|^{2}_{R}+\|\bar{y}_{i}(t_{\ell+1})\|^{2}_{Q}\big{)}$
$\displaystyle\quad+\lambda_{g}\bar{n}\|g(t_{\ell+1})\|^{2}+(\lambda_{h}/\bar{n})\|h(t_{\ell+1})\|^{2}+\|\bar{\xi}_{L}(t_{\ell+1})\|^{2}_{P}$
$\displaystyle\quad-\Big{[}\sum_{i=0}^{L-1}\big{(}\|\bar{u}^{*}_{i}(t_{\ell})\|^{2}_{R}+\|\bar{y}^{*}_{i}(t_{\ell})\|^{2}_{Q}\big{)}+\lambda_{g}\bar{n}\|g^{*}(t_{\ell+1})\|^{2}$
$\displaystyle\quad+(\lambda_{h}/\bar{n})\|h^{*}(t_{\ell+1})\|^{2}+\|\bar{\xi}^{*}_{L}(t_{\ell+1})\|^{2}_{P}\Big{]}.$

In addition, since $\bar{u}_{i}(t_{\ell+1})=K\bar{\xi}(t_{\ell+1})$, it can be
deduced recursively from (11) that

$\displaystyle\sum_{i=0}^{L-1}\big{(}\|\bar{u}_{i}(t_{\ell+1})\|^{2}_{R}+\|\bar{y}_{i}(t_{\ell+1})\|^{2}_{Q}\big{)}\leq\|{\xi}_{t_{\ell+1}}\|_{P_{r}}^{2}\\!-\\!\|\bar{\xi}_{L}(t_{\ell})\|_{P_{r}}^{2}$
$\displaystyle\leq
2\bar{\lambda}_{P_{r}}\big{(}\|\xi_{e}(t_{\ell+1})\|^{2}+\|\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|^{2}\big{)}.$

Notice from (21) that

$\displaystyle\quad~{}\|g(t_{\ell+1})\|^{2}\leq\|H^{{\ddagger}}_{u\xi}\|^{2}(\|\bar{u}_{[0,L-1](t_{\ell})}\|^{2}+\|\xi_{t_{\ell+1}}\|^{2})$
$\displaystyle\leq\|H^{{\ddagger}}_{u\xi}\|^{2}(1+\bar{\lambda}_{P_{r}}/\underline{\lambda}_{R})\|\xi_{t_{\ell+1}}\|^{2}$
$\displaystyle\leq
2\|H^{{\ddagger}}_{u\xi}\|^{2}\Big{(}1+\frac{\bar{\lambda}_{P_{r}}}{\underline{\lambda}_{R}}\Big{)}\big{(}\|\xi_{e}(t_{\ell+1})\|^{2}+\|\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|^{2}\big{)}.$

Since $h_{[-\eta,-1]}(t_{\ell+1})=n_{[t_{\ell+1}-\eta,t_{\ell+1}-1]}$ and
$h_{i}(t_{\ell+1})=0$ for $i\in\\{0,\cdots,L-1\\}$, one has that

$(\lambda_{h}/\bar{n})\|h(t_{\ell+1})\|^{2}\leq\lambda_{h}\eta\bar{n}.$

Part b: Upper-bounding $\mathcal{T}_{2}$. It follows recursively from (23)
that

$\displaystyle\quad W(\xi_{t_{\ell+1}})-W(\xi_{t_{\ell}})$
$\displaystyle\leq-\frac{1}{2}\|\xi_{[t_{\ell},t_{\ell+1}-1]}\|^{2}+c_{1}\|\bar{u}^{*}_{[0,\tau_{\ell}-1]}(t_{\ell})\|^{2}+c_{2}\|y_{[t_{\ell},t_{\ell+1}]}\|^{2}$
$\displaystyle\leq-\frac{1}{2}\|\xi_{[t_{\ell},t_{\ell+1}-1]}\|^{2}\\!+\\!c_{1}\|\bar{u}^{*}_{[0,\tau_{\ell}-1]}(t_{\ell})\|^{2}\\!+\\!2c_{2}\|\bar{y}^{*}_{[0,\tau_{\ell}-1]}(t_{\ell})\|^{2}$
$\displaystyle\quad+2c_{2}\|y_{[t_{\ell},t_{\ell+1}]}-\bar{y}^{*}_{[0,\tau_{\ell}-1]}(t_{\ell})\|^{2}.$
(27)

Since
$\gamma={\min\\{\underline{\lambda}_{Q},\underline{\lambda}_{R}\\}}/{\max\\{c_{1},2c_{2}\\}}$,
we have that

$\displaystyle\quad\gamma(c_{1}\|\bar{u}^{*}_{[0,\tau_{\ell}-1]}(t_{\ell})\|^{2}\\!+\\!2c_{2}\|\bar{y}^{*}_{[0,\tau_{\ell}-1]}(t_{\ell})\|^{2})$
$\displaystyle\leq\sum_{i=0}^{\tau_{\ell}-1}\big{(}\|\bar{u}_{i}^{*}(t_{\ell})\|_{R}^{2}+\|\bar{y}_{i}^{*}(t_{\ell})\|_{Q}^{2}\big{)}.$
(28)

Substituting (14) and (III-B) into (III-B), we arrive at

$\displaystyle\mathcal{T}_{2}$
$\displaystyle\leq-\frac{\gamma}{2}\|\xi_{[t_{\ell},t_{\ell+1}-1]}\|^{2}+2\bar{\lambda}_{Q}\sum_{i=0}^{\tau_{\ell}-1}\|h^{*}_{i}(t_{\ell})\|^{2}$
$\displaystyle\quad+\sum_{i=0}^{\tau_{\ell}-1}\big{(}\|\bar{u}_{i}^{*}(t_{\ell})\|_{R}^{2}+\|\bar{y}_{i}^{*}(t_{\ell})\|_{Q}^{2}\big{)}$
$\displaystyle\quad+\bar{\lambda}_{Q}\big{(}2\eta\bar{n}^{2}+2\|h^{*}_{[-\eta,-1]}(t_{\ell})\|^{2}\big{)}\sum_{i=0}^{\tau_{\ell}-1}(\rho^{i+\eta})^{2}.$
(29)

Combining the results in Parts a and b, it holds that

$\displaystyle
V(\xi_{t_{\ell+1}})-V(\xi_{t_{\ell}})\leq-\frac{\gamma}{2}\|\xi_{[t_{\ell},t_{\ell+1}-1]}\|^{2}+2\bigg{[}\bar{\lambda}_{Q}+2\bar{\lambda}_{P_{r}}$
$\displaystyle+\lambda_{g}\bar{n}\|H_{u\xi}^{\ddagger}\|^{2}\Big{(}1+\frac{\bar{\lambda}_{P_{r}}}{\underline{\lambda}_{R}}\Big{)}\bigg{]}\bigg{[}\big{(}\eta\bar{n}^{2}+\|h^{*}_{[-\eta,-1]}(t_{\ell})\|^{2}\big{)}$
$\displaystyle\times\sum_{i=0}^{\tau_{\ell}-1}(\rho^{\eta+i})^{2}+\sum_{i=0}^{\tau_{\ell}-1}\|h^{*}_{i}(t_{\ell})\|^{2}\bigg{]}+\lambda_{h}\eta\bar{n}^{2}+2\bigg{[}\bar{\lambda}_{P_{r}}$
$\displaystyle+\lambda_{g}\bar{n}\|H_{u\xi}^{\ddagger}\|^{2}\Big{(}1+\frac{\bar{\lambda}_{P_{r}}}{\underline{\lambda}_{R}}\Big{)}\bigg{]}\|\bar{\xi}^{*}_{\tau_{\ell}}(t_{\ell})\|^{2}.$
(30)

In addition, (17) indicates that $\tau_{\ell}\leq\check{\tau}_{\ell}$, and
hence (III-B) obeys

$\displaystyle\quad V(\xi_{t_{\ell+1}})-V(\xi_{t_{\ell}})$
$\displaystyle\leq-\Big{(}\frac{\gamma}{2}\\!-\\!\sigma\Big{)}\|\xi_{t_{\ell}}\|^{2}+2\sigma\eta\bar{n}^{2}\\!+\\!\\!\\!\sum_{i=1}^{\tau_{\ell}-1}\\!\\!\Big{(}\sigma\|\bar{\xi}_{i}(t_{\ell})\|^{2}\\!-\\!\frac{\gamma}{2}\|\xi_{t_{\ell}+i}\|^{2}\Big{)}$
$\displaystyle\leq-\Big{(}\frac{\gamma}{2}-\sigma\Big{)}\|\xi_{t_{\ell}}\|^{2}+2\sigma\eta\bar{n}^{2}+\frac{\gamma}{2}\sum_{i=1}^{\tau_{\ell}-1}\|\bar{\xi}_{e}(t_{\ell}+i)\|^{2}$
$\displaystyle\leq-(\gamma/2-2\sigma)\|\xi_{t_{\ell}}\|^{2}+\alpha_{1}(\bar{n})$

where
$\alpha_{1}(\bar{n}):=\bar{n}[J^{*}_{L}(t_{\ell})/\lambda_{h}+\sum_{i=1}^{L-1}(\rho^{\eta+i})(\eta+J^{*}_{L}(t_{\ell})/\lambda_{h})]+2\sigma\eta\bar{n}^{2}$
is a $\mathcal{K}_{\infty}$ function, and $\bar{\sigma}\leq\gamma/4$.
Moreover, according to [15, Lemma 1], the Lyapunov function satisfies
$\underline{\lambda}_{P_{\xi}}\|\xi_{t_{\ell}}\|^{2}\leq
V(\xi_{t_{\ell}})\leq\gamma_{1}\|\xi_{t_{\ell}}\|^{2}+\alpha_{2}(\bar{n})$ for
every $\ell\in\mathbb{N}_{0}$ and some constant $\gamma_{1}>0$ with function
$\alpha_{2}(\cdot)\in\mathcal{K}_{\infty}$. Hence, we have that

$\displaystyle V(\xi_{t_{\ell+1}})$
$\displaystyle\leq\underbrace{\Big{(}1-\frac{\gamma-4\sigma}{2\gamma_{1}}\Big{)}}_{\stackrel{{\scriptstyle\triangle}}{{=}}\gamma_{2}<1}V(\xi_{t_{\ell}})+\underbrace{\frac{\gamma-4\sigma}{2\gamma_{1}}\alpha_{2}(\bar{n})+\alpha_{1}(\bar{n})}_{\stackrel{{\scriptstyle\triangle}}{{=}}\alpha_{3}(\bar{n})\in\mathcal{K}_{\infty}}.$

In conclusion, the Lyapunov function decreases at every self-triggered time.
Since problem (10) is feasible at $t_{0}$, it is feasible at $t_{\ell}$ for
all $\ell\in\mathbb{N}_{0}$. ∎

######  Remark 2.

Condition (18) ensures that problem (10) is feasible at $t_{\ell+1}$ if it is
feasible at $t_{\ell}$, which is also known as _recursive feasibility_. In
addition, condition (19) guarantees the Lyapunov function of system (3)
decreases along the self-triggered times. Therefore, the self-triggering law
reduces the communication load while ensuring the recursive feasibility of
(10). Similar to model-based self-triggered control, the choice of the
contraction factor $\sigma\in(0,1)$ determines the trade-off between having a
better control performance and incurring less transmissions. The smaller
$\sigma$, the better the system performance.

Merging the data-driven MPC and self-triggering schemes in Lemma 3, our data-
driven self-triggered controller is presented in Alg. 1, with stability
guarantees provided below.

Algorithm 1 Data-driven Self-triggered Predictive Control.

1: Input: Prediction window $L\geq{\eta}+1$; coefficients $R\succ 0$, $Q\succ
0$, $P\succ 0$, $\lambda_{g}>0$, $\lambda_{h}>0$, $\sigma\in(0,1)$,
$\epsilon>0$ and $r>0$ satisfying (16); noise bound $\bar{n}$; data
$\\{u^{p}_{[0,N-1]},y^{p}_{[0,N-1]}\\}$ of (3) from initial condition
$x_{0}^{p}$, with $u^{p}_{[0,N-1]}$ persistently exciting of order $L+\eta+1$.

2: Construct Hankel matrices for the input, output, and the extended state
trajectories, i.e., $H_{L+\eta}(u^{p})$, $H_{L+\eta}(y^{p})$, and
$H_{u\xi}:=[H_{L+\eta}^{\top}(u^{p})~{}H_{1}^{\top}(\xi^{p}_{N-L-\eta+1})]^{\top}$.

3: Compute $\rho^{i}$ for $i\in\\{1,2,\ldots,L-1\\}$ from (15).

4: If $t=t_{\ell}$, do

5:

  * Use the past $\eta$ measurements, i.e., $u_{[t-\eta,t-1]}$ and $\zeta_{[t-\eta,t-1]}$, to solve problem (10).

Determine the next self-triggered time $t_{\ell+1}$ such that the inter-
triggering time obeys (17)–(19).

Set $u_{t}=\bar{u}^{*}_{0}(t)$.

Set $t_{\ell}=t_{\ell+1}$

6: Else if $t\neq t_{\ell}$

7: Set $u_{t}=\bar{u}^{*}_{t-t_{\ell}}(t_{\ell})$.

8: End if

9: Set $t=t+1$ and go back to 4.

######  Theorem 1.

Let Assumptions 1—5 hold. Consider system (3) whose control input is generated
by solving problem (10) at each self-triggered time $t_{\ell}$. If the
sequence $\\{t_{\ell}\\}_{\ell\in\mathbb{N}_{0}}$ satisfies the conditions in
Lemma 3, then system (3) achieves practical exponential stability under the
controller (10).

###### Proof.

Recall from Lemma 3 that it holds for any $t_{\ell}$

$V(\xi_{t_{\ell+1}})\leq\gamma_{2}V(\xi_{t_{\ell}})+\alpha_{3}(\bar{n}).$ (31)

Recursively, it can be deduced that

$V(\xi_{t_{\ell}})\leq\gamma_{2}^{\ell}V(\xi_{t_{0}})+\alpha_{3}(\bar{n})\sum_{i=0}^{\ell-1}\gamma_{2}^{i}.$
(32)

In addition, noticing that
$\underline{\lambda}_{P_{\xi}}\|\xi_{t_{\ell}}\|^{2}\leq
V(\xi_{t_{\ell}})\leq\gamma_{1}\|\xi_{t_{\ell}}\|^{2}+\alpha_{2}(\bar{n})$, we
arrive at

$\|\xi_{t_{\ell}}\|^{2}\leq\frac{\gamma_{1}}{\underline{\lambda}_{P_{\xi}}}\gamma_{2}^{\ell}\|\xi_{t_{0}}\|^{2}+\Big{(}\alpha_{2}(\bar{n})\gamma_{2}^{\ell}+\alpha_{3}(\bar{n})\sum_{i=0}^{\ell-1}\gamma_{2}^{i}\Big{)}$

which verifies the system’s practical exponential stability. ∎

##  IV Data-driven Self-triggered Predictive Control for State Feedback
Systems

In this section, we consider a special case of the system in (3) having $C=I$
and $D=0$; that is, when input-state data become available. In this setup, a
state feedback controller is employed. Although simpler than the MPC scheme,
the state feedback controller does not provide a predicted system trajectory,
which makes design of the self-triggering mechanism difficult. To address this
issue, a data-based norm maximization problem is put forth such that a future
system trajectory can be predicted from the input-state data. The key idea
behind constructing a predicted trajectory for our self-triggered state
feedback controller is combining the results of Lemmas 1 and 2. Therefore, the
trajectory prediction along with the self-triggering mechanism here builds on
the idea of (10).

Consider the following state feedback controller

$\displaystyle x_{t+1}$ $\displaystyle=Ax_{t}+Bu_{t},\quad$ $\displaystyle t$
$\displaystyle\in\mathbb{N}_{0}$ (33a) $\displaystyle u_{t}$
$\displaystyle=K{\zeta}_{t_{\ell}},\quad$ $\displaystyle t_{\ell}$
$\displaystyle\leq t<t_{\ell+1}$ (33b) $\displaystyle\zeta_{t_{\ell}}$
$\displaystyle=x_{t_{\ell}}+n_{t_{\ell}}$ (33c)

where the gain matrix $K\in\mathbb{R}^{n_{u}\times n_{x}}$ can be obtained via
a data-based linear matrix inequality method in e.g. [2]. At each self-
triggered time $t_{\ell}$, the state $x_{t_{\ell}}$ is sampled and sent to the
controller. Similar as in Section III, the transmission is corrupted by noise
$n_{t}$, and hence the controller receives $\zeta_{t_{\ell}}$. The control
input $u_{t_{\ell}}$ is sent to the plant, and kept the same within the
interval $[t_{\ell},t_{\ell+1})$. See Fig. 2 for an illustration.

Figure 2: Data-driven self-triggered state feedback controller for (33).

Instead of input-output data in Section III, input-state data are available in
this section, and hence Assumptions 2 and 3 are replaced as follows.

######  Assumption 6 (_Unknown system model_).

The system matrices $(A,B)$ in (33) are unknown, but some pre-collected input-
state data $\\{u^{p}_{t},x^{p}_{t}\\}_{t=0}^{N-1}$ are available.

######  Assumption 7 (_Pre-collected data_).

Let $\\{u^{p}_{t}\\}_{t=0}^{N-1}$ be any sequence persistently exciting of
order $L+n_{x}+1$, where $L\geq 2$. The state sequence
$\\{x^{p}_{t}\\}_{i=0}^{N-1}$ is generated by system (33) offline using the
input sequence $\\{u^{p}_{t}\\}_{t=0}^{N-1}$, from any initial condition
$x^{p}_{0}$.

Different from the (data-driven) MPC, (data-driven) state-feedback control
does not produce a system trajectory. Hence, we resort to the following
problem to predict the states for ensuing $L$ times at each triggered time
$t_{\ell}$

$\displaystyle
J^{*}_{L}(x_{t_{\ell}}):=\underset{g(t_{\ell}),h(t_{\ell})\atop\bar{x}_{i}(t_{\ell})}{\max}$
$\displaystyle\sum_{i=0}^{L-1}\|\bar{x}_{i}(t_{\ell})\|_{\infty}$ (34a)
$\displaystyle{\rm s.t.}\quad$
$\displaystyle\left[\begin{matrix}u(t_{\ell})\\\
\bar{x}(t_{\ell})\end{matrix}\right]=\left[\begin{matrix}H_{L}(u^{p})\\\
H_{L}(x^{p})\end{matrix}\right]g(t_{\ell})$ (34b)
$\displaystyle\,\bar{x}_{0}(t_{\ell})+h(t_{\ell})=\zeta_{t_{\ell}}$ (34c)
$\displaystyle\,\|h(t_{\ell})\|\leq\bar{n}$ (34d)

where
$\bar{x}(t_{\ell})=[\bar{x}_{0}^{\top}(t_{\ell})~{}\cdots~{}\bar{x}^{\top}_{L-1}(t_{\ell})]^{\top}\in\mathbb{R}^{n_{x}L}$
collects the predicted states for times from $t_{\ell}$ to $t_{\ell}+L-1$,
$u(t_{\ell})=[u_{t_{\ell}}^{\top}~{}\cdots~{}u_{t_{\ell}}^{\top}]^{\top}\in\mathbb{R}^{n_{u}L}$
repeats the control input $u_{t_{\ell}}$ for all times in
$[t_{\ell},t_{\ell}+L-1]$, and $h(t_{\ell})\in\mathbb{R}^{n_{x}}$ is a slack
variable mitigating the noise $n_{t}$. Upon solving (34), the optimal state
$\bar{x}^{*}(t_{\ell})$ predicts the worst (in terms of norm) state trajectory
under noise $n_{t}$ and control inputs $u(t_{\ell})$, which is subsequently
used to characterize the actual future state following the same procedure as
in Section III. Specifically, capitalizing on $\bar{x}^{*}(t_{\ell})$, the
self-triggering mechanism is constructed as follows

$\left\\{\begin{aligned}
t_{\ell+1}&:=t_{\ell}+\min{\\{L-1,\tau_{\ell}\\}},~{}~{}t\in\mathbb{N}_{0};~{}~{}t_{0}:=0\\\
\tau_{\ell}&:=\min\\{\tau\in\mathbb{N}:\varphi(\bar{x}_{\tau}(t_{\ell}),\,\zeta_{t_{\ell}},\,\rho^{\tau}\\!,\,\bar{n},\,\tau)>\sigma\|x_{t_{\ell}}\|_{\infty}\\}\end{aligned}\right.$
(35)

where $\sigma\in(0,1)$ is a threshold parameter,
$\rho^{\tau_{\ell}}:=\min\\{\rho|\rho>=A^{\tau_{\ell}}\\}$ for
$\tau=1,\cdots,L-1$ is defined in Lemma 2 with $C=I$, and
$\varphi:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{x}}\times\mathbb{R}\times\mathbb{R}\times\mathbb{N}\rightarrow\mathbb{R}_{>0}$
is a certain function to quantify the error induced by the self-triggered
sampling. The error between the states at the current triggered time
$t_{\ell}$ and its next one satisfies

$\|\zeta_{t_{\ell}}-x_{t_{\ell+1}}\|_{\infty}\leq\|x_{t_{\ell}}-\bar{x}^{*}_{\tau_{\ell}}(t_{\ell})\|_{\infty}+\|x_{t_{\ell+1}}-\bar{x}^{*}_{\tau_{\ell}}(t_{\ell})\|_{\infty}+\bar{n}.$
(36)

It follows from Lemma 2 that the error between the predicted state and the
actual state can be bounded by

$\displaystyle\|x_{t_{\ell+1}}-\bar{x}^{*}_{\tau_{\ell}}(t_{\ell})\|_{\infty}\leq\rho^{\tau_{\ell}}(\bar{n}+\|h^{*}(t_{\ell})\|_{\infty}).$
(37)

Substituting (37) into (36), the self-triggering function $\varphi(\cdot)$ in
(35) can be given by

$\displaystyle\varphi(\bar{x}(t_{\ell}),$
$\displaystyle~{}h(t_{\ell}),\zeta_{t_{\ell}},\rho^{\tau},\bar{n},\tau)$ (38)
$\displaystyle:=\|x_{t_{\ell}}-\bar{x}_{\tau}(t_{\ell})\|_{\infty}+\rho^{\tau_{\ell}}(\bar{n}+\|h(t_{\ell})\|_{\infty})+\bar{n}.$

Based on (35), (37), and (38), the following stability result comes ready
whose proof is similar to that of the model-based self-triggered control in
[38] and is thus omitted here.

######  Lemma 4.

Let Assumptions 1, 6, and 7 hold. Consider the state feedback system in (33).
If the self-triggered times satisfy (35) with function $\varphi(\cdot)$
defined by (38), and the parameters $h^{*}(t_{\ell})$, $\bar{x}^{*}(t_{\ell})$
obtained by solving problem (34), then system (33) achieves input-to-state
stability under the proposed self-triggering mechanism.

######  Remark 3 (_Novelty_).

The closest references to the results in this section are [8] and [10], yet
they distinguish from the present work considerably. First, the work [8] is
devoted to estimating the maximum sampling interval for state feedback systems
under different levels of process noise and controller gains. Therefore, they
only use the pre-collected data to estimate upper bounds of $\|A^{i}\|$ for
all $i\in\\{1,\cdots,\bar{h}\\}$. This plays a similar role as the norm
maximization problem (15). However, the key idea of our self-triggering
mechanism is to construct a future state trajectory using the most recently
received state, and (15) offers a way for estimating one of the parameters
required for constructing the trajectory. In addition, although the work [10]
considers the self-triggered control problem by generalizing the estimation
method in [8], its system model as well as associated self-triggered
controller design are completely different from ours.

##  V Numerical Examples

To validate the effectiveness of our proposed data-driven self-triggering
mechanism as well as predictive controllers, several numerical examples are
provided in this section.

###  V-A Output feedback system

Consider a linearized version of the four-tank system in, e.g., [39, 15]. The
system matrices are given by

$\displaystyle A=\left[\begin{matrix}0.927&0&0.041&0\\\ 0&0.918&0&0.033\\\
0&0&0.924&0\\\
0&0&0&0.937\end{matrix}\right],B=\left[\begin{matrix}0.017&0.001\\\
0.001&0.023\\\ 0&0.061\\\ 0.072&0\end{matrix}\right]$ $\displaystyle
C=\left[\begin{matrix}1&0&0&0\\\ 0&1&0&0\end{matrix}\right],\quad D=0.$

with the observability index $\eta=2$. To this end, an input-output trajectory
of length $N=800$ was obtained by means of simulating the open-loop system
off-line using persistently exciting input sequence. Let the network-induced
noise $n_{t}$ satisfies $\bar{n}=0.0015$. Parameters of Alg. 1 were set as
follows: the prediction horizon $L=11$, cost matrices $Q=I$, $R=8\cdot
10^{-3}I$, coefficients $\lambda_{g}\bar{n}=10^{-6}$,
$\lambda_{h}/\bar{n}=500$ the input constraint set $\mathbb{U}=[-2,2]^{2}$,
and the input-output equilibrium
$(u^{e},y^{e})=([1,1]^{\top},[0.65,0.77])^{\top}$. Over a simulation horizon
of $200$ time steps, Fig. 3 depicts the output trajectory of system using
data-driven self-triggering MPC scheme with $\sigma=0.88$. Evidently, the
system converges to the setpoint with $37$ measurement packets transmitted to
the controller. In addition, since $\eta=2$, at most $2\times 37=74$ outputs
are sent based on the packetized transmission protocol in Section II-A, which
is also much lesser than $200$, and illustrates the effectiveness of the self-
triggering mechanism. Figs. 4 and 5 compare the cost function
$J_{L}^{*}(t_{\ell})$ and the self-triggered times for different $\sigma$
values.

Figure 3: Output trajectory (top) and inter-sampling time
$t_{\ell+1}-t_{\ell}$ (bottom) with $\sigma=0.88$.
Figure 4: Optimal cost $J^{*}_{L}$ with different $\sigma$.
Figure 5: Inter-sampling time $t_{\ell+1}-t_{\ell}$ with different $\sigma$.

###  V-B State feedback systems

In this part, we consider state feedback systems and illustrate the
effectiveness of the self-triggering mechanism in Section IV. To guarantee the
system performance as done in the model-based event-triggered case e.g., [40],
the following $\mathcal{L}_{2}$-like condition is incorporated and examined

$\kappa\sum_{i=1}^{\tau_{\ell}}\|\bar{x}^{*}_{i}(t_{\ell})\|\leq\mu\tau_{\ell}\bar{n}$
(39)

where constants $\kappa>0$ and $\mu>0$ balance between achieving a better
system performance and incurring less transmissions.

####  V-B1 Example 1

Consider the system [10]

$\dot{x}(t)=\left[\begin{matrix}0&1\\\
0&-0.1\end{matrix}\right]x(t)+\left[\begin{matrix}0\\\
0.1\end{matrix}\right]u(t),~{}~{}t\geq 0.$ (40)

Here, we consider a discrete-time version of this system with a sampling
period of $0.1s$. Let $\sigma=0.27$, $\kappa=0.1$, $\mu=200$, and use the
initial state $x_{0}=[3~{}-2]^{\top}$. System performance is depicted in Fig.
6, where only $14$ out of $200$ samples are needed. This outperforms the
result in the previous work [10], which uses $15$ measurements under noise
bounded by $10^{-4}$.

Figure 6: State trajectory of the order-$2$ system under (34) with
$\sigma=0.27$.

####  V-B2 Example 2

Consider the inverted pendulum control problem in [40]. The linearized system
is given by

$\displaystyle\dot{x}(t)=\left[\begin{matrix}0&1&0&0\\\
0&0&\frac{m_{1}g}{-m_{2}}&0\\\ 0&0&0&1\\\
0&0&\frac{g}{\ell}&0\end{matrix}\right]x(t)+\left[\begin{matrix}0\\\
\frac{1}{m_{2}}\\\ 0\\\ \frac{-1}{m_{2}\ell}\end{matrix}\right]u(t)$ (41)

where $m_{1}=1$, $m_{2}=10$, $\ell=3$, and $g=10$. Let $\sigma=0.27$,
$\kappa=0.1$, $\mu=200$, and the initial state
$x_{0}=[0.98~{}0~{}0.2~{}0]^{\top}$. It can be seen from Fig. 7 that the
system converges to zero with $62$ samples over $[0,200]$, which is less than
that of [10].

Figure 7: State trajectory of the order-$4$ system under (34) with
$\sigma=0.3$.

##  VI Conclusions

This paper puts forth a data-driven self-triggering control framework for
unknown linear systems through trajectory prediction. For output feedback
systems, a data-driven MPC scheme is developed, which generates a sequence of
control inputs once an event (transmission) is triggered. The predicted output
trajectory from the MPC is further used to design a self-triggering law that
is purely data-based. Both feasibility as well as practical exponential
stability were established for the resulting data-driven self-triggered MPC.
Moreover, when a state feedback control law is considered, a norm maximization
problem was developed to predict future system states, thus enabling the data-
driven self-triggered control. Finally, numerical examples are presented to
validate the practical merits of the proposed data-driven control methods and
theory.

## References

  * [1] Z. Hou and Z. Wang, “From model-based control to data-driven control: Survey, classification and perspective,” _IEEE Trans. Ind. Inform._ , vol. 235, pp. 3–35, Jun. 2013.
  * [2] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” _IEEE Trans. Autom. Control_ , vol. 65, no. 3, pp. 909–924, Mar. 2020.
  * [3] K. J. Åström and B. Wittenmark, “On self-tuning regulators,” _Automatica_ , vol. 9, no. 2, pp. 185–199, Mar. 1973.
  * [4] J. Chen, J. Sun, and G. Wang, “From unmanned systems to autonomous intelligent systems,” _Eng._ , vol. 8, no. 8, pp. 2–6, 2022.
  * [5] G. Liu, “Coordination of networked nonlinear multi-agents using a high-order fully actuated predictive control strategy,” _IEEE/CAA J. Autom. Sinica_ , Jan. 2022, doi: 10.1109/JAS.2022.105449.
  * [6] J. C. Willems, I. Markovsky, P. Rapisarda, and B. L. M. De Moor, “A note on persistency of excitation,” _Syst. Control Lett._ , vol. 56, no. 4, pp. 325–329, May, 2005.
  * [7] H. J. van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: A new perspective on data-driven analysis and control,” _IEEE Trans. Autom. Control_ , vol. 65, no. 11, pp. 4753–4768, Jan. 2020.
  * [8] S. Wildhagen, J. Berberich, M. Hertneck, and F. Allgöwer, “Data-driven estimation of the maximum sampling interval: analysis and controller design for discrete-time systems,” _arXiv:2101.00903_ , Jan. 2021\.
  * [9] J. G. Rueda-Escobedo, E. Fridman, and J. Schiffer, “Data-driven control for linear discrete-time delay systems,” _IEEE Trans. Autom. Control_ , pp. 1–16, Jul. 2021, doi: 10.1109/TAC.2021.3096896.
  * [10] X. Wang, J. Berberich, J. Sun, G. Wang, F. Allgöwer, and J. Chen, “Data-driven control of event- and self-triggered discrete-time systems,” _arXiv:2202.08019_ , 2021.
  * [11] G. Bianchin, M. Vaquero, J. Cortes, and E. Dall’Anese, “Online stochastic optimization for unknown linear systems: Data-driven synthesis and controller analysis,” _arXiv:2108.13040_ , Aug. 2021.
  * [12] C. De Persis and P. Tesi, “Low-complexity learning of linear quadratic regulators from noisy data,” _Automatica_ , vol. 128, p. 109548, Jun. 2021\.
  * [13] A. Bisoffi, C. De Persis, and P. Tesi, “Trade-offs in learning controllers from noisy data,” _Syst. Control Lett._ , vol. 154, p. 104985, Aug. 2021.
  * [14] F. Zhao, X. Li, and K. You, “Data-driven control of unknown linear systems via quantized feedback,” _Proc. Mach. Learn. Res._ , vol. 144, pp. 1–12, 2022\.
  * [15] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Data-driven model predictive control with stability and robustness guarantees,” _IEEE Trans. Autom. Control_ , vol. 66, no. 4, pp. 1702–1717, Jun. 2021.
  * [16] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in _Proc. Eur. Control Conf._ , Naples, Italy, Jun. 25-28, 2019, pp. 307–312.
  * [17] F. Dörfler, J. Coulson, and I. Markovsky, “Bridging direct & indirect data-driven control formulations via regularizations and relaxations,” _IEEE Trans. Autom. Control_ , pp. 1–15, Feb. 2022, doi: 10.1109/TAC.2022.3148374.
  * [18] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Robust constraint satisfaction in data-driven MPC,” in _Proc. IEEE Conf. Decis. Control_ , Jeju, Korea, Dec. 14-18, 2020, pp. 1260–1267.
  * [19] W. Liu, J. Sun, G. Wang, F. Bullo, and J. Chen, “Data-driven resilient predictive control under Denial-of-Service,” _arXiv:2110.12766_ , Oct. 2021\.
  * [20] G. Baggio, D. S. Bassett, and F. Pasqualetti, “Data-driven control of complex networks,” _Nat. Commun._ , no. 1429, pp. 1–13, Mar. 2021.
  * [21] F. Pasqualetti, F. Dörfler, and F. Bullo, “Attack detection and identification in cyber-physical systems,” _IEEE Trans. Autom. Control_ , vol. 58, no. 11, pp. 2715–2729, Jun. 2013.
  * [22] Z. Pang, C. Bai, G. Liu, Q. Han, and X. Zhang, “A novel networked predictive control method for systems with random communication constraints,” _J. Syst. Sci. Complex._ , vol. 34, no. 4, pp. 1364–1378, Aug. 2021.
  * [23] W. P. M. H. Heemels, K. H. Johansson, and P. Tabuada, “An introduction to event-triggered and self-triggered control,” in _Proc. IEEE Conf. Decis. Control_ , Maui, HI, USA, Dec. 10-13 2012, pp. 3270–3285.
  * [24] M. Lv, D. Wang, Z. Peng, L. Lu, and H. Wang, “Event-triggered neural network control of autonomous surface vehicles over wireless network,” _Sci. China Inf. Sci._ , vol. 63, Mar. 2020, doi: 10.1007/s11432-019-2679-5.
  * [25] Z. Chen, Q. L. Han, Y. Yan, and Z. G. Wu, “How often should one update control and estimation: Review of networked triggering techniques,” _Sci. China Inf. Sci._ , vol. 63, no. 5, pp. 1–18, Feb. 2020.
  * [26] D. Antunes and W. P. M. H. Heemels, “Rollout event-triggered control: Beyond periodic control performance,” _IEEE Trans. Autom. Control_ , vol. 59, no. 12, pp. 3296–3311, Aug. 2014.
  * [27] M. Wakaiki and H. Sano, “Event-triggered control of infinite-dimensional systems,” _SIAM J. Control Optim._ , vol. 58, no. 2, pp. 605–635, Mar. 2020\.
  * [28] W. P. M. H. Heemels and M. C. F. Donkers, “Model-based periodic event-triggered control for linear systems,” _Automatica_ , vol. 49, no. 3, pp. 698–711, Mar. 2013.
  * [29] A. Anta and P. Tabuada, “To sample or not to sample: Self-triggered control for nonlinear systems,” _IEEE Trans. Autom. Control_ , vol. 55, no. 9, pp. 2030–2042, Sept. 2010.
  * [30] T. M. P. Gommans and W. P. M. H. Heemels, “Resource-aware MPC for constrained nonlinear systems: A self-triggered control approach,” _Syst. Control Lett._ , vol. 79, pp. 59–67, May, 2015.
  * [31] Z. Sun, L. Dai, K. Liu, D. V. Dimarogonas, and Y. Xia, “Robust self-triggered MPC with adaptive prediction horizon for perturbed nonlinear systems,” _IEEE Trans. Autom. Control_ , vol. 64, no. 11, pp. 4780–4787, Mar. 2019\.
  * [32] J. Almeida, C. Silvestre, and A. M. Pascoal, “Self-triggered output feedback control of linear plants in the presence of unknown disturbances,” _IEEE Trans. Autom. Control_ , vol. 59, no. 11, pp. 3040–3045, Nov. 2014\.
  * [33] J. O’Reilly, _Mathematics in Science and Engineering, Observers for Linear Systems_. London: Academic Press, 1983\.
  * [34] F. M. Callier and C. A. Desoer, _Linear System Theory_. Springer New York, 1991.
  * [35] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “On the design of terminal ingredients for data-driven MPC,” _arXiv:2101.05573_ , May, 2021.
  * [36] J. B. Rawlings, D. Q. Mayne, and M. M. Diehl, _Model Predictive Control: Theory and Design_ , 2nd ed. Nob Hill Publishing, 2019.
  * [37] C. Cai and A. R. Teel, “Input–output-to-state stability for discrete-time systems,” _Automatica_ , vol. 44, no. 2, pp. 326–336, Sept. 2008.
  * [38] M. Wakaiki, “Self-triggered stabilization of discrete-time linear systems with quantized state measurements,” _arXiv:2103.15362_ , Jul. 2021.
  * [39] T. Raff, S. Huber, Z. K. Nagy, and F. Allgöwer, “Nonlinear model predictive control of a four tank system: An experimental stability study,” in _Proc. IEEE Int. Conf. Control Appl._ , Munich, Germany, Oct. 4-6, 2006, pp. 237–242.
  * [40] C. Peng and Q.-L. Han, “A novel event-triggered transmission scheme and $\mathcal{L}_{2}$ control co-design for sampled-data control systems,” _IEEE Trans. Autom. Control_ , vol. 58, no. 10, pp. 2620–2626, Mar. 2013\.