# Optimal Sensor Precision for Multi-Rate Sensing for Bounded Estimation Error

Niladri Das and Raktim Bhattacharya Niladri Das (graduate student) and Raktim
Bhattacharya (associate professor) are with the Department of Aerospace
Engineering, Texas A& M University, College Station, TX, 77845 USA. E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We address the problem of determining optimal sensor precisions for estimating
the states of linear time-varying discrete-time stochastic dynamical systems,
with guaranteed bounds on the estimation errors. This is performed in the
Kalman filtering framework, where the sensor precisions are treated as
variables. They are determined by solving a constrained convex optimization
problem, which guarantees the specified upper bound on the posterior error
variance. Optimal sensor precisions are determined by minimizing the $l_{1}$
norm, which promotes sparseness in the solution and indirectly addresses the
sensor selection problem. The theory is applied to realistic flight mechanics
and astrodynamics problems to highlight its engineering value. These examples
demonstrate the application of the presented theory to a) determine redundant
sensing architectures for linear time invariant systems, b) accurately
estimate states with low-cost sensors, and c) optimally schedule sensors for
linear time-varying systems.

###### Index Terms:

Sensor precision, Kalman filtering, Convex optimization.

##  I Introduction

Kalman filtering provides the state estimate with minimum error covariance for
a given sensor precision. Sensor precision is defined as the inverse of sensor
noise. In this paper, we look at the inverse problem. Given a dictionary of
sensors, what is the most noise in these sensors such that the error
covariance is below a given error bound. This problem is of significant
engineering value since it allows system designers to determine the optimal
accuracy of sensing components that satisfy a given system-level error budget.
Since cost of sensors are proportional to the accuracy, the problem has
economical implications as well. In this paper we present convex optimization
formulations to solve the optimal sensing precision problem in a Kalman
filtering framework, for linear discrete-time time-varying systems with multi-
rate sensing.

The proposed framework also impacts other important problems in sensing, such
as sensor scheduling and selection. Existing algorithms for sensor scheduling
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10] and selection [11, 12, 13, 14, 15] assume that
the sensor’s noise variance is given. The framework in this paper can be used
to determine them optimally.

While this paper focuses on determining the optimal sensor precisions, the
presented framework can also be applied to schedule and select sensors.
Starting with a dictionary of sensors with unknown precisions, minimizing the
$l_{1}$ norm of the sensor precisions will promote sparseness in the solution.
Sensors with zero precisions can be eliminated from the system, and sensors
with non zero precisions will guarantee the required estimation accuracy.
Therefore, it is possible to simultaneously determine the optimal sensor
precisions and prune out unnecessary sensors using the proposed framework.
This is appealing because the problem can be solved in a convex optimization
framework in a general setting. This is in contrast with the NP-hard
formulations and heuristic methods to solve them.

###  I-A Key Contribution

The problem of determining optimal sensor precision was first introduced in
[16], in the context of output-feedback controller design for continuous-time
systems without uncertainty and steady-state performance guarantees. Recently,
we extended that work to determine optimal sensing precisions for continuous-
time robust output-feedback control, with guaranteed $\mathcal{H}_{2}$
performance [17].

In this paper, we look at the problem of determining the optimal sensor
precision for state estimation of linear time-varying, discrete-time,
stochastic systems with multi-rate sensing. To the best of our knowledge, this
is the first paper that determines optimal sensor precision for these systems.
The main results are presented as two theorems, deriving the convex
optimization formulations, which achieve the objective mentioned above. The
first theorem addresses the problem of determining optimal sensor precisions
to bound the estimation error over one time step. The second theorem
determines the optimal sensor precision to bound the steady-state error for
periodic time-varying systems.

###  I-B Layout

Section II presents the preliminaries for this paper, where we define the
model of the dynamical system and the sensor model assumed in the paper. We
also define the problem to determine the optimal precisions for a given upper
bound in the estimation error. Section III derives the optimization problem
for determining the sensing precision if the objective is to bound the error
after a one-time step. Section IV, derives the optimization problem for
determining the sensor precisions to bound the steady-state estimation error
for periodic systems. Section V highlights the engineering relevance of the
proposed framework, where the theoretical results are applied to state-
estimation problems from flight mechanics and astrodynamics. In that section,
we present three examples that demonstrate how the proposed theory can be
applied to identify redundancy in sensing, accurately estimate states with
cheap sensors, and optimally schedule sensors for time-varying periodic
systems.

##  II Preliminaries

We focus on determining the optimal sensor precision for linear time-varying
discrete-time stochastic systems described by the model of the form:

$\displaystyle\boldsymbol{x}_{k+1}$
$\displaystyle=\boldsymbol{A}_{k}\boldsymbol{x}_{k}+\boldsymbol{B}_{k}\boldsymbol{w}_{k},$
(1a) $\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=\boldsymbol{C}_{k}\boldsymbol{x}_{k}+\boldsymbol{n}_{k},$ (1b)

where $k=0,1,2,...$ are the time indices, $\boldsymbol{A}_{k}\in^{n_{x}\times
n_{x}},\boldsymbol{B}_{k}\in^{n_{x}\times n_{w}}$ and
$\boldsymbol{C}_{k}\in^{n_{y_{k}}\times\boldsymbol{n}_{x}}$ are the system
matrices, $\boldsymbol{x}_{k}\in^{n_{x}}$ is the $n_{x}$ dimensional state of
the model at time instant $k$, $\boldsymbol{w}_{k}\in^{n_{w}}$ is the $n_{w}$
dimensional zero-mean Gaussian additive process noise variable with
$\mathbb{E}[\boldsymbol{w}_{k}\boldsymbol{w}_{l}^{T}]=\boldsymbol{Q}_{k}$,
where $\mathbb{E}[.]$ denotes the expected value. The $n_{y_{k}}$ dimensional
observations at time $k$ is denoted by $\boldsymbol{y}_{k}\in^{n_{y_{k}}}$,
which is corrupted by an $n_{y_{k}}$ dimensional additive observation noise
$\boldsymbol{n}_{k}\in^{n_{y_{k}}}$ at time instant $k$. The sensor noise at
each time instant is a zero mean Gaussian random variable with
$\mathbb{E}[\boldsymbol{n}_{k}\boldsymbol{n}_{l}^{T}]=\boldsymbol{R}_{k}\delta_{kl}$.
The initial conditions are
$\mathbb{E}\left[\boldsymbol{x}_{0}\right]=\boldsymbol{\mu}_{0}$ and
$\mathbb{E}\left[(\boldsymbol{x}_{0}-\boldsymbol{\mu}_{0}){(\boldsymbol{x}_{0}-\boldsymbol{\mu}_{0})}^{T}\right]=\boldsymbol{\Sigma}_{0}$.
The process noise $\boldsymbol{w}_{k}$, observation noise
$\boldsymbol{n}_{k}$, and initial state variable $\boldsymbol{x}_{0}$ are
assumed to be independent.

In Kalman filtering, the propagation equations for the mean and the error
covariance are given by [18]

$\displaystyle\boldsymbol{\mu}^{-}_{k}=\boldsymbol{A}_{k}\boldsymbol{\mu}^{+}_{k-1}$
$\displaystyle\text{ and
}\boldsymbol{\Sigma}_{k}^{-}=\boldsymbol{A}_{k}\boldsymbol{\Sigma}^{+}_{k-1}\boldsymbol{A}_{k}^{T}+\boldsymbol{B}_{k}\boldsymbol{Q}_{k}\boldsymbol{B}_{k}^{T},$

and update equation for mean is
$\boldsymbol{\mu}^{+}_{k}=\boldsymbol{\mu}^{-}_{k}+\boldsymbol{K}_{k}(\boldsymbol{y}_{k}-\boldsymbol{C}_{k}\boldsymbol{\mu}^{-}_{k})$,
and for variance is

$\boldsymbol{\Sigma}^{+}_{k}=(\boldsymbol{I}_{n_{x}}-\boldsymbol{K}_{k}\boldsymbol{C}_{k})\boldsymbol{\Sigma}_{k}^{-}(\boldsymbol{I}_{n_{x}}-\boldsymbol{K}_{k}\boldsymbol{C}_{k})^{T}+\boldsymbol{K}_{k}\boldsymbol{R}_{k}\boldsymbol{K}_{k}^{T},$
(2)

where $\boldsymbol{K}_{k}$ is the Kalman gain determined by minimizing
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}_{k}\right)$ for a given
$\boldsymbol{R}_{k}$.

In this paper, we treat $\boldsymbol{R}_{k}$ as a variable, and for a given
error bound $\gamma_{d}\in^{+}$ (positive real number) we determine the
maximum $\boldsymbol{R}_{k}$ such that
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}_{k}\right)\leq\gamma_{d}$. Assuming
$\boldsymbol{R}_{k}$ to be diagonal, i.e.
$\boldsymbol{R}_{k}:=\mathbf{diag}\left(\boldsymbol{r}_{k}\right)$, where
$\boldsymbol{r}_{k}:=\begin{bmatrix}r_{1}&r_{2}&\cdots&r_{n_{y_{k}}}\end{bmatrix}^{T}$,
with $r_{i}>0$. As shown later, it is convenient to formulate the problem in
terms of sensor precisions, defined by $\boldsymbol{S}_{k}$, which is the
inverse of sensor noise $\boldsymbol{R}_{k}$, i.e.
$\boldsymbol{S}_{k}:=\boldsymbol{R}_{k}^{-1}$, resulting in $s_{i}:=1/r_{i}$.
Defining
$\boldsymbol{s}:=\begin{bmatrix}s_{1}&s_{2}&\cdots&s_{n_{y}}\end{bmatrix}^{T}$,
we can determine the least precise sensors, i.e. minimize
$\mathbf{tr}\left(\boldsymbol{S}_{k}\right)$, for which
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}_{k}\right)\leq\gamma_{d}$ is
guaranteed. Since more precise sensors are more expensive, satisfying required
accuracy with least precise sensors has favorable economic implications.

Minimization of $\mathbf{tr}\left(\boldsymbol{S}_{k}\right)$ also has
sparseness implications as $\mathbf{tr}\left(\boldsymbol{S}_{k}\right)$, for
$\boldsymbol{S}_{k}\geq 0$, is equivalent to the $\|\boldsymbol{s}\|_{1}$.
Since it is well-known that the $l_{1}$ norm is sparseness promoting [19],
minimizing $\mathbf{tr}\left(\boldsymbol{S}_{k}\right)$ will result in a
sparse solution that satisfies
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}_{k}\right)\leq\gamma_{d}$, if a
sparse solution exists for the problem. Consequently, sensors with zero
precisions would not contribute to achieving
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}_{k}\right)\leq\gamma_{d}$, and thus
can be removed from the system. It should also me mentioned that for a given
system parameters, there exists a range of $\gamma_{d}$, for which the
minimization problem is infeasible.

