Abderrahman Iggidr
email: abderrahman.iggidr@inria.fr

CONTROLLABILITY, OBSERVABILITY AND STABILITY OF MATHEMATICAL MODELS

Keywords: accessibility, asymptotic stability, attractivity, chemostat, closed-loop, controllability, differential equation, finite dimensional systems, input, Lie algebra, Lotka-Volterra systems, Lyapunov functions, nonlinear systems, observability, observer, output, stabilization, state. Controllability, Observability

Bounded set: A subset B of IR n is bounded if there exists a positive real number M such that x ≤ M for all x ∈ B.

Closure of a set: The closure Ā of a set A is the smallest closed set containing A.

Connected set: A subset Ω of IR n is said to be connected if it can not be represented as the union of two disjoint non-empty open sets of IR n .

State vector: Vector whose elements are the state variables of a dynamical system.

State space:

The space to which the state vector belongs.

Continuous time systems: Systems whose evolution is described by differential equations.

Discrete time systems: Systems whose evolution is described by difference equations.

Time invariant systems: Dynamical systems for which the time derivative of the state vector at a time t depends only on the state vector at the current time t and does not depend explicitly on time. For instance, the parameters of the model of a time-invariant system are constants.

Vector field: A vector field X on IR n is a function that associates a vector X(x) ∈ IR n to every point x ∈ IR n . A vector field X is always associated to a differential equation ẋ = X(x).

Characteristic polynomial:

The characteristic polynomial of a square n × n matrix A is defined by P A (X) = det(A -XI n ) where I n is the n × n identity matrix and det(M ) is the determinant of the square matrix M .

Eigenvalue: A scalar λ is an eigenvalue of a square matrix A if there exists a vector v = 0 such that A.v = λv. The eigenvalues of a square matrix A are the roots of its characteristic polynomial.

Introduction

The main goal of this article is to develop in some details some notions of Control Theory introduced in the article (Basic principles of mathematical modeling). It concerns more specifically three structural properties of control systems: controllability, stability and observability. Based on the references listed in the end of this article, we give a survey of these three properties with applications to various LSS examples. By a control system we mean a dynamical system evolving in some state space and that can be controlled by the user. More precisely we are interested in the study of systems that can be modeled by differential (respectively difference) equations of the form ẋ(t) = dx dt = X(x(t), u(t)),

y(t) = h(x(t)), (1)
where the variable t represents the time, the vector x(t) is the state of the system at time t, the vector u(t) is the input or the control, i.e., the action of the user or of the environment and the vector y(t) is the output of the system that is, the available information that can be measured or observed by the user. The dynamics function X indicates how the system changes over time.

We shall address the following problems:

• When the whole state x(t) is not available for measurement, how is it possible to use the information provided by y(t) together with the dynamics (1) in order to get a "good" estimation of the real state of the system? This turns out to be an observability problem. • How to use the control u(t) in order to meet some specified needs? This is the problem of controllability and stabilizability. To illustrate these various concepts, we explain them through a simple LSS system: an epidemic model for the transmission of an infectious disease. An homogeneous population is divided into four classes S, E, I and T according to the heath of its individuals. Let S(t) denote the number of individuals who are susceptible to the disease, i.e., who are not yet infected at time t. E(t) denotes the number of members at time t who are exposed but not yet infected. I(t) denotes the number of infected individuals, that is, who are infectious and able to spread the disease by contact with individuals who are susceptible. T (t) is the number at time t of treated individuals. The total population size is denoted N = S + E + I + T . The dynamic of the disease can be described by the following differential system:

                           dS dt = bN -µS -β SI N , dE dt = β SI N -(µ + )E, dI dt = E -(r + d + µ)I, dT dt = rI -µT, dN dt = (b -µ)N -dI, (2)
where the parameter b is the rate for natural birth and µ that of natural death. The parameter Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 4

β is the transmission rate, d is the rate for disease-related death, is the rate at which the exposed individuals become infective and r is the per-capita treatment rate. The parameters b, µ, β, d and are assumed to be constant. The per-capita treatment rate may vary with time. The system (2) can be written by using the fractions s = S/N , e = E/N , i = I/n and τ = T /N of the classes S, E, I and T in the population. These fractions satisfy the system of differential equations: (3) Since s + e + i + τ = 1, it is sufficient to consider the following system

                  
             ds dt = b -bs -(β -d)si, de dt = βsi -(b + )e + dei, di dt = e -(r + d + b)i + di 2 .
(4)

If we suppose that the proportion of infected individuals can be measured and that we can control the population by choosing the treatment rate then the above system can be seen as a control system of the form (1) with: the state of the system is x(t) = (s(t), e(t), i(t)), the control is u(t) = r(t) and the measurable output is y(t) = i(t). The state space is

Ω = {x ∈ IR 3 : 0 ≤ s ≤ 1, 0 ≤ e ≤ 1, 0 ≤ i ≤ 1, s + e + i ≤ 1}.
The aim of Control Theory is to give answers to the following questions:

1. Given two configurations x 1 and x 2 of the system (3), is it possible to steer the system from x 1 to x 2 by choosing an appropriate treatment strategy r(t) ? For a given initial state x 0 = (s 0 , e 0 , i 0 ) and a given time T , what are all the possible states that can be reached in time T by using all the possible control functions r(t)? 2. Does the system (3) possess equilibrium points? Are they stable or unstable? How can the control be chosen in order to make a given equilibrium (for instance, the disease-free equilibrium) globally asymptotically stable? 3. Since only some variables can be measured (for instance, the proportion i(t) of infected individuals in the above epidemic example), is it possible to use them together with the dynamics (3) of the system in order to estimate the non measurable variables s(t) and e(t)? and how this can be done? The questions 1 involve the controllability of control systems. We shall give in Section 2 a general definition of this notion as well as some simple useful criteria that allow to study the controllability of some nonlinear systems. The questions 2 invoke the stability and the stabilization of nonlinear system, this will be the subject of Section 3 where some classical stability results are exposed. The questions 3 are connected to the observability and the construction of observers for dynamical systems. These problems will be addressed in Section 4 and Section 5.

In each section we shall apply the different tools to various models such as predator-prey systems, fisheries and bioreactors. Concerning the epidemic model (2-4), the divers problems mentioned above are still under investigation. We can give here just some partial answers. It can be shown that if the per-capita treatment rate satisfies

r ≥ r 0 = β + b -d -b (5)
then the system has a unique equilibrium point in the domain Ω. This equilibrium is the disease free equilibrium s = 1, e = i = τ = 0 and it is globally asymptotically stable, that is, if (s(0), e(0), i(0)) ∈ Ω is any initial condition then the solution of the system (4) starting from this state will converge to the disease free equilibrium, i.e., s(t) → 1, e(t) → 0 and i(t) → 0 as t → ∞. Therefore, if the treatment satisfy the condition (5) then the disease will die out. If r < r 0 then the disease free equilibrium becomes unstable and in this case, there is another equilibrium which belongs to the interior of the domain Ω. This equilibrium is the endemic equilibrium and it is globally asymptotically stable within the interior of Ω provided that r < r 0 . The disease will be endemic in this case.

Controllability

In this section, we present an elementary overview of an important property of a system, namely that of controllability. This concept has been briefly introduced and explained in (Basic principles of mathematical modeling) . Here, we deal with the problem of testing controllability of systems of the type ẋ(t) = X(x(t), u(t)).

(6)

The state vector x(t) belongs to the state space M which will always be here IR n or an open connected subset of IR n . The control functions t → u(t) are defined on [0, ∞) and take values in a connected subset U of IR m . These control functions are assumed to belong to an admissible control set U ad . This admissible control set is generally specified by the problem considered.

Here it is assumed that U ad contains the set of piecewise constant functions as a dense subset, i.e., any admissible control function can be "approximated" by piecewise constant functions.

The dynamics function X : M ×IR n -→ IR n is assumed to be analytic on (x, u). For each control value u ∈ U , we denote by X u the vector field ("vector function") defined by X u (x) = X(x, u) for all x ∈ M .

A real-valued function f defined on some open set O ⊂ IR is said to be analytic if for each x 0 ∈ O, there exists a positive real number r > 0 such that f (x) is the sum of a power series for all x with |x -x 0 | < r:

f (x) = ∞ n=0
a n (x -x 0 ) n for all x satisfying |x -x 0 | < r.

For an analytic function, the coefficients a n can be computed as

a n = f (n) (a) n!
In a similar way a function f : V ⊂ IR k → IR of several variables is said to be analytic if it is locally given by power series. A vector function

X : V ⊂ IR k → IR n x = (x 1 , . . . , x k ) → X(x) =    X 1 (x 1 , . . . , x k ) . . . X n (x 1 , . . . , x k )   
is analytic if all its components X i are analytic.

For each x ∈ M and each control function u(.) ∈ U ad , we denote by X u(.) t (x) the solution of (6) satisfying X u(.) 0 (x) = x. We assume moreover that X u(.) t (x) is defined for all t ∈ [0, ∞). For instance, if X is a linear vector function, that is, X(x, u) = Ax + Bu, with A and B being matrices with appropriate dimensions, then X u(.) t (x) = e tA x + t 0 e (t-s)A Bu(s)ds.

What is controllability

A state x 1 ∈ M is reachable from x 0 ∈ M at time T ≥ 0 if there exists a control u(.) ∈ U ad such that x 1 = X u(.)
T (x 0 ). The set of reachable states from x 0 at time T will be denoted A T (x 0 ) and the reachable set from x 0 will be denoted

A(x 0 ) = T ≥0 A T (x 0 ).
For the class of piecewise constant controls, the reachable set at time T from x 0 is the set of points of M of the form:

x = X up tp o . . . oX u 2 t 2 oX u 1 t 1 (x 0 ), u 1 , . . . , u p ∈ U, t 1 + . . . + t p = T.
This formula means that the constant control u i is applied during time t i , i.e., u(t

) = u 1 for t ∈ [0, t 1 ), u(t) = u 2 for t ∈ [t 1 , t 1 + t 2 ), . . . , u(t) = u i for t ∈ [t 1 + . . . t i-1 , t 1 + . . . t i-1 + t i ).
The system (6) will be said to be accessible from x 0 if the reachable set A(x 0 ) has a nonempty interior in the state space M . When the set A(x 0 ) is equal to the whole state space M then the system is controllable from x 0 and it is controllable (or completely controllable) if this property holds for any x 0 ∈ M . More precisely, the system (6) is controllable if any point of M is reachable from any other point of M , i.e., for any two states x 0 , x 1 ∈ M there exist a finite time T (that may depend on (x 0 , x 1 )) and an admissible control function u(.)

: [0, T ] → U such that x 1 = X u(.)
T (x 0 ). The system (6) is strongly controllable if for any given time T > 0 any point of the state space M is reachable from any other point of M in T or fewer units of time.

Before giving some criteria for accessibility and controllability of control systems, we give a simple LSS example to illustrate these notions.

Example: Let us consider a simple Lotka-Volterra system that models a population of a harvested prey in the presence of a predator Ṅ = aN -bN P -q uN Ṗ = -cP + ebN P

In these equations N and P represent, respectively, the prey and the predator populations, a is the prey growth rate, b is the predator attack rate, c is the predator mortality rate, and e is the conversion efficiency of predators. The parameter q represents the catchability coefficient of the prey and u is the harvesting effort. The term quN is the rate of harvest when the effort is u and the prey stock is N . All the variables and the constants are positive. Here the state is x = (N, P ), the state space is the positive quadrant M = IR 2 + = {(N, P ) ∈ IR 2 : N > 0, P > 0}, the effort u represents the harvesting policy and can be seen as an input or a control term. We suppose that u can take two values 0 or 1. Therefore our admissible control set U ad is simply the set of piecewise constant functions defined on [0, ∞) with values in U = {0, 1}.

The dynamics vector function is X(x, u) = aN -bN P -q uN -cP + ebN P . The values of the control u generate two vector fields X 0 (x) = aN -bN P -cP + ebN P and X 1 (x) = aN -bN P -qN -cP + ebN P .

The trajectory of a vector field X through x ∈ M is the curve

{X t (x), t ∈ [0, ∞)} where X t (x)
is the solution of ẋ = X(x) with initial condition X 0 (x) = x (we suppose that the solution is defined for all t ≥ 0). Here, the trajectories of the vector fields X 0 and X 1 in the interior of the positive quadrant are closed curves (see Figure 2). More precisely the trajectories are the level curves of, respectively, the real-valued functions

V 0 (N, P ) = e b N -c log( e b c N ) + b P -a log( b a P ) -(a + c) V 1 (N, P ) = e b N -c log( e b c N ) + b P -(a -q) log( b a-q P ) -(a -q + c).
The horizontal and the vertical axes are particular trajectories for both X 0 and X 1 . Now, given an initial state x 0 , what is the reachable set from x 0 ? if x 0 is of the form x 0 = (N 0 , 0) then the reachable set is A(x 0 ) = {(N, 0) : N ≥ N 0 } and if x 0 = (0, P 0 ) then A(x 0 ) = {(0, P ) : P ≤ P 0 } since the vertical and the horizontal axes are positively invariant sets for the system (7). The interior M of the positive quadrant is positively invariant. Therefore if x 0 belongs to M then the whole trajectory starting from x 0 is contained in M and so A(x 0 ) ⊂ M . The reachable set of such state x 0 is actually the whole set M . For, let x 1 be any point of M , we shall show that x 1 ∈ A(x 0 ). There exist two positive numbers r 0 and r 1 such that the point x 0 belongs to the level curve C r 0 of V 0 defined by V 0 (x) = r 0 and the point x 1 belongs to the set C r 1 : V 0 (x) = r 1 . It can be proved that along each level curve C r of X 0 , there is a point where the vector field X 1 points inward C r and a point where the vector field X 1 points outward C r . This implies that the control u = 1 allows to pass from a given level set C r to a "smaller" one as well as to a "larger" one. Therefore one can go from the level curve C r 0 to the level curve C r 1 by using a piecewise constant control u = 0 or u = 1. Hence the system can be driven from the state x 0 to the state x 1 . The system is then completely controllable on the state space M . It must be noticed that there are many other possible strategies to achieve the same goal, for instance, the same reasoning can be done by considering the level curves of X 1 instead of those of X 0 . The figure 3 illustrates possible control schemes to go from x 0 to x 1 .

