# GRM: Gradient Rectification Module for Visual Place Retrieval

Boshu Lei1 Wenjie Ding2 Limeng Qiao2 Xi Qiu2 *This work was done in MEGVII
Inc. 1Boshu Lei is with Faculty of Automation, Xi’an Jiaotong Univerisity, Xi’
an, China2 Wenjie Ding, Limeng Qiao and Xi Qiu are with MEGVII Inc.

###### Abstract

Visual place retrieval aims to search images in the database that depict
similar places as the query image. However, global descriptors encoded by the
network usually fall into a low dimensional principal space, which is harmful
to the retrieval performance. We first analyze the cause of this phenomenon,
pointing out that it is due to degraded distribution of the gradients of
descriptors. Then, we propose Gradient Rectification Module (GRM) to alleviate
this issue. GRM is appended after the final pooling layer and can rectify
gradients to the complementary space of the principal space. With GRM, the
network is encouraged to generate descriptors more uniformly in the whole
space. At last, we conduct experiments on multiple datasets and generalize our
method to classification task under prototype learning framework.

##  I INTRODUCTION

Visual place retrieval (VPR) has gained much popularity in recent years. Given
a query image, VPR aims to find images in the database that depict similar
places. The typical pipeline of this task is to perform rough matching using
global features and refine similarity using local features [1, 2, 3]. Most of
the methods aim to obtain distinctive global and local features by improving
the optimization function [4], aggregation methods [5, 6] and local branch
structures [1, 2]. However, few methods consider the feature distribution,
which is very important for VPR task because it directly affects the retrieval
results.

Early work to address the significance of feature distribution for retrieval
is [7]. The author claims that whitening can balance out the co-occurrence of
features, which is beneficial to retrieval. Later, many works try to eliminate
the unbalanced distribution of descriptors, like loss regularizor [8, 9] and
normalization technique[10, 11, 12]. However, no one goes deeper into the
question, what makes the model generate such unbalanced features. After
thorough analysis, we attribute it to the gradients of descriptors in training
process. During training, we find that both the descriptors and the
corresponding gradients are trapped in the same low dimensional principal
space. As updates to the descriptors, gradients trapped in principal space
won’t push descriptors to other dimensions, resulting in the network
generating degraded features and thus harmful to the retrieval performance. We
term the span of the eigenvectors with large eigenvalues from feature
covariance matrix as principal space and the span of remaining vectors as
complementary space.

Figure 1: Gradient Rectification. The circles in the figure represent the
distribution of descriptors and gradients respectively. Without GRM, network
generates degraded features when gradients fall into the same principal space
of descriptors (Up). When GRM is applied, gradients are pulled out of this
space (in orthogonal direction) and the network will generate descriptors more
uniformly (Down).

We now propose to modify the distribution of descriptors through a Gradient
Rectification Module (GRM) as illustrated in Fig. 3. First, we propose a new
memory queue method to estimate the co-variance matrix of descriptors. This
covariance matrix is used to generate a projection matrix to rectify
gradients. Second, we propose our rectification criterion. We inhibit the
component of gradients along the principal space and stimulate the component
along directions orthogonal to the principal space. After doing so, gradients
will be deviated from the principal space and descriptors are pulled to the
complementary space of the principal space. In the forward pass, GRM stores
descriptors in the memory queue for covariance estimation. In the backward
pass, GRM modifies gradients based on the covariance matrix estimated from
queue. It’s worth noting that GRM is plug-and-play and parallel to those
methods that aim to improve the loss function or local branches.

In conclusion, our contributions are threefold:

1) We point out that gradients fall into the same space of descriptors in VPR
task, hindering the network from exploring other dimensions of the feature
space.

2) We propose a simple yet effective Gradient Rectification Module (GRM). By
modifying the gradients of descriptors, GRM pulls descriptors to the
complement space of the principal space and alleviates the problem above.

3) We conduct extensive experiments in various datasets in VPR and
classification tasks. Experimental results show that the proposed method
significantly outperforms state of the arts.

##  II Related Works

###  II-A Global Feature Descriptor

Early approaches to extract global feature descriptor are based on aggregation
of hand crafted local features. Methods are like Bag of Words (BoW) [13] and
Vector of Locally Aggregated Descriptors (VLAD) [14]. Sparse keypoint location
or dense sampling on image are applied to aggregate local features. When deep
learning technology is introduced into these frameworks, new methods like
NetVLAD [5], NetBoW [15], NetFV [6] emerge. Other methods apply simple pooling
layer, eg. max pooling or GeM pooling [16], on top of the CNN architecture
[12, 4]. These methods are trained with Triplet Loss [5], or ranking-based
loss [17, 18]. Most recently, [4] proposes to use a continuous measure to
calculate the similarity between images. Instead of depending on the binary
relation, they use the calculated similarity score to train the neural network
and achieve huge success in MSLS challenge [19].

Figure 3: The learning paradigm of the proposed method.  Forward descriptors
are accumulated in the memory queue and later used to generate the projection
matrix. Backward gradients are multiplied by this matrix when flowing through
the GRM.

###  II-B Local Feature Refinement

Early methods for refinement are separated from global feature matching
scheme. Local features like SIFT[20], SURF[21] are extracted and matched based
on nearest neighbor criterion. Then model-based methods, such as RANSAC[22]
and PROSAC[23] are used to generate transformation hypotheses based on
correspondences. The final retrieval result is re-ranked by the number of
inliers under the hypothesis. DELF[3] proposes to use the feature map before
pooling for later refinement process. PatchNetVLAD[1] proposes to extract
NetVLAD[5] descriptors on image patches for later refinement procedure.[2] is
the first to jointly learn local and global features for retrieval task.

###  II-C Whitening and Spreading Methods

According to [24], when the features become white, the Fisher information
matrix reduces to an identity matrix. In that case, the simple stochastic
gradient descent optimization coincides with the natural gradients descent
optimization which is more suitable for difficult optimzation tasks. Whitening
methods can be classified into two categories. The first category is to add
extra penalty terms in loss function. Zhang et al. [8] proposed to add a
regularizor to spread vectors uniformly in the feature space. In [9], the
author provides a thorough analysis of the uniformity and alignment metric of
the feature distribution space. Additional loss penalties are added to adjust
the uniformity and alignment of feature space. Another method is to process
the features directly, either in the forward pass or the backward pass. Batch
Normalization (BN) [10] is the most widely used method in the forward pass of
deep learning model. Later, [11] constructs the whitening matrix based on
zero-phase component analysis to mitigate the stability issue in BN and PCA
[25]. [12] appends a PCA projection layer behind the network to spread
descriptors. [26] finds that the whitening matrix exerted on features can be
instead applied on model weights. They prove that their methods are
mathematically equivalent to ZCA [11] but the computation complexity is
largely reduced. WADAM [27] extends the idea in [26] with momentum, adaptive
dampening and proposes new optimizers for network training.

In VPR, [7] points out that whitening can balance out the co-occurrence of
features, which is generally beneficial to retrieval tasks. Whitened features
can boost the performance of the network on the current datasets and improve
the generalization ability of the network across datasets [4]. However,
Whitening adds a strong prior on the distribution of the descriptors and is
not involved in the training process.

Our work focuses on the global descriptor extraction procedure, since the
quality of the global descriptor is the bottleneck of later local descriptor
matching process. Unlike previous works in VPR [4, 17, 28, 5], we do not
propose new aggregation methods or new loss functions, but instead focus on
the feature distribution that lacks attention. We propose a new method to
optimize the feature distribution, which is parallel to most VPR methods
before-head.

##  III Method

###  III-A Gradient Distribution Problem

We claim that pair-wise loss function in VPR is one of the reasons for the
degraded distribution of gradients. A general pair-wise loss function can be
written as

$L=\sum_{i}f(s_{i,1},...,s_{i,N})$ (1)

where $s_{i,j}$ is the similarity between query $i$ and sample $j$. It is
usually defined as cosine similarity or L2 distance. These two are equivalent
when the descriptors are normalized. For simplicity, we will use the L2
distance.

$s_{i,j}=||\vec{p}_{i}-\vec{p}_{j}||^{2},\quad\vec{p}_{i}\in R^{C}$ (2)

$\vec{p}_{i}$ and $\vec{p}_{j}$ are the descriptors encoded by the network.
The gradient of the descriptor is

$\frac{\partial L}{\partial\vec{p}_{i}}=\sum_{j=1}^{N}\frac{\partial
f}{\partial s_{i,j}}2(\vec{p}_{i}-\vec{p}_{j})$ (3)

Since $2\frac{\partial f}{\partial s_{i,j}}$ is a scalar, Eq. 3 can be written
as

$\frac{\partial
L}{\partial\vec{p}_{i}}=\sum_{j=1}^{N}\alpha_{j}(\vec{p}_{i}-\vec{p}_{j})$ (4)

Here, we use contrastive loss as an example. The contrastive loss (CL) is
shown in Eq. 5.

$L_{CL}(\vec{p_{i}})=\sum_{j=1}^{N}\phi_{i,j}s_{i,j}+(1-\phi_{i,j})max(\tau-
s_{i,j},0)$ (5)

