Robineau Fabien
email: fabien.robineau@imag.fr

Vuillerme Nicolas

Jean-Pierre Orliaguet

Yohan Payan

Tongue Liminary Threshold Identification to Electrotactile Stimulation

Many applications use electrostimulation of the human skin to provide tactile sensation. The effect of electrotactile stimulations were studied on a 6x6 matrix of tactile electrodes placed on the anterior part of the tongue. The liminary threshold with continuous or discontinuous waveform and patterns with 2 or 4 electrodes was investigated.

The result suggest that for energy saving and to improve the yield, it would probably be better to use discontinuous stimulation with two electrode patterns.

Introduction

Electrotactile stimulation can provide sensory information not available to damaged classical sensory system (vision, hearing, vestibular or somatosensory systems). Many electrotactile human machines interface (HMI) have been developed and applied to various surface areas [START_REF] Bach-Y-Rita | Seeing with the brain[END_REF]. The first development of electrotactile visual substitution systems (TVSS) was designed to provide distal spatial information to blind people via a 20x20 electrodes matrix [START_REF] Bach-Y-Rita | Vision substitution by tactile image projection[END_REF]. Later Bach-y-Rita and his collaborators converged towards the electro-stimulation of the tongue surface, called "Tongue Display Unit" (TDU), which provides a practical HMI [START_REF] Bach-Y-Rita | Form perception with a 49-point electrotactile stimulus array on the tongue[END_REF]. The human tongue is very sensitive, highly mobile and discriminative (spatial threshold < 2 mm). Moreover, the environment of the mouth offers a protected volume.

As concerns the tongue sensitivity, most authors agree that the tactile sensitivity threshold of the anterior part of the tongue, and particularly the tip, is lower than the other part [START_REF] Marlow | General sensory innervation of the human tongue[END_REF][START_REF] Maeyaman | Minimal two-point discrimination in human tongue and palate[END_REF]. Moreover, it was showed that the mechano-receptive innervation is denser on the tip than in the other regions. However, few research has been investigated about the asymmetry of tongue sensitivity with respect to the median sulcus (right and left tongue sides) [START_REF] Ringel | Oral perception. I. Two-point discrimination[END_REF].

TIMC-IMAG laboratory recently developed medical electrotactile human machine interfaces using the TDU of Bach-y-Rita with a 6x6 matrix. These devices, which feasibility was evidenced, are used for disable persons assistance and for computer-aided surgery. Either the TDU provides electrotactile signal to supply punctual information, or the functioning is almost continuous to give orientation information. An application to prevent pressure ulcer formation in paraplegics, developed as an "alarm", provides via the TDU, tactile information about excess of pressure at the skin/seat interface [START_REF] Moreau-Gaudry | A New Health Strategy to Prevent Pressure Ulcer Formation in Paraplegics using Computer and Sensory Substitution via the Tongue[END_REF]. In addition, a device to improve human ankle joint position sense was investigated [START_REF] Vuillerme | Improving human ankle joint position sense using an artificial tongue-placed tactile biofeedback[END_REF]. On the other hand, electrotactile information comparable to continuous information is investigating. In this way, a biofeedback system aims at improving human balance control for people with visual or hearing impairments [START_REF] Vuillerme | Controlling posture using a plantar pressure-based, tongue-placed tactile biofeedback system[END_REF]. Also, is studied a system to provide to a surgeon, via the TDU device, orientation information to accuracy guide a needle until a target inside a body [START_REF] Robineau | Guiding the surgical gesture using an electro-tactile stimulus array on the tongue: A feasibility study[END_REF].

In these studies, a minimum of four electrode patterns were used because the subjects were able to discriminate without problem the stimuli. In parallel, a more ergonomic wireless 6x6 TDU device was recently developed in our laboratory. It is inserted in a dental retainer including microelectronics, antenna and power supply [START_REF] Robineau | Guiding the surgical gesture using an electro-tactile stimulus array on the tongue: A feasibility study[END_REF]. To increase the comfort of such an embedded device, the size and the energy consumption have to be even more reduced. The intensity on the matrix surface has to be as low as possible, but of course sufficient to be perceived by the subject (i.e. higher than the "liminary threshold" of the human tongue). Another option to decrease such intensity would be to reduce the number of activated electrodes. However, such a reduction should not decrease the discrimination of the electrotactile patterns.

The present study aims at evaluating the liminary thresholds of the tongue with continuous or discontinuous electrotactile stimulations among subjects, and comparing such thresholds with 2 or 4 electrode patterns.

Methods

Subjects

Two groups of ten subjects (age: 28,3 ± 3,9 years; body weight: 68,8 ± 10,2 kg; height: 174,6 ± 9,9 cm) voluntarily participated in this experiment. None of the

Proceedings of ENACTIVE/07

4th International Conference on Enactive Interfaces Grenoble, France, November 19 th -22 nd , 2007 subjects presented any history of sensory/motor or cognitive problem.

Apparatus

Electrotactile stimuli were delivered to the dorsum of the tongue via a ribbon TDU derived from the one developed by Bach-y-Rita [START_REF] Bach-Y-Rita | Seeing with the brain[END_REF][START_REF] Bach-Y-Rita | Form perception with a 49-point electrotactile stimulus array on the tongue[END_REF]. This electrotactile device consisted of an array of 36 tactile electrodes (6×6 matrix, radius: 0.7mm each), embedded in a 1,5×1,5 cm plastic strip (Fig. 1, left panel). The tip of the TDU was inserted in the oral cavity and held lightly between the lips, maintaining the array in close and permanent contact with the surface of the tongue (Fig. 1, right panel). A flexible cable, made of a thin (100 µm) strip of polyester material, connected the matrix to an external electronic device. This device delivered the electrical signals that activated the tactile receptors on the anterior superior part of the tongue. As Kaczmarek suggested [START_REF] Kaczmarek | Electrotactile adaptation on the abdomen: Preliminary results[END_REF], the frequency of the DC stimulating pulses was fixed to 50Hz for all trials and subjects. Because of the conductive properties of the saliva and the epidermis thickness, the TDU only required a 5-15V input voltage and a 0.4-4.0mA current. Each goldplated circular electrode (Ø = 1.4mm, inter-centre distance = 2.3mm) received monophasic pulses. Considering that each electrode could be independently activated, four tactile patterns couples were designed. As illustrated in Fig. 2, these patterns consist in the activation of either 2 or 4 electrodes located on the front, the left, the back and the right side of the matrix.

Task and procedure

Before the experiment, a set of instructions describing the different kind of patterns used in the study and the required task was orally presented to each participant. The TDU matrix was placed inside the mouth and in contact with the anterior superior surface of the tongue. The participants were asked to keep the TDU in the same position on her/his tongue as well as possible during the course of the experiment.

Two groups carried out the experiment. The first group performed the test with continuous waveform (CT condition) and the second group with discontinuous waveform (DCT condition). For each group, the patterns with two electrodes were first tested, followed by the patterns with four electrodes.

A psychophysical technique referred to as the "method of limits" [START_REF] Gescheider | Psychophysics: The Fundamentals[END_REF] was employed in determining each pattern's liminary thresholds for each subject. An ascending threshold was established by beginning with a very low, previously tested stimulation threshold, gradually increasing by steps of 150 mV until the subject perceived a stimulation. A descending threshold was established by beginning with a stimulation clearly detectable (previously tested) and gradually decreasing by steps of 150 mV until the subject could no longer perceive any stimulation. A dedicated software administered the test of the liminary threshold detection. A total of 5 ascending and 5 descending series were alternatively run to each of the eight patterns for each subject. During the continuous stimulation test, the stimuli intensity was maintained 3s then increased or decreased as a function of the phase (ascending or descending). In the case of the discontinuous stimulation, the stimuli were sent during 1 s every 3 s in order to avoid a sensory adaptation and evoke a good recovery.

The subject was required to signal when she/he was absolutely certain of the appearance (ascending test) or the disappearance (descending test) of the lingual stimulation, by pressing a user-defined key on the laptop placed in front of her/him. To indicate the direction of the test, an up or down arrow was shown on the computer screen. A red point was also shown to specify when there was stimulation. Each testing session lasted about 2 hours.

Data analysis

The experiment yielded a total of 100 thresholds per pattern tested. The identified threshold for each pattern and for each subject was defined as the average of the median value for the 5 ascending responses and the median value for the 5 descending responses. The threshold values were calculated as the value at 50% of the 5 ascending or descending responses, with a linear interpolation supposing that the probability variation was linear between two values which surrounded 50% .

Proceedings of ENACTIVE/07

4th International Conference on Enactive Interfaces Grenoble, France, November 19 th -22 nd , 2007 Data were computed and treated with Statistica statistics software. They were then submitted to four separated 2 Groups (Continuous vs Discontinuous) x 2 Electrode Patterns (2 electrodes vs 4 electrodes) analyses of variances (ANOVAs). The level of significance was set at 0.05. Results show a main effect for the group, for the front position (F (1;18)=15.01, p < 0.01, figure 3A). Liminary threshold is higher for the continuous than for the discontinuous condition.

Results

Moreover, the results show a significant difference between 2 and 4 electrode patterns. For all positions, the patterns with 2 electrodes require an intensity threshold higher than the patterns with 4 electrodes (front, left, back and right, respectively F (1;18)=20.27, p < 0.001; F (1;18)=13.33, p < 0,01; F (1;18)=18.68, p < 0.001 and F (1;18)=16.13, p < 0.001).

