# Maximum likelihood identification of uncontrollable linear time-invariant
models for offset-free control††thanks: This report is an extended version of
a submitted paper. ††thanks: This work was supported by the National Science
Foundation (NSF) under Grant 2138985.

Steven J. Kuntz and James B. Rawlings The authors are with the Department of
Chemical Engineering, University of California, Santa Barbara, CA 93106 USA
(e-mail<EMAIL_ADDRESS>jbraw@ucsb.edu). _Corresponding author: Steven J.
Kuntz_.

###### Abstract

Maximum likelihood identification of linear time-invariant models is a
difficult problem because it is, in general, a nonlinear semidefinite program,
with semidefinite covariance matrix arguments and semidefinite filter
stability constraints. To enforce filter stability, we establish a general
theory of closed constraints on the system eigenvalues using LMI regions. To
solve the identification problem, we employ a Cholesky factorization method
that reduces the semidefinite program to a standard nonlinear program.
Finally, we apply the identification algorithm to a class of linear plant and
disturbance models commonly used in offset-free model predictive control
applications. Specifically, we consider models that are structured with
uncontrollable, integrating disturbance states. We solve this disturbance
modeling problem, and validate the resulting controller and estimator
performance, in two real-world case studies: first, a low-cost benchmark
temperature control laboratory, and second, an industrial-scale chemical
reactor at Eastman Chemical’s Kingsport plant.

##  I Introduction

Linear system identification is an important problem in control applications
and theory, with a longstanding history of applied use and a large body of
literature on its theory [1, 2, 3]. In particular, stochastic linear time-
invariant (LTI) state-space models are used in a variety of control contexts
to represent dynamics with process and measurement uncertainty. Maximum
likelihood (ML) identification is the preferable method in parametric
identification for its desirable statistical properties (consistency,
asymptotic efficiency) and ability to handle general parameterizations,
constraints, and stochastic noise models [4, 5, 6].

The main computational challenge to ML identification of LTI models is that,
in general, the problem is a nonlinear semidefinite program (SDP), with
semidefinite matrix arguments and semidefinite filter stability constraints.
In the ML identification literature, nonlinear SDPs are avoided either by
using the expectation maximization (EM) algorithm [7, 8, 9, 10], by invoking
simplified covariance matrix parameterizations (positive diagonal matrices,
scaled rotation matrices) [11, 12, 13, 14], or by minimizing the determinant
of the sample covariance [4, 5, 6, 15, 16, 17]. Each of these strategies
impose a specific structure on the estimates. The chief advantage of the EM
algorithm is that, for black-box models, there are closed-form solutions for
the iterates. When further structure on the model is imposed, it may invoke an
optimization problem within the EM algorithm iterates, significantly slowing
down the computation [18]. Diagonal covariances are a highly constrained
structure, and the type of simple covariance structures available in the
literature do not scale to high-dimensional systems. Finally, the minimum
determinant approach requires a fully parameterized covariance matrix for the
Kalman filter innovations. Moreover, these last two approaches introduce a
filter stability constraint that is not explicitly enforced in the current
literature.

EM does not have strong convergence guarantees even in the best case scenario.
While it can be shown that the EM iterates produce, almost surely, an
increasing sequence of likelihood values [7, 9], slow convergence at low noise
levels has been reported on a range of problems [18, 19, 20, 21, 22, 23], with
hundreds or thousands of iterates being common for small ML problems in even
the latest EM works [18]. Interior point, and even gradient methods [23], are
therefore preferable to the standard EM approach.

##### Linear identification of nonlinear systems

In a wide variety of control applications, including chemical processes [24,
25, 26], aerospace vehicles [17, 27], combustion engines [28], nautical
vehicles [11, 29], and speech recognition [8], linear approximations of the
nonlinear plant are beneficial for the convenience of linear identification
relative to that of nonlinear identification and the ability to meet strict
computational constraints, e.g., for on-line optimal control. Linear black-box
models are particularly useful when first-principles knowledge of the plant
dynamics is not available.

The main difficulty of linear identification of nonlinear systems is plant-
model mismatch. With ML identification, properties of the estimates are
dependent on the plant’s stochastic behavior [30, 31]. For stationary, input-
free models, the solution to the mismatched problem can be interpreted as
(asymptotically) minimizing the Kullback-Leibler divergence between the power
spectral densities of the model and plant [32]. However, there are still gaps
in the treatment of inputs, state-space models, and arbitrary nonlinear
plants. Rather than address this theory, we turn to control and estimation
applications that specifically address plant-model mismatch.

##### Application to offset-free control

Model predictive control (MPC) is a widely-used advanced control method in
which an optimal control problem is solved on-line, based on the current state
or state estimate, and the first input in the solution trajectory is injected
into the plant [33, 34]. Model quality is the main contributor to the
performance of an MPC implementation [35, 36], and therefore high-quality
identification algorithms are of relevance to MPC implementations.

Inherent to the use of identified plant models in MPC is the problem of plant-
model mismatch. Moreover, many applications require the rejection of a
stochastic, possibly nonstationary, disturbance process with unknown or un-
modeled components (e.g., environmental or upstream disturbances, demand
changes). In linear offset-free MPC, the stochastic LTI state-space model is
augmented with _uncontrollable_ integrating disturbance modes to achieve
offset-free control in the presence of plant-model mismatch and persistent
disturbances [37, 38]. We refer to these models as _linear augmented
disturbance models_.

Linear augmented disturbance models (or the resulting observer) can either be
tuned or identified. Tuning of disturbance models can be roughly divided into
three categories: pole placement [39, 40, 41, 42], diagonal covariance matrix
tuning [25, 43, 44], and direct filter gain tuning [45, 46, 47]. On the other
hand, disturbance models have been identified only via autocovariance least
squares (ALS) estimation [48] and (approximate) ML estimation [49, 50, 51].
Only our prior work in [50, 51] integrates the plant and disturbance
identification in a single step. However, the method in [50, 51] uses a nested
ML estimation approach in which unstructured stochastic LTI models are
augmented with integrating disturbance modes. The nested ML method can be
improved by consolidation to a single step ML step.

##### Direct data-driven control

The approach discussed so far is an indirect data-driven control design of
offset-free MPC. A potential alternative is the direct data-driven control
approach, where the control law is designed according to data [52, 53, 54,
55]. The drawback of this approach is its reliance on Willem’s Fundamental
Lemma [56], which does not admit the required linear augmented disturbance
model structure.

##### Contributions and outline

The main contributions of this work are (i) a method for directly solving
constrained ML identification problems as NLPs, more efficiently and on a
wider class of systems than the state-of-the-art, and (ii) real-world case
studies of the application of this algorithm to offset-free control. The ML
problem is stated in Section II. In Section III, stability and other
eigenvalue constraints are formulated. Inspired by the smoothed spectral
radius and abscissa formulations of [57, 58], and we present a novel barrier
function theory of constraints on the system eigenvalues via the linear matrix
inequality (LMI) regions of [59]. In Section IV, we present a novel
modification of the Burer-Monteiro-Zhang (BMZ) method [60], reformulating a
class of nonlinear SDPs as NLPs via Cholesky factorization. In identification
problems, Cholesky factor subsitution has only been used for replacing
semidefinite covariance matrices [49], not for general matrix inequalities. In
Section V, we used the ML identification algorithm for two real-world
applications of offset-free MPC: first, a benchmark temperature
microcontroller [61], and second, an industrial-scale chemical reactor at
Eastman Chemical’s Kingsport plant [51]. The advantage of ML identification
over other disturbance modeling techniques is demonstrated. Finally, in
Section VI, we conclude with a discussion of future work on applying the
identification problem to other parts of the control architecture (e.g.,
performance monitoring, steady-state optimization, and automated MPC upkeep)
and extensions to direct data-driven control.

This report is an extended version of a submitted work, and contains proofs of
minor results and details on the industrial-scale reactor case study that were
omitted from the journal version due to page limitations. Compared to the
journal version, this report additional contains the following additions:

  * •

the proof of Proposition 6 in Appendix A;

  * •

a longer version of Lemma 13;

  * •

an additional discussion of LMI region properties in Section III;

  * •

an explicit counterexample of [62, Thm. 1] in 18;

  * •

the proof of Proposition 20 in Appendix B;

  * •

the proof of Proposition 21(b,c) in Appendix C; and

  * •

additional figures, data, and discussion of the industrial-scale reactor study
in Section V.

##### Notation

Denote the set of $n\times n$ symmetric, positive definite, positive
semidefinite matrices, lower triangular, and positive lower triangular by
$\mathbb{S}^{n}$, $\mathbb{S}^{n}_{++}$, $\mathbb{S}^{n}_{+}$,
$\mathbb{L}^{n}$, and $\mathbb{L}^{n}_{++}$, respectively. Recall
$M\in{}^{n\times n}$ is positive definite if and only if there exists a unique
$L\in\mathbb{L}^{n}_{++}$, called the _Cholesky factor_ , such that
$M=LL^{\top}$. Denote the matrix direct sum and the Kronecker product by
$\oplus$ and $\otimes$, respectively, defined as in [63]. Define the set of
eigenvalues of a matrix $A\in{}^{n\times n}$ by $\lambda(A)\subset\mathbb{C}$.
The spectral radius and spectral abscissa are defined as
$\rho(A)\coloneqq\max_{\lambda\in\lambda(A)}|\lambda|$ and
$\alpha(A)\coloneqq\max_{\lambda\in\lambda(A)}\mathrm{Re}(\lambda)$,
respectively. We say a matrix $A$ is Schur (Hurwitz) stable if $\rho(A)<1$
($\alpha(A)<0$). We use $\sim$ as a shorthand for “distributed as” and
$\overset{\textnormal{iid}}{\sim}$ as a shorthand for “independent and
identically distributed as.” The complement, interior, closure, and boundary
of a set $S$ are denoted $S^{c}$, $\textnormal{int}(S)$, $\textnormal{cl}(S)$,
and $\partial S$, respectively.

##  II Maximum likelihood identification

In this section, we formulate a maximum likelihood problem for identifying
models of the following form:

$\displaystyle x_{k+1}$ $\displaystyle=A(\theta)x_{k}+B(\theta)u_{k}+w_{k}$
(1a) $\displaystyle y_{k}$ $\displaystyle=C(\theta)x_{k}+D(\theta)u_{k}+v_{k}$
(1b) $\displaystyle x_{0}$
$\displaystyle\sim\mathcal{N}(\hat{x}_{0}(\theta),\hat{P}_{0}(\theta))$ (1c)
$\displaystyle\begin{bmatrix}w_{k}\\\ v_{k}\end{bmatrix}$
$\displaystyle\overset{\textnormal{iid}}{\sim}\mathcal{N}(0,S(\theta))$ (1d)

where $x\in{}^{n}$ are the model states, $u\in{}^{m}$ are the inputs,
$y\in{}^{p}$ are the outputs, $w\in{}^{n}$ and $v\in{}^{p}$ are the process
and measurement noises, and
$\mathcal{M}\coloneqq(A,B,C,D,\hat{x}_{0},\hat{P}_{0},S)$ are functions, to be
defined, that map the model parameters $\theta\in\Theta$ to system matrices or
vectors of appropriate dimensions. The ML estimate $\hat{\theta}_{N}$ is
defined a maximizer of $p(\mathbf{y}_{N-1}|\mathbf{u}_{N-1},\theta)$, or
equivalently, a solution to

$\min_{\theta\in\Theta}f_{N}(\theta)\coloneqq-\sum_{k=0}^{N-1}\ln
p(y_{k}|\mathbf{u}_{k-1},\mathbf{y}_{k-1},\theta).$ (2)

The noise covariance matrix $S(\theta)$ may be partitioned as

$S(\theta)=\begin{bmatrix}Q_{w}(\theta)&S_{wv}(\theta)\\\
[S_{wv}(\theta)]^{\top}&R_{v}(\theta)\end{bmatrix}$ (3)

where $Q_{w}(\theta)\in\mathbb{S}^{n}_{+}$ is the process noise covariance,
$S_{wv}(\theta)$ is the cross-covariance, and
$R_{v}(\theta)\in\mathbb{S}^{p}_{+}$ is the measurement noise covariance.
Throughout, we impose the stronger requirement $R_{v}(\theta)\succ 0$ on the
measurement noise covariance.

###  II-A The constraint set

The main difficulty of the ML problem eq. 2 is that the parameter constraint
set $\Theta$ must necessarily contain matrix inequalities that make eq. 2 a
nonlinear SDP. For example, we expect consistency (i.e., positive
semidefiniteness) of the covariance matrices and nondegeneracy of the
measurement noise distribution.

######  Assumption 1.

For all $\theta\in\Theta$, we have $\hat{P}_{0}(\theta)\succeq 0$,
$S(\theta)\succeq 0$, and $R_{v}(\theta)\succ 0$.

Other matrix inequalities may arise as stability, eigenvalue, or other system-
level constraints. Stability and other eigenvalue constraints are explicitly
covered in Section III.

We structure the constraint set to facilitate the Cholesky factor-based
reformulation in Section IV. To define this structure, we first need to define
some additional notation. Consider the index sets
$\mathcal{L}^{n}\coloneqq\set{(i,j)\in\mathbb{N}^{2}}{1\leq i\leq j\leq n}$
and $\mathcal{D}^{n}\coloneqq\set{(i,i)\in\mathbb{N}^{2}}$ corresponding to
the sparsity patterns of $n\times n$ lower triangular and diagonal matrices,
respectively. With a slight abuse of notation, we define the direct sum of
index sets $\mathcal{I}\subseteq\mathcal{L}^{n}$ and
$\mathcal{J}\subseteq\mathcal{L}^{m}$ by

$\mathcal{I}\oplus\mathcal{J}\coloneqq\mathcal{I}\cup\set{(i+n,j+n)}{(i,j)\in\mathcal{L}^{m}}\subseteq\mathcal{L}^{n+m}.$

For each $\mathcal{I}\subseteq\mathcal{L}^{n}$, define the sets

$\displaystyle\mathbb{S}^{n}[\mathcal{I}]$
$\displaystyle\coloneqq\set{S\in\mathbb{S}^{n}}{S_{ij}=0\;\forall\;(i,j)\in\mathcal{L}^{n}\setminus\mathcal{I}}$
$\displaystyle\mathbb{L}^{n}[\mathcal{I}]$
$\displaystyle\coloneqq\set{L\in\mathbb{L}^{n}}{L_{ij}=0\;\forall\;(i,j)\not\in\mathcal{I}}$
$\displaystyle\mathbb{L}^{n}_{++}[\mathcal{I}]$
$\displaystyle\coloneqq\set{L\in\mathbb{L}^{n}_{++}}{L_{ij}=0\;\forall\;(i,j)\not\in\mathcal{I}}.$

Using this notation, we make the following assumptions about the structure of
the parameter constraint set $\Theta$.

######  Assumption 2.

There exist index sets
$\mathcal{D}^{n_{\Sigma}}\subseteq\mathcal{I}_{\Sigma}\subseteq\mathcal{L}^{n_{\Sigma}}$
and
$\mathcal{D}^{n_{\mathcal{A}}}\subseteq\mathcal{I}_{\mathcal{A}}\subseteq\mathcal{L}^{n_{\mathcal{A}}}$
and differentiable functions
$g:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow{}^{n_{g}}$,
$h:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow{}^{n_{h}}$,
$H:{}^{n_{\beta}}\rightarrow\mathbb{S}^{n_{\Sigma}}$, and
$\mathcal{A}:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow\mathbb{S}^{n_{\mathcal{A}}}[\mathcal{I}_{\mathcal{A}}]$
such that

$\Theta=\\{\;(\beta,\Sigma)\in{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]\;|\;g(\beta,\Sigma)=0,\\\
h(\beta,\Sigma)\leq 0,\;\Sigma\succeq
H(\beta),\;\mathcal{A}(\beta,\Sigma)\succeq 0\;\\}.$

Moreover, $\textnormal{cl}(\Theta_{++})=\Theta$ where

$\Theta_{++}\coloneqq\\{\;(\beta,\Sigma)\in{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]\;|\;g(\beta,\Sigma)=0,\\\
h(\beta,\Sigma)\leq 0,\;\Sigma\succ H(\beta),\;\mathcal{A}(\beta,\Sigma)\succ
0\;\\}.$

Under 2, the parameter set $\Theta$ contains four constraints: two standard
vector equality and inequality constraints and two matrix inequalities. The
purpose of the elimination algorithm is to transform the matrix inequalities
into vector equalities while introducing as few new variables as possible. The
first matrix inequality enforces a lower bound $H(\beta)$ on the sparse
symmetric matrix argument $\Sigma$. In many cases this lower bound is simply
zero or a small diagonal matrix, but we have left the bound general for
illustrative purposes. This inequality is in the form used by [60], where
variables in the Cholesky factorization

$\displaystyle\Sigma$ $\displaystyle=L_{\Sigma}L_{\Sigma}^{\top}+H(\beta),$
$\displaystyle L_{\Sigma}$ $\displaystyle\in\mathbb{L}^{n_{\Sigma}}_{++}$ (4)

are algorithmically eliminated to write $\Sigma$ in terms of just $H(\beta)$
and a sparse lower triangular matrix
$L^{\mathcal{I}_{\Sigma}}\in\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]$.
Similarly, in Section IV, we use the sparsity structure in the second matrix
inequality $\mathcal{A}(\beta,\Sigma)$ to algorithmically eliminate variables
in the squared slack variable transformation

$\displaystyle\mathcal{A}(\beta,\Sigma)$
$\displaystyle=L_{\mathcal{A}}L_{\mathcal{A}}^{\top},$ $\displaystyle
L_{\mathcal{A}}$ $\displaystyle\in\mathbb{L}^{n_{\mathcal{A}}}_{++}$ (5)

by writing $L_{\mathcal{A}}$ in terms of a sparse lower triangular matrix
$L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]$.
The second part of 2 guarantees the existence and uniqueness of these Cholesky
factors $(L_{\Sigma},L_{\mathcal{A}})$, helps to avoid divisions by zero
during the variable elimination procedure, and allows taking limits.

######  Remark 3.

2 rules out direct use of the strict inequality $R_{v}(\theta)\succ 0$. To
satisfy 1, we use the closed constraint $R_{v}(\theta)\succeq\delta I_{p}$
with a small backoff $\delta>0$.

######  Remark 4.

The index set $\mathcal{I}_{\Sigma}$ provides enough sparsity in $\Sigma$ to
reduce the problem to $n_{\beta}+|\mathcal{I}_{\Sigma}|$ variables (excluding
slack terms). This is important when $\Sigma$ has a particular structure. The
most common example is the block diagonal structure $\Sigma=\hat{P}_{0}\oplus
Q_{w}\oplus R_{v}\in\mathbb{S}^{2n+p}[\mathcal{I}_{\Sigma}]$ where
$\mathcal{I}_{\Sigma}\coloneqq\mathcal{L}^{n}\oplus\mathcal{L}^{n}\oplus\mathcal{L}^{p}$.
We may further restrict $Q_{w}$ and $R_{v}$ to take block tridiagonal and
diagonal structures, e.g.,

$\displaystyle Q_{w}$ $\displaystyle=\begin{bmatrix}Q_{1,1}&Q_{1,2}\\\
Q_{1,2}^{\top}&Q_{2,2}&\ddots\\\ &\ddots&\ddots&Q_{\tilde{n}-1,\tilde{n}}\\\
&&Q_{\tilde{n}-1,\tilde{n}}^{\top}&Q_{\tilde{n},\tilde{n}}\end{bmatrix}$
$\displaystyle R_{v}$ $\displaystyle=R_{1}\oplus\ldots\oplus R_{\tilde{n}}$

that arise in sequentially interconnected processes such as chemical plants.
Adding a $Q_{1,\tilde{n}}$ block can account for an overall recycle loop. Note
that if we parameterize the block tridiagonal $Q_{w}$ via a sparse shaping
matrix (i.e., $Q_{w}=G_{w}G_{w}^{\top}$), then there are more parameters than
if the sparsity of $Q_{w}$ is known _unless_ the rank of $G_{w}$ is known to
be low.

######  Remark 5.

