Finding cones for K-cooperative systems

TL;DR

This paper introduces a numerical cone finding algorithm to analyze nonlinear systems via differential positivity, offering a new tool for studying multi-stability beyond traditional Lyapunov methods.

Contribution

The paper presents a novel cone finding algorithm that enhances the analysis of nonlinear systems through differential positivity, demonstrated on complex multi-stable and controlled oscillator examples.

Findings

Effective in identifying cones for nonlinear systems

Applicable to multi-stable and oscillatory dynamics

Provides a numerical alternative to Lyapunov analysis

Abstract

We design and test a cone finding algorithm to robustly address nonlinear system analysis through differential positivity. The approach provides a numerical tool to study multi-stable systems, beyond Lyapunov analysis. The theory is illustrated on two examples: a consensus problem with some repulsive interactions and second order agent dynamics, and a controlled duffing oscillator.

Figures

Equations

$$x(0)\in \mathcal{K}⟹x(t)\in \mathcal{K}\text{ for all }t\ge 0.$$

$$x_1(0)⪯x_2(0)⟹x_1(t)⪯x_2(t),$$

$$x_1{⪯}_{\mathcal{K}}{x}_2⟺x_2-x_1\in \mathcal{K}$$

$$h^T\delta x=0⟹h^T\partial f(x)\delta x>0$$

$$\dot{x}=f(x)\dot{\delta x}=\partial f(x)\delta x,$$

$$\delta x(0)\in \mathcal{K}\setminus \{0\}⟹\delta x(t)\in \text{mbox}int(\mathcal{K})\text{mbox}fort>0$$

$$\mathcal{K}_{\mathcal{H}}(H)=\{r:Hr\geq 0\}\$$

$$\mathcal{K}_{\mathcal{R}}(R)=\{r:r=Rp,p\geq 0\}\$$

$$(\alpha I+\partial f(x))R=RP$$

$$\partial f(x)\in conic-hull(\mathcal{A})\setminus \{0\},$$

$$conic-hull(\mathcal{A}):=\left\{A:A=\sum _{i=1}^k{p}_i{A}_i,p_i\ge 0\right\}$$

$$({\alpha }_{i}I+A_i)R=R{P}_i$$

$$\partial f(x)=\sum _{i=1}^k{p}_i^x{A}_i,$$

$$(\alpha I+\partial f(x))R=\left(\sum _{i=1}^k{\alpha }_i{p}_i^{x}I+\sum _{i=1}^k{p}_i^x{A}_i\right)R=$$

$$=\sum _{i=1}^k{p}_i^x({\alpha }_{i}I+A_i)R=\sum _{i=1}^k{p}_i^{x}R{P}_i=R\sum _{i=1}^k{p}_i^x{P}_i=RP$$

$$W_{A_i,{\tau }_i,w_i}:={\tilde{A}}_{A_i,{\tau }_i}-w_i{\tilde{r}}_i{{\tilde{h}}_i}^T,$$

$$\dot{x_i}=\sum _{j=1}^N{f}_{ij}(x_j-x_i)0<i\le N,$$

$$\mathcal{A}=\{A_{(1,1)},A_{(-1,1)},A_{(1,-1)},A_{(-0.9,-0.9)}\}$$

$${\dot{x}}_i=v_i,\tau {\dot{v}}_i=-v_i+\sum _{j=1}^N{f}_{ij}(x_j-x_i)0<i\le N,$$

$$\begin{array}{rl}\dot{x_p}& =x_v\\ \dot{x_v}& =-\alpha (x_p)-c{x}_v+k_f{x}_i\\ L\dot{x_i}& =k_p(x_{ref}-x_p)-k_e{x}_v-R{x}_i\end{array}$$

Full text

I INTRODUCTION

Despite the ubiquity of multi-stable systems (including, for example, bi-stable switches) in diverse areas of engineering, there is a lack of general tools for their analysis. Analyzing their behavior often involves finding ways of using Lyapunov theory confined to the basin of attraction of each equilibrium. This makes the process cumbersome and leads to local/regional stability results. A similar situation is encountered when trying to apply traditional Lyapunov stability theory to mono-stable scenarios where the equilibrium state depends on the parameters. Differential analysis [10] overcomes these limitations by using the linearized dynamics to characterize the behavior of a system independently of the location of the system attractor. In contraction theory, for example, convergence to a unique fixed point is established from the stability of the linearized dynamics computed along any system trajectory, using tools such as matrix measures or Lyapunov inequalities applied to the linearized vector field [17, 24, 26, 11].

