Open AccessArticle

Determining Asymptotic Stability and Robustness of Networked Systems

Elias August

Department of Engineering, Reykjavik University, Menntavegur 1, 102 Reykjavik, Iceland

Systems 2020, 8(4), 39; https://doi.org/10.3390/systems8040039

Submission received: 27 August 2020 / Revised: 5 October 2020 / Accepted: 26 October 2020 / Published: 29 October 2020

(This article belongs to the Section Complex Systems and Cybernetics)

Download

Browse Figures

Versions Notes

Abstract

This paper is motivated by the notion that coupling systems allows for mitigating the failure of individual ones. We present a novel approach to determining asymptotic stability and robustness of a network consisting of coupled dynamical systems, where individual system dynamics are represented through polynomial or rational functions. The analysis relies on a local analysis; thus, making it computationally implementable. We present an efficient computational method that relies on semidefinite programming. Importantly, for cases where multiple equilibrium points exist, we show how to determine regions around an asymptotically stable equilibrium point that bounds solutions. These regions increase when systems are coupled as we observe when applying the presented analysis framework to a mathematical model of a continuous stirred tank reactor. Importantly, the presented work has implications to other fields as well.

Keywords:

coupled nonlinear systems; stability; robustness; semidefinite programming; sum of squares decomposition; continuous stirred tank reactor

1. Introduction

We are currently experiencing rapid digitalisation and automatisation in industry and society, which are culminating in the notions of smart factories, smart cities, the industry of things, and the internet of things. This makes networked systems and their control increasingly gaining in importance [1,2,3]. Significantly, networked sensor and actuator systems are also already pervasive in many traditional industries, for instance, the petroleum industry [4]. Additionally, networks of coupled dynamical systems appear also in other areas such as modelling of the heart [5] and of flocking behaviour [2]. An important question in network science, which links the dynamics of individual nodes with the emergent properties of macroscopic networks [6], is how much failure a networked system can tolerate. From the perspective of control engineering, a networked system can be viewed as a system of coupled dynamical systems and the question can be formulated in terms of asymptotic stability, the ability to return to a steady state, and robustness, the ability to withstand failures and perturbations.

In recent years, tools have been developed to determine whether a given complex network, consisting of connected systems, is robustly controllable, that is, whether any desired state can be reached from any initial state for almost all values of network edge-weights. For linear systems, conditions based on structural controllability, which identifies the minimal number of inputs or controllable nodes for system controllability of a complex network, were developed in [7,8]. For nonlinear systems, an equivalent approach based on feedback vertex set control, which is an attractor-based control method and requires thorough system knowledge, has been presented in [9]. However, these works and the majority of related ones rather consider systems that do not represent dynamical processes, as highlighted in [10], where the authors studied networks with node behaviour, which better represents dynamical behaviour that is typical for many systems, and examined how this influences the properties of the network. Using their analysis framework, they generated more realistic, albeit still simplified, results.

In this work, we provide certificates for asymptotic stability and robustness of complex networked systems consisting of coupled realistic system models through local analysis, which makes the approach computationally implementable and does not require simplifying the individual node dynamics. For this, we use so-called barrier certificates [11] and semidefinite programming [12,13]. For dynamical systems, the authors of [11] introduced a safety verification method that relies on state functions termed barrier certificates. Such functions form contour levels and, thus, create different regions in the state space. For a given set of initial conditions, the zero-level set of a barrier certificate bounds corresponding system trajectories and, by doing so, separates safe from unsafe regions. This means that one obtains certificates or proofs of system safety without the need to compute system behaviour explicitly. The use of different types of barrier certificates is popular in fields as diverse as biological model-testing, using the noise inherent in such systems [14], and avoiding collisions in multi-robot systems [15]. For our results, barrier certificates provide means to assess the set of initial conditions that will converge towards a desired operating point.

Related work has been presented in [16], where global stability for coupled differential equations was investigated assuming a Lyapunov function [17] of the form

V (x) = x - x^{*} ln x

, where

x^{*}

is a positive equilibrium point. This approach was later extended to coupled systems with time-varying coupling [18]. As a novelty, the approach presented in this paper provides not only means to search for a Lyapunov function that guarantees asymptotic stability of the coupled system, for systems with multiple equilibria, it also uses Lyapunov theory to provide regions of attraction, even when model parameters are uncertain (robustness).

In many areas of research, realistic dynamical system models are based on ordinary differential equations described by polynomial or rational functions. Assuming such models, in Section 2.1 and Section 2.2, we provide results on testing for asymptotic stability and robustness of a networked system through a local analysis of the system. Locality makes the analysis feasible and allows for its computational implementation, as we show in Section 2.3. Meanwhile, for simplicity, we first assume all-to-all coupling, in Section 2.4, we extend the results to arbitrary coupling topologies. Section 3 provides an application to process engineering, where we analyse the stability and robustness of a high biomass concentration steady state for connected continuous stirred tank reactors. Finally, we discuss our results and conclude the paper in Section 4.

2. Asymptotic Stability of Coupled Dynamical Systems

In this paper, we consider bidirectional coupling only. Furthermore, let us consider N dynamical systems that are linearly all-to-all coupled. Then, we represent the i-th system through