The index set $\mathcal{I}_{\mathcal{A}}$ provides sparsity on the range of
$\mathcal{A}$ to eliminate all but $|\mathcal{I}_{\mathcal{A}}|$ entries of
the squared slack term $L_{\mathcal{A}}$. Typically $\mathcal{A}$ is a block
diagonal matrix of eigenvalue constraints, so the dense matrix
$L_{\mathcal{A}}$ is reduced to a sparse matrix
$L_{\mathcal{A}}=L_{\mathcal{A},1}\oplus\ldots\oplus
L_{\mathcal{A},{\tilde{n}_{\mathcal{A}}}}$.

###  II-B Kalman filtering and the log-likelihood

In this subsection, we derive exact and approximate expressions for the
negative log-likelihood. For brevity, we drop the dependence on $\theta$ where
appropriate and write

$\displaystyle x_{k+1}$ $\displaystyle=Ax_{k}+Bu_{k}+w_{k}$ $\displaystyle
y_{k}$ $\displaystyle=Cx_{k}+Du_{k}+v_{k}$ $\displaystyle x_{0}$
$\displaystyle\sim\mathcal{N}(\hat{x}_{0},\hat{P}_{0})$
$\displaystyle\begin{bmatrix}w_{k}\\\ v_{k}\end{bmatrix}$
$\displaystyle\overset{\textnormal{iid}}{\sim}\mathcal{N}(0,S)$

and $\mathcal{M}=(A,B,C,D,\hat{x}_{0},\hat{P}_{0},S)$.

####  II-B1 Time-varying Kalman filter formulation

Consider the Kalman filter in innovations form

$\displaystyle\hat{x}_{k+1}$
$\displaystyle=A\hat{x}_{k}+Bu_{k}+\mathcal{K}_{k}e_{k}$ (6a) $\displaystyle
y_{k}$ $\displaystyle=C\hat{x}_{k}+Du_{k}+e_{k}$ (6b) $\displaystyle e_{k}$
$\displaystyle\sim\mathcal{N}(0,\mathcal{R}_{k})\quad\textnormal{(indep.)}$
(6c) where $\displaystyle\hat{P}_{k+1}$ $\displaystyle\coloneqq
A\hat{P}_{k}A^{\top}+Q_{w}-\mathcal{K}_{k}\mathcal{R}_{k}\mathcal{K}_{k}^{\top}$
(6d) $\displaystyle\mathcal{K}_{k}$
$\displaystyle\coloneqq(A\hat{P}_{k}C^{\top}+S_{wv})\mathcal{R}_{k}^{-1}$ (6e)
$\displaystyle\mathcal{R}_{k}$ $\displaystyle\coloneqq
C\hat{P}_{k}C^{\top}+R_{v}.$ (6f)

Since the $e_{k}$ in eq. 6 are mutually independent, we have the negative log-
likelihood

$f_{N}(\theta)\propto\frac{1}{2}\sum_{k=0}^{N-1}\ln\det\mathcal{R}_{k}(\theta)+|e_{k}(\theta)|_{[\mathcal{R}_{k}(\theta)]^{-1}}^{2}$

and we can write eq. 2 equivalently as

$\min_{\theta\in\Theta}\frac{1}{2}\sum_{k=0}^{N-1}\ln\det\mathcal{R}_{k}(\theta)+|e_{k}(\theta)|_{[\mathcal{R}_{k}(\theta)]^{-1}}^{2}$
(7)

where the $e_{k}(\theta)$ and $\mathcal{R}_{k}(\theta)$ are given by the
recursion eq. 6. For black-box covariance models

$\mathcal{M}(\beta,\hat{P}_{0}\oplus
S)=(A(\beta),B(\beta),C(\beta),D(\beta),\hat{x}_{0}(\beta),\hat{P}_{0},S)$

the covariance consistency constraint

$\Theta\coloneqq\set{(\beta,\Sigma)\in{}^{n_{\beta}}\times\mathbb{S}^{n+p}[\mathcal{L}^{n}\oplus\mathcal{L}^{n+p}]}{\Sigma\succeq\begin{bmatrix}0\\\
&\varepsilon I_{p}\end{bmatrix}}$

suffices, without any vector constraints or general matrix inequalities. These
constraints say nothing about stability of the filter eq. 6, so unstable
filters may be realized during optimization, producing numerically infinite or
undefined values.

####  II-B2 Steady-state Kalman filter formulation

In most situations, the state error covariance matrix converges exponentially
fast to a steady-state solution $\hat{P}_{k}\rightarrow\hat{P}$, so it
suffices to consider the following steady-state filter:

$\displaystyle\hat{x}_{k+1}$ $\displaystyle=A\hat{x}_{k}+Bu_{k}+Ke_{k}$ (8a)
$\displaystyle y_{k}$ $\displaystyle=C\hat{x}_{k}+Du_{k}+e_{k}$ (8b)
$\displaystyle e_{k}$
$\displaystyle\overset{\textnormal{iid}}{\sim}\mathcal{N}(0,R_{e})$ (8c)

where $K\coloneqq(A\hat{P}C^{\top}+S_{wv})R_{e}^{-1}$, $R_{e}\coloneqq
C\hat{P}C^{\top}+R_{v}$, and $\hat{P}$ is the unique, stabilizing solution to
the discrete algebraic Riccati equation (DARE),

$\hat{P}=A\hat{P}A^{\top}+Q_{w}-(A\hat{P}C^{\top}+S_{wv})\\\
\times(C\hat{P}C^{\top}+R_{v})^{-1}(A\hat{P}C^{\top}+S_{wv})^{\top}.$ (9)

Recall a solution to the DARE eq. 9 is stabilizing if the resulting
$A_{K}\coloneqq A-KC$ is stable.

Convergence of $\hat{P}_{k}$ to $\hat{P}$ is equivalent to the solution to the
DARE eq. 9 being unique and stabilizing. We generally assume such a solution
exists, but for completeness, we state the following proposition, adapted from
[64, Thm. 18(iii)] (see Appendix A for proof).

######  Proposition 6.

Assume $R_{v}\succ 0$ and consider the full rank factorization

$\begin{bmatrix}Q_{w}&S_{wv}\\\
S_{wv}^{\top}&R_{v}\end{bmatrix}=\begin{bmatrix}\tilde{B}\\\
\tilde{D}\end{bmatrix}\begin{bmatrix}\tilde{B}^{\top}&\tilde{D}^{\top}\end{bmatrix}$

Then the following statements are equivalent:

  1. 1.

The DARE eq. 9 has a unique, stabilizing solution $\hat{P}\succeq 0$.

  2. 2.

The error covariance converges exponentially fast
$\hat{P}_{k}\rightarrow\hat{P}$ for any $\hat{P}_{0}\succeq 0$.

  3. 3.

$(A,C)$ is detectable and $(A-FC,\tilde{B}-F\tilde{D})$ is stabilizable for
all $F\in{}^{n\times p}$.

Under the steady-state approximation, the likelihood is

$\tilde{f}_{N}(\theta)\coloneqq\frac{N}{2}\ln\det
R_{e}(\theta)+\frac{1}{2}\sum_{k=0}^{N-1}|e_{k}(\theta)|_{[R_{e}(\theta)]^{-1}}^{2}$

and we can approximate solutions to eq. 2 by solving

$\min_{\theta\in\Theta}\frac{N}{2}\ln\det
R_{e}(\theta)+\frac{1}{2}\sum_{k=0}^{N-1}|e_{k}(\theta)|_{[R_{e}(\theta)]^{-1}}^{2}$
(10)

where the $e_{k}(\theta)$ are given by the recursion eq. 8 and $R_{e}(\theta)$
is found by solving the DARE eq. 9.

While $R_{e}(\theta)$ and $K(\theta)$ could be defined via $\hat{P}(\theta)$,
taken as the function that returns solutions to the DARE eq. 9 and therefore
enforcing filters stability, it is more convenient to directly parameterize
these matrices. In fact, it is equivalent to consider the model structure

$\mathcal{M}=\left(A,B,C,D,\hat{x}_{0},0,\begin{bmatrix}KR_{e}K^{\top}&KR_{e}\\\
R_{e}K^{\top}&R_{e}\end{bmatrix}\right)$

where $K:\Theta\rightarrow{}^{n\times p}$ and
$R_{e}:\Theta\rightarrow\mathbb{S}^{p}_{++}$ are now given functions that are
not explicitly related to the DARE eq. 9. We write this model parameterization
as

$\mathcal{M}_{\textnormal{KF}}=(A,B,C,D,\hat{x}_{0},K,R_{e})$ (11)

using the subscript KF to denote that $\mathcal{M}_{\textnormal{KF}}$
represents a Kalman filter in innovation form eq. 11. Since the model
structure eq. 11 no longer enforces filter stability, the following stability
assumption is required.

######  Assumption 7.

Given the model structure eq. 11, we have $\rho(A_{K}(\theta))<1$ for all
$\theta\in\Theta$.

Since $\rho$ is continuous but not differentiable, 7 cannot be directly
implemented in a form satisfying 2, even with a backoff term to make the
inequality nonstrict, i.e., $\rho(A_{K}(\theta))\leq 1-\delta$ where
$\delta>0$. We deal with enforcing the stability constraint in Section III.

####  II-B3 Minimum determinant formulation

Suppose, in the model structure eq. 11, that $R_{e}$ is parameterized fully,
and separately from the other terms, i.e.,

$\mathcal{M}_{\textnormal{KF}}(\beta,\tilde{\Sigma}\oplus
R_{e})=(A(\beta,\tilde{\Sigma}),B(\beta,\tilde{\Sigma}),C(\beta,\tilde{\Sigma}),\\\
D(\beta,\tilde{\Sigma}),\hat{x}_{0}(\beta,\tilde{\Sigma}),K(\beta,\tilde{\Sigma}),R_{e}).$

Moreover, assume $R_{e}$ is constrained separately as well, i.e.,

$\Theta=\\{\;(\beta,\tilde{\Sigma}\oplus
R_{e})\in{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]\;|\;\tilde{g}(\beta,\tilde{\Sigma})=0,\\\
\tilde{h}(\beta,\tilde{\Sigma})\leq
0,\;\tilde{\Sigma}\succeq\tilde{H}(\beta),\;\;R_{e}\succeq\varepsilon
I_{p},\;\tilde{\mathcal{A}}(\beta,\tilde{\Sigma})\succeq 0\;\\}$

for some sufficiently small $\varepsilon>0$, differentiable functions
$(\tilde{g},\tilde{h},\tilde{H},\tilde{\mathcal{A}})$, and index set
$\mathcal{I}_{\Sigma}\coloneqq\mathcal{I}_{\tilde{\Sigma}}\oplus\mathcal{L}^{p}$
where
$\mathcal{D}^{\tilde{n}_{\Sigma}}\subseteq\mathcal{I}_{\tilde{\Sigma}}\subseteq\mathcal{L}^{\tilde{n}_{\Sigma}}$
and $n_{\tilde{\Sigma}}\coloneqq n_{\Sigma}-p$. Then we can always solve eq.
10 stagewise, first in $R_{e}$, and then in the remaining variables
$(\beta,\tilde{\Sigma})$. Solving the inner problem gives the solution

$\hat{R}_{e}(\beta,\tilde{\Sigma})\coloneqq\frac{1}{N}\sum_{k=0}^{N-1}e_{k}(\beta,\tilde{\Sigma})[e_{k}(\beta,\tilde{\Sigma})]^{\top}$

where we use the fact that $e_{k}$ is only dependent on
$(\beta,\tilde{\Sigma})$, and we assume
$\hat{R}_{e}(\beta,\tilde{\Sigma})\succeq\varepsilon I_{p}$ for all
$(\beta,\tilde{\Sigma})\in\tilde{\Theta}$ where

$\tilde{\Theta}\coloneqq\\{\;(\beta,\tilde{\Sigma})\in{}^{n_{\beta}}\times\mathbb{S}^{\tilde{n}_{\Sigma}}[\mathcal{I}_{\tilde{\Sigma}}]\;|\;\tilde{g}(\beta,\tilde{\Sigma})=0,\\\
\tilde{h}(\beta,\tilde{\Sigma})\leq
0,\;\tilde{\Sigma}\succeq\tilde{H}(\beta),\;\tilde{\mathcal{A}}(\beta,\tilde{\Sigma})\succeq
0\;\\}.$

The outer problem can be written

$\min_{(\beta,\tilde{\Sigma})\in\tilde{\Theta}}\det\hat{R}_{e}(\beta,\tilde{\Sigma}).$
(12)

The problem eq. 12 is of relevance for avoiding the nonlinear SDP formulation
of eq. 10, both in the early ML identification literature [4, 5, 6] and in
recent works [15, 16, 17]. None of these works consider filter stability
constraints. To the best of our knowledge, only [18] consider the ML problem
eq. 7 with stability constraints, but they consider open-loop stability (i.e.,
$\rho(A)<1$) and use the EM algorithm. To satisfy 7, we must consider filter
stability constraints (i.e., $\rho(A_{K})<1$).

######  Remark 8.

For real-world data, $\det\hat{R}_{e}(\beta,\tilde{\Sigma})=0$ is not
attainable because that would imply some direction of $y_{k}$ were perfectly
modeled. Therefore, a constant $\varepsilon>0$ exists such that the lower
bound $\hat{R}_{e}(\beta,\tilde{\Sigma})\succeq\varepsilon I_{p}$ is satisfied
for all $(\beta,\tilde{\Sigma})\in\tilde{\Theta}$. Moreover, we need not
explicitly choose $\varepsilon>0$ to satisfy the constraint $R_{e}\succ 0$ of
1.

##  III Eigenvalue constraints

In this section, we describe stability and other eigenvalue constraints
conforming to 2. Notably, we consider the smoothed spectral radius and
abscissa formulations of [57, 58], the linear matrix inequality (LMI) region
approach of [59], and a novel barrier function method that combines the two
approaches. While we consider the matrix $A$ in this section, constraints can
be added to any square matrix of interest, such as $A$, $A_{K}$, or some
submatrix thereof.

###  III-A Spectral radius and abscissa bounds

Recall that, for any $Q\in\mathbb{S}^{n}_{++}$, a matrix $A$ is Schur stable
if and only if there exists $P\in\mathbb{S}^{n}_{++}$ such that

$P-APA^{\top}=Q.$ (13)

Likewise, for any $Q\in\mathbb{S}^{n}_{++}$, a matrix $A$ is Hurwitz stable if
and only if there eixsts $P\in\mathbb{S}^{n}_{++}$ such that

$AP+PA^{\top}=-Q.$ (14)

Moreover, solutions to eqs. 13 and 14, when they exist, are uniquely given by
$P_{d}(A,Q)\coloneqq\sum_{k=0}^{\infty}A^{k}Q(A^{\top})^{k}$ and
$P_{c}(A,Q)\coloneqq\int_{t=0}^{\infty}e^{At}Qe^{A^{\top}t}dt$, respectively.
It is a well-known observation that discrete- and continuous-time stability of
$A$ are equivalent to finiteness of the matrices $P_{d}(A,Q)$ and
$P_{c}(A,Q)$, respectively.

######  Lemma 9 ([57, Lem. 5.1] and [58, Lem. 2.1]).

Let $Q\in\mathbb{S}^{n}_{++}$ and $\|\cdot\|:{}^{n\times
n}\rightarrow\mathbb{R}_{\geq 0}$ be a submultiplicative norm. A matrix
$A\in{}^{n\times n}$ is Schur stable (Hurwitz stable) if and only if
$\|P_{d}(A,Q)\|$ ($\|P_{c}(A,Q)\|$) is finite.

Inspired by this observation, [57, 58] impose an upper bound on the norm of
$P$ satisfying either eq. 13 or eq. 14.

####  III-A1 Smoothed spectral radius

Let $W,V\in\mathbb{S}^{n}_{++}$ and consider the function
$\phi_{d}(A,s)\coloneqq\textnormal{tr}(VP_{d}(A/s,W))$. In [57], the implicit
function theorem is used to show the existence of a _smoothed spectral radius_
$\rho_{\varepsilon}(A)$ satisfying

$\phi_{d}(A,\rho_{\varepsilon}(A))=\varepsilon^{-1}.$ (15)

Properties of $\rho_{\varepsilon}(A)$ are reiterated in the following theorem.

######  Theorem 10 ([57, Thms. 5.4, 5.6]).

There exists a function
$\rho_{(\cdot)}(\cdot):\mathbb{R}_{>0}\times{}^{n\times n}\rightarrow\real$
such that, for each $A\neq 0$ and $\varepsilon>0$, $\rho_{\varepsilon}(A)$
uniquely solves eq. 15, and $\rho_{\varepsilon}(0)=0$. Moreover,
$\rho_{(\cdot)}(\cdot)$ has the following properties:

  1. 1.

$\rho_{(\cdot)}(\cdot)$ is analytic on $\mathbb{R}_{>0}\times{}^{n\times
n}\setminus\set{0}$ and continuous on $\mathbb{R}_{>0}\times\set{0}$;

  2. 2.

$\rho_{\varepsilon}(A)>\rho(A)$ for all $A\neq 0$ and $\varepsilon>0$;

  3. 3.

$\rho_{\varepsilon}(A)\searrow\rho(A)$ as $\varepsilon\searrow 0$ for all
$A\in{}^{n\times n}$;

  4. 4.

for each $\varepsilon,s>0$, $\rho_{\varepsilon}(A)\leq s$ if and only if there
exists $P\succeq 0$ such that $s^{2}P-APA^{\top}=W$ and
$\textnormal{tr}(VP)\leq\varepsilon^{-1}$.

The first property of Theorem 10 establishes the smoothness of
$\rho_{\varepsilon}$. The second and fourth properties of Theorem 10 let us
construct constraint sets satisfying 2 and 7. Finally, the third property
demonstrates that $\rho_{\varepsilon}$ is an approximation of $\rho$ in the
sense that $\varepsilon$ can be made arbitrarily small to keep
$\rho_{\varepsilon}(A)-\rho(A)$ arbitrarily small.

####  III-A2 Smoothed spectral abscissa

Let $W,V\in\mathbb{S}^{n}_{++}$ and consider the function
$\phi_{c}(A,s)\coloneqq\textnormal{tr}(VP_{c}(A-sI,W))$. In [58], the implicit
function theorem is used to show the existence of a _smoothed spectral
abscissa_ $\alpha_{\varepsilon}(A)$ satisfying

$\phi_{c}(A,\alpha_{\varepsilon}(A))=\varepsilon^{-1}.$ (16)

Properties of $\alpha_{\varepsilon}(A)$ are reiterated in the following
theorem.

######  Theorem 11 ([58, Thms. 2.5, 2.6]).

There exists a function
$\alpha_{(\cdot)}(\cdot):\mathbb{R}_{>0}\times{}^{n\times n}\rightarrow\real$
such that, for each $A\in{}^{n\times n}$ and $\varepsilon>0$,
$\alpha_{\varepsilon}(A)$ uniquely solves eq. 15. Moreover,
$\alpha_{(\cdot)}(\cdot)$ has the following properties:

  1. 1.

$\alpha_{(\cdot)}(\cdot)$ is analytic on $\mathbb{R}_{>0}\times{}^{n\times
n}$;

  2. 2.

$\alpha_{\varepsilon}(A)>\alpha(A)$ for all $A\in{}^{n\times n}$ and
$\varepsilon>0$;

  3. 3.

$\alpha_{\varepsilon}(A)\searrow\alpha(A)$ as $\varepsilon\searrow 0$ for all
$A\in{}^{n\times n}$; and

  4. 4.

for each $\varepsilon,s>0$, $\alpha_{\varepsilon}(A)\leq s$ if and only if
there exists $P\succeq 0$ such that $(A-sI)P+P(A-sI)^{\top}=-W$ and
$\textnormal{tr}(VP)\leq\varepsilon^{-1}$.

Similarly to the properties of Theorem 10, the first property of Theorem 11
gives $\alpha_{\varepsilon}$ its smoothness property, the second and fourth
properties allow us to construct constraints that satisfy 2, and the third
property demonstrates that $\alpha_{\varepsilon}$ is an approximation of
$\alpha$.

###  III-B $\mathcal{D}$-stability constraints