With this background, we next present the formulation to determine the optimal
precision for multi-rate information fusion in the Kalman filtering framework.

##  III Optimal Sensing Precision for Update after One Time Step

Here we present the formulation that determines the optimal sensing precision
that guarantees a bounded estimation error after one update. The problem is
formulated in a batch processing framework, where $m$ measurements are
collected before the state is updated. Specifically, given the state
uncertainty at time $t_{km}$, the objective is to determine the precisions of
these $m$ measurements, such that the state estimate at time $t_{(k+1)m}$ has
bounded error.

Consider sensing over $m$ time steps, as shown below.

$t_{km}$$\boldsymbol{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}x}}_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}km}}$$t_{km+1}$$\boldsymbol{y}_{km+1}$$t_{km+2}$$\boldsymbol{y}_{km+2}$$t_{(k+1)m}$$\begin{matrix}\boldsymbol{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}x}}_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}(k+1)m}}\\\\[5.69054pt]
\boldsymbol{y}_{(k+1)m}\end{matrix}$ Figure 1: Multi-rate measurements over
$m$ time steps.

This allows us to model multi-rate sensing, with $m$ being the least-common-
multiple of the various sensing intervals. With each measurement
$\boldsymbol{y}_{km+j}$, for $j=1,\cdots,m$, we associate sensor noises
$\boldsymbol{r}_{km+j}\in^{n_{y_{km+j}}}$.

In conventional Kalman filtering, the sensor noises $\boldsymbol{r}_{km+j}$
are known and the objective is to estimate the state at time $t_{(k+1)m}$,
given the posterior at time $t_{km}$ and measurements
$\boldsymbol{y}_{km+1},\cdots,\boldsymbol{y}_{(k+1)m}$, at times
$t_{km+1},\cdots,t_{(k+1)m}$. This scenario is common in control system
applications where the control-loop is band limited (to prevent excitation of
high-frequency dynamics), but the sensing loop can be faster. In such a
scenario, the Kalman filter determines the state estimate every $m$ time
steps, by batch-processing $m$ measurements. In this paper, we are interested
in maximizing $\boldsymbol{r}_{km+j}$ for which the posterior state estimation
error at time $t_{(k+1)m}$ satisfies a given upper bound.

###  III-A Augmented Dynamical System

To determine the posterior at time $t_{(k+1)m}$, we need to propagate the
state uncertainty from $t_{km}$ to $t_{(k+1)m}$, to obtain the prior at
$t_{(k+1)m}$. This is done by lifting discrete-time signals defined over times
$t_{k}$ to signals defined over times $t_{km}$, for $k=0,1,\cdots,\infty$. We
define lifted signals $\boldsymbol{X}_{k}\in^{mn_{x}}$ and
$\boldsymbol{Y}_{k}\in^{n_{y_{k,m}}}$ as a vector with $m$ consecutive state
vectors and measurements stacked vertically, respectively, i.e.

$\displaystyle\boldsymbol{X}_{k}:=\begin{pmatrix}\boldsymbol{x}_{km+1}\\\
\vdots\\\ \boldsymbol{x}_{(k+1)m}\end{pmatrix},$
$\displaystyle\boldsymbol{Y}_{k}:=\begin{pmatrix}\boldsymbol{y}_{km+1}\\\
\vdots\\\ \boldsymbol{y}_{(k+1)m}\end{pmatrix},$ (3)

where $m\geq 1$, and $n_{y_{k,m}}=\sum_{j=1}^{m}n_{y_{km+j}}$.

Using the system defined in equation (1a), the state $\boldsymbol{X}_{k}$ can
be expressed as [20]

$\boldsymbol{X}_{k}=\boldsymbol{\mathcal{A}}_{k}\boldsymbol{x}_{km}+\boldsymbol{\mathcal{B}}_{k}\boldsymbol{W}_{k},$
(4)

where

$\displaystyle\boldsymbol{\mathcal{A}}_{k}$
$\displaystyle:=\begin{bmatrix}\boldsymbol{A}_{km}\\\
\boldsymbol{A}_{km+1}\boldsymbol{A}_{km}\\\ \vdots\\\
\displaystyle\prod_{i=0}^{m-1}\boldsymbol{A}_{km+i}\end{bmatrix},$ (5)
$\displaystyle\boldsymbol{\mathcal{B}}_{k}$
$\displaystyle:=\begin{bmatrix}\boldsymbol{B}_{km}&0&...&0\\\
\boldsymbol{A}_{km+1}\boldsymbol{B}_{km}&\boldsymbol{B}_{km+1}&...&0\\\
\vdots&\vdots&&\vdots\\\
\displaystyle\prod_{i=1}^{m-1}\boldsymbol{A}_{km+i}\boldsymbol{B}_{km}&...&...&\boldsymbol{B}_{km+m-1}\end{bmatrix},$
(6) $\displaystyle\boldsymbol{W}_{k}$
$\displaystyle:=\begin{pmatrix}\boldsymbol{w}_{km}\\\ \vdots\\\
\boldsymbol{w}_{km+m-1}\end{pmatrix},$ (7)

with the product terms in $\boldsymbol{\mathcal{A}}_{k}$ and
$\boldsymbol{\mathcal{B}}_{k}$ defined as

$\displaystyle\prod_{i=i_{1}}^{i_{2}}\boldsymbol{A}_{k+i}:=\boldsymbol{A}_{k+i_{2}}\times\cdots\times\boldsymbol{A}_{k+i_{1}+1}\boldsymbol{A}_{k+i_{1}}\text{if
$i_{2}\geq i_{1}$}.$ (8)

The augment measurement model from (1b) is

$\displaystyle\boldsymbol{Y}_{k}$
$\displaystyle=\boldsymbol{\mathcal{C}}_{k}\boldsymbol{X}_{k}+\boldsymbol{N}_{k},$
(9)

where

$\displaystyle\boldsymbol{\mathcal{C}}_{k}$
$\displaystyle:=\mathbf{diag}\left(\boldsymbol{C}_{km+1},...,\boldsymbol{C}_{(k+1)m}\right),$
(10)

and

$\displaystyle\boldsymbol{N}_{k}$
$\displaystyle:=\begin{pmatrix}\boldsymbol{n}_{km+1}\\\ \vdots\\\
\boldsymbol{n}_{(k+1)m}\end{pmatrix}.$ (11)

We define the augmented process noise and observation noise variances as:

$\displaystyle\boldsymbol{\mathcal{Q}}_{k}:=\mathbb{E}\left[\boldsymbol{W}_{k}{\boldsymbol{W}_{k}}^{T}\right]=\mathbf{diag}\left(\boldsymbol{Q}_{km},\cdots,\boldsymbol{Q}_{km+m-1}\right),$

and

$\displaystyle\boldsymbol{\mathcal{R}}_{k}:=\mathbb{E}\left[\boldsymbol{N}_{k}{\boldsymbol{N}_{k}}^{T}\right]=\mathbf{diag}\left(\boldsymbol{R}_{km},\cdots,\boldsymbol{R}_{km+m-1}\right),$
(12)

where

$\displaystyle\boldsymbol{Q}_{km+j}$
$\displaystyle:=\mathbb{E}\left[\boldsymbol{w}_{km+j}\boldsymbol{w}^{T}_{km+j}\right],$
$\displaystyle\boldsymbol{R}_{km+j}$
$\displaystyle:=\mathbb{E}\left[\boldsymbol{n}_{km+j}\boldsymbol{n}^{T}_{km+j}\right].$

###  III-B Uncertainty Propagation and Measurement Update

The prior statistics of $\boldsymbol{X}_{k}$ are related to the posterior
statistics of $\boldsymbol{x}_{km}$ i.e.
$(\boldsymbol{\mu}_{km}^{+},\boldsymbol{\Sigma}^{+}_{km})$ as

$\displaystyle\boldsymbol{\bar{X}}_{k}^{-}$
$\displaystyle:=\mathbb{E}\left[\boldsymbol{X}_{k}\right]=\boldsymbol{\mathcal{A}}_{k}\boldsymbol{\mu}_{km}^{+},$
(13) $\displaystyle\boldsymbol{P}^{-}_{k}$
$\displaystyle:=\mathbb{E}\left[(\boldsymbol{X}_{k}^{-}-\boldsymbol{\bar{X}}_{k}^{-})(\boldsymbol{X}_{k}^{-}-\boldsymbol{\bar{X}}_{k}^{-})^{T}\right],$
$\displaystyle=\boldsymbol{\mathcal{A}}_{k}\boldsymbol{\Sigma}^{+}_{km}\boldsymbol{\mathcal{A}}_{k}^{T}+\boldsymbol{\mathcal{B}}_{k}\boldsymbol{\mathcal{Q}}_{k}\boldsymbol{\mathcal{B}}_{k}^{T}.$
(14)

Prior statistics of the augmented state, i.e. $\boldsymbol{\bar{X}}_{k}^{-}$
and $\boldsymbol{P}^{-}_{k}$, can be updated using the augmented measurements
$\boldsymbol{Y}_{k}$ to obtain posterior
$(\boldsymbol{\bar{X}}_{k}^{+},\boldsymbol{P}^{+}_{k})$, using similar steps
as in standard Kalman filtering, i.e.

$\displaystyle\boldsymbol{\bar{X}}_{k}^{+}$
$\displaystyle:=\boldsymbol{\mathcal{A}}_{k}\boldsymbol{\mu}^{+}_{km}+\boldsymbol{\mathcal{K}}_{k}(\boldsymbol{Y}_{k}-\boldsymbol{\mathcal{C}}_{k}\boldsymbol{\mathcal{A}}_{k}\boldsymbol{\mu}^{+}_{km}),$
(15) $\displaystyle\boldsymbol{P}^{+}_{k}$
$\displaystyle:=(\mathbf{I}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{k})\boldsymbol{P}^{-}_{k},$
(16)

where

$\displaystyle\boldsymbol{\mathcal{K}}_{k}:=\boldsymbol{P}^{-}_{k}\boldsymbol{\mathcal{C}}_{k}^{T}\Big{[}\boldsymbol{\mathcal{C}}_{k}\boldsymbol{P}^{-}_{k}\boldsymbol{\mathcal{C}}_{k}^{T}+\boldsymbol{\mathcal{R}}_{k}\Big{]}^{-1}.$
(17)