N P u=0 0 t 2 3 o X t 1 o X t 0 1 x 0 ( ) X x 1 0 2 3 o X 1 o X 1 x 0 ( ) X x 1 τ τ τ 1 -- --
Figure 3: Two control strategies to drive the system from the state x 0 to the state x 1 .

Controllability of linear systems

We have seen that the above system is controllable on the positive quadrant. To show this we have used the geometrical properties of the trajectories of the vector fields generated by the different values of the control. Unfortunately, it is not often possible to know explicitly the phase portrait of a vector field. So there is a need to have calculable criteria for testing controllability or at least accessibility of a given system. For linear systems, there is a simple controllability criterion known as Kalman's controllability rank condition. A linear system is a system governed by

ẋ = X(x, u) = Ax + Bu, x ∈ IR n , u ∈ U ⊂ IR m , (8)
where A and B are respectively n × n and n × m matrices. We start by exploring the controllability properties when the control set is U = IR m , i.e., there are no restrictions on the size of controls. It can be proved that the following are equivalent: the linear system ( 8) is (i) controllable, (ii) accessible from any point x 0 ∈ IR n , (iii) accessible from the origin x 0 = 0. The reachable set from the origin A(0) is a linear subspace of IR n . More precisely, it is the image of the linear map:

IR m+n -→ IR n (u 1 , . . . , u n ) -→ (B, AB, A 2 B, . . . , A n-1 B)    u 1 . . . u n    = n i=1 A i-1 Bu i .
The n × nm block matrix R(A, B) = (B, AB, A 2 B, . . . , A n-1 B) whose columns are the columns of B, AB, . . . , A n-1 B is called the Kalman controllability matrix. The controllability of the linear system ( 8) is related to the rank of this matrix as follows Theorem 2.1 The system (8) is controllable if and only if the n × nm Kalman controllability matrix

R(A, B) = B, AB, A 2 B, . . . , A n-1 B
is of rank n (the dimension of the state space).

In this case we also say that the pair (A, B) is controllable. It must be noticed that for a controllable linear system, every state x 1 can be reached from any state x 0 in a time interval of arbitrary length (i.e., as small as one wants) provided that the control values set U is unbounded. For linear systems with unbounded control values set, there is an equivalence between controllability and strong controllability.

When the linear system (8) is not controllable, we have rank(R(A, B)) = r < n. Let c 1 , c 2 , . . . , c r be r linearly independent columns of the matrix R(A, B) and let e 1 , . . . , e n-r be any vectors of IR n such that {c 1 , c 2 , . . . , c r , e 1 , . . . , e n-r } is a basis of IR n . Then the matrix T = (c 1 , c 2 , . . . , c r , e 1 , . . . , e n-r ) is an n × n invertible matrix and with the linear change of coordinates x = T z, the system ( 8) is transformed into

   ż1 = A 11 z 1 + A 12 z 2 + B 1 u, ż2 = A 22 z 2 , z 1 ∈ IR r , z 2 ∈ IR n-r (9)
where the pair (A 11 , B 1 ) is a controllable pair of matrices. We can remark that the control cannot act on the z 2 component of the state. The eigenvalues (or the eigenvectors) of the matrix A 22 are called the "uncontrollable modes" of the system.

Controllability with bounded controls: The Kalman's rank condition is a necessary condition for controllability of linear systems. It is a sufficient condition when the control values set U is unbounded. When there is a bound on the magnitude of the control u, the rank condition is no more sufficient to get controllability as it is illustrated by the following linear system

ẋ = x + u, x ∈ IR, u ∈ IR. ( 10
)
The rank controllability condition is obviously satisfied for this system. hence it is controllable and even strongly controllable if the control u can take any real value. Now, suppose that the admissible controls are bounded by a positive constant K, that is, for all u(.) ∈ U ad , one has |u(t)| ≤ K for all positive time t. The solution of (10) is

x(t) = e t x(0) + t 0 e -s u(s)ds .
Since |u(s)| ≤ K, we have t 0 e -s u(s)ds ≤ K. It follows that x(t) > K for all positive time t if x(0) > 2K, and therefore the system is not controllable on IR.

The above example shows that controllability of a linear system cannot in general be expected if there are restrictions on the size of controls. There is actually an additional necessary condition for controllability when the set U is bounded: Theorem 2.2 a. Suppose that the control set U is bounded. Then a necessary condition for controllability of the linear system (8) is that every eigenvalue of A have its real part equal to zero. b. If this necessary condition is satisfied and the set U is a neighborhood of the origin in IR m , then the linear control system with controls in U is controllable whenever the Kalman's rank condition holds, i.e., rank (B, AB, A 2 B, . . . , A n-1 B) = n.

It is interesting to remark that • the controllability of linear systems with unbounded controls is a "generic" property of linear systems. This means, roughly speaking, that almost all the pairs (A, B) are controllable when U = IR m . More precisely, if a given pair (A, B) does not satisfy the rank condition then by a "small perturbation" of the entries of the matrices, the pair becomes controllable; • linear systems are generically not controllable with bounded control because the necessary condition is easily violated when the entries of the matrix A are perturbed.

Accessibility criteria for nonlinear systems

For each control value u ∈ U , let X u = X(., u) be the vector function (vector field) corresponding to the constant control u. We then have a family F = {X u , u ∈ U } of vector fields parametrized by the values of the control. To get an accessibility criterion for the system (6) analogous to the Kalman's rank condition for linear systems we need to introduce some new vector fields generated by the elements of F. For any vector functions X and Y , it is possible to define a new vector function [X, Y ], called the Lie bracket of X and Y , and defined by

∀x ∈ IR n , [X, Y ](x) = ∂Y ∂x (x)X(x) - ∂X ∂x (x)Y (x)
where ∂Y ∂x (x) is the Jacobian of the vector function evaluated at the point x. In a coordinates system (x 1 , x 2 , . . . , x n ), if

X(x) =      X 1 (x 1 , x 2 , . . . , x n ) X 2 (x 1 , x 2 , . . . , x n ) . . . X n (x 1 , x 2 , . . . , x n )      and Y (x) =      Y 1 (x 1 , x 2 , . . . , x n ) Y 2 (x 1 , x 2 , . . . , x n ) . . . Y n (x 1 , x 2 , . . . , x n )      then, [X, Y ](x) =      Z 1 (x 1 , x 2 , . . . , x n ) Z 2 (x 1 , x 2 , . . . , x n ) . . . Z n (x 1 , x 2 , . . . , x n )      with Z i (x) = j=n j=1 ∂Y i ∂x j (x)X j (x) - ∂X i ∂x j (x)Y j (x).
For example, let X be a linear vector function, say, X(x) = Ax with A an n × n matrix and let Y be a constant vector function, that is,

Y (x) = b ∈ IR n . Then [X, Y ] is a constant vector function and [X, Y ](x) = -A b, for all x ∈ IR n .
A geometric interpretation of the Lie bracket: Let us consider a control system

ẋ = u 1 X 1 (x) + u 2 X 2 (x) x ∈ IR n , (u 1 , u 2 ) ∈ IR 2 . ( 11
)
Let x 0 ∈ IR n be a given state such that the vectors X 1 (x 0 ) and X 2 (x 0 ) are linearly independent, and let E(x 0 ) = span{X 1 (x 0 ), X 2 (x 0 )} be the vector subspace of IR n spanned by the vectors X 1 (x 0 ) and X 2 (x 0 ). From the point x 0 , we can steer in all directions contained in E(x 0 ) by using appropriate constant controls. Is it possible to steer in another direction that does not belong to E(x 0 )? If the answer is no then the system can never be accessible from x 0 . It turns out that one also can steer in the direction of [X 1 , X 2 ](x 0 ) by choosing appropriate piecewise constant control functions. So, if [X 1 , X 2 ](x 0 ) does not belong to E(x 0 ), then one has a third steering direction. This shows that the Lie bracket plays a prominent role in determining the reachable set from a given state x 0 .

To the control system (6) and the generated family of vector functions F we associate a set of vector functions called the accessibility Lie algebra of the system (6) and denoted Lie(F). The set Lie(F) is the Lie algebra of vector functions generated by the family F: every element of Lie(F) is a combination of repeated Lie brackets of the form

[X u k , [X u k-1 , [. . . , [X u 2 , X u 1 ] . . .]]], with u i ∈ U, X u i ∈ F.
For instance, if u 1 , u 2 are any two control values then the vector functions

[X u 1 , X u 2 ], [X u 1 , [X u 1 , X u 2 ]], [X u 2 , [X u 1 , [X u 1 , X u 2 ]]],
are in Lie(F).

The control system (6) satisfies the accessibility rank condition at the point x 0 if the set

Lie(F)(x 0 ) = {X(x 0 ), X ∈ Lie(F)}
is a vector space of dimension n (the dimension of the state space). This rank condition is analogous to the Kalman's rank condition. It allows to get the following accessibility criterion for analytic nonlinear systems:

Theorem 2.3 The analytic control system (6) is accessible from the point x 0 if and only if the accessibility rank condition holds at x 0 .

Remarks: 1. If the control system is not analytic but only smooth (C ∞ ), then the rank condition is still sufficient but not necessary. 2. For the linear system (8) with U = IR m , it is easy to see that the elements of Lie(F)(x 0 ) are generated by elements of the form A k Bu, k = 0, . . . , n -1, or Ax 0 + Bu, where u ∈ IR m . Therefore, for the linear system (8), the accessibility rank condition at the origin reduces to the Kalman controllability rank condition.

Example: Consider the following system defined on IR 3 :

ẋ = u 1   -2x 3 e x 1 e x 1 0   X 1 (x) +u 2   -2(x 2 + 1) 0 1   X 2 (x) , x = (x 1 , x 2 , x 3 ), (u 1 , u 2 ) ∈ IR 2 . ( 12
)
For any state x ∈ IR 3 , the vectors X 1 (x) and

X 2 (x) are linearly independent but [X 1 , X 2 ](x) = 2(x 2 + 1)X 1 (x). Since [X 1 , X 2 ] is collinear with X 1 , Lie(X 1 , X 2
) has a constant dimension equal to 2. Therefore, since the system ( 12) is analytic and the rank condition fails, it is is not accessible from any point x ∈ IR 3 , that is, the reachable set A(x) has an empty interior whatever what the control is.

This fact can be seen directly by performing a change of coordinates that allows to rewrite the system (12) in a simpler form. Let Φ : IR 3 -→ IR 3 the map defined by

x = (x 1 , x 2 , x 3 ) → Φ(x) = y = (y 1 , y 2 , y 3 ) with y 1 = x 1 + 2x 2 x 3 + 2x 3 , y 2 =
x 2 and y 3 = x 3 . The map Φ is bijective and its inverse is given by: Φ -1 (y 1 , y 2 , y 3 ) = (y 1 -2y 2 y 3 -2y 3 ). It is clear that the maps Φ and Φ -1 are continuously differentiable. Therefore Φ is a global diffeomorphism and allows to define a new coordinate system y = Φ(x). We also make the inputs transformation

v 1 = exp(y 1 -2y 2 y 3 -2y 3 )u 1 , v 2 = u 2 .
Hence in the new coordinates, the system ( 12) is expressed by:

   ẏ1 = 0 ẏ2 = v 1 ẏ3 = v 2
It is clear that the reachable set from any state x ∈ IR 3 is contained in a vector space of dimension 2 and, hence, it is of empty interior.

Control affine systems: A control affine system is a control system of the form

ẋ = X 0 (x) + u 1 X 1 (x) + . . . + u m X m (x), x ∈ M, u = (u 1 , . . . , u m ) ∈ U ⊂ IR m , (13)
where X 0 , X 1 , . . . , X m are analytic vector fields on the state space M , and u 1 , . . . , u m are the control functions. The vector field X 0 is called the drift and the vector fields X 1 , . . . , X m are called controlled vector fields. We assume that the control values set U contains m linearly independent points of IR m . For example if m = 2 then this assumption is fulfilled if U contains the two points (1, 0) and (0, 1). The interest of this assumption is that the Lie algebra generated by the family F corresponding to the control system (13) is independent of U and it is equal to the Lie algebra generated by the vector fields X 0 , X 1 , . . . , X m , that is,

Lie(F) = Lie{X 0 , X 1 , . . . , X m }.
For control affine systems (13), with U satisfying the above assumption, a controllability criterion is available:

Theorem 2.4 Assume that dim(Lie{X 1 , . . . , X m }(x)) = n = dimM for all x ∈ M . Then
(i) the control affine system (13) is controllable whenever there are no restrictions on the size of controls, (ii) the driftless (X 0 = 0) system ( 13) is controllable even if there are restrictions on the size of controls, provided that the convex hull of the constraint set U is a neighborhood of the origin in IR m .

Example: Consider the following simplified model of maneuvering a car:

  ẋ1 ẋ2 ẋ3   = u 1   0 0 1   X 1 (x) +u 2   sin x 3 cos x 3 0   X 2 (x) , x = (x 1 , x 2 , x 3 ) ∈ IR 3 , (u 1 , u 2 ) ∈ U ⊂ IR 2 . ( 14
)
The center of the front axle has coordinates (x 1 , x 2 ) ∈ IR 2 , while the rotation of this axle is given by the angle x 3 . The controls are the steering wheel moves u 1 and the engine speed u 2 . This system is a driftless control affine system. It is reasonable to suppose that the controls take values in a bounded neighborhood of the origin in IR 2 . We have dim

(Lie{X 1 , X 2 }(x)) = 3 for all x ∈ IR 3 since [X 1 , X 2 ](x) =   -cos x 3 sin x 3 0   and for all x ∈ IR 3 , rank{X 1 (x), X 2 (x), [X 1 , X 2 ](x)} = 3.
Therefore, Theorem 2.4 applies and hence the system is controllable.