We may wish to place the eigenvalues in a specified region. In [59] the
problem of placing eigenvalues in any convex, open region of the complex plane
is posed in terms of solving a linear matrix inequality (LMI). These so-called
“LMI regions” are defined as follows.

######  Definition 12.

A subset $\mathcal{D}\subseteq\mathbb{C}$ that takes the form

$\mathcal{D}=\set{z\in\mathbb{C}}{f_{\mathcal{D}}(z)\succ 0}$

is called an _LMI region_ with the _characteristic function_
$f_{\mathcal{D}}:\mathbb{C}\rightarrow\mathbb{C}^{m\times m}$, defined as

$f_{\mathcal{D}}(z)\coloneqq M_{0}+M_{1}z+M_{1}^{\top}\overline{z}$

and _generating matrices_ $(M_{0},M_{1})\in\mathbb{S}^{n}\times{}^{n\times
n}$.

The following lemma defines the four basic LMI regions: shifted half-planes,
circles centered on the real axis, conic sections, and horizontal bands.

######  Lemma 13.

For each $s,x_{0}\in\real$, the subsets

$\displaystyle\mathcal{D}_{1}$
$\displaystyle\coloneqq\set{z\in\mathbb{C}}{\textnormal{Re}(z)<s}$
$\displaystyle\mathcal{D}_{2}$
$\displaystyle\coloneqq\set{z\in\mathbb{C}}{|z-x_{0}|<s}$
$\displaystyle\mathcal{D}_{3}$
$\displaystyle\coloneqq\set{z\in\mathbb{C}}{|\textnormal{Im}(z)|<s(\textnormal{Re}(z)-x_{0})}$
$\displaystyle\mathcal{D}_{4}$
$\displaystyle\coloneqq\set{z\in\mathbb{C}}{|\textnormal{Im}(z)|<s}$

are LMI regions with characteristic functions

$\displaystyle f_{\mathcal{D}_{1}}(z)$ $\displaystyle\coloneqq
2s-z-\overline{z}$ $\displaystyle f_{\mathcal{D}_{2}}(z)$
$\displaystyle\coloneqq\begin{bmatrix}s&-x_{0}\\\
-x_{0}&s\end{bmatrix}+\begin{bmatrix}0&1\\\
0&0\end{bmatrix}z+\begin{bmatrix}0&0\\\ 1&0\end{bmatrix}\overline{z}$
$\displaystyle f_{\mathcal{D}_{3}}(z)$
$\displaystyle\coloneqq-2sx_{0}I_{2}+\begin{bmatrix}s&1\\\
-1&s\end{bmatrix}z+\begin{bmatrix}s&-1\\\ 1&s\end{bmatrix}\overline{z}$
$\displaystyle f_{\mathcal{D}_{4}}(z)$
$\displaystyle\coloneqq-2sI_{2}+\begin{bmatrix}0&1\\\
-1&0\end{bmatrix}z+\begin{bmatrix}0&-1\\\ 1&0\end{bmatrix}\overline{z}.$

###### Proof.

The first identity follows from the formula
$2\textnormal{Re}(z)=z+\overline{z}$. For the second identity, we have
$f_{\mathcal{D}_{2}}(z)=\begin{bmatrix}s&z-x_{0}\\\
\overline{z}-x_{0}&s\end{bmatrix}\succ 0$ if and only if $s>0$ and
$s^{2}>|z-x_{0}|^{2}$, or equivalently, $|z-x_{0}|<s$. For the third identity,
we have
$f_{\mathcal{D}_{3}}(z)=\begin{bmatrix}2s(\textnormal{Re}(z)-x_{0})&2\iota\textnormal{Im}(z)\\\
-2\iota\textnormal{Im}(z)&2s(\textnormal{Re}(z)-x_{0})\end{bmatrix}\succ 0$ if
and only if $2s(\textnormal{Re}(z)-x_{0})>0$ and
$4s^{2}(\textnormal{Re}(z)-x_{0})^{2}>4|\textnormal{Im}(z)|^{2}$, or
equivalently, $|\textnormal{Im}(z)|<s(\textnormal{Re}(z)-x_{0})$. For the
fourth identity, we have
$f_{\mathcal{D}_{4}}(z)=\begin{bmatrix}2s&2\iota\textnormal{Im}(z)\\\
-2\iota\textnormal{Im}(z)&2s\end{bmatrix}\succ 0$ if and only if $2s>0$ and
$4s^{2}>4|\textnormal{Im}(z)|^{2}$, or equivalently, $|\textnormal{Im}(z)|<s$.
∎

######  Remark 14.

For continuous-time systems, $\mathcal{D}_{1}$ corresponds to a minimum decay
rate of $s>0$, $\mathcal{D}_{3}$ corresponds to a minimum damping ratio
$-\cos(\theta)$, and $\mathcal{D}_{2}\cap\mathcal{D}_{3}$ implies to a maximum
natural frequency $r\sin(\theta)$, where $\theta=\tan^{-1}(s)$ and $s<0$ [59].
For discrete-time systems, $\mathcal{D}_{2}$ corresponds to a minimum decay
rate of $-\ln r$, and $\mathcal{D}_{2}\cap\mathcal{D}_{3}$ implies a minimum
damping ratio $-\cos(\tan^{-1}(\theta/\ln r))$ and maximum natural frequency
$(\ln(r)^{2}+\theta^{2})/\Delta$, where $\theta=\tan^{-1}(s)$, $s>0$, and
$\Delta$ is the sample time.

######  Remark 15.

An LMI region $\mathcal{D}$ is convex, open, and symmetric about the imaginary
axis. The intersection of two LMI regions
$\mathcal{D}\coloneqq\mathcal{D}_{1}\cap\mathcal{D}_{2}$ is an LMI region with
the characteristic function $f_{\mathcal{D}}(z)=f_{\mathcal{D}_{1}}(z)\oplus
f_{\mathcal{D}_{2}}(z)$. An LMI region $\mathcal{D}$ with characteristic
function $f_{\mathcal{D}}$ also has characteristic function
$Mf_{\mathcal{D}}(\cdot)M^{\top}$ for any nonsingular $M\in{}^{m\times m}$. By
this property, we can construct any convex polyhedron that is symmetric about
the real axis by intersecting left and right half-planes, horizontal strips,
and conic sections. Moreover, since any convex region can be approximated, to
any desired accuracy, by a convex polyhedron, the set of LMI regions is dense
in the space of convex subsets of $\mathbb{C}$ that are symmetric about the
real axis. For an in-depth discussion of LMI region geometry and other
properties, see [65].

Throughout, assume the LMI region $\mathcal{D}$ is nonempty, not equal to
$\mathbb{C}$, and its characteristic function $f_{\mathcal{D}}$ and generating
matrices $(M_{0},M_{1})$ are fixed. We seek LMI conditions under which the
eigenvalues of $A\in{}^{n\times n}$ lie in $\mathcal{D}$ or
$\textnormal{cl}(\mathcal{D})$.

####  III-B1 $\mathcal{D}$-stability

In the following definition, we generalize the notion of asymptotic stability
to include pole placement within a given LMI region $\mathcal{D}$.

######  Definition 16.

We say the matrix $A\in{}^{n\times n}$ is _$\mathcal{D}$ -stable_ if
$\lambda(A)\subset\mathcal{D}$.

In [59], $\mathcal{D}$-stability of a matrix $A\in{}^{n\times n}$ is shown to
be equivalent to the strict feasibility of the following system of matrix
inequalities:

$\displaystyle M_{\mathcal{D}}(A,P)$ $\displaystyle\succ 0,$ $\displaystyle P$
$\displaystyle\succ 0$ (17)

where $M_{\mathcal{D}}:{}^{n\times
n}\times\mathbb{S}^{n}\rightarrow\mathbb{S}^{nm}$ is defined by

$M_{\mathcal{D}}(A,P)\coloneqq M_{0}\otimes
P+M_{1}\otimes(AP)+M_{1}^{\top}\otimes(AP)^{\top}.$ (18)

This fact is restated in the following theorem.

######  Theorem 17 ([59, Thm. 2.2]).

The matrix $A\in{}^{n\times n}$ is $\mathcal{D}$-stable if and only if eq. 17
holds for some $P\in\mathbb{S}^{n}$.

####  III-B2 Marginal $\mathcal{D}$-stability

The drawback of Theorem 17 is strictness of the matrix inequalities eq. 17,
making them inadmissible in constraint sets satisfying 2. Suppose we relax the
first inequality,

$\displaystyle M_{\mathcal{D}}(A,P)$ $\displaystyle\succeq 0,$ $\displaystyle
P$ $\displaystyle\succ 0.$ (19)

Since $M_{\mathcal{D}}(A,P)$ is linear in $P$, feasibility of eq. 19 is
equivalent to feasibility of

$\displaystyle M_{\mathcal{D}}(A,P)$ $\displaystyle\succeq 0,$ $\displaystyle
P$ $\displaystyle\succeq P_{0}$ (20)

for some fixed $P_{0}\in\mathbb{S}^{n}_{++}$.111For any $P_{0}\succ 0$ and $P$
satisfying eq. 19, define the scaling factor
$\gamma\coloneqq\|P_{0}\|_{2}\|P^{-1}\|_{2}$ and a rescaled solution
$P^{*}\coloneqq\gamma P$. Then $P^{*}\succeq P_{0}$ and
$M_{\mathcal{D}}(A,P^{*})=\gamma M_{\mathcal{D}}(A,P)\succeq 0$. Therefore the
system eq. 19 is admissible to constraint sets satisfying 2 with minor
alteration.

An attempt was made in [62, Thm. 1] to characterize the class of matrices
$A\in{}^{n\times n}$ for which eq. 19, but this theorem does not correctly
treat eigenvalues on the boundary $\partial\mathcal{D}$. We restate [62, Thm.
1] below as a conjecture and disprove it with a simple counterexample.

######  Conjecture 18 ([62, Thm. 1]).

The matrix $A\in{}^{n\times n}$ satisfies
$\lambda(A)\subset\textnormal{cl}(\mathcal{D})$ if and only if eq. 19 holds
for some $P\in\mathbb{S}^{n}$.

###### Counterexample.

Let $\mathcal{D}$ be the left half-plane, consider the Jordan block
$A=\begin{bmatrix}0&1\\\ 0&0\end{bmatrix}$, and suppose
$P=\begin{bmatrix}p_{11}&p_{12}\\\
p_{12}&p_{22}\end{bmatrix}\in\mathbb{S}^{2}$ such that eq. 19 holds. Then
$\lambda(A)\subset\textnormal{cl}(\mathcal{D})$ and

$0\preceq M_{\mathcal{D}}(A,P)=-\begin{bmatrix}2p_{12}&p_{22}\\\
p_{22}&0\end{bmatrix}$

which implies $p_{12}=p_{22}=0$, a contradiction of eq. 19. ∎

The correction to 18 requires a more careful treatment of eigenvalues lying on
the the LMI region’s boundary $\partial\mathcal{D}$. Recall that a matrix $A$
is marginally stable (in continuous- or discrete-time) if all its eigenvalues
lie in the closure of the stability region (left half-plane or unit disc)
_and_ the non-simple eigenvalues lie strictly in the interior of that region.
Replacing the stability region with $\mathcal{D}$, we propose the following
definition of _marginal $\mathcal{D}$-stability_.

######  Definition 19.

We say a matrix $A\in{}^{n\times n}$ is _marginally $\mathcal{D}$-stable_ if
$\lambda(A)\subseteq\textnormal{cl}(\mathcal{D})$ and $\lambda\in\mathcal{D}$
for all non-simple eigenvalues $\lambda\in\lambda(A)$.

In the following proposition, we show $\mathcal{D}$-stability of
$A\in{}^{n\times n}$ is equivalent to feasibility of eq. 19 (see Appendix B
for proof).

######  Proposition 20.

The matrix $A\in{}^{n\times n}$ is marginally $\mathcal{D}$-stable if and only
if eq. 19 holds for some $P\in\mathbb{S}^{n}$.

####  III-B3 $\mathcal{D}$-stability barrier functions

Thus far, we have proposed LMI equivalences without regard to the topology of
the set of $\mathcal{D}$-stable and marginally $\mathcal{D}$-stable matrices

$\displaystyle\mathbb{A}_{\mathcal{D}}^{n}$
$\displaystyle\coloneqq\set{A\in{}^{n\times n}}{\textnormal{$A$ is
$\mathcal{D}$-stable}}$ $\displaystyle\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$
$\displaystyle\coloneqq\set{A\in{}^{n\times n}}{\textnormal{$A$ is marginally
$\mathcal{D}$-stable}}.$

The following proposition characterizes the topology of
$\mathbb{A}_{\mathcal{D}}^{n}$ and $\mathbb{A}_{\mathcal{D}}^{n}$ (see
Appendix C for proof).

######  Proposition 21.

  1. (a)

$\mathbb{A}_{\mathcal{D}}^{n}$ is open.

  2. (b)

$\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$ is not open if (i) $n\geq 2$ or (ii)
$\partial\mathcal{D}\cap\real$ is nonempty.

  3. (c)

$\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$ is not closed if (i) $n\geq 4$ or (ii)
$\partial\mathcal{D}\cap\real$ is nonempty and $n\geq 2$.

  4. (d)

$\textnormal{cl}(\mathbb{A}_{\mathcal{D}}^{n})=\set{A\in{}^{n\times
n}}{\lambda(A)\subset\textnormal{cl}(\mathcal{D})}$.

Proposition 21 reveals a weakness of the marginal $\mathcal{D}$-stability
constraints eqs. 19 and 20. Since $\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$ is
not closed and $P\in\mathbb{S}^{n}_{+}$ is neither bounded nor regularized,
the optimizer may fail to converge in $P$ as it seeks an $A$ that is not a
limit point in $\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$. While the optimizer can
approach $A$ in this case, $P$ grows unbounded along the path of iterates, and
the optimizer itself does not converge.

Motivated by the smoothed spectral radius and abscissa, we upper bound
$\textnormal{tr}(VP)$ for some $V\in\mathbb{S}^{n}_{++}$. This still allows
rescaling of $P$, however, so we also lower bound $M_{\mathcal{D}}(A,P)$,

$\displaystyle M_{\mathcal{D}}(A,P)$ $\displaystyle\succeq M,$ $\displaystyle
P$ $\displaystyle\succeq 0,$ $\displaystyle\textnormal{tr}(VP)$
$\displaystyle\leq\varepsilon^{-1}$ (21)

and assume $M\in\mathbb{S}^{nm}_{+}$ is chosen in a way that implies eq. 17.
Consider the parameterized linear SDP,

$\phi_{\mathcal{D}}(A)\coloneqq\inf_{P\in\mathbb{S}^{n}_{+}}\textnormal{tr}(VP)\;\textnormal{subject
to}\;M_{\mathcal{D}}(A,P)\succeq M.$ (22)

The optimal value function $\phi_{\mathcal{D}}:{}^{n\times
n}\rightarrow\mathbb{R}_{\geq 0}\cup\set{\infty}$ is a barrier function for
the constraint $A\in\mathbb{A}_{\mathcal{D}}^{n}$. Proposition 22 establishes
properties of $\phi_{\mathcal{D}}$ and its $\varepsilon^{-1}$-sublevel sets
(see Appendix D for proof).

######  Proposition 22.

Let $V\in\mathbb{S}^{n}_{++}$ and $M\in\mathbb{S}^{n}_{+}$ such that
$M_{\mathcal{D}}(A,P)\succeq M$ implies $M_{\mathcal{D}}(A,P)\succ 0$. Then

  1. (a)

$\phi_{\mathcal{D}}$ is continuous on $\mathbb{A}_{\mathcal{D}}$;

  2. (b)

for each $\varepsilon>0$, the $\varepsilon^{-1}$-sublevel set of
$\phi_{\mathcal{D}}$,

$\displaystyle\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)\coloneqq$
$\displaystyle\set{A\in{}^{n\times
n}}{\phi_{\mathcal{D}}(A)\leq\varepsilon^{-1}}$ (23) $\displaystyle=$
$\displaystyle\set{A\in{}^{n\times n}}{\exists P\succeq
0:\textnormal{\lx@cref{creftype~refnum}{eq:lmi:cont} holds}}$ (24)

is closed; and

  3. (c)

$\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)\nearrow\mathbb{A}_{\mathcal{D}}^{n}$
as $\varepsilon\searrow 0$.

######  Remark 23.

It is sufficient, but not necessary, to choose a positive definite lower bound
$M\succ 0$.

######  Remark 24.

To reconstruct the $s$-sublevel sets of the smoothed spectral radius via
Proposition 22, we set $M=sW\oplus 0_{n\times n}$ for any $W,V\succ 0$ and
$s>0$ and apply the Schur complement lemma to
$M_{\mathcal{D}_{2}}(A,P)/s-M/s$, where $\mathcal{D}_{2}$ is the circle
defined in Lemma 13 with $x_{0}=0$, and $M_{\mathcal{D}_{2}}$ is defined by
the generating matrices used in Lemma 13. Then the $\varepsilon^{-1}$-sublevel
set of $\phi_{\mathcal{D}_{2}}$ equals the $s$-sublevel set of
$\rho_{\varepsilon}$.

######  Remark 25.

The $\varepsilon^{-1}$-sublevel sets of $\phi_{\mathcal{D}_{1}}$ equal the
$s$-sublevel sets of $\alpha_{\varepsilon}$, where $\mathcal{D}_{1}$ is the
shifted half-plane defined in Lemma 13, and $M=W$ for any $W,V\succ 0$.

###  III-C Discussion

We conclude this section with a discussion on implementing eigenvalue
constraints. Discrete LTI models with eigenvalues having positive real parts
have a one-to-one correspondence with continuous LTI models [66]. As such, it
is important to satisfy not only the filter stability constraint
$\rho(A_{K})<1$ (and, if desired, open-loop stability constraint $\rho(A)<1$),
but it is also desirable to satisfy differentiability constraints
$\alpha(-A)<0$ and $\alpha(-A_{K})<0$. In practice, we find
$\rho_{\varepsilon}(A_{K})\leq 1-\delta$ and $\alpha_{\varepsilon}(-A_{K})\leq
0$ to be important constraints in ML identification of linear augmented
disturbance models. Without the former constraint, the optimizer will
frequently pick unstable filters that evaluate to infinite objective values,
and without the latter constraint, the optimal filter equations degenerate to
aphysical solutions that incorrectly estimate $u\mapsto y$ gain matrices.
Last, we refer the reader to [67] for other examples of matrix inequalities
(e.g., detectability, minimum phase) that may be useful in system
identification.

##  IV Cholesky reparameterizations

In this section, we seek to transform the ML problems eqs. 2, 7 and 10 from
nonlinear SDPs to standard NLPs, while introducing as few new variables as
possible. The main idea is to combine a squared slack variable substitution
eqs. 4 and 5 with an elimination scheme that adds only
$|\mathcal{I}_{\mathcal{A}}|$ variables to the optimization problem.

For this section, we define the following notation. For each
$\mathcal{I}\subseteq\mathcal{L}^{n}$, let $\pi_{\mathcal{I}}^{L}:{}^{n\times
n}\rightarrow\mathbb{L}^{n}[\mathcal{I}]$ and $\pi_{\mathcal{I}}:{}^{n\times
n}\rightarrow\mathbb{S}^{n}[\mathcal{I}]$ denote the orthogonal projections
(in the Frobenius norm) from n×n onto the subspaces
$\mathbb{L}^{n}[\mathcal{I}]$ and $\mathbb{S}^{n}[\mathcal{I}]$, respectively.
Let $\textnormal{chol}:\mathbb{S}^{n}_{++}\rightarrow\mathbb{L}^{n}_{++}$
denote the invertible function that maps a positive definite matrix to its
Cholesky factor.

###  IV-A Burer-Monteiro-Zhang method

In the Burer-Monteiro-Zhang (BMZ) method [60], an invertible map is
constructed between the interior of

$\Theta_{\textnormal{BMZ}}\coloneqq\set{(\beta,\Sigma)\in{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]}{\Sigma\succeq
H(\beta)}$ (25)

and

$\Theta_{\textnormal{Chol}}\coloneqq{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]$
(26)

