# Towards Understanding Sensor and Control Nodes Selection in Nonlinear
Dynamic Systems: Lyapunov Theory Meets Branch-and-Bound

Sebastian A. Nugroho<EMAIL_ADDRESS>Ahmad F. Taha
<EMAIL_ADDRESS>Department of Electrical and Computer Engineering, The
University of Texas at San Antonio, 1 UTSA Circle, San Antonio, TX 78249

###### Abstract

Sensor and actuator selection problems (SASPs) are some of the core problems
in dynamic systems design and control. These problems correspond to
determining the optimal selection of sensors (measurements) or actuators
(control nodes) such that certain estimation/control objectives can be
achieved. While the literature on SASPs are indeed inveterate, the vast
majority of the work focuses on linear(ized) representation of the network
dynamics, resulting in the placements of sensors or actuators (SAs) that are
valid for confined operating regions. As an alternative, herein we propose a
new general framework for addressing SASPs in nonlinear dynamic systems
(NDSs), assuming that the inputs and outputs are linearly coupled with the
nonlinear dynamics. This is investigated through (i) classifying and
parameterizing the NDSs into various nonlinear function sets, (ii) utilizing
rich Lyapunov theoretic formulations, and (iii) designing a new customized
branch-and-bound (BnB) algorithm that exploits problem structure of the SASPs.
The newly designed BnB routines are computationally more attractive than the
standard one and also directly applicable to solve SASPs for linear systems.
In contrast with contemporary approaches from the literature, our approach is
suitable for finding the optimal SAs combination for stable/unstable NDSs that
ensures stabilization of estimation error and closed-loop dynamics through a
simple linear feedback control policy.

###### keywords:

Sensor selection , actuator selection , nonlinear dynamic systems , Lipschitz
continuous , mixed-integer semidefinite programming , branch-and-bound
algorithm

††journal: Automatica

##  1 Introduction and Paper Contributions

Sensor and actuator selection problems (SASPs)—which fall into (i) finding
optimal geographic placements or time-varying selection of sensors or
actuators (SAs) (but not necessarily simultaneously), while (ii) optimizing
state estimation or control metrics—have become prevalent research topics in
numerous science and engineering fields. For instance, in power networks
application, sensor selection corresponds to the placement of phasor
measurement units for the purpose of power systems monitoring [1] while
actuator selection corresponds to the placement of energy storage devices to
ensure system controllability [2]. Both problems are crucial to achieve a more
economical and robust power systems operations while ensuring that certain
minimal requirements, such as allowable tolerance on dynamic estimation error
and voltage/frequency deviation, are met.

Extensive studies have been carried out in the literature to address SASPs for
particularly linear dynamic systems—see [3, 4, 5, 6] for some notable
references. Unlike SASPs in linear(ized) systems which yield feasible or
optimal SAs configurations for specific operating points, SASPs for nonlinear
dynamic systems (NDSs) offer SA selections that are applicable for a much
larger operating regions or even globally—this is demonstrated in our prior
work through a simple example [7]. Indeed, SASPs for NDSs have not received
comparable attention in the literature as the SASPs for linear systems, until
recently. In [1], the concept of empirical Gramian to quantify the
observability of NDSs is considered jointly with the sensor selection problem
(SSP) to determine the combination of sensors that maximizes the logarithmic
determinant of the empirical Gramian matrix. The joint problem of sensor
selection and state observation is investigated in [8]. In particular, a new
method for reconstructing the initial states of NDSs while simultaneously
selecting sensors for a given observation window is presented. A different
approach is pursued in [9] where the authors introduce a randomized algorithm
for dealing with SSP and develop theoretical bounds for eigenvalue and
condition number of observability Gramian. Recently, a method for selecting
the control nodes and designing control actions for NDSs, which is based on an
open-loop predictive control framework, is proposed in [10]. In spite of these
efforts, their potential applicability to address SASPs in unstable NDSs along
with their ability to incorporate SAs selection with some estimation/control-
theoretic metrics, unfortunately, remains unclear.

The objective of this paper is to find the minimal SAs combination for
stable/unstable NDSs that ensures stabilization for estimation error (via
minimal sensor selection) and closed loop system through feedback control (via
minimal actuator selection). The proposed algorithms developed herein are
based on the following observations. First, numerous observer and controller
designs for NDSs have been developed in the past two decades. By using
Lyapunov stability theory and given that the nonlinearities in NDSs satisfy
some function sets such as bounded Jacobian [12], Lipschitz continuous [13],
one-sided Lipschitz [14], quadratically inner-bounded [14], and quadratically
bounded [15], the procedure to compute stabilizing observer/controller gain
matrix can be posed as convex semidefinite programmings (SDPs). Secondly,
physical states in many practical NDS models are almost always bounded. Hence,
it can be shown or safely assumed that the corresponding nonlinearities
satisfy the properties of the aforementioned function sets in confined state-
space regions—regions that are much larger than operating regions of
linearization points.

Without loss of generality, this paper focuses on SASPs for Lipschitz NDSs in
which the vector-valued nonlinearity in the system dynamics satisfies the
Lipschitz continuity assumption. The proposed method can be easily extended
for other types of non-Lipschitz nonlinearities. It is assumed throughout the
paper that the inputs and outputs of the NDSs form linear relations with the
nonlinear dynamics. It is also worth noting that we do not address the
simultaneous selection of SAs. That is, we only consider the minimal selection
of sensors with observer design or actuators with controller design. A
preliminary version of this work appeared in [16], where the SSP for Lipschitz
NDSs is studied and solved via a general purpose solver [17]. The key
contributions of this paper are summarized as follows:

  * •

We present a novel generalized framework to address SASPs in NDSs by: (i)
parameterizing NDSs’ nonlinearities based on the corresponding function sets,
(ii) formulating the SASPs using a plethora of observer and/or controller
designs, and (iii) reformulating the resulting SASPs into convex mixed-integer
SDPs (MISDPs).

  * •

Particularly for Lipschitz NDSs with Lipschitz constant $\gamma_{l}$, we
formalize NDSs observability and controllability (coined as
$\gamma_{l}$-observability and $\gamma_{l}$-controllability, respectively) as
SDP feasibility problems and show that it is crucial to obtain the smallest
Lipschitz constant in order to obtain less number of SAs in SASPs. In
particular, we theoretically investigate the relationship between the
Lipschitz constant of the system and the number of activated SAs.

  * •

A new customized branch-and-bound (BnB) algorithm for MISDPs to solve SASPs is
proposed. This BnB algorithm, referred to as structure-exploiting BnB (SE-
BnB), utilizes heuristics and exploits problem structure to efficiently find
optimal and suboptimal solutions for SASPs, owing to the reduced number of
constraints and/or variables of the corresponding SDP problems when finding
upper an lower bounds on each node of the BnB tree.

  * •

We showcase the computational advantages of the proposed BnB routines for
solving SSP in comparison with our implementation of the standard BnB
algorithm and the $\ell_{1}$-norm relaxation technique [18].

The distinction of our approach to address SASPs for NDSs in comparison with
the other contemporary methods includes (a) its ability to find optimal
placement of SAs with stabilizing estimation/control gain matrix, and (b) its
flexibility to incorporate estimation/control metrics for continuous-
time/discrete-time NDSs having different classes of nonlinearity. The paper’s
organization and notation are as follows. Section 2 formalizes the NDSs model.
Section 3 investigates the relation between Lipschitz constant and
$\gamma_{l}$-observability while Section 4 constructs the SSP for Lipschitz
NDSs. Next, Section 5 focuses on transforming SSP into a convex MISDP form and
provides a detailed description for the SE-BnB algorithm. In Section 6, some
discussions related to actuator selection problem, solving SASPs beyond
Lipschitz NDSs, and robust SSP are provided. Finally, Section 7 presents some
numerical results and Section 8 concludes the paper.

Notation. The notation $\bm{1}$ represents a matrix of appropriate dimension
with elements of 1. The notations $\mathbb{R}^{n}$ and $\mathbb{R}^{p\times
q}$ denote the sets of row vectors with $n$ elements and matrices with size
$p$-by-$q$ with elements in $\mathbb{R}$. Given $\bm{A}$ as a matrix, its
$i$-th and $j$-th element is denoted by $A_{ij}$. The operators
$\mathrm{Blkdiag}(\cdot)$ constructs a block diagonal matrix,
$\mathrm{Diag}(\cdot)$ constructs a diagonal matrix from a vector,
$\mathrm{vec}(\cdot)$ constructs a vector by stacking each column of a matrix,
$\otimes$ denotes the Kronecker product, and $\odot$ denotes the Hadamard
product. The symbol $*$ is used to represent symmetric entries in symmetric
matrices. The set $\mathbb{I}(n)$ is defined as
$\mathbb{I}(n):=\\{i\in\mathbb{N}\,|\,1\leq i\leq n\\}$, which is usually used
to represent the set of indices.

Table 1: Dynamic models of some nonlinear dynamic systems in electric power systems, traffic, combustion networks, …. Type of Systems | Nonlinear Dynamics Model | Description
---|---|---
Electric power grids[19] |  $\dot{\delta_{i}}=\omega_{i}-\omega_{0}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\;$ $\dot{\omega}_{i}=\frac{\omega_{0}}{2H_{i}}\bigg{(}P_{\textrm{m}i}-\sum\limits_{j=1}^{N}E_{i}E_{j}\left(\bar{G}_{ij}\cos(\delta_{i}-\delta_{j})+\bar{B}_{ij}\sin(\delta_{i}-\delta_{j})\right)-\frac{K_{\textrm{D}i}}{\omega_{0}}(\omega_{i}-\omega_{0})\bigg{)}$ |  2-order swing equation $\bm{x_{i}}:=[\delta_{i}\;\;\omega_{i}]^{\top}$ $\delta_{i}$: rotor angle, $\omega_{i}$: rotor speed
Highway traffic [20] | $\dot{\rho}_{i}=\frac{v_{f}}{l}\left(\rho_{i-1}-\rho_{i}+\hat{\rho}_{j}-\alpha(k)\check{\rho}_{k}\right)-\delta\left(\rho_{i-1}^{2}-\rho_{i}^{2}+\hat{\rho}_{j}^{2}-\alpha(k)\check{\rho}_{k}^{2}\right)$ |  Free-flow condition $\bm{x_{i}}:=[\rho_{i}]$, $\rho_{i}$: traffic density
Combustion networks[21] | $\dot{\chi}_{i}=\sum\limits_{j=1}^{n_{r}}(\beta_{ji}-\alpha_{ji})\Bigg{(}d^{(f)}_{j}\prod\limits^{n}_{k=1}\chi_{k}^{\alpha_{jk}}-d^{(b)}_{j}\prod\limits^{n}_{k=1}\chi_{k}^{\beta_{jk}}\Bigg{)}$ |  $\bm{x_{i}}:=[\chi_{i}]$ $\chi_{i}$: chemical concentration
Oscillator synchronization [22] | $\dot{\theta}_{i}=\omega_{i}+\sum\limits_{j=1}^{N}K_{ij}\sin(\theta_{j}-\theta_{i})$ |  Kuramoto model $\bm{x_{i}}:=[\theta_{i}]$, $\theta_{i}$: oscillator phase
Epidemic outbreaks[23] | $\dot{p}_{i}=-\delta p_{i}+\sum\limits_{j=1}^{N}a_{ij}\beta p_{j}(1-p_{i})$ |  $\bm{x_{i}}:=[p_{i}]$, $p_{i}\in[0,1]$ $p_{i}$: infection probability
Mass-spring-damper systems |  $m_{i}\ddot{x}_{i}=d_{i}\dot{x}_{i-1}-(d_{i}+d_{i+1})\dot{x}_{i}+d_{i+1}\dot{x}_{i+1}+k_{i}{x}_{i-1}-(k_{i}+k_{i+1}){x}_{i}+k_{i+1}{x}_{i+1}$ $+a_{i}{x}_{i-1}^{3}-(a_{i}+a_{i+1}){x}_{i}^{3}+a_{i+1}{x}_{i+1}^{3}\qquad\qquad\qquad\qquad\qquad\qquad\quad$ |  $\bm{x_{i}}:=[x_{i}]$ $x_{i}$: mass relative position

##  2 System Description and Preliminaries

Consider a continuous-time networked NDS comprised of $N$ subsystems (also
referred to as nodes) modeled as

$\displaystyle\dot{\bm{x}}(t)$
$\displaystyle=\bm{A}\bm{x}(t)+\bm{G}\bm{f}(\bm{x})+\bm{B}\bm{u}(t),\quad\bm{y}(t)=\bm{C}\bm{x}(t).$
(1)

Numerous NDSs—see Table 1 for some notable examples—can be represented in the
form of (1). In this model, the global state
$\bm{x}\in\mathbfcal{X}\subset\mathbb{R}^{n_{x}}$ consists of $N$ numbers of
nodal state
$\bm{x_{i}}\in\mathbfcal{X}_{\bm{i}}\subset\mathbb{R}^{n_{x_{i}}}$. Likewise,
each subsystem $\bm{i}$ is also comprised of $n_{u_{i}}$inputs and $n_{y_{i}}$
measurements such that
$\bm{u_{i}}\in\mathbfcal{U}_{\bm{i}}\subset\mathbb{R}^{n_{u_{i}}}$ and
$\bm{y_{i}}\in\mathbb{R}^{n_{y_{i}}}$. The notations $\bm{u_{i}}$ and
$\bm{y_{i}}$ denote the nodal input and measurements vectors whereas
$\mathbfcal{X}$ and $\mathbfcal{U}$ represent the operating region and
admissible inputs for NDS (1). The global input and output vectors are
$\bm{u}\in\mathbfcal{U}\subset\mathbb{R}^{n_{u}}$ and
$\bm{y}\in\mathbb{R}^{n_{y}}$. Matrices $\bm{A}$, $\bm{B}$, $\bm{C}$, and
$\bm{G}$ are all assumed to have appropriate dimensions and in particular,
$\bm{B}$ and $\bm{C}$ are representative of input-to-state and state-to-output
mappings and supposed to posses the following structure:
$\bm{B}=\mathrm{Blkdiag}\left(\bm{B}_{1},\ldots,\bm{B}_{N}\right)$,
$\bm{C}=\mathrm{Blkdiag}\left(\bm{C}_{1},\ldots,\bm{C}_{N}\right)$. Matrix
$\bm{A}$ expresses the linear dynamics of the system as well as interaction
between subsystems while matrix $\bm{G}$ depicts the distribution of nonlinear
mapping $\bm{f}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{g}}$ that may
capture any linear and nonlinear phenomena in the system.

The ultimate objective of this work is to search for both optimal and
suboptimal SAs configurations or locations for NDSs (or equivalently, the
number of nonzero columns of $\bm{B}$ and nonzero rows of $\bm{C}$) while
satisfying some user-defined constraints, which generally may include closed-
loop system stability, convergence of estimation error dynamics, and
actuator/sensor constraints. It is worth mentioning that this paper emphasizes
solving SASPs for continuous-time NDSs with Lipschitz nonlinearities since (i)
Lipschitz is one of widely used function sets both in observer and controller
designs—for example, see [13, 20]—and (ii) the proposed methodology to solve
SASPs can be directly extended for other classes of nonlinearities mentioned
previously as well as for discrete-time NDSs with additional
estimation/control objectives To proceed, the subsequent assumption is
considered in the paper.

######  Assumption 1.

The mapping $\bm{f}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{g}}$ is
locally Lipschitz continuous in $\mathbfcal{X}$ such that for any
$\bm{x},\hat{\bm{x}}\in\mathbfcal{X}$

$\displaystyle\left\lVert\bm{f}(\bm{x})-\bm{f}(\hat{\bm{x}})\right\rVert_{2}\leq\gamma_{l}\left\lVert\bm{x}-\hat{\bm{x}}\right\rVert_{2},$
(2)