Remark: It must be noticed that the sufficient controllability rank condition of the above theorem 2.4 invokes the Lie algebra generated only by the controlled vector fields X 1 , . . . , X m . This is a stronger condition than the accessibility rank condition that invokes the Lie algebra Lie{X 0 , X 1 , . . . , X m }. A control affine system may satisfy the accessibility rank condition without being controllable as illustrated by the following two-dimensional example:

   ẋ1 = u, ẋ2 = x 2 1 , (x 1 , x 2 ) ∈ M = IR 2 , u ∈ U = IR. (15)
Here, The drift vector field is

X 0 = 0 x 2 1
, and the controlled vector field is

X 1 = 1 0 . It is obvious that for all x ∈ IR 2 , dim(Lie{X 1 }(x)) = 1 . We have [X 0 , X 1 ] = 0 2x 1 , and [X 1 , [X 0 , X 1 ]] = 0 2 . Therefore, dim(Lie{X 0 , X 1 }(x)) = 2 = dimM
everywhere, and hence, according to Theorem 2.3, the system is accessible. However, it is not controllable since the reachable set from a given point x 0 = (x 1 0 , x 2 0 ) is contained in the half-space x 2 ≥ x 2 0 (see Figure 4). Theorem 2.4 says that controllability with bounded controls is possible for driftless control affine systems, i.e., X 0 = 0. There also exists a controllability criterion with bounded controls for the system (13) in the presence of a drift vector field X 0 which is not equal to zero but which is positively Poisson stable. To simplify matters, we assume that the state space M is an open connected subset of IR n . Let X be an analytic vector field on M and let for each point x ∈ M , X t (x) be the solution of the differential equation ẋ = X(x) that satisfies X 0 (x) = x. We assume that the vector field X is complete, which means that for any point x ∈ M , the solution

X t (x) is defined for all t ∈ [0, ∞). A point x ∈ M is said to be positively Poisson stable for X if ∀ > 0, ∀T > 0, ∃t ≥ T : X t (x) -x < .
The vector field X is called positively Poisson stable if the set of positively Poisson stable points for X is dense in M . A set N ⊂ M is dense in M if for any point p ∈ M and any number > 0, there exists q ∈ N such that p -q < . This means "near any point of M , there are points of N".

The linear vector field X(x) = Ax with A = 0 1 -1 0 is a positively Poisson stable vector field on IR 2 . Another example of positively Poisson stable vector field is given by the vector field X 0 corresponding to u = 0 in the predator-prey system (7).

Theorem 2.5 Assume that the drift vector field X 0 is positively Poisson stable. Then the system (13) with controls constrained by

