C Poussot-Vassal

C Spelta

O Sename

S M Savaresi

L Dugard

Survey on Some Automotive Semi-Active Suspension Control Methods: a Comparative Study on a Single-Corner Model

Keywords: Semi-active suspension, Nonlinear control, LPV control, Performance evaluation, Benchmark

), here, a complete benchmark to evaluate the main semi-active suspension control design methods is proposed. The present paper gives a quick picture of the present state of the art in the semi-active suspension control field in terms of comfort and road-holding trade-off.

INTRODUCTION

In the automotive field, semi-active suspensions have received a lot of attention since they seem to provide the best compromise between cost (energy-consumption, actuators/sensors hardware) and performance. The research in this field follows two mainstreams: the study of new technologies of semi-active actuation of damping (like electro-hydraulic, electro-rheological and magnetorheological damper), and the design of semi-active control strategies (see e.g. [START_REF] Karnopp | Vibration Control Using Semi-Active Force Generators[END_REF]; [START_REF] Tseng | Semi-Active Control Laws -Optimal and Sub-Optimal[END_REF]; [START_REF] Hrovat | Survey of advanced suspension developments and related optimal control application[END_REF]; [START_REF] Kawabe | New semi-active suspension controller design using quasi-linearization and frequency shaping[END_REF]; Ahmadian and Reichert (2001); [START_REF] Guardabassi | Approximate Linearization via Feedback -an Overview[END_REF]; [START_REF] Hong | Modified skyhook control of semi-active suspensions: A new model, gain scheduling, and hardware-in-the-loop tuning[END_REF]; [START_REF] Giua | Design of a Predictive Semiactive Suspension System[END_REF]; [START_REF] Geurts | Extremely randomized trees[END_REF]; Giorgetti et al. (2006a); [START_REF] Canale | Semiactive suspension control using fast model-predictive techniques[END_REF]; [START_REF] Savaresi | A single-sensor control strategy for semi-active suspensions[END_REF] and references therein).

In a previous work, [START_REF] Poussot-Vassal | A Methodology for Optimal Semi-Active Suspension Systems Performance Evaluation[END_REF] proposed a methodology to evaluate the best performance trade-off, in term of comfort and road-holding, a semi-active suspension system can achieve. Based on this recent development, the main contribution of this paper is to provide a full analysis and benchmark of some of the most common and dedicated semi-active suspension control law performances, both for the comfort and for the road-holding objectives.

This work is carried out on the basis of the well known single-corner vehicle model including a semi-active suspension system (see Figure 1). Performances are then thoroughly evaluated trough dedicated frequency domain simulations.

The paper is organized as follows: Section 2 presents the semi-active suspension models and recalls the performance criteria used in [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF]. In Section 3 and 4, some classical dedicated semi-active control approach for comfort and road-holding performance enhancement are recalled. More generic or complex control algorithms are given in Section 5. All the proposed algorithms are evaluated with numerical simulations and compared to the optimal and passive performances (obtained in [START_REF] Poussot-Vassal | A Methodology for Optimal Semi-Active Suspension Systems Performance Evaluation[END_REF]), using both frequency responses and performance index, providing then a(n almost) complete semi-active suspension control benchmark. Conclusions are gathered in Section 7.

SINGLE-CORNER MODEL AND PERFORMANCE DEFINITIONS

Single-corner model

The single-corner car model is the basic model used for suspension analysis. It consists in double mass / spring / damper system describing the dynamics of the chassis and the center of the wheel. The suspension system is modeled by a spring and a damping system, and the tire is often reduced to a spring element. Figure 1 illustrates this singlecorner model, defined through relation (1).

         M z = -k(z -z t -∆ s ) -c( ż -żt ) -M g mz t = k(z -z t -∆ s ) + c( ż -żt ) -k t (z t -z r -∆ t ) -mg ċ = -β(c -c in ) z t -z r < ∆ t (1)
where z, z t , and z r are the vertical positions of the body, of the unsprung mass, and of the road profile respectively. M is the single-corner body mass; m is the unsprung mass (tire, wheel, brake calliper, suspension links, etc.). g is the gravitational constant. k ∈ R + and k t ∈ R + are the stiffness of the suspension spring and of the tire, respectively; ∆ s ∈ R + and ∆ t ∈ R + are the length of the unloaded suspension spring and tire, respectively. c ∈ R + and c in ∈ R + are the actual and the requested damping coefficients of the shock-absorber, respectively. The damping-coefficient variation is ruled by a 1st-order dynamic, where β ∈ R + is the bandwidth. The actual damping coefficient c always remains in that interval: c min ≤ c ≤ c max , where c min ∈ R + and c max ∈ R + are the shock-absorber technological limitations. This last inequality is the so-called "passivity-constraint" of a semiactive suspension -guaranteeing that the actuator only dissipates energy (see e.g. [START_REF] Savaresi | A single-sensor control strategy for semi-active suspensions[END_REF]; [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF]). Since the control signal c in modifies the damping coefficient c, a state variable, model (1) is nonlinear. When a passive suspension is considered, (1) is reduced to a 4th-order linear system (by simply setting ċ = 0 and c = constant -e.g. nominal damping).

For nonlinear control design, the following model is commonly used:

M z = -k(z -z t ) -c( ż -żt ) mz t = k(z -z t ) + c( ż -żt ) -k t (z t -z r ) ċ = β(c in -c) (2)
Similarly, in order to apply linear control techniques, a linear-like semi-active suspension model is (commonly) defined as follows:

   M z = -k(z -z t ) -c 0 ( ż -żt ) -F d mz t = k(z -z t ) + c 0 ( ż -żt ) -k t (z t -z r ) + F d Ḟd = -β(F d -u) (3)
With reference to models (1) and ( 2), the linearised model (3) includes symbols with the following meaning: c 0 denotes the nominal damping (which can be used as a design parameter when synthesizing the controller); F d is an additional damping force commanded by the control variable u, according to a 1st order actuation dynamic.

In order to be fully equivalent to (1) and (2), the control signal u should respect the passivity-constraint. In models (1) and (2) this constraint is described as c min ≤ c in ≤ c max . It is straightforward to see that in model (3) this constraint is recast as follows: the additional force reference u should satisfies u ∈ D(c min , c max , c 0 ) ⊆ R, where the dissipative D(c min , c max , c 0 ) set is defined as follows (see Figure 2):

∀ (U, V ) ∈ R × R | U -(c max -c 0 )V (c min -c 0 )V -U ≥ 0 (4)
where c min and c max are the minimal and maximal damping factors of the considered controlled damper, normalized around c 0 ∈ R + . Note that in formulation (3) of the single-corner model, the system with no control (u = 0) is damped thanks to c 0 (hence stable and damped). This remark is practically important for numerical reasons when some control approaches (such as LMI-based, MPC, . . . ) are applied.

Performance indexes definition

Since the paper aims at evaluating the trade-off of the semi-active suspension laws in terms of comfort and roadholding, the performance metric proposed in Poussot-Vassal et al. ( 2010) is used here. This criteria focuses on two specific signals, representing either the comfort or the road-holding performance (for further detail, refer to [START_REF] Hrovat | Survey of advanced suspension developments and related optimal control application[END_REF]; [START_REF] Kiencke | Automotive Control Systems[END_REF] and book of the authors [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF]), namely:

• The vertical chassis acceleration z (or displacement z) response to road disturbances z r , between 0 and 20Hz, representing the acceleration felt by the driver, i.e. the comfort specification. • The vertical wheel deflection z t -z r response to road disturbances z r , between 0 and 30Hz, representing the ability of the wheel to stay in contact with the road, i.e. the road-holding specification.

The common objective is then to minimize either the transfer z r to z -identically z -(comfort) or the transfer z r to z t -z r (road-holding) or a combination of these two transfer. More formally, let define the function C :

R × R × R → R, as C(x, f , f ) = f f |x(f )| 2 df
, where x(f ) represents the frequency dependent signal of interest, obtained by a nonlinear frequency response algorithm (similar to the one given in [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF]; [START_REF] Poussot-Vassal | A Methodology for Optimal Semi-Active Suspension Systems Performance Evaluation[END_REF]); f and f represent the interval limits of interest. Then, the comfort and road-holding criteria are respectively defined as:

• J z , Comfort criteria:

J z = C(F z , 0, 20) C(F nom z , 0, 20) (5)
• J z def t , Road-holding criteria:

J z def t = C(F z def t , 0, 30) C(F nom z def t , 0, 30) (6)
where F z and F z def t are the frequency responses of the controlled suspension; while F nom z and F nom z def t are the frequency responses of the passive uncontrolled reference suspension system with a nominal damping c = 1500N m/s (i.e. model ( 1) with c = 1500 and c in = 1500).

COMFORT-ORIENTED SEMI-ACTIVE SUSPENSION CONTROL

This section focusses on the main comfort-oriented dedicated semi-active suspension control strategies. This section gathers, in an increasing complexity order, the main semi-active control laws for comfort objective.

Skyhook two-state damper control (SH-2)

The two-state Skyhook control is an on/off strategy that switches between high and low damping coefficients in order to achieve body comfort specifications. This control law is recalled as follows.

Proposition 1. (SH-2 states control). Given system (2), the SH-2 law is defined as:

c in = c min if ż żdef ≤ 0 c max if ż żdef > 0 (7)
Basically, this control law consists in a switching controller which deactivates the controlled damper when the body speed and suspension deflection speed have opposite signs. The controlled damper technology only needs to have two damping coefficient states. This control strategy presents the advantage to be simple but requires two sensors.

Many studies have concerned the Skyhook control strategy since it represents a simple but efficient way to achieve good comfort requirement (see e.g. [START_REF] Simon | An Investigation of the Effectiveness of Skyhook Suspensions for Controlling Roll Dynamics of Sport Utility Vehicles Using Magneto-Rheological Dampers[END_REF]; Ahmadian et al. ( 2004)). Some extended versions of the Skyhook control have been also developed, such as the adaptive one in [START_REF] Song | Analysis and strategy for superharmonics with semiactive suspension control systems[END_REF] or the gain-scheduled one in [START_REF] Hong | Modified skyhook control of semi-active suspensions: A new model, gain scheduling, and hardware-in-the-loop tuning[END_REF].

Skyhook linear approximation damper control (SH-L)

