High-order maximum principles for the stability analysis of positive bilinear control systems

# High-order maximum principles for the stability analysis of positive bilinear control systems

Gal Hochma and Michael Margaliot Research supported in part by the Israel Science Foundation (ISF). GH is with the School of Electrical Engineering-Systems, Tel Aviv University, Israel 69978. Email: gal.hoc@gmail.comMM (corresponding author) is with the School of Electrical Engineering-Systems and with the Sagol School of Neuroscience, Tel Aviv University, Israel 69978. Email: michaelm@eng.tau.ac.ilAn abridged version of this paper was presented at the nd IEEE Conference on Decision and Control [1].
###### Abstract

We consider a continuous-time positive bilinear control system (PBCS), i.e. a bilinear control system with Metzler matrices. The positive orthant is an invariant set of such a system, and the corresponding transition matrix  is entrywise nonnegative for all time . Motivated by the stability analysis of positive linear switched systems (PLSSs) under arbitrary switching laws, we fix a final time  and define a control as optimal if it maximizes the spectral radius of . A recent paper [2] developed a first-order necessary condition for optimality in the form of a maximum principle (MP). In this paper, we derive higher-order necessary conditions for optimality for both singular and bang-bang controls. Our approach is based on combining results on the second-order derivative of the spectral radius of a nonnegative matrix with the generalized Legendre-Clebsch condition and the Agrachev-Gamkrelidze second-order optimality condition.

Positive switched systems, stability under arbitrary switching laws, variational approach, high-order maximum principles, Perron-Frobenius theory.

## I Introduction

Consider the continuous-time linear switched system

 ˙x(t) =Aσ(t)x(t), x(0) =x0, (1)

where  is the state vector, and  is a piecewise constant function referred to as the switching signal. This models a system that can switch between the two linear subsystems

 ˙x=A0x and ˙x=A1x.

Recall that (I) is said to be globally uniformly asymptotically stable (GUAS) if there exists a class  function111A continuous function belongs to the class if it is strictly increasing and . A continuous function belongs to the class if for each fixed , belongs to , and for each fixed , the mapping is decreasing and as .  such that for any initial condition  and any switching law , the corresponding solution of (I) satisfies

 |x(t)|≤β(|x0|,t), for all t≥0.

This implies in particular that

 limt→∞x(t)=0, for all σ and all x0∈Rn, (2)

and or linear switched systems, (2) is in fact equivalent to GUAS (see, e.g., [3]). Switched systems and, in particular, their stability analysis are attracting considerable interest in the last two decades; see e.g. the survey papers [4, 5, 6, 7] and the monographs [8, 9, 10, 11, 12, 13].

It is well-known that a necessary (but not sufficient) condition for GUAS of (I) is the following.

###### Assumption 1

The matrix  is Hurwitz for all .

Recall that a linear system

 ˙x=Ax, (3)

with , is called positive if the positive orthant

 Rn+:={x∈Rn∣xi≥0,i=1,…,n}

is an invariant set of the dynamics, i.e.,  implies that  for all .

Positive systems play an important role in systems and control theory because in many physical systems the state-variables represent quantities that can never attain negative values (e.g. population sizes, probabilities, concentrations, buffer loads) [14, 15, 16]. A necessary and sufficient condition for (3) to be positive is that is a Metzler matrix, that is, for all . If  is Metzler then  is (entrywise) nonnegative for all . By the Perron–Frobenius theory, the spectral radius of  (i.e., the eigenvalue with maximal absolute value) is real and nonnegative, and since  is non-singular, it is in fact positive.

If both  and  are Metzler and  then (I) is called a positive linear switched system (PLSS). Mason and Shorten [17], and independently David Angeli, posed the following.

###### Conjecture 1

If (I) is a PLSS, then Assumption 1 provides a sufficient condition for GUAS.

Had this conjecture been true, it would have implied that determining GUAS for a PLSS is relatively simple. (See [18] for analysis of the computational complexity of determining whether any matrix in a convex set of matrices is Hurwitz.) Gurvits, Shorten, and Mason [19] proved that Conjecture 1 is in general false (see also [20]), but that it does hold when  (even when the number of subsystems is arbitrary). Their proof in the planar case is based on showing that the PLSS admits a common quadratic Lyapunov function (CQLF). (For more on the analysis of switched systems using CQLFs, see [5, 4, 21, 22, 23].) Margaliot and Branicky [24] derived a reachability–with–nice–controls–type result for planar bilinear control systems, and showed that the proof of Conjecture 1 when  follows as a special case. Fainshil, Margaliot, and Chigansky [25] showed that Conjecture 1 is false already for the case . In general, it seems that as far as the GUAS problem is concerned, analyzing PLSSs is not simpler than analyzing linear switched systems.

There is a rich literature on sufficient conditions for GUAS, see, e.g., [5, 6, 8, 7, 12]. A more challenging problem is to determine a necessary and sufficient condition for GUAS. What makes this problem difficult is that the set of all possible switching laws is huge, so exhaustively checking the solution for each switching law is impossible.

A natural idea is to try and characterize a “most destabilizing” switching law  of the switched system, and then analyze the behavior of the corresponding trajectory . If  converges to the origin, then so does any trajectory of the switched system and this establishes GUAS. This idea was pioneered by E. S. Pyatntisky [26, 27], who studied the celebrated absolute stability problem (ASP). This variational approach was further developed by several scholars including N. E. Barabanov and L. B. Rapoport, and proved to be highly successful; see the survey papers [28, 29, 30], the related work in [31, 32], and the recent extensions to the stability analysis of discrete–time linear switched systems in [33, 34].

A first attempt to extend the variational approach to the stability analysis of PLSSs was taken in [35] using the classical Pontryagin maximum principle (PMP). Recently, Fainshil and Margaliot [2] developed an alternative approach that combines the Perron-Frobenius theory of nonnegative matrices with the standard needle variation used in the PMP.

The goal of this paper is to derive stronger, higher-order necessary conditions for optimality. We thus begin by reviewing the first-order MP in [2].

### I-a Stability analysis of PLSSs: a Variational Approach

The variational approach to the stability analysis of a linear switched system includes several steps. The first step is relaxing (I) to the bilinear control system (BCS)

 ˙x =(A+uB)x,u∈U, (4) x(0) =x0,

where  , , and  is the set of measurable controls taking values in . Note that for  [], Eq. (4) yields  [, i.e., trajectories of the BCS corresponding to piecewise constant bang-bang controls are also trajectories of the original switched system.

The BCS (4) is said to be globally asymptotically stable (GAS) if  for all  and all . Since every trajectory of the switched system (I) is also a trajectory of (4), GAS of (4) implies GUAS of the linear switched system. It is not difficult to show that the converse implication also holds, so the BCS is GAS if and only if the linear switched system is GUAS. Thus, the GUAS problem for the switched linear system (I) is equivalent to the GAS problem for the BCS (4).

From here on we assume that the switched system is positive, i.e.  is Metzler for all . For the BCS, this implies that if , then  for all  and all . Thus (4) becomes a positive bilinear control system (PBCS).

For , and , let  denote the solution at time  of the matrix differential equation

 ddtC(t,a,u) =(A+Bu(t))C(t,a,u), C(a,a,u) =I. (5)

It is straightforward to verify that the solution of (4) satisfies  for all  and all . In other words,  is the transition matrix from time  to time  of (4) corresponding to the control . To simplify the notation, we will sometimes omit the dependence on  and just write .

When the initial time is  we write (I-A) as

 ˙C(t) =(A+Bu(t))C(t), C(0) =I. (6)

For a PBCS,  is a non-negative matrix for all  and all . Since it is also non-singular, the spectral radius  is a real and positive eigenvalue of , called the Perron root. If this eigenvalue is simple then the corresponding eigenvector , called the Perron eigenvector, is unique (up to multiplication by a scalar). The next step in the variational approach is to relate  to GAS of the PBCS.

Define the generalized spectral radius of the PBCS (4) by

 ρ(A,B):=limsupt→∞ρt(A,B),

where

 ρt(A,B):=maxu∈U(ρ(C(t,u)))1/t. (7)

Note that the maximum here is well-defined, as the reachable set of (I-A) corresponding to  is compact [36]. In fact, this is why we consider a bilinear control system with controls in  rather than the original linear switched system with piecewise constant switching laws.

The next result relates the GAS of the PBCS to .

###### Theorem 1

The PBCS (4) is GAS if and only if

 ρ(A,B)<1.

Thm. 1 already appeared in [2], but without a proof. For the sake of completeness we include its proof in the Appendix.

###### Remark 1

It follows from (7) and Thm. 1 that if  for some  and , then the PBCS is not GAS. Indeed, for any integer , define  via the periodic extension of , and let  denote the corresponding solution of (I-A) at time . Then

 ρ(¯C(kT))=(ρ(C(T)))k,

so (7) yields

 ρkT(A,B) ≥(ρ(¯C(kT)))1/(kT) =(ρ(C(T)))1/T ≥1,

and this implies that

Thm. 1 motivates the following optimal control problem.

###### Problem 1

Consider the PBCS (I-A). Fix an arbitrary final time . Find a control  that maximizes .

The main result in [2] is a first-order necessary condition for optimality. Let  denote the transpose of the matrix .

###### Theorem 2

[2] Consider the PBCS (I-A). Suppose that  is an optimal control for Problem 1. Let  denote the corresponding solution of (I-A) at time , and let . Suppose that  is a simple eigenvalue of . Let  [] be an eigenvector of  [] corresponding to , normalized such that . Let  be the solution of

 ˙q =−(A+Bu∗)′q, (8) q(T) =w∗,

and let  be the solution of

 ˙p =(A+Bu∗)p, (9) p(0) =v∗.

Define the switching function  by

 m(t):=q′(t)Bp(t). (10)

Then for almost all ,

 u∗(t)={1,m(t)>0,−1,m(t)<0. (11)

This MP has some special properties.

###### Remark 2

First, note that (8) implies that

 q′(0)=q′(t)C∗(t), for all t∈[0,T].

In particular, substituting  yields

 q′(0) =q′(T)C∗(T) =(w∗)′C∗(T) =ρ∗(w∗)′,

as  is an eigenvector of  corresponding to the eigenvalue . Since scaling  by a positive constant has no effect on the sign of , this means that the final condition  in (8) can be replaced by the initial condition . This leads to an MP in the form of a one-point boundary value problem (with the unknown  as the initial conditions at time ).

###### Remark 3

Note that

 m(T) =q′(T)Bp(T) =(w∗)′BC∗(T)p(0) =(w∗)′Bρ∗v∗ =q′(0)Bp(0) =m(0). (12)

Thus, the switching function is “periodic” in the sense that

One difficulty in applying Theorem 2 is that both  and  are unknown. There are cases where this difficulty may be alleviated somewhat berceuse  can be expressed in terms of . The next example demonstrates this.

###### Example 1

Consider an optimal bang-bang control in the form

 u∗(t)={1,t∈(0,t1),−1,t∈(t1,T),

where . The corresponding transition matrix is

 C∗(T)=exp((A−B)τ2)exp((A+B)τ1),

where  and . Thus,  and  satisfy

 exp((A−B)τ2)exp((A+B)τ1)v∗ =ρ∗v∗,

and

 exp((A+B)′τ1)exp((A−B)′τ2)w∗ =ρ∗w∗. (13)

Suppose that  and  are symmetric matrices. Then (13) becomes

 exp((A+B)τ1)exp((A−B)τ2)w∗=ρ∗w∗,

and multiplying this on the left by  yields

 C∗(T)exp((A−B)τ2)w∗=ρ∗exp((A−B)τ2)w∗.

Since the Perron eigenvector of  is unique (up to multiplication by a constant) this means that

 exp((A−B)τ2)w∗=rv∗,

for some

The MP in Theorem 2 is a necessary, but not sufficient, condition for optimality and it is possible of course that a control satisfying this MP is not an optimal control. The next example demonstrates this.

###### Example 2

Consider a PBCS satisfying the following properties:

• The matrix is symmetric. Its maximal eigenvalue  is simple with corresponding eigenvector , and

 z′Bz=0. (14)
• The matrices  and  are Metzler;

• .

(A specific example is , and . Indeed, here , and it is straightforward to verify that all the properties above hold.)

Consider the possibility that the singular control  is optimal. Then

 ρ(u):=ρ(exp(AT))=exp(μT).

Since  is symmetric, the corresponding right and left eigenvector is , so in the MP . Thus, (9) and (8) yield

 p(t) =exp(At)z =exp(μt)z,

and

 q(t) =exp(A′(T−t))z =exp(μ(T−t))z.

Substituting this in (10) yields

 m(t) =exp(μT)z′Bz≡0.

Thus,  (vacuously) satisfies Thm. 2. However, since  the control  yields

 ρ(~u):=exp((A+B)T)>ρ(u),

so clearly  is not an optimal control.

The reason that  in Example 2 cannot be ruled out is that Thm. 2 is a first-order MP. More specifically, its derivation is based the following idea. Suppose that  is a candidate for an optimal control. Introduce a new control  by adding a needle variation to , i.e.

 ~u(t):={a,t∈[τ,τ+ϵ),u(t),otherwise,

where  is a Lebesgue point of , and  is sufficiently small, and analyze the difference  to first-order in . For ,

 C(T,~u)=exp(A(T−τ−ϵ))exp((A+aB)ϵ)exp(Aτ),

so

 ddϵC(T,~u)∣∣∣ϵ=0 =aexp(A(T−τ))Bexp(Aτ).

Combining this with known results on the derivative of a simple eigenvalue of a matrix (see, e.g. [37, Chapter 6]) yields

 ρ(C (T,~u))=ρ(C(T,u))+ϵaw′exp(A(T−τ))Bexp(Aτ)v+o(ϵ). (15)

If  then  for all sufficiently small  and thus  is not optimal. However, in Example 2 the term multiplying  in (15) is zero for all , , and , so a first-order analysis cannot rule out the possibility that  is optimal.

Summarizing, Example 2 suggests that there is a need for a higher-order MP, i.e., an MP that takes into account higher-order terms in the Taylor expansion of  with respect to , and can thus be used to rule out the optimality of a larger set of controls.

In the next section, we apply the generalized Legendre-Clebsch condition to derive a high-order necessary condition for a singular control to be optimal. We also combine known results on the second-order derivative of the Perron root [38] and the Agrachev-Gamkrelidze second-order variation for bang-bang controls (see, e.g., [39]) to derive a second-order MP for bang-bang controls. The proofs of these results are given in Section III.

## Ii Main results

Our first result is a high-order necessary condition for singular optimal controls for Problem 1. Without loss of generality (see [40]), we assume that the singular control is . Let  denote the Lie-bracket of the matrices .

### Ii-a High-order MP for singular controls

###### Theorem 3

Consider the PBCS (4). Suppose that the conditions of Thm. 2 hold, and that  is an optimal control. Then

 (w∗)′[B,[B,A]]v∗≤0. (16)
###### Example 3

Consider the specific PBCS with  given in Example 2. In this case,

 [B,[B,A]]=[6.818.421.4−6.8],

and , so

 (w∗)′[B,[B,A]]v∗=20.

It follows from (16) that  is not an optimal control. Note that we were not able to derive this conclusion using the first-order MP in Thm. 2

### Ii-B Second-order MP for bang-bang controls

In this section, we derive an Agrachev-Gamkrelidze-type second-order MP for optimal bang-bang controls for Problem 1. Note that for an optimal bang-bang  we have

 C∗(T)=exp((A+B)τk)…exp((A+B)τ2)exp((A−B)τ1),

with  and . Any cyclic shift of , e.g.,

 exp((A−B)τ1)exp((A+B)τk)…exp((A+B)τ2)

also corresponds to an optimal control (as a product of matrices and its cyclic shift have the same spectral radius). This means that we can always assume that  is a switching point of , and then (3) implies that  is also a switching point of .

Let  denote the set of all vectors  satisfying

 α1+⋯+αk=0. (17)

We can now state the main result in this section.

###### Theorem 4

Suppose that  is an optimal control for Problem 1, that the conditions of Thm. 2 hold, and that the switching function (10) admits a finite number of zeros at , with , so that  for , for , for , and so on, with . Denote , and . Define matrices , , by

 H1 :=P, (18) H2 :=Q, H3 :=exp(−τ2Q)Pexp(τ2Q), H4 :=exp(−τ2Q)exp(−τ3P)Qexp(τ3P)exp(τ2Q), H5 :=exp(−τ2Q)exp(−τ3P)exp(−τ4Q)Pexp(τ4Q)exp(τ3P)exp(τ2Q), ⋮

Then

 q′(t1)k∑i=1αiHip(t1)=0,for all α∈Pk. (19)

Furthermore,

 rk(α):=q′(t1)∑1≤i

satisfies

 rk(α)≤0,for all α∈Qk, (20)

where

 Qk:={α∈Pk:k∑i=1αiHip(t1)=0}. (21)

We refer to the control  defined above as a control with  bang arcs. As will be shown in the proof, condition (19) is a first-order condition (that can also be derived using the first-order MP). Condition (20) however is a second-order condition, and it is meaningful for values  that make a certain first-order variation vanish, i.e. that belong to .

Note that the conditions in Thm. 4 are given in terms of  and . It is possible of course to state them in terms of  and , but this leads to slightly more cumbersome expressions.

The next example demonstrates the calculations for a control with two bang arcs.

###### Example 4

Consider an optimal control in the form

 u∗(t)={1,t∈(0,t1),−1,t∈(t1,T),

where . In this case, (19) becomes

 q′(t1)(α1(A+B)+α2(A−B))p(t1)=0,for % all α∈P2,

and the definition of  yields

 α1q′(t1)((A+B)−(A−B))p(t1)=0,for all α1∈R.

Of course, this is just the conclusion that we can get from the first-order MP, as at the switching point  we must have

 0=m(t1)=q′(t1)Bp(t1).

The second-order term is

 r2(α) =α1α2q′(t1)[H1,H2]p(t1) =−α21q′(t1)[A+B,A−B]p(t1) =2α21q′(t1)[A,B]p(t1),

so (20) becomes

 r2(α)≤0,for all α∈Q2, (22)

where

 Q2={α1∈R:α1Bp(t1)=0}.

Again, this provides information that can also be derived from the first-order MP, as the fact that  and  implies that

 ˙m(t1)≤0.

and differentiating (10) yields

 ˙m(t1)=q′(t1)[A,B]p(t1).

Thus, , so (22) actually holds for all

However, for a control with more than two bang arcs the second-order condition does provide new information. The next simple example demonstrates this.

###### Example 5

Consider the PBCS (I-A) with

 A=[−5/23/23−5/2],B=[3/2−1/21−3/2].

Note that  is Metzler for all . Consider the control

 u(t)={1,t∈(0,t1)∪(t2,t3),−1,t∈(t1,t2)∪(t3,T), (23)

with , , , and . The corresponding transition matrix is

 C(T)=exp(A−B)exp(A+B)exp(A−B)exp(A+B).

Let . The spectral radius of  is

 ρ=(9+7exp(5)+9exp(10)+3s(exp(5)−1)25exp(10))2,

and it is a simple eigenvalue. The Perron right and left eigenvectors of  are

 v=[exp(5)−1+s2+8exp(5)]′

and

 w=[exp(5)−1+s4+exp(5)]′.

Calculating the switching function  defined in (10) yields the behavior depicted in Fig. 1. Note that  for , and  for , so the control  satisfies the first-order MP.

We now show that the second-order MP implies that  is not an optimal control. Eq. (18) yields

 H1 =A+B, H2 =A−B, H3 =exp(−(A−B))(A+B)exp(A−B), H4 =exp(−(A−B))exp(−(A+B))(A−B)exp(A+B)exp(A−B).

Note that

 H3 =exp(−(A−B))(A−B+2B)exp(A−B) =A−B+2exp(−(A−B))Bexp(A−B). (24)

Our goal is to find  such that . Indeed, this will imply that  is not optimal. It turns out that we can find such an  satisfying  and . Since  must be zero, . Then

 4∑i=1¯αiHi =A+B+¯α2(A−B)−(1+¯α2)(A−B+2exp(−(A−B))Bexp(A−B)) =2B−2(1+¯α2)exp(−(A−B))Bexp(A−B),

so

 4∑i=1¯αiHip(t1) =2(B−(1+¯α2)exp(−(A−B))Bexp(A−B))exp(A+B)v =2Bexp(A+B)v−2√ρ(1+¯α2)exp(−(A−B))Bv.

A tedious but straightforward calculation shows that

 Bexp(A+B)v=exp(−5)exp(−(A−B))Bv,

so  for

 ¯α2=(√ρexp(5))−1−1.

Summarizing, . The second-order term is

 r4(¯α) =q′(t1)∑1≤i

and a calculation yields

 r4(¯α)=1050exp(5)(exp(5)−1)(12(s−3)+exp(5)(−67+s+exp(5)(67+36exp(5)+12s)))(4+exp(5))(1+4exp(5))(9−3s+exp(5)(7+9exp(5)+3s))2.

Clearly, , so the second-order MP implies that  in (23) is not an optimal control. The reason that  here satisfies the conditions of the first-order MP is that it actually minimizes the spectral radius at time . Thus, the second-order MP plays here a similar role to the second-derivative of a function: it allows to distinguish between a maximum point and a minimum point.

## Iii Proofs

### Iii-a Proof of Thm. 3

Assume that  is an optimal control. The corresponding solution of (I-A) is . For , consider the control

 ~u(t):=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩0,t∈[0,T−4ϵ1/3),−1,t∈[T−4ϵ1/3,T−3ϵ1/3),1,t∈[T−3ϵ1/3,T−ϵ1/3),−1,t∈[T−ϵ1/3,T). (25)

Then

 C(T,~u) =exp((A−B)ϵ1/3)exp((A+B)2ϵ1/3)exp((A−B)ϵ1/3)exp(A(T−4ϵ1/3)) =exp((A−B)ϵ1/3)exp((A+B)2ϵ1/3)exp((A−B)ϵ1/3)exp(−4ϵ1/3A)C∗(T),

and it follows from the computation in [40, p. 719] (see also [41]) that

 C(T,~u)=exp(2ϵ3[B,[B,A]])C∗(T)+o(ϵ). (26)

Note that this implies that any result derived using  will be a high-order MP, as the width of the needle variations in (25) is of order  yet the perturbation in  with respect to  is of order . By (26),

 ddϵC(T,~u)|ϵ=0=(2/3)[B,[B,A]]C∗(T),

so

 ρ(C(T,~u) )−ρ(C∗(T)) =(2ϵ/3)(w∗)′[B,[B,A]]C∗(T)v∗+o(ϵ) =(2ϵ/3)(w∗)′[B,[B,A]]ρ(C∗(T))v∗+o(ϵ).

If  then  for all sufficiently small , and this contradicts the optimality of . This proves (16).

### Iii-B Proof of Theorem 4

The proof is based on introducing a new control defined by a perturbation of the switching times  to , , …, . Here, and . Define  by for , for , and so on. Note that (17) implies that the time length of the perturbed control is

 ~tk−~t0=tk+s(α1+…αk)−t0=tk−t0=T.

Denote the corresponding transition matrix by . Note also that  for any , so . Our goal is to derive an expression for the difference  in the form

 e(s,α)=sz1(α)+s22z2(α)+o(s2), (27)

where  denotes a function  that satisfies .

Suppose for a moment that  [] for some