We approach the problem from the perspective of differential positivity [12], which extends differential analysis approaches to multi-stable systems. This theory enables us to reduce the analysis of the system to the problem of finding a cone that is forward invariant under some linear operators. The existence of a suitable cone guarantees that almost every bounded trajectory converges to a fixed point. In this paper we address this problem algorithmically, proposing a numerical method to build a cone that is invariant (and contracting) for the linearized dynamics of the system.

On linear systems, differential positivity corresponds to classical positivity [18, 7, 22]. In the nonlinear setting, differential positivity shows tight connections with monotone systems [27, 15, 2], whose salient feature is that their trajectories preserve a partial order on the system state space. The connection builds on the fact that every cone induces a partial order on the system state space and that positivity of the linearization entails order preservation among systems trajectories. In this sense, monotonicity is an integral version of differential positivity.

This connection with monotone systems is particularly relevant because, traditionally, most of the literature takes monotonicity as an intrinsic property of a system. However, certifying monotonicity for a system is hard, with basic tests known for the case of orthant positivity (Metzler Jacobian, conditions based on graph connectivity) [1, 3, 5]. In this sense, our algorithm provides a way to establish monotonicity for systems that so far have remained elusive to the property.

At the most general level, differential positivity is defined with respect to a cone field on manifolds and can also capture limit cycle behavior. In this paper, however, we restrict ourselves to strict variants of differential positivity with constant cone fields. On vector spaces, finding a conic set that is forward invariant and contracting for the linearized dynamics has the practical interpretation that rays on the boundary of the cone are mapped into its interior. This form of projective contraction, connected to Perron-Frobenius theory [6], guarantees that the asymptotic behavior of the system is essentially one dimensional and, for systems on a vector space, leads to trajectories that almost globally converge to some fixed point.

To obtain conditions we can test more easily and to make an explicit distinction with the general theory of differential positivity, while maintaining the strong convergence results, we introduce the notion of strict K-cooperativity. The name is picked to highlight the connection with K-cooperative systems [15], a large subclass of monotone systems.

Finding a cone that is forward invariant under some linear operators is a difficult task computationally [21]. In this paper we derive basic necessary conditions for the existence of such a cone, then we leverage the classical theory of positive matrices [4] to algorithmically address its construction. The algorithm is illustrated and tested on a consensus dynamics problem with nonlinear and repulsive interactions, and on a bistable electro-mechanical system. We find cones for different parameter values and study their robustness to static uncertainties. The examples show the potential of our approach in nonlinear analysis.

**Notation: A matrix $P$ is non-negative, $P\geq 0$ (positive, $P>0$), if all its elements are non-negative (positive). The interior of a set $\mathcal{S}$ is denoted by $\operatorname{int}(\mathcal{S})$. A proper cone is a set $\mathcal{K}$ such that: (i) if $r_{1},\,r_{2}\in\mathcal{K}\text{ and }0\leq p_{1},\,p_{2}\in\mathbb{R}$, then $p_{1}r_{1}+p_{2}r_{2}\in\mathcal{K}$; (ii) $\operatorname{int}(\mathcal{K})\neq\{0\}$; and (iii) if $r\in\mathcal{K}$, then $-r\notin\mathcal{K}$. The dual cone of a cone $\mathcal{K}$ is denoted by $\mathcal{K}^{*}=\{h:{h}^{T}r\geq 0,\,\forall r\in\mathcal{K}\}$. Considering two cones $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$, $\mathcal{K}_{1}\subseteq\mathcal{K}_{2}$ denotes the usual set inclusion. We will use $\mathcal{K}_{1}\subset\mathcal{K}_{2}$ to denote $r\in\mathcal{K}_{1}\setminus\{0\}\implies r\in\operatorname{int}(\mathcal{K}_{2})$. **

II POSITIVITY AND MONOTONICITY

II-A Positivity

A linear dynamical system is positive if trajectories starting in a cone remain in the cone. Namely, a continuous time system $\dot{x}=Ax$, $x\in\mathbb{R}^{n}$, is positive with respect to the proper cone $\mathcal{K}$ if all its trajectories $x(\cdot)$ satisfy

[TABLE]

The classical example is given by systems whose state vector is constrained to be element-wise positive. The cone corresponds to the positive orthant $\mathcal{K}=\mathbb{R}^{n}_{+}$.

The trajectories of strictly positive systems ($x(0)\in\mathcal{K}\setminus\{0\}\implies x(t)\in\operatorname{int}(\mathcal{K})\text{ for all }t>0$) converge to a ray. This qualitative behavior makes positive systems central in engineering [18, 7, 22]. Convergence to a ray is a consequence of the Perron-Frobenius theorem which states that every strictly positive map $A$ admits a strictly dominant eigenvalue, that is, in the continuous time case, a right-most simple eigenvalue [18, Ch. 6]. The associated eigenvector is the dominant eigenvector. The Perron-Frobenius theorem can be seen as a consequence of the Banach contraction mapping theorem: for any positive $t>0$, the semiflow $e^{At}$ maps any ray in $\mathcal{K}$ into the interior of $\mathcal{K}$, which guarantees contraction of Hilbert’s metric [6]. Thus, $e^{At}$ has a fixed point, the dominant eigenvector, lying within $\mathcal{K}$ [4, 28].

There are two main approaches to certify strict positivity for closed linear systems: one requires finding a right-most isolated eigenvalue in $A$; the other searches for a cone that satisfies (1). The algorithm we present in Section IV-B approaches the latter numerically and, in contrast to other methods, allows for extensions to the nonlinear setting.

II-B Monotonicity and K-cooperativity

A dynamical system is monotone if its trajectories preserve some partial ordering $\preceq$ on the system state space. Namely, for any pair of trajectories $x_{1}(\cdot)$ and $x_{2}(\cdot)$, monotonicity requires:

[TABLE]

for any $t\geq 0$. Under mild conditions, almost all bounded trajectories of a monotone system converge to fixed points [15], a feature that makes them central in system theory and feedback control [27, 2, 19, 16].

Our interest in monotone systems stems from the fact that any proper cone $\mathcal{K}$ induces a partial order:

[TABLE]

Equation (3) suggests a new characterization for monotonicity, based on system linearization [12]:

Definition 1** (Strict K-cooperativity)**

Let $\mathcal{K}\subseteq\mathbb{R}^{n}$ be a proper cone. A system $\dot{x}=f(x)$ is strictly K-cooperative with respect to $\mathcal{K}$ if for all $x\in\mathcal{X},\,\delta x\in\mathcal{K}\setminus\{0\},\,h\in\mathcal{K}^{*}\setminus\{0\},$ the strict sub-tangentiality condition is met:

[TABLE]

This implies that, for a strictly K-cooperative system, any trajectory $(x(\cdot),\delta x(\cdot))$ of the prolonged system

[TABLE]

satisfies

[TABLE]

Proposition 1

Any strictly K-cooperative system is strictly differentially positive.

Proof:

From (5) and (6), [12, Definition 2] holds for all $T>0$. ∎

Strictly K-cooperative systems enjoy strong convergence properties, as shown by the theory of differentially positive systems [12, Corollary 5]:

Proposition 2** (Steady state of K-cooperative systems)**

Consider a strictly K-cooperative dynamical system $\dot{x}=f(x)$, $x\in\mathbb{R}^{n}$. Suppose that all trajectories are bounded. Then, for almost all $x\in\mathbb{R}^{n}$, the $\omega$-limit set $\omega(x)$ is a fixed point.

Proving the existence of a cone $\mathcal{K}$ to certify K-cooperativity is a difficult task. In this paper we propose an algorithm that, given a system, finds a cone $\mathcal{K}$ that satisfies (4). Indeed, by Proposition 2, the algorithm guarantees that the system’s trajectories almost globally converge to some fixed point.

III POLYHEDRAL CONES

III-A Representations, Membership, and Tests

We work with proper polyhedral cones adopting the following representations: {LaTeXdescription}

for any matrix $H$,

[TABLE]

for any matrix $R$,

[TABLE]

Strict inequalities characterize the interior of the cones. In the H-representation, the rows of $H$ define half-spaces through the origin, with the overall cone given by their intersection. In the R-representation, the columns of $R$ are generator rays and the overall cone is their conical hull. For a proper cone we need at least $n$ independent columns in $R$ and $n$ independent rows in $H$, where $n$ is the dimension of the space.

By the Minkowski-Weyl theorem [14], any polyhedral cone admits both representations. However, in dimensions above three, it can be computationally expensive to convert between the two representations. From this point onward, we focus on R-representation cones, noting that the analysis can be easily extended to H-representation cones using duality.

Proposition 3

Suppose $\mathcal{K}=\mathcal{K}_{\mathcal{R}}(R)$ for some matrix $R$. An $n$ dimensional system $\dot{x}=f(x),\,x\in\mathcal{X}$, is strictly K-cooperative with respect to $\mathcal{K}$ if for all $x\in\mathcal{X}$, there exists a matrix $P>0$ such that, for some $\alpha$,

[TABLE]

Proof:

Let $\delta x\in\mathcal{K}\setminus\{0\},\,h\in\mathcal{K}^{*}\setminus\{0\},\,{h}^{T}\delta x=0$, and $\mathcal{K}_{\mathcal{R}}(R)$ is proper. By definition, $\delta x=Rp$ for some element-wise positive vector $0\neq p\geq 0$. Note that at least one element of $p$ is strictly greater than zero.

Then, for any $\alpha\in\mathbb{R}$, ${h}^{T}\partial f(x)\delta x={h}^{T}(\alpha I+\partial f(x))\delta x={h}^{T}(\alpha I+\partial f(x))Rp={h}^{T}RPp>0$ and (4) holds. The strict inequality follows from the fact that $0\neq p\geq 0$, $0\neq h^{T}R\geq 0$, and $P>0$. ∎

Proposition 3 provides simple geometric conditions for K-cooperativity with respect to any polyhedral cone. These conditions need to hold for every $x\in\mathcal{X}$, which can be tested numerically through conical relaxations, as shown in Section III-B.

III-B Conical Relaxation

The tests in Proposition 3 can be made numerically tractable by finding a finite family of matrices $\mathcal{A}:=\{A_{1},A_{2},\dots,A_{k}\}$ such that, for all $x\in\mathcal{X}$:

[TABLE]

where, for finite $k$:

[TABLE]

A set $\mathcal{A}$ that satisfies (10) is referred to as a conical relaxation of the dynamics $\dot{x}=f(x)$.

Proposition 4

Suppose $\mathcal{K}=\mathcal{K}_{\mathcal{R}}(R)$ for some matrix $R$. An $n$ dimensional system $\dot{x}=f(x),\,x\in\mathcal{X}$, is strictly K-cooperative with respect to $\mathcal{K}$ if (10) holds and, for all $A_{i}\in\mathcal{A}$, there exists a matrix $P_{i}>0$ such that

[TABLE]

for some $\alpha_{i}\in\mathbb{R}$.

Proof:

Since (10) holds, $\partial f(x)$ can be written as:

[TABLE]

for some $p_{i}^{x}\geq 0$ and $\sum_{i=1}^{k}p_{i}^{x}\neq 0$.

Define $\alpha:=\sum_{i=1}^{k}\alpha_{i}p_{i}^{x}$ and $P:=\sum_{i=1}^{k}p_{i}^{x}P_{i}$, and note that $P>0$. Then:

[TABLE]

Thus, strict K-cooperativity follows from Proposition 3. ∎

Proposition 4 shows that K-cooperativity with respect to a given polyhedral cone can be determined by solving $k$ Linear Programming (LP) problems. However, some methods of building the family of matrices $\mathcal{A}$ lead to a combinatorial explosion in the number of matrices to test. General methods for producing tight conical relaxations are left as future work.

Proposition 4 also provide ways to test K-cooperativity for systems subject to uncertainties. Suppose that (11) holds for $\dot{x}=f(x)$. Then it also holds for any perturbed system $\dot{x}=f(x)+g(x)$ where $\partial g(x)\in\operatorname{conic-hull}(\mathcal{A})$. Moreover, in a situation where the uncertainties in the system can be represented in the linearizations by $\partial f(x)+\mathcal{Q}$, where $\mathcal{Q}$ is a given family of perturbations, the family of matrices $\mathcal{A}_{\mathcal{Q}}$ can be suitably adapted to the perturbations. If $\partial f(x)+Q\in\operatorname{conic-hull}(\mathcal{A}_{\mathcal{Q}})$, then (11) with $\mathcal{A}_{\mathcal{Q}}$ guarantees strict K-cooperativity of the perturbed system.

Given a conical relaxation, testing K-cooperativity with respect to a given cone is a tractable problem. However, deciding whether or not a system is K-cooperative with respect to some cone is a much harder question, even when a conical relaxation is used [21]. The work in Section IV is a first attempt in this direction.

IV FINDING CONES

IV-A Necessary Conditions

We present two necessary conditions for the existence of a cone that satisfies (11). The first is a spectral condition that can be verified for every $A_{i}\in\mathcal{A}$ individually, summarized in the following Proposition:

Proposition 5** (Spectral Condition)**

(11) can hold for a given $\mathcal{A}$ only if every $A_{i}\in\mathcal{A}$ has a strictly dominant eigenvalue.

Proof:

As in the proof of Proposition 3, (11) implies that every $A_{i}\in\mathcal{A}$ satisfies the strict sub-tangentiality condition with respect to $\mathcal{K}$.

This implies that a system $\dot{z}=A_{i}z$ is strictly positive with respect to $\mathcal{K}$, which in turn requires $A_{i}$ to have a strictly dominant eigenvalue (see [4, Theorems 4.3.37 and 4.3.41] for detailed proofs). ∎

The second is a geometric condition on the compatibility between the invariant spaces of each $A_{i}\in\mathcal{A}$ based on the following property:

Proposition 6

For any cone $\mathcal{K}$ that satisfies (11), every $A_{i}\in\mathcal{A}$ must have its dominant right eigenvector $\bar{r}_{i}$ in $\operatorname{int}(\mathcal{K})$ and its dominant left eigenvector $\bar{h}_{i}$ in $\operatorname{int}(\mathcal{K}^{*})$.

Proof:

Continuing from the proof of Proposition 5, if a system $\dot{z}=A_{i}z$ is strictly positive with respect to $\mathcal{K}$, a dominant right eigenvector of $A_{i}$ must be in $\operatorname{int}(\mathcal{K})$ [4, Theorem 4.3.34]. Additionally, the dual system $\dot{\eta}=A_{i}^{T}\eta$ must be strictly positive with respect to $\mathcal{K}^{*}$ [4, Theorem 4.3.45]. As such, a dominant right eigenvector of $A_{i}^{T}$ must be in $\operatorname{int}(\mathcal{K}^{*})$. This is a dominant left eigenvector of $A_{i}$. ∎

We observe that the property outlined in Proposition 6 refers to a common cone $\mathcal{K}$, therefore, it can be used to characterize a necessary condition based on the construction of two useful cones $\mathcal{K}_{inner}(\mathcal{A})$ and $\mathcal{K}_{outer}(\mathcal{A})$. These cones will be used to initialize and terminate our cone finding algorithm (Section IV-B). The construction of $\mathcal{K}_{inner}(\mathcal{A})$ and $\mathcal{K}_{outer}(\mathcal{A})$ is detailed below, building on the ‘Orientation Trick’ of [9].

Proposition 7** (Geometric Condition)**

(11) can hold for a given $\mathcal{A}$ only if $\mathcal{K}_{inner}(\mathcal{A})\subset\mathcal{K}_{outer}(\mathcal{A})$

Proof:

Assume (11). This implies that Proposition 6 holds for a given family of matrices $\mathcal{A}$ and cone $\mathcal{K}$. Consider the left and right dominant eigenvectors of Proposition 6. Then ${\bar{h}_{j}}^{T}\bar{r}_{k}>0$ for all $(j,k)$.

If $\tilde{h}_{1}=\bar{h}_{1}$, by construction $\tilde{r}_{k}=\bar{r}_{k}$ and $\tilde{h}_{j}=\bar{h}_{j}$. If $\tilde{h}_{1}=-\bar{h}_{1}$, by construction $\tilde{r}_{k}=-\bar{r}_{k}$ and $\tilde{h}_{j}=-\bar{h}_{j}$. In both cases ${\bar{h}_{j}}^{T}\bar{r}_{k}>0\iff{\tilde{h}_{j}}^{T}\tilde{r}_{k}>0$. But if $\mathcal{K}_{inner}(\mathcal{A})\not\subset\mathcal{K}_{outer}(\mathcal{A})$ there exists some $(j,k)$ such that ${\tilde{h}_{j}}^{T}\tilde{r}_{k}\leq 0$; which leads to a contradiction. ∎

Note that checking whether $\mathcal{K}_{inner}(\mathcal{A})\subset\mathcal{K}_{outer}(\mathcal{A})$ amounts to computing $\tilde{H}$ and $\tilde{R}$ and verifying that $\tilde{H}\tilde{R}>0$. As such, it is a very scalable test.

Remark 1

Because of the complexity of producing tight conical relaxations and finding cones, these two conditions can be used to quickly rule out K-cooperativity by testing sets of randomly sampled $\partial f(x)$ of a candidate system.

IV-B Cone Finding Algorithm

A simple procedure to find a cone that satisfies (11) for a given $\mathcal{A}$ is described below. We begin by illustrating the main ideas on the easier case of a single matrix $A$. Given an initial cone $\mathcal{K}_{\mathcal{R}}(R^{(0)})$ and a linear operator $W$, consider the matrix $R^{(1)}$ formed by combining the columns of $R^{(0)}$ and $WR^{(0)}$, $R^{(1)}=[R^{(0)},WR^{(0)}]$. This can be repeated recursively such that $R^{(k+1)}=[R^{(k)},WR^{(k)}]$. If $W:=(\alpha I+A)$, $R^{(k+1)}$ will trivially satisfy $(\alpha I+A)R^{(k)}=R^{(k+1)}P$ with $P\geq 0$. Moreover, the projective contraction property of positive systems tells us that if $\{A\}$ satisfies Proposition 5 and $\mathcal{K}_{inner}(\{A\})\subset\mathcal{K}_{\mathcal{R}}(R^{(0)})\subset\mathcal{K}_{outer}(\{A\})$, there will be an $\alpha$ for which the sequence of cones $R^{(k)}$ converges to a cone $\mathcal{K}_{\mathcal{R}}(R^{(k)})\to\mathcal{K}(R)$ in a finite number of steps. It follows that $\mathcal{K}_{\mathcal{R}}(R)$ will satisfy (11) for $\{A\}$ and $P\geq 0$.

This motivates an algorithm for finding cones for a family of systems. We will first define a time-step based parametrization. Then, the technical requirement of satisfying (11) with $P_{i}>0$ is tackled by adding a widening operation to $W_{i}$. A procedure for initializing $R^{(0)}$ is also needed. Finally, unlike the single matrix case, no guarantees about convergence can be made and additional termination conditions need to be provided.

Time-step based parametrization: We denote the dominant eigenvalue of a matrix $\{A_{i}\}$ satisfying Proposition 5 as $\lambda_{i}$. The values of $\alpha_{i}$ for which the single matrix algorithm described above is guaranteed to converge are the values for which $(\alpha_{i}I+A_{i})$ has a real simple positive eigenvalue with a larger absolute value than any other eigenvalue. Such an eigenvalue will be referred to as an absolutely dominant eigenvalue of $(\alpha_{i}I+A_{i})$. We adopt the parametrization $\tilde{A}_{A_{i},\tau_{i}}:=I+\tau_{i}(A_{i}-\lambda_{i}I)$ which corresponds to taking $\alpha_{i}=1/\tau_{i}-\lambda_{i}$. The only free parameter in this formulation is $\tau_{i}$, which acts like a discretization time-step. For any $A_{i}$ satisfying Proposition 5, there is a non-empty open interval $\tau_{i}\in(0,T_{i})$ for which $\tilde{A}_{A_{i},\tau_{i}}$has a dominant eigenvalue.

Widening operation: Setting $W_{A_{i},\tau_{i}}:=\tilde{A}_{A_{i},\tau_{i}}$ and using $W_{A_{i},\tau_{i}}$ to product the sequence of matrices $R^{(k)}$, we have that $(\alpha_{i}I+A_{i})R^{(k)}=R^{(k+1)}P_{i}$ with $P_{i}\geq 0$. However, for the strong convergence properties of strict K-cooperativity we require $P_{i}>0$. With this aim, define $W_{A_{i},\tau_{i},w_{i}}$as:

[TABLE]

where $\tau_{i},w_{i}>0$ and fixed, and $\tilde{r}_{i}$ and $\tilde{h}_{i}$ are the right and left dominant eigenvectors of $A_{i}$ in $\mathcal{K}_{inner}(\mathcal{A})$ and $\mathcal{K}_{outer}(\mathcal{A})$ respectively. Now, if a cone $\mathcal{K}_{\mathcal{R}}(R^{(k+1)})$ is formed as $R^{(k+1)}=[R^{(k)},W_{A_{i},\tau_{i},w_{i}}R^{(k)}]$, then $W_{A_{i},\tau_{i},w_{i}}R^{(k)}=R^{(k+1)}\bar{P}_{i}$ for some $\bar{P}_{i}\geq 0$. But, $\mathcal{K}_{inner}(\mathcal{A})\subset\mathcal{K}_{\mathcal{R}}(R^{(k)})\subset\mathcal{K}_{outer}(\mathcal{A})$, so ${\tilde{h}_{i}}^{T}R^{(k)}>0$, and $\tilde{r}_{i}=R^{(k+1)}\tilde{P}_{i}$ for some $\tilde{P}_{i}>0$. Hence, we get $\tilde{A}_{A_{i},\tau_{i}}R^{(k)}=w_{i}\tilde{r}_{i}{\tilde{h}_{i}}^{T}+R^{(k+1)}\bar{P}_{i}=R^{(k+1)}P_{i}$ for some $P_{i}>0$. The additional widening provided by $-w_{i}\tilde{r}_{i}{\tilde{h}_{i}}^{T}$ can be interpreted as moving new rays slightly away from $\tilde{r}_{i}$, and $w_{i}$ can be thought of as a ‘widening coefficient’. Given $\tau_{i}\in(0,T_{i})$, there will be a maximum $w_{i}$ for which $W_{A_{i},\tau_{i},w_{i}}$ still has an absolutely dominant eigenvalue. In what follows, the set of $W_{A_{i},\tau_{i},w_{i}}$ operators is denoted by $\mathcal{W}_{\mathcal{A},\tau,w}$.