Finally, the ANOVA shows a significant interaction between the groups and electrode numbers conditions for the front position (Fig. 3A), (F (1;18)=6.22, p < 0.03). Planned comparison analyses evidence a significant difference between two and four electrodes for continuous stimulation (F (1;18)=24.48, p < 0.001), whereas no significant difference was observed for the discontinuous stimulation (F (1;18)=2.01, p > 0.17).

Discussion

Tactile vibratory sensations are transmitted on the tongue surface via electrotactile stimulation. The electric current passing through the skin stimulates cutaneous afferent fibers at the location of the tactile electrodes. Our laboratory addressed various medical applications using the electrotactile feedback to supplement the damaged sensory information or provide accurate orientation information not easily accessible.

The front part of the tongue shows a significant difference for the continuous vs discontinuous conditions. As previously reported in the literature [START_REF] Maeyaman | Minimal two-point discrimination in human tongue and palate[END_REF] [10] [START_REF] Pleasonton | Sensitivity of the tongue to electrical stimulation[END_REF] [13], the sensitivity of the tongue tip is higher than all the others parts. Interestingly, our result further demonstrated that on the tip tongue, the liminary threshold is higher using continuous stimulations than discontinuous. This high tip tongue sensitivity of the tongue anterior part can explain that the accommodation effect is faster with the continuous than with discontinuous stimulations. Indeed, after accommodation, the induced sensation becomes less clear and less discriminable. Therefore, the performances of the subjects are reduced. An option to limit this adaptation could be to find the best stimulation frequency to improve the subjects performance (alarm or continuous guidance) and to optimize the tongue tactile neuronal adaptation/recovery in order to allow a regular use.

Proceedings of ENACTIVE/07

4th International Conference on Enactive Interfaces Grenoble, France, November 19 th -22 nd , 2007

In this experiment, all the electrode number main effects were significant. For each position (front, left, back and right), the electrotactile liminary threshold was higher for two electrodes than for four electrodes. This result is quite natural and consistent with the literature [START_REF] Ajdukovic | The relationship between electrode area and sensory qualities in lelectrical tongue stimulation[END_REF]. However, it is interesting to compare the reduced intensity saved with respect to the electrode number activated to optimize the energy consumption. Although the two and four electrode patterns intensity liminary thresholds are very closed (only a 0,9 ratio between both), there are twice more activated electrodes between both patterns. Consequently, patterns with two electrodes should be preferred to patterns with four electrodes, assuming that the subjects are still able to efficiently perceive such patterns.

Moreover, the observation of a significant interaction between both conditions (CT/DCT vs Electrode number) on the front part shows that the liminary threshold difference between two and four electrode patterns observed for the continuous stimulation is not present under discontinuous stimulation. These results suggest that, for energy saving and to improve the yield, it would probably be better to use discontinuous stimulation with two electrode patterns.

Obviously, complementary studies should be performed to investigate more quantitatively lingual electrotactile stimulation. Psychophysical experiments have shown that it possible to provide electrotactile sensation on the tongue with a very light current. Also, further experiments could allow optimising the comfortability threshold of the wireless TDU for a daily use.

Fig. 1

 1 Fig. 1 On the left panel, the Tongue Display Unit. It constitutes of 2D electrode array (1,5 cm × 1,5 cm) including 36 gold-plated contacts each with 1.4 mm diameter, arranged in a 6×6matrix. On the right panel, the TDU placed on the tongue.

Fig. 2

 2 Fig. 2 Electrode Schemes of the four patterns couples tested on the anterior part of the tongue. Left panel: patterns with 2 electrodes; right panel: patterns with 4 electrodes.

Fig. 3

 3 Fig. 3 Mean intensity liminary threshold for CT and DCT condition as a function of electrode number, 2 or 4, for A) the front, B) the left, C) the back, D) the right positions.

Figure 3

 3 Figure 3 illustrates the mean intensity liminary threshold (± standard deviations) obtained for CT and DCT conditions as a function of number of activated electrodes (2 or 4), for the four front (A), left (B), back (C) and right (D) positions.Results show a main effect for the group, for the front position (F (1;18)=15.01, p < 0.01, figure3A). Liminary threshold is higher for the continuous than for the discontinuous condition.Moreover, the results show a significant difference between 2 and 4 electrode patterns. For all positions, the patterns with 2 electrodes require an intensity threshold higher than the patterns with 4 electrodes (front, left, back and right, respectively F (1;18)=20.27, p < 0.001; F (1;18)=13.33, p < 0,01; F (1;18)=18.68, p < 0.001 and F (1;18)=16.13, p < 0.001).Finally, the ANOVA shows a significant interaction between the groups and electrode numbers conditions for the front position (Fig.3A), (F (1;18)=6.22, p < 0.03). Planned comparison analyses evidence a significant difference between two and four electrodes for continuous stimulation (F (1;18)=24.48, p < 0.001), whereas no significant difference was observed for the discontinuous stimulation (F (1;18)=2.01, p > 0.17).

Acknowledgement

This research was supported by the Gravit ("Grenoble Valorisation") consortium.


Louis Guillot
email: guillot@lma.cnrs-mrs.fr

Christophe Vergez
email: vergez@lma.cnrs-mrs.fr

Bruno Cochelin
email: bruno.cochelin@centrale-marseille.fr

Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method

Keywords: Nonlinear dynamics, Delay equations, Harmonic balance method, Numerical continuation, Periodic solutions, Quadratic recast, Asymptotic Numerical Method

This article presents an extension of the Asymptotic Numerical Method combined with the Harmonic Balance Method to the continuation of periodic orbits of Delay Differential Equations. The equations can be forced or autonomous and possibly of neutral type. The approach developed in this paper requires the system of equations to be written in a quadratic formalism which is detailed. The method is applied to various systems, from Van der Pol and Duffing oscillators to toy models of clarinet and saxophone. The Harmonic Balance Method is ascertained from a comparison to standards time-integration solvers. Bifurcation diagrams are drawn which are sometimes intricate, showing the robustness of this method.

Introduction

The aim of this article is to extend the framework of the numerical continuation of periodic solutions of differential systems using the Harmonic Balance Method (HBM) coupled with the Asymptotic Numerical Method (ANM) to time-delayed systems.

For dynamical systems without delay, the efficiency of this method has already been shown in several works, see [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF] and [START_REF] Vigué | Regularized friction and continuation: Comparison with coulomb's law[END_REF] for some examples. The work [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF] treats the case of quasiperiodic solutions using the same framework. To achieve better efficiency, the choice to recast the original system into a quadratic version has been done. The procedure to obtain such a recast has been explained in great details in [START_REF] Guillot | A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] and is briefly recalled in this paper. The way to recast nonpolynomial non-linearities was explained for the first time in [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF] and a way to recast fractional order derivatives is already known from [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF]. As opposed to already existing methods of continuation of delay differential equations (DDE), see [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using dde-biftool[END_REF] and its adaptation for neutral systems [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF], the method developed in this paper is based on the frequency representation of the periodic solutions.

This article contains three sections. In the first section the ANM and the HBM are recalled. It is admitted that the method can be applied to virtually any explicit or implicit differentialalgebraic equations. It is extended in a second section to the case of explicit or implicit delay algebraic-differential equations which can be forced or autonomous. Neutral systems can also be treated. The effect of time-delay on Fourier series is derived and the use of auxiliary variables to write the delayed variable in the frequency domain is detailed. This makes the quadratic recast of a delayed variable clear. In the last section, the method is applied to several examples, showing its versatility and its robustness particularly when applied to neutral systems. The convergence of the Fourier series with the order of truncation is shown and the solutions are compared to standard time-integration solvers. The continuation with respect to the delay is also explored in the first example.

2 Continuation of periodic solutions

Harmonic Balance Method

As described in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], the periodic solutions t → x(t) of a dynamical systems are represented by their truncated Fourier decomposition which is found through the Harmonic Balance Method :

x(t) = H k=-H x k exp(ikωt) = a 0 + H k=1 a k cos(kωt) + b k sin(kωt) (1)
where ω is the pulsation of the oscillations and H is the number of harmonics i.e the order of truncation of the Fourier series. Given that x -k is the complex conjugate of x k , the complex form and the real form of the Fourier series are equivalent, hence both can be used depending on which is more convenient. For example, let's seek periodic solutions of

ẍ + 2µ ẋ + x + x 3 = F cos(ωt) (2)
The system is rewritten quadratically as in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] :

ẍ + 2µ ẋ + x + xy = F cos(ωt) y = x 2 (3)
The next step is to compute products and derivatives directly on the complex form of the Fourier series [START_REF] Atay | Van der pol's oscillator under delayed feedback[END_REF]. The products can be written as a convolution :

H k=-H x k exp(ikωt) H k=-H y k exp(ikωt) = H k=-H exp(ikωt) h1+h2=k |h 1 |≤H,|h 2 |≤H x h1 y h2 (4)
The derivatives can be written :

ẋ(t) = H k=-H ikωx k exp(ikωt) ẍ(t) = H k=-H -k 2 ω 2 x k exp(ikωt) (5)
Now it is possible to write (3) in the frequency domain by replacing with the formula (4) and [START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF] and by balancing the harmonics. This leads to a set of algebraic equations.

Continuation through Asymptotic Numerical Method