$\phi_{i,j}=1$ if sample i and j are positive pair and 0 otherwise. $\tau$ is
the threshold. The derivative is

$\frac{\partial
L}{\partial\vec{p}_{i}}=\sum_{j=1}^{N}[\phi_{i,j}-(1-\phi_{i,j})sign(\tau-
s_{i,j})](\vec{p}_{i}-\vec{p}_{j})$ (6)

In Eq. 6, since $\phi_{i,j}-(1-\phi_{i,j})sign(\tau-s_{i,j})$ is a scalar, the
gradient is a linear combination of the query and database descriptors. Since
the projection of descriptors in the principal space is much larger than the
projection in the complement space, a linear combination of these principal
components is generally larger than other components. Therefore, gradients
will also have large components in the principal space.

###  III-B Gradient Rectification Module

The whole process can be summarized in Fig. 3. During training, we feed images
into the network to get descriptors. These descriptors are used to update the
estimated covariance matrix $P$. Loss is computed according to the learning
task, for instance, GCL [4] or Triplet Loss. Then we backward the loss and
capture gradients at the descriptor level. The projection matrix $P^{*}$ is
applied here and the modified gradients flow into the network.

####  III-B1 Covariance Matrix Estimation

A brutal way to estimate the covariance matrix is to feed all training
examples into the network and calculate the covariance matrix with them after
each epoch. This can be applied when the size of the dataset is relatively
small. However, when the size of the dataset is large, this operation is time
consuming. What is worse, since the descriptors are usually designed to be
high dimensional, the number of examples inside one batch is inadequate to
provide an accurate estimation of the covariance matrix, if using the regular
formula below:

$P=\frac{1}{B-1}\sum_{i=1}^{B}(\boldsymbol{x}-\bar{\boldsymbol{x}})(\boldsymbol{x}-\bar{\boldsymbol{x}})^{T}.$
(7)

Here, $P$ is the estimated covariance matrix. $B$ is the batch size.
$\boldsymbol{x}$ is the descriptor and $\bar{\boldsymbol{x}}$ is the mean
descriptor within the current batch.

Therefore, non-singularity and positive-definite property of the covariance
matrix may not be preserved. To this end, we borrow the idea, memory queue,
from unsupervised learning domain [29]. We manage a queue of size K. The
current batch is enqueued to the memory queue and the oldest batch in the
queue is removed. The queue represents a sampled subset of all the instances.
Covariance matrix are estimated using the standard form from instances in
memory queue.  In order to prevent numerical instability, we also add a small
number ($10^{-3}$) to each of the diagonal elements of P.

####  III-B2 Projection Matrix Generation

We apply eigenvalue decomposition on $P$:

$P=Udiag(\lambda_{1},...,\lambda_{n})U^{T}$ (8)

$\lambda_{1},...,\lambda_{n}$ are eigenvalues in descending order. Then we
form the projection matrix $P^{*}$ to modify gradients.

$P^{*}=U\mathop{diag}(\frac{\bar{\lambda}}{\lambda_{1}},...,\frac{\bar{\lambda}}{\lambda_{n}})^{s}U^{T}$
(9)

$\bar{\lambda}$ is the mean of the eigenvalues. To eliminate the scale
problem, we put $\bar{\lambda}$ instead of a constant on the numerator. Since
at the initial stage of the training process, the elements in the descriptors
are all near zero, causing all the eigenvalues to be very small. A constant in
the numerator may cause a large scaling value along this dimension, leading to
unstable training at the beginning. $s$ is a hyper-parameter controlling the
rectification rate.  Large $s$($>1$) causes big rectifying value and leads to
unstable training while small $s$($<0.5$) is unable to push the network to
generate descriptors uniformly. We manually set this parameter to 1.

####  III-B3 Gradient Rectification

We rectify gradients using projection matrix $P^{*}$ as follows

$g^{*}=P^{*}g$ (10)

Modified gradients $g^{*}$ flows into the backbone network. Next, optimizers
like SGD or Adam[30] are adopted to update network’s parameters. In the ideal
case, when the covariance matrix of gradients is only different from the
covariance matrix of the descriptors by a constant multiplier, it is trivial
to show that the variance of gradients is identical along each dimension.

The space complexity of GRM is $(K+C)C$, where $C$ is the dimension of
descriptors. $KC$ is the size of the memory queue and $C^{2}$ is the size of
the covariance matrix. The time complexity for covariance matrix estimation is
$O(KC^{2})$ and for eigenvalue decomposition is $O(C^{3})$. GRM is not
involved in the inference stage, it doesn’t impact inference speed of the
model.

###  III-C Datasets

####  III-C1 MSLS

MSLS is a large scale place recognition dataset that contains images taken in
30 different cities [19]. The training set contains over 500k query images and
900k map images. The validation set consists of 19k map images and 11k query
images and the test set has 39k map images and 27k query images. Two images
are considered as similar if they are taken by cameras located within 25m of
distance and with less than 40∘ of viewpoint variation.

####  III-C2 Pittsburg

This dataset contains images of urban environments gathered from google street
view in the city of Pittsburgh, Pennsylvania, USA [5]. We use the test set of
Pitts30K to test our model. The test set of Pitts30K contains 6k query images
and 10k map images.

####  III-C3 Tokyo 24/7

It consists of images taken in Tokyo, Japan, with large variations of
illumination, as they are taken during day and night [31]. The query set
consists of 315 images and the map contains 76k photos.

####  III-C4 Nordland

This dataset is recorded on a train at four different seasons [32]. Each image
is taken at the same place as corresponding ones in other seasons. We use the
test subset, which consists of 3.5k images in each season, to test our model.
The tolerance for positive match is set to 10 frames and we use images in
summer as queries and images in winter as databases, as PatchNetVLAD[1].

###  III-D Implementation Details

We use GCL Loss [4] to train our model. We select ResNet50 [33] and ResNext101
[34] as backbone to extract features. GeM pooling is applied to get the global
descriptor. The memory queue size (K) is 10,240. We use Adam optimizer[30] and
set the learning rate to 1e-4. We train our model for 200 epochs (10k pairs in
MSLS per epoch). Following GCL [4], we use the weights pretrained on ImageNet
and only train the last two blocks of the backbone. For GRM, the rectification
rate is set to 1. We train our network on MSLS training set and test them on
other datasets. Top-N recall (R@N) is adopted as evaluation metric.

TABLE I: Comparison with State-of-the-art Methods.

Recall (%) | MSLS (Test. set) | MSLS (Val. set) | Pittsburgh 30k | Nordland | Tokyo247
---|---|---|---|---|---
|  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10
NetVLAD 64 [5] |  45.1 |  58.8 |  63.7 |  70.1 |  80.8 |  84.9 |  68.6 |  84.7 |  88.9 |  - |  - |  - |  34.0 |  47.6 |  57.1
NetVLAD 16 [5] |  39.4 |  53.0 |  57.5 |  70.5 |  81.1 |  84.3 |  70.3 |  84.1 |  89.1 |  - |  - |  - |  37.8 |  53.3 |  61.0
Patch-NetVLAD [1] |  48.1 |  57.6 |  60.5 |  79.5 |  86.2 |  87.7 |  87.7 |  94.5 |  95.9 |  71.3 |  80.9 |  82.2 |  86.0 |  88.6 |  90.5
GCL (ResNet50) [4] |  52.9 |  65.7 |  71.9 |  74.6 |  84.7 |  88.1 |  79.9 |  90.0 |  92.8 |  44.7 |  58.8 |  65.3 |  58.7 |  71.7 |  76.8
GCL (ResNext101) [4] |  62.3 |  76.2 |  81.1 |  80.9 |  90.7 |  92.6 |  79.2 |  90.4 |  93.2 |  69.9 |  85.6 |  90.1 |  58.1 |  74.3 |  78.1
TransVPR w/o. rerank [28] |  48.0 |  67.1 |  73.6 |  70.8 |  85.1 |  89.6 |  73.8 |  88.1 |  91.9 |  - |  - |  - |  - |  - |  -
Ours (ResNet50) |  62.3 |  75.5 |  78.7 |  82.4 |  90.2 |  91.7 |  82.8 |  91.9 |  93.9 |  66.2 |  89.9 |  93.0 |  55.9 |  70.8 |  76.2
Ours (ResNext101) |  64.9 |  77.2 |  81.1 |  82.6 |  90.7 |  92.6 |  81.3 |  91.4 |  93.5 |  82.3 |  92.4 |  95.2 |  62.5 |  75.2 |  82.5

###  III-E Comparison with State-of-the-art Methods.

For a fair comparison, we only compare the first stage of all methods without
the re-ranking process. We compare the performance of GRM with state of the
arts in Table I. On MSLS dataset, even though the baseline methods, GCL-
ResNet50 and GCL-ResNext101, are superior than all the other methods, our
methods can still boost the performance. ResNext101 trained by our method can
achieve the best R@5 in the official MSLS challenge, at 77.2%. For ResNet50,
the R@5 increases nearly 10%, and is comparable to the larger ResNext101
model. Our methods can also bring huge improvement on Nordland dataset. This
is due to the fact that season variation consists a large portion of the MSLS
dataset. Therefore, our methods trained on MSLS generalize well on Nordland
dataset, increasing over 5% on ResNext101 and over 30% on ResNet50. In urban
environment, Pittsburg 30K and Tokyo24/7, our method can bring moderate
improvement on the baseline method, between 1~2%, compared to the baseline
method. Our best R@5 on Pittsburg30K (91.9%) and Tokyo24/7 (75.2%) is lower
than PatchNetVLAD (94.5% on Pittsburg30K and 88.6% on Tokyo24/7). This is
because Patch-NetVLAD is fine-tuned on Pittsburg30K training dataset, with
images including drastic variations in illumination. However, in order to keep
aligned with other methods, we only train our method on MSLS dataset.