The state at time $t_{(k+1)m}$ can be determined from $\boldsymbol{X}_{k}$ as

$\boldsymbol{x}_{(k+1)m}:=\boldsymbol{M}_{m}\boldsymbol{X}_{k},$

where

$\displaystyle\boldsymbol{M}_{m}:=\begin{bmatrix}\boldsymbol{0}_{n_{x}\times
n_{x}(m-1)}&\boldsymbol{I}_{n_{x}}\end{bmatrix}.$ (18)

The posterior statistics of $\boldsymbol{x}_{(k+1)m}$ can then be determined
from the posterior statistics of $\boldsymbol{X}_{k}$ using

$\displaystyle\boldsymbol{\mu}_{(k+1)m}^{+}:=\boldsymbol{M}_{m}\boldsymbol{\bar{X}}_{k}^{+},\text{
and
}\boldsymbol{\Sigma}_{(k+1)m}^{+}:=\boldsymbol{M}_{m}\boldsymbol{P}^{+}_{k}\boldsymbol{M}_{m}^{T}.$
(19)

###  III-C Optimal Sensor Precision for a Single Measurement Update

Here we present a convex optimization framework for determining the nosiest
sensors, for which the estimation errors are below a given upper bound after
one measurement update. That is, maximize
$\mathbf{tr}\left(\boldsymbol{\mathcal{R}}_{k}\right)$ or minimize
$\mathbf{tr}\left(\boldsymbol{\mathcal{S}}_{k}\right)$ where
$\boldsymbol{\mathcal{S}}_{k}:=\boldsymbol{\mathcal{R}}_{k}^{-1}$, for which
$\mathbf{tr}\left(\boldsymbol{\Sigma}_{(k+1)m}^{+}\right)\leq\gamma_{d}$,
given $\boldsymbol{\Sigma}_{km}^{+}$. This is achieved by solving the
following optimization problem.

######  Theorem 1.

Optimal sensor precision $\boldsymbol{s}_{k}\in^{n_{y_{k,m}}}\geq 0$, which
satisfies
$\mathbf{tr}\left(\boldsymbol{\Sigma}_{(k+1)m}^{+}\right)\leq\gamma_{d}$, is
given by the solution of the following optimization problem,

$\left.\begin{aligned}
&\min_{\boldsymbol{s}_{k},\boldsymbol{\mathcal{K}}_{k},\boldsymbol{F}}{\mathbf{tr}\left(\boldsymbol{W}\boldsymbol{\mathcal{S}}_{k}\right)},\text{
subject to }\\\
&\begin{bmatrix}\boldsymbol{F}&\boldsymbol{M}_{12}\sqrt{\boldsymbol{P}^{-}_{k}}&\boldsymbol{\mathcal{K}}_{k}\\\
(\ast)^{T}&\boldsymbol{I}_{mn_{x}}&\boldsymbol{0}_{mn_{x}\times
n_{y_{k,m}}}\\\
(\ast)^{T}&(\ast)^{T}&\boldsymbol{\mathcal{S}}_{k}\end{bmatrix}\geq 0\\\
&0\leq\boldsymbol{s}_{k}\leq\boldsymbol{s}_{k}^{\text{max}},\end{aligned}\right\\}$
(20)

where
$\boldsymbol{M}_{12}:=\boldsymbol{M}_{m}(\boldsymbol{I}_{mn_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{k})$
, $\boldsymbol{\mathcal{K}}_{k}$ and $\boldsymbol{F}$ are the design variables
with $\mathbf{tr}\left(\boldsymbol{F}\right)\leq\gamma_{d}$,
$\boldsymbol{M}_{n}$ and system parameter $\boldsymbol{\mathcal{C}}_{k}$ is
defined in (18) and (10). The variable $\boldsymbol{W}$ is a diagonal matrix,
which is user defined and serves as a normalizing weight on
$\boldsymbol{\mathcal{S}}_{k}:=\boldsymbol{\mathcal{R}}_{k}^{-1}$, where
$\boldsymbol{\mathcal{R}}_{k}$ is defined in (12).

###### Proof.

We can write the posterior error variance as

$\displaystyle\boldsymbol{P}^{+}_{k}=(\boldsymbol{I}_{mn_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{k})\boldsymbol{P}^{-}_{k}(\boldsymbol{I}_{mn_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{k})^{T}+\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{R}}_{k}\boldsymbol{\mathcal{K}}_{k}^{T}.$
(21)

The optimal $\boldsymbol{\mathcal{K}}_{k}$ is determined by minimizing
$\mathbf{tr}\left(\boldsymbol{P}^{+}_{k}\right)$ and is given by (17).
However, in this formulation, we leave
$\boldsymbol{\mathcal{K}}_{k}\in^{mn_{x}\times n_{y_{k,m}}}$ as a variable,
and write
$\mathbf{tr}\left(\boldsymbol{M}_{m}\boldsymbol{P}^{+}_{k}\boldsymbol{M}_{m}^{T}\right)\leq\gamma_{d}$
equivalently as

$\displaystyle\boldsymbol{F}-\boldsymbol{M}_{12}\boldsymbol{P}^{-}_{k}\boldsymbol{M}_{12}^{T}-\boldsymbol{M}_{m}\boldsymbol{\mathcal{K}}_{k}^{T}\boldsymbol{\mathcal{R}}_{k}\boldsymbol{\mathcal{K}}_{k}^{T}\boldsymbol{M}^{T}_{m}\geq
0,$ and $\displaystyle\mathbf{tr}\left(\boldsymbol{F}\right)\leq\gamma_{d},$

where
$\boldsymbol{M}_{12}:=\boldsymbol{M}_{m}(\boldsymbol{I}_{mn_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{k})$
and $\boldsymbol{F}\in\mathbb{S}^{n_{x}}_{+}$. Representing
$\sqrt{\boldsymbol{P}^{-}_{k}}$ as the principal matrix square-root of
$\boldsymbol{P}^{-}_{k}$, substituting
$\boldsymbol{\mathcal{S}}_{k}:=\boldsymbol{\mathcal{R}}_{k}^{-1}$ which is
defined in (12) , and using Schur complement we get the following LMI [21],

$\displaystyle\begin{bmatrix}\boldsymbol{F}&\boldsymbol{M}_{12}\sqrt{\boldsymbol{P}^{-}_{k}}&\boldsymbol{M}_{m}\boldsymbol{\mathcal{K}}_{k}\\\
(\ast)^{T}&\boldsymbol{I}_{mn_{x}}&\boldsymbol{0}_{mn_{x}\times
n_{y_{k,m}}}\\\
(\ast)^{T}&(\ast)^{T}&\boldsymbol{\mathcal{S}}_{k}\end{bmatrix}\geq 0.$

Combining the inequalities and minimizing
$\mathbf{tr}\left(\boldsymbol{W}\boldsymbol{\mathcal{S}}_{k}\right)$ we get
(20). ∎

######  Remark 1.

In theorem 1 the sparseness of the solution can be improved by iteratively
solving the optimization problem (20) with weights
$\boldsymbol{W}_{j+1}:=(\boldsymbol{\mathcal{S}}_{k,j}^{\ast}+\epsilon\boldsymbol{I})^{-1}$,
with $\boldsymbol{W}_{1}:=\boldsymbol{I}_{n_{y_{k,m}}}$ and $\boldsymbol{I}$
is the identity matrix of appropriate dimension, where subscript $j$ denotes
the iteration index [22, 23] and $\boldsymbol{\mathcal{S}}_{k,j}^{\ast}$ is
the solution to (20) in the $j^{\text{th}}$ iteration of improving sparsity in
the solution.

##  IV Optimal Sensor Precision for Bounded Steady-State Errors

In this section, we present the result that determines the optimal sensor
precision for bounded steady-state error, assuming the system to be
$m$-periodic. If the system in (1a) is $m$-periodic, i.e.
$\boldsymbol{A}_{km+j}=\boldsymbol{A}_{(k+1)m+j}$,
$\boldsymbol{B}_{km+j}=\boldsymbol{B}_{(k+1)m+j}$, and
$\boldsymbol{C}_{km+j}=\boldsymbol{C}_{(k+1)m+j}$ for $j=1,\cdots,m$; it will
be of interest to determine the sensing precisions that bound the steady-state
errors, assuming it exists.

###  IV-A Augmented Dynamical System

From (4), the augmented dynamics of the $m$-periodic system is given by,

$\displaystyle\boldsymbol{x}_{(k+1)m}$
$\displaystyle=\boldsymbol{M}_{m}\boldsymbol{\mathcal{A}}_{k}\boldsymbol{x}_{km}+\boldsymbol{M}_{m}\boldsymbol{\mathcal{B}}_{k}\boldsymbol{W}_{k}.$
(22)

In this section we generalize the sensor model in (1b) by including the
process noise in the measurement[24]. This scenario, for example, occurs in
measurements from accelerometers where the disturbance forces algebraically
impact the accelerations. The new measurement model is therefore,

$\boldsymbol{y}_{k}=\boldsymbol{C}_{k}\boldsymbol{x}_{k}+\boldsymbol{D}_{k}\boldsymbol{w}_{k}+\boldsymbol{n}_{k}.$
(23)

Consequently, the augmented sensor model is

$\displaystyle\boldsymbol{Y}_{k}$
$\displaystyle=\boldsymbol{\mathcal{C}}_{k}\boldsymbol{X}_{k}+\boldsymbol{\mathcal{D}}_{k}\boldsymbol{W}_{k}+\boldsymbol{N}_{k},$
$\displaystyle=\boldsymbol{\mathcal{C}}_{k}(\boldsymbol{\mathcal{A}}_{k}\boldsymbol{x}_{km}+\boldsymbol{\mathcal{B}}_{k}\boldsymbol{W}_{k})+\boldsymbol{\mathcal{D}}_{k}\boldsymbol{W}_{k}+\boldsymbol{N}_{k},$
$\displaystyle=\boldsymbol{\mathcal{C}}_{k}\boldsymbol{\mathcal{A}}_{k}\boldsymbol{x}_{km}+\left(\boldsymbol{\mathcal{C}}_{k}\boldsymbol{\mathcal{B}}_{k}+\boldsymbol{\mathcal{D}}_{k}\right)\boldsymbol{W}_{k}+\boldsymbol{N}_{k},$
(24)

where
$\boldsymbol{\mathcal{D}}_{k}:=\mathbf{diag}\left(\boldsymbol{D}_{km+1},...,\boldsymbol{D}_{(k+1)m}\right)$.

The presence of $\boldsymbol{W}_{k}$ in (24) makes derivation of the Kalman
filter complicated. This is circumvented by assuming the process noise to be
colored, or filtered white noise. That is, we model the process noise as

$\displaystyle\boldsymbol{Z}_{k+1}=\boldsymbol{\mathcal{G}}\boldsymbol{Z}_{k}+\boldsymbol{\mathcal{H}}\boldsymbol{\Lambda}_{k},\;\boldsymbol{W}_{k}=\boldsymbol{Z}_{k},$
(25)

where $\boldsymbol{\Lambda}_{k}$ is white noise, $\boldsymbol{Z}_{k}$ is the
filter state, and the pair
$(\boldsymbol{\mathcal{G}}_{k},\boldsymbol{\mathcal{H}}_{k})$ defines the
filter.

If white noise $\boldsymbol{\lambda}_{k}\in^{n_{w}}$ is filtered via

$\boldsymbol{z}_{km+1}=\boldsymbol{G}\boldsymbol{z}_{km}+\boldsymbol{H}\boldsymbol{\lambda}_{km},$
(26)

then for the augmented system,

$\displaystyle\boldsymbol{Z}_{k}:=\begin{pmatrix}\boldsymbol{z}_{km}\\\
\vdots\\\
\boldsymbol{z}_{(k+1)m-1}\end{pmatrix},\;\boldsymbol{\Lambda}_{k}:=\begin{pmatrix}\boldsymbol{\lambda}_{km}\\\
\vdots\\\ \boldsymbol{\lambda}_{(k+1)m-1}\end{pmatrix},$ (27)

$\boldsymbol{\mathcal{G}}:=\boldsymbol{I}_{q}\otimes\boldsymbol{G}$, and
$\boldsymbol{\mathcal{H}}:=\boldsymbol{I}_{q}\otimes\boldsymbol{H}$, where
$\boldsymbol{G}\in^{n_{w}\times n_{w}}$ and $\boldsymbol{H}\in^{n_{w}\times
n_{w}}$ define the filter in (26).

Introducing a new state variable

$\displaystyle\boldsymbol{\Gamma}_{k}:=\begin{bmatrix}\boldsymbol{x}_{km}\\\
\boldsymbol{Z}_{k}\end{bmatrix}\in^{N_{x}},$ (28)

where $N_{x}:=n_{x}+mn_{w}$, we can write the dynamics of
$\boldsymbol{\Gamma}_{k}$ and measurement $\boldsymbol{Y}_{k}$ as

$\displaystyle\boldsymbol{\Gamma}_{k+1}=\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\Gamma}_{k}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\Lambda}_{k},$
(29a)
$\displaystyle\boldsymbol{Y}_{k}=\boldsymbol{\mathcal{C}}_{m}\boldsymbol{\Gamma}_{k}+\boldsymbol{N}_{k}.$
(29b)

where

$\displaystyle\boldsymbol{\mathcal{A}}_{m}$
$\displaystyle:=\begin{bmatrix}\boldsymbol{M}_{m}\boldsymbol{\mathcal{A}}_{k}&\boldsymbol{M}_{m}\boldsymbol{\mathcal{B}}_{k}\\\
\boldsymbol{0}_{qn_{w}\times n_{x}}&\boldsymbol{\mathcal{G}}\end{bmatrix},$
(30a) $\displaystyle\boldsymbol{\mathcal{B}}_{m}$
$\displaystyle:=\begin{bmatrix}\boldsymbol{0}_{n_{x}\times qn_{w}}\\\
\boldsymbol{\mathcal{H}}\end{bmatrix},$ (30b)
$\displaystyle\boldsymbol{\mathcal{C}}_{m}$
$\displaystyle:=\begin{bmatrix}\boldsymbol{\mathcal{C}}_{k}\boldsymbol{\mathcal{A}}_{k}&(\boldsymbol{\mathcal{C}}_{k}\boldsymbol{\mathcal{B}}_{k}+\boldsymbol{\mathcal{D}}_{k})\end{bmatrix}.$
(30c)

Note that for the $m$-periodic system, the matrices
$\boldsymbol{\mathcal{A}}_{m}$, $\boldsymbol{\mathcal{B}}_{m}$, and
$\boldsymbol{\mathcal{C}}_{m}$ are time invariant.

States $\boldsymbol{x}_{km}$ can be recovered from $\boldsymbol{\Gamma}_{k}$
as

$\displaystyle\boldsymbol{x}_{km}=\boldsymbol{M}_{x}\boldsymbol{\Gamma}_{k},$
(31)

where $\boldsymbol{M}_{x}$ is a mask-matrix defined by

$\displaystyle\boldsymbol{M}_{x}:=\begin{bmatrix}\boldsymbol{I}_{n_{x}}&\boldsymbol{0}_{n_{x}\times
qn_{w}}\end{bmatrix}.$ (32)

We next define

$\displaystyle\boldsymbol{\mathcal{Q}}_{m}:=\mathbb{E}\left[\boldsymbol{\Lambda}_{k}\boldsymbol{\Lambda}_{k}^{T}\right],\text{
and
}\boldsymbol{\mathcal{R}}_{m}:=\mathbb{E}\left[\boldsymbol{N}_{k}\boldsymbol{N}^{T}_{k}\right].$
(33)

For steady-state analysis, we assume
$(\boldsymbol{\mathcal{A}}_{m},\boldsymbol{\mathcal{R}}_{m}^{1/2}\boldsymbol{\mathcal{C}}_{m})$
is detectable and
$(\boldsymbol{\mathcal{A}}_{m},(\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}^{T}_{m})^{1/2})$
is stabilizable [18].

