A robust identification method for stochastic nonlinear parameter varying systems

1.   Introduction

Linear Parameter Varying (LPV) models have become progressively popular over the last two decades ^[20,28,43]. LPV models retain the design synthesis advantages of the Linear Time Invariant (LTI) setting, while being able to represent nonlinear dynamics with truthfulness. The LPV framework has been applied for control ^[37], observation ^[29] and identification ^[3,5] purposes. In this paper, we address the two latter topics.

Recent works have shown the particular interest of the class of Nonlinear Parameter Varying (NLPV) models ^[6,33]. In contrast to the "regular" LPV paradigm, which encapsulates the nonlinearities into time-varying scheduling parameters $\rho$, NLPV models retain some explicit nonlinearities that are "easily handled" outside of the scheduling term $\rho$, such as local Lipschitz nonlinear terms. On one hand, the LPV class exhibits linearity w.r.t. the state space, on the other hand, the NLPV class is nonlinear on both parameter and state spaces. The main interest in using the explicit nonlinearities coupled to the LPV structure is that they can be taken into account by the control design method, reducing the over-bounding of the nonlinearities by the scheduling parameter and, thus, enabling less conservative control performances. Anyhow, the development of identification and state estimation methods specifically conceived for NLPV models is yet unseen.

Over the past few decades, particular attention has been devoted to nonlinear system identification under non-Gaussian measurement noises^[8,9,36,40]. These stochastic tools are generic and could be applied to the case of NLPV models, yet they could benefit from the availability of the scheduling parameter data if properly designed.

Motivated by this research gap, the aim of this paper is to extend and generalise the Masreliez-Martin filter ^[40] joint state and parameter estimation method for the class of NLPV models with explicit Lipchitz nonlinearities. We stress that the main interest regarding the considered method is that it offers a simple recursive implementation and convergence properties verifiable through an LMI problem. The main novelty is that, by exploiting the local Lipschitz condition of the NLPV model, better state and parameter estimates are obtained. Furthermore, the use of the scheduling parameter data also enhances the results, since more information regarding the system dynamics is available. Accordingly, the major contributions are:

● The proposition of a robust identification procedure for joint parameter and state estimation for the class of NLPV models. Assuming that the explicit nonlinear term is structurally known and Lipschitz, the algorithm is is based on an adaptation of the modified Masreliez-Martin filter. The method considers coloured noise upon the measurement data and a random walk behaviour for the parameter variations.

● A sufficient LMI condition is provided for the verification of the local convergence of the proposed parameter and state estimation tool. This LMI remedy is derived from the direct Lyapunov method applied over an Ordinary Differential Equation (ODE) associated to the estimation error. The solution of the LMI also serves as an estimate for the zone of attraction of the algorithm. We also discuss how the the Lipschitz constant directly affects the size of this basin of attraction.

● The application of the proposed NLPV identification procedure to a mechatronic test-rig of an automotive semi-active suspension system. The effectiveness of the proposed method is discussed in terms of state estimation error and parameter identification precision, measured through standard indexes.

Remark 1. We stress that in the previous works by the Authors ^[40], the joint parameter-state estimation of a nonlinear stochastic system is considered and verified through simulations on a mathematical model. In this paper, the joint estimation algorithm of a stochastic NLPV system is considered, for which the nonlinear part is explicit, a priori known, and obeys a local Lipschitz condition. The proposed algorithm is verified through experimental results on a real vehicle system. Due to LPV parameter-dependent model matrices and to the nonlinear Lipschitz condition of the explicit nonlinear term, the proposed algorithm is completely different from the prior ^[40] in terms of matrix structures, gains and recursive estimation laws, but also conceptually, exploiting the availability of scheduling parameters to refine the joint estimation results.

The rest of this paper is organised as follows. The background state-of-the-art regarding nonlinear and LPV identification is discussed is Section 2, where the formal problem setup is also given. The proposed NLPV state and parameter estimation algorithm is presented in Section 3. The LMI problem used to verify the local asymptotic convergence of the method is also presented therein. The experimental validation results are given in Section 4. This paper ends with conclusions in Section 5.

1.1.   Notation

In the sequel, the set of nonnegative real number is denoted by $\mathbb{R}_+$, whist the set of nonnegative integers including zero is denoted by $\mathbb{N}$. $\mathbb{R}^q$ represents a $q$-dimensional real space. The index set $\mathbb{N}_{[a, b]}$ represents $\{i \, \in \, \mathbb{N} \, | \, a \, \leq \, i \, \leq \, b\}$, with $0 \, \leq \, a \, \leq \, b$. The identity matrix of size $j$ is denoted as $\mathbb{I}_{j}$; col$\{a \, , \, b \, , \, c\}$ denotes the vectorization (collection) of the entries and diag$\{v\}$ denotes the diagonal matrix generated with the line vector $v$.

2.   Background and problem formulation

System identification ^[2] is a topic of great importance for coherent control synthesis and implementation. The vast majority of control strategies are model-driven, which means that their success relies on the accurate knowledge of the model parameters. Accordingly, the issue of system identification under stochastic process signals deserves special attention ^[10,27,39], since real processes are not purely deterministic and stochastic noises are always present.

