Determining Asymptotic Stability and Robustness of Networked Systems

: 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 efﬁcient computational method that relies on semideﬁnite 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 ﬁelds as well.


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 Sections 2.1 and 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.

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 ẋi where f i : R n → R n describes the system's dynamics for all i, i ∈ {1, 2, . . ., N}, x i (t) ∈ R n , and D ∈ R n×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) = 0 if and only if states x i (k) and x j (l) are coupled, for all i, j and i = j.Positive constant k corresponds to the coupling strength.Next, we provide results on asymptotic stability of the coupled system based on Lyapunov theory.

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 ẋ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, where The second equality in the first line of Equation ( 3) follows from the symmetry of M. Thus, if holds for all i unless x 1 = 0 and X i = 0 for all i, and, thus, unless x = 0, then 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 ∈ Q ⊂ R r .Then, we can rewrite condition (4) by letting g(q 1 , x 1 ) ≡ f 1 (x 1 ) and g(q i , x i ) ≡ f i (x i ).In this case, if condition (4) holds, for all x 1 , x i ∈ R n and all q 1 , q i ∈ Q, then asymptotic stability of the coupled system is guaranteed.

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 Ṽ(x 1 , X i ) < c, c ∈ R, c > 0, and all i then Ṽ(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.

Coupled Linear Systems
Consider N coupled linear dynamical systems of the form ẋ1 = Ax and ẋi = Bx for all i > 1, where A, B ∈ R n×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 V(X i ) = X T i PX i , where P = P T 0, inequality (4) becomes 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 given A, B, D, N, k search for M, P subject to M = M T 0, P = P T 0, (5).
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.

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 2d, x ∈ 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]: where f i are polynomial functions.Now, F(x) is a SOS if and only if there exists a positive semidefinite matrix R and The entries of vector χ consist of all monomial combinations of the elements of vector x up to degree d (including x 0 (i) = 1) and, thus, its length is = ( 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].
Thus, j = 1, . . ., 5, and, for example, c 1 = 2, c 3 = −1, In general, in order to find R, we solve the optimisation problem associated with the following semidefinite program: 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 Ṽ(x 1 , X i ), given by Equation ( 2), such that (4) holds can be cast as follows Here, some additional remarks: • Consider a rational function F(x), F(x) = f (x) g(x) , where f (x) and g(x) are polynomial functions.Then, F(x) ≥ 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 , 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 = ( n+d d ).

Different Coupling Configurations
The presented results can be extended to any symmetric linear coupling topology with corresponding Laplacian matrix −C by letting V(X i ) = X T i PX 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 − 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 ∈ R N×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 ∈ R N , where λ min (−C) denotes the smallest positive eigenvalue of the matrix −C.
By the Spectral Theorem, there exists a unitary matrix S such that S −1 = S T and −SCS T = Λ −C , where Λ −C is the diagonal matrix containing the eigenvalues of −C and, without loss of generality, its first entry is 0.Moreover, implies that SUS 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 ∈ R N and Λ C being the diagonal matrix containing the eigenvalues of C, which completes the proof.
Consider, again, difference X i ≡ x i − x 1 , where X i ∈ R n and X 1 = 0 [25].Let matrix F ∈ R N×N be defined as follows, Then, Moreover, note that for symmetric coupling determined by Laplacian matrix −C, C = CF and where the inequality follows from Lemma 1, αN = λ min (−(FC) gL ), and (FC) gL is the matrix obtained by deleting the first row and the first column of FC.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(AB) = eig(BA), has been proven in [28] and, thus, eig(FC) = eig(CF) = eig(C).
Since the first row of matrix FC is zero, it follows that λ min (−(FC) gL ) = λ min (−C).
The different coupling configurations are illustrated in Figure 1.

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 → ∅ s.The reaction describing the growth of product by consuming substrate is given by xs k → 2x, where reaction constant k is described below.

𝑥𝑠 → 2𝑥
,  The corresponding mathematical model based on ordinary differential equations is given by 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 kxs, where k = 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 x = x + x * and s = s + s * , first, solve the following problem given Equation (16), u = 2, β ∈ R search for M subject to p 1 ( x, s), p 2 ( x, s), p 3 ( x, s) are SOS functions of degree 4, −0.001 for different values of β, and, then, solve 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.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.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), given Equation (16), subject to p j ( x1 , s1 , xi , si ) are SOS functions of degree 4, j ∈ {1, . . ., 5}, By, first, solving a problem analogous to problem (18), with 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 ≤ u ≤ 2.153, which translates to 4.618 ≤ x * ≤ 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 where p 6 ( x1 , s1 , xi , si ) and p 7 ( x1 , s1 , xi , si ) are SOS functions of degree 4. We, again, observe that the regions that bound solutions are given by the black solid lines.

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.

R, R m×n
real numbers, m × n real matrices A (i,j) (i, j)th entry of matrix A ∈ R m×n I n the identity matrix, I n ∈ R n×n A T transpose of matrix A ∈ R m×n ẋ 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 ⊗ B, A ∈ R m×n , B ∈ R p×q The Kronecker product: 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 π N .(d) Ring of 2k-nearest neighbour coupled systems [29]:

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

Figure 3 .
Figure3.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.