An improved version of Skyhook control has been used to handle variable damping, either with discrete damping coefficients, or with continuously variable damper, as illustrated in [START_REF] Sohn | Semi-active control of the macpherson suspension system: Harwarein-the-loop simulations[END_REF]; [START_REF] Sammier | Skyhook and H ∞ control of active vehicle suspensions: some practical aspects[END_REF]. The linear approximation of the Skyhook control algorithm, adapted to semi-active suspension actuators, is given as:

Proposition 2. (SH-L). Given system (2), the SH-L law is defined as:

c in =    c min if ż żdef ≤ 0 sat αc max żdef + (1 -α)c max ż żdef if ż żdef > 0
(8) where α ∈ [0; 1] is a tuning parameter that modifies the closed-loop performances and sat denotes that c in ∈ [c min ; c max ].

When α = 1, this control law is equivalent to the Skyhook two-state one. As the two-state control, the linear approximation consists in a switching controller which modifies the damping factor according to the body speed and suspension deflection speed. The innovation rely in the fact that, according to the second expression (when ż żdef > 0), such a control provides an infinite number of damping coefficients. As a matter of fact, this control law requires a continuously variable controlled damper (e.g. an MR dampers). From the computational point of view, this control law also requires two measurements and is simple to implement, but suffers of żdef zero crossing as well (which is practically complex to measure).

Acceleration Driven Damper control (ADD)

The ADD control is a semi-active control law described in [START_REF] Savaresi | Acceleration driven damper (ADD): an optimal control algorithm for comfort oriented semi-active suspensions[END_REF], which consists in changing the damping factor using the acceleration knowledge. Proposition 3. (ADD). Given system (2), the ADD law is defined as:

c in = c min if z żdef ≤ 0 c max if z żdef > 0 (9)
This strategy shows to be optimal in the sense that it minimizes the vertical body acceleration when no road information is available. Since it requires the same number of sensors as the Skyhook two-state and the linear approximations control law, this control law is simple from the implementation point of view. Note that the control law is very similar to the two-state approximation of the Skyhook algorithm, with the difference that the switching law depends on the body acceleration (z), instead of the body speed (which is easier to practically measure). It is also to note that the ADD design is well adapted to comfort improvement but not to road-holding. Nevertheless, as the previous approach, the "switching dynamic" may influences the closed-loop performances, implying chattering effects.

Power Driven Damper (PDD)

In [START_REF] Morselli | Control of a port hamiltonian systems by dissipative devices and its application to improve the semi-active suspension behavior[END_REF], the authors propose a semiactive suspension control strategy using the port Hamiltonian techniques, which provide powerful tools for modeling mechatronics systems with dissipative components. Based on this observation, it is straightforward that this framework perfectly fits the semi-active suspension problem. Proposition 4. (PDD). Given system (2), the PDD control approach is described by:

c in =            c min if kz def żdef + c min żdef ≥ 0 c max if kz def żdef + c max żdef < 0 c min + c max 2 if z def = 0 and żdef = 0 - kz def żdef otherwise (10)
where k is the stiffness of the considered suspension.

In [START_REF] Morselli | Control of a port hamiltonian systems by dissipative devices and its application to improve the semi-active suspension behavior[END_REF] (see also Figures 3,6 and 7), the authors show that this strategy provides results comparable to those of the ADD control law, while avoiding the chattering effect of the damping control value. The additional cost is the need for the knowledge of the spring stiffness k and a more complex rule.

Mixed Skyhook-Acceleration Driven Damper (SH-ADD)

The Mixed Skyhook-ADD rational mixes the the best behaviour of SH and ADD, without an increasing of either the computational effort or the hardware complexity. The key idea exploits a very simple but effective frequency range selector, which is able to distinguish the instantaneous dynamical behaviour of the suspension: in the case of low frequency dynamics the SH is selected, while the ADD is selected otherwise. The resulting control law is incredibly simple and requires the same apparatus as SH.

Proposition 5. (SH-ADD). Given system (2), the mixed SH-ADD control approach is described by:

c in =    c max if (z 2 -α 2 ż2 ) ≤ 0 AND ż żdef > 0 OR (z 2 -α 2 ż2 ) > 0 AND ż żdef > 0 c min otherwise (11)
where α ∈ R + is the tuning parameter allowing for frequency range selector, i.e. adjusts the "switch" between the SH and the ADD.

The amount (z 2 -α 2 ż2 ) hence can be considered as a simple "frequency-range selector". The parameter α represents the frequency limit between the low and the high frequency ranges, and it is the only tuning knob of the control strategy (11). Specifically the value of α is set at the cross-over frequency (in rad/s) between SH and ADD. For a standard motorcycle suspension it has to be selected around 19rad/s (3Hz). A simplified version of this algorithm, employing one single sensor has also been developed, leading to very satisfactory results both in simulation and experimental benchmark (see [START_REF] Savaresi | A single-sensor control strategy for semi-active suspensions[END_REF]).

ROAD-HOLDING ORIENTED SEMI-ACTIVE SUSPENSION CONTROL

Complementary to comfort-oriented control strategies, a very few studies have been devoted to the possible improvement of road-holding, using suspension actuator. Since few years, the studies on Global Chassis Control have emphasized that the suspension system may also help getting better road-holding and even handling, encouraging the community to reach this new objective.

Ground-hook 2 states (GH-2)

In a dual way to the Skyhook case, the 2-states Groundhook control (see [START_REF] Valasek | Development of semi-active roadfriendly truck suspensions[END_REF]) consists in a switching control law depending now on the sign of the product between the suspension deflection velocity żdef and the velocity of the unsprung mass żt , as: Proposition 6. (GH-2). Given system (2), the GH-2 control approach is:

c in = c min if -żt żdef ≤ 0 c max if -żt żdef > 0 (12)
This control has globally the same properties as the SH-2 one, but focussing around the unsprung mass.

Ground-hook linear (GH-L)

In this case, the semi-active damper allows to continuously change the damping coefficient, according to: Proposition 7. (GH-L). Given system (2), the GH-L control approach is defined through:

c in =    c min if -żt żdef ≤ 0 sat αc max żdef + (1 -α)c max żt żdef if -żt żdef > 0
(13) where α ∈ [0; 1] is a tuning parameter that modifies the closed-loop performances and sat denotes that c in ∈ [c min ; c max ].

ADVANCED SEMI-ACTIVE SUSPENSION CONTROL

Until now, the presented control approaches where completely dedicated to semi-active suspension systems and based on the nonlinear model (2). In this section we describe some of the main semi-active suspension control method, based on more classical (or usual) control tools, and relying on model (3).

Clipped approaches

Many works have concerned the application of classical control methods (e.g. H ∞ , H 2 , pole placement, disturbance rejection, optimal, active Skyhook . . . ). However, most of the results were obtained for active suspensions, as in [START_REF] Zin | Robust LPV -H ∞ control for active suspensions with performance adaptation in view of global chassis control[END_REF]. When applied to semi-active dampers, the dissipative constraint of the damper is usually handled using a simple projection (i.e. saturation, as shown in Figure 2). Even if it is not always referred to as the "clipped approach", the latter is very widespread in control strategies for semi-active suspension (see e.g. [START_REF] Karnopp | Active damping in road vehicle suspension systems[END_REF] and [START_REF] Margolis | Semi-active control of wheel hop in ground vehicles[END_REF]). In the control step, the force applied by the semi-active damper is then chosen to be as close as to the force required by the controller for a given suspension deflection speed and for the possible range of forces the damper can deliver. This simple strategy has been then applied in many cases (see e.g. [START_REF] Rossi | H ∞ control of automotive semi-active suspensions[END_REF]; [START_REF] Du | Semi-active H ∞ control with magneto-rheological dampers[END_REF]; [START_REF] Sename | Robust H ∞ control of quarter-car semi-active suspensions[END_REF]).

The question that arises is: is optimal the clipped-optimal? If not, how far is it from the real optimal one? How would look like the optimal semi-active one? Clipped approaches lead to unpredictable behaviours and ensure neither closed-loop internal stability nor performances any longer. As a matter of fact, active control applied on a semi-active damper results in a "synthesize and try" method.

To cope with this last drawback, some modern control techniques have been applied to the specific semi-active suspension problem.

Hybrid MPC control approaches

In Giorgetti et al. (2006a), the authors introduce an hybrid model predictive optimal controller (using receding horizon). They solve an off-line optimization process which is a finite horizon optimal regulation problem s.t.:

J * (ξ, x(k)) = min ξ x T (N )Q N x(N )+ N -1 k=1 x T (k)Qx(k)+y 2 (k) (14) subject to,      x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k) 0 ≤ u(k) żdef (k) |u(k)| ≤ Λ (15)
where Q is a performance index and Q N is the final weight, as in the optimal control theory. Matrices A, B, C and D in (15) define the LTI single-corner model (3), Λ is the maximal force allowed by the considered controlled damper and u(k) żdef (k) ≥ 0 guarantees the passivity constraint. ξ is a vector composed by the sequence of control signals (from 0 to N -1) to be applied, where, N is the prediction horizon. Giorgetti et al. (2006a) show that choosing N = 1 leads to performances that are identical to those of the clipped-optimal approach, and by increasing N (e.g. until 40), the performances can be significantly improved. The implemented control law does not involve any optimization procedure since the control algorithm provides a collection of affine gains over a polyhedral partition of the system states x (e.g. [START_REF] Borrelli | An efficient algorithm for computing the state feedback optimal control law for discrete time hybrid systems[END_REF]). By the way, this approach exhibits notable drawbacks, such as high complexity and switching between control regions and requires a complete full-state measurement.

LPV semi-active control

In an other philosophy, the "LPV semi-active" control adjustment, introduced in Poussot-Vassal et al. (2008), is a robust semi-active suspension control design using a Linear Parameter Varying (LPV) approach. More specifically, the proposed semi-active suspension control strategy is designed so that it minimizes the H ∞ performance criteria while guaranteeing the dissipative constraint thanks to a specific parameter dependent structure and a scheduling strategy design. The "LPV semi-active" controller design is summarized in the following algorithm (for H ∞ performances).

In this algorithm, ρ is a parameter function of the type of semi-active actuator. It varies as a function of the deflection velocity and the required damping force (for more details refer to [START_REF] Poussot-Vassal | A New Semi-active Suspension Control Strategy Through LPV Technique[END_REF]; [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF]). The main interest of such approach is that it presents a large flexibility concerning the performances and the type and number of involved sensors. Additionally, it provides a robustness certificate thanks to the polytopic Algorithm 1 H ∞ "LPV semi-active" (1) Based on initial system (3), and ρ-dependent performance weighting functions, construct the LPV polytopic system. Note that to complete this step, weighting filters should be defined and parametrized (refer to [START_REF] Poussot-Vassal | A New Semi-active Suspension Control Strategy Through LPV Technique[END_REF] and Chapter 8 of [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF], for details). ( 2) Compute the H ∞ LMI based "LPV semi-active" feasibility conditions. (3) Reconstruct the LPV controller and obtain two dynamical full order controllers are obtained which are solution of the robust H ∞ "LPV semi-active" control problem.

C(ρ) = A c (ρ) B c (ρ) C c (ρ) 0 and C(ρ) = A c (ρ) B c (ρ) C c (ρ) 0 (16)
(4) Apply control law (where ρ is varying):

u = |ρ -ρ| ρ -ρ C(ρ) + |ρ -ρ| ρ -ρ C(ρ) z def (17)
approach. The inherent counterpart of this flexibility is an important complexity in the design step.

NUMERICAL SIMULATIONS

In order to benchmark the performances of all these algorithms, the criteria presented in [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF] and recalled in Section 2 is applied. All simulation are carried out using a motorcycle parameter set (available in [START_REF] Poussot-Vassal | A Methodology for Optimal Semi-Active Suspension Systems Performance Evaluation[END_REF] and [START_REF] Savaresi | Semi-Active Suspension Control Design for Vehicles[END_REF]).

Note that in the following we do not simulate neither the Clipped approach (because of the nearly infinite way of adjusting) nor the MPC approach (complex to adjust due to multiple parameters, as well), but interested reader should refer to [START_REF] Tseng | Semi-Active Control Laws -Optimal and Sub-Optimal[END_REF]; [START_REF] Giorgetti | Hybrid model predictive control application toward optimal semi-active suspension[END_REF].

Frequency domain analysis

First, we compare the frequency responses of F z , the transfer between z r and z (comfort), and F z def t , the transfer from z r to z def t = z t -z r . On Figure 3, the SH 2-states, SH linear (with α = 0), ADD and PDD are compared with the passive ones with either c min or c max , illustrating to attenuate well the F z transfer while degrading the F z def t one.

Similarly, Figure 4 compares the performance attenuation of the two road-holding control laws, namely, GH 2-states and GH linear (with α = 0) with the passive ones with either c min or c max .

It clearly emphasizes the improvement on F z def t while degrading the F z one. Moreover, it confirms the observation made in [START_REF] Morselli | Control of a port hamiltonian systems by dissipative devices and its application to improve the semi-active suspension behavior[END_REF] concerning the performances of the ADD and the PDD approaches.

Finally, the Mixed SH-ADD (whith α adjusted to 20Hz) and "LPV semi-active" (with either a comfort or roadholding adjustment) frequency responses are compared to the passive cases on Figure 5. the two configurations tested shows the large flexibility of the approach. In the next subsection, all these approaches are gathered and compared to the passive and optimal performances, using the index defined in Section 2.

On this last

Performance indexes & trade-off analysis

Since improvement in one criteria implies a deterioration on the other, the performance index introduced in Section 2 is also evaluated for each control strategy. On Figure 6, all the performance indexes are compared to the passive with three configurations: c = {c min , 1500, c max } (note that c = 1500 has been selected as the nominal reference value since it provides a nice compromise between comfort and road-holding).

In order to be even more precise, on Figure 7, the trade-off between comfort and road-holding is illustrated for all the previous strategies using a diagram with Comfort in the x-axis and Road-holding in the y-axis. This figure also includes the criteria evaluation for a passive suspension with damping varying from c = 100N/m/s to c = 10000N/m/s and the optimal bound numerically computed using the optimization framework presented in [START_REF] Poussot-Vassal | A Methodology for Optimal Semi-Active Suspension Systems Performance Evaluation[END_REF].

The interesting point of this diagram is that (i) it illustrates the interest of the control to enhance the passive performances, (ii) shows the optimal performances of the SH-ADD approach in term of comfort (blue cross close to the optimal bound, in dashed red) and (iii) emphasizes the interest of the flexible LPV semi-active rule, allowing to achieve good comfort performances (red left-oriented triangle) or road-holding ones (blue right-oriented triangle).

CONCLUSIONS AND COMPLEXITY REMARKS

In this paper we provide a benchmark of the main semiactive suspension control strategies. The proposed benchmark is of course not exhaustive but still quite large. We illustrate the interest of the recent dedicated semi-active suspension approach (such as SH-ADD), presenting a low complexity, while providing great comfort performances. Nevertheless, other approaches (such as LPV semi-active, or Hybrid MPC) shows some interests in the fact that they are quite flexible, but far more complex to adjust and implement. The proposed benchmark also include a complete evaluation of the passive suspension and the optimal performances such a system can achieve, and thus extends the results obtained in [START_REF] Poussot-Vassal | A Methodology for Optimal Semi-Active Suspension Systems Performance Evaluation[END_REF].

REFERENCES Ahmadian, M. and Reichert, B. (2001). System Nonlinearities Induced by Skyhook Dampers. Shock and Vibration, 8(2), 95-104. Ahmadian, M., Song, X., and Southward, S. (2004). Nojerk skyhook control methods for semiactive suspensions. Transactions of the ASME, 126, 580-584.

Fig. 1 .

 1 Fig. 1. Single-corner vehicle model.

Fig. 2 .

 2 Fig. 2. Illustration of the dissipative D(c min , c max , c 0 ) set as a function of c 0 . Left: c 0 = 0, right: c 0 = cmin+cmax 2

Fig. 3 .

 3 Fig. 3. Comfort oriented control law frequency response F z (top) and F z def t (bottom).

Fig. 4 .

 4 Fig. 4. Road-holding oriented control law frequency response F z (top) and F z def t (bottom).

Fig. 5 .

 5 Fig. 5. Road-holding oriented control law frequency response F z (top) and F z def t (bottom).

Fig. 6 .

 6 Fig. 6. Normalized performance criteria comparison for different comfort (resp. road-holding) oriented control strategies. Comfort criteria -when cost function is J c (left histogram set) and Road-Holding criteria -when cost function is J rh (right histogram set).

Fig. 7 .

 7 Fig. 7. Normalized performance criteria trade-off for the presented control algorithms, compared to the passive suspension system, with damping value c ∈ [c min ; c max ] (solid line with varying color), optimal comfort and road-holding bounds (red dash dotted).


Bertrand Rivet

Mohammad Niknazar

Christian Jutten

Nonparametric modelling of ECG: Applications to denoising and to single sensor fetal ECG extraction

Keywords: non-parametric modelling, source extraction, denoising, fetal ECG extraction

In this work, we tackle the problem of fetal electrocardiogram (ECG) extraction from a single sensor. The proposed method is based on non-parametric modelling of the ECG signal described thanks to its second order statistics. Each assumed source in the mixture is thus modelled as a second order process thanks to its covariance function. This modelling allows to reconstruct each source by maximizing the related posterior distribution. The proposed method is tested on synthetic data to evaluate its performance behavior to denoise ECG. It is then applied on real data to extract fetal ECG from a single maternal abdominal sensor.

Introduction

Fetal electrocardiogram (f-ECG) extraction from maternal abdominal ECG sensors is an old problem. Since the first works of Cremer [START_REF] Cremer | Über die direkte Ableitung der Aktionsströme des menschlichen Herzens vom Ösophagus und über das Elektrokardiogramm des fötus[END_REF] who produced a very primitive record of fetal rate activity, this problem is still of interest nowadays since it is a fascinating issue due to the characteristics of the involved signals. Indeed, the f-ECG is definitively less powerful than the mother's ECG (m-ECG), moreover the recorded signals are also contaminated by noise due to electromyogram (EMG) or to power line interference and they are also influenced by the fluctuation of the baseline. Among the several approaches used to tackle the extraction of f-ECG, one can quote methods which require several sensors e.g., adaptive filtering [START_REF] Widrow | Adaptive noise cancelling: Principles and applications[END_REF], blind source separation [START_REF] Lathauwer | Fetal electrocardiogram extraction by source subspace separation[END_REF][START_REF] Zarzoso | Noninvasive fetal electrocardiogram extraction : blind source separation versus adaptative noise cancellation[END_REF] or quasi-periodic analysis [START_REF] Tsalaile | Sequential blind source extraction for quasi-periodic signals with time-varying period[END_REF].

In this paper, we consider the same issue but assuming that only a single sensor is available. In this case, one can extract f-ECG by singular value decomposition [START_REF] Kanjilal | Fetal ecg extraction from single-channel maternal ecg using singular value decomposition[END_REF] or by nonlinear decomposition such as shrinkage wavelet denoising [START_REF] Mochimaru | Detecting the fetal electrocardiogram by wavelet theory-based methods[END_REF] or nonlinear projections [START_REF] Richter | Fetal ecg extraction with nonlinear state-space projections[END_REF]. Moreover, state modelling as Kalman filtering [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF] has been applied to overcome the lack of information provided by a single sensor. Among these methods, the latter has been shown to be the most efficient [START_REF] Sameni | A review of fetal ECG signal processing; issues and promising directions[END_REF]. However, Kalman filtering relies on a very strong assumption: the state equation, which models the dynamical evolution of the unobserved state. As a consequence, Kalman filtering needs reliable prior about the state to perform accurately. To overcome the potential lack of prior information about the system, we propose in this study to model the second order statistics of the signal instead of the signal itself. This article is organized as follows. Section 2 presents the proposed approach to model a signal thanks to its second order statistics. The proposed algorithm to extract f-ECG is then introduced in Section 3. Numerical experiments and results are given in Section 4 before conclusion and perspectives in Section 5.

Nonparametric modelling of ECG

As already proposed in [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF], one can choose a parametric model of ECG: each beat of an ECG signal is a summation of 5 Gaussian functions, each of them modelling the P, Q, R, S and T waves as illustrated in Figure 1:

z(θ) = i∈{P,Q,R,S,T } a i exp - (θ -θ i ) 2 2σ 2 i , (1)
where a i , θ i and σ i are the amplitude, the position and the width of each wave, respectively. Note that in this model, the beats are defined in phase θ ∈ [-π, π], so that each beat is assumed to have a linear variation of phase with respect to the time, even if each beat has not the same duration. This model can then be used in an extended Kalman filtering to denoise a single ECG or extract f-ECG from a mixture of m-ECG and f-ECG [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF][START_REF] Sameni | Extraction of Fetal Cardiac Signals from an Array of Maternal Abdominal Recordings[END_REF]. This method is thus a parametric method since the unknown amplitude z(θ) is explicitly parameterized. On the other hand, nonparametric methods perform estimation, prediction or denoising without explicitly parameterizing the unknown amplitude z(θ). For instance, a well known approach is the spline smoothing [START_REF] Kimeldorf | A correspondence between bayesian estimation on stochastic processes and smoothing by splines[END_REF]. If one considers the ECG z(θ) as a statistical process, it can be fully described at the second order by Fig. 2. Two functions drawn at random from a zero-mean GP with covariance function [START_REF] Lathauwer | Fetal electrocardiogram extraction by source subspace separation[END_REF]. The shaded area represents plus and minus two times the standard deviation for the prior. On the right, the related σ(θ) and l d (θ) functions.

its mean function m(θ) = E[z(θ)] and covariance function k(θ 1 , θ 2 ) △ = E[(z(θ 1 ) - m(θ 1 ))(z(θ 2 ) -m(θ 2 ))] [8]
. Obviously, the ECG signal is almost surely a more complex statistical process than a simple second order one. As a consequence, considering only its second order statistics, it relies among the Gaussian process (GP) framework which is widely used in machine learning e.g., [START_REF] Mackay | Introduction to Gaussian processes[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. A GP z(θ) is a distribution over functions denoted as GP(m(θ), k(θ 1 , θ 2 )). In this case, the statistical latent process z(θ) is not directly parameterized as in parametric model, but its statistics are it thanks to hyper-parameters. This means that one has to choose a class of semidefinite positive functions k(θ 1 , θ 2 ) which describes the expected second order properties of the latent process.

In this study, we propose to use the following non-stationary covariance function

k(θ 1 , θ 2 ) = σ(θ 1 )σ(θ 2 ) 2l d (θ 1 )l d (θ 2 ) l d (θ 1 ) 2 + l d (θ 2 ) 2 exp - θ 1 -θ 2 2 l d (θ 1 ) 2 + l d (θ 2 ) 2 , (2) with σ(θ) = a m + (a M -a m ) exp - (θ -θ 0 ) 2 2σ 2 T , l d (θ) = l M -(l M -l m ) exp - (θ -θ 0 ) 2 2σ 2 l ,
where σ(θ) and l d (θ) allow to have a time-varying power (between a m and a M ) and a time-varying length scale correlation (between l m and l M ), respectively. Indeed, as shown in Fig. 1, an ECG beat can be decomposed into three parts: the P wave, the QRS complex and the T wave. The P and T waves share the same kind of second order statistics: a larger length scale and a lower power than the QRS complex. Fig. 2 shows two functions drawn at random from the zero-mean GP prior with covariance function [START_REF] Lathauwer | Fetal electrocardiogram extraction by source subspace separation[END_REF]. This figure illustrates the flexibility of such a representation compared to model (1) since with the same prior GP 0, k(θ 1 , θ 2 ) , it can generate a multitude of different shapes with the same prior.

3 Denoising of ECG and extraction of fetal ECG from a single sensor Suppose that the observed values x(t) differ from the ECG, s(t), by additive noise n(t):

x(t) = s(t) + n(t), (3)
and that this noise is uncorrelated with s(t). The aim of this study is to infer the values of s(t) from x(t), i.e. to denoise or extract the ECG from the observations. Moreover, it is assumed that the ECG signal, s(t), is a succession of beats, each of them following a zero-mean GP defining by [START_REF] Lathauwer | Fetal electrocardiogram extraction by source subspace separation[END_REF] and that the additive noise follows a zero-mean GP whose covariance function k n (t, t ′ ) is given by

k n (t, t ′ ) = σ 2 n exp - (t -t ′ ) 2 2 l 2 n + σ 2 w δ(t -t ′ ), (4)
where δ(•) is the delta Dirac function. In this expression, the first term is useful for instance to model a baseline variation as a stationary process for which the correlation is almost unity between close samples and decreases as their distance increases compared to the length scale l n . The second term models a white Gaussian noise of power σ 2 w . From ( 3) and ( 4), the covariance function of observation x(t) is thus expressed as

k x t, t ′ = k s t, t ′ + k n t, t ′ , (5)
where

k s t, t ′ = N n=1 N n ′ =1 k t -τ n , t ′ -τ n ′
and {τ n } 1≤n≤N is the set of R peak instants that can be estimated easily from the raw signals. From this modelling and assuming that the observed process x(t) has been recorded at times {T m } 1≤m≤M , the covariance matrix of this process is thus given by K x , whose (i, j)th entry is

K x (i, j) = k x (T i , T j ). (6)
One can then infer on the value s(t) thanks to the maximization of the a posteriori distribution of s(t)

given x = [x(T 1 ), • • • , x(T M )] T by ŝ(t) = k T K -1 x x, (7)
where

k = [k(t, T 1 ), • • • , k(t, T M )] T .
It is interesting to note that as soon as σ 2 w = 0, matrix K x is invertible as the summation of definite positive matrices and a diagonal matrix σ 2 w I. This algorithm needs some comments. First of all, the recorded signal, x(t), does not need to be regularly sampled and one can observe from ( 7) that the value of the latent process, s(t), can be predicted at any time t even for

t = T i , ∀i ∈ {1, • • • , T M }. Moreover, the hyper-parameters θ = {a m , a M , σ 2 T , l m , l M , σ 2 l , θ 0 , {τ k } k , σ 2 n , l 2 n , σ 2 w } defining k(•, •) and k n (•,
•) need to be estimated. This can be done by maximizing the evidence (or log marginal likelihood) given by log

p x|{T i } i , θ = - 1 2 x T K s + K n -1 x - 1 2 log K s + K n - M 2 log(2π). ( 8
)
The optimization of the latter equation is obtained thanks to a gradient ascent method, assuming that the initial parameter values are not so far from its actual values.

Fetal ECG extraction from a single abdominal sensor is then a direct extension of the proposed method by modelling the recorded signal x(t) as

x(t) = s m (t) + s f (t) + n(t), (9)
where s m (t) is the signal related to the mother, s f (t) is related to the fetus and n(t) is the additive noise. All these signals are modelled by zero-mean GPs with covariance functions k m (•, •) and k f (•, •) defined by ( 2) and k n (•, •) obtained from (4), respectively. In this case, the estimation of s f (t) is given by

ŝf (t) = k T f K m + K f + K n -1 x, (10)
where

k f = [k f (t, T 1 ), • • • , k f (t, T M )] T .

Numerical experiments

In this section we first investigate the performance of the proposed method on synthetic data to denoise ECG (Section 4.1). An illustration of f-ECG extraction is then provided on real data (Section 4.2).

Synthetic data: ECG denoising

The performance of the proposed algorithm to denoise ECG is assessed. In the first experiment, each beat of the ECG signal is generated by model (1). To mimic the variability presented in a real ECG, the waves amplitudes and P-R and R-T intervals are randomly changed (3%) around their average values. The ECG signal is then obtained as the summation of several beats with random global amplitudes and random R-R intervals. To ensure the consistency of the results, the whole procedure has been repeated one thousand times by regenerating all random parameters of the signal and noise samples. In this experiment, 1500 samples are used with 15 heart beats simulated at 100Hz sampled frequency. It is worth noting that the proposed method does not assume that the maxima of the R peaks are located at observed samples but can also appear in between samples. The proposed method is compared to the extended Kalman filtering (EKF) and smoothing (EKS) [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF]. The state model is chosen equal to (1) (i.e. the same model than the one used to generate data) whose parameters are equal to average values.

Quantitative results are shown in Fig. 3 which compares the output signal-tonoise ratio (SNR) achieved after denoising versus different input SNRs. As one can see (Fig. 3(b)), the proposed method increases the SNR with a gain between 3dB to 18dB. Contrary to extended Kalman filtering, the proposed method always improves the SNR. Indeed, in the case of high input SNR, EKS and EKF deteriorate the SNR: this can be explained by the variability of the simulated ECG as confirmed by Fig. 3(a), since this phenomenon is not observed without variability. Moreover, one can see that the variability decreases the overall performance, but the proposed method keeps the best performance by a smaller decrease than EKS or EKF.

Real data: f-ECG extraction

In this section, we illustrate (Fig. 4) the proposed method to extract f-ECG from a single sensor on the well-known DaISy fetal ECG database [START_REF] De Moor | Daisy: A database for identification of systems[END_REF]. As one can see, the proposed method provides suitable estimations of both maternal and fetal ECG even when mother's and fetus's R peaks are concomitant (e.g., the fourth, seventh and tenth mother's beats). Moreover, a visual inspection of the residual noise n(t) = x(t) -ŝm (t) -ŝf (t) confirms the validity of the assumed modelling [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. Indeed, this residual noise is effectively composed of a smooth varying baseline (dark curve) related to the first term of covariance function (4) plus a quasi white noise (validated by its covariance function empirical estimation). Moreover, both contributions are decorrelated with the mother's and fetus's ECG signals.

Conclusions and perspectives

In this paper, a non-parametric model of ECG signals is derived. By considering them as second order processes, which are fully defined by their mean and covariance functions, one can model a large class of signals with few hyper-parameters. From this modelling, denoising or extraction methods are directly obtained as the maximization of the posterior distribution. Numerical experiments show that the proposed method outperforms an extended Kalman filtering especially in presence of slightly random state parameters. Indeed, Gaussian processes realize a tradeoff between the suitable description of the signal by its second order statistics and its intrinsic variabilities. Finally, the main advantage of the proposed method is its flexibility and it provides a mix between purely data based methods as principal component analysis and parametric model based methods as Kalman filtering. Future work will deal with a computationally efficient implementation of hyper-parameters estimation of the proposed method as well as an online implementation.

Fig. 1 .

 1 Fig. 1. Modelling (1) of the amplitude of one beat.

Fig. 3 .

 3 Fig. 3. ECG denoising: output SNRs vs. the input SNR without Fig. 3(a) and with Fig. 3(b) parameters variability. In the two figures, the black line corresponds to the same input and output SNRs. In each case, the median are plotted, as well as the first and last quartiles as error bars.

  Fig.4. Fetal ECG extraction. Top to bottom: recorded signal x(t), estimated mother's ECG ŝm(t), estimated fetal's ECG ŝf (t) and residual noise n(t) (light gray curve) with estimated baseline (dark curve), respectively.

  4. Fetal ECG extraction. Top to bottom: recorded signal x(t), estimated mother's ECG ŝm(t), estimated fetal's ECG ŝf (t) and residual noise n(t) (light gray curve) with estimated baseline (dark curve), respectively.

⋆ Christian Jutten is also with Institut Universitaire de France.


Abdelmounime El Magri
email: a.elmagri@yahoo.fr

Fouad Giri

Abderrahim Elfadili

Luc Dugard

El Fadili

Wind Sensorless Control of Wind Energy Conversion System with PMS Generator

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

ost wind turbines currently in operation are constrained to work at constant speeds (with a tolerable variation of 1-2%) despite wind speed variations. In such operation conditions, the aerodynamic performance is optimal (in the sense of maximum wind energy extraction) only for one wind speed value depending on the considered value of turbine speed. Achieving maximum wind energy extraction in presence of varying wind speed conditions needs a varying turbine speed operation mode. Specifically, the turbine rotor velocity must be controlled so that its power-speed working point is constantly maintained near the optimal position (Fig. 4). This control objective is commonly referred to 'maximum power point tracking (MPPT)' and its achievement guarantees optimal aerodynamic efficiency.

In this paper, the focus is made on the wind energy conversion system Fig 1. This includes a permanent magnet synchronous generator (PMSG) that converts wind turbine power to output voltage whose amplitude and frequency vary with wind speed. PMSG generators offer several benefits, in wind power applications due to their high power density, high efficiency (as the copper losses in the rotor disappear), and absence of gearbox and reduced active weight. These features make it possible to achieve with PMSG's high varying speed control performance and highly reliable operation conditions (reduced need for maintenance). In the light of the above discussion, the global system including the wind turbine, the PMSG and the AC/DC/AC power converter has to be controlled to achieve accurate referencespeed tracking. The rotor speed reference ( ref  ) must be updated online, tacking into account the current wind velocity ( win v ), so that the MPPT requirement is realized.

Existing MPPT methods can be divided in two classes. The first one includes those methods using explicitly the turbine power characteristic, e.g. [START_REF] Senjyu | Sensor-less maximum power point tracking control for wind generation system with squirrel cage induction generator[END_REF]. The second class includes methods based on extremum-seeking techniques not necessitating the knowledge of turbine characteristics, e.g. [START_REF] Hong | MPPT for PM wind generator using gradient approximation[END_REF] [5] [START_REF] Kesraoui | Maximum power point tracker of wind energy conversion system[END_REF]. In fact, the required wind speed measurement is a sort of average value of wind speed along the turbine blade. Besides this is not easy to measure, the use of wind sensors entails additional cost and inevitably reduce reliability.

Designing rotor speed reference update laws that are able to meet the MPPT requirement without necessitating wind velocity sensors is a crucial issue that has yet to be solved. A solution is presently developed based on the wind turbine aerofoil characteristic curve. The obtained online speedreference optimizer represents a key feature of our control strategy which, besides, involves turbine rotor speed regulation, DC link voltage regulation and adaptive state observation. Speed regulation is performed using a nonlinear regulator which, in addition, regulates the d-component of the stator current to zero optimizing thus the delivered stator current. The voltage at the DC link between the AC/DC rectifier and inverter is regulated to a constant reference value (Fig. 2). In fact, the DC link voltage reference value should equal the nominal value of the PMS Generator. This regulation loop controls the reactive power control delivered to the grid. All previous regulation loops are designed using the backstepping technique [START_REF] Krstic | Nonlinear and adaptive control design[END_REF]. The nonlinear multiloop (speed, current, voltage) controller thus obtained (Fig. 3) is formally shown, using Lyapunov stability, to meet its control objectives.

II. MODELING THE ASSOCIATION SYNCHRONOUS GENERATOR-AC/DC/AC CONVERTER A. Modeling of the combination 'PMSG-DC/AC converter'

The controlled system is illustrated by Fig 1 . It includes a combination 'synchronous Generator-rectifier, on one hand, and a tri-phase DC/AC inverter, the other hand. The rectifier is a AC/DC converter operating, like the DC/AC rectifier, according to the known Pulse Wide Modulation (PWM) principle.

Wind Sensorless Control of Wind Energy Conversion System with PMS Generator

A. El Magri, F. Giri, A. Elfadili, L. Dugard M Such modeling is generally performed in the d-q rotating reference frame because the resulting components sd i and sq i turn out to be DC currents. It is shown in [START_REF] Magri | Nonlinear Control of Associations Including Wind Turbine, PMSG and AC/DC/AC Converters -Speed regulation and power factor correction[END_REF] that the association PMSG-PWM DC/AC converter model, expressed in the d-q coordinates, can be given by the following state space form:

J T x J K x J F dt dx g M     2 1 1 (3a) dc s s M s s v u L x L K x x p x L R dt dx 1 1 3 1 2 2 1      (3b) dc s s s v u L x x p x L R dt dx 2 2 1 3 3 1     (

,

represent the average d-and q-axis of the three-phase duty ratio system ) , , ( . with ) (z denote the average value on the modulation(MLI) period of ) (z .

B. DC/AC inverter modeling

The inverter circuit (DC/AC) is presented in Fig. 1. The power supply net is connected to a converter which consists of a three phase converter that has 6 semiconductors (IGBTs with anti-parallel diodes for bidirectional current flow mode) displayed in three legs 1, 2 and 3. The 6 semiconductors are considered as ideal switches. Only one switch on the same leg can be conducting at the same time.

Applying Kirchhoff's laws, this subsystem is described by the following set of differential equations: To simplify the three-phase representation (4a) for the synthesis of control laws, the Park transformation is invoked again.

e k v i dt d L   (4a) i r dc i i dt dv C   2 (4b)     123 123 n T i i k i  (4c) where     T n n n n i i i i
dc nd nq n nd nd v u L i E L dt di 0 0 1 1      (6a) dc nq nq n nq nq v u L i E L dt di 0 0 1 1      (6b) i R dc i i dt dv C   2 (6c)
where (

nq nd E E ,
), (

nq nd i i , ) and ( nq nd u u ,
) denote the averaged network voltage and current and input control of the inverter in dq-coordinate (Park's transformation).

The power absorbed by the DC/AC converter is given by the well known expression

dc i Load v i P 
. On the other hand, the power released by the network is given by

    nq nq nd nd n T n OUT i E i E i e P    123 123
. Using the power conservation principle, one has

OUT Load P P 
or, equivalently

nq nq nd nd dc i i E i E v i  
. Also, from (6a-b) one gets that:

  R dc nq nq nd nd dc dc i v i E C i E C dt dv v     1 1 2 (7a) dc nd nq n nd nd v u L i E L dt di 0 1 1      (7b) dc nq nd n nq nq v u L i E L dt di 0 0 1 1      (7c)
Let us introduce the state variables

2 4 dc v x  , nd i x  5 , nq i x  6 and nq nd u u u u   4 3
, represent the average d-and q-axis components of the three-phase duty ratio system i k . The state space equations obtained up to now are put together to get a state-space model of the whole system including the AC/DC/AC converters combined with the synchronous generator. For convenience, the whole model is rewritten here for future reference: CO3. The inverter output currents (

J T x J K x J F dt dx g M     2 1 1 (9a) dc s s M s s v u L x L K x x p x L R dt dx 1 1 3 1 2 2 1      (9b) dc s s s v u L x x p x L R dt dx 2 2 1 3 3 1     (9c)   R dc nq nd i v x E C x E C dt dx     6 5 4 1 1 (9d) cd n nd v u L x E L dt dx 3 0 6 0 5 1 1      (9e) cd n nq v u L x E L dt dx 4 0 5 0 6 1 1      (9f)
3 2 1 , , n n n i i i
) must be sinusoidal with the same frequency as the supplied power grid and the reactive power in the AC grid must be well regulated.

CO4. Controlling the continuous voltage

dc v making it track a given reference signal dcref v
. This generally is set to a constant value equal to the nominal voltage entering the converter and machine.

CO5. Regulate the current sd i to a reference value dref i , preferably equal to zero in order to optimize the stator current (see e.g. [START_REF] Muhammad | Power electronics handbook[END_REF]).

IV. STATE FEEDBACK CONTROLLER DESIGN

A. Wind sensorless rotor speed reference optimization

The optimization strategy is to search a reference speed opt  , in order to achieve optimal speed ratio working conditions of the wind turbine to capture the maximum energy from the wind. The wind power acting on the swept area of the blade  is a function of the air density  ( The summits of these curves give the maximum 'extractable' power opt P and so represent the optimal points. Each one of these points is characterized by the optimal speed opt  . It is readily seen from Fig. 4 that for any wind velocity value, say i win v , there is a unique couple ( i  , i P ) that involves the largest extractable power. The set of all such optimal couples ( i  , i P ) is represented by the blue curve in Fig. 4. A number of such couples have been collected from Fig. 4 and interpolated to get a polynomial function ) (

opt opt P F   . Let
the obtained polynomial be denoted:

0 1 1 1 ... ) ( h P h P h P h P F n n n n        (11)
The function (.) F defines the speed-reference optimizer (Fig. 5). Indeed, suppose that, at some instant, the wind velocity is s m v win / 14 0  and the rotor speed is 0  . One can see on Fig. 4 which transmitted wind power 0 P corresponds to this couple ) , (

0 0  win v
. The point is that 0 P is easily online computed. Given only the value of 0 P , the speedreference optimizer gives a new rotor speed reference value ) ( 0

1 P F  
. This corresponds to point A on Fig. 4. Assume that, a speed regulator (that has to be determined) is available that makes the machine rotor rotate at the new speed reference i.e. 1    . Then, according to Fig. 4, the wind turbine provides a new extractable power value equal to 1 P . Then, the speed-reference optimizer will suggest (to the speed regulator) a new speed reference 2  (point B). This process will continue until the achievement of the optimal point ( opt opt P 

,

) (point C on Fig. 4).

Speed regulator design for synchronous generator

The regulator design is based on equations (9a-b) where the input signal 1 u stands as the actual input, in order to guarantee speed reference tracking. Following the backstepping technique [START_REF] Krstic | Nonlinear and adaptive control design[END_REF], let 1 z denotes the speed tracking error:

* 1 1 1 1 x x x z ref      (12)
In view of (9a), the above error undergoes the following equation:

* 1 2 1 1 / ) ( x J T x K x F z g M        (
  * 1 1 1 1 * ) / ( ) / ( x J T x J F z c g        (with 0 1  c ) (14) Indeed, if *    one have 1 1 1 z c z   
which clearly is asymptotically stable with respect the Lyapunov function:

In effect, one then has:

0 2 1 1 1 1 1     z c z z V   (16) As 2 ) / ( x J K M   
, is just a virtual control input, one cannot set *    . Nevertheless, the above expression of *  retained as a first stabilization function and a new error is introduced:

* 2     z (17)
Using ( 14)-( 17), it follows from (13) that the 1 z -dynamics undergoes the following equation:

2 1 1 1 z z c z     ( 18
)
The next step consists in determining the control input 1 u so that the ( 21 , z z

) error system is asymptotically stable. First, let us obtain the trajectory of the error 2 z . Deriving this with respect to time and using (17) gives:

* 2 2 ) / (        x J K z M (19)
Using ( 14) and (9a-b) in (19), one gets:

dc s M v u JL K z c z c t x z 1 2 1 1 2 1 2 ) / ( ) , (       (20) where     * 1 2 2 2 1 2 1 3 1 2 / / / . ) / ( ) , ( x J T J FT J x FK x F JL K x K x x L p x R t x g g M s M M s s             (21)
The error equations ( 18) and (20) are given in more compact form:

2 1 1 1 z z c z     (22a) dc s M v u JL K z c z c t x z 1 2 1 1 2 1 2 ) / ( ) , (       (22b)
To determine a stabilizing control law for (22b), let us consider the quadratic Lyapunov function candidate:

2 2 2 1 2 2 1 2 5 . 0 5 . 0 5 . 0 z z z V V     (23)
Using (18), the time derivative of 2 V can be rewritten as:

2 2 2 1 2 1 1 2 z z z z z c V       (24a)
This shows that, for the ( 21 , z z )-system to be globally asymptotically stable, it is sufficient to choose the control 1 u so that:

2 2 2 2 1 1 2 z c z c V     (24b) where 0 2  c
is a new design parameter. In view of (24a), equation (24b) is ensured if:

1 2 2 2 z z c z     (25)
Comparing ( 25) and (22b) yields the following backstepping control law:

    dc M s v K t x z c z c c L J u 3 / ) , ( ) 1 ( ) ( 1 2 1 2 2 1 1        (26)
C. d-axis current regulation The d-axis current 3

x undergoes equation (9c) in which the following quantity:

s dc L v u x x p v / 2 1 2   (27)
As the reference signal dref i is null, it follows that the tracking error

dref i x z   3 3 undergoes the equation: v z L R z s s    3 3 ) / (  (28)
To get a stabilizing control signal for this first-order system, consider the quadratic Lyapunov function

2 3 3 5 . 0 z V  . It is easily checked that, if the virtual control is let to be 3 3 ) / ( z c L R v s s     where 0 3  c
, then:

2 3 3 3 z c V    ( 30
)
which is negative definite. Furthermore, substituting (29) in (28), one gets the closed-loop equation:

3 3 3 z c z    (31)
Now, it is readily observed that the actual control input is obtained substituting (29) into ( 27) and solving the resulting equation for 2 u . Doing so, one gets:

    dc s s s v L x x p z L R z c u / / 1 2 3 3 3 2    (32)
From ( 22a), ( 25) and (31),it is clear that the subsystem is globally exponentially stable.

D. Reactive power and DC voltage controller

In controlling a PFC, the main objective is to obtain a sinusoidal output current and the injection or extraction a desired reactive power in the electric network. The continuous voltage x and the second ensures the injection of the desired reactive power.

1) DC voltage loop

Based on equations (9d-e), first equation involving the control input 3 u will now be designed, using the backstepping technique [START_REF] Krstic | Nonlinear and adaptive control design[END_REF], that the squared DC-link voltage . As the system (9d-e) is of relative degree 2, the design towards that equation is performed in two steps.

Step 1. Let 4 z denote the speed tracking error:

* 4 4 4 x x z   (34)
In view of (9d), the above error undergoes the following equation:

* 4 3 1 3 1 5 4 ) , ( / x z x C x E z i i nd            (35)       C x x x p z L R z c L x t x z c z c c K L J x E s s s M s nq / ) / ( ) , ( ) ( ) 3 / ( ) ( 3 1 2 3 3 3 2 1 2 1 2 2 1 6             In (35), the quantity C x E nd / 5 1    stands up as a (virtual)
control input for the 4 z -dynamics because the actual control input 3 u act on 4 z indirectly through 1  . Following the backstepping design technique, the Lyapunov function candidate is considered as: . Deriving 4 V along the trajectory of (35) yields:

  * 4 5 4 4 4 4 ) , ( / x z x C x E z z z V nd          (36)
This suggests for the (virtual control) 1

 the following control law:

* 4 4 4 * 1 ) , ( x z x z c        (37) with 0 4  c a design parameter. Indeed, substituting * 1  to C x E ) , ( / x z x z c C x E z nd             (40)
Using (39) and (9d-e) in (40) yields:

cd nd i i v u L C E z x z 3 0 5 1 6 1 1 5 ) , (         (41) With 6 0 2 * 4 4 4 1 ) , ( ) , ( x C E L C E x z x z c z x n nd nd            
To determine a stabilizing control law for (9d-e), let us consider the quadratic Lyapunov function candidate:

2 5 2 4 5 5 . 0 5 . 0 z z V   (42a)
Using ( 39)-(41), one gets from (42a) that:

  cd nd i i v u L C E z x z z z c V 3 0 5 1 6 1 1 4 5 2 4 4 5 ) / ( ) , (            (42b)
This suggests for the control variable 3 u the following choice: ;

  cd nd i i v E L C z x z z c u / ) , ( 0
z z c z z z c z         (45)
2) Reactive Power loop Here, the focus is made on the control objective CO3 that involves the reactive power n Q which is required to track its reference * n Q . The electrical reactive power injected in the grid is given by ) , (

n nq nd n Q E x E x x        .
As the equation ( 46) is a first order, it can be (globally asymptotically) stabilized using a simple proportional control law:

  ) , ( / 6 5 2 6 6 0 3 4 x x z c L u E u E v nq nd cd      (47) with 0 6  c
. Then the control law 4 is given as:

    nd nq cd E u E v x x z c L u / / ) , ( 3 6 5 2 6 6 0 4      (48)
It can be easily checked that the dynamic of 6 z undergoes the following equation:

6 6 6 z c z    (49)
Then, from (45) and (49), it is clear that the closed loop subsystem is globally exponentially stable and the error vector ( ) converges exponentially fast to zero, whatever the initial conditions.

V. SIMULATION

The closed loop system is simulated, within the Matlab/Simulink environment, using the electro-mechanical characteristics of Table 1. The controller, mathematically defined by (26), (32), (43) and ( 48) is simulated. The used numerical values for the design parameters are chosen as follows: that proved to be suitable. The controller performances will be evaluated in presence of (time-varying) wind velocity . According to the control design, the remaining closed loop inputs are kept constant:

V V A i dcref ref sd 700 , 0   and VAR Q ref n 0  .

B. Illustration of the controller performances

In this subsection, the controller performances are illustrated by the curves in Fig. 9. The tracking quality is quite satisfactory as the response time (after each change in the wind speed) is less than s 1 . Curves (10a) and (10b) respectively show the reactive power injected in the three-phase network n Q (equals zero) and electric power n P , produced by the machine and transferred to the grid through the tri-phase inverter.

Fig. 11a shows that the DC-link voltage dc v is tightly regulated: it quickly settles down after each change in the wind speed. The wave frame of the line current 1 n i is showed in fig 11b. It is seen that the current amplitude changes whenever the wind velocity varies. The current remains (almost) all time sinusoidal and in phase with the network voltage complying with the PFC requirement.

VI. CONCLUSION

We have addressed the problem of sensor less control of the wind energy conversion system. Maximum wind energy extraction is achieved by running the wind turbine generator in variable-speed mode without using the sensor wind velocity. The controlled system is an association including wind turbine, permanent magnet synchronous aero-generator and AC/DC/AC converter connected with a tri-phase net work. The system dynamics have been described by the averaged sixth order nonlinear state-space model (9a-f). The multi loops nonlinear controller defined by the control laws (26), (32), ( 43) and (48) has been designed. The backstepping design technique is used. The controller has been designed to: (i) satisfactory rotor speed reference tracking for extracting maximum power; (ii) tight regulation of the stator d-axis; (iii) power factor correction; (iv) well regulated DC-link voltage (vdc). These results have been confirmed by a simulation study.

  is the sinusoidal threephase net voltages (with known constant frequency n  ), dc v denotes the voltage in capacitor 2C, i i designates the input current inverter, and i k is the switch position function taking values in the discrete set { 1 , 0 }.

Fig. 1 .

 1 Fig.1. The AC/DC/AC Converter power circuit in WECS

C

  is a non-linear function of the tip speed-ratiowin v R /   (where R is the turbine radius), which depends on the wind velocity and the rotation speed of the generator rotor  ( s rad / ). Fig.4represents the transmitted power according to the rotor PMSG speed for various values of the wind speed.

Fig . 4 :

 4 Fig .4 : Turbine Power Characteristics (Pitch angle=0°).

  lead to two control loops. The first loop ensures the regulation of the DC voltage 4

  is a new design parameter. Indeed, substituting (43) in (42b) yields:

  Figs 9a and 9b show that the machine speed, their respective references.

Fig. 9 .

 9 Fig. 9. a) Reference and real rotor speed (m/s); b) d-axis current sd i (A).

  Fig. 10. a) Reactive power n Q

Fig. 11 .

 11 Fig. 11. a) DC link voltage; b) wave frame of the line input current.

Fig. 8 .

 8 Fig. 8. a) Wind speed; b) Optimized Rotor Speed reference


Miroslav Šiman
email: siman@utia.cas.cz

SOME STATISTICS BASED ON PROJECTION PURSUIT IN A GENERAL REGRESSION CONTEXT

Keywords: exact computation, projection pursuit, quantile regression, multivariate quantile, projection depth, regression depth exact computation, projection pursuit, quantile regression, multivariate quantile, projection depth, regression depth

whether they are published or not. The documents may come

Recently, a new theory of directional (regression) quantiles has been developed, see Paindaveine andŠiman (2010b, 2010c) for the algorithms and Hallin et al. (2010a), Paindaveine and Šiman (2010a), and [START_REF] Kong | Quantile tomography: using quantiles with multivariate 11 Page 12 of 13[END_REF] for the theory. It gives rise to many promising directional statistics but their directional dependence makes it difficult to exactly compute their suprema or infima over infinitely many unit directions, which would be very useful for overall statistical inference. This paper addresses the problem and solves it in some important cases.

Besides, the projectional quantile regression of Paindaveine andŠiman (2010a, 2010c) 1 appears very useful for the exact computation of many statistics based on projection pursuit methodology where the same problem arises as well; see e.g. [START_REF] Friedman | Exploratory projection pursuit[END_REF] and [START_REF] Huber | Projection pursuit[END_REF]. We illustrate this on a few important examples including multiple-output regression generalizations of location projection depth, skewness and kurtosis. These concrete multivariate regression extensions are obtained quite naturally by means of regression quantiles, introduced in Koenker and Bassett (1978) and masterly reviewed in [START_REF] Koenker | Quantile Regression[END_REF], but they seem, quite surprisingly, not to have yet been explicitly considered in the literature.

The projection pursuit plays a crucial role in the definitions of various statistical depths including halfspace depth, weighted halfspace depth and projection depth; see [START_REF] Rousseeuw | The depth function of a population distribution[END_REF], [START_REF] Hlubinka | Weighted halfspace depth[END_REF] and [START_REF] Zuo | Projection-based depth functions and associated medians[END_REF] with references therein, respectively.

The projection depth leads the way and beats many of its location depth competitors in numerous respects, see Zuo andSerfling (2000a, 2000b). It has already been generalized even to regression fits (i.e. hyperplanes) in Zuo and Cui (2004). Finally, it is here extended to arbitrary points in a general multiple-output regression context, which seems more intuitive, more useful and more promising for future development; see also [START_REF] Liu | Linear fitting by simplicial intercept depth (SID): reflection invariance and robustness[END_REF], Paindaveine and Šiman (2010a), and [START_REF] Rousseeuw | Regression depth[END_REF] for other notions of depth in a regression framework.

Exact computation of some special forms of location projection depth has already been briefly outlined in [START_REF] Zuo | Projection-based affine equivariant multivariate location estimators with the best possible finite sample breakdown point[END_REF] for odd numbers of bivariate data points. Besides, there is another paper dealing with exact computation of bivariate projection depth that is accompanied with software implementation in R, see [START_REF] Zuo | Exact computation of the bivariate projection depth and Stahel-Donoho estimator[END_REF] and Zuo and Ye (2009). Although that paper has not yet been published and made available, its title indicates that it deals only with the bivariate and location case. In other words, the present paper almost surely brings up some new ideas regarding this computational issue, at least in its treatment of projection regression depth in spaces of general dimension. Needless to say that the results presented here can lead to a real breakthrough in the dissemination of projection depth and related methods. of various statistics based on univariate (regression) quantiles, see e.g. [START_REF] Kim | On more robust estimation of skewness and kurtosis[END_REF] and [START_REF] White | Modeling autoregressive conditional skewness and kurtosis with multi-quantile CAViaR[END_REF]. Two of these univariate quantile-based shape characteristics are considered even here and extended by means of projection pursuit to the multiple-output context in the same way as their moment-based competitors in [START_REF] Malkovich | On tests of multivariate normality[END_REF].

These extensions are then shown exactly computable even when the output dimension exceeds two or three, which may lead to many applications, including new tests for (conditional) multivariate symmetry or normality. This paper proceeds as follows. Section 2 introduces one fundamental lemma that is then heavily applied to the multiple-output quantile regression in Section 3 and to the projectional quantile regression in Section 4. Its subsections discuss the exact computation of the multivariate shape statistics and projection depth in a general regression context.

Section 5 concludes with some final remarks.

This brief communication blends so many various theoretical concepts together that their joint detail presentation in the Introduction would only confuse the reader. This is why their definitions are left only for the sections where they are employed. Sections 3 and 4 are independent, which explains why they sometimes use the same notation for similar but different entities.

THEORY

Hereafter, we adopt basic terminology of [START_REF] Boyd | Convex optimization[END_REF].

Definition. Function F on R m is directional (or, scale invariant) if F (cu) = F (u) ∀c > 0 ∀u ∈ R m , and is sign-directional if sign(F (cu)) = sign(F (u)) ∀c > 0 ∀u ∈ R m .
Lemma 1. Let us assume that the whole space R m , m > 1, can be partitioned into a finite number of non-degenerate closed convex polyhedral cones C q with disjunct interiors, say R m = N q=1 C q , and let us write Γ := {C q : q = 1, . . . , N } for this conic segmentation, determined by its finite vertex set V := N q=1 V q where V q := {u ∈ R m : u = 1, u lies in an edge of C q }. 2. If F is directional and quasiconcave on each

C q ∈ Γ, then inf u =1 F (u) = min u∈V F (u).
3. If F is sign-directional and quasiconcave on each

C q ∈ Γ, then F (u) > 0 for all u ∈ S m-1 := {u ∈ R m : u = 1} if and only if min u∈V F (u) > 0.
All these statements hold even when the assumptions are not met by F itself, but only by

a function u → F (u)g( u )
where g is an arbitrary positive function on (0, ∞) such as

g(x) = x.
Proof. The surface D of the convex hull of V consists of q closed convex polyhedral facets

D q := { i λ i v i : i λ i = 1, λ i ≥ 0, v i ∈ V q } determined by C q , D = N q=1 D q . If F is directional, then sup u =1 F (u) = sup u∈D F (u) = max q=1,...,N sup u∈Dq F (u) and also inf u =1 F (u) = inf u∈D F (u) = min q=1,...,N inf u∈Dq F (u). If F is quasiconvex on each C q ,
then it is also quasiconvex on each D q and sup u∈Dq F (u) = max u∈Vq F (u), see [START_REF] Boyd | Convex optimization[END_REF]. This proves the first statement.

If F is quasiconcave on each C q , then -F is quasiconvex there by definition, which implies

inf u∈Dq F (u) = min u∈Vq F (u) and so inf u =1 F (u) = min u∈V F (u).
It also holds for any sign-directional function F that F (u) > 0 on S m-1 if and only if

F (u) > 0 on D, which is therefore equivalent to inf u∈D F (u) = min q=1,...,N min u∈Vq F (u) > 0 thanks to the quasiconcavity of F on each C q .
Finally, F (u) and F (u)g( u ) differ on S m-1 only by a positive multiplicative constant, which concludes the proof.

Lemma 1 will often be invoked below, but only with either g( u ) = 1 or g( u ) = u (convex in u), which should always be kept in mind. The notation introduced in Lemma 1 will also be used hereinafter, together with the following definition.

Definition. Let us consider some indexed conic segmentations Γ(i), i ∈ I, with their vertex sets V(i), i ∈ I. Then their common refinement Γ(I) is defined as the conic segmentation determined by V(I) := i∈I V(i).

All applications of Lemma 1 are likely to use the theory of quasiconvex functions, especially the following composition rules from [START_REF] Boyd | Convex optimization[END_REF].

. if g : R m → R is quasiconvex and h : R → R is nondecreasing, then f = h • g is quasiconvex 2. if g is quasiconvex, then f (u) := g(Au + b) is quasiconvex 3. if g is quasiconvex, then f (u) := g((Au + b)/(c u + d)) is quasiconvex on the set {u ∈ R m : c u + d > 0} 4. if g is a convex function, h is a concave function, and both g(u) > 0 and h(u) > 0 on a convex set C, then f (u) = g(u)/h(u) is quasiconvex on C
where all the compositions are well-defined by assumption.

For the sake of simplicity, only suprema of directional statistics will be discussed hereinafter (without any significant loss of generality). This is why these composition rules are formulated only for quasiconvex functions here. They can be used to check or extend all the claims made below.

Let us consider m-dimensional responses Y i associated with (p+1)-dimensional regressors

X i = (1, Z i ) , i = 1, . . . , n > p + 2,
and assume that the random sample (Y i , Z i ) ∈ R m+p , i = 1, . . . , n, comes from a continuous distribution. Section 3 and 4 discuss the consequences of Lemma 1 in the empirical regression context with a suitable conic segmentation. where w i > 0, i = 1, . . . , n, are some positive weights, τ ∈ (0, 1) stands for any fixed quantile level, ρ τ (x) = x(τ -I(x < 0)) is the well-known quantile check function and

MULTIPLE-OUTPUT QUANTILE REGRESSION

Ψ τ u = n i=1 w i ρ τ (b τ u Y i -a τ u X i ) (2)
denotes the optimal value of the objective function. Typically, w i = 1, i = 1, . . . , n, as in Hallin et al. (2010a).

It follows from [START_REF] Paindaveine | Computing multiple-output regression quantile regions[END_REF] that there exists a finite conic segmentation Γ(τ ) = {C q (τ ) : q = 1, . . . N τ } of R m almost surely for all but a finite number of τ 's such that it meets the assumptions of Lemma 1 and (a τ u , b τ u ) = (a τ q , b τ q ) /t τ q u, Ψ τ u = ψ τ q /t τ q u and t τ q u > 0 for any u in any C q (τ ) where a τ q , b τ q , t τ q and ψ τ q are some constant vectors or scalars possibly different on each C q (τ ). Paindaveine and Šiman (2010b) also provide an algorithm (and its Matlab implementation by means of parametric programming) that can find such conic segmentation Γ(τ ), its vertex set V(τ ) and corresponding characteristics a τ q , b τ q , t τ q and ψ τ q , q = 1, . . . , N τ . Clearly, (a

τ u , b τ u ) = -(a (1-τ )(-u) , b (1-τ )(-u)
) and thus V(τ ) = -V(1 -τ ). Note as well that Γ(T ) is easy to construct for any finite index set T of τ 's, at least for m = 2, and that a τ u , b τ u and Ψ τ u can be computed for any u ∈ R m \ {0} by means of available a τ q , b τ q , t τ q and ψ τ q , which may be more effective then their straightforward computation from scratch.

Simple Inference in Multiple-Output Quantile Regression

Hallin et al. (2010a) already showed that both a τ u , b τ u and Ψ τ u can be very useful for statistical inference and that they can be combined into directional statistics that are constant on each C q (τ ) and therefore easy to maximize or minimize over S m-1 .

Lemma 1 can be applied in this context to extend these results considerably. For example, it implies sup u =1 G(u) = max u∈V(τ ) G(u) where G(u) can stand e.g. for a τ u , b τ u and Ψ τ u . Similarly, it guarantees for any suitable τ 1 , τ 2 ∈ (0, 1) that Ψ τ 1 u /Ψ τ 2 u has its supremum over S m-1 respectively equal to its maximum over the vertex set V(T ) corresponding to

Γ(T ), T = {τ 1 , τ 2 }.
w i ρ τ (u Y i -a X i ) (3)
where w i > 0, i = 1, . . . , n, are some positive weights, τ ∈ (0, 1) denotes the quantile level,

ρ τ (x) = x(τ -I(x < 0)
) is the well-known quantile check function and

Ψ τ u := n i=1 w i ρ τ (u Y i -a τ u X i ) (4)
stands for the optimal value of the objective function. Often w i = 1, i = 1, . . . , n, and then a τ u is nothing but the standard sample regression τ -quantile of projections u Y i 's.

It was shown in Paindaveine andŠiman (2010a, 2010c) that there exists a finite conic segmentation Γ(τ ) = {C q : q = 1, . . . , N τ } of R m almost surely for all but a finite number of τ 's such that each C q is a non-degenerate closed convex polyhedral cone where a τ u and Ψ τ u are linear in u, i.e. a τ u = A τ q u, A τ q ∈ R (p+1)×m , and Ψ τ u = ψ τ q u, ψ τ q ∈ R m , for any u ∈ C q (τ ). Paindaveine and Šiman (2010c) also provide an algorithm (and its Matlab implementation by means of parametric programming) that can find such a conic segmentation Γ(τ ), its vertex set V(τ ) and corresponding matrices A τ q and vectors ψ τ q , q = 1, . . . , N τ . Clearly, a τ u = -a (1-τ )(-u) and thus V(τ ) = -V(1 -τ ). Note as well that Γ(T ) is easy to construct for any index set T of τ 's, at least for m = 2, and that a τ u and Ψ τ u can be computed for any u ∈ R m \ {0} by means of available A τ q and ψ τ q , which may be more effective then their straightforward computation from scratch.

Simple Inference in Projectional Quantile Regression

Paindaveine and Šiman (2010a) already showed that both a τ u and Ψ τ u can be very useful for statistical inference. Now Lemma 1 makes such inference possible and simple in many cases. For example, its direct application guarantees for any suitable τ 1 , τ 2 ∈ (0, 1) that

sup u =1 G(u) = max u∈V(τ 1 )∪V(τ 2 ) G(u) where G(u) can stand for Ψ τ 1 u /Ψ τ 2 u , a τ 1 u /Ψ τ 2 u , or 1/Ψ τ 1 u , among others.

Ratios and Inverses of Regression L-Statistics

Lemma 1 can also be applied to ratios R u (x) of directional regression L-statistics:

R u (x) = τ ∈T N c τ a τ u x τ ∈T D d τ a τ u x := N u (x) D u (x) .
The denominator D u (x) is linear in u on each cone from Γ(T D ) and therefore it is positive for all u ∈ S m-1 if and only if its minimum over V(T D ) is greater than zero, which can be easily checked. If this is the case, then sup u =1 R u (x) can be computed as max u∈V(

T D ∪T N ) R u (x) because R u (x) is then even quasilinear in u on each cone from Γ(T D ∪ T N ). Lemma 1 also implies sup u∈S 1/D u (x) = max u∈V(T D ) 1/D u (x) if D u (x) > 0 ∀u ∈ S m-1 .
For example, given x ∈ 1 × R p , we can consider two successful quantile-based shape measures from Kim and White (2006) adjusted to the general directional regression framework considered here, namely the skewness coefficient ξ u (x) and the kurtosis measure κ u (x):

ξ u (x) = (a (3/4)u -a (2/4)u )x -(a (2/4)u -a (1/4)u )x (a (3/4)u -a (1/4)u )x , κ u (x) = (a (7/8)u -a (5/8)u )x + (a (3/8)u -a (1/8)u )x (a (6/8)u -a (2/8)u )x .
Both these statistics can be rewritten in the form of R u (x). Indeed, the former would Let us generalize location projection depth, thoroughly investigated in [START_REF] Zuo | Projection-based depth functions and associated medians[END_REF], to the points in the general regression setup. In the empirical case discussed here, the projection regression depth P D n (y, x) of y ∈ R m given x = (1, z ) ∈ 1 × R p may be defined as

P D n (y, x) = 1 1 + sup u =1 g n (y, x, u) with g n (y, x, u) = u y -µ u (x) σ u (x)
where µ u (x) and σ u (x) > 0 are some sample suitable regression measures of location and scatter of u Y given X, respectively.

In this example, let us consider only µ u (x) and σ u (x) in the form of regression L-statistics,

µ u (x) := τ ∈Tµ c τ a τ u x and σ u (x) := τ ∈Tσ d τ a τ u x,
where c τ ∈ R, τ ∈ T µ , and d τ ∈ R, τ ∈ T σ . This restriction does not seem too severe after considering the close asymptotic relationship between R-, M-, and L-statistics thoroughly described in [START_REF] Jurečková | Robust statistical procedures: asymptotics and interrelations[END_REF].

Generally speaking, the projection (regression) depth is difficult to compute fast and exactly for m > 1 because the supremum in its definition then involves infinitely many unit directions. Besides, its regression definition seems reasonable only if σ u (x) > 0 for any u ∈ S m-1 , which need not be true for p > 0 due to the quantile crossing phenomenon.

Nevertheless, Lemma 1 implies that σ u (x) > 0 for any u ∈ S m-1 if and only if min u∈V(Tσ) σ u (x) > 0, which can be easily checked. And if this condition holds, then g n (y, x, u) is directional and quasiconvex on each C q ∈ Γ(T µ ∪ T σ ), which further results in P D n (y, x) = 1 1 + max u∈V(Tµ∪Tσ) g n (y, x, u) .

Note that this formula can be used for fast and exact computation of projection depth even in the general regression context and for m > 2. The choice of T µ and T σ allows for a 9 compromise between robustness and speed as the computational time of finding Γ(τ ) and related quantities grows with decreasing |τ -0.5|, see [START_REF] Paindaveine | Computing multiple-output regression quantile regions from projection quantiles[END_REF].

FINAL REMARKS

The minimizers (a τ u , b τ u ) in Section 3 and a τ u in Section 4 need not be defined uniquely, not even in the purely location setup for integer values of nτ , but there is only a finite number of such exceptional quantile levels with probability one. Besides, this potential ambiguity is not too worrying because all such minimizers always form a convex polyhedral set that shrinks towards their unique population limits under very mild conditions, see Hallin et al. (2010a) and Chapter 4 in [START_REF] Koenker | Quantile Regression[END_REF] for a detailed discussion of these asymptotic properties. Anyway, such an ambiguity can often be fixed easily by a tiny perturbation of τ or the data points, which can have hardly any impact on practical applications. This is why these minimizers are considered uniquely defined throughout this paper.

Clearly, the theory presented here is especially useful when the number of cones in the conic segmentations is quite low. In Sections 3 and 4, this often happens for small n or extreme τ . Then even the search for approximate extremes on a sampled dense grid of directions would be much more complicated, in particular for m > 2.

It should also be pointed out that the weights w i 's used in Sections 3 and 4 can be really useful, both for handling multiple observations and for computing local constant counterparts to the regression quantiles of Sections 3 and 4, see [START_REF] Hallin | Local bilinear multiple-output quantile regression and regression depth[END_REF]. They would also lead to interesting local versions of the shape measures and projection regression depth considered in Section 4.

F

  Skewness and kurtosis of univariate (conditional) distributions can be described by means 2

4

  Page 5 of 13 URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.It holds for arbitrary parameters A, b, c and d of the right dimensions that 1

  Hallin et al. (2010a) introduced a new concept of directional multiple-output regression quantiles. In the general empirical case, they can be defined for any u ∈ R m \ {0} by means of the optimization problem(a τ u , b τ u ) = argmin (a ,b ) n i=1 w i ρ τ (b Y i -a X i ) subject to u b =

6

  Page can consider projectional regression quantiles from[START_REF] Kong | Quantile tomography: using quantiles with multivariate 11 Page 12 of 13[END_REF] andPaindaveine and Šiman (2010a). In the general empirical case, they can be defined for any u ∈ R m \ {0} by means of the optimization problem a τ u := argmin a∈R p+1 n i=1

7

  Page

  lead to T D (ξ) = {1/4, 3/4} and T N (ξ) = {1/4, 2/4, 3/4} while the latter would correspond to T D (κ) = {2/8, 6/8} and T N (κ) = {1/8, 3/8, 5/8, 7/8}. If their denominators are always positive, which holds almost trivially in the location case with p = 0, then we can meaningfully define overall skewness ξ(x) := sup u∈S ξ u (x) and overall kurtosis κ(x) := sup u∈S κ u (x) by means of projection pursuit and compute them as max u∈V(T D (ξ)∪T N (ξ)) ξ u (x) and max u∈V(T D (κ)∪T N (κ)) κ u (x), respectively.

ca Communications in Statistics -Simulation and Computation

	Page 10 of 13
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.

ca Communications in Statistics -Simulation and Computation

  Zuo, Y., and Cui, H. (2004). Statistical depth functions and some applications. Advances inMathematics, 33, 1-25. Zuo, Y., and Serfling, R. (2000a). On the performance of some robust nonparametric location measures relative to a general notion of multivariate symmetry. Journal of Statistical Planning and Inference, 84, 55-79. Zuo, Y., and Serfling, R. (2000b). General notions of statistical depth functions. Annals of Statistics, 28, 461-482. Zuo, Y., and Ye, X. (2009). ExPD2D: Exact computation of bivariate projection depth based on Fortran code. R package version 1.0.1, http://cran.r-project.org/package=ExPD2D.

	Page 14 of 13
	o F o F o
	r r r
	P P P
	e e r e e r e e r
	R e v i e w R e v i e w R e v i e w
	O n O n O n
	l y l y l y
	13

URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.

URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.caCommunications in Statistics -Simulation and Computation

Davy Paindaveine, Marc Hallin, Claude Adan, Nancy de Munck and Romy Genin for all the good they did for him (and for all the good he could learn from them) during his recent stay at Université Libre de Bruxelles.

ACKNOWLEDGEMENTS

This research work of Miroslav Šiman was supported by Project 1M06047 of the Ministry of Education, Youth and Sports of the Czech Republic. The author would also like to thank 10