V Ciarla

C Canudas De Wit

F Quaine

V Cahouet

Oscillation annealing in Electronic Power Steering (EPS) systems

This document is a working report concerning the Task 5.1 of the project VolHand 09 VTT 14 and refers to the paper [1]. This task is competence of the GipsaLab. This report presents several aspects of modelling, observation and control towards a new generation of Electrical Power Steering(EPS) systems. In particular authors design an optimal control to reject oscillations of the steering column. Authors also revisited the LuGre tire dynamic friction model by improving the transient behaviour between the sticking phases and the dynamic ones. Simulation of the proposed control are shown at the end of the document.

Introduction

According to the document ANR 09 VTT VOLHAND Doc B, presented on October 2011, the research team in Gipsa-Lab is involved into the task 5.1 (T5.1) for the design of the control law to eliminate the high frequency oscillations of the steering column. To attempt to this job, the following steps have been done:

1. The study of a mechanical model of the Electronic Power Steering (EPS) system and design of an LQR regulator to reject the typical oscillations, due to the torsion of the steering wheel (Paragraph 2)

2. The design of an estimation of the full set of state variables as well as of the exogenous torques is carried out.

3. Modelling of the tyre/road contact friction, in order to test the mathematical model uder the most realistic conditions (Paragraph 3).

4. Study of the amplification curves, in order to provide the correct steering assistance to the driver.

The general architecture that authors propose to study the EPS system is shown in Figure 1. This report concerns the block 1 (observer/control), and the block 2 (tire-friction model) of Figure 1, with the additional contribution that the observer also estimate the driver's and load torque. The design of the specific reference model (block 3) is under current investigation (tasks 5.2 -5.3 and 5.4).

Mechanical Model

Oscillation annealing 2.1 Column model

To begin with, it is essential to define the mechanical system to study. Figure 2 shows an explanatory schema of the system. From left to right, it is possible to find:

• the driver's exerted torque, τ v , as an exogenous input;

• the steering wheel, represented as an inertia J v with a given local viscosity B v ;

• the stiffness of the steering column represented by the coefficient k; • the assistance motor, similarly to the steering wheel, modelled as an inertia J m and a local damping B m ;

τ v J v B v θ v k θ s = θ m /N 2 J m B m N 2 N 1 u τ a θ m
• gear ratio N 2 that connects the assistance motor to the steering column;

• gear ratio N 1 that connects the the steering column to the wheels;

• the control inputs u and the exogenous input τ a , given by the road/tyre friction torque.

The mechanical equations governing the system explained above are inspired from [START_REF] Canudas De Wit | Control design for an electro power steering system: part I The reference model[END_REF] and describe the system as follows:

J v θv = τ v -k(θ v -θ s ) -B v θv (1)
J T θs = -k(θ s -θ v ) -N 2 2 B m θs + τ a N 1 + N 2 u (2)
with

J T = J c + N 2 2 J m + J w N 2 1
. The constants of the model are reported in Table 1, while θ m is the motor angle [r ad], θ v is the steering wheel angle [r ad] and θ s is the driving shaft angle after the gear ratio, i.e.:

θ s = θ m N 2 (3)

State-space representation

Being the first objective to compensate the oscillations due to the torsion force

F t = k(θ v -θ s )
of the steering wheel, the synthesis of an optimal Linear Quadratic Regulator (LQR) controller is

Symbol Description Value

J v

Steering wheel inertia 0.025 kgm

x =    x 1 x 2 x 3    =    θv θs θ v -θ s    (4)
The model can be formulated into the following state-space form:

ẋ = Ax + Bu + Gw (5)
with the disturbances vector w T = τ v , τ a T and following state matrices

A =      - B v J v 0 -k J v 0 - N 2 2 B m J T k J T 1 -1 0      (6) B =     0 N 2 J T 0     (7) G =     1 J v 0 0 1 N 1 J T 0 0     = g 1 g 2 (8)

Open loop model

Once the space state representation of the mechanical model has been completely defined, the interest is to observe how the model responds to a variation of the torque exerted by the driver on the steering wheel. For that to be done, authors compute the transfer function of the steering wheel's angular velocity θv in relation to the driver's exerted torque τ v . Such transfer function is easily calculated from the state-space matrices as follows:

G ol = θv τ v = C(sI -A) -1 G 1 (9)
where C = 1, 0, 0 is the output vector. The frequency response of the open loop system, described in Eq. (9), is shown in Figure 3. The system has a significant peak for a frequency of about 12 Hz, that might cause important oscillations on the steering wheel. The goal of the control design is to compensate these oscillations, on order to improve the driver's comfort and safety.

Optimal control based on full state-feedback

To compensate the oscillations that might appear in the steering wheel, an optimal LQR controller is implemented. This controller is calculated so that the state-feedback law u = -K x (with K the state-feedback gain) minimizes a cost function given by the expression:

J (x, u) = ∞ 0 (x T Q x + u T Ru) dt ( 10
)
where the constant matrices Q and R are chosen in order to penalize, respectively, the states and the control inputs in the cost function. To solve this optimization problem it is necessary to solve the following algebraic Riccati equation: whose solution P is used to calculate the state-feedback gain: K = R -1 B T P. To synthesize this controller the first step to do is to select the adequate penalization matrices for our specific problem. Three possible choices are encountered:

A T P + PA -PBR -1 B T P + Q = 0 (11)
1. penalization on the variation of the torsion and the torsion rate;

2. penalization on the torsion rate;

3. penalization on the torsion.

Penalization on the variation of the torsion and the torsion rate

In this case, the cost function to minimize has the form:

J (x, u) = ∞ 0 (q 1 (x 1 -x 2 ) 2 + q 2 x 2 3 + Ru)dt (12)
To respect the restriction on the matrix Q (it has to be positive or semi-definite positive and symmetric) an easy way to build it is to take an "output" matrix representing the state variables to penalize. In this case:

C q = q 1 -q 1 0 0 0 q 2 ( 13
)
and the matrix penalizing the states is computed as follows:

Q = C T q C q =    q 1 -q 1 0 -q 1 q 1 0 0 0 q 2    (14)
Thus, the value R of the penalization on the control input is normalized and set to 1. Consequently, the values of the parameters q 1 and q 2 that define the matrix Q must be sufficiently high, such that the states are correctly penalized in the cost function. The chosen values are, respectively, 3 and 12. Therefore:

Q =    3 -3 0 -3 3 0 0 0 12    ; R = 1 (15)
Once the state-feedback gain has been calculated, the closed loop transfer function is given by the following equation:

G cl = θv τ v = C(sI -(A -BK)) -1 G 1 (16)
with (A -BK) the closed-loop state space matrix. To prove the performances of the state-feedback 4 the resonance peaks that might cause undesirable oscillations have been eliminated thanks to the optimal linear control computed and implemented. However, the gain for low frequencies has been sensibly increased. This fact might affect the performances of the system.

Penalization on the variation of the torsion rate

Next, it is also possible to penalize only the torsion rate. This implies to minimize the following cost function:

J (x, u) = ∞ 0 (q 1 (x 1 -x 2 ) 2 + Ru)dt (17)
For this reason, the matrix used to penalize the states in the cost function is:

Q =    q 1 -q 1 0 -q 1 q 1 0 0 0 0    (18)
with q 1 = 7 and R = 1. It can be observed in Figure 5 that the frequency response in this case is very close to that obtained in the case of penalizing the variation of the torsion and the torsion rate.

Penalization on the torsion

Finally, the objective is to penalize only the torsion; that is, to have a cost function of the form:

J (x, u) = ∞ 0 (q 2 x 2 3 + Ru)dt (19)

Full State Penalization Penalization on x1-x2

Figure 5: Full state penalization (solid) vs torsion velocity penalization Thus, following the same reasoning as in the previous cases,

C q = 0 0 q 2 (20)
the corresponding covariance matrix is

Q =    0 0 0 0 0 0 0 0 q 2    (21)
with q 2 = 200 and R = 1. Figure 6 shows that the fact of only penalizing the torsion affects slightly

Full State Penalization Penalization on x3

Figure 6: Full state penalization (solid) vs torsion penalization on the frequency response since there is a small peak where the oscillations occurred for the open loop case. However, the magnitude of such a peak is not at all significant for the performance of the closed-loop system.

Optimal control based on acceleration

To evaluate the performances of the controller in terms of comfort improvement, it is common to consider the transfer function between the steering wheel acceleration θv and the driver's exerted torque. Such a transfer function is calculated from the state-space matrices as follows:

G ol = θv τ v = T (sI -A) -1 G 1 ( 22
)
where T = -

B v J v , 0, -k J v
. The corresponding frequency response is shown in Figure 7 the system has a significant peak for a frequency of about 12 Hz that might cause important oscillations on the steering wheel that would rather be avoided so as to improve the driving comfort. In Figure 7, it is possible to observe the benefits of the proposed controller. The resonance peak that might cause undesirable oscillations have been eliminated thanks to the proposed optimal linear control. Note also that the low-frequency gain has kept unchanged.

Optimal "output" feedback controller

All the calculations that have been previously carried out assumed that the whole set of state-space variables was measurable, that is, it was possible to obtain exact values of those variables. On the contrary, the actual mechanical model is not provided with sensors to measure all the variables, because of the important extra cost for the installation of these devices. On a common power steering problem, the variables that are measured are those concerning the assistance motor, either θ s or θs . We assume then θs as the output of the system, that is, C = (0, 1, 0). The other state variables are not directly measurable so it is necessary to calculate a state observer to reconstruct the full state so as to implement the full state-feedback computed in previous sections.