###  IV-B Steady-state Variance

Let $\bar{\boldsymbol{\Gamma}}^{-}_{k}$ and $\boldsymbol{P}^{-}_{k}$ be the
prior mean and variance of $\boldsymbol{\Gamma}_{k}$ at time $k$. This defines
the prior random variable
$\boldsymbol{\Gamma}_{k}^{-}\sim\mathcal{N}(\bar{\boldsymbol{\Gamma}}^{-}_{k},\boldsymbol{P}^{-}_{k})$,
where $\mathcal{N}(\cdot,\cdot)$ defines a Gaussian distribution.

In Kalman filtering we assume the posterior is a linear function of the prior
and the measurement, i.e.

$\displaystyle\boldsymbol{\Gamma}^{+}_{k}:=(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{\Gamma}^{-}_{k}+\boldsymbol{\mathcal{K}}_{k}\boldsymbol{Y}_{k},$
(34)

where $\boldsymbol{\mathcal{K}}_{k}\in^{N_{x}\times n_{y_{k,m}}}$ is the
unknown gain.

The coefficient $\boldsymbol{\mathcal{K}}_{k}$ is determined by minimizing the
posterior variance. However, in this formulation, we leave it as a free
variable along with $\boldsymbol{\mathcal{R}}_{m}$. Both these variables will
be jointly determined in a single optimization problem, presented in theorem
2.

Using (34), the posterior variance is given by

$\displaystyle\boldsymbol{P}^{+}_{k}=(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{P}^{-}_{k}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{m})^{T}+\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{R}}_{m}\boldsymbol{\mathcal{K}}^{T}_{k}.$
(35)

Using (29a), the prior mean and variance of $\boldsymbol{\Gamma}_{k}$ at time
$k+1$ is given by

$\displaystyle\bar{\boldsymbol{\Gamma}}^{-}_{k+1}$
$\displaystyle=\boldsymbol{\mathcal{A}}_{m}\bar{\boldsymbol{\Gamma}}^{+}_{k},$
(36a) $\displaystyle\boldsymbol{P}^{-}_{k+1}$
$\displaystyle=\boldsymbol{\mathcal{A}}_{m}\boldsymbol{P}^{+}_{k}\boldsymbol{\mathcal{A}}^{T}_{m}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T},$
(36b)

which defines the random variable
$\boldsymbol{\Gamma}_{k+1}^{-}\sim\mathcal{N}(\bar{\boldsymbol{\Gamma}}^{-}_{k+1},\boldsymbol{P}^{-}_{k+1})$.

Replacing $\boldsymbol{P}^{-}_{k}$ from (35) in (36b), we get the propagation
equation for the prior variance

$\displaystyle\boldsymbol{P}^{-}_{k+1}=$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{P}^{-}_{k}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{C}}_{m})^{T}\boldsymbol{\mathcal{A}}^{T}_{m}$
$\displaystyle+\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{k}\boldsymbol{\mathcal{R}}_{m}\boldsymbol{\mathcal{K}}^{T}_{k}\boldsymbol{\mathcal{A}}^{T}_{m}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T}.$
(37)

Steady-state variance $\boldsymbol{P}_{\infty}$ is determined by solving

$\displaystyle\boldsymbol{P}_{\infty}=$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{P}_{\infty}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})^{T}\boldsymbol{\mathcal{A}}^{T}_{m}$
$\displaystyle+\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{R}}_{m}\boldsymbol{\mathcal{K}}^{T}_{\infty}\boldsymbol{\mathcal{A}}^{T}_{m}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T},$
(38)

where $\boldsymbol{\mathcal{K}}_{\infty}$ is the steady-state gain. The
steady-state variance of $\boldsymbol{x}_{(k+1)m}$ is then given by
$\boldsymbol{\Sigma}_{\infty}:=\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}$.

######  Remark 2.

Equations (37) becomes the Riccati difference equation (RDE) if
$\boldsymbol{\mathcal{K}}_{k}$ is determined by minimizing the posterior
variance. Consequently, it transforms (38) to the algebraic Riccati equation
(ARE). That is, for $\boldsymbol{\mathcal{K}}_{k}$ given by (17), (37)
transforms to

$\displaystyle\boldsymbol{P}_{k+1}=$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{P}_{k}-\boldsymbol{P}_{k}\boldsymbol{\mathcal{C}}_{m}^{T}[\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}_{k}\boldsymbol{\mathcal{C}}_{m}^{T}+\boldsymbol{\mathcal{R}}_{m}]^{-1}\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}_{k})\boldsymbol{\mathcal{A}}_{m}^{T}$
$\displaystyle+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}^{T}_{m},$
(39)

and (38) transforms to

$\displaystyle\boldsymbol{P}_{\infty}=$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{P}_{\infty}-\boldsymbol{P}_{\infty}\boldsymbol{\mathcal{C}}_{m}^{T}[\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}_{\infty}\boldsymbol{\mathcal{C}}_{m}^{T}+\boldsymbol{\mathcal{R}}_{m}]^{-1}\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}_{\infty})$
$\displaystyle\times\boldsymbol{\mathcal{A}}_{m}^{T}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}^{T}_{m},$
(40)