Besides the identification problem regarding model parameters, state estimation is a crucial aspect for the control of processes modelled through state-space representations, for which only output data is available. The classic state-feedback design requires online information of the state variables, which can be provided by observers/estimators. Due to these reasons, we note that joint state and parameter estimation is very useful and much facilitates the operation of adaptive and fault-tolerant state-feedback methods, since state observation and parameter variation estimates are provided by a single operator.

The Extended Kalman Filter (EKF) is a very well established method ^[25] for the joint estimation of model parameters and system states, as seen in different applications ^{[17,47,49,50]}. The EKF is popular for many reasons, but especially because it is not a complex algorithm, it has a recursive implementation law, and ensures asymptotic convergence. Its main drawback is that it is based on the assumption that the measurement noise has a Gaussian distribution. Many empirical and practical essays have shown that this assumption is rather weak ^[40]. It is a fact that measurement data from real processes may present inconsistent values at some points, i.e. "outliers". These corruptions are typically non-Gaussian and significantly reduce the effectiveness of algorithms such as the EKF. Research has progressively focused on overcoming the EKF restrictions; some of the recent approaches are: maximum likelihood methods ^[24,45], prediction error techniques ^[42] and robust identification algorithms ^[40]. These are all computationally less expensive than the EKF and require less restrictive assumptions on the measurement data.

The known fact is that the measurements have inconsistent observations with the largest part of the observation population. Justification of the proposed approach, which considers the presence of outliers, was confirmed in practice ^[21,32]. The presence of outliers can destroy the good features of linearly recursive algorithms, which are typically designed for the estimation in the presence of Gaussian noises. Therefore, it is very important to design a robust algorithm which is insensitive to outliers, since they exist in all real instrumented systems.

We note that many of the popular identification methods, such as the KF and the EKF, have been extended for the case of nonlinear ^[11] and LPV models ^[4]. For the latter class of models, the known parameter-varying data is used to schedule the estimation problem, which has input-output/input-state linearity. The state-of-the-art regarding LPV identification includes subspace techniques ^[12,15], successive LTI approximations ^[13], robust filtering ^[5], multi-step canonical variate analysis ^[22,35], prediction error methods ^[51] and nonlinear programs ^[23]. The review book ^[14] thoroughly discusses these recent developments and the future trends on LPV identification. A clear investigation gap is to consider NLPV models, for which there exists an explicit nonlinear function which must be taken into account by the estimation procedure, since input-output/input-state linearity no longer holds. Up to the Authors' knowledge, there has been no technique specifically devoted for the class of NLPV models, albeit the fact that many of NL identification methods can be applied to the NLPV case. The focus of this paper is given w.r.t. to this matter, pursuing the generalisation of a nonlinear joint parameter and state estimation method to the NLPV class, by exploiting the availability of the scheduling map data and the properties of the explicit nonlinearity. The purpose of exploiting the explicit nonlinearity is to seek better estimation results than just applying the nonlinear method itself.

Nonlinearity is a key feature present in almost all physical system controlled over large operating conditions. The extended Masreliez-Martin filter (EMMF) ^[26] is a an interesting option for the robust identification of nonlinear systems under stochastic, non-Gaussian corruptions. This robust method is based on Huber's statistics foundations ^[21]. Instead of using a Gaussian distribution assumption as the EKF, the EMMF is synthesised considering that the noise is coloured and belongs to some known distribution class. Furthermore, this algorithm ensures robustness and insensitivity w.r.t. the differences between this distribution class and the real outlier distribution. This paper considers the adaptation of the EMMF to the NLPV case due to its rather simple implementation and easily verifiable convergence property, as previously discussed ^[40,41]. Furthermore, we stress that the EMMF framework offers offset-free joint parameter and state estimation, with a zone of attraction that can be estimated from an LMI problem.

Remark 2. System identification theory has been set up since the $70$'s as the problem of generating parameter estimates in a stochastic framework ^[25,38], since input data samples are consider as deterministic signals corrupted by noise. As demonstrated by then ^[16,18,38], the input signal needs to be independent from the process noise for the estimation of the parameter to be unbiased, in the majority of cases. This property remains true in the context of this work. We note that in the experimental validation results (Section 4), we consider a PRBS input signal which is independent from this process noise.

2.1.   Problem formulation

In this paper, we generalise and apply the EMMF filter filter for the joint state and parameter estimation problem of NLPV models with strcuturally known Lipschitz nonlinearities. The task of estimation both states and parameters is done by using a joint formulation. For this goal, the following first-order affine Lipschitz NLPV system model is considered:

with model matrices given as follows:

being $x(k) \, \in \, \mathbb{R}^{n_x}$ the state vector (to be estimated), $\theta (k) \, \in \, \mathbb{R}^{n_p}$ the parameter vector (also to be estimated), $y(k) \, \in \, \mathbb{R}^{n_y}$ the output measurement vector, $u(k) \, \in \, \mathbb{R}^{n_u}$ the input variable and $\rho_k \, \in \, \mathbb{R}^{n_\rho}$ the scheduling vector. The explicit nonlinear term $\Phi \, : \, \mathbb{R}^{n_x} \, \to \, \mathbb{R}$ is structurally known and Lipschitz, as details Assumption 1. We note that $u$ and $\rho$ are, a priori, deterministic and known signals. Furthermore, we stress that $x(k)$ and $\theta(k)$ are independent from $\rho(k)$ and $u(k)$.

Assumption 1. The nonlinear term $\Phi (x(k))$ in the considered NLPV model of Eq (2.1) obeys a local Lipschitz condition around $x$, this is:

Remark 3. Assuming that the structure of the explicit nonlinearity of the NLPV is known is not at all absurd. Note that, in the majority of LPV cases, these nonlinearities are embedded to the scheduling parameters $\rho$. On the contrary, for the NLPV case, they are willingly made explicit (outside of $\rho$) because the designer knows how they are physically derived. In many processes, when two physical systems are cascaded (actuator and the actual process, for instance), and one of them has a known nonlinearity, this nonlinearity re-appears in the input/output relationship for the complete cascaded block. As an example, take a vehicle suspension system: the dynamics of the suspended vehicle chassis depend on the suspension control system (damper and spring), which inherently have known nonlinear behaviours; for the case of electro-rheological suspensions, this behaviour is a hyperbolic tangent function ^[33].

Assumption 2. The NLPV system in Eq (2.1) is corrupted by process noise $w_x(k) \, \in \, \mathbb{R}^{n_x}$ and measurement noise $w_y \, \in \, \mathbb{R}^{n_y}$. We consider that the process noise is Gaussian, whereas the measurement noise is a Gaussian mixture. This is: the process noise is a zero-mean white noise, i.e. $p(w_x) \, : \, \mathcal{N}(0, Q_x(k))$, where $p(\cdot)$ stands for the probability density function and $Q_x(k)$ for the co-variance matrix related to $w_x$, whilst the measurement noise $w_y$ is non-Gaussian, including the presence of outliers. This noise is distributed as follows:

which is a mixture of a nominal normal distribution $p_1(w_y) \, : \, \mathcal{N}(0, Q_y^1(k))$ and a contaminating probability density $p_2(w_y) \, : \, \mathcal{N}(0, Q_y^2(k))$, being $Q_y^1(k)$ and $Q_y^2(k)$ the co-variance matrices for the nominal and the contamination term. The degree of contamination parameter abides to: $0 \, < \, \epsilon \, < \, 1$.

Remark 4. The output measurement noise $w_y$, for $\epsilon \, = \, 0$, would have a Gaussian distribution $p_1(w_y)$, whereas, for $\epsilon \, = \, 1$, it would also present a normal distribution, of $p_2(w_y)$. In both cases, there would be a complete absence of outliers, which is not envied if the stochasticity of the process is to be taken into account.

Assumption 3. The change in the model parameters has a random walk behaviour, this is:

or, more generally,

where $p(w_\theta) \, : \, \mathcal{N}(0, Q_\theta(k))$, being $Q_\theta (k)$ the co-variance matrix for the random walk distribution. Matrix $C_\theta$ is a prior known and non-singular; it can be used to plug information on the model parameter change phenomenon that occurs.

Despite the fact that many system parameters can be available with some reasonable accuracy, given within a certain range, some model parameters may be entirely unknown because manufacturers consider these data as proprietary information. For example, the precise determination of leakage and friction coefficients causes great difficulty in the control of electro-rheological actuators ^[29]; accordingly loss of effectiveness faults must be included to reflect the variation of the physical system parameters. Due to the obvious presence of non-Gaussian noises, as well as the impossibility of determining the exact values of the system parameters, the observed system is considered as a stochastic system whose parameters are successfully determined by the proposed identification procedure.

Bearing in mind the previous Assumptions, the goal of this paper is reiterated: propose an algorithm to jointly estimate the states $x(k)$ and parameters $\theta _k$ of the NLPV system described by Eq (2.1), for all $k$. Moreover, asymptotic convergence of the estimation is required, despite the presence of non-Gaussian measurement noise (outliers) $w_y$. Output data $y(k)$, input data $u(k)$ and scheduling data $\rho (k)$ is considered available. Formally, this Problem is described:

Problem 1. For the NLPV system in Eq (2.1), consider an augmented states vector given by:

which compiles not only the states, but also the varying process model parameters. Consider $q_k(\cdot)$ as a nonlinear map at instant $k$. Then, find an estimation law of the following form:

with $\xi_k \, : \, \mathcal{N}\left(0, Q_\xi(k)\right)$, with $Q_\xi(k) \, = \, \mathit{\text{diag}}\{Q_x(k) \, , \, (1-\epsilon)Q_y^1(k) + \epsilon Q_y^2(k) \, , \, Q_\theta(k)\}$, encompassing all disturbances:

This estimation should minimise the sensitivity of $\hat{z}$ to $w_y$ and ensure that the expectation of estimation error $E((e(k)) \, = \, E\left(z(k) - \hat{z}(k)\right)$ has local asymptotic convergence for any state error starting condition $e_x(0)\, = \, x(0)-\hat{x}(0)$ given within a basin of attraction $\Omega$, which is equivalent to:

In the sequel, Problem 1 is solved by the use of a modified EMMF. More specifically, the Fischer information of the regular EMMF in the a posterior co-variance matrix is replaced by an approximation given by the derivative of Huber's function, which increases the flexibility of the algorithm and speeds up its convergence rate. The benefits of this modified version of the EMMF is discussed and illustrated in ^[40,41].

Remark 5. We note that the EMMF has been shown to possess offset-free convergent robust filtering qualities ^[41], which can be checked through a set of linear inequalities. If such inequalities are verified, the estimates from the EMMF, given in an augmented state form alike of $z(k)$, retain low sensitivity to measurement outliers.

Remark 6. The considered NLPV model presented in Eqs (2.1)-(2.4) has affine parameter-dependence w.r.t. the scheduling variable $\rho_k$. Anyhow, the method could directly consider polynomial, polytopic and fractional parameter dependencies, for which only the form of matrices $A(\cdot)$, $B(\cdot)$ and $C(\cdot)$ have to be altered. Moreover, we note that the proposed joint estimation problem, given in the recursive form of Eq (2.10), is set to benefit from the availability of the scheduling parameter data $\rho_k$, which coordinates the estimation function $q_k(\cdot)$, as shown in the sequel. In contrast to the nonlinear EMMF, the proposed NLPV EMMF uses an explicit $\rho_k$-dependent time-varying function $q_k(\cdot)$. As done in other LPV identification methods ^[14], the scheduling variable $\rho_k$ data is used to provide more data to the filter, which is potentially better than using a nonlinear EMMF that disregards the knowledge of $\rho_k$.

Lemma 1. Consider a joint state and parameters estimation law of the form given by Eq (2.10), in Problem 1. The convergence of the estimation error $e(k)$ can be analysed with an associated ODE system. If this associated system is locally asymptotically stable, which can be verified through the direct Lyapunov method, then, the expectation of the estimation error $E\left(e(k)\right)$ converges.

Proof. Provided in the original paper by L. Ljung ^[25].

3.   Robust estimation procedure for NLPV systems

In order to develop the proposed robust method for joint estimation of states and parameters of affine NLPV systems with explicit Lipchitz nonlinearities, given in the form of Eq (2.1), the full expression for $z(k)$, in terms of the nonlinear maps $q_k(\cdot)$ and $h_k(\cdot)$ are detailed:

Then, based on these expressions, the extended robust filter (EMMF) ^[40] is proposed for the NLPV case. The estimation of $z(k)$, namely $\hat{z}(k)$, addressing Problem 1, has the following form:

being implied that:

For the outputs contaminated with outliers, the nonlinear function $\psi (\cdot)$ issues a transformation upon $v(k)$, with is associated to the output prediction error. This nonlinearity is a Huber's influence function, as gives:

with its derivative with respect to $v_j(k)$ given by:

Remark 7. The sign function is defined as usual, being equal to identity for positive entries, null for null entries and minus identity for negative entries.

Remark 8. The use of the nonlinear Huber's influence function, given in Eq (3.13), is set to assign less weight to small quantities of larger residuals, such that the presence of outliers does not have significant influence on the final estimate $\hat{z}(k)$, while giving unit weight to the main part of the data population with moderate residuals.

The following block structured form of the EMMF estimation gain and co-variance matrices are given:

Finally, the estimation law from Eq (3.6) is given in a compact, recursive form, as follows:

which is directly equivalent to:

We note how Eq (3.18) explicitly expresses the state and parameter estimates w.r.t. previous scheduling parameters data, benefiting from its availability. In a full nonlinear formulation, as in the original EMMF, with unknown scheduling parameter data $\rho_k$, the nonlinear map $f_{k-1}(\cdot)$ is only written in terms of $\hat{x}(k-1)$, $u(k-1)$, and $\hat{\theta}_{k-1}$.

Remark 9. The initial conditions for the estimation algorithm are:

These values ($\hat{x}(0)$ and $\theta_0$) represent some a priori knowledge about the NLPV system variables. For the case when these are not known, they can be chosen, for simplicity, as null entries.

Remark 10. The main difference from the proposed joint state and parameter estimation algorithm is the use of Huber's nonlinearity $\Psi(\cdot)$. Its insertion suitably transforms the output prediction errors to penalize the presence of outliers and, thus, eliminate their influence upon the resulting estimation.

3.1.   Convergence analysis

In order verify that the proposed joint state and parameter estimation algorithm converges, i.e. $\lim_{k\to \infty} e(k) \to 0$, an associated ODE system is presented and its stability is verified using a direct Lyapunov method over the associated ODE from Lemma 1.

Lemma 2. The expectation of the joint parameter and state estimation error derived from proposed method in Eqs (3.18)-(3.19) is associated to an ODE system^*, which is given by:

^*With $e_x(\tau) \, = \, x(\tau) - \hat{x}(\tau)$ and $e_\theta (\tau) \, = \, \theta (\tau) - \hat{\theta}(\tau)$.

The nonlinear functions $f_x^e(\cdot)$ and $f_\theta^e(\cdot)$ derive from the difference between the Lipschitz affine NLPV system dynamics and the estimation from the algorithm, as follows:

Proof. The proof follows from ^[25]. Take a continuous version of the discrete time system where $\tau \, = \, k'T_s$, where $T_s$ is the sampling period of the NLPV system in Eq (2.1). Then, taking $k \, = \, \frac{k'}{T_s}$, the following continuous-time approximation is valid: $\frac{dr(\tau)}{d\tau} \, = \, r(k)-r(k-1)$.

Firstly, take the dynamics of the parameter estimation error, for which the parameter estimates are given by Eq (3.19), considering $w_\theta$ null:

It directly follows that:

which asymptotically converges to some constant $\overline{e_\theta}$ due to the construction of $M(\tau)$.

Then, take the dynamics of the state estimation error, for which the state estimates are given by Eq (3.18):

which, by exploiting the affine nature of $A$ and $B$, leads to:

Finally, since $e_\theta(\tau) \, = \, \overline{e_\theta}$ and taking $u(k-1) \, = \, \overline{u}$ as constant, one obtains:

This concludes the proof.

Lemma 3. The ODE system of Eq (3.21) is locally asymptotically stable for any real Lipschitz constant $\Gamma$ if there exists a positive definite matrix $P_L\, = \, P_L^T \, \succ 0$ and a scalar $\epsilon _\gamma \, \in \, \mathbb{R}$ such that the following LMI holds:

Then, for any starting condition $x(0)\, \in\, \Omega$, which is the sub-level ellipsoid associated to $P_L$, the proposed algorithm in Eq (2.10) yields local asymptotically stability.

Proof. In order verify the stability property of the ODE system of Eq (3.21), the direct Lyapunov method is used. Assuming there exists some positive definite matrix $P_L$, consider the following positive definite candidate Lyapunov function:

The local basin of attraction is given w.r.t. the state error estimation $e_x(k)$, since the parameter error $e_\theta(k)$ asymptotically converges to some constant $\overline{e_\theta}$ by construction (refer to the proof of Lemma 2). This basin of attraction can be estimated as sub-level set ellipsoid related to $P_L$, i.e.:

Differentiating $V$ w.r.t. $\tau$ along the solution of the ODE system, we get:

Note that the latter term $\frac{\partial \Psi(v(\tau))}{\partial\tau}$ can be expanded as follows:

Moreover, note that $v(\tau) \, = \, T(\tau)[C(\rho_k, \theta(\tau))x(\tau) - C(\rho_k, \hat{\theta}(\tau))\hat{x}(\tau)]$. Then, it holds for $\frac{dT(\tau)}{d\tau} \, = \, 0$ that:

Recall from previous development that $\frac{\partial \Psi (v(\tau))}{\partial v}$ is naturally bounded, i.e.:

which means that $\left(\frac{d\Psi (v(\tau))}{d\tau} \, \leq \, \frac{dv(\tau)}{d\tau}\right)$ holds. Thus:

Defining $\lambda (\tau) \, = \, \left(e_x(\tau)^T \, \left[\left(\Phi(x(\tau))-\Phi(\hat{x}(\tau))\right)\overline{u}\right]^T \, \Psi (v(\tau))^T \right)^T$, it follows that:

with:

The Lipschitz condition from Eq (2.5) holds as:

which is analogous to:

Since the direct Lyapunov method requires $\frac{dV(\tau)}{d\tau}$ to be nonpositive, the explicit nonlinearities can be overlapped with the aid of the Lipschitz condition, taking:

for $\epsilon _\gamma$ being a (scalar) slack variable.

Finally, the $\mathcal{S}$-procedure ^[7] can be applied to Eq (3.42), considering $\lambda (\tau)^TM_\Gamma \lambda (\tau) \, = E\left\{\lambda (\tau)^TM_\Gamma \lambda (\tau)\right\}$. By doing so, condition:

is equivalently expressed in the form of the LMI in Lemma 3. Therefore, if the positive definite $P_L$ and a scalar $\epsilon _\gamma$ indeed exists, the Lyapunov condition is verified and the proposed algorithm guarantees offset-free joint state and parameter estimation. This concludes the proof.

Remark 11. We note that inequality (3.43) is a potentially conservative sufficient condition that ensures the convergence of the NLPV joint estimation algorithm. Therefore, if there exists some positive definite $P_L$ which verifies LMI in Lemma 3, the local zone of attraction within which the ODE from Eq (3.21) displays asymptotic convergence may be small. The size of this zone depends on a series of matters, such as the model of the process and the Lipschitz constant $\Gamma$. In Section 4, we show an illustration of this attraction region for a suspension system example.

4.   Experimental validation results

In this Section, we present experimental validation results concerning the proposed joint parameter and state estimation, considering Lipschitz NLPV in form of Eq (2.1).

4.1.   Semi-Active suspension test-rig

The identified process is that of a semi-active suspension system. The used data is provided by the sensors of a $1/5$-scaled vehicle mechatronic test-bench, equipped with $4$ semi-active Electro-rheological (ER) dampers with the force range of $\pm \, 50 \, \rm{N}$. This vehicle test-rig, shown in Figure 1, allows testing under different configurations and use-cases (for full details, refer to www.gipsa-lab.fr/projet/inove). Each damper is controlled through a PWM signal given within the $[0 \, 0.35]$ range^†. Below each wheel of the vehicle lies an OMRO$N^{TM}$ linear servomotor which is able to mimic various road type conditions, with a maximal velocity threshold of $1.5 \, \rm{m/s}$.

^†These PWM modules operate at $25 \, \rm{kHz}$ rate duty-cycles and vary the electric field that is applied over each ER damper chambers and, thus, change the ER fluid flows viscosity and the delivered damping force.

4.2.   Lipschitz NLPV model

In order to describe the vertical dynamics of this experimental platform, a decoupled, quarter-car modelling framework is followed: each corner of the vehicle is represented by an individual set of equations. Each damper is controlled seeking to minimise the effects of the road upon the safety and comfort of the passengers on-board. The dynamic coupling effects between the corners are neglected. Through the sequel, the presented results are those from the front-left corner test-bench.

Figure 2 shows a schematic diagram of the quarter-car representation of a vehicle with four suspension units. Each semi-active suspension system comprises a spring with a stiffness parameter $k_s$ and a controlled damper of variable damping coefficient $c(\cdot)$. The chassis body at each side is represented by a sprung mass $m_s$, which is connected to the wheel-link, represented by the unsprung mass $m_{us}$. The wheel is represented by a spring with stiffness $k_t$. The vertical dynamics are described regarding the displacements of the sprung and unsprung masses: $z_s(t)$ and $z_{us}(t)$, respectively. The road profile is $z_r(t)$ a disturbance signal to this system.

The force provided by each ER damper is described through a dynamic version ^[33] of the widely-used nonlinear hyperbolic tangent model ^[19], as follows:

where $z_{def}(t) \, = \, \left(z_s(t) - z_{us}(t)\right)$ represents the suspension deflection and $u(t)$ stands for duty cycle of a PWM signal that regulates the ER damper. Regarding control purposes, this PWM variable $u(t)$ is the control input to the suspension system.

The spring and tire forces are given by:

The dynamics of the sprung and unsprung masses are obtained using regular laws of motion around an origin equilibrium:

Table 1 presents the model parameter values and their descriptions. The nominal values for spring stiffness, passive damper stiffness and viscous damping coefficient are provided by the manufacturers of these models. The remaining values were previously identified using regular state-of-the-art procedures ^[48].

We must stress that the damping characteristics of the ER damper vary over time. According to the lifespan of these components, the passive stiffness $k_0$ and viscous damping coefficient $c_0$ have variations, for multiple reasons ^[29,31], such as small oil leakages, air pressure inside the damper chamber, influence of external (high) temperatures, etc. Therefore, we assume that these parameters should be identified online and are subject to a random walk behaviour and additive noise, in the fashion of Eq (2.8). Accordingly, it is implied that:

where $\theta_0$ comprises the values for $k_0$ and $c_0$ given in Table 1 (nominal values).

The two major control objectives ^[34] of semi-active suspensions systems are: vehicle body isolation and passenger comfort enhancement. These two goals are physically conflicting: while stiff/high damping enhances passenger comfort, smooth/low damping enables easier road holding. An accurate knowledge of the actuator dynamics (damping force delivered by the ER damper) is necessary to correctly design control strategies regarding these two objectives. Accordingly, the online estimation of parameters $\theta$ is essential.

We note that the active hysteresis coefficients $k_1$, $c_1$ and the dynamic yield force of the ER damper $f_c$ appear only in the controlled part of the damper force $F_{ER}(t)$, see Eq (4.1). This means that possible variations upon these parameters (due to external factors, for instance) can be treated directly as actuator loss of effectiveness, as done in previous papers ^[29]. Therefore, their online estimation is not as important as the one regarding the passive terms, present in $F_d(t) - F_{ER}(t)$, which directly appear in the dynamics and cannot be indirectly accounted for.

Furthermore, as illustrates Table 1, we note that the tire stiffness and spring stiffness values are much larger that those for the damping coefficient, which means that the model is not so sensible to their variations. Changes upon the sprung mass and unsprung mass are usually accounted for by robustness, since any suspension control requires the tolerance of these parameters within a nominal interval ^[30], since these parameters vary according to the amount of vehicle passengers.

In order to cast this quarter-car suspension system dynamics into a Lipschitz NLPV state-space model, the system states are selected as follows:

The available measurements of this system are the vertical acceleration variables, given through on-boards sensors (inertial units/accelerometers), as follows:

Then, considering $w_y$ as coloured measurement noise and $\dot{z}_r(t)$ as the load disturbance variable, an NLPV state-space formulation is found with $\rho = u(t)$ as the scheduling variable:

where $\omega(t) \, = \, \left[\dot{z}_r(t) \, , \, w_y (t)\right]^T$ is a concise disturbance vector. The nonlinearity is:

with $\Gamma_{in} \, = \, \left[k_1, c_1, 0, -c_1, 0\right]$. Notice that Eq. (4.8) verifies a local Lipschitz condition in $x$, as expected by Assumption 1:

with $\Gamma \, = \, ||\Gamma_{in}||$ is the smallest Lipchitz constant that verifies Eq (4.9).

The state-space matrices $A_c$, $B_c$, $C_c$, $D_{c_1}$ and $D_{c_2}$ are given by:

This model is Euler-discretized with a sampling period of $T_s \, = \, 5 \, \rm{ms}$, which is an operational constraint of the considered test-bench. This leads to discrete-time Lipschitz NLPV system of the form in Eq (2.1). W.r.t. to the original problem setup model, the discretization implies $A(\cdot) \, = \, \mathbb{I}_{n_x} + T_sA_c(\cdot)$, $B(\cdot) \, = \, T_sB_c(\cdot)$ and $C(\cdot) \, = \, T_s C_c(\cdot)$.

Regarding the previous notation, the road profile derivative disturbance $\dot{z}_r(t)$ stands for $w_x(t)$. We note that in many modern cars, cameras (and adaptive estimation algorithms) are used to preview the future road profiles, e.g. ^[30,44,46]. This is trivial and widely seen in the automotive suspension literature ^[30]. Thus, from the viewpoint of the joint estimation algorithm, $w_x$ is a known variable.

4.3.   Convergence verification

It must be remarked that the considered experimental test-bench has physical limits to the motion variables ($z_s$, $\dot{z}_s$, etc.). Therefore, it is implied that the states $x$ have upper and lower limits; mathematically, this is expressed as a box-type set constraint $x \, \in \, \mathcal{X} \, \forall \, t$, where:

The numerical values for these upper and lower (physical) limits of the states are given below, in Table 2.

Before presenting the actual estimation results obtained with the proposed algorithm, its convergence property regarding state estimation must be analysed. In fact, state estimation convergence can be tested with the LMI provided in Lemma 3. For the considered system, the corresponding LMI is indeed verified, which ensures that the Lyapunov convergence condition holds. The symmetric positive definite matrix $P_L$ for which this LMI holds is the following, with $\epsilon _\gamma \ = \, 8.02 \, 10^{-8}$:

The sub-level ellipsoid set of $P_L$ determines the zone of attraction estimate for the proposed NLPV EMMF. This is, for any starting condition $x_0 \, \in \, \Omega$, convergence is verified. We note that the LMI problem from Lemma 3 is evaluated offline within $2.10\, \rm{s}$, using Matlab software and SDPT$3$ solver on an i5 CPU@2.4 GHz Macintosh with 8 GB of RAM.

Regarding this matter, Figure 3 shows this zone of attraction in three different $2$D cuts: on the left-side, the region of attraction for states $x_1$ and $x_2$ is given; in the center, the region for states $x_3$ and $x_4$ is shown; while on the right-side, the region for states $x_4$ and $x_5$ is exhibited. In this Figure, the (minimal and maximal) physical limits of each state are also displayed (those in Table 2). Clearly, for any physically possible $x_0$, the algorithm will ensure convergence, since all possible $x_0$ are indeed contained inside $\Omega$: $\mathcal{X} \subset \Omega$, which means that $\forall \, x_0 \, \left(x_0 \, \in \, \mathcal{X} \, \to \, x_0 \, \in \, \Omega \right)$ and, thus, the proposed solution is recursively feasible and ensures convergence when applied to the considered suspension system process.

4.4.   Joint state and parameter estimation

For the application of the proposed algorithm, a real input set comprising $t\, = \, 5 \, \rm{s}$ of vertical vehicle motions is considered (i.e. $k\, = \, 1000$ iterations). The road profile $z_r(t)$ represents a vehicle running in a straight line on a dry road when it encounters a sequence of $13 \, \rm{mm}$ bumps along its four wheels. This real road profile $z_r(t)$ applied to the test-bench through the controlled servomotors under the wheels and its derivative (input of the NLPV model) $w_x(t) = \dot{z}_r(t)$ are presented in Figure 4.

The test-bed model is then fed with a PRBS signal bounded within $[0 \, , \, 0.3]$, which is given in Figure 5. This input signal is of PRBS type in order to provide frequency-rich output signals $y$ (chassis body and wheel-link accelerations). These outputs, which are corrupted with coloured noise, are presented in Figure 6. As shown in this Figure, we note that the measurement noise is quite significant: the average noise-to-signal ratios for the chassis body and wheel-link acceleration outputs are of roughly $10\, \%$ and $5\, \%$, respectively.

The achieved estimation results are now presented. In order to show the advantages of the proposed method conceived for Lipschitz NLPV systems, it is compared against: i) an extended Kalman filter ^[47], which is tuned based on a tangent linerised model of the process around the origin equilibrium; and ii) a robust EMMF ^[41], which considers nonlinear process dynamics and does not make use of the scheduling variable and does not take into account the Lipschitz behaviour along the input trajectory. Through the sequel, NLPV-MEMMF denotes the results achieved with the proposed algorithm, EKF those obtained with a the extended Kalman filter and RMMF those for the robust nonlinear filter.

We note that the results comprise the application of these methods for different degrees of contamination $\epsilon$, as gives Eq (2.6). For presentation simplicity, the estimation results are shown only for $\epsilon = 0.05$. Anyhow, we assess the results considering other contamination levels through performance indexes.

The following results, considering the joint estimation of the NLPV system states $x$ and the variations of the damper force passive parameters $\theta_k$, were elaborated on a i5 CPU@2.4 GHz (2 Cores) Macintosh with 8 GB of RAM, with the aid of Matlab. The average execution time required to evaluate the recursive law given in Eq (2.10) is of $3.43\, \rm{ms}$ per iteration, which means that the proposed method can be used in online vehicle suspension applications, running on embedded on-board micro-controllers.

Considering a $\epsilon = 0.05$, Figure 7 shows the identification of the varying parameters over the $5 \, \rm{s}$ of data, while Figure 8 shows the estimation of the five system states. Clearly, convergent state estimation results and adequate parameter identification curves are obtained with the proposed method. These results could certainly serve for the computation of online state-feedback control laws, such as robust MPCs, for instance.

In order to quantify the performance of the proposed tool against the other methods, Figures 9, 10 and 11 show, respectively the Mean Square Error (MSE) index (in logarithmic scale), computed as the squared mean $2$-norm of the deviance between real variables $z$ and their estimations $\hat{z}$ (along the dataset), for different leves of contamination: $\epsilon = 0.05$, $0.1$ and $0.2$, respectively. Note that $z$ comprises states and parameters. Clearly, the smallest MSE results are obtained with the proposed NLPV-EMMF method, for all degrees of contamination, indicating the effectiveness of the proposed tool. Recall that, in the logarithmic scale, results closer to $-\infty$ indicate smaller estimation errors.

Complementary, to quantify the identification performance of the three methods regarding state and parameter estimations separately, the Robust Akaike criterion (RAC) ^[1] is considered. This criterion measures the inconsistency between the probability density of the estimation results against the exact probability density of the parameter variations (random walk). This criterion has a dynamic solution, as follows: RAC$(k) = \ln \left(\frac{1}{k-1}\left(\text{RAC}(k-1) + ||\theta_k - \hat{\theta}_k ||^2\right)\right)$.

Table 3 shows the means values obtained for the parameter estimation considering the RAC index (RAC) and the logarithmic Mean Square Error (MSE) for the state estimation; both indexes are computed for the $\epsilon = 0.05$ case. Evidently, the proposed NLPV-MEMMF achieves the best results. This is very interesting and means that the proposed solution can indeed serve for an online parameter and state estimation tool for controlled Lipschitz NLPV processes subject to inherent stochasticity, benefiting from the availability of the scheduling parameter data and the Lipschitz characteristic of $\Phi(x)$.

The presented results have shown that the widely-used extended Kalman filter (EKF) is very sensitive to the presence of non-Gaussian noises, as opposed to the proposed robust joint estimation algorithm (NLPV-MEMMF). Also, it is clearly shown that the included modifications to extended Masreliez-Martin filer (RMMF) for joint estimation increase the practical usability and convergence rate of the algorithm, which is now specially tendered for Lipschitz NLPV systems. In order to show robustness of the proposed robust joint estimation algorithm for systems with parameter faults with respect to these other conventional (and widely-used) joint estimation algorithms, the algorithms are tested over $1000$ random and independent simulations, for different contamination degrees. Regarding these tests, it is particularly important to note that NLPV-MEMMF maintains its high performances regarding the other methods, for all contamination degrees.

To conclude this Section, Table 4 provides the statistical data based on these random independent tests. Regarding this Table, for the cases of levels of high contamination degree ($\epsilon = 0.2$), it can be seen that the worst results obtained by NLPV-MEMMF are even better than the best result obtained by other methods. It can be clearly seen that the superiority of the proposed robust algorithm is greater in higher degrees of contamination, since the proposed method yields the least amount of variance in the considered performance index. Also, it is important to notice that, in presence of non-Gaussian noises, the proposed NLPV-MEMMF is an attractive alternative solution which outperforms other algorithms and, at the same time, has reduced complexity in comparison with them.

5.   Conclusions

This paper proposed a robust system identification algorithm specific conceived for Nonlinear Parameter-Varying models subject to stochasticity. The explicit nonlinear term obeys a Lipschitz condition and the matrices are affine on the scheduling parameters. The proposed approach, based on a modified extended Masreliez-Martin filter, jointly estimates the NLPV system states together with the physical parameters with precision. The estimation error asymptotically converges, despite the non-Gaussian disturbances. The effectiveness of the proposed algorithm is verified using real data from vehicle suspension mechatronic test-bed. As illustrated, the method can effectively estimate estates and stochastic parameter variations, becoming an interesting option for online state-feedback control policies, since it only uses output data. The major drawback of the proposed method is that it requires the explicit time-varying nonlinearity to be structurally known and Lipschitz, as well as needing a linear parametrisation of a priori evaluated basis functions. For future works, the Authors plan on comparing the method against artificial neural networks and genetic algorithms set for the structural discovery of the model dynamics, which do not require such assumptions.

Acknowledgments

The Authors thank Prof. Olivier Sename and his team at gipsa-lab for the application example (dataset and models). M. M. Morato is partially supported by ITEA3 European project $15016$ EMPHYSIS (Embedded Systems With Physical Models in the Production Code Software) and CNPq project $304032/2019-0$. V. Stojanovic has been supported by the Serbian Ministry of Education, Science and Technological Development ($451-03-68/2020-14/200108$).

Conflict of interest

The authors declare no potential conflict of interests.