{\dot{x}}_{i} = f_{i} (x_{i}) + \sum_{j = 1}^{N} k D (x_{j} - x_{i}),

(1)

where

f_{i} : R^{n} \to R^{n}

describes the system’s dynamics for all i,

i \in {1, 2, \dots, N}

x = {[x_{1} \dots x_{N}]}^{T}

x_{i} (t) \in R^{n}

, and

D \in R^{n \times n}

. Vector x is the state vector. Matrix D is the matrix that denotes the variables that are used in the coupling; that is,

D_{(k, l)} \neq 0

if and only if states

x_{i_{(k)}}

and

x_{j_{(l)}}

are coupled, for all

i, j

and

i \neq j

. Positive constant k corresponds to the coupling strength. Next, we provide results on asymptotic stability of the coupled system based on Lyapunov theory.

2.1. Unique Equilibrium Point

First, let us assume that the networked system described by N coupled systems given by Equation (1) has a unique global equilibrium point. Without loss of generality, let it be the origin; otherwise, the equilibrium point can be shifted to the origin through a change of variables. Moreover, we assume that

{\dot{x}}_{1} = f_{1} (x_{1})

is asymptotically stable. We do this, because, without proof, the networked system can only be asymptotically stable if at least one individual system is asymptotically stable, whose state we denote by

x_{1}

, without loss of generality. If we consider all other systems to be similar to this system but to possess a certain amount of bounded parameter uncertainty then condition (4), which we derive in the following, guarantees that the coupled system is asymptotically stable. This allows us to assess the robustness of a networked system.

Now, consider the following Lyapunov function,

V (x) = \sum_{i = 2}^{N} (\frac{1}{2} \frac{1}{N - 1} x_{1}^{T} M x_{1} + \bar{V} (X_{i})) = \sum_{i = 2}^{N} \tilde{V} (x_{1}, X_{i}),

(2)

where

X_{i} = x_{i} - x_{1}

M = M^{T} ≻ 0

, and

\tilde{V} (x_{1}, X_{i}) \geq 0

. Moreover,

\tilde{V} (x_{1}, X_{i}) > 0

if neither

x_{1} = 0

nor

X_{i} = 0

, for all i. Thus,

V (x) > 0

unless

x = 0

. Then,

\begin{matrix} \dot{V} (x) & = & \sum_{i = 2}^{N} (\frac{1}{2} \frac{1}{N - 1} ({\dot{x}}_{1}^{T} M x_{1} + x_{1}^{T} M {\dot{x}}_{1}) + \dot{\bar{V}} (X_{i})) = \sum_{i = 2}^{N} (\frac{1}{N - 1} x_{1}^{T} M {\dot{x}}_{1} + \dot{\bar{V}} (X_{i})) \\ = & \sum_{i = 2}^{N} (\frac{1}{N - 1} x_{1}^{T} M (f_{1} (x_{1}) + k D \sum_{j = 2}^{N} X_{j}) + \frac{\partial \bar{V} (X_{i})}{\partial X_{i}} (f_{i} (x_{i}) - f_{1} (x_{1}) - k N D X_{i})) \\ = & \sum_{i = 2}^{N} (\frac{1}{N - 1} x_{1}^{T} M f_{1} (x_{1}) + x_{1}^{T} k M D X_{i} + \frac{\partial \bar{V} (X_{i})}{\partial X_{i}} (f_{i} (x_{i}) - f_{1} (x_{1}) - k N D X_{i})) . \end{matrix}

(3)

The second equality in the first line of Equation (3) follows from the symmetry of M. Thus, if

\frac{1}{N - 1} x_{1}^{T} M f_{1} (x_{1}) + x_{1}^{T} k M D X_{i} + \frac{\partial \bar{V} (X_{i})}{\partial X_{i}} (f_{i} (x_{i}) - f_{1} (x_{1}) - k N D X_{i}) < 0

(4)

holds for all i unless

x_{1} = 0

and

X_{i} = 0

for all i, and, thus, unless

x = 0

, then

\dot{V} < 0

unless

x = 0

, which guarantees asymptotic stability of the coupled system. Finally, for all i, we consider the differences in function

f_{i} (x_{i})

to depend on r real parameters that are uncertain but whose values are bounded. Let these parameters be elements of vector

q \in Q \subset R^{r}

. Then, we can rewrite condition (4) by letting

g (q_{1}, x_{1}) \equiv f_{1} (x_{1})

and

g (q_{i}, x_{i}) \equiv f_{i} (x_{i})

. In this case, if condition (4) holds, for all

x_{1}, x_{i} \in R^{n}

and all

q_{1}, q_{i} \in Q

, then asymptotic stability of the coupled system is guaranteed.

2.2. Multiple Equilibrium Points

Second, if the networked system possesses multiple equilibrium points, then inequality (4) can only hold in a certain region of the state space. In order to estimate this region, we assume, again, that the equilibrium point, whose stability properties we wish to investigate, is the origin. Then, if inequality (4) holds for

\tilde{V} (x_{1}, X_{i}) < c

c \in R

c > 0

, and all i then

\tilde{V} (x_{1}, X_{i})

decreases with time for all i and, thus, so does

V (x)

, which means that the origin is asymptotically stable in the region bounded by

