Particle filtering of dynamical networks: Highlighting observability issues

TL;DR

This paper introduces a particle filtering framework to evaluate observability in high-dimensional dynamical networks, revealing how dynamics and topology influence sensor placement and network state reconstruction.

Contribution

It proposes a novel particle filtering approach to assess network observability, incorporating the effects of dynamics, topology, and sensor variable choices.

Findings

Heterogeneous nodal dynamics improve observability with dynamical affinity.

Sensor variable selection impacts particle filtering performance.

PF framework aligns with established observability metrics in benchmarks.

Abstract

In a network of high-dimensionality, it is not feasible to measure every single node. Thus, an important goal in the literature is to define the optimal choice of sensor nodes that provides a reliable state reconstruction of the network system state-space. This is an observability problem. In this paper, we propose a particle filtering (PF) framework as a way to assess observability properties of a dynamical network, where each node is composed by an individual dynamical system. The PF framework is applied on two benchmarks, networks of Kuramoto and R\"ossler oscillators, to investigate how the interplay between dynamics and topology impacts the network observability. Based on the numerical results, we conjecture that, when the network nodal dynamics are heterogeneous, better observability is conveyed for sets of sensor nodes that share some dynamical affinity to its neighbourhood. Moreover, we also investigate how the choice of an internal measured variable of a multidimensional sensor node affects the PF performance. The PF framework effectiveness as an observability measure is compared to a well-consolidated nonlinear observability metric for a small network case and some chaotic systems benchmarks.

Figures

Tables

Table 1. Table 1: Comparison of the observability degree and the PF estimation error of chaotic systems benchmarks.

Chaotic system		${\delta }_o$11footnotemark: 1	$\eta$ [Q1 Q2 Q3]22footnotemark: 2	${\Sigma }_w$	${\Sigma }_v$
	x	0.022	[0.68 1.02 1.45]	0.01	0.5
Rössler system	y	0.133	[0.62 0.96 1.39]	0.01	0.5
	z	0.0063	[1.29 1.91 2.92]	0.01	0.5
		$y⊳x⊳z$	$y⊳x⊳z$
	x	0.0104	[0.12 0.20 0.40]	0.001	0.05
Cord system	y	0.0005	[0.27 0.39 0.59]	0.01	0.5
	z	0.0005	[0.28 0.42 0.62]	0.01	0.5
		$x⊳y\approx z$	$x⊳y\approx z$
	x	–	[0.24 0.37 0.54]	0.05	0.6
Double-scroll33footnotemark: 3	y	–	[0.43 0.64 0.99]	0.01	0.2
	z	–	[0.22 0.35 0.56]	0.05	0.6
		$x⊳z⊳y$	$z\approx x⊳y$
	x	$6.5\times {10}^{-6}$	[2.56 3.64 5.14]	1	2
Lorenz system	y	$8.2\times {10}^{-6}$	[2.13 2.85 3.88]	1	2
	z	$2.2\times {10}^{-5}$	[8.73 17.0 33.6]	1	2
		$z⊳x⊳y$	$y⊳x⊳z$
	x	$0.0142$	[4.07 6.29 8.25]	0.005	0.25
Lorenz’84	y	$0.0009$	[1.99 2.68 3.83]	0.01	0.5
	z	$0.0014$	[2.05 2.80 4.07]	0.01	0.5
		$x⊳z\approx y$	$y\approx z⊳x$

Equations

$${\mathbf{x}}_k$$

$${\mathbf{y}}_k$$

$${\mathbf{x}}_k^{(p)}\sim q({\mathbf{x}}_k^{(p)}|{\mathbf{x}}_{k-1}^{(p)},{\mathbf{y}}_k),\forall p=1,2,\dots ,N_p,$$

$$W_k^{(p)}=W_{k-1}^{(p)}\cdot \frac{p({\mathbf{y}}_k|{\mathbf{x}}_k^{(p)})\cdot p({\mathbf{x}}_k^{(p)}|{\mathbf{x}}_{k-1}^{(p)})}{q({\mathbf{x}}_k^{(p)}|{\mathbf{x}}_{k-1}^{(p)},{\mathbf{y}}_k)},\forall p.$$

$${\overline{W}}_k^{(p)}=\frac{W_k^{(p)}}{{\sum }_{p=1}^{N_p}{W}_k^{(p)}},\forall p.$$

$${\hat{N}}_{\mathrm{eff}}=\frac{1}{{\sum }_{i=1}^{N_p}({\overline{W}}_k^{(p)}{)}^2}.$$

$${\hat{\mathbf{x}}}_k=\frac{1}{N_p}\sum _{p=1}^{N_p}{\mathbf{x}}_k^{(p)}{\overline{W}}_k^{(p)},$$

$$\dot{\mathbf{x}}$$

$$\mathbf{y}$$

$${\mathcal{O}}_h(\mathbf{x})=\left[\begin{array}{c}\text{partialderivative}{\mathcal{L}}_f^{0}h(\mathbf{x})\mathbf{x}\\ \text{partialderivative}{\mathcal{L}}_f^{1}h(\mathbf{x})\mathbf{x}\\ ⋮\\ \text{partialderivative}{\mathcal{L}}_f^{n-1}h(\mathbf{x})\mathbf{x}\end{array}\right],$$

$${\delta }_o(\mathbf{x})=\left|\frac{{\lambda }_{\min }\left({\mathcal{O}}_h(\mathbf{x}{)}^{\mathsf{T}}{\mathcal{O}}_h(\mathbf{x})\right)}{{\lambda }_{\max }\left({\mathcal{O}}_h(\mathbf{x}{)}^{\mathsf{T}}{\mathcal{O}}_h(\mathbf{x})\right)}\right|,$$

$${\delta }_o=\frac{1}{t_1-t_0}{\int }_{t_0}^{t_1}{\delta }_o(\mathbf{x}(\tau ))\mathrm{d}\tau .$$