u i ∈ {-a i , a i }, a i > 0, i = 1 . . . m is controllable on M if for each point x ∈ M , one has dim (Lie{X 0 , X 1 , . . . , X m }(x)) = n = dimM .
Remark: The control constraint U can actually be any subset of IR m whose convex hull contains the origin in its interior. For example, if m = 2 then the above theorem can be applied with U = {(- 1 2 ,0),(3,0),(0,(0,1 4 )}.

We go back to the harvested prey system (7) which can be written ẋ = X 0 (x) + uX 1 (x), with x = (N, P ) ∈ M = {(N, P ) ∈ IR 2 : N > 0, P > 0}, X 0 = aN -bN P -cP + ebN P , and

X 1 = -qN 0 .
Let us evaluate the rank of the Lie algebra generated by the vector fields X 0 and X 1 : The vectors X 0 (x) and X 1 (x) are linearly independent everywhere (in M ) except at the points of the vertical line N = c eb so the vector space spanned by those two vectors is not of dimension 2 at every point x ∈ M . Hence, we have to compute their Lie bracket

[X 0 , X 1 ](x) = 0 ebN P
. It is clear that the vector space spanned by the three vectors

X 0 (x), X 1 (x) and [X 0 , X 1 ](x) is of dimension 2 at any point x ∈ M . Therefore, we have dim (Lie{X 0 , X 1 }(x))
= 2 = dimM and since the vector field X 0 is positively Poisson stable (its trajectories are closed curves in M ), we can apply Theorem 2.5, that is, the system ( 7) is controllable on M with control constrained by, for instance, u ∈ {-1, 1}. However, we have seen that the system (7) can be controlled with a piecewise control that can take only the values 0 and 1. This fact has been established by using a qualitative geometrical reasoning. In this situation, the control values set is U = {0, 1} and its convex hull is co(U ) = [0, 1] which does not contain zero in its interior so Theorem 2.5 can not apply. Below, we give another result that can be helpful to study controllability even if the system is not control affine nor the origin is an interior point of the convex hull of the control values set U .

Theorem 2.6 Consider the analytic control system (6). Assume that (i) the system is symmetric, i.e., X(x, -u) = -X(x, u) for each x ∈ M and each u ∈ U , or (ii) each vector field of the family F = {X u , u ∈ U } is positively Poisson stable.

Then the system ( 6) is controllable if and only if the accessibility rank condition holds:

rank(Lie(F)(x)) = n = dimM, for all x ∈ M.
Let us consider again the predator-prey system (7):

   Ṅ = aN -bN P -q uN, Ṗ = -cP + ebN P, x = (N, P ), N > 0, P > 0, u ∈ {0, 1}.
(16)

For this system we have U = {0, 1} and

F = {X u , u ∈ U } = {X 0 , X 1 } with X 0 (x) = aN -bN P -cP + ebN P , X 1 (x) = aN -bN P -qN -cP + ebN P .
It is important to note that X 1 is not the controlled vector field X 1 (here, X 0 = X 0 and X 1 = X 0 + X 1 ). The two vector fields X 0 and X 1 are positively Poisson stable on M because, as we have seen before, their trajectories are the closed level curves of, respectively, the realvalued functions

V 0 (N, P ) = e b N -c log( e b c N ) + b P -a log( b a P ) -(a + c) V 1 (N, P ) = e b N -c log( e b c N ) + b P -(a -q) log( b a-q P ) -(a -q + c).
Therefore, by Theorem 2.6, system ( 16) is controllable if and only if rank(Lie{X 0 , X 1 }(x)) = 2 for all x ∈ M . This rank condition holds since

[X 0 , X 1 ](x) = 0 ebqN P and rank{X 0 (x), X 1 (x), [X 0 , X 1 ](x)} = 2 for all x ∈ M .
It is worthwhile noticing that the system (16) remains controllable if U = {u 1 , u 2 } for any non-negative real numbers u 1 = u 2 . This results from the fact that the generated vector fields X u are positively Poisson stable for any value of u and that the generated Lie algebra satisfies the rank condition provided u 1 = u 2 .

Stability

General definitions

Let us consider a system of ordinary differential equations

ẋ = X(x) , x ∈ Ω , ( 17
)
where Ω is an open connected subset of R n and X is a locally Lipschitz continuous map from Ω to R n . For each x ∈ Ω, let us denote by X t (x) the solution of (17) satisfying X 0 (x) = x.

An equilibrium point or a steady state is a state x e ∈ Ω satisfying X(x e ) = 0. Corresponding to each equilibrium point x e , we have a constant solution X t (x e ) ≡ x e of (17).

Let x e ∈ Ω be an equilibrium point. The system ( 17) is stable (we also say Lyapunov stable) at x e or x e is a stable equilibrium position for (17), if for each > 0 there exists a positive real number δ such that for each x with x -x e < δ, the solution X t (x) is defined for all t ≥ 0 and satisfies X t (x) -x e < for all t > 0. When ( 17) is not Lyapunov stable at x e , we say that it is unstable at x e , or that x e is an unstable equilibrium for the system (17).

Lyapunov stability of an equilibrium x e means that all solutions starting at nearby points stay nearby. The Lyapunov stability is an important property. For, let x e be a desired steady-state of our system. Unpredictable perturbations may cause the system to deviate from x e . Lyapunov stability guarantees that every state value taken by the system in its future evolution is not too far from the desired one if the perturbations are small.

The steady state x e is said to be attractive (We also say that ( 17) is attractive at x e ) if there exists a neighborhood U ⊂ Ω of x e such that for any initial condition x belonging to U, the corresponding solution X t (x) of ( 17) is defined for all t ≥ 0 and tends to x e as t tends to infinity, i.e., lim t→+∞ X t (x) = x e .

Let A be the set of points x ∈ Ω such that lim t→+∞

X t (x) = x e
holds for all solutions X t (x) starting from x. Then A is connected and it is called the region of attraction of x e .

The equilibrium x e is globally attractive if A = Ω, i.e., lim t→+∞ X t (x) = x e , ∀x ∈ Ω.

Remark: An equilibrium can be attractive without being stable. This can be illustrated by the following classical two-dimensional system written in polar cordinates (r, θ):

ṙ = r(1 -r) θ = sin 2 θ 2 .
This system has two equilibria: the origin 0 of IR 2 and the point P defined by r = 1 and θ = 0. It can be proved that P is attractive and that its attraction region is A(P ) = IR 2 \{0}. However it is not stable. (see figure 5)

-2 -1 0 1 2 -2 -1 1 2
. P S( x e , ε)

x e x e S( x e , ε)

Figure 6: The phase portrait of a stable and an asymptotically stable equilibrium. S(x e , r) is a sphere centered at x e with radius r.

The steady state x e is said to be marginally stable if it is stable but not asymptotically stable. For instance, the origin is a marginally stable equilibrium point for the following planar system ẋ = -y ẏ = x.

The system (17) is exponentially stable (respectively globally exponentially stable) at x e if there exist two positive constants K an λ such that

X t (x) -x e < K x -x e e -λt
for all x in a neighborhood of x e (respectively for all x ∈ Ω) and all positive time t.

An equilibrium x e is an isolated equilibrium point for the system (17) if it has a neighborhood U such that X(x) = 0 for all x ∈ U. In the sequel, only isolated equilibrium points will be considered and since any isolated equilibrium can be shifted to the origin (of IR n ) by a change of variables, x = x -x e , we shall assume that x e = 0.

Remark:

The above definitions mean that the solution of the differential system (17) has to be explicitly known before stability conclusions can be drawn. Except for a linear system, it is often not possible to compute the analytic expression of the solution of (17). To overcome this difficulty, some tools have been developed to explore the stability properties of nonlinear systems. The most effective are the linearization technique and the Lyapunov's second method. The linearization allows to derive only local conclusions while Lyapunov's second method allows to get global results.

Stability of linear systems

For linear time-invariant system

ẋ = Ax, x ∈ IR n , ( 18
)
the solution is given by X t (x) = e tA x. Therefore, the previous stability conditions can be expressed in terms of the eigenvalues of the matrix A. Let σ(A) = {λ i = α i + jβ i , i = 1 . . . k ≤ n} be the spectrum of A, i.e., the set of all the distinct eigenvalues of A, α i is the real part of the eigenvalue λ i and β i the imaginary part. Let Ψ A be the minimal polynomial of A, that is, the monic polynomial π(λ) of least degree such that π(A) = 0. We have

Ψ A (λ) = (λ -λ 1 ) n 1 (λ -λ 2 ) n 2 . . . (λ -λ k ) n k , with n 1 + n 2 + . . . + n k ≤ n.
The stability properties of the system (18) can then be summarized as follows:

Theorem 3.1 The system ( 18) is asymptotically stable if and only if α i < 0 for all eigenvalues.

In this case, we say that A is a stable matrix or a Hurwitz matrix. If there is an eigenvalue λ i for which α i > 0, then the system (18) is unstable. In this case the matrix A is called unstable.

A matrix which is neither stable nor unstable is called critical. The eigenvalues of a critical matrix are all with non-positive real parts and at least one of them is with zero real part. The eigenvalues with vanishing real parts are called critical eigenvalues or critical characteristic roots.

Theorem 3.2 If the matrix A is critical then the equilibrium of the equation ( 18) is marginally stable (stable but not attractive) if all the critical characteristic roots are simple roots of the minimal polynomial Ψ A . Otherwise the equilibrium is unstable.

For discrete-time system x(k + 1) = Ax(k), we have the following analogous criteria:

• the discrete-time system is unstable if there is an eigenvalue λ i such that |α i | > 1 or if there is an eigenvalue λ i which is not a simple root of the minimal polynomial Ψ A (n i ≥ 2) for which |α i | ≥ 1; • the discrete-time system is Lyapunov stable if |α i | ≤ 1 for all eigenvalues that are simple roots of Ψ A and |α i | < 1 for all repeated eigenvalues; • the discrete-time system is asymptotically stable if |α i | < 1 for all eigenvalues. It must be noted that for linear systems (continuous or discrete-time), we have the following equivalences: asymptotic stability ⇐⇒ global asymptotic stability ⇐⇒ exponential stability.

Linearization and stability of nonlinear systems

Let x = 0 be an equilibrium point for the nonlinear system (17) where the vector function X : Ω -→ IR n is assumed to be continuously differentiable and Ω is a neighborhood of the origin (i.e., the origin is an interior point of Ω). In a coordinate system (x 1 , x 2 , . . . , x n ), the system (17) is given by the following system of differential equations:

     ẋ1 = X 1 (x) = X 1 (x 1 , x 2 , . . . , x n ) . . . ẋn = X n (x) = X n (x 1 , x 2 , . . . , x n ). (19)

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 21

Let A be the Jacobian matrix of the vector function X at the origin, that is,

A = ∂X ∂x (0) =         ∂X 1 ∂x 1 (0) ∂X 1 ∂x 2 (0) . . . . . . ∂X 1 ∂xn (0) ∂X 2 ∂x 1 (0) ∂X 2 ∂x 2 (0) . . . . . . ∂X 2 ∂xn (0) . . . . . . . . . ∂Xn ∂x 1 (0) ∂Xn ∂x 2 (0) . . . . . . ∂Xn ∂xn (0)         (20)
Then we can write

ẋ = X(x) = Ax + Y (x)
where the vector function Y satisfies lim

x →0

Y (x) x = 0. The linear system ẋ = Ax is called the linearized system of the nonlinear system (19) at the origin. This linear approximation is valid for small deviation from the equilibrium point x = 0. Often it is sufficient to draw conclusions about the local stability properties of the nonlinear system thanks to the following theorem known as Lyapunov's first method.

Theorem 3.3 If all the eigenvalues of the Jacobian matrix A have negative real parts, then the origin is an asymptotically stable equilibrium point for the nonlinear system ( 19). If at least one eigenvalue of A has a positive real part, then the origin is unstable for the nonlinear system (19).

The linearization technique or the Lyapunov's first method is a simple technique, and is usually the first method used in the stability analysis of an equilibrium point for a given nonlinear system. However, it does not say how large is the region of attraction of the considered equilibrium when all the eigenvalues of the Jacobian matrix have negative real parts. It also must be noted that the Lyapunov's first method is inconclusive when all the eigenvalues of the Jacobian matrix have non-positive real parts and some of them actually have a zero real part.

Lyapunov functions

The Lyapunov's first method uses the linearized system to reveal the stability properties of the nonlinear system (19). The Lyapunov's second method (or the direct method) works explicitly with the nonlinear system (19). This method can often be used to determine the stability of the equilibrium when the information obtained from the linearization is inconclusive. It also has the advantage of enabling the analysis to extend beyond only a small neighborhood of the equilibrium. The Second Method of Lyapunov is based on the use of auxiliary functions called Lyapunov functions.

Let V be a real-valued function defined and continuous in Ω. We say that the function V is a Lyapunov function for the system (17) on U ⊂ Ω if it is non-increasing along the solutions of the system (17), that is, V (X t (x)) ≤ V (X t (x)) for all x ∈ U and all 0 ≤ t ≤ t . When the function V is of class C 1 then its value never increases along the trajectories of the system if the time derivative of t → V (X t (x)) is non-positive, that is, for all x ∈ U,

V (x) = X.V (x) ≤ 0.
The function X.V is called the Lie derivative of V along the vector function X, it is defined by

X.V (x) = d dt V X t (x) t=0 .
Let (x 1 , . . . , x n ) be a coordinate system. If X(x) = (X 1 (x), . . . , X n (x)) T , . , . is a scalar product and ∇V is the gradient of V in these coordinates, then

X.V (x) = X(x) , ∇V (x) = n i=1 X i (x) ∂V ∂x i (x).
This formula shows that the time derivative of the function V along the solutions of the system ( 17) can be computed without explicitly integrating the differential equation.

A real-valued function V is said to be positive definite at a point x e ∈ Ω if there is a neighborhood U of x e such that V (x) > 0 for each x ∈ U \ {x e }. Similarly, the function V is said to be negative definite at x e on U if V (x) < 0 for all x ∈ U \ {x e }.

We go back to the system (17) and we suppose that the origin is an equilibrium point, i.e., X(0) = 0. The stability properties of the system (17) can be studied with the help of Lyapunov functions thanks to the following result known as Lyapunov's second method. This method can be seen as a generalization of the energy method. It has been inspired by the fact that a stable equilibrium state for a mechanical system corresponds to a local minimum of the total energy.

Theorem 3.4 (Lyapunov 1892) Suppose there exists a C 1 function V defined on some neighborhood U of the origin such that i. The function V is positive definite at the origin, i.e., V (x) > 0 for all x ∈ U \ {0} and V (0) = 0. ii. V (x) ≤ 0 for all x ∈ U. Then the origin is a Lyapunov stable equilibrium point for the system (17). If moreover the function x → V (x) is negative definite, that is V (x) < 0 for all x ∈ U \ {0} and that V (0) = 0, then the origin is an asymptotically stable equilibrium point for the system (17).

There also exists a global version of this theorem. To this end we need to introduce the notion of a proper function. A continuous real valued function V : Ω -→ V (Ω) ⊂ IR is said to be proper if, for all α ∈ V (Ω) ⊂ IR, the set {x ∈ Ω : V (x) ≤ α} is a compact subset of Ω. Proper functions are useful for proving that the solutions of the system (17) are bounded: Suppose that the system (17) has a Lyapunov function V on Ω, that is V (x) ≤ 0, for all x ∈ Ω. If the function V is proper then all the solutions are bounded, i.e., ∀x ∈ Ω, ∃M ≥ 0 such that X t (x) ≤ M, ∀t ≥ 0.

Since V is assumed to be continuous and Ω is a connected open subset of IR n , we have

V (Ω) = (a, b) (or V (Ω) = [a, b)
), b can be a real or +∞. So, another characterization for the function V to be proper is that

lim x→∂Ω V (x) = b, lim x →+∞ x∈Ω V (x) = b For instance, V : x → arctan( x ) is proper on IR n . The function W : (x 1 , x 2 ) → x 2 1 + x 2 2 1 + x 2 2
is not proper on IR 2 because lim

x 2 →+∞ W (0, x 2 ) = 1.
Now, we give the global version of the Lyapunov's theorem:

Theorem 3.5 If there exists a C 1 function V defined on the whole state space Ω and satisfying i. The function V is positive definite at the origin, i.e, V (x) > 0 for all x ∈ Ω \ {0} and

V (0) = 0. ii. V (x) < 0 for all x ∈ Ω \ {0}.
iii. The function V is proper. Then the origin is a globally asymptotically stable equilibrium for the system (17).

Example: The Predator-Prey Model

Let us consider again the simple Lotka-Volterra predator-prey model ( 7) with u = 0:

Ṅ = aN -bN P Ṗ = -cP + ebN P (21)
All the variables and constants are positive. This system has two equilibria: the trivial equilibrium (0, 0) and the other equilibrium (N * , P

* ) = ( c eb , a b
).

The linearization of the system (21) around the point (0, 0) leads to Ṅ Ṗ = a 0 0 -c A 0

N P

The matrix A 0 has a positive eigenvalue which implies that the point (0, 0) is an unstable equilibrium point.

The linearization of the system (21) around the point (N * , P * ) leads to

Ṅ Ṗ = 0 -c/e ea 0 A 1 N P .
The matrix A 1 has two eigenvalues whose real parts are zero, so the linearization method does not give any information about the stability of the nonlinear system (21).

In this situation, the Lyapunov's theorem can be helpful as we shall see. Let Ω be the interior of the positive quadrant, i.e., Ω = IR 2 +> = {(N, P ) ∈ IR 2 : N > 0, P > 0}. Let V be the function defined on Ω by

V (N, P ) = e b N -c log( e b c N ) + b P -a log( b a P ) -(a + c).
This function is positive definite at (N * , P * ): V (N * , P * ) = 0 and V (N, P ) > 0 for all N = N * and P = P * . The function V is proper on Ω because we have

V (Ω) = [0, +∞), lim (N,P ) →+∞ V (N, P ) = +∞, lim N →0
V (N, P ) = +∞ and lim P →0

V (N, P ) = +∞.

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 24

Moreover its time derivative is V = 0. Therefore, by Lyapunov's theorem, the nontrivial equilibrium (N * , P * ) is stable and, since V is proper on Ω, all the solutions are bounded. Actually, the above function V , known as the ecological Lyapunov function, is a constant of motion. Hence, the solutions of (21) lie on fixed curves defined by V = k. For a given initial population distribution (N 0 , P 0 ), the corresponding solution (N (t), P (t)) of ( 21) satisfies V (N (t), P (t)) = V (N 0 , P 0 ), for all positive time t. Now, if we add to the previous predator-prey system terms representing decay due to crowding, we get

Ṅ = aN -bN P -f N 2 Ṗ = -cP + ebN P -gP 2 (22)
The nontrivial equilibrium for this system is N * = ag + bc eb 2 + f g and P * = abe -cf eb 2 + f g . This equilibrium belongs to Ω if and only if abe -cf > 0.

A candidate Lyapunov function is

V (N, P ) = eN * N N * -log N N * + P * P P * -log P P * -eN * -P * .
This function is positive definite at (N * , P * ) and it is proper on Ω. Let us evaluate its time Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 25 derivative along the solutions of the system (22):

V (N, P ) = e 1 - N * N aN -bN P -f N 2 + 1 - P * P -cP + ebN P -gP 2 = e (N -N * )   a -bP -f N -(a -bP * -f N * ) =0   + (P -P * )   -c + ebN -gP -(-c + ebN * -gP * ) =0   = e (N -N * ) (-f (N -N * ) -b(P -P * )) + (P -P * )(eb(N -N * ) -g(P -P * )) = -ef (N -N * ) 2 -g(P -P * ) 2 .
Hence V (N, P ) < 0 for all (N, P ) ∈ Ω \ {(N * , P * )} and V (N * , P * ) = 0. Thus, by the global version of Lyapunov's theorem, we conclude that the nontrivial equilibrium (N * , P * ) is globally asymptotically stable.

In order to use the original Lyapunov's theorems for the purpose of proving the asymptotic stability of a given system, we have to find a positive definite function V whose time-derivative V is negative definite. This is a difficult task in general. The definiteness of the derivative V can actually be relaxed by using the LaSalle's Invariance Principle that we expose hereafter.

A set G ⊂ Ω is an invariant set for the system (17) if whenever a solution belongs to G at some time, then it belongs to G for all future and past time. That is, if x ∈ G then X t (x) ∈ G, for all t ∈ IR. For example, the following sets are invariant sets for the predator-prey system ( 21 A set G is positively invariant (negatively invariant) if:

x ∈ G =⇒ X t (x) ∈ G, ∀t ≥ 0 (∀t ≤ 0).
The LaSalle's Invariance Principle can be stated as follows Theorem 3.6 Let Ω be a subset of IR n . Assume that Ω is positively invariant for system (17). Let V : Ω → IR be a C 1 scalar function such that V (x) ≤ 0 in Ω. Let E be the set of points within Ω where V (x) = 0, and let L be the largest invariant set within E. Then every bounded solution starting in Ω tends to the set L as time goes to infinity.

This theorem is a very useful tool for system analysis. Unlike Lyapunov's theorem, it does require neither the function V to be positive definite nor the function V to be negative definite.

However, it gives only information about the attraction. Therefore it can be used to show that the solutions tend to an equilibrium if the set L is reduced to this equilibrium but it does not allow to say if this equilibrium is stable or not. Actually, when we are interested in establishing asymptotic stability of an equilibrium point (assumed to be the origin of IR n ), we use the following corollary of the LaSalle's Invariance Principle:

Corollary 3.1 Let Ω be an open connected subset of IR n such that x = 0 ∈ Ω, and x = 0 is an equilibrium state for the system (17). Let U be a neighborhood of the origin in Ω and let V : U → IR be a C 1 positive definite function, such that V (x) ≤ 0 in U. Let E = {x ∈ U : V (x) = 0}, and assume that largest invariant set within E is reduced to the origin. Then, the origin is asymptotically stable.

If the above conditions hold for U = Ω and the function V is proper on Ω, then, the origin is a globally asymptotically stable equilibrium point for the system (17).

This result contains the original Lyapunov's theorem as a special case. It must be noted that the set Ω needs not to be bounded.

Example: Consider again the predator-prey system (21) and suppose that the crowding affects only the prey growth:

Ṅ = aN -bN P -f N 2 , Ṗ = -cP + ebN P. ( 23
)
This system has two equilibria: the trivial equilibrium N = P = 0 (which is unstable) and a nontrivial equilibrium N * = c eb and P * = abe -cf eb 2 . This equilibrium belongs to Ω (The interior of the positive quadrant) if and only if abe -cf > 0. We take again the following candidate Lyapunov function Its time-derivative along the solutions of ( 23) is

V (N, P ) = -ef (N -N * ) 2 ≤ 0.
In this case V is not negative definite so Lyapunov's theorem does not apply. Hence, we have to apply the LaSalle's Invariance principle. Here, E = {(N, P ) ∈ Ω : V (N, P ) = 0} = {(N * , P ), P > 0}. On E, the vector field is Ṅ Ṗ = (a -bP -f N * )N * 0 (cf Figure 8 ).

It is then easy to see that the only invariant set within E is the equilibrium point (N * , P * ) and this proves that this equilibrium is globally asymptotically stable.

When the vector field X and the Lyapunov function are analytic, the set L can be computed according to the following formula:

L = {x ∈ Ω : X k .V (x) = 0, k = 1, 2, . . .}.
Hence, we just have to solve a system of equations. For the previous system (23), we have, with x = (N, P ) and X = (aN -bN P -f N 2 , -cP

+ ebN P ) T , X.V (x) = 0 =⇒ N -N * = 0 (i) X 2 .V (x) = -2ef (N -N * )(aN -bN P -f N 2 ) = 0 (ii) X 3 .V (x) = -2ef (aN -bN P -f N 2 ) -2ef (N -N * )(-bN )(-cP + ebN P ) +(a -bP -2f N )(-2ef (N -N * )(aN -bN P -f N 2 ) = 0 (iii)
By combining (i) and (iii), we get N = N * and P = P * .

Limit cycle

Thus far the only type of attractors encountered have been equilibrium points. For nonlinear systems, an other type of attractor is possible, the limit cycle.

Here we consider two-dimensional systems of the form

   ẋ1 = X 1 (x 1 , x 2 ), ẋ2 = X 2 (x 1 , x 2 ), x = (x 1 , x 2 ) ∈ Ω ⊂ IR 2 , X(x) = (X 1 (x), X 2 (x)), (24)
where X 1 and X 2 are continuously differentiable.

A solution of (24) through the point x ∈ Ω is said to be periodic if there exists T > 0 such that

X t+T (x) = X t (x), ∀t ∈ IR. ( 25
)
The trajectory corresponding to a periodic solution is called a periodic orbit or a closed orbit.

The period of a periodic solution is the smallest T > 0 such that (25) holds. A constant solution (equilibrium position) is a trivial periodic solution. Here we are interested in nontrivial periodic solutions. A system oscillates when it has a (nontrivial) periodic solution.

Consider the linear two-dimensional system ẋ1 = -x 2 , ẋ2 = x 1 . All the solutions are periodic and the corresponding orbits are circles centered at the origin. One can remark that there is a continuum of periodic orbits and that the amplitude of the oscillations depend on the initial condition. The above system is an harmonic oscillator.

Now consider the following nonlinear two-dimensional system

ẋ1 = -x 2 + x 1 (1 -x 2 1 -x 2 2 ), ẋ2 = x 1 + x 2 (1 -x 2 1 -x 2 2 ). ( 26
)
The origin is an equilibrium point and there is no other equilibrium point. The linearization of the system (26) at the origin is given by the matrix

A = 1 -1 1 1
This matrix has two eigenvalues whose real parts are positive. Therefore the origin is an unstable equilibrium point for (26).

The unit circle is an invariant set of the system (26). This fact can be easily seen by representing the system in the polar coordinates (r, θ) defined by:

x 1 = r cos θ, x 2 = r sin θ, which yields ṙ = r(1 -r 2 ), θ = 1.
This shows that the sets defined respectively by r = 0 and r = 1 are invariant sets of the system. The first set r = 0 is just the origin and the second r = 1 corresponds to the unit circle. The polar form also shows that the unit circle is the unique periodic orbit for the system (26). One also can remark that if r < 1 then ṙ > 0 and if r > 1 then ṙ < 0. Therefore all the trajectories (except the trivial solution X t (0) ≡ 0) spiral toward the unit circle from inside or outside (see Figure 11). The unit circle is a stable limit cycle.

Remark: It is possible to show that the unit circle attract all the solutions (except the trivial one) of the system (26) by considering the real-valued function

V (x 1 , x 2 ) = (1 -x 2 1 -x 2
2 ) 2 and applying the LaSalle's Invariance Principle (Theorem 3.6).

A closed orbit γ is called a limit cycle if it is an isolated closed orbit, that is, there exists a neighborhood of γ which contains no other closed orbits of the system (24). A limit cycle is a special periodic solution of the planar system ( 24) that attract all other nearby solutions as t → ∞ or t → -∞. Geometrically this means that the nearby non-closed trajectories spiral toward it, either from inside or outside as t → ∞ or as t → -∞.

-2

-1 1 2 -2 -1 1 2 x 1
x 2

Figure 11:

The trajectories of (26) corresponding to different initial conditions: the unit circle is a limit cycle.

A limit cycle is called stable when nearby trajectories approach it as time t approaches infinity. It is unstable if the all the trajectories starting arbitrarily close to it will tend away from it as t → ∞.

The unit circle is a stable limit cycle for the system (26). It is an unstable limit cycle for

ẋ1 = x 2 -x 1 (1 -x 2 1 -x 2 2 ), ẋ2 = -x 1 -x 2 (1 -x 2 1 -x 2 2 ). ( 27
)
This system is the same as (26) in reversed time. Therefore, its phase portrait is the same as the one of (26) but the arrow heads are reversed.

We present a celebrated result known as Poincaré-Bendixson theorem which is very useful for establishing the existence of periodic orbits, and for exploring the stability properties of nonlinear two dimensional systems. It must be noted however that it has no generalization to higher dimensional systems.

Theorem 3.7 (Poincaré-Bendixson) Let D be a positively invariant compact set for the planar system (24) containing a finite number of equilibria. Let x be a point of D, and consider the corresponding solution X t (x) of ( 24). Then one of the following possibilities holds.

i. The solution X t (x) is a periodic solution;

ii. the solution X t (x) approaches a periodic solution as t → ∞;

iii. the solution X t (x) approaches a cycle graph as t → ∞.

A cycle graph is the union of equilibria and trajectories X t (y) such that lim t→-∞ X t (y) and lim t→∞ X t (y) are equilibria. A cycle graph can be seen as the union of equilibria and orbits connecting them.

The Poincaré-Bendixson theorem says that a bounded trajectory that does not approach any equilibrium point is either a closed periodic orbit or approaches a closed periodic orbit as t → ∞.

An implication of this theorem is that a nonempty compact set that is positively or negatively invariant contains either an equilibrium point or a limit cycle.

The following criterion is practical in ruling out the presence of closed orbits in a region of the plane.

Theorem 3.8 (Bendixson's negative criterion) If on a simply connected domain D ⊂ Ω (that is D is a connected set without holes) the expression

∂X 1 ∂x 1 + ∂X 2 ∂x 2
is not identically zero and does not change sign, then the system (24) has no closed orbits lying entirely in D.

A consequence of the above criterion is that limit cycles can only be obtained with nonlinear systems.

Another criterion that can be used to preclude the existence of periodic solutions in a region of the plane when the Bendixsons fails is the following result due to Dulac.

Theorem 3.9 (Dulac's negative criterion) Suppose there exists a real-valued function ρ, continuously differentiable in a simply connected domain D ⊂ Ω, such that the expression

∂(ρX 1 ) ∂x 1 + ∂(ρX 2 ) ∂x 2
is not identically zero and does not change sign, then the system (24) has no closed orbits lying entirely in D.

Example: Consider the following simple epidemic model which is a simplified version of the system (2) introduced in Section 1 without the class E and without a disease-related death (d = 0).

               Ṡ = bN -µS -β SI N , İ = β SI N -(r + µ)I, Ṫ = rI -µT, Ṅ = (b -µ)N.
(28)

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 32

The proportions s = S/N , i = I/N , and τ = T /N satisfy the following system of differential equations

     ṡ = b -bs -βsi, i = βsi -(r + b)i, τ = ri -bτ. (29)
Since s + i + τ = 1 it is sufficient to study the planar system

ṡ = b -bs -βsi = X 1 (s, i), i = βsi -(r + b)i = X 2 (s, i). ( 30
)
The system evolves in the set

D = {(s, i) ∈ IR 2 : 0 ≤ s ≤ 1, 0 ≤ i ≤ 1, 0 ≤ s + i ≤ 1}
and has two equilibria:

The disease-free equilibrium:

s * = 1, i * = 0.
The endemic equilibrium:

s * = b + r β = 1 R 0 , i * = b (β -(b + r)) (b + r) β
. This equilibrium belongs to

D only if b + r < β.
The linearization technique allows to get the following local stability properties: The diseasefree equilibrium (s * = 1, i * = 0) is asymptotically stable if b + r > β and it is unstable if b + r < β. The endemic equilibrium is asymptotically stable when it exists, that is, if b + r < β.

To rule out the existence of periodic orbits we first apply the Bendixson's criterion on D which is simply connected:

∂X 1 ∂s + ∂X 2 ∂i = -b -βi + βs -(r + b) ≤ -b -βi + β -(r + b) since s ≤ 1.
Therefore the planar system (30) has no periodic solutions if b + r > β. When b + r ≤ β, Bendixson's criterion does not allow to conclude. We apply the Dulac's criterion with the function ρ(s, i)

= 1 si which is continuously differentiable in D 1 = {(s, i) ∈ D : s > 0, i > 0}: ∂(ρX 1 ) ∂s + ∂(ρX 2 ) ∂i = - b s 2 i < 0, ∀(s, i) ∈ D 1 .
Thus, for any values of the positive parameters, the system has no periodic orbits lying entirely in D 1 . The set D \ D 1 can not contain a periodic orbit since it is just the union of two segments. Therefore there are no periodic orbits in D.

Now the Poincaré-Bendixson theorem allows to get the following global stability properties:

• For any values of the parameters, the set D 0 = {(s, i) ∈ D : i = 0} is positively invariant and the disease-free equilibrium attracts all the trajectories emanating from this set. • If b + r ≥ β then the disease-free equilibrium (s * = 1, i * = 0) is globally attractive in D:

the disease dies out. • If b + r < β then the endemic equilibrium is globally asymptotically stable in the region D \ D 0 = {(s, i) ∈ D : i > 0}: the disease will spread.

possible to stabilize the system (31) by a continuous feedback law u(x); however, it is possible to stabilize it by a time-varying feedback u(x, t) if some technical conditions hold. Moreover, for a given positive number T > 0, it is possible to choose a periodic stabilizer u(x, t) that satisfies: u(x, t + T ) = u(x), for all t ≥ 0.

Examples: We present some examples to illustrate the above obstructions to stabilizability:

1. The first example (due to R. Brockett) concerns a system that satisfies the conditions (i) and (iii) but does not satisfy the second condition (ii):

ẋ = X(x, u), with (x 1 , x 2 , x 3 ) ∈ IR 3 , (u 1 , u 2 ) ∈ IR 2 , and X(x, u) =   u 1 u 2 x 2 u 1 -x 1 u 2   (33)
On the one hand, this driftless three-dimensional system is controllable (see Theorem 2.4).

Hence the condition (i) is satisfied. On the other hand, the condition (iii) is obviously fulfilled since the linearized system at the origin

   ẋ1 = u 1 , ẋ2 = u 2 , ẋ3 = 0 ,
has two controllable modes and one uncontrollable mode associated with a vanishing eigenvalue. However the condition (ii) is not satisfied since it is not possible to find (x, u) ∈ IR 3 × IR 2 such that X(x, u) = (0, 0, ) with = 0. Therefore, the system (33) can not be stabilized at the origin by means of a continuous feedback law u(x).

Consider the system

   ẋ1 = x 1 (x 2 1 + x 2 2 ) , ẋ2 = u , (x 1 , x 2 ) ∈ IR 2 , u ∈ IR . (34)
The conditions (ii) and (iii) are not violated but the condition (i) is not satisfied because for any α > 0, the origin can not be asymptotically reached from any point of the open half space x 1 > α. Hence there exists no continuous feedback law u(x 1 , x 2 ) which makes the origin asymptotically stable for the system (34).

The two-dimensional system

   ẋ1 = x 1 -x 3 2 , ẋ2 = u , (x 1 , x 2 ) ∈ IR 2 , u ∈ IR , (35)
satisfies the conditions (i) and (ii) but does not satisfy the condition (iii) because the linearized system at the origin ẋ1 = x 1 , ẋ2 = u , has an uncontrollable mode associated with a positive eigenvalue. Therefore, the system (35) can not be stabilized at the origin by means of a C 1 feedback law u(x 1 , x 2 ).

Sufficient stabilizability conditions

For linear systems, the stabilization problem has been completely solved and the construction of static feedback stabilizers can be done in a systematic way. For nonlinear systems, the problem is a hard task and still under intensive investigation. No general method is available but many powerful techniques have been developed, each one allows to compute the stabilizers for a class of nonlinear systems. Here, we shall expose just one of these techniques because it is one of the most general and very simple to apply. We consider an affine control system defined by smooth (C ∞ ) vector fields X, Y 1 , . . . , Y m , with an equilibrium point at the origin of IR n :

   ẋ = X(x) + m i=1 u i Y i (x). X(0) = 0, Y i (0) = 0. (36)
If there exists a positive definite and proper smooth function V : IR n → IR such that: (i) The Lie-derivative of V with respect to the vector field X satisfies

X •V (x) ≤ 0, ∀x ∈ IR n , (ii) The set W = {x ∈ IR n | X k+1 •V (x) = X k •Y i •V (x) = 0, k ∈ IN, i = 1, . . . , m} is reduced
to the set {0}. Then the affine system (36) is globally stabilizable by the smooth state feedback control law

u(x) = -    Y 1 •V (x)
. . .

Y m •V (x)    . (37)
Remark: The above feedback is known as the Jurdjevic-Quinn feedback. It must be noted that if the assumptions (i) and (ii) are satisfied then, for any smooth real-valued function β satisfying β(x) > 0 for all x ∈ IR n \ {0}, the feedback law

u(x) = -β(x)    Y 1 •V (x) . . . Y m •V (x)    (38)
also stabilizes the system (36) at the origin because the time derivative of the function V along the solutions of the closed-loop system (36-38) is

V (x) = X •V (x) -β(x) m i=1 (Y i •V (x)) 2 ≤ 0
and the assumption (ii) implies that the largest invariant set contained in the set { V = 0} is reduced to the origin. Hence the LaSalle's Invariance Principle allows to conclude that the origin is a globally asymptotically stable equilibrium for the closed loop-system (36-38).

Example: A harvested fish population. Consider the following stage-structured model of a fish population that has been built by Touzeau-Gouzé. It includes (n + 1) stages represented by their abundance x i (t), stage 0 being the pre-recruits stage. Each stage is characterized by its fecundity, mortality and predation rates:

         ẋ0 (t) = -αx 0 (t) -m 0 x 0 (t) + n i=1 f i l i x i (t) -n i=0 p i x i (t)x 0 (t), ẋ1 (t) = αx 0 (t) -αx 1 (t) -m 1 x 1 (t), . . . ẋn (t) = αx n-1 (t) -αx n (t) -m n x n (t), (39)
where • m i : linear mortality rate.

• α: linear aging coefficient.

• p 0 : juvenile competition parameter.

• p i : predation rate of class i on class 0.

• f i : fecundity rate of class i.

• l i : reproduction efficiency of class i.

The mortality coefficient can be written as a sum of the natural mortality rate M i and the fishing mortality coefficient F i . Hence one can write: m i = M i + q i E, where q i is the catchability of stage i and E is the fishing effort that can be seen as a control term. Denoting

α i = α + M i , we get          ẋ0 (t) = -α 0 x 0 (t) + n i=1 f i l i x i (t) -n i=0 p i x i (t)x 0 (t) ẋ1 (t) = αx 0 (t) -(α 1 + q 1 E(t))x 1 (t) . . . ẋn (t) = αx n-1 (t) -(α n + q n E(t))x n (t) (40)
The origin is an equilibrium point which corresponds to an extinct population and is therefore not very interesting. Under some conditions and for a constant fishing effort Ē, there exists another nontrivial equilibrium x * whose coordinates are

             x * 0 = n i=1 f i l i π i -α 0 p 0 + n i=1 p i π i x * i = π i x * 0 π i = α i i j=1 (α i + q i Ē) (41)
The goal is to compute the fishing effort as feedback control E(x) = Ē + u(x), in order to maintain the fish population around its steady state x * . The state x * becomes a globally asymptotically stable equilibrium for the closed-loop system within Ω, the positive quadrant.

The system can be rewritten as follows

ẋ =      -α 0 x 0 + n i=1 f i l i x i -n i=0 p i x i x 0 αx 0 -(α 1 + q 1 Ē)x 1 . . . αx n-1 -(α n + q n Ē)x n      X(x) +u      0 -q 1 x 1 . . . -q n x n      Y (x) . ( 42
)
Let V be the following candidate Lyapunov function

V (x) = 1 2 (x 0 -x * 0 ) 2 + n i=1 n j=i k j π j α i + q i Ē x i -x * i π i 2 , k i = f i l i -p i x * 0 i = 1, . . . , n.
V is a positive definite function on Ω provided that

x * 0 < min i=1,...,n n j=i f j l j π j n j=i p j π j
Its derivative along the drift vector field X satisfies

X.V (x) = X(x), ∇V (x) ≤ - 1 2 n i=1 k i π i (x 0 -x * 0 ) -( x i -x * i π i ) 2 ≤ 0
A candidate stabilizer given by ( 37)

is u = Φ(x) = -Y (x), ∇V (x) = n i=1 γ i q i x i (x i -x * i ), with γ i = n j=i k j π j (α i + q i Ē)π 2 i
. However, the function Φ is unbounded and takes positive as well as negative values, hence E(x) = Ē + u(x) can take negative values which is not possible in practice (E is a fishing effort). Thus, instead of using the feedback given by the formula (37), we use the feedback law given by the formula (38) u(x) = β(x)Φ(x) and we choose the function β such that β(x)Φ(x) ≥ -Ē. A simple computation shows that Φ(x) ≥ -n i=1 γ i q i

x * i 4 . Therefore we take β

(x) = 4 Ē n i=1 γ i q i x * i and u(x) = 4 Ē n i=1 γ i q i x * i Φ(x) . (43)
We then have u(x) ≥ -Ē which ensures that E(x) ≥ 0 for all x ∈ Ω. The figure 12 presents a simulation for a five stages system controlled by the feedback law given by the formula (43). x 3

Stage i i = 0 i = 1 i = 2 i = 3 i = 4 p i 0.2 0 0.1 0.1 0.1 f i 0 0.5 0.5 0.5 l i 0 10 20 15
q i 0 0 0 0.1 0.15 α 0.8 α i 1.3 1 1 0.9 0.85 Ē 1
x 4 E

x 0 * x 1 * x 3 * x 2 *
x *

E -

Figure 12: The time evolution of the different stages for the closed-loop system (42-43) with the parameter values given in Table 1 and with initial conditions:

x 0 (0) = 20, x 1 (0) = 20, x 2 (0) = 45, x 3 (0) = 10 and x 4 (0) = 3.
where x ∈ M = IR, u ∈ IR, y ∈ IR. Every pair of distinct states x 1 and x 2 are distinguishable: it is sufficient to take u = c = 0. However, the input u(t) ≡ 0 is not an universal input: it does not distinguish the states x and -x.

The system (44) is observable if any pair of distinct initial states (x 1 , x 2 ) are distinguishable.

The system (44) is uniformly input observable (or observable for any input) if for any input u(.) ∈ U ad and for any pair of distinct initial states (x 1 , x 2 ), there exists a time t ≥ 0 such that y(x 1 , u(.), t) = y(x 2 , u(.), t).

Example: Bacterial growth in a chemostat.

Consider the following model of the chemostat:

       ẋ = µ(s) x -ux ṡ = -kµ(s) x -u(s -s in ) y = x (45)
where x(t) and s(t) are respectively the concentration in micro-organisms and substrate, the function µ(s) is the absorbing rate of the substrate by the micro-organisms and k is a constant. For this system, the input is the flow rate u and the output is usually the concentration x(t). If the map s → µ(s) is injective then the chemostat ( 45) is uniformly input observable.

Observability also means that the data of the output y(t) and the input u(t) on any finite time interval [t 0 , t 1 ] allow to recover the initial state x 0 and therefore the trajectory starting from this initial state.

Observability of linear systems

In this section, we study the observability properties of linear systems

             ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t), x(t) ∈ IR n , u(t) ∈ U ⊂ IR m , y(t) ∈ IR q ,
A, B, and C are respectively n × n, n × p and m × n matrices.

(

) 46
For the linear time-invariant system (46), we have y (x 1 , u(.), t) = Ce tA x 1 + C t 0 e (t-s)A Bu(s)ds .

Hence, y (x 1 , u(.), t) -y (x 2 , u(.), t) = Ce tA (x 1 -x 2 ).

Thus, x 1 and x 2 are indistinguishable if and only if Ce tA (x 1 -x 2 ) = 0 for all t ≥ 0. By analyticity of Ce tA x, this is equivalent to say that all the derivatives of t → Ce tA (x 1 -x 2 ) vanish at t = 0, i.e,

d k dt k Ce tA (x 1 -x 2 ) t=0 = 0, for all k = 0, 1, 2, . . . ,
Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 41 which can be written :

CA k (x 1 -x 2 ) = 0, for all k = 0, 1, 2, . . . .

Thanks to Cayley-Hamilton theorem, this is equivalent to :

C(x 1 -x 2 ) = CA(x 1 -x 2 ) = . . . = CA n-1 (x 1 -x 2 ) = 0.
Therefore x 1 and x 2 are indistinguishable if and only if

x 1 -x 2 ∈ ker C ∩ ker CA ∩ ker CA 2 ∩ . . . ∩ ker CA n-1
We remark that this condition is independent of the input. So the linear system ( 46) is uniformly input observable if and only if it is observable if and only if the matrix:

O (C,A) =         C CA CA 2 . . . CA n-1        
is of rank n. In this case we say that the system (46), or the pair (C, A), satisfies the Kalman rank condition for observability.

Another interpretation of observability for linear systems: Let us consider now the system (46) with a single output, that is, y(t) ∈ IR (q = 1). If the system considered is observable (we also say the pair (C, A) is observable) then the linear change of coordinates z = O (C,A) x allows us to write system (46) in the following observability canonical form : ż(t) = Ãz(t) + Bu(t),

y(t) = Cz(t), (47) with
: Ã =        0 1 0 0 . . . 0 0 0 1 0 . . . 0 . . . . . . 0 0 0 0 . . . 1 -a 0 -a 1 -a 2 -a 3 . . . -a n-1        and C = 1 0 0 . . . 0 .
Let us compute the successive time derivatives of the output:

                   y(t) = Cz(t) ẏ(t) = C ż(t) = C Ãz(t) + C Bu(t) ÿ(t) = C Ã2 z(t) + C Ã Bu(t) + C B u(t) . . . y (n-1) (t) = C Ãn-1 z(t) + C Ãn-2 Bu(t) + C Ãn-3 B u(t) + . . . + C Bu (n-2) (t)

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 42

This can be written in a condensed form :

O ( C, Ã) z(t) =           y(t) ẏ(t) -C Bu(t) ÿ(t) -C Ã Bu(t) -C B u(t) . . . y (n-1) (t) -C Ãn-2 Bu(t) -C Ãn-3 B u(t) -. . . -C Bu (n-2) (t)           .
From this relation we see that, if the functions t → y(t) and t → u(t) are known then one can compute the state vector z(t) uniquely since the matrix O ( C, Ã) is invertible.

Observability of nonlinear systems

We go back to the nonlinear system (44). To simplify the notations, we shall consider systems with single output, i.e., q = 1 and y(t) = h(x(t)) ∈ IR.

In order to derive an observability condition for (44), we need to recall that a smooth (C ∞ ) vector field X defined on M operates on C ∞ (M ), the set of C ∞ functions Φ : M -→ IR, by Lie differentiation in the following way:

C ∞ (M ) -→ C ∞ (M ) Φ -→ X.Φ with X.Φ (x) = d dt Φ X t (x) t=0
In the coordinate system (x 1 , . . . , x n ), let . , . denote a scalar product and ∇Φ the gradient of Φ in these coordinates. If

X(x) = n i=1 X i ∂/∂x i (x) =       X 1 (x 1 , . . . , x n ) X 2 (x 1 , . . . , x n ) . . . X n (x 1 , . . . , x n )       , then X.Φ (x) = ∇Φ(x) , X(x) = n i=1 X i (x) ∂Φ ∂x i (x)
The function X.Φ is called the Lie derivative of Φ along the vector field X. It is also denoted L X Φ. For a given positive integer k > 0, the Lie derivative of order k of Φ along X is defined by induction as follows:

X k .Φ = X.(X k-1 .Φ) .
Example: ẋ = X(x) = Ax, y = h(x) = Cx. Here, we have X t (x 0 ) = e tA x 0 . Hence X.h(x 0 ) = d dt Ce tA x 0 t=0

= CAx 0 and it is easy to see that X k .h(x 0 ) = CA k x 0 .

The observation space of (44) O is the linear space (over IR) of functions on M containing the observation function h and which is closed under Lie differentiation by all elements of F = {X u , u ∈ U }. (F is just the set of vector fields corresponding to constant controls).

It can be proved that O can be defined as the set of all the linear combination of all repeated Lie derivatives of functions of the form X u k . . . X u 2 X u 1 .h. i.e.:

O = span IR {(X u l ) k l . . . (X u 2 ) k 2 .(X u 1 ) k 1 .h : l ≥ 0, u 1 , . . . , u l ∈ U, k i = 0, 1, 2, . . .}.
For analytic systems, the observability is equivalent to the fact that the observability space O separates the points of M .

Remark:

The observation space O contains the output function and all derivatives of the output function along the system trajectories. In particular, for a system without input

ẋ(t) = X(x(t)), y(t) = h(x(t)),
O is constructed by taking y = h(x) together with all repeated time derivatives ẏ = X.h(x), ÿ = X 2 .h(x), . . .. For the linear system (46), the observation space is generated by the functions Cx, CAx, . . . , CA n-1 x.

Remark: We have seen that the observability of linear systems does not depend on the input. This is no more true for nonlinear systems as it can be shown by the following example:

   ẋ = 1 0 1 -1 x + u 0 1 1 0 x, y = x 1 . (48)
Here

x = (x 1 , x 2 ) ∈ M = IR 2 , u ∈ U = {0, 1}
. This is a very simple nonlinear system : it is a bilinear system. This system is observable, because with the input u(t) ≡ 1, we obtain an observable linear system ẋ = Ax , y = Cx, with A = 1 1 2 -1 and C = 1 0 . It is easy to see that the Kalman rank condition for observability is satisfied. One also can remark that it is possible to reconstruct the state (x 1 (t), x 2 (t)) from y(t) = x 1 (t) and ẏ(t) = x 1 (t)+x 2 (t). Hence, the input u(t) ≡ 1 distinguish every pair of distinct initial states x 0 and x0 and this proves that the system (48) is observable. However this system is not uniformly input observable because the input u(t) ≡ 0 does not distinguish states x 0 and x0 satisfying: A common form for the observers is the following ẋ(t) = g(x(t), y(t), u(t)),

where g is smooth. In practice, an on-line estimation of the state is obtained by choosing an initial condition x(0) and integrating equation ( 53) on a computer. We must notice here that it is not possible to do the same thing with the system (51) (which is supposed to model a real system) because we do not know the initial condition x(0). Let e(t) = x(t) -x(t) be the estimation error. The error equation is then given by ė = f (x, u) -g(x, y, u).

(54)

If the system ( 53) is an observer for the system (51), then e = 0 must be an equilibrium for (54) which implies that g(x, h(x), u) = f (x, u) for all x and all admissible inputs u. Furthermore the null solution e(t) ≡ 0 of (54) has to be globally asymptotically stable. If it is only locally asymptotically stable then the dynamical system (53) is a local observer for the system (51). System ( 53) is said to be an exponential observer if the estimation error decreases at an exponential rate, that is, there exist positive constants c and a in such a way that the solutions x(t) and x(t) of ( 51) and ( 53) satisfy for any initial conditions x(0) and x(0)

x(t) -x(t) ≤ c x(0) -x(0) e -a t , ∀t > 0.

To construct an observer for a given system, one has to find the "good" function g that satisfy the above conditions. For linear systems this problem has been solved by the Luenberger observer. For nonlinear systems, there is no "universal" solution but several methods have been developed for some classes of systems. Usually, the simplest and the most natural way to construct an observer for the system (51) is to take a copy of it and to add a corrective term that depends on the difference h(x(t)) -y(t), for instance g(x, y, u) = f (x, u) + K(x) (h(x) -y) .

Observers for linear systems

We consider an observable linear system

             ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t), x(t) ∈ IR n , u(t) ∈ U ⊂ IR m , y(t) ∈ IR q ,
A, B, and C are respectively n × n, n × m and q × n matrices.

An observer (called Luenberger Observer) for this system is

ẋ(t) = Ax(t) + Bu(t) + K y(t) -C x(t) (56)
where the n × q matrix K is to be computed. The observer is an n-dimensional system with state vector x(t). The inputs of the observer consist of a copy of the inputs of the original system (55) and the measurements y(t) available from the original system (55). The estimation error e(t) = x(t) -x(t) is governed by the differential equation

ė(t) = Ae(t) -K y(t) -C x(t) = Ae(t) -K Cx(t) -C x(t) = (A -KC)e(t) .
System ( 55) is assumed to be observable, so the pair (A, C) is observable and this is equivalent to say that the pair (A T , C T ) ( T = transpose) is controllable. For a pair of controllable matrices, we can apply the Pole-Shifting Theorem which says that given any nth-order real polynomial p(λ) = λ n + a n-1 λ n-1 + . . . + a 0 there is a real matrix F such that the matrix A T + C T F has p(λ) as its characteristic polynomial. In other words, for any given set S = {α 1 , . . . α n } of n complex numbers satisfying z ∈ S ⇒ z ∈ S, it is possible to find a matrix F in such a way that the spectrum of

A T + C T F is σ(A T + C T F ) = S.
Since the spectrum of a real matrix M and the one of its transpose M T are equal, we have

σ(A T + C T F ) = σ (A T + C T F ) T = σ(A + F T C).
In particular, there exists a matrix F such that all the eigenvalues of (A + F T C) are with negative real part. Therefore, if we take K = -F T then the estimation error satisfies e(t) ≤ c e -αt , where c > 0, α > 0, and α = max

λ∈σ(A-KC) |Re(λ)|.
It follows that the Luenberger observer ( 56) is an exponential observer for the system (55). Moreover the rate of convergence can be arbitrary chosen.

Example: Consider a Leslie type system for the dynamic of a three stage structured population

    ẋ1 (t) ẋ2 (t) ẋ3 (t)     =     -α 1 0 µ 3 β 1 -α 2 0 0 β 2 -α 3     A     x 1 (t)
x 2 (t)

x 3 (t)     .
The entries of the Leslie matrix A are assumed to be independent of the environment and we suppose that we observe only the last stage, that is, the output of the system is y(t) = x 3 (t).

Here, C = 0 0 1 . The transfer from stage i -1 to stage i being positive, i.e, β 1 > 0 and β 2 > 0, the Kalman observability rank condition is fulfilled. Therefore the following system is an exponential observer (or software sensor) for the above Leslie system:

     ẋ1 (t) ẋ2 (t) ẋ3 (t)      =     -α 1 0 µ 3 β 1 -α 2 0 0 β 2 -α 3         x1 (t) x2 (t) x3 (t)     + (y(t) -x3 (t))     k 1 k 2 k 3     .
If the above map Φ defines a global diffeomorphism then the system (60) can be rewritten:

                         ż =       0 1 0 ... 0 0 0 1 ... 0 . . . . . 0 0 ... 0 1 0 0 ... 0 0       A z +        0 0 . . . 0 ψ(z)        + u        g 1 (z) g 2 (z) . . . g n-1 (z) g n (z)        = F (z) + uG(z) y = z 1 = (1, 0, . . . , O) C z. (61)
If moreover the system is observable for any input then the map G satisfies g i (z) = g i (z 1 , . . . , z i ). This means that ∂ ∂z j Y X i-1 h(φ -1 (z)) = 0, for i = 1, . . . , n -1 and j = i + 1, . . . , n

If ψ and g i are globally Lipschitz then an exponential observer for the system (61) is, for θ large enough, given by:

ż = F (ẑ) + uG(ẑ) -S -1 θ C T (C ẑ -y). (62)
For the original system (60), the observer is:

ẋ = X(x) + uY (x) - ∂Φ ∂x -1 x=x S -1 θ C T (C x -y). (63)
For the system (60), it is possible to give another type of observer called the extended Kalman filter:

ż = F (ẑ) + uG(ẑ) -1 r S(t) -1 C T (C ẑ -y), Ṡ = -SQ θ S -[A (ẑ, u)] T S -SA (ẑ, u) + 1 r C T C, (64)
where r, θ are positive real numbers,

Q θ = ∆ θ Q∆ θ , Q is a given symmetric positive definite n × n matrix, A (ẑ, u) is the Jacobian matrix of F (z) + uG(z) evaluated at z = ẑ, i.e, A (ẑ, u) = ∂ ∂z (F (z) + uG(z)) z=ẑ .
Once again, for θ ≥ 1 and large enough, the system (64) is an exponential observer for the system (61). The difference with the Luenberger-like observer (62) is that the gain S -1 θ C T used in the correction term is not constant but it is dynamically computed as a solution of a Ricatti matrix differential equation and hence, takes into account the information appearing at the current time t. Therefore, the Kalman observer is more robust with respect to noise that may affect the measurement output.

Example: The chemostat. We consider the following chemostat model

       ẋ = µ(s) x -ux, ṡ = -kµ(s) x -u(s -s in ), y = x. (65)
Where µ(s) = ms a + s is the Monod absorption response of the substrate by the micro-organisms.

The domain Ω = {x = (x, s) ∈ IR 2 : x > 0, s > 0, s + kx < s in } is positively invariant under the flow of ( 65). This can be proved by checking that the vector field is always tangent or pointing inside the boundary of Ω. However, the system is not observable on the closure Ω of Ω because all the initial conditions (0, s) produce the same output y(t) ≡ 0 and so, for s 1 = s 2 , the states (0, s 1 ) and (0, s 2 ) are indistinguishable whatever the input is. Therefore, we shall choose an open subset Ω 1 of Ω such that the system is observable on Ω 1 . To this end we assume that 0 < u min ≤ u ≤ u max < m. Let d be the positive number defined by d = au min m -u max and c

be any positive number satisfying kc + d < s in . We, then, define Ω 1 by

Ω 1 = Ω ∩ {(x, s) : x > c, s > d}.
It is easy to see that Ω 1 is positively invariant and that the system is observable on Ω 1 for all input satisfying the condition 0 < u min ≤ u ≤ u max < m.

We perform the change of coordinates given by Φ : z 1 = y = x, z 2 = µ(s)x = ms a + s x.

Conversely:

Φ -1 : x = z 1 , s = az 2 mz 1 -z 2 .
The domain Ω is transformed by Φ in

Φ(Ω) = {(z 1 , z 2 ) ∈ IR 2 : z 1 > 0, 0 < z 2 < mz 1 , mkz 2 1 -kz 1 z 2 + (a + s in )z 2 -ms in z 1 < 0}, and the domain Ω 1 is transformed by Φ in Φ(Ω 1 ) = {(z 1 , z 2 ) ∈ IR 2 : z 1 > c, dm a + d z 1 < z 2 < mz 1 , mkz 2 1 -kz 1 z 2 +(a+s in )z 2 -ms in z 1 < 0}
The sets Ω and Φ(Ω) (corresponding to the parameters: a = 10, m = 1, s in = 30, k = 1) are drawn in Figure 14, and the sets Ω 1 and Φ(Ω 1 ) are drawn in Figure 15.

The map Φ is bijective from Ω 1 to Φ(Ω 1 ). The Jacobian matrices of Φ and Φ -1 are which are invertible for all (x, s) ∈ Ω 1 and all (z 1 , z 2 ) ∈ Φ(Ω 1 ). Therefore Φ is a diffeomorphism from Ω 1 to Φ(Ω 1 ) (it is also a diffeomorphism from Ω to Φ(Ω)). With the new coordinates, the system (65) is defined, for (z 1 , z 2 ) ∈ Φ(Ω 1 ), by

dΦ dx =   1 0 ms a + s max (a + s) 2   , dΦ -1 dz =    1 0 -amz 2 (mz 1 -z 2 ) 2 amz 1 (mz 1 -z 2 ) 2    ,
         ż1 = z 2 -uz 1 , ż2 = z 2 2 z 1 - kz 2 (mz 1 -z 2 ) 2 amz 1 + u -z 2 - (mz 1 -z 2 )z 2 mz 1 + (mz 1 -z 2 ) 2 amz 1 s in , y = z 1 . (66)
We extend this system to the whole IR 2 as follows:

             ż1 = z 2 -uz 1 , ż2 = ψ(z 1 , z 2 ) + ug(z 1 , z 2 ), y = z 1 , (z 1 , z 2 ) ∈ IR 2 , ( 67
)
where ψ and g are any C ∞ (or at least continuous) functions that are globally Lipschitz on IR 2 and such that their restrictions to Φ(Ω 1 ) are:

ψ/Φ(Ω 1 ) : (z 1 , z 2 ) → z 2 2 z 1 - kz 2 (mz 1 -z 2 ) 2 amz 1 , g/Φ(Ω 1 ) : (z 1 , z 2 ) → -z 2 - (mz 1 -z 2 )z 2 mz 1 + (mz 1 -z 2 ) 2 amz 1 s in
The system (67) has the same form as the system (61) and it satisfies all the required conditions on IR 2 . Therefore a Luenberger-like observer can be built. To this end, we compute the following matrices

S -1 θ = 2θ θ 2 θ 2 θ 3 , S -1 θ C T = 2θ θ 2 .
And then an exponential observer for ( 67) is ż1 = ẑ2 -uẑ 1 -2θ(ẑ 1 -y), ż2 = ψ(ẑ) + ug(ẑ) -θ 2 (ẑ 1 -y).

Therefore, for all α > 0, it is possible to find θ ≥ 1 such that the solutions of (68-67) satisfy:

∃ C 1 > 0, ∀t ≥ 0, ∀(ẑ(0), z(0)) ∈ IR 2 × IR 2 , ẑ(t) -z(t) ≤ C 1 ẑ(0) -z(0) exp(-αt). Now, let Φ : IR 2 → IR 2 be a global diffeomorphism such that Φ and its inverse Φ-1 are globally Lipschitz and such that the restriction of Φ to the compact set Ω 1 is equal to Φ, that is, Φ(x) = Φ(x) for all x ∈ Ω 1 . An estimation of the state vector x(t) = (x(t), s(t)) is then given by:

x(t) = (x(t), ŝ(t)) = Φ-1 (ẑ 1 (t), ẑ2 (t)) = Φ-1 (ẑ(t)),

where ẑ(t) is given by the differential system (68). This estimation converges exponentially to the real state x(t). Indeed, let K 1 and K 2 be the Lipschitz constants of respectively Φ and Φ-1 , then we have ∀t ≥ 0, x(t) -x(t) = Φ-1 (ẑ(t)) -Φ-1 z(t)

≤ K 2 ẑ(t) -z(t) ≤ K 2 C 1 ẑ(0) -z(0) exp(-αt)

≤ K 2 K 1 C 1 x(0) -x(0) exp(-αt).

Remark: Another way (which is less rigorous than the above one and which is often used for the simulations purpose) to compute the estimations (x(t), ŝ(t)) is just to apply the formula (63) without extending the system to the whole space IR 2 . To this end we need to compute the following matrices The observer for the chemostat is then given by:

∂Φ

     ẋ = µ(ŝ)
x -ux -2θ(x -x), ṡ = -kµ(ŝ) x -u(ŝ -s in ) --2 ŝ(a + ŝ)θ ax + (a + ŝ) 2 θ 2 amx (x -x).

(70)

In Figure 16, we present a simulation with the parameter values: a = 10, m = 1, s in = 30, k = 1 and with a constant input u = 0.1. We have chosen θ = 7. The first figure represents the evolution of, respectively, the biomass x and its estimation x. The second represents the evolution of, respectively, the substrate s and its estimation ŝ. The curves show the convergence of the estimations provided by the observer (70) to the states of the system (65) and it can be noted that the convergence is quite fast. x ^(t)

Figure 16: The convergence of the estimations x(t) and ŝ(t) to the real states x(t) and s(t).

Observers for a class of non-affine control systems

Consider the following analytic system defined on IR n

                       ẋ(t) =        ẋ1 (t) . . . ẋi (t) . . . ẋn (t)        =        f 1 (x 1 , x 2 , u)
. . . f i (x 1 , x 2 , . . . , x i , x i+1 , u) . . .

f n (x 1 , x 2 , . . . , x n , u)        = f (x, u) y = h(x 1 , u) (71)
where:

• Each of the map f i is globally Lipschitz with respect to (x 1 , . . . , x i ) uniformly with respect to x i+1 and u, i.e, there exists a constant M i that does not depend neither on x i+1 nor on u such that f i (x 1 , . . . , x i , x i+1 , u) -f i (z 1 , . . . , z i , x i+1 , u) ≤ M i (x 1 , . . . , x i ) -(z 1 , . . . , z i ) .

• There exist two real numbers α, β, 0 < α < β, such that

α ≤ | ∂h ∂x 1 | ≤ β, α ≤ | ∂f i ∂x i+1 | ≤ β, for i = 1, . . . n -1.
Let A(t) and C(t) be respectively n × n and 1 × n time-dependent real matrices: , C(t) = φ 1 (t), 0, . . . , 0 .

A(t) =        0 φ 2 (
If the functions φ i (t) satisfy α ≤ φ i (t) ≤ β, for i = 1, . . . n -1.

Then, there is a real λ > 0, a vector K ∈ IR n , and a symmetric positive definite n × n matrix S, depending only on α and β, such that:

A(t) -KC(t) T S + S A(t) -KC(t) ≤ -λId.

The above vector K allows to construct a Luenberger-like observer for the system (71) as follows: This observer is an exponential one. Moreover, it is possible to choose the convergence rate, that is, for any α > 0, one can find θ large enough in such a way that for any initial conditions (x 0 , x0 ), the corresponding solutions x(t) of the system (71) and x(t) of the observer (72) satisfy ∀t ≥ 0, x(t) -x(t) ≤ P (α)e -αt x0 -x 0 , where P is some polynomial of degree n.

ẋ = f (x, u) -∆ θ K(h(x 1 , u) -y), (72)
Remark: We have seen that the above observer constructions are made under the assumption that some functions have to be globally Lipschitz. This is a very restrictive condition. However, for real systems, the state space (or at least the set of interest) is often a bounded connected open subset Ω of IR n and so the global Lipschitz condition is met on Ω the closure of Ω. Therefore the functions considered can be extended to the whole IR n by smooth and globally Lipschitzian functions on IR n .

Asymptotic observers

For all the above observers, it is possible to assign an arbitrarily (fast) exponential rate of convergence of the error estimation. However, they require the full knowledge of the structure of the model, i.e, the function X(x, u) has to be exactly known and the system has to satisfy some strong observability conditions. When one of the two above conditions is not satisfied, it is still possible in some situations to construct an estimator of the state of the given system but in general, it will not be possible to choose the convergence rate of the estimate to the real state. To illustrate this, we give just two examples. The first one concerns linear systems that are not completely observable, i.e, the rank of the observability matrix is strictly less than the dimension of the state space. The second example concerns biotechnological processes whose dynamical models are partially known, especially the mathematical structure of the biological kinetics is unknown.

Observers for detectable linear systems: Consider again the linear system

             ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t),
x(t) ∈ IR n , u(t) ∈ U ⊂ IR m , y(t) ∈ IR q , A, B, and C are respectively n × n, n × m and q × n matrices,

and assume that the pair (C, A) is not observable, that is, the rank of the observability matrix O (C,A) is less than the dimension n of the state space. In this case the Luenberger observer does not work but one can still built an observer for the system (73) if it is detectable which means that the unobservable space corresponds to the stable modes of the matrix A. More precisely let r = rank O (C,A) , there exists a linear change of coordinates x = T z under which the equations ( 73) are transformed into

       ż1 = A 11 z 1 + B 1 u, ż2 = A 12 z 1 + A 22 z 2 + B 2 u, y(t) = C 1 z 1 , (74)
where z 1 ∈ IR r , z 2 ∈ IR n-r , A 11 is r × r, A 12 is (n -r) × r, A 22 is (n -r) × (n -r), B 1 is r × m, B 2 is (n -r) × m and C 1 is q × r. The pair (C 1 , A 11 ) is observable.

The system is detectable if A 22 is a stable matrix, that is, all its eigenvalues are with negative real part. An observer for (74) is then given by

Figure 1 :

 1 Figure 1: This diagram illustrates the dynamical transfer of the population.

  ds dt = b -bs -(β -d)si, de dt = βsi -(b + )e + dei, di dt = e -(r + d + b)i + di 2 , dτ dt = ri -bτ + diτ.

Figure 2 :

 2 Figure 2: Trajectories of X 0 and X 1 for a = 0.4, b = 0.01, c = 0.3, e = 0.5, q = 0.15.

Figure 4 :

 4 Figure 4: Reachable set from the state x 0 for the system (15).

Figure 5 :

 5 Figure 5: Attractivity without stability.

Figure 7 :

 7 Figure 7: The trajectories of (21) for a = 0.4, b = 0.01, c = 0.3, e = 0.5.

  ): • The positive quadrant, • {(N, 0), N ∈ IR}, • {(0, P ), P ∈ IR}, • {(N, P ) : e b N -c log( e b c N ) + b P -a log( b a P ) = k}, k being a positive constant.

  eN * -P * .

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: The set E where V = 0.

Figure 13 :

 13 Figure 13: An observer.

Figure 14 :

 14 Figure 14: The invariant domain Ω and its image by Φ.

Figure 15 :

 15 Figure 15: The invariant domain Ω 1 and its image by Φ.

  , and Stability of Mathematical Models -Abderrahman Iggidr 54

  ż1 = A 11 ẑ1 + B 1 u + K(y -C 1 ẑ1 ), ż2 = A 12 ẑ1 + A 22 ẑ2 + B 2 u.(75)Indeed, the estimation error satisfies the differential equation ė1 = (A 11 -KC 1 )e 1 , ė2 = A 12 e 1 + A 22 e 2 .

  The steady state x e is said to be asymptotically stable (We also say that (17) is asymptotically stable at x e ) if it is Lyapunov stable and attractive. It is globally asymptotically stable if it is Lyapunov stable and globally attractive.

						S(x	e , δ)
	. x	. x e	S(x	e , δ)	. x	x e	.
	X (x) t
						X (x) t
	is a stable equilibrium	is asymptotically stable

Asymptotic stability of the equilibrium x e means that all solutions starting near x e not only stay nearby, but also tend to the equilibrium x e as time goes to infinity.

Table 1 :

 1 The parameter values used in the simulation.

	Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 38
					23
					22
					21
					20
	2	4	6	8	19	0	2	4	6	8

Acknowledgements

The author is grateful to Prof. C. Lobry for the opportunity to contribute to the EOLSS and to Dr. P. Adda and Prof. G. Sallet for valuable discussions on the subject of the manuscript. A special thank to the anonymous referee for valuable comments and suggestions which helped the author to improve the presentation of this article.

Stabilization

Now, we consider a controlled system with an equilibrium state at the origin of IR n    ẋ(t) = X(x(t), u(t)), X(0, 0) = 0, x(t) ∈ IR n , u(t) ∈ U ⊂ IR m .

(31)

We assume that the function X is sufficiently smooth and we want to know if it is possible to control this system in order to stabilize it at its equilibrium state. There are two ways to achieve this goal. The first one consists in finding the control u as a function of the time variable t and corresponds to what is called an open-loop control. The second strategy is to build the stabilizing control as a function u(x) of the state x, this is a feedback control or a closed-loop control. Here, we develop the basis of the second strategy.

We shall say that the system (31) is stabilizable if there exists a state feedback control u(x), which is at least continuous ( as a function of the state x), and such that the origin is an asymptotically stable equilibrium point for the closed-loop system ẋ(t) = X(x(t), u(x(t))).

The function u(x) is called a static stabilizing feedback law. The system (31) can also be dynamically stabilizable by use of a dynamic feedback. This means that we control the system by using other systems. So, we add to the system (31) another dynamical system whose input is x. Hence, we obtain an interconnected system:

We shall then say that the system (31) is dynamically stabilizable if there exists a feedback law u(x, y) that stabilizes the system (32) at (x, y) = (0, 0). In this article, emphasis will be put on the static stabilization.

A necessary condition for the existence of a C 1 static feedback u(x) that stabilizes locally the system (31) around its equilibrium point 0 is that (i) there exists a neighborhood N of the origin such that for each state ξ ∈ N there is a control u ξ steering the system from x = ξ at t = 0 to x = 0 at t = ∞, i.e., lim

contains a neighborhood of the origin; (iii) the linearized system has no uncontrollable modes (see Controllability of linear systems 2.2) associated with eigenvalues whose real part is positive. A system ẋ = X(x, u) satisfying the condition (i) is said to be asymptotically controllable to the origin or open-loop stabilizable.

The conditions (i) and (ii) are necessary even if one requires that the stabilizer u(x) is only continuous. The condition (iii) is not necessary if the stabilizer is required to be only continuous. When the condition (i) holds but the condition (ii) is not satisfied, then it is not

Observability

We consider the standard input-output finite dimensional system:

Here M is an open connected subset of IR n . We assume that the vector field X and the map h are C ∞ or analytic (when needed). The point x(t) in M is called the state of the system, u(t) is the input and y(t) is the measurable output of the system. The introduction of y is due to the fact that usually we do not have access to the whole state: we can observe or measure only a part of the actual state of the system. The class of admissible inputs U ad = {u(.) : t ∈ IR + → u(t) ∈ U } is contained in the set of measurable controls with values in U and contains the class of piecewise constant controls with values in U . We shall use the notation X u(.) t (x 0 ) to denote the solution of the differential equation ( 44) corresponding to the admissible control u(.) ∈ U ad and with initial condition x 0 : X u(.)

Two states x 1 , x 2 ∈ M are said to be indistinguishable if y (x 1 , u(.), t) = y (x 2 , u(.), t) for any admissible control (or input) u(.) ∈ U ad and any t for which both sides are defined. Roughly speaking, this means that the information provided by the measurable output is not enough to tell us if the evolution of the system is given by the solution of (44) emanating from the state x 1 or by the one emanating from the state x 2 .

Example: Consider the single input system with two outputs

We have, for any real x, any integer k and any input u

Hence the states x and x + 2kπ are indistinguishable.

Two states x 1 , x 2 ∈ M (x 1 = x 2 ) are said to be distinguishable if there exists an admissible control (or input) u(.) ∈ U ad and a time t ≥ 0 such that y (x 1 , u(.), t) = y (x 2 , u(.), t). An admissible input which distinguishes every pair of states is called an universal input.

Consider the single input single output system

Examples from life support systems

Leslie type systems: Models describing the growth of a stage structured population are usually represented by Leslie-type systems :

In these models, all the variables are positive, x 1 represents the youngest stage (often eggs). The function F 1 describes the so-called recruitment function, i.e., eggs laying from the other stages. The function x i-1 → F i (x i-1 , ., .) corresponds to the transfer from stage i -1 to stage i, and it is often an increasing function. The fishery model (40) is a particular example of Leslie-type systems (49) where the input is the fishing effort. Suppose we measure only the last stage (the oldest), that is the output of the system (49) is y = x n . This system is observable for any input. Indeed, suppose that for two initial states x and x we have the same output, i.e., x n (t) = xn (t), for all t ≥ 0 then by differentiation we get:

Since ∂F n /∂x n-1 is positive, the map x n-1 → F n (x n-1 , x n , u) is injective. Hence, the above equality (together with x n (t) ≡ xn (t)) implies that x n-1 (t) ≡ xn-1 (t). Once again, by differentiation of this equality we get:

Thanks to the same argument (x i-1 → F i (x i-1 , x i , u) is monotone), we deduce that x n-2 (t) ≡ xn-2 (t). And so on, we show that for all t ≥ 0,

Hence, y(x, u(.), t) ≡ y(x, u(.), t) implies that x = x which proves that the Leslie-type system ( 49) is uniformly input observable.

Trophic chains: Models describing trophic chains represent the dynamics of an ecosystem from nutrient (x 1 ), phytoplankton (x 2 ), ... to higher levels such as fish (x n ). They can be written as

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 45

The functions F i are generally monotone functions of the variable x i+1 , i.e., x i+1 → F i (., ., x i+1 , .) is monotone. These systems are uniformly input observable if one can measure on-line the first compartment of the chain or the last one, i.e., y = x 1 or y = x n .

Observers

We consider an input-output nonlinear system described by

where x(t) ∈ IR n is the state of the system at time t, u(t) ∈ U ⊂ IR m is the input and y(t) ∈ IR q is the measurable output of the system.

Here, f : IR n × IR m → IR n and h : IR n × IR m → IR q are smooth functions.

We have seen that the observability of the system (51) means that it is theoretically possible to recover the state x(t) from the input u(t) and the output y(t) together with their repeated time derivative u, ü, . . . , ẏ, ÿ, . . . , y (k) , . . . along the solution of the system. However, in practice, the use of the derivatives may not give sufficiently accurate performance, especially in the case of noisy measurements as it can be illustrated by the following example:

This system is observable for any input. Moreover, the state can be completely recovered because x 1 is measured and one can deduce x 2 and x 3 by x 2 = ẏ and x 2 = ÿ. However, if the measurement of x 1 is corrupted by a sinusoidal disturbance with a small amplitude then the error on x 2 and x 3 (computed according to the above rules) can become very large. For instance, if y(t) = x 1 (t) + sin(100t) then ẏ(t) = x 2 (t) + 100 cos(100t) and ÿ(t) = x 3 (t) -10000 sin(100t).

To avoid the use of the time derivatives of the output and at the same time to counter the induced effects of disturbances, another tool has been developed to reconstruct the state from available outputs. This tool is called an observer or a software sensor. An observer for the the system ( 51) is a dynamical system whose inputs are the inputs and outputs of the system (51), which produces an estimate x(t) of the state x(t) such that x(t) -x(t) → 0 as t → +∞ and the estimation error must remain small if it starts small.

The gain matrix K = k 1 k 2 k 3 T can be selected in order to force the eigenvalues of A-KC

to have desired values. For instance, if one wants to have σ(A -KC) = {-3, -2, -1}, then the coefficients of the matrix K are:

The corresponding estimation error satisfies for all initial conditions x(0) and x(0) (recall that x(0) is unknown but x(0) can be chosen by the user) and any positive time t e(t) ≤ c x(0) -x(0) e -3t .

Discrete-time systems: Now we consider an n-dimensional discrete-time linear system

If the pair (C, A) is observable then, by the Pole-Shifting Theorem, it is possible to find a matrix K in such a way that A -KC is a nilpotent matrix, i.e, (A -KC) n = 0. The corresponding observer x(k + 1) = Ax(k) + Bu(k) + K y(k) -C x(k) is called a deadbeat observer. It has the particularity that the estimate becomes equal to the state after n steps.

Some nonlinear observers

The construction of nonlinear observers attracted much attention during the last decades.

Efficient methods have been designed by many authors. For lack of space we shall only expose the high-gain constructions developed recently by J.P. Gauthier, H. Hammouri, I. Kupka and S. Othman. To avoid complex calculus, we consider only single output nonlinear system.

System with no input

Let us consider an observable single output system

The output function h together with its first n -1 derivatives along the vector field X allow to define the following map

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 49

We assume that Φ is a global diffeomorphism. This assumption implies that the state can be recovered from the output and its first n -1 time-derivatives. In the coordinates defined by z = Φ(x), the system (57) is governed by the following differential equation:

Assume that ψ is globally Lipschitz. Then, for θ ≥ 1 large enough, an exponential observer (a Luenberger type observer) for the system (58) is given by the following dynamical system:

where

, and K is chosen such that the matrix A -KC has all its eigenvalues with negative real part. This is possible since the pair (C, A) is observable.

Remark: In the equation of the observer (59), the term ∆ θ K can be replaced by S -1 θ C T with S θ being the solution of

The matrix S θ can be analytically computed by

Affine control systems

Consider a single input single output (SISO) nonlinear system that can be modeled by

On the one hand, since the pair (C 1 , A 11 ) is observable, it is possible to choose the matrix K in such a way that the matrix A 11 -KC 1 has all its eigenvalues with negative real part. Therefore e 1 (t) converges to zero. On the other hand, the matrix A 22 is stable so e 1 (t) converges to zero as well. However, we must notice that one can choose the convergence speed of e 1 (t) whereas the convergence rate of e 2 (t) is completely determined by the eigenvalues of A 22 .

An observer for a biotechnological process: Let us consider a more general model of the chemostat:

(76)

Here we assume that the concentration of micro-organisms x(t) is measured on-line. It is reasonable to assume that the dilution rate (or the flow rate) u(t), the substrate concentration s in in the inflow and the constant k are known. However the function µ(x, s) is is often a very complex function of the state of the process which depends on the operating conditions. Its analytical expression is still the subject of intensive investigation. For the same process there exist several possible models in the literature, so, which one to choose? For this reason, it is not possible to use one of the above exponential observers. Therefore another type of estimators has been developed for this kind of systems. These estimators are often called asymptotic observers. For the chemostat (76) an asymptotic estimator for the unmeasured substrate concentration s(t) can be built as follows. We perform a linear change of coordinates:

With these new coordinates, the system (76) becomes: ż1 = µ(z 1 , z 2 -kz 1 ) -u(t)z 1 , ż2 = -u(t)z 2 + u(t)s in .

The advantage of the above change of coordinates is that the dynamic of the unmeasured variable z 2 does not depend on the unknown function µ(x, s). Therefore a candidate estimator for z 2 can be given by ż(t) = -u(t)ẑ(t) + u(t)s in .

The dynamic of the estimation error e(t) = ẑ2 (t) -z 2 (t) is ė = -u(t)e.

If the dilution rate u(t) satisfies the condition 0 < u min ≤ u(t), for all positive time t, then |e(t)| ≤ |e(0)|e -u min t , ∀t ≥ 0.

Controllability, Observability, and Stability of Mathematical Models -Abderrahman Iggidr 58

Hence an estimate of the unmeasured substrate concentration s(t) is given by ŝ(t) which is computed by ż(t) = -u(t)ẑ(t) + u(t)s in ŝ(t) = ẑ(t) -kx(t).

The convergence of the estimator is exponential but its rate of convergence can not be chosen by the user, it is completely determined by the dilution rate u(t).