Equation (40) has a unique positive semi-definite solution if
$(\boldsymbol{\mathcal{A}}_{m},\boldsymbol{\mathcal{R}}_{m}^{1/2}\boldsymbol{\mathcal{C}}_{m})$
is detectable and
$(\boldsymbol{\mathcal{A}}_{m},(\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}^{T}_{m})^{1/2})$
is stabilizable.

###  IV-C Optimal Sensor Precision for Bounded Steady-State Estimation Error

Let $\boldsymbol{P}^{d}_{\infty}$ be the desired steady-state error, and let
us assume that it is the solution of (40) for some
$\boldsymbol{\mathcal{R}}^{d}_{m}$, i.e.

$\displaystyle\boldsymbol{P}^{d}_{\infty}=$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{P}^{d}_{\infty}-\boldsymbol{P}^{d}_{\infty}\boldsymbol{\mathcal{C}}_{m}^{T}\left[\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}^{d}_{\infty}\boldsymbol{\mathcal{C}}_{m}^{T}+\boldsymbol{\mathcal{R}}^{d}_{m}\right]^{-1}\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}^{d}_{\infty})$
$\displaystyle\times\boldsymbol{\mathcal{A}}_{m}^{T}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}^{T}_{m}.$

Therefore, for any
$\boldsymbol{\mathcal{R}}_{m}-\boldsymbol{\mathcal{R}}^{d}_{m}\leq 0$

$\displaystyle\boldsymbol{P}^{d}_{\infty}\geq$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{P}^{d}_{\infty}-\boldsymbol{P}^{d}_{\infty}\boldsymbol{\mathcal{C}}_{m}^{T}[\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}^{d}_{\infty}\boldsymbol{\mathcal{C}}_{m}^{T}+\boldsymbol{\mathcal{R}}_{m}]^{-1}\boldsymbol{\mathcal{C}}_{m}\boldsymbol{P}^{d}_{\infty})$
$\displaystyle\times\boldsymbol{\mathcal{A}}_{m}^{T}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}^{T}_{m},$
(41)

which makes the solution of the RDE monotonic.

According to Lemma 2 in [25], if for some $k$ the solution of the RDE in (39)
is monotonic, i.e. $\boldsymbol{P}_{k}\geq\boldsymbol{P}_{k+1}$, then
$\boldsymbol{P}_{k+i}\geq\boldsymbol{P}_{k+i+1}$ for all $i\geq 1$. Therefore,
(41), guarantees that the evolution of $\boldsymbol{P}^{d}_{\infty}$ is
monotonic, and $\boldsymbol{P}_{\infty}$ is guaranteed to satisfy
$\boldsymbol{P}^{d}_{\infty}\geq\boldsymbol{P}_{\infty}$. The gap in the
inequality can be minimized by maximizing
$\mathbf{tr}\left(\boldsymbol{\mathcal{R}}_{m}\right)$.

The optimization problem that determines the optimal sensor precision to
guarantee
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}^{d}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$,
for a given $\gamma_{d}$ is presented as the following theorem.

######  Theorem 2.

Optimal sensor precision $\boldsymbol{s}\in^{n_{y_{k,m}}}\geq 0$, which
satisfies
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}^{d}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$
is given by the solution of the following optimization problem,

$\displaystyle\min_{\boldsymbol{s}_{m},\boldsymbol{Z},\boldsymbol{P}^{d}_{\infty},\boldsymbol{\mathcal{K}}_{\infty}}{\mathbf{tr}\left(\boldsymbol{W}\boldsymbol{\mathcal{S}}_{m}\right)}\text{
subject to }$ (42a)
$\displaystyle\begin{bmatrix}\boldsymbol{M}_{11}&\boldsymbol{M}_{x}\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})&\boldsymbol{M}_{x}\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{\infty}\\\
(\ast)^{T}&\boldsymbol{Z}&\boldsymbol{0}_{N_{x}\times n_{y_{k,m}}}\\\
(\ast)^{T}&(\ast)^{T}&\boldsymbol{\mathcal{S}}_{m}\end{bmatrix}\geq 0,$ (42b)
$\displaystyle\begin{bmatrix}2\boldsymbol{I}_{N_{x}}&\boldsymbol{P}^{d}_{\infty}&\boldsymbol{Z}\\\
\boldsymbol{P}^{d}_{\infty}&\frac{1}{\delta}\boldsymbol{I}_{N_{x}}&\boldsymbol{0}_{N_{x}\times
N_{x}}\\\ \boldsymbol{Z}&\boldsymbol{0}_{N_{x}\times
N_{x}}&\delta\boldsymbol{I}_{N_{x}}\end{bmatrix}\geq 0,$ (42c)
$\displaystyle\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}^{d}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d},$
(42d) $\displaystyle 0\leq\boldsymbol{s}_{m}\leq\boldsymbol{s}_{\text{max}},$
(42e)

$\boldsymbol{M}_{11}:=\boldsymbol{M}_{x}\left(\boldsymbol{P}^{d}_{\infty}-\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T}\right)\boldsymbol{M}^{T}_{x}$,
$\boldsymbol{M}_{x}$ defined in (32), $\boldsymbol{\mathcal{B}}_{m}$ and
$\boldsymbol{\mathcal{Q}}_{m}$ are system parameters,
$\boldsymbol{\mathcal{S}}_{m}:=\mathbf{diag}\left(\boldsymbol{s}_{m}\right)$,
$n_{s}:=\sum_{j=1}^{m}p_{j}$, $\boldsymbol{Z}\in\mathbb{S}_{+}^{N_{x}}$ is a
design variable, and $p_{j}$ is the dimension of the j${}^{\text{th}}$ sensor.
Variables $\gamma_{d}>0$ and $\delta>0$ are user specified. The variable
$\boldsymbol{W}$ is a diagonal matrix, which is also user defined, and serves
as a normalizing weight on $\boldsymbol{\mathcal{S}}_{m}$.

###### Proof.

From (37), monotonicity of $\boldsymbol{P}_{\infty}^{d}$ is guaranteed if

$\displaystyle\boldsymbol{P}^{d}_{\infty}\geq$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{P}^{d}_{\infty}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})^{T}\boldsymbol{\mathcal{A}}^{T}_{m}$
$\displaystyle+\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{R}}_{m}\boldsymbol{\mathcal{K}}^{T}_{\infty}\boldsymbol{\mathcal{A}}^{T}_{m}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T}.$

Introducing a new variable $\boldsymbol{Z}\in\mathbb{S}_{+}^{N_{x}}$, and the
relaxation

$\boldsymbol{Z}^{-1}\geq\boldsymbol{P}^{d}_{\infty},$

the condition for monotonicity can then be written as

$\displaystyle\boldsymbol{P}^{d}_{\infty}\geq$
$\displaystyle\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{Z}^{-1}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})^{T}\boldsymbol{\mathcal{A}}^{T}_{m}$
$\displaystyle+\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{R}}_{m}\boldsymbol{\mathcal{K}}^{T}_{\infty}\boldsymbol{\mathcal{A}}^{T}_{m}+\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T}.$

However, we want to enforce monotonicity of
$\boldsymbol{M}_{x}\boldsymbol{P}^{d}_{\infty}\boldsymbol{M}_{x}^{T}$, i.e.

$\displaystyle\boldsymbol{M}_{x}\boldsymbol{P}^{d}_{\infty}\boldsymbol{M}^{T}_{x}\geq$
$\displaystyle\boldsymbol{M}_{x}\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{M}_{x}\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})\boldsymbol{Z}^{-1}$
$\displaystyle\times(\boldsymbol{I}_{N_{x}}-\boldsymbol{M}_{x}\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})^{T}\boldsymbol{\mathcal{A}}^{T}_{m}\boldsymbol{M}^{T}_{x}$
$\displaystyle+\boldsymbol{M}_{x}\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{R}}_{m}\boldsymbol{\mathcal{K}}^{T}_{\infty}\boldsymbol{\mathcal{A}}^{T}_{m}\boldsymbol{M}^{T}_{x}$
$\displaystyle+\boldsymbol{M}_{x}\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T}\boldsymbol{M}^{T}_{x}.$

Using Schur complement [26], and substituting
$\boldsymbol{\mathcal{S}}_{m}:=\boldsymbol{\mathcal{R}}_{m}^{-1}$, we get

$\begin{bmatrix}\boldsymbol{M}_{11}&\boldsymbol{M}_{x}\boldsymbol{\mathcal{A}}_{m}(\boldsymbol{I}_{N_{x}}-\boldsymbol{\mathcal{K}}_{\infty}\boldsymbol{\mathcal{C}}_{m})&\boldsymbol{M}_{x}\boldsymbol{\mathcal{A}}_{m}\boldsymbol{\mathcal{K}}_{\infty}\\\
(\ast)^{T}&\boldsymbol{Z}&\boldsymbol{0}_{N_{x}\times n_{y_{k,m}}}\\\
(\ast)^{T}&(\ast)^{T}&\boldsymbol{\mathcal{S}}_{m}\end{bmatrix}\geq 0,$

where
$\boldsymbol{M}_{11}:=\boldsymbol{M}_{x}\left(\boldsymbol{P}^{d}_{\infty}-\boldsymbol{\mathcal{B}}_{m}\boldsymbol{\mathcal{Q}}_{m}\boldsymbol{\mathcal{B}}_{m}^{T}\right)\boldsymbol{M}^{T}_{x}$.

The relaxation $\boldsymbol{Z}^{-1}\geq\boldsymbol{P}^{d}_{\infty}$ can be
written as
$\boldsymbol{P}^{d}_{\infty}\boldsymbol{Z}\leq\boldsymbol{I}_{N_{x}}$, which
is non convex. However, we know that

$\boldsymbol{P}^{d}_{\infty}\boldsymbol{Z}+\boldsymbol{Z}\boldsymbol{P}^{d}_{\infty}\leq\delta\boldsymbol{P}^{d}_{\infty}\boldsymbol{P}^{d}_{\infty}+\frac{1}{\delta}\boldsymbol{Z}\boldsymbol{Z},$

for a given $\delta$ following from a special case of Young’s relation [26].
Therefore,

$\delta\boldsymbol{P}^{d}_{\infty}\boldsymbol{P}^{d}_{\infty}+\frac{1}{\delta}\boldsymbol{Z}\boldsymbol{Z}\leq
2\boldsymbol{I}_{N_{x}},$ (43)