The Asymptotic Numerical Method (ANM) that has first been described in [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF][START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF], relies on a high-order Taylor series representation of the solution-branch. This technique has already proven its efficiency for a lot of applications in engineering [START_REF] Cottanceau | A finite element/quaternion/asymptotic numerical method for the 3d simulation of flexible cables[END_REF], mechanics [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF] or acoustics for example. While some implementations relying on automatic differentiation do exist, see [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF] for example, the choice is made here to work with a quadratic framework. A generic implementation of this latest approach which minimizes problem-dependent implementation has been developed [START_REF] Guillot | A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF]. It is based on the numerical continuation of algebraic systems of the form

R(V ) = 0 (6)
where V ∈ R n+1 and R(V ) ∈ R n is an analytic function of its argument. This system is always written in a quadratic format as a prerequisite of the method. Again, this formalism is not a constraint that we suffer but a choice that allows to treat a very wide range of problems as shown in [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF]. The system (6) can then be written

R( Ṽ ) = C + L( Ṽ ) + Q( Ṽ , Ṽ ) = 0 (7)
where Ṽ ∈ R N +1 , N ≥ n and C, L and Q are respectively a constant, a linear and a quadratic operator. For the small example (3), the introduction of an auxiliary variable y = x 2 in the time domain allows to have almost directly a quadratic algebraic system. The only non-quadratic terms are due to the second order derivative (5) that leads to cubic terms x k ω 2 . The introduction of the auxiliary variable Ω = ω 2 solves this problem. The coupling between ANM and HBM in the quadratic framework [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] can be automatized and is detailed in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. Again, the key idea of the ANM is to compute a high-order Taylor series development of the solution-branch. Let consider that Ṽ0 is a regular solution of the system [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF]. Let Ṽ1 be a tangent vector at Ṽ0 . The classical arc-length parameter a = Ṽ t 1 ( Ṽ0 -Ṽ ) is introduced and the solution-branch around Ṽ = Ṽ0 is written as a power series with respect to a

Ṽ (a) = Ṽ0 + a Ṽ1 + a 2 Ṽ2 + a 3 Ṽ3 + • • • (8)
More on the implementation of continuation in the case of algebraic systems coming from the harmonic balance method is available in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. The next section generalizes this approach for systems with time-delayed variables in their formulation.

Continuation of periodic solution of systems with time delay

In this section a system of delay differential equations, forced or autonomous and possibly neutral, of the following general form is considered:

f (t, x(t), ẋ(t), x(t -τ ), ẋ(t -τ ), λ) = 0 (9)
Where f : R × R N 4 × R → R N is an analytic function that is periodic with respect to its first variable, x : R → R N is an unknown periodic function of unknown pulsation ω, τ is a real number called the delay and λ is a control parameter of interest called the continuation parameter. The continuation parameter λ may be equal to the delay τ as shown in section 4.1.

A systematic way to recast the equations [START_REF] Colinot | Formulation analytique de la pression interne d'un modèle idéalisé d'instrument conique à anche[END_REF] quadratically is first to introduce three auxiliary variables :

       f (t, x(t), y(t), x τ (t), y τ (t), λ) = 0 y(t) -ẋ(t) = 0 x τ (t) -x(t -τ ) = 0 y τ (t) -y(t -τ ) = 0 (10)
The system has now 4N +1 unknowns for 4N equations. The equations on delayed variables appear isolated and a classical quadratic recast [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF][START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF][START_REF] Guillot | A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] can be applied to the function f to write the system (10) quadratically :

       g(t, x(t), y(t), x τ (t), y τ (t), λ, v(t)) = 0 y(t) -ẋ(t) = 0 x τ (t) -x(t -τ ) = 0 y τ (t) -y(t -τ ) = 0 (11)
where v : R → R Naux is a function of the time that can be expressed as quadratic expressions of all the variables. The system has now 4N +1+N aux unknowns for 4N +N aux quadratic equations. The two last equations are the only one that are not treated in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF][START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF][START_REF] Guillot | A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF], thus, the Fourier coefficients of x τ are now derived from the Fourier coefficients of x with their quadratic recast.

Harmonic Balance Method applied to delayed variables

The order of truncation of the Fourier series is noted H. The truncated Fourier series of t → x(t) is:

x(t) = a 0 + H k=1 a k cos(kωt) + b k sin(kωt) (12)
The Fourier series of t → x(tτ ) is:

x(t -τ ) = a 0 + H k=1 a k cos(kω(t -τ )) + H k=1 b k sin(kω(t -τ )) = a 0 + H k=1 a k (cos(kωt) cos(kωτ ) + sin(kωt) sin(kωτ )) + H k=1 b k (sin(kωt) cos(kωτ ) -cos(kωt) sin(kωτ )) = a 0 + H k=1 [a k cos(kωτ ) -b k sin(kωτ )] cos(kωt) + H k=1 [b k cos(kωτ ) + a k sin(kωτ )] sin(kωt) (13)
Therefore, if a τ k and b τ k are the Fourier coefficients of t → x(tτ ), their expression is given by:

a τ 0 = a 0 a τ k = a k cos(kωτ ) -b k sin(kωτ ) , k ≥ 1 b τ k = b k cos(kωτ ) + a k sin(kωτ ) , k ≥ 1 (14)

Asymptotic Numerical Method applied to delayed variables

For k ≥ 1, the auxiliary variables c τ k = cos(kωτ ) and s τ k = sin(kωτ ) are defined to get a quadratic recast for ( 14):

τ = ωτ a τ 0 = a 0 a τ k = a k c τ k -b k s τ k b τ k = b k c τ k + a k s τ k s τ k = sin(kτ ) and d(s τ k ) = kc τ k × d(τ ) c τ k = cos(kτ ) and d(c τ k ) = -ks τ k × d(τ ) (15)
where in the last two equations, the differentiated forms are used as explained in [START_REF] Guillot | A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] to compute the Taylor series expansions (8) of s τ k and c τ k . The last equations are defined as "initial conditions" plus a "differential system". For the computation of the Taylor series (8), the "initial conditions" are used to compute the zero-th order coefficient Ṽ0 and the differential system to compute higher order coefficients Ṽp , p ≥ 1. This type of computation of Taylor series development is explained in details in the book [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]. The main ideas are recalled below.

Let x = x 0 + ax 1 + a 2 x 2 + . . . be a given Taylor series and let s = sin(x) and c = cos(x) be its sine and its cosine. To compute the coefficients of the Taylor series of s (= s 0 + as 1 + a 2 s 2 + . . . ) and c (= c 0 + ac 1 + a 2 c 2 + . . . ), the following procedure is applied:

s 0 = sin(x 0 ) and c 0 = cos(x 0 ).

-For p ≥ 1, identify the term of order p -1 in the differentiated equations ds = cdx and dc = -sdx.

For the second step, the development in Taylor series of the differentiated equations are written:

s 1 + 2as 2 + 3a 2 s 3 + . . . = c 0 + ac 1 + a 2 c 2 + . . . x 1 + 2ax 2 + 3a 2 x 3 + . . . c 1 + 2ac 2 + 3a 2 c 3 + . . . = -s 0 + as 1 + a 2 s 2 + . . . x 1 + 2ax 2 + 3a 2 x 3 + . . . (16)
For p = 1, the identification of the constant term of the system [START_REF] Hu | Resonances of a harmonically forced duffing oscillator with time delay state feedback[END_REF] gives s 1 = c 0 x 1 and c 1 = -s 0 x 1 . For p ≥ 1, the identification of the terms of order p -1 in the system [START_REF] Hu | Resonances of a harmonically forced duffing oscillator with time delay state feedback[END_REF] gives the following recursive formula:

s p = 1 p [c 0 × px p + c 1 × (p -1)x p-1 + • • • + c p-2 × 2x 2 + c p-1 × x 1 ] c p = -1 p [s 0 × px p + s 1 × (p -1)x p-1 + • • • + s p-2 × 2x 2 + s p-1 × x 1 ] (17)

Applications

In this section, four examples are given which have been explored with our method. The first one, a delayed Van der Pol oscillator [START_REF] Atay | Van der pol's oscillator under delayed feedback[END_REF], is autonomous. It is used to show the convergence of the Fourier series to a reference solution of the system, computed with a standard solver [START_REF] Shampine | Solving ddes in matlab[END_REF]. Moreover, a bifurcation diagram with respect to the delay is drawn and the location of the many Hopf bifurcations is ascertained with the results of a linear stability analysis. The second example, a forced delayed Duffing oscillator [START_REF] Hu | Resonances of a harmonically forced duffing oscillator with time delay state feedback[END_REF], exhibits a lot of resonances (both super-harmonic and sub-harmonic). Its bifurcation diagram is represented with a focus on two symmetry-breaking bifurcations. The third example is a toy model of clarinet [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] that is purely algebraic. Its periodicsolutions are squared shaped and their Fourier series converge rather slowly. This example shows the robustness of our approach. The fourth and last example is a neutral model of saxophone [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF][START_REF] Colinot | Formulation analytique de la pression interne d'un modèle idéalisé d'instrument conique à anche[END_REF] using transcendental functions in its writing. A solution is compared to a finite-difference scheme from the literature [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] and a bifurcation diagram is drawn, showing an expected behaviour for this type of system according to [START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the helmholtz motion analogy[END_REF].

A delayed Van der Pol oscillator

In this part we consider the Van der Pol oscillator (dimensionless) with delayed feedback that has been studied in [START_REF] Atay | Van der pol's oscillator under delayed feedback[END_REF]:

ü(t) -ε(1 -u(t) 2 ) u(t) + u(t) = εku(t -τ ) ( 18
)
The continuation is performed with respect to τ . Introducing the auxiliary variables

v(t) = u(t), r(t) = 1 -u(t) 2 and F nl (t) = v(t)r(t) -ku(t -τ ).
In order to minimize the number of auxiliary variables, the variable F nl is the sum of the delayed variable ku(tτ ) (with k constant) and the product v × r. Hence the delayed variable is not isolated like in ( 10) and ( 11), but as it appears linearly, this allows to use the approach explained in the previous section. The following quadratic recast is used:

v(t) + u(t) -εF nl (t) = 0 v(t) -u(t) = 0 r(t) -1 + u(t) 2 = 0 F nl (t) -v(t)r(t) -ku(t -τ ) = 0 (19)
To validate the method, the figure 1 shows the difference between time integration and harmonic balance for the periodic solution of equation ( 18) obtained with ε = 2, τ = 0.2 and k = 0.2. The number of harmonics is varied to show the convergence of the harmonic balance towards the time integration solution. The reference solution is a time integration performed with dde23 built-in solver of MATLAB. Its absolute and relative tolerance thresholds have been set to 10 -8 .

With a low number of harmonics, i.e for H < 10 , the solution is not accurate and even the period of the oscillation is not computed correctly. Indeed, the plot H = 5 stops before t = 10 s, conversely to the other plots. From H = 25 to H = 50, the solution gets much closer to the reference. Solutions for H = 100 and H = 200 look superimposed. In fact, the tolerance at which the reference has been computed, i.e 10 -8 , is reached for H = 50. When computing the solution with a higher number of harmonics, the time-integration reference is probably farther from the true solution than the solution obtained with the harmonic balance method. In addition, the timeintegration solver used takes approximately ten times the time needed to compute the solution with H = 50. Moreover, once a solution has been computed the computation of the bifurcation diagram through ANM is very fast : for the same computation time as one time-integration over twenty periods (to get rid of the transient), all the steady state solutions between ε = 2 and ε = 5 are available with our method.

The figure 2 shows a part of the bifurcation diagram of the system [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] with the delay τ as a continuation parameter. There are infinitely many Hopf bifurcations in this system. The analytic position of the Hopf bifurcations is determined by the stability analysis of the equilibrium u = 0 and is given by the solutions of the system [START_REF] Atay | Van der pol's oscillator under delayed feedback[END_REF]: Fig. 1 Differences over one period between several solutions of equation ( 18) obtained with the harmonic balance method with increasing H and a reference solution obtained with time integration. Here ε = 2, τ = 0.2 and k = 0.2.

1 -ω 2 = kε cos(τ ω) εω = kε sin(τ ω) (20)
where ω is the pulsation of the arising periodic solution at the corresponding Hopf point. The continuation of the periodic solution branches has been initialized with an analytic approximation with H = 1 and the pulsation ω solution of the system (20) for the corresponding value of τ . In this particular case, the robustness of the ANM allows to follow all the solution branches from just one starting point. In the general case, it would probably necessitate to initialize with an analytic approximation at every Hopf bifurcations. On the figure 2 some parts of the diagram are of null amplitude : these are just some parts where the solution found is the equilibrium u = 0.

Fig. 2 Bifurcation diagram of equation ( 18) with respect to the delay τ , ε = 0.1 and k = 2. There are as many Hopf bifurcations as branches starting from a null amplitude. The red dots are the analytical position of Hopf bifurcations given by the solutions of system (20).

A forced delayed Duffing oscillator

The example presented in this section is taken from [START_REF] Hu | Resonances of a harmonically forced duffing oscillator with time delay state feedback[END_REF]. It is an example of a forced delayed differential equation. The dimensionless equation of motion (Eq.( 3) in the reference) is :

ẍ(t) + 2ζ ẋ(t) + x(t) + µx 3 (t) = 2ux(t -τ ) + 2v ẋ(t -τ ) + f cos(λt) (21)
The continuation is performed with respect to λ. The following quadratic recast is used :

ẋ(t) -y(t) = 0 ẏ(t) + 2ζy(t) + x(t) + µr(t)x(t) -2ux(t -τ ) -2vy(t -τ ) = f cos(λt) r(t) -u(t) 2 = 0 ( 22
)
where λ is the continuation parameter and the set of parameters used is : ζ = ε, µ = ε, u = ε, v = -ε, τ = π and f = 5 with ε = 0.05. Again, the time-delayed terms are not isolated but as they appear linearly in the equations, the recast of section 3 is applicable. In [START_REF] Hu | Resonances of a harmonically forced duffing oscillator with time delay state feedback[END_REF], the primary and the 1/3 subharmonic resonances were studied with the method of multiple scales. In this present article a whole bifurcation diagram is drawn, showing both primary and 1/3 subharmonic resonances. To find the starting point of the continuation, the approximations computed in [START_REF] Hu | Resonances of a harmonically forced duffing oscillator with time delay state feedback[END_REF] are used. The forcing amplitude has been chosen high to have a rich dynamics as can be seen in the bifurcation diagram which is represented in figure 3. To draw this diagram, the order of truncation of the Fourier series has been set to H = 25 which is sufficient away from superharmonic resonances. The 1/3 subharmonic resonance is disconnected from the main branch. It corresponds to a solutionbranch where the pulsation of the oscillation verifies ω = λ 3 .

Fig. 3 Bifurcation diagram of the amplitude of the solution with respect to the forcing pulsation λ. On the primary resonance solution-branch, there are several super-harmonic resonances appearing for λ = 1 3 , 1 5 , 1 7 , 1 9 , . . . . There are also two symmetry-breaking bifurcation, noted SB on the zoom of the diagram, detected with the method explained in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] for which even harmonics appear.

The zoom on the top of figure 3 shows the many superharmonic resonances and the symmetrybreaking (SB) bifurcations. The superharmonic resonances appear for λ 1 2p+1 , p ∈ N and in these areas the solutions become more complex as the harmonics 2p + 1 are close to the natural pulsation of the system. In addition, two symmetry-breaking bifurcations are detected and the emerging solution is then continued with the method explained in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF]. On the bifurcated branches, the solutions have strong even harmonics appearing and their phase diagrams (x, ẋ) are no longer symmetric.

A clarinet toy model

The following example is a simplified model of clarinet-like instrument, taken from [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] [Eq.( 4)-( 5)-( 6)-( 7)-( 10)]. This system has been chosen to highlight that our approach copes even with a degenerate form of the general system [START_REF] Colinot | Formulation analytique de la pression interne d'un modèle idéalisé d'instrument conique à anche[END_REF]. In this example, there are no time derivatives in the equations of the motion. The unknowns are functions of the time and the dynamics comes from the coupling between a non-linearity and an equation with time delayed variable. As in the reference, the reed dynamics is neglected and the dimensionless equations can be written :

u(t) = 0 if p(t) ≤ γ -1 ζ(p(t) -γ + 1) |γ -p(t)|sign(γ -p(t)) if p(t) > γ -1 p -(t) = -λp + (t -2ℓ c ) (23)
where p and u are respectively the acoustic pressure and the volume flow in the mouthpiece of the instrument. p + and p -are respectively the forward and the backward progressive waves, linked to p and u through p = p -+ p + and u = p +p -. ℓ is the length of the instrument and c is the speed of sound. The delay 2ℓ c in equation ( 23) represents the time for a pressure wave to go back and forth the cylindrical bore of the clarinet. λ is the reflection coefficient at the end of the instrument which can be seen as a damping term. γ is the static blowing pressure in the mouth of the player and ζ is the reed parameter. All the variables but the time are dimensionless in this model. In the following, the time is normalized by 2ℓ c so that the value of the delay τ is 1. For simplicity, the nonlinear characteristics (the first equation of the system ( 23)) is written under the assumption that -1 < pγ < 0 which means according to [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] that the reed channel is open and air is allowed to enter the clarinet. The system [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF] becomes

u(t) = ζ(p(t) -γ + 1) γ -p(t) p -(t) = -λp + (t -1) (24)
Its quadratic recast is

p -(t) = -λp + (t -1) v(t) 2 = γ -p(t) u(t) = ζ(p(t) -γ + 1)v(t) (25)
The variable v(t) is defined implicitly by the second equation. To ensure that the results are correct, we have to check that v(t) is non-negative. The parameters chosen are λ = 0.95 and ζ = 0.8. The continuation is performed with respect to γ, the dimensionless blowing pressure. The bifurcation diagram is plotted in figure 4. The values for which the system encounters a Hopf bifurcation γ H = 0.35 or a period doubling bifurcation γ P D = 0.43 are the same as the one given in [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] (table 1) up to the third significant digit. On the figure 4, the Hopf bifurcation occurs at the lowest point of the solution-branch. As the mean value of p is not zero, the L 2 norm of p is not 0 at the bifurcation point. The exact solution of this system in the time-domain is a square and the decomposition in Fourier series is obviously not very well suited for this discontinuous type of signal. However, the numerical method is so robust that it is able to continue this solution on the whole domain which is remarkable. This also explains why the diagram represents the L 2 norm of p and not its amplitude : Gibbs phenomenon leads to an overshoot at each of the two discontinuous point of the signal. This gives wrong amplitude for the acoustic pressure. The analytical solutions given in [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF] were used as starting points for the continuation process, using their truncated decomposition in Fourier series.

A neutral saxophone model

In this section, a neutral system that is a dimensionless model of saxophone taken from [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] is studied. It has also been the subject of the conference talk [START_REF] Colinot | Formulation analytique de la pression interne d'un modèle idéalisé d'instrument conique à anche[END_REF]. The derivative of Eq.( 40) in [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] is taken to get rid of the integrals. This will necessitate to be taken into account when considering the Fourier developments of the unknowns. The following equation is obtained :