where $\gamma_{l}\geq 0$ is the corresponding Lipschitz constant.

Readers are referred to [24] for scalable numerical methods to compute
Lipschitz constants—including the ones corresponding to the other function
sets. In the next section, we formalize the notion of observability for
Lipschitz NDSs which is crucial in the context of SSP.

##  3 Observability for Lipschitz NDS

In contrast to linear systems, quantifying observability for nonlinear systems
is, without doubt, much more difficult and less straightforward. To that end,
in this section we focus on the concept of observability for NDS (1) in regard
to solvability of the corresponding SDPs and Lipschitz constant. To proceed,
consider a NDS in the form of (1). A standard Luenberger-like observer for NDS
(1) can then be constructed as follows

$\displaystyle\dot{\hat{\bm{x}}}(t)$
$\displaystyle=\bm{A}\hat{\bm{x}}(t)+\bm{G}\bm{f}(\hat{\bm{x}})+\bm{B}\bm{u}(t)+\bm{L}(\bm{y}(t)-\hat{\bm{y}}(t))$
(3a) $\displaystyle\hat{\bm{y}}(t)$ $\displaystyle=\bm{C}\hat{\bm{x}}(t),$
(3b)

where in (3), $\bm{L}\in\mathbb{R}^{n_{x}\times n_{y}}$ is the observer gain
matrix, $\hat{\bm{x}}$ is the estimated states, and $\hat{\bm{y}}$ is the
estimated outputs. The objective here is to compute $\bm{L}$ that makes the
estimation error dynamics, where estimation error is constructed as
$\bm{e}(t):=\bm{x}(t)-\hat{\bm{x}}(t)$, to converge asymptotically towards
zero. The next definition quantifies observability for NDS (1), posed as a SDP
feasibility problem, for a given $\gamma_{l}$.

######  Definition 1.

NDS (1) is said to be $\gamma_{l}$-observable if and only if there exist
$\bm{P}\in\mathbb{S}^{n_{x}}_{++}$, $\bm{Y}\in\mathbb{R}^{n_{x}\times n_{y}}$,
and $\epsilon\in\mathbb{R}_{++}$ such that the following linear matrix
inequality (LMI) is feasible

$\displaystyle\begin{bmatrix}\bm{A}^{\top}\bm{P}+\bm{P}\bm{A}-\bm{C}^{\top}\bm{Y}^{\top}-\bm{Y}\bm{C}+\epsilon\gamma_{l}^{2}\bm{I}&*\\\
\bm{G}^{\top}\bm{P}&-\epsilon\bm{I}\end{bmatrix}$ $\displaystyle\prec 0.$ (4)

LMI (4) originates from [13] with special case $\bm{G}=\bm{I}$ and it can be
proven that (4) is indeed sufficient and necessary for the existence of
stabilizing matrix $\bm{L}$ for the estimation error dynamics, provided that
$V(\bm{e})=\bm{e}(t)^{\top}\bm{P}\bm{e}(t)$ is the Lyapunov function candidate
and $\epsilon\in\mathbb{R}_{+}$. Once the LMI is solved, $\bm{L}$ can be
recovered as $\bm{P}^{-1}\bm{Y}$. The advantages of utilizing Definition 1 to
quantify observability are twofold. First, if there exists any solution for
(4) then one can immediately obtain stabilizing gain matrix $\bm{L}$ and
second, the observability prevails for various operating points in
$\mathbfcal{X}$. These are in contrast to other types of observability for NDS
such as empirical observability Gramian, as it only quantifies observability
based on local behaviors of the systems around certain operating points and,
upon determining the observability, one is still required to design a state
estimator algorithm to estimate the actual states.

At this point, we are now compelled to consider the following question
regarding the feasibility of SDP that corresponds to
$\gamma_{l}$-observability for NDS (1):

Q1: What role does constant $\gamma_{l}$ have on $\gamma_{l}$-observability?

To answer Q1, we need to further study the LMI given in (4). For the sake of
simplicity, let us assume that $\epsilon>0$ is a constant. Then, realize that
$\bm{Y}\bm{C}$ can be expressed as

$\displaystyle\bm{Y}\bm{C}=\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}Y_{ij}\check{\bm{C}}_{ij},$

where $\check{\bm{C}}_{ij}\in\mathbb{R}^{n_{y}\times n_{x}}$ is constructed as

$\displaystyle\mathrm{row}\left(\check{\bm{C}}_{ij}\right)_{k}=\begin{cases}\left\\{[C_{kl}]\right\\}^{n_{x}}_{l=1}\in\mathbb{R}^{1\times
n_{x}},&\;\text{if}\;j=k\\\ \bm{0}_{1\times n_{x}},&\;\text{otherwise}.\\\
\end{cases}$

Now let
$\check{\bm{C}}_{S,ij}:=\check{\bm{C}}_{ij}+\check{\bm{C}}_{ij}^{\top}$ for
each $i,j$. By using the above expression, the term
$\bm{Y}\bm{C}+\bm{C}^{\top}\bm{Y}^{\top}$ can be expressed as

$\displaystyle\bm{Y}\bm{C}+\bm{C}^{\top}\bm{Y}^{\top}$
$\displaystyle=\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}Y_{ij}\check{\bm{C}}_{ij}+\left(\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}Y_{ij}\check{\bm{C}}_{ij}\right)^{\top}$
$\displaystyle=\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}Y_{ij}\check{\bm{C}}_{S,ij}=\sum_{k=1}^{n_{x}n_{y}}y_{v,k}\bar{\bm{C}}_{k},$
(5)

in which
$\bm{y}_{v}:=[y_{v,1}\,\,y_{v,2}\,\,\cdots\,\,y_{v,n_{x}n_{y}}]^{\top}=\mathrm{vec}(\bm{Y})$
and $\bar{\bm{C}}_{k}\in\mathbb{S}^{n_{x}}$ for each
$k\in\mathbb{I}(n_{x}n_{y})$ is associated with the corresponding
$\check{\bm{C}}_{S,ij}$. Using (5), LMI (4) along with $\bm{P}\succ 0$ can be
written in the following compact form

$\displaystyle\mathcal{A}(\bm{P},\bm{y}_{v})+\mathcal{A}_{0}\succ 0,$ (6a) in
which
$\mathcal{A}:\mathbb{S}^{n_{x}}\times\mathbb{R}^{n_{x}n_{y}}\rightarrow\mathbb{S}^{2n_{x}}\times\mathbb{S}^{n_{x}}$
is constructed as
$\mathcal{A}(\bm{P},\bm{y}_{v}):=\mathcal{A}_{P}(\bm{P})+\mathcal{A}_{Y}({\bm{y}_{v}})$
and $\displaystyle\mathcal{A}_{P}(\bm{P})$
$\displaystyle:=\mathrm{Blkdiag}\left(-\begin{bmatrix}\bm{A}^{\top}\bm{P}+\bm{P}\bm{A}&*\\\
\bm{G}^{\top}\bm{P}&\bm{O}\end{bmatrix},\bm{P}\right)$ (6b)
$\displaystyle\mathcal{A}_{Y}({\bm{y}_{v}})$
$\displaystyle:=\mathrm{Blkdiag}\left(\begin{bmatrix}\sum_{i=1}^{n_{x}n_{y}}y_{v,i}\bar{\bm{C}}_{i}&*\\\
\bm{O}&\bm{O}\end{bmatrix},\bm{O}\right)$ (6c) $\displaystyle\mathcal{A}_{0}$
$\displaystyle:=\mathrm{Blkdiag}\left(\begin{bmatrix}-\epsilon\gamma_{l}^{2}\bm{I}&*\\\
\bm{O}&\epsilon\bm{I}\end{bmatrix},\bm{O}\right).$ (6d)

Using the expression provided in (6), we now present the first theoretical
result of the paper.

######  Theorem 1 (Larger Lipschitz Constant Reduces Observability).

There exist $\bm{P}\in\mathbb{S}^{n_{x}}$ and
$\bm{y}_{v}\in\mathbb{R}^{n_{x}n_{y}}$ such that
$\mathcal{A}(\bm{P},\bm{y}_{v})+\mathcal{A}_{0}\succ 0$ if and only if there
is no $\bm{Z}\in\mathbb{S}^{3n_{x}}$ in the form of

$\displaystyle\bm{Z}=\mathrm{Blkdiag}\left(\begin{bmatrix}\bm{Z}_{1}&*\\\
\bm{Z}_{2}^{\top}&\bm{Z}_{3}\end{bmatrix},\bm{Z}_{4}\right),$ (7a) such that
$\bm{Z}\succeq 0$ and $\bm{Z}\neq 0$ satisfying
$\displaystyle\bm{Z}_{4}=\bm{Z}_{1}\bm{A}^{\top}+\bm{A}\bm{Z}_{1}+\bm{G}\bm{Z}_{2}^{\top}+\bm{Z}_{2}\bm{G}^{\top}$
(7b)
$\displaystyle\sum_{j=1}^{n_{x}}\sum_{k=1}^{n_{x}}\bar{C}_{i,jk}Z_{1,jk}=0,\;\;\forall
i\in\mathbb{I}(n_{x}n_{y})$ (7c)
$\displaystyle\mathrm{tr}(\bm{Z}_{3})\leq\gamma_{l}^{2}\,\mathrm{tr}(\bm{Z}_{1}).$
(7d)

The proof of the above theorem is presented in B. Problem (7) is referred to
throughout the section as the alternative problem for (6a). To answer Q1, we
need to look into inequality (7d). Realize that the term
$\mathrm{tr}(\bm{Z}_{3})$ is lower bounded by zero since $\bm{Z}\succeq 0$.
Consequently, we get

$\displaystyle 0\leq\mathrm{tr}(\bm{Z}_{3})$
$\displaystyle\leq\gamma_{l}^{2}\,\mathrm{tr}(\bm{Z}_{1}).$ (8)

It is seen from (8) that the feasible set of (7) expands as the value of
$\gamma_{l}$ increases. In a particular case when $\gamma_{l}=0$, the term
$\mathrm{tr}(\bm{Z}_{3})$ is enforced to be equal to zero. Since eigenvalues
of $\bm{Z}_{3}$ cannot be negative, then they must be equal to zero. On the
other hand, if $\gamma_{l}$ is considerably large, then there is not much
restriction on the upper bound of $\mathrm{tr}(\bm{Z}_{3})$. Hence, it is
important to get smaller Lipschitz constant $\gamma_{l}$ since this will make
the problem described in (7) to have narrower feasible space. This finding is
in accordance with empirical observations from the literature: SDP (4) has no
solution if $\gamma_{l}$ is selected to be sufficiently large—such large
Lipschitz constant is labeled as conservative. This provides an answer for Q1.
The next section introduces the SSP for Lipschitz NDSs and presents some
properties useful for the development of a more efficient BnB algorithm.

##  4 Sensor Selection and NDS Observability

In this section, we wish to answer the question below:

Q2: What is the impact of the number of placed sensors on
$\gamma_{l}$-observability?

The above question may seem trivial since it is expected that dynamic systems
have higher degree of observability as more sensors are placed or utilized.
The objective of this section is to demonstrate the above conjecture from
$\gamma_{l}$-observability stand point. To begin with, we formally construct
the SSP for NDS (1) as follows. Let $\gamma_{i}\in\\{0,1\\}$ be a binary
variable that ascertains the activation or deactivation of sensor on each
subsystem $i$. That is, $\gamma_{i}=1$ if the sensor measuring subsystem $i$
is activated and $\gamma_{i}=0$ otherwise. These variables can be combined
into $\bm{\gamma}$, giving
$\bm{\gamma}:=\left[\gamma_{1}\,\,\gamma_{2}\,\,\cdots\,\,\gamma_{N}\right]^{\top}$.
The following conventions are adopted in this section.

######  Definition 2.

Let $\mathcal{S}_{\gamma}$ be an $N$-tuple representing the selection of
sensors, i.e.,
$\mathcal{S}_{\gamma}:=(\gamma_{1},\gamma_{2},\ldots,\gamma_{N})$. The set of
active sensor is constructed as
$\mathcal{S}_{\gamma}^{(a)}:=\left\\{\gamma\in\mathcal{S}_{\gamma}\,|\,\gamma=1\right\\}$.
Let
$\bar{\mathcal{S}}_{\gamma}:=(\gamma_{1}\bm{1}_{n_{y_{1}}},\gamma_{2}\bm{1}_{n_{y_{2}}},\ldots,\gamma_{N}\bm{1}_{n_{y_{N}}})$
where $\mathrm{card}(\bar{\mathcal{S}}_{\gamma})=n_{y}$. The selection of
sensors in matrix form can be written as
$\bm{\Gamma}:=\mathrm{Blkdiag}\big{(}\gamma_{1}\bm{I}_{n_{y_{1}}},\gamma_{2}\bm{I}_{n_{y_{2}}},\ldots,\gamma_{N}\bm{I}_{n_{y_{N}}}\big{)}$.

By using the above definition, NDS (1) together with sensor selection can be
conveniently expressed as

$\displaystyle\dot{\bm{x}}(t)$
$\displaystyle=\bm{A}\bm{x}(t)+\bm{G}\bm{f}(\bm{x})+\bm{B}\bm{u}(t)$ (9a)
$\displaystyle\bm{y}(t)$ $\displaystyle=\bm{\Gamma}\bm{C}\bm{x}(t).$ (9b)

Additionally, it is also beneficial to consider a set
$\mathcal{G}_{\gamma}\subseteq\\{0,1\\}^{N}$ representing logistic constraints
and availability of sensors such that $\bm{\gamma}\in\mathcal{G}_{\gamma}$
might be imposed. This in turn allows a simplified high level formulation of a
general SSP

$\displaystyle\mathbf{(P1)}\;\;\;\operatorname*{minimize}\;\;\;$
$\displaystyle\bm{c}^{\top}\bm{\gamma}+\mathrm{EstimationObjective}$
$\displaystyle\operatorname*{subject\ to}\;\;\;$
$\displaystyle\eqref{eq:gen_dynamic_systems_sen_sel},\;\bm{\gamma}\in\mathcal{G}_{\gamma},\;\mathrm{EstimationConstraints}.$

The objectives of P1 are threefold: (i) performing state estimation for NDS
(9) while (ii) utilizing smallest number of sensors as possible (or satisfying
a given constraint over the collections of library of sensors) and (iii)
optimizing a specific estimation metric. Vector $\bm{c}\in\mathbb{R}^{N}_{+}$
in the objective function of P1 assigns weights for each sensor $\gamma_{i}$.
Given P1, SSP for Lipschitz NDSs assuming no estimation objective (this is
considered in Section 6), can be constructed by incorporating (9) into LMI
(4). The resulting problem is given below.

$\displaystyle\mathbf{(P2)}\;\operatorname*{minimize}_{\bm{P},\bm{Y},\epsilon,\bm{\gamma}}\;\;\bm{c}^{\top}\bm{\gamma}$
(10a) $\displaystyle\operatorname*{subject\
to}\;\;\begin{bmatrix}\bm{A}^{\top}\bm{P}+\bm{P}\bm{A}+\epsilon\gamma_{l}^{2}\bm{I}&\\\
-\bm{Y}\bm{\Gamma}\bm{C}-\bm{C}^{\top}\bm{\Gamma}\bm{Y}^{\top}&*\\\
\bm{G}^{\top}\bm{P}&-\epsilon\bm{I}\end{bmatrix}\preceq 0$ (10b)
$\displaystyle\quad\quad\;\;\bm{P}\succ
0,\;\epsilon>0,\;\bm{\gamma}\in\mathcal{G}_{\gamma},\;\bm{\gamma}\in\\{0,1\\}^{N}.\vspace*{-0.05cm}$
(10c)