guarantees
$\boldsymbol{P}^{d}_{\infty}\boldsymbol{Z}\leq\boldsymbol{I}_{N_{x}}$. The
inequality in (43), can be written as the following linear matrix inequality

$\begin{bmatrix}2\boldsymbol{I}_{N_{x}}&\boldsymbol{P}^{d}_{\infty}&\boldsymbol{Z}\\\
\boldsymbol{P}^{d}_{\infty}&\frac{1}{\delta}\boldsymbol{I}_{N_{x}}&\boldsymbol{0}_{N_{x}\times
N_{x}}\\\ \boldsymbol{Z}&\boldsymbol{0}_{N_{x}\times
N_{x}}&\delta\boldsymbol{I}_{N_{x}}\end{bmatrix}\geq 0.$

The optimal precision is given by minimizing
$\mathbf{tr}\left(\boldsymbol{\mathcal{S}}_{m}\right)$. ∎

######  Remark 3.

The parameter $\delta$ can be tweaked to improve the solution, using
techniques from successive convex over-bounding techniques described in [27].

######  Remark 4.

Like in theorem 1, the sparseness of the solution can be improved by
iteratively solving the optimization problem in theorems 2 with weights
$\boldsymbol{W}_{j+1}:=(\boldsymbol{\mathcal{S}}_{m}^{\ast})^{-1}_{j}$, with
$\boldsymbol{W}_{1}:=\boldsymbol{I}_{n_{s}}$, where subscript $j$ denotes the
iteration index and $(\boldsymbol{\mathcal{S}}_{m}^{\ast})_{j}$ denotes the
optimal $\boldsymbol{\mathcal{S}}_{m}$ calculated at the $j^{\text{th}}$
iteration.

######  Remark 5.

The optimal $\boldsymbol{\mathcal{S}}^{\ast}_{m}$ can be conservative. The
actual steady-state variance, denoted by $\boldsymbol{P}_{\infty}$, can be
much smaller than $\boldsymbol{P}^{d}_{\infty}$. At the same time, it is
possible that $\boldsymbol{P}_{\infty}$ is smaller than
$\boldsymbol{P}^{d}_{\infty}$ without requiring $\boldsymbol{P}^{d}_{\infty}$
to be monotonic. Thus, theorem 2 is conservative and this results in more
precision than that required to achieve $\boldsymbol{P}^{d}_{\infty}$.

We remove the conservativeness by applying the following strategy. The
optimization problem in theorem 2 is solved to determine the sparse solution
$\boldsymbol{\mathcal{S}}^{\ast}_{m}$ for which $\boldsymbol{P}^{d}_{\infty}$
is contractive. It is then scaled by $\xi$ to reduced the gap between
$\boldsymbol{P}^{d}_{\infty}$ and $\boldsymbol{P}_{\infty}$, where optimal
$\xi^{\ast}$ is obtained using the bisection algorithm described in Algorithm
1.

Algorithm 1 Bisection algorithm for optimal scaling of sensor precision.

Define: $\xi_{\text{min}}:=0$

Define: $\xi_{\text{max}}:=10^{3}$ # Something large

Solve optimization problem (42) to get
$\boldsymbol{\mathcal{R}}^{\ast}_{m}:=(\boldsymbol{\mathcal{S}}^{\ast}_{m})^{-1}$.

Define: MAXITER = $100$ # Something large

for i = 1:MAXITER

$\xi=\frac{1}{2}(\xi_{\text{min}}+\xi_{\text{max}})$

$\boldsymbol{P}_{\text{ss}}:=$ solution of ARE in (38) with noise
$\xi\boldsymbol{\mathcal{R}}_{m}^{\ast}$

if
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\text{ss}}\boldsymbol{M}^{T}_{x}\right)<\gamma_{d}$

$\xi_{\text{min}}:=\xi$

else

$\xi_{\text{max}}:=\xi$

end

end

##  V Examples

Next, we apply theorems 1 and 2 to three estimation problems related to
aerospace engineering, highlighting their engineering value. The first example
demonstrates the application of theorem 2 to determine redundant sensing
architectures for linear time invariant (LTI) systems. The second example
demonstrates accurate state estimation with low-cost sensors. Finally, the
third example demonstrates sensor scheduling for a linear time-varying system.

###  V-A Time-Invariant System: Flight Control Example

Here we demonstrate practical applications of the result presented in theorem
2 to a linear time-invariant discrete-time system. It is applied to a steady-
state estimation problem for an aircraft model. We first present the details
of the aircraft model. We then present two examples, which highlight different
applications of theorem 2.

####  V-A1 Model

Let us consider the longitudinal motion model of an aircraft, where the states
of the system are velocity $V$ in $ft/s$, angle of attack $\alpha$ in $rad$,
pitch angle $\theta$ in $rad$, and pitch rate $q$ in $rad/s$, i.e.

$\boldsymbol{x}:=\begin{bmatrix}V&\alpha&\theta&q\end{bmatrix}^{T}.$ (44)

We consider onboard sensors that measure body acceleration
$\dot{u}\,(ft/s^{2})$ along roll axis, body acceleration $\dot{w}\,(ft/s^{2})$
along yaw axis, angle of attack $\alpha(rad)$, pitch rate $q(rad/s)$, and
dynamic pressure $\bar{q}:=\frac{1}{2}\rho V^{2}\,(lb/ft^{2})$, where $\rho$
is the atmospheric density. Variables $u$ and $w$ are defined as
$u:=V\cos(\alpha)$, and $w:=V\sin(\alpha)$.

Therefore, the vector of measured outputs is

$\boldsymbol{y}:=\left[\dot{u},\dot{w},\alpha,q,\bar{q}\right]^{T}.$ (45)

In a typical aircraft, these measurements are available from the
accelerometers, angle-of-attack sensors, gyro sensors, and pitot tube,
respectively.

The dynamics and measurement model is given by the following equations

$\displaystyle\frac{d\boldsymbol{x}}{dt}$
$\displaystyle=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{B}d,$ (46a)
$\displaystyle\boldsymbol{y}$
$\displaystyle=\boldsymbol{C}\boldsymbol{x}+\boldsymbol{D}d+\boldsymbol{n},$
(46b)

where 111More accurate data is available upon request.

$\displaystyle\boldsymbol{A}=\begin{bmatrix}-0.0179&33.2244&-32.1700&0.6728\\\
-0.0001&-1.4528&0&0.9323\\\ 0&0&0&1.0000\\\
-0.0000&-4.1970&0&-1.8836\end{bmatrix},$
$\displaystyle\boldsymbol{B}=\begin{bmatrix}0.5697&-0.0029&0&-0.4670\end{bmatrix}^{T},$
$\displaystyle\boldsymbol{C}=10^{3}\times\begin{bmatrix}-0.0000&0.0332&-0.0322&0.0007\\\
-0.0001&-1.3544&0&0.8692\\\ 0&0.0010&0&0\\\ 0&0&0&0.0010\\\
0.0017&0&0&0\end{bmatrix},$
$\displaystyle\boldsymbol{D}=\begin{bmatrix}0.5697&-2.7345&0&0&0\end{bmatrix}^{T},$

$d$ is the disturbance acting on the system, and $\boldsymbol{n}$ is the
sensor noise. Note that there is a direct feed-through term because the
disturbance directly impacts the acceleration measurements. In this example,
we model the disturbance as a filtered white noise, filtered by
$\frac{1}{s/\omega_{c}+1}$, where $\omega_{c}$ is the cutoff frequency. The
disturbance is generated by a vibrating control surface in the aircraft.

Let the filter state be $x_{d}$, and the filter dynamics be given by

$\dot{x}_{d}=\omega_{c}(-x_{d}+w),\text{ and }d=x_{d},$

where $w$ is white noise with a given variance. In this example, we choose the
variance to be $5$ $\text{deg}^{2}$. This corresponds to small angular
deflections in the control surface, which is measured in degrees.

####  V-A2 Example 1: Redundancy in Sensing

Here we discretize the continuous time system with sampling time $\Delta
t:=0.01\ s$ and determine the least precision needed to achieve a steady-error
that satisfies
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$
for $\gamma_{d}:=0.1$, where $\boldsymbol{P}_{\infty}$ quantifies the actual
steady-state error. In this example, theorem 2 is applied with $m=1$. For the
filter, the cutoff frequency $\omega_{c}$ is chosen to be $10$ rad/s. Finally,
$\delta=200$ was chosen to implement the constraint in (42c).

Figure 2: Sensor precisions for the five sensors from different algorithms.
The respective units are $(ft/s^{2})^{2}$, $(ft/s^{2})^{2}$, $rad^{2}$,
$(rad/s)^{2}$, and $(lb/ft^{2})^{2}$

Fig.(2) shows the sensor precisions from the unweighted optimization
(indicated by legend “$1$”), the precisions from iteratively weighted
optimization to improve sparseness (indicated by legend “$1/\textbf{s}$”), and
finally, the scaled precision to remove the conservativeness in the optimal
solution (indicated by legend “$\textbf{s}/\xi^{\ast}$”). The top-panel in
fig.(2) shows the sparse solution, and the bottom panel shows the same data on
the logarithmic scale. We observe that out of the five sensors chosen in the
design, only two significantly contribute to the required estimation accuracy.

Iteratively weighted optimization significantly improves the sparseness in the
solution by several orders of magnitude. We also observe that the weighted
optimization solution is conservative, and the precisions can be further
reduced to get closer to the boundary of
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$.

Fig.(3(a)) shows the scaled optimal precisions for $\xi^{\ast}=64.1106$. The
precision values in fig.(3(a)) indicate that only angular velocity measurement
$q$ and dynamic pressure data $\bar{q}$ are needed in higher precision to
estimate all four states of the system with the required accuracy.

(a) Optimal scaled sensor precisions satisfying
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$
for $\gamma_{d}:=0.1$.

(b) State estimation with only $q$ and $\bar{q}$: required precisions to
achieve
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$
for $\gamma_{d}:=0.1$.

To completely remove the sensor for $\dot{u},\dot{v}$, and $\alpha$, we can
set their corresponding precisions to exactly zero, prior to $\xi$ scaling.
Since detectability and stabilizability conditions are verified to be
satisfied for this sensor configuration, (38) has a unique solution.
Therefore, using algorithm 1 we can determine the optimal $\xi$ scaling that
guarantees
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$.
For this example, we get $\xi^{\ast}=3.375$, and the scaled sensor precisions
are shown in fig.(3(b)). We see from fig.(3(b)) that the required precisions
for $q$ and $\bar{q}$ in this case are much higher than those in fig.(3(a)).