$${\dot{x}}_i={\omega }_i+\rho \sum _{j=1}^N{a}_{ij}\sin (x_j-x_i),i=1,\dots ,N,$$

$${\mathbf{y}}_k=C{\mathbf{x}}_k+{\mathbf{v}}_k,$$

$${\mathbf{x}}_k^{(p)}\sim \mathcal{N}\left(f({\mathbf{x}}_{k-1}^{(p)}),{\Sigma }_w\right),p=1,\dots ,N_p.$$

$$W_k^{(p)}$$

$$=W_{k-1}^{(p)}\cdot \frac{1}{(2\pi {)}^{q/2}{\sigma }_v^q}\cdot \exp -\frac{1}{2}\text{norm}{\frac{C{\mathbf{x}}_k^{(p)}-{\mathbf{y}}_k}{{\sigma }_v}}^2,\forall p.$$

$$\{\begin{array}{ll}{\dot{x}}_i& =-y_i-z_i\\ {\dot{y}}_i& =x_i+a_i{y}_i+\rho {\sum }_{j=1}^N{a}_{ij}(y_j-y_i)\\ {\dot{z}}_i& =b_i+z_i(x_i-c),\end{array},i=1,\dots ,N$$

$${\mathbf{y}}_k=C\cdot {\left[\begin{array}{cccc}{\mathbf{x}}_{1,k}^{\mathsf{T}}& {\mathbf{x}}_{2,k}^{\mathsf{T}}& \dots & {\mathbf{x}}_{N,k}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}+{\mathbf{v}}_k$$

$${\eta }_i=\text{median}\left\{\frac{1}{\Delta {\mathbf{x}}_i}\sqrt{\frac{1}{N}\sum _k{\left({\mathbf{x}}_{i,k}-{\hat{\mathbf{x}}}_{i,k|q}\right)}^2}\right\},$$

$$C=\left[\begin{array}{ccccc}1& 0& \dots & 0& 0\\ 0& 0& \dots & 0& 1\end{array}\right].$$

Full text

**Particle filtering of dynamical networks: Highlighting observability issues **

Arthur N. Montanari,

Graduate Program in Electrical Engineering of the Universidade Federal de Minas Gerais (UFMG), Avenida Antônio Carlos 6627, 31270-901, Belo Horizonte, Minas Gerais, Brazil.

Luis A. Aguirre

Departamento de Engenharia Eletrônica, UFMG.

Abstract

In a network of high-dimensionality, it is not feasible to measure every single node. Thus, an important goal in the literature is to define the optimal choice of sensor nodes that provides a reliable state reconstruction of the network system state-space. This is an observability problem. In this paper, we propose a particle filtering (PF) framework as a way to assess observability properties of a dynamical network, where each node is composed by an individual dynamical system. The PF framework is applied on two benchmarks, networks of Kuramoto and Rössler oscillators, to investigate how the interplay between dynamics and topology impacts the network observability. Based on the numerical results, we conjecture that, when the network nodal dynamics are heterogeneous, better observability is conveyed for sets of sensor nodes that share some dynamical affinity to its neighbourhood. Moreover, we also investigate how the choice of an internal measured variable of a multidimensional sensor node affects the PF performance. The PF framework effectiveness as an observability measure is compared to a well-consolidated nonlinear observability metric for a small network case and some chaotic systems benchmarks.

The high-dimensionality of network systems makes it unfeasible to measure every single node. Physically, it is not (yet) possible to measure every single neuron of the one hundred billion neurons present in a brain network. Likewise, it is not of economic interest to place a phasor measurement unit in every single electrical substation of a power system. Indeed, the problem of observability is recurrent in many network systems. Thus, it is only natural that in the past ten years, the literature has been flooded by innumerable methods that provide a minimum selection of sensor nodes that claim to be sufficient and/or necessary to render a network observable. However, a framework that allows one to compare the effectiveness and practicability of these methods is still missing. Moreover, the problem becomes intrinsically more complex when dealing with dynamical networks—that is, networks whose nodes are individual dynamical systems. In this paper, we develop a framework based on Bayesian filtering as a way to assess the degree of observability conveyed by a given set of sensor nodes. This allows a further investigation of the interplay between observability, nodal dynamics and topology on a network system, as well as a means of comparison of different methods provided in the literature.

1 Introduction

Introduced by Rudolf Kalman [26], observability is a property that determines if the initial state of a dynamical system can be reconstructed solely based on input and output signals [10]. The classic concept of observability, and—by duality—controllability, is based on a crisp definition, i.e. the system is either observable (controllable) or not.

Perhaps, a more practical question might be whether an observable system is almost unobservable. To that end, this crisp concept was extended by metrics that quantify observability in a continuous manner, measuring how well the system trajectory can be reconstructed from its observations [18, 1, 45]. This approach is useful to investigate which variables are more relevant to render a system observable for practical purposes, e.g. state reconstruction or system embedding[28, 30].

Since, in network systems, it is practically and numerically unfeasible to perform measurements on every single node, a selection of sensor nodes, in which both quantity and positioning must be optimal, is required. This is an observability problem. However, due to high dimensionality, classical observability metrics face serious numerical and scalability issues. To circumvent these challenges, Liu et al. pioneering works in controllability [35] and observability [36] took a graph-inspired approach, grounded on Lin’s controllability (observability) definition [33, 9], to develop an algorithm that select the “minimum” set of driver (sensor) nodes under which a complex network is controllable (observable). Because of the intuitive representation of network systems by graph models and the high scalability of graph metrics, many following works [27, 42, 25, 37, 19, 51] embraced this graph-theoretical approach to study the controllability and observability in networks.

Faced with several different proposals of observability (controllability) methods, a recurrent goal in the literature revolves around the question whether a novel method really estimates a “better” set of sensor (driver) nodes than other methods, or if the estimated set at least renders a system observable (controllable) for practical purposes. However, a benchmark framework to compare and validate the efficacy of novel metrics in the literature is still absent, as pointed out in Ref. [8].

For instance, it is arguable that works grounded on crisp definitions of observability (controllability), such as Lin’s structural definition[33, 9] and following works [35, 36], are not really apt to classify if a certain set of sensor nodes is really the best one or even if its feasible. Indeed, the pioneering methods of Liu et al.[35, 36] have been shown to underestimate the required set of sensor[23, 32] and driver[20, 27, 47] nodes. A natural approach to circumvent this problem is to extend metrics that gradually quantify system observability (controllability), under a specific set of sensor (driver) nodes, to a network context [39, 21, 45].

On the other hand, most network systems are only investigated exclusively from its network topology, i.e. the corresponding adjacency matrix [27], which is insufficient from an observability point of view [5]. However, most networks used as benchmarks in [35, 52], including food webs, regulatory networks, power grids, electronic circuits and neuronal networks, are in fact dynamical networks (see Table 1 of Ref. [27, Supplementary information]), i.e. a connected graph in which each node has an individual (multidimensional) dynamical system. Cowan et al. [14] show that significantly different results could be derived by considering the presence of self-edges that dictate the independent dynamic behaviour of any given node. This includes a new level of depth to the observability problem of network systems, specially when nonlinear dynamics are considered[14, 39, 21, 45, 44, 27, 23, 5, 32].

As pointed out in Ref. [5], a dynamical network can be decomposed in three levels: i) the nodal dynamics, described by an individual dynamical system in each node; ii) the network topology, described by a graph; and iii) the full network, which is a combination of the aforementioned levels. Figure 1 illustrates the corresponding graph of a given network topology when nodal dynamics are considered. Despite recent investigations using small motifs [38, 50, 44], the interplay between observability, nodal dynamics, especially with nonlinearities, and network topology is still not clear [23, 5].

1.1 Paper contributions

In this paper, we present a Bayesian filter formulation, based on the particle filter111Also known as sequential importance resampling (SIR) filter. (PF), for dynamical network applications as a way to assess its observability properties. This is motivated by the fact that if the measured signals do not provide relevant information to the filter, i.e. they convey poor observability, then the update stage of the filter is impaired and, consequently, the estimates show poor performance. Thus, we argue that the provided PF formulation can be used as benchmark for observability studies, that is, as a means to validate and compare, on a practical level, different observability metrics and sensor selection methods proposed in the literature. Likewise, a moving horizon estimation technique [23] as well as the fitting error from training a Gauss-Newton algorithm [22] or a reservoir computer [8] have also been proposed in the literature as alternatives to assess the “reconstruction quality” of a dynamical system trajectory.

Having established a numerical framework based on the PF, this paper investigates a number of issues related to the interplay between the network topology and nodal dynamics from an observability point of view. To this end two benchmark systems were used: a network of Kuramoto oscillators and a network of Rössler systems. In the former, every node is composed of a 1-dimensional model, the focus is exclusively on the impact of the choice of sensor nodes. In the latter, each node has a 3-dimensional model and this allows the study of the influence of sensor node selection and the choice of dynamical variables at each sensor node  in terms of observability.

The paper is organized as follows. Section 2 provides a background on observability metrics and PF. Section 3 describes the PF methodology and numerical setup applied in each network benchmark, while Section 4 presents the numerical results. Section 5 concludes the work.

2 Background

2.1 Particle filtering

PF is a suboptimal state estimator based on sequential Monte Carlo methods to solve statistical problems. This method uses an independent number of stochastic samples, so-called particles, extracted directly from the state-space, that are recursively located, propagated and weighted accordingly to Bayes’ theorem. The method describes the states probability density function (PDF) by sampling approximation, not relying on any linearisation method or function approximation. Thus, it is applicable to any nonlinear and non-Gaussian problem, such as systems with several nonlinearities (e.g. chaotic systems) or multi-modal PDFs. The drawback is the computational burden, although the ever-increasing computational power already allows PF to be used in online applications [17], specially multi-target tracking [48].

This section provides background on the PF framework. The reader is referred to Refs. [16, 12, 17, 7] for details. Let us consider the following discrete-time autonomous nonlinear dynamical system:

[TABLE]

where $k=1,2,\dots$ are the time indices, $\bm{x}_{k}\in\mathbb{R}^{n}$ is the $n$-dimensional state at time instant $k$, $\bm{y}_{k}\in\mathbb{R}^{q}$ is the $q$-dimensional observation at time $k$, $f:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ and $h:\mathbb{R}^{n}\mapsto\mathbb{R}^{q}$ are nonlinear functions, and $\bm{w}_{k}\in\mathbb{R}^{n}$ and $\bm{v}_{k}\in\mathbb{R}^{q}$ are independent additive process and observation noises, respectively with known distributions.

Given the system model (1) and a sequence of observations $\bm{y}_{1:k}=\left\{\bm{y}_{1},\bm{y}_{2},\ldots,\bm{y}_{k}\right\}$, the goal is to estimate the posterior distribution $p(\bm{x}_{k}|\bm{x}_{1:k-1},\bm{y}_{1:k})$ sequentially, that is, from a prior distribution $p(\bm{x}_{k-1}|\bm{x}_{1:k-2},\bm{y}_{1:k-1})$. It is assumed that model (1) is a Markovian process, i.e. future states depend only upon the present state, yielding $p(\cdot|\bm{x}_{1:k-1})=p(\cdot|\bm{x}_{k-1})$ and $p(\cdot|\bm{y}_{1:k})=p(\cdot|\bm{y}_{k})$. Recall that the PF approximates a posterior distribution $p(\bm{x}_{k}|\bm{x}_{k-1},\bm{y}_{k})$ with random samples, which are propagated from time instant $k-1$ to $k$ by a mathematical model. The particles are assigned to weights based on the likelihood distribution $p(\bm{y}_{k}|\bm{x}_{k}^{(p)})$ of the observations at time instant $k$, where $\bm{x}_{k}^{(p)}$ is a realization of the $p$-th particle. Finally, the filter output is the weighted average of all particle realizations. Algorithm 1 summarizes these proceedings.

Some remarks are in order: (i) for resampling methods details, e.g. multinomial, residual and systematic, see Ref. [7]; (ii) a common threshold[12] is $N_{\rm t}=N_{p}/2$ or $N_{p}/3$; and (iii) weights are usually calculated in logarithmic scale to improve numerical precision. Noise process must be carefully designed: a smaller noise process causes sample impoverishment in very few iterations, while higher noise process can provide estimates with larger deviations [12].

2.2 Observability

This section provides background on observability metrics for nonlinear systems. The reader is referred to Refs. [24, 30, 3] for details. Let us consider the continuous-time autonomous nonlinear dynamical system

[TABLE]

where $\bm{x}\in\mathbb{R}^{n}$, $\bm{y}\in\mathbb{R}^{q}$, $f:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ and $h:\mathbb{R}^{n}\mapsto\mathbb{R}^{q}$.

The observability matrix of (7) can be formulated as

[TABLE]

where $\mathcal{L}^{j}_{f}h(\bm{x})$ is the $j$th-order Lie derivative of $h$ along the vector field $f$. System (7) is said to be observable if (8) is full rank $\forall\bm{x}$, i.e. ${\rm rank}(\mathcal{O}_{h}(\bm{x}))=n,~{}\forall\bm{x}$. Otherwise, it is unobservable.

In order to quantify the system observability in a continuous manner, it is helpful to assess how far $\mathcal{O}_{h}(\bm{x})$ is from being rank-deficient. Following the ideas in Refs. [18, 1], adapted to the nonlinear case[31, 28], this can be achieved by computing a coefficient $\delta_{o}(\bm{x})$ that quantifies the numerical ill-conditioning of $\mathcal{O}_{h}(\bm{x})$ at a given state $\bm{x}$:

[TABLE]

where $0\leq\delta_{o}(\bm{x})\leq 1$ and $\lambda_{\rm max}(\cdot)$ indicates the maximum eigenvalue (likewise for $\lambda_{\rm min}$). In principle, the higher $\delta_{o}(\bm{x})$, the more observable (7) is. If $\delta_{o}(\bm{x})=0$, the system is unobservable. Averaging $\delta_{o}$ along a trajectory over interval $t\in\left[t_{0},t_{1}\right]$ yields a “global” observability coefficient

[TABLE]

We note that the global coefficient $\delta_{o}$ is not absolute, but relative to the dynamical system. In other words, for a given system, two or more measuring functions can be investigated on the basis of their respective coefficient $\delta_{o}$. However, such observability coefficients of different systems are not necessarily comparable. Check Ref. [2] for details.

3 Particle Filtering of Dynamical Networks

A graph is defined as $\mathcal{G}=\{\mathcal{V},\mathcal{E},A\}$, where $\mathcal{V}=\{v_{1},v_{2},\dots,v_{N}\}$ and $\mathcal{E}\subseteq\mathcal{V}=\mathcal{V}\times\mathcal{V}=\{e_{1},e_{2},\dots,e_{m}\}$ are finite sets of $N$ nodes and $m$ edges, respectively. The adjacency matrix $A=[a_{ij}]$ is a mapping which associates elements (edges) of $\mathcal{E}$ to a pair of elements (nodes) of $\mathcal{V}$. If $(v_{i},v_{j})$ are unidirectionally linked, then $a_{ij}=a_{ji}=1$ (unweighted). Otherwise, $a_{ij}=0$. The cardinality of a set is referred to as $\left|\mathcal{V}\right|=m$.

In the following experiments, there is a dynamical system at each node. A benchmark network is used to generate the data that will feed a particle filter. A hat, e.g. $\hat{\bm{x}}$, is used to indicate estimated values.

The computational framework is based on the PF for convenience. Despite the fact that PFs suffer with the inconvenient curse of dimensionality [43, 46], we deal only with networks with at most 45 states. Furthermore, as our second benchmark involves chaotic behaviour, the PF is seen to be a better alternative due to several nonlinearities and non-Gaussian uncertainties in the system [13, 34]. We argue that Kalman filter-based estimators are not recommended in such cases since its covariance matrix is not computationally feasible in high-dimensional systems, especially those with high linear dependence among its states, such as network systems.

3.1 Network of Kuramoto oscillators

Consider a network of Kuramoto oscillators where each oscillator $v_{i}\in\mathcal{V}$ is represented by a phase angle $x_{i}\in\mathbb{R}$. The dynamical benchmark network can be represented by the following continuous state-space model [49, 11]:

[TABLE]

where $N$ is the network size, $\rho$ is the coupling strength, and $\omega_{i}>0$ is the natural frequency of the $i$th oscillator. Hence the state vector of the full network is $\bm{x}=\begin{bmatrix}x_{1}&x_{2}&\dots&x_{N}\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{N}$. Note that the coupling among the oscillators is additive, diffusive and proportional to $\rho$.

In the simulations $\rho=0.1$, $\omega_{i}\sim\mathcal{N}\left(1,0.03\right)$ and initial conditions $x_{i}(t_{0})\sim\mathcal{U}(-\pi,+\pi)$, $\forall i$, where $\mathcal{N}\left(\mu,\sigma^{2}\right)$ denotes a Gaussian distribution with average $\mu$ and variance $\sigma^{2}$, and $\mathcal{U}\left(a,b\right)$ denotes a uniform distribution within limits $[a,\,b]$. Numerical integration is performed using a fourth-order Runge Kutta algorithm for simulation time of $500$ s, with time step $0.1$ s.

The output vector is obtained as follows:

[TABLE]

where $C\in\mathbb{R}^{q\times n}$ is the output matrix, $q$ is the number of sensor nodes, $\bm{v}_{k}\sim\mathcal{N}(0,\Sigma_{v})$, $\Sigma_{v}=\sigma_{v}^{2}\cdot I_{q}$ is the covariance matrix of $\bm{v}_{k}$, $\sigma_{v}=0.1$ and $I_{q}$ is an identity matrix of dimension $q$. It is assumed that measures are taken independently. The set of sensor nodes is $\mathcal{S}=\{s_{1},s_{2},\ldots,s_{q}\}\subseteq\mathcal{V}$. A node $v_{j}$ will be listed as a sensor node $s_{i}$ if $c_{ij}=1$, where $c_{ij}$ is an element of $C$. If measures are taken on all nodes, then $\mathcal{S}\equiv\mathcal{V}$.

The goal is to investigate how the selection of sensor nodes affects the PF performance. Other important aspects such as robustness to model parameter variation, noise level and dimensionality have been investigated elsewhere [12, 43]. Thus, in the simulations the true values of $\omega_{i}$ will be used.

The bootstrap PF is implemented as follows. Particles are propagated through a discrete numerical approximation $f(\bm{x}_{k-1}^{(p)},\omega)$ of the continuous model (11), obtained analytically via the backward Euler method with time step $0.1$ s. We introduce an additive process noise $\bm{w}_{k}\sim\mathcal{N}\left(0,\Sigma_{w}\right)$—where $\Sigma_{w}=\sigma_{w}^{2}\cdot I_{n}$ is the covariance matrix of $\bm{w}_{k}$ and $\sigma_{w}=0.1$—so that the filter particles have a better coverage of the state-space, as well as avoiding sample impoverishment in few steps. We set $\hat{\bm{x}}(t_{0})=\bm{x}(t_{0})$ for fast convergence and $N_{p}=500$.

In this paper, we implement exclusively the bootstrap PF variant, in which the proposal distribution is assumed to be equal to the transition distribution, yielding $q(\bm{x}_{k}^{(p)}|\bm{x}_{k-1}^{(p)},\bm{y}_{k})=p(\bm{x}_{k}^{(p)}|\bm{x}_{k-1}^{(p)})$. Based on Algorithm 1, we summarize the implemented PF framework for the Kuramoto network in Algorithm 2.

3.2 Network of Rössler systems

Consider a network of Rössler oscillators, in a chaotic regime, coupled by the $y$ variable [40, 6]. Hence, at each node $v_{i}\in\mathcal{V}$ there is a system with three state variables $\bm{x}_{i}=\begin{bmatrix}x_{i}&y_{i}&z_{i}\end{bmatrix}^{\mathsf{T}}$. The dynamical network can be represented by the following state-space model:

[TABLE]

where $N$ is the number of nodes, $n=3N$ is the number of state variables of the network and $(a_{i},b_{i},c_{i})$ are the parameters of the $i$-th Rössler system. Numerical integration is performed using fourth-order Runge Kutta algorithm for a total simulation time of $500$ s, with time step $0.01$ s. Simulation parameters were set to $\rho=0.1$, initial conditions $\bm{x}_{i}(t_{0})\sim\mathcal{N}(0,I_{3})$, $(b_{i},c_{i})=(2,4)$, and $N$ equidistant values within interval $a_{i}\in[0.388,0.408],\forall i$. This corresponds to spiral dynamics for each node in the case of no coupling.

Measures are taken linearly on the network of oscillators as follows:

[TABLE]

where $C\in\mathbb{R}^{q\times n}$, $\bm{v}_{k}\sim\mathcal{N}(0,\Sigma_{v})$ and $\sigma_{v}=0.1$. Note that the observation matrix $C$ defines the network sensor nodes as well as which variables are measured at each oscillator. Hence, given a sensor node, all variables could be measured $s_{i}=(x_{i},y_{i},z_{i})$, or just one (e.g. $s_{i}=y_{i}$) or a combination of two variables.

Algorithm 2 summarizes the implemented PF for the Rössler network. The true values of $(a_{i},b_{i},c_{i})$ are used. We set $N_{p}=500$, $\bm{w}_{k}\sim\mathcal{N}(0,\Sigma_{w})$ and $\sigma_{w}=0.01$ for the reasons mentioned in Section 3.1. Function $f(\bm{x}_{k-1}^{(p)},a,b,c)$ is a discrete numerical approximation of the continuous model (15) via the backward Euler method, with time step $0.01$ s.

4 Numerical Results

The dynamical networks and PFs detailed in Section 3 are considered in this section. Measures are taken independently at each node and the network graph of interest is a chain graph of $N$ nodes.

4.1 Particle filtering framework as a validation benchmark

The quality of a Bayesian filter estimate depends on the observability conveyed by the measured variables. In other words, the higher the observability rank of a network system observed using a particular set of sensor nodes, the better the PF performance, that is, the smaller is the difference between its estimates and the true value. In this section, we analyze the PF framework effectiveness as an alternative mean to assess observability properties of a network system by comparing its performance with the results yielded by the well-consolidated observability rank in (10). Since $\delta_{o}$ becomes unfeasible for high-order systems, we investigate the PF framework on a Kuramoto dynamical network of $N=4$ nodes.

Moreover, it is important to notice that, unlike the observability metrics presented in Section 2.2, the PF is not a deterministic tool since its results heavily depend on the particle realizations. A statistical analysis is required. Thus, in the following numerical simulations, the applied performance index $\eta_{i}$ is the median of the normalized root-mean-square error222Simulations were carried out for different performance metrics, including the mean absolute error, mean absolute percentage error and mean relative error, yielded similar results. (NRMSE) per node $v_{i}$ over 100 Monte Carlo simulations:

[TABLE]

where $\hat{\bm{x}}_{i,k|q}$ refers to estimate $\hat{\bm{x}}_{i,k}$ of the $q$th Monte Carlo simulation (with different process noise realizations), and $\Delta\bm{x}_{i}=\max{\bm{x}_{i}}-\min{\bm{x}_{i}}$ is the measured data range. Note that median is preferred to the mean since the PDF of $\hat{\bm{x}}_{i,k}$ is non-Gaussian. The PF overall performance for a dynamical network is taken as $\eta=\frac{1}{n}\sum_{i=1}^{n}\eta_{i}$.

Fig. 2 compares the computed observability coefficients $\delta_{o}$ and the PF framework performance index $\eta$ for every possible combination of sensor nodes in this network. Firstly, it is interesting to note that all choices of sensor nodes are observable (i.e. $\delta_{o}>0$ for all sets $\mathcal{S}$). This is in line with Liu et al. structural approach [33, 36] since the whole graph is a root-strongly connected component (SCC). However, from a practical point-of-view, as pointed out by $\delta_{o}$ and $\eta$ metrics, some choices of sensor nodes (specially when $\left|S\right|=1$) are almost unobservable. Thus, although Lin’s structural definition of controllability (observability)[33, 9] holds for almost every choice of the parameter space (since it does not depend on the nonzero specific entries of the adjacency matrix), it might fail—for practical purposes—under specific entries[23, 50].

It is evident that the PF performance $\eta$ is not exactly an equivalent metric to $\delta_{o}$. Specially for $|\mathcal{S}|=1$, where the correlation between $\eta$ and $\delta_{o}$ is the lowest. A possible justification is that, since systems of poor observability compromises the PF update stage efficacy, the PF realizations exhibits large interquartile ranges in its estimates which, in turn, also compromises the statistical relevance of $\eta$ for $|\mathcal{S}|$. In other words, if the error bars are taken into account the perfornace of the PF for all options with $|\mathcal{S}|=1$, it is seen that performance is almost equally poor from a statistical point of view. However, when $|\mathcal{S}|\geq 2$, we note that the standard deviation of the PF estimates $\eta$ reduces considerably and results become statistically relevant. Indeed, for $|\mathcal{S}|\geq 2$, there is a high correlation between $\eta$ and $\delta_{o}$, except in two occasions ($\mathcal{S}_{8}$ and $\mathcal{S}_{9}$ are switched in Figs. 2a,b). We argue that although the PF is not effective to quantify or rank the observability of dynamical systems of very poor observability indexes, it is quite robust when the observability index is a little bit better.

Appendix A further investigates the PF efficacy as an observability quantification in chaotic systems benchmarks explored in Ref. [4].

4.2 Network of Kuramoto oscillators

Consider a Kuramoto dynamical network of $N=15$ nodes where measures are taken independently on the chain extremities ($\mathcal{S}_{1}=\{v_{1},v_{15}\}$, see Fig. 3). Thus:

[TABLE]

where $C\in\mathbb{R}^{2\times 15}$. Following the results in Section 4.1, it is implied that a node (state) with the worst PF estimate (highest $\eta_{i}$) is the least observable node from a particular set of sensor nodes. Hence, the optimal selection of sensor nodes is actually the one that provides minimum overall estimation error $\eta$ when compared to all other options with the same number of sensor nodes.

Guan et al. show, through their proposed metric, that the observability level of a particular node seems to (intuitively) decrease with the increase of the path length between this node and a particular set of sensor nodes [22]. Such analysis is directly related to the network topology and, considering the chain network setup in Fig. 3, one might argue that the expected “least observable node” from $\mathcal{S}=\{v_{1},v_{15}\}$ would be the central node $v_{8}$. However, the numerical setup from which the authors of Ref. [22] draw this particular conclusion assume that all nodal dynamics are equal (i.e. every node is governed by the same ordinary differential equations with the same parametric values).

But if the nodal dynamics were heterogeneous (i.e. a network with Kuramoto oscillators of different natural frequencies $\omega_{i}$ along the network structure), would $v_{8}$ still be the least observable one? It is known that the location of oscillators (with different $\omega_{i}$) in a network affects the synchronization process[6, 15], but how does it affect observability? For instance, Refs. [50, 23] showed through “small motifs” that the optimal $\mathcal{S}$ depends not only on the ordinary differential equations that govern the system behaviour but also on its specific parameters.

It is here argued that the observability that results from the interplay between dynamics and topology can be understood assessing the performance of the PF estimates for different networks that share the same topology, but different dynamics. Figure 4 presents the attained filter estimates $\eta_{i}$ for different allocations of $\omega_{i}$ over the chain network in Fig. 3.

Despite all networks in Figs. 4a,c,d being composed by the same Kuramoto oscillators, we see that their placing order along the chain graph affects the PF performance in different manners. Depending on the oscillators location, the PF worst estimate shifts away from $v_{8}$, presenting different $\eta_{i}$ distributions over the network. For instance, in Figs. 4a,b we see a bimodal distribution of $\eta_{i}$, while in Fig. 4c,d there is a somewhat unimodal distribution. There is an implication between the degree of observability and the graph topology and dynamic interactions that seems to be, likewise in synchronization analysis [6], related to the degree of mismatch between the dynamics of an oscillator and its neighbours. Adopting the convention in Ref. [6], we refer as soft transition if there is a relatively small frequency mismatch in a neighbourhood, or as hard transition otherwise.

In the following analyses, we conjecture that the oscillators nodes chosen as sensor nodes not only have to be topologically central to the network, but also must have a dynamical behaviour that is representative to (similar to) the dynamics of the other oscillators. In a sense, node oscillators chosen as sensor nodes in a dynamical network should be central both topologically and dynamically.

Consider Fig. 4a. Note that, although the central node $v_{8}$ is the farthest from a topological sense, it does not convey the worst observability, i.e. the higher estimate error $\eta_{i}$. Following the conjecture, despite its topological distance, $\omega_{8}$ is very close to the sensor nodes natural frequencies $(\omega_{1},\omega_{15})$. On the other hand, estimates $(\eta_{4},\eta_{10})$ and their respective neighbourhood show worse performance than $\eta_{8}$, despite being relatively closer than $v_{8}$ to the sensor nodes. The dynamical behaviour, however, is not quite central, with a high frequency mismatch between the oscillators $(v_{4},v_{10})$ and their respective nearest sensor node $(v_{1},v_{15})$. Interestingly, $\eta_{4}<\eta_{10}$ can be justified by the fact that the degree of frequency mismatch between the neighbourhood of $v_{4}$ and its closer sensor node $s_{1}$ is smaller compared to $v_{10}$ and $s_{2}$. From the point-of-view of observability, in some way, the presence of dynamical centrality compensates for the lack of topological centrality, and vice-versa. Nevertheless, $v_{2}$ still provides good observability (low $\eta_{2}$) despite having a high frequency mismatch between $(v_{1},v_{2})$. This might be attributed to the close topological proximity of these nodes on the network, relieving the dynamic difference effects. In Fig. 4b, by only adjusting the oscillators frequencies $\omega_{1}$ and $\omega_{15}$ at the chain ends to respectively closer frequencies of $\omega_{4}$ and $\omega_{10}$ (and neighbourhood), the estimate error $\eta_{i}$ at this sections greatly reduced, with $\eta$ being approximately $30\%$ smaller than in Fig. 4a.

Albeit the estimation error $\eta$ is certainly related to the degree of observability conveyed by the respective choice of $\mathcal{S}$, it is important to note that results are also impacted by the intrinsic dynamical relations. For instance, certain placing orders of oscillators along the network chain can give rise to formation of synchronous cluster manifolds, which effectively reduce the system dimension when complete synchronization holds, leading to a better convergence of the estimator. Moreover, the numerical results are also affected by the tunned levels of process and measurement noise as well as the discretization model of the PF.

Figures 4c,d support this analysis with different arrangements. In Fig. 4c, the frequency mismatch over the network chain is relatively low, pointing to soft transitions that favours global synchronization behaviour in dynamical networks [6]. Thus, the tendency to form a single synchronous manifold for the whole network, rather than isolated clustered manifolds in Fig. 4a, might lead to observations of better quality for the sensor nodes and, hence, a lower overall $\eta$ as seen in Fig. 4c. Moreover, being the natural frequencies arranged in a homogeneous manner, the central node might have the worst estimates due to its large topological distance from both sensor nodes as well as the high dynamic behaviour differences between $v_{8}$ and $s_{1}$ or $s_{2}$. The $\eta_{i}$ is mostly unimodal, although some slight peaks in $\eta_{6}$ and $\eta_{10}$ might be attributed, respectively, to relative soft transitions on pairs $(v_{6},v_{7})$ and $(v_{10},v_{11})$. Likewise, in Fig. 4d, the worst estimates are central, due to high topological distance to sensor nodes and high dynamic behaviour differences.

Using the oscillators location along the chain graph of Fig. 3, we show in Fig. 5 how the PF performance changes with the sensor locations. For instance, in Fig. 5a, we set a pair of sensor nodes which provide the smallest path length from all other nodes in the network, while, in Fig. 5b, we set $(v_{4},v_{10})$ as the sensor nodes. Since $\eta$ is smaller in Fig. 5b, we argue that choosing a set of sensor nodes that is not only central to the network topology, but also “dynamically central” (representative of its neighbourhood dynamics) seems to convey better degrees of observability. At least for the Kuramoto dynamical network, presence of synchronous manifolds, dynamical correlations and topological proximity seem to affect the network system observability. Nevertheless, the interplay between synchronization, observability and dynamics has already been addressed in Ref. [29], with systems of lower dimensionality. The work pointed to a not exclusive dependence between synchronization and observability, a remark that seems to be also relevant in high-dimensional systems.

Finally, in Fig. 6, we investigate the improvement in the PF performance $\eta$ by gradually increasing the number of sensor nodes. Sensor nodes are allocated in such a way that the sum of the path lengths between every single node and the closest sensor node is minimum. As the number of sensors increase, the PF usually benefits from more valuable information and, consequently, $\eta$ has a tendency to decrease in an exponential manner. This is related to Haber et al. [23] conclusion that there is an intrinsic obstacle when trying to infer the network states from a small number of sensors. Interestingly, in some instances, $\eta$ increases after adding a sensor node. Due to the sensor nodes allocation method, from an even number of nodes to an odd number, not only the number of sensors increase but their positions change as well. Thus, increasing the number of sensor nodes does not guarantee better average observability. This result is also present in Fig. 2, where some sensor sets of cardinality $2$ provide better observability than a set of cardinality $3$. For instance, Fig. 6b shows clearly that the PF performance, and the network observability, are better when not including the central node as a sensor node, which always happens when there is an odd number of sensors. Interestingly, the higher the observability, the lower the median and the interquartile range of $\eta$—indicating that the PF estimates are more reliable.

4.3 Network of Rössler systems

Consider a Rössler dynamical network of $N=15$ nodes ($n=3N$ states) where measures are taken independently on the chain extremities ($\mathcal{S}=\{v_{1},v_{15}\}$, see Fig. 3). In this network system, each node is composed by three state variables whose observability rank is $y\rhd x\rhd z$, with $y$ and $x$ being good variables to reconstruct the dynamics, and $z$ a poor one [30]. In this section, based on the PF performance, we investigate how the observability of nodal dynamics affects the network observability. Thus, we compare the numerical results for the following sets of sensor nodes: (i) all the variables of a sensor node are measured, i.e. $\mathcal{S}_{1}=\{(x_{1},y_{1},z_{1}),(x_{15},y_{15},z_{15})\}$; and (ii) only one variable of a sensor node is measured i.e. $\mathcal{S}_{2}=\{x_{1},x_{15}\}$, $\mathcal{S}_{3}=\{y_{1},y_{15}\}$ and $\mathcal{S}_{4}=\{z_{1},z_{15}\}$.

Fig. 7 presents the attained filter estimates $\eta_{i}$ for each of the described sets of sensor nodes. As described in Section 3.2, the implemented Rössler system network has a soft transition in parameter $a_{i}$ along the network chain. Although it is not our goal to evaluate the effects of topology in network observability in this section, we note that the central node $v_{8}$ has the highest $\eta_{i}$ in every occasion, just as seen in Fig. 4c.

