Identification of linear systems with nonlinear Luenberger Observers

The design of a nonlinear Luenberger observer for a parametrized linear system is studied. From an observability assumption of the system, the existence of such an observer is concluded. In a second step, a constructive novel algorithm for the identification of multi-input multi-output linear systems is suggested and is implemented on a second order system.


I. INTRODUCTION
In this note, a parametrized linear system described by the following equation is considered: where θ in Θ ⊂ R q is a vector of unknown parameters and Θ is a known set, u in R m is a controlled input (in L ∞ loc (R + )).The state x is in R n and y is the measured outputs in R p .Mappings A : Θ → R n×n , B : Θ → R n×m and C : Θ → R p×n are known C 1 matrix valued functions.
In the following, an estimation problem is considered: The design of an observer to estimate the state and the unknown parameters of the system from the knowledge of y is seek.In other word, an asymptotic observer for the extended (nonlinear) n + q dimensional system ẋ = A(θ )x + B(θ )u , θ = 0 , y = C(θ )x (2) has to be designed.The nonlinear Luenberger methodology introduced in [8] (see also [13], [5], [7], [2], [1]) is a method which permits to design an observer based on weak observability assumptions.Following this approach, the first step is to design a C 1 function (x, θ , w) → T (x, θ , w) such that the following partial differential equation (PDE) is satisfied: where Λ is a Hurwitz squared matrix, L a matrix and g is a controlled vector field.Dimensions of the matrices and of the vector field g must be taken consistently.This will be precisely defined in the following.The interest on this mapping is highlighted if (z(•), w(•)), the solution initiated from (z 0 , w 0 ) of the dynamical system, is considered: ż = Λz + Ly , ẇ = g(w, u) .
In other words, z provides an estimate of the function T .The second step of the Luenberger design is to left invert the function T in order to reconstruct the extended state (x, θ ) from the estimate of T .Hence, a mapping T * has to be constructed such that Of course, this property requires the mapping T to be injective.The final observer is simply A particularly interesting feature of this observer is that its convergence rate can be made as large as requested (see [1]).In this paper, this strategy is adopted to suggest a solution to the state and parameters estimation for system (2).This paper can be seen as an extension of the result of [12] in which a nonlinear Luenberger observer is constructed for a harmonic oscillator which fits in the class of system studied.
As will be seen in Section IV, our approach can be employed to identify the state and all parameters describing an unknown linear model.This topic has been widely studied in the literature (see the books [4], [10], [11]).Adaptive identifiers can be traced back to G. Kreisselmeier in [6].In this paper three algorithms were given for identification of a SISO system in a specific controllable and observable representation.This work has then been extended in many directions to allow time varying matrices and MIMO systems (see for instance [3], [9], [14]).
The algorithm given in Section IV gives a new approach to address the same problem.It has the advantage to allow for a prescribed convergence rate for MIMO systems.Moreover, its structure is relatively simple (see ( 27)) and its dimension significantly smaller than that of available algorithms.One of the main difference is that, in contrast to all other available approaches, no optimization step is required.
The paper is organized as follows: In Section II, the existence of a mapping T which solves the partial differential equation ( 3) is discussed.Section III is devoted to the study of the injectivity of the mapping T assuming some observability properties.In Section IV, a left inverse of the mapping T is constructed to get the observer when considering a specific structure for the matrices A, B and C.This leads to a novel solution for the identification of MIMO linear time invariant systems.

Notations:
• Given a matrix A in R n×n , σ (A) denotes its spectrum.
• 1 n denotes the n dimensional real vector composed of 1.
• I n denotes the n dimensional identity matrix.
• ⊗ stands for the Kronecker product between two matrices.
• For a vector or a matrix | • | denotes the usual 2-norm.
• Given a set C, Cl(C) is its closure.

II. EXISTENCE OF THE MAPPING T
In [2], it is shown that, in the autonomous case the existence of the mapping T , solution of the PDE (3), is obtained for almost all Hurwitz matrices Λ.In the context of the controlled system (2), the same type of result still holds.Moreover, an explicit solution of the PDE (3) may be given.
Theorem 1 (Existence of T ): Let r be a positive integer.For all r-uplet of negative real numbers (λ 1 , . . ., λ r ) such that Proof: Keeping in mind that the spectrum of Λ and A(θ ) are disjoint (see (6)), the matrix M i (θ ) in R p×n can be defined for all i in {1, . . ., r}: For all i = 1, . . ., r, let T i : R n × Θ × R m → R p and the vector field g i : R m × R m → R m be defined as: It can be noticed that T i is solution to the PDE Hence, the solution of the PDE (3) is simply taken as This ends the proof.
Remark 1: Note that in the particular case in which the system is autonomous and there is only one output, the mapping T is given as This matrix is solution to the following parametrized Sylvester equation Hence, taking r = n, the well known Luenberger observer introduced in [8] in the case of autonomous systems is recovered.Note however that, here, the injectivity is more involved than in the context of [8] since θ is unknown.
Remark 2: Note that if the set Θ is bounded, then it is ensured that there exists (λ i )'s which satisfy equation (6).Indeed, if Θ is bounded, then the set ( θ ∈Θ σ {A(θ )}) is a bounded set.This can be obtained from the fact that each eigenvalues λ in σ {A(θ )} is a zero of the characteristic polynomial: where µ i (θ ) are continuous function of θ .Boundedness of Θ together with the continuity of the µ i 's imply that there

III. INJECTIVITY
As seen in the previous section, it is known that if the following dynamical extension is considered: with z in R rp and w in R rm , then it yields that along the solution of the system (2)-( 12), equation ( 4) is true.Consequently, T (x, θ , w) defined in ( 7)-( 8) is asymptotically estimated.The question that arises is whether this information is sufficient to get the knowledge of x and θ .This is related to the injectivity property of this mapping.As shown in [2], in the autonomous case this property is related to the observability of the extended system (2).With observability, it is sufficient to take r large enough to get injectivity.Here, the same type of result holds if it is assumed an observability uniform with respect to input in a specific set.
In order to simplify the presentation, the SISO case is considered.The following strong observability assumption is made: Assumption 1 (Uniform differential injectivity): There exist two bounded open subsets C θ ⊂ Θ, C x ⊂ R n , an integer r and U r a bounded subset of R r such that the mapping and S is the shift matrix operator such that for all s = (s 1 , . . ., s r ), S × s = (0, s 1 , . . .s r−1 ) is injective in (x, θ ), uniformly in (x, θ , v) ∈ Cl(C θ )×Cl(C x )×U r and full rank.More precisely, there exists a positive real number L such that for all (x, θ ) and The following result establishes an injectivity property for large eigenvalues of the observer.
Theorem 2: Assume Assumption 1 holds.Let u(•) be a bounded C r−1 ([0, +∞]) function with bounded r − 1 first derivatives, i.e. there exists a positive real number u such that For all r-uplet of distinct negative real numbers ( λ1 , . . ., λr ), for all positive time τ and for all w 0 in R r , there exist two positive real numbers k * and L T such that for all k > k * the mapping defined in ( 7)-( 8) with λ i = k λi , i = 1, . . ., r satisfies the following injectivity property in ) is in U r then for all (x, θ ) and (x * , θ * ) in C x × C θ the following inequality holds: where w(•) is the solution of the w dynamics in (12) initiated from w 0 .
Proof: First of all, picking k sufficiently large implies that the matrix M i which satisfies is well defined.On another hand, assume that k is sufficiently large such that for all i in [1, r], This implies that, for θ in C θ : Note moreover that R i is a C 1 function which satisfies for θ in C θ : Moreover, M i can be stated in the form which will be useful in the following: Let K be the matrix in R r×r defined as Note that: On another hand, for all t, since u is C r and w(•) being solution of (12) (in the SISO case), one get: which implies that w i satisfies: where: R wi (t) = w and with: Hence, with (13), it yields that for all t : where C is a positive real number which depends on w i (0) and (u(0), . . ., u (r−2) (0)).Keeping in mind that λi is negative, when t is larger than τ > 0, the previous inequality satisfies where R wi0 depends on k but not on t.This implies that by collecting terms of higher order in 1 k in a function denoted R MBi : and with (16) and using the fact that C θ and u ( j) (t) are bounded, for all t ≥ τ: Finally: and for all (x, θ ) in C x × C θ which is bounded and all t ≥ τ, it yields: This implies: where Ṽ in R r×r is the Vandermonde matrix defined as: and R T is a C 1 function which satisfies for all (x, θ Hence the mapping R T is globally Lipschitz with a Lipschitz constant in o 1 k r .Hence, it is possible to find k 0 such that for all k ≥ k 0 and all quadruples (x, x * , θ , θ * ) in C 2 x × C 2 θ , for all t ≥ τ: It can be shown that the result holds with this value of k 0 .Employing (17), it yields that, for all t: Consider now t 1 ≥ τ, the last term of the previous inequality can be lower-bounded by (18).Moreover, if u(t 1 ) . . ., u (r−1) (t 1 ) is in U r , the other term can be lowerbounded based on Assumption 1 and the result follows.

A. State and parameter estimation algorithm
A concrete algorithm for computing simultaneously the state and the parameter vector θ of the linear system ( 1) is now derived.Such estimates can be obtained as: where z(t) is the value obtained by integrating the dynamical system (12) starting from an arbitrary user-defined initial value (z(0), w(0)).The formulation (19) can be seen as one way of inverting the mapping T (see ( 5)).A distinctive feature of the optimization problem (19) which is worth pointing out is that all data involved in it are signals at time t.
Hence there is no explicit use of the data relevant to the past.Solving the problem (19) can be very challenging, because the mapping T is in general nonlinear and can be highly complex.As a consequence, the optimization problem ( 19) is non-convex so that solving it directly might require an iterative searching scheme.One simple way to tackle it is as follows.Let im T ⊂ R n z , n z = pr, denote the range of the map T .Assume that z(t) ∈ im T for some t.Such an assumption seems quite reasonable since z(t) converges exponentially to im T as stated by (4).Under this assumption, the estimates x(t) and θ (t) can be searched such that: But explicit determination of the expression of those estimates is still missing.For this purpose, let us consider particular canonical structures for the matrix-valued functions A, B, C in (1) : Note that assuming such structures (20) for A, B,C is without loss of generality: any input-output behavior of a linear MIMO system can be described with a model of this structure.Such a realization is however not necessarily minimal but it is observable for any vector θ .The problem of identifying canonical representations of linear systems has been deeply investigated in the last three decades employing mainly standard identification techniques (see for instance [4] and references therein).Our approach offers a new algorithm to address it.The interest of this one comes from the fact that the dimension of the proposed observer may not be too high and that the convergence rate is tunable (see [1]).Let the mapping T be as in ( 7)- (8).
To carry out the calculations, it is convenient to observe that the matrices A(θ a ) and C in (20) can be written as: where ⊗ refers to the Kronecker product, Exploiting these expressions and doing a little manipulation through Kronecker algebra lead to Let Applying the Sherman-Morrison-Woodbury formula, one gets: where it is assumed that 1 − e T 1 J −1 i θ a = 0. Indeed it can be easily checked that this is true provided the assumption (6) holds.Combining ( 22) and (23) gives us: , which, together with (21), reveals that: Rearranging the previous expression, one gets: Finally, the following linear equation for ( θa , θb , x) must be solved: where P(z, w) = P 1 (z, w) whenever the matrix P is full rank.However, the solution (26) can be efficiently computed only if the matrix P(z(t), w(t)) has full column rank at any time t .This is indeed related to the injectivity of the map T and more precisely to its rank.Indeed, let us give the following proposition: Proposition 1: Assume that the matrix-valued functions A, B, C in system (1) have the canonical structure (20).Assume moreover that condition (6) holds.Then, for all (z, x, θ , w) such that z = T (x, θ , w) and such that: then the matrix P defined in (24)-(25) has full column rank.
Proof: The proof follows from equation (25) in the case in which z i = T i (x, θ , w).Indeed, differentiating with respect to (x, θ ) this equation yields to: This implies that Hence: Again, since condition (6) holds, it yields that Diag(1 − V T 1 θ a , ..., 1 − V T r θ a ) is invertible and consequently P is full rank if ∂ T ∂ (x,θ ) (x, θ , w) is.Consequently, since it is known that z converges asymptotically to im T , if T is injective and full rank, the observer given in (26) is well defined after a certain time.However, in the transient period, there is no guarantee that this matrix is well defined.This problem has been deeply investigated in [12] considering an autonomous second order system with only one parameter.

B. Numerical illustration
In order to show the performances of the observer ( 12)-(26), a second order linear system in the form (1) is simulated in Matlab where: 0 ) θ i 's, the initial value x(0) and the input signal are given in Table I.
u(t) = sin(t) + sin(2.2t)+ sin(4.5t)+ sin(6.7t)+ sin(7t) It is assumed that the state (x 1 , x 2 ) and the parameters (θ 1 , θ 2 , θ 3 , θ 4 ) are unknown and are estimated with: where, P is obtained from (24).Λ, L and the initial configuration of z(0) and w(0) are given in Table II.The ordinary differential equations in (27) are solved with the Matlab lsim function.As is proved in Theorem 2 and Proposition 1, the rank of the matrix P is related to observability property and to the number of λ i 's considered.In this simulation, there are r = 7 eigenvalues which is more than the least that can be taken, since there are 6 unknown parameters.The simulation results are depicted in Figures 1-2: the estimated states and parameters converge to their true values after a transition phase due to the arbitrary choice of z(0)  and w(0).The convergence speed of this observer can be tuned arbitrarily by selecting larger eigenvalues λ i 's (see [1] for more details).From Theorem 2 and Proposition 1 it is inferred that if the input is properly selected, and if one can take a sufficiently large number of eigenvalues, the matrix P is asymptotically full rank (however, in the transient period, there is no such guarantee).This is the reason why in the algorithm, when the matrix P is too badly conditioned, the parameters are kept constant and the state estimate is obtained by integrating the model.For instance, this occurs at the beginning of the simulation (up to the time 4s).To see the observer robustness, a pseudo random noise with standard uniform distribution on the interval ±10%y(t) is added to the output signal.The figures 3 and 4 show the efficiency of the estimation algorithm with such output noise.

V. CONCLUSIONS
The design of a nonlinear Luenberger observer to estimate the state and the unknown parameters of a parametrized linear system was studied here.In a first part of the study, a Luenberger observer was shown to exist.This result is obtained from the injectivity property of a certain mapping.This result has been used on a second order linear system example and seems to offer very promising results.

TABLE I SYSTEM
CONFIGURATION.

TABLE II OBSERVER
CONFIGURATION.