Therefore, from a sensor pruning perspective, an ad-hoc approach would be to
start with a dictionary of sensors and determine the optimal sensor precisions
using theorem 2, then assign zero precisions to those sensors with small
precisions, and finally apply algorithm 1 to arrive at the optimal precisions
of the reduced number of sensors. More sophisticated algorithms [12, 28, 14]
for sensor pruning or sensor selection can also be applied. These algorithms
assume sensor precisions are known, which can be determined from theorems 1 or
2 .

####  V-A3 Example 2: Accurate State Estimation with Low Precision Sensors

In this section, we apply theorem 2 to explore the tradeoff between sensing
rate and sensing precision. We use the same F16 example described in §V-A1. In
this example, the continuous-time model is discretized with $dt=1/1000\ s$,
using Tustin’s method [29]. The augmented system is created with pitch rate
$q=10(rad/s)$. This formulation captures a scenario where the sensor data is
available at $1$ Khz, but the state estimates are needed at $100$ Hz. We
assume that the state’s estimates are used by some control law executing at
$100$ Hz.

Figure 4: Optimal scaled precisions for high sensing rate, guaranteeing
$\mathbf{tr}\left(\boldsymbol{M}_{x}\boldsymbol{P}_{\infty}\boldsymbol{M}_{x}^{T}\right)\leq\gamma_{d}$
for $\gamma_{d}:=0.1$.

In this example, the $5$ physical sensors in (45) are treated as virtual
sensors over $10$ time steps and are assigned an unknown precision. Thus,
there are $50$ virtual sensors. Optimization in (42), results in the sparse
precisions shown in fig.(4). The $y$-axis are the five physical sensors, and
the $x$-axis are times steps from $1$ to $q$. The heat-map shows the
precisions of the $5$ sensors across the $10$ time steps. The three panels in
fig.(4) are solutions with $\boldsymbol{s}_{\text{max}}=5,2.5,1$ respectively,
and scaled by $\xi^{\ast}$ determined by Algorithm 1. They have the same
required accuracy, defined by $\gamma_{d}=0.1$. In the top panel, we see that
the required precisions for sensors with $1000$ Hz are much lower than the
results shown in fig.(3(b)), which is for $100$ Hz sensing. Therefore, less
precise data at a higher rate can achieve the same accuracy. It is also
interesting to note that the sensors with nonzero precision are different in
the two cases. In the $100$ Hz example, the angular velocity sensor is the
most precise, followed by the dynamic pressure sensor. Other sensors have very
low precision. In the $1000$ Hz example, the angular velocity sensor and angle
of attack sensors have very low precisions. Still, acceleration measurements
have relatively higher precisions, with the dynamic pressure sensor the most
precise. We also observe that the dynamic pressure sensor plays an important
role in both cases, and the angle-of-attack sensor has very low precision in
both cases.

It is also interesting to note that the sensor values at the beginning and the
end of the 10 time-step window have higher precision. However, as
$\boldsymbol{s}_{\text{max}}$ is reduced, we observe that intermediate values
of $\bar{q}$ are needed to achieve the same accuracy. We can infer from this
observation that if available sensor precision is low, we can achieve a higher
estimation accuracy by fusing data at a higher rate. Theorem 2, determines
data from which sensors are needed at a higher rate and the corresponding
precisions to achieve this accuracy. This is very useful from a real-time
scheduling perspective because from fig.(4) we can determine exactly when to
poll the sensors. This optimizes sensor polling, reduces the associated
delays, and improves real-time schedulability.

###  V-B Time-Varying System: Satellite Tracking Problem

This example applies theorem 1 to a linear time-varying discrete-time system.
Here we consider the problem of determining the optimal sensor precision for
tracking a space object with the required accuracy.

####  V-B1 Model

Here we consider a simple satellite dynamics model [30] with $J_{2}$
perturbation given by

$\displaystyle\ddot{r}$
$\displaystyle=-\frac{\mu_{E}}{r^{2}}+\dot{\theta}^{2}r+\frac{3J_{2}}{2r^{4}}\left(3\sin(\theta)^{2}-1\right),$
(47a) $\displaystyle\ddot{\theta}$
$\displaystyle=-\frac{2\dot{\theta}\dot{r}}{r}-\frac{3J_{2}}{r^{4}}\cos(\theta)\sin(\theta).$
(47b)

where $r$ is the distance of the satellite from the centre of the orbit and
$\theta$ is the angular position of the satellite in the orbit. Length and
time in the dynamics are normalized using $R_{E}$ (radius of Earth), and
$T_{p}$ (time for one orbit) respectively. The nominal trajectory is the
solution of (47) with normalized initial condition

$\displaystyle r_{0}$ $\displaystyle=\frac{R_{E}+h}{R_{E}}=1.0533,$ (48a)
$\displaystyle\dot{r}_{0}$ $\displaystyle=0,$ (48b) $\displaystyle\theta_{0}$
$\displaystyle=0,$ (48c) $\displaystyle\dot{\theta}_{0}$
$\displaystyle=\left(\frac{V_{\theta}T_{p}}{R_{E}}\right)\frac{1}{r_{0}}=6.2832.$
(48d)

The parameters necessary to simulate the system are provided in table I.

$R_{E}=6378.1363$ km $\mu_{E}=398600.4415$ km3/s2 $T_{p}=5.48\times 10^{3}$ s
$J_{2}=1.7555\times 10^{10}$ km5/s2 $V_{\theta}=7.7027$ km/s $h$ = 340 km

TABLE I: Parameters in the satellite dynamics model.

Equation (47) is linearized about the nominal trajectory to obtain a
continuous-time periodic system. We augment the linear model with process
noise, to account for the effects of sporadic thrusts that are necessary for
orbital station keeping. The augmented model is given by,

$\displaystyle\dot{\boldsymbol{x}}=\boldsymbol{A}(t)\boldsymbol{x}+\boldsymbol{B}\boldsymbol{w}(t),$
(49)

where
$\boldsymbol{x}:=\begin{bmatrix}r&\dot{r}&\theta&\dot{\theta}\end{bmatrix}^{T}$
is the state vector,
$\boldsymbol{w}(t):=\begin{bmatrix}w_{r}(t)&w_{\theta}(t)\end{bmatrix}^{T}$ is
a zero-mean Gaussian random process,

$\displaystyle\boldsymbol{A}(t):=\begin{bmatrix}0&1.0&0&0\\\
a_{21}(t)&0&a_{23}(t)&12.59\\\ 0&0&0&1.0\\\
a_{41}(t)&-12.21&a_{43}(t)&0\end{bmatrix},$
$\displaystyle\boldsymbol{B}:=\begin{bmatrix}0&0\\\ 1&0\\\ 0&0\\\
0&1\end{bmatrix},$ (50a)

with

$\displaystyle a_{21}(t)$ $\displaystyle=0.416\cos(12.4t)+126.4,$
$\displaystyle a_{23}(t)$ $\displaystyle=0.2113\sin(12.4t),$ $\displaystyle
a_{41}(t)$ $\displaystyle=0.2774\sin(12.4t),$ $\displaystyle a_{43}(t)$
$\displaystyle=-0.1408\cos(12.4t).$

In this example, we assume the mass of the satellite is $100$ kg, and the
satellite sporadically applies maximum of 1mN of thrust for orbital station
keeping. The normalized accelerations due to these thrusts are modeled as
zero-mean Gaussian random processes $w_{r}(t)$ and $w_{\theta}(t)$, with
$\mathbb{E}\left[w_{r}w_{r}^{T}\right]=\mathbb{E}\left[w_{\theta}w^{T}_{\theta}\right]=0.0471^{2}$.
Fig.(5) shows the propagation of mean $\boldsymbol{\mu}(t)$ and variance
$\boldsymbol{\Sigma}(t)$ for the time-varying linear system, with

$\displaystyle\boldsymbol{\mu}(t_{0})$
$\displaystyle:=\begin{bmatrix}50/R_{E}\\\ 0\\\ 0\\\ 0\end{bmatrix},\text{ and
}$ (51a) $\displaystyle\boldsymbol{\Sigma}(t_{0})$
$\displaystyle:=0.01\times\mathbf{diag}\left(\boldsymbol{\mu}(t_{0})\right).$
(51b)

The evolution equation for $\boldsymbol{\mu}(t)$ and $\boldsymbol{\Sigma}(t)$
are given by

$\displaystyle\dot{\boldsymbol{\mu}}(t)$
$\displaystyle=\boldsymbol{A}(t)\boldsymbol{\mu}(t),$ (52)
$\displaystyle\dot{\boldsymbol{\Sigma}}(t)$
$\displaystyle=\boldsymbol{A}(t)\boldsymbol{\Sigma}(t)+\boldsymbol{\Sigma}(t)\boldsymbol{A}^{T}(t)+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T},$
(53)

where $\boldsymbol{Q}:=0.0471^{2}\times\boldsymbol{I}_{2}$.

Figure 5: Uncertainty propagation with the linear time-varying periodic
system. Solid line shows the evolution of the mean perturbation, and the
shaded region shows $\boldsymbol{\mu}_{i}\pm\sqrt{\boldsymbol{\Sigma}_{ii}}$
for $i=1,2,3$ and $4$.

We discretize the normalized time interval $[0,1]$ with $dt=0.1$, resulting in
the temporal grid $\\{t_{k}\\}$, where $t_{k}:=kdt$. We assume that
measurements are available at these times. The dynamics in (49) is discretized
over $\\{t_{k}\\}$, and is given by

$\displaystyle\boldsymbol{x}_{k+1}=\boldsymbol{A}_{d}(t_{k})\boldsymbol{x}_{k}+\boldsymbol{w}_{k},$
(54)

where

$\displaystyle\boldsymbol{A}_{d}(t_{k})$
$\displaystyle:=\boldsymbol{\Phi}(t_{k+1},t_{k}),$ (55)
$\displaystyle\boldsymbol{w}_{k}$
$\displaystyle:=\int_{t_{k}}^{t_{k+1}}\boldsymbol{\Phi}(\tau,t_{k})\boldsymbol{B}\boldsymbol{w}(\tau)d\tau,$
(56)

and $\boldsymbol{\Phi}(\cdot,\cdot)$ is the state-transition matrix, which is
obtained by numerical integration of the fundamental matrix.