###  III-F Comparison with Other Whitening Methods

We compare the performance of GRM and other feature whitening methods. We
select GCL [4] as our baseline and use multiple strategies to expand the
feature space. For Spread Loss [8]. We set the weight for the regularization
term to 1. For AP Loss [17]. We use the GCL Label in [4] to sample 2 positive
images and 7 negative images. Negative images are half hard negative
(similarity score greater than 0 but less than 0.5) and half negative
(similarity score equals 0). We compare our method to [27]. Standard Adam [30]
optimizer is replaced by WADAM [27] in the GCL [4] framework.

The training procedure follows the previous section and we show the R@N index
on MSLS val and test sets in Table II. GRM exceeds all methods in Table II.
The performance drops when adding the extra regularization in [8]. We
conjecture that the inclusion of regularization might break the original
distribution and leads to a sub-optimal distribution. As we can see, AP Loss
[17], as the representation of multi-pair form, increase R@5 by 4% on MSLS
test set but is still lower than our method. Meanwhile, AP Loss suffers from
large memory consumption and long training time. As for the WADAM [27], it can
bring slight improvement on MSLS val., but when evaluated on test set, it
under-performs the baseline method.

The result corroborates that it is neither the weak repulsive force nor the
limited batch size that causes the unbalanced distribution of descriptors but
the phenomenon that gradients and descriptors fall into the same low
dimensional principal space. Since gradients won’t push the descriptors out of
this principal space, descriptors will dwell in this principal space and never
be uniformly distributed.

TABLE II: Comparison with Whitening Method.

Recall (%) | MSLS (Test. set) | MSLS (Val. set)
---|---|---
|  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10
ResNet50 [4] |  52.9 |  65.7 |  71.9 |  74.6 |  84.7 |  88.1
Spread Loss [8] |  30.8 |  47.2 |  53.6 |  49.1 |  63.4 |  68.6
AP Loss [17] |  53.4 |  69.7 |  75.0 |  72.3 |  82.4 |  84.9
WADAM [27] |  48.4 |  62.6 |  68.7 |  72.9 |  85.1 |  87.9
ResNet50 + GRM |  62.3 |  75.5 |  78.7 |  82.4 |  90.2 |  91.7

Figure 4: Gradient and descriptor distribution. (A) and (C) are alignment of
eigenvectors between different training stage w/o. GRM and with GRM. (B) and
(D) are the eigenvalue distribution of descriptors and gradients with and w/o.
GRM.

###  III-G Ablation Study

####  III-G1 Gradient distributions analysis

In order to illustrate the problem that the descriptors encoded by the network
fall into a low dimensional space and do not change during the training
process, we plot the correlation between the principal space of descriptors in
different stages during training process in Fig. 4. In the upper figure of
(A), pixel value at (i, j) is the dot product between the i-th eigenvetor of
the feature space encoded by model after training 30 epochs and the j-th
eigenvector of the feature space encoded by model after training 90 epochs,
without GRM. Here, all the eigenvectors are sorted by their corresponding
eigenvalue. The i-th eigenvector is the eigenvector corresponds to the i-th
largest eigenvalue. Therefore, the brighter the pixel is, the more aligned two
eigenvectors are. In the lower figure of (A), pixel value at (i, j) is the dot
product between the i-th eigenvetor of the feature space and the j-th
eigenvector of the gradient space. Fig. 4 (C) is the same as (A) but the model
is trained with GRM.

In Fig. 4 (A), We notice that brightest pixels are on the diagonal line of the
image, indicating that eigenvectors are aligned respectively during training
process and the principle space does not change. The same situation remains
between the principle space of descriptors and gradients. Given that
descriptors and gradients all concentrate in a low principal space (the green
curve in Fig. 4 B and D), we claim that descriptors trapped in a low
dimensional space is the result of the gradients which fall into the same
dimensional space.

When GRM is applied, bright pixels disperse in Fig. 4 (C). Hence the principal
space between descriptors and gradients are not aligned. Gradients are pulled
out of the principal space of descriptors and influence the distribution of
descriptors in later steps. We also plot the eigenvalue of the descriptor
space and the gradient space in Fig. 4 (B) and (D). We can find that the
descriptors and gradients distribute more uniformly when GRM (red curve) is
applied.

####  III-G2 Experiments in Classification Tasks.

To evaluate the performance of the proposed method on classification tasks, we
test our method on CIFAR-10, CIFAR-100, Caltech 101 and Caltech 256 datasets
[35]. We adopt the prototype learning method in [36]. We use ResNet18 as
backbone to extract features and use average pooling to get the descriptor.
The feature dimension is set to 512. We use the GCPL loss in [36] to train our
network. The gradient rectification module is applied to descriptors and
prototypes simultaneously. The model is optimized using SGD with momentum of
0.9. The initial learning rate is set to 0.05 and is multiplied by 0.7 every
20 epochs. We initialize all prototypes as zero [36] to avoid noise. The whole
training process lasts for 200 epochs. We use the top-1 accuracy to evaluate
all the models. Results are summarized in Table III.

GRM can bring improvement in prototype learning framework. It brings around
0.5% improvement on CIFAR-10 and CIFAR-100 using the GCPL framework. On
Caltech dataset, even though the GCPL framework [36] underperforms the
baseline method, our method can still achieve better results, about 1~2% than
the baseline method.

TABLE III: Performance on classification task.

Accuracy (%) |  CIFAR-10 |  CIFAR-100 |  Caltech101 |  Caltech256
---|---|---|---|---
ResNet18 (Baseline) |  93.2 |  73.3 |  74.3 |  73.5
ResNet18 + GCPL |  94.4 |  75.5 |  71.5 |  71.4
ResNet18 + GCPL + GRM |  94.9 |  75.8 |  75.2 |  75.3

####  III-G3 Comparison with Other Estimation Method

Another way to estimate the covariance matrix is running average. The
procedure is given below:

$\displaystyle N_{k+1}$ $\displaystyle=N_{k}+b$ (11)
$\displaystyle\bar{x}_{k+1}$
$\displaystyle=\frac{N_{k+1}-b}{N_{k+1}}\bar{x}_{k}+\frac{b}{N_{k+1}}\sum_{i=1}^{b}x_{i}$
$\displaystyle P_{k+1}$
$\displaystyle=\frac{N_{k+1}-b}{N_{k+1}}P_{k}+\frac{b}{N_{k+1}}\sum_{i}^{b}(x_{i}-\bar{x}_{k+1})(x_{i}-\bar{x}_{k+1})^{T}$

$b$ is the batch size. $k$ is the number of step. $x$ is the descriptor. $P$
is the co-variance matrix. $\bar{x}$ is the running mean. $P$ is initialized
as identity and $\bar{x}$ is initialized as zero. In each step, the projection
matrix is calculated from $P_{k}$ and multiplies the backward gradients. We
use the same setting to test the running average method and compare it with
baseline and memory queue method.

The running average method can still outperform the baseline method on MSLS.
It gets a top-5 recall rate of 73.2% on MSLS test set, higher than the
baseline (65.5%), which is still significant. On MSLS validation set, it gets
top-5 recall at 88.4%, 4% higher than the baseline method. It is even slightly
better on Pittsburg30K than the memory queue method. Moreover, running average
method doesn’t need to maintain a large memory queue to store descriptors,
cutting off the memory consumption.

TABLE IV: Comparison with Running Average Method.

method | MSLS (Test. set) | MSLS (Val. set)
---|---|---
|  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10
GCL (ResNet50) [4] |  52.9 |  65.7 |  71.9 |  74.6 |  84.7 |  88.1
ResNet50 + Average (Sqrt) |  61.5 |  73.2 |  77.2 |  84.0 |  88.7 |  90.7
ResNet50 + Bank (Linear) |  62.3 |  75.5 |  78.7 |  82.4 |  90.2 |  91.7

####  III-G4 Memory Queue Size

Here we conduct experiments to analyze the effect of memory queue size K. The
result is summarized in Table V. We find that small queue will harm the
performance of the model. This is because the estimated covariance matrix is
not accurate from few samples. Queue size larger than 10,000 suffices an
accurate estimation.

TABLE V: Memory Queue Size Comparison

Memory Queue Size | MSLS (Test. set) | MSLS (Val. set)
---|---|---
|  R@1 |  R@5 |  R@10 |  R@1 |  R@5 |  R@10
GCL (ResNext101) [4] |  62.3 |  76.2 |  81.1 |  80.9 |  90.7 |  92.6
K=5120 |  59.1 |  71.4 |  75.7 |  80.4 |  87.6 |  89.9
K=7680 |  62.4 |  73.8 |  76.8 |  81.9 |  89.1 |  90.6
K=10240 |  64.9 |  77.2 |  81.1 |  82.6 |  90.7 |  92.6