V (x) = c

. As we will show in Section 3, in addition to above robustness analysis, determining this region allows us to assess how far a system can deviate from a particular steady state and still remain attracted by the latter.

2.3. Computational Implementation

2.3.1. Coupled Linear Systems

Consider N coupled linear dynamical systems of the form

{\dot{x}}_{1} = A x

and

{\dot{x}}_{i} = B x

for all

i > 1

, where

A, B \in R^{n \times n}

, matrix A is stable, and matrix B is a matrix that is similar to A, however, with some uncertain but bounded entries. Then, for

\bar{V} (X_{i}) = X_{i}^{T} P X_{i}

, where

P = P^{T} ≻ 0

, inequality (4) becomes

\begin{matrix} x_{1}^{T} (M (\frac{1}{N - 1} A - k D) + P (A - k N D)) x_{1} \\ + x_{1}^{T} (k M D + P (k N D - B)) x_{i} + x_{i}^{T} P (k N D - A) x_{1} + x_{i}^{T} P (B - k N D) x_{i} < 0 . \end{matrix}

(5)

Given k, the problem of determining matrices M and P such that Equation (2) and inequality (5) hold can be cast as a semidefinite programme of the following form

\begin{matrix} given & A, B, D, N, k \\ search for & M, P \\ subject to & M = M^{T} ≻ 0, P = P^{T} ≻ 0, (5) . \end{matrix}

(6)

Problem (6) is a classic case of linear matrix inequalities [19] that can be solved efficiently using, for example, YALMIP, a free, third-party MATLAB toolbox [20], which can handle bounded uncertainties in B [21] to obtain robustness certificates.

2.3.2. Coupled Nonlinear Systems

For nonlinear systems, we investigate those ones whose dynamics are modelled using polynomial or rational functions. Consider the real-valued polynomial function

F (x)

of degree

2 d

x \in R^{n}

. Testing for non-negativity is NP-hard [22]. However, a sufficient condition for

F (x)

to be nonnegative is that it can be decomposed into a sum of squares (SOS) [23]:

F (x) = \sum_{i} f_{i}^{2} (x) \geq 0,

(7)

where

f_{i}

are polynomial functions. Now,

F (x)

is a SOS if and only if there exists a positive semidefinite matrix R and

F (x) = \sum_{i} f_{i}^{2} (x) = χ^{T} R χ, χ = [1, x_{(1)}, x_{(2)}, \dots, x_{(n)}, x_{(1)} x_{(2)}, \dots, x_{(n)}^{d}] .

(8)

The entries of vector

χ

consist of all monomial combinations of the elements of vector x up to degree d (including

x_{(i)}^{0} = 1

) and, thus, its length is

ℓ = (\binom{n + d}{d})

. Note that R is not necessarily unique. However, Equation (8) poses certain constraints on R of the form trace

(A_{j} R) = c_{j}

, where

A_{j}

and

c_{j}

are appropriate matrices and constants respectively. As an illustration, in the following, we reproduce Example 3.5 of [23].

Consider the following polynomial function, where

z_{1} = x_{1}^{2}

z_{2} = x_{2}^{2}

z_{3} = x_{1} x_{2}

\begin{matrix} F (x) & = & 2 x_{1}^{4} + 2 x_{1}^{3} x_{2} - x_{1}^{2} x_{2}^{2} + 5 x_{2}^{4} \\ = & {[\begin{matrix} x_{1}^{2} \\ x_{2}^{2} \\ x_{1} x_{2} \end{matrix}]}^{T} [\begin{matrix} q_{11} & q_{12} & q_{13} \\ q_{12} & q_{22} & q_{23} \\ q_{13} & q_{23} & q_{33} \end{matrix}] [\begin{matrix} x_{1}^{2} \\ x_{2}^{2} \\ x_{1} x_{2} \end{matrix}] \\ = & q_{11} x_{1}^{4} + q_{22} x_{2}^{4} + (q_{33} + 2 q_{12}) x_{1}^{2} x_{2}^{2} + 2 q_{13} x_{1}^{3} x_{2} + 2 q_{23} x_{1} x_{2}^{3} . \end{matrix}

This leads to the following set of linear equalities:

q_{11} = 2, q_{22} = 5, q_{33} + 2 q_{12} = - 1, 2 q_{13} = 2, 2 q_{23} = 0 .

Thus,

j = 1, \dots, 5

, and, for example,

c_{1} = 2

c_{3} = - 1

