Electronic Journal of Qualitative Theory of Differential Equations Proc. 8th Coll. QTDE, 2008, No. 16 1-14; http://www.math.u-szeged.hu/ejqtde/

About stability on bounded domains of the state space

Vladimir Răsvan

Department of Automatic Control,

University of Craiova,

A.I.Cuza, 13, Craiova, RO-200585, Romania

             e-mail: vrasvan@automation.ucv.ro

Abstract

Modern aircraft control systems contain the so-called rate limiter elements which are in fact incorporating the saturation function which is a bounded nonlinearity. In certain critical cases the so-called sector rotation - a standard procedure in absolute stability - leads to the fact that the sector conditions are broken outside a bounded domain. In order to apply standard results of the absolute stability theory there will be combined the results of the hyperstability theory with those arising from Liapunov function theory since existence of a Liapunov function is up to now the best way to estimate stability domains. At the same time the stability conditions will be expressed in the language of some frequency domain inequality as required by the conditions of the practical problem that generated the mathematical one.

AMS Classification: 93D10, 34D20, 34D23 Keywords: absolute stability, bounded state space domain, frequency

domain inequality

This paper is in final form, has been submitted to the Proceedings of 8th Colloquium on the Qualitative Theory of Differential Equations and no version of it will be submitted for publication elsewhere.

1 A motivating problem

A genuine contemporary challenge for both engineers and mathematicians is the problem of the so called PIO - P(ilot) I(n the loop) O(scillations) of modern combat (but also civil) aircraft. Their mechanism which is better and better understood shows a self-sustained oscillation proneness of the feedback system pilot-aircraft: there exist situations when, paradoxically, the pilot’s efforts to control the aircraft result in some kind of de-stabilizing that generates uncontrollable self-sustained oscillations. Today the aircraft specialists consider three kinds of PIO that may be present in aircraft dynamics. Without discussing here their significance, we just mention that our paper deals with the theory of the so-called PIO-II: their defining conditions are characterized by the modeling assumption that all elements of the system are linear except the so-called rate limiters of the actuators.

A. The rate limiter is an engineering structure that has to be modelled. The common structure is by now that of fig.1

Figure 1: Rate-limiter modeling

Let us consider a standard classical example arising apparently from aircraft control [2]

ÿ + 2ζω_{n}ẏ + ω_{n}^{2}y = −K_{p}y  (1)

with 0 < ζ < 1, K_{p} > 0. Obviously this linear system is exponentially stable for all K_{p} > 0; here K_{p} could be the equivalent pilot gain. But the pilot command is usually applied through a rate limiter having the structure of fig.1 and this leads to a system of augmented dimension

\begin{matrix}ẅ + 2ζω_{n}ẇ + ω_{n}^{2}w = K_{p}z , z = Rξ \\ ξ̇ = sat (K_{s} (r/R − ξ)) , r = −w + v(t) \\ sat(e) = \frac{e}{max \{1, |e|\}} \end{matrix}  (2)

The dimension augmentation limits K_{p} > 0 even in the linear case: application of the Routh Hurwitz criterion for the linearized system

\begin{matrix}ẅ + 2ζω_{n}ẇ + ω_{n}^{2}w = K_{p}Rξ \\ ξ̇ = −K_{s}ξ − \frac{K_{s}}{R} (w − v(t)) \end{matrix}  (3)

will give

K_{p} < 2ζω_{n}(K_{s} + 2ζω_{n} + \frac{ω_{n}^{2}}{K_{s}})  (4)

In the nonlinear case we have to remember that saturation is a nonlinear function confined to a sector (fig.2) namely (0,1]

0 ≤ \frac{sat(e)}{e} = \frac{1}{max\{1, |e|\}} ≤ 1  (5)

Figure 2: Sector restricted nonlinear function

For v(t) ≡ 0 we obtain a system that is at equilibrium in 0

\begin{matrix}ẅ + 2ζω_{n}ẇ + ω_{n}^{2}w = K_{p}Rξ \\ ξ̇ = −φ(σ) \\ σ = \frac{K_{s}}{R} w + K_{s}ξ \end{matrix}  (6)

This system corresponds to the critical case of a simple zero root in the absolute stability problem. If the Popov frequency domain inequality is used

1 + ℜe (1 + ıωβ)γ(ıω) ≥ 0  (7)

where the transfer function γ(λ) reads here

γ(λ) = K_{s}(\frac{1}{λ} + \frac{K_{p}}{λ(λ^{2} + 2ζω_{n}λ + ω_{n}^{2})})  (8)

A tedious but straightforward computation will give the same inequality (4) as absolute stability condition; this means that (4) is a necessary and sufficient condition of absolute stability hence if K_{p} - the pilot gain in the aircraft dynamics interpretation - is larger than the value prescribed by the RHS of (4), self sustained oscillations will occur - the PIO-II type.

B. We shall consider now a more recent version of aircraft dynamics - the longitudinal short period motion controlled by two actuators - the elevon control and the canard control; we write down the model in deviations with the notations of the field

\frac{d}{dt}Δα = q

\frac{d}{dt} q = M_{α}Δα + M_{q}q + M_{δe}Δ_{e} + M_{δe}Δ_{c}

\frac{d}{dt}Δ_{e} = ω_{a}ψ(−k_{α}Δ_{α} − k_{q}q − Δ_{e})

\frac{d}{dt}Δ_{c} = ω_{a}ψ(−k_{α}Δ_{α} − k_{q}Δ_{q} − Δ_{c})

(9)

where ψ(σ) is the saturation function as previously. The characteristic equation of the uncontrolled motion is

λ^{2} − M_{q}λ − M_{α} = 0  (10)

with M_{q} < 0, M_{α} > 0: we are in the case of the unstable aircraft which displays a saddle point at the equilibrium - one of the eigenvalues is strictly positive.

Considering first the linearized case with ψ(σ) ≡ γσ, 0 < γ < 1 which has the characteristic equation

P_{γ}(λ) ≡ (λ + ω_{a}γ)[λ^{3} + (ω_{a}γ − M_{q})λ^{2} + (ω_{a}γA_{1} − M_{α})λ + ω_{a}γA_{0}] (11)

where we denoted: M_{δ} = M_{δe} + M_{δc}, A_{0} = M_{δ}k_{α} − M_{α}, A_{1} = M_{δ}k_{q} − M_{q}. The (gain) Hurwitz condition will be

ω_{a}γ > ξ_{+} = \frac{A_{0} + M_{α} + A_{1}M_{q} + \sqrt{(A_{0} + M_{α} + A_{1}M_{q})^{2} − 4A_{1}M_{q}M_{α}}}{2A_{1}} (12)

This shows that the true stability sector (in the nonlinear case) has to be expected at most

\frac{ξ_{+}}{ω_{a}} < \frac{ψ(σ)}{σ} < 1

In order to see the significance of this fact, introduce a new state variable - the elevon/canard unsymmetry

ζ = Δ_{e} − Δ_{c}

to obtain the new system

\frac{d}{dt}Δα = q

\frac{d}{dt} q = M_{α}Δα + M_{q}q + (M_{δe} + M_{δc})Δ_{e} − M_{δc}ζ

\frac{d}{dt}Δ_{e} = ω_{a}ψ(−k_{α}Δ_{α} − k_{q}q − Δ_{e})

\frac{dζ}{dt} = ω_{a}[ψ(−k_{α}Δ_{α} − k_{q}q − Δ_{e}) − ψ(−k_{αΔ_{α}} − k_{q}q − Δ_{e} + ζ)]

(13)

which has the invariant set ζ ≡ 0. If we consider our system confined to the invariant set, it will have a reduced (-1) order

\begin{matrix}\frac{d}{dt}Δα = q \\ \frac{d}{dt} q = M_{α}Δα + M_{q}q + (M_{δe} + M_{δc})Δ_{e} \\ \frac{d}{dt}Δ_{e} = ω_{a}ψ(−k_{α}Δ_{α} − k_{q}q − Δ_{e}) \end{matrix}  (14)

and will contain a single nonlinear element; the sector rotation

φ(σ) = −ξ_{+}σ − ω_{a}ψ(−σ)

will give the system

\frac{d}{dt}Δα = q

\frac{d}{dt} q = M_{α}Δα + M_{q}q + M_{δ}Δ_{e}

\frac{d}{dt}Δ_{e} = −ξ_{+}(k_{α}Δ_{α} + k_{q}q + Δ_{e}) − φ(k_{α}Δ_{α} + k_{q}q + Δ_{e})

0 < \frac{φ(σ)}{σ} < ω_{a} − ξ_{+}

(15)

One may see that the nonlinear function was confined to the sector with 0 as lower limit.

In fact since the initial nonlinear function - the saturation - had already 0 as lower sector limit, the sector of the rotated nonlinearity should have been (−ξ_{+}, ω_{a} − ξ_{+}) but we already know that (−ξ_{+}, 0) is not a stability sector. For this reason we have restricted ourselves to the positive sector (0, ω_{a} − ξ_{+}). Unfortunately the saturation nonlinearity breaks the positive sector for large deviations (fig.3) This implies the application of the absolute stability methods for bounded domains of the state space.

Figure 3: Sector rotation for saturation.

2 The mathematical problem and the main result

We shall consider the following system

ẋ = Ax − bφ(c^{∗}x)  (16)

with the usual notations of the absolute stability theory, where φ(σ) satisfies the sector condition

\underline{φ} ≤ \frac{φ(σ)}{σ} ≤ φ̄  (17)

(see also fig.2 where the notations are slightly different) and assume that there exists some φ_{0} ∈ [\underline{φ}, φ̄] such that A − φ_{0}bc^{∗} should be a Hurwitz matrix; this is some kind of minimal stability since if we want system (16) to be absolutely stable i.e. asymptotically stable for all functions satisfying (17) then we have to assume minimally this stability for at least a single linear characteristic of the sector.

A. In the following we shall define an auxiliary problem for system (16). If (17) hold then φ(σ) verifies

(φ(σ) − \underline{φ}σ)(φ̄σ − φ(σ)) ≥ 0  (18)

Moreover, if σ(t) is some differentiable function on some interval, the derivative being integrable on that interval, then we shall have

\begin{matrix}∫_{0}^{t}(\underline{φ}σ(τ) − φ(σ(τ)))\frac{dσ}{dt} (τ)dτ = \underline{Ψ}(σ(0)) − \underline{Ψ}(σ(t)) \\ ∫_{0}^{t}(φ(σ(τ)) − φ̄σ(τ))\frac{dσ}{dt} (τ)dτ = \overline{Ψ}(σ(0)) − \overline{Ψ}(σ(t)) \end{matrix}  (19)

where

\underline{Ψ}(σ) = ∫_{0}^{σ}(φ(λ) − \underline{φ}λ)dλ ≥ 0 , \overline{Ψ}(σ) = ∫_{0}^{σ}(φ̄λ − φ(λ))dλ ≥ 0 (20)

Consider now the linear controlled system

ẋ = Ax + bu(t)  (21)

and associate to it the following integral index

η(0, t) = ∫_{0}^{t} 𝔉(u(τ), x(τ))dτ  (22)

where 𝔉(u, x) is the following quadratic form of n + 1 variables

𝔉(u, x) = α_{0}(u + \underline{φ}c^{∗}x)(u + φ̄c^{∗}x) +

(23)

+ (α_{1}(u + \underline{φ}c^{∗}x) − α_{2}(u + φ̄c^{∗}x))(c^{∗}Ax + c^{∗}bu)

and the integral index is defined for any pair of integrable vector valued functions. If, additionally x(t) and u(t) satisfy (21) then we may integrate by parts in (22) to obtain

η(0, t) = \frac{1}{2}(α_{1}\underline{φ} − α_{2}φ̄)(c^{∗}x(τ))^{2}|_{0}^{t} + ∫_{0}^{t} 𝔊(u(τ), x(τ))dτ  (24)

where 𝔊(u, x) is the following quadratic form

𝔊(u, x) = α_{0}(u^{2} + (\underline{φ} + φ̄)uc^{∗}x + \underline{φ}φ̄(c^{∗}x)^{2}) + (α_{1} − α_{2})c^{∗}(Ax + bu) (25)

Following V.M. Popov [4] we associate to the system defined by (21) and (24) the so-called system’s characteristic function χ : ℂ × ℂ ↦ ℂ

χ(λ, σ) = \frac{1}{2}(α_{1}\underline{φ} − α_{2}φ̄)(λ + σ)γ(λ)γ(σ) + α_{0}χ_{1}(λ, σ) + (α_{1} − α_{2})χ_{2}(λ, σ)

χ_{1}(λ, σ) = 1 + \frac{1}{2}(\underline{φ} + φ̄)(γ(λ) + γ(σ)) + \underline{φ}φ̄γ(λ)γ(σ)

χ_{2}(λ, σ) = \frac{1}{2}(λγ(λ) + σγ(σ)) , γ(σ) = c^{∗}(σI − A)^{−1}b

(26)

All functions are rational since they depend on the strictly proper rational function γ(σ) - the transfer function of the linear system defined by the controlled system (21) with the output ν = c^{∗}x.

We shall consider now the positivity theory [4] - a generalization of the Yakubovich-Kalman-Popov Lemma - applied to the system defined by (21) and (24) with the characteristic function (26). If the frequency domain inequality below holds

χ(−ıω, ıω) = α_{0}(1 + (\underline{φ} + φ̄)ℜe γ(ıω) + \underline{φ}φ̄|γ(ıω)|^{2}) + βℜe ıωγ(ıω) ≥ 0 (27)

for some α_{0} ≥ 0, β ∈ ℝ, then there exist a scalar γ_{0}, a n-dimensional vector w and a n×n Hermitian matrix H such that the integral index η(0, t) should be given the form

η(0, t) = [\frac{1}{2}(α_{1}\underline{φ} − α_{2}φ̄)(c^{∗}x(τ))^{2} + x^{∗}(τ)Hx(τ)]|_{0}^{t} +

(28)

+ ∫_{0}^{t} |γ_{0}u(τ) + w^{∗}x(τ)|^{2}dτ

A small comment is necessary: if (27) holds for β ≥ 0 then we may take α_{1} ≥ 0, α_{2} = 0 and if (27) holds for β ≤ 0 then we may take α_{1} = 0, α_{2} ≥ 0.

The form (28) will turn to be useful in the construction of the suitable Liapunov function.

We shall make now use of the property of minimal stability in the sense of V.M. Popov [4]. Due to the existence of φ_{0} ∈ [\underline{φ}, φ̄] such that A − φ_{0}bc^{∗} is a Hurwitz matrix, we may chose in the system (16),(22)-(23) the control u as

u(t) = −φ_{0}c^{∗}x(t) + ρ(t),

where x(t) = e^{(A−φ_{0}bc^{∗})t}x_{0} and ρ(t) chosen according to [4], Chapter 5, in function of the various combinations of the free parameters α_{i} ≥ 0 (as they may follow from the fulfilment of the frequency domain inequality (27) for some α_{0} ≥ 0 and real β) to obtain η(0, t) ≤ 0. Therefore the matrix

P = H + \frac{1}{2}(α_{1}\underline{φ} − α_{2}φ̄)cc^{∗}  (29)

results nonnegative definite. If additionally, (A, b) is a controllable pair then, as shown in (op. cit.), Chapter 4, we have even P a positive definite matrix for all combinations of the free parameters α_{i} ≥ 0 (again, as they may follow from the fulfilment of the frequency domain inequality (27) for some α_{0} ≥ 0 and real β). This fact will turn extremely useful in the following.

B. We shall return now to the nonlinear system (16) and make firstly a rather obvious but useful remark: let z(t) be some solution of (16) corresponding to some initial condition z(0) = x_{0}; if we define u(t) = −φ(c^{∗}z(t)) using this solution and apply it to (21), the solution x(t) of (21) defined by x_{0} and u(t) will coincide with z(t).

Therefore we may consider the solutions of (16) as the solutions of (21) with the control input u(t) = −φ(c^{∗}x(t)). We thus consider the integral index (24) along the solutions of (21) with the control input defined as above. Taking into account (19) and (23) we find

η(0, t) = α_{1}\underline{Ψ}(ν(0)) − α_{1}\underline{Ψ}(ν(t)) + α_{2}\overline{Ψ}(ν(0)) − α_{2}\overline{Ψ}(ν(t)) +

(30)

+ α_{0} ∫_{0}^{t}(\underline{φ}ν(τ) − φ(ν(τ)))(\overline{φ}ν(τ) − φ(ν(τ)))dτ

with ν = c^{∗}x. We then equate (28) and (30) and use the notation (29) to find the so-called basic stability equality

x^{∗}(t)P x(t) + α_{1}\underline{Ψ}(c^{∗}x(t)) + α_{2}\overline{Ψ}(c^{∗}x(t)) =
= x_{0}^{∗}P x_{0} + α_{1}\underline{Ψ}(c^{∗}x_{0}) + α_{2}\overline{Ψ}(c^{∗}x_{0})−
−α_{0} ∫_{0}^{t}(\underline{φ}c^{∗}x(τ) − φ(c^{∗}x(τ)))(\overline{φ}c^{∗}x(τ) − φ(c^{∗}x(τ)))dτ−
−∫_{0}^{t} | − γ_{0}φ(c^{∗}x(τ)) + w^{∗}x(τ)|^{2}dτ

Introduce now the following state function that may be considered as a candidate Liapunov function

V (x) = x^{∗}P x + α_{1}\underline{Ψ}(c^{∗}x) + α_{2}\overline{Ψ}(c^{∗}x)  (31)

If (20) and (29) are taken into account then one can see that

V (x) = x^{∗}Hx + (α_{1} − α_{2})∫_{0}^{c^{∗}x} φ(λ)dλ  (32)

which is the standard Liapunov function “quadratic form plus integral of the nonlinear function” occurring in the absolute stability theory. The basic stability equality, which is written along the solutions of (16) becomes

V (x(t)) = V (x_{0}) − α_{0} ∫_{0}^{t}(\underline{φ}c^{∗}x(τ) − φ(c^{∗}x(τ)))(\overline{φ}c^{∗}x(τ) − φ(c^{∗}x(τ)))dτ −

(33)

− ∫_{0}^{t} | − γ_{0}φ(c^{∗}x(τ)) + w^{∗}x(τ)|^{2}dτ

what shows that V is at least non-increasing along the trajectories of (16).

We need now some information about the sign of V (x) itself. But V (x) is clearly positive definite since it has the form (31) hence V (x) ≥ x^{∗}P x due to the fact that α_{i} ≥ 0, \underline{Ψ}(ν) ≥ 0, \overline{Ψ}(ν) ≥ 0. Since P > 0 as proved we deduce V to be positive definite, being already at least non-increasing along the solutions of system (16). Therefore stability follows. For the asymptotic stability we need some additional information on the RHS of (33): if, for instance, the sector inequalities are strict, then V (x(t)) is strictly decreasing and we may apply the theorem of Liapunov on asymptotic stability to obtain this property. Since the inequality V (x) ≥ x^{∗}P x with P > 0 shows the Liapunov function to be radially unbounded, the asymptotic stability will be global.

We may summarize the results proved above in the following

Theorem 1. Consider system (16) with φ : ℝ ↦ ℝ subject to inequalities (17) assumed to be strict. If there exists φ_{0} ∈ (\underline{φ}, φ̄) such that A − φ_{0}bc^{∗} is a Hurwitz matrix and also the numbers α_{0} > 0 and β ∈ ℝ such that the frequency domain inequality (27) holds, then the equilibrium of (16) at the origin is globally asymptotically stable for all functions satisfying (17) with strict inequalities.

We observe, for the completeness of the results, that the sector conditions may be allowed to be non-strict provided the frequency domain inequality (27) is strict or system (16) has A dichotomic i.e. without eigenvalues on the imaginary axis. Details may be found in [4].

3 Application to the case of a bounded domain in the state space