where $H:{}^{n_{\beta}}\rightarrow\mathbb{S}^{n_{\Sigma}}$ and
$\mathcal{I}_{\Sigma}$ is some index set satisfying
$\mathcal{D}^{n_{\Sigma}}\subseteq\mathcal{I}_{\Sigma}\subseteq\mathcal{L}^{n_{\Sigma}}$.
The set $\Theta_{\textnormal{BMZ}}$ satisfies 2 so long as $H$ is
differentiable since

$\textnormal{int}(\Theta_{\textnormal{BMZ}})=\set{(\beta,\Sigma)\in{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]}{\Sigma\succ
H(\beta)}.$

Recall $\Sigma\succ H(\beta)$ if and only if
$\Sigma=L_{\Sigma}L_{\Sigma}^{\top}+H(\beta)$ for some
$L_{\Sigma}\in\mathbb{L}^{n_{\Sigma}}_{++}$. With
$\mathcal{J}_{\Sigma}\coloneqq\mathcal{L}^{n_{\Sigma}}\setminus\mathcal{I}_{\Sigma}$,
we split $L_{\Sigma}$ into the sum of
$L^{\mathcal{I}_{\Sigma}}\in\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]$
and
$L^{\mathcal{J}_{\Sigma}}\in\mathbb{L}^{n_{\Sigma}}[\mathcal{J}_{\Sigma}]$,

$\Sigma=(L^{\mathcal{I}_{\Sigma}}+L^{\mathcal{J}_{\Sigma}})(L^{\mathcal{I}_{\Sigma}}+L^{\mathcal{J}_{\Sigma}})^{\top}+H.$
(27)

But $\Sigma\in\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]$, so we can take

$\pi_{\mathcal{J}_{\Sigma}}[(L^{\mathcal{I}_{\Sigma}}+L^{\mathcal{J}_{\Sigma}})(L^{\mathcal{I}_{\Sigma}}+L^{\mathcal{J}_{\Sigma}})^{\top}+H]=0$
(28)

to produce $|\mathcal{J}_{\Sigma}|$ equations to eliminate the
$|\mathcal{J}_{\Sigma}|$ variables of $L^{\mathcal{J}_{\Sigma}}$. For each
$(H,L^{\mathcal{I}_{\Sigma}})\in\mathbb{S}^{n_{\Sigma}}\times\mathbb{L}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]$,
define

$L^{\mathcal{J}_{\Sigma}}_{ij}\coloneqq-\frac{1}{L_{jj}^{\mathcal{I}_{\Sigma}}}\left(H_{ij}+\sum_{k=1}^{j-1}(L^{\mathcal{I}_{\Sigma}}_{ik}+L^{\mathcal{J}_{\Sigma}}_{ik})(L^{\mathcal{I}_{\Sigma}}_{jk}+L^{\mathcal{J}_{\Sigma}}_{jk})\right)$
(29)

for each $(i,j)\in\mathcal{J}_{\Sigma}$, in a top-to-bottom and left-to-right
order. So long as we never divide by zero, each $L^{\mathcal{J}_{\Sigma}}$ is
fully defined by $H$ and $L^{\mathcal{I}_{\Sigma}}$ via eq. 29, and we have
the following lemma.

######  Lemma 26 ([60, Lem. 1]).

For each
$(\beta,L^{\mathcal{I}_{\Sigma}})\in{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]$
such that $L^{\mathcal{I}_{\Sigma}}_{ii}\neq 0$ for each
$i\in\mathbb{I}_{1:n_{\Sigma}}$, there is a unique
$L^{\mathcal{J}_{\Sigma}}\in\mathbb{L}^{n_{\Sigma}}[\mathcal{J}_{\Sigma}]$
satisfying eq. 28.

Let
$L^{\mathcal{J}_{\Sigma}}:{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]\rightarrow\mathbb{L}^{n_{\Sigma}}[\mathcal{J}_{\Sigma}]$
be the function that maps each
$(\beta,L^{\mathcal{I}_{\Sigma}})\in\Theta_{\textnormal{Chol}}$ to the
$L^{\mathcal{J}_{\Sigma}}\in\mathbb{L}^{n_{\Sigma}}[\mathcal{J}_{\Sigma}]$
defined by eq. 29 with $H=H(\beta)$. Similarly, let
$\Sigma:{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]\rightarrow\mathbb{S}^{n_{\Sigma}}[\mathcal{I}_{\Sigma}]$
be the map defined by eq. 27 for each
$(\beta,L^{\mathcal{I}_{\Sigma}})\in\Theta_{\textnormal{Chol}}$ with
$H=H(\beta)$ and
$L^{\mathcal{J}_{\Sigma}}=L^{\mathcal{J}_{\Sigma}}(\beta,L^{\mathcal{I}_{\Sigma}})$.
Finally, we let

$\Psi(\beta,L^{\mathcal{I}_{\Sigma}})\coloneqq\left(\beta,\Sigma(\beta,L^{\mathcal{I}_{\Sigma}})\right)$
(30)

which has the inverse

$\Psi^{-1}(\beta,\Sigma)\coloneqq\left(\beta,\pi_{\mathcal{I}_{\Sigma}}^{L}[\textnormal{chol}(\Sigma-H(\beta))]\right)$
(31)

and we have the following lemma.

######  Lemma 27 ([60, Lem. 2]).

The function $\Psi$ defined by eq. 30 is a bijection between
$\Theta_{\textnormal{Chol}}$ and
$\textnormal{int}(\Theta_{\textnormal{BMZ}})$.

Differentiability of $\Psi$ and $\Psi^{-1}$ follow from differentiability of
$H$ and the algorithm eq. 29. In fact, these functions are as smooth as $H$
is, so if $H$ is analytic, so are $\Psi$ and $\Psi^{-1}$. More importantly,
the bijection $\Psi$ allows us to transform the minimum of a continuous
function in $\Theta_{\textnormal{BMZ}}$ to an infimum of a function in
$\Theta_{\textnormal{Chol}}$, given by the following theorem.

######  Theorem 28 ([60, Thm. 1]).

For any continuous function
$f:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow\real$ that attains a
minimum in $\Theta_{\textnormal{BMZ}}$,

$\min_{(\beta,\Sigma)\in\Theta_{\textnormal{BMZ}}}f(\beta,\Sigma)=\inf_{(\beta,L^{\mathcal{I}_{\Sigma}})\in\Theta_{\textnormal{Chol}}}f_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})$
(32)

where $f_{\Psi}\coloneqq f\circ\Psi$.

We reiterate the proof of Theorem 28 for illustrative purposes.

###### Proof.

Continuity of $f$ implies its minimum over $\Theta_{\textnormal{BMZ}}$ equals
its infimum over $\textnormal{int}(\Theta_{\textnormal{BMZ}})$, i.e.,

$\min_{(\beta,\Sigma)\in\Theta_{\textnormal{BMZ}}}f(\beta,\Sigma)=\inf_{(\beta,\Sigma)\in\textnormal{int}(\Theta_{\textnormal{BMZ}})}f(\beta,\Sigma)$

Since $\Psi$ is a bijection, we can transform the optimization variables as
follows:

$\inf_{(\beta,\Sigma)\in\textnormal{int}(\Theta_{\textnormal{BMZ}})}f(\beta,\Sigma)=\inf_{(\beta,L^{\mathcal{I}_{\Sigma}})\in\Psi^{-1}(\textnormal{int}(\Theta_{\textnormal{BMZ}}))}f(\Psi(\beta,L^{\mathcal{I}_{\Sigma}}))$

and
$\Theta_{\textnormal{Chol}}=\Psi^{-1}(\textnormal{int}(\Theta_{\textnormal{BMZ}}))$
and $f_{\Psi}=f\circ\Psi$ imply eq. 32. ∎

###  IV-B Modified Burer-Monteiro-Zhang method

For constraint sets $\Theta$ satisfying 2, the BMZ method only suffices to
eliminate the variables $\Sigma$ and the matrix inequality $\Sigma\succeq
H(\beta)$. The addition of the general matrix inequality
$\mathcal{A}(\beta,\Sigma)\succeq 0$ requires a similar procedure to introduce
squared slack variables, although the elimination procedure will no longer
eliminate as many variables as are introduced. We further complicate the
generalization of Lemma 27 and Theorem 28 by requiring general vector equality
and inequality constraints $g(\beta,\Sigma)=0$ and $h(\beta,\Sigma)\leq 0$.

Lemma 27 gives the invertible map
$\Psi:\Theta_{\textnormal{Chol}}\rightarrow\Theta_{\textnormal{BMZ}}$ with
which we can define

$\displaystyle g_{\Psi}$ $\displaystyle\coloneqq g\circ\Psi,$ $\displaystyle
h_{\Psi}$ $\displaystyle\coloneqq h\circ\Psi,$
$\displaystyle\mathcal{A}_{\Psi}$
$\displaystyle\coloneqq\mathcal{A}\circ\Psi.$

In this notation, consider the set

$\Theta_{\textnormal{NLP}}\coloneqq\\{\,(\beta,L^{\mathcal{I}_{\Sigma}})\in{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]\;|\;g_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})=0,\\\
h_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})=0,\;\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\succ
0\,\\}.$ (33)

Recall $\Sigma(\beta,L^{\mathcal{I}_{\Sigma}})\succ H(\beta)$ holds
automatically for all
$(\beta,L^{\mathcal{I}_{\Sigma}})\in\Theta_{\textnormal{Chol}}\coloneqq{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]$.
Then by a simple change of variables we have
$\Psi^{-1}(\Theta_{++})=\Theta_{\textnormal{NLP}}$. Since
$\Theta_{\textnormal{NLP}}$ and $\Theta_{++}$ are subsets of
$\Theta_{\textnormal{Chol}}$ and
$\textnormal{int}(\Theta_{\textnormal{BMZ}})$, respectively, we have proven
the following lemma.

######  Lemma 29.

The function $\Psi$ defined by eq. 30 is a bijection between
$\Theta_{\textnormal{NLP}}$ and $\Theta_{++}$.

Next, we remove the matrix inequality
$\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\succ 0$ from our
representation of the constraint set $\Theta_{\textnormal{NLP}}$. Recall
$\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\succ 0$ if and only if
$\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})=L_{\mathcal{A}}L_{\mathcal{A}}^{\top}$
for some $L_{\mathcal{A}}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}$. With
$\mathcal{J}_{\mathcal{A}}\coloneqq\mathcal{L}^{n_{\mathcal{A}}}\setminus\mathcal{I}_{\mathcal{A}}$,
we split $L_{\mathcal{A}}$ into the sum of
$L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]$
and
$L^{\mathcal{J}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}[\mathcal{J}_{\mathcal{A}}]$,

$\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})=(L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}})(L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}})^{\top}.$
(34)

But
$\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\in\mathbb{S}^{n_{\mathcal{A}}}[\mathcal{I}_{\mathcal{A}}]$,
so we can take

$\pi_{\mathcal{J}_{\mathcal{A}}}[(L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}})(L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}})^{\top}]=0$

to produce $|\mathcal{J}_{\mathcal{A}}|$ equalities with which to eliminate
the $|\mathcal{J}_{\mathcal{A}}|$ free variables in
$L^{\mathcal{J}_{\mathcal{A}}}$. For each
$L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]$,
we define

$L^{\mathcal{J}_{\mathcal{A}}}_{ij}\coloneqq-\frac{1}{L_{jj}^{\mathcal{I}_{\mathcal{A}}}}\sum_{k=1}^{j-1}(L^{\mathcal{I}_{\mathcal{A}}}_{ik}+L^{\mathcal{J}_{\mathcal{A}}}_{ik})(L^{\mathcal{I}_{\mathcal{A}}}_{jk}+L^{\mathcal{J}_{\mathcal{A}}}_{jk})$
(35)

for each $(i,j)\in\mathcal{J}_{\mathcal{A}}$, in a top-to-bottom and left-to-
right order. Each $L^{\mathcal{J}_{\mathcal{A}}}$ satisfying eq. 34 is thus
fully defined by $L^{\mathcal{I}_{\mathcal{A}}}$ via eq. 35, and we have the
following corollary of Lemma 26.

######  Corollary 30.

For each
$L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}[\mathcal{I}_{\mathcal{A}}]$
such that $L^{\mathcal{I}_{\mathcal{A}}}_{ii}\neq 0$ for each
$i\in\mathbb{I}_{1:n_{\mathcal{A}}}$, there exists a unique
$L^{\mathcal{J}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}[\mathcal{J}_{\mathcal{A}}]$
satisfying eq. 34.

Let
$L^{\mathcal{J}_{\mathcal{A}}}:\mathbb{L}^{n_{\mathcal{A}}}_{++}\rightarrow\mathbb{L}^{n_{\mathcal{A}}}[\mathcal{I}_{\mathcal{A}}]$
denote the map described by Corollary 30, and let
$\tilde{\mathcal{A}}(L^{\mathcal{I}_{\mathcal{A}}})\coloneqq(L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}}(L^{\mathcal{I}_{\mathcal{A}}}))(L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}}(L^{\mathcal{I}_{\mathcal{A}}}))^{\top}$.
Then $\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\succ 0$ if and only
if
$\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})=\tilde{\mathcal{A}}(L^{\mathcal{I}_{\mathcal{A}}})$
for some
$L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]$.
Subsuming this constraint $g_{\Psi}$, we have

$\tilde{g}(\beta,L^{\mathcal{I}_{\Sigma}},L^{\mathcal{I}_{\mathcal{A}}})\coloneqq\begin{bmatrix}g_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\\\
\textnormal{vecs}_{\mathcal{I}_{\mathcal{A}}}\left(\mathcal{A}_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})-\tilde{\mathcal{A}}(L^{\mathcal{I}_{\mathcal{A}}})\right)\end{bmatrix}$

for all
$(\beta,L^{\mathcal{I}_{\Sigma}},L^{\mathcal{I}_{\mathcal{A}}})\in{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]\times\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]$,
where
$\textnormal{vecs}_{\mathcal{I}_{\mathcal{A}}}:\mathbb{S}^{n_{\mathcal{A}}}\rightarrow{}^{|\mathcal{I}_{\mathcal{A}}|}$
vectorizes the $|\mathcal{I}_{\mathcal{A}}|$ entries of the argument
corresponding to the index set $\mathcal{I}_{\mathcal{A}}$. We have proven the
following lemma.

######  Lemma 31.

The set $\Theta_{\textnormal{NLP}}$, defined by eq. 33, equals

$\\{\,(\beta,L^{\mathcal{I}_{\Sigma}})\in{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]\;|\;\exists
L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]:\\\
\tilde{g}(\beta,L^{\mathcal{I}_{\Sigma}},L^{\mathcal{I}_{\mathcal{A}}})=0,\;h_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}}))\leq
0\,\\}.$

Finally, we have the following equivalence between minimization problems over
the parameter sets $\Theta$ and $\Theta_{\textnormal{NLP}}$.

######  Proposition 32.

For any continuous function
$f:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow\real$ that attains a
minimum in $\Theta$,

$\min_{(\beta,\Sigma)\in\Theta}f(\beta,\Sigma)=\inf_{(\beta,L^{\mathcal{I}_{\Sigma}})\in\Theta_{\textnormal{NLP}}}f_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})$

where $f_{\Psi}\coloneqq f\circ\Psi$.

###### Proof.

The proof follows that of Theorem 28, noting that 2 gives
$\textnormal{cl}(\Theta_{++})=\Theta$ and therefore the minimum of $f$ over
$\Theta$ equals the infimum of $f$ over $\Theta_{++}$. ∎

###  IV-C $\varepsilon$-approximate solutions

The strict inequalities implied by constraints
$\mathbb{L}^{n_{\Sigma}}_{\varepsilon}[\mathcal{I}_{\Sigma}]$ and
$\mathbb{L}^{n_{\mathcal{A}}}_{\varepsilon}[\mathcal{I}_{\mathcal{A}}]$ are
not amenable to implementation in standard NLP software. In [60], a log-
barrier approach is used to achieve global convergence for a class of linear
SDPs. Instead, we consider a constant backoff on the inequalities, providing a
small but nonzero lower bound on the diagonal elements. For each
$\mathcal{D}^{n}\subseteq\mathcal{I}\subseteq\mathcal{L}^{n}$, we define

$\mathbb{L}^{n}_{\varepsilon}[\mathcal{I}]\coloneqq\set{L\in\mathbb{L}^{n}_{++}[\mathcal{I}]}{L_{ii}\geq\varepsilon\;\forall
i\in\mathbb{I}_{1:n}}.$

These sets imply an additional $n$ inequality constraints, but do not have any
ill-posed strictness requirements. Therefore we can optimize over the
restricted set

$\Theta_{\varepsilon}\coloneqq\\{\,(\beta,L^{\mathcal{I}_{\Sigma}})\in{}^{n_{\beta}}\times\mathbb{L}^{n_{\Sigma}}_{\varepsilon}[\mathcal{I}_{\Sigma}]\;|\;\exists
L^{\mathcal{I}_{\mathcal{A}}}\in\mathbb{L}^{n_{\mathcal{A}}}_{\varepsilon}[\mathcal{I}_{\mathcal{A}}]:\\\
\tilde{g}(\beta,L^{\mathcal{I}_{\Sigma}},L^{\mathcal{I}_{\mathcal{A}}})=0,\;h_{\Psi}(\beta,L^{\mathcal{I}_{\Sigma}})\leq
0\,\\}.$ (36)

In the following proposition we show, for any continuous function $f$, the
infimum of $f$ over $\Theta_{\varepsilon}$ converges to the minimum of $f$
over $\Theta$, so long as 2 is satisfied (see Appendix E for proof).

######  Proposition 33.

For any continuous
$f:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow\real$ that attains a
minimum in $\Theta$, let
$\mu_{0}\coloneqq\min_{(\beta,\Sigma)\in\Theta}f(\beta,\Sigma)$ and

$\mu_{\varepsilon}\coloneqq\inf_{(\beta,L^{\mathcal{I}_{\Sigma}})\in\Theta_{\varepsilon}}f(\beta,L^{\mathcal{I}_{\Sigma}})$
(37)

where $f(\beta,L^{\mathcal{I}_{\Sigma}})\coloneqq
f(\Psi(\beta,L^{\mathcal{I}_{\Sigma}}))$. If 2 is satisfied, then
$\mu_{\varepsilon}\searrow\mu$ as $\varepsilon\searrow 0$.

With requirements on the objective $f$, convergence of
$\varepsilon$-approximate solutions to the exact solution is guaranteed by the
following proposition (see Appendix E for proof).

######  Proposition 34.

For any continuous
$f:{}^{n_{\beta}}\times\mathbb{S}^{n_{\Sigma}}\rightarrow\real$, consider the
set-valued function $\hat{\theta}:\mathbb{R}_{\geq
0}\rightarrow\mathcal{P}(\Theta)$, defined as

$\hat{\theta}_{\varepsilon}\coloneqq\operatorname*{argmin}_{(\beta,\Sigma)\in\Psi(\Theta_{\varepsilon})}f(\beta,\Sigma)$

for all $\varepsilon>0$, and

$\hat{\theta}_{0}\coloneqq\operatorname*{argmin}_{(\beta,\Sigma)\in\Theta}f(\beta,\Sigma).$

If there exists $\alpha\in\real$ and compact $C\subseteq\Theta$ such that

$\Theta_{f\leq\alpha}\coloneqq\set{(\beta,\Sigma)\in\Theta}{f(\beta,\Sigma)\leq\alpha}$

is contained in $C$ and $\Theta_{f\leq\alpha}\cap\Theta_{++}$ is nonempty,
then there exists $\overline{\varepsilon}>0$ such that, for all
$\varepsilon_{0}\in[0,\overline{\varepsilon})$,

  1. (a)

$f$ achieves a minimum in $\Theta$ and $\hat{\theta}_{0}$ is nonempty;

  2. (b)

if $\varepsilon_{0}>0$, then $f$ achieves a minimum in
$\Psi(\Theta_{\varepsilon_{0}})$ and $\hat{\theta}_{\varepsilon_{0}}$ is
nonempty;

  3. (c)

$\mu_{\varepsilon}$ is continuous and $\hat{\theta}_{\varepsilon}$ is outer
semicontinuous at $\varepsilon=\varepsilon_{0}$; and

  4. (d)