p + √ 3(1 -F ′ (p)) ṗ + p = pτ + √ 3(-1 -F ′ (p τ )) ṗτ + p τ ( 26
)
where p is the acoustic pressure in the mouthpiece of the instrument; γ and ζ have the same meaning as in section 4.3. Everytime there is a τ index, it means that the variable is delayed, i.e p τ (t) = p(tτ ). In order to have a smooth function, the nonlinear characteristics F (p) is an exponential regularized law [START_REF] Colinot | Formulation analytique de la pression interne d'un modèle idéalisé d'instrument conique à anche[END_REF] :

F (p) = ζ exp(1) √ 3 exp(3(p -γ))(γ -p) (27)
To recast the equation (26) quadratically, the idea is to write it on the form I + √ 3 ṗ = I τ -√ 3 ṗτ with I = p -√ 3F ′ (p) ṗ + p. It is then necessary to write F ′ (p) quadratically. This is done using the same idea as for the recast in equation [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF] ; setting ∆p = pγ, the exponential is written using A = exp(3∆p) and its differentiated form dA = 3Ad(∆p). These remarks lead to introduce the following auxiliary variables :

∆p = p -γ A = exp(3∆p) F ′ (p) = ζ exp(1) √ 3 A(-1 -3∆p) q = ṗ I = p + p - √ 3F ′ (p)q (28)
which allows to write the equation (26) in a quadratic format :

I -I τ + √ 3(q + q τ ) = 0 ∆p -(p -γ) = 0 A -exp(3∆p) = 0 and dA -3A d(∆p) = 0 F ′ (p) -ζ exp(1) √ 3 A(-1 -3∆p) = 0 q -ṗ = 0 I -(p + p - √ 3F ′ (p)q) = 0 (29)
The definition of A (third line, left hand side in (29)) is used only to compute the mean value (or the constant term) of its Fourier series and its differentiated form (third line, right hand side in (29)) is used to compute higher order terms. This has been explained in [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF]. The dimensionless delay τ can be expressed in term of the length of the conical resonator ℓ of the saxophone, truncated at its input to plug the mouthpiece, and the length x 1 of the missing part of the cone (from the apex to the input of the truncated cone). This is not our purpose here to understand the refinements of this model, the interested reader can find the details in [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF]. Below, the value of τ is set to 2 equation. This implies that after the Harmonic Balance Method (HBM) has been performed, the constant coefficient equation (i.e its mean value) states 0 = 0. As this leads to a singular system, this equation is replaced by p 0 = 0 where p 0 is the constant term [the a 0 of the Fourier decomposition ( 12)] of p.

The waveform of each unknown over one period is a by-product of the method and is available for each value of γ, i.e for each point of the bifurcation diagram 5. The waveform of p is compared to the one obtained with the numerical integration algorithm presented in [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] for γ = 0.5 and is in very good agreement as can be seen in figure 6. On the same figure, the pressure obtained by integration with the built-in solver ddensd [START_REF] Shampine | Dissipative approximations to neutral ddes[END_REF] of MATLAB is presented. The three methods give very close results. However, as expected in [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF], their scheme is unstable and the Hopf bifurcation Fig. 6 Comparison of the periodic solutions of equations (26) over one period for ζ = 1.2 and γ = 0.5. In blue, the solution obtained with Harmonic Balance Method with H = 50 and extracted from the bifurcation diagram shown in figure 5. In red dashed, the solution obtained with the time-integration algorithm of [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] with 5000 samples per period. In yellow dotted, the solution obtained with ddensd Matlab's solver with absolute and relative thresholds set to 10 -3 . A zoom on the zone of beating reed is shown to enhance the differences between the waveforms. is detected around γ = 0.285, that is below the theoretical value of γ = 1 3 which is found with our method.

The scheme used in [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] is very fast but requires a high sampling frequency (more than 1000 samples per period) to converge. The ddensd solver takes a longer time especially when the tolerance chosen is decreased to the minimum value possible of 10 -5 . In our approach, the transient is not available but the shape of the steady state is obtained on all the solution-branch instantly through the Fourier series construction, as well as the unstable periodic regimes. In addition, the Hopf bifurcation and the limit point can be deduced from the bifurcation diagram.

The model of saxophone presented should have a zero mean over one period by physical consideration [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF]. For the signals represented in figure 6, the acoustic pressure obtained by Harmonic Balance has a mean value around the numerical precision 10 -16 while the time integration solver from [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] gives a mean around 10 -4 and ddensd gives 10 -3 .

Conclusion

An extension of the HBM, coupled with the ANM continuation, to time-delayed systems is developed. Various types of systems exhibiting constant time-delays is explored. Especially, delayed differential equations of neutral type do not require a specific setting in the framework presented. The available tools of the ANM continuation allow to compute bifurcated branches without the need of a specific implementation.

The coupling between HBM and ANM shows its efficiency and its robustness compared to time-integration methods. Especially in some cases where there are not many numerical methods available (neutral systems) and even when the Fourier series of the periodic solution have bad convergence properties. It is ascertained by a comparison to time-domain solvers. The algorithms used for this article are freely available online at http://manlab.cnrs-mrs.fr.

Fig. 4

 4 Fig.4Bifurcation diagram of the toy model of clarinet[START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF] with respect to the dimensionless blowing pressure γ, with H = 200. The period doubling bifurcation is visible as the connection of the branches. The red dotted lines correspond to section of branches where the variable v(t) is negative, i.e where the simplified model is not valid anymore to represent the clarinet functioning.

√Fig. 5

 5 Fig.5Bifurcation diagram of a model of saxophone with respect to the dimensionless blowing pressure γ and for ζ = 1.2. The sub-critical Hopf bifurcations can be seen clearly. It is not surprising for the saxophone to show such sub-critical bifurcations[START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the helmholtz motion analogy[END_REF]. This diagram has been computed using H = 50 harmonics.

Acknowledgements The authors want to thank Tom Colinot for the help with the saxophone model and its associated time-integration method. This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR).


Pascal Andre

Christian Attiogbé

Arnaud Lanoix

Pascal André

Systematic Construction of Critical Embedded Systems using Event-B

Keywords: Embedded control systems, Specification method, Event-B patterns

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

Introduction

Modelling and analysis of complex systems without dedicated methods is painful, inefficient and time-consuming. Methods and tools are required for efficient system engineering; this is particularly true for formal software engineering.

Unlike many other types of software, embedded systems are often developed for specific target environments (processors, vehicles, medical devices, etc.) and very often they should run for long times (even years), once they have been implemented in their so called critical environments. Therefore, embedded systems and their construction have stringent robustness requirements; one have to develop them accordingly to get them reliable at runtime. The target environments for the development of each embedded system do not help the advent or the expansion of tools and methods dedicated to this type of software. But, there are numerous models for embedded real-time systems [START_REF]Communicating Embedded Systems: Software and Design[END_REF],

Considering that i) the requirements for reliability and correct construction of the models and the derived embedded systems are of great importance, and that ii) the development of these systems lacks of methods to guide the developers, we are motivated to contribute to fill the gap between these needs and the state of the art.

In this work we propose a correct-by-construction method dedicated to critical embedded control systems. This method, based on Event-B, is intended to guide step by step the specifier or the engineer to drive its development from requirements to concrete software, defining abstract models, and refining them in a systematic way.

In Section 2 we introduce the proposed method with the details of each step. Section 3 illustrates the application of the method on a common case study, called the landing gear system. In Section 4 we evaluate the application of the method and the case study and comment related studies. Section 5 concludes this work.

A Method to Construct Correct Embedded Systems

We present a stepwise and systematic method (named Heñcher) to construct critical embedded control systems using Event-B. We reuse the approach already established and demonstrated in several case studies [START_REF] Damchoom | Applying Event and Machine Decomposition to a Flash-Based Filestore in Event-B[END_REF][START_REF] Damchoom | Modelling and Proof of a Tree-Structured File System in Event-B and Rodin[END_REF][START_REF] Abrial | Modeling in Event-B -System and Software Engineering[END_REF], following which complex systems can be constructed by combining 1) horizontal refinement with feature augmentation where we have to build a global abstract model of a the whole system (a controller and its physical environment) (Sect. 2.1) and 2) structural refinement (making the abstract structures more and more concrete (Sect. 2.2). But we extend it and provide dedicated guidelines at different steps, which help the developers to reach quickly a correct control system. Fig. 1 illustrates the Event-B patterns from the most abstract model which describes only the interface of the controller, to the systematic decomposition into two parts: the Controller and the physical Environment. The high level state space of any control system can be described by the elicitation of the interface variables between the digital part (the controller) and the physical part (the controlled environment) of the considered system.

Step 1: Characterise the abstract model of a considered system Fig. 2. A generic shape for event-based model of a control system Fig. 2 depicts a general principle that may govern the organisation of event-based models of control systems. The dashed ovals are representative of the parametric events families; They should be replaced by the effective events related to the logic of a specific case study. Besides, the identified physical devices to be controlled should be precisely listed. The behaviour of each one will be specified later.

Step 1.1: Elicit the interface We distinguish in our method, a first step which consists in the description of the interface variables of the controller. There are three categories of variables at the interface of a controller with the controlled environment.

the input variables: they give the sensed state of the environment (the value given by sensors); they are read by the controller; the state variables: they are set and modified by the controller; they are used for monitoring the whole system; the order output variables: they are those used to send the orders that stimulate the physical environment; their values will be used by actuators.

These categories of variables will be used at different levels of the modelling and refinement of the system at hand: to introduce the events of the first abstract model and, for refining gradually the first model. Additionally, we have the category of internal variables, which are only used inside the controller.

Step 1.2: Elicit the global properties of the system The required system properties, including safety, liveness and non-functional properties should be explicitly named and listed in their informal form. These properties will be formalised and gradually introduced with their formal form during the various construction steps of our method.

Step 1.3: Start with a first abstract model Use the interface variables resulting from

Step 1.1 to build a first Event-B abstract model. This abstract model comprises in the related B machine clauses, the interface variables together with the appropriate abstract sets and properties which characterise them. This Event-B model will be enriched to obtain a global abstract model of the system including its control part, its physical part and their related properties. Global properties of the system (among those already listed in Step 1.2) should be formalised and introduced according to the available variables. Notice that the enrichment incorporates gradually the details of the physical environment (sensors and actuators) and the corresponding properties.

Step 1.4: List the events of the abstract model The Event-B abstract model resulting from the previous step will be enriched with a series of events built by defining a family of events related to each category of the interface variables: sense events, monitor events and stimulate events.

-Sense events family This family gathers the events used to set and to modify the values of the input variables. For each variable of this category, define an event named after the variable, with the prefix sense_. The link with the physical (state variables of the) environment is done later by refining these events in Step 2.2. -Monitor events family These events modify the appropriate state variables. For each variable of this category, define an event named after the variable with the prefix monitor_, to set the variable according to the current state of the controller and the input data. -Stimulate events family The events of this family modify the order output variables; each variable of this category is set with an event named after the variable with the prefix stmlt_. These events use the internal variables, the input variables and the output variables. Associated with these events to stimulate the physical devices, we may have as many events to stop the stimulation of the devices; accordingly these events have their name prefixed with stop_.

These three families of events, together with the reaction events family introduced later, are compliant with the standard sense-decision-control of the control cycle.

Step 2: Extend the abstract global model

Use feature augmentation [START_REF] Damchoom | Applying Event and Machine Decomposition to a Flash-Based Filestore in Event-B[END_REF][START_REF] Damchoom | Modelling and Proof of a Tree-Structured File System in Event-B and Rodin[END_REF][START_REF] Abrial | Modeling in Event-B -System and Software Engineering[END_REF] to integrate the controlled environment. This is precisely achieved on the basis of the sense events family, which in turn need the description of the controlled environment. The global properties listed before are also gradually formalised in the model, as invariants, as soon as the appropriate variables are available. Two sub-steps are distinguished but no matter their order during the development.

Step 2.1: Introduce the physical environment and the reaction events family It consists in adding successively to the model, events to propagate the values of interface output variables inside the physical environment. These events simulate or stop the behaviours of physical devices via the actuators. The feature augmentation is used to introduce the physical state variables, invariants and appropriate behaviours. Depending on the cases, either one simulates the behaviour of the physical devices with an abstract model, or the values of the output variables (from the interface) result in signals sent to the environment. In this last case we do not have dedicated events in the abstract model. Accordingly, the behaviours of the physical devices should be formally described. These behaviours, systematically guarded by the values of the output variables, may impact the state of the environment and finally they may impact the sensors.

State automata can often be used to capture the behaviour of a physical device; describing the automata with B events is then straightforward. The description of the physical part behaviour results in the family of reaction events. These events should be named using as prefix the identifier of the physical part that they describe.

Step 2.2: Detail the sense events family Each event of the sense family, updates an input variable according to the state of the sensors; for this purpose the event needs the model and the behaviour of its related sensor. Therefore the feature augmentation consists in introducing the model of the sensors and their related behaviour, as variables, invariants and related events. The behaviour of the sensors should consider the possible failures (anomalies, malfunctioning or physical defects); specific events should be described for each such possible failure. Practically during all the refinement steps, it is recommended to proceed incrementally with several small refinement steps dedicated to variables and events. This is necessary to master the proof complexity.

Step 3: Introduce the specific properties According to the system one has to build, besides the global properties gradually introduced with the variables, additional specific properties should be integrated at the abstract level to constrain the functioning of the system.

1. Reachability property with partial ordering: specific events (not at the same granularity with the Event-B events) with timestamps may be systematically used to order and to reason on reachability properties. 2. Non functional properties: specific properties related to nonfunctional requirements should be gradually introduced here. No matter the way they are described, provided that the mathematical support of Event-B is very large, and that external modules may be used to analyse them.

Vertical Process: Building the Concrete Parts of the System

The aim of this second process is to build the digital part and possibly the physical part of the system. The global Event-B abstract model resulting from the horizontal process should then be decomposed into various parts leading to specific components. At least we have a decomposition into a control software and a physical part. The decomposition can be performed as soon as one want to go into the details of one of the specific part by considering that the other part will stay as it is; that means no modification of the other part cannot be considered when we are refining a given one. Typically, from a decomposition step, the digital part will be refined until code by considering the events and variables of the physical part as they was at the decomposition step.

Step 4: Refine the global abstract model

We recommend to perform structural refinements as needed by the specific model to be refined. New internal B events may be added to refine the events of each family of events (sense, monitor, or stimulate). The state space variables of the global abstract model may be refined with more details in the invariant. At the end of this step, be sure that, the events of both parts are all in place, that the global required properties of the system are all in place (they cannot be introduced later after the decomposition).

Like with other formal models, an Event-B model can be animated, i.e. when appropriate values are provided for the variables in the model, its behaviour can be observed step by step according to the semantics of the model. Animation capabilities are helpful during all the refinement steps where we still have all the events of the global system together; it will not be possible to animate the whole system after its decomposition.

Step 5: Decompose into software and physical parts A decomposition paradigm is already supported by the Event-B method. It consists in splitting a given machine into several ones which will be refined independently. The decomposition splits the variables of the state space and/or the behaviour of the machine; however resulting machines cannot contradict each other by modifying the variables and their related properties once they have been separated. Two approaches exist for this purpose: the Abrial'style decomposition (called the A-style decomposition) [START_REF] Abrial | Refinement, Decomposition, and Instantiation of Discrete Models: Application to Event-B[END_REF] based on shared variables, and the Butler'style decomposition (called the B-style decomposition) [START_REF] Butler | Decomposition Structures for Event-B[END_REF][START_REF] Silva | Shared Event Composition/Decomposition in Event-B. In: 9th International Symposium onFormal Methods for Components and Objects[END_REF] based on shared events. In the A-style decomposition, events are first split between Event-B sub-components and then shared variables of the sub-components are used to introduce external events in the sub-components; these external events should be refined in the same way. In the B-style decomposition, variables are first partitioned between the sub-components and then shared events (which use the variables of both sub-components) are split between the sub-components according to the used variables.

We adopt the A-style decomposition which is more relevant when considering a list of specific events to be split relatively to a control part and a physical part. The methodological guide to achieve the decomposition is as follows: the digital part is made with all the events defined in the sense events, the monitor events and the stimulate events families whereas the physical environment gathers all the events defined in the reaction events families. Moreover, each part must have an abstract view of the other through external variables and events.

Step 6: Refine the control software and the physical environment separately

Step 6.1: Refining the control software Use Structural refinements based on the stimulate events family to refine the controller. The involved categories of variables are the input variables, the state variables and the output variables. Typically, the values of the output variables are synthesised from the other ones. This can be done through simple control functions or through sub-modules. When there are sub-modules, the input variables may be spawned inside the sub-modules; in the same way output variables may be updated by promotion from the sub-modules if any. Therefore one have to incorporate successively in the Event-B model the events to set and modify the output variables; they describe the result of the behaviour of the control part. State automata help to catch these behaviours; then the events of the B models encode the automata. We give now some recurrent patterns to help in modelling control part. i). Composition of several redundant sub-modules: when a controller is made of several redundant modules, it is straightforward to describe a generic module and use an indexing function to compose several instances of such modules (see Fig. 3).

-Encasing variables inside modules: the values coming from outside one or several modules can be systematically encased inside the modules with a dedicated event that spawn the events. -Promoting variables outside a module: in a symmetric way, the values going outside a module or several modules can be systematically described using a promotion pattern (with a dedicated event) for merging the output variables of the internal computing modules. ii). When the modules are not redundant, each one should be refined separately, but the treatment we have described for the inputs and outputs variables is the same.

Step 6.2: Refining the controlled (physical) environment Many cases can be considered depending on the system to be studied; either the physical devices are already available, or one has to build the physical devices from the formal models, or one has to build a part of the physical devices. Nevertheless, the exchange of signal with actuators is the standard way to act on physical devices.

i). In the case where the physical devices are available, with the related actuators, the refinement is straightforward; it consists for the events of the physical part to output the correct signal values (for example on/off values are encoded as voltage) as input of the actuators. But the physical devices may be emulated in preliminary studies before implementing the control part on the real devices. Mathematical models and dedicated system engineering tools are available as explained hereafter. ii). When the output variables of the controller cannot be immediately encoded as signal values, transformation can be achieved via appropriate mathematical functions and models. This should be done, starting from the requirements and the properties of the previous model, for example by external modules or functions written with tools like Matlab or Simulink1 or SciLab2 ; related works can be found in [START_REF] Su | Complementary Methodologies for Developing Hybrid Systems with Event-B[END_REF]. These tools generate executable codes dedicated to the target devices. They are also equipped with specific functions to handle time requirements.

A Running Case Study

The proposed method is applied on the Landing Gear (LG) case study, a benchmarking example proposed at the ABZ'2014 conference to compare different formal methods in terms of expressivity, performance, and ease of use. A prerequisite for reading this section is the detailed specification of this critical embedded system given in [START_REF] Boniol | The landing gear system case study[END_REF]. A summary of the LG system is depicted in Fig. 4. The LG system is in charge of manoeuvring 3 landing boxes: front, left and right. Each landing box contains a landing gear, an associated door and the corresponding hydraulic cylinders in charge to move gears and doors. The system is made of a controller (the digital part) and the controlled physical environment (i.e. the 3 landing gear boxes and a pilot interface) which interact via sensors and actuators; the sensors provide to the digital part the information on the state of its physical part; the actuators engage the orders of the controller on the physical part. The physical devices already exist, we will not build them; the challenge deals with the digital control part only (see page 2 of [START_REF] Boniol | The landing gear system case study[END_REF]).

Horizontal Process: Building an Abstract Global Model of the System

We give the main elements resulting from the successive application of the steps proposed in the method (Sect. 2).

Step 1.1: Elicitation and modelling of the interface variables The requirement document listed several triplicated input variables: handle, analogical switch, gear states, doors, • • •. We model (Step 1.3) them with a type T RIPLE = {1, 2, 3} used as an index of the function variables:

GEAR = {FG, LG, RG} analogical_switch ∈ T RIPLE → AnalSW STAT E DOOR = {FD, RD, LD} handle ∈ T RIPLE → HSTAT E HSTAT E = {hDown, hU p} gear_extended ∈ (T RIPLE × GEAR) → BOOL AnalSW STAT E = {openSW, closedSW } door_closed ∈ (T RIPLE × DOOR) → BOOL • • • • • •
The function variable handle ∈ T RIPLE → HSTAT E captures precisely the requirement handle i ∈ {hDown, hU p} with i ∈ {1, 2, 3}. The state variables are the states of the gears, doors, anomalies, etc. They are modelled as follows:

gears_locked_down ∈ BOOL ∧ gears_maneuvering ∈ BOOL ∧ anomaly ∈ BOOL ∧ • • •
The output variables hold the values computed for various electro-valves:

general_EV ∈ BOOL ∧ close_EV ∈ BOOL ∧ open_EV ∈ BOOL ∧ • • •
The lights which indicate the position of the gears and doors to the pilot are described as internal variables: greenLight, orangeLight, redLight. These variables are bound to the output state variable gears_locked_down with an invariant predicate. Another internal variable order is used to record the action of the pilot on the handle.

The LG system is controlled digitally in the normal mode until an anomaly is detected. A permanent failure leads to an emergency mode where the system is controlled analogically. Accordingly the internal boolean variable anomaly is used to denote that an anomaly has been detected or not.

Step 1.2: Elicitation of the global properties of the LG system Most of the normal mode requirements are safety properties. Some identified ones are the following: The first Event-B abstract model resulting from Step 1.3, gathers all the variables of the interface, their related invariants and initialisations. Event-B contexts are used to model the static part with the various sets and definitions that we have introduced.

Step 1.4: The families of events of the abstract model A thorough analysis of the two action sequences (outgoing sequence and retraction sequence) described in the LG system helps us to capture the behaviour of the digital part and to derive the events. We use here state automata to make it clear the interaction between the different components (actions of the pilot, the controller, the orders received by the environment).

In the sense event family we have listed for example the event sense_gear to modify the input variable gear_extended listed above. In the same way we have listed the other events sense_door, etc. Examples of events we have identified for the control events family are: stmlt_general_EV to stimulate the general electro valve, stmlt_door_opening, stmlt_gear_outgoing, stop_stmlt_general_EV, stop_ stmlt_gear_outgoing, etc. Each one modifies its related variable, for instance the event stop_stmlt_gear_outgoing sets the variable extend_EV to FALSE. Examples of event we have classified in the monitor events family are: monitor_ anomaly, monitor_gears_locked_Down, monitor_ gears_maneuvering.

Step 2: Extension of the abstract global model with the event families We achieve many refinement steps, by feature augmentation, to integrate gradually the variables and events related to the physical devices: the sensors, the doors and the gears.

Following Step 2.1, we define the behaviours of physical devices. For instance, the door behaviour is first captured with a state automata; the transitions of the automata are then described as events. For this purpose we use a transition function doorState ∈ DOOR → DSTAT E where DSTAT E = {ClosedLocked,ClosedUnlocked, OpenUnloc-ked} is the enumerated set of the identified door states. The set DOOR contains the three door identifiers. The function doorState is a total function; this captures the requirement that all the three doors are controlled via the state transition.

The starting transition of the door behaviour is enabled by the open_EV order given by the digital part. Therefore there is a synchronisation between the digital part and the motion of the doors. We only give below the description of the starting event Door_openDoor_cl2cu; the other necessary events are similar. where @g0 general_EV = T RUE @g1 order = hDown @g2 ran(handle) = {hDown} @g3 ran(door_closed) = {FALSE} // the three doors are in the open position @g4 ran(door_open) = {T RUE} @next nextOGseq = 3 @gano anomaly = FALSE // no anomaly detected @notretract retract_EV = FALSE then @a1 extend_EV := T RUE @a2 nextOGseq := nextOGseq + sequenceStep end The variable nextOGseq controls the evolution of the outgoing sequence; it indicates in the event guards, the next step in the outgoing sequence. We note that the events in the sense event family anticipate their real future specifications, which are related to the physical part introduced later. When we have introduced the various events families and the related variables, it becomes clear for us that we have the complete control loop.

Following Step 2 the properties (listed in Step 1.2 above) are formalised as first order predicates, integrated into the invariant of the abstract model and, proved along the horizontal refinement. As an example, the requirement R 51 is described as follows.

((open_EV = T RUE ∨ close_EV = T RUE ∨ extend_EV = T RUE ∨ retract_EV = T RUE) ⇒ general_EV = T RUE)
To sum up, the global Event-B abstract model results from a series of refinement of contexts and machines.

Step 3: Dealing with specific properties In this case study, reachability is one of the specific properties. Based on the idea of Lamport's logical clocks [START_REF] Lamport | Time, Clocks, and the Ordering of Events in a Distributed System[END_REF], we implement a technique that captures the reachability requirement R 1 given in page 13 of [START_REF] Boniol | The landing gear system case study[END_REF]. For that purpose, we introduce the notion of control cycle, a period of time during which one can observe several events, especially a chain of events denoting an outgoing sequence or a retraction sequence; a typical control cycle is one starting with an event (downH) which denotes the hDown order and terminating by an event (dcge) which denotes the fact that "the gears are locked down and the doors are seen closed"; similarly, another control cycle is started when the handle triggers an order hU p. A dedicated variable endCycle is used to control the start and the end of each control cycle.

Given a set obsEvents of events (for instance the starting of an outgoing sequence, a door closed, a gear locked in a position, etc.) and a logical clock modelled as a natural number, the occurrences of the events can be ordered by the timestamps given by the clock. In our case two events cannot happen at the same time. We use a partial function ldate ∈ obsEvents → N to record the timestamps of the events. We can compare and reason on the timestamps of any events happening during a sequence and specifically within the specific event sequence called control cycle. An example is as follows.

∀d j.(((d j ∈ N) ∧ (dcge ∈ dom(ldate)) ∧ (d j = ldate(dcge)) ∧(endCycle = T RUE) ∧ d j < llc) ⇒ ∃di.((di ∈ N) ∧ (downH ∈ dom(ldate)) ∧ (di = ldate(downH)) ∧ (di < d j)∧ ∀ii.(ii ∈ N ∧ di ≤ ii ∧ ii < d j ⇒ ldate ∼ [{ii}] = {upH})))
To put into practice in Event-B with Rodin, we defined the set obsEvents in the context of our machines, and the above property is included in the invariants of the abstract model.

Vertical Process: Building the Concrete Parts of the LG System

The vertical process includes several refinements (in Step 4) described below following the proposed method. We end our process with the Step 5. The Step 6 was not performed for the LG case study because only the digital part will be refined with the objective to build the software part. The variables and events which are specific to the behaviour of the physical part are not refined but we keep them in the model in order to preserve animation capabilities. This approach is very pragmatic.

Step 4: Structural refinements of the global abstract model The requirement document details the inner structure of the digital part; it is made of two redundant computing modules. We achieve structural refinement steps to overcome the details of the behaviour of the digital part. a) Introducing the two computing modules with refinements Both modules have the same interface (input and output variables) inherited from the abstract model of the digital part. Each interface variable of a module k (where k ∈ {1, 2}) is inherited from a variable (for instance gear_extended) of the digital part of the abstract model and it is denoted by k_gear_extended(k) where k is an index. An enumerated set CompModule = {1, 2} is used for the indexes. Therefore each interface variable of the computing modules is specified with the following shape:

k_gear_extended ∈ CompModule → ((T RIPLE × GEAR) → BOOL)
The binding between the two modules interface variables and those of the abstract module is achieved via refinements where new variables and related events are introduced. b) Spawning the inputs inside the computing modules with refinements We introduced new events (prefixed with spawn_) to push the value of each input variable (for example handle) at the abstract level, in the corresponding variable (for example k_handle) of each computing module. As the inputs of the modules should be the same, an invariant is defined in each case of variable spawning in order to guarantee the correctness of the binding between the input variable of the digital part and the same input of the computing modules. The following event pattern spawns the variables at the interfaces of the computing modules. We have identified a reusable specification rule: a new event is introduced along with each new k-indexed variable. This event should copy the variable at high level (the digital part) into the indexed variables at the low level. Furthermore, the existing events, whose guards or actions involve the spawned variables, should be refined by extending their guards and actions in order to satisfy the binding between the variables and the associated k-indexed variables. One noticeable feature in this case is that when we have a non-deterministic event of abstract level (as for the value of the sensors), then in the refinement the event should be refined (not extended). This is another reusable specification rule we have identified. c) Merging the outputs of the computing modules with refinements As depicted in Fig. 3, the k-indexed output variables (for example k_extend_EV (1) and k_extend_EV (2)) are merged using a logical OR to set the corresponding variable (for example extend_EV ) at the output of the digital part. Therefore the event that sets the variable should be guarded by the availability of the merged value. As explained before, a binding invariant should be provided for each variable and the related k-indexed variable. Several refinements are used to introduce the appropriate events. d) Specifying the behaviour of the computing modules The two computing modules have the same behaviour. It is made of the events that monitor the system and set accordingly the state output variables and the input variables of the digital part, the events that give orders (control decision) to the physical part through the order output variables. It results in the k-indexed form of the events related to the three categories of the interface variables and the internal variables.

We stopped our construction at this stage. However following the guidelines provided in the method, it remains to perform the decomposition step in the basis of the sense, monitor, control events families (Step 6).

Assessment and Discussion

Coverage Applying the proposed method helped us a lot in mastering the case study. The resulting Event-B model presented in this article covers the main aspects of the landing system: the digital part with modules redundancy, its physical part and their interactions. The model covers mainly the safety properties of the LG system; liveness properties are treated by adapting Lamport's logical clocks [START_REF] Lamport | Time, Clocks, and the Ordering of Events in a Distributed System[END_REF]; but we have not deal with time constraints. Code generation was out of the scope of the current work. Nevertheless the management of huge B models is still tedious, since modifying the models already equipped with a lot of variables an events, at more abstract level requires redoing several steps of modelling, refinements and proving.

Experimentation with Rodin and statistics The Rodin tool is very efficient Statistics on Proof Obligations are given in Tab. 1. From a total of 619 POs, 547 of them were automatically discharged by Rodin and 6 of them were interactively discharged. Most of the POs at the abstract levels were proved. The undischarged POs are related to the structural refinement and specifically they are related to the binding invariants.

Managing very large models requires a rigorous slicing and several small steps of refinements. This is the reason why we have introduced many refinements, but it is still not enough, the slicing should be of finer grain. Moreover a good naming discipline is necessary at each level of the modelling. As far as the ProB animation tool (integrated in Rodin) is concerned, it is very helpful to tune the Event-B models.

Related Works

The state of the art lacks of assistance methods. The four-variable model of software-controlled embedded systems originally proposed by Parnas and Madey has been used successfully in the development of safety-critical applications in various industries. But as mentioned by [START_REF] Patcas | From system requirements to software requirements in the four-variable model[END_REF], the model does not explicitly specify the software requirements, but rather bounds them by specifying the system requirements and the input and output hardware interfaces of the system.

We share the same the motivations with [START_REF] Hudon | Development of Control Systems Guided by Models of Their Environment[END_REF]. However the authors propose a method to synthesise the controller from the environment behaviour. Indeed they introduced the controller and its interface as a solution to the problem of maintaining a desired behaviour of an autonomous system. In our approach the controller is not synthesised to maintain a specific behaviour; it is built simultaneously with the environment according to the given control requirements, but the environment behaviour is less constrained by the controller.

As far as method is concerned in the treatment of the LG system benchmark, all the related B specifications of the LG System are based on refinements. They do not describe a precise methodological process. Often, the authors need about ten refinements to include properties and requirements. The distinction between them is rather the way the refinements are organised rather than on methodological assumptions.

Su and Abrial [START_REF] Su | Aircraft landing gear system: Approaches with event-b to the modeling of an industrial system[END_REF] mentioned that there is no definitive answer for applying some recipe since the question vary from one project to another. They propose a light methodology including three steps: informal requirements, refinement strategy and formal model. They excluded features like redundancy and simplified time constraints. The systematic refinement strategy integrates progressively the devices, which is rather specific to the case study. Accordingly, our method focuses on a more general refinement strategy.

Mammar and Laleau adopted the four-variable model of software-controlled embedded systems originally proposed by Parnas and Madey. They used a series of refinements [START_REF] Mammar | Modeling a landing gear system in event-b[END_REF] first according to a variable classification first (monitor, control, output) then including timing aspects, failure cases and last properties. Mammar and Laleau focused on the control part only. Since it appears to be a logical organisation, a separation of concerns, this ordering delays most of the proof work to the last refinements. It lost modularity and extensibility. R. Banach used Hybrid Event-B to lead his study [START_REF] Banach | The landing gear case study in hybrid event-b[END_REF]; this extension enables one to carry continuous varying behaviours. R. Banach proposed a proof of concept of the language extension rather than a method or a full answer to the case study. However, hybrid-B seems adequate to refine the physical part of our current specification and especially to model time requirements.

Hansen et al. focused on the validation of the case using ProB rather than on the methodology of specifying with B [START_REF] Hansen | Validation of the abz landing gear system using prob[END_REF]. As a matter of fact, the temporal properties are naturally introduced using LTL expressions. Another interesting feature is the ability to visualise the system execution. The counterpart is a simplified specification (no redundancy, no physical part, no failure). The refinements start with physical devices (door, gears, electro-valves), then the output, sensor and controllers are introduced as refinements and finally the general control (switch, valves, lights).

Conclusion

We proposed a method (named Heñcher) to guide step by step the construction of embedded control systems with Event-B. We build on the well-known structure of control systems and on the experiments of several case studies where the Event-B was used and where some methodological guidelines was provided [START_REF] Damchoom | Modelling and Proof of a Tree-Structured File System in Event-B and Rodin[END_REF][START_REF] Damchoom | Applying Event and Machine Decomposition to a Flash-Based Filestore in Event-B[END_REF]. We provide a systematic use of the interface of the controller to build the components of the abstract model of the control system and, also how the features of the control system should be used to guide the successive refinements of the abstract model. A non trivial case study served as illustration and assessment of the proposal.

One flaw of the Event-B top-down approach is the constraint imposed by the evolution of the global abstract model defined before its refinement to the concrete models. This constraint prevents for an incremental model evolution. Indeed, if we miss some features in the abstract state, we will have to reconsider completely the structural refinements. It would be interesting to be able to mix both horizontal and vertical refinements in an incremental view of the design method. In [START_REF] Back | Software Construction by Stepwise Feature Introduction[END_REF] Back, have proposed guidelines for this purpose; an adaptation of this work to Event-B is likely to be interesting.

The reuse of existing independent models, with a bottom-up approach, would be interesting for managing large Event-B models. A typical example is the composition of existing models to build a given abstract model where each part can be modified and refined separately.

We plan to develop an assistance tool to help the user with various patterns, in the form of Event-B machines derived from the interface variables which will be extracted from a sketched graphical view of its control system (as in Fig. 2).

Fig. 1 .

 1 Fig. 1. Synoptic structure of the Event-B models of the construction

Fig. 3 .

 3 Fig. 3. Modules redundancy

Fig. 4 .

 4 Fig. 4. Global architecture of the LG system

R

  21 We can not observe a retraction sequence (consequence of the order hU p) if the handle is down. Using the enumerated set HSTAT E which permits only one value from two for the variable order. R 31 The gears outgoing event occurs if doors are open locked. R 41 Opening and closing doors electro-valve are not stimulated simultaneously. R 51 It is not possible to stimulate the manoeuvring EV (opening, closure, outgoing or retraction) without stimulating the general EV.

  event Door_openDoor_cl2cu /* Door's Behaviour (for the three doors). The first transition of the Door Automata */ where @g1 open_EV = T RUE // all the doors Electro Valves are on @g2 ran(doorState) = {notOpenLocked} then @a1 doorState := DOOR × {notOpenNotLocked} // door is being opened end The following event describes an event of the control event family. event stmlt_gear_outgoing /* stimulate gear outgoing electro valve once the three doors are in the open position */

  event spawn_handleDown // spawn handleDown within the k CompModules where @g1 ran(handle) = {hDown} then @a1 k_handle := {1 → (T RIPLE × (ran(handle))), 2 → (T RIPLE × (ran(handle)))} end

Table 1 .

 1 Statistics of PO generated and proved with Rodin

	for proving the Event-			Total Auto Manual Review. Undis.
	B models; a very high percentage (∼ 90%) of proof obligations was	LandingSys5 Landing_DP_Ctx LandingSysDP_A	619 547 6 Abstract model 0 0 0 115 114 1	0 0 0	66 0 0
	automatically discharged. The specifications are	LandingSysDP_SWITCH_A LandingSysDP_DOOR_A LandingSysDP_DOOR_GEAR_A	5 42 42 0 3 0 79 79 0	0 0 0	2 0 0
	available online 3 . The current version is par-tial as we focus on rep-	LandingSysDP_DOOR_GEAR_TIME_A Models of the vertical refinement 2 2 LandingSysDP_DGT_R1_In 52 50 0 0 LandingSysDP_DGT_R2_INOUT 56 56 0	0 0 0	0 2 0
	resentative events in-stead of being exhaus-	LandingSysDP_DGT_R3_INOUTDOOR LandingSysDP_DGT_R3INOUTDOORGEAR 140 120 0 128 81 5	0 0	42 20
	tive.

http://uk.mathworks.com/products/control/

http://www.scilab.org