In P2 the objective is to minimize the number of the activated sensors while
(a) finding a stabilizing observer gain matrix $\bm{L}$ and (b) satisfying
sensor’s logistic constraints. As $\bm{\Gamma}$ is a binary matrix variable,
P2 is a nonconvex optimization problem with mixed-integer bilinear matrix
inequalities (MIBMIs) because of the $\bm{Y}\bm{\Gamma}$ term. Our approach to
address this problem is given in the next section. Following (5), the term
$\bm{Y}\bm{\Gamma}\bm{C}+\bm{C}^{\top}\bm{\Gamma}\bm{Y}^{\top}$ is equal to

$\displaystyle\sum_{j=1}^{n_{y}}\gamma_{j}\left(\sum_{i=1}^{n_{x}}Y_{ij}\check{\bm{C}}_{S,ij}\right)=\sum_{k=1}^{n_{x}n_{y}}\bar{\gamma}_{k}y_{v,k}\bar{\bm{C}}_{k},$
(11)

where each $\bar{\gamma}_{k}$ corresponds to
${\gamma}_{j}\in\bar{\mathcal{S}}_{\gamma}$. As such we may write
$\bar{\gamma}_{k}\triangleleft{\gamma}_{j}$. For the sake of analysis, let us
assume that $\bm{\gamma}$ is fixed. From (11), define
$\mathcal{A}_{Y,\Gamma}({\bm{y}_{v}})$ as follows

$\displaystyle\mathcal{A}_{Y,\Gamma}({\bm{y}_{v}})$
$\displaystyle:=\mathrm{Blkdiag}\left(\begin{bmatrix}\sum_{k=1}^{n_{x}n_{y}}\bar{\gamma}_{k}y_{v,k}\bar{\bm{C}}_{k}&*\\\
\bm{O}&\bm{O}\end{bmatrix},\bm{O}\right),$ (12a) such that, by having
$\mathcal{A}_{\Gamma}(\bm{P},\bm{y}_{v}):=\mathcal{A}_{P}(\bm{P})+\mathcal{A}_{Y,\Gamma}({\bm{y}_{v}})$,
matrix inequality (6a) now becomes
$\displaystyle\mathcal{A}_{\Gamma}(\bm{P},\bm{y}_{v})+\mathcal{A}_{0}=\mathcal{A}_{P}(\bm{P})+\mathcal{A}_{Y,\Gamma}({\bm{y}_{v}})+\mathcal{A}_{0}\succ
0.$ (12b)

The next result is established due to Theorem 1.

######  Theorem 2.

Consider a certain sensor configuration $\mathcal{S}_{\gamma}$. There exist
$\bm{P}\in\mathbb{S}^{n_{x}}$ and $\bm{y}_{v}\in\mathbb{R}^{n_{x}n_{y}}$
satisfying $\mathcal{A}_{\Gamma}(\bm{P},\bm{y}_{v})+\mathcal{A}_{0}\succ 0$ if
and only if there is no $\bm{Z}\in\mathbb{S}^{3n_{x}}$ provided in the form of
(7a) such that $\bm{Z}\succeq 0$ and $\bm{Z}\neq 0$ satisfying (7b), (7d), and
$\bm{Z}_{1}\in\mathcal{Z}\big{(}\mathcal{S}_{\gamma}^{(a)}\big{)}$ where

$\displaystyle\mathcal{Z}\big{(}\mathcal{S}_{\gamma}^{(a)}\big{)}:=$
$\displaystyle\Bigg{\\{}\vphantom{\Bigg{|}}\bm{Z}_{1}\in\mathbb{S}^{n_{x}}\,\Bigg{|}\,\sum_{j=1}^{n_{x}}\sum_{k=1}^{n_{x}}\bar{\gamma}_{i}\bar{C}_{i,jk}Z_{1,jk}=0,$
$\displaystyle\quad\;\;\forall
i\in\mathbb{I}(n_{x}n_{y})\;\wedge\;\bar{\gamma}_{i}\triangleleft{\gamma}_{m},\;\gamma_{m}\in\mathcal{S}_{\gamma}^{(a)}\vphantom{\Bigg{|}}\Bigg{\\}}.$
(13)

It follows from Theorem 2 (see C for the proof) that including more sensors
will reduce the feasible set of the alternative problem and thus, in general,
utilizing more sensors reduces the feasible set of the alternative problem. By
reducing this set, one can expect that the constraint will be inconsistent
thus making the alternative problem infeasible and consequently, giving a
feasible solution for (12b). This provides an answer for Q2. From Theorem 2,
the following additional result is obtained.

######  Corollary 1.

Suppose $\mathcal{S}_{\gamma,1}$ and $\mathcal{S}_{\gamma,2}$ are two distinct
sensor’s configurations so that
$\mathcal{S}_{\gamma,1}^{(a)}\subset\mathcal{S}_{\gamma,2}^{(a)}$. If (12b)
has no solution for $\mathcal{S}_{\gamma,2}$, then it also has no solution for
$\mathcal{S}_{\gamma,1}$.

The above corollary essentially states that if the system is not
$\gamma_{l}$-observable for a particular sensor configuration
$\mathcal{S}_{\gamma,2}$, then the system is also not $\gamma_{l}$-observable
for any sensor configuration, represented by $\mathcal{S}_{\gamma,1},$ having
fewer number of sensors if
$\mathcal{S}_{\gamma,1}^{(a)}\subset\mathcal{S}_{\gamma,2}^{(a)}$. The proof
of Corollary 1 follows directly from the fact that
$\mathcal{S}_{\gamma,1}^{(a)}\subset\mathcal{S}_{\gamma,2}^{(a)}$ implies
$\mathcal{Z}\big{(}\mathcal{S}_{\gamma,1}^{(a)}\big{)}\subset\mathcal{Z}\big{(}\mathcal{S}_{\gamma,2}^{(a)}\big{)}$.
That is, since the alternative problem with set of constraints
$\mathcal{Z}\big{(}\mathcal{S}_{\gamma,2}^{(a)}\big{)}$ has nonempty feasible
set, then this set is never empty for any $\mathcal{S}_{\gamma,1}$ since
$\mathcal{Z}\big{(}\mathcal{S}_{\gamma,1}^{(a)}\big{)}$ is less constrained.
By contrapositive, then it is also true that the system is
$\gamma_{l}$-observable for any sensor configuration $\mathcal{S}_{\gamma,2}$
provided that it is also $\gamma_{l}$-observable for a certain sensor
configuration $\mathcal{S}_{\gamma,1}$ such that
$\mathcal{S}_{\gamma,1}\subset\mathcal{S}_{\gamma,2}$.

Now, a natural question to ask is whether there is a relation between the
value of Lipschitz constant $\gamma_{l}$ and the number of activated sensors
to achieve $\gamma_{l}$-observability for SSP. Consider a SSP where the
objective is to achieve $\gamma_{l}$-observability with minimum number of
activated sensors. From previous results, we know that minimizing $\gamma_{l}$
while incorporating more sensors reduces the feasible set of the alternative
problem. However, as we want to incorporate the least number of sensors as
possible, then we can only rely on using the smallest Lipschitz constant
$\gamma_{l}$. The rationale here is, if the alternative problem can be made
infeasible from (8) by reducing the Lipschitz constant, then we can utilize
less number of activated sensors to achieve $\gamma_{l}$-observability. In the
next section, we discuss our approach to solve SSP formulated in P2.

##  5 Solving the SSP Through Customized BnB

###  5.1 From MIBMI to Convex MISDP

Our approach involves transforming the SSP to a convex MISDP. Specifically,
our prior work [16] reformulates the SSP P2 into a covex MISDP using
McCormick’s relaxation [25]. This is not the only method to convert MIBMI into
MISDP: the big-M method, which is popular in disjunctive programming, can also
be employed—see [6]. The McCormick’s reformulation is performed by defining a
new matrix variable $\bm{M}:=\bm{Y}\bm{\Gamma}$ where
$\bm{M}\in\mathbb{R}^{n_{x}\times n_{y}}$ given that $\bm{Y}$ is bounded such
that
$\stackunder[1.2pt]{$\bm{Y}$}{\rule{3.44444pt}{0.32289pt}}\leq\bm{Y}\leq\bar{\bm{Y}}$,
while the big-M assumes that $-L\bm{1}\leq\bm{Y}\leq L\bm{1}$ for a
sufficiently large constant $L>0$. For convenience, we consider $\bm{Y}\in
B(\bm{Y})$ where $B(\bm{Y})$ is defined as

$\displaystyle
B(\bm{Y}):=\begin{cases}\stackunder[1.2pt]{$\bm{Y}$}{\rule{3.44444pt}{0.32289pt}}\leq\bm{Y}\leq\bar{\bm{Y}},&\;\text{for
McCormick},\\\ -L\bm{1}\leq\bm{Y}\leq L\bm{1},&\;\text{for Big-M}.\\\
\end{cases}$

It is apparent here that the bounds on $\bm{Y}$ are determined by a single
constant $L$ for big-M and matrices
$\stackunder[1.2pt]{$\bm{Y}$}{\rule{3.44444pt}{0.32289pt}},\bar{\bm{Y}}$ for
McCormick, making big-M to be a special case of McCormick’s relaxation. The
resulting problem can now be written as

$\displaystyle(\mathbf{P3})\;\operatorname*{minimize}_{\bm{P},\bm{Y},\bm{M},\epsilon,\bm{\gamma}}\;\;\bm{c}^{\top}\bm{\gamma}$
(14a) $\displaystyle\operatorname*{subject\
to}\;\;\begin{bmatrix}\bm{A}^{\top}\bm{P}+\bm{P}\bm{A}+\epsilon\gamma_{l}^{2}\bm{I}&\\\
-\bm{M}\bm{C}-\bm{C}^{\top}\bm{M}^{\top}&*\\\
\bm{G}^{\top}\bm{P}&-\epsilon\bm{I}\end{bmatrix}\preceq 0$ (14b)
$\displaystyle\quad\quad\quad\quad\;\;\bm{M}=\bm{Y}\bm{\Gamma},\;\bm{Y}\in
B(\bm{Y}),\;\eqref{eq:sen_sel_obs_2}.\vspace*{-0.05cm}$ (14c)

Without loss of generality, it is assumed throughout the section that the
number of variables in $\bm{\gamma}$ are equal to the number of rows in
$\bm{C}$. By applying big-M or McCormick’s relaxation, the nonconvex MISDP P3
can be converted to an equivalent convex MISDP given in the next theorem.

######  Theorem 3.

Problem P3 is equivalent to

$\displaystyle\mathbf{(P4)}\;\operatorname*{minimize}_{\bm{v}}\;\;$
$\displaystyle\bm{c}^{\top}_{v}\bm{v}$ (15a)
$\displaystyle\operatorname*{subject\ to}\;\;$
$\displaystyle\mathcal{L}_{v}(\bm{v})+\mathcal{C}_{v}\succeq 0,$ (15b)
$\displaystyle v_{i}\in\\{0,1\\},\;\;\forall i\in\mathbb{I}_{v}(\bm{\gamma}),$
(15c) where $\mathcal{L}_{v}(\cdot)$ denotes LMI terms, $\mathcal{C}_{v}$
denotes the constant terms, and vector $\bm{v}\in\mathbb{R}^{n_{v}}$ populates
all decision variables in the following fashion
$\displaystyle\bm{v}:=\begin{bmatrix}\mathrm{vec}(\bm{P})^{\top}\;\;\epsilon\;\;\mathrm{vec}(\bm{M})^{\top}\;\;\mathrm{vec}(\bm{Y})^{\top}\;\;\bm{\gamma}^{\top}\end{bmatrix}^{\top},$
(15d)

with
$n_{v}={\bar{n}_{x}+1+2n_{x}n_{y}+n_{y}},\bar{n}_{x}:=\tfrac{1}{2}n_{x}(n_{x}+1)$.

The proof of Theorem 3 is detailed in D. The notation
$\mathbb{I}_{v}(\bm{\gamma})$ in (15c) is useful for indicating the index $i$
of $\bm{v}$ such that $v_{i}=\gamma_{j}$ for each corresponding
$j\in\mathbb{I}(\bm{\gamma})$; similar notations are also used to label any
variable in $\bm{v}$ that corresponds to $\bm{P}$, $\bm{Y}$, $\bm{M}$, and
$\epsilon$. Specifically, the inequality $\bm{\Xi}\,\bm{\nu}\leq\bm{\varphi}$
which is embedded in LMI (15b) essentially encompasses the McCormick’s
envelope on the bilinear equality $M_{ij}=Y_{ij}\gamma_{j}$ expressed as

$\displaystyle\begin{bmatrix}1&-1&-\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\\\
-1&1&\bar{Y}_{ij}\\\ 1&0&-\bar{Y}_{ij}\\\
-1&0&\bar{Y}_{ij}\end{bmatrix}\begin{bmatrix}M_{ij}\\\ Y_{ij}\\\
\gamma_{j}\end{bmatrix}\leq\begin{bmatrix}-\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\\\
\bar{Y}_{ij}\\\ 0\\\ 0\end{bmatrix},$ (16)

which is equivalent to the previous bilinear constraint (similarly for big-M).
Now that we have a convex MISDP, problem P4 can be solved using one of the
aforementioned techniques. In the next section, we develop a customized BnB
algorithm that can be used to compute both optimal and suboptimal solutions
for P4.

###  5.2 Structure-Exploiting BnB Algorithm

A natural way to solve such MISDP is through BnB algorithm. In a general
purpose BnB algorithm, an SDP is solved in every node, in which the SDP is
obtained from relaxing the integer variables to their corresponding convex
relaxations. Thus, the complexity of BnB algorithm stems from the difficulties
in solving the SDPs and the number of SDPs that need to be solved. In a
practical situation, finding a relatively good solution in a much shorter time
is more desirable than finding an optimal one. Despite that some general
purpose BnB algorithms for MISDPs are available, these conventional solvers
cannot be configured to make them amenable for exploiting the structure of the
problem at hand in order to reduce the complexity of solving SSP.

1 input: $\mathcal{F}$, $\epsilon_{f}$, $\mathrm{maxBranch}$

2 initialize: $\mathcal{C}=\emptyset$, $\mathcal{W}=\emptyset$,
$\mathrm{LB},\mathrm{UB}=\infty$, $k=0$, $\bm{v}^{*}=\bm{0}$

3 process root node: $\mathcal{T}\leftarrow\mathrm{RelaxSol}(\mathcal{F})$

4 update: $\mathrm{LB}\leftarrow\bm{c}_{v}^{\top}\mathcal{T}(2)$,
$\mathcal{C}\leftarrow\mathcal{C}\cup\\{\mathcal{T}\\}$

5 while  _$\mathrm{UB}-\mathrm{LB} >\epsilon_{f}$ and $k<\mathrm{maxBranch}$_
do

6 $\tilde{\mathcal{T}}=\operatorname*{arg\
min}_{\mathcal{T}_{k}\in\mathcal{C}}\bm{c}_{v}^{\top}\mathcal{T}_{k}(2)$,
$\mathcal{C}\leftarrow\mathcal{C}\setminus\\{\tilde{\mathcal{T}}\\}$,
$\mathrm{LB}\leftarrow\bm{c}_{v}^{\top}\tilde{\mathcal{T}}(2)$

7 if  _$\tilde{\mathcal{T}}(2)_{i}\in\\{0,1\\},\;\forall
i\in\mathbb{I}_{v}(\bm{\gamma})$_ then

8 $\tilde{\mathcal{T}}(3)\leftarrow\tilde{\mathcal{T}}(2)$

9

10 else

11
$\\{\tilde{\mathcal{T}}(3),\mathcal{W}\\}\leftarrow\mathrm{FeasSol}(\tilde{\mathcal{T}},\mathcal{W})$