It is easy to verify that if $\mathbb{E}\left[\boldsymbol{w}(t)\right]=0$,
then $\mathbb{E}\left[\boldsymbol{w}_{k}\right]=0$. Therefore,
$\boldsymbol{w}_{k}$ is a zero-mean random process. To determine the optimal
precision, we need to quantify
$\boldsymbol{Q}_{k}:=\mathbb{E}\left[\boldsymbol{w}_{k}\boldsymbol{w}^{T}_{k}\right]$,
which is difficult to determine from (56). Instead, we use the covariance
$\boldsymbol{\Sigma}(t)$, determined by solving (53), and the discrete-time
covariance propagation equation, to determine the time-varying
$\boldsymbol{Q}_{k}$. It is given by

$\displaystyle\boldsymbol{Q}_{k}:=\boldsymbol{\Sigma}(t_{k+1})-\boldsymbol{A}_{d}(t_{k})\boldsymbol{\Sigma}(t_{k})\boldsymbol{A}^{T}_{d}(t_{k}).$
(57)

####  V-B2 Example 3: Optimal Sensor Scheduling

In this example, we consider a set of 10 laser-ranging sensors located on the
surface of the Earth, at angular positions $\theta(t_{k})$. For the periodic
system described above, the objective is to determine the optimal sensor
precisions such that
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}(t_{k}=1)\right)\leq\gamma_{d}$,
given $\boldsymbol{\Sigma}^{-}(t_{k}=0)$. Here we apply theorem 1 to determine
the optimal sensor precisions, which are shown in fig.(6) for various values
of $\boldsymbol{s}_{\text{max}}$.

The optimization is done with
$\boldsymbol{\Sigma}(t_{k}=0)=\boldsymbol{\Sigma}(t_{0})$, and
$\gamma_{d}=0.1\times\mathbf{tr}\left(\boldsymbol{\Sigma}^{-}(t_{k}=1)\right)$,
where $\boldsymbol{\Sigma}^{-}(t_{k}=1)$ is the prior obtained at $t_{k}=1$.
It is obtained by propagating $\boldsymbol{\Sigma}(t_{k}=0)$ using (53).
Variance $\boldsymbol{\Sigma}(t_{0})$ is defined in (51b). The value of
$\gamma_{d}$ specifies that we want the trace of the posterior to be 10% of
the prior at $t_{k}=1$.

Figure 6: Optimal precisions for
$\boldsymbol{s}_{\text{max}}=2500,1500,819.60.$

From fig.(6), we see that there is a tradeoff between sensor precision and
sensing frequency. With low precision sensors
($\boldsymbol{s}_{\text{max}}=819.60$), we need to sense the satellite at all
the sites to guarantee
$\mathbf{tr}\left(\boldsymbol{\Sigma}^{+}(t_{k}=1)\right)<\gamma_{d}$. As
$\boldsymbol{s}_{\text{max}}$ is increased, the sensing becomes more sparse.
With $\boldsymbol{s}_{\text{max}}=2500$, we only need to get two range
measurements, in order to estimate the state at $t_{k}=1$ with the required
accuracy.

This is significant for large-scale spatio-temporal sensing, and especially
for tracking space objects. Currently, there are about $500K$ space objects,
and only $30K$ are tracked. The number of sensing sites is significantly lower
than that. Therefore, by applying theorem 1, it is possible to determine the
sparse sensing schedule for each object, thereby increasing the ability to
track more objects with guaranteed accuracy.

Precision in laser-ranging sensors is determined by the energy of the laser
beam and the reflectivity of the object being sensed. A lower value of
precision implies lower energy requirements for sensing the object, which
results in optimal sensor designs. Modern satellites have reflectors, which
allows precise sensing with low powered lasers. Some space objects (e.g.,
asteroids, etc.) have poor reflectivity, which fundamentally limits how
accurately it can be sensed. This limit can be incorporated using the variable
$\boldsymbol{s}_{\text{max}}$.

##  VI Summary & Conclusion

This paper presented convex optimization problem-formulations to determine
optimal sensor precisions that guarantee a specified estimation accuracy. The
formulation is presented in a general multi-rate sensing framework, with
linear time-varying discrete-time system dynamics. Optimality is achieved by
minimizing sensor precisions, subject to the upper bound on the estimation
error, as defined in the discrete-time Kalman filtering framework. Since the
minimization of precisions is done with respect to the $l_{1}$ norm, the
proposed optimization framework can also be used to determine sparse sensing
architectures. This will be valuable in the design of large-scale sensor
networks. We have shown the proposed theory’s engineering value by applying it
to realistic flight mechanics and astrodynamics problems.

##  VII Acknowledgements

This research has been supported by the National Science Foundation grant
#1762825.

## References

  * [1] Carlos Guestrin, Andreas Krause, and Ajit Paul Singh.  Near-optimal sensor placements in gaussian processes.  In Proceedings of the 22nd International Conference on Machine learning. ACM Press, 2005.
  * [2] Vijay Gupta, Timothy H. Chung, Babak Hassibi, and Richard M. Murray.  On a stochastic sensor selection algorithm with applications in sensor scheduling and sensor coverage.  Automatica, 42(2):251–260, feb 2006.
  * [3] Lei Zhang and Dimitrios Hristu-Varsakelis.  Communication and control co-design for networked control systems.  Automatica, 42(6):953–958, jun 2006.
  * [4] Ying He and Edwin K.P. Chong.  Sensor scheduling for target tracking: A monte carlo sampling approach.  Digital Signal Processing, 16(5):533–545, sep 2006.
  * [5] Dawei Shi and Tongwen Chen.  Optimal periodic scheduling of sensor networks: A branch and bound approach.  Systems & Control Letters, 62(9):732–738, sep 2013.
  * [6] Sergio Pequito, Soummya Kar, and A. Pedro Aguiar.  A structured systems approach for optimal actuator-sensor placement in linear time-invariant systems.  In 2013 American Control Conference. IEEE, jun 2013.
  * [7] Syed Talha Jawaid and Stephen L. Smith.  Submodularity and greedy algorithms in sensor scheduling for linear dynamical systems.  Automatica, 61:282–288, nov 2015.
  * [8] Tyler H. Summers, Fabrizio L. Cortesi, and John Lygeros.  On submodularity and controllability in complex dynamical networks.  IEEE Transactions on Control of Network Systems, 3(1):91–101, mar 2016.
  * [9] Duo Han, Junfeng Wu, Huanshui Zhang, and Ling Shi.  Optimal sensor scheduling for multiple linear dynamical systems.  Automatica, 75:260–270, jan 2017.
  * [10] Xudong Chen, Mohamed-Ali Belabbas, and Tamer Başar.  Optimal capacity allocation for sampled networked systems.  Automatica, 85:100–112, nov 2017.
  * [11] S. Joshi and S. Boyd.  Sensor selection via convex optimization.  IEEE Transactions on Signal Processing, 57(2):451–462, feb 2009\.
  * [12] Neil K. Dhingra, Mihailo R. Jovanovic, and Zhi-Quan Luo.  An ADMM algorithm for optimal sensor and actuator selection.  In 53rd IEEE Conference on Decision and Control. IEEE, dec 2014\.
  * [13] Sundeep Prabhakar Chepuri and Geert Leus.  Sparsity-promoting sensor selection for non-linear measurement models.  IEEE Transactions on Signal Processing, 63(3):684–698, feb 2015\.
  * [14] V. Tzoumas, A. Jadbabaie, and G. J. Pappas.  Sensor placement for optimal kalman filtering: Fundamental limits, submodularity, and algorithms.  In 2016 American Control Conference (ACC). IEEE, jul 2016.
  * [15] Haotian Zhang, Raid Ayoub, and Shreyas Sundaram.  Sensor selection for kalman filtering of linear dynamical systems: Complexity, limitations and greedy algorithms.  Automatica, 78:202–210, 2017.
  * [16] Faming LI, Maurcio C. de Oliveira, and Robert E. Skelton.  Integrating information architecture and control or estimation design.  SICE Journal of Control, Measurement, and System Integration, 1(2):120–128, 2008.
  * [17] Radhika Saraf, Raktim Bhattacharya, and Robert Skelton.  $\mathcal{H}_{2}$ optimal sensing architecture with model uncertainty.  In 2017 American Control Conference (ACC), pages 2429–2434. IEEE, 2017.
  * [18] Brian DO Anderson and John B Moore.  Optimal filtering.  Courier Corporation, 2012.
  * [19] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learning: data mining, inference, and prediction.  Springer Science & Business Media, 2009.
  * [20] Wolfgang Koch and Felix Govaers.  Accumulated state densities and their use in decorrelated track-to-track fusion.  In 2015 International Conference on Military Communications and Information Systems (ICMCIS), pages 1–9. IEEE, 2015.
  * [21] Stephen Boyd, Laurent El Ghaoui, Eric Feron, and Venkataramanan Balakrishnan.  Linear matrix inequalities in system and control theory.  SIAM, 1994.
  * [22] Bhaskar D Rao and Kenneth Kreutz-Delgado.  An affine scaling methodology for best basis selection.  IEEE Transactions on signal processing, 47(1):187–200, 1999.
  * [23] Emmanuel J Candes, Michael B Wakin, and Stephen P Boyd.  Enhancing sparsity by reweighted $l_{1}$ minimization.  Journal of Fourier analysis and applications, 14(5-6):877–905, 2008\.
  * [24] Ameet S. Deshpande.  Bridging a gap in applied kalman filtering: Estimating outputs when measurements are correlated with the process noise [focus on education].  IEEE Control Systems, 37(3):87–93, jun 2017.
  * [25] Robert R Bitmead, Michel R Gevers, Ian R Petersen, and R John Kaye.  Monotonicity and stabilizability-properties of solutions of the riccati difference equation: Propositions, lemmas, theorems, fallacious conjectures and counterexamples.  Systems & Control Letters, 5(5):309–315, 1985.
  * [26] Ryan James Caverly and James Richard Forbes.  Lmi properties and applications in systems, stability, and control theory, 2021.
  * [27] EC Warner and JT Scruggs.  Iterative convex overbounding algorithms for bmi optimization problems.  IFAC-PapersOnLine, 50(1):10449–10455, 2017.
  * [28] Fu Lin, Makan Fardad, and Mihailo R. Jovanovic.  Design of optimal sparse feedback gains via the alternating direction method of multipliers.  IEEE Transactions on Automatic Control, 58(9):2426–2431, sep 2013\.
  * [29] Gene F Franklin, J David Powell, Michael L Workman, et al.  Digital control of dynamic systems, volume 3.  Addison-wesley Reading, MA, 1998.
  * [30] John L Junkins and Hanspeter Schaub.  Analytical mechanics of space systems.  American Institute of Aeronautics and Astronautics, 2009.