Now if a cone $\mathcal{K}_{\mathcal{R}}(R^{(k+1)})$ is formed with $R^{(k+1)}=[R^{(k)},W_{1}R^{(k)},\dots,W_{l}R^{(k)}]$ with $W_{i}\in\mathcal{W}_{\mathcal{A},\tau,w}$, there must then exist $\alpha_{i}$ such that $(\alpha_{i}I+A_{i})R^{(k)}=R^{(k+1)}P_{i}$ with $P_{i}>0$ for all $A_{i}\in\mathcal{A}$. As such, if $R^{(k)}$ converges to some matrix $R$ as $k\to\infty$, then $\mathcal{K}_{\mathcal{R}}(R)$ satisfies (11).

Initializing $R^{(0)}$: An initial cone $\mathcal{K}_{\mathcal{R}}(R^{(0)})$ that satisfies $\mathcal{K}_{inner}(\mathcal{A})\subset\mathcal{K}_{\mathcal{R}}(R^{(0)})\subset\mathcal{K}_{outer}(\mathcal{A})$ can be generated by starting with $\mathcal{K}_{inner}(\mathcal{A})$ and adding random vectors of small magnitude until $\mathcal{K}_{inner}(\mathcal{A})\subset\mathcal{K}_{\mathcal{R}}(R^{(0)})$, restarting with a smaller noise magnitude if at any point $\mathcal{K}_{\mathcal{R}}(R^{(0)})\not\subset\mathcal{K}_{outer}(\mathcal{A})$. We call this operation Initialize$(\mathcal{K}_{inner},\mathcal{K}_{outer})$.

Termination conditions: The algorithm terminates if (11) is satisfied, which implies strict K-cooperativity. This is tested at every iteration. If after an iteration of the algorithm the cone is no longer in the interior of $\mathcal{K}_{outer}$, or if the maximum allowed number of iterations has been exceeded, the algorithm terminates without finding a cone. Algorithm 2 below summarizes the overall process:

Remark 2

In the implementation of the algorithm we also remove redundant rays in $\mathcal{K}_{\mathcal{R}}(R^{(k)})$ and test (11) every $N$ iterations for improved speed.

V EXAMPLES

V-A Consensus

We consider K-cooperativity of standard consensus dynamics [20, 23] represented by:

[TABLE]

where each agent $x_{i}$ is modeled by a simple integrator driven by the weighted differences with its neighboring agents, characterized by functions $f_{ij}:\mathbb{R}\to\mathbb{R}$ such that $f_{ij}(0)=0$. The system has a continuum of equilibria $x_{e}=\alpha 1_{N}$ where $\alpha\in\mathbb{R}$ and $1_{N}$ is the vector of $N$ ones. The agents reach consensus when $x_{1}=\dots=x_{N}$.

For linear positive weights convergence to consensus is guaranteed by Perron-Frobenius theory [18, 25]. This result extends to nonlinear strictly increasing weights: their slope is strictly positive, the system is cooperative, and bounded trajectories converge to the consensus equilibria [15]. Convergence is also guaranteed for the case of unconstrained linear weights, whenever $1_{N}$ is a dominant eigenvector of the system.

The presence of uncertainties or nonlinear weights that are not strictly increasing makes the problem more challenging. However, from Proposition 2, we can use K-cooperativity to study consensus problems: the trajectories of strictly K-cooperative consensus dynamics will converge to consensus whenever $1_{N}$ belongs to the interior of $\mathcal{K}$ [8, Section V]. This demonstrates how Algorithm 2 provides a numerical tool to study nonlinear and robust consensus problems away from positive weights.

For illustration, we apply Algorithm 2 to the network of agents in Fig. 1, where the black edges represent linear weights normalized to $1$, and blue and red edges represents the nonlinear weights $f_{15}$ and $f_{42}$, respectively. Using $A_{a,b}$ to denote the Jacobian of the consensus dynamics in Fig. 1, for $f_{15}^{\prime}=a$ and $f_{42}^{\prime}=b$, we used Algorithm 2 to find a common cone $\mathcal{K}$ for the polytope of linearizations:

[TABLE]

This polytope is visualized in Fig. 1 and includes, for example: $-1\leq f_{15}^{\prime}\leq 1$ and $f_{42}^{\prime}=1$, or $-0.9\leq f_{15}^{\prime}\leq 1$ and $-0.9\leq f_{42}^{\prime}\leq 1$.

Consensus is thus reached when uncertainties of magnitude less than one affect the weight between nodes $1$ and $5$ or between nodes $4$ and $2$. Consensus is also reached for any nonlinear weight between these nodes whose slope is bounded between $-1$ and $1$. Indeed, $f_{15}(x_{5}-x_{1})=\sin(x_{5}-x_{1}$) and $f_{42}(x_{2}-x_{4})=x_{2}-x_{4}$ is compatible with consensus. Finally, consensus is reached when both $f_{15}$ and $f_{42}$ are nonlinear/uncertain provided that they are constrained to a smaller negative range.

We complete the example by considering agent dynamics extended to second order networks of the form:

[TABLE]

where $\tau$ is a homogeneous time constant and $f_{ij}$ captures the coupling between agents. We adopt the topology in Fig. 1. For the second order case, consensus equilibria are given by the subspace $(x,v)_{e}=(1_{N},0)$.

The necessary conditions in Propositions 5 and 7 are not satisfied for $\tau\geq 1$. A cone $\mathcal{K}$ with $(1_{N},0)\in\mathcal{K}$ was successfully found for $\tau=0.3$. Thus, the trajectories of (14) with $\tau=0.3$ asymptotically converge to consensus for any perturbed/nonlinear $f_{15}$ and $f_{42}$ constrained to the intervals outlined above. An example trajectory with $\tau=0.3$, $f_{15}=0.9\sin(x_{5}-x_{1})$, and $f_{42}=-0.9\sin(x_{2}-x_{4})$ is shown in Fig. 2.

V-B Controlled Duffing Oscillator

We study the three dimensional system given by:

[TABLE]

This is based on the Duffing oscillator example in [13]. The first two states capture the mechanics of the planar mechanical system and the DC motor inertia, while the third state is the electrical equation of the DC motor. We use $\alpha(\cdot)$ to characterize a nonlinear spring. Although simple, this system captures a rich range of qualitative behaviors including mono-stability, bi-stability, and limit cycles. We consider parameters $c=5,k_{f}=1,L=0.1,k_{e}=1,R=1$.

Algorithm 2 shows that (15) is strictly K-cooperative with respect to the cone shown in Fig. 3 for $-2\leq\alpha^{\prime}\leq 5$ and $0\leq k_{p}^{\prime}\leq 3$. This includes spring characteristics that can lead to bi-stability, as shown in Fig. 3 (right).

Following the discussion on robustness at the end of Section III-B, the cone found can also be leveraged to explore K-cooperativity for different sets of parameters. For example, taking $c=8$, $0\leq k_{p}^{\prime}\leq 3$, and all other parameters as above, one can immediately show (using LP) that system (15) is strictly K-cooperative with respect to the same $\mathcal{K}_{\mathcal{R}}(R)$ for the wider range $-2.6<\alpha^{\prime}<6.1$.

VI CONCLUSIONS

K-cooperativity combined with a cone finding algorithm make the analysis of multi-stable nonlinear systems accessible. The approach was illustrated on two examples, showing the potential of the theory in application. The performance of Algorithm 2 strongly depends on the analyzed dynamics. For example, the cone found for (13) had 52 rays but the one for (14) had 11738 when $\tau=0.3$ (10 hours of computation), and a cone with 4278 rays was found when $\tau=0.1$ (20 minutes of computation). Generally, matrices with a small gap between their dominant eigenvalue and their complex eigenvalues (in the conical relaxation (10)) lead to cones with a large number of rays. Future work will focus on scaling up the approach to systems of large dimension. Following [9], we will also study the completeness of the algorithm and we will characterize conditions for its convergence.