Clearly, the observability order of Rössler system is carried over to its high-dimensional network, as $\eta_{y}<\eta_{x}<\eta_{z}$. The poor performance of $\eta_{z}$ is likely associated with the little information (poor observability) that the $z$ variable has regarding the system cyclic behaviour, which hampers state reconstruction, especially when $z_{i}\approx 0$. Curiously, as the path length from a node estimate to a sensor node increases, $\eta_{z_{i}}$ reduces rather than increasing as in the other examples. This counter-intuitive phenomenon is yet to be explained. Despite this, all estimates $\eta_{z_{i}}$ in Fig. 7d are always inferior compared to the others. Moreover, the result that $\eta_{y}<\eta_{xyz}$, although only slightly, brings us to a recurring conclusion found in Ref. [41] that, in some cases, to measure all variables is not necessarily better than to measure the best one in terms of observability.

The attained effects of a multidimensional node over the network observability in this study can be studied under a graph point-of-view [5]. A single multidimensional node can be decomposed in multiple additional and interconnected nodes (states) as illustrated in Fig. 1. In the case of a Rössler system, based on system model (15), a single node $v_{i}$ can be “unfolded” into three nodes: $(v_{x_{i}},v_{y_{i}},v_{z_{i}})$, where their interconnections can be deduced from the system equation following Lin’s structural controllability approach [33]. The reader is referred to Refs. [32, 5] for further details. Nevertheless, notice that $v_{z_{i}}$ has a direct path to $v_{x_{i}}$ that is weighted by the state $z_{i}$. Thus, when $z_{i}\rightarrow 0$, this path vanishes and $v_{z_{i}}$ loses information of $v_{x_{i}}$ and, consequently, of all network. This behaviour is probably another way of understanding the loss of performance in Fig. 7d.

5 Conclusion

In this paper, we develop a particle filtering framework as a way to assess observability properties of a dynamical network, where each node is composed by an individual dynamical system. This is motivated by the fact that the quality of the filter estimation depends on the “level” of observability provided by the signal. Firstly, we show the PF effectiveness to quantify the observability level in a network by comparing this framework to a traditional observability metric in a system of lower dimensional order. Secondly, we show, in a network of $15$ and $45$ states (composed by Kuramoto oscillators and Rössler systems, respectively), that a good choice of sensors is related not only to the topological position of a sensor node in a graph, but also to its dynamic behaviour—that is, the precise values of the model parameters. These results corroborate with Ref. [23]. We conjecture that the sensor node dynamics must share some dynamical affinity to its neighbours in order to be representative of their behaviour and, thus, carry relevant information to the measurement signals. Thirdly, we investigate how the choice of an internal measured variable of a sensor node can affect the overall network observability. As shown in a Rössler network, the PF performance, and thus the observability, is better when measures are taken on variables of higher observability of the uncoupled node dynamics.

This work might provide insight for future researches on dynamical networks and the interplay between topology and dynamics. The detailed framework of Bayesian filtering as a performance measure of observability can be used as a validation platform over different benchmarks of networks in order to compare novel observability metrics in literature. This framework, however, is dependent on an appropriate tuning of the filter internal parameters, such as the process and measurement noise levels. This feature is also present in other observability metrics that relies on the fitting error of a state reconstructor [22, 8].

Appendix A Observability and particle filtering

of chaotic systems benchmarks

This section aims to validate, along with Section 4.1 results, the PF framework effectiveness—when well tuned—as an observability metric. To this purpose, we investigate the observability degree $\delta_{o}$ and the overall PF estimate error $\eta$ on five 3-dimensional chaotic systems benchmarks explored in [4]. Table 1 presents the observability rank derived by $\delta_{o}$ and $\eta$ (for a 100 Monte Carlo iterations) when variables $x$, $y$ and $z$ are recorded one at a time. For details on the chaotic systems ordinary differential equations as well as chosen parameters, we refer the reader to Ref. [4].

All numerical integrations were performed using a fourth-order Runge Kutta algorithm for a simulation time of $200$ s, with time step $0.01$ s (except for the cord system, which we set $0.001$ s). For all chaotic benchmarks, the PF framework follows Algorithm 2. Since in some chaotic systems the state variables magnitude range are significantly different, random noise $\bm{w}_{k}\in\mathbb{R}^{3}$ and $\bm{v}_{k}\in\mathbb{R}^{3}$ are independent, but not identical, random variables following a normal distribution $\mathcal{N}(0,\Sigma)$. The diagonal entries of $\Sigma_{w}$ and $\Sigma_{v}$ are shown in Table 1. In order to tune the PF for each chaotic system benchmark, $\Sigma_{w}$ and $\Sigma_{v}$ were adjusted following a model validation testing described in Ref. [7, Section 7.6.4], which is grounded on a whiteness test analysis for a PF framework. This model validation test was not applied in Section 3 since it is not quite robust for high-dimensional systems.

All results are statistically consistent, as shown by the quartile ranges. In the chaotic systems of uttermost observability properties, the Rössler and Cord systems, both metrics $\delta_{o}$ and $\eta$ ranked the variables from best to worst degree of observability in the same order. In the double-scroll attractor, both metrics indicate $y$ as the variable that conveys worse observability—which is expected since this variable does not distinguish the three equilibrium points. Despite the disagreement regarding the $x$ and $z$ variables, both convey good observability. The Lorenz system provides an interesting result. It is known that the $z$ variable does not distinguish between the two different equilibrium points on each of the attractor wings. Since this information is not detectable by a local measure, such as $\delta_{o}$, it ranks $z$ as variable of good observability[28]. The PF estimation error, on the other hand, detects this global property, inferring that $z$ is the variable of worst observability. This result has also been detected using the training error of a reservoir computer[8]. Finally, albeit the degree of observability conveyed by $x$ in the Lorenz’84 is different for each method, both recognizes that the $y$ and $z$ variables provides similar results.

Acknowledgements

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001, and by the Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brasil (CNPq).