12 find: $i\in\mathbb{I}_{v}(\bm{\gamma})$ such that
$\tilde{\mathcal{T}}(2)_{i}\notin\\{0,1\\}$ and
$\tilde{\mathcal{T}}(2)_{i}=\min_{j\in\mathbb{I}_{v}(\bm{\gamma})}\left\lvert\tilde{\mathcal{T}}(2)_{j}-\tfrac{1}{2}\right\rvert$

13 $\mathcal{F}_{L}:=\\{\bm{v}\in\tilde{\mathcal{T}}(1)\,|\,v_{i}\leq 0\\}$

14 if
_$\forall\,\bm{\gamma}\in\\{0,1\\}^{n_{y}},\bm{\gamma}\in\mathbfcal{I}(\mathcal{F}_{L})$
, then $\bm{\gamma}\notin\mathcal{W}$_ then

15 $\mathcal{T}_{L}\leftarrow\mathrm{RelaxSol}(\mathcal{F}_{L})$,
$\mathcal{C}\leftarrow\mathcal{C}\cup\\{\mathcal{T}_{L}\\}$

16 $\mathcal{F}_{R}:=\\{\bm{v}\in\tilde{\mathcal{T}}(1)\,|\,v_{i}\geq 1\\}$

17 if
_$\forall\,\bm{\gamma}\in\\{0,1\\}^{n_{y}},\bm{\gamma}\in\mathbfcal{I}(\mathcal{F}_{R})$
, then $\bm{\gamma}\notin\mathcal{W}$_ then

18 $\mathcal{T}_{R}\leftarrow\mathrm{RelaxSol}(\mathcal{F}_{R})$,
$\mathcal{C}\leftarrow\mathcal{C}\cup\\{\mathcal{T}_{R}\\}$

19 $k\leftarrow k+1$

20

21 if  _$\bm{c}_{v}^{\top}\tilde{\mathcal{T}}(3) <\mathrm{UB}$_ then

22 $\mathrm{UB}\leftarrow\bm{c}_{v}^{\top}\tilde{\mathcal{T}}(3)$,
$\bm{v}^{*}\leftarrow\tilde{\mathcal{T}}(3)$

23 foreach  _$\mathcal{T}_{j}\in C$_ do

24 if  _$\bm{c}_{v}^{\top}\mathcal{T}_{j}(2) >\mathrm{UB}$_ then

25 $\mathcal{C}\leftarrow\mathcal{C}\setminus\\{\mathcal{T}_{j}\\}$

26

27

28

29output: $\mathrm{LB}$, $\mathrm{UB}$, $\bm{v}^{*}$

Algorithm 1 SE-BnB Algorithm

To that end, we present a customized BnB algorithm that can be utilized to
compute optimal and suboptimal solutions for P4 in a more efficient manner,
owing it to (a) its ability to exploit problem structure and (b) heuristics to
quickly find good suboptimal solutions. To proceed, define

$\displaystyle\mathcal{F}$
$\displaystyle:=\\{\bm{v}\,|\,\mathcal{L}_{v}(\bm{v})+\mathcal{C}_{v}\succeq
0,\;v_{i}\in\\{0,1\\},\;\forall i\in\mathbb{I}_{v}(\bm{\gamma})\\}$
$\displaystyle\mathcal{R}$
$\displaystyle:=\\{\bm{v}\,|\,\mathcal{L}_{v}(\bm{v})+\mathcal{C}_{v}\succeq
0,\;v_{i}\in[0,1],\;\forall i\in\mathbb{I}_{v}(\bm{\gamma})\\},$

i.e., $\mathcal{F}$ denotes the feasible set for P4 whereas $\mathcal{R}$ is
its convex relaxation. Note that we have $\mathcal{F}\subseteq\mathcal{R}$
because in $\mathcal{R}$, the integer constraints are relaxed. The proposed
BnB algorithm solves two convex problems on each BnB node in order for it to
be able to find a good suboptimal solution immediately: the convex relaxation
of P4 to get a local lower bound, and P4 with fixed SAs combination to get a
global upper bound. Within this section, it is supposed that each node $i$ is
representable by a $3$-tuple
$\mathcal{T}_{i}:=\left(\mathcal{F}_{i},\stackunder[1.2pt]{$\bm{v}$}{\rule{3.44444pt}{0.32289pt}},\bar{\bm{v}}\right)$,
where

$\displaystyle\stackunder[1.2pt]{$\bm{v}$}{\rule{3.44444pt}{0.32289pt}}=\operatorname*{arg\
min}_{\bm{v}\in\mathcal{R}_{i}}\bm{c}_{v}^{\top}{\bm{v}},\;\quad\;\bar{\bm{v}}\in\\{\bm{v}\;|\;\bm{v}\in\mathcal{F}_{i}\\},$

such that $\mathcal{F}_{i}=\mathcal{T}_{i}(1)$,
$\stackunder[1.2pt]{$\bm{v}$}{\rule{3.44444pt}{0.32289pt}}=\mathcal{T}_{i}(2)$,
and $\bar{\bm{v}}=\mathcal{T}_{i}(3)$. We also use an abstract set
$\mathcal{I}(\cdot)\subseteq\mathbb{R}^{n_{y}}$ where
$\mathcal{I}(\mathcal{T}_{i})\subset\mathcal{T}_{i}(1)$ to represent the
restrictions (or integer cut) on integer variable $v_{i}$, where
$i\in\mathbb{I}_{v}(\bm{\gamma})$, which is added in the branching process of
BnB algorithm. To stop the algorithm when a suboptimal solution is sought, two
stopping criteria characterize by $\epsilon_{f}\in\mathbb{R}_{+}$ and
$\mathrm{maxBranch}\in\mathbb{N}_{+}$ are employed. Let $\mathrm{LB}$ and
$\mathrm{UB}$ be the known best lower and upper bounds at the end of the
algorithm. When the optimality gap $\mathrm{UB}-\mathrm{LB}$ is zero, the
solution is optimal. However, if $\mathrm{UB}-\mathrm{LB}\leq\epsilon_{f}$ the
solution is regarded to be $\epsilon_{f}$-optimal. Alternatively, the solution
is only suboptimal. The notation $\mathrm{maxBranch}$, sometimes abbreviated
as $\mathrm{mBr}$, determines the number of maximum allowable branching. These
BnB routines, referred to as structure-exploiting BnB (SE-BnB), are detailed
in Algorithm 1. The algorithmic function $\mathrm{RelaxSol}(\cdot)$ solves the
relaxed problem
$\stackunder[1.2pt]{$f$}{\rule{3.44444pt}{0.32289pt}}=\min_{\bm{v}\in\mathcal{R}}c_{v}^{\top}{\bm{v}}$
for a given $\mathcal{F}$ where
$\stackunder[1.2pt]{$\bm{v}$}{\rule{3.44444pt}{0.32289pt}}$ is the minimizer
and wrap the results into $\mathcal{T}$. In addition to this, Algorithm 1 uses
a heuristic procedure for efficiently finding a feasible solution for P4 on
each node, if such solution exists. This heuristic generates a random
combination of sensors $\bm{\gamma}_{c}$. To do so, we define the set
$\mathcal{W}$ as

$\displaystyle\mathcal{W}:=\left\\{\bm{\gamma}\in\\{0,1\\}^{n_{y}}\,\Big{|}\,\begin{bmatrix}\bm{w}^{\top}\,\,\bm{\gamma}^{\top}\end{bmatrix}^{\top}\notin\mathcal{F},\;\forall\bm{w}\right\\},$

i.e., the set of combinations of sensors that are infeasible for P4. By
exploiting the property proposed in Corollary 1, we are able to effectively
generate a candidate of sensor combination $\bm{\gamma}_{c}$ that is not
positively infeasible for P4. This heuristic, embedded in the function
$\mathrm{FeasSol}(\cdot)$, is presented in Algorithm 2. This idea is
implemented in Steps 1 and 1 of Algorithm 1. Now, some modifications for
making the SE-BnB algorithm to be more efficient are discussed.

1 input: $\tilde{\mathcal{T}}$, $\mathcal{W}$

2 candidate generation: generate a random sensor combination
$\bm{\gamma}_{c}\in\\{0,1\\}^{n_{y}}$ such that
$\bm{\gamma}_{c}\notin\mathcal{W}$ and
$\bm{\gamma}_{c}\in\mathcal{I}(\tilde{\mathcal{T}}(1))$

3 solve:
$\bar{f}=\min_{\bm{v}\in\tilde{\mathcal{T}}(1),\,v_{i}=\gamma_{c,i}}c_{v}^{\top}{\bm{v}}$
(problem P6) with minimizer $\bar{\bm{v}}$

4 if  _$\bar{f}=\infty$_ then

5 $\mathcal{W}\leftarrow\mathcal{W}\cup\\{\bm{\gamma}_{c}\\}$

output: $\bar{\bm{v}}$, $\mathcal{W}$

Algorithm 2 $\mathrm{FeasSol}$

  * •

In $\mathrm{RelaxSol}(\cdot)$, the relaxed P4 is solved. Realize that an
integer cut is appended in every node (apart from the root node) that
corresponds to whether a particular relaxed integer variable, say $v_{i}$ for
$i\in\mathbb{I}_{v}(\bm{\gamma})$, is constrained to be $v_{i}\leq 0$ or
$v_{i}\geq 1$. As $v_{i}\in[0,1]$, then the only feasible solution is
$v_{i}=0$ or $v_{i}=1$. Either case, the resulting McCormick’s envelope (or
big-M) yields the best relaxation which makes (16) redundant. At this
instance, constraint (16) can be discarded and replaced with $M_{ij}=0$ or
$M_{ij}=Y_{ij}$, depending on the value of $v_{i}$.

  * •

Recall that the SE-BnB algorithm attempts to find a feasible solution on each
node by fixing $\bm{\gamma}$ and solve P4—performed in
$\mathrm{FeasSol}(\cdot)$. Since P4 is equivalent to P3 and P3 becomes a SDP
when $\bm{\gamma}$ is fixed, then one can solve P3 to find a feasible
solution. Furthermore, given a fixed $\bm{\gamma}$, the equality constraint
$\bm{M}=\bm{Y}\bm{\Gamma}$ becomes redundant and hence, this constraint can be
removed as well as the variable $\bm{M}$ from P3. Moreover, when a certain
variable $\bm{\gamma}_{i}=0$, due to redundancy, then the corresponding row of
$\bm{C}$ and column of $\bm{Y}$ can also be removed to obtain smaller number
of variables (essentially, only consider the nonzero row of $\bm{C}$ and
column of $\bm{Y}$).

The aforementioned procedures enable the SE-BnB to exploit the structure of
SSP and as a result, the complexity of solving the associated SDPs required to
find bounds on each node of the BnB tree can be lowered.

##  6 Some Important Extensions

###  6.1 Controllability and Actuator Selection

In this section, NDS controllability and actuator selection problem (ASP) are
discussed, in which we consider a similar concept to quantify controllability
for a given $\gamma_{l}$.

######  Definition 3.

NDS (1) is said to be $\gamma_{l}$-controllable if and only if there exist
$\bm{Q}\in\mathbb{S}^{n_{x}}_{++}$, $\bm{X}\in\mathbb{R}^{n_{u}\times n_{x}}$,
and $\sigma\in\mathbb{R}_{++}$ such that the LMI below is feasible

$\displaystyle\small\begin{bmatrix}\bm{Q}\bm{A}^{\top}+\bm{A}\bm{Q}-\bm{X}^{\top}\bm{B}^{\top}-\bm{B}\bm{X}+\sigma\bm{G}\bm{G}^{\top}&*\\\
\bm{Q}&-\frac{\sigma}{\gamma_{l}^{2}}\bm{I}\end{bmatrix}$ $\displaystyle\prec
0.$ (17)

If the NDS is $\gamma_{l}$-controllable, then from the above LMI, the
stabilizing controller gain matrix is given as $\bm{K}=\bm{X}\bm{Q}^{-1}$ with
control action $\bm{u}(t)=-\bm{K}\bm{x}(t)$ [26]. It is worth noting that,
unlike (4), LMI (17) is a sufficient condition. A similar question to Q1 now
emerges: what role does $\gamma_{l}$ have on $\gamma_{l}$-controllability? By
performing a similar analysis as in Section 3 (not shown here due to space
constraint), it is revealed that if $\gamma_{l}$ is made to be sufficiently
small, then there exists at least a solution for LMI (17) (if $\sigma$ is
fixed).

Next, the relation between actuators and NDS controllability is discussed. Let
$\pi_{i}\in\\{0,1\\}$ be a binary variable representing the activation or
deactivation of actuator/control node on each subsystem $i$ such that
$\pi_{i}=1$ if control node at $i$ is activated and $\pi_{i}=0$ otherwise. For
compactness, let us define
$\bm{\pi}:=\left[\pi_{1}\,\,\pi_{2}\,\,\cdots\,\,\pi_{N}\right]^{\top}$ and
$\bm{\Pi}$ as

$\bm{\Pi}:=\mathrm{Blkdiag}\big{(}\pi_{1}\bm{I}_{n_{u_{1}}},\pi_{2}\bm{I}_{n_{u_{2}}},\ldots,\pi_{N}\bm{I}_{n_{u_{N}}}\big{)}.$

As such, ASP for NDS (1) can be expressed as

$\displaystyle\dot{\bm{x}}(t)=\bm{A}\bm{x}(t)+\bm{G}\bm{f}(\bm{x})+\bm{B}\bm{\Pi}\bm{u}(t),\quad\bm{y}(t)=\bm{C}\bm{x}(t),$
(18)

from which a simplified high level formulation for ASP can be formulated as
follows

$\displaystyle(\mathbf{P5})\;\;\;\operatorname*{minimize}\;\;\;$
$\displaystyle\bm{c}^{\top}\bm{\pi}+\mathrm{ControlObjective}$
$\displaystyle\operatorname*{subject\ to}\;\;\;$
$\displaystyle\eqref{eq:gen_dynamic_systems_act_sel},\;\bm{\pi}\in\mathcal{G}_{\pi},\;\mathrm{ControlConstraints}.$

The goals in the above are threefold: (i) performing feedback stabilizing
control on NDS (18) while (ii) utilizing smallest number of actuators as
possible (or satisfying a given constraint over the collections of library of
actuators) and (iii) optimizing a specific control metric. In P5, vector
$\bm{c}\in\mathbb{R}^{N}_{+}$ assigns weights for each actuator $\pi_{i}$
whereas $\mathcal{G}_{\pi}\subseteq\\{0,1\\}^{N}$ represents actuator’s
logistic constraints. The ASP for Lipschitz NDSs can be formulated as

$\displaystyle(\mathbf{P6})\;\operatorname*{minimize}_{\bm{Q},\bm{X},\epsilon,\bm{\pi}}\;\;\bm{c}^{\top}\bm{\pi}$
(19a) $\displaystyle\operatorname*{subject\
to}\;\;\begin{bmatrix}\bm{Q}\bm{A}^{\top}+\bm{A}\bm{Q}+\sigma\bm{G}\bm{G}^{\top}&\\\
-\bm{X}^{\top}\bm{\Pi}\bm{B}^{\top}-\bm{B}\bm{\Pi}\bm{X}&*\\\
\bm{Q}&-\frac{\sigma}{\gamma_{l}^{2}}\bm{I}\end{bmatrix}\preceq 0$ (19b)
$\displaystyle\quad\quad\;\;\bm{Q}\succ
0,\;\sigma>0,\;\bm{\pi}\in\mathcal{G}_{\pi},\;\bm{\pi}\in\\{0,1\\}^{N}.\vspace*{-0.05cm}$
(19c)

In P6 the objective is to minimize the number of the activated actuators while
(a) finding a stabilizing controller gain matrix $\bm{K}$ and (b) satisfying
actuator’s logistic constraints. Notice that P6 is also a nonconvex
optimization problem with a MIBMI constraint. Further analysis shows that, in
general, in order to employ less number of actuators, one should obtain the
least possible Lipschitz constant.