##  IV Conclusions

In this article, we point out the gradient problem in training deep neural
network for global feature extraction in VPR. The degraded distribution of
gradients leads to the unbalanced distribution of descriptors. We propose our
new GRM and memory queue method to rectify gradients. This method can improve
the existing state of the art methods on global descriptor extraction. We also
demonstrate that our method could be generalized to classification task. For
future work, the GRM Module can be applied to other metric learning tasks like
unsupervised learning or re-ID in tracking task.

## References

  * [1] S. Hausler, S. Garg, M. Xu, M. Milford, and T. Fischer, “Patch-netvlad: Multi-scale fusion of locally-global descriptors for place recognition,” in _IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2021, virtual, June 19-25, 2021_. Computer Vision Foundation / IEEE, 2021, pp. 14 141–14 152.
  * [2] B. Cao, A. Araujo, and J. Sim, “Unifying deep local and global features for image search,” in _Computer Vision - ECCV 2020 - 16th European Conference, Glasgow, UK, August 23-28, 2020, Proceedings, Part XX_ , ser. Lecture Notes in Computer Science, A. Vedaldi, H. Bischof, T. Brox, and J. Frahm, Eds., vol. 12365. Springer, 2020, pp. 726–743. [Online]. Available: https://doi.org/10.1007/978-3-030-58565-5_43
  * [3] H. Noh, A. Araujo, J. Sim, T. Weyand, and B. Han, “Large-scale image retrieval with attentive deep local features,” in _IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017_. IEEE Computer Society, 2017, pp. 3476–3485. [Online]. Available: https://doi.org/10.1109/ICCV.2017.374
  * [4] M. Leyva-Vallina, N. Strisciuglio, and N. Petkov, “Generalized contrastive optimization of siamese networks for place recognition,” _CoRR_ , vol. abs/2103.06638, 2021. [Online]. Available: https://arxiv.org/abs/2103.06638
  * [5] R. Arandjelovic, P. Gronát, A. Torii, T. Pajdla, and J. Sivic, “Netvlad: CNN architecture for weakly supervised place recognition,” _CoRR_ , vol. abs/1511.07247, 2015. [Online]. Available: http://arxiv.org/abs/1511.07247
  * [6] A. Miech, I. Laptev, and J. Sivic, “Learnable pooling with context gating for video classification,” _CoRR_ , vol. abs/1706.06905, 2017. [Online]. Available: http://arxiv.org/abs/1706.06905
  * [7] H. Jégou and O. Chum, “Negative evidences and co-occurences in image retrieval: The benefit of pca and whitening,” in _ECCV_ , 2012.
  * [8] X. Zhang, F. X. Yu, S. Kumar, and S. Chang, “Learning spread-out local feature descriptors,” in _IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017_. IEEE Computer Society, 2017, pp. 4605–4613. [Online]. Available: https://doi.org/10.1109/ICCV.2017.492
  * [9] F. Wang and H. Liu, “Understanding the behaviour of contrastive loss,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2021, pp. 2495–2504.
  * [10] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in _Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015_ , ser. JMLR Workshop and Conference Proceedings, F. R. Bach and D. M. Blei, Eds., vol. 37. JMLR.org, 2015, pp. 448–456. [Online]. Available: http://proceedings.mlr.press/v37/ioffe15.html
  * [11] A. Kessy, A. Lewin, and K. Strimmer, “Optimal whitening and decorrelation,” _The American Statistician_ , vol. 72, no. 4, pp. 309–314, 2018.
  * [12] A. Gordo, J. Almazán, J. Revaud, and D. Larlus, “Deep image retrieval: Learning global representations for image search,” in _Computer Vision \- ECCV 2016 - 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part VI_ , ser. Lecture Notes in Computer Science, B. Leibe, J. Matas, N. Sebe, and M. Welling, Eds., vol. 9910. Springer, 2016, pp. 241–257. [Online]. Available: https://doi.org/10.1007/978-3-319-46466-4_15
  * [13] G. Csurka, C. Dance, L. Fan, J. Willamowski, and C. Bray, “Visual categorization with bags of keypoints,” in _Workshop on statistical learning in computer vision, ECCV_ , vol. 1, no. 1-22. Prague, 2004, pp. 1–2.
  * [14] H. Jégou, M. Douze, C. Schmid, and P. Pérez, “Aggregating local descriptors into a compact image representation,” in _The Twenty-Third IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2010, San Francisco, CA, USA, 13-18 June 2010_. IEEE Computer Society, 2010, pp. 3304–3311. [Online]. Available: https://doi.org/10.1109/CVPR.2010.5540039
  * [15] E. Ong, S. Husain, M. Bober-Irizar, and M. Bober, “Deep architectures and ensembles for semantic video classification,” _CoRR_ , vol. abs/1807.01026, 2018. [Online]. Available: http://arxiv.org/abs/1807.01026
  * [16] F. Radenović, G. Tolias, and O. Chum, “Fine-tuning cnn image retrieval with no human annotation,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 41, no. 7, pp. 1655–1668, 2018.
  * [17] J. Revaud, J. Almazán, R. S. Rezende, and C. R. de Souza, “Learning with average precision: Training image retrieval with a listwise loss,” in _2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019_. IEEE, 2019, pp. 5106–5115. [Online]. Available: https://doi.org/10.1109/ICCV.2019.00521
  * [18] F. Çakir, K. He, X. Xia, B. Kulis, and S. Sclaroff, “Deep metric learning to rank,” in _IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2019, Long Beach, CA, USA, June 16-20, 2019_. Computer Vision Foundation / IEEE, 2019, pp. 1861–1870.
  * [19] F. Warburg, S. Hauberg, M. López-Antequera, P. Gargallo, Y. Kuang, and J. Civera, “Mapillary street-level sequences: A dataset for lifelong place recognition,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020_. Computer Vision Foundation / IEEE, 2020, pp. 2623–2632.
  * [20] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _Int. J. Comput. Vis._ , vol. 60, no. 2, pp. 91–110, 2004. [Online]. Available: https://doi.org/10.1023/B:VISI.0000029664.99615.94
  * [21] H. Bay, T. Tuytelaars, and L. Van Gool, “Surf: Speeded up robust features,” in _Computer Vision – ECCV 2006_ , A. Leonardis, H. Bischof, and A. Pinz, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2006, pp. 404–417.
  * [22] J. Philbin, O. Chum, M. Isard, J. Sivic, and A. Zisserman, “Object retrieval with large vocabularies and fast spatial matching,” in _2007 IEEE Conference on Computer Vision and Pattern Recognition_ , 2007, pp. 1–8.
  * [23] O. Chum and J. Matas, “Matching with prosac - progressive sample consensus,” in _2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)_ , vol. 1, 2005, pp. 220–226 vol. 1.
  * [24] G. Desjardins, K. Simonyan, R. Pascanu, _et al._ , “Natural neural networks,” _Advances in neural information processing systems_ , vol. 28, 2015.
  * [25] L. Huang, D. Yang, B. Lang, and J. Deng, “Decorrelated batch normalization,” in _2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018_. Computer Vision Foundation / IEEE Computer Society, 2018, pp. 791–800.
  * [26] S. Markovich-Golan, B. Battach, and A. Bleiweiss, “Feature whitening via gradient transformation for improved convergence,” _CoRR_ , vol. abs/2010.01546, 2020. [Online]. Available: https://arxiv.org/abs/2010.01546
  * [27] H. Yong and L. Zhang, “An embedded feature whitening approach to deep neural network optimization,” in _the European Conference on Conputer Vision_ , 2022\.
  * [28] R. Wang, Y. Shen, W. Zuo, S. Zhou, and N. Zheng, “Transvpr: Transformer-based place recognition with multi-level attention aggregation,” _CoRR_ , vol. abs/2201.02001, 2022. [Online]. Available: https://arxiv.org/abs/2201.02001
  * [29] K. He, H. Fan, Y. Wu, S. Xie, and R. B. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020_. Computer Vision Foundation / IEEE, 2020, pp. 9726–9735. [Online]. Available: https://doi.org/10.1109/CVPR42600.2020.00975
  * [30] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2015. [Online]. Available: http://arxiv.org/abs/1412.6980
  * [31] A. Torii, R. Arandjelovic, J. Sivic, M. Okutomi, and T. Pajdla, “24/7 place recognition by view synthesis,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 40, no. 2, pp. 257–271, 2018. [Online]. Available: https://doi.org/10.1109/TPAMI.2017.2667665
  * [32] D. Olid, J. M. Fácil, and J. Civera, “Single-view place recognition under seasonal changes,” _CoRR_ , vol. abs/1808.06516, 2018. [Online]. Available: http://arxiv.org/abs/1808.06516
  * [33] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016_. IEEE Computer Society, 2016, pp. 770–778. [Online]. Available: https://doi.org/10.1109/CVPR.2016.90
  * [34] S. Xie, R. B. Girshick, P. Dollár, Z. Tu, and K. He, “Aggregated residual transformations for deep neural networks,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017_. IEEE Computer Society, 2017, pp. 5987–5995. [Online]. Available: https://doi.org/10.1109/CVPR.2017.634
  * [35] L. Fei-Fei, R. Fergus, and P. Perona, “One-shot learning of object categories,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 28, no. 4, pp. 594–611, 2006. [Online]. Available: https://doi.org/10.1109/TPAMI.2006.79
  * [36] H.-M. Yang, X.-Y. Zhang, F. Yin, and C.-L. Liu, “Robust classification with convolutional prototype learning,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 3474–3482.