if $\hat{\theta}_{0}$ is a singleton, then $\limsup_{\varepsilon\searrow
0}\hat{\theta}_{\varepsilon}=\hat{\theta}_{0}$.

##  V Case Studies

In this section, we apply the identification methods outlined in the previous
sections to design Kalman filters for the linear augmented disturbance models
used in offset-free MPC:

$\displaystyle\begin{bmatrix}\hat{x}_{k+1}\\\ \hat{d}_{k+1}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}A&B_{d}\\\
0&I_{p}\end{bmatrix}\begin{bmatrix}\hat{x}_{k}\\\
\hat{d}_{k}\end{bmatrix}+\begin{bmatrix}B\\\
0\end{bmatrix}u_{k}+\begin{bmatrix}K_{x}\\\ K_{d}\end{bmatrix}e_{k}$ (38a)
$\displaystyle y_{k}$
$\displaystyle=\begin{bmatrix}C&C_{d}\end{bmatrix}\begin{bmatrix}\hat{x}_{k}\\\
\hat{d}_{k}\end{bmatrix}+e_{k}$ (38b) $\displaystyle e_{k}$
$\displaystyle\overset{\textnormal{iid}}{\sim}\mathcal{N}(0,R_{e})$ (38c)

where $\hat{x}\in{}^{n}$ denote plant state estimates, $\hat{d}\in{}^{n_{d}}$
denote disturbance state estimates, $(B_{d},C_{d})\in{}^{n\times
n_{d}}\times{}^{p\times n_{d}}$ denote disturbance shaping matrices, and
$(K_{x},K_{d})\in{}^{n\times p}\times{}^{n_{d}\times p}$ are plant and
disturbance state filter gains.

Figure 1: Benchmark temperature Control Laboratory (TCLab) [61].

In the first case study, we consider the TCLab (Figure 1), an Arduino-based
temperature control laboratory that serves as a low-cost222The TCLab is
available for under $40 from https://apmonitor.com/heat.htm and
https://www.amazon.com/gp/product/B07GMFWMRY. benchmark for linear MIMO
control [61]. We identify the TCLab from open-loop data and use the resulting
model to design an offset-free MPC. We compare closed-loop control and
estimation performance of these models to that of offset-free MPCs designed
with the identification methods from [50, 51]. In the second case study, data
from an industrial-scale chemical reactor is used to design Kalman filters for
the linear augmented disturbance model, and the closed-loop estimation
performance is compared to that of the designs proposed in [51].

To aid in optimizer convergence and guarantee uniqueness of the solution, to
any likelihood function $f_{N}(\theta)$, we add a regularization term or prior
distribution,

$\min_{\theta\in\Theta}f_{N}(\theta)-\ln p(\theta)$

where $p(\theta)$ is the prior density. Throughout, we use the following
prior,

$-\ln
p(\beta,\Sigma)\propto\frac{\rho}{2}\left(\|\beta-\beta_{0}\|_{2}^{2}+\textnormal{tr}(\Sigma-\Sigma_{0})\right)$
(39)

where $\rho>0$ and $(\beta_{0},\Sigma_{0})\in\Theta_{++}$ is the initial guess
given to the optimizer. The prior eq. 39 is equivalent to a Frobenius norm
regularizer in the Cholesky factor space,

$-\ln
p(\Psi(\beta,L^{\mathcal{I}_{\Sigma}}))\propto\frac{\rho}{2}\left(\|\beta-\beta_{0}\|_{2}^{2}+\|L^{\mathcal{I}_{\Sigma}}-L^{\mathcal{I}_{\Sigma}}_{0}\|_{\mathrm{F}}^{2}\right.\\\
\left.+\|L^{\mathcal{J}_{\Sigma}}(\beta,L^{\mathcal{I}_{\Sigma}})-L^{\mathcal{J}_{\Sigma}}_{0}\|_{\mathrm{F}}^{2}\right)$
(40)

where
$L^{\mathcal{I}_{\Sigma}}_{0}\coloneqq\pi_{\mathcal{I}_{\Sigma}}^{L}[\textnormal{chol}(\Sigma_{0}-H(\beta_{0}))]$
and $L^{\mathcal{J}_{\Sigma}}_{0}\coloneqq
L^{\mathcal{J}_{\Sigma}}(\beta_{0},L^{\mathcal{I}_{\Sigma}}_{0})$. When
$\Sigma$ is block diagonal,
$L^{\mathcal{J}_{\Sigma}}(\beta,L^{\mathcal{I}_{\Sigma}})=0$ and the last term
of eq. 40 vanishes. The squared slack variables
$L_{\mathcal{A}}=L^{\mathcal{I}_{\mathcal{A}}}+L^{\mathcal{J}_{\mathcal{A}}}$
are not regularized.

Throughout these experiments, we use the steady-state filter likelihood eq. 10
and a Cholesky factor diagonal backoff $\varepsilon$ of $10^{-6}$. When
constraints are considered, they either take the form of smoothed spectral
radius and abscissa constraints,

$\displaystyle P_{d}-\frac{A_{K}P_{d}A_{K}^{\top}}{(1-\delta)^{2}}$
$\displaystyle=\varepsilon_{d}I_{n+p},$ $\displaystyle P_{d}$
$\displaystyle\succeq 0,$ $\displaystyle\textnormal{tr}(P_{d})$
$\displaystyle\leq\varepsilon_{d}^{-1}$ (41a) $\displaystyle
A_{K}P_{c}+P_{c}A_{K}^{\top}$ $\displaystyle=\varepsilon_{c}I_{n+p},$
$\displaystyle P_{c}$ $\displaystyle\succeq 0,$
$\displaystyle\textnormal{tr}(P_{c})$ $\displaystyle\leq\varepsilon_{c}^{-1}$
(41b)

or the LMI region barrier function constraints,

$\displaystyle P_{d}-\frac{A_{K}P_{d}A_{K}^{\top}}{(1-\delta)^{2}}$
$\displaystyle\succeq\varepsilon_{d}I,$ $\displaystyle P_{d}$
$\displaystyle\succeq 0,$ $\displaystyle\textnormal{tr}(P_{d})$
$\displaystyle\leq\varepsilon_{d}^{-1}$ (42a) $\displaystyle
A_{K}P_{c}+P_{c}A_{K}^{\top}$ $\displaystyle\succeq\varepsilon_{c}I,$
$\displaystyle P_{c}$ $\displaystyle\succeq 0,$
$\displaystyle\textnormal{tr}(P_{c})$ $\displaystyle\leq\varepsilon_{c}^{-1}.$
(42b)

Each optimization problem is formulated in CasADi as a
$\varepsilon$-approximate problem in Cholesky factor form eqs. 36 and 37, and
solved with IPOPT. Wall times for a single-thread run on an Intel Core
i9-10850K processor are reported. The initial guesses for these models are
based on a nested ML identification approach described in [50, 51]. This
approach effectively augments a standard identification method (e.g., PCA, Ho-
Kalman, canonical correlation analysis algorithms), so we refer to the initial
guess models as “augmented” versions of the standard method being used.

###  V-A TCLab

Figure 2: TCLab identification data and noise-free responses $\hat{y}_{k}=\sum_{j=1}^{k}\hat{C}\hat{A}^{j-1}\hat{B}u_{k-j}$ of a few selected models. TABLE I: TCLab model fitting results. ∗ The augmented PCA/CCA identification methods are not iterative. ∗∗ The maximum number of iterations was set at 500. Model | Time (s) | Iterations | Log-likelihood
---|---|---|---
Augmented PCA | 0.01 | N/A∗ | 3823.4
Augmented CCA | 0.04 | N/A∗ | 2415.8
Unregularized ML | 121.9∗∗ | 500∗∗ | -9431.4
Regularized ML | 7.5 | 17 | -9411.7
Constrained ML 1 | 36.3 | 61 | -9407.5
Constrained ML 2 | 45.1 | 78 | -9412.5
Constrained ML 3 | 14.2 | 23 | -9407.5
Constrained ML 4 | 25.1 | 54 | -9412.5
Figure 3: TCLab models open-loop and closed-loop (filter) eigenvalues.

Unless otherwise specified, the TCLab is modeled as a two-state system of the
form eq. 1, with internal temperatures as plant states
$x=\begin{bmatrix}T_{1}&T_{2}\end{bmatrix}^{\top}$, heater voltages as inputs
$u=\begin{bmatrix}V_{1}&V_{2}\end{bmatrix}^{\top}$, and measured temperatures
$y=\begin{bmatrix}T_{m,1}&T_{m_{2}}\end{bmatrix}^{\top}$ as outputs.
Throughout, we choose $n_{d}=p$ to satisfy the offset-free necessary
conditions in [37, 38], and we consider output disturbance models
$(B_{d},C_{d})=(0_{2\times 2},I_{2})$. We use an observability canonical form
[68] with $(A,B)$ fully parameterized and $C=I_{2}$ to guarantee model
identifiability and make the states interpretable as internal temperatures.
The remaining model terms $(K_{x},K_{d},R_{e})$ are fully parameterized. As in
eq. 38, the models do not include a passthrough term, i.e., $D=0$.

Eight TCLab models were considered:

  1. 1.

Augmented PCA: the 6-state TCLab model used in [50], where principle component
analysis on a $400\times 5100$ data Hankel matrix is used to determine the
states in the disturbance-free model.

  2. 2.

Augmented CCA: a 2-state augmented canonical correlation analysis (CCA) model,
based on the canonical method of [69].

  3. 3.

Unregularized ML: a 2-state model, fit directly to eq. 10 without
regularization or constraints, using Augmented CCA as the initial guess.

  4. 4.

Regularized ML: the same as Unregularized ML but with an added regularizer eq.
39 with $\rho=10^{-2}N$.

  5. 5.

Constrained ML 1: the same as Unregularized ML but with the smoothed LMI
region constraints eq. 42 with $\delta=\varepsilon_{d}=10^{-3}$ and
$\varepsilon_{c}=10$.

  6. 6.

Constrained ML 2: the same as Constrained ML 1 but with $\delta=2\times
10^{-3}$.

  7. 7.

Constrained ML 3: the same as Constrained ML 1, but using the smoothed
spectral radius and abscissa constraints eq. 41 (and the same constants).

  8. 8.

Constrained ML 4: the same as Constrained ML 3 (i.e., using eq. 41) but with
$\delta=2\times 10^{-3}$.

Each ML model uses Augmented CCA as the initial guess as it has the smallest
number of states. The augmented PCA model is, in effect, an unsupervised
learner of the state estimates, and therefore does not produce a parsimonious
state description.

In Figure 2, the identification data is presented along with the noise-free
responses $\hat{y}_{k}=\sum_{j=1}^{k}\hat{C}\hat{A}^{j-1}\hat{B}u_{k-j}$ of a
few selected models. Computation time, number of IPOPT iterations, and the
_unregularized_ log-likelihood value are reported in Table I. The open-loop
$A$ and closed-loop $A_{K}$ eigenvalues of each model are plotted in Figure 3.

Except for the augmented PCA model, all of the open-loop eigenvalues cluster
around the same region of the complex plane (fig. 3). The closed-loop filter
eigenvalues are placed similarly, with the exception of the unregularized ML
model, which has a slightly unstable filter. Despite these differences, the ML
models all have about the same unregularized log-likelihood value (Table I)
and appear to have identical noise-free responses (Figure 2). Here, the
constraints eqs. 41 and 42 appear to enforce filter stability and aid in
convergence with only a small computational penalty.

As reported in Table I, the unregularized ML model fitting did not converge.
It is our experience that encountering iterates with unstable filters and
filters with eigenvalues having negative real parts can cause convergence
issues. This is because the likelihood function becomes sensitive to small
changes in the parameter values due to filter instability or rapid
oscillations in the filter predictions.

Figure 4: TCLab setpoint tracking tests. Figure 5: TCLab disturbance rejection
tests.

Figure 6: TCLab setpoint tracking test performance.

Figure 7: TCLab disturbance rejection test performance.

To test offset-free control performance, we performed two sets of closed-loop
experiments on offset-free MPCs designed with the models. In Figure 4,
identical setpoint changes were applied to a TCLab running at a steady-state
power output of 50%. The setpoint changes were tracked with the offset-free
MPC design described in [50]. In Figure 5, step disturbances in the output
$p_{i}$ and the input $m_{i}$ are injected into a plant trying to maintain a
given steady-state temperature. The setpoints are tracked with the offset-free
MPC design described in [50].

Control performance indexed by squared distance from the setpoint
$\ell_{k}\coloneqq\|y_{k}-y_{\textnormal{sp},k}\|_{2}^{2}$. Estimation
performance is indexed by squared filter errors $e_{k}^{\top}e_{k}$. For any
signal $a_{k}$, we define a $T$-sample moving average by $\langle
a_{k}\rangle_{T}\coloneqq T^{-1}\sum_{j=0}^{T-1}a_{k-j}$. Setpoint tracking
performance is reported in Figure 6, and disturbance rejection performance is
reported in Figure 7. From Figures 6 and 7, it is clear that the regularized
ML model delivers the best overall performance.

Figure 8: TCLab identification index data for (left) setpoint tracking and
(right) disturbance rejection tests.

To investigate the accuracy of the stochastic models, we consider the inverse-
covariance weighted squared norm of the filter errors $q\coloneqq
e^{\top}R_{e}^{-1}e$ as an identification index. Recall the signal $e_{k}$ is
an i.i.d., zero-mean Gaussian process, i.e.,
$e_{k}\overset{\textnormal{iid}}{\sim}\mathcal{N}(0,R_{e})$, and therefore the
index $q_{k}$ is i.i.d. with a $\chi^{2}_{p}$ distribution. Moreover, the
moving average $\langle q_{k}\rangle_{T}$ is distributed as $\chi^{2}_{pT}/T$,
although it is no longer independent in time. In Figure 8, histograms of
$\langle q\rangle_{T}$ are plotted against their expected distribution for
$T\in\set{1,10,100}$ and the augmented PCA, augmented CCA, and regularized ML
models. The extreme discrepancies between the augmented PCA and CCA models’
performance index $\langle q\rangle_{T}$ and the reference distribution
$\chi^{2}_{pT}/T$ are primarily due to the augmented PCA/CCA models
significantly overestimating $R_{e}$,

$\displaystyle\hat{R}_{e}^{\textnormal{aug-PCA}}$
$\displaystyle=\begin{bmatrix}0.5871&0.3365\\\ 0.3365&0.2878\end{bmatrix},$
$\displaystyle\hat{R}_{e}^{\textnormal{aug-CCA}}$
$\displaystyle=\begin{bmatrix}0.5889&0.0918\\\ 0.1791&0.3152\end{bmatrix},$
$\displaystyle\hat{R}_{e}^{\textnormal{ML}}$
$\displaystyle=\begin{bmatrix}0.0107&0.0006\\\ 0.0006&0.008\end{bmatrix}.$

The reference distribution and the ML model’s $\langle q\rangle_{T}$
distribution diverge at large $T$ since, due to plant-model mismatch, the
filter’s innovation errors are slightly autocorrelated.

###  V-B Eastman reactor

$\overline{F}_{2}$
---
$\overline{T}_{H}$
---
$\overline{F}_{1}$
---
$F_{1,\textrm{sp}}$
---
$T_{\textrm{sp}}$
---
$F_{2,\textrm{sp}}$
---
MeOH flow rate, $F_{1}$
---
Reactor temp, $T$
---
Slurry flow rate, $F_{2}$
---
Hot oil TC.op, $T_{H}$
---
Quality param 1, $r_{1}$
---
Quality param 2, $r_{2}$
---
MPC
---
Hot oil return
---
Hot oil supply
---
Recycle
---
Sludge
---
Xylene
---
Crude DMT, H2O, MeOH
---
Slurry
---
MeOH
---
Figure 9: Schematic of the DMT reactor and MPC control strategy. TABLE II: Eastman reactor model fitting results. ∗ The augmented HK/CCA identification methods are not iterative. Model | Time (s) | Iterations | Log-likelihood
---|---|---|---
Augmented HK | 0.08 | N/A∗ | -7143.1
Augmented CCA | 0.07 | N/A∗ | -11288.4
Unregularized ML 1 | 23.7 | 93 | -14345.5
Unregularized ML 2 | 19.4 | 73 | -14345.5
Regularized ML 1 | 10.0 | 28 | -13053.6
Regularized ML 2 | 8.8 | 22 | -13587.7
Constrained ML 1 | 44.6 | 115 | -13034.1
Constrained ML 2 | 37.1 | 91 | -13587.7
Figure 10: Training data and noise-free responses for the Eastman reactor
models (Augmented HK and ML models using Augmented HK as the initial guess).
Figure 11: Training data and noise-free responses for the Eastman reactor
models (Augmented CCA and ML models using Augmented CCA as the initial guess).
Figure 12: Eastman reactor models open-loop and closed-loop (filter)
eigenvalues.

A schematic of the chemical reactor considered in the next case study is
presented in Figure 9.The control objective of the chemical reactor considered
in the next case study is to pick three inputs (the reactant flow rates and
utility temperatures $u=\begin{bmatrix}F_{1}&T_{H}&F_{2}\end{bmatrix}^{\top}$)
that steer the system to three setpoints (the output, a specified reactor
temperature $y=T$, and the flowrates
$\begin{bmatrix}u_{1}&u_{2}\end{bmatrix}^{\top}=\begin{bmatrix}F_{1}&F_{2}\end{bmatrix}^{\top}$)
without offset. See [51] for more details about the reactor operation. As in
Section V-A, we choose $n_{d}=p$, consider output disturbance models
$(B_{d},C_{d})=(0_{3\times 1},1)$, and use an observability canonical form
[68] this time with the parameterization $A=\begin{bmatrix}0&1&0\\\ 0&0&0\\\
a_{1}&a_{2}&a_{3}\end{bmatrix}$ and $C=\begin{bmatrix}1&0&0\end{bmatrix}$.
Again, the remaining model terms $(B,K_{x},K_{d},R_{e})$ are fully
parameterized. As in eq. 38, the models do not include a passthrough term,
i.e., $D=0$. We fit eight models (two augmented, six ML) of this form to
closed-loop data of this core 3-input, 1-output system:

  1. 1.

Augmented HK: a Ho-Kalman-based subspace model is augmented with a disturbance
model, as detailed in [50].

  2. 2.

Augmented CCA: a CCA subspace model, based on the method of [69], is augmented
with a disturbance model, as detailed in [51].

  3. 3,4)

Unregularized ML 1 and 2: a model is fit directly to eq. 10 without
regularization or constraints, using Augmented HK and CCA (resp.) as the
initial guesses.

  4. 5,6)

Regularized ML 1 and 2: the same as Unregularized ML 1 and 2, but with an
added regularizer eq. 39 with $\rho=4N$.

  5. 7,8)

Constrained ML 1 and 2: the same as Regularized ML 1 and 2, but with the
smoothed spectral radius and abscissa constraints eq. 41 with $\delta=0$,
$\varepsilon_{d}=10^{-3}$, and $\varepsilon_{c}=10$.

Computation time, number of IPOPT iterations, and the _unregularized_ log-
likelihood value are reported in Table I.In Figures 10 and 11, the
identification data is presented along with the noise-free responses
$\hat{y}_{k}=\sum_{j=1}^{k}\hat{C}\hat{A}^{j-1}\hat{B}u_{k-j}$ of the Ho-
Kalman-based and CCA-based models, respectively. The open-loop $A$ and closed-
loop $A_{K}$ eigenvalues of each model are plotted in Figure 12.

The eigenvalue clustering of the reactor models is more variable than that of
the TCLab models. Unregularized models converge to the same solution despite
different initial guesses, and regularization pulls the eigenvalues towards
the subspace model clustering around $z=0.98$. A notable feature of the
unregularized models is their strong oscillating modes at $\lambda=0.7246\pm
0.3141\iota$. While these oscillations are substantially dampened by
regularization and constraints, they persist for most models. These
oscillations are likely true plant dynamics caused by oscillations in the
lower-level PID loops (i.e., level or flow controls). This effect could be
seen as an advantage or disadvantage, with the model either correctly fitting
desired oscillations so they can be corrected in a supervisory MPC layer, or
the model incorrectly overfitting to unwanted process dynamics. Removing
unwanted process dynamics may be desirable when experimentation is costly, as
it allows the user to “fix” the data rather than run another experiment.

Figure 13: Test performance for the Eastman reactor models (Augmented HK and
ML models using Augmented HK as the initial guess).

Figure 14: Test performance for the Eastman reactor models (Augmented CCA and
ML models using Augmented CCA as the initial guess).

The models’ estimation performances are compared against each other in Figures
13 and 14 using one of the test datasets from in [51]. Control performance
could not be compared without a costly redesign of the existing MPC strategy.
The unregularized models perform the best here, providing over 70% reduction
in average filter error for the HK-based models. While the regularized models
still capture much of this improvement, they do not have the same level of
filter accuracy. This is again evidence of the oscillating modes being a real
part of the reactor dynamics rather than a feature of model overfitting.

##  VI Conclusion

We conclude with a discussion of possible future directions of research for
applying our ML identification scheme. An advantage of ML identification is
the large body of literature on the asymptotic distribution and statistical
efficiency of the parameter estimates. From the asymptotic distribution,
decision functions can be constructed to map the on-line MPC performance to a
re-identification signal. These decision functions can be constructed so as to
not alarm unless sufficiently exciting data is available. Decision-theoretic
re-identification therefore has lower cost and risk compared to classic
adaptive control or online reinforcement learning methods that require a
persistently exciting identification signal. This approach could bring
statistical data efficiency to the fields of adaptive control and online
reinforcement learning. Recent work on direct data-driven control has
incorporated likelihood functions with measurement noise models into the
control design [70]. To the best of our knowledge, no current work has
considered process noise, Kalman filter forms, or structuring the model with
uncontrollable integrators for offset-free MPC. There is a future possibility
of direct data-driven offset-free MPC design with both optimal control and
estimation performance.

##  Appendix A Proof of Proposition 6

Silverman [64] contains a more complete characterization of the DARE solutions
for regulation problems with cross terms. However, this admits additional
nullspace terms into the gain matrix which the Kalman filtering problem does
not allow. We avoid nullspace terms through the assumption $R_{v}\succ 0$ and
therefore streamline the proof of Proposition 6.

For the following definitions and lemmas, consider the system matrices
$\mathcal{W}\coloneqq(A,B,C,D)$ corresponding to a noise-free system.

######  Definition 35.

The system $\mathcal{W}$ is _left invertible on $\mathbb{I}_{0:k-1}$_ if

$0=\begin{bmatrix}D\\\ CB&D\\\ \vdots&\ddots&\ddots\\\
CA^{k-2}B&\ldots&CB&D\end{bmatrix}\begin{bmatrix}u_{0}\\\ \vdots\\\
u_{k-1}\end{bmatrix}$

implies $u_{0}=0$. The system $\mathcal{W}$ is _left invertible_ if there is
some $j\in\mathbb{N}$ such that $\mathcal{W}$ is left invertible on
$\mathbb{I}_{0:k-1}$ for all $k\geq j$.

######  Definition 36.

The system $\mathcal{W}$ is _strongly detectable_ if $y_{k}\rightarrow 0$
implies $x_{k}\rightarrow 0$.

The following lemmas are taken directly from [64, Thms. 8, 18(iii)], but the
proofs are omitted for the sake of brevity.

######  Lemma 37 ([64, Thm. 8]).

If $\mathcal{W}$ is left invertible, then $\mathcal{W}$ is strongly detectable
if and only if $(A-BF,C-DF)$ is detectable for all $F$ of appropriate
dimension.

######  Lemma 38 ([64, Thm. 18(iii)]).

If $\mathcal{W}$ is left invertible, then the DARE

$P=A^{\top}PA-(A^{\top}PB+C^{\top}D)(B^{\top}PB+D^{\top}D)^{-1}\\\
\times(B^{\top}PA+D^{\top}C)$

has a unique, stabilizing solution333Contrary to in Section II, here we mean
the solution $P$ is stabilizing when $A-BK(P)$ is stable, where
$K(P)\coloneqq(B^{\top}PB+D^{\top}D)^{-1}B^{\top}P$. if and only if
$\mathcal{W}$ is stabilizable and semistrongly detectable.

For the remainder of this section, we consider the full rank factorization

$\begin{bmatrix}Q_{w}&S_{wv}\\\
S_{wv}^{\top}&R_{v}\end{bmatrix}=\begin{bmatrix}\tilde{B}\\\
\tilde{D}\end{bmatrix}\begin{bmatrix}\tilde{B}^{\top}&\tilde{D}^{\top}\end{bmatrix}$

and the dual system
$\tilde{\mathcal{W}}\coloneqq(A^{\top},C^{\top},\tilde{B}^{\top},\tilde{D}^{\top})$
to analyze the properties of the original system eq. 1. The following lemma
relates the properties $R_{v}\succ 0$ and left invertability of
$\tilde{\mathcal{W}}$.

######  Lemma 39.

If $R_{v}\succ 0$ then $\tilde{\mathcal{W}}$ is left invertible.

###### Proof.

Left invertability on $\mathbb{I}_{0:k-1}$ is equivalent to

$0=\begin{bmatrix}\tilde{D}^{\top}\\\
\tilde{B}^{\top}C^{\top}&\tilde{D}^{\top}\\\ \vdots&\ddots&\ddots\\\
\tilde{B}^{\top}(A^{\top})^{k-2}C^{\top}&\ldots&\tilde{B}^{\top}C^{\top}&\tilde{D}^{\top}\end{bmatrix}\begin{bmatrix}u_{0}\\\
\vdots\\\ u_{k-1}\end{bmatrix}$ (43)

implying $u_{0}=0$. But $R_{v}=\tilde{D}\tilde{D}^{\top}\succ 0$, so
$\tilde{D}^{\top}$ has a zero nullspace. For each $k\in\mathbb{N}$, the
coefficient matrix of eq. 43 has a zero nullspace. Thus, $u_{0}=0$ and
$\tilde{\mathcal{W}}$ is left invertible. ∎

Finally, we can prove Proposition 6.

###### Proof of Proposition 6.

By Lemma 39, we have that $\tilde{\mathcal{W}}$ is left invertible. Therefore,
by Lemma 38, the DARE eq. 9 has a unique, stabilizing solution if and only if
$\tilde{\mathcal{W}}$ is stabilizable and strongly detectable. But by Lemma 37
and duality, the latter statement is true if and only if $(A,C)$ is detectable
and $(A-FC,\tilde{B}-F\tilde{D})$ is stabilizable for all $F\in{}^{n\times
p}$. ∎

##  Appendix B Proof of Proposition 20

Throughout this appendix, we define the set of $n\times n$ Hermitian,
Hermitian positive definite, and Hermitian positive semidefinite matrices as
$\mathbb{H}^{n}$, $\mathbb{H}^{n}_{++}$, and $\mathbb{H}^{n}_{+}$. Notice that
$f_{\mathcal{D}}$ maps to Hermitian matrices so we can write it as
$f:\mathbb{C}\rightarrow\mathbb{H}^{m}$. We define the extension of
$M_{\mathcal{D}}$ to complex arguments $M_{\mathcal{D}}:\mathbb{C}^{n\times
n}\times\mathbb{H}^{n}_{+}\rightarrow\mathbb{H}^{nm}$ as

$M_{\mathcal{D}}(A,P)\coloneqq M_{0}\otimes
P+M_{1}\otimes(AP)+M_{1}^{\top}\otimes(AP)^{\mathrm{H}}.$

To show Proposition 20, we need a preliminary result about Hermitian positive
semidefinite matrices, generalized from Lemma A.1 in [59].

######  Lemma 40.

For any $M\in\mathbb{H}^{n}$, if $M\succeq 0$ ($M\succ 0$) then
$\textnormal{Re}(M)\succeq 0$ ($\textnormal{Re}(M)\succ 0$).

###### Proof.

With $M=\textnormal{Re}(M)+\iota\textnormal{Im}(M)$, it is clear $M$ Hermitian
implies $\textnormal{Re}(M)$ is symmetric and $\textnormal{Im}(M)$ is skew-
symmetric. Thus $v^{\top}Mv=v^{\top}\textnormal{Re}(M)v$ for all $v\in{}^{n}$,
and positive (semi)definiteness of $M$ implies positive (semi)definiteness of
$\textnormal{Re}(M)$. ∎

In proving Proposition 20, we take the approach of [59] but are careful to
distinguish eigenvalues on the interior $\mathcal{D}$ from those on the
boundary $\partial\mathcal{D}$.

###### Proof of Proposition 20.

($\Leftarrow$) Suppose that $M_{\mathcal{D}}(A,P)\succeq 0$ for some $P\succ
0$ and let $\lambda\in\lambda(A)$. Then there exists a nonzero
$v\in\mathbb{C}^{n}$ for which $v^{\mathrm{H}}A=\lambda v^{\mathrm{H}}$.
Consider the identity

$\displaystyle(I_{m}$ $\displaystyle\otimes
v)^{\mathrm{H}}M_{\mathcal{D}}(A,P)(I_{m}\otimes v)$
$\displaystyle=M_{0}\otimes
v^{\mathrm{H}}Pv+M_{1}\otimes(v^{\mathrm{H}}APv)+M_{1}^{\top}\otimes(v^{\mathrm{H}}PA^{\top}v)$
$\displaystyle=M_{0}\otimes
v^{\mathrm{H}}Pv+M_{1}\otimes(\overline{\lambda}v^{\mathrm{H}}Pv)+M_{1}^{\top}\otimes(\lambda
v^{\mathrm{H}}Pv)$
$\displaystyle=v^{\mathrm{H}}Pv(M_{0}+M_{1}\lambda+M_{1}^{\top}\overline{\lambda})$
$\displaystyle=v^{\mathrm{H}}Pvf_{\mathcal{D}}(\lambda).$

The assumption $P\succ 0$ implies $v^{\mathrm{H}}Pv>0$, and
$M_{\mathcal{D}}(A,P)\succeq 0$ further implies
$f_{\mathcal{D}}(\lambda)\succeq 0$. Therefore
$\lambda\in\textnormal{cl}(\mathcal{D})$.

Next suppose $\lambda\in\lambda(A)$ is non-simple and
$\lambda\in\partial\mathcal{D}$. Then there exists nonzero
$v_{1},v_{2}\in\mathbb{C}^{n}$ (linearly independent) such that
$v^{\mathrm{H}}f_{\mathcal{D}}(\lambda)v=0$, $v_{1}^{\mathrm{H}}A=\lambda
v_{1}^{\mathrm{H}}$, and $v_{2}^{\mathrm{H}}A=\lambda
v_{2}^{\mathrm{H}}+v_{1}$. Because $\mathcal{D}$ is open,
$\lambda\in\partial\mathcal{D}=\textnormal{cl}(\mathcal{D})\setminus\mathcal{D}$
must satisfy both $f_{\mathcal{D}}(\lambda)\succeq 0$ and
$f_{\mathcal{D}}(\lambda)\not\succ 0$. Therefore $f_{\mathcal{D}}(\lambda)$ is
singular, and there exists a nonzero vector $v\in\mathbb{C}^{m}$ such that
$v^{\mathrm{H}}f_{\mathcal{D}}(\lambda)v=0$. With the $2\times 2$ matrices

$\displaystyle\tilde{P}$ $\displaystyle=\begin{bmatrix}p_{11}&p_{12}\\\
p_{12}&p_{22}\end{bmatrix}\coloneqq\begin{bmatrix}v_{1}^{\mathrm{H}}\\\
v_{2}^{\mathrm{H}}\end{bmatrix}P\begin{bmatrix}v_{1}&v_{2}\end{bmatrix}\succ
0$ $\displaystyle\tilde{J}$ $\displaystyle\coloneqq\lambda
I_{2}+\begin{bmatrix}0&1\\\ 0&0\end{bmatrix}$

we have
$\begin{bmatrix}v_{1}&v_{2}\end{bmatrix}^{\mathrm{H}}A=\tilde{J}\begin{bmatrix}v_{1}&v_{2}\end{bmatrix}^{\mathrm{H}}$
and therefore

$\displaystyle(I_{m}$
$\displaystyle\otimes\begin{bmatrix}v_{1}&v_{2}\end{bmatrix})^{\mathrm{H}}M_{\mathcal{D}}(A,P)(I_{m}\otimes\begin{bmatrix}v_{1}&v_{2}\end{bmatrix})$
$\displaystyle=M_{0}\otimes\tilde{P}+M_{1}\otimes\tilde{J}\tilde{P}+M_{1}^{\top}\otimes(\overline{\tilde{J}}\tilde{P})^{\top}$
$\displaystyle=M_{\mathcal{D}}(\tilde{J},\tilde{P})\succeq 0.$

Next, we have

$\displaystyle\tilde{M}$ $\displaystyle\coloneqq
K_{2,m}M_{\mathcal{D}}(\tilde{J},\tilde{P})K_{2,m}^{\top}$
$\displaystyle=\tilde{P}\otimes M_{0}+\tilde{J}\tilde{P}\otimes
M_{1}+(\overline{\tilde{J}}\tilde{P})^{\top}\otimes M_{1}^{\top}$
$\displaystyle=\tilde{P}\otimes
f_{\mathcal{D}}(\lambda)+\begin{bmatrix}p_{12}(M_{1}+M_{1}^{\top})&p_{22}M_{1}\\\
p_{22}M_{1}^{\top}&0\end{bmatrix}\succeq 0.$

Finally,

$(I_{2}\otimes v)^{\mathrm{H}}\tilde{M}(I_{2}\otimes v)\\\
=\begin{bmatrix}p_{12}v^{\mathrm{H}}(M_{1}+M_{1}^{\top})v&p_{22}v^{\mathrm{H}}M_{1}v\\\
p_{22}v^{\mathrm{H}}M_{1}^{\top}v&0\end{bmatrix}\succeq 0.$

But $\tilde{P}\succ 0$ implies $p_{22}>0$, so the above matrix inequality
implies $v^{\mathrm{H}}M_{1}v=0$. Moreover, with
$v^{\mathrm{H}}f_{\mathcal{D}}(\lambda)v=0$, we also have
$v^{\mathrm{H}}M_{0}v=0$ and therefore $f(z)\equiv 0$ and $\mathcal{D}$ is
empty, a contradiction. Therefore each $\lambda\in\lambda(A)$ non-simple
implies $\lambda\in\mathcal{D}$.

($\Rightarrow$) Suppose $\lambda(A)\subset\textnormal{cl}(\mathcal{D})$ and
$\lambda\in\lambda(A)$ non-simple implies $\lambda\in\mathcal{D}$.

If $A=\lambda$ is a (possibly complex) scalar, then it lies in
$\textnormal{cl}(\mathcal{D})$ by assumption, with
$M_{\mathcal{D}}(\lambda,p)=pf_{\mathcal{D}}(\lambda)\succeq 0$ for all $p>0$.

If $A=\lambda I_{n}+N$ is a (possibly complex) Jordan block, where
$N\in{}^{n\times n}$ is a shift matrix and $n>1$, then $\lambda\in\mathcal{D}$
and $f_{\mathcal{D}}(\lambda)\succ 0$. Let
$T_{k}\coloneqq\textnormal{diag}(k^{n-1},\ldots,k,1)$ for each
$k\in\mathbb{N}$. Then $T_{k}^{-1}AT_{k}=\lambda
I_{n}+k^{-1}N\rightarrow\lambda I_{n}$ as $k\rightarrow\infty$. Moreover,
because $M_{\mathcal{D}}$ is continuous, we have

$M_{\mathcal{D}}(T_{k}^{-1}AT_{k},I_{n})\rightarrow M_{\mathcal{D}}(\lambda
I_{n},I_{n})=f_{\mathcal{D}}(\lambda)\otimes I_{n}\succ 0.$

Therefore there exists some $k_{0}\in\mathbb{N}$ such that
$M_{\mathcal{D}}(T_{k}^{-1}AT_{k},I_{n})\succ 0$ for all $k\geq k_{0}$. With
$P\coloneqq T_{k}T_{k}^{\top}$, we have

$\displaystyle M_{\mathcal{D}}(A,P)$ $\displaystyle=M_{0}\otimes
T_{k}T_{k}^{\top}+M_{1}\otimes(AT_{k}T_{k}^{\top})$
$\displaystyle\qquad+M_{1}^{\top}\otimes(AT_{k}T_{k}^{\top})^{\top}$
$\displaystyle=(I_{m}\otimes T_{k})(M_{0}\otimes I_{n}+M_{1}\otimes
T_{k}^{-1}AT_{k}$
$\displaystyle\qquad+M_{1}^{\top}\otimes(T_{k}^{-1}AT_{k})^{\top})(I_{m}\otimes
T_{k})^{\top}$ $\displaystyle=(I_{m}\otimes
T_{k})M_{\mathcal{D}}(T_{k}^{-1}AT_{k},I_{n})(I_{m}\otimes T_{k})^{\top}\succ
0.$

Finally, for any $A\in{}^{n\times n}$, let
$A=V(\bigoplus_{i=1}^{p}J_{i})V^{-1}$ denote the Jordan decomposition of $A$,
where $J_{i}=\lambda_{i}I_{n_{i}}+N_{i}$, $\lambda_{i}\in\lambda(A)$, $N_{i}$
are shift matrices, and $n=\sum_{i=1}^{p}n_{i}$. We have already shown that
for each $i\in\mathbb{I}_{1:p}$, there exists $P_{i}\succ 0$ such that
$M_{\mathcal{D}}(J_{i},P_{i})\succeq 0$. Then with $\tilde{P}\coloneqq
V(\bigoplus_{i=1}^{p}P_{i})V^{-1}$, we have

$\displaystyle(I_{m}$ $\displaystyle\otimes
V^{-1})M_{\mathcal{D}}(A,\tilde{P})(I_{m}\otimes V^{-1})^{\mathrm{H}}$
$\displaystyle=M_{0}\otimes\left(\bigoplus_{i=1}^{p}P_{i}\right)+M_{1}\otimes\left(\bigoplus_{i=1}^{p}J_{i}P_{i}\right)$
$\displaystyle\qquad+M_{1}\otimes\left(\bigoplus_{i=1}^{p}J_{i}P_{i}\right)^{\top}$
$\displaystyle=K_{n,m}\left(\bigoplus_{i=1}^{p}K_{m,n_{i}}M_{\mathcal{D}}(J_{i},P_{i})K_{m,n_{i}}^{\top}\right)K_{n,m}^{\top}\succeq
0$

and therefore $M_{\mathcal{D}}(A,\tilde{P})\succeq 0$. Last, Lemma 40 gives
$M_{\mathcal{D}}(A,P)\succeq 0$ with $P\coloneqq\textrm{Re}(\tilde{P})$ since

$M_{\mathcal{D}}(A,P)=M_{\mathcal{D}}(A,\textrm{Re}(\tilde{P}))=\textrm{Re}(M_{\mathcal{D}}(A,\tilde{P})).\qed$

##  Appendix C Proof of Proposition 21

To show Proposition 21(a), we first require the following eigenvalue
sensitivity result due to [63, Thm. 7.2.3].

######  Theorem 41 ([63, Thm. 7.2.3]).

For any $A\in\mathbb{C}^{n\times n}$, denote its Schur decomposition by
$A=Q(D+N)Q^{\mathrm{H}}$, where $Q\in\mathbb{C}^{n\times n}$ is unitary,
$D\in\mathbb{C}^{n\times n}$ is diagonal, and $N\in\mathbb{C}^{n\times n}$ is
strictly upper triangular.444A matrix $U$ is strictly upper triangular if
$U_{ij}=0$ for all $i\geq j$. Let $p$ be the smallest integer for which
$M^{p}=0$ where $M_{ij}\coloneqq|N_{ij}|$. Then

$\min_{\lambda\in\lambda(A)}|\mu-\lambda|\leq\max\set{c}{E\|,(c\|E\|)^{1/p}}$

where $c\coloneqq\sum_{k=0}^{p-1}\|N\|^{k}$.

###### Proof of Proposition 21.

Throughout this proof, we show a set $S$ is not open (or not closed) by
demonstrating that $S^{c}$ (or $S$) does not contain all its limit points.

(a)—For any $A\in\mathbb{A}_{\mathcal{D}}^{n}$, continuity of
$f_{\mathcal{D}}$ gives the existence of a function $\delta(\lambda)>0$ such
that $f_{\mathcal{D}}(z)\succ 0$ for all $|z-\lambda|<\delta(\lambda)$ and
$\lambda\in\lambda(A)$. Let
$\delta\coloneqq\min_{\lambda\in\lambda(A)}\delta(\lambda)$. By Theorem 41 and
norm equivalence, there exist $c>0$ and $p\in\mathbb{I}_{1:n}$ such that

$\max_{\mu\in\lambda(A+E)}\min_{\lambda\in\lambda(A)}|\lambda-\mu|\leq\max\set{c}{E\|_{\textnormal{F}},(c\|E\|_{\textnormal{F}})^{1/p}}$

for all $E\in{}^{n\times n}$. Therefore there exists a $\varepsilon>0$ such
that

$\max_{\mu\in\lambda(A+E)}\min_{\lambda\in\lambda(A)}|\lambda-\mu|<\delta$

for all $E\in\mathcal{B}\coloneqq\set{E^{\prime}\in{}^{n\times
n}}{\|E^{\prime}\|_{\textnormal{F}}<\varepsilon}$. Finally, $A+\mathcal{B}$ is
a neighborhood of $A$ contained in $\mathbb{A}_{\mathcal{D}}^{n}$, and, since
$A\in\mathbb{A}_{\mathcal{D}}^{n}$ was chosen arbitrarily,
$\mathbb{A}_{\mathcal{D}}^{n}$ is open.

(b)(i)—Because $\mathcal{D}$ is open, nonempty, and not equal to
$\mathcal{D}$, $\partial\mathcal{D}$ is nonempty. Let
$\lambda\in\partial\mathcal{D}$ and $\lambda_{k}\in\mathcal{D}^{c}$ be a
sequence for which $\lambda_{k}\rightarrow\lambda$. By symmetry, we also have
$\overline{\lambda}\in\mathcal{D}$ and
$\overline{\lambda}_{k}\in\mathcal{D}^{c}$.

For $n=2$, we have
$A\coloneqq\begin{bmatrix}\textnormal{Re}(\lambda)&-\textnormal{Im}(\lambda)\\\
\textnormal{Im}(\lambda)&\textnormal{Re}(\lambda)\end{bmatrix}\in{}^{2\times
2}$ has eigenvalues $\lambda,\overline{\lambda}\in\mathcal{D}$, and
$A_{k}\coloneqq\begin{bmatrix}\textnormal{Re}(\lambda_{k})&-\textnormal{Im}(\lambda_{k})\\\
\textnormal{Im}(\lambda_{k})&\textnormal{Re}(\lambda_{k})\end{bmatrix}\in{}^{2\times
2}$ has eigenvalues $\lambda_{k},\overline{\lambda_{k}}\in\mathcal{D}^{c}$ for
each $k\in\mathbb{N}$. The corresponding eigenvectors are
$\begin{bmatrix}\pm\iota\\\ 1\end{bmatrix}\in\mathbb{C}^{2}$. Therefore
$A\in\tilde{\mathbb{A}}_{\mathcal{D}}^{2}$ but
$A_{k}\in(\tilde{\mathbb{A}}_{\mathcal{D}}^{2})^{c}$ for each
$k\in\mathbb{N}$, and the limit $A_{k}\rightarrow A$ gives us that
$(\tilde{\mathbb{A}}_{\mathcal{D}}^{2})^{c}$ does not contain all its limit
points.

For $n>2$, let $A_{0}\in\tilde{\mathbb{A}}_{\mathcal{D}}^{n-2}$, and we can
extend the prior argument with the sequence $B_{k}\coloneqq A_{k}\oplus
A_{0}\in(\tilde{\mathbb{A}}_{\mathcal{D}}^{n})^{c},k\in\mathbb{N}$ that
converges to $B\coloneqq A\oplus
A_{0}\in\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$.

(b)(ii)—By part (b)(i), it suffices to consider the case $n=1$. By closure and
convexity of $\mathcal{D}$, $\mathcal{D}\cap\real$ is either a closed line
segment, a closed ray, or itself. In other words, $\mathcal{D}\cap\real$ is
open if and only if it has no endpoints. Moreover, since
$\partial\mathcal{D}\cap\real$ is the set of the endpoints of
$\mathcal{D}\cap\real$, $\mathcal{D}\cap\real$ is open if and only if
$\partial\mathcal{D}\cap\real$ is empty. Finally, since
$\tilde{\mathbb{A}}_{\mathcal{D}}^{1}=\mathcal{D}\cap\real$,
$\tilde{\mathbb{A}}_{\mathcal{D}}^{1}$ is open if and only if
$\partial\mathcal{D}\cap\real$ is empty.

(c)(i)—Let $\lambda\in\partial\mathcal{D}$. Suppose $n=4$. Then
$\overline{\lambda}\in\partial\mathcal{D}$ by symmetry. Because $\mathcal{D}$
is open, there exists a sequence $\lambda_{k}\in\mathcal{D}$ such that
$\lambda_{k}\rightarrow\lambda$, and by symmetry, we also have
$\overline{\lambda}_{k}\in\mathcal{D}$ and
$\overline{\lambda}_{k}\rightarrow\overline{\lambda}$. Consider again the
$2\times 2$ matrices $A$ and $A_{k}$ from part (b)(i), which have eigenvalues
$\lambda,\overline{\lambda}\in\mathcal{D}$ and
$\lambda_{k},\overline{\lambda_{k}}\in\mathcal{D}^{c}$, respectively. Then the
block matrices $B\coloneqq\begin{bmatrix}A&I_{2}\\\
0&A\end{bmatrix}\in{}^{4\times 4}$ and
$B_{k}\coloneqq\begin{bmatrix}A_{k}&I_{2}\\\
0&A_{k}\end{bmatrix}\in{}^{4\times 4}$ have the same eigenvalues, but this
time the eigenvectors are $\begin{bmatrix}\pm\iota\\\ 1\\\ 0\\\
0\end{bmatrix},\begin{bmatrix}0\\\ 0\\\ \pm\iota\\\
1\end{bmatrix}\in\mathbb{C}^{4}$ and the eigenvalues are non-simple. Since
$\lambda$ is a non-simple eigenvalue on the boundary of $\mathcal{D}$, we have
$B\not\in\tilde{\mathbb{A}}_{\mathcal{D}}^{4}$. However, $\lambda_{k}$ are all
in the interior of $\mathcal{D}$, so
$B_{k}\in\tilde{\mathbb{A}}_{\mathcal{D}}^{4}$. Since $B_{k}\rightarrow B$,
the set $\tilde{\mathbb{A}}_{\mathcal{D}}^{4}$ does not contain all its limit
points.

On the other hand, let $\lambda\in\partial\mathcal{D}$ and suppose $n>4$.
Similarly to part (b)(i), with any
$\tilde{A}_{0}\in\tilde{\mathbb{A}}_{\mathcal{D}}^{n-4}$, we can extend the
argument for the $n=4$ case with the sequence $\tilde{A}_{k}\coloneqq
B_{k}\oplus\tilde{A}_{0}\in\tilde{\mathbb{A}}_{\mathcal{D}}^{n},k\in\mathbb{N}$
that converges to $\tilde{A}\coloneqq
B\oplus\tilde{A}_{0}\in(\tilde{\mathbb{A}}_{\mathcal{D}}^{n})^{c}$.

(c)(ii)—Let $\lambda\in\partial\mathcal{D}\cap\real$ and $n\geq 2$. Because
$\mathcal{D}$ is convex, open, and nonempty, there exists $\varepsilon>0$ such
that exactly one of the real intervals $(\lambda,\lambda+\varepsilon)$ or
$(\lambda-\varepsilon,\lambda)$ is contained in $\mathcal{D}$, whereas the
other is contained in $\textnormal{int}(\mathcal{D}^{c})$. Without loss of
generality, assume
$(\lambda-\varepsilon,\lambda)\subseteq\mathcal{D}$.555Otherwise, take the
reflection about the imaginary axis $-\mathcal{D}$ and
$-\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$. Then
$A_{k}\coloneqq(\lambda-\varepsilon/k)I_{n}+N_{n}\in\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$
for each $k\in\mathbb{N}$, but $A_{k}\rightarrow\lambda
I_{n}+N_{n}\in(\tilde{\mathbb{A}}_{\mathcal{D}}^{n})^{c}$ and therefore
$\tilde{\mathbb{A}}_{\mathcal{D}}^{n}$ does not contain all its limit points.
(d)—Since
$\overline{\mathbb{A}}_{\mathcal{D}}^{n}\coloneqq\set{A\in{}^{n\times
n}}{\lambda(A)\subset\textnormal{cl}(\mathcal{D})}$ contains
$\mathbb{A}_{\mathcal{D}}^{n}$, it suffices to show any
$A\in\overline{\mathbb{A}}_{\mathcal{D}}^{n}$ is a limit point of
$\mathbb{A}_{\mathcal{D}}^{n}$. Denote the Jordan form by
$A=V\left(\bigoplus_{i=1}^{p}\mu_{i}I_{n_{i}}+N_{n_{i}}\right)V^{-1}$, where
$V\in{}^{n\times n}$ is invertible, $\mu_{i}\in\lambda(A)$,
$n=\sum_{i=1}^{p}n_{i}$, and $N_{i}\in{}^{n_{i}\times n_{i}}$ is a shift
matrix. Because $\mu_{i}\in\textnormal{cl}(\mathcal{D})$, there exists a
sequence $\mu_{i,k}\in\mathcal{D}$ such that $\mu_{i,k}\rightarrow\mu_{i}$.
Then $A_{k}\coloneqq
V\left(\bigoplus_{i=1}^{p}\mu_{i,k}I_{n_{i}}+N_{i}\right)V^{-1}\in\mathbb{A}_{\mathcal{D}}^{n}$
and $A_{k}\rightarrow A$. ∎

##  Appendix D Proof of Proposition 22

To prove Proposition 22(a,b), we use sensitivity results on the value
functions of parameterized nonlinear SDPs,

$V(y)\coloneqq\inf_{x\in\mathbb{X}(y)}F(x,y)$ (44)

where the set-valued function
$\mathbb{X}:{}^{m}\rightarrow\mathcal{P}({}^{n})$ is defined by

$\mathbb{X}(y)\coloneqq\set{x\in{}^{n}}{G(x,y)\succeq 0}.$

Consider also the graph of the set-valued function $\mathbb{X}$,

$\mathbb{Z}\coloneqq\set{(x,y)\in{}^{n+m}}{G(x,y)\succeq 0}.$

Notice that $\mathbb{Z}$ is closed if $G$ is continuous. We say Slater’s
condition holds at $y\in{}^{m}$ if there exists $x\in{}^{n}$ such that
$x\in\textnormal{int}(\mathbb{X}(y))$, or equivalently, $G(x,y)\succ 0$.

In [71, Prop. 4.4], continuity of a general class of optimization problems is
considered. In the following proposition, we state the specialization to
nonlinear SDPs.

######  Proposition 42 ([71, Prop. 4.4]).

Let $y_{0}\in{}^{m}$ and suppose

  1. (i)

$F$ and $G$ are continuous on n+m;

  2. (ii)

there exist $\alpha\in\real$ and compact $C\subset{}^{n}$ such that, for each
$y$ in a neighborhood of $y_{0}$, the level set

$\textnormal{lev}_{\leq\alpha}F(\cdot,y)\coloneqq\set{x\in\mathbb{X}(y)}{F(x,y)\leq\alpha}$

is nonempty and contained in $C$; and

  3. (iii)

Slater’s condition holds at $y_{0}$.

Then $F(\cdot,y)$ attains a minimum on $\mathbb{X}(y)$ for all $y\in N_{y}$,
and $V(y)$ is continuous at $y=y_{0}$.

###### Proof.

See [71, Prop. 4.4] and the discussions in [71, pp. 264, 483–484, 491–492]. ∎

Finally, we prove Proposition 22.

###### Proof of Proposition 22.

Let $\textnormal{vec}:{}^{n\times n}\rightarrow{}^{n^{2}}$ and
$\textnormal{vecs}:{}^{n\times n}\rightarrow{}^{(1/2)(n+1)n}$ denote the
vectorization and symmetric vectorization operators, respectively.

(a)—With $x\coloneqq\textnormal{vecs}(P)$, $y\coloneqq\textnormal{vec}(A)$,
$F(x,y)\coloneqq\textnormal{tr}(VP)$, and $G(x,y)\coloneqq
P\oplus(M_{\mathcal{D}}(A,P)-M)$, we can use Proposition 42 to show the
continuity of $\phi_{\mathcal{D}}$ on $\mathbb{A}_{\mathcal{D}}^{n}$. Let
$A_{0}\in\mathbb{A}_{\mathcal{D}}^{n}$. Condition (i) of Proposition 42 holds
by assumption. Slater’s condition (iii) holds because for any $P\succ 0$ such
that $M_{\mathcal{D}}(A_{0},P)\succ 0$, we can define $P_{0}\coloneqq\gamma
P\succ 0$ for some
$\gamma>\gamma_{0}\coloneqq\|M\|\times\|[M_{\mathcal{D}}(A_{0},P)]^{-1}\|$ to
give

$M_{\mathcal{D}}(A_{0},P_{0})=\gamma
M_{\mathcal{D}}(A_{0},P)\succ\gamma_{0}M_{\mathcal{D}}(A_{0},P)\succeq M.$

Moreover, by continuity of $M_{\mathcal{D}}$, there exists a neighborhood
$N_{A}$ of $A_{0}$ such that $M_{\mathcal{D}}(A,P_{0})\succ M$ for all $A\in
N_{A}$. Letting $\alpha\coloneqq\textnormal{tr}(VP_{0})>0$, we have that the
set

$\set{P\in\mathbb{S}^{n}_{+}}{\textnormal{tr}(VP)\leq\alpha}$

is compact and contains the nonempty level set

$\set{P\in\mathbb{P}(A)}{\textnormal{tr}(VP)\leq\alpha}$

for all $A\in N_{A}$. Taking the image of each of the above sets under the
vecs operation gives condition (ii) of Proposition 42. All the conditions of
Proposition 42 are thus satisfied for each
$A_{0}\in\mathbb{A}_{\mathcal{D}}^{n}$, and we have $\phi_{\mathcal{D}}$ is
continuous on $\mathbb{A}_{\mathcal{D}}^{n}$.

(b)—Continuity of $\phi_{\mathcal{D}}$ on $\mathbb{A}_{\mathcal{D}}^{n}$
implies closure of the sublevel sets of $\phi_{\mathcal{D}}$, and eq. 24
follows by definition of $\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)$.

(c)—First, $M_{\mathcal{D}}(A,P)\succ 0$ implies $P\succ 0$ since, if
$M_{\mathcal{D}}(A,P)\succ 0$ and $P\succeq 0$ but $P\not\succ 0$, there
exists a nonzero $v\in{}^{n}$ such that $Pv=0$ and

$(I_{m}\otimes v)^{\top}(M_{\mathcal{D}}(A,P))(I_{m}\otimes
v)=M\otimes(v^{\top}Pv)\\\
+M_{1}\otimes(v^{\top}APv)+M_{1}^{\top}\otimes(v^{\top}PA^{\top}v)=0$

a contradiction of the assumption $M_{\mathcal{D}}(A,P)\succ 0$. Moreover, for
any $P\succ 0$ such that $M_{\mathcal{D}}(A,P)\succ 0$, we have
$M_{\mathcal{D}}(A,P)\succeq\gamma M_{\mathcal{D}}(A,P)\succeq M$ with
$P\coloneqq\gamma P$ and
$\gamma\coloneqq\|M\|\times\|[M_{\mathcal{D}}(A,P)]^{-1}\|$, so feasibility of
eq. 17 is equivalent to feasibility of

$\displaystyle M_{\mathcal{D}}(A,P)$ $\displaystyle\succ M,$ $\displaystyle P$
$\displaystyle\succeq 0$

and therefore
$\bigcup_{\varepsilon>0}\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)=\mathbb{A}_{\mathcal{D}}^{n}$.
But $\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)$ is monotonically
decreasing,666By “monotonically decreasing” we mean
$\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)\supseteq\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon^{\prime})$
for all $\varepsilon\leq\varepsilon^{\prime}$. so
$\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)\nearrow\bigcup_{\varepsilon>0}\mathbb{A}_{\mathcal{D}}^{n}(\varepsilon)=\mathbb{A}_{\mathcal{D}}^{n}$
as $\varepsilon\searrow 0$. ∎

##  Appendix E Proof of Propositions 33 and 34

Starting with Proposition 33:

###### Proof of Proposition 33.

Since $\mu_{\varepsilon}$ is nondecreasing and bounded from below by $\mu$, it
suffices to show that for each $\delta>0$, there exists a
$\overline{\varepsilon}>0$ such that
$\mu_{\overline{\varepsilon}}-\mu<\delta$.

Let $(\beta^{*},\Sigma^{*})\in\Theta$ denote a point for which
$\mu=f(\beta^{*},\Sigma^{*})$. If $(\beta^{*},\Sigma^{*})\in\Theta_{++}$, we
could simply choose $\overline{\varepsilon}>0$ small enough to put
$(\beta^{*},\Sigma^{*})$ in $\Theta_{\overline{\varepsilon}}$ and achieve
$\mu_{\overline{\varepsilon}}-\mu=0<\delta$.

Instead, we assume $(\beta^{*},\Sigma^{*})\not\in\Theta_{++}$. By 2, there
exists a sequence $(\beta_{k},\Sigma_{k})\in\Theta_{++},k\in\mathbb{N}$ such
that $\beta_{k}\rightarrow\beta$ and $\Sigma_{k}\rightarrow\Sigma$ as
$k\rightarrow\infty$. Defining $\nu_{k}\coloneqq f(\beta_{k},\Sigma_{k})$, we
have $\nu_{k}\rightarrow\mu$ by continuity of $f$. Therefore, there exists
some $k_{0}\in\mathbb{N}$ such that $\nu_{k}-\mu<\delta$ for all $k\geq
k_{0}$.

For each $(\beta_{k},\Sigma_{k})\in\Theta_{++}$, there exist unique
$L^{\mathcal{I}_{\Sigma}}_{k}\in\mathbb{L}^{n_{\Sigma}}_{++}[\mathcal{I}_{\Sigma}]$
and
$L^{\mathcal{I}_{\mathcal{A}}}_{k}\in\mathbb{L}^{n_{\mathcal{A}}}_{++}[\mathcal{I}_{\mathcal{A}}]$
such that the constraints

$\displaystyle\tilde{g}(\beta_{k},L^{\mathcal{I}_{\Sigma}}_{k},L^{\mathcal{I}_{\mathcal{A}}}_{k})$
$\displaystyle=0$ $\displaystyle h(\beta_{k},L^{\mathcal{I}_{\Sigma}}_{k})$
$\displaystyle\leq 0$

are satisfied (by Lemmas 29 and 31). Let $\overline{\varepsilon}$ be the
minimum value over all the diagonal elements of
$L^{\mathcal{I}_{\Sigma}}_{k_{0}}$ and
$L^{\mathcal{I}_{\mathcal{A}}}_{k_{0}}$. Then
$(\beta_{k_{0}},L^{\mathcal{I}_{\Sigma}}_{k_{0}})\in\Theta_{\overline{\varepsilon}}$
by construction, so $\nu_{k_{0}}\geq\mu_{\overline{\varepsilon}}$ by
optimality, and therefore
$\mu_{\overline{\varepsilon}}-\mu\leq\nu_{k_{0}}-\mu<\delta$. ∎

As in Appendix D, we use sensitivity results of [71] on optimization problems
to to prove Proposition 34. This time, however, we consider the continuity of
the value function for parameterized NLPs on Banach spaces. Let $\mathcal{X}$,
$\mathcal{Y}$, and $\mathcal{K}$ be Banach spaces and consider the
parameterized NLP,

$V(y)\coloneqq\inf_{x\in\mathbb{X}(y)}F(x,y)$ (45)

where the set-valued function
$\mathbb{X}:\mathcal{Y}\rightarrow\mathcal{P}(\mathcal{X})$ is defined by

$\mathbb{X}(y)\coloneqq\set{x\in\mathcal{X}}{G(x,y)\in K}$

for some $G:\mathcal{X}\times\mathcal{Y}\rightarrow\mathcal{K}$ and
$K\subseteq\mathcal{K}$ is closed. Let $X^{0}(y)$ denote the (possibly empty)
set of solutions to eq. 45. Define the graph of the set-valued function
$\mathbb{X}(\cdot)$ by

$\mathbb{Z}\coloneqq\set{(x,y)\in\mathcal{X}\times\mathcal{Y}}{G(x,y)\in K}.$

Notice that $\mathbb{Z}$ is closed if $G$ is continuous and $K$ is closed.

######  Proposition 43 ([71, Prop. 4.4]).

Let $y_{0}\in\mathcal{Y}$ and assume:

  1. (i)

$F$ and $G$ are continuous on $\mathcal{X}\times\mathcal{Y}$ and $K$ is
closed;

  2. (ii)

there exist $\alpha\in\real$ and a compact set $C\subseteq\mathcal{X}$ such
that, for every $y$ in a neighborhood of $y_{0}$, the level set

$\set{x\in\mathbb{X}(y)}{f(x,y)\leq\alpha}$

is nonempty and contained in $C$; and

  3. (iii)

for any neighborhood $N_{x}$ of the solution set $X^{0}(y_{0})$, there exists
a neighborhood $N_{y}$ of $y_{0}$ such that $N_{x}\cap\mathbb{X}(y)$ is
nonempty for all $y\in N_{y}$;

then $V(y)$ is continuous and $X^{0}(y)$ is outer semicontinuous at $y=y_{0}$.

###### Proof of Proposition 34.

First, we must specify $\overline{\varepsilon}$. For each
$\theta\in\Theta_{++}$, let

$\varepsilon(\theta)\coloneqq\max\set{\varepsilon>0}{\theta\in\Psi(\Theta_{\varepsilon})}$

where the maximum is achieved since there is a finite number of diagonal
elements of the Cholesky factors that must be lower bounded. Now we specify
$\overline{\varepsilon}$ as the supremum of $\varepsilon(\theta)$ over all
$\theta\in\Theta_{f\leq\alpha}\cap\Theta_{++}$,

$\overline{\varepsilon}\coloneqq\sup\set{\varepsilon(\theta)}{\theta\in\Theta_{f\leq\alpha}\cap\Theta_{++}}$

so that, for any $\varepsilon\in(0,\overline{\varepsilon})$,
$\Theta_{f\leq\alpha}\cap\Psi(\Theta_{\varepsilon})$ is nonempty and is
contained in the compact set $C$.

(a)—Following the proof of [71, Prop. 4.4], we have (i) $F$ is continuous and
(ii) the level set $\Theta_{f\leq\alpha}$ is nonempty and contained in the
compact set $C$, which implies $\Theta_{f\leq\alpha}$ is a compact level set
and therefore the minimum of $f$ over $\Theta_{f\leq\alpha}$ is achieved and
equals the minimum over $\Theta$. Moreover, $\hat{\theta}_{0}$ must be
nonempty.

(b)—Similarly to part (a), we have, for each
$\varepsilon\in(0,\overline{\varepsilon})$, that the level set
$\Theta_{f\leq\alpha}\cap\Psi(\Theta_{\varepsilon})$ is nonempty and contained
in the compact set $C$, so $f$ achieves its minimum over
$\Psi(\Theta_{\varepsilon})$ and $\hat{\theta}_{\varepsilon}$ is nonempty.

(c)—Consider the graph of the constraint function,

$\mathbb{Z}\coloneqq\set{(\theta,\varepsilon)\in\Theta\times\mathbb{R}_{\geq
0}}{\theta\in\Psi(\Theta_{\varepsilon})\textnormal{ if }\varepsilon>0}.$

Consider a sequence $(\theta_{k},\varepsilon_{k})\in\mathbb{Z},k\in\mathbb{N}$
that is convergent
$(\theta_{k},\varepsilon_{k})\rightarrow(\theta,\varepsilon)$. Then
$\varepsilon\geq 0$, otherwise the sequence would not converge. Moreover,
$\theta\in\Theta$ since
$\theta_{k}\in\Psi(\Theta_{\varepsilon_{k}})\subseteq\Theta$ for all
$k\in\mathbb{N}$ and $\Theta$ contains all its limit points. If
$\varepsilon=0$, then $(\theta,\varepsilon)\in\mathbb{Z}$ trivially. On the
other hand, if $\varepsilon>0$, then $\varepsilon(\theta_{k})$ converges to
$\varepsilon(\theta)$ because $\Psi$ is continuous and the max can be taken
over a finite number of elements of $\Psi^{-1}(\theta_{k})$. Moreover,
$\varepsilon(\theta_{k})$ and upper bounds $\varepsilon_{k}$ because
$\theta_{k}\in\Psi(\Theta_{\varepsilon_{k}})$, so
$\varepsilon(\theta)\geq\varepsilon$. Finally, we have
$\theta\in\Psi(\Theta_{\varepsilon})$, $(\theta,\varepsilon)\in\mathbb{Z}$,
and $\mathbb{Z}$ is closed.

Let $\varepsilon_{0}\geq 0$ and $N_{\theta}$ be a neighborhood of
$\hat{\theta}_{\varepsilon_{0}}$. With

$\delta\coloneqq\sup\set{\varepsilon(\theta)}{\theta\in N_{\theta}}>0$

we have $N_{\theta}\cap\Theta$ and $N_{\theta}\cap\Psi(\Theta_{\varepsilon})$
are nonempty for all $\varepsilon\in(0,\varepsilon_{0}+\delta)$.

Finally, the requirements of [71, Prop. 4.4] are satisfied for all
$\varepsilon_{0}\in[0,\overline{\varepsilon})$, so $\mu_{\varepsilon}$ is
continuous and $\hat{\theta}_{\varepsilon}$ is outer semicontinuous at
$\varepsilon=\varepsilon_{0}$.

(d)—The last statement follows by the definition of outer semicontinuity and
the fact that the $\limsup$ is nonempty. ∎

## References

  * [1] L. Ljung, _System Identification: Theory for the User_ , 2nd ed. New Jersey: Prentice Hall, 1999.
  * [2] E. J. Hannan and M. Deistler, _The Statistical Theory of Linear Systems_ , ser. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 2012.
  * [3] R. H. Shumway and D. S. Stoffer, _Time Series Analysis and its Applications_ , 4th ed. Cham, Switzerland: Springer International Publishing, 2017.
  * [4] K. Åström and P. Eykhoff, “System identification—A survey,” _Automatica_ , vol. 7, no. 2, pp. 123–162, 1971.
  * [5] L. Ljung, “Convergence analysis of parametric identification methods,” _IEEE Trans. Auto. Cont._ , vol. 23, no. 5, pp. 770–783, 1978.
  * [6] K. J. Åström, “Special section system identification tutorial: Maximum likelihood and prediction error methods,” _Automatica_ , vol. 16, no. 5, pp. 551–574, 1980.
  * [7] R. H. Shumway and D. S. Stoffer, “An approach to time series smoothing and forecasting using the EM algorithm,” _J. Time Series Anal._ , vol. 3, pp. 253–264, 1982.
  * [8] V. Digalakis, J. Rohlicek, and M. Ostendorf, “ML estimation of a stochastic linear system with the EM algorithm and its application to speech recognition,” _IEEE Trans. Speech and Audio Proc._ , vol. 1, no. 4, pp. 431–442, 1993.
  * [9] S. Gibson and B. Ninness, “Robust maximum-likelihood estimation of multivariable dynamic systems,” _Automatica_ , vol. 41, no. 10, pp. 1667–1682, 2005.
  * [10] M. Dewar and V. Kadirkamanathan, “A Canonical Space-Time State Space Model: State and Parameter Estimation,” _IEEE Trans. Signal Process._ , vol. 55, no. 10, pp. 4862–4870, 2007.
  * [11] C. G. Källström and K. J. Åström, “Experiences of System Identification Applied to Ship Steering Dynamics,” _IFAC Proceedings Volumes_ , vol. 12, no. 8, pp. 173–184, 1979.
  * [12] N. R. Kristensen, H. Madsen, and S. B. Jørgensen, “A method for systematic improvement of stochastic grey-box models,” _Comput. Chem. Eng._ , vol. 28, no. 8, pp. 1431–1449, 2004.
  * [13] E. Shahinfar, M. Bozorg, and M. Bidoky, “Parameter Estimation of an AUV Using the Maximum Likelihood method and a Kalman Filter with Fading Memory,” _IFAC Proceedings Volumes_ , vol. 43, no. 16, pp. 1–6, 2010.
  * [14] L. Simpson, A. Ghezzi, J. Asprion, and M. Diehl, “An Efficient Method for the Joint Estimation of System Parameters and Noise Covariances for Linear Time-Variant Systems,” in _2023 62nd IEEE Conference on Decision and Control (CDC)_ , Dec 2023, pp. 4524–4529.
  * [15] T. McKelvey, A. Helmersson, and T. Ribarits, “Data driven local coordinates for multivariable linear systems and their application to system identification,” _Automatica_ , vol. 40, no. 9, pp. 1629–1635, 2004.
  * [16] T. Ribarits, M. Deistler, and B. Hanzon, “An analysis of separable least squares data driven local coordinates for maximum likelihood estimation of linear systems,” _Automatica_ , vol. 41, no. 3, pp. 531–544, 2005.
  * [17] P. Li, I. Postlethwaite, and M. Turner, “Parameter Estimation Techniques for Helicopter Dynamic Modelling,” in _2007 American Control Conference_ , 2007, pp. 2938–2943.
  * [18] J. Umenberger, J. Wågberg, I. R. Manchester, and T. B. Schön, “Maximum likelihood identification of stable linear dynamical systems,” _Automatica_ , vol. 96, pp. 280–292, 2018.
  * [19] R. A. Redner and H. F. Walker, “Mixture Densities, Maximum Likelihood and the EM Algorithm,” _SIAM Rev._ , vol. 26, no. 2, pp. 195–239, 1984\.
  * [20] O. Bermond and J.-F. Cardoso, “Approximate likelihood for noisy mixtures,” in _International Workshop on Independent Component Analysis (ICA ’99)_ , vol. 99, Aussois, France, 1999, pp. 325–330.
  * [21] K. B. Petersen, O. Winther, and L. K. Hansen, “On the Slow Convergence of EM and VBEM in Low-Noise Linear Models,” _Neural Comput._ , vol. 17, no. 9, pp. 1921–1926, 2005.
  * [22] K. B. Petersen and O. Winther, “Explaining slow convergence of EM in low noise linear mixtures,” Informatics and Mathematical Modelling, Technical University of Denmark, Technical Report 2005-2, 2005. [Online]. Available: https://core.ac.uk/download/pdf/13755227.pdf
  * [23] R. K. Olsson, K. B. Petersen, and T. Lehn-Schiøler, “State-Space Models: From the EM Algorithm to a Gradient Approach,” _Neural Comput._ , vol. 19, no. 4, pp. 1097–1111, 2007.
  * [24] T. Westerlund, “A digital quality control system for an industrial dry process rotary cement kiln,” _IEEE Trans. Auto. Cont._ , vol. 26, no. 4, pp. 885–890, 1981.
  * [25] M. H. Caveness and J. J. Downs, “Reactor control using infinite horizon model predictive control,” AIChE Spring Meeting, Atlanta, GA, April 2005.
  * [26] H. Raghavan, A. K. Tangirala, B. R. Gopaluni, and S. L. Shah, “Identification of chemical processes with irregular output sampling,” _Control Eng. Pract._ , vol. 14, no. 5, pp. 467–480, 2006.
  * [27] L. W. Taylor, “On-Orbit Systems Identification of Flexible Spacecraft,” _IFAC Proceedings Volumes_ , vol. 18, no. 5, pp. 511–516, 1985.
  * [28] H. Melgaard, E. Hendricks, and H. Madsen, “Continuous Identification of a Four-Stroke SI Engine,” in _1990 American Control Conference_ , 1990, pp. 1876–1881.
  * [29] K. Åström and C. Källström, “Identification of ship steering dynamics,” _Automatica_ , vol. 12, no. 1, pp. 9–22, 1976.
  * [30] A. B. W. Jr., “On the numerical solution of initial/boundary-value problems in one space dimension,” _SIAM J. Numer. Anal._ , vol. 19, no. 4, pp. 683–697, 1982.
  * [31] H. White, “Maximum Likelihood Estimation of Misspecified Dynamic Models,” in _Misspecification Analysis_ , ser. Lecture Notes in Economics and Mathematical Systems, T. K. Dijkstra, Ed. Berlin, Heidelberg: Springer, 1984, pp. 1–19.
  * [32] B. D. O. Anderson, J. B. Moore, and R. M. Hawkes, “Model approximations via prediction error identification,” _Automatica_ , vol. 14, no. 6, pp. 615–622, 1978.
  * [33] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” _Control Eng. Pract._ , vol. 11, no. 7, pp. 733–764, 2003.
  * [34] J. B. Rawlings, D. Q. Mayne, and M. M. Diehl, _Model Predictive Control: Theory, Design, and Computation_ , 2nd ed. Santa Barbara, CA: Nob Hill Publishing, 2020, 770 pages, ISBN 978-0-9759377-5-4.
  * [35] W. M. Canney, “The future of advanced process control promises more benefits and sustained value,” _Oil & Gas J._, vol. 101, no. 16, pp. 48–54, Apr 2003.
  * [36] M. L. Darby and M. Nikolaou, “MPC: Current practice and challenges,” _Control Eng. Pract._ , vol. 20, no. 4, pp. 328 – 342, 2012.
  * [37] K. R. Muske and T. A. Badgwell, “Disturbance modeling for offset-free linear model predictive control,” _J. Proc. Cont._ , vol. 12, no. 5, pp. 617–632, 2002.
  * [38] G. Pannocchia and J. B. Rawlings, “Disturbance models for offset-free MPC control,” _AIChE J._ , vol. 49, no. 2, pp. 426–437, 2003.
  * [39] M. Wallace, B. Das, P. Mhaskar, J. House, and T. Salsbury, “Offset-free model predictive control of a vapor compression cycle,” _J. Proc. Cont._ , vol. 22, no. 7, pp. 1374–1386, 2012.
  * [40] M. Wallace, P. Mhaskar, J. House, and T. I. Salsbury, “Offset-Free Model Predictive Control of a Heat Pump,” _Ind. Eng. Chem. Res._ , vol. 54, no. 3, pp. 994–1005, 2015.
  * [41] P. Schmid and P. Eberhard, “Offset-free Nonlinear Model Predictive Control by the Example of Maglev Vehicles,” _IFAC–P. Online_ , vol. 54, no. 6, pp. 83–90, 2021.
  * [42] F. Xu, Y. Shi, K. Zhang, and X. Xu, “Real-Time Application of Robust Offset-Free MPC in Maglev Planar Machine,” _IEEE Trans. Indus. Elec._ , vol. 70, no. 6, pp. 6121–6130, 2023.
  * [43] R. Huang, S. C. Patwardhan, and L. T. Biegler, “Offset-free Nonlinear Model Predictive Control Based on Moving Horizon Estimation for an Air Separation Unit,” _IFAC Proceedings Volumes_ , vol. 43, no. 5, pp. 631–636, 2010.
  * [44] L. N. Petersen, N. K. Poulsen, H. H. Niemann, C. Utzen, and J. B. Jørgensen, “Comparison of three control strategies for optimization of spray dryer operation,” _J. Proc. Cont._ , vol. 57, 2017.
  * [45] S.-K. Kim, D.-K. Choi, K.-B. Lee, and Y. I. Lee, “Offset-Free Model Predictive Control for the Power Control of Three-Phase AC/DC Converters,” _IEEE Trans. Indus. Elec._ , vol. 62, no. 11, pp. 7114–7126, 2015.
  * [46] F. A. Bender, S. Goltz, T. Braunl, and O. Sawodny, “Modeling and Offset-Free Model Predictive Control of a Hydraulic Mini Excavator,” _IEEE Trans. Auto. Sci. Eng._ , vol. 14, no. 4, pp. 1682–1694, 2017.
  * [47] D. Deenen, E. Maljaars, L. Sebeke, B. De Jager, E. Heijman, H. Grüll, and W. P. M. H. Heemels, “Offset-free model predictive control for enhancing MR-HIFU hyperthermia in cancer treatment,” _IFAC–P. Online_ , vol. 51, no. 20, pp. 191–196, 2018.
  * [48] B. J. Odelson, A. Lutz, and J. B. Rawlings, “Application of autocovariance least-squares methods to laboratory data,” TWMCC, Department of Chemical Engineering, University of Wisconsin-Madison, Tech. Rep. 2003–03, September 2003\. [Online]. Available: https://engineering.ucsb.edu/~jbraw/jbrweb-archives/tech-reports/twmcc-2003-03.pdf
  * [49] M. A. Zagrobelny and J. B. Rawlings, “Identifying the uncertainty structure using maximum likelihood estimation,” in _American Control Conference_ , Chicago, IL, July 1–3, 2015, pp. 422–427.
  * [50] S. J. Kuntz and J. B. Rawlings, “Maximum likelihood estimation of linear disturbance models for offset-free model predictive control,” in _American Control Conference_ , Atlanta, GA, June 8–10, 2022, pp. 3961–3966.
  * [51] S. J. Kuntz, J. J. Downs, S. M. Miller, and J. B. Rawlings, “An industrial case study on the combined identification and offset-free control of a chemical process,” _Comput. Chem. Eng._ , vol. 179, 2023.
  * [52] F. Dorfler, J. Coulson, and I. Markovsky, “Bridging Direct & Indirect Data-Driven Control Formulations via Regularizations and Relaxations,” _IEEE Trans. Auto. Cont._ , vol. 68, no. 2, pp. 883–897, 2022.
  * [53] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Linear Tracking MPC for Nonlinear Systems—Part II: The Data-Driven Case,” _IEEE Trans. Auto. Cont._ , vol. 67, no. 9, pp. 4406–4421, 2022.
  * [54] Z. Yuan and J. Cortés, “Data-Driven Optimal Control of Bilinear Systems,” _IEEE Ctl. Sys. Let._ , vol. 6, pp. 2479–2484, 2022.
  * [55] G. Bianchin, M. Vaquero, J. Cortés, and E. Dall’Anese, “Online Stochastic Optimization for Unknown Linear Systems: Data-Driven Controller Synthesis and Analysis,” _IEEE Trans. Auto. Cont._ , pp. 1–15, 2023.
  * [56] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. M. De Moor, “A note on persistency of excitation,” _Sys. Cont. Let._ , vol. 54, no. 4, pp. 325–329, 2005.
  * [57] M. Diehl, K. Mombaur, and D. Noll, “Stability Optimization of Hybrid Periodic Systems via a Smooth Criterion,” _IEEE Trans. Auto. Cont._ , vol. 54, no. 8, pp. 1875–1880, 2009.
  * [58] J. Vanbiervliet, B. Vandereycken, W. Michiels, S. Vandewalle, and M. Diehl, “The Smoothed Spectral Abscissa for Robust Stability Optimization,” _SIAM J. Optim._ , vol. 20, no. 1, pp. 156–171, 2009\.
  * [59] M. Chilali and P. Gahinet, “$H_{\infty}$ design with pole placement constraints: an LMI approach,” _IEEE Trans. Auto. Cont._ , vol. 41, no. 3, pp. 358–367, 1996.
  * [60] S. Burer, R. D. Monteiro, and Y. Zhang, “Solving a class of semidefinite programs via nonlinear programming,” _Math. Prog._ , vol. 93, no. 1, pp. 97–122, Jun 2002.
  * [61] J. Park, R. A. Martin, J. D. Kelly, and J. D. Hedengren, “Benchmark temperature microcontroller for process dynamics and control,” _Comput. Chem. Eng._ , vol. 135, Apr 2020.
  * [62] D. N. Miller and R. A. De Callafon, “Subspace identification with eigenvalue constraints,” _Automatica_ , vol. 49, no. 8, pp. 2468–2473, 2013.
  * [63] G. H. Golub and C. F. Van Loan, _Matrix Computations_ , 4th ed. Baltimore, MD: The Johns Hopkins University Press, 2013.