###  6.2 NDS Parameterization for Non-Lipschitz Systems

Note that many observer/controller designs for function sets other than
Lipschitz continuous such as one-sided Lipschitz, quadratically inner-bounded,
bounded Jacobian, and quadratically bounded—see [12, 14, 15] for notable
examples—can also be considered and utilized to construct SASPs. For example,
the nonlinearities in a NDS may satisfy one-sided Lipschitz (OSL) and
quadratically inner-bounded (QIB) conditions—see Definition 4.

######  Definition 4 (OSL & QIB).

The mapping $\bm{f}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{g}}$ in (1) is
locally OSL in $\mathbfcal{X}$ if for any
$\bm{x},\hat{\bm{x}}\in\mathbfcal{X}$ then

$\displaystyle\langle\bm{G}(\bm{f}(\bm{x})-\bm{f}(\hat{\bm{x}})),\bm{x}-\hat{\bm{x}}\rangle\leq\gamma_{s}\left\lVert\bm{x}-\hat{\bm{x}}\right\rVert_{2}^{2},$
for $\gamma_{s}\in\mathbb{R}$ and QIB in $\mathbfcal{X}$ if for any
$\bm{x},\hat{\bm{x}}\in\mathbfcal{X}$ it holds that
$\displaystyle\langle\bm{G}(\bm{f}(\bm{x})-\bm{f}(\hat{\bm{x}})),\bm{G}(\bm{f}(\bm{x})-\bm{f}(\hat{\bm{x}}))\rangle\leq$
$\displaystyle\qquad\gamma_{q1}\left\lVert\bm{x}-\hat{\bm{x}}\right\rVert_{2}^{2}+\gamma_{q2}\langle\bm{G}(\bm{f}(\bm{x})-\bm{f}(\hat{\bm{x}})),\bm{x}-\hat{\bm{x}}\rangle,$
for $\gamma_{q1},\gamma_{q2}\in\mathbb{R}$.

In Definition 4, the notation $\langle\cdot,\cdot\rangle$ denotes the standard
inner product. By considering OSL and QIB, the parameterization of NDSs
entails to finding the constants $\gamma_{s}$, $\gamma_{q1}$, and
$\gamma_{q2}$. These constants can be computed either analytically or
numerically via stochastic point-based method or deterministic interval-based
method. If numerical approach is pursued, then the computation of these
constants reduces to solving global optimization problems. Our particular work
in [24] investigates NDSs parameterization using interval-based global
maximization.

###  6.3 From OSL and QIB to Robust Selection

Now we demonstrate how the proposed methodology can be easily extend to solve
SASPs for other class of nonlinearities. Now suppose that $\bm{f}(\cdot)$ in
NDS (1) is both OSL and QIB with constants $\gamma_{s}$, $\gamma_{q1}$, and
$\gamma_{q2}$. Using observer design presented in [27], the corresponding SSP
can be formed as

$\displaystyle\;\operatorname*{minimize}_{\bm{P},\bm{Y},\epsilon_{1},\epsilon_{2},\bm{\gamma}}\;\;\bm{c}^{\top}\bm{\gamma}$
(21a) $\displaystyle\operatorname*{subject\
to}\;\;\begin{bmatrix}\bm{A}^{\top}\bm{P}+\bm{P}\bm{A}+\epsilon_{1}\gamma_{s}\bm{I}\\\
+\epsilon_{2}\gamma_{q1}\bm{I}-\bm{Y}\bm{\Gamma}\bm{C}-\bm{C}^{\top}\bm{\Gamma}\bm{Y}^{\top}&*\\\
\bm{G}^{\top}\bm{P}+\dfrac{\gamma_{q2}\epsilon_{2}-\epsilon_{1}}{2}\bm{I}&-\epsilon_{2}\bm{I}\end{bmatrix}\preceq
0$ (21b) $\displaystyle\quad\quad\quad\quad\;\;\bm{P}\succ
0,\;\epsilon_{1},\epsilon_{2}>0,\;\bm{\gamma}\in\mathcal{G}_{\gamma},\;\bm{\gamma}\in\\{0,1\\}^{N},\vspace*{-0.05cm}$
(21c)

while the SSP for discrete-time NDSs is formulated as [28]

$\displaystyle\;\operatorname*{minimize}_{\bm{P},\bm{Y},\epsilon_{1},\epsilon_{2},\bm{\gamma}}\;\;\bm{c}^{\top}\bm{\gamma}$
(22a) $\displaystyle\operatorname*{subject\ to}\;\;$ (22b)
$\displaystyle\begin{bmatrix}-\bm{P}+\epsilon_{1}\gamma_{s}\bm{I}+\epsilon_{2}\gamma_{l1}\bm{I}&*&*\\\
\bm{P}\bm{A}-\bm{Y}\bm{\Gamma}\bm{C}+\dfrac{\gamma_{q2}\epsilon_{2}-\epsilon_{1}}{2}\bm{I}&\bm{P}-\epsilon_{2}\bm{I}&*\\\
\bm{P}\bm{A}-\bm{Y}\bm{\Gamma}\bm{C}&\bm{O}&-\bm{P}\end{bmatrix}\preceq 0$
(22c) $\displaystyle\;\;\bm{P}\succ
0,\;\epsilon_{1},\epsilon_{2}>0,\;\bm{\gamma}\in\mathcal{G}_{\gamma},\;\bm{\gamma}\in\\{0,1\\}^{N}.\vspace*{-0.05cm}$
(22d)

Notice that the problem described in (21) possesses similar nonconvexity as in
P2. In addition to the above, the proposed sensor selection framework can also
facilitate to include robustness in the objective. For instance, by
considering the Lipschitz condition, robust SSP with $\mathcal{L}_{\infty}$
metric can be constructed as [20]

$\displaystyle\operatorname*{minimize}_{\bm{P},\bm{Y},\epsilon,\alpha,\mu_{0,1,2}}\quad\mu_{0}\mu_{1}+\mu_{2}$
(23a) $\displaystyle\operatorname*{subject\ to}$
$\displaystyle\begin{bmatrix}\bm{A}^{\top}\bm{P}+\bm{P}\bm{A}-\bm{C}^{\top}\bm{\Gamma}\bm{Y}^{\top}\\\
-\bm{Y}\bm{\Gamma}\bm{C}+\alpha\bm{P}+\epsilon\gamma_{l}^{2}\bm{I}&*&*\\\
\bm{P}&-\epsilon\bm{I}&*\\\
\bm{B_{\mathrm{w}}}^{\top}\bm{P}-\bm{D_{\mathrm{w}}}^{\top}\bm{Y}^{\top}&\bm{O}&-\alpha\mu_{0}\bm{I}\end{bmatrix}\preceq
0$ (23b) $\displaystyle\begin{bmatrix}-\bm{P}&*&*\\\
\bm{O}&-\mu_{2}\bm{I}&*\\\ \bm{Z}&\bm{O}&-\mu_{1}\bm{I}\end{bmatrix}\preceq
0,$ (23c) $\displaystyle\;\;\bm{P}\succ
0,\;\epsilon,\mu_{0},\mu_{1},\mu_{2}\geq
0,\;\alpha>0,\;\bm{\gamma}\in\mathcal{G}_{\gamma},\;\bm{\gamma}\in\\{0,1\\}^{N},\vspace*{-0.05cm}$
(23d)

where matrices $\bm{B_{\mathrm{w}}}$ and $\bm{D_{\mathrm{w}}}$ represent how
the disturbances are distributed, $\bm{Z}$ is the performance matrix with
$\bm{z}(t)=\bm{Z}\bm{e}(t)$, and $\bm{e}(t)$ is the estimation error. The
objective herein is to select the best sensor combination, given a fixed
number of activated sensors, that minimizes the worst case disturbance
attenuation $\mu=\sqrt{\mu_{0}\mu_{1}+\mu_{2}}$. If the problem described in
(23) is solved, then it is ensured that
$\lim_{t\rightarrow\infty}\sup\left\lVert\bm{z}(t)\right\rVert_{2}\leq\mu\left\lVert\bm{w}(t)\right\rVert_{\mathcal{L}_{\infty}}$
where $\bm{w}(t)$ denotes the disturbance vector [20]. By fixing some
constants, (21) and (23) share similar nonconvexity with P2 and thus the
proposed methodology can be used to solve the SSP. Provided that the variables
$\alpha$ and $\mu_{0}$ or $\mu_{1}$ are fixed, this problem along with the
ones described in (21) and (22) share similar nonconvexity as P2 and
therefore, by using the proposed methodology, they can be solved using the SE-
BnB algorithm. Next, we focus on numerically testing the SE-BnB algorithm in
solving SSP for a NDS satisfying the Lipschitz condition.

(a)  (b)

Figure 1: Numerical experiment results for network of nonlinear unstable nodes
to find suboptimal solutions via the SE-BnB and S-BnB: (a) computational time
for different number of maximum allowable branching $\mathrm{mBr}$ and (b) the
resulting sensor costs. The solid lines correspond to SE-BnB whereas dashed
lines correspond to S-BnB. To compensate the heuristics on both SE-BnB and
S-BnB, the numerical test is performed three times for each value of $N$ and
the results shown correspond to the average values.

(a)  (b)

Figure 2: Comparison between SE-BnB and $\ell_{1}$-norm relaxation to search
for suboptimal solutions: (a) computational time and (b) the resulting sensor
costs. To compensate the heuristics on SE-BnB, the numerical test is performed
three times and the results shown correspond to the average values.

##  7 Numerical Assessment

In this section, we demonstrate the proposed methodology to solve SSP for a
NDS on network of nonlinear unstable nodes. The simulations are performed
using MATLAB R2019a, where YALMIP’s [17] optimization package with MOSEK [29]
solver are used to solve the corresponding SDPs problem. The network of
nonlinear unstable nodes considered here is adapted from [30], which is
initially a network of linear systems. In order to obtain nonlinear systems, a
sinusoidal function is introduced on each node such that each node $i$ has the
following dynamics

$\displaystyle\dot{\bm{x}}_{i}=\begin{bmatrix}\zeta_{1i}&1\\\
1&\zeta_{2i}\end{bmatrix}\bm{x}_{i}+\begin{bmatrix}0\\\
\beta_{i}\end{bmatrix}\sin(x_{2i})+\sum_{j\neq
i}e^{\alpha(i,j)}\bm{x}_{j}+\begin{bmatrix}0\\\ 1\end{bmatrix}u_{i},$ (24)

where the constants $\zeta_{1i}$, $\zeta_{2i}$, and $\beta_{i}$ are randomly
generated within $[-2,2]$, $[-2,2]$, and $[-1,1]$ respectively for each $i$.
The coupling between nodes $i$ and $j$ is determined by the function of
Euclidean distance $\alpha(i,j)$. This distance is randomly generated inside a
box of size $5\times 5$. It can be verified that the Lipschitz constant for
this system is $\gamma_{l}=\sqrt{N}$ where $N$ is the number of nodes. For the
sake of simplicity, the weighting vector is set to be $\bm{c}=\bm{1}$ where
$\bm{Y}$ is bounded such that $-10^{3}\leq Y_{i,j}\leq 10^{3}$ for each $i,j$.
From (24), each node $i$ can have at most two measurements—hence for $N$
subsystems, there are at most $2N$ measurements.

###  7.1 SE-BnB and Standard BnB

In the first case of this numerical assessment, we evaluate the performance of
the SE-BnB algorithm to find suboptimal solutions for SSP. The tolerance
defining optimality gap is set to be $10^{-4}$ with integer tolerance
$10^{-6}$ while restricting the maximum number of allowable branching. We also
invoke a logistic constraint such that at least $20\%$ of the available
measurements are equipped with sensor. We test the new SE-BnB algorithm with
various number of nodes and compare the results with different values of
$\mathrm{maxBranch}$ (or $\mathrm{mBr}$) against our own implementation of the
standard BnB algorithm (S-BnB)—that is, only takes MISDP P4—while still
utilizes the heuristic discussed in Section 5.2. The results of this numerical
test are illustrated in Fig. 1. It can be seen that, in overall, the SE-BnB
requires less computational time than the S-BnB. This is due to the fact that,
in the SE-BnB, the corresponding SDPs solved to obtain local lower bound and
global upper bound on each BnB node have reduced complexities as opposed to
the S-BnB. It is also showcased in Fig. 1b that this reduction in
computational time does not contribute to the resulting sensor costs since
sensor costs are merely determined by the number of explored nodes as well as
the heuristics used to find global upper bounds. Both SE-BnB and S-BnB show
similar sensor costs. These costs are reduced, despite of some irregularities,
as more branches are explored—these irregularities are attributed to the
heuristics used in finding global upper bounds.

###  7.2 SE-BnB and $\ell_{1}$-Norm Relaxation

Next, we assess the ability of our SE-BnB algorithm in finding good suboptimal
solutions relative to $\ell_{1}$-norm relaxation, which is developed in [18],
to solve a class of combinatorial optimization problems. This technique finds
its popularity in SASPs for linear systems—for example, see [31]. This
particular approach works as follows. First, note that each $\gamma_{i}$
determines the values of $\mathrm{col}\left(\bm{Y}\right)_{i}$, i.e.,
$\mathrm{col}\left(\bm{Y}\right)_{i}$ can be any value if and only if
$\gamma_{i}=1$ and zero otherwise. As such, the integer variable $\bm{\gamma}$
can be removed from optimization variables and consequently, the SSP boils
down into minimizing the number of nonzero columns of $\bm{Y}$, i.e.,
$\sum_{i=1}^{n_{y}}\left\lVert\mathrm{col}\left(\bm{Y}\right)_{i}\right\rVert_{\ell_{0}}$.
Since $\ell_{0}$-norm is nonconvex, it is replaced by $\ell_{1}$-norm,
yielding the following relaxed problem

$\displaystyle(\mathbf{P7})\;\operatorname*{minimize}_{\bm{P},\bm{Y},\epsilon,\bm{\gamma}}\;\;$
$\displaystyle\sum_{i=1}^{n_{y}}c_{i}w_{i}^{(k)}\left\lVert\mathrm{col}\left(\bm{Y}\right)_{i}\right\rVert_{\ell_{1}}$
(25a) $\displaystyle\operatorname*{subject\ to}\;\;$
$\displaystyle\eqref{eq:LMI_lip_obs},\;\bm{P}\succ
0,\;\epsilon>0,\vspace*{-0.25cm}$ (25b)

where $w_{i}^{(k)}$ is the corresponding weighting factor for column $i$ at
iteration $k$. When $\bm{w}^{(k)}$ is fixed, P7 becomes a convex SDP. This
weighting factor is updated at each iteration $k$ using the solution from
previous iteration through the following update rule [18]

$\displaystyle
w_{i}^{(k+1)}:=\frac{1}{\left\lVert\mathrm{col}\left(\bm{Y}\right)_{i}^{(k)}\right\rVert_{\ell_{1}}+\kappa},\;\forall
i=1,2,\ldots,n_{y},$ (26)

where $\kappa>0$ is a relatively small constant. Algorithm 3, adapted from
[18], provides the detailed steps.

In this numerical experiment, both algorithms are configured such that the
objective is to find at least $N-1$ activated sensors. The $\ell_{1}$-norm
relaxation is used with $\epsilon_{w}=0.1$ and $500$ maximum iterations. Two
different values for $\kappa$ are considered which correspond to two scenarios
for this approach: $\ell_{1}$-norm relaxation-1 uses $\kappa=10^{-4}$ while
$\ell_{1}$-norm relaxation-2 uses $\kappa=10^{-5}$. The corresponding results
are illustrated in Fig. 2. In particular, Fig. 2b indicates that our SE-BnB
algorithm is able to find suboptimal sensor combinations with $N-1$ activated
sensors while the $\ell_{1}$-norm relaxation fails to find sensor combinations
under such prescribed limitation. On the matter of computational time, as
shown in Fig. 2a, the SE-BnB algorithm requires less computational time
compared to the $\ell_{1}$-norm relaxation. This demonstrates an advantage of
the proposed SE-BnB algorithm in finding suboptimal solutions for SSP over the
$\ell_{1}$-norm relaxation technique, at least for the network of nonlinear
unstable nodes described in (24).

1 input: $\epsilon_{w}$, $\kappa$, $\mathrm{maxIter}$

2 initialize: $k=0$, $\bm{w}^{(1)}=\bm{1}$

3 do

4 $k\leftarrow k+1$

5 solve P7 in (25), obtain $\bm{Y}^{(k)}$

6 foreach  _$i\in\mathbb{I}(n_{y})$_ do

7 update $w_{i}^{(k+1)}$ from $\bm{Y}^{(k)}$ using (26)

8

9

10while
_$\left\lvert\left\lVert\bm{w}^{(k+1)}\right\rVert_{2}-\left\lVert\bm{w}^{(k)}\right\rVert_{2}\right\rvert
>\epsilon_{w}$ and $k<\mathrm{maxIter}$ _

11find: a feasible solution of LMI (25b) where $\bm{Y}$ corresponds to the
nonzero column of $\bm{Y}^{(k)}$ from the last iteration

12 output: $\bm{P}$, $\bm{Y}$

Algorithm 3 Iterative $\ell_{1}$-Norm Relaxation for SSP

###  7.3 Corroborating the Theory

Finally, we employ SE-BnB algorithm to find optimal solutions for SSP while
using different values of Lipschitz constant. This numerical test is intended
to find empirical relation between the value of Lipschitz constant and the
number of activated sensors. For this purpose, we consider the so called
Lipschitz multiplier $\eta>0$ such that the Lipschitz constant for $N$ number
of subsystems is formulated as $\gamma_{l}=\eta\sqrt{N}$. Fig. 3 illustrates
the corresponding results. It can be observed that increasing Lipschitz
constant (that is, by equivalently increasing $\eta$) raises the number of
activated sensors—this is indicated by the nondecreasing fashion in the number
of sensor costs as the value of $\eta$ increases. This result corroborates the
theoretical foundation given in Section 4 and Theorem 1. That is, in the
context of SSP, increasing Lipschitz constant shrinks the feasible space of
the SSP problem, hence requiring more sensors to stabilize the estimation
error dynamics.

##  8 Paper Summary and Limitations

A novel general framework for dealing with SASPs for NDSs is presented. Our
approach is built upon SDP formulations for observer/controller designs for
various classes of NDSs developed in the literature. Specifically, this paper
focuses on addressing SASPs for Lipschitz NDSs and our investigations show
that smaller Lipschitz constant allows less number SAs are required to achieve
stabilization in control/estimation purpose. A customized BnB algorithm,
referred to as SE-BnB, which exploits problem structure and utilizing new
heuristics to efficiently obtain optimal and suboptimal solutions for SASPs,
is proposed. The main advantage of our framework as opposed to other
approaches from the literature is its capability to obtain the optimal SAs
combination by means of a BnB algorithm for stable/unstable NDS that ensures
stability for either estimation and/or stabilization/tracking purposes.

The algorithms presented in this paper come with their limitations. First,
this paper is not concerned with scaling the problem for extremely large-scale
systems. This is a result of using SDPs and LMIs that still have serious
scalability issues—in addition to the worst case complexity of BnB routines.
Second, this paper also does not address the simultaneous selection of SAs and
also assumes that all SAs are operating in ideal conditions (e.g. no
saturation effect). To that end, our future work will focus on (i)
investigating different approaches to solve SASPs other than using BnB
algorithm such as the generalized Benders decomposition [32] and branch-and-
cut algorithm [33], (ii) extending the proposed approach for addressing the
robust SASPs as well as the simultaneous selection of SAs through output
feedback control and observer-based control policies, and (iii) considering
the effect of SAs saturation in SASPs.

Figure 3: Optimal sensor costs computed through the SE-BnB for different
Lipschitz multiplier $\eta$ where the Lipschitz constant is given by
$\gamma_{l}=\eta\sqrt{N}$.

## Acknowledgments

This material is based upon work supported by the National Science Foundation
(NSF) under Grants 1728629, 1917164, 2151571, and 2152450. We gratefully
acknowledge NSF’s support and the thorough comments and suggestions made by
the editor and the reviewers.

## References

  * [1] J. Qi, K. Sun, W. Kang, Optimal pmu placement for power system dynamic state estimation by using empirical observability gramian, IEEE Transactions on power Systems 30 (4) (2015) 2041–2054.
  * [2] S. Pequito, N. Popli, S. Kar, M. D. Ilić, A. P. Aguiar, A framework for actuator placement in large scale power systems: Minimal strong structural controllability, in: 2013 5th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 2013, pp. 416–419.
  * [3] T. H. Summers, J. Lygeros, Optimal sensor and actuator placement in complex dynamical networks, IFAC Proceedings Volumes 47 (3) (2014) 3784 – 3789, 19th IFAC World Congress.  doi:https://doi.org/10.3182/20140824-6-ZA-1003.00226.
  * [4] N. K. Dhingra, M. R. Jovanović, Z.-Q. Luo, An admm algorithm for optimal sensor and actuator selection, in: IEEE Conference on Decision and Control, IEEE, 2014, pp. 4039–4044.
  * [5] A. F. Taha, N. Gatsis, T. Summers, S. A. Nugroho, Time-varying sensor and actuator selection for uncertain cyber-physical systems, IEEE Transactions on Control of Network Systems 6 (2) (2019) 750–762.  doi:10.1109/TCNS.2018.2873229.
  * [6] S. A. Nugroho, A. F. Taha, N. Gatsis, T. H. Summers, R. Krishnan, Algorithms for joint sensor and control nodes selection in dynamic networks, Automatica 106 (2019) 124 – 133.  doi:https://doi.org/10.1016/j.automatica.2019.04.047.
  * [7] S. A. Nugroho, A. F. Taha, On the need for sensor and actuator placement algorithms in nonlinear systems: Wip abstract, in: Proceedings of the 10th ACM/IEEE International Conference on Cyber-Physical Systems, ICCPS ’19, ACM, New York, NY, USA, 2019, pp. 304–305.  doi:10.1145/3302509.3313313.
  * [8] A. Haber, F. Molnar, A. E. Motter, State observation and sensor selection for nonlinear networks, IEEE Transactions on Control of Network Systems 5 (2) (2018) 694–708.  doi:10.1109/TCNS.2017.2728201.
  * [9] S. D. Bopardikar, O. Ennasr, X. Tan, Randomized sensor selection for nonlinear systems with application to target localization, IEEE Robotics and Automation Letters 4 (4) (2019) 3553–3560.  doi:10.1109/LRA.2019.2928208.
  * [10] A. Haber, S. A. Nugroho, P. Torres, A. F. Taha, Control node selection algorithm for nonlinear dynamic networks, IEEE Control Systems Letters 5 (4) (2021) 1195–1200.
  * [11] A. Haber, Joint sensor node selection and state estimation for nonlinear networks and systems, IEEE Transactions on Network Science and Engineering (2021) 1–1doi:10.1109/TNSE.2021.3069890.
  * [12] G. Phanomchoeng, R. Rajamani, The bounded jacobian approach to nonlinear observer design, in: Proceedings of the 2010 American Control Conference, 2010, pp. 6083–6088.  doi:10.1109/ACC.2010.5531543.
  * [13] G. Phanomchoeng, R. Rajamani, Observer design for lipschitz nonlinear systems using riccati equations, in: Proceedings of the 2010 American Control Conference, 2010, pp. 6060–6065.  doi:10.1109/ACC.2010.5531294.
  * [14] M. Abbaszadeh, H. J. Marquez, Nonlinear observer design for one-sided lipschitz systems, in: Proceedings of the 2010 American Control Conference, 2010, pp. 5284–5289.  doi:10.1109/ACC.2010.5530715.
  * [15] Y. Guo, Y. Chen, C. Zhang, Decentralized state-observer-based traffic density estimation of large-scale urban freeway network by dynamic model, Information 8 (3) (2017) 95.
  * [16] S. A. Nugroho, A. F. Taha, Sensor placement strategies for some classes of nonlinear dynamic systems via lyapunov theory, in: 2019 Conference on Decision and Control (CDC), 2019.
  * [17] J. Löfberg, YALMIP: A toolbox for modeling and optimization in matlab, in: Proc. IEEE Int. Symp. Computer Aided Control Syst. Design, 2004, pp. 284–289.
  * [18] E. J. Candès, M. B. Wakin, S. P. Boyd, Enhancing sparsity by reweighted $\ell-1$ minimization, Journal of Fourier Analysis and Applications 14 (5) (2008) 877–905.  doi:10.1007/s00041-008-9045-x.
  * [19] A. Bergen, V. Vittal, Power Systems Analysis, Prentice Hall, 2000.
  * [20] S. A. Nugroho, A. F. Taha, C. G. Claudel, A control-theoretic approach for scalable and robust traffic density estimation using convex optimization, IEEE Transactions on Intelligent Transportation Systems 22 (1) (2021) 64–78.  doi:10.1109/TITS.2019.2953023.
  * [21] F. Perini, E. Galligani, R. D. Reitz, An analytical jacobian approach to sparse reaction kinetics for computationally efficient combustion modeling with large reaction mechanisms, Energy & Fuels 26 (8) (2012) 4804–4822.  doi:10.1021/ef300747n.
  * [22] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, Lect. Notes Phys. 39 (1975) 420–422.  doi:10.1007/BFb0013365.
  * [23] A. Lajmanovich, J. A. Yorke, A deterministic model for gonorrhea in a nonhomogeneous population, Mathematical Biosciences 28 (3) (1976) 221 – 236.
  * [24] S. A. Nugroho, A. F. Taha, , V. Hoang, Nonlinear dynamic systems parameterization using interval-based global optimization: Computing lipschitz constants and beyond (2020).  arXiv:2004.12061.
  * [25] G. P. McCormick, Computability of global solutions to factorable nonconvex programs: Part i—convex underestimating problems, Mathematical programming 10 (1) (1976) 147–175.
  * [26] M. Yadegar, A. Afshar, M. Davoodi, Observer-based tracking controller design for a class of lipschitz nonlinear systems, Journal of Vibration and Control 24 (11) (2018) 2112–2119.  doi:10.1177/1077546317721597.
  * [27] W. Zhang, H. Su, H. Wang, Z. Han, Full-order and reduced-order observers for one-sided lipschitz nonlinear systems using riccati equations, Communications in Nonlinear Science and Numerical Simulation 17 (12) (2012) 4968–4977.
  * [28] W. Zhang, H. Su, F. Zhu, D. Yue, A note on observers for discrete-time lipschitz nonlinear systems, IEEE Transactions on Circuits and Systems II: Express Briefs 59 (2) (2012) 123–127.  doi:10.1109/TCSII.2011.2174671.
  * [29] E. D. Andersen, K. D. Andersen, The Mosek Interior Point Optimizer for Linear Programming: An Implementation of the Homogeneous Algorithm, Springer US, Boston, MA, 2000, pp. 197–232.
  * [30] N. Motee, A. Jadbabaie, Optimal control of spatially distributed systems, IEEE Transactions on Automatic Control 53 (7) (2008) 1616–1629.  doi:10.1109/TAC.2008.929366.
  * [31] A. Argha, S. W. Su, A. Savkin, B. Celler, A framework for optimal actuator/sensor selection in a control system, International Journal of Control 92 (2) (2019) 242–260.  doi:10.1080/00207179.2017.1350755.
  * [32] J. Zhang, X. Wang, D. J. Chmielewski, Covariance-based hardware selection part iv: Solution using the generalized benders decomposition, AIChE Journal 62 (10) (2016) 3628–3638.  doi:10.1002/aic.15285.
  * [33] K. Kobayashi, Y. Takano, A branch-and-cut algorithm for solving mixed-integer semidefinite optimization problems, Computational Optimization and Applications 75 (2) (2020) 493–513.
  * [34] V. Balakrishnan, L. Vandenberghe, Semidefinite programming duality and linear time-invariant systems, IEEE Transactions on Automatic Control 48 (1) (2003) 30–41.

##  Appendix A Adjoint Mapping for Symmetric Matrices

Consider $T:\mathcal{V}\rightarrow\mathcal{W}$ where $\mathcal{V}$ and
$\mathcal{W}$ are finite dimensional vector spaces with the associated inner
product given as $\langle\cdot,\cdot\rangle_{\mathcal{V}}$ and
$\langle\cdot,\cdot\rangle_{\mathcal{W}}$, respectively. For every
$\bm{w}\in\mathcal{W}$, the adjoint of $T$, which is denoted by $T^{\star}$,
is a mapping such that $\langle
T(\bm{v}),\bm{w}\rangle_{\mathcal{V}}=\langle\bm{v},T^{\star}(\bm{w})\rangle_{\mathcal{W}}$
for a unique $\bm{v}\in\mathcal{V}$ that satisfies $T(\bm{v})=\bm{w}$. For
LMIs represented by linear mapping
$\mathcal{A}:\mathcal{V}\rightarrow\mathbb{S}^{n}$, its adjoint
$\mathcal{A}^{\star}:\mathbb{S}^{n}\rightarrow\mathcal{V}$ is therefore a
mapping such that for every $\bm{v}\in\mathcal{V}$ and
$\bm{S}\in\mathbb{S}^{n}$, it holds that
$\langle\mathcal{A}(\bm{v}),\bm{S}\rangle_{\mathcal{V}}=\langle\bm{v},\mathcal{A}^{\star}(\bm{S})\rangle_{\mathbb{S}^{n}}$;
see reference [34].

##  Appendix B Proof of Theorem 1

The proof relies on the notion of adjoint mapping for both finite-dimensional
vector spaces and symmetric matrices, as discussed in A. Let
$\mathcal{A}:\mathbb{S}^{n_{x}}\times\mathbb{R}^{n_{x}n_{y}}\rightarrow\mathbb{S}^{2n_{x}}\times\mathbb{S}^{n_{x}}$
be expressed as
$\mathcal{A}(\bm{P},\bm{y}_{v})=\mathcal{A}_{P}(\bm{P})+\mathcal{A}_{Y}({\bm{y}_{v}})$
for $\bm{P}\in\mathbb{S}^{n_{x}}$ and $\bm{y}_{v}\in\mathbb{R}^{n_{x}n_{y}}$.
First, it can be proven that for $\mathcal{A}_{P}(\bm{P})$ described in (6b)
and matrix $\bm{Z}$ given in (7a),
$\langle\mathcal{A}_{P}(\bm{P}),\bm{Z}\rangle_{\mathrm{dom}(\mathcal{A})}$ is
equal to

$\displaystyle\langle\bm{P},-\bm{Z}_{1}\bm{A}^{\top}-\bm{A}\bm{Z}_{1}-\bm{G}\bm{Z}_{2}^{\top}-\bm{Z}_{2}\bm{G}^{\top}+\bm{Z}_{4}\rangle_{\mathrm{im}(\mathcal{A})},$
yielding
$\mathcal{A}_{P}^{\star}(\bm{Z})=-\bm{Z}_{1}\bm{A}^{\top}-\bm{A}\bm{Z}_{1}-\bm{G}\bm{Z}_{2}^{\top}-\bm{Z}_{2}\bm{G}^{\top}+\bm{Z}_{4}$.
Second, let us reformulate $\mathcal{A}_{Y}({\bm{y}_{v}})$ described in (6c)
into $\displaystyle\mathcal{A}_{Y}({\bm{y}_{v}})$
$\displaystyle=\sum_{i=1}^{n_{x}n_{y}}y_{v,i}\bar{\bm{C}}_{D,i},$ (27a) where
$\bar{\bm{C}}_{D,i}$ in (27a) is nothing but
$\displaystyle\bar{\bm{C}}_{D,i}=\mathrm{Blkdiag}\left(\begin{bmatrix}\bar{\bm{C}}_{i}&*\\\
\bm{O}&\bm{O}\end{bmatrix},\bm{O}\right).$ Then, it is straightforward to show
that
$\displaystyle\langle\mathcal{A}_{Y}({\bm{y}_{v}}),\bm{Z}\rangle_{\mathrm{dom}(\mathcal{A})}$
$\displaystyle\;\;\;=\langle{\bm{y}_{v}},\begin{bmatrix}\mathrm{tr}(\bar{\bm{C}}_{D,1}\bm{Z})\,\,\mathrm{tr}(\bar{\bm{C}}_{D,2}\bm{Z})\,\,\cdots\,\,\mathrm{tr}(\bar{\bm{C}}_{D,{n_{x}n_{y}}}\bm{Z})\end{bmatrix}\rangle_{\mathrm{im}(\mathcal{A})}$
$\displaystyle\;\;\;=\langle{\bm{y}_{v}},\begin{bmatrix}\mathrm{tr}(\bar{\bm{C}}_{1}\bm{Z}_{1})\,\,\mathrm{tr}(\bar{\bm{C}}_{2}\bm{Z}_{1})\,\,\cdots\,\,\mathrm{tr}(\bar{\bm{C}}_{n_{x}n_{y}}\bm{Z}_{1})\end{bmatrix}\rangle_{\mathrm{im}(\mathcal{A})}.$
The above indicates that
$\mathcal{A}_{Y}^{\star}({\bm{Z}})=\left\\{\begin{bmatrix}\mathrm{tr}(\bar{\bm{C}}_{i}\bm{Z}_{1})\end{bmatrix}\right\\}^{n_{x}n_{y}}_{i=1}$.
From these results, it can be confirmed that the adjoint mapping of
$\mathcal{A}(\bm{P},\bm{y}_{v})$ is equivalent to
$\mathcal{A}^{\star}(\bm{Z})=\left(\mathcal{A}_{P}^{\star}(\bm{Z}),\mathcal{A}_{Y}^{\star}(\bm{Z})\right)$.
The theorem of alternatives [34, Theorem 1] declares that whether there exist
$\bm{P}\in\mathbb{S}^{n_{x}}$ and $\bm{y}_{v}\in\mathbb{R}^{n_{x}n_{y}}$ such
that $\mathcal{A}(\bm{P},\bm{y}_{v})+\mathcal{A}_{0}\succ 0$ or there exists
$\bm{Z}$ as given in (7a) such that $\bm{Z}\succeq 0$ and $\bm{Z}\neq 0$
satisfying $\mathcal{A}^{\star}(\bm{Z})=0$ and
$\mathrm{tr}(\mathcal{A}_{0}\bm{Z})\leq 0$. Notice that
$\mathcal{A}^{\star}(\bm{Z})=0$ is indeed equivalent to
$\displaystyle\mathcal{A}^{\star}(\bm{Z})=0\Leftrightarrow\begin{cases}\mathcal{A}_{P}^{\star}(\bm{Z})=0\\\
\mathcal{A}_{Y}^{\star}(\bm{Z})=0\end{cases}.$ (27b) The first right-hand side
of (27b) implies
$\displaystyle-\bm{Z}_{1}\bm{A}^{\top}-\bm{A}\bm{Z}_{1}-\bm{G}\bm{Z}_{2}^{\top}-\bm{Z}_{2}\bm{G}^{\top}+\bm{Z}_{4}=0,$
(27c) whereas the second condition implies that
$\displaystyle\mathrm{tr}(\bar{\bm{C}}_{i}\bm{Z}_{1})=\sum_{j=1}^{n_{x}}\sum_{k=1}^{n_{x}}\bar{C}_{i,jk}Z_{1,jk}=0,\;\forall
i\in\mathbb{I}(n_{x}n_{y}).$ (27d) The results from (27c) and (27d) establish
(7b) and (7c), respectively. Next, note that
$\mathrm{tr}(\mathcal{A}_{0}\bm{Z})\leq 0$ implies
$\displaystyle\mathrm{tr}(\mathcal{A}_{0}\bm{Z})=\mathrm{tr}(-\epsilon\gamma_{l}^{2}\bm{Z}_{1})+\mathrm{tr}(\epsilon\bm{Z}_{3})\leq
0.$ (27e) Since $\epsilon>0$, then (27e) establishes (7d). This concludes the
proof.

$\blacksquare$

##  Appendix C Proof of Theorem 2

The proof of the proposition is similar to that of Theorem 1. The only
difference lies in the adjoint of $\mathcal{A}_{Y,\Gamma}({\bm{y}_{v}})$. From
(12a), it is not difficult to show that

$\displaystyle\langle\mathcal{A}_{Y}({\bm{y}_{v}}),\bm{Z}\rangle_{\mathrm{dom}(\mathcal{A})}$
$\displaystyle\;\;\;=\langle{\bm{y}_{v}},\begin{bmatrix}\mathrm{tr}(\bar{\gamma}_{1}\bar{\bm{C}}_{1}\bm{Z}_{1})\,\,\cdots\,\,\mathrm{tr}(\bar{\gamma}_{n_{x}n_{y}}\bar{\bm{C}}_{n_{x}n_{y}}\bm{Z}_{1})\end{bmatrix}\rangle_{\mathrm{im}(\mathcal{A})}.$

Since $\bar{\gamma}_{i}=0$ where $\bar{\gamma}_{i}\triangleleft{\gamma}_{j}$
for which ${\gamma}_{j}\notin\mathcal{S}_{\gamma}^{(a)}$, then

$\displaystyle\mathcal{A}_{Y,\Gamma}^{\star}({\bm{Z}})=\left\\{\begin{bmatrix}\mathrm{tr}(\bar{\gamma}_{i}\bar{\bm{C}}_{i}\bm{Z}_{1})\end{bmatrix}\right\\}^{n_{x}n_{y}}_{i=1},\;\forall\bar{\gamma}_{i}\triangleleft{\gamma}_{j},\;\bar{\gamma}_{j}\in\mathcal{S}_{\gamma}^{(a)},$

which is equivalent to having
$\bm{Z}_{1}\in\mathcal{Z}\big{(}\mathcal{S}_{\gamma}^{(a)}\big{)}$ where
$\mathcal{Z}\big{(}\mathcal{S}_{\gamma}^{(a)}\big{)}$ is defined in (13), thus
completing the proof.  $\blacksquare$

##  Appendix D Proof of Theorem 3

The equality constraint $\bm{M}=\bm{Y}\bm{\Gamma}$ in P3 is equivalent to

$\displaystyle M_{ij}=Y_{ij}\gamma_{j}\Leftrightarrow
M_{ij}=\begin{cases}Y_{ij},&\text{if}\;\;\gamma_{j}=1\\\
0,&\text{if}\;\;\gamma_{j}=0,\end{cases}\vspace{-0.0cm}$ (28)

where
$Y_{ij}\in\left[\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij},\bar{Y}_{ij}\right]$
for all $i,j$ for McCormick’s reformulation and $Y_{ij}\in\left[-L,L\right]$
for all $i,j$ for big-M. In what follows, we apply McCormick’s relaxation to
transform the consequence in (28) into convex MISDP as it generalizes the
big-M method. To that end, by realizing that $\bar{Y}_{ij}-Y_{ij}$,
$Y_{ij}-\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}$,
$1-\gamma_{j}$, and $\gamma_{j}$ are all nonnegative and noting that
$M_{ij}=Y_{ij}\gamma_{j}$, we get

$\displaystyle\left(\bar{Y}_{ij}-Y_{ij}\right)\left(1-\gamma_{j}\right)$
$\displaystyle\geq 0$ $\displaystyle\Leftrightarrow M_{ij}$ $\displaystyle\geq
Y_{ij}+\bar{Y}_{ij}\left(\gamma_{j}-1\right)$ (29a)
$\displaystyle\left(\bar{Y}_{ij}-Y_{ij}\right)\gamma_{j}$ $\displaystyle\geq
0\Leftrightarrow\bar{Y}_{ij}\gamma_{j}\geq M_{ij}$ (29b)
$\displaystyle\left(Y_{ij}-\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\right)\left(1-\gamma_{j}\right)$
$\displaystyle\geq 0$ $\displaystyle\Leftrightarrow-M_{ij}$
$\displaystyle\geq-
Y_{ij}+\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\left(1-\gamma_{j}\right)$
(29c)
$\displaystyle\left(Y_{ij}-\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\right)\gamma_{j}$
$\displaystyle\geq 0\Leftrightarrow
M_{ij}\geq\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\gamma_{j}.$
(29d)

Next, notice that there are only two possible values for $\gamma_{j}$. That
is, for $\gamma_{j}=1$, then substituting $\gamma_{j}=1$ to (29) yields

$\left.\begin{aligned} Y_{ij}\geq&\,M_{ij}\geq Y_{ij}\\\
\bar{Y}_{ij}\geq&\,M_{ij}\geq\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\end{aligned}\right\\}\Rightarrow
M_{ij}=Y_{ij}.$

On the other hand, for $\gamma_{j}=0$, substituting $\gamma_{j}=0$ to (29)
yields

$\left.\begin{aligned} 0\geq&\,M_{ij}\geq 0\\\
\bar{Y}_{ij}\geq&\,Y_{ij}\geq\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\end{aligned}\right\\}\Rightarrow
M_{ij}=0.$

Since this holds for all $i,j$, then (28) and (29) are equivalent. For big-M
method, one can obtain a similar equivalency by substituting
$\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}=-L$ and
$\bar{Y}_{ij}=L$ in (29). Next, define $\bm{\sigma}_{1}$, $\bm{\sigma}_{2}$,
$\bm{\Omega}^{\prime}$ and $\bm{\Omega}$ as follows

$\displaystyle\bm{\sigma}_{1}$
$\displaystyle:=\begin{bmatrix}1&-1&1&-1\end{bmatrix}^{\top},\;$
$\displaystyle\bm{\sigma}_{2}$
$\displaystyle:=\begin{bmatrix}-1&1&0&0\end{bmatrix}^{\top},$
$\displaystyle\bm{\Omega}^{\prime}$
$\displaystyle:=\underbrace{\begin{bmatrix}\mathrm{vec}(\bm{\Omega})&\ldots&\mathrm{vec}(\bm{\Omega})\end{bmatrix}}_{n_{y}\,\text{times}},$
$\displaystyle\bm{\Omega}$
$\displaystyle:=\begin{bmatrix}\bm{\omega}_{11}&\bm{\omega}_{12}&\cdots&\bm{\omega}_{1n_{y}}\\\
\bm{\omega}_{21}&\bm{\omega}_{22}&\cdots&\bm{\omega}_{2n_{y}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{\omega}_{n_{x}1}&\bm{\omega}_{n_{x}2}&\cdots&\bm{\omega}_{n_{x}n_{y}}\end{bmatrix},\,$

where $\bm{\omega}_{ij}$ is given by

$\displaystyle\bm{\omega}_{ij}:=\begin{cases}\begin{bmatrix}-\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}&\bar{Y}_{ij}&-\bar{Y}_{ij}&\stackunder[1.2pt]{$Y$}{\rule{3.44444pt}{0.32289pt}}_{ij}\end{bmatrix}^{\top},&\;\text{for
McCormick},\\\ \begin{bmatrix}L&L&-L&-L\end{bmatrix}^{\top},&\;\text{for
Big-M}.\\\ \end{cases}$

By using the above notations, (29) can be written as
$\bm{\Xi}\,\bm{\nu}\leq\bm{\varphi}$ where where matrix
$\bm{\Xi}\in\mathbb{R}^{(2n_{x}n_{y}+n_{y})\times(4n_{x}n_{y})}$ is given as

$\displaystyle\bm{\Xi}:=\begin{bmatrix}\bm{I}_{n_{x}n_{y}}\otimes\bm{\sigma}_{1}\;\vline\;\bm{I}_{n_{x}n_{y}}\otimes\bm{\sigma}_{2}\;\vline\;\bm{\Omega}^{\prime}\odot\bm{I}_{n_{y}}\otimes\bm{I}_{4n_{x}}\end{bmatrix},$
(30a) vector $\bm{\nu}\in\mathbb{R}^{2n_{x}n_{y}+n_{y}}$ is given as
$\displaystyle\bm{\nu}:=\begin{bmatrix}\mathrm{vec}(\bm{M})^{\top}\;\;\mathrm{vec}(\bm{Y})^{\top}\;\;\bm{\gamma}^{\top}\end{bmatrix}^{\top},$
(30b) and vector $\varphi\in\mathbb{R}^{4n_{x}n_{y}}$ is given as
$\displaystyle\bm{\varphi}:=\mathrm{vec}(\bm{\Omega})\odot\bm{I}_{n_{x}n_{y}}\otimes\begin{bmatrix}1&1&0&0\end{bmatrix}^{\top}.$
(30c)

This yields the following problem

$\displaystyle\;\operatorname*{minimize}_{\bm{P},\bm{Y},\bm{M},\epsilon,\bm{\gamma}}\;\;$
$\displaystyle\bm{c}^{\top}\bm{\gamma}$ (31a)
$\displaystyle\operatorname*{subject\ to}\;\;$
$\displaystyle\eqref{eq:sen_sel_obs_const_1},\,\eqref{eq:sen_sel_obs_2},\,\bm{Y}\in
B(\bm{Y}),\,\bm{\Xi}\,\bm{\nu}\leq\bm{\varphi},$ (31b)

Since other constraints except $\bm{M}=\bm{Y}\bm{\Gamma}$ in P3 are convex,
then P3 and the problem described in (31) are equivalent. The next step is to
transform problem (31) into standard inequality form. To do so, notice that
the logistic constraint in (10c) can be implicitly represented by
$\bm{H}^{\top}\bm{\gamma}\leq\bm{h}$ in which
$\bm{H}\in\mathbb{R}^{n_{y}\times n_{h}}$ and $\bm{h}\in\mathbb{R}^{n_{h}}$
whereas $\bm{Y}\in B(\bm{Y})$ can be represented by
$\bm{\Phi}^{\top}\mathrm{vec}(\bm{Y})\leq\bm{\psi}$ where
$\bm{\Phi}=\begin{bmatrix}\bm{I}&-\bm{I}\end{bmatrix}\in\mathbb{R}^{n_{x}n_{y}\times
2n_{x}n_{y}}$ and
$\bm{\psi}=\begin{bmatrix}\mathrm{vec}\left(\bar{\bm{Y}}\right)^{\top}&\mathrm{vec}\left(\stackunder[1.2pt]{$\bm{Y}$}{\rule{3.44444pt}{0.32289pt}}\right)^{\top}\end{bmatrix}$
for McCormick’s relaxation—similarly for big-M. Next, define
$\mathbb{I}_{v}(\bm{P})$, $\mathbb{I}_{v}(\epsilon)$,
$\mathbb{I}_{v}(\bm{M})$, $\mathbb{I}_{v}(\bm{Y})$, and
$\mathbb{I}_{v}(\bm{\gamma})$ as the set of indices for $\bm{v}$—given in
(15d)—that correspond to $\bm{P}$, $\epsilon$, $\bm{M}$, $\bm{Y}$, and
$\bm{\gamma}$, which can be combined into index set $\mathbb{I}(\bm{v})$.
Next, define

$\mathcal{E}_{P}:=\\{\bm{E}_{P,1},\bm{E}_{P,2},\ldots,\bm{E}_{P,\bar{n}_{x}}\\}$

as the set of standard basis in $\mathbb{S}^{n_{x}}$. Likewise,

$\mathcal{E}_{M}:=\\{\bm{E}_{M,1},\bm{E}_{M,2},\ldots,\bm{E}_{M,n_{x}n_{y}}\\}$

is defined as the set of standard basis in $\mathbb{R}^{n_{x}\times n_{y}}$.
As such, the objective function (31a) can be written as
$\bm{c}^{\top}_{v}\bm{v}$ where vector $\bm{c}_{v}$ is constructed as
$\begin{bmatrix}\bm{0}^{1\times(\bar{n}_{x}+1+2n_{x}n_{y})}\;\;\bm{c}^{\top}\end{bmatrix}^{\top}$
while the linear terms in problem (31) can be written as

$\displaystyle\mathcal{L}_{v}(\bm{v}):=\sum_{i\in\mathbb{I}(\bm{v})}v_{i}\,\mathrm{Blkdiag}\left(\bm{A}^{(1)}_{i},\bm{A}^{(2)}_{i},\bm{A}^{(3)}_{i},\bm{A}^{(4)}_{i}\right),$

where $\bm{A}^{(1)}_{i}$, $\bm{A}^{(2)}_{i}$, $\bm{A}^{(3)}_{i}$,
$\bm{A}^{(4)}_{i}$ are detailed as follows

$\displaystyle\bm{A}^{(1)}_{i}$
$\displaystyle=\begin{cases}\mathrm{Blkdiag}\left(\begin{bmatrix}-\bm{A}^{\top}\bm{E}_{P,i}-\bm{E}_{P,i}\bm{A}&*\\\
-\bm{G}^{\top}\bm{E}_{P,i}&\bm{O}\end{bmatrix},\bm{E}_{P,i},\bm{O}\right),\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\;\;\,\;\text{if}\;i\in\mathbb{I}_{v}(\bm{P})&\\\
\mathrm{Blkdiag}\left(-\gamma_{l}^{2}\bm{I},\bm{I},\bm{O},\bm{I}\right),\qquad\qquad\,\text{if}\;i\in\mathbb{I}_{v}(\epsilon)&\\\
\mathrm{Blkdiag}\left(\bm{E}_{M,i}\bm{C}+\bm{C}^{\top}\bm{E}_{M,i}^{\top},\bm{O}\right),\;\;\text{if}\;i\in\mathbb{I}_{v}(\bm{M})&\\\
\bm{O},\qquad\qquad\qquad\qquad\qquad\qquad\quad\;\;\text{otherwise}\end{cases},$
$\displaystyle\bm{A}^{(2)}_{i}$
$\displaystyle=\begin{cases}\mathrm{Diag}\left(-H_{1i},-H_{2i},\ldots,-H_{n_{h}i}\right),&\,\text{if}\;i\in\mathbb{I}_{v}(\bm{Y})\\\
\bm{O},&\,\text{otherwise}\end{cases},$ $\displaystyle\bm{A}^{(3)}_{i}$
$\displaystyle=\begin{cases}\mathrm{Diag}\left(-\Phi_{1i},-\Phi_{2i},\ldots,-\Phi_{2n_{x}n_{y}i}\right),&\text{if}\;i\in\mathbb{I}_{v}(\bm{Y})\\\
\bm{O},&\text{otherwise}\end{cases},$ $\displaystyle\bm{A}^{(4)}_{i}$
$\displaystyle=\begin{cases}\mathrm{Diag}\left(-\Xi_{1i},-\Xi_{2i},\ldots,-\Xi_{4n_{x}n_{y}i}\right),&\text{if}\;i\in\mathbb{I}_{v}(\bm{z})\\\
\bm{O},&\text{otherwise}\end{cases},$

where
$\mathbb{I}_{v}(\bm{z}):=\left(\mathbb{I}_{v}(\bm{M}),\mathbb{I}_{v}(\bm{Y}),\mathbb{I}_{v}(\bm{\gamma})\right)$.
Accordingly, the constant terms in problem (31) can be written as

$\mathcal{C}_{v}:=\mathrm{Blkdiag}\left(\bm{C}^{(1)},\bm{C}^{(2)},\bm{C}^{(3)},\bm{C}^{(4)}\right)$

where $\bm{C}^{(1)}$, $\bm{C}^{(2)}$ $\bm{C}^{(3)}$, $\bm{C}^{(4)}$ are
constructed as

$\displaystyle\bm{C}^{(1)}$
$\displaystyle=-\kappa\bm{I},\quad\bm{C}^{(2)}=\mathrm{Diag}\left(h_{1},h_{2},\ldots,h_{n_{h}}\right),$
$\displaystyle\bm{C}^{(3)}$
$\displaystyle=\mathrm{Diag}\left(\psi_{1},\psi_{2},\ldots,\psi_{2n_{x}n_{y}}\right),\;\;$
$\displaystyle\bm{C}^{(4)}$
$\displaystyle=\mathrm{Diag}\left(\varphi_{1},\varphi_{2},\ldots,\varphi_{4n_{x}n_{y}}\right),$

in which $\kappa>0$ is a relatively small constant to ensure positive
definiteness. The resulting problem can now be presented as

$\displaystyle\operatorname*{minimize}_{\bm{v}}\;\;$
$\displaystyle\bm{c}^{\top}_{v}\bm{v}$ $\displaystyle\operatorname*{subject\
to}\;\;$ $\displaystyle\mathcal{L}_{v}(\bm{v})+\mathcal{C}_{v}\succeq 0,$
$\displaystyle v_{i}\in\\{0,1\\},\;\;\forall i\in\mathbb{I}_{v}(\bm{\gamma}),$

which is P4. This ends the proof.  $\blacksquare$


# Exploring Multiple Neighborhood Neural Cellular Automata (MNNCA) for
Enhanced Texture Learning

Magnus Petersen
FIAS
Goethe University Frankfurt
Ruth-Moufang-Straße 1, 60438 Frankfurt am Main
<EMAIL_ADDRESS>

###### Abstract

Cellular Automata (CA) have long been foundational in simulating dynamical
systems computationally. With recent innovations, this model class has been
brought into the realm of deep learning by parameterizing the CA’s update rule
using an artificial neural network, termed Neural Cellular Automata (NCA).
This allows NCAs to be trained via gradient descent, enabling them to evolve
into specific shapes, generate textures, and mimic behaviors such as swarming.
However, a limitation of traditional NCAs is their inability to exhibit
sufficiently complex behaviors, restricting their potential in creative and
modeling tasks. Our research explores enhancing the NCA framework by
incorporating multiple neighborhoods and introducing structured noise for seed
states. This approach is inspired by techniques that have historically
amplified the expressiveness of classical continuous CA. All code and example
videos are publicly available on Github.

##  1 Introduction

Cellular Automata (CA) are mathematical models composed of a grid of cells of
either continuous or discrete states. The state of each cell evolves over time
based on a set of rules that consider the states of neighboring cells.
Historically, John von Neumann first conceptualized CA in the 1950s with the
aim of creating a self-replicating system von Neumann [1966]. Later, in the
1970s, John Conway popularized CA with his "Game of Life,"Conway [1970] a
simple yet powerful demonstration of how basic rules can lead to complex and
unpredictable patterns.

These models have since been recognized as a potent tool in the computational
modeling of dynamical systems. Originating from simple rules and interactions,
CA can produce intricate patterns and behaviors, making them an attractive
model for a myriad of applications ranging from physics to biology. With the
advent of deep learning, there has been a growing interest in merging
traditional CA principles with neural networks, leading to the development of
Neural Cellular Automata (NCA) Gilpin [2019]. By parameterizing the update
rules of CA with artificial neural networks, NCAs can be trained, adapted, and
refined using gradient descent, opening up new avenues for modeling and
simulation. This approach has been used to train the NCAs to grow into shapes
Randazzo et al. [2023] form textures Mordvintsev et al. [2021] and has even
been extended to graphs Grattarola et al. [2021] to allow the modeling of an
even wider range of dynamics.

However, while NCAs offer promising capabilities, they are not without
limitations. Their potential in emulating complex behaviors is often limited.
Recognizing this gap, our study seeks to augment the expressiveness of NCAs.
Drawing inspiration from techniques developed by Slackermanz for continuous CA
SLACKERMANZ [2021] that have enhanced CAs capabilities, we delve into the
integration of multiple neighborhoods and the introduction of structured noise
as NCA seed state leading to more complex and creatively interesting dynamics.

##  2 Method

Our approach to enhancing Neural Cellular Automata (NCA) involves two primary
modifications. First, we integrate multiple neighborhoods into the NCA
framework, allowing for richer interactions and dependencies between cells and
multiscale dynamics if different neighborhood sizes are used. Second, we
introduce structured noise, like Perlin noise Perlin [2002], as the initial
condition making the generation of structured textures easier for the
automata. To test the improvement in complex behavior we trained the model on
a texture loss Gatys et al. [2016]Heitz et al. [2021].

Figure 1: Illustration of the MNNCA algorithm.

##  3 Results

(a) Target Texture hip [2023]

(b) NCA after 50 time steps

(c) NCA after 200 time steps

(d) MNNCA after 50 time steps

(e) MNNCA after 200 time steps

(f) MNNCA after 600 time steps

Figure 2: Comparison of textures generated by NCA and MNNCA at different time
steps. Both models were trained with a 256x256 resolution and sampled at
512x512.

##  4 Ethical Implications

Given the nature of our study, which primarily focuses on computational
modeling and its technical advancements for topics like texture generation, we
have not identified any direct ethical issues or implications. However, we
acknowledge the importance of ongoing ethical reflection as this technology
evolves and is applied in broader contexts.

## References

  * dre [2023] dreamstime, 2023.  URL thumbs.dreamstime.com/b/ trama-verniciata-con-arcobaleno- senza-saldatura-neretto-psichedelico-neon- artistico-scarabocchio-di-matite- lavato-imperfetto-190793728.jpg.
  * hip [2023] Digital multicolor curves image, 2023.  URL i0.hippopx.com/photos/521/378/436/digital-multicolor-colorful-curves-preview.jpg.
  * Conway [1970] John Conway.  The game of life.  _Scientific American_ , October 1970.
  * FractalCaleidoscope [2023] FractalCaleidoscope.  Stock psychedelic texture, 2023.  URL images-wixmp-ed30a86b8c4ca887773594c2.wixmp.com/f/056fb59c-3e38- 4de4-813f-a3cc869f33e7/dees4mu-1b402149-a698-4cd3-b384-483ba01ebd9e.jpg/v1/fill/w_894,h_894,q_70,strp/stock__psychedelic_texture_by_fractalcaleidoscope_dees4mu-pre.jpg.
  * Gatys et al. [2016] Leon A Gatys, Alexander S Ecker, and Matthias Bethge.  Image style transfer using convolutional neural networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2414–2423, 2016.
  * Gilpin [2019] William Gilpin.  Cellular automata as convolutional neural networks.  _Physical Review E_ , 100(3):032402, 2019.
  * Grattarola et al. [2021] Daniele Grattarola, Lorenzo Livi, and Cesare Alippi.  Learning graph cellular automata.  _Advances in Neural Information Processing Systems_ , 34:20983–20994, 2021.
  * Heitz et al. [2021] Eric Heitz, Kenneth Vanhoey, Thomas Chambon, and Laurent Belcour.  A sliced wasserstein loss for neural texture synthesis.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 9412–9420, 2021.
  * Mordvintsev et al. [2021] Alexander Mordvintsev, Eyvind Niklasson, and Ettore Randazzo.  Texture generation with neural cellular automata.  _arXiv preprint arXiv:2105.07299_ , 2021.
  * Perlin [2002] Ken Perlin.  Improving noise.  In _Proceedings of the 29th annual conference on Computer graphics and interactive techniques_ , pages 681–682, 2002.
  * Randazzo et al. [2023] Ettore Randazzo, Alexander Mordvintsev, and Craig Fouts.  Growing steerable neural cellular automata.  _arXiv preprint arXiv:2302.10197_ , 2023.
  * SLACKERMANZ [2021] SLACKERMANZ.  Understanding multiple neighborhood cellular automata, 2021.  https://slackermanz.com/understanding-multiple-neighborhood-cellular-automata/.
  * von Neumann [1966] John von Neumann.  _Theory of Self-Reproducing Automata_.  University of Illinois Press, Urbana, IL, 1966.

##  5 Supplementary

##  6 Update Rules Explanation

Our method can be used in different ways:

  1. 1.

Sum Update: This rule directly aggregates all MNNCA output channel values
without any weighting or selection. Mathematically, it is represented as:

$y=\sum_{i}y_{\text{MNNCA output channel}[i]}$

  2. 2.

Random Update: This rule employs a stochastic approach by generating a random
mask to weight and mix the MNNCA output channel values. The mask is normalized
using a softmax function to ensure the sum of weights equals one. The output
is given by:

$y=\sum_{i}y_{\text{MNNCA output channel}[i]}\times\text{choice\\_mask}[i]$

  3. 3.

Based on Environment Channel Value: This rule derives weights from one of the
channels of the environment to mix the MNNCA output channel values. For two
outputs, the rule can be represented as:

$y=y_{\text{MNNCA output channel}[0]}\times\text{choice\\_mask}+y_{\text{MNNCA
output channel}[1]}\times(1-\text{choice\\_mask})$

For more than two outputs, the weights are derived from multiple channels of
the environment channel value and normalized using a softmax function.

  4. 4.

Based on MNNCA Output Channel Value: This rule uses the MNNCA output channel
values themselves to determine the weights for mixing. For two outputs, the
rule is:

$y=y_{\text{MNNCA output channel}[0]}\times\text{choice\\_mask}+y_{\text{MNNCA
output channel}[1]}\times(1-\text{choice\\_mask})$

For multiple outputs, specific channels from all MNNCA output channel values
are aggregated to determine the weights, which are then normalized using a
softmax function.

The demonstrations in the paper are based on 3. However, the code offers the
ability to experiment with all four.

###  6.0.1 Experimental Details

Table 1: Experimental Settings Model | Parameter Count | Rule Count | Batch Count
---|---|---|---
NCA | 2806 | 1 | 3000
NCA XL | 37576 | 1 | 3000
MNNCA | 37536 | 2 | 3000

###  6.0.2 Ablation Experiments

NCA & Uniform Noise

NCA & Perlin Noise

MNNCA & Uniform Noise

MNNCA & Perlin Noise

Figure 3: Illustrations of various timesteps, 50, 200, and 600, of classical
NCA and MNNCA, with uniform and Perlin noise.Trained on texture hip [2023].

NCA & Uniform Noise

NCA & Perlin Noise

MNNCA & Uniform Noise

MNNCA & Perlin Noise

Figure 4: Illustrations of various timesteps, 50, 200, and 600, of classical
NCA and MNNCA, with uniform and Perlin noise. Trained on texture dre [2023].

NCA & Uniform Noise

NCA & Perlin Noise

MNNCA & Uniform Noise

MNNCA & Perlin Noise

Figure 5: Illustrations of various timesteps, 50, 200, and 600, of classical
NCA and MNNCA, with uniform and Perlin noise. Trained on texture
FractalCaleidoscope [2023].

###  6.0.3 Comparison with Fixed Parameter Count

To ensure a fair comparison, given that MNNCAs inherently possess more
parameters, we trained an NCA with a parameter size comparable to that of
MNNCA. Details of this can be found in Section 6.0.1. The observed results
suggest that the enhanced capability of MNNCAs to generate complex textures is
not solely attributed to their increased parameter count as the NCA still
generates less structured textures.

NCA

NCA XL

Figure 6: Illustrations of various timesteps, 50, 200, and 600, of classical
NCA with uniform noise and different sizes.Trained on texture hip [2023].