A_{1} = [\begin{matrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}], and A_{3} = [\begin{matrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{matrix}] .

In general, in order to find R, we solve the optimisation problem associated with the following semidefinite program:

\begin{matrix} minimise & trace (A_{0} R) \\ subject to & trace (A_{j} R) = c_{j}, j = 1, \dots, m \\ R = R^{T} ⪰ 0 . \end{matrix}

(9)

In this paper, to solve SOS programmes, we use SOSTOOLS [24], a free, third-party MATLAB toolbox. For instance, the problem of determining Lyapunov function

\tilde{V} (x_{1}, X_{i})

, given by Equation (2), such that (4) holds can be cast as follows

\begin{matrix} given & f_{1} (x_{1}), f_{i} (x_{i}), D, N, k \\ search for & \tilde{V} (x_{1}, X_{i}) \\ subject to & - ε (x_{1}^{T} x_{1} + x_{i}^{T} x_{i}) + \tilde{V} (x_{1}, X_{i}) is SOS, \\ - ε (x_{1}^{T} x_{1} + x_{i}^{T} x_{i}) - \dot{\tilde{V}} (x_{1}, X_{i}) is SOS, 0 < ε ≪ 1 . \end{matrix}

(10)

Here, some additional remarks:

Consider a rational function $F (x)$ , $F (x) = \frac{f (x)}{g (x)}$ , where $f (x)$ and $g (x)$ are polynomial functions. Then, $F (x) \geq 0$ if (9) is feasible with $χ^{T} R χ = F (x) g^{2} (x)$ or with $χ^{T} R χ = F (x) g (x)$ if $g (x) > 0$ .
If

$\begin{matrix} F (x) + p (x) h (x) = \sum_{i} g_{i}^{2} (x) \geq 0, p (x) \geq 0, \\ h (x) = \{\begin{matrix} \leq 0 if a_{i} \leq x_{(i)} \leq b_{i} \forall i \\ > 0 otherwise \end{matrix}, \end{matrix}$

then $F (x) \geq 0$ if $a_{i} \leq x_{(i)} \leq b_{i}$ for all i, where $a_{i}, b_{i}$ are constants. This can be used to show that $F (x)$ is nonnegative in a specific region of the state and/or parameter space.
Finally, the computational time necessary to solve SOS decomposition problems scales badly with the size of the problem, as the length of vector is $ℓ = (\binom{n + d}{d})$ .

2.4. Different Coupling Configurations

The presented results can be extended to any symmetric linear coupling topology with corresponding Laplacian matrix

- C

by letting

\bar{V} (X_{i}) = X_{i}^{T} P X_{i}

, for all i, in Equation (2), where

P = P^{T} ≻ 0

. To show this, first, we define the Laplacian matrix of the completely connected graph given by

- U \in R^{N \times N}

through

- U = N I_{N} - E ⪯ 0, E = e e^{T}, e = {[1 \dots 1]}^{T} .

(11)

Note that the eigenvalues of

- U

are

λ_{1} = 0

and

λ_{i} = N

for all

i > 1 .

Then, the following lemma establishes a relation between C and U.

Lemma 1.

Let

U, C \in R^{N \times N}

, where

- U

is the Laplacian matrix of the completely connected graph defined in Equation (11) and

- C

is a symmetric Laplacian matrix of rank

(N - 1)

. Then, for all

y \in R^{N}

y^{T} C y \leq \frac{λ_{min} (- C)}{N} y^{T} U y,

(12)

where

λ_{min} (- C)

denotes the smallest positive eigenvalue of the matrix

- C

Proof.

First, note that

C_{(i, j \neq i)} \geq 0

. Now,

U E = 0

and

C E = 0

. Hence,

U U = - N U

and

C U = - N C

. By the Spectral Theorem, there exists a unitary matrix S such that

S^{- 1} = S^{T}

and

- S C S^{T} = Λ_{- C}

, where

Λ_{- C}

is the diagonal matrix containing the eigenvalues of

- C

and, without loss of generality, its first entry is 0. Moreover,

S C U S^{T} = S C S^{T} S U S^{T} = - S C S^{T} N = S C S^{T} Λ_{U}

implies that

S U S^{T} = Λ_{U}

, where

Λ_{U}

is the diagonal matrix containing the eigenvalues of U, which are

- N

but for the first one, which is 0. Now, for all

y \in R^{N}

and

Λ_{C}

being the diagonal matrix containing the eigenvalues of C,

y^{T} C y = y^{T} S^{T} Λ_{C} S y = v^{T} Λ_{C} v \leq \frac{λ_{min} (- C)}{N} v^{T} Λ_{U} v = \frac{λ_{min} (- C)}{N} y^{T} S^{T} Λ_{U} S y = \frac{λ_{min} (- C)}{N} y^{T} U y,

which completes the proof. □

Consider, again, difference

X_{i} \equiv x_{i} - x_{1}

, where

X_{i} \in R^{n}

and

X_{1} = 0

[25]. Let matrix

F \in R^{N \times N}

be defined as follows,

F = [\begin{matrix} 0 & 0 & 0 & \dots & 0 \\ - 1 & 1 & 0 & \dots & 0 \\ - 1 & 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ - 1 & 0 & 0 & \dots & 1 \end{matrix}] .

(13)

Then,

X = (F \otimes I_{n}) x, X, x \in R^{N n},

(14)

and

\dot{X} = (F \otimes I_{n}) \dot{x} .

(15)

Moreover, note that for symmetric coupling determined by Laplacian matrix

- C

C = C F

and

F U = - N F

. Let

\tilde{f} (x) = {[f_{1} (x_{1}) \dots f_{N} (x_{N})]}^{T}

. Then,

\begin{matrix} \sum_{j = 2}^{N} \dot{\bar{V}} (X_{i}) & = & X^{T} (I_{N} \otimes P) \dot{X} = X^{T} (F \otimes P) \dot{x} \\ = & X^{T} (F \otimes P) (\tilde{f} (x) + k (C \otimes D) x) \\ = & X^{T} (F \otimes P) (\tilde{f} (x) + k (C F \otimes D) x) \\ = & X^{T} (F \otimes P) (\tilde{f} (x) + k (C \otimes D) X) \\ \leq & X^{T} (F \otimes P) (\tilde{f} (x) + α k (U \otimes D) X) \\ = & \sum_{i = 2}^{N} X_{i}^{T} P (f_{i} (x_{i}) - f_{1} (x_{1})) + α k X^{T} (F U \otimes P D) X \\ = & \sum_{i = 2}^{N} X_{i}^{T} P (f_{i} (x_{i}) - f_{1} (x_{1})) - α N k X^{T} (F \otimes P D) X \\ = & \sum_{i = 2}^{N} X_{i}^{T} P (f_{i} (x_{i}) - f_{1} (x_{1}) - α k N D X_{i}), \end{matrix}

where the inequality follows from Lemma 1,

α N = λ_{min} (- {(F C)}_{gL})

, and

{(F C)}_{gL}

is the matrix obtained by deleting the first row and the first column of

F C

. It is the so-called grounded Laplacian [26,27], where the deleted row and deleted column correspond to the grounded variable—in our case it is

x_{1}

. A grounded Laplacian matrix is a principal sub-matrix of the Laplacian matrix. The fact that, for square matrices

A, B

eig (A B) = eig (B A)

, has been proven in [28] and, thus,

eig (F C) = eig (C F) = eig (C)

. Since the first row of matrix

F C

is zero, it follows that

λ_{min} (- {(F C)}_{gL}) = λ_{min} (- C)

. For illustration, we provide examples of

λ_{min} (- C)

of different unweighted coupling configurations.

(a): All-to-all coupling ( $C = U$ ): $λ_{min} (- C) = N$ .
(b): Star-configuration: $λ_{min} (- C) = 1$ .
(c): Ring of diffusively coupled systems: $λ_{min} (- C) = 4 {sin}^{2} \frac{π}{N}$ .
(d): Ring of $2 k$ -nearest neighbour coupled systems [29]: $λ_{min} (- C) ≃ 2 π^{2} k (k + 1) (2 k + 1) / 3 N^{2}$ if $k ≪ N$ .

The different coupling configurations are illustrated in Figure 1.

3. Coupled Continuous Stirred Tank Reactors

Consider a continuous stirred tank reactor for biomass production held at a constant temperature [30]. Concentrations of feed or substrate are denoted by s and those of the product by x (Figure 2). The reaction describing the inflow and outflow of substrate and outflow/harvesting of product is given by

x \to \emptyset ⇄ s

. The reaction describing the growth of product by consuming substrate is given by

x s \overset{k}{\to} 2 x

, where reaction constant k is described below.

The corresponding mathematical model based on ordinary differential equations is given by

\begin{matrix} \dot{x} & = & \frac{2 x s}{{(s + 1)}^{2}} - 0.4 x, \\ \dot{s} & = & u - \frac{2 x s}{{(s + 1)}^{2}} - 0.4 s, \end{matrix}

(16)

where we denote the input rate of substrate by u. Moreover, we assume an outflow rate of

0.4

(unitless) and a reaction rate of

k x s

, where

k = \frac{2}{{(s + 1)}^{2}}

implies inhibition for high concentrations of substrate (poisoning). We denote equilibrium points by

x^{*}

and

s^{*}

. Figure 3 shows the phase plot with the vector field for

u = 2

, which depicts the behaviour of the system, determined through the direction set by the vector arrows, in dependence of its state. For

u = 2

, there are three equilibrium points, which are clearly visible. The equilibrium point that corresponds to higher concentrations of product is given by

x^{*} = 4.618

and

s^{*} = 0.382

. Note that the system described by Equation (16) is a nonnegative system such that for nonnegative initial conditions, solutions remain in the nonnegative orthant [31]. To obtain a set that bounds solutions, we let

\bar{x} = x + x^{*}

and

\bar{s} = s + s^{*}

, first, solve the following problem

\begin{matrix} given & Equation (16), u = 2, β \in R \\ search for & M \\ subject to & p_{1} (\bar{x}, \bar{s}), p_{2} (\bar{x}, \bar{s}), p_{3} (\bar{x}, \bar{s}) are SOS functions of degree 4, \\ - 0.001 ({\bar{x}}^{2} + {\bar{s}}^{2}) + V_{1} is SOS, V_{1} = \frac{1}{2} {[\begin{matrix} \bar{x} \\ \bar{s} \end{matrix}]}^{T} M [\begin{matrix} \bar{x} \\ \bar{s} \end{matrix}], \\ - 0.001 ({\bar{x}}^{2} + {\bar{s}}^{2}) - {\dot{V}}_{1} \cdot {(s + 1)}^{2} + p_{1} (\bar{x}, \bar{s}) ({\bar{x}}^{2} + {\bar{s}}^{2} - β) - p_{2} (\bar{x}, \bar{s}) x - p_{3} (\bar{x}, \bar{s}) s is SOS, \end{matrix}

(17)

for different values of

β

, and, then, solve

\begin{matrix} given & V_{1} \\ minimise & - γ, γ \in R \\ subject to & p_{1} (\bar{x}, \bar{s}), p_{2} (\bar{x}, \bar{s}), p_{3} (\bar{x}, \bar{s}) are SOS functions of degree 2, \\ V_{1} - γ - p_{1} (\bar{x}, \bar{s}) ({\bar{x}}^{2} + {\bar{s}}^{2} - β) - p_{2} (\bar{x}, \bar{s}) x - p_{3} (\bar{x}, \bar{s}) s is SOS . \end{matrix}

(18)

The boundary of the set is given by

V_{1} = γ

. Figure 3 depicts the sets defined by

β = 7

and

β = 9.95

. Their boundaries are given by the dashed red lines. For

β > 9.95

, programme (17) is infeasible and remains so even after increasing the degree of polynomials

p_{1}

p_{2}

, and

p_{3}

. Note that outside the regions, whose boundaries are given by the dashed red line in Figure 3, one might risk washout or depletion of product by being attracted to the equilibrium point given by

x^{*} = 0

and

s^{*} = 5

Now, by interconnecting reactors as in Equation (1), we can mitigate the danger of washout that might result, for example, from disturbances. Note that in our modelling, we assume full control over all flow rates, that is, in all inlets, where we let

u = 2

, and between all reactors, where

k = 5

. For

N > 1

, we solve the following problem, which is analogous to problem (17),

\begin{matrix} given & Equation (16), Equation (1), u = 2, k = 5, D = I_{2}, N > 1, β \in R \\ search for & \tilde{V} ({\bar{x}}_{1}, {\bar{s}}_{1}, {\bar{x}}_{i}, {\bar{s}}_{i}) \\ subject to & p_{j} ({\bar{x}}_{1}, {\bar{s}}_{1}, {\bar{x}}_{i}, {\bar{s}}_{i}) are SOS functions of degree 4, j \in {1, \dots, 5}, \\ \tilde{V} = \frac{1}{2} \frac{1}{N - 1} x_{1}^{T} M x_{1} + \bar{V} (X_{i}) is a SOS function of degree 2, \\ - 0.001 ({\bar{x}}_{1}^{2} + {\bar{s}}_{1}^{2} + {\bar{x}}_{i}^{2} + {\bar{s}}_{i}^{2}) + \tilde{V} is SOS, \\ - 0.001 ({\bar{x}}_{1}^{2} + {\bar{s}}_{1}^{2} + {\bar{x}}_{i}^{2} + {\bar{s}}_{i}^{2}) - \dot{\tilde{V}} \cdot {({\bar{s}}_{1} + 1)}^{2} {({\bar{s}}_{i} + 1)}^{2} \\ + p_{1} \cdot ({\bar{x}}_{1}^{2} + {\bar{s}}_{1}^{2} + {\bar{x}}_{i}^{2} + {\bar{s}}_{i}^{2} - β) - p_{2} \cdot x_{1} - p_{3} \cdot s_{1} - p_{4} \cdot x_{i} - p_{5} \cdot s_{i} is SOS . \end{matrix}

(19)

By, first, solving a problem analogous to problem (18), with

V_{1} = \frac{1}{2} \frac{1}{N - 1} {[\begin{matrix} {\bar{x}}_{1} \\ {\bar{s}}_{1} \end{matrix}]}^{T} M [\begin{matrix} {\bar{x}}_{1} \\ {\bar{s}}_{1} \end{matrix}]

, for

β = 12

and

β = 17.5

, for different

N > 1

, and, then, assuming that all systems are at equilibrium but the first one, we find that solutions remain bounded and return to equilibrium if

x_{1}

and

s_{1}

are initially within the regions bounded by the black solid lines in Figure 3. Finally, for investigating the robustness of our previous results, we repeat above computation with additional free variables that account for uncertainty in model parameters and corresponding constraints. For instance, for an uncertain input, such as

2 \leq u \leq 2.153

, which translates to

4.618 \leq x^{*} \leq 5

x^{*}

becomes an additional free variable. We add the additional constraints to the last line of programme (19). They have the following form

+ p_{6} \cdot (x_{1}^{*} - 4.168) (x_{1}^{*} - 5) + p_{7} \cdot (x_{i}^{*} - 4.168) (x_{i}^{*} - 5),

where

p_{6} ({\bar{x}}_{1}, {\bar{s}}_{1}, {\bar{x}}_{i}, {\bar{s}}_{i})

and

p_{7} ({\bar{x}}_{1}, {\bar{s}}_{1}, {\bar{x}}_{i}, {\bar{s}}_{i})

are SOS functions of degree 4. We, again, observe that the regions that bound solutions are given by the black solid lines.

4. Discussion and Conclusions

In this paper, we presented a novel approach for determining asymptotic stability and robustness of a network consisting of coupled dynamical systems, where individual system dynamics are represented through polynomial or rational functions. The analysis framework relies on local analysis; thus, making it computationally implementable. With respect to the latter, we presented an efficient method that relies on semidefinite programming. Importantly, for cases, where multiple equilibrium points exist, we show how to determine regions around an asymptotically stable equilibrium point that bounds solutions. The novelty of the presented research lies in the analysis of complex networks that consist of realistic nonlinear systems, which are of relevance to science and industry, by employing barrier certificates to provide guarantees for asymptotic stability and robustness, as opposed to previous research, which does rather not consider realistic dynamical systems as network nodes.

For example, for biomass production in continuous stirred tank reactors, on the one hand, one wishes to increase yield by providing large concentrations of substrate; on the other hand, too much of it can lead to a decrease or even washout of product. Intuitively, coupling a few reactors should allow for greater deviations of individual reactors from the optimal operating point without risking washout. Using a model of such a networked system, we showed how to determine just how much deviation is (surely) safe. Importantly, assuming all reactors at the optimal operating point but one, we showed that deviations deemed unsafe for a single reactor are now still within the safe zone. Finally, wind farms, the industrial internet of things, and the spread of epidemics, have all in common complex dynamical behaviour because of the interconnectedness of an increasingly large number of individual systems or individuals. However, the fact that they consist of complex networks makes achieving the goal of robustly controlling such systems particularly challenging. The research presented in this paper provides an approach to overcoming this challenge.

Funding

This research received no external funding.

Conflicts of Interest

The author declares no conflict of interest.

Notation

$R, R^{m \times n}$	real numbers, $m \times n$ real matrices
$A_{(i, j)}$	$(i, j)$ th entry of matrix $A \in R^{m \times n}$
$I_{n}$	the identity matrix, $I_{n} \in R^{n \times n}$
$A^{T}$	transpose of matrix $A \in R^{m \times n}$
$\dot{x}$	derivative of x with respect to time variable t
$A ≻ 0$ , $B ⪰ 0$	matrix A is positive definite, matrix B is positive semidefinite
$A ≺ 0$ , $B ⪯ 0$	matrix A is negative definite, matrix B is negative semidefinite
$A \otimes B$ , $A \in R^{m \times n}$ , $B \in R^{p \times q}$	The Kronecker product: $[\begin{matrix} A_{(1, 1)} B & \dots & A_{(1, n)} B \\ A_{(2, 1)} B & \dots & A_{(2, n)} B \\ ⋮ & ⋱ & ⋮ \\ A_{(m, 1)} B & \dots & A_{(m, n)} B \end{matrix}]$

References

Bulatov, Y.N.; Kryukov, A.V. A multi-agent control system of distributed generation plants. In Proceedings of the International Conference on Industrial Engineering, Applications and Manufacturing (ICIEAM), St. Petersburg, Russia, 16–19 May 2017. [Google Scholar]
Olfati-Saber, R.; Fax, J.A.; Murray, R.M. Consensus and Cooperation in Networked Multi-Agent Systems. Proc. IEEE 2007, 95, 215–233. [Google Scholar] [CrossRef] [Green Version]
Ge, X.; Yang, F.; Han, Q.L. Distributed networked control systems: A brief overview. Inf. Sci. 2017, 380, 117–131. [Google Scholar] [CrossRef]
Pajic, M.; Sundaram, S.; Pappas, G.J.; Mangharam, R. The wireless control network: A new approach for control over networks. IEEE Trans. Autom. Control 2011, 56, 2305–2318. [Google Scholar] [CrossRef] [Green Version]
dos Santos, A.M.; Lopes, S.R.; Viana, R.L. Rhythm synchronization and chaotic modulation of coupled Van der Pol oscillators in a model for the heartbeat. Physica A 2004, 338, 335–355. [Google Scholar] [CrossRef]
Vespignani, A. Twenty years of network science. Nature 2018, 558, 528–529. [Google Scholar] [CrossRef]
Liu, Y.Y.; Slotine, J.J.; Barabási, A.L. Controllability of complex networks. Nature 2011, 473, 167–173. [Google Scholar] [CrossRef]
Liu, Y.Y.; Barabási, A.L. Control principles of complex systems. Rev. Mod. Phys. 2016, 88, 035006. [Google Scholar] [CrossRef] [Green Version]
Zanudo, J.G.T.; Yang, G.; Albert, R. Structure-based control of complex networks with nonlinear dynamics. Proc. Natl. Acad. Sci. USA 2017, 114, 7234–7239. [Google Scholar] [CrossRef] [Green Version]
Leitold, D.; Vathy-Fogarassy, Á.; Abonyi, J. Controllability and observability in complex networks—The effect of connection types. Sci. Rep. 2017, 7, 151. [Google Scholar] [CrossRef]
Prajna, S.; Jadbabaie, A.; Pappas, G. Stochastic Safety Verification Using Barrier Certificates. In Proceedings of the IEEE Conference on Decision and Control (CDC) (IEEE Cat. No.04CH37601), Nassau, Bahamas, 14–17 December 2005; pp. 929–934. [Google Scholar]
Vandenberghe, L.; Boyd, S. Semidefinite programming. SIAM Rev. 1996, 38, 49–95. [Google Scholar] [CrossRef] [Green Version]
Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, UK, 2004. [Google Scholar]
August, E. Using Noise for Model-Testing. J. Comput. Biol. A J. Comput. Mol. Cell Biol. 2012, 19, 968–977. [Google Scholar] [CrossRef] [Green Version]
Wang, L.; Ames, A.; Egerstedt, M. Safety Barrier Certificates for Collisions-Free Multirobot Systems. IEEE Trans. Robot. 2017, 33, 661–674. [Google Scholar] [CrossRef]
Li, M.; Shuai, Z. Global-stability problem for coupled systems of differential equations on networks. J. Differ. Equ. 2010, 248, 1–20. [Google Scholar] [CrossRef] [Green Version]
Khalil, H.K. Nonlinear Systems, 3rd ed.; Prentice-Hall: Upper Saddle River, NJ, USA, 2000. [Google Scholar]
Feng, J.; Luan, X. Stability analysis for coupled systems with time-varying coupling structure. Adv. Differ. Equ. 2017, 2017, 191. [Google Scholar] [CrossRef] [Green Version]
Boyd, S.; El Ghaoui, L.; Feron, E.; Balakrishnan, V. Studies in Applied Mathematics. In Linear Matrix Inequalities in System and Control Theory; SIAM: Philadelphia, PA, USA, 1994; Volume 15. [Google Scholar]
Löfberg, J. YALMIP: A Toolbox for Modeling and Optimization in MATLAB. In Proceedings of the CACSD Conference, New Orleans, LA, USA, 2–4 September 2004. [Google Scholar]
Löfberg, J. Automatic robust convex programming. Optim. Methods Softw. 2012, 27, 115–129. [Google Scholar] [CrossRef] [Green Version]
Murty, K.G.; Kabadi, S.N. Some NP-complete problems in quadratic and nonlinear programming. Math. Program. 1987, 39, 117–129. [Google Scholar] [CrossRef]
Parrilo, P.A. Semidefinite programming relaxations for semialgebraic problems. Math. Program. Ser. B 2003, 96, 293–320. [Google Scholar] [CrossRef]
Papachristodoulou, A.; Anderson, J.; Valmorbida, G.; Prajna, S.; Seiler, P.; Parrilo, P.A. SOSTOOLS—Sum of Squares Optimization Toolbox, User’s Guide. 2018. Available online: http://sysos.eng.ox.ac.uk/sostools/ (accessed on 29 October 2020).
Belykh, V.N.; Belykh, I.V.; Hasler, M. Connection graph stability method for synchronized coupled chaotic systems. Physica D 2004, 195, 159–187. [Google Scholar] [CrossRef]
Miekkala, U. Graph properties for splitting with grounded laplacian matrices. BIT 1993, 33, 485–495. [Google Scholar] [CrossRef]
Jadbabaie, A.; Motee, N.; Barahona, M. On the Stability of the Kuramoto Model of Coupled Nonlinear Oscillators. In Proceedings of the 2004 American Control Conference, Boston, MA, USA, 30 June–2 July 2004; pp. 4296–4301. [Google Scholar]
Thurston, H.S. On the Characteristic Equations of Products of Square Matrices. Am. Math. Mon. 1931, 38, 322–324. [Google Scholar] [CrossRef]
Barahona, M.; Pecora, L.M. Synchronization in Small-World Systems. Phys. Rev. Lett. 2002, 89, 054101. [Google Scholar] [CrossRef] [Green Version]
Chmiel, H. Bioprozesstechnik, 3rd ed.; Springer: Heidelberg, Germany, 2001. [Google Scholar]
Haddad, W.M.; Chellaboina, V.; August, E. Stability and Dissipativity Theory for Nonnegative Dynamical Systems: A Thermodynamic Framework for Biological and Physiological Systems. In Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No.01CH37228), Orlando, FL, USA, 4–7 December 2001; pp. 442–458. [Google Scholar]

Figure 1. Four different coupling configurations. Here,

N = 10

in (a,c,d),

N = 11

in (b), and

k = 2

in (d).

Figure 1. Four different coupling configurations. Here,

N = 10

in (a,c,d),

N = 11

in (b), and

k = 2

in (d).

Figure 2. Schematic view of a continuous stirred tank reactor.

Figure 3. Phase plot with vector field and corresponding contour lines, for the system described by Equation (16), and regions that bound solutions (red dashed lines). For coupling of a few such systems as in Equation (1) and for

N > 1

, solutions are bounded by the black solid lines.

N > 1

, solutions are bounded by the black solid lines.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

August, E. Determining Asymptotic Stability and Robustness of Networked Systems. Systems 2020, 8, 39. https://doi.org/10.3390/systems8040039

AMA Style

August E. Determining Asymptotic Stability and Robustness of Networked Systems. Systems. 2020; 8(4):39. https://doi.org/10.3390/systems8040039

Chicago/Turabian Style

August, Elias. 2020. "Determining Asymptotic Stability and Robustness of Networked Systems" Systems 8, no. 4: 39. https://doi.org/10.3390/systems8040039

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Determining Asymptotic Stability and Robustness of Networked Systems

Abstract

1. Introduction

2. Asymptotic Stability of Coupled Dynamical Systems

2.1. Unique Equilibrium Point

2.2. Multiple Equilibrium Points

2.3. Computational Implementation

2.3.1. Coupled Linear Systems

2.3.2. Coupled Nonlinear Systems

2.4. Different Coupling Configurations

3. Coupled Continuous Stirred Tank Reactors

4. Discussion and Conclusions

Funding

Conflicts of Interest

Notation

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI