Thibaut Raharijaona

Lorris Dola

Bruno Boisseau

John-Jairo Martinez-Molina

Nicolas Marchand

Franck Ruffier

Event-based speed control on a sensor-less miniature thruster

The main contribution of this paper is the experimental evaluation of the performance of different event-based control strategies in the case of a sensor-less miniature thruster.

While most control engineering applications considers periodic control with equidistant sample intervals, the event-based approach updates the control signal only when required. The performances in terms of thruster power consumption and control update number are evaluated by comparing the eventbased control system with respect to the time-based loop for the angular speed regulation of a sensor-less miniature thruster. The experimental results clearly show an important reduction of the control updates while the power consumptions are similar.

I. INTRODUCTION

Due to the necessity of saving the energy and increasing the autonomy of embedded systems, event-based control becomes an interesting research topic. While the classical socalled time-based framework of controlled systems consists in sampling the system uniformly in time with a constant sampling period and in computing and updating the control law every time instants, the event-based approach aims at updating asynchronously the control signal. The control signal is updated only when required. Some works addressed more recently event-based sampling where the control law is event-driven [START_REF] Astrom | Comparison of riemann and lebesgue sampling for first order stochastic systems[END_REF], [START_REF] Sandee | Eventdriven control as an opportunity in the multidisciplinary development of embedded controllers[END_REF], [START_REF] Lunze | A state-feedback approach to event-based control[END_REF], [START_REF] Anta | To sample or not to sample: Self-triggered control for nonlinear systems[END_REF].

For an application point of view, there are several possible ways to save the energy and improve the autonomy. The first way aims at reducing the computation load and the second way tends to reduce the actuation load. Low-power electronic components in all embedded and miniaturized mechatronic applications encourages companies to develop asynchronous versions of the existing time-triggered components. A significant power consumption reduction can be achieved by decreasing the samplings and consequently the CPU utilization. In [START_REF] Årzén | A simple event-based pid controller[END_REF], an event-based proportional integral derivative (PID) controller is presented with the simulation results applied to a double tank process. For stability reason, a maximal sampling period is introduced in the initial approach. In [START_REF] Durand | Further results on event-based pid controller[END_REF], this maximal period is removed and different event-based PID algorithms without safety limit condition were also developed. In [START_REF] Durand | Asyncar, a radio-controlled vehicle for asynchronous experiments implementation of an event-based cruise control[END_REF], the efficiency of an approach without safety limit condition is highlighted by implementing such a controller in a real-time testbed. A radio-controlled vehicle embeds an event-based cruise control mechanism. In , the reduction of computational load is investigated in [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF] and an event-driven system theory is drawn for disturbed linear systems.

In this paper, using miniature actuators which equipped the bio-inspired robot LORA III [START_REF] Roubieu | Two-directional 1-g visual motion sensor inspired by the fly's eye[END_REF], an experimental setup was developed to test some event-based techniques for speed tracking and disturbance rejection. This work aims at comparing the performances of time-based and event-based control laws in terms of power consumption and control updates in particular an event-based state feedback approach that uses the instantaneous current information to limit the starting peak current and the power consumption. The paper is structured as follows. Firstly, the experimental platform is presented in section II. The event-based angular speed control laws of the miniature thruster are designed in section III. Experimental results are provided in section IV.

II. THE MINIATURE SENSOR-LESS THRUSTER TESTBED

The experimental platform is depicted in Fig. 1A. It is based on a miniature direct current thruster GWS EDF-50 (see Fig. 1A, power: 24.41W, mass: 29g) which was previously embedded on the bio-inspired robot LORA III [START_REF] Roubieu | Two-directional 1-g visual motion sensor inspired by the fly's eye[END_REF]. The maximum allowable current drawn is 3.4A under 7.2V. Some extra components were added to make possible the sensor-less thruster speed regulation.

A. The testbed

The aim of the experimental set-up is to implement different kinds of time-based and event-based control laws from the literature for thruster's speed regulation and disturbance rejection without any angular speed measurement. The performances in term of power consumption and control updates are considered. As presented in the right side of Fig. 1B, a second thruster generates wind disturbances.

B. The data acquisition and control

The data acquisition and control are performed in a realtime framework using a dSPACE electronic board which is connected to the experimental platform. The control laws are implemented using Matlab/Simulink. The measurements and performances are monitored online and one can analyze them off line.

C. Estimation of the rotor speed

In this section, we propose to estimate the rotor speed. It simplifies mechanical assemblies and reduces the size of electronic circuits. Moreover, one can significantly reduce the sampling frequency F s = 20kHz required if a reflective sensor would measure the thruster's angular speed. Therefore estimating the thruster's angular speed allows to save energy and computational resources. The speed estimation simply requires that the armature voltage and the motor current be measured. The estimation principle is presented in [START_REF] Viollet | A 1-gram dual sensorless speed governor for micro-air vehicles[END_REF]. One can recall the main steps.

The dynamics of a DC motor is described by the following equations:

di(t) dt = - R L .i(t) - K e .M L .ω h + 1 L .u(t) (1)
dω h (t) dt = K t J.M .i(t) - D J.M 2 .ω h (t) - 1 J.M .c 0 (t) (2)
• ω h (t) : the thruster angular speed (rad/s) • c 0 (t): an additive bounded disturbance • u(t) : the armature voltage (V)

• i(t) : the armature current (A) The parameters are defined as:

• D : the drag coefficient of the thruster (N.m.s/rad) • J : the inertia of the motor (kg.m²) • K t : the torque constant (N.m/A)

• K e : the back-EMF voltage constant (V.s/rad) • L : the inductance (H) • M : the reduction ratio between the thruster and the rotor

• R : the armature resistance (Ω)

In the following, we consider that the mechanical time constant is greater than the electrical time constant. Therefore, we assume that the armature current i(t) is constant during the thruster's speed ω h (t) transient phase. As shown in Fig. 2A, the measurements of the armature voltage u(t) = v bat -v 1 and the armature current i(t) = v1-v2 Rs measured thanks to a current shunt resistor R s , gives an estimated thruster's speed ω est (t) = ωh (t) in Volt such that : ω est (t) = u(t) -R s i(t). In Fig. 2B, the measured current i(t) is noisy and it is not possible to use it for estimation. The current measurement is filtered using a Kalman filter whose the process noise covariance is set 0.0034 and the measurement noise covariance is set to 20 after experiments. The measured current and the estimated current using the Kalman filter are shown in Fig. 2B.

In the Laplace domain, one writes: Ω est (s) = U (s) -R s I(s) and from (1) in steady state, U (s) = K e M Ω h (s) + RI(s). From (2), assuming D is small and c 0 (t) = 0, in the Laplace domain

I(s) = sJM
Kt Ω h (s). Finally, the estimated velocity is written such as:

Ω est (s) = K e M Ω h (s) + RI(s) -R s I(s) Ω est (s) = K e M 1 + sJ(R -R s ) K e K t Ω h (s) (3)
We define τ = J(R-Rs) KeKt . In practice, one can remark that τ ≪ 1, therefore Ω est (s) ≃ K e M Ω h (s). In Fig. 3, the estimated speed in dashed red line is compared to the measured speed in solid blue line using a reflective sensor at the sampling frequency F s = 20kHz. One can remark that the speed estimation is fair. The gain K e M from (3) is experimentally tuned and fixed for the following work. Finally, we estimate the thruster's speed without any velocity sensor using only the following measurements u(t), v 1 (t) and v 2 (t). The estimated speed will feed the controllers for regulation and disturbance rejection at a lower sampling frequency F s = 5kHz.

D. Discrete-time model identification

Under the assumption of homogenous magnetic field, the miniature direct current thruster is supposed to be a linear plant. Equations ( 1) and ( 2) are considered. The variable to be controlled is the angular speed ω h (t) which is estimated as detailed in sec. II-C. In addition, state-space oriented control is investigated in order to limit the starting peak current. Therefore, we define as state variables the armature current i(t) and the thruster speed ω h (t) such as:

x(k + 1) = A.x(k) + B.u(k) y(k) = C.x(k) (4)
where

x(k) = [i(k), ω h (k)] T , u(k) is the control input signal, y(k) is the measurement, A is a 2 × 2 matrix, B is a 2 × 1 matrix and C = 1 0 0 1 .
For the discrete time state-space model identification, data are recorded at the sampling frequency F s = 5kHz. The measured current varies from 0.2A to 0.5A and the estimated thruster's speed varies from approximately less than 900 rad/s to 1100 rad/s. The identified state-space model for the miniature thruster is given by the representation:

               x(k + 1) = 0.9873 0 -4.262 × 10 -5 0.9541 .x(k) + 49.18 0.2022 .u(k) y(k) = 1 0 0 1 .x(k) (5)
For the state space model identification, the best fits for the speed and the current are respectively equal to 87.7% and 62.2%.

The model validation is performed using a new set of data. Figures 4A and4B exhibit the results. The fits for the speed and the current are respectively equal to 88.7% and 61.6%.

III. EVENT-BASED SPEED CONTROL

The thruster's speed is usually controlled with a timebased proportional-integral (PI) strategy or time-based state feedback. In this section, we propose to design some eventbased approaches introduced in [START_REF] Durand | Further results on event-based pid controller[END_REF] and [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF] for both regulation and disturbance rejection and to compare the results to the ones obtained with the time-based controllers.

A. Time-based PI control

The proportional-integral controller in frequency domain is:

U (s) = K p E(s) + 1 T i s E(s) (6)
where U (.) is the control signal and E(.) the error between the estimated output of the controller system and a given setpoint. K p and T i are some tunable parameters. A discrete time controller is finally obtained, that is:

   u P = K p .(ω ref (k) -ω h (k)) u I (k) = u I (k -1) + Kp Ti .T s .(ω ref (k) -ω h (k)) u(k) = u p (k) + u I (k) (7)
where T s = 2.10 -4 s (ie: F s = 5kHz). In order to limit the amplitude of the control input and the peaks of current, we choose the same settling time at 95% for the closed-loop and the open-loop systems i.e. τ settling = 0.173s. Then T i = 0.0442s and K p = 2.9377.10 -4 are the parameters of the PI controller which gives a gain margin of 64.9dB and a phase margin of 90°.

B. Time-based state feedback control

A discrete-time state feedback controller with integral action for trajectory tracking is also designed such that:

u(k) = -K i K ω h K z . i(k) ω h (k) z(k)
T . Using a backward discretization, the discrete-time control signal is computed as follows:

   u P (k) = -[K i K ω h ].[i(k) ω h (k)] T u I (k) = u I (k -1) -K z .T s .(ω ref (k) -ω h (k)) u(k) = u P (k) + u I (k) (8)
Where K i = -0.2270, K ω h = -0.0001 and K z = -0.0015. The settling time at 95% is the same as for the PI controller.

C. Event-based state feedback control

The event-based state feedback controller is designed using [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF]. A geometric interpretation is proposed in Fig. 5. The control update u(k) is computed at a constant sampling rate T s when the state of the system lies outside a pre-defined set β. When the state is inside the set β, the control keeps its previous value. If we consider the system x(k + 1) = A.x(k) + B.u(k), we design the set β where the control is held by β := {xǫR n | |x| < e T }. The threshold e T ensures a certain level of performance. The discrete-time state feedback control signal is computed as follows: Let introduce the instant τ k . τ k denotes the event times by using a zero-order hold such that: u(t) = u(k) for all t ∈ [τ k , τ k+1 ). Therefore, x(k) = x(τ k ) and u(k) = u(τ k ). In this paper, the controller checks whether the state is inside or outside the set β at a constant rate. This rate is chosen to be equal to the sampling rate of the controller T s such as τ k+1 = τ k + T s . Then as long as the state is outside the set β, the control is updated at a constant sampling time. When the state reaches the set β the control is held to the value just before entering the set β. The gains K i , K ω h and K z are the same as for the time-based state feedback controller.

u(k) = -Kx(k) if |x(k)| e T u(k -1) if |x(k)| < e T (9)

D. Event-based PI control without safety limit condition

The event-based control described in [START_REF] Durand | Further results on event-based pid controller[END_REF] is designed. The event occurs when |e(k) -e(k -1)| > e lim where e is the error e(k) = ω ref (k) -ω h (k) and e lim is a pre-defined threshold. As proposed in [START_REF] Durand | Further results on event-based pid controller[END_REF], we implement an algorithm without the safety limit condition firstly introduced in [START_REF] Årzén | A simple event-based pid controller[END_REF] in order to prevent the output from high overshoot. The hybrid algorithm is a mix between the approach with an exponential forgetting factor and the approach with a saturation of the product h act .e. The resulting algorithm is: The power consumptions are slightly identical (0.45W). In Fig. 7C, the duty cycle of the PWM is only updated during the changes of reference and keeps its previous value during the steady-state. Figure 7D

B. Disturbance rejection

In this section the disturbances are gusts of wind generated by a second thruster. The latter is controlled in open-loop with different duty cycles and activated at different moments. The disturbances occurs at each interval t ∈ [3s; 4s], t ∈ [5s; 6s] and t ∈ [7s; 8s] with a duty cycle respectively equals to 0.5, 0.6 and 0.6 otherwise the motor is turned off with a duty cycle equal to 0. The speed of the main thruster is regulated at 1100 rad/s. Figure 9A shows that the eventbased PI controller rejects the disturbances despite larger overshoots than the ones obtained with the time-based PI controller. The control updates are displayed in Fig. 9B, the reduction of the control updates is about 92%. The power consumptions are the same with the time-based and the event-based controllers (0.725W). The comparison of the control input signals is presented in Fig. 9C and Fig. 9D. From Fig. 10A, the event-based state feedback controller rejects the disturbances. The control updates are reduced up to 96% (see Fig. 10B). The duty cycle is only updated when the disturbances are added or removed (see Fig. 10C). Figure 10D plots the currents. The power consumptions are slightly the same with the time-based controller (0.70W) and with the event-based controller (0.72W). Figure 11 shows the trajectories of the states under the disturbance rejection with the time-based and event-based state feedback controllers. As soon as the states are outside Set(A), the control is updated such that the states reach again the set. Once the control is held, the states converge inside the set. One can say that this strategy succeeds in limiting the power consumption and starting peak current by using the instantaneous current information.

V. CONCLUSION An experimental platform was developed and aimed at testing some event-based control schemes using a miniature thruster. The event-based control laws performed a reference tracking and a disturbance rejection. The thruster's speed estimation was implemented thanks to the measurements of only three voltages. Therefore, the speed estimation was fed into the designed controllers. The event-based approach was performed with a PI controller and a state feedback controller. The event-based state feedback approach was based on a set defined by the controlled output and a inner state, the current, which is somehow an image of the instantaneous energy consumption. The event-based approaches allowed to significantly reduce the number of control updates while having a similar energy consumption compared to a timebased approaches while applying output disturbances.

Fig. 1 .

 1 Fig. 1. A: Picture of the experimental platform which is equipped with a miniature thruster used to compare the event-based algorithms and a additional thruster which generates wind disturbance B: Sketch of the sensor-less miniature thruster testbed using a dSPACE electronic board

Fig. 2 .

 2 Fig. 2. A: Electric diagram of the DC thruster. The current i(t) is obtained with the measured voltages v 1 and v 2 . The voltages are filtered with a lowpass filter. The voltage U P W M controls the thruster speed through the transistor B: The current i(t) is filtered with a Kalman filter. The top plot is the duty-cycle applied in open-loop and the bottom plot shows the measured current in blue line and the filtered current in red dashed line.

Fig. 3 .

 3 Fig. 3. Comparison between the measured speed and the estimated speed: the top plot is the duty-cycle applied in open-loop and the bottom plot shows the measured speed in blue line and the estimated speed in red dashed line.

Fig. 4 .

 4 Fig. 4. Comparison between the measured outputs and the simulated model outputs with the same input. A: Comparison between the estimated current in dashed blue line and the simulated current in red solid line. B: Comparison between the estimated speed in dashed blue line and the simulated speed in red solid line.

Fig. 5 .

 5 Fig. 5. Geometric interpretatrion of the event-based state feedback control proposed in [8]. As long as the state x(k) is inside the set β, the control input signal is u(k) = u(k -1).

Fig. 6 .

 6 Fig. 6. A: Trajectory tracking: comparison between the time-based PI in dashed blue line and the event-based PI from [6] in red line. The figure shows the evolution of the thruster speed with both types of control. B: The event-based control update is also displayed. The update is performed when the indicator is at 1. C: The duty cycles of the PWM in % with the time-based PI in dashed blue line and the event-based PI in red line. D: The current with the time-based PI in dashed blue line and the event-based PI in red line.

Fig. 7 .

 7 Fig. 7. A: Trajectory tracking: comparison between the time-based state feedback in dashed blue line and the event-based state feedback from [8] in red line. B: The event-based control update is also displayed. The update is performed when the indicator is at 1. C: The duty cycles with the time-based state feedback in dashed blue line and the event-based state feedback in red solid line. D: The current with the time-based state feedback in dashed blue line and the event-based state feedback in red solid line.

  displays the currents. The trajectories of the state ω h versus i with the time-based and event-based controllers are plotted in Fig. 8. One can see that the trajectory of the states converges to the different performance sets Set(A), Set(B) and Set(C).

Fig. 8 .

 8 Fig. 8. The figure shows the phase portrait in closed loop: the current i(t) and the thruster speed ω h (t). The trajectories of the time-based controller and the event-based one are respectively plotted in dashed blue line and in red solid line. The performance sets are represented by the grey rectangles.

Fig. 9 .

 9 Fig. 9. A: Disturbance rejection: comparison between the time-based PI in dashed blue line and the event-based PI from [6] in red line. B: The event-based control update is also displayed. The update is performed when the indicator is at 1. C: The duty cycles of the PWM in % with the time-based PI in dashed blue line and the event-based PI in red solid line. D: The current with the time-based PI in dashed blue line and the event-based PI in red solid line.

Fig. 11 .

 11 Fig. 11. The figure shows the phase portrait in closed loop: the current i(t) and the thruster speed ω h (t). The trajectories of the time-based controller and the event-based one are respectively plotted in dashed blue line and in red solid line. The performance set called Set(A) is represented by the grey rectangle.

ACKNOWLEDGMENTS

We thank M. Boyron for his assistance with the electrical design and J. Diperi for his involvement in the mechanical design of the testbed. This work is part of the project PEPS SpikeRob, supported by The French National Center for Scientific Research (CNRS).

Algorithm 1 Hybrid algorithm from [START_REF] Durand | Further results on event-based pid controller[END_REF] % Inputs ysp = u(1);% SetPoint y = u( 2

IV. EXPERIMENTAL RESULTS

In this section, we propose to test the different eventbased control strategies and to compare their performances to the time-based approaches on the miniature thruster for trajectory tracking and disturbance rejection.

A. Trajectory tracking

In Fig. 6A, the thruster's speed with the event-based PI controller from Sec. III-D is plotted and is compared to the time-based PI controller. The event-based approach performs a trajectory tracking without static error. One can remark increases of the overshoots during the changes of reference with the event-based approach. In Fig. 6A, the trajectory tracking is performed with only 1778 updates which is a reduction of 96.5% of the control updates (see Fig. 6B). The power consumption is smaller with the event-based controller (0.52W) than with the time-based one (0.54W). Figures 6C and6D show the comparison between the duty cycles of the PWM in % and the currents. In Fig. 7A, the thruster's speed with the event-based state feedback controller from Sec. III-C is plotted and is compared to the time-based state feedback controller.

We define the performance set β such that:

) with ω tolerance = 40rad/s, i set = 0.3A and i tolerance = 0.1A. Figure 7A shows that the event-based state feedback performs the reference tracking with a static error lower than ω tolerance = 40rad/s.

The reduction of the control updates is about 91% with the event-based state feedback (see Fig. 7B).


Mohamed Boumediene

Christophe Cudel
email: christophe.cudel@uha.fr

Michel Basset
email: michel.basset@uha.fr

Mohammed Boumediene
email: mohammed.boumediene@uha.fr

Abdelaziz Ouamri
email: abdelaziz.ouamri@univ-usto.dz

M Boumediene

Triangular traffic signs

Keywords: Traffic sign detection, Advanced driver assistance systems, Computer vision

come

Introduction

Traffic signs allow the regulation and control of the traffic flow. These signs are important for road users-drivers or pedestrians. Although traffic signs are clearly visible along the roads, the driver can be distracted and miss signs, which can cause dangerous situations. Therefore, a Traffic Sign Recognition (TSR) system assists the driver and gives him information about potential dangerous situations.

Two main approaches to TSR systems are proposed in the literature: the GPS (Global Positioning System) based approach and the vision based approach. The first one recognizes the current traffic signs by extracting information from the GPS navigation data. The disadvantage of this approach is that the GPS map is neither complete nor up-to-date and there are numerous situations where the GPS signal is not available, e.g, in a tunnel. Therefore, the vision based approach is necessary. This second approach uses an optical sensor (vehicle-mounted camera) and locates and recognizes the traffic signs appearing in the captured frames. Recently, the MIPS (Modélisation, Intelligence, Processus et Systémes) laboratory has tried to merge both systems in order to design more robust TSR systems [START_REF] Daniel | Conflict management in multi-sensor dempster-shafer fusion for speed limit determination[END_REF].

This paper focuses on a computer vision approach. Two distinct steps usually compose TSR system by vision, which are detection and recognition. First, Regions Of Interest (ROI) containing the candidate traffic signs are detected in the current frame. This detection is based on the sign characteristics, such as colour and shape, which help to distinguish traffic signs from other objects appearing in the road environment. Then, recognition or classification of the ROIs consists in identifying sign pictograms figuring in each ROI.

Each traffic sign is characterized by its colour, shape and pictogram. In the literature numerous methods use colour information to extract the ROIs containing road signs. There are different colour spaces for segmentation, the one which is most often used is the Hue-Saturation-Value (HSV) space [START_REF] Piccioli | A robust method for road sign detection and recognition[END_REF], [START_REF] Paclka | Road Sign Classification using the Laplace Kernel Classifier[END_REF] and [START_REF] Liu | Scale and skew-invariant road sign recognition[END_REF]. In [START_REF] Arlicot | Circular Road Sign Extraction from Street Level Images using Colour, Shape and Texture Database MAPS[END_REF], adapted thresholds for saturation and hue are defined for more accurate segmentation. A comparison between various segmentation techniques using different colour spaces for traffic sign detection is presented in [START_REF] Gómez-Moreno | Goal Evaluation of Segmentation Algorithms for Traffic Sign Recognition[END_REF]. The results show that normalization improves performance and represents a low cost operation.

Advanced Driver Assistance Systems (ADAS) can contain up to three different applications: TSR, pedestrian and lane detection. Usually, pedestrian and lane detection use only a grey-scale camera and only TSR can use colour information. But colour information is sometimes perturbed by traffic sign surface reflectance which depends on the weather, daylight conditions and the age of signs. This is one of the reasons why a greyscale method is proposed here, which offers the same performance as colour methods.

The aim is to detect triangular traffic signs in a context of grey-scale images, using the original operator RSLD which exploits the symmetry of a triangle. This paper is organized in three sections. Section 2 presents a state of the art of related studies. Section 3 details the proposed traffic sign detection and the results are presented in Sec. 4, followed by a conclusion and discussion on future research work.

Related Studies

The methods proposed in the literature can be classified into two categories: the Hough-like approaches and the template matching approaches. The Hough transform is a popular tool in computing vision, e.g, for line detection and circle detection. For the detection of traffic signs, a voting procedure is carried out in a parameter space, from which the ROIs are obtained by computing the local-maxima in the parameter space. For the template matching approach, a set of shape sign templates is used for detection. Matching can be performed through correlation or by a classifier, such as the Support Vector Machine (SVM). The next paragraphs describe recent algorithms, as well as colour or grey-scale images.

[13] uses the classic Hough algorithm for searching straight lines in order to detect triangular signs. The purpose is to detect three straight lines intersecting each other, forming a 60 degree-angle. The advantage of the method proposed in [START_REF] Garcia-Garrido | Fast Traffic Sign Detection and Recognition Under Changing Lighting Conditions[END_REF] is that it can work day and night with no change of the algorithm. This is possible because the edge points used by the Hough transform are obtained with an extended Canny filter. The two Canny-thresholds are adapted by the image histogram.

In [START_REF] Cyganek | Real-Time Detection of the Triangular and Rectangular Shape Road Signs[END_REF], a colour segmentation is performed to obtain a binary image with selected regions of the candi-Fig. 1 The voting scheme of VBT and STUVT. a Scheme for the voting method of VBT. b Scheme for the voting method of STUVT [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF]. c example of V accu and B accu [START_REF] Belaroussi | Angle Vertex and Bisector Geometric Model for Triangular Road Sign Detection[END_REF] date signs. A simple threshold method helps to perform a segmentation in the Hue-Saturation-Intensity (HSI) space. The salient points related to the shape of the traffic sign are detected in the binary image by a Detector of Local Binary Features (DLBF) [START_REF] Cyganek | Real-Time Detection of the Triangular and Rectangular Shape Road Signs[END_REF]. The set of the salient points are clustered for relating single points, one at each corner. The author does not well explain how the corners are connected for the recognition of triangular shapes.

The symmetric nature of triangular shapes is used in [START_REF] Loy | Fast Shape-based Road Sign Detection for a Driver Assistance System[END_REF]. An extended fast radial symmetry transform is proposed, the Regular Polygon Detector (RPD), which computes the possible shape centroid locations into the grey-scale image. For regular polygons with n corners, the radius is defined as the perpendicular distance from an edge to the centroid. As for detecting the circle shape, each edge point votes for a potential radius along the line of the gradient vector. All votes are accumulated into a matrix and are used with the equiangular information of the regular polygons for computing the centroid of the shape.

The Vertex and Bisector Transformation (VBT) is defined in [START_REF] Belaroussi | Angle Vertex and Bisector Geometric Model for Triangular Road Sign Detection[END_REF]. This transformation is based on the image gradient and produces two arrays: V accu and B accu , accumulating evidence of respectively the angle vertex and the angle bisector. The pair edge points (P i , P j ) with non-parallel orientation gradient vote for the vertex, (A) in Fig. 1a, to which they belong. This pair also votes for the bisector (AB) and helps to obtain two accumulators (V accu and B accu ). The centre of the triangular sign is represented by the local maxima of B accu and its corners are determined by the local maxima of V accu (see Fig. 1c).

Similar to VBT, a Single Target Vote for Upright Triangles (STVUT) is proposed in [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF]. The difference from VBT is that the STUVT computes one accumulator only and the vote relates triple edge points with pairwise non-parallel orientations (see Fig. 1b). To decrease the computational time of the voting process, a heuristic method is proposed to choose the appropriate triple edge points. Pair edge points (P i , P j ) are chosen in the same row image with mirror orientations. The position of the third edge point P k is between P i and P j .

Several preprocessing steps are proposed and evaluated in [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF] to increase the performance of the Houghlike methods. The preprocessing steps transform the colour image into a gradient image which assigns magnitude and orientation to significant pixels, e.g, red pixels which are related to warning traffic signs. Learned Colour Gradient (LCG) and Expected Orientations (EO) provide best results with RPD, VBT and STVUT [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF]. LCG products the gradient image according to local colour distribution. EO is a similar step, but orientation gradient is introduced to enhance triangles detection.

The colour and shape characteristics are used in [START_REF] Gao | Recognition of traffic signs based on their colour and shape features extracted using human vision models[END_REF]. Colour segmentation helps to extract the ROIs containing possible traffic signs from the rest of the image. A quad-tree histogram is applied to segmentation which, uses recursively, a division of the image into quadrants until all elements are homogeneous. To classify the ROIs into a particular group of signs, e.g, triangular signs, the Histograms of Oriented Gradients (HOG) is used as descriptor.

A genetic algorithm is used in [START_REF] De La Escalera | Traffic sign recognition and analysis for intelligent vehicles[END_REF] to search for possible signs in the image. First, an analysis of the hue and saturation components of the image helps to detect the regions which satisfy some colour restrictions. Based on the perimeters of these regions, the genetic algorithm performs a global search to detect the traffic signs where each individual is represented by a sign model with a certain position, a scale and a rotation.

Two steps are proposed in [START_REF] Zaklouta | Warning traffic sign recognition using a HOG-based K-d tree[END_REF] for the detection of warning traffic signs. In the first step, a linear SVM classifier is used, and it is trained with the HOG features of all type of triangular traffic signs. To remove false detections, a blackhat filter [START_REF] Zaklouta | Warning traffic sign recognition using a HOG-based K-d tree[END_REF] is applied to the image by emphasizing the dark pixels with high contrast in their local environment. Therefore, the border and the pictogram of the traffic sign are emphasized and the result is used to eliminate false detections.

[1] proposes a matching method based on an Edge Orientation Histogram (EOH). Using the grey-scale image, a Gaussian pyramidal of three levels is created. At each level, a scanning window is used and the matching measurements, between a set of templates and the region defined by the scanning window are analyzed in order to determine the presence or absence of traffic signs. The matching measurements are based on edge orien-tation histograms defined on different sub-parts which represent a template.

The Local Contour Pattern operator (LCP) is defined in [START_REF] Parada-Loira | Local Contour Patterns for Fast Traffic Sign Detection[END_REF] and [START_REF] Landesa-Vzquez | Fast Real-Time Multiclass Traffic Sign Detection based on Novel Shape and Texture Descriptors[END_REF] as a measurement over contour images to compute the local geometrical structures. For each contour point, obtained with the Canny filter, the local geometrical structure is represented by a LCP code-word. In [START_REF] Parada-Loira | Local Contour Patterns for Fast Traffic Sign Detection[END_REF] and [START_REF] Landesa-Vzquez | Fast Real-Time Multiclass Traffic Sign Detection based on Novel Shape and Texture Descriptors[END_REF], 30 LCP code-words are used to represent 4 sets of orientations. The scanning windows are subdivided into subparts according to geometrical structures defining the triangular signs and circular signs. The 4-bin histogram, related to the 4 sets of orientations, is computed on each subpart and used to detect the traffic signs. In [START_REF] Landesa-Vzquez | Fast Real-Time Multiclass Traffic Sign Detection based on Novel Shape and Texture Descriptors[END_REF], the checking of the candidate signs is based on quantum features which are learned by AdaBoost algorithm. Like Haarlike features, quantum features are a simple intensity level comparison between pixels which are not adjacent.

The method proposed here cannot be classified as a Hough-like approach or template matching approach. It maybe be close to Hough-like approach because it is based on line fitting, but there is no voting process. The main idea is to use the principle of coordinates inversion combined with line detection to detect the symmetry of triangles. The interest of the method is to work on grey-scale images with a final algorithm easily implementable in real-time, and with performance similar to that of the latest latest methods published.

The method proposed

The authors propose a method to detect triangular traffic signs in four steps. This method is based on the Harris corner detector and an original algorithm, RANSAC for Symmetric Lines Detection (RSLD). Figure 2 shows the four steps proposed in the algorithm. The method starts by computing the gradient image following x and y directions (Eq. 1). The first step consists in computing the corners inside the images, using the well-known Harris detector to ensure the extraction of triangular traffic sign corners. The second step uses corner coding in order to select up and down corner candidates related to triangular shapes oriented upwards and downwards respectively. The third step proposes the RSLD method for traffic sign triangular shape detection. The final step checks the presence of triangle baseline to confirm the triangular traffic sign detection.

I x = ∂I ∂x , I y = ∂I ∂y . (1)
Fig. 2 Diagram of the algorithm proposed

Corner detection

For the detection of triangle corners in images, the Harris detector [START_REF] Harris | A combined corner and edge detector[END_REF] is one of the most efficient methods with a low computation time. The Harris detector is based on the covariance matrix M (Eq. 2). This matrix describes the gradient distribution in a local neighbourhood for each pixel. The neighbourhood depends on the Gaussian window W used. The eigenvalues, λ 1 and λ 2 , of this matrix represent the principal gradient components in two orthogonal directions. Both significant eigenvalues indicate the presence of two main directions in the neighbourhood, which can describe a corner at the corresponding position. [START_REF] Harris | A combined corner and edge detector[END_REF] proposes an alternative to eigenvalue computation with the criterion R described in Eq. [START_REF] Belaroussi | Angle Vertex and Bisector Geometric Model for Triangular Road Sign Detection[END_REF] where k is a scalar fixed at 0.04.

I 2 x = I 2 x ⊗ W, I 2 y = I 2 y ⊗ W, I x I y = I x I y ⊗ W, M = [ I 2 x I x I y I x I y I 2 y ]
.

(2)

λ 1 + λ 2 = T r 2 (M ), λ 1 * λ 2 = Det(M ), R = Det(M ) -k * T r 2 (M ). (3)
The set of corners are computed by finding the localmaxima in the R matrix. An example of detection is shown in Fig. 9a.

Corner Coding

The Harris detector is usually employed to compute geometric relations between images or in image registration. These methods associate feature vector descriptors to describe each corner in order to perform a comparison with representative vectors contained in a database [START_REF] Lin | Image registration based on corner detection and affine transformation[END_REF]. One major limitation of descriptors for traffic sign corners is the symmetric neighbourhood considered to create them. For example, the Harris Feature Vector Descriptor (HFVD) proposed in [START_REF] Wang | Harris Feature Vector Descriptor (HFVD)[END_REF] is based on a circular neighbourhood partition. This approach is not adapted for traffic sign detection because the feature descriptor will not only include information of road signs but also of the background scene. This is the reason why an alternative to centered descriptors is proposed: coding the main directions attached to a corner in an appropriate neighbourhood.

The aim is to detect up and down corners related to triangular shapes oriented upwards or downwards respectively, using gradient images (I x , I y ) previously computed to identify both types of corners. A coded image is generated where each pixel is coding according to its orientation. Coded image is obtained with Eq. 4. Intermediate calculations available in the Harris computation are used to define A = (I 2

x > T ), B = (I 2 y > T ), and C = (I x I y < 0). T is a threshold used to select a significant gradient, its value is fixed at 1000. Five classes {0, 1, 2, 3, 4} are defined, in which the non-edge pixels are coded by class 0, and the edge pixels are coded following the four directions depicted in Fig. 3. The diagonal edges are represented by classes 1 and 2. The horizontal and vertical edges are represented by classes 3 and 4 respectively.

Class = A • B • (2 -C) + 3 • A • (1 -B) + 4 • B • (1 -A)(4)
The aim of corner coding is to select corner candidates which represent triangular traffic sign positions. This filtering is based on coded images. For each corner position reported in coded images, an up or down corner is detected according to class distribution. Two class distributions are proposed in adapted subregions. Figure 4a shows subregions for up corner coding. In this case, a subregion concerns the neighbourhood under the corner position where the shape information of road signs must be found. For each corner, if class 2 is dominant in the South-West (SW) subregion and class 1 in the South-East (SE) subregion, then the corner is classified as an up corner. A similar process is applied for down corners in agreement with the subregions shown in Fig. 4b. Figure 5 illustrates the result of the corner coding step applied to a synthetic image, and Fig. 9b on a road scene image. Figure 6 shows the number of remaining triangular traffic sign candidates compared to the Harris corner number. The corner coding is able to select candidate corners only, and to significantly reduce the number of corners to be analyzed in the next steps, which should confirm that a corner represents a triangular traffic sign or not.

RSLD for triangle identification

Principle The methods described in [START_REF] Loy | Fast Shape-based Road Sign Detection for a Driver Assistance System[END_REF], [START_REF] Liu | An MRF Model Based Algorithm of Triangular Shape Object Detection in Color Images[END_REF], [START_REF] Maldonado-Bascón | Road-Sign Detection and Recognition Based on Support Vector Machines[END_REF], and [START_REF] Cardarelli | Road signs shapes detection based on Sobel phase analysis[END_REF] detect the borders of traffic signs with a line detector and analyze the relative positions of these lines with different criteria such as symmetry, the distance to the border or the position of the intersection points. To reduce the computation time for line detection, the RSLD method based on RANSAC algorithm has been designed for this application.

RANSAC algorithm is a well-known, powerful method to fit a model in a set of N points containing data matching with the model (inliers) and out matching data (outliers). It is an iterative algorithm, composed of two mains steps. The first step consists in generating The idea is to straighten the two diagonal sides of triangular shapes. This helps to detect only one segment instead of two with a RANSAC algorithm. This principle is illustrated in Figs. 7a and7b. For up corner candidates, a simple inversion of the y-coordinate applied on coded pixels of class 2 transforms the two diagonal segments into a single segment. For down corner candidates, the same inversion is performed with coded pixels of class 1. This approach converts the symmetry detection into a simple line detection while preserving the symmetry criteria. Indeed, if the corner candidate is really positioned at the top or bottom of a triangular shape, this method will produce a long segment centered around this corner. If the corner is not a true candidate, the segment will be too short and will not detected.

Interest of the proposal There are two advantages which justify this method: accuracy is better, depending on the outliers, and the number of computations is reduced compared to other methods.

Concerning the accuracy, this approach was compared with the separate detection of the two segments, using RANSAC to retrieve one or both of the segments. This comparison was performed by adding outliers to a synthetic image containing two segments representing the diagonal sides of a triangle. For each measurement, the RMSE (root-error-mean-square) geometric distance was computed between the segment(s) model and the segment(s) detected according to the level of the outliers added to the synthetic image. The comparison is performed between two approaches:

without the y-coordinate inversion: detection of two segments using the RANSAC method for each segment. The adapted subregions were defined around the corner candidate (see Fig. 7) and RANSAC was applied to each region separately. -with the y-coordinate inversion: detection of a single segment with the RANSAC method applied to appropriate subregions after the y-coordinate inversion process. Figure 8 shows an example of corresponding synthetic images with and without the y-coordinate inversion method. Inliers were generated from equations: tan( π 3 ) x + y = 0 and -tan( π 3 ) x + y = 0 with Gaussian noise (σ = 3) added to their coordinates. This noise represents perturbations, which generates realistic simulations. The measurements were performed according to outliers percentage added. For a given outlier percentage, the experiment was repeated 1000 times to obtain representative results.

The results of these experiments are reported in Fig. 10. The possibility of easier detection of a long single segment instead of two short, ones is thus con- firmed. The curve presenting the lower RMSE is the one where the y-coordinate inversion was performed.

Concerning the cost of computation, the number of iterations with and without the coordinate inversion was compared. The number of RANSAC iterations, K (Eq. 5), depends on the percentage of the inliers w, complexity of the model n and probability p that RANSAC gives a suitable model [START_REF] Fischler | Random Sample Consensus: A Paradigm for Model Fitting with Apphcatlons to Image Analysis and Automated Cartography[END_REF].

K = log(1 -p) log(1 -w n ) (5)
The RANSAC computation time T can be estimated using Eq. 6 [START_REF] Choi | Performance Evaluation of RANSAC Family[END_REF] where t g and t e represent, respectively, the calculation time for the two RANSAC steps: generation and evaluation.

T = K (t g + N t e ) (6)
The detection of diagonal triangle sides need the RANSAC algorithm to be applied two times in the subregions described in Fig. 7a. First in subregion C and secondly in D. Here, A and B are unused subregions. C and D are considered to contain N 2 samples each. The global computation cost T CD of this detection is then defined as the sum of the computation time for C and D subregions:

T CD = T C + T D , = 2 K (t g + N 2 t e ), = K t g + T. (7)
The RANSAC with y-coordinate inversion used here, is applied to subregions A and D and the number of points in the dataset is then N . The duration T AD of the method can be obtained by adding the y-coordinate inversion time t i to T .

T AD = t i + T (8)
As the y-coordinate inversion is performed once, and not at each iteration, it can be concluded that t i < K t g . In this case T AD < T CD , which shows that the proposed method is faster than two successive RANSAC segment detections. Figure 11 compares the computation cost, with Matlab scripts, for both methods according to the outlier percentage. These curves confirm that the method needs low computation time. In practice, for real images, the computation cost is reduced by 40%.

Multi-scale Detection

The traffic sign scale in images depends on its distance from the camera. To use the present RSLD method efficiently, the region around the candidate corner where the inversion is applied must be defined precisely. The region size must correspond to the road sign's dimensions and without depth information, it is impossible to estimate this size. A multi-scale detection process helps to solve this problem. First, a small region r1 is centered on the selected corner where the corresponding points for RANSAC algorithm are extracted, as shown in Fig. 12. If no line is detected, there is no triangle related to this candidate corner. Otherwise, there remains the length of the line detected in r1. This process is repeated in regions r2 to r5 until the length line L is stabilized, as shown in Fig. 12. This multi-scale detection helps to estimate the resolution of potential triangles detected and gives the diagonal side length which is L 2 . At this stage, only diagonal triangle sides are detected and it is necessary to confirm that this shape is a triangular traffic sign with a baseline detection. A typical example is given in Fig. 9c.

Baseline Detection

In road scene images, numerous objects present shapes which are similar to triangular traffic signs. With two diagonal sides and a symmetrical appearance, they can be detected thanks to the previous step (see Fig. 9c). These objects are then considered as false positive detections and increase the computation time during the traffic sign recognition phase. The aim is to eliminate such false detections in a final step which consists in triangular shape baseline detection. This step exploits the length L 2 obtained with multi-scale detection. First, a horizontal subregion at L 2 sin( π 3 ) distance from the corner position is defined, as illustrated in Fig. 13. Secondly, as the sides of an equilateral triangle are equal, this subregion contains a baseline if the number of pixels from class 3 is L 2 ± ϵ.

Evaluation & Discussion

This section, presents the evaluation of the method proposed by comparing its performance with that of methods presented in recent publications [START_REF] Loy | Fast Shape-based Road Sign Detection for a Driver Assistance System[END_REF], [START_REF] Belaroussi | Angle Vertex and Bisector Geometric Model for Triangular Road Sign Detection[END_REF] and [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF].

The comparisons are performed with two public image databases and on our experimental images.

A first comparison concerns results [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF] published in 2011. The author of [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF] uses the public available Stereopolis database [START_REF]Stereopolis database[END_REF] to compute detection rates of the RPD, VBT and STVUT methods combined with the LCG and EO preprocessing steps. The database is constituted of 847 images from Paris, and there are 27 triangular traffic signs. In table 1, the detection rates of RPD, VBT, and STVUT algorithms are compared with the proposed method. The results show that the detection rate obtained is equivalent to that obtained with STVUT, the best algorithm presented in [START_REF] Houben | A single target voting scheme for traffic sign detection[END_REF]. These methods are Hough-like approaches using colour infor- mation, while the proposed method presents the advantage of similar performance on grey-scale images. The detection rate may even be improved with the introduction of colour information.

Secondly, the authors of [START_REF] Belaroussi | Angle Vertex and Bisector Geometric Model for Triangular Road Sign Detection[END_REF] propose to evaluate the robustness of the VBT algorithm [START_REF] Belaroussi | Détection des panneaux de signalisation routière par accumulation bivariée[END_REF] in the presence of Gaussian noise. This evaluation has been presented in Grigorescu's database of traffic sign [START_REF] Grigorescu | Distance Sets for Shape Filters and Shape Recognition[END_REF] which contains 48 traffic scene images with 40 triangular signs. The experimental conditions described in [START_REF] Belaroussi | Détection des panneaux de signalisation routière par accumulation bivariée[END_REF] have been reproduced. The VBT algorithm and the proposed method are compared by adding different levels of Gaussian noise to images. Figure 15 shows the detection rate for σ varying from 0 to 40. The proposed method provides the best detection rate, for example for σ = 5; the proposed algorithm scores at 100% while the VBT algorithm at 90% only. With σ = 15, the detection rate stays around 90% (see Fig. 16). So, the present approach is more robust to noise than the VBT algorithm because the RSLD algorithm is based on RANSAC which is robust in a noisy environment and the coded image which helps to only select the pixels related to the diagonal sides of a triangular shape.

Thirdly, robustness to rotation must be evaluated. The idea is to reproduce rotation that can appear in images when the car takes a tight turn. As the proposed algorithm is based on the pixel classification according to the orientation gradient, the rotation angle of the sign can decrease the correct detection rate. Figure 17 illustrates the detection in traffic scenes taken from the Grigorescu's database after adding an angular rotation. The conclusion is that the proposed method keeps the same detection rate if the rotation remains under 10 degrees.

Finally, the following results have been obtained with images of the experimental system proposed. This system consists of an Evaluation Camera Kit (ECK1000, Cypress Semiconductor Corporation), designed for automotive applications. The algorithm is implemented in an embedded computer with the Real Time Multisensor Advanced Prototyping Software (RTMaps) [START_REF]Itempora, Real Time Multisensor Advanced Prototyping Software[END_REF], coupled with the Open Computer Vision (OpenCV) [START_REF][END_REF] li- brary. The present method runs in real-time at the camera video rate (30 frames per second). Different traffic situations, have been tests in urban, rural and highway scenes. More than 14,000 frames sizing 640x480 pixels, have been used for the evaluation. The road sections contain 10 triangular warning signs and 13 yield signs. Detection is considered correct if a sign is detected in more than three consecutive frames. The results obtained show a detection rate of 91% (see Fig. 18). The undetected traffic signs are those located on highway exits; they are far from the vehicle and thus not clearly visible. As they are not on the main driving road, their non-detection does not affect the performance of ADAS.

Conclusion

This paper, has presented a novel approach to triangular traffic sign detection. The idea is to detect the corners and the symmetric sides of a triangle. The corners are detected using a well-known Harris detector followed by a corner coding process. This process identifies the up and down corners related to the warning and yield traffic signs. The aim is to detect the symmetric and diagonal sides of a triangle with an original algorithm named RANSAC for Symmetric Lines Detection (RSLD). This algorithm is based on the y-coordinate inversion and transforms a symmetric detection problem into a simple line detection. The multi-scale approach developed helps to locate an adapted subregion for a baseline triangle detection.

The experiments and evaluation show that the proposed method gives a detection rate equivalent to that of the latest methods published. This method is easy to implement in real-time. A future aim is to combine this triangular traffic sign detection with an eye tracking system developed in the MIPS laboratory [START_REF] Bernet | Autocalibration-based partioning relationship and parallax relation for head-mounted eye trackers[END_REF]. For this kind of application, where videos are available, the classification of traffic signs by type is sufficient and will give interesting statistics for driver evaluation.

Future studies will focus on the RSLD method. Instead of performing a simple y-coordinate inversion,

Fig. 3 2 Fig. 4 Fig. 5

 3245 Fig. 3 The edge angle space divided into 4 bins. Classes 3 and 4 represent the horizontal and the vertical pixel edges respectively. The diagonal pixel edges are represented by classes 1 and 2

Fig. 6

 6 Fig.6Results of corner coding on 32 images of Grigorescu's database of traffic signs[START_REF] Grigorescu | Distance Sets for Shape Filters and Shape Recognition[END_REF]. The top curve represents the number of detected corners with the Harris detector. The bottom curve represents the number of selected corners with the proposed corner coding method

Fig. 7

 7 Fig. 7 The inversion around the identified corners. a represents a vertical y-coordinate inversion applied around the selected corners. b represents a vertical y-coordinate inversion applied around the down corners.

Fig. 8

 8 Fig. 8 Synthetic images. The (o) represent inliers, real data related to the triangle's sides with Gaussian noise, and the ( * ) points represent outliers which are uniformly distributed. a represents two sides of a triangle and b represents these sides after applying a y-coordinate inversion.

Fig. 9

 9 Fig. 9 Example of triangular traffic sign detection. a Results of Harris Detector (step 1). b The corner candidates (step 2). c Results of the RSLD proposed (step 3). d Results of Baseline Detection (step 4).

Fig. 10

 10 Fig. 10 Comparison between RANSAC detection with and without the y-coordinate inversion.

  Fig. 11 Computing time of RANSAC detection with and without y-coordinate inversion.

Fig. 12

 12 Fig. 12 Multi-scale Detection.

Fig. 13

 13 Fig. 13 Baseline Detection. a the horizontal subregion and its histogram for the triangle oriented upwards. b the horizontal subregion and its histogram for the triangle oriented downwards.

Fig. 14 Fig. 15

 1415 Fig. 14 Detection Results: warning and yield traffic signs

Fig. 16

 16 Fig. 16 Detection with an additive noise in the images

Fig. 17

 17 Fig. 17 Detection with rotation distortion. a An image rotated by an angle of +10 degrees (left) and -10 degrees (right). b An image rotated by an angle of +14 degrees (left) and -14 degrees (right)

  Fig. 17 Detection with rotation distortion. a An image rotated by an angle of +10 degrees (left) and -10 degrees (right). b An image rotated by an angle of +14 degrees (left) and -14 degrees (right)

Table 1

 1 Comparison of detection rates on Stereopolis database

	Method	Detection rate (%)
	RPD (LCG)	52
	VBT (LCG)	74
	STVUT (LCG)	78
	RPD (EO)	48
	VBT (EO)	81
	STVUT (EO)	81
	Proposed method	81


Benjamin Mourllion
email: benjamin.mourllion@uha.fr

Abderazik Birouche

A Mourllion

Birouche

Modal Truncation for Linear Hamiltonian Systems: a Physical Energy Approach

Keywords: Modal Analysis, Model Order Reduction, Hamiltonian system, H ∞ -norm

This paper proposes a new criterion to select the eigenmodes to be left out in a modal truncation procedure. This criterion takes account of the dynamics of the system on the one hand and, on the other hand, the input/output matrices of the system. Therefore, this approach is satisfactory from a system and control theory point-of-view. Moreover, this criterion is physically meaningful because its computation is directly linked with the physical energy supplied to each eigenmode.

MCS/CCS/AMS Classification/CR Category numbers: 70S05, 78M34

Notations

• M = M t ≻ 0 (resp. 0) means that M is a symmetric positive-definite matrix (resp. positive-semidefinite).

Introduction

Realistic systems are often complex and characterized by models of a very high order. These models are usually difficult to study in a context of analysis. The interest of model order-reduction is that high order models generate high order controllers/observers leading to numerous complexities and difficulties in synthesis, simulation and implementation. Therefore, it is necessary to reduce the order of these systems. The problem of model order-reduction has become the focus of different research areas (mathematics, mechanics, control theory and computer engineering). Several studies have been devoted to this issue following the class of system considered.

In the control theory domain, state-of-the-art methods for reducing the order of linear models can be found in [START_REF] Antoulas | Approximation of Large-scale Dynamical Systems[END_REF][START_REF] Antoulas | A survey of model reduction methods for large-scale systems[END_REF][START_REF]Model Order Reduction: Theory, Research Aspects and Applications[END_REF]. One important contribution to model reduction is the balanced truncation method [START_REF] Glover | All optimal Hankel-norm approximation of linear multiariable systems and their L ∞ -error bounds[END_REF][START_REF] Moore | Principal component analysis in linear systems: controllability, observability and model reduction[END_REF]. In this approach, each state is equally controllable and observable and the reduced order model is obtained by truncating the least controllable and observable states. Balanced truncation is a common method to reduce the order of Linear Time Invariant (LTI) systems thanks to its simplicity of implementation; the stability of the reduced system is guaranteed and the approximation error is bounded. Another important solution to model reduction problems is the optimization approach which minimizes the H ∞ -norm (or H 2 -norm) between the original full order model and the reduced order model [START_REF] Lam | H ∞ and L ∞ /L 2 model reduction for system input with sector nonlinearities[END_REF][START_REF] Yan | An approximate approach to H 2 optimal model reduction[END_REF]. The main drawback of these solutions is that the reduced models obtained do not necessarily have obvious physical meanings. To overcome this problem, some methods propose several adaptations of the balanced truncation in order to preserve the second-order form of the systems [START_REF] Chahlaoui | Second-order balanced truncation[END_REF][START_REF] Meyer | Balancing and model reduction for second-order form linear systems[END_REF][START_REF] Stykel | Balanced truncation model reduction of second order systems[END_REF]. All these methods cannot be applied to conservative systems where H ∞ -norm is infinite. Indeed, in this case, the reachability Gramian and observability Gramian cannot be computed (since the Hankel singular values are infinite), and so, the truncation is not possible. It can be noticed that van der Schaft [START_REF] Van Der Schaft | Model reduction of linear conservative mechanical systems[END_REF] proposed an alternative approach by associating a gradient (or reciprocal) system with the conservative system [START_REF] Van Der Schaft | Time-reversible Hamiltonian systems[END_REF][START_REF] Willems | Dissipative dynamical systems II: Linear systems with quadratic supply rates[END_REF] but this ad hoc method is not straightforwardly interpretable.

In the mechanical and physical domain, several other methods have been developed, such as Guyan Condensation [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF], Component Mode Synthesis [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF], Modal Truncation, for example. State-of-the-art and comparative studies of these methods can be found in [START_REF] Koutsovasilis | Comparison of model reduction techniques for large mechanical systems[END_REF].

It is well known that modal truncation helps to keep stronger physical meanings of states but the selection of eigenmodes to be removed is not easy to perform and remains an open issue. As underlined by Wortelboer [27, p. 60], all sorts of criteria for selection have been proposed in literature. Besides, selection criteria are often not made explicit. It is difficult to trace the origins of most modal reduction methods. This problem will be the topic of the present paper.

The main contribution of this paper is the definition of a new criterion to select the eigenmodes to be left out in the modal truncation step. This criterion is based of the physical energy supplied by the control command and the passivity [START_REF] Willems | Dissipative dynamical systems I: General theory[END_REF][START_REF] Willems | Dissipative dynamical systems II: Linear systems with quadratic supply rates[END_REF] of physical systems. In order to preserve as much as possible the connection between this new criterion and the physical energy, the Hamiltonian formulation of physical systems will be used.

This paper is organized as follows. Section 2 describes the problems and drawbacks of conventional modal truncation. Section 3, considering classic (i.e. undamped [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF][START_REF] Goldstein | Classical Mechanics[END_REF][START_REF] Landau | Mechanics -Third edition[END_REF]) Hamiltonian systems, presents a strong relation between the energy supplied and the Fourier transform of the input is introduced. This relation leads to the formulation of a new criterion. Two examples of modal energy computation are then presented. The first one deals with a finite-time input signal, while the second one deals with a infinite-time input signal. Section 4 shows how this criterion can be used in the eigenmode selection in modal truncation. One example is studied, with two different modal truncation procedures. Finally, the conclusion will be followed by a few perspectives.

Motivation and Problem Positioning

Introduction

Generally, modal analysis deals with physical second-order form systems given by:

Mq(t) + C q(t) + Kq(t) = bu(t) y(t) = c q(t) with dim(q) = n (1)
where M, C and K are respectively the mass (or inertia) matrix, damping matrix and stiffness matrix. To ensure the physical feasibility of the second order form system equation (1), these three n × n matrices have to verify the so-called structural conditions:

• M = M t ≻ 0 • C = C t 0 • K = K t ≻ 0 B. Mourllion and A. Birouche

Modal basis

As M = M t ≻ 0 ⇒ det(M) = 0, the inertia matrix is invertible and system (1) has n eigenmodes, solution of: The mass and stiffness matrices are symmetric and positive definite, then the modal matrix is diagonal and the eigenfrequencies are real and positive and the eigenmodes are real. Moreover, the eigenmodes are orthogonal

M -1 K Φ = ΦΩ
φ t i Mφ j = 0 and φ t i Kφ j = 0 if ω i = ω j .
Therefore, the modal transformation gives the following diagonal matrices

M = Φ t MΦ K = Φ t KΦ
Due to the damping matrix, orthogonal modes can be coupled through the following quantity:

C = Φ t CΦ.
If the damping matrix can be written as C = a m M + a k K (known as Rayleigh damping condition [START_REF] Rayleigh | Theory of Sound (2 volumes)[END_REF]) or as C = M ∞ j=-∞ a j M -1 K j (known as Caughey damping condition [START_REF] Caughey | Classical normal modes in damped linear dynamic systems[END_REF]), the modal damping matrix can be written as

C = diag{2ξ i ω i }
then the damped eigenmodes remain uncoupled.

Drawbacks of modal truncation

One can rewrite the system presented in equation (1) in order to obtain the following state-space realization:

January 31, 2013 16:58

Dynamical Systems FinalRevisedMBPaper

Modal Truncation for Linear Hamiltonian Systems 5

                   q(t) q(t) = A -M -1 C -M -1 K I 0 q(t) q(t) + B M -1 b 0 u(t) y(t) = c 0 C q(t) q(t)
Assuming that the eigenmodes are uncoupled, the 2n eigenvalues λ i of the matrix A are:

eig(A) = λ 1,2 = -ξ 1 ω 1 ± ω 1 1 -ξ 2 1 , . . . , λ 2n-1,2n = -ξ n ω n ± ω n 1 -ξ 2 n .
Usually, the aim of modal truncation is to reduce the order of a system by leaving out the eigenmodes corresponding to some criteria enumerated below.

• In the damped case, the eigenmodes to be neglected can be selected according to (1) the eigenfrequencies (usually the higher ones), (2) the dynamics (usually the more damped ones), (3) the H ∞ -norm (the smaller ones). Indeed, the transfer function can be written as:

G(s) = c 1 . . . c 2n    s -λ 1 . . . s -λ 2n    -1    b 1 . . . b 2n    = 2n i=1 cibi s-λi .
The H ∞ -norm of each term of the sum is given by cibi s-λi ∞ = |cibi| |Reλi| . So, in order to obtain a small error upper bound in the H ∞ -norm, the terms with the smaller H ∞ -norm have to be removed.

• In the undamped case (the dissipation matrix is zero: C = 0), the poles of the system are located on the imaginary axis, therefore the H ∞ -norm (or L ∞ -norm) and the Hankel singular values cannot be computed. Currently, the only way to perform modal truncation is to leave out the eigenmodes corresponding to some of the eigenfrequencies (usually the higher ones).

As underlined in [START_REF] Van Der Schaft | Model reduction of weakly damped mechanical systems[END_REF], from a system and control theory point-of-view, in the undamped case, the reduced model may have disadvantages since the analysis is only concerned with the system matrix and not with the input and output matrices. In the following Sections, a new criterion is defined in order to take account of both the system dynamics and the input and output matrices. B. Mourllion and A. Birouche

Energy-Based Criterion

The aim of this Section is to define a new criterion based on the physical energy. Dealing with passive systems [START_REF] Willems | Dissipative dynamical systems I: General theory[END_REF][START_REF] Willems | Dissipative dynamical systems II: Linear systems with quadratic supply rates[END_REF] and focusing on the energy supplied, a natural approach is to use the Hamiltonian framework.

Hamiltonian formulation

Generally, a dissipative Hamiltonian system is defined [START_REF] Van Der Schaft | Port-Hamiltonian systems: an introductory survey[END_REF] as

ẋ(t) = (J(x(t)) -R(x(t)))∇ x H(x(t)) + B(x(t))u(t) y(t) = B t (x(t))∇ x H(x(t)) (2)
with:

• x = q t p t t the state vector.

• q = q 1 . . . q n t ∈ M is the vector of the generalized coordinates, belonging to the configuration manifold M with dim(M ) = n.

• p = p 1 . . . p n t ∈ T *
q M is the vector of the generalized momentum, belonging to the cotangent space to M at the configuration q. Therefore, the vector x belongs to T * M the cotangent bundle(also known as phase space) of M and dim(x) = 2n.

• H(x) is the Hamiltonian function and represents the total energy of the system, • ∇H represents the column gradient vector.

• J(x) = -J t (x), the symplectic matrix. It is smooth with respect to x.

• R(x) = R t (x) 0, the dissipation matrix. It is a smooth function with respect to

x. • B is a column vector.

A convenient property of this kind of systems is that the energy balance can be easily obtained:

d dt H(x(t)) = ∇ t x H(x(t)) ẋ(t) = ∇ t x H(x(t))[(J(x(t)) -R(x(t)))∇ x H(x(t)) + B(x(t))u(t)] = ∇ t x H(x(t))J∇ x H(x(t)) =0 -∇ t x H(x(t))R(x(t))∇ x H(x(t)) + ∇ t x H(x(t))B(x(t))u(t) = y t (t)u(t) -∇ t x H(x(t))R(x(t))∇ x H(x(t))
Since R 0, the Hamiltonian systems are passive [START_REF] Van Der Schaft | Port-Hamiltonian systems: an introductory survey[END_REF].

In this paper, the focus will only be on linear invariant Hamiltonian systems. In this case, matrices J, R and B are independent of x. Moreover, the Hamiltonian function is defined by the following quadratic form:

H(x) = 1 2 x t Ex
where E = E t 0 is called the energy matrix.

The system presented in equation ( 2) has the following state-space realization (called energy variable representation [START_REF] Polyuga | Structure preserving model reduction of port-hamiltonian systems[END_REF]):

ẋ(t) = (J -R)Ex(t) + Bu(t) y(t) = B t Ex(t)
The Hamiltonian function is defined over the cotangent bundle T * M . The generalized velocities q = q1 . . . qn t belong to the tangent space to M at the configuration q: T q M . In the linear case, the value of the kinetic energy E k (p) = p t M -1 p is equal to the value of the kinetic co-energy E * k (q) = qt M q (Legendre transform of E k ). More properties of energy and co-energy can be found in [START_REF] Polyuga | Structure preserving model reduction of port-hamiltonian systems[END_REF][START_REF] Woodson | Electromechanical Dynamics-Part I[END_REF] for instance.

Relation Between the Input Fourier Transform and the Energy Supplied

Considering the linear Hamiltonian systems, an important relation between the Fourier Transform of the input system with the energy supplied by this input is presented in this Section. To the author's knowledge, this relation has never been presented and used in the related literature. This is the key point on which the new criterion is based. Consider the undamped harmonic oscillator of the form.

d 2 q dt 2 (t) + ω 2 0 q(t) = 1 m u(t) with ω 2 0 = k m . (3)
Defining p = m q, the Hamiltonian function is H(q, p)

= 1 2m p 2 + 1 2 kq 2 = 1 2 m q2 + 1 2 kq 2 ,
this system can be rewritten as equation ( 2) with J = 0 1 -1 0 , R = 0 and B = 0 1 .

The energy variable representation is immediately obtained:

A B C D =     0 1 m -k 0 0 1 0 1 m 0     B.
H(x(t)) = 1 2m |U (ω 0 )| 2 = CB 2 |U (ω 0 )| 2 ∀t ≥ T ( 4
)
where U is the Fourier transform of u.

Proof :

The oscillator is assumed to be initially in equilibrium with a nil mechanical energy. Therefore, the energy supplied to the system is its mechanical energy at the time T : H T and remains constant after this time:

H(t) = H T ∀t ≥ T .
Defining y(t) = q(t), the energy supplied by the input force u is given by its work:

W = ∞ -∞ y(t)u(t)dt = 1 2π ∞ -∞ Y (ω)U * (ω)dω
where the second equality is obtained by the application of the Parseval theorem and the uppercase variables are the Fourier transforms of the respective lowercase variables.

The derivative property of the Fourier transform leads to Y (ω) = ωQ(ω) where Q (the Fourier transform of q) is obtained from the equation (3):

Q(ω) = 1 m U (ω) ω 2 0 -ω 2
The singular points (ω = ±ω 0 ) lead to a problem in the evaluation of W .

To obtain a better definition of the work, the singular points are slightly shifted out of the path of integration through the introduction of an infinitesimally small friction coefficient ε > 0 (ε 2 negligible compared with ω 2 0 ). The previous Fourier transforms becomes:

Q ε (ω) = 1 m U (ω) ω 2 0 -(ω -ε) 2 and Y ε (ω) = ωQ ε (ω) = ω 1 m U (ω) ω 2 0 -(ω -ε) 2 ,
Moreover, as y(t) and u(t) are real, their Fourier transforms have an Hermitian symmetry (i.e. F (-ω) = F * (ω)), so:

W ε = 1 π Re ∞ 0 Y ε (ω)U * (ω)dω = 1 mπ Re ∞ 0 ω ω 2 0 -(ω -ε) 2 |U (ω)| 2 dω
Now, using the so-called Dirac-Plemelj formula (also known as the Sokhotsky-Plemelj January 31, 2013 16:58

Dynamical Systems FinalRevisedMBPaper

Modal Truncation for Linear Hamiltonian Systems 9 formula, see [START_REF] Gakhov | Boundary Value Problems. Adiwes International Series in Mathematics[END_REF] for instance):

lim η→0 + 1 x ± η = P 1 x ∓ πδ(x) (5)
the following equation is obtained letting ε to zero:

W = 1 m ∞ 0 ω -δ(ω 0 + ω)P 1 ω 0 -ω + δ(ω 0 -ω)P 1 ω 0 + ω |U (ω)| 2 dω
where the first part of the integral is always zero. Therefore, this leads to:

W = 1 2m |U (ω 0 )| 2 (6)
Note that T is the time from which the input is zero. This time can be either finite or infinite.

Several remarks can be made about equation ( 6):

(1) The work W is always positive or nil, this is consistent with the fact that Hamiltonian systems are passive. (2) If the control signal is a sinusoidal u(t) = sin(ω u t), its Fourier transform is U (ω) = δ(ω u -ω). Therefore, if ω u = ω 0 , the supplied energy tends to infinity due to the resonance effect or else the supplied energy is zero. (3) If the control signal is a unitary impulse u(t) = δ(t), its Fourier transform is the unitary constant U (ω) = 1. The supplied energy is therefore W = 1 2m . This is consistent with the shock theory (conservation of the impulse). Indeed, by definition, the impulse conservation is given by: u average ∆t = m∆v = 1 (since u is a unitary impulse). Moreover, the velocity before the impact is assumed to be zero so the following relation v = 1 m is obtained. Just after the impact, the supplied energy is only the kinetic energy:

E k (v) = 1 2 mv 2 = 1 2m = W .

Energy Supplied to Eigenmodes

Consider the following Hamiltonian system: Let Φ be the modal matrix transformation defined in Section 2.2 and assume (with no loss of generality) that Φ t Φ = ΦΦ t = I n . Define Θ = diag{Φ, Φ} and consider the following transformation x = Θx. By underlining that the gradient is a covariant vector (∇ x = Θ∇ Θ t x = Θ∇ x) and that the modal transformation is a canonical transformation (Θ t JΘ = J), the Hamiltonian system in its modal basis can be written as follows:

ẋ(t) = J∇ x H(x(t)) + Bu(t) y(t) = B t ∇ x H(x(t)) (7
ẋ(t) = J∇ x H(x(t)) + Bu(t) y(t) = Bt ∇ x H(x(t)) (8)
with:

x(t) = q p and B = 0 b = Θ t B H(x(t)) = 1 2 xt (t) Ẽ x(t) and Ẽ = Θ t EΘ = K 0 0 M-1 M = diag{ m1 , .
. . , mn } and K = diag{ k1 , . . . , kn }

As eigenmodes are orthogonal, the Hamiltonian function H is the sum of the n Hamiltonian sub-functions Hi of each following Hamiltonian sub-system (written in the energy variable representation):

Ãi Bi Ci Di =     0 1 mi -ki 0 0 bi 0 ci 0    
where bi is the ith entry of the vector b and ci = bi mi . Now, using equation ( 6), ∀t ≥ T , the n Hamiltonian sub-functions Hi can be written as follows: Hi (x i (t)) = 1 2 C i B i |U (ω i )| 2 (with ω 2 i = ki mi and xi the associated modal statevector).

H(x(t))

= n i=1 Hi (x i (t)) = n i=1 1 2 C i B i |U (ω i )| 2 (9)

Example of Modal Energy Computation

In this section, consider the following system composed of 3 masses chained by 3 springs:

m 1 m 2 m 3 k 1 k 2 k 3 u y
The input is a force u applied to the third mass and the output is the velocity y of the same mass. This system can be described by the following matrices:

M = diag(m 1 , m 2 , m 3 ); K =    k 1 + k 2 -k 2 0 -k 2 k 2 + k 3 -k 3 0 -k 3 k 3    ; b =    0 0 1    ; and c = b t
The energy supplied to the 3 modes of the system is computed considering the 2 following cases:

(1) Rectangular input: u(t) = A ∀0 ≤ t ≤ T (u = 0 elsewhere) (2) Exponential input: u(t) = A exp(-t/t0) ∀t ≥ 0 (u = 0 elsewhere)

The numerical results were obtained with: m 1 = m 2 = m 3 = 1 kg and k 1 = k 2 = k 3 = 4 N/m. The computed eigenfrequencies (in rad/s) are: ω 1 = 0.8901, ω 2 = 2.4940 and ω 3 = 3.6039.

• In case 1, the amplitude of the input was set with A = 1 and the length T of the rectangular input was set at the specific value T = 8π/ω 2 ≃ 10.0774 s. Figure 1 shows the outputs of the three modes and the output of the system on the left. On the right, one can observe the energy supplied to the three eigenmodes and to the whole system. The dashed lines refer to the computed values (using equations ( 6) and ( 9)) of the 4 Hamiltonian functions of each eigenmode and the system. Note that the energy related to the second eigenmode H 2 is zero. This is due to the fact that the length of the rectangular input is a multiple of 2π/ω 2 . Figure 2 shows the modulus of the Bode system transfer function in the upper part and in the lower part one can observe the modulus of the input Fourier transform. The 3 squares on the curve show the values at the 3 eigenfrequencies. For the second eigenfrequency, this value is zero which explains why H 2 is equal to zero. • In case 2, the amplitude of the input was set with A = 1 and t 0 was set at the specific value t 0 = 4 s. Figure 3 shows the outputs of the three modes and the output of the system on the left. On the right, one can observe the energy supplied to the three eigenmodes and to the whole system. The dashed lines refer to the computed values (using equations ( 6) and ( 9)) of the 4 Hamiltonian functions of each eigenmode and the system. B. Mourllion and A. Birouche

Truncation Step

Physically, it makes sense to keep the most energetic modes, i.e. the modes which have stored most of the energy. This criterion, based on both the eigenfrequency and input/output matrices, is dependent of the type of input signal. If the type of input signal is previously known, its Fourier transform can be directly used to define the criterion and perform a modal truncation in accordingly. In contrast, if there is no knowledge on the type of input signal, the two following cases can be considered for the definition of the criterion:

(1) Dirac (impulse) case: u(t) = δ(t) ⇋ U (ω) = 1 which leads to

H(x(t)) = n i=1 Hi (x i (t)) = n i=1 1 2 C i B i ∀t > 0 (10)
(2) Heaviside (step) case:

u(t) = Γ(t) ⇋ U (ω) = P 1 ω + 1 2 δ(ω) which leads to H(x(∞)) = n i=1 Hi (x i (∞)) = n i=1 1 2 C i B i 1 ω 2 i (11)
These signals are suggested for two main reasons:

(1) Firstly, these functions are frequently used in system analysis (impulse for the study of the impulse response and step for the study of second order structured system) and reduction. When a balanced truncation is performed on an asymptotic stable system, the reachability gramian W A considered is solution of the Lyapunov equation AW A + W A A t + BB t = 0 but can also be understood as W A = ∞ 0 x δ (τ )x t δ (τ )dτ where x δ is the state impulse response. Indeed,

W A = ∞ 0 t 0 exp(A(t -τ ))Bδ(τ )dτ t 0 exp(A(t -τ ))Bδ(τ )dτ t dt = ∞ 0 exp(At)BB t exp(A t t)dt.
(2) Secondly, when the impulse function is used, the criterion is only based on the input output matrices (independently of the eigenfrequencies since its Fourier transform is a constant). When the primitive of the impulse (namely the step) is used, the Fourier transform can be seen as a penalty function which tends to favorise eigenmodes with small eigenfrequencies. Which can be a trade-off between the proposed approach and the classical approach (remove the higher eigenfrequencies).

Although there is no theoretical link between them, it is worth underlining that there is a strong similarity between the second case (equation [START_REF] Koutsovasilis | Comparison of model reduction techniques for large mechanical systems[END_REF]) and the truncation criterion based on the H ∞ -norm in the damped case. In fact, if the damping matrix is not equal B. Mourllion and A. Birouche to zero, then the modal matrices A i can be written as:

Ãi = 0 1 mi -ki -2ξ i ω i
The H ∞ -norm of the ith Hamiltonian sub-system can be approximated by: Gi ∞ ≃

1 2 C i B i 1
ξiωi . Indeed, in the SISO (single input single output) case, the H ∞ -norm corresponds to the peak gain of the frequency response. This peak occurs at the pseudoeigenfrequency which is close to the eigenfrequency. Therefore, the following approximation can be made: G i ∞ ≃ σ max (G i (ω i )). Using the triangle inequality, the following relation, quite similar to the relation [START_REF] Koutsovasilis | Comparison of model reduction techniques for large mechanical systems[END_REF], is then obtained:

G ∞ ≤ n i=1 Gi ∞ = n i=1 1 2 C i B i 1 ξ i ω i (12)

Example of Reduction

Consider a system of the form (7) with The eigenfrequencies of the system: ω 2 1 = 1, ω 2 2 = 2 and ω 2 3 = 5 (in rad/s) were set arbitrary. The matrices K and b were generated pseudo-randomly.

x = q p B = 0 b H(x) = 1 2 x t K 0 0 M -1 x M = I 3 K =    1
The corresponding state-space representation of this system is composed of 6 states (3 generalized positions and 3 generalized velocities). In order to reduce it by removing one eigenmode, two energy-based modal reductions were performed:

(1) The first one was computed using a unitary impulse function (denoted by the superscript ↑). The energy supplied to the 3 eigenmodes is given in the following As H ↑ 1 > H ↑ 3 > H ↑ 2 , the eigenmode to be removed is the second one (H ↑ 2 = 4.5002 J).

(2) The second one was computed using a Heaviside (step) function with a unitary amplitude (denoted by the superscript ). The energy supplied to the 3 eigenmodes is given in the following table :   Index mode i 1 2 3 Eigenfrequency (rad/s) ω i 1.0000 1.4142 2.2361 Energy (J) H i 12.5003 2.2501 1.6000 As H 1 > H 2 > H 3 , the eigenmode to be removed is the third one (H 3 = 2.2361 J).

The difference between H ↑

1 and H 1 is only due to numerical approximation. Figure 5 shows the Bode diagram of the original system and the 2 reduced systems. In this example, the truncation with respect to a step function gives the same results as the conventional modal truncation, since the Hamiltonian sub-functions are ordered as their corresponding eigenfrequencies. On the other hand, the truncation with respect to an impulse function does not remove the same eigenmodes.

As the H ∞ -norm (or L ∞ -norm) cannot be computed, the responses to an impulse and a step input are plotted in figure 6 in order to compare the reduced models. 16:58 Dynamical Systems FinalRevisedMBPaper

Conclusion

This paper has discussed the limitations of traditional modal truncation and proposed a new method for the selection of the eigenmodes to be left out. This approach is based on an energetic criterion defined in the same way as the ones for damped system (modal truncation based on the H ∞ -norm or balanced truncation based on the Hankel singular values). The new criterions considers the dynamics of the system on the one hand, and, on the other hand the input/output matrices. So, the approach is satisfactory from the point-of-view of a system and control theory. This study focused on conservative systems but the proposed method can be applied on dissipative systems. Further work will propose a trade-off between the energy supplied and the energy dissipated with each eigenmode in order to select the states to be neglected for dissipative systems.

2

 2

  1 , . . . , ω n } the eigenfrequencies matrix, Φ = [φ 1 , . . . , φ n ] the modal matrix.

Figure 1 .

 1 Figure 1. Modal response and modal energy: case 1 (rectangular input).

Figure 2 .Figure 3 .

 23 Figure 2. Bode diagram and input Fourier transform: case 1 (rectangular input).

Figure 4 .

 4 Figure 4. Bode diagram and input Fourier transform: case 2 (exponential input).

Figure 5 .

 5 Figure 5. Bode Diagram obtained for impulse-based reduction and step-based reduction.

  Mourllion and A. Birouche Proposition 3.1 Considering a forcing function (or distribution) u with the support [0, T ] (where T > 0 can be finite or infinite). The energy supplied to an harmonic undamped oscillator given by equation (3) initially at equilibrium with nil mechanical energy is

January 31, 2013 16:58

Acknowledgement

This work was supported by the French national project INOVE/ ANR 2010 BLAN 0308.