Being our system ẋ = Ax + Bu + Gw the observer will be:

ẋ = (A -L C)x + Bu + Gw + L y (23)
In Eq. ( 23), L is the observer gain and is chosen via pole placement method or via the loop recovery strategy. Thus L is chosen so as to have a fast error dynamics and hence a fast convergence of the observer towards the real states.

Simulation results

In order to see the performances of the proposed controller, following simulations are carried out using as input to the system the driver's exerted torque profile shown in Fig. 8-a.

Open loop behaviour: the first simulation concerns the system in open loop. As shown at the Figure 8-b, once the steering wheel is released by the driver at the time instant 16 s, the steering wheel suffers significant oscillations that would cause an undesirable driving feeling and might even be dangerous when we are in driving situations at fast speeds.

Closed-loop behaviour: the optimal linear controller computed before has the aim to eliminate the oscillations found in open loop thanks to the assistance motor. As shown at the Figure 8-c, the oscillations derived from the release of the steering wheel have been satisfactorily compensated by the optimal controller. The behavior of the observer is shown by Figure 8-d. From this figure we see that the performance of the observer is completely satisfactory.

Road/tyre contact friction torque

The model used until now takes into account all the elements that go from the steering wheel to the steering shaft. Load torques exerted by the tyre/road contact need to be considered explicitly in the model. The contact friction force is an element of the driving that is very important since it is present in all the different driving situations that a driver may find. To be able to reproduce this contact friction force it is necessary to have a model of the vehicle. However, the fact of developing and implementing a complete vehicle model entails a huge complexity. Thus, and as it is not essential for the development of our simulations, the model chosen is a simplified model of the tyres assuming that the normal force on them is equal to a quarter of the total weight of the vehicle. Figure 9 shows a tyre as well as both the normal and the lateral velocities (v x and v y ) that might appear when driving. This lateral force that takes part of the model when turning is the reason why the so-called slip angle β appears. Thus, the slip angle has to be computed. It depends from

v y v v x β

Figure 9: Forces action on the tire when turning. Slip angle β show the direction of the vehicle resulting velocity the speed of the vehicle and the radius of curvature given by the curve traced by the vehicle in a given moment, that is: β(t) = f (v, ρ). The radius of curvature can be considered as being inversely proportional to the steering wheel angle: an infinite radius of curvature means the vehicle is going straight and hence the steering wheel angle is zero. The slip angle is considered to be zero in the cases when the vehicle is stopped or it is not in a cornering situation. Besides, the slip angle is proportional to the speed of the vehicle and the steering wheel angle. Furthermore, the values of the slip angle will not be very high as the situations that will be handled are mainly low speed manoeuvres that would not possibly exceed 30 km/h. Thus, the slip angle does not normally go above 5 d e g.

Dynamic tyre friction model

The dynamic model used so as to represent the forces and torques caused by the contact of the tyres and the road is an extension to the LuGre dynamic friction model presented in the paper [START_REF] Velenis | Dynamic Tire Friction Models for Combined Longitudinal and Lateral Vehicle Motion[END_REF]. The objective of this friction model is to capture the different phenomena presents on a real driving situation so as to obtain simulations as trustworthy as possible. The most important phenomena are described below:

• Self-alignement torque: it is the torque that appears as the vehicle rolls along that tends to rotate the wheels around its vertical axis in order to make them return to the straight position.

• Sticking torque: it is the torque that opposes the movement of the tyres when turning them around the vertical axis. It is especially important for low speeds of the vehicle

The dynamic LuGre Friction model proposed in [START_REF] Velenis | Dynamic Tire Friction Models for Combined Longitudinal and Lateral Vehicle Motion[END_REF][START_REF] Canudas De Wit | Dynamic friction models for road/tire longitudinal interaction[END_REF] derives from a distributed friction model. It is described by four differential equations, which capture the average behaviour of the internal friction states. To introduce this model, let zi (t) (i = x, y), ẑy (t) denote the internal states of the system:

żi (t) = v r i -C 0i (v r )z i (t) -κ ss i |ωr|z i (t) (i = x, y) (24) ży (t) = G F n L v r y -C 0 y (v r )ẑ y (t) -ν ss |ωr|ẑ y (t) + |ωr| L zy (t) (25)
where v r = [ v r x v r y ] T is the vector of the relative speeds:

v r x = ωr -vcos(β) (26) v r y = -vsin(β) (27)
with ω that corresponds to the angular velocity of the wheel of radius r, while v is the velocity of the vehicle and β is the slip angle.

The scalar function C 0i (v r ) (with i = x, y) is peculiar of the LuGre model and is given by:

C 0i (v r ) = λ(v r )σ 0i µ 2 ki (28) with λ(v r ) = ||M 2 k v r || g(v r ) (29)
In Eq. ( 29) it is possible to recognize the the matrix of the kinetic friction coefficients M k = µ k x 0 0 µ k y > 0 for the motion along the x and the y directions, respectively; note also that each parameter µ ki in Eq. ( 28) corresponds to one element on the main diagonal of matrix M k . We find also the function g v r , given by:

g v r = M 2 k v r M k v r + M 2 s v r M s v r - M 2 k v r M k v r e - v r v s γ ( 30
)
with M s = µ s x 0 0 µ s y > 0 that is the matrix of static friction coefficients.

To evaluate the constants κ ss i and ν ss , we do the hypothesis that the steady-state solution of the lumped model is the same with the steady-state solution of the distributed one. The complete expression of these parameters is given from:

κ ss i = 1 |ωr| v r i zss i -C 0i (v r ) (31)
ν ss = 1 |ωr| 1 ẑss y Gv r y F n L + |ωr|z ss y L -C 0 y (32)
with zss i and ẑss y , that are the steady-states of the system; the explicit expression for these parameters is given in the Appendix. The auto-aligning torque can be written in terms of the mean states zy (t) and ẑy (t) as follows:

M sel f -al i g = F n L σ 0 y 1 2 zy (t) -ẑy (t) + σ 1 y 1 2 ży (t) -ży (t) + σ 2 y 1 2 v r y - G F n L (33)
To calculate the sticking torque we introduce the last equation of the LuGre model: where the angular velocity of the wheel rim is φ = θs /N 1 , and and the function

żz = φ - σ 0z | φ| g z ( φ) z z (t) (34)
g z φ = µ kz + µ sz -µ kz e -φ φs 2 (35)
where µ kz and µ sz are, respectively, the kinetic and static friction coefficients across the z-axis, while φs is the Stribeck velocity.

The sticking torque can be evaluated as follows:

M s t icking = -LF n (σ 0z z z (t) + σ 1z żz (t) + σ 2z φ) (36)
Hence, the total torque generated by the contact of the tyres and the road is:

τ a = M sel f -al i g + M s t icking (37)
The parameters used in this model are reported in Table 2.

Friction model improvements

During the analysis of equations of the LuGre model, there are some parameters that are crucial to the behaviour of the sticking torque and can be set in order to improve the realistic behaviour of the model.

The approximate torque exerted on the wheels

This parameter F n L can be chosen experimentally by measuring the necessary torque to be exerted in order to start moving the wheels without assistance when being at a standstill. As by now these measures are not available, it can be assumed that a torque of about 50 Nm appears on the wheels before starting to move. Since the sticking torque is introduced directly on the steering column on the model presented, the motor/steering colon gear ratio has to be taken into account. Hence, F n L should have an approximate value of F t = 50 N 2

Coefficients of the sticking torque

This coefficient σ 0z affects directly on the settling time of the sticking torque. It is important to notice that very fast dynamics make the sticking torque vary in steps. This includes very high frequency components that are likely to cause oscillations and instabilities. Thus, it is advisable to choose a value for σ 0z low enough so as to avoid too fast dynamics. In this case the value chosen was σ 0z = 0.0001. The value σ 1z is chosen conveniently to avoid overshoot in the dynamics of the sticking torque.

Dependency on vehicle's speed

One possible improvement of the previous model concerns the the dependence of the sticking torque to the vehicle velocity, v. In fact, it can be proved experimentally that, as the speed of the vehicle grows, the self-alignment torque dominates over the sticking torque, and inversely. The previous model in its actual form, does not respect this observation. It is then necessary to weight the sticking torque as a function the velocity of the vehicle. Thus, a possible modification along these observations is:

M s t icking = -LF n (σ 0z z z + σ 1z żz + σ 2z φ)e -|v|/v k (38)
where v k is a positive constant.

Simulation results

It is also important to check if the results follow the expected logic when we vary the vehicle's speed.

In order to do this, several simulations were carried out for velocities going from 0 to 30 km h -1 and over a range time of several seconds, in order to see the effects at steady-state. The choice of this speed range is due to the fact that an EPSs operate in this range. The results obtained from those simulations are shown in Figure 10 Results are as expected: the sticking torque decreases exponentially as the velocity increases, so as the contribution of the self-alignment is more important as the speed increases.

A Appendix

A.1 Constants of the LuGre model

To evaluate the constants κ ss i (for i = x, y)in Eq. ( 31) and ν ss in Eq. (32), we have to know following steady-state equations of the system:

zss i = 1 F n L 0 z ss i (ζ)f n (ζ)dζ = = 1 F n C 1y α 1 ζ L 2 2 -C 1y α 1 C 2y 2 -C 2y C 2y + ζ L e
z ss i (ζ) = C 1i (1 -e - ζ C 2i ) (41)
f n (ζ) is the trapezoidal load distribution:

f n (ζ) =    α 1 ζ, for 0 ≤ ζ ≤ ζ L f ma x , for ζ L ≤ ζ ≤ ζ R α 2 ζ + β 2 , for ζ R ≤ ζ ≤ L (42)
G is the following function:

G = L 0 f n (ζ)ζdζ = = α 1 3 ζ 3 L + f ma x 2 (ζ 3 R -ζ 3 L ) + α 2 3 (L 3 -ζ 3 R ) + β 2 2 (L 2 -ζ 2 R ) (43

Figure 1 :

 1 Figure 1: General architecture of the EPS system: block 1 concerns the observer and the control; block 2 includes the tyre-friction model and block 3 includes the reference model.

Figure 2 :

 2 Figure 2: Mechanical model of the EPSs

Figure 3 :

 3 Figure 3: Frequency response of the transfer function G ol in Eq. (9)

Figure 4 :

 4 Figure 4: Open loop vs closed loop frequency response

Figure 7 :

 7 Figure 7: Open loop (blue) vs closed loop frequency response. Weighting parameters are: q 1 = 3 and q 2 = 12, r = 1.

Figure 8 :

 8 Figure 8: (a) Real driver's steering torque used for the simulation. (b) Profile of the steering wheel speed in open loop. (c) Profile of the steering wheel speed in closed loop. (d) Estimation of the steering wheel speed.

:•

  Figure 10(a) shows the different curves obtained for the sticking torque for different driving speeds, while Figure 10(b) shows the different curves obtained for the self-alignment torque.

+-+++

  C 1y f max ζ L -ζ R + C 2y (e -C 1y α 2 C 2y C 2y + L e -L C 2y -C 2y C 2y + ζ R e -C 1y β 2 (L -ζ R ) + C 1y α 2 2 (L 2 -ζ R 2 ) + C 1y C 2y β 2 e (σ 0 y zss y + σ 2 v r y ) -M ss z F n Lσ 0 y -Gσ 2 y v r y F n Lσ 0 y (40)z ss i (ζ)is the steady-state function of the distributed model:

Table 1 :

 1 Constant parameters of the EPSs proposed. At this purpose, it is necessary to use a state-space representation of the system. The following choice of the state variables is made:

Table 2 :

 2 Constant parameters of the LuGre model

	Symbol	Value	Description
	r	wheel radius	0.38 m
	v	vehicle speed	range from 0 to 30 km h -1
	ω	angular velocity of the wheel	ω = v/r rad s -1
	µ k x	kinetic friction coeff. x-axis	0.75
	µ k y	kinetic friction coeff. y-axis	0.75
	µ kz	kinetic friction coeff. z-axis	0.76
	µ s x	static friction coeff. x-axis	1.35
	µ s y	static friction coeff. y-axis	1.40
	µ sz v s φs	static friction coeff. z-axis Stribeck relative velocity Stribeck relative velocity	0.91 3.96 m s -1 74 rad s -1
	γ	steady state constant	1
	σ 0 y σ 1 y σ 2 y	normalized rubber stiffness normalized rubber damping normalized viscous relative damping	6000 N m -1 0.3568 N m -1 s -1 0.0001 N m -1 s -1
	L	patch length	0.15 m
	ζ L	left patch length	0.0030 m
	ζ R	right patch length	0.1155 m
	F ma x	max value of normal load distribution	1900 N
	F n α 1 α 2	normal value of normal load distribution coeff. for Self-Align. torque coeff. for Self-Align. torque	249.37 N 63000 N m -1 -55000 N m -1
	β 2	coeff. for Self-Align. torque	8260 N
	σ 0z	coeff. for Stick. torque	20
	σ 1z	coeff. for Stick. torque	0.0023
	σ 2z G	coeff. for Stick. torque load distribution function	0.0001 16.83 Nm 2
	κ ss i ν ss function used to approx. steady behavior function used to approx. steady behavior	11.9 -0.8

  ) assuming f n (ζ) as in Eq. (42).M ss z is the steady-state auto-aligning moment of the distributed model: L α 1 σ 0y C 2y -C 2y + ζ L e C 1y f max σ 0y C 2y C 2y + ζ L e -C 1y β 2 σ 0y C 2y C 2y + L e

		L
	M ss z =		σ 0y z ss y (ζ) + σ 2y v ry f n (ζ)	L 2	-ζ dζ =
		0
	= C 1y α 1 σ 0y 2 C 2y	3 -C 2y	3	2 ζ L C 2y	+	ζ L C 2y 2 2 + 2 e	-	ζ L C 2y	-	1 3	α 1 ζ L	3 (C 1y σ 0y + σ 2y v ry )+
	-	1 2	C 1y C 2y , -	ζ L C 2y	+	L α 1 4	(C 1y σ 0y ζ L	2 + v ry σ 2y ζ L	2 )+
	+ -	ζ L C 2y -C 2y C 2y + ζ R e	-	ζ R C 2y	+
	+	C 1y f max σ 0y 2	(ζ L	2 -ζ R	2 ) +	f max σ 2y v ry 2	(ζ L	2 -ζ R	2 )+
	-	C 1y L f max σ 0y 2	(ζ L + ζ R ) -	L f max σ 2y v ry 2	(ζ L -ζ R )+
	-	1 2	C 1y C 2y L f max σ 0y e	-	ζ L C 2y -e	-	ζ R C 2y	+	C 1y α 2 σ 0y ζ R 3	3	+
	-	C 1y L 3 α 2 σ 0y 12	+	C 1y β 2 σ 0y ζ R 2	2	-	L 3 α 2 σ 2y v ry 12	+
	-C 1y α 2 σ 0y C 2y	3	2 L C 2y	+	L 2 C 2y	2 + 2 e	-L C 2y -C 2y	3	2 ζ R C 2y	+	ζ R C 2y 2 2 + 2 e	-	ζ R C 2y	+
	+ σ 2y v ry ζ R	2 α 2 3	ζ R +	β 2 2	-L C 2y -C 2y C 2y + ζ R e	-	ζ R C 2y	+
	+	1 2	C 1y L α 2 σ 0y C 2y C 2y + L e	-L C 2y -C 2y C 2y + ζ R e	-	ζ R C 2y	+
	-	C 1y L β 2 σ 0y ζ R 2		-	L β 2 σ 2y v ry ζ R 2	-	C 1y L α 2 σ 0y ζ R 4	2	-	L α 2 σ 2y v ry ζ R 4	2	+
	+	1 2	C 1y C 2y L β 2 σ 0y e	-L C 2y -e	-	ζ R C 2y
															(44)

V. Ciarla, C. Canudas de Wit, F. Quaine, V. Cahouet-GIPSA-Lab -Report VOLHAND Task 5.1

Conclusions

This document presented simulation results of a detailed model of an EPS system including a dynamic model that reproduces the physical phenomena involved in driving. From an adequate mechanical model of the steering column, the natural oscillations that exist were eliminated. Furthermore, a satisfactory estimator of the exogenous torques was developed. Next, a dynamic friction model based on the LuGre model was successfully implemented in the simulation obtaining physically logical results.


V Ciarla

C Canudas De Wit

F Quaine

V Cahouet

Exogenous input estimation in Electronic Power Steering (EPS) systems

This document is a working report concerning the Task 5.1 of the project VolHand 09 VTT 14 and refers to the paper [1]. This task is competence of the GipsaLab. This report presents several aspects of modelling, observation and control towards a new generation of Electrical Power Steering(EPS) systems. In particular authors design an observer to estimate the driver applied torque (steering wheel torque) and the load torque (tire/ground contact friction). Simulation of the proposed control and observer are shown at the end of the document.

Introduction

According to the document ANR 09 VTT VOLHAND Doc B, presented on October 2011, the research team in Gipsa-Lab is involved into the task 5.1 (T5.1) for the design of the control law to eliminate the high frequency oscillations of the steering column. To attempt to this job, the following steps have been done:

1. The study of a mechanical model of the Electronic Power Steering (EPS) system and design of an LQR regulator to reject the typical oscillations, due to the torsion of the steering wheel.

2. The design of an estimation of the full set of state variables as well as of the exogenous torques is carried out (Paragraph 2)

3. Modelling of the tyre/road contact friction, in order to test the mathematical model uder the most realistic conditions.

4. Study of the amplification curves, in order to provide the correct steering assistance to the driver.

The general architecture that authors propose to study the EPS system is shown in Figure 1. This report concerns the block 1 (observer/control) to estimate the driver's and load torque.

Mechanical Model

Tyre/Road Friction Model Reference Model

φ β V τ v τ a Observer -K x y u ass

Extended state-space representation

Let consider that both exogenous torques are slowly time-varying: τa = τv ≈ 0. The extended state-space representation of the EPS system, shown in [START_REF] Illàn | Oscillation annealing and driver/tire load torque estimation in electric power steering systems[END_REF], with z is now defined as follows:

z = θv , θs , θ v -θ s , τ v , τ a T ( 1
)
in which both exogenous torques are added as state variables, is give by

ż = Ãe z + Be u (2)
The state matrices for this system are shown below:

Ãe = A G 2×3 2×2 ; Be =    B 0 0    ; C e = C 0 0 (3)
If only the motor velocity θs is used as a available output, then the observability matrix for this new extended system has rank 4 in spite of dealing with a system of dimension 5. To avoid this problem, it is possible to measure a second signal: the torsion force,

y 2 = k(θ v -θ s ) = k x 3 (4)
In this case, the output matrix C e of the extended system becomes:

Ce = 0 1 0 0 0 0 0 k 0 0 . ( 5
)
Computing the new resulting observability matrix, it can be easily check that it has now full rank.

The system is then observable and it is concluded that both exogenous torques can be correctly estimated. The poles of the observer's dynamic (matrix Ãe -L e Ce ) should be fast enough compared to those resulting from the linear controller. The performances of the observer have been tested in simulations. The results obtained for the estimation of τ v and τ a are shown in Fig. 2(a) and Fig. 2(b). The maximum value reached by the observation error can be reduced at the price of increasing the observer gain. Measurement noise will limit at certain point the maximum possible value of the observer gain.

Conclusions

This document presented simulation results of a satisfactory estimator of the exogenous torques acting on the EPS system.

Figure 1 :

 1 Figure 1: General architecture of the EPS system: block 1 concerns the observer and the control; block 2 includes the tyre-friction model and block 3 includes the reference model.

Figure 2 :

 2 Figure 2: (a) Real driver's torque (solid) vs. estimated driver's torque (dotted). (b) Real load torque (solid) vs. estimated load torque (dotted).

V. Ciarla, C. Canudas de Wit, F. Quaine, V. Cahouet-GIPSA-Lab -Report VOLHAND Task 5.1


V Ciarla

C Canudas De Wit

F Quaine

V Cahouet

Preliminary study on Electronic Power Assistance Steering (EPAS) systems

. The basic concepts of the Electric Power Steering (EPS) systems with a realistic model for the friction contact, that acts on the wheels are discussed. A relation between the assistance and the driver's torque is provided, under the hypothesis of a position-oriented control of the movement and the Stevens' power law. Finally, the simulation results proposed at the end of this paper validate the shape of the booster curves and are in accord with the initial hypothesis.

Introduction

The assistance provided by the electric motor of EPS systems reproduces the hysteresis cycle of the hydraulic valve, as shown in Fig. 1: the amplification torque U grows, according to the value of the driver's torque τ v and the vehicle speed V , then saturates. If the driver's torque changes in direction, we have the downhill phase of the amplification torque. During this phase, the amplification torque decreases, until reaching the lower saturation values. In literature, we can find many stationary models for the EPS system ( [START_REF] De | Fun-to-drive by feedback[END_REF], [START_REF] Canudas De Wit | Stability analysis via passivity of the lateral actuator dynamics of a heavy vehicle[END_REF], [START_REF] Claeys | Modeling and control of steering actuator for heavy duty vehicle[END_REF], [START_REF] Zhang | Modeling and characteristic curves of electric power steering system[END_REF], [START_REF] Chabaan | Optimal control and gain scheduling of electrical power steering systems[END_REF], [START_REF] Jiang | Low-pass filter based automotive eps controller and comparative full-vehicle tests[END_REF], [START_REF] Tsung Hsien | Design of control logic and compensation strategy for electric power steering systems[END_REF], [START_REF] Burton | Innovation drivers for electric power-assisted steering[END_REF], [START_REF] Hoon | Control logic for an electric power steering system using assist motor[END_REF]), where the hysteresis cycle is approximated by a single curve. Nevertheless, these investigations do not provide ground foundations (other than these curves reproduce the behaviour of the older power steering systems, based on hydraulic valves), neither evidence that they are optimal for the driver in any sense. To justify the optimality of these curves (or to develop new ones), we propose a model for the steering column and the assistance, by introducing a relationship between the assistance and the driver's exerted torque. This relationship is based on the Steven's power law. The simulations at the end of this report shows some assistance curves, that are very similar to those shown in literature.

Torques acting on the EPS system

Works in [START_REF] Illán | Oscillation annealing and driver/tire load torque estimation in electric power steering systems[END_REF] show that it is possible to compensate the torsion torque in the steering column with a inner control law. Under this hypothesis, it is possible to describe the dynamics of the EPS system, as follows:

J θ + B θ + τ a N 1 = τ v + U (1)
J = J v + J T = total inertia at the steering wheel level

B = B v + N 2 2 B m + σ 2 N 2 1
= the total viscosity of the system at the steering wheel level

τ a = friction torque τ v = driver's exerted torque
At stand still, the main force opposing to the driver's effort is the contact friction. In this situation, it is possible to assume that the assistance U acts to compensate, first of all, a part of the contribution of the load torque τ a in Eq. ( 1). Indeed at steady-state conditions:

τ a N 1 = τ v + U (2)
In order to design the assistance amplification, the hypothesis that the driver's perception of the load torque obeys to the Stevens' power law ([14, 8, 1]) is done:

δ τ a N 1 n = τ v ( 3
)
where:

δ = the intensity magnitude of the external stimulus n = the rate of growth of the sensation of the stimulus and depends on the sensory modality of the driver (e.g. perception of force or perception of position).

The value of δ is tuned to obtain a realistic value of the amplification torque, while to choose the value of n, it is mandatory to refer to the literature, concerning some neuro-muscular studies.

During the steering movement, some authors ( [START_REF] Hoult | A neuromuscular model featuring co-activation for use in driver simulation[END_REF][START_REF] Pick | A mathematical model of driver steering control including neuromuscular dynamics[END_REF]), show that drivers prefer exerting their control, basing on the transfer function between the friction contact load and the steering angle (position-oriented control), rather than the transfer function between the friction torque and the steering torque (torque-oriented control). This can be explained by the fact that the neuro-muscular system is better suited to control the steering angle than the exerted torque, thanks to the gamma (γ) motor neurons, that adapt the length of the spindles according to the hand wheel angle [START_REF] Pick | A mathematical model of driver steering control including neuromuscular dynamics[END_REF].

Under the hypothesis of a control-oriented control, we will use n < 1.

Relation with the assistance torque

To obtain a relation between the assistance torque and the driver's torque in steady-state, Eq. ( 2) is substituted in Eq. ( 3), as follows:

δ τ v + U n = τ v (4)
By solving for U, it is possible to obtain:

U = n τ v δ -τ v (5)
When the dynamic terms of inertia and viscosity are included, the same closed-form solution for the assistance U is not reachable, but it follows:

J θ + B θ + δ τ a N 1 n = τ v (6)
Let consider Eq. ( 1), we can express the assistance torque U, as the following difference:

U = τ a N 1 -δ τ a N 1 n (7)
In this case, it is not possible to retrieve the corresponding value for the assistance torque in closed form, because it is necessary to solve Eq. ( 6) to obtain the corresponding values in Eq. [START_REF] Claeys | Modeling and control of steering actuator for heavy duty vehicle[END_REF]. For this reason, it is necessary to introduce the optimization problem, shown in Section 4.

Problem formulation

Under the hypothesis of a position-oriented control of the movement, we propose an optimization problem, based on the minimization of the jerk [START_REF] Hogan | On rhythmic and discrete movements: reflections, definitions and implications for motor control[END_REF] to find the optimal trajectories to substitute in Eq. ( 1) and Eq. ( 6), to obtain the corresponding driver's torque for different values of δ to substitute in Eq. ( 7) to obtain the amplification curves.

The considered movement consists into bringing the steering wheel from 0 rad to π/2 rad and returning it to its initial position. The problem is to find the optimal trajectory θ , that minimizes the cost index I:

I = T 0 d 3 θ dt 3 2 dt (8)
under the constraint to have a driver's torque τ v = 0 at the beginning and at the end of the movement.

The optimization problem can be written as follows:

min θ 1 0 d 3 θ dt 3 2 dt (9)
with the boundary conditions on the trajectory

θ (0) = 0 θ (1) = 0 θ (0.5) = π/2 θ (0) = 0 θ (1) = 0 θ (0.5) = 0 θ (0) = 0 θ (1) = 0 (10)
and on the driver's torque

τ v (0) = 0 τ v (1) = 0 (11)
Details of this procedure can be found in the paper [START_REF] Ciarla | Genesis of booster curves in electric power assistance steering systems[END_REF].

Simulation results

Fig. 2 shows the optimal trajectories obtained from the simulation. As expected from an optimization that minimizes the jerk, the optimal trajectories are smooths and regulars. Note that the angular trajectory is not symmetric in the last phase. This feature allows to obtain both the load and the driver's torques to be null at the final time. In Fig. 3, at steady-state the curve has an exponential behaviour. This characteristic is more accentuated, as the coefficient n diminishes. The amplification is almost negligible for driver's torques in the interval -1 ≤ τ v ≤ +1 N m, while it grows for more important values. The hysteresis effect is present when the optimization in transient state is carried out, due to the influence of the inertial and viscosity torques. Nevertheless, it is possible to conclude that the commonly used amplification curves correspond to the profile obtained with n < 1. This means that the Stevens' coefficient for the rate of growth of the stimulus corresponds to the hypothesis of the position-oriented control of the movement.

Figure 1 :

 1 Figure1: Example of the booster curve, provided by the EPAS system, reported in[START_REF] De | Fun-to-drive by feedback[END_REF].

Figure 2 :

 2 Figure 2: Evolution of the angular and speed profile used to evaluate the analytic solution of the contact friction model.

Figure 3 :

 3 Figure 3: Amplification torque obtained to compensate the friction torque with rate of the sensation stimulus n = 0.5 in steady-state case (dotted) and dynamical one (full line)

V. Ciarla, C. Canudas de Wit, F. Quaine, V. Cahouet-GIPSA-Lab -Report VOLHAND Task 5.2


Sylvain Prigent

Xavier Descombes

Didier Zugaj

Laurent Petit

Anne-Sophie Dugaret

Philippe Martel

Josiane Zerubia

Morphème / Ayin

† Morphème

Skin lesion evaluation from multispectral images

Keywords: skin, hyperpigmentation, change detection, statistical inferences, multispectral imaging peau, hyperpigmentation, détection de changements, inférences statistiques, images multispectral

During evaluation of skin disease treatments, dermatologists have to clinically measure the pathology severity of each patient during treatment periods. Such a process is sensitive to intra-and inter-dermatologist diagnosis. To make this severity measurement more robust we propose to use image processing to quantify the pathology severity. We focus on a hyperpigmentation disorder called melasma. During a treatment period, multispectral images are taken on patients receiving the same treatment. After co-registration and classication steps, we propose algorithms to measure the intensity, the size and the homogeneity evolution of the pathological area. Obtained results are compared with a dermatologist diagnosis using statistical tests on a full clinical study.

Introduction

In dermatology, specic protocols are established to quantify the evolution of lesions severity along a treatment period. In this report, we focus on a hyperpigmentation disorder called melasma. This disease shows darker and irregular spots on the face caused by an abnormal melanocytes activity in response to a hormonal reaction. Melasma has the specicity to be symmetric in the two cheeks. Thereby, for melasma, the protocol consists in applying a studied treatment in one cheek, and a comparator treatment in the second cheek [START_REF] Pandya | Guidelines for clinical trials in melasma[END_REF]. The comparator is often chosen as the vehicle used for the studied treatment without active product. To make a clinical study reliable, patients are randomized [START_REF] Pandya | Guidelines for clinical trials in melasma[END_REF]. The ecacy of a treatment is clinically measured by dermatologists. To this end, a clinical score is evaluated by one or several dermatologists in each cheek of each patient at several time moments during the treatment. Finally, a statistical test such as the Student test [START_REF] Mood | Introduction to the theory of statistics[END_REF] or the Wilcoxon test [START_REF] Mood | Introduction to the theory of statistics[END_REF] is performed on these measurements to obtain a quantication of the studied treatment ecacy relatively to the comparator treatment [START_REF] Pandya | Reliability assessment and validation of the melasma area and severity index (masi) and a new modied masi scoring method[END_REF]. The most widely used score for melasma severity evaluation is the Melasma Area and Severity Index (MASI) [4,[START_REF] Manaloto | Erbium:yag laser resurfacing for refractory melasma[END_REF][START_REF] Hurley | Ecacy of glycolic acid peels in the treatment of melasma[END_REF]. It has been introduced in [START_REF] Ck Kimbrough-Green | Topical retinoic acid (tretinoin) for melasma in black patients. a vehiclecontrolled clinical trial[END_REF]. The MASI split the face into 4 regions: forehead, left cheek, right cheek, and chin. For each region, three scores are evaluated: Area (A): Measures the area aected by the disease. It is estimated on a scale ranging from 0 to 6 as: 0 = 0%; 1 ≤ 10%; 2 = 10 -29%; 3 = 30 -49%; 4 = 50 -69%; 5 = 70 -89% and 6 = 90 -100%.

RR n°8136

Darkness (D): Color of the pathological area compared with normal skin. It is graded on a scale from 0 to 4 as follows: 0=normal skin color without evidence of hyperpigmentation; 1=barely visible hyperpigmentation; 2=mild hyperpigmentation; 3=moderate hyperpigmentation; 4=severe hyperpigmentation. Homogeneity (H): Homogeneity of the hyperpigmentation, graded on a scale from 0 to 4 as follows: 0=normal skin color without evidence of hyperpigmentation; 1=specks of involvement; 2=small patchy areas of involvement < 1.5 cm diameter; 3=patches of involvement > 2 cm diameter; 4=uniform skin involvement without any clear areas).

To integrate these three criteria measured on the four regions: forehead (F), left cheek (LC), right cheek (RC) and chin (C), the following equation is proposed in [START_REF] Ck Kimbrough-Green | Topical retinoic acid (tretinoin) for melasma in black patients. a vehiclecontrolled clinical trial[END_REF]:

M ASI = 0.3(D F + H F )A F + 0.3(D LC + H LC )A LC + 0.3(D RC + H RC )A RC + 0.1(D C + H C )A C , (1)
where D F , D LC , D RC , D C are darknesses measured respectively on the forehead, left cheek, right cheek and chin, H F , H LC , H RC , H C homogeneities measured respectively on the forehead, left cheek, right cheek and chin and A F , A LC , A RC , A C areas measured respectively on the forehead, left cheek, right cheek and chin. This MASI denition has been discussed in [START_REF] Pandya | Reliability assessment and validation of the melasma area and severity index (masi) and a new modied masi scoring method[END_REF]. The authors claim that the H criterion has signicant variations from one dermatologist to another.

They then propose remove the H criterion from Equation (1):

M ASI = 0.3D F A F + 0.3D LC A LC + 0.3D RC A RC + 0.1D C A C . (2)
In this report, we propose a electronic score called dierential score. The idea is to use images with an associated image processing algorithm in order to automatically evaluate a severity score containing area (A), darkness (D), and homogeneity (H) measurements. Such a computer based score avoids the variability on the measures obtained by human inspection. To describe the proposed protocol, we introduce multispectral imaging for skin analysis, present the clinical study used for experimental results, and describe the proposed scheme.

Multispectral imaging

Multispectral images are data cubes where each slice corresponds to a monochrome image taken at a specied wavelength. Such images allow the combination of spectral and spatial information of the imaged area. A multispectral image can be seen as a 2D spectrocolorimeter. Several algorithms have been proposed to analyse the spectrum of the skin, and showed that the skin spectra contain useful information. The most popular sprectral analysis algorithm for the skin is probably the one proposed by Stamatas et al. in [START_REF] Stamatas | Non-invasive measurements of skin pigmentation in situ[END_REF][START_REF] Stamatas | In vivo measurement of skin erythema and pigmentation: new means of implementation of diuse reectance spectroscopy with a commercial instrument[END_REF]. The authors model the skin absorbency as a mixture of the melanin, deoxy-haemoglobin and oxy-haemoglobin components. A linear regression associated with the Beer-Lamber law provides the proportions of the components mentioned above for each pixel of a multispectral image. The images we use are 1200*960 pixels, and contain 18 spectral bands with wavelengths from 405 nm to 970 nm.

Clinical study

To evaluate the proposed dierential score, we use a clinical study that contains two treatments denoted T and A. Each treatment is tested on a population of 22 patients. For each patient, one cheek receives one of the treatments T or A and the second cheek receives a comparator that is the vehicle used for treatments T and A. The study lasts 3 months and two multispectral images Inria are taken for each patient (one per cheek) every month. Thus, we have 2 * 22 * 2 = 88 series of four images corresponding to times t 0 = baseline, t 1 = after one month, t 2 = after two months and t 3 = after three months.

Proposed methodology

To evaluate the ecacy of a treatment, we need to measure the severity of the melasma on each image taken along the clinical study. To quantify the severity for an image, we rst classify the pathological area in the image, and then extract features describing the disease.

As we are interested in the pathology evolution, we measure change features from t 0 to t k according to the clinical protocol. Thus, the proposed scheme is drawn is Fig. 1. We rst classify each image taken at baseline (t 0 ). The classication step is not described in this report. In [10] we compare classication methods for skin hyperpigmentation. In [START_REF] Prigent | Classication of skin hyper-pigmentation lesions with multi-spectral images[END_REF] we proposed a method based on support vector machines with a global training and an interactive software to perform the classication on all the baseline images. In the second step, we extract a spectral feature describing the studied disease. We dene this feature as the linear combination of bands that maximizes the darkness of the disease. To compute the change detection, time series of images need to be co-registered. This is done with the algorithm described in [START_REF] Myronenko | Intensity-based image registration by minimizing residual complexity[END_REF] initialized by an ane registration. To evaluate the spatial changes of the disease, we propose two change detection algorithms between monochrome images obtained by the spectral feature. These two algorithms aim at quantifying criterion A and H, respectively. Finally, criteria calculated on the spectral feature map and the change map give an dierential score for each patient at each measurement time. To detail the presented steps, the report is organized as follows. Section 2 describes the spectral feature extraction method and the darkness criterion. Section 3 describes the proposed method to quantify the changed area between t 0 and t k . Section 4 proposes to use a multi-scale change detection analysis method based on statistical inferences to get the homogeneity criterion.

Section 5 proposes an dierential score from the criteria D, A and H computed in Section 2, 3 and 4. The last section concludes and discusses the proposed dierential score. To extract a spectral feature of the studied disease, let's assume that we have extracted a region of healthy skin and a region of pathological skin for each series of co-registered multispectral images. For a given patient, at a given time t k , we thus have four regions: I pS We propose to use these four regions to compute the weights of a spectral feature M dened by:

M (I) = N b b=1 m b I(b), (3)
where m b is the weight of the b st spectral band I(b) of the image I and N b the number of spectral bands in the initial multispectral image. We look for features that give the highest contrast between the healthy areas and the pathological areas.

RR n°8136

In the literature, feature extraction for spectral images has been widely studied. Methods of features extraction considering a linear combination of spectral bands are PCA [START_REF] Fauvel | Spectral and spatial methods for the classication of urban remote sensing data[END_REF], NAPCA [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF], SVD [START_REF] Ball | Hyperspectral pixel unmixing via spectral band selection and dc-insensitive singular value decomposition[END_REF], or ICA [START_REF] Cardoso | High-order contrasts for independent component analysis[END_REF][START_REF] Hyvarinen | Fast and robust xed-point algorithms for independent component analysis[END_REF] that tend to extract features by minimizing the correlation or higher order statistics. The MNF approach [START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[END_REF] searches features by factorizing the signal matrix into two non negative matrices. This factorization can be done by minimizing the mean square error or the Kullback-Leibler divergence for example [START_REF] Dhillon | Generalized nonnegative matrix approximations with bregman divergences[END_REF][START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. Algorithms such as N-FINDR [START_REF] Winter | Fast autonomous spectral end-member determination in hyperspectral data[END_REF][START_REF] Dobigeon | Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery[END_REF][START_REF] Chang | Random n-nder (n-ndr) endmember extraction algorithms for hyperspectral imagery[END_REF], PPI [START_REF] Boardman | Mapping target signatures via partial unmixing of aviris data[END_REF][START_REF] Plaza | A fast iterative algorithm for implementation of pixel purity index[END_REF] or VCA [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] analyse the geometry of the data set in the spectral N b dimensional space to extract features dened as extreme points in this set. We can also cite supervised methods like DAFE [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF], NWFE [START_REF] Lee | Feature extraction based on decision boundaries[END_REF][START_REF] Lee | Decision boundary feature extraction for neural networks[END_REF] or maximum likelihood based approach [START_REF] Riedmann | Supervised band selection for optimal use of data from airborne hyperspectral sensors[END_REF].

Proposed methods

For skin spectral feature extraction, we focus on two approaches: one involving an objective function optimization, and one based on independent component analysis (ICA). Before going to these two approaches we rst introduce a normalization step that combines the four regions

I pS t k , I hS t k , I pC t k , I hC t k
. The rst normalization combines the healthy and the pathological areas. We can either ignore the healthy area or compute the dierence between healthy and pathological areas.

We do not consider dividing healthy and pathological areas since division is more sensitive to noise. The second normalization combines the studied treatment S and the comparator treatment C. We use the method proposed in clinical protocols that calculate the dierence S -C. Finally, depending on the use of the healthy area we have two possible normalizations:

D e t k = (µ M (I pS t k ) ) -(µ M (I pC t k ) ) (4)
and

D e t k = (µ M (I hS t k ) -µ M (I pS t k ) ) -(µ M (I hC t k ) -µ M (I pC t k ) ), (5)
where D e t k denotes the darkness for the patient e at t k and µ M (I hS t k

) (respectively µ M (I pS t k ) , µ M (I hC t k ) , µ M (I pC t k
) ) is the spatial average of the spectral feature M on the region I hS t k (respectively

I pS t k , I hC t k , I pC t k
).

Objective function

In this section, we propose an objective function that maximizes both the distance between pathological and healthy areas, and the deviation between measurement times. The distance between healthy and pathological areas can be measured by a dierence to avoid noise introduced by divisions. Thus, the objective function is:

f = N b b=1 m b tn t=t1 Ne e=1 µ I e,hS t 0 -µ I e,pS t 0 -µ I e,hC t 0 + µ I e,pC t 0 -µ I e,hS t + µ I e,pS t +µ I e,hC t -µ I epC t . (6)
The function f is linear with respect to the parameters m b of M . The optimization is then trivial. The solution is the spectral band with the highest contrast if we consider only positive m b . If we consider also negative m b , the solution is given by two bands. These solutions are optimal with respect to the darkness criterion.

Inria

Figure 2: Feature obtained by ICA on all the images at t 0 on a clinical study of melasma. X-axis:

indexes of spectral bands, Y-axis: coecient given to spectral bands.

ICA-based method

In this section we propose to use ICA to get a spectral feature of melasma. Performing ICA on a multispectral image gives, at most, as many features as the number of bands in the image. We then select the most representative feature of the pathological area. The selected component is the one that maximizes the deviation between a healthy reference area and the pathological area of interest. If we calculate ICA separately on all the images we get slightly dierent features from one images to another. On Fig. 2 one can see that features obtained in all the images of a clinical study have similar shape. Several features are negative due to the ICA model indetermination.

To get a single feature for all the images, we propose to calculate the average feature with a -1 coecient for the negatives features. It can be noticed that the coecients in the near infra-red (IR) bands are negative. During experiments we observe that this near infra-red information biased the contrast measurements. The IR coecients for the ICA based feature are then set to zero.

Obtained results

Proposed feature extraction methods have been tested in a clinical study involving three active treatments that we denote S t , A d1 and A d2 . Clinical conclusions state that S t has no ecacy whereas A d1 and A d2 have ecacy. Furthermore, A d2 has a stronger eect than A d1 . The protocol of this study is the same as the one presented in section 1.2.

We calculate a spectral feature for each treatment S t , A d1 and A d2 individually. Then we compute the features in all the images of the study on the four regions

I pS t k , I hS t k , I pC t k , I hC t k manually segmented. It results in a darkness measurement D e t k
for each patient. We nally apply a Wilcoxon hypothesis test for each treatment between t k and t 0 to get a quantication of the treatment ecacy. Table 1 shows the p-values obtained by the Wilcoxon test with the proposed methods and with the luminance feature that we use as a reference feature. Equation ( 5) is used for normalization. As one can see, all the features agree with the clinical conclusions.

The dierence is in how fast features detect an ecacy of the treatment. The method based on the objective function optimization gives the fastest ecacy detection compared to the ICAbased method and the luminance. To understand these dierences, we plot the features obtained with all the methods in Fig. 3. As the luminance coecient is not linear, we use its linear approximation Y . We observe that the feature obtained by f selects a specic band at 590 nm RR n°8136 for the treatments S t , and A d2 and the band 645 nm for the treatment A d1 . ICA selects spectral bands mainly around 600 nm, but sets signicant coecients to wavelengths between 400 and 600 nm. The Y feature selects bands in the 500-600 nm wavelengths. Due to the Stamatas et al. model for pigmentation we should focus on the wavelength between 600 and 700 nm. Thus, heigh coecients in the spectral ranges out of 600-700 nm explain why the ICA based method and the Y features detect later the treatments ecacy than the feature estimated by f . Table 2 shows the results obtained with the Wilcoxon tests on the same studies with the normalization of Equation ( 4). The results are comparable to the ones with the normalization of Equation ( 5).

This means that the calibration of the camera is robust and that the normalization with the healthy area is not mandatory. As this healthy normalization can introduce noise due to the variation of the healthy area, we prefer to use the normalization of Equation (4). To validate the spectral features calculated on the test clinical study, we apply them on the validation study described in section 1. Inria

RR n°8136 3 Area criterion

To calculate an area criterion, we propose to compute the spatial changes between t k and t 0 . Lets denote C e t0 a classication mask at t 0 for the patient e obtained by the method proposed in [START_REF] Prigent | Classication of skin hyper-pigmentation lesions with multi-spectral images[END_REF]. We denote by B e t k/0 a binary change map between t 0 and t k for the patient e. Then, the area criterion can be calculated by:

A e t k/0 = card{x : B e t k/0 (x) = 1} card{x : C e t0 (x) = 1} . ( 7
)
In the literature, numerous methods have been proposed to quantify the changes between two images. Two state of the art articles are [START_REF] Lu | Change detection techniques[END_REF][START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF]. Change detection methods consider the image either globally by subtraction or division [START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF][START_REF] Skifstad | Illumination independent change detection for real world image sequences[END_REF][START_REF] Liu | Statistical change detection with moments under time-varying illumination[END_REF][START_REF] Li | Integrating intensity and texture dierences for robust change detection[END_REF] or locally with a window that browses the image [START_REF] Poor | An Introduction to Signal Detection and Estimation[END_REF][START_REF] Aach | Statistical model-based change detection in moving video[END_REF][START_REF] Rignot | Change detection techniques for ers-1 sar data[END_REF][START_REF] Bruzzone | An adaptive semiparametric and context-based approach to unsupervised change detection in multitemporal remote-sensing images[END_REF][START_REF] Inglada | A new statistical similarity measure for change detection in multitemporal sar images and its extension to multiscale change analysis[END_REF][START_REF] Yakimovsky | Boundary and object detection in real world images[END_REF].

Using a local window tends to a loss of resolution. Furthermore, the procedure of dividing images is sensitive to noise. That is why, we choose to perform change detection by subtraction:

F e t k/0 = M (I e t k ) -M (I e t0 ), (8)
where M (I

B e t k/0 (x) = 1 if F e t k/0 (x) > τ and C e t0 (x) = 1 0 otherwise (9)
where x denotes the position of any pixel in the image. The choice of τ is crucial. We propose a method that computes one threshold τ per image series estimated on the histogram of F e t 1/0 . First, we assume that the histogram of F e t 1/0 is monomodal symmetric and centred in zero in the case of no changes F e t k/0 . This assumption is justied since F e t 1/0 is a centred noise in zero in the case of no changes. Then, we assume that there are no signicant changes in the negative part (i.e. the disease does not become worst). Thus, in the case of no disease decreasing, the number of positive and negative pixels in F e t 1/0 should be the same. With this model, the area of changes in F e t 1/0 is A = ( >0)-( ≤0) ( >0)+( ≤0) , where ( > 0) (respectively ( ≤ 0)) denote the number of positive (respectively negative) pixels in F e t 1/0 . To calculate τ from F e t 1/0 , we sort the positive pixels of F e t 1/0 and take the threshold τ such as: ( > 0 and < τ ) = ( ≤ 0).

An inconvenience of thresholding methods for change detection is their sensitivity to noise. To decrease noise, we lter the change map B e t k/0 with opening and closing mathematical morphology operators using a disk of two pixels as a structural element. Fig. 9 d-f ) shows obtained results for an image time series. For this patient, dermatologists clinically measure a decrease of the pathology. As one can see, the area measured by the algorithm shows a signicant evolution at time t 1 and then a gradual evolution of the pathological area. To validate this area measurement, we calculate change areas A e t k/0 for all the patients of the validation clinical studies. As for the darkness measurement, we compute the Wilcoxon test for each treatment between t k and t 0 in the areas measurements. Results are shown in Table 4. Results show that the treatment T has an eect whereas the treatment A does not. This is in agreement with the dermatologists diagnosis.

Inria

Homogeneity criterion

With the proposed darkness (section 2) and area (section 3) criteria, we obtained a denition of a dierential score with two components D and A like in the clinical MASI dened in [START_REF] Pandya | Reliability assessment and validation of the melasma area and severity index (masi) and a new modied masi scoring method[END_REF]. In this section we study a dierential and numerical equivalent of the homogeneity clinical criterion proposed in [START_REF] Ck Kimbrough-Green | Topical retinoic acid (tretinoin) for melasma in black patients. a vehiclecontrolled clinical trial[END_REF]. This criterion jointly measures the distribution and the size of brown areas constituting the pathology area. To get such information, we adapt the methodology of Statistical Parametric Mapping (SPM) [START_REF] Friston | Statistical Parametric Mapping: The Analysis of Functional Brain Images[END_REF] developed for functional MRI. This approach models a change map as a Gaussian random eld. Then a multi-scale analysis partition this map into regions and a probability is assigned to each region. These probabilities measure the likelihood for the regions to be a realization of the random eld. For our application, we apply this methodology on each change map F e t k/0

. To model F e t k/0 as a Gaussian random eld, we need to normalize it.

Then we apply the multi-scale analysis to get regions and calculate probabilities.

Change map normalization

Change map thresholding

F N e t k/0
satises the conditions to apply the Gaussian eld theory [START_REF] Alder | The Geometry of Random Fields[END_REF]. An interesting characteristic of a Gaussian eld F of size S is the DT (Dierential Topology) characteristic dened by:

χ(A u ) = (-1) (d-1) d-1 k=0 (-1) k χ k (A u ) (10)
with

A u = {s ∈ S, F (s) ≥ u}, ( 11
)
where u is a given threshold, d the dimension of the eld and χ k (A u ) the number of points s ∈ S satisfying the conditions: (a)

F (s) = u, (b) F (i) (s) = 0, ∀i ∈ [1, d -1], (c) F (d) (s) > 0 and (d) the (d -1) × (d -1
) matrix of second order partial derivatives of F(s) has exactly k negative eigenvalues. In the above expression, F (i) denotes the i th derivative of F . In the 2D case (i.e. d = 2), for a given threshold u the expectation of the DT characteristic of a centred and reduced Gaussian eld can be written as follows [START_REF] Alder | The Geometry of Random Fields[END_REF]:

E χ(Au) = S(2π) 3 2 |Λ| 1 2 σ -3 ue -u 2 2σ 2 (12)

RR n°8136

with S the surface of the eld, Λ the 2 × 2 correlation matrix of the eld derivatives, and σ 2 = E(F 2 ).

In order to estimate the likelihood of a pixel of F N e t k/0 and its neighbourhood under the Gaussian eld assumption, we compare its characteristic at a threshold u with the DT characteristic. Two parameters can be used to compute this likelihood: the maximum intensity and the spatial extent of the considered neighbourhood above u. Let x 0 denotes the maximum intensity of R x0 u a region above the threshold u. The likelihood of this region under the Gaussian eld assumption is [START_REF] Poline | Combining spatial extent and peak intensity to test for activations in functional imaging[END_REF]:

P (R x0 u ∈ F ) = E χ(Ax 0 ) E χ(Au) = x 0 u e u 2 -x 2 0 2 . ( 13
)
Let S 0 denotes the spatial extent of R S0 u a region above the threshold u. The expectation for this region under the Gaussian eld assumption is [START_REF] Poline | Combining spatial extent and peak intensity to test for activations in functional imaging[END_REF]:

E(R S0 u ∈ F ) = E(N u ) E χ(Au) , (14)
where E(N u ) is the expectation to have N u pixels above the threshold u in a Gaussian random eld F . As F follows a standard normal distribution,

E(N u ) = S ∞ u (2π) 1 2 e -x 2 2 = SΦ(-u) (15)
then,

E(R S0 u ∈ F ) = Φ(-u) (2π) 3 2 |Λ| 1 2 ue -u 2 2 . (16)
Thus, the number of pixels in a region above the threshold u has an exponential distribution with the parameter λ u = 1/E(R S0 u ∈ F ). We can then write the probability of a region R S0 u to appear in a realization of the Gaussian eld by:

P (R S0 u ∈ F ) = e -λuS0 = e (2π) 3 2 |Λ| 1 2 S 0 ue -u 2 2 Φ(-u) (17)
Finally, for a spatial region above u with the maximum intensity x 0 and the spatial extend S 0 , its probability to be in F is [START_REF] Poline | Combining spatial extent and peak intensity to test for activations in functional imaging[END_REF]:

P (R x0,S0 u ∈ F ) = min(P (R x0 u ∈ F ), P (R S0 u ∈ F )). (18)
Thereby, to detect changes between t 0 and t k we dene a family of n u thresholds U = {u 1 , ..., u nu }. For each u i , F e t k/0 is segmented into regions denoted R xo,S0 ui . This segmentation is obtained by a connected components analysis of the class dened by the pixels of F e t k/0 above the threshold u i . Then we assign a probability P (R xo,S0 ui ∈ F ) to each region R xo,S0 ui . This process gives a statistical map of the changes denoted SM e t k between t 0 and t k . The algorithm 1 summarizes the procedure for obtaining SM t k between I 0 and I k for a given patient.

Final normalization

The statistical approach proposed in Section 4.2 makes the hypothesis that changes are rare events. In our images, a large part of the pathological area can change. Thus, if we apply directly the methodology describes above, we will detect only a small part of the changes (see 3: Choose a family U = u 1 , ..., u nu 4: for i = 1 to n u do 5: Compute the connected components above u i 6: for each connected component do 7:

Calculate P (R x0,S0 u ∈ F ) with eq. 18 8: end for 9: end for 10: Merge the obtained connected clusters obtained with the family U to get SM e t k gure 9 g-i)). Homogeneity variations cannot be quantied in these maps. We propose then to introduce the normalization:

F N d e t k/0 = F N e t k/0 + 1 N r t k r=t0 µ Ct 0 (F N e t r/0 ) -µ Ct 0 (F N e t 1/0 ) (19)
where N r is the number of time measurements between t 0 and t k . µ Ct 0 (F N e t r/0 ) (respectively µ Ct 0 (F N e t 1/0 )) denote the average of the pixels intensities in the interest area of F N e t r/0

(respectively

F N e t 1/0
) delimited by C t0 . If the normalized data with Equation ( 19) do not respect the rare event hypothesis, the calculated statistics will be biased. Nevertheless, we do not use the absolute values of these statistics but their spatial repartition. Fig. 9 To validate the proposed homogeneity criterion, we apply it on the validation study. The homogeneity measure makes sense only when some changes appear in the image between t 0 and . To qualitatively evaluate the homogeneity criterion, we plot its evolution for the population of patients whose pathology evolve. Due to the darkness and the area criteria, the considered patients are those who receive the treatment T in the validation study. Fig. 4 shows time plots of the homogeneity measurement for each patient of the group receiving the treatment T . If the treatment has an eect, we expect to see the homogeneity measurement growing. In this population of 22 patients, we observed that for only 7 patients the H criterion is not interpretable. Indeed, for ve patients less than 10% of the pathological area changes at time t 2 = 8 or t 3 = 12 and for two patients there were less than 10% of the pathological area that change in all the time sequence. For the others fteen patients we globally observe the expected growth of the homogeneity with time. For most patients (in blue) main changes occurred in the second measurement time t 1 = 4, and for 4 patients (in red), main changes occurred in time t 2 = 8. These results show that the homogeneity criterion gives relevant interpretation of disease progression during treatment. This criterion is interpretable only on patients whose pathology has evolved. This is why we recommend to use H e t k as an extra criterion to analyse the evolution of a patient pathology.

Graphical representation

The clinical MASI as described in [START_REF] Pandya | Reliability assessment and validation of the melasma area and severity index (masi) and a new modied masi scoring method[END_REF] and [START_REF] Ck Kimbrough-Green | Topical retinoic acid (tretinoin) for melasma in black patients. a vehiclecontrolled clinical trial[END_REF] calculates a linear combination of darkness, area and homogeneity measurements (see Equations ( 1) and ( 2)). For the dierential score, such linear combination does not have any meaning. Indeed, the three criteria have dierent dynamics, and combining them can bias the nal measurement. That is why we prefer to keep the three measurements independent. To make the dierential score easily readable for an individual patient we propose the following graphical representation using colored bars.

Darkness representation

The proposed darkness criterion is a dierential measure between two monochrome images. Its theoretical magnitude is then [-255, 255] for 8 bits images. In the experiments, we observe that the interval is [-12, 12]. To make a graphical representation of D, we use a function to transform Inria

D G = D √ 1 + D 2 (22)
where D G is the darkness value used for the graphical representation. Fig. 5 shows a plot of the function of Equation [START_REF] Dobigeon | Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery[END_REF]. The darkness bar contains a chip with a color from red (-12) to blue (+12). This chip position on the bar depends on the D G value from -12 to 12. Fig. 6 shows three examples of darkness representation bars.

Area representation

The area criterion we proposed has values in the interval [0, 1]. The graphical representation we choose is then a percentage bar. Fig. 7 shows examples of graphical representation of the area criterion.

(a) (b) (c)

Homogeneity representation

For the graphical representation of the homogeneity criterion, we propose to use f pe t k (s). red, that means there is no homogeneity variation. The slower the gradient changes from blue to red, the more heterogeneous the variation is. Fig. 8 shows examples of the homogeneity bar.

Dierential score representation

Finally, the dierential score representation is the combination of the three bars of darkness, area and homogeneity. An example for a patient whose pathology decrease is shown in Fig. 9 and an example for a patient whose pathology does not decrease is shown in Fig. 10. As one can see, with this graphical representation, we can both quantify and qualify the patient pathology evolution.

Conclusion

In this report, we proposed a strategy to calculate the severity of skin hyperpigmentation using

Figure 1 :

 1 Figure 1: Proposed scheme to evaluate a treatment ecacy on a clinical study

  regions on the cheek treated with the studied treatment, and I pC t k and I hC t k the pathological and healthy regions on the cheek treated with the comparator treatment.

Figure 3 :

 3 Figure 3: Spectral features obtained by the studied methods. X-axis: wavelengths in nanometers, Y-axis: coecients m b .

Table 1 : 1 A d1 4 .

 114 Results obtained by the Wilcoxon test on a clinical study with the normalization of Equation (5) for the darkness criterion. Feature Traitement t1 -t0 t2 -t0 t3 -t0 L St 7.959 10 -1 5.014 10 -1 7.173 10 -1 A d1 3.793 10 -1 1.128 10 -2 8.793 10 -2 A d2 9.798 10 -2 3.204 10 -3 5.312 10 -4 f St 1.476 10 -1 1.089 10 -1 7.563 10 -1 A d1 5.571 10 -2 3.400 10 -2 4.373 10 -2 A d2 8.360 10 -3 3.204 10 -3 9.350 10 -4 ICA St 5.349 10 -1 8.767 10 -1 9.176 10 -080 10 -1 2.987 10 -2 9.798 10 -2 A d2 7.873 10 -2 5.233 10 -3 9.350 10 -4

2 .

 2 The p-values of the Wilcoxon tests on the validation study are shown in Table 3. For f we choose the feature b 590 -b 405 . Obtained results are coherent with the clinical analysis.

First, we smooth the image F e t k/ 0

 0 with a Gaussian convolution lter to increase its spatial correlation. Second, as the interest area of F e t k/0 comes from a unique class, its histogram can be modelled as monomodal in the hypothesis of no changes. We use a histogram specication to transform it into a Gaussian distribution. We center and reduce the specied histogram to standardize it. In the following we will denote F N e t k/0 the normalized change map.

Inria Algorithm 1 2 :

 12 Multi-scale analysis Require: I 0 , I k , M p0 1: F e t k/0 = I k -I 0 Smooth and normalize F e t k/0 to get F N e t k/0 (see section 4.1)

Figure 4 :

 4 Figure 4: Time sequences of the homogeneity measurements for the patients of the group where an evolution has been measured by the area and darkness criteria. X-axis: time in weeks (t 0 = 1, t 1 = 4, t 2 = 8, t 3 = 12), Y-axis: homogeneity value.

Figure 5 :Figure 6 :

 56 Figure 5: Plot of the darkness transformation of Equation (22). X-axis: D in the scale values [-255, 255], Y-axis: D G values in the scale [0,1] .

Figure 7 :Figure 8 :

 78 Figure 7: Area graphical representation. a) for a patient with stationary pathological area, b) for a patient whose pathological area decrease by 26%, c) for a patient whose pathological area decrease by 86%.

  multispectral images. The proposed criterion called dierential score has been developed to be a computer score related to the clinical MASI criterion. The proposed dierential score has been tested on a full clinical study and gave relevant results with the clinical analysis. The proposed methodology does not take into account any specicity of the studied pathology. We can thus adapt it to other pathologies like vitiligo, scars or rosacea.

Figure 9 :Figure 10 :

 910 Figure 9: Changes map and graphical representation of the dierential score for a patient whose pathology decreases during the treatment period. a-c) equalized dierence image on the spectral feature map between each time t 1 , t 2 et t 3 and t 0 , d-f ) change detection map used to measure the area for each times t 1 , t 2 et t 3 , g-i) homogeneity maps SM e t kfor the three time t 1 , t 2 et t 3 without using the normalization of Equation (19), j-l) homogeneity maps SM e t k for the three time t 1 , t 2 et t 3 using the normalization of Equation (19), m-o) graphical representation of the dierential score for the three times t 1 , t 2 et t 3 .

Table 2 :

 2 Results obtained by the Wilcoxon test on a clinical study with the normalization of

	Equation (4) for the darkness criterion.
	Feature Traitement t1 -t0	t2 -t0	t3 -t0
		St	7.959 10 -1	6.791 10 -1 6.416 10 -1
	L	A d1	1.306 10 -2 2.289 10 -2 9.798 10 -2
		A d2	1.788 10 -1 7.169 10 -3 1.918 10 -3
		St	7.173 10 -1	5.349 10 -1 5.014 10 -1
	f	A d2	2.289 10 -2 1.737 10 -2 3.400 10 -2
		A d3	1.997 10 -2 1.609 10 -3 9.725 10 -3
		St	6.416 10 -1	5.694 10 -1 5.014 10 -1
	ICA	A d1	1.997 10 -2 3.400 10 -2 8.793 10 -2