A. The motivating application showed that it is possible for the sector conditions to hold only in a bounded domain of the state space, of the form |c^{∗}x| ≤ ξ̄. Under these circumstances all development of the previous section keeps its validity provided we remain in the state space domain defined by the above inequality. Consequently we need finding invariant subsets contained in |c^{∗}x| ≤ ξ̄. This is the reason why we took the Liapunov approach|c^{∗}x| ≤ ξ̄: indeed, the easiest to recognize invariant sets of the system are those of the form V (x) ≤ c where c > 0 is some constant. Therefore the set of interest to our application would be

𝔐 = sup_{c} {{x ∈ ℝ^{n} : V (x) < c} ⊂ {x ∈ ℝ^{n} : |c^{∗}x| < σ_{0}}}  (34)

with V the Liapunov function of (32), H being the Hermitian matrix whose existence is ensured by the frequency domain inequality (27) and which may be determined by solving Linear Matrix Inequalities of Lurie type [1] and the available MATLAB software. Observe that the supremum problem of (34) also can be solved by adequate software.

The algorithm described above solves in a satisfactory way the mathematical problem issued from the practical one. The practical problem itself requires some other additional features among which the necessity to express the result in the language of those system parameters which may be measured and evaluated by the customer (in the aircraft case - by the so called pilot ratings).

B. We turn to the equations of the aircraft application given by (15). The transfer function that is associated to the linear part of (15) is

ϑ(λ) = \frac{λ^{2} + A_{1}λ + A_{0}}{(λ + ξ_{+} − M_{q})(λ^{2} + A_{1}ξ_{+} − M_{α})}  (35)

with the poles ω_{1} = ξ_{+} − M_{q} > 0, ±ıω_{0}, ω_{0}^{2} = A_{1}ξ_{+} − M_{α} > 0 and A_{0} = M_{δ}k_{α} − M_{α} > 0 (from the numerical data: in fact k_{α} which is a design parameter is always chosen as such), A_{1} = M_{δ}k_{q} − M_{q} > 0.

The sector of interest being now (0, ω_{a} − ξ_{+}) the frequency domain condition (27) will be read as a standard Popov inequality

\frac{1}{ω_{a} − ξ_{+}} + ℜe(1 + ıωβ)ϑ(ıω) ≥ 0  (36)

corresponding to the critical case of a pair of purely imaginary poles. The unique choice for β is as follows

β̄ = (\frac{A_{0}ω_{1}}{ω_{0}^{2}} + A_{1} − ω_{1})(A_{1}ω_{1} + ω_{0}^{2} − A_{0})^{−1} > 0

and (36) reads

\frac{1}{ω_{a} − ξ_{+}} + \frac{β̄ω^{2} + A_{0}ω_{1}/ω_{0}^{2}}{ω_{1}^{2} + ω^{2}} ≥ 0  (37)

We may see that, in principle we may accept even infinitely large stability sectors since the second term in (37) is strictly positive for all ω ∈ ℝ and, therefore, the first term given by a finite sector condition is not necessary for the fulfilment of (37). This shows in fact that the system could bear infinite linear and nonlinear gain without loosing its global asymptotic stability. The problem is here the sector breaking as a consequence of the saturation type of the nonlinear function. It is this fact that may induce instability and requires estimate of the invariant sets enclosed in the domain where the sector restrictions are observed.

And we thus arrive to the problem of the best Liapunov function. A rather general way of constructing it has been indicated above, involving theory and practice of Linear Matrix Inequalities. But, since the system is of low order, revisiting a classical book containing analytical methods for Liapunov functions construction in the problem of the absolute stability [3] might be rewarding.

References

[1] S. Boyd, L. El Ghaoui, E. Feron and V. Balakrishnan, Linear matrix inequalities in system and control theory, SIAM, Philadelphia, 1994.

[2] D. Graham and D. T. McRuer, Analysis of nonlinear control systems, J. Wiley,& Sons Inc., 1961.

[3] A.M. Letov, Stability of Nonlinear Controls (in Russian), Gostekhizdat, Moscow, 1955 (English version by Princeton Univ. Press, 1961); second edition, 1962 (English version, 1963)

[4] V. M. Popov, Hyperstability of Control Systems (in Romanian), Editura Academiei, Bucharest, 1966 (English version by Springer Verlag, 1973).

(Received August 15, 2007)