- Email: [email protected]
road friction force estimation
ARTICLE IN PRESS
Engineering Applications of Artificial Intelligence 21 (2008) 442–456 www.elsevier.com/locate/engappai
Neural network based tire/road friction force estimation Jadranko Matusˇ ko, Ivan Petrovic´, Nedjeljko Peric´ Faculty of Electrical Engineering and Computing, University of Zagreb, Croatia Received 14 May 2006; received in revised form 1 March 2007; accepted 9 May 2007 Available online 2 July 2007
Abstract This paper deals with the problem of robust tire/road friction force estimation. Availability of actual value of the friction force generated in contact between the tire and the road has significant importance for active safety systems in modern cars, e.g. anti-lock brake systems, traction control systems, vehicle dynamic systems, etc. Since state estimators are usually based on the process model, they are sensitive to model inaccuracy. In this paper we propose a new neural network based estimation scheme, which makes friction force estimation insensitive to modelling inaccuracies. The neural network is added to the estimator in order to compensate effects of the friction model uncertainties to the estimation quality. An adaptation law for the neural network parameters is derived using Lyapunov stability analysis. The proposed state estimator provides accurate estimation of the tire/road friction force when friction characteristic is only approximately known or even completely unknown. Quality of the estimation is examined through simulation using one wheel friction model. Simulation results suggest very fast friction force estimation and compensation of the changes of the model parameters even when they vary in wide range. r 2007 Elsevier Ltd. All rights reserved. Keywords: Tire/road friction; Neural network; Estimation; Lyapunov stability; Antilock braking systems; Traction control systems; Lugre friction model; Model uncertainty estimation; Radial basis function networks
1. Introduction Active safety and comfortability of modern cars has become a major concern of many research activities. Above all, this is caused by the market demands for more sophisticated equipment in nowadays cars with no significant effect to the overall car price. As a consequence, these enhancements are primarily directed to better usability of the existing hardware. Almost every new car is equipped with sophisticated electronic systems such as the anti-lock brake system (ABS), traction control system (TCS) (also called ASR—anti slip regulation), electronic stability program (ESP), etc. ABS works by sensing slippage at the wheels during braking and adjusting braking pressure to ensure maximum contact force between the tires and the road with the purpose of not only shortening the braking distance but also of ensuring the driveability in critical situations (Kiencke and Nielsen, Corresponding author.
E-mail addresses: [email protected] (J. Matusˇ ko), [email protected] (I. Petrovic´), [email protected] (N. Peric´). 0952-1976/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.engappai.2007.05.001
2000). Similar to ABS, TCS also tries to maximize friction force between the tires and the road but during acceleration. This is done by controlling the engine power and the brakes. Both ABS and TCS are primarily concerned with control of the longitudinal (front-rear) motion of the vehicle. On the other hand an ESP system deals with the problem of the uncontrolled lateral vehicle motion (van Zanten et al., 1998). It uses brake actuators and engine torque to control yaw rate of the vehicle but also to limit the side slip angle, ensuring vehicle to follow the path intended by the driver as best as possible. All of the above-mentioned control systems have the purpose of efficient transmission of the forces from vehicle wheels to the road surface. Amount of the forces transferred in this contact is determined by the tire/road friction characteristics. Typical friction characteristics for different road conditions are shown in Fig. 1. Obviously, the goal is to maintain the operating point near the peak values of the tire/road friction characteristics. However, finding the optimal operating point is not a trivial task if the value of the tire/road friction force is not available. Therefore, friction estimation task plays important role in
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
443
presented in (Burckhardt, 1993): mðsÞ ¼ c1 ð1 ec2 s Þ c3 s.
Braking/traction force
Dry asphalt
The main drawback of this model is its nonlinearity in parameters ci , which makes the parameter estimation process difficult. To overcome this problem a rational function is used in Kiencke and Daiss (1994) to describe m2s relationship:
Wet asphalt
mðsÞ ¼ Snow
Slip Fig. 1. Typical tire/road static curves for different road conditions.
the design of the automotive control systems. Beside estimation accuracy, it is important to ensure high estimation availability. This means that estimator should provide useful information not only for various road conditions, but also in situations when friction modelling uncertainties are present. In literature one can find many various approaches to the tire/road friction modelling and estimation. The short overviews are given in Sections 2 and 3, respectively. Neural network based friction force estimator, which we propose as a valuable alternative to other approaches, is presented in Section 4. The key idea is to use a neural network, as the universal approximator, to compensate friction model inaccuracies. It has been added to the state estimator, which is based on nominal friction model, to provide robust estimation. Neural network parameters are adapted online with the training algorithm based on Lyapunov technique. Results of simulation experiments with our estimator is given in Section 5. The paper ends with the conclusion remarks in Section 6. 2. Overview of the tire/road friction modelling approaches Building an accurate tire/road friction model is important for the following reasons: (i) simulation of the overall vehicle behavior, (ii) design of the automotive control systems and (iii) estimation of the friction force. Historically, tire/road friction effects were described by the static models, either physical or empirical. The static friction models usually describe relationship between friction force F f (or more often actual friction coefficient m ¼ F f =F n , where F n is normal force in the tire/road contact) and dimensionless quantity s called slip using algebraic equations. The quantity slip s is usually defined as ro v , (1) s¼ maxðro; vÞ where o is wheel angular speed, v is vehicle speed and r is wheel radius. Simple static tire/road friction model is
(2)
ks . as2 þ bs þ 1
(3)
This model is primarily used for estimation purposes since it is linear in parameters ðk; a; bÞ and thus classical recursive least square method can be used as estimation algorithm. However, the most widely used static model is Pacejka model, also called ‘‘Magic formula’’ (Bakker et al., 1987). The model uses trigonometric functions to describe relationship between slip and friction force generated in the contact between the tire and the road: F f ðsÞ ¼ D sinðC arc tanðB s EðB s arc tanðB sÞÞÞÞ, (4) where B; C; D; E are model parameters. It is proven that Pacejka model has ability to approximate experimental data (Bakker et al., 1987). It can be used for both simulation and estimation tasks. However, as all other static models it has significant disadvantage related to inability to describe low slip effects. To resolve this problem a deeper insight is needed into physical effects taking place in the tire/road contact surface. The best insight in what is happening in the tire/ road contact surface is given by so called ‘‘brush’’ models (van Zanten et al., 1989). These models assume that contact between the tire and the road is realized through many tiny and massless elastic elements, called bristles. Bristle deformations are direct cause of the generated friction force. When a bristle enters into contact with the road, it starts to deform with deformation rate proportional to the relative velocity (adhesion region in Fig. 2). After reaching deformation limit, there is no additional deformation and bristle starts to slide over the contact patch (sliding region in Fig. 2). The original brush model describes the deformation of a bristle as follows: ( vr ; zðzÞozm ; z_ðzÞ ¼ (5) 0; zðzÞXzm ; where z is the bristle deformation, vr ¼ ro v is relative velocity and z is a spatial variable related to the bristle position while passing through the contact patch. Friction force generated by a bristle is given by F ðzÞ ¼ s0 zðzÞ þ s1 z_ðzÞ,
(6)
where s0 and s1 are the rubber longitudinal lumped stiffness and damping, respectively. Overall friction force generated in the tire/road contact is obtained by integration of the bristle forces over the whole contact patch. While saturation function (5) is used to describe deformation of a bristle in the original brush model, the
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
444
=
Ff Fn
confirmed that overshoot exists in the friction force at transition from the adhesion regime to the sliding regime of the bristles in the contact patch, but LuGre model cannot capture it. Also, LuGre model predicts the same initial slope of the static characteristics (at s ¼ 0) for all road conditions (all values of coefficient y), which is in conflict with experimental results in Dieckmann (1992). In spite of these possible shortcomings LuGre model is still the most widely used dynamic tire/road friction model and therefore is used in our work as the basis for development of the tire/ road friction force estimators. We expect that it will be straightforward to adapt our estimator to possible future changes of the LuGre model.
Point of maximum friction
Sliding region
Adhesion region
Slip s Fig. 2. Regions on the tire/road static friction curve.
3. Overview of the tire/road friction estimation approaches
first-order nonlinear function is used in LuGre model (Canudas de Wit et al., 1995; Canudas de Wit and Tsiotras, 1999): s0 jvr j z_ðzÞ ¼ vr zðzÞ, (7) gðvr Þ where gðvr Þ is Stribeck function defined as pffiffiffiffiffiffiffi gðvr Þ ¼ mC þ ðmS mC Þ e vr =vs .
(8)
In (8) mC denotes normalized static friction, mS normalized maximum friction and vs Stribeck velocity. Although LuGre model (6)–(8) has compact form, it is primarily used for simulation purposes, since spatial integration over the contact patch is computationally very demanding. For estimation as well as control system design purposes a less complex model needs to be used. Therefore, instead of describing dynamics of each particular bristle deformation, an average bristle deflection is introduced in (Canudas de Wit and Tsiotras, 1999): Z 1 L z¯ ¼ zðzÞ dz, (9) L 0 where L is the length of the contact patch. Usage of the average value of the bristles deformation led to the lumped LuGre model given by s jv j _z¯ ¼ vr y 0 r z¯ , (10) gðvr Þ where the coefficient of road adhesion y is introduced in order to model friction characteristic changes caused by the road condition changes. Lumped LuGre model (10) has a steady-state error, which is eliminated in Deur (2001) and Deur et al. (2004) by introducing an additional term s0 jvr j k (11) þ rjoj z¯ , z_¯ ¼ vr y gðvr Þ L where k is a constant. Some researchers criticized LuGre model as not being able to capture all dynamic effects in the contact patch. For example, in Bozˇic´ et al. (2002) it is experimentally
Estimation of the friction force generated in the contact patch between the tire and the road is a quite demanding task. Special hardware sensors described in Hozwarth and Eichhorn (1993) are considered as inadequate due to their high prices. More appropriate solution is to use readily accessible, easily measurable, signals to estimate the tire/ road friction in an indirect way. In literature many different approaches to the problem of tire/road friction estimation have been proposed. Generally, all approaches are directed either to the estimation of actual friction force (or friction coefficient) or just to the qualitative estimation of road conditions. While later approaches are simpler to develop and implement, they do not enable applications of the advanced ABS and TCS that search for optimal values of the slip s (point of maximum friction, see Fig. 2). Instead, the designer has to predefine referent slip for all road conditions based on the experience gained through experimental tests. On the other hand, if the estimator estimates actual friction force, and if the slip value is known, it is possible to implement advanced ABS and TCS. A road condition estimator based on Kalman filter is proposed in Gustafsson (1997). The approach supposes that the slope of the friction force characteristic at low slip depends on the road conditions (see Fig. 3), which is experimentally confirmed in Dieckmann (1992). The appropriate model which describes the tire/road friction in low slip region is given by the following equation: 1=k sðtÞ ¼ ðmðtÞ1Þ þ eðtÞ, (12) d where k is low-slip slope of the friction characteristic, d is the slip offset and eðtÞ is the white noise signal (see Fig. 3). In order to ensure good tracking performance this estimator has an additional subsystem for detection of abrupt changes of the road conditions. In Canudas de Wit and Horowitz (1999), a state estimator based on passivity theory is presented. It uses standard dynamic LuGre model for description of the tire/ road friction. This estimator is able to estimate all state
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
sph a
lt
Ff FN
Dry a
=
ow
Sn
1
2
Slip s
Fig. 3. Slopes of the m2s curves near zero slip for different road conditions.
variables as well as the coefficient of the road adhesion y, which represents road conditions (dry asphalt, wet asphalt, snow, ice, etc.). Convergence of the estimation error to zero is guaranteed by the adaptation law derived applying Lyapunov stability theory. However, proposed estimator is sensitive to the model inaccuracies and measurement noise. Another approach to the simultaneous estimation of the state variables and the coefficient of the road adhesion is presented in Matusˇ ko et al. (2002). The estimator is based on modified LuGre dynamic friction model and uses extended Kalman filter, additionally augmented by an integral term (PI Kalman filter) in order to capture variations of the road adhesion coefficient. Simulation results suggest that the estimator is appropriate when significant level of the measurement noise is present. However the approach relies on exact tire/road friction model and thus uncertainties in the friction model parameters may significantly degrade estimation accuracy. In Pal et al. (1994) a neural network has been applied for the road condition estimation. The network parameters were tuned by an off-line training algorithm applied on the field test data acquired for different road conditions. The approach suffers from two big drawbacks: (1) the output of the estimator can get only several discrete values meaning dry, wet, snow or icy road, but cannot estimate friction force or state variables, and (2) long series of field test data are required to train the neural network. A somewhat different approach is presented in Umeno et al. (2002), where spectral analysis of the wheel speed signal is used to determine the road conditions. The approach is based on the tire vibration model, instead of the friction model, and relies on the assumption that the road excites the tire/road vibrational system with the strength dependant on the road conditions. As a result of this excitation higher frequency harmonic (around 40 Hz) occurs in signals measured by wheel speed sensors. The main
445
advantage of this estimator is its ability to estimate the road condition even at zero slip (free rolling motion). The necessity of high-resolution wheel speed sensors and problems due to weak excitations on icy roads are significant drawbacks of the method. In the introductory section, it was pointed out that the estimator should posses high availability. The majority of the estimators are based on some type of the tire/road friction model, where usual assumption is that model accurately describes effects in the tire/road contact patch. However, this is only partially fulfilled, since models used for estimation purposes are rather simple, and thus modelling errors are inherent to them. Besides, parameters of the processes in the contact patch may vary with time and initially accurate model could drift away from the true one. Therefore, it is important to ensure robustness of the estimation to the model uncertainties. 4. Neural network estimator design As aforementioned, vehicle behavior on the road is dominantly determined by the amount of forces transferred between vehicle wheels and road. Generally, tire/road friction forces at the vehicle wheels may mutually differ and therefore it is important to estimate them separately for each wheel. This is the reason why one-wheel friction model (Fig. 4) is used in this paper as the basis for design of the tire/road friction estimator. One-wheel friction model is described by the following set of equations: m_v ¼ F f þ F n s2 vr ,
(13)
_ ¼ rF f þ ur so o, Jo
(14)
where m is the wheel mass, J is the moment of inertia of the wheel, ur is the torque acting on wheel axis, s2 is the viscous relative damping and so is the coefficient of viscous friction. Dynamics of the friction force F f in Eqs. (13) and (14) is modelled using lumped LuGre dynamic friction model described in Section 2. For the sake of clarity the model equations are again written here, but in somewhat different form (further on variable z is used instead of variable z¯ to denote average bristle deflection): F f ¼ F n ðs0 z þ s1 z_Þ,
(15)
s0 jvr j k þ rjoj z ¼ vr f ðvr Þ z f o ðoÞ z, z_ ¼ vr y gðvr Þ L (16) ω
v r
Fn
Ff
Fig. 4. One wheel friction model.
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
446
gðvr Þ ¼ mC þ ðmS mC Þ e
pffiffiffiffiffiffiffiffiffi
jvr =vs j
,
(17)
where f ðvr Þ ¼ yðs0 jvr j=gðvr ÞÞ and f o ðoÞ ¼ ðk=LÞrjoj. Although LuGre dynamic friction model can satisfactorily describe dynamics of the friction at the tire/road contact patch, many field test data are usually required to parameterize it accurately. In this paper, it is assumed that functions f ðvr Þ and f o ðoÞ in Eq. (16) are only partially known and the aim is to design friction force estimator which provides accurate estimation of the friction force under these assumptions. The nonlinear function f ðvr Þ can be written as the sum of a known function f n ðvr Þ (nominal value of function f ðvr Þ), and an unknown function Df ðvr Þ: s0 jvr j s0 jvr j f ðvr Þ ¼ y þD y ¼ f n ðvr Þ þ Df ðvr Þ. (18) gðvr Þ n gðvr Þ In order to compensate the unknown nonlinear function Df ðvr Þ in Eq. (18) a neural network is employed. Neural networks commonly consist of big number of highly interconnected neurons organized in two or more layers (McCulloch and Pitts, 1943). These interconnections give them universal approximation capability, i.e. a neural network can approximate arbitrary continuous nonlinear function with desired accuracy (Chen et al., 1995; Chen and Chen, 1995; Park and Sandberg, 1991) if consisted of at least two layers. Generally, a two-layer neural network with linear output layer can be described by the following equation: yNN ¼ W2 wðW1 ; uÞ, n
can be represented with a scaling factor b. If its optimal value bn is found, then f o ðoÞ can be written as k (21) f o ðoÞ ¼ bn rjoj ¼ bn f on ðoÞ, L where f on ðoÞ is the function f o ðoÞ with the nominal value of parameters k and L. The task of the proposed estimator is to estimate optimal values of the network parameters Wn2 and parameter bn . As the first step in the estimator design procedure the following state variable transformation was performed (Canudas de Wit and Horowitz, 1999): x1 ¼ rF n s1 z þ Jo, x2 ¼ z.
One-wheel friction model (13)–(17) in the space of new variables can be written as x_ ¼ Ax þ Bf ðy; uÞDx þ Bf o ðyÞDx þ Ru þ Ey, y ¼ CT x, " x¼
x1 x2
2 s0 A ¼ 4 s1
# ;
" D ¼ ½ 0 1 ;
R¼
3 1 6 J 7 7 C¼6 4 rF n s1 5; J
3
1
0
0
1
u¼
v ur
"
5;
0
"
# ;
0
#
; 1 3 2 s0 J so 5 ; E ¼ 4 s1 r
B¼
# ;
y ¼ ½ o :
For the system model given by Eq. (23) the state estimator in standard predictor–corrector form (Fig. 5) is given by +
y(t)
E
u(t)
−
K
A
1 s
R
B
ˆ x(t)
C
D
X
(20)
Uncertainty in function f o ðoÞ is the consequence of variations of the parameters k and L. Parameter k is lumped coefficient of the LuGre model and varies between 1.1 and 1.4 depending on the tire type (Deur et al., 2004). On the other hand the length of the contact patch L varies with variations of the pressure in the tire. These variations
0
0
2
m
f ðvr Þ ¼ f n ðvr Þ þ Df ðvr Þ f n ðvr Þ þ Wn2 wðvr Þ.
ð23Þ
where
(19)
where u 2 R is network input, y 2 R is network output, w: Rn 7!Rm is nonlinear kernel function and W1 and W2 are parameter matrices of hidden and output network layer, respectively. Various neural network structures are obtained by choosing different kernel functions. As mentioned above, friction characteristic may vary with time. This implies that neural network parameters (weights) should be updated in an on-line manner. Hereafter we suppose the usage of a two-layer neural network with hidden layer weights W1 set off-line, and linear output layer weights W2 updated on-line. The reasons for choosing such a network architecture are computational simplicity of the on-line training algorithm and satisfactory approximation of function Df ðvr Þ. If parameters W2 are adjusted in such a way that their optimal values Wn2 are found, then neural network of proper structure gives good approximation of the function Df ðvr Þ, and Eq. (18) can be written as
ð22Þ
fn (y,u) ˆ
fωn (y) ˆ ) Δf (y,u,W 2
Fig. 5. State estimator principle scheme.
ˆ y(t)
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
^ ^ bn Bf on ðyÞDx bBf on ðyÞDx ~ ^ ¼ bn Bf ðyÞDx~ þ bBf ðyÞDx,
the following equations: ^ ^ 2 wðy; uÞDx^ ^ þ BW x_^ ¼ Ax^ þ Bf n ðy; uÞDx^ þ Bbf on ðyÞDx ^ þ Ru þ Ey þ Kðy yÞ, T
^ y^ ¼ C x.
on
ð24Þ
The designer has the ability to influence performances of the estimator (24) by tuning the value of the gain vector ^ 2 of K ¼ ½k1 k2 T and adapting the output layer weights W ^ In order to guarantee the neural network and parameter b. convergence of the estimation, the adaptation algorithm of ^ 2 and parameter b^ as well as rules the network weights W for the vector K selection are derived applying Lyapunov stability analysis. The results of the analysis are summarized in Assumption 1 and Theorem 1. Assumption 1. Vehicle velocity v and wheel angular velocity o are bounded, i.e.
~ 2 wðy; uÞDx^ x_~ ¼ ðA KCT Þx~ þ B½Wn2 wðy; uÞDx~ þ W n ~ ^ þ f ðy; uÞDx~ þ b f ðyÞDx~ þ bf ðyÞDx, n
on
on
~ y~ ¼ CT x,
ð31Þ
^ ^ 2 and b~ ¼ b b. ~ 2 ¼W W where W 2 Let Lyapunov function candidate be given as n
n
~ T PW W ~ 2 g þ gb~ 2 , V ¼ x~ T Px~ þ trfW 2
(32)
where P is a symmetric positive definite matrix, and trfg is matrix trace operator defined as (Horn and Johnson, 1999) n X
xi;i .
(33)
i¼1
ð25Þ
It is important to note that Assumption 1 is introduced only due to theoretical reasons, since in reality quantities v and o are always bounded. Theorem 1. Let the system and state estimator be given by Eqs. (23) and (24), respectively, and let the following conditions for the gain vector K be satisfied: k1 4rF N s1 k2 , k2 o0.
According to the Lyapunov theorem of stability (Slotine and Li, 1990), the time derivative of the Lyapunov function should be negative semidefinite in order to ensure system stability: _~ ~_ 2 g þ 2gb~ bo0. ~ T PW W V_ ¼ x_~ T Px~ þ x~ T Px_~ þ 2 trfW 2 V_ ¼ x~ T ððA KCT ÞT P þ PðA KCT ÞÞx~ þ 2x~ T PBWn2 wðy; uÞDx~ ~ 2 wðy; uÞDx^ þ 2x~ T PBW
ð26Þ
_^ 1 ~ x^ T DT wT ðy; uÞ, W 2 ¼ PW y
(27)
_ ^ b^ ¼ g1 y~ f on ðyÞDx,
(28)
þ 2x~ T PBf n ðy; uÞDx~ þ 2bn x~ T PBf on ðyÞDx~ _^ ~ _^ ~ T PW W þ 2b~ x~ T PBf on ðyÞDx^ 2trfW 2 g 2gbb. 2
Gain vector K ¼ ½k1 k2 is chosen in such a way that estimation error dynamics (31) is passive (Slotine and Li, 1990), i.e. transfer function GðsÞ ¼ CT ðsI A þ KCT Þ1 B
Proof of the Theorem 1. Estimation error dynamics is obtained by subtracting Eqs. (24) from Eqs. (23):
GðsÞ ¼
on
T
~ y~ ¼ C x,
ð29Þ
^ Taking into account where x~ ¼ x x^ and y~ ¼ y y. Eqs. (18), (20) and (21), and following equalities: ^ 2 wðy; uÞDx^ ¼ BW wðy; uÞDx~ BW2 wðy; uÞDx BW 2 ~ 2 wðy; uÞDx, ^ þ BW n
ð35Þ
T
where PW is a symmetric positive definite matrix and g is a positive scalar coefficient. With the adaptation laws (27) and (28), and condition (26), the convergence of the state estimates toward their true values is ensured.
x_~ ¼ ðA KCT Þx~ þ Bf ðy; uÞDx Bf n ðy; uÞDx^ þ Bf o ðyÞDx ^ ^ 2 wðy; uÞDx, ^ Bbf ðyÞDx^ BW
(34)
By inserting Eq. (31) into (34) it is obtained:
Let the adaptation laws of the neural network output-layer ^ 2 and parameter b^ be given by the following parameters W equations:
n
ð30Þ
on
equations describing estimation error dynamics can be written as
trfXg ¼
jvjovmax , jojoomax .
447
(36)
is strictly positive real function. By inserting (23) into (36) the transfer function of the estimation error dynamics becomes rF N s1 =J s þ rF N s0 =J þ sðs0 =s1 þ k1 =J rF N s1 k2 =JÞ rF N s0 k2 =J b1 s þ b0 ¼ 2 . ð37Þ s þ a1 s þ a0 s2
Sufficient conditions for strict positiveness of the transfer function (37) are a0 ; a1 40 ) k2 o0 ^ k1 4rF N s1 k2 Js0 =s1 , b0 a0 þ ðb1 a1 b0 Þo2 40 ) k2 o0 ^ k1 4rF N s1 k2 ,
ð38Þ
i.e. k2 o0 ^ k1 4rF N s1 k2 .
(39)
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
448
If the system is passive, then, from Kalman–Yakubovich lemma (Slotine and Li, 1990), it follows: 8Q; 9P: ðA KCT ÞT P þ PT ðA KCT Þ ¼ Q, PB ¼ C,
ð40Þ
where Q is a symmetric positive definite matrix. By taking into account Eq. (40) time derivative of the Lyapunov function (35) becomes ~ 2 wðy; uÞDx^ V_ ¼ x~ T Qx~ þ 2x~ T CWn2 wðy; uÞDx~ þ 2~yT W þ 2x~ T Cf n ðy; uÞDx~ þ 2bn x~ T Cf on ðyÞDx~ þ 2b~ x~ T Cf ðyÞDx^ on
_^ ~ T PW W ^_ 2 g 2gb~ b. 2 trfW 2
ð41Þ
Using the following feature of the matrix trace operator trfg: xT y ¼ trfxyT g;
x; y 2 Rn ,
(42)
Eq. (41) can be written in the form V_ ¼ x~ T Qx~ þ 2x~ T CWn2 wðy; uÞDx~ þ 2x~ T Cf n ðy; uÞDx~ _^ ~ y~ T f ðyÞDx^ gbÞ þ 2bn x~ T Cf on ðyÞDx~ þ 2bð on _^ ~ T ðy~ x^ T DT wT ðy; uÞ PW W þ 2 trfW ð43Þ 2 Þg. 2 If the adaptation laws of the neural network output _^ _^ weights W 2 and parameter b are chosen as _^ 1 ~ x^ T DT wT ðy; uÞ, W 2 ¼ PW y
(44)
_ ^ b^ ¼ g1 y~ T f on ðyÞDx,
(45)
then time derivative of the Lyapunov function becomes V_ ¼ x~ T Qx~ þ 2x~ T CWn2 wðy; uÞDx~ þ 2x~ T Cf n ðy; uÞDx~ ~ þ 2bn x~ T Cf on ðyÞDx.
ð46Þ
The following assessment is valid for Eq. (46): pffiffiffiffiffiffi ~ 2 þ 2 N h kCkF kWn2 kF kDkF kxk ~ 2 V_ p lmin ðQÞ kxk ~ 2 þ 2 maxff n ðy; uÞgkCkF kDkF kxk ~ 2, þ 2jbn j maxff on ðyÞgkCkF kDkF kxk
ð47Þ
where lmin ðQÞ is the least eigenvalue of matrix Q, k k is the Euclidian vector norm, k kF is the Frobenius matrix norm, N h is the number of the neurons in network hidden layer and maxff n ðy; uÞg and maxff o ðyÞg are maximum values of the nonlinear functions f n ðy; uÞ and f o ðyÞ given by maxff n ðy; uÞg ¼ y maxff o ðyÞg ¼
s0 jvr max j , gðvr max Þ
k rjomax j. L
ð48Þ
Kalman–Yakubovich lemma, given by Eq. (40), enables that matrix Q is chosen in such a way that following
inequality is satisfied: pffiffiffiffiffiffiffi lmin ðQÞ42 N h kCkF kWn2 kF kDkF þ 2 maxff n ðy; uÞgkCkF kDkF þ 2jbn j maxff on ðyÞgkCkF kDkF ,
ð49Þ
and that is the condition which ensures negative semidefiniteness of the time derivative of the Lyapunov function V_ p0. Consequently, convergence of the state estimates toward their true values is guaranteed. & Remark 1. It is obvious from Eq. (18) that function f ðy; uÞ actually depends only on relative velocity vr ¼ ro v, and therefore vr is the only needed input to the neural network. In that way the choice of the hidden layer parameters becomes very simple, and the needed number of hidden neurons becomes relatively small, which are both very beneficial for practical application. Remark 2. The estimator design procedure could also be performed applying a neural network to approximate the whole nonlinearity in (16), i.e. without splitting it in two separate nonlinear functions f and f o . In that case functions f and f o become f ¼ ðys0 jvr j=gðvr ÞÞ þ ðk=LÞrjoj and f o ¼ 0, respectively. Derivation of the selection rules for the vector K and the neural network weights adaptation law remains the same, and final results of the Lyapunov analysis summarized in Theorem 1 stay valid too. One needs only to set the parameter b equal to zero in all equations in which it appears and to exclude its adaptation equation (28) from Theorem 1.However, if a neural network is used to approximate the whole nonlinearity in (16), the input vector of the network x1 must be twodimensional, either x1 ¼ ½o; v or x1 ¼ ½o; vr . Of course, in that case an advanced algorithm for off-line setting of the hidden layer parameters, and increased number of hidden neurons, need to be used. Remark 3. Although we prove convergence of the state estimates for neural networks with fixed hidden layer parameters, it can also be proven for cases when hidden network parameters are trained on-line. For example, using Lyapunov function as in Baric´ et al. (2005) instead of (32), the convergence can be proven for multi-layer perceptron networks with on-line training of all its parameters. However, one has to be aware of increased computational complexity of the training algorithm. 5. Simulation tests Proposed friction force estimator given by Eqs. (24), (27) and (28) and the one-wheel friction model given by Eqs. (13)–(17) have been implemented in Matlab/Simulink environment and numerous simulation experiments have been performed in order to validate estimator performances. The model parameters used in simulations are given in Table 1. The architecture of the neural network used in simulations, selection of the gain vector K,
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456 Table 1 Model parameters used in simulations Values
s0 r s1 m s2 J Fn mC mS k vS L
200 (N/m) 0.25 (m) 0.4948 (N s/m) 350 (kg) 0.0018 (N s/m) 0.2344 ðkg m2 Þ 1000 ðkg m2 =s2 Þ 0.7 1.9 1.1 3.5 (m/s) 0.2 (m)
x1
W1 − w11
+
− w12 • • •
+
•
v1
y1
x2
2
v11
y11
x21
•
2
v12
y12
x22
− w1Nh
• • • 2
v1Nh
• • •
y1Nh x2Nh
0.2
−15
−10
−5 0 5 10 Relative velocity vr [m/s]
15
20
Then, spreads of the neurons were chosen using P-nearest ðP ¼ 2Þ neighbor heuristics (Moody and Darken, 1989): vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u P u1 X kwj wi k, (51) s1;i ¼ t P j¼1
W2
+
w22 w2Nh
y2 = yNN
giving values: s1 ¼ ½2:291; 1:732; 1:225; 1; 0:922; 0:742; 0:592; 0:5; . . . , 0:5; 0:592; 0:742; 0:922; 1; 1:225; 1:732; 2:291.
w2(N
h+1)
1 matrix P1 W and coefficient g , and simulation results are presented in succession.
5.1. Neural network used in simulations Neural network used in simulations was a RBF network (Fig. 6) described by the following set of equations:
y1i ¼ eðv1i =s1i Þ ;
0.4
Fig. 7. Set of RBF hidden layer activation functions.
Fig. 6. Two-layer RBF neural network.
2
0.6
−20
1
v1i ¼ kx1 w1i k2 ;
0.8
0
w21
• • •
•
Activation functions
1
Parameters
+
449
i ¼ 1 . . . N h,
i ¼ 1 . . . N h,
T
x2 ¼ ½y1 1 , y2 ¼ W2 x2 , yNN ¼ y2 ,
ð50Þ
where x1 ¼ vr is input to the network, v1 is the input vector to Gaussian activation functions of the hidden layer neurons, N h is the number of hidden layer neurons, W1 and s1 are vectors of centers and spreads of hidden layer neurons, respectively, and x2 and W2 are input and weight vectors of the output layer, respectively. Network hidden layer parameters were set off-line and fixed. N h ¼ 16 neurons were used with centers set to the following values: W1 ¼ ½10; 6; 3:5; 2:5; 1:5; 0:8; 0:4; 0:1; . . . , 0:1; 0:4; 0:8; 1:5; 2:5; 3:5; 6; 10.
With these values of W1 and s1 , the RBF network input space is well covered with activation functions in the range of vr 2 ð10; 10Þ (see Fig. 7), which is satisfactory for very wide range of vehicle velocities (up to 180 km/h). Of course, input range of the network can be easily extended in order to cover wider range if necessary. Hidden layer neurons are set more densely for smaller values of relative velocity because nonlinear effects in approximated function Df ðvr Þ are much stronger for small than for large values of vr , what can be easily checked analyzing expressions (17) and (18). Moreover, ABS and TCS usually keep the tire slip at values jsjo0:1, which means that relative velocities are in the range ð5; 5Þ. Higher slip can rarely occur, e.g. when road conditions abruptly change in large amount. Initial values of the neural network output layer weights W2 can be conveniently set to zero, and parameter b to 1. These parameters are then adapted on-line by adaptation laws (27) and (28). 5.2. Selection of the gain vector K, matrix P1 W and coefficient g1 According to Theorem 1, passivity of the estimation error dynamics is ensured if the gain vector K ¼ ½k1 k2 T is selected respecting conditions (26). Passivity area of the estimation error for model parameters given in Table 1 is shown in Fig. 8. By choosing proper values of gains k1 and k2 inside the passivity area desired dynamics of the estimation error with respect to the friction model uncertainties can be obtained. This dynamics is given by transfer function (37), which defines behavior of the ~ and also output error y, ~ when errors estimation error x,
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
450
in friction model occur. Since compensation of the friction model uncertainties is performed by neural network training according to (27) and (28), which use y~ as their input, it is desirable that dynamics given by (37) is as fast as possible, what can be achieved by proper selection of gains k1 and k2 . However, too fast dynamics can result in oversensitivity to measurement noise. Therefore, a tradeoff is necessary. It is suggested that dynamics of (37) is chosen several times faster than desired dynamics of estimation error compensation.
k1
Area of passivity
2000
Dynamics of the transfer function (37) can be adjusted by proper placement of its poles:
p1;2 ¼
a1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a21 4a0 2
.
(52)
Both poles must satisfy characteristic equation of transfer function (37): p2i þ pi a1 þ a0 ¼ 0;
i ¼ 1; 2,
(53)
wherefrom, by substitution of a0 and a1 with corresponding physical variables, the following equations are obtained: p21 þ p1
s0 k1 rF N k2 ðp1 s1 þ s0 Þ ¼ 0, þ p1 s1 J J
(54)
p22 þ p2
s0 k1 rF N k2 ðp2 s1 þ s0 Þ ¼ 0. þ p2 s1 J J
(55)
1000 {k1,k2} = {0,−0.4} k2
1
F
k1
=r
−1000
k2 σ1 n
Solution of Eqs. (54) and (55) yields the following expressions for k1 and k2 :
−2000
Fig. 8. Passivity area of the estimation in k1 , k2 plane.
k1 ¼
c2 d 0 c0 d 2 , d 2 c1 d 1 c2
(56)
k2 ¼
c1 d 0 c0 d 1 , d 1 c2 c1 d 2
(57)
7
2
6
1.9
5 1.8
S
4 3
1.7
2 1.6
1 0
1.5 0
10
20 30 Time t [s]
40
50
0
10
20 30 Time t [s]
40
50
0
10
20 30 Time t [s]
40
50
0.8 0.7 0.3 L [m]
C
0.6 0.5
0.25
0.4 0.2 0.3 0.2
0.15 0
10
20 30 Time t [s]
40
50
Fig. 9. Model parameters during simulation.
ARTICLE IN PRESS Input torque to wheel axis ur[Nm]
J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
451
60 50 40 30 20 10 0 −10 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 10. Input torque to the wheel axis ur ðN mÞ.
Friction force Ff [N]
250 True values Estimated values
200 150 100 50 0 -50 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 11. Friction force estimation.
Friction force Ff [N] 165
175
160
170
155
165
150
160
145
155
140
150
135 9.95
10
10.05
145 14.95
210
200
200
180
190
15
15.05
160
180 140
170
120
160 150 24.9
24.95
25
25.05
100 25.1 35.95 Time t [s]
36
36.05
36.1
Fig. 12. Friction force estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
452
where s0 ; s1 s0 d 0 ¼ p22 þ p2 ; s1
c0 ¼ p21 þ p1
(37) becomes p1 rF N ; c2 ¼ ðp1 s1 þ s0 Þ, J J p rF N ðp2 s1 þ s0 Þ. d1 ¼ 2 ; d2 ¼ J J
c1 ¼
GðsÞ ¼
1 1 1 1 ¼ . k2 1 sJ=rF N s1 k2 k2 1 þ sT e
(59)
Obviously, estimation error dynamics is now of the first order and depends only on the gain k2 , which can be set to get desired value of time constant T e . For example, for k2 ¼ 0:4 and parameter values from Table 1, time constant T e gets the value of T e 5 ms. Therefore, with K ¼ ½0 0:4T (see Fig. 8) one can expect response time of y~ to the friction model errors of about 15 ms (three time
ð58Þ
Apart from poles p1 and p2 , dynamics of the transfer function (37) is also influenced by its zero z1 ¼ s0 =s1 , whose location does not depend on the selected gain vector K. However, if the gain k1 is set to value k1 ¼ 0 perfect zero/pole cancellation takes place and transfer function
x 10−4 Bristle deflection z[m]
12 True values Estimated values
10 8 6 4 2 0 -2 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 13. Bristle deflection estimation.
Bristle deflection z [m] x 10−4
x 10−4
8.05
8.65 8.6
8 8.55 7.95 8.5 7.9 9.95
10
10.05
8.45 14.95
x 10−3
15
15.05
x 10−4
1.03
5.8
1.025 5.7
1.02 1.015
5.6 1.01 1.005
5.5
1 24.9
24.95
25
25.05
5.4 25.1 35.95 Time t [s]
36
36.05
36.1
Fig. 14. Bristle deflection estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
Wheel angular velocity ω [rad/s]
constants), and response time of the estimation error compensation several time slower. 1 Choices of matrix P1 in (28) W in (27) and coefficient g are less critical than choice of the gain vector K as they cannot cause violence of the passivity condition. However, 1 P1 determine the convergence rate of the network W and g parameters W2 and parameter b, and consequently convergence rate of the estimation error compensation. If they are greater the convergence is faster, but the estimator is more sensitive to the measurement noise. By simulations we found out that it is appropriate to set g1 ¼ 50 and P1 W ¼ 50 I, where I is identity matrix of size ½N h þ 1 N h þ 1.
453
5.3. Results of simulation experiments Results of a simulation experiment showing robust behavior of the proposed estimator are given in Figs. 9–18. The main purpose of the simulation experiment was validation of the estimator robustness to changes of the friction model parameters. Four parameters of the friction model were changed (see Fig. 9): (i) the road condition parameter y changed from value 1 that corresponds to dry asphalt to value 5 that corresponds to the snow road condition at t ¼ 10 s; (ii) model parameter mS changed its value from 1.9 to 1.55 at t ¼ 15 s; (iii) model parameter mC changed its value from 0.68 to 0.3 at t ¼ 25 s; and (iv) the
70 True values Estimated values
60 50 40 30 20 10 0 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 15. Wheel rotational speed estimation.
Wheel angular velocity ω[rad/s] 20.85
11.2
20.8
11.15
20.75 11.1
20.7
11.05
20.65 20.6
11
20.55
10.95 9.95
10
10.05
20.5 14.95
15
15.05
60.8
43.8
60.78 43.6 60.76 43.4
60.74 60.72
43.2 60.7 43 24.9
24.95
25
25.05
60.68 25.1 35.95 Time t [s]
36
36.05
36.1
Fig. 16. Wheel rotational speed estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
454
Δ(f(vr) + fω(ω))[1/s]
250 True values Estimated values
200 150 100 50 0 −50 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 17. Estimation of the model uncertainty Dðf ðvr Þ þ f o ðoÞÞ:
Model uncertainty Δ(f(vr) + fω(ω))[1/s] 40
12 10
35 8 30
6 4
25 2 0 9.95
10
10.05
20 14.95
200
50
190
40
180
30
170
20
160
10
150 24.9
24.95
25
25.05
0 25.1 35.95 Time t [s]
15
36
15.05
36.05
36.1
Fig. 18. Estimation of the model uncertainty Dðf ðvr Þ þ f o ðoÞÞ—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
contact patch length L changed its value from 0.2 to 0.3 (m) at t ¼ 36 s. The one-wheel friction model and estimator were excited with the input torque ur of the waveform shown in Fig. 10. Responses of actual and estimated values of the most important process variables to the imposed input torque are shown in Figs. 11–18. Fig. 11 shows waveforms of actual and estimated friction force during the whole experiment and Fig. 12 detailed views around time instances of the friction model parameter changes. Figs. 13 and 14 show waveforms of the bristle deflection and Figs. 15 and 16 waveforms of the wheel angular velocity. Finally, Figs. 17 and 18 show estimation of the
friction model uncertainty Dðf ðvr Þ þ f o ðoÞÞ. It is obvious that the proposed estimator provides fast and robust tracking of the actual friction force either it changes due to road condition changes or due to model parameter changes. The estimator response time is below 50 ms in all cases, which is very fast considering wide ranges of imposed changes. In order to check noise influence on the estimator behavior, random noise of variance 0:5 (higher than expected in real situations) was added to the wheel angular speed signal and the same experiment was repeated without changing the estimator parameters. Comparing obtained results for friction force estimation, shown in Figs. 19
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
455
Friction force Ff [N]
250 True values Estimated values
200 150 100 50 0 −50 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 19. Friction force estimation.
Friction force Ff [N] 170
180 175
160
170 165
150 160 155
140
150 130 9.8
9.9
10
10.1
10.2
145 14.8
210
170
200
160
14.9
15
15.1
15.2
35.9
36
36.1
36.2
150
190
140 180 130 170
120
160 150 24.8
110 24.9
25
25.1
100 25.2 35.8 Time t [s]
Fig. 20. Friction force estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
and 20, with noiseless results, shown in Figs. 11 and 12, it can be concluded that proposed estimation scheme successfully copes with noise. Fast and robust friction force estimation in highly changing road conditions, and particularly robustness to the friction model uncertainties, make possible usage of the proposed estimation scheme very attractive. However, it should be noted that currently in most middle class cars only road conditions are qualitatively estimated, for which low-resolution sensors and low-performance microcontrollers are used. In order to implement our estimator more advanced sensors and microprocessor systems are needed. Usage of such systems can be expected in the near future,
since the trend in automotive industry is increasing safety and autonomy of the cars, which cannot be achieved without them. We expect that our estimation scheme will find its way to real applications in that framework. 6. Conclusion A new scheme for robust tire/road friction force estimation is proposed. An RBF neural network is used to compensate negative influence of the friction model uncertainties to the accuracy of the force estimation. By splitting the nonlinearity of the friction model in two parts it is achieved that RBF network input comprises only one
ARTICLE IN PRESS 456
J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
variable, which enables that the centers and the spreads of the hidden layer neurons can be set off-line in a very simple manner. The output layer weights are adapted on-line by an adaptation law derived applying Lyapunov stability analysis. Such an approach ensures fast and robust estimation of the tire/road friction force even in cases when wide and abrupt changes of the friction parameters occur, which is verified by simulation experiments on the one-wheel friction model. The authors are hopeful to do experimental testing in their future work. Acknowledgments Support from the Ford Motor Company and the Ministry of Science, Education and Sports of the Republic of Croatia is gratefully acknowledged. References Bakker, E., Pacejka, H.B., Lidner, L., 1987. A new tire model with an application in vehicle dynamics studies. SAE Paper. Baric´, M., Petrovic´, I., Peric´, N., 2005. Neural network based sliding mode control of electronic throttle. Engineering Applications of Artificial Intelligence 18 (8), 951–961. Bozˇic´, A., Petrovic´, I., Peric´, N., Matusˇ ko, J., 2002. Experimental investigations of the rubber-asphalt sliding pair. In: ASME International Congress and Exp.—Symposium on Advanced Automotive Technologies, New Orleans, USA. Burckhardt, M., 1993. Fahrwerktechnik: Radschlupf-regelsysteme. Vogel, Wuertzburg, Germany. Canudas de Wit, C., Horowitz, R., 1999. Observers for tire/road contact friction using only wheel angular velocity information. In: 38th IEEE Conference on Decision and Control, Phoenix, USA. Canudas de Wit, C., Tsiotras, P., 1999. Dynamic tire friction models for vehicle traction. In: 38th IEEE Conference on Decision and Control, Phoenix, USA. Canudas de Wit, C., Olsson, H., Astro¨m, K.J., Lischinsky, P., 1995. A new model for control of systems with friction. IEEE Transactions on Automatic Control 40 (3), 419–425. Chen, T., Chen, H., 1995. Approximation capability to functions of several variables, nonlinear functionals, and operators by radial basis function neural networks. IEEE Transactions on Neural Networks 6 (4), 904–910.
Chen, T., Chen, H., Liu, H., 1995. Approximation capability in c(r) n by multilayer feedforward networks and related problems. IEEE Transactions on Neural Networks 6, 25–30. Deur, J., 2001. Modeling and analysis of longitudinal tire dynamics based on the lugre friction model. In: 3rd IFAC Workshop Advances in Automotive Control, Karlsruhe, Germany. Deur, J., Asgari, J., Hrovat, D., 2004. A 3d brush-type dynamic tire friction model. Vehicle System Dynamics 42, 133–173. Dieckmann, T., 1992. Assessment of road grip by way of measured wheel. In: XXIV FISITA Congress, London. Gustafsson, F., 1997. Slip-based tire-road friction estimation. Automatica 3 (6), 1087–1099. Horn, R.A., Johnson, C.R., 1999. Matrix Analysis. Cambridge University Press, Cambridge. Hozwarth, F., Eichhorn, U., 1993. Non-contact sensors for road conditions. Sensors and Actuators (37–38), 121–127. Kiencke, U., Daiss, A., 1994. Estimation of tyre friction for enhanced abssystems. In: International Symposium on Advanced Vehicle Control, AVEC 1994, Tokyo, Japan. Kiencke, U., Nielsen, L., 2000. Automotive Control Systems. Springer, Berlin. Matusˇ ko, J., Petrovic´, I., Peric´, N., 2002. Application of extended Kalman filter for road condition estimation. In: 10th International Power Electronics and Motion Control Conference, EPE-PEMC 2002, Cavtat and Dubrovnik, Croatia. McCulloch, W.W., Pitts, W., 1943. A logical calculus of the ideas imminent in nervous activity. Bulletin of Mathematical Biophysics 5, 115–133. Moody, J., Darken, C., 1989. Fast-learning in networks of locally-tuned processing units. Neural Computations 1, 281–294. Pal, C., Higiwara, I., Morishita, S., Inoue, H., 1994. Application of neural networks in real time identification of dynamic structural response and prediction of road friction coefficient m from steady state automobile response. In: International Symposium on Advanced Vehicle Control, AVEC 1994, Hiroshima, Japan. Park, J., Sandberg, I.W., 1991. Universal approximation using radialbasis-function networks. Neural Computation 3 (2), 246–257. Slotine, J.J., Li, W., 1990. Applied Nonlinear Control. Prentice-Hall, Englewood Cliffs, NJ. Umeno, T., Ono, E., Asano, K., Ito, S., Yamamoto, M., Yasui, Y., Sawada, M., 2002. Estimation of tire–road friction using wheel speed vibration. In: International Symposium on Advanced Vehicle Control, AVEC 2002, Hiroshima, Japan. van Zanten, A., Ruf, W.D., Lutz, A., 1989. Measurement and simulation of transient tire forces. SAE Technical Paper Series, Paper 890640. van Zanten, A.T., Erhardt, R., Leindesfeind, K., Pfaff, G., 1998. Vdc systems development and perspective, In: SAE International Congress and Exposition, Detroit, USA.
Engineering Applications of Artificial Intelligence 21 (2008) 442–456 www.elsevier.com/locate/engappai
Neural network based tire/road friction force estimation Jadranko Matusˇ ko, Ivan Petrovic´, Nedjeljko Peric´ Faculty of Electrical Engineering and Computing, University of Zagreb, Croatia Received 14 May 2006; received in revised form 1 March 2007; accepted 9 May 2007 Available online 2 July 2007
Abstract This paper deals with the problem of robust tire/road friction force estimation. Availability of actual value of the friction force generated in contact between the tire and the road has significant importance for active safety systems in modern cars, e.g. anti-lock brake systems, traction control systems, vehicle dynamic systems, etc. Since state estimators are usually based on the process model, they are sensitive to model inaccuracy. In this paper we propose a new neural network based estimation scheme, which makes friction force estimation insensitive to modelling inaccuracies. The neural network is added to the estimator in order to compensate effects of the friction model uncertainties to the estimation quality. An adaptation law for the neural network parameters is derived using Lyapunov stability analysis. The proposed state estimator provides accurate estimation of the tire/road friction force when friction characteristic is only approximately known or even completely unknown. Quality of the estimation is examined through simulation using one wheel friction model. Simulation results suggest very fast friction force estimation and compensation of the changes of the model parameters even when they vary in wide range. r 2007 Elsevier Ltd. All rights reserved. Keywords: Tire/road friction; Neural network; Estimation; Lyapunov stability; Antilock braking systems; Traction control systems; Lugre friction model; Model uncertainty estimation; Radial basis function networks
1. Introduction Active safety and comfortability of modern cars has become a major concern of many research activities. Above all, this is caused by the market demands for more sophisticated equipment in nowadays cars with no significant effect to the overall car price. As a consequence, these enhancements are primarily directed to better usability of the existing hardware. Almost every new car is equipped with sophisticated electronic systems such as the anti-lock brake system (ABS), traction control system (TCS) (also called ASR—anti slip regulation), electronic stability program (ESP), etc. ABS works by sensing slippage at the wheels during braking and adjusting braking pressure to ensure maximum contact force between the tires and the road with the purpose of not only shortening the braking distance but also of ensuring the driveability in critical situations (Kiencke and Nielsen, Corresponding author.
E-mail addresses: [email protected] (J. Matusˇ ko), [email protected] (I. Petrovic´), [email protected] (N. Peric´). 0952-1976/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.engappai.2007.05.001
2000). Similar to ABS, TCS also tries to maximize friction force between the tires and the road but during acceleration. This is done by controlling the engine power and the brakes. Both ABS and TCS are primarily concerned with control of the longitudinal (front-rear) motion of the vehicle. On the other hand an ESP system deals with the problem of the uncontrolled lateral vehicle motion (van Zanten et al., 1998). It uses brake actuators and engine torque to control yaw rate of the vehicle but also to limit the side slip angle, ensuring vehicle to follow the path intended by the driver as best as possible. All of the above-mentioned control systems have the purpose of efficient transmission of the forces from vehicle wheels to the road surface. Amount of the forces transferred in this contact is determined by the tire/road friction characteristics. Typical friction characteristics for different road conditions are shown in Fig. 1. Obviously, the goal is to maintain the operating point near the peak values of the tire/road friction characteristics. However, finding the optimal operating point is not a trivial task if the value of the tire/road friction force is not available. Therefore, friction estimation task plays important role in
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
443
presented in (Burckhardt, 1993): mðsÞ ¼ c1 ð1 ec2 s Þ c3 s.
Braking/traction force
Dry asphalt
The main drawback of this model is its nonlinearity in parameters ci , which makes the parameter estimation process difficult. To overcome this problem a rational function is used in Kiencke and Daiss (1994) to describe m2s relationship:
Wet asphalt
mðsÞ ¼ Snow
Slip Fig. 1. Typical tire/road static curves for different road conditions.
the design of the automotive control systems. Beside estimation accuracy, it is important to ensure high estimation availability. This means that estimator should provide useful information not only for various road conditions, but also in situations when friction modelling uncertainties are present. In literature one can find many various approaches to the tire/road friction modelling and estimation. The short overviews are given in Sections 2 and 3, respectively. Neural network based friction force estimator, which we propose as a valuable alternative to other approaches, is presented in Section 4. The key idea is to use a neural network, as the universal approximator, to compensate friction model inaccuracies. It has been added to the state estimator, which is based on nominal friction model, to provide robust estimation. Neural network parameters are adapted online with the training algorithm based on Lyapunov technique. Results of simulation experiments with our estimator is given in Section 5. The paper ends with the conclusion remarks in Section 6. 2. Overview of the tire/road friction modelling approaches Building an accurate tire/road friction model is important for the following reasons: (i) simulation of the overall vehicle behavior, (ii) design of the automotive control systems and (iii) estimation of the friction force. Historically, tire/road friction effects were described by the static models, either physical or empirical. The static friction models usually describe relationship between friction force F f (or more often actual friction coefficient m ¼ F f =F n , where F n is normal force in the tire/road contact) and dimensionless quantity s called slip using algebraic equations. The quantity slip s is usually defined as ro v , (1) s¼ maxðro; vÞ where o is wheel angular speed, v is vehicle speed and r is wheel radius. Simple static tire/road friction model is
(2)
ks . as2 þ bs þ 1
(3)
This model is primarily used for estimation purposes since it is linear in parameters ðk; a; bÞ and thus classical recursive least square method can be used as estimation algorithm. However, the most widely used static model is Pacejka model, also called ‘‘Magic formula’’ (Bakker et al., 1987). The model uses trigonometric functions to describe relationship between slip and friction force generated in the contact between the tire and the road: F f ðsÞ ¼ D sinðC arc tanðB s EðB s arc tanðB sÞÞÞÞ, (4) where B; C; D; E are model parameters. It is proven that Pacejka model has ability to approximate experimental data (Bakker et al., 1987). It can be used for both simulation and estimation tasks. However, as all other static models it has significant disadvantage related to inability to describe low slip effects. To resolve this problem a deeper insight is needed into physical effects taking place in the tire/road contact surface. The best insight in what is happening in the tire/ road contact surface is given by so called ‘‘brush’’ models (van Zanten et al., 1989). These models assume that contact between the tire and the road is realized through many tiny and massless elastic elements, called bristles. Bristle deformations are direct cause of the generated friction force. When a bristle enters into contact with the road, it starts to deform with deformation rate proportional to the relative velocity (adhesion region in Fig. 2). After reaching deformation limit, there is no additional deformation and bristle starts to slide over the contact patch (sliding region in Fig. 2). The original brush model describes the deformation of a bristle as follows: ( vr ; zðzÞozm ; z_ðzÞ ¼ (5) 0; zðzÞXzm ; where z is the bristle deformation, vr ¼ ro v is relative velocity and z is a spatial variable related to the bristle position while passing through the contact patch. Friction force generated by a bristle is given by F ðzÞ ¼ s0 zðzÞ þ s1 z_ðzÞ,
(6)
where s0 and s1 are the rubber longitudinal lumped stiffness and damping, respectively. Overall friction force generated in the tire/road contact is obtained by integration of the bristle forces over the whole contact patch. While saturation function (5) is used to describe deformation of a bristle in the original brush model, the
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
444
=
Ff Fn
confirmed that overshoot exists in the friction force at transition from the adhesion regime to the sliding regime of the bristles in the contact patch, but LuGre model cannot capture it. Also, LuGre model predicts the same initial slope of the static characteristics (at s ¼ 0) for all road conditions (all values of coefficient y), which is in conflict with experimental results in Dieckmann (1992). In spite of these possible shortcomings LuGre model is still the most widely used dynamic tire/road friction model and therefore is used in our work as the basis for development of the tire/ road friction force estimators. We expect that it will be straightforward to adapt our estimator to possible future changes of the LuGre model.
Point of maximum friction
Sliding region
Adhesion region
Slip s Fig. 2. Regions on the tire/road static friction curve.
3. Overview of the tire/road friction estimation approaches
first-order nonlinear function is used in LuGre model (Canudas de Wit et al., 1995; Canudas de Wit and Tsiotras, 1999): s0 jvr j z_ðzÞ ¼ vr zðzÞ, (7) gðvr Þ where gðvr Þ is Stribeck function defined as pffiffiffiffiffiffiffi gðvr Þ ¼ mC þ ðmS mC Þ e vr =vs .
(8)
In (8) mC denotes normalized static friction, mS normalized maximum friction and vs Stribeck velocity. Although LuGre model (6)–(8) has compact form, it is primarily used for simulation purposes, since spatial integration over the contact patch is computationally very demanding. For estimation as well as control system design purposes a less complex model needs to be used. Therefore, instead of describing dynamics of each particular bristle deformation, an average bristle deflection is introduced in (Canudas de Wit and Tsiotras, 1999): Z 1 L z¯ ¼ zðzÞ dz, (9) L 0 where L is the length of the contact patch. Usage of the average value of the bristles deformation led to the lumped LuGre model given by s jv j _z¯ ¼ vr y 0 r z¯ , (10) gðvr Þ where the coefficient of road adhesion y is introduced in order to model friction characteristic changes caused by the road condition changes. Lumped LuGre model (10) has a steady-state error, which is eliminated in Deur (2001) and Deur et al. (2004) by introducing an additional term s0 jvr j k (11) þ rjoj z¯ , z_¯ ¼ vr y gðvr Þ L where k is a constant. Some researchers criticized LuGre model as not being able to capture all dynamic effects in the contact patch. For example, in Bozˇic´ et al. (2002) it is experimentally
Estimation of the friction force generated in the contact patch between the tire and the road is a quite demanding task. Special hardware sensors described in Hozwarth and Eichhorn (1993) are considered as inadequate due to their high prices. More appropriate solution is to use readily accessible, easily measurable, signals to estimate the tire/ road friction in an indirect way. In literature many different approaches to the problem of tire/road friction estimation have been proposed. Generally, all approaches are directed either to the estimation of actual friction force (or friction coefficient) or just to the qualitative estimation of road conditions. While later approaches are simpler to develop and implement, they do not enable applications of the advanced ABS and TCS that search for optimal values of the slip s (point of maximum friction, see Fig. 2). Instead, the designer has to predefine referent slip for all road conditions based on the experience gained through experimental tests. On the other hand, if the estimator estimates actual friction force, and if the slip value is known, it is possible to implement advanced ABS and TCS. A road condition estimator based on Kalman filter is proposed in Gustafsson (1997). The approach supposes that the slope of the friction force characteristic at low slip depends on the road conditions (see Fig. 3), which is experimentally confirmed in Dieckmann (1992). The appropriate model which describes the tire/road friction in low slip region is given by the following equation: 1=k sðtÞ ¼ ðmðtÞ1Þ þ eðtÞ, (12) d where k is low-slip slope of the friction characteristic, d is the slip offset and eðtÞ is the white noise signal (see Fig. 3). In order to ensure good tracking performance this estimator has an additional subsystem for detection of abrupt changes of the road conditions. In Canudas de Wit and Horowitz (1999), a state estimator based on passivity theory is presented. It uses standard dynamic LuGre model for description of the tire/ road friction. This estimator is able to estimate all state
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
sph a
lt
Ff FN
Dry a
=
ow
Sn
1
2
Slip s
Fig. 3. Slopes of the m2s curves near zero slip for different road conditions.
variables as well as the coefficient of the road adhesion y, which represents road conditions (dry asphalt, wet asphalt, snow, ice, etc.). Convergence of the estimation error to zero is guaranteed by the adaptation law derived applying Lyapunov stability theory. However, proposed estimator is sensitive to the model inaccuracies and measurement noise. Another approach to the simultaneous estimation of the state variables and the coefficient of the road adhesion is presented in Matusˇ ko et al. (2002). The estimator is based on modified LuGre dynamic friction model and uses extended Kalman filter, additionally augmented by an integral term (PI Kalman filter) in order to capture variations of the road adhesion coefficient. Simulation results suggest that the estimator is appropriate when significant level of the measurement noise is present. However the approach relies on exact tire/road friction model and thus uncertainties in the friction model parameters may significantly degrade estimation accuracy. In Pal et al. (1994) a neural network has been applied for the road condition estimation. The network parameters were tuned by an off-line training algorithm applied on the field test data acquired for different road conditions. The approach suffers from two big drawbacks: (1) the output of the estimator can get only several discrete values meaning dry, wet, snow or icy road, but cannot estimate friction force or state variables, and (2) long series of field test data are required to train the neural network. A somewhat different approach is presented in Umeno et al. (2002), where spectral analysis of the wheel speed signal is used to determine the road conditions. The approach is based on the tire vibration model, instead of the friction model, and relies on the assumption that the road excites the tire/road vibrational system with the strength dependant on the road conditions. As a result of this excitation higher frequency harmonic (around 40 Hz) occurs in signals measured by wheel speed sensors. The main
445
advantage of this estimator is its ability to estimate the road condition even at zero slip (free rolling motion). The necessity of high-resolution wheel speed sensors and problems due to weak excitations on icy roads are significant drawbacks of the method. In the introductory section, it was pointed out that the estimator should posses high availability. The majority of the estimators are based on some type of the tire/road friction model, where usual assumption is that model accurately describes effects in the tire/road contact patch. However, this is only partially fulfilled, since models used for estimation purposes are rather simple, and thus modelling errors are inherent to them. Besides, parameters of the processes in the contact patch may vary with time and initially accurate model could drift away from the true one. Therefore, it is important to ensure robustness of the estimation to the model uncertainties. 4. Neural network estimator design As aforementioned, vehicle behavior on the road is dominantly determined by the amount of forces transferred between vehicle wheels and road. Generally, tire/road friction forces at the vehicle wheels may mutually differ and therefore it is important to estimate them separately for each wheel. This is the reason why one-wheel friction model (Fig. 4) is used in this paper as the basis for design of the tire/road friction estimator. One-wheel friction model is described by the following set of equations: m_v ¼ F f þ F n s2 vr ,
(13)
_ ¼ rF f þ ur so o, Jo
(14)
where m is the wheel mass, J is the moment of inertia of the wheel, ur is the torque acting on wheel axis, s2 is the viscous relative damping and so is the coefficient of viscous friction. Dynamics of the friction force F f in Eqs. (13) and (14) is modelled using lumped LuGre dynamic friction model described in Section 2. For the sake of clarity the model equations are again written here, but in somewhat different form (further on variable z is used instead of variable z¯ to denote average bristle deflection): F f ¼ F n ðs0 z þ s1 z_Þ,
(15)
s0 jvr j k þ rjoj z ¼ vr f ðvr Þ z f o ðoÞ z, z_ ¼ vr y gðvr Þ L (16) ω
v r
Fn
Ff
Fig. 4. One wheel friction model.
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
446
gðvr Þ ¼ mC þ ðmS mC Þ e
pffiffiffiffiffiffiffiffiffi
jvr =vs j
,
(17)
where f ðvr Þ ¼ yðs0 jvr j=gðvr ÞÞ and f o ðoÞ ¼ ðk=LÞrjoj. Although LuGre dynamic friction model can satisfactorily describe dynamics of the friction at the tire/road contact patch, many field test data are usually required to parameterize it accurately. In this paper, it is assumed that functions f ðvr Þ and f o ðoÞ in Eq. (16) are only partially known and the aim is to design friction force estimator which provides accurate estimation of the friction force under these assumptions. The nonlinear function f ðvr Þ can be written as the sum of a known function f n ðvr Þ (nominal value of function f ðvr Þ), and an unknown function Df ðvr Þ: s0 jvr j s0 jvr j f ðvr Þ ¼ y þD y ¼ f n ðvr Þ þ Df ðvr Þ. (18) gðvr Þ n gðvr Þ In order to compensate the unknown nonlinear function Df ðvr Þ in Eq. (18) a neural network is employed. Neural networks commonly consist of big number of highly interconnected neurons organized in two or more layers (McCulloch and Pitts, 1943). These interconnections give them universal approximation capability, i.e. a neural network can approximate arbitrary continuous nonlinear function with desired accuracy (Chen et al., 1995; Chen and Chen, 1995; Park and Sandberg, 1991) if consisted of at least two layers. Generally, a two-layer neural network with linear output layer can be described by the following equation: yNN ¼ W2 wðW1 ; uÞ, n
can be represented with a scaling factor b. If its optimal value bn is found, then f o ðoÞ can be written as k (21) f o ðoÞ ¼ bn rjoj ¼ bn f on ðoÞ, L where f on ðoÞ is the function f o ðoÞ with the nominal value of parameters k and L. The task of the proposed estimator is to estimate optimal values of the network parameters Wn2 and parameter bn . As the first step in the estimator design procedure the following state variable transformation was performed (Canudas de Wit and Horowitz, 1999): x1 ¼ rF n s1 z þ Jo, x2 ¼ z.
One-wheel friction model (13)–(17) in the space of new variables can be written as x_ ¼ Ax þ Bf ðy; uÞDx þ Bf o ðyÞDx þ Ru þ Ey, y ¼ CT x, " x¼
x1 x2
2 s0 A ¼ 4 s1
# ;
" D ¼ ½ 0 1 ;
R¼
3 1 6 J 7 7 C¼6 4 rF n s1 5; J
3
1
0
0
1
u¼
v ur
"
5;
0
"
# ;
0
#
; 1 3 2 s0 J so 5 ; E ¼ 4 s1 r
B¼
# ;
y ¼ ½ o :
For the system model given by Eq. (23) the state estimator in standard predictor–corrector form (Fig. 5) is given by +
y(t)
E
u(t)
−
K
A
1 s
R
B
ˆ x(t)
C
D
X
(20)
Uncertainty in function f o ðoÞ is the consequence of variations of the parameters k and L. Parameter k is lumped coefficient of the LuGre model and varies between 1.1 and 1.4 depending on the tire type (Deur et al., 2004). On the other hand the length of the contact patch L varies with variations of the pressure in the tire. These variations
0
0
2
m
f ðvr Þ ¼ f n ðvr Þ þ Df ðvr Þ f n ðvr Þ þ Wn2 wðvr Þ.
ð23Þ
where
(19)
where u 2 R is network input, y 2 R is network output, w: Rn 7!Rm is nonlinear kernel function and W1 and W2 are parameter matrices of hidden and output network layer, respectively. Various neural network structures are obtained by choosing different kernel functions. As mentioned above, friction characteristic may vary with time. This implies that neural network parameters (weights) should be updated in an on-line manner. Hereafter we suppose the usage of a two-layer neural network with hidden layer weights W1 set off-line, and linear output layer weights W2 updated on-line. The reasons for choosing such a network architecture are computational simplicity of the on-line training algorithm and satisfactory approximation of function Df ðvr Þ. If parameters W2 are adjusted in such a way that their optimal values Wn2 are found, then neural network of proper structure gives good approximation of the function Df ðvr Þ, and Eq. (18) can be written as
ð22Þ
fn (y,u) ˆ
fωn (y) ˆ ) Δf (y,u,W 2
Fig. 5. State estimator principle scheme.
ˆ y(t)
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
^ ^ bn Bf on ðyÞDx bBf on ðyÞDx ~ ^ ¼ bn Bf ðyÞDx~ þ bBf ðyÞDx,
the following equations: ^ ^ 2 wðy; uÞDx^ ^ þ BW x_^ ¼ Ax^ þ Bf n ðy; uÞDx^ þ Bbf on ðyÞDx ^ þ Ru þ Ey þ Kðy yÞ, T
^ y^ ¼ C x.
on
ð24Þ
The designer has the ability to influence performances of the estimator (24) by tuning the value of the gain vector ^ 2 of K ¼ ½k1 k2 T and adapting the output layer weights W ^ In order to guarantee the neural network and parameter b. convergence of the estimation, the adaptation algorithm of ^ 2 and parameter b^ as well as rules the network weights W for the vector K selection are derived applying Lyapunov stability analysis. The results of the analysis are summarized in Assumption 1 and Theorem 1. Assumption 1. Vehicle velocity v and wheel angular velocity o are bounded, i.e.
~ 2 wðy; uÞDx^ x_~ ¼ ðA KCT Þx~ þ B½Wn2 wðy; uÞDx~ þ W n ~ ^ þ f ðy; uÞDx~ þ b f ðyÞDx~ þ bf ðyÞDx, n
on
on
~ y~ ¼ CT x,
ð31Þ
^ ^ 2 and b~ ¼ b b. ~ 2 ¼W W where W 2 Let Lyapunov function candidate be given as n
n
~ T PW W ~ 2 g þ gb~ 2 , V ¼ x~ T Px~ þ trfW 2
(32)
where P is a symmetric positive definite matrix, and trfg is matrix trace operator defined as (Horn and Johnson, 1999) n X
xi;i .
(33)
i¼1
ð25Þ
It is important to note that Assumption 1 is introduced only due to theoretical reasons, since in reality quantities v and o are always bounded. Theorem 1. Let the system and state estimator be given by Eqs. (23) and (24), respectively, and let the following conditions for the gain vector K be satisfied: k1 4rF N s1 k2 , k2 o0.
According to the Lyapunov theorem of stability (Slotine and Li, 1990), the time derivative of the Lyapunov function should be negative semidefinite in order to ensure system stability: _~ ~_ 2 g þ 2gb~ bo0. ~ T PW W V_ ¼ x_~ T Px~ þ x~ T Px_~ þ 2 trfW 2 V_ ¼ x~ T ððA KCT ÞT P þ PðA KCT ÞÞx~ þ 2x~ T PBWn2 wðy; uÞDx~ ~ 2 wðy; uÞDx^ þ 2x~ T PBW
ð26Þ
_^ 1 ~ x^ T DT wT ðy; uÞ, W 2 ¼ PW y
(27)
_ ^ b^ ¼ g1 y~ f on ðyÞDx,
(28)
þ 2x~ T PBf n ðy; uÞDx~ þ 2bn x~ T PBf on ðyÞDx~ _^ ~ _^ ~ T PW W þ 2b~ x~ T PBf on ðyÞDx^ 2trfW 2 g 2gbb. 2
Gain vector K ¼ ½k1 k2 is chosen in such a way that estimation error dynamics (31) is passive (Slotine and Li, 1990), i.e. transfer function GðsÞ ¼ CT ðsI A þ KCT Þ1 B
Proof of the Theorem 1. Estimation error dynamics is obtained by subtracting Eqs. (24) from Eqs. (23):
GðsÞ ¼
on
T
~ y~ ¼ C x,
ð29Þ
^ Taking into account where x~ ¼ x x^ and y~ ¼ y y. Eqs. (18), (20) and (21), and following equalities: ^ 2 wðy; uÞDx^ ¼ BW wðy; uÞDx~ BW2 wðy; uÞDx BW 2 ~ 2 wðy; uÞDx, ^ þ BW n
ð35Þ
T
where PW is a symmetric positive definite matrix and g is a positive scalar coefficient. With the adaptation laws (27) and (28), and condition (26), the convergence of the state estimates toward their true values is ensured.
x_~ ¼ ðA KCT Þx~ þ Bf ðy; uÞDx Bf n ðy; uÞDx^ þ Bf o ðyÞDx ^ ^ 2 wðy; uÞDx, ^ Bbf ðyÞDx^ BW
(34)
By inserting Eq. (31) into (34) it is obtained:
Let the adaptation laws of the neural network output-layer ^ 2 and parameter b^ be given by the following parameters W equations:
n
ð30Þ
on
equations describing estimation error dynamics can be written as
trfXg ¼
jvjovmax , jojoomax .
447
(36)
is strictly positive real function. By inserting (23) into (36) the transfer function of the estimation error dynamics becomes rF N s1 =J s þ rF N s0 =J þ sðs0 =s1 þ k1 =J rF N s1 k2 =JÞ rF N s0 k2 =J b1 s þ b0 ¼ 2 . ð37Þ s þ a1 s þ a0 s2
Sufficient conditions for strict positiveness of the transfer function (37) are a0 ; a1 40 ) k2 o0 ^ k1 4rF N s1 k2 Js0 =s1 , b0 a0 þ ðb1 a1 b0 Þo2 40 ) k2 o0 ^ k1 4rF N s1 k2 ,
ð38Þ
i.e. k2 o0 ^ k1 4rF N s1 k2 .
(39)
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
448
If the system is passive, then, from Kalman–Yakubovich lemma (Slotine and Li, 1990), it follows: 8Q; 9P: ðA KCT ÞT P þ PT ðA KCT Þ ¼ Q, PB ¼ C,
ð40Þ
where Q is a symmetric positive definite matrix. By taking into account Eq. (40) time derivative of the Lyapunov function (35) becomes ~ 2 wðy; uÞDx^ V_ ¼ x~ T Qx~ þ 2x~ T CWn2 wðy; uÞDx~ þ 2~yT W þ 2x~ T Cf n ðy; uÞDx~ þ 2bn x~ T Cf on ðyÞDx~ þ 2b~ x~ T Cf ðyÞDx^ on
_^ ~ T PW W ^_ 2 g 2gb~ b. 2 trfW 2
ð41Þ
Using the following feature of the matrix trace operator trfg: xT y ¼ trfxyT g;
x; y 2 Rn ,
(42)
Eq. (41) can be written in the form V_ ¼ x~ T Qx~ þ 2x~ T CWn2 wðy; uÞDx~ þ 2x~ T Cf n ðy; uÞDx~ _^ ~ y~ T f ðyÞDx^ gbÞ þ 2bn x~ T Cf on ðyÞDx~ þ 2bð on _^ ~ T ðy~ x^ T DT wT ðy; uÞ PW W þ 2 trfW ð43Þ 2 Þg. 2 If the adaptation laws of the neural network output _^ _^ weights W 2 and parameter b are chosen as _^ 1 ~ x^ T DT wT ðy; uÞ, W 2 ¼ PW y
(44)
_ ^ b^ ¼ g1 y~ T f on ðyÞDx,
(45)
then time derivative of the Lyapunov function becomes V_ ¼ x~ T Qx~ þ 2x~ T CWn2 wðy; uÞDx~ þ 2x~ T Cf n ðy; uÞDx~ ~ þ 2bn x~ T Cf on ðyÞDx.
ð46Þ
The following assessment is valid for Eq. (46): pffiffiffiffiffiffi ~ 2 þ 2 N h kCkF kWn2 kF kDkF kxk ~ 2 V_ p lmin ðQÞ kxk ~ 2 þ 2 maxff n ðy; uÞgkCkF kDkF kxk ~ 2, þ 2jbn j maxff on ðyÞgkCkF kDkF kxk
ð47Þ
where lmin ðQÞ is the least eigenvalue of matrix Q, k k is the Euclidian vector norm, k kF is the Frobenius matrix norm, N h is the number of the neurons in network hidden layer and maxff n ðy; uÞg and maxff o ðyÞg are maximum values of the nonlinear functions f n ðy; uÞ and f o ðyÞ given by maxff n ðy; uÞg ¼ y maxff o ðyÞg ¼
s0 jvr max j , gðvr max Þ
k rjomax j. L
ð48Þ
Kalman–Yakubovich lemma, given by Eq. (40), enables that matrix Q is chosen in such a way that following
inequality is satisfied: pffiffiffiffiffiffiffi lmin ðQÞ42 N h kCkF kWn2 kF kDkF þ 2 maxff n ðy; uÞgkCkF kDkF þ 2jbn j maxff on ðyÞgkCkF kDkF ,
ð49Þ
and that is the condition which ensures negative semidefiniteness of the time derivative of the Lyapunov function V_ p0. Consequently, convergence of the state estimates toward their true values is guaranteed. & Remark 1. It is obvious from Eq. (18) that function f ðy; uÞ actually depends only on relative velocity vr ¼ ro v, and therefore vr is the only needed input to the neural network. In that way the choice of the hidden layer parameters becomes very simple, and the needed number of hidden neurons becomes relatively small, which are both very beneficial for practical application. Remark 2. The estimator design procedure could also be performed applying a neural network to approximate the whole nonlinearity in (16), i.e. without splitting it in two separate nonlinear functions f and f o . In that case functions f and f o become f ¼ ðys0 jvr j=gðvr ÞÞ þ ðk=LÞrjoj and f o ¼ 0, respectively. Derivation of the selection rules for the vector K and the neural network weights adaptation law remains the same, and final results of the Lyapunov analysis summarized in Theorem 1 stay valid too. One needs only to set the parameter b equal to zero in all equations in which it appears and to exclude its adaptation equation (28) from Theorem 1.However, if a neural network is used to approximate the whole nonlinearity in (16), the input vector of the network x1 must be twodimensional, either x1 ¼ ½o; v or x1 ¼ ½o; vr . Of course, in that case an advanced algorithm for off-line setting of the hidden layer parameters, and increased number of hidden neurons, need to be used. Remark 3. Although we prove convergence of the state estimates for neural networks with fixed hidden layer parameters, it can also be proven for cases when hidden network parameters are trained on-line. For example, using Lyapunov function as in Baric´ et al. (2005) instead of (32), the convergence can be proven for multi-layer perceptron networks with on-line training of all its parameters. However, one has to be aware of increased computational complexity of the training algorithm. 5. Simulation tests Proposed friction force estimator given by Eqs. (24), (27) and (28) and the one-wheel friction model given by Eqs. (13)–(17) have been implemented in Matlab/Simulink environment and numerous simulation experiments have been performed in order to validate estimator performances. The model parameters used in simulations are given in Table 1. The architecture of the neural network used in simulations, selection of the gain vector K,
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456 Table 1 Model parameters used in simulations Values
s0 r s1 m s2 J Fn mC mS k vS L
200 (N/m) 0.25 (m) 0.4948 (N s/m) 350 (kg) 0.0018 (N s/m) 0.2344 ðkg m2 Þ 1000 ðkg m2 =s2 Þ 0.7 1.9 1.1 3.5 (m/s) 0.2 (m)
x1
W1 − w11
+
− w12 • • •
+
•
v1
y1
x2
2
v11
y11
x21
•
2
v12
y12
x22
− w1Nh
• • • 2
v1Nh
• • •
y1Nh x2Nh
0.2
−15
−10
−5 0 5 10 Relative velocity vr [m/s]
15
20
Then, spreads of the neurons were chosen using P-nearest ðP ¼ 2Þ neighbor heuristics (Moody and Darken, 1989): vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u P u1 X kwj wi k, (51) s1;i ¼ t P j¼1
W2
+
w22 w2Nh
y2 = yNN
giving values: s1 ¼ ½2:291; 1:732; 1:225; 1; 0:922; 0:742; 0:592; 0:5; . . . , 0:5; 0:592; 0:742; 0:922; 1; 1:225; 1:732; 2:291.
w2(N
h+1)
1 matrix P1 W and coefficient g , and simulation results are presented in succession.
5.1. Neural network used in simulations Neural network used in simulations was a RBF network (Fig. 6) described by the following set of equations:
y1i ¼ eðv1i =s1i Þ ;
0.4
Fig. 7. Set of RBF hidden layer activation functions.
Fig. 6. Two-layer RBF neural network.
2
0.6
−20
1
v1i ¼ kx1 w1i k2 ;
0.8
0
w21
• • •
•
Activation functions
1
Parameters
+
449
i ¼ 1 . . . N h,
i ¼ 1 . . . N h,
T
x2 ¼ ½y1 1 , y2 ¼ W2 x2 , yNN ¼ y2 ,
ð50Þ
where x1 ¼ vr is input to the network, v1 is the input vector to Gaussian activation functions of the hidden layer neurons, N h is the number of hidden layer neurons, W1 and s1 are vectors of centers and spreads of hidden layer neurons, respectively, and x2 and W2 are input and weight vectors of the output layer, respectively. Network hidden layer parameters were set off-line and fixed. N h ¼ 16 neurons were used with centers set to the following values: W1 ¼ ½10; 6; 3:5; 2:5; 1:5; 0:8; 0:4; 0:1; . . . , 0:1; 0:4; 0:8; 1:5; 2:5; 3:5; 6; 10.
With these values of W1 and s1 , the RBF network input space is well covered with activation functions in the range of vr 2 ð10; 10Þ (see Fig. 7), which is satisfactory for very wide range of vehicle velocities (up to 180 km/h). Of course, input range of the network can be easily extended in order to cover wider range if necessary. Hidden layer neurons are set more densely for smaller values of relative velocity because nonlinear effects in approximated function Df ðvr Þ are much stronger for small than for large values of vr , what can be easily checked analyzing expressions (17) and (18). Moreover, ABS and TCS usually keep the tire slip at values jsjo0:1, which means that relative velocities are in the range ð5; 5Þ. Higher slip can rarely occur, e.g. when road conditions abruptly change in large amount. Initial values of the neural network output layer weights W2 can be conveniently set to zero, and parameter b to 1. These parameters are then adapted on-line by adaptation laws (27) and (28). 5.2. Selection of the gain vector K, matrix P1 W and coefficient g1 According to Theorem 1, passivity of the estimation error dynamics is ensured if the gain vector K ¼ ½k1 k2 T is selected respecting conditions (26). Passivity area of the estimation error for model parameters given in Table 1 is shown in Fig. 8. By choosing proper values of gains k1 and k2 inside the passivity area desired dynamics of the estimation error with respect to the friction model uncertainties can be obtained. This dynamics is given by transfer function (37), which defines behavior of the ~ and also output error y, ~ when errors estimation error x,
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
450
in friction model occur. Since compensation of the friction model uncertainties is performed by neural network training according to (27) and (28), which use y~ as their input, it is desirable that dynamics given by (37) is as fast as possible, what can be achieved by proper selection of gains k1 and k2 . However, too fast dynamics can result in oversensitivity to measurement noise. Therefore, a tradeoff is necessary. It is suggested that dynamics of (37) is chosen several times faster than desired dynamics of estimation error compensation.
k1
Area of passivity
2000
Dynamics of the transfer function (37) can be adjusted by proper placement of its poles:
p1;2 ¼
a1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a21 4a0 2
.
(52)
Both poles must satisfy characteristic equation of transfer function (37): p2i þ pi a1 þ a0 ¼ 0;
i ¼ 1; 2,
(53)
wherefrom, by substitution of a0 and a1 with corresponding physical variables, the following equations are obtained: p21 þ p1
s0 k1 rF N k2 ðp1 s1 þ s0 Þ ¼ 0, þ p1 s1 J J
(54)
p22 þ p2
s0 k1 rF N k2 ðp2 s1 þ s0 Þ ¼ 0. þ p2 s1 J J
(55)
1000 {k1,k2} = {0,−0.4} k2
1
F
k1
=r
−1000
k2 σ1 n
Solution of Eqs. (54) and (55) yields the following expressions for k1 and k2 :
−2000
Fig. 8. Passivity area of the estimation in k1 , k2 plane.
k1 ¼
c2 d 0 c0 d 2 , d 2 c1 d 1 c2
(56)
k2 ¼
c1 d 0 c0 d 1 , d 1 c2 c1 d 2
(57)
7
2
6
1.9
5 1.8
S
4 3
1.7
2 1.6
1 0
1.5 0
10
20 30 Time t [s]
40
50
0
10
20 30 Time t [s]
40
50
0
10
20 30 Time t [s]
40
50
0.8 0.7 0.3 L [m]
C
0.6 0.5
0.25
0.4 0.2 0.3 0.2
0.15 0
10
20 30 Time t [s]
40
50
Fig. 9. Model parameters during simulation.
ARTICLE IN PRESS Input torque to wheel axis ur[Nm]
J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
451
60 50 40 30 20 10 0 −10 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 10. Input torque to the wheel axis ur ðN mÞ.
Friction force Ff [N]
250 True values Estimated values
200 150 100 50 0 -50 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 11. Friction force estimation.
Friction force Ff [N] 165
175
160
170
155
165
150
160
145
155
140
150
135 9.95
10
10.05
145 14.95
210
200
200
180
190
15
15.05
160
180 140
170
120
160 150 24.9
24.95
25
25.05
100 25.1 35.95 Time t [s]
36
36.05
36.1
Fig. 12. Friction force estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
452
where s0 ; s1 s0 d 0 ¼ p22 þ p2 ; s1
c0 ¼ p21 þ p1
(37) becomes p1 rF N ; c2 ¼ ðp1 s1 þ s0 Þ, J J p rF N ðp2 s1 þ s0 Þ. d1 ¼ 2 ; d2 ¼ J J
c1 ¼
GðsÞ ¼
1 1 1 1 ¼ . k2 1 sJ=rF N s1 k2 k2 1 þ sT e
(59)
Obviously, estimation error dynamics is now of the first order and depends only on the gain k2 , which can be set to get desired value of time constant T e . For example, for k2 ¼ 0:4 and parameter values from Table 1, time constant T e gets the value of T e 5 ms. Therefore, with K ¼ ½0 0:4T (see Fig. 8) one can expect response time of y~ to the friction model errors of about 15 ms (three time
ð58Þ
Apart from poles p1 and p2 , dynamics of the transfer function (37) is also influenced by its zero z1 ¼ s0 =s1 , whose location does not depend on the selected gain vector K. However, if the gain k1 is set to value k1 ¼ 0 perfect zero/pole cancellation takes place and transfer function
x 10−4 Bristle deflection z[m]
12 True values Estimated values
10 8 6 4 2 0 -2 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 13. Bristle deflection estimation.
Bristle deflection z [m] x 10−4
x 10−4
8.05
8.65 8.6
8 8.55 7.95 8.5 7.9 9.95
10
10.05
8.45 14.95
x 10−3
15
15.05
x 10−4
1.03
5.8
1.025 5.7
1.02 1.015
5.6 1.01 1.005
5.5
1 24.9
24.95
25
25.05
5.4 25.1 35.95 Time t [s]
36
36.05
36.1
Fig. 14. Bristle deflection estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
Wheel angular velocity ω [rad/s]
constants), and response time of the estimation error compensation several time slower. 1 Choices of matrix P1 in (28) W in (27) and coefficient g are less critical than choice of the gain vector K as they cannot cause violence of the passivity condition. However, 1 P1 determine the convergence rate of the network W and g parameters W2 and parameter b, and consequently convergence rate of the estimation error compensation. If they are greater the convergence is faster, but the estimator is more sensitive to the measurement noise. By simulations we found out that it is appropriate to set g1 ¼ 50 and P1 W ¼ 50 I, where I is identity matrix of size ½N h þ 1 N h þ 1.
453
5.3. Results of simulation experiments Results of a simulation experiment showing robust behavior of the proposed estimator are given in Figs. 9–18. The main purpose of the simulation experiment was validation of the estimator robustness to changes of the friction model parameters. Four parameters of the friction model were changed (see Fig. 9): (i) the road condition parameter y changed from value 1 that corresponds to dry asphalt to value 5 that corresponds to the snow road condition at t ¼ 10 s; (ii) model parameter mS changed its value from 1.9 to 1.55 at t ¼ 15 s; (iii) model parameter mC changed its value from 0.68 to 0.3 at t ¼ 25 s; and (iv) the
70 True values Estimated values
60 50 40 30 20 10 0 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 15. Wheel rotational speed estimation.
Wheel angular velocity ω[rad/s] 20.85
11.2
20.8
11.15
20.75 11.1
20.7
11.05
20.65 20.6
11
20.55
10.95 9.95
10
10.05
20.5 14.95
15
15.05
60.8
43.8
60.78 43.6 60.76 43.4
60.74 60.72
43.2 60.7 43 24.9
24.95
25
25.05
60.68 25.1 35.95 Time t [s]
36
36.05
36.1
Fig. 16. Wheel rotational speed estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
454
Δ(f(vr) + fω(ω))[1/s]
250 True values Estimated values
200 150 100 50 0 −50 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 17. Estimation of the model uncertainty Dðf ðvr Þ þ f o ðoÞÞ:
Model uncertainty Δ(f(vr) + fω(ω))[1/s] 40
12 10
35 8 30
6 4
25 2 0 9.95
10
10.05
20 14.95
200
50
190
40
180
30
170
20
160
10
150 24.9
24.95
25
25.05
0 25.1 35.95 Time t [s]
15
36
15.05
36.05
36.1
Fig. 18. Estimation of the model uncertainty Dðf ðvr Þ þ f o ðoÞÞ—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
contact patch length L changed its value from 0.2 to 0.3 (m) at t ¼ 36 s. The one-wheel friction model and estimator were excited with the input torque ur of the waveform shown in Fig. 10. Responses of actual and estimated values of the most important process variables to the imposed input torque are shown in Figs. 11–18. Fig. 11 shows waveforms of actual and estimated friction force during the whole experiment and Fig. 12 detailed views around time instances of the friction model parameter changes. Figs. 13 and 14 show waveforms of the bristle deflection and Figs. 15 and 16 waveforms of the wheel angular velocity. Finally, Figs. 17 and 18 show estimation of the
friction model uncertainty Dðf ðvr Þ þ f o ðoÞÞ. It is obvious that the proposed estimator provides fast and robust tracking of the actual friction force either it changes due to road condition changes or due to model parameter changes. The estimator response time is below 50 ms in all cases, which is very fast considering wide ranges of imposed changes. In order to check noise influence on the estimator behavior, random noise of variance 0:5 (higher than expected in real situations) was added to the wheel angular speed signal and the same experiment was repeated without changing the estimator parameters. Comparing obtained results for friction force estimation, shown in Figs. 19
ARTICLE IN PRESS J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
455
Friction force Ff [N]
250 True values Estimated values
200 150 100 50 0 −50 0
5
10
15
20
25 Time t [s]
30
35
40
45
50
Fig. 19. Friction force estimation.
Friction force Ff [N] 170
180 175
160
170 165
150 160 155
140
150 130 9.8
9.9
10
10.1
10.2
145 14.8
210
170
200
160
14.9
15
15.1
15.2
35.9
36
36.1
36.2
150
190
140 180 130 170
120
160 150 24.8
110 24.9
25
25.1
100 25.2 35.8 Time t [s]
Fig. 20. Friction force estimation—detailed view around time instances of the model parameter changes (dashed—true values, solid—estimated values).
and 20, with noiseless results, shown in Figs. 11 and 12, it can be concluded that proposed estimation scheme successfully copes with noise. Fast and robust friction force estimation in highly changing road conditions, and particularly robustness to the friction model uncertainties, make possible usage of the proposed estimation scheme very attractive. However, it should be noted that currently in most middle class cars only road conditions are qualitatively estimated, for which low-resolution sensors and low-performance microcontrollers are used. In order to implement our estimator more advanced sensors and microprocessor systems are needed. Usage of such systems can be expected in the near future,
since the trend in automotive industry is increasing safety and autonomy of the cars, which cannot be achieved without them. We expect that our estimation scheme will find its way to real applications in that framework. 6. Conclusion A new scheme for robust tire/road friction force estimation is proposed. An RBF neural network is used to compensate negative influence of the friction model uncertainties to the accuracy of the force estimation. By splitting the nonlinearity of the friction model in two parts it is achieved that RBF network input comprises only one
ARTICLE IN PRESS 456
J. Matusˇko et al. / Engineering Applications of Artificial Intelligence 21 (2008) 442–456
variable, which enables that the centers and the spreads of the hidden layer neurons can be set off-line in a very simple manner. The output layer weights are adapted on-line by an adaptation law derived applying Lyapunov stability analysis. Such an approach ensures fast and robust estimation of the tire/road friction force even in cases when wide and abrupt changes of the friction parameters occur, which is verified by simulation experiments on the one-wheel friction model. The authors are hopeful to do experimental testing in their future work. Acknowledgments Support from the Ford Motor Company and the Ministry of Science, Education and Sports of the Republic of Croatia is gratefully acknowledged. References Bakker, E., Pacejka, H.B., Lidner, L., 1987. A new tire model with an application in vehicle dynamics studies. SAE Paper. Baric´, M., Petrovic´, I., Peric´, N., 2005. Neural network based sliding mode control of electronic throttle. Engineering Applications of Artificial Intelligence 18 (8), 951–961. Bozˇic´, A., Petrovic´, I., Peric´, N., Matusˇ ko, J., 2002. Experimental investigations of the rubber-asphalt sliding pair. In: ASME International Congress and Exp.—Symposium on Advanced Automotive Technologies, New Orleans, USA. Burckhardt, M., 1993. Fahrwerktechnik: Radschlupf-regelsysteme. Vogel, Wuertzburg, Germany. Canudas de Wit, C., Horowitz, R., 1999. Observers for tire/road contact friction using only wheel angular velocity information. In: 38th IEEE Conference on Decision and Control, Phoenix, USA. Canudas de Wit, C., Tsiotras, P., 1999. Dynamic tire friction models for vehicle traction. In: 38th IEEE Conference on Decision and Control, Phoenix, USA. Canudas de Wit, C., Olsson, H., Astro¨m, K.J., Lischinsky, P., 1995. A new model for control of systems with friction. IEEE Transactions on Automatic Control 40 (3), 419–425. Chen, T., Chen, H., 1995. Approximation capability to functions of several variables, nonlinear functionals, and operators by radial basis function neural networks. IEEE Transactions on Neural Networks 6 (4), 904–910.
Chen, T., Chen, H., Liu, H., 1995. Approximation capability in c(r) n by multilayer feedforward networks and related problems. IEEE Transactions on Neural Networks 6, 25–30. Deur, J., 2001. Modeling and analysis of longitudinal tire dynamics based on the lugre friction model. In: 3rd IFAC Workshop Advances in Automotive Control, Karlsruhe, Germany. Deur, J., Asgari, J., Hrovat, D., 2004. A 3d brush-type dynamic tire friction model. Vehicle System Dynamics 42, 133–173. Dieckmann, T., 1992. Assessment of road grip by way of measured wheel. In: XXIV FISITA Congress, London. Gustafsson, F., 1997. Slip-based tire-road friction estimation. Automatica 3 (6), 1087–1099. Horn, R.A., Johnson, C.R., 1999. Matrix Analysis. Cambridge University Press, Cambridge. Hozwarth, F., Eichhorn, U., 1993. Non-contact sensors for road conditions. Sensors and Actuators (37–38), 121–127. Kiencke, U., Daiss, A., 1994. Estimation of tyre friction for enhanced abssystems. In: International Symposium on Advanced Vehicle Control, AVEC 1994, Tokyo, Japan. Kiencke, U., Nielsen, L., 2000. Automotive Control Systems. Springer, Berlin. Matusˇ ko, J., Petrovic´, I., Peric´, N., 2002. Application of extended Kalman filter for road condition estimation. In: 10th International Power Electronics and Motion Control Conference, EPE-PEMC 2002, Cavtat and Dubrovnik, Croatia. McCulloch, W.W., Pitts, W., 1943. A logical calculus of the ideas imminent in nervous activity. Bulletin of Mathematical Biophysics 5, 115–133. Moody, J., Darken, C., 1989. Fast-learning in networks of locally-tuned processing units. Neural Computations 1, 281–294. Pal, C., Higiwara, I., Morishita, S., Inoue, H., 1994. Application of neural networks in real time identification of dynamic structural response and prediction of road friction coefficient m from steady state automobile response. In: International Symposium on Advanced Vehicle Control, AVEC 1994, Hiroshima, Japan. Park, J., Sandberg, I.W., 1991. Universal approximation using radialbasis-function networks. Neural Computation 3 (2), 246–257. Slotine, J.J., Li, W., 1990. Applied Nonlinear Control. Prentice-Hall, Englewood Cliffs, NJ. Umeno, T., Ono, E., Asano, K., Ito, S., Yamamoto, M., Yasui, Y., Sawada, M., 2002. Estimation of tire–road friction using wheel speed vibration. In: International Symposium on Advanced Vehicle Control, AVEC 2002, Hiroshima, Japan. van Zanten, A., Ruf, W.D., Lutz, A., 1989. Measurement and simulation of transient tire forces. SAE Technical Paper Series, Paper 890640. van Zanten, A.T., Erhardt, R., Leindesfeind, K., Pfaff, G., 1998. Vdc systems development and perspective, In: SAE International Congress and Exposition, Detroit, USA.