A d2

1.476 10 -1 7.169 10 -3 1.123 10 -3

Table 3 :

 3 Results obtained by the Wilcoxon test on the validation clinical study for the darkness criterion.

	Test Treatment t1 -t0	t2 -t0	t3 -t0
	L	A T	8.455 10 -1 7.210 10 -1 1.270 10 -1 9.741 10 -1 9.397 10 -3 1.637 10 -3
	f	A T	1.116 10 -1 5.373 10 -1 9.777 10 -2 4.750 10 -1 3.859 10 -3 9.337 10 -3
	ICA	A T	7.210 10 -1 5.589 10 -1 1.627 10 -1 8.966 10 -1 6.389 10 -3 3.478 10 -3

Table 4 :

 4 p-values of the Wilcoxon tests on the validation clinical study for the area criterion.

	Criterion Treatment	t 1 -t 0	t 2 -t 0	t 3 -t 0
		A	2.840 10 -1 3.986 10 -1 1.045 10 -1
	Area	T	4.552 10 -1 3.478 10 -3 2.438 10 -4

  is an increasing function from [0, 1] to [0, 1]. As the interpretation of f e s (p) is not straightforward, we prefer to work with f pe t k (s) that is the inverse function of f e s (p). A linear interpolation is used to calculate f pe t k (s) with a regular sampling along s. f pe t k (s) is also an increasing function from [0, 1] to [0, 1]. In the following, we use f pe t k (s) for both calculating a graphical representation of the homogeneity criterion (see section 5), and a scalar H e

					j-l) shows SM e t k	cartographies
	obtained with the normalization of Equation (19). The pathology evolution is now visible during
	the time sequence.
	To quantify a treatment ecacy from homogeneity maps, we integrate them into a scalar
	criterion. To do so, we dene a function f e s (p) that represents the spatial repartition of the
	changes by:
		f e s (p) =	card{x : SM e t k (x) ≤ p} card{SM e t k }	.	(20)
	f e s (p) t k	describing the homogeneity of the
	patient e at time t k . To integrate f pe t k (s) in a scalar H e t k ∈ [0, 1], we can compute:
					1
		H e t k = 2	0	|f pe t k (s) -s|ds.	(21)
	H e t k	represents the area between the function f pe t k (s) and the identity function Id(s) = s. Id
	corresponds to the most heterogeneous repartition that can be expected. Then, more homoge-
	neous changes yield higher values H e t k	. To calculate H e t k	with discrete data, we use the Riemann
	integral.

  In fact this information is richer than H e t k . The graphical bar of homogeneity is then a plot of f pe t k (s) with a gradient color. Values of f pe t k (s) close to zero are in blue and values of f pe t k (s) close to one are in red. Intermediate values are represented with the jet colormap. If the gradient is all

Prigent & Descombes & others

RR n°8136

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

S. Prigent, X. Descombes and J. Zerubia would like to thank Galderma R&D for partial funding this research work and for providing the data.