# Model Predictive Control with adaptive resilience for Denial-of-Service
Attacks mitigation on a Regulated Dam

Raffaele G. Cestari1, Stefano Longari1, Stefano Zanero1, Simone Formentin1
1Department of Electronics, Information, and Bioengineering, Politecnico di
Milano, Milano, Italy. Email to<EMAIL_ADDRESS>paper
is partially supported by FAIR (Future Artificial Intelligence Research)
project, funded by the NextGenerationEU program within the PNRR-PE-AI scheme
(M4C2, Investment 1.3, Line on Artificial Intelligence), by the Italian
Ministry of Enterprises and Made in Italy in the framework of the project 4DDS
(4D Drone Swarms) under grant no. F/310097/01-04/X56 and by the PRIN PNRR
project P2022NB77E “A data-driven cooperative framework for the management of
distributed energy and water resources” (CUP: D53D23016100001), funded by the
NextGeneration EU program.

###### Abstract

In recent years, SCADA (Supervisory Control and Data Acquisition) systems have
increasingly become the target of cyber attacks. SCADAs are no longer
isolated, as web-based applications expose strategic infrastructures to the
outside world connection. In a cyber-warfare context, we propose a Model
Predictive Control (MPC) architecture with adaptive resilience, capable of
guaranteeing control performance in normal operating conditions and driving
towards resilience against DoS (controller-actuator) attacks when needed.
Since the attackers’ goal is typically to maximize the system damage, we
assume they solve an adversarial optimal control problem. An adaptive
resilience factor is then designed as a function of the intensity function of
a Hawkes process, a point process model estimating the occurrence of random
events in time, trained on a moving window to estimate the return time of the
next attack. We demonstrate the resulting MPC strategy’s effectiveness in 2
attack scenarios on a real system with actual data, the regulated Olginate dam
of Lake Como.

##  I INTRODUCTION

In recent years, the number of attacks on SCADA (Supervisory Control and Data
Acquisition) systems has soared, leading to 2021 as the year in which, for the
first time, the manufacturing and critical infrastructure industries suffered
the most tremendous economic and structural damage, according to the IBM 2021
report. The most famous is the Stuxnet malware (2010), which targeted a
Nuclear Power Plant in Iran. Night Dragon (oil and gas, 2010), Shamoon (oil
and gas, 2012), New York dam (2013), Ukraine power grid (2015-2016), up to the
most recent ones (OMRON, Schneider PLCs, 2022) are just some of the malwares
and targets of the last few years. Citing Drias et al. [1], avoiding attack is
quite impossible. Industrial networks should perform profiling, network
traffic monitoring, and attack detections. Industrial network infrastructure
must be able to ensure control operations under attacks due to the criticality
of the systems that they operate. As they point out, it is impossible to
prevent the attacks. This depends on the disparity between the resources
needed by attackers and defenders: developing malware is relatively simple and
easy to deploy, and conversely, repairing the damage caused by it is expensive
in terms of money and time. For this reason, it is crucial to study the
development of control schemes from a performance maximization perspective and
consider the possibility that they are victims of cyber threats. SCADA systems
often constitute national and international safety critical infrustructures
such as power grids and dams. For this reason, they represents privileged
targets of cyber-warfare and might cause enormous social and economic damages
on national and international scale (see [2], [3], [4]).

The primary attacks in the literature are Deception, False Data injection,
Eavesdropping, and Denial of Service (DoS). See [5]. The first two consist of
the injection of corrupted data into the system. Deception, more insidious
than false data injection, introduces corrupted data with continuity
properties, making detecting the attack more complex. See [6]. Eavesdropping
consists of violating the confidentiality of the communication, sniffing
packets on the Local Area Network (LAN), or intercepting wireless
transmissions.

DoS attacks are the most common (an easy to deploy) and lead to the
interruption of communications between the controller and actuator (or
analogously between controller and sensors) and, in the worst case, to the
application of malicious control actions. See [7], [8]. Given its ease of
development and applicability, in this work we focus on developing a model
predictive control architecture capable of mitigating the impact of DoS
attacks on the communication between controller and actuator.

In the literature, several formulations of Resilient MPC aim to respond
through an appropriate control scheme design to one or more of these threats.
The authors of [9] propose an observer that, combined with MPC, is capable of
detecting deception attacks. In [10] they propose an artificial neural network
(ANN) based MPC to identify occurrence of false data injection attacks in DC
microgrids. The work in [11] proposes a resilient MPC architecture against DoS
attacks jamming the communication between controller and actuator. Their main
idea is to store actuator-side the whole vector of optimal control actions
computed by the MPC over the prediction horizon in such a way that, even if
the data transmission between controller and actuator is compromised, the
actuator is still capable of executing optimal control actions.

Our problem statement takes a stronger assumption: we assume that the actuator
side is hijacked by the attacker. Communication between controller and
actuator is compromised and the actuator is no longer in control, hence the
proposed solution of [11] would not be feasible.

Suppose that the attacker can access the actuator at any time and tamper with
the optimal control action for that instant. Our architecture will mitigate
the damage caused by fiddling with the optimal control. In particular, the
proposed solution employs two MPC controllers which simultaneously solve two
different optimization problems. The master controller, hereafter called
ResMPC, is responsible for calculating the optimal control action applied to
the system. This controller determines the optimum by balancing the
performances obtainable on the prediction horizon by including a
”regularizing” term that weighs the distance of its optimal control action
from that provided by the second controller, from now on called SafeMPC, which
instead assumes that from the next instant and after that, the actuator will
be out of control and applying the optimal control action on the system will
no longer be possible.

From this perspective, it is clear how impactful it is to correctly balance
performance tracking with respect to the safe control action that assumes the
presence of an attack in the immediate future. Choosing the conservative path
in the absence of attacks will lead to a loss of performance, while in the
presence of frequent attacks, pursuing only the quality of the control action
on the prediction horizon will cause a blind approach to the onset of attacks
with a consequent damage overhead compared to the precautionary solution.

For this reason, we propose an online adaptive solution that can estimate the
onset of attacks in the near future and increase or decrease the resilience
factor responsible for the balance between performance and safety.

This is possible thanks to Hawkes processes ([12]), a class of point processes
that allows us to estimate the return time of an event whose occurrence is
apparently random in time. Hawkes processes find wide use in the literature in
all those applications where predicting events without uniform temporal
sampling is necessary. Some examples are finance ([13], [14]), neuroscience
([15]), seismology ([16]), and social media ([17]). We have references also in
the field of cybersecurity. The authors of [18] used multivariate Hawkes
processes to predict attack frequency. In our architecture, at each simulation
step, we train the Hawkes process on a fixed-length moving window (expanding
at the beginning of the simulation as we assume that we have no data before
that collected online) and calculate its intensity function as the return time
of the attack varies. Since the intensity function measures the probability of
event occurrence, it will increase as the expected return time decreases. By
defining a fixed look-ahead horizon, we integrate the intensity function to
correspondingly define the value of the resilience factor. In this way, we
algebraically associate the intensity function with the resilience factor,
representing the estimated attack condition in the near future.

We apply the proposed methodology on a SCADA system, linear with saturation,
represented by the water release from the Olginate dam responsible for
regulating the Lake Como basin, Lombardy, Italy. The system to be regulated
(with application-related problems) and the ”performance-oriented” reference
MPC controller are described in [19] and [20]. The data used in the simulation
are authentic and come from the dam’s historical regulation archive. This case
study is well representative of a critical infrastructure for the management
of water resources in the Lombardy region, in which the lack of ”resilience”
to cyber attacks could lead to direct (flooding of roads, houses, commercial
and industrial activities in the city of Como) and indirect (failure to
release water to agricultural districts in periods of drought leads to
significant economic losses) economic damage. To focus exclusively on the
attack impact, we assume the perfect knowledge of the system disturbances
(e.g., the water inflow in the lake).

Given the application, we assume that we are working in a cyber warfare
context in which the attacker is interested in undermining the territorial and
infrastructural stability of the country. From this perspective, we assume
that it is capable of exercising through a DoS attack the interruption of
communications between the dam controller and the actuation of the release
valve and that it takes its place, applying the policy of ”maximum damage”
that is, given the conditions (states) of the system, evaluate which control
action could cause the greatest possible damage. To make this possible, we
assume that the attacker has access, at every time instant, to the measured
states of the system and that to determine the ”optimal attack”, it solves an
optimization problem equal and opposite to the one solved, at the same
instant, by the MPC controller, in order to maximize the damage caused by the
malicious control action.

We evaluate performance on actual data of year 2005, comparing 3 control
architectures:

  * •

MPC: the standard ”performance-oriented” MPC, described in [19] and [20].

  * •

SafeMPC: the ”safety-oriented” MPC, applying to the actual system exclusively
the control action computed by the SafeMPC controller.

  * •

Adaptive Resilient MPC (ARMPC): the proposed architecture, adapting online
resilience by balancing performance and safety given the expected return time
of the attacks.

We compare the controllers in 2 attack scenarios:

  * •

A time-varying periodic attack pattern, in which the attacks happen
periodically with changing period.

  * •

A pseudo-random binary attack pattern, in which the attacks happen pseudo-
randomly given the firing of a pseudo-random binary periodic triggering
signal.

We show that adaptive resilience is mandatory for adaptation to changing
attack conditions and that the proposed architecture, Adaptive Resilient MPC,
outperforms the competitors in both the attack scenarios.

The main contributions of the paper are the following:

  * •

As far as we are aware, for the first time, we propose a model predictive
control architecture capable of mitigating DoS attacks where the attacker
overrides the optimal control action and replaces it with the optimal
adversarial control action.

  * •

We propose an online adapting solution capable of predicting the attack
occurrence and correspondingly increasing the safety level, mitigating
damages, and, simultaneously, keeping a satisfactory performance tracking in
periods of absence of attacks.

  * •

For the first time, as far as we are aware, we exploit Hawkes processes in the
context of model predictive control to forecast cyber attack occurrence.

The remainder of the paper is as follows. In Section II we describe the
Adaptive Resilient MPC architecture and the optimal attack problem. In
particular, Section II-A describes the optimal control problem solved by the
main ResMPC controller. Section II-B describes the optimal control problem
solved by the auxiliary ”safety” controller. Section II-C describes the
adversary optimal control problem solved at each time step by the attacker.
Section II-D shows why an adapting resilience factor is needed. Section II-E
describes the Hawkes process and the resilience factor computation. In Section
III we show the numerical simulations and results in the 2 attack scenarios.
Section IV summarizes the main contributions of the paper and the future work.

##  II Mathematical Formulation

The Resilient MPC architecture consists in 2 MPC controllers. The main one,
ResMPC, is responsible for computing the optimal control action $u_{res}$, at
each time step, considering simultaneously control performance and safety. The
auxiliary controller, SafeMPC, computes, at each time step, the safe control
action $u_{safe}$ which is fed to the ResMPC. The attacker interrupts the
communication between controller and actuator and replaces the optimal control
action $u_{res}$ with its optimal malicious action $u_{DoS}$. Attack
occurrence is saved as a binary variable in a database. The attack detection
module provides the resilience factor $\lambda$ to the ResMPC after training
on historical archive and current time attack occurrence. Figure 1 summarizes
the Adaptive Resilient MPC architecture.

Figure 1: Adaptive Resilient MPC Architecture

###  II-A Resilient MPC Control Problem

ResMPC controller solves, at each time step, the following constrained, convex
optimization problem:

$\operatorname{min}_{u}\\{c_{f}||\textbf{s}-s_{max}||_{2}+c_{d}||\textbf{s}-s_{min}||_{2}+c||\epsilon||_{2}+\lambda||\textbf{u}-u_{safe}||\\}$
(1a) $s_{h+1}=s_{h}+3600\cdot 24(\hat{q_{h}}-u_{h})$ (1b) $u_{min}\leq
u_{h}\leq u_{max}$ (1c) $u_{h}\geq d_{h}+\epsilon_{h}\\\ $ (1d)
$u_{h}-u_{h-1}\leq 240\,\,m^{3}/s\,\,\,\,if\,\,l_{h}\leq 1.1\,\,m\\\ $ (1e)
$u_{h}-u_{h-1}\leq 360\,\,m^{3}/s\,\,\,\,if\,\,l_{h}>1.1\,\,m\\\ $ (1f)
$\textbf{s}\in R^{H},\textbf{u}\in R^{H},\epsilon\in R^{H}$ (1g)
$c_{f}=10^{-4},c_{d}=10^{-3},c=100$ (1h)

Equation (1a) shows the minimization problem. u $[m^{3}/s]$ is the control
action (e.g., the water release going out from the lake) vector over the
prediction horizon $H$. $c_{f}$, $c_{d}$ and $c$ are numerical coefficients
tuned to have the same order of magnitude for the 3 objectives whose value is
collected in Equations (1h), for further details, see [19], [20]. s $[m^{3}]$
is the lake volume vector. $s_{min}$ and $s_{max}$ are the minimum and maximum
lake volumes. The control goal is, simultaneously, maintaining the lake volume
between the minimum and maximum thresholds (dry and flood levels, first and
second cost components respectively) and satisfying the agricultural water
demand (third component). The fourth cost component is the key of resilience:
we want to minimize the distance of the following optimal control action that
will be applied to the system from the optimally computed safe control action
$u_{safe}$ $[m^{3}/s]$. $\lambda$ is the resilience factor, increasing it
leads the controller to safer decisions at the cost of losing performance
quality. Equation (1b) shows the daily mass balance equation for the lake
volume. $\hat{q}_{h}$ is the predicted water inflow entering the lake (in this
work we assume it coincides with the actual one) in the $h$ time-step of the
prediction horizon. Equation (1c) shows the control action constraints,
minimum $u_{min}$ and maximum $u_{max}$ respectively. These constraints are
non-linear functions of the lake level. For additional details see [19].
Equation (1d) shows an additional soft constraint on control action. Water
release must be close to the agricultural water demand $d$. The constraint
must be soft to include cases where the release is lower than the request (dry
periods) and where it must be above the request (flood periods). This is
ensured through the slack variable $\epsilon$, an additional penalty term in
the cost function (1a). Equations (1e) and (1f) ensure that the variation of
control action is always below the daily limits imposed by the Olginate dam
operation. These boundaries depend on $l\,\,[m]$, lake level.

###  II-B Safe MPC Control Problem

The safe controller solves, at each time step, a constrained convex
optimization problem whose role is computing the next safe control action
$u_{safe}$:

$\operatorname{min}_{u_{safe}}\\{c_{f}||\textbf{s}-s_{max}||_{2}+c_{d}||\textbf{s}-s_{min}||_{2}+c||\epsilon||_{2}\\}$
(2a) $s_{h+1}=s_{h}+3600\cdot 24(\hat{q_{h}}-u_{h})$ (2b) $u_{min}\leq
u_{safe}\leq u_{max}$ (2c) $u_{safe}\geq d_{h}+\epsilon_{h}\\\ $ (2d)
$u_{safe}-u_{t-1}\leq 240\,\,m^{3}/s\,\,\,\,if\,\,l_{t}\leq 1.1\,\,m\\\ $ (2e)
$u_{safe}-u_{t-1}\leq 360\,\,m^{3}/s\,\,\,\,if\,\,l_{t}>1.1\,\,m\\\ $ (2f)
$\textbf{s}\in R^{H},\textbf{u}\in R^{H},\epsilon\in R^{H}$ (2g)
$u_{1}=u_{safe}$ (2h) $u_{2}=...=u_{h}=...=u_{DoS}\,\,\,\forall\,\,2\leq h\leq
H$ (2i)

The structure is similar to the ResMPC, but there are substantial differences:

  * •

The optimally computed control action is exclusively the one calculated for
the next time instant, as Equation (2h) specifies. The controller is designed
assuming that from the following time instant, a DoS attack might occur. If it
happens, the controller will lose the link with the corresponding actuator
responsible for the execution of the computed control action. The system will
go open loop until the DoS attack finishes.

  * •

The controller assumes that the corrupted control action is equal to $u_{DoS}$
and that remains constant for the whole prediction horizon after the current
instant.

The other equations follow the same structure of the ResMPC. The optimal
control action $u_{safe}$ is computed minimizing the cost obtained after the
evolution of the system trajectory as a consequence of the corrupted actions.
Remark: we assume reasonably that we can anticipate which is the value of the
malicious control action $u_{DoS}$. This assumption derives from the fact
that, we assume the attacker will hit with the purpose of causing the highest
possible damages, thus applying exactly $u_{DoS}$ whose computation is
described in the next section.

Both the ResMPC and the SafeMPC optimization problems are solved using the
convex optimization solver CVX. See [21].

###  II-C Optimal Attacker

The attacker solves an ”adversary” optimization problem at each time instant.
Given the same information on the system (states and inputs), it calculates
the optimal control action, which, if applied at that instant on the system,
maximizes the objective function (which in the MPC optimization problem was
instead minimized). We note that the solved problem is no longer convex but
rather concave (maximizing a convex cost function is equivalent to minimizing
a concave one). This constrains the selection of solvers for the optimization
problem. We choose to replace CVX (used in MPC controllers) with Gurobi
([22]), a commercial solver capable of dealing with non-convex functions.

$\operatorname{max}_{u_{DoS}}\\{c_{f}||\textbf{s}-s_{max}||_{2}+c_{d}||\textbf{s}-s_{min}||_{2}+c||\epsilon||_{2}\\}$
(3a) $s_{h+1}=s_{h}+3600\cdot 24(\hat{q_{h}}-u_{h})$ (3b) $0\leq u_{DoS}\leq
u_{max}$ (3c) $u_{DoS}\geq d_{h}+\epsilon_{h}\\\ $ (3d) $\textbf{s}\in
R^{H},u_{DoS}\in R^{H},\epsilon\in R^{H}$ (3e)

The Equations (3) describe the constrained concave problem solved at each time
step by the attacker. We highlight the maximization problem in Equation (3a).
We also highlight the difference (with respect to Equations (1c) and (2c)) in
the lower bound in Equation (3c) on attacker control action which is $0$
$m^{3}/s$ and not $u_{min}$ (as it is the value given by the required minimum
environmental flow (see [19]), not respected by the attacker.). We assume the
attack duration is fixed and equal to $1$ day.

Figure 2: Maximum Damage control (upper panel), Lake Como level (lower panel)

Figure 2 shows in the lower panel the Lake Como level and in the upper panel
the correspondingly computed optimal ”maximum damage” control action, obtained
as solution of Equations (3). In the dry periods, the best choice for the
attacker is applying the maximum water release (compliant with the upper
bound, given by the dam physical structure), wasting the resource. In flood
periods, the best choice is the opposite, closing the dam and amplyifing the
flood effect. The optimal attack adapts online, monitoring the lake condition
at each time step.

###  II-D Adaptive Resilience Motivations

In this section, we illustrate why adopting a resilience factor $\lambda$ that
can adapt online is necessary. Given how much we want to weigh the
conservative control action of ”safety,” the correct level of resilience
depends on the attack condition. If the attack frequency is high, a
conservative solution is best. However, in sporadic attacks, the performance-
oriented solution is convenient. Table I shows the control cost variation
($\Delta J$ [-]) with respect to the minimum, between the MPC, Resilient MPC
(with fixed resilience factor $\lambda=\bar{\lambda})$ and Safe MPC control
strategies, for different attack frequencies. The results confirm the previous
discussion. Resilient MPC outperforms the competitors only in a ”in-between
scenario” where the attack frequency is of 1 attack every 2 weeks. This
motivates the need of resorting to an adaptive resilience factor
$\lambda=\lambda(t)$. Remark: $\bar{\lambda}=\lambda^{*}=50$ is computed after
sensitivity analysis, finding the best $\bar{\lambda}$ minimizing the
Resilient MPC control cost, as Figure 3 shows as best trade off solution for
the Pareto front. Remark: Control cost $J$ is the MPC control cost described
in Equation (1a), without the resilience term (e.g.,
$J=c_{f}||\textbf{s}-s_{max}||_{2}+c_{d}||\textbf{s}-s_{min}||_{2}+c||\epsilon||_{2}$)
evaluated a posteriori on system trajectories, for each control strategy.
Remark: $J_{demand}=c||u-d||_{2}$, $J_{dry}=c_{d}||\textbf{s}-s_{min}||_{2}$
are the dominating cost components used in Figure 3. Flood cost component is
equal for all the tested $\bar{\lambda}$ hence is not shown.

Figure 3: $\bar{\lambda}$ sensitivity analysis for fixed resilience, Resilient MPC. $\Delta J$ [-] | MPC | Resilient MPC | Safe MPC
---|---|---|---
No attacks | $0$ | $4943$ | $41845$
1 attack per week | $14645$ | $12516$ | $0$
1 attack every 2 weeks | $1488$ | $0$ | $9853$
TABLE I:  $\Delta$ Control Cost - MPC, Resilient MPC, Safe MPC - No attacks, 1
attack per week, 1 attack every 2 weeks.

###  II-E Hawkes-based Adaptive Resilience

In this section we show how we derived a time-varying resilience factor
adapting to varying attack conditions. At each time step, we train a Hawkes
process, we compute the estimated intensity function for different attack
return times and at last, we integrate the intensity function to retrieve the
resilience factor. The higher the probability of having attacks shortly, the
higher will be the resilience factor.

Each data (attack) vector on which the process is trained is a binary vector
($1$ if attack occurs, $0$ otherwise) of length equal to the simulated time
steps up to the current time. At each step, the attack vector is updated and
stored into database.

Hawkes processes are a class of point processes. The Hawkes choice depends on
their wide application in predicting the occurrence of random events in time.
In Hawkes, the intensity function is defined (according to [23]) as:

$\gamma(t)=\mu+\sum_{k}^{k(t_{k}<t)}\alpha\cdot e^{-\beta\cdot(t-t_{k})},$ (4)

where $t$ is current time, $t_{k}$ is the event time (in our case, attack
time), $\mu$ is the baseline, (e.g., the event rate in the absence of past
events), $\alpha$ is the weight associated with each event, $\beta$ is the
decay rate and controls the exponential decay over time. The collection
$\theta=[\mu,\alpha,\beta]$ constitutes the parameter vector that must be
identified during Hawkes training. $\theta^{*}=[\mu^{*},\alpha^{*},\beta^{*}]$
is identified through maximum likelihood estimation (MLE) on a moving window
$W=W^{*}$ of historical data. MLE identifies $\theta^{*}$ maximizing the joint
probability of the observed data, e.g., the likelihood function
$\mathcal{L}_{n}(\theta)$,

$\theta^{*}=\underset{\theta\in\Theta}{\operatorname{arg\;max}}\,\mathcal{L}_{n}(\theta;\mathbf{y}),$
(5)

where y is the collection of observed data samples and $\Theta$ is the
parameter space. The window size affects the quality of predicted intensity
function, however:

  * •

a larger window mandates a higher warm-up time in which the attack detection
process is not active.

  * •

a larger window slows down the adaptive capacity of the Hawkes process to
sudden attack pattern changes in time.

  * •

a shorter window leads to less data available for Hawkes process training,
increasing the probability of overfitting.

For the above reasons, our choice is a window size of $W^{*}=30$, fulfilling
the rule of thumb of a ratio of $10$ between the data numerousness and the
number of parameters, guaranteeing overfitting avoidance, and at the same
time, minimizing the warm-up time. After Hawkes process training, we can
compute the intensity function for different event (attack) return times. The
higher will be the intensity function on the corresponding horizon, the higher
the probability of having an attack within the given return time. Figure 4
shows the Hawkes intensity function (after training) for different estimated
return times, under different attack vectors (with different attack return
periods). The intensity function stays at its maximum for the length of the
expected return time, then it decreases. At each time step, the intensity
function is ”discretely integrated” to retrieve the resilience factor.

$\lambda(t)=\sum_{l}^{L_{MAX}}\gamma_{l}(t)$ (6)

where $l$ is the varying return time and $L_{MAX}=100$ is the maximum return
time considered. The choice of this hyperparameter does not affect Hawkes
prediction quality as it represents exclusively the number of lags considered
for the intensity function evaluation. The resilience factor (e.g., intensity
function integral) is higher as the estimated attack return time is shorter.
The Hawkes process fitting procedure is performed in an expanding window
fashion on simulation data if the current simulation time is lower than the
training window size (hence, the warm-up period is required). When the
simulation time is equal to or greater than the training size, only the last
$W^{*}$ elements of the window (including the current one) are read from
database and used for training. For perspective at a glance see the flow chart
in Figure 1 (bottom left for the attack detection components). The algorithm
is implemented in Python, thanks to the Machine Learning open-source library
Tick, which implements the Hawkes training procedure. See [24].

Algorithm 1 Adaptive Resilient MPC Algorithm

1:for t = 1,…,$T_{sim}$ do

2: Solve Attacker Problem with (3), return $u_{DoS}$.

3: Train Hawkes Process with (4), return $\gamma$.

4: Compute the resilience factor $\lambda$ with (6).

5: Solve SafeMPC with (2), return $u_{safe}$.

6: Solve ResMPC with (1), return $u_{res}$.

7: If attack occurs, apply $u_{DoS}$ to the system.

8: Else, apply $u_{res}$.

9: Update attack occurrence database. End for.

Algorithm 1 describes the Adaptive Resilient Model Predictive Control
procedure. $T_{sim}$ is the simulation length, in our case of $1$ year ($365$
days).

Figure 4: Hawkes intensity function as function of return time, for different
attack return times.

##  III Numerical Simulations and Results

In this section we show the numerical simulations and results obtained with
the 3 compared control strategies: MPC, Adaptive Resilient MPC, and Safe MPC
in 2 different attack scenarios.

In the first case, we consider a time-varying periodic attack return time. In
the second we tackle an attack occurrence prescribed by a pseudo random binary
signal (PRBS).

We perform the numerical simulations with real data of year 2005. All
computations are carried out on an Intel Core i7-8750H with 6 cores, at 2.20
GHz (maximum single core frequency: 4.10 GHz), with 16 GB RAM and running
Matlab R2021b. The control system runs in Matlab-Simulink, the MPC
optimization steps are implemented through the CVX solver. The optimal
attacker runs in Python 3.7 through the Gurobi solver (Gurobipy extension).
Hawkes algorithm runs in Python 3.7 through Tick statistical learning library.

The first attack scenario is a time varying periodic attack defined by the
following piece-wise conditions. This scenario is relatively easy to detect
but still, is representative for the effectiveness of the introduction of the
degree of freedom given by the time-varying resilience factor, and,
simultaneously, the capability of Hawkes process to adapt to changing
patterns.

$T=\begin{cases}7&\text{if }t<50\text{ or }t>250\\\ 14&\text{if }50\leq
t<150\\\ 21&\text{if }150\leq t\leq 250\end{cases}$ (7)

Equations (7) prescribe a periodic attack pattern with changing return period
$T$ of $1$, $2$, and $3$ weeks. Remark: simulation horizon is of $1$ year
($365$ days), daily resolution.

The second attack scenario is defined by a PRBS (pseudo random binary signal)
with clock time $\tau=10$ (e.g., the number of instants for which the signal
must stay constant) with the assumption of not having attacks in consecutive
days. This scenario is closer to reality as the attacker would like to be as
unpredictable as possible in the definition of the attack times. The choice of
a PRBS as triggering signal allows us to define a random signal with still a
pattern that could be tracked by the Hawkes process, given by the choice of
the clock period.

Figure 5 shows the time-varying resilience factor obtained through Equations
(4) and (6) in time-varying return period scenario (upper panel) and in PRBS
scenario (lower panel). In the upper panel, we highlight the required warm-up
window $W$ in which even when the first attack occurs the resilience factor
increases but takes a whole window to detect the first pattern ($7$ days).
When $T$ return period changes the process tracks the changing pattern with a
delay. This is the proof that, the Hawkes process can detect time-varying
periodic patterns.

Figure 5, lower panel, highlights how the resilience factor adapts also in
case of pseudo-random signal occurrence. Hawkes as soon as detects the
occurrence of the attacks increases the resilience level leading to a more
conservative control action. However, it modulates the resilience on the
estimated attack return period ”charging and discarging the integral” given by
Equation (6).

As Table II shows, in both the scenarios the proposed ARMPC architecture
outperforms the competitors ”performance-oriented” and ”safety-oriented”
strategies. Cost margin is higher in the time-varying return period scenario
as Hawkes detects with accuracy the changing return time of the attack. In the
pseudo-random case the advantage is still evident even if with lower margin.

Figure 5: Resilience Factor $\lambda$ and attack occurrence. $\Delta J$ Adaptive Resilient MPC [-] | MPC | Safe MPC
---|---|---
Time-varying Return Period | $6720$ | $1349$
PRBS | $7532$ | $119$
TABLE II:  $\Delta$ Control Cost with respect to Adaptive Resilient MPC - MPC,
Safe MPC - Time-varying Return Period, PRBS.

Figure 6 shows the Lake Como level in 2005 (upper panel) and water release
(lower panel) with the 3 control strategies under optimal attack with time-
varying return period attack pattern. SafeMPC is conservative. It always
increases lake level (reducing water release). Standard ”performance-oriented”
MPC works as if attacks do not happen. Adaptive Resilient MPC instead,
estimates the attack occurrence time and increases resilience level as the
return time gets closer. This is evident in correspondence of correctly
estimated return times where the lake level at first sticks to the lower
bound, then increases, for the resilience ”activation” and then decreases,
given the occurrence of the attack. In this way, it avoids lake level negative
peaks and, at the same time, tracks at best the agricultural demand. Under
attack conditions, all the strategies see $u_{DoS}$ control action (as the
attack ”spikes” underline in Figure 6, lower panel), given the current system
condition.

Figure 6: Upper panel: Lake Como level, lower panel: water release - 2005 -
Time-varying Attack Return Period.

##  IV CONCLUSIONS

The proposed adaptive resilience architecture (ARMPC) allows us to maximize
the performance of an MPC control scheme even in the presence of DoS attacks,
which inhibit the policy the controller decides and replace it with an optimal
attack, maximizing system damages. The use of Hawkes processes to estimate the
return time of the attack and the linking of its intensity function
(estimation of the probability of event occurrence) with the resilience factor
introduces a degree of freedom to the architecture that allows us to achieve
the same level of performance in normal conditions and to robustify the
control policy in attack conditions. We demonstrate the approach’s validity
with the attack’s occurrence given by a periodic return time with a variant
period and given by a pseudo-random signal. In future work, we will use the
same architecture in different applications, even on smaller-scale hardware
systems. We will also extensively validate the strategy over multiple years of
real system data.

## References

  * [1] Z. Drias, A. Serhrouchni, and O. Vogel, “Analysis of cyber security for industrial control systems,” in _2015 international conference on cyber security of smart cities, industrial control system and communications (ssic)_. IEEE, 2015, pp. 1–8.
  * [2] R. Gandhi, A. Sharma, W. Mahoney, W. Sousan, Q. Zhu, and P. Laplante, “Dimensions of cyber-attacks: Cultural, social, economic, and political,” _IEEE Technology and Society Magazine_ , vol. 30, no. 1, pp. 28–38, 2011\.
  * [3] M. Rudner, “Cyber-threats to critical national infrastructure: An intelligence challenge,” _International Journal of Intelligence and CounterIntelligence_ , vol. 26, no. 3, pp. 453–481, 2013.
  * [4] P. Lis and J. Mendel, “Cyberattacks on critical infrastructure: An economic perspective,” _Economics and Business Review_ , vol. 5, no. 2, pp. 24–47, 2019.
  * [5] E. Irmak and İ. Erkek, “An overview of cyber-attack vectors on scada systems,” in _2018 6th international symposium on digital forensic and security (ISDFS)_. IEEE, 2018, pp. 1–5.
  * [6] D. Du, R. Chen, X. Li, L. Wu, P. Zhou, and M. Fei, “Malicious data deception attacks against power systems: A new case and its detection method,” _Transactions of the Institute of Measurement and Control_ , vol. 41, no. 6, pp. 1590–1599, 2019.
  * [7] J. D. Markovic-Petrovic and M. D. Stojanovic, “Analysis of scada system vulnerabilities to ddos attacks,” in _2013 11th international conference on telecommunications in modern satellite, cable and broadcasting services (telsiks)_ , vol. 2. IEEE, 2013, pp. 591–594.
  * [8] R. Kalluri, L. Mahendra, R. S. Kumar, and G. G. Prasad, “Simulation and impact analysis of denial-of-service attacks on power scada,” in _2016 national power systems conference (NPSC)_. IEEE, 2016, pp. 1–5.
  * [9] Y. Liu, Y. Chen, M. Li, and Z. Wan, “Mpc for the cyber-physical system with deception attacks,” in _2020 chinese control and decision conference (ccdc)_. IEEE, 2020, pp. 3847–3852.
  * [10] M. R. Habibi, H. R. Baghaee, F. Blaabjerg, and T. Dragičević, “Secure mpc/ann-based false data injection cyber-attack detection and mitigation in dc microgrids,” _IEEE Systems Journal_ , vol. 16, no. 1, pp. 1487–1498, 2021.
  * [11] Q. Sun, K. Zhang, and Y. Shi, “Resilient model predictive control of cyber–physical systems under dos attacks,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 7, pp. 4920–4927, 2019.
  * [12] P. J. Laub, T. Taimre, and P. K. Pollett, “Hawkes processes,” _arXiv preprint arXiv:1507.02822_ , 2015.
  * [13] R. G. Cestari, F. Barchi, R. Busetto, D. Marazzina, and S. Formentin, “Hawkes-based cryptocurrency forecasting via limit order book data,” _arXiv preprint arXiv:2312.16190_ , 2023.
  * [14] E. Bacry, I. Mastromatteo, and J.-F. Muzy, “Hawkes processes in finance,” _Market Microstructure and Liquidity_ , vol. 1, no. 01, p. 1550005, 2015.
  * [15] P. Reynaud-Bouret, V. Rivoirard, and C. Tuleau-Malot, “Inference of functional connectivity in neurosciences via hawkes processes,” in _2013 IEEE global conference on signal and information processing_. IEEE, 2013, pp. 317–320.
  * [16] J. Kwon, Y. Zheng, and M. Jun, “Flexible spatio-temporal hawkes process models for earthquake occurrences,” _Spatial Statistics_ , vol. 54, p. 100728, 2023\.
  * [17] M.-A. Rizoiu, Y. Lee, S. Mishra, and L. Xie, “Hawkes processes for events in social media,” in _Frontiers of multimedia research_ , 2017, pp. 191–218.
  * [18] Y. Bessy-Roland, A. Boumezoued, and C. Hillairet, “Multivariate hawkes process for cyber insurance,” _Annals of Actuarial Science_ , vol. 15, no. 1, pp. 14–39, 2021.
  * [19] R. G. Cestari, A. Castelletti, and S. Formentin, “Hourly operation of a regulated lake via model predictive control,” _IFAC-PapersOnLine_ , vol. 55, no. 33, pp. 7–12, 2022.
  * [20] ——, “Scenario-based model predictive control of water reservoir systems,” _IFAC-PapersOnLine_ , vol. 56, no. 3, pp. 313–318, 2023.
  * [21] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
  * [22] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023. [Online]. Available: https://www.gurobi.com
  * [23] A. G. Hawkes, “Spectra of some self-exciting and mutually exciting point processes,” _Biometrika_ , vol. 58, no. 1, pp. 83–90, 1971.
  * [24] E. Bacry, M. Bompaire, S. Gaïffas, and S. Poulsen, “Tick: a python library for statistical learning, with a particular emphasis on time-dependent modelling,” 2018.