B Tolbi

H Tebbikh

H Alla

Fault-tolerant continuous flow systems modelling

Keywords: Continuous flow systems, Elementary Hybrid Petri Nets, faults modelling, Hybrid Automata, Timed Transition Systems, Timed bisimilarity

Fault-tolerant continuous flow systems modelling

B Tolbi, H Tebbikh, Hassane Alla

Introduction

Hybrid dynamic systems (HDSs) are dynamic systems integrating explicitly and simultaneously continuous and discrete event systems, which require for their description the use of a continuous model, and a discrete event model [START_REF] Pettersson | Hybrid modelling focused on hybrid petri nets[END_REF]. In this paper, a particular class of hybrid dynamic systems, which are continuous flow systems is considered, for example transportation systems, production systems, communications systems...etc. This class involves hybrid systems, which are said positive, i.e. all the state variables take only positive values. In the literature, a particular effort was given to study this class of hybrid dynamic systems. Indeed, it is sufficiently rich to allow a realistic modelling of many actual problems. This assumption allows is introduced to represent independently faults. Faults modelling on discrete part will not be presented in this paper, only on continuous part. Secondly, a simple modelling technique will be done to describe system without faults. Finally, a composition of independent models will be presented to get the global model. A powerful HPN model is necessary to show the faults consequences on the continuous part, and for behaviour analysis, HA model is also necessary. For this, formal semantics of EHPNs and HA will be given in terms of TTS and the analysis power of HA will be coupled to the modelling power of EHPNs using a structural translation from EHPN to HA that preserves the semantics (in the sense of timed bisimilarity) of both models. The technique of bisimilarity will be used to prove the correctness of our translation procedure by giving the proof of this one.

Outline of the paper. Section 2 introduces faults modelling using the general model of HPNs, in this section the easy way of faults modelling will be presented. Section 3 explains our proposal procedure to translate the sub-HPN named EHPNs into HA. In this Section, also the way how to prove the bisimilarity between EHPNs and HA as TTS formalism will be showed. Section 4 provides a case study of our translation method. Finally, our ongoing work and perspectives will be given in Section 5.

Faults modelling

Let us note that a system fault corresponds to a dysfunction state, while a failure or a fault source of failure is an event that can lead to a state of dysfunction. In the context of HDSs, the occurrence of fault is, also, the passage towards a dysfunction state. This passage can be modelled by a transition on a fault, where one considers an event-based modelling [START_REF] Sampath | Diagnosability of discrete-event systems[END_REF]. Other authors consider a state-based modelling [START_REF] Zad | Fault diagnosis in finite-state automata and timed discreteevent systems[END_REF], an approach combining the advantages of state and event based modelling is proposed in [START_REF] Sayed Mouchaweh | Decentralized diagnosis by boolean discrete event system model: Application on manufacturing systems[END_REF]. A fault can also be represented as an execution of a given supervision pattern, which is a temporal property related to the occurrence of a set of trajectories/events that must be diagnosed [START_REF] Jeron | Supervision patterns in discrete event systems diagnosis[END_REF]. These techniques use models including faulty behaviours. There are also approaches, which use fault-free models; they are based on comparing the systems outputs with the models nominal outputs. The fault-free modelling approach proposed by [START_REF] Pandalai | Template languages for fault monitoring of timed discrete event processes[END_REF] uses condition templates to determine if the system generates events in the right order or within the given time delays. In (Sayed-Mouchaweh, 2012) expert knowledge is associated with condition templates to identify the faults related to missing or unexpected events, and progressive monitoring is used to reduce the set of fault candidates after the occurrence of new observable events. Another practical fault-free modelling approach for fault diagnosis of manufacturing systems has been proposed in [START_REF] Roth | The concept of residuals for fault localization in discrete event systems[END_REF].

Hybrid Petri nets

In this formalism, the firing of a continuous transition describes the material flow, while the firing of a discrete transition models the occurrence of an event that can for example change firing speeds of the continuous transitions.

The work presented in this paper is interested, and restricted, only to EHPNs formalism, since there is a decoupling between the discrete and the continuous parts (one part may influence the behaviour of the other one, but there is no transformation of discrete marking into continuous marking or vice-versa). They combine a time PN and a constant speed continuous PN (CCPN)(R. David & Alla, 2010). A formal presentation of EHPNs will be given in Definition 3.

This paper is based on some of techniques which are presented at the beginning of this section, but here is not needed to rebuild automaton model to new or build both automaton model (model without faults and model with faults) [START_REF] Derbel | Online Diagnosis of Systems with Rectangular Hybrid Automata Models[END_REF]. This work is aimed to build the system complete model (model with faults), by modelling the system in normal functioning, and an independent modelling of faults. The model that will be used is HPN. In this model, the faults are considered as unobservable events, this leads to decompose the sets of events in two subsets: the set of observable events and the set of unobservable events (Σ o and Σ u ). If the unobservable event is associated with a discrete transition, then this transition is said to be unobservable. This work can be extended into general HPN, but the translation algorithm will be difficult, mostly the continuous-discrete decoupling will be losed and therefore the structural character of translation. However, this problem could be studied as a future perspective.

Definition 1: A hybrid PN model in normal functioning is an 8-uplet N =(P , T , P re, P ost, h, Σ, I, V , M 0 ), such that: 1. P ={P 1 , P 2 , ..., P n } is a finite set of n places.

P = P C ∪ P D ; 2. T ={T 1 , T 2 , ..., T m } is a finite set of m transitions. T = T C ∪ T D ;
• T D is a finite set of discrete transitions; • T C is a finite set of continuous transitions; 3. P re : P × T →N and P ost : P × T →N are the backward and forward incidence mappings; These mappings are such that: ∀(P i , T j ) ∈ P D × P C , pre(P i , T j ) = post(P i , T j ); 4. h : P ∪ T → {C, D }, defines the set of continuous nodes (h(x) = C) and discrete nodes (h(x) = D);

5. Σ is a set of events; Σ=Σ o ∪ Σ u • Σ o is a sub-set of observable events; • Σ u is a sub-set of unobservable events; • Σ → Q + × Q + associates to each event σ j , an occurrence interval d j = [α j β j ].
6. I : T D → Σ associates an event to each D-transition T j ; 7. V :T C → R + associates a maximal firing speed V j to each C-transition T j . 8. M 0 is the initial marking, such that M =(m C , m D ) T m C and m D are the discrete and the continuous markings.

In order to illustrate the HPN modelling let us consider the following example.

Example 1: Let consider a road section that can tolerate a maximum number of cars 150 (assuming that average distance between two cars is L=4m). Assuming that the section has an entry ramp, the two section entries are regulated by a traffic light, the green light turns on in the first and the second entry, respectively after d 1 = [20 25] and d 2 = [25 30]. Cars enter the section with an average speed of 30 km/h from both entries and leave it with an average speed of 48 km/h. The example is given in figure 1(a). gives the HPN in normal functioning of the road section. The continuous part is represented by a double line, the continuous place P 5 represents the section and it is associated by its maximum number of cars (150) which are represented by tokens. The continuous transitions are associated by their maximal firing speeds "for example: V 5 = (48km/h × 1000)/(3600 × 4m) ". The discrete part is represented by simple line where the transitions are associated by occurrence durations (d 1 and d 2 ) of events (σ 1 and σ 2 ) which correspond to the occurrence durations of green lights on the two section entries.

Our objective is to describe independently faults behaviour, representing just the faults. Faults modelling on the discrete part is not discussed here as it is widely studied. Only faults that have a consequence on the continuous part are presented. Fault occurring on the continuous aspect may be modelled by a continuous transition controlled by the occurrence of an unobservable event associated to a discrete transition. This is defined bellow.

Definition 2: A single fault model is a sub-hybrid PN F which contains: 1. P 1 ,P 2 ,two discrete places; 2. T 1 a discrete transition associated with fault σ f (it is an unobservable event); 3. T 2 a continuous transition with constant firing speed V 2 ; 4. P re f : P 2 × T 2 →N and P ost f : P 2 × T 2 →N are the backward and forward incidence mappings of fault such that:

P re f (P 2 , T 2 ) = P ost f (P 2 , T 2 ) = 1; 5. M (P 1 ) = 1, M (P 2 ) = 0 is the initial marking.
Let us now complete example 1. A faulty situation of traffic can be observed if the green light is turned on, in both entries; the passage to this faulty state is due to the occurrence of the fault event (σ f ). The HPN model of this fault is given in Figure 2(a). It models the fault and its influence on the continuous behaviour. After occurrence of fault event, P 2 is marked and continuous transition T 2 is enabled at its maximal speed V 2 . This gives an additive flow to place P 3 , this behaviour is shown by the evolution graph defined in Figure 2(b). This evolution graph contains two states. The change from one state to another is assured at the fault occurrence, and that is reflected with firing speed modification of the continuous transition that becomes strictly positive. This change in dynamics will be useful in the diagnoser synthesis in our future work. Remark 1: It is possible to have a subtractive flow at the occurrence of the fault, for that we must reverse the arc that links transition T 2 and place P 3 .

To get the global model, it is needed to perform a composition of both models (normal functioning and fault models), this composition is structural and independent from the initial marking. And then the model with faults is translated directly to get an automaton model ready for analysis, even in the presence of faults. The global model is defined bellow. Definition 3: A global HPN model G is a structure obtained from the composition of the models N and F (defined respectively in Definition 1 and 2). It is an 8-uplet G =(P G , T G , P re g , P ost g , h, Σ g , I, V , M 0 g), such that: 1. P G ={P 1 , P 2 , ..., P n+2 } is a global finite set of n + 2 places (P 1 and P 2 are the two discrete places of F); 2. T G ={T 1 , T 2 , ..., T m+2 } is a global finite set of m + 2 transitions (T 1 and T 2 is a fault discrete and continuous transition respectively of F); 3. P re g , P ost g are the global backward and global forward incidence mappings; 4. h : P ∪ T → {C, D }, defines the set of continuous nodes (h(x) = C) and discrete nodes (h(x) = D); 5. Σ g is the set of all events; Σ g =Σ o ∪ Σ u + {σ f }; 6. I : T D → Σ g associates an event to each D-transition T j ; 7. V :T C → R + associates a maximal firing speed V j to each C-transition T j . 8.

M 0g =[1 0 M 0 ] T is the global initial marking, (M 0 is the initial marking of F).
Remark 2:

• For multiplicative faults, the structure that corresponds to each type of fault must be duplicated; • A fault that corresponds to a structure change leads to a new HPN (this is not formalized in our paper). Figure 3 gives the global model of the road section.

Hybrid automata

To get the powerful model, it is necessary to couple the analysis power of HA to the modelling power of EHPN by preceding a structural translation from EHPN to HA which is defined bellow.

Definition 4: A hybrid automata is a 6-uplet HA=(Loc, x, E, δ, F , inv), such that: 1. Loc is a finite set of locations; 2. x is the continuous state space. x C is the vector of real-valued variables modelling the continuous places marking and x D is the vector of clocks corresponding to enabled transitions. A valuation is a function that assigns a real-valued v(x) ∈ R to each variable x i ∈ x.

3. E is a set of events; 4. δ is a finite set of transitions, each transition is a quintuple T =(q, a, g, init, q ) such that:

• q ∈ Loc is the source location;

• a ∈ E is the event associated to the firing of T ;

• g is the transition guard, it is a linear predicate on x; a transition can be fired whenever its guard is satisfied; • init is a reset function that affects a linear expression to variables of x when taking the corresponding transition; • q ∈ Loc is the target location; 5. F is a function that assigns to each location a continuous linear vector field on x. While in discrete location q, the continuous variables m i ∈ m C evolve according to a differential equation of the form ṁi = B i , where B i ∈ R + is a the dynamic balance of the continuous place P i , and the clocks t j ∈ x D evolves according to the differential equation ṫj = 1; 6. inv is a function that affects to each location q, a linear predicate inv(q) that must be satisfied by the continuous variables in order to stay in the location q. Most of works developed on HDSs return to a generalization of results coming from continuous systems (component approach), where few switching are considered. These approaches extend results of the continuous domain, such as unified models [START_REF] Branicky | Studies in hybrid systems: Modeling, analysis and control[END_REF], bond-graphs with switches [START_REF] Buisson | Descriptor systems for the knowledge modelling and simulation of hybrid physical systems[END_REF]. The most considered faults are sensors faults, such as sensor offset or sensor stuck-off/stuck-on, which represent discrepancies between the measured and real values of system variables.

The approaches coming from DESs (Discrete Event Systems) consider the global behaviour of a process (system approach), where breaking dynamics are modelled. The main models are time Petri nets [START_REF] Berthomieu | Ifip congress series[END_REF], timed automata [START_REF] Alur | A theory of timed automata[END_REF], hybrid automata [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF], hybrid Petri nets [START_REF] Bail | Hybrid petri nets[END_REF], batch Petri nets [START_REF] Demongodin | Generalised batches petri net: Hybrid model for high speed systems with variable delays[END_REF]. In these models the discrete part is complex while the continuous one is often simpler. The faults are boolean and lead to on/off behaviours. The major interest of our work is to consider the faults in the continuous part in an elegant way. Indeed, by modifying the balance marking of a continuous place in a HPN, it is possible to model a breakdown of a sensor or an additional flow due to a leakage or a feeding. This physical aspect of the fault description highlights the practical interest of this modelling.

From elementary hybrid Petri nets to hybrid automata

In this Section, hybrid automata are built from EHPNs so that the behaviours of both are in a one-to-one correspondence using the proposed algorithm. The HA resulting from the translation and the EHPN model will be proved to be timed bisimilar. To achieve this goal, formal semantics of HA and EHPNs are given in terms of TTS. Therefore, the bisimilarity between a TTS of an EHPN and a TTS of a HA is demonstrated. A bisimilarity between two transition systems ensures that any action of one of the systems can be simulated by the other (equivalent behaviours).

Structural translation

Our objective is to realize a structural translation from the sub-model EHPNs to HA. An EHPN is characterized by the following structure linking a continuous place to a discrete transition or a continuous transition to a discrete place (there is no transformation of tokens). Each marking in time PN will correspond to a CCPN configuration. An EHPN model can therefore be studied in a hierarchical manner, firstly the discrete part is considered, and then for each reachable discrete marking, the continuous configuration is considered. The translation process of this work is given in a informal way by the following steps:

Step 1: Build the reachable marking graph of subjacent autonomous PN by isolating the time PN part from continuous part (CCPN) of hybrid model;

Step 2: Associate respectively the dynamics1 and the guards to each location and transition according to the continuous part configuration of hybrid model, where each reachable discrete marking correspond to a CCPN configuration and null flow vectors must be suppressed.

Step 3: Build the final HA, where the value of the entry guard determines the target location whose downstream configuration of each location. The unreachable locations must be deleted.

The main part of our work concerns the structural modelling of faults by EHPN tool, which is a powerful tool for modelling. This contribution is not addressed in the paper [START_REF] Ghomri | Synthése de contrleur de systèmes hybrides á flux continu par réseaux de petri hybrides[END_REF], only the translation part is similar. Our translation procedure shows the effect of continuous place on discrete one and separates the dynamics of the discrete part from the continuous part at a certain level (there are no marking transformation, from discrete to continuous or vice-versa). In our case, the structure given in Figure 5 is added and for the first time, a mathematic proof for our proposed procedure is presented. This last will be detailed at the end of Section 3.3.

Semantics of the two basic models

TTS is used to describe the set of states and transitions between these states [START_REF] Larsen | Model-checking for real-time systems[END_REF]. Two types of transitions are possible in a TTS: a continuous transition or time transition that describes the flow of time or a continuous evolution and a discrete transition describing the evolution following a discrete action (event). If an initial state is isolated in a TTS, the different paths starting from this state represent the different possible evolutions for modelled system. In this work, the emptying of a C-place (C1-event (R. David & Alla, 2010)) is considered like discrete transition. The semantics of HA and EHPNs are defined bellow.

Definition 5: Hybrid Automaton Semantics. The semantics of a hybrid automaton HA is defined as timed transition system S HA = (S A , S A 0 , E, -→ A ) so that: • S A = (q, x) is the set of states; • S A 0 = (q 0 , x 0 ) is the initial state; • E is the set of actions (events);

• -→ A ∈ S A × (E ∪ R + ) × S A
is the set of transitions that can be discrete or continuous, such that: Discrete transition: firing transition T j in the hybrid automaton.

(q, x)

Tj -→ A (q , x ) if:   
∃(q, a, g, init, q ) ∈ δ so that: g(x) = true, inv(q)(x) = true and inv(q )(x ) = true Continuous transition: Time elapsing. (q, x)

d -→ A (q , x ) if:    q = q x = F (d) ∀ 0 ≤ d ≤ d, inv(q)(x + d ) = true
Definition 6: Elementary Hybrid PN Semantics. The semantic of an EHPN is a timed transition system S EHP N = (S E , S E 0 , Σ g , -→ E ) so that:

• S E = (m D , m C , t, v) is the set of states; • S E 0 = (m D0 , m C0 , 0, v 0 ) is the initial state; • Σ g is the set of actions (events), Σ g = T D ∪ {m i = 0} • -→ E ∈ S E × (Σ g ∪ R + ) × S E
is a discrete or a continuous transitions, such that:

The discrete transition (m D , m C , t, v)

a -→ E (m D , m C , t , v
) with a ∈ Σ g is subdivided into two forms:

• a ∈ T D (the firing of a transition T j in the T-time PN) is the first form of discrete transitions: t j is the time associated to the transition T j ; t k is the new time associated to the transitions newly enabled (T k ) after firing of

T j . (m D , m C , t, v) a -→ E (m D , m C , t , v ) if:            M ≥ (., T j ), T j ∈ T D2 and m D = m D -P re(., T j ) + P ost(., T j ), m C = m C α(T j ) ≤ t j ≤ β(T j ) if ↑ enabled(T k , M, T j ) : t k = 0, else t k = t k ∀(P k ∈ T • i ), T j ∈ T D2 {T i , T i+1 } ∈ T C , and v i = 0 ∀(P j ∈ T • j ), T j ∈ T D1 , and v i+1 = V i+1 2 • a ∈ {m i = 0} (the emptying of a continuous place P i ) is the he second form of discrete transitions: (m D , m C , t, v) a -→ E (m D , m C , t , v ) if: v = f (m C , t) m i = 0 • The continuous transition is the time elapsing: (m D , m C , t, v) d -→ E (m D , m C , t , v ) if:        m C = m C + B.d t = t + d ∀T k ∈ T D , if m ≥ P re(., T k ) then t k ≤ β(T k ) ∀P i ∈ P C , if m i > 0 then m i > 0
After construction of timed transition systems of HA and EHPN, the bisimilarity between semantics of the two models will be demonstrated subsequently.

Timed similarity and bisimilarity

In this sub-section, the proof of the correctness of the translation algorithm is given. This will guaranty that the translation is correct and terminates. It will be proved that the semantics of an EHPN and its translated HA are timed bisimilar. The timed similarity and bisimilarity relations between both semantics are given bellow to introduce the proof.

Definition 7: Timed Similarity and Bisimilarity relation between EHPN and HA.

• Let S EHP N = (S E , S E 0 , Σ g , -→ E ) be the TTS of EHPN and S HA = (S A , S A 0 , E, -→ A ) be the TTS of HA. Let Φ be a binary relation on S E × S A . We write sΦs for (s, s ) ∈ Φ. Φ is a timed relation from S EHP N by S HA if the following assertions are checked:

     if s 1 ∈ S E 0 , then ∃s 2 ∈ S A 0 such that s 1 Φs 2 ; if s 1 δ -→ E s 1 , with δ ∈ R + and s 1 Φs 2 then ∃s 2 δ -→ A s 2 such that s 1 Φs 2 ; if s 1 a -→ E s 1 , with a ∈ Σ g and s 1 Φs 2 then ∃s 2 a -→ A s 2 with a ∈ E such that s 1 Φs 2 . • ∀(S EHP N ΦS HA )∧(S HA Φ -1 S EHP N ), then S EHP N ≈ S HA . ≈ is a bisimilarity relation between S EHP N and S HA . • ∀S 1 = (q 1 , x 1 ) ∈ S E is a state of S EHP N , S 2 = (q 2 , x 2 ) ∈ S A is a state of S HA and S 1 ΦS 2 ⇔ m D = (q 2 ) and x 1 = x 2
, where is a function associating to each discrete marking a given location.

Remark 3: Relation Φ expresses the time evolution given by the two identical dynamics of the EHPN and the HA. They are calculated from the firing conditions of both models: balance marking for the continuous evolution and time interval for the discrete one.

Theorem 1: The TTS of EHPN (S EHP N ) and The TTS of HA (S HA ) are timed bisimilar, S EHP N ≈ S HA . (≈ is a timed bisimilarity relation).

Proof. ∀S 1 = (q 1 , x 1 ) ∈ S E and S 2 = (q 2 , x 2 ) ∈ S A such that S 1 ΦS 2 , We write S 1 ΦS 1 for (S 1 , S 1 ) ∈ Φ. S EHP N ΦS HA if the following assertions are checked :

if S 1 d∈R + -→ E S 1 , then ∃S 2 d∈R + -→ A S 2 such that S 1 ΦS 2 ; if S 1 a∈Σg -→ E S 1 , then ∃S 2 a∈Σg -→ A S 2 such that S 1 ΦS 2 ;
Continuous transitions: According to the semantics of EHPN (Definition 6) the continuous transition d∈R + -→ E is the time elapsing. t and m C check the following:

t = t + d; m C = m C + B.d and m D = m D ⇒ q 1 = q 1 , x 1 = F (d).
According to our algorithm the transition d∈R + -→ E is translated to a transition d = [α j β j ] between sub-location (Algorithm1, line 9) or between location (Algorithm1, line 29), so only dynamics that will be changed during this time, then

x 2 = F (d) and ∀m D = m D ⇔ (q 2 ) = (q 2 ) ⇒ q 2 = q 2 . S 1 ΦS 2 , (q 1 = q 1 ) ∧ (q 2 = q 2 ) ⇒ q 1 = q 2 x 1 = F (d) and

x 2 = F (d) ⇒ x 1 = x 2 ⇒ S 1 ΦS 2 (1)
Discrete transitions: According to the semantics of EHPN, the discrete transition -→ E : if M ≥ P re(., T j ) ⇒ M = g(x 1 ) = true, and if T j is fired, then m C = m C ; m D = m D -P re(., T j ) + P ost(., T j ) ⇒ inv(q 2 )(x 1 ) = true and inv(q 2 )(x 1 ) = true According to our algorithm (Algorithm 1, line 5-9), T j correspond to the transition with guard m i = S ⇒ m i = g(x 2 ) if g(x 2 ) = true and inv(q 2 )(x 2 ) = true, then inv(q 2 )(x 2 ) = true. inv(q 2 )(x 2 ) = true and g(q 2 )(x 2 ) = true such that S 1 ΦS 2 then

x 1 = x 2 ⇒ S 1 ΦS 2 (2)
b) the emptying of a continuous place P i is the discrete transition a∈{mi=0}

-→ E : It is considered as a discrete event of type (C1-event), If m i = 0, g(x 1 ) = 0 ⇒ m i = g(x 1 ) = true ⇒ m D = m D , and m C = m C -P re(., P i ) then g(q 2 )(x 2 ) = true this emptying is translated in this work (Algorithm1, line 38) to a transition with a guard

m i = 0 ⇒ g(x 2 ) = 0 ⇒ m i = g(x 2 ). If g(x 2 ) = true and inv(q 2 )(x 2 ) = true, then g(q 2 )(x 2 = true), such that S 1 ΦS 2 then x 1 = x 2 ⇒ S 1 ΦS 2 (3)
From ( 1), ( 2) and (3) we write: for all S 1 and S 2 on relation,

S EHP N ΦS HA (4)
With the same method, according to the semantics of HA (Definition 5) and according to our algorithm, it has been proved that:

S HA Φ -1 S EHP N (5)
From ( 4) and ( 5), it has been deduced that S EHP N and S HA are timed bisimilar for all S 1 and S 2 on relation. We write:

S EHP N ≈ S HA (6)
From ( 6) it has been concluded that the semantics of EHPN are timed bisimilar to semantics of HA, this means that a behavioral equivalence between EHPN and HA has been provided. This last ensures the correctness of our translation procedure presented in this work.

Case study

The liquid heating system [START_REF] Derbel | Diagnostic à base de modèles des systèmes temporises et d'une sous-classe de systèmes dynamiques hybride[END_REF], illustrated in Figure 6(a), is composed of two valves V 1 and V 2 , a tank, a resistor, a thermostat and two level sensors: c 1 and c 2 monitor respectively the maximal and the minimal levels. The liquid to be heated is introduced through the valve V 1 with a flow rate 5 v.u./t.u. (volume units per time unit). When the liquid level reaches the maximum level 500 v.u, the sensor c 1 generates a notification event to controller for the closing of valve V 1 . Thereafter, the liquid is heated for a period of 40t.u. Then the liquid is discharged through the valve V 2 with a flow rate 8 v.u./t.u. until the tank is empty. At this time, the sensor c 2 generates a notification event; therefore the controller starts a new heating cycle. Assuming that the functioning of this system can be affected by two faults:

• The first fault corresponds to existence of a leak in the tank with a flow rate 0.5 v.u./t.u.

This fault is represented by the unobservable event σ 1 ; • The second fault corresponds to the blocking valve V 1 in the closed position. This fault is represented by the unobservable event σ 2 . s 1 and s 2 are observable events and they correspond to notifications generated, respectively, by the sensors c 1 and c 2 . The event u is introduced to allow the representation of an autonomous transition. At the initial time t = 0, the tank is empty.

The elementary HPN model in normal functioning is given in Figure 6(b). With the transitions T 4 and T 5 are associated maximal firing speeds V 4 = 5 and V 5 = 8; and with the transitions T 1 , T 2 and T 3 are associated durations d 1 = 0, d 2 = 40 and d 3 = 0 of events s 1 , u and s 2 successively, knowing that d 1 and d 3 are a synchronization durations only. The occurring of leak fault of continuous aspect is modelled by the continuous transition T 2 attached to the continuous place P 3 with an arc from the place to the transition T 2 . This transition is controlled by the occurrence of the unobservable event σ 1 associated to the discrete transition T 1 (Figure 6(c)).

The blocking valve fault is modelled by the discrete transition T 1 attached to the discrete place P 2 with an arc from the place to the transition T 1 . This transition is controlled by the occurrence of the unobservable event σ 2 associated to the discrete transition T 1 , this fault can occurred at any time during filling process (Figure 6(d)). To get the global model given in Figure 7, it is necessary to perform a structural composition of normal functioning model and fault models. The model with faults is translated directly using the presented algorithm, the final HA is obtained as follows:

Step 1: Build the reachable marking graph of subjacent autonomous PN (Figure 8) by isolating the time PN part from continuous part (CCPN) of hybrid model;

Step 2: Associate respectively the dynamics and the guards to each location and transition (Figure 10) according to the continuous part configuration, where each reachable discrete marking correspond to a CCPN configuration (Figure 9) and null flow vectors are suppressed;

Step 3: The value of the entry guard determines the target location whose downstream configuration of each location, the unreachable locations must be deleted (this is not the case in our example). The final HA is given in Figure 11.

By applying this translation to the case study, a HA model bisimilar to EHPN model is obtained, then all the properties on the HA remain true for the EHPN. Moreover, the ease of faults modelling that affect the continuous part of heating system, and the change of system dynamics caused by these faults are shown. A structural modelling has been done and the same model presented in [START_REF] Derbel | Diagnostic à base de modèles des systèmes temporises et d'une sous-classe de systèmes dynamiques hybride[END_REF] is obtained, with an automatic and systematic way without imagining faults behaviour. In addition, the dynamics gap appears clearly on locations of the faulty behaviour. This gap can be observed on locations L 5 , L 6 and L 7 and it could be remarked on the dynamic value which becomes respectively 4.5 instead of 5, -0.5 instead of 0 and -8.5 instead of -8 (subtractive flow). This corresponds to the basic idea which will be used for the diagnoser synthesis.

Conclusion

In this paper, EHPNs is used to describe the evolution in normal operation of a type of HDSs called continuous flow systems. A new faults modelling method based on fault structure is proposed to model faults that can be occurred in this type of systems. This has allowed us to get the global model by performing a structural composition of both models, the physical system model in normal functioning with the faults model. No dynamics synchronization is used. Then, the model with faults is translated directly to get an automaton model ready for analysis using a new translation approach. Timed bisimilarity property between EHPN and HA as TTS models has been proved to ensure the correctness of our translation algorithm. A liquid heating system has illustrated the proposed approach. As future work, a diagnoser synthesis of this type of systems modelled by EHPNs could be developed.

  Figure 1. (a). Road section. (b) Hybrid PN of the normal functioning.

Figure

  Figure1(b) gives the HPN in normal functioning of the road section. The continuous part is represented by a double line, the continuous place P 5 represents the section and it is associated by its maximum number of cars (150) which are represented by tokens. The continuous transitions are associated by their maximal firing speeds "for example: V 5 = (48km/h × 1000)/(3600 × 4m) ". The discrete part is represented by simple line where the transitions are associated by occurrence durations (d 1 and d 2 ) of events (σ 1 and σ 2 ) which correspond to the occurrence durations of green lights on the two section entries.Our objective is to describe independently faults behaviour, representing just the faults. Faults modelling on the discrete part is not discussed here as it is widely studied. Only faults that have a consequence on the continuous part are presented. Fault occurring on the continuous aspect may be modelled by a continuous transition controlled by the occurrence of an unobservable event associated to a discrete transition. This is defined bellow.

  Figure 2. (a) Hybrid PN of the fault model of the road section. (b) Corresponding evolution graph.

Figure 3 .

 3 Figure 3. Global HPN model of the road section.

Figure 4 .

 4 Figure 4. HA of the road section.

Figure 4

 4 Figure4gives the HA of the road section, it is obtained thanks to Algorithm 1 given in Appendix 1. In each node, the time derivative of the continuous place corresponds to its balance marking (input flow minus output flow).

Figure 5 .

 5 Figure 5. Mandatory structure in an EHPN.

  firing of a transition T j in the T-time PN is the discrete transition a∈T D

  Figure 6. (a). Liquid heating system. (b) Elementary HPN model in normal functioning. (c) First fault model.(d) Second fault model.

Figure 7 .

 7 Figure 7. Global elementary HPN model.

Figure 8 .

 8 Figure 8. The reachable marking graph of subjacent autonomous PN.

Figure 9 .

 9 Figure 9. CCPN configuration of each reachable discrete marking.

Figure 10 .

 10 Figure 10. Initial hybrid automaton.

Figure 11 .

 11 Figure 11. The final hybrid automaton of the liquid heating system.

the marking balance corresponds to the difference between the input and the output flow for each continuous place Pi of a CCPN.

To calculate instantaneous firing speeds see (R.David & Alla, 2010).

The instantaneous firing speeds are calculated according to an algorithm presented in details in (R.David & Alla, 2010).

Acknowledgements

We warmly thank Latéfa Ghomri Doctor at University of Tlemcen, Algeria, and Yamen El Touati Doctor at university of Manouba, Tunisia for their support and guidance in Section 3.3 of this paper.

Appendix 1: Translation algorithm.

The translation algorithm of an EHPN into a HA of this work presented in Section 3.1 is detailed on the following algorithm:

Algorithm 1 Translation algorithm Initialization ; Isolate the timed PN part from continuous part (CCPN) of hybrid model; Build the reachable marking graph of Subjacent autonomous PN;

for

Burst each location which corresponds to T j ∈ T D2 into two sub-locations L k1 and L k2 ; Associate the time clock to location L k2 with initialization at his input; Connect L k1 and L k2 for each location with a transition and associate him the guard m i = S ; Create an output transition for each L k2 and associate him the guard


Hervé Leblond
email: herve.leblond@univ-angers.fr

Philippe Grelu

Dumitru Mihalache

Houria Triki

Few-cycle solitons in supercontinuum generation dynamics

We review several propagation models that do not rely on the slowly-varying-envelope approximation (SVEA), and can thus be considered as fundamental models addressing the formation and propagation of few-cycle pulsed field structures and solitary waves arising in the course of intense ultrashort optical pulse evolution in nonlinear media and beyond octave-bandwidth optical spectrum broadening. These generic models are: the modified-Korteweg-de-Vries (mKdV), the sine-Gordon (sG), and the mixed mKdV-sG equations. To include wave polarization dynamics, the vector extensions of both mKdV and sG equations are introduced. Multi-octave-spanning supercontinuum generation and few-cycle soliton structures are highlighted from numerical simulations.

Introduction

At the turning of the last century, the improvement of titanium-sapphire laser cavities mode locked through the Kerr-lensing effect, including for instance chirped mirrors for precise dispersion compensation, allowed an efficient exploitation of the very large gain medium bandwidth with the direct generation of amazingly short laser pulses, with a duration down to a few optical cycles. As one of the first realizations, Ref. [1] reported in 2001 the generation of 5 fs-long pulses with a spectrum extending from 600 to 1200 nm, namely, pulses about two optical cycles long spanning over an octave across the central 800-nm wavelength, directly from the laser oscillator. For such intense ultrashort pulses, even diluted media such as the surrounding air can manifest as strongly dispersive and nonlinear ones.

The question was then, how the propagation of such pulses could be modeled. Popular ultrashort propagation models had been centered around the slowly-varyingenvelope approximation (SVEA) that is fundamentally a perturbative approach valid when the pulse duration is much longer than the optical cycle or, equivalently, when the spectral bandwidth is much smaller than the central frequency, which can thus be considered as a carrier frequency. At the lowest orders of dispersion and nonlinearity, propagation in cubic nonlinear media can be effectively described by the nonlinear Schrödinger (NLS) equation, which can be accurate for pulsewidths down to about ten optical cycles at moderate optical intensities. In addition, the NLS equation has served as a central model in the investigation of envelope soliton formation and propagation. With an increase in the optical intensity, additional nonlinear effects need to be included, such as Raman and self-steepening effects. As the optical bandwidth increases, higher-order dispersion terms cannot be neglected. Thus, within the frame of the SVEA, the answer to shorter and higher-intensity pulse modeling implies an increasing number of terms added to the NLS equation, which is the so-called generalized NLS equation (GNLSE) approach [2]. The GNLSE implies numerical integration, and has indeed provided successful modeling of supercontinuum generation experiments in well-characterized optical media such as silica glass fibers [3]. Despite its actual success in a range of applications for pulse duration down to a few optical cycles, the GNLSE relies on a large number of coefficients, and it is not integrable in general.

In the present contribution, our general motivation is to unveil stable few-cycle pulse propagation effects in the context of supercontinuum generation spanning over several octaves. This investigation, from the theoretical point of view, can hardly be undertaken in the frame of an arbitrarily high-order GNLSE since, beyond the mathematical complexity, the validity of the latter approach remains questionable. Instead, we propose to consider field models that do not use the SVEA, and rely on few equation parameters, allowing the investigation of general classes of few-cycle pulse dynamics.

At this point, let us emphasize on two main different types of solitons. This classification is decisive in the frame of the mathematical soliton theory of integrable equations and the inverse-scattering-transform method. One type corresponds to the envelope soliton, which is the usual description of optical solitons in the picosecond range (Fig. 1(a)). As discussed above, the definition of the envelope assumes the SVEA, which means that the pulse spatial extension L is much larger than the wavelength λ. The other type, called a solitary wave, consists in a single oscillation, without any carrier wave (Fig. 1(b)). It is best portayed by the hydrodynamical soliton that was first observed by John Scott Russel in 1834 [START_REF] Russell | Meeting of the British Association for the Advancement of Science[END_REF], and was modeled by the Korteweg-de Vries (KdV) nonlinear partial differential equation, which was in fact discovered by Boussinesq [START_REF] Boussinesq | [END_REF]6]. The case of few-cycle pulses (FCPs) may appear as an intermediate situation, where we could consider an envelope and a carrier. However, since the pulse duration L has the same order of magnitude as the wavelength λ, the SVEA, which is intrinsically a perturbative approach, is not valid in general, and uselessly complicates the mathematical approach of solitary waves. Therefore a completely different approach to few-optical-cycle pulse propagation has been developed, based on KdV-type models [7,8]; see also the reviews [9][10][11]. The first objective of this approach is to highlight a few classes of few-cycle pulse formation and solitary wave dynamics, which can be relevant in supercontinuum generation and other ultrashort pulse propagation in highly nonlinear waveguides. Note that in the following, the transverse spatial dimensions are excluded. At this early stage, since the parameters of the KdV-type equations are not easily connected with the experimental characterization parameters of materials, we leave a quantitative comparison between model predictions and experiments for future investigations, while focusing on the major physical mechanisms at play.

Models for few-cycle solitons 2.1 Assumptions

With the goal to build simple but effective models, some approximations were required, while taking into account the very large bandwidth required for FCP soliton propagation. The absorption spectrum of most transparent media contains two absorption bands in general: one below, and one above the transparency window. In a first approximation, we can reduce both absorption bands to single transitions as sketched in Fig. 2, so that the general material becomes modeled by the very simple model of a two-component medium, each component being described by density matrix equations for two-level models.

The approximation used to simplify the mathematical model is that the bandwidth of the transparency window is very large, i.e., we assume that the FCP duration τ p is such that ω 1 ≪ (1/τ p ) ≪ ω 2 , where ω 1 and ω 2 are the central frequencies of the infrared (IR) and ultraviolet (UV) transitions, respectively.

Models

The mKdV equation. In a first stage, the two components are treated separately. We first consider the UV transition only, with the assumption that (1/τ p ) ≪ ω 2 . This assumption corresponds to a long-wave approximation, which is performed according to the standard procedure of the reductive perturbation method, or multiscale expansion method [12]. It results in an asymptotic model, which is the modified Korteweg-de Vries (mKdV) equation, as [8] ∂E ∂ζ = 1 6

d 3 k dω 3 ω=0 ∂ 3 E ∂τ 3 - 6π nc χ (3) (ω; ω, ω, -ω) ω=0 ∂E 3 ∂τ . ( 1
)
Here E is the electric field, ζ is the propagation distance and τ is the retarded time in a frame moving at the speed of light c/n, where c is the light speed in vacuum, and n is the medium optical index in the limit of long wavelengths and low amplitudes. The coefficients of Eq. (1) reflect in a simple general form the dispersion relation k(ω) of the medium and its cubic nonlinear susceptibility χ (3) . It should be emphasized that the mKdV equation is completely integrable by means of the inverse scattering transform (IST) method [13,[START_REF] Dodd | Solitons and Nonlinear Wave Equations[END_REF].

The sine-Gordon equation. In a second stage, we consider the IR transition only, with the assumption ω 1 ≪ (1/τ p ). This corresponds to a short-wave approximation, which is another situation well-known within the frame of the reductive perturbation method, although a bit less standard than the previous one. It results in the sine-Gordon (sG) equation, as

[8] ∂ 2 ψ ∂z∂t = c 1 sin ψ. (2)
In Eq. ( 2), c 1 = w ∞ /w r is the normalized initial population difference, and the field ψ is defined so that ∂ψ/∂t = E/E r is the normalized electric field. The propagation distance z and retarded time t are also normalized, see [8].

The mKdV-sG equation. Then the two approximations are brought together to yield a general model, which is known as the mKdV-sG equation. It is expressed in terms of the field ψ as [START_REF] Leblond | [END_REF]:

∂ 2 ψ ∂z∂t + c 1 sin ψ + c 2 ∂ ∂t ∂ψ ∂t 3 + c 3 ∂ 4 ψ ∂t 4 = 0, (3)
and in terms of the normalized electric field u = ∂ψ/∂t = E/E r as

∂u ∂z + c 1 sin t u + c 2 ∂u 3 ∂t + c 3 ∂ 3 u ∂t 3 = 0. ( 4
)
Equation ( 3) or ( 4) is integrable by the IST method in some special cases: when it reduces to the mKdV equation (c 1 = 0), or to the sG equation (c 2 = c 3 = 0) [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF], or if a specific relation between the coefficients is satisfied, namely c 3 = 2c 2 [START_REF] Konno | [END_REF]. In each integrable case it possesses an analytical solution, called a breather, which is a localized oscillatory pulse that propagates without deformation except for periodic oscillations (Fig. 3). The breather is hence a stable short pulse that propagates without being spread out by dispersion, and closely corresponds to the intended shape of a few-cycle optical soliton.

A few generalizations 3.1 General Hamiltonian

We assumed that each of the IR and UV absorption bands could be reduced to a single transition. How does the models described above generalize to a higher number of atomic levels? In the case of the UV transitions, it has been shown that the mKdV model as given by Eq. ( 1) generalizes without modification to a general Hamiltonian, with an arbitrary number of atomic levels [18]. In the case of the IR transitions, the generalization is not so simple.

Indeed, the sG model (Eq. ( 2)) originates from the following coupled equations:

∂E ∂z = -N ε 0 c ΩQ, (5)
∂w ∂t = -EQ, (6)
∂Q ∂t = |µ| 2 Ew, (7)
in which Ω is the angular frequency of the transition, N is the density of atoms, µ is the dipolar momentum matrix element, w is the population difference, and Q is some auxiliary field, while obviously ε 0 is the dielectric permittivity of vacuum, and is Planck's constant. The system ( 5)-( 7) is formally identical to the equations of the self-induced transparency (SIT) equations, with the major difference that here E and Q are real fields, and Q is not the polarization density, but is proportional to its τ -derivative. We recall that the SIT model concerns the envelope of the field close to resonance, while Eqs. ( 5)-( 7) concern the electric field itself, far from the resonance.

Since the population difference is explicitly involved, the generalization to a multilevel model requires one population difference w j variable for each transition. As an example, for a four-level system with two allowed transitions, the system ( 5)- (7) generalizes to [19]

∂E ∂z = -N Ω ε 0 c (Ω 1 Q 1 + Ω 2 Q 2 ), ( 8
)
∂w j ∂t = -EQ j , j = 1, 2, (9)
∂Q j ∂t = |µ j | 2 Ew j , j = 1, 2, (10)
which in turn reduces to

∂u ∂z + c 1 sin τ udτ ′ + qc 1 sin ν τ udτ ′ = 0, (11)
in normalized form, with some real coefficients c 1 , q, and ν depending on the Ω j , |µ j |, and initial values of w j . Equation ( 11) is a generalized double sine-Gordon equation. Under certain conditions, it becomes precisely the so-called double sine-Gordon equation, when the argument of the second sine function is the second-harmonic of argument of the first one (ν = 2). Both original and generalized double sG equations admit breather solutions, examples of which are shown on Fig. 4.

Circular polarization -long wave case

Importantly, both mKdV and sG models can be generalized in order to incorporate wave polarization dynamics. In the case of the mKdV model, i.e., when only the UV transition is considered, assuming as above that the wave central frequency is very far from the resonance, the normalized vector electric field (U, V ) obeys the following system of two coupled mKdV equations [20]:

U Z = U T T T + U 2 + V 2 U T , (12)
V Z = V T T T + U 2 + V 2 V T . (13)
Setting

f = U + iV, (14)
the coupled mKdV equations Eqs. ( 12) and ( 13) reduce to

f Z = f T T T + |f | 2 f T , (15)
which is the complex modified Korteweg-de Vries (cmKdV) equation. However, confusion must be avoided between equation [START_REF] Leblond | [END_REF] and the other cmKdV equation

f Z = f T T T + |f | 2 f T . (16)
Indeed, Eq. ( 16) is completely integrable [21] while Eq. ( 15) is not. Equation ( 15) is usually referred to as cmKdV I, while Eq. ( 16) is called cmKdV II.

Linearly polarized solitons can be described by the solutions of Eq. ( 15) of the form

f = f 0 (Z, T )e iφ , (17)
where f 0 is a solution of the real mKdV equation and φ is a constant. More interesting would be a circularly polarized soliton, of the form

f = u(T -wZ)e i(ωT -kZ) . ( 18
)
However no exact, even numerical, steady state solution of this type do exist. It is however shown that very robust approximate ones do exist. An example is shown on Fig. 5.

The circularly polarized FCP is unstable in the sub-cycle range (for ωτ p < 1.4), as appears on Fig. 6. The initial circularly polarized sub-cycle pulse evolves into a single-humped, linearly polarized pulse, which is nothing but a fundamental soliton of the real mKdV equation.

Circular polarization -short wave case

The case of one IR transition, i.e., when we intend to generalize the sG model to a vectorial normalized field (U, V ), is far from being identical, although analogous features are observed. Since the excited state for x-polarization cannot be identified with the excited state for y-polarization, two distinct population differences w 1 and w 2 are required. The asymptotic equations cannot be reduced to coupled sG-type equations, but yield a complicated system instead, which can be written in complex form as [22]: where f is the complex electric field defined as in Eq. ( 14) above, P = P x + iP y is the polarization density, w = (w 1 + w 2 )/2, s = (w 2 -w 1 )/2 -iS are the population differences (S is the coherence between the two excited states), and the real quantity K and the complex one M are auxiliary fields.

f Z = -P T , (19)
P T = -M -iU K, ( 20
) M T = wf -sf * , ( 21
) K T = Im (f * P ) , ( 22
) s T = f M, ( 23
) w T = -3Re (f * M ) , (24)
As in the case of the cmKdV model, we can retrieve the scalar model form for linear polarizations by setting f as in Eq. ( 17); then the system ( 19)-( 24) reduces to the sG equation. More interesting are the solutions corresponding to circularly polarized FCP solitons. As in the case of the cmKdV, it is seen that no exact solutions of this type exist, whereas robust approximate ones do exist.

An approximate solution valid in the SVEA can be derived, as

f = be i k-b 2 k ω 2 Z-ωT sech b T + k ω Z . ( 25
)
The solution given by ( 25) can be used as initial data for the numerical resolution of the system ( 19)- (24). It is shown that the circular pulse is stable in the SVEA limit, but it is unstable in the sub-cycle range. Figure 8 shows the initial circular polarization, and the final state with two pulses having orthogonal linear polarizations.

4 Supercontinuum generation

Supercontinuum generation in the mKdV and sG models

It is known that an intense laser pulse launched in an optical fiber experiences spectral broadening, which can extend to several octaves, generating white light flashes, whose coherence largely depends on the input pulse parameters [23][24][25][26]. Experiments have been performed in particular in photonic crystal fibers and in tapered fibers, where confinement and nonlinearity can be exacerbated. A GNLSE model is used in general to confront the experimental data. However, the GNLSE model requires numerous parameters, whereas fundamentally, it remains based on the SVEA, whose validity is restricted to less than one octave spectral width. In spite of the good agreement between some theoretical works using GNLSE and experimental results, the principle of a SVEA is highly questionable from the fundamental point of view. We show in the present paper that the non-SVEA models presented above are consistent with supercontinuum generation with a reduced number of model parameters.

The mKdV equation ( 4) with c 1 = 0, in a normalized (dimensionless) form with c 2 = c 3 = 1, is solved numerically using a standard fourth-order Runge-Kutta scheme, and finite differences for the t-derivatives. The initial condition is a Gaussian pulse, with carrier angular frequency ω and duration τ :

u(0, t) = A sin(ωt)e -t 2 /τ 2 . ( 26
)
The ratio between the optical period 2π/ω and the pulse duration τ is chosen to be the same as in a 100 fs long pulse with wavelength λ = 1 µm. Therefore, it is more convenient to do as if t was expressed in fs and ω = 2πc/λ, with c = 0.3 µm/fs, and τ = FWHM/ √ 2 ln 2, where FWHM is the full width at half maximum of the pulse. Numerical resolution shows that a very broad spectrum quickly develops, see Fig. 9.

Generation of long wavelengths

Among the questions which naturally arise, let us first ask ourselves how frequencies lower than the initial pulse frequency ω can be generated? The low-frequency spectrum is more easily seen when the spectrum is plotted against the wavelength λ. Figure 10(a) shows this spectrum for the mKdV model, whereas Fig. 10(b) is the analogous one for the sG model, starting with the same input (26) as above, corresponding to a 100 fs pulse, λ = 1 µm. By comparing Fig. 10(a) to Fig. 10(b) it is seen that analogous spectrum profiles are obtained when either the pulse duration or its amplitude is varied. In the case of the mKdV model, a few specific features are observed, see Fig. 10(a). A supercontinuum is generated in both cases; there is a strong asymmetry between Stokes and anti-Stokes sides. The asymmetry is more pronounced in the case of the sG model (Fig. 10(b)), which in addition shows more broadening towards large wavelengths than the mKdV model (Fig. 10(a)), for which the spectrum first extends towards the short wavelength, the anti-Stokes side.

Recalling that the mKdV model accounts for the UV transitions, while the sG model accounts for the IR transitions, a complete description of a general medium naturally requires both models. The extension of the spectrum towards low frequencies at the beginning of the supercontinuum generation process is frequently attributed to the effect of Raman scattering, where Raman scattering is due to the interaction of light with the IR transitions of the medium. The spectral broadening towards low frequencies we observe is thus nothing but the one due to Raman scattering. Hence, although we cannot pretend that the sG model describes the Raman effect alone, we see that it gives an account of the corresponding spectral broadening.

Self-phase modulation

Let us consider a quasi-monochromatic wave with angular frequency ω and wave vector k, as u = U (z, t)e i(kz-ωt)

+ cc + u 1 (z, t), ( 27
)
where U is the amplitude of the fundamental wave, cc holds for 'complex conjugate', and u 1 is some small correction. Equation ( 27) is reported into Eq. ( 4), in which we neglect dispersion, and disregard third harmonic generation, which yields

U = A 2 exp i c 1 2ω 3 + 3ωc 2 A 2 4 z . (28)
Equation ( 28) accounts for the self-phase modulation (SPM), a well-known phenomenon which leads to a broadening of the spectrum with typical oscillations of the spectral density. We can compare the analytical expression of the self-phase modulation given by Eq. ( 28) to the numerical solution of the mKdV equation, see Fig. 11. It is seen that the analytic envelope approximation coincides with numerics until z ≃ 4, then actual broadening becomes asymmetric, while analytic formula remains symmetric. This shows that the very beginning of the broadening process is actually due to self-phase modulation.

High harmonics generation

Lots of high harmonics are created and involved in the process of supercontinuum generation. Such effect cannot be accounted for within the SVEA unless an amplitude is explicitly included for each harmonic in the model. Up to 15 harmonics can be seen in the example of Fig. 12. Oscillations of their amplitudes can be observed, and are easily retrieved by a simple analytical approach. We can gain further understanding from a simple analytic approach of the third harmonic generation. Let us indeed consider the approximate solution of the mKdV-sG equation given by ( 27) and (28). The correction term u 1 in ( 27) is a third harmonic

u 1 (z, t) = V (z, t)e i(k ′ z-3ωt) + cc + u 2 (z, t), (29)
with k ′ = k(3ω) and u 2 a higher-order correction term. Its amplitude V can be computed with the same approximations as above, as

∂V ∂z = i -c 1 6ω 3 + 3ωc 2 U 3 e i(3k-k ′ )z . (30)
Assuming that U is a constant, and no incident third harmonic (V = 0 at z = 0), this expression leads to well-known oscillations of the third harmonic amplitude due to the phase mismatch ∆k = 3k -k ′ , with V ∝ sin ∆kz 2 . Notice that, if the SPM effect is important, the phase of U in Eq. ( 30) evolves according to Eq. ( 28), and consequently the phase mismatch ∆k should be replaced by

∆k ′ = 3k -k ′ + 3 c 1 2ω 3 + 3ωc 2 . ( 31
)
These simple analytical results are found to be in quantitative agreement with numerics. Let us now determine the spectral width of harmonics. Assuming a Gaussian profile

u = Ae -t 2 τ 2 e -iω0t , (32)
its Fourier transform has the well-known expression

û = 2 √ π Aτ e -(ω -ω 0 ) 2 τ 2 4 . (33)
The nth harmonic evolves essentially as u n does, and is straightforwardly seen that its duration is τ / √ n. Consequently its Fourier transform is given by replacing τ by τ / √ n (and A by A n ) in (33), and hence the spectral width of the nth harmonic is

2 √ n τ . (34)
It increases as √ n. Figure 13 shows the comparison between the actual spectrum obtained by resolution of the mKdV equation (blue), and Gaussian profiles with width increasing as √ n (red). The agreement allows to conclude that the spectral width of the harmonics is actually governed by this law, at least at the considered stage of the process.

We expect also some spectral broadening due to parametric interaction between the sidebands of the harmonics and the ones of the fundamental. We denote by ω the fundamental angular frequency, and assume some sideband ω + δω. Then the third harmonics contains the sideband 3 (ω + δω). It may interact with the fundamental with frequency ±ω, according to 3 (ω + δω) -ω -ω = ω + 3δω. Such a parametric process relates a sideband at ω + δω to another one at ω + 3δω, and is expected to result in spectral broadening. However, it was difficult to get convincing numerical evidence of the occurring of this process, due to the concurrence of another one, which will be discussed in the next section. 5 Few cycle solitons in supercontinuum generation

The sG model

We solve numerically the sG equation, Eq. ( 4) with c 2 = c 3 = 0, and c 1 = 50 (a change in c 1 results only in a change of the unit along the z-axis), with an initial data in a form of a Gaussian pulse as Eq. (26). For moderate amplitudes and high frequencies, a few FCP solitons form and tend to separate (Fig. 14). However, due to the dispersion properties of the sG equation and the characteristics of the input pulse, their velocity is quite small; they don't separate completely, but interact forming some interference pattern, see panel (a) in Fig. 15. Nonlinear interaction between FCP solitons produces very sharp peaks, which induce a huge broadening of the spectrum, see see panel (b) in Fig. 15. Finally, the evolution of the supercontinuum generation process according to the sG model can be summarized as illustrated by Fig. 16. The input spectrum (z = 0) is the blue curve. First, it occurs the generation of harmonics (dark green, z = 4), then the harmonics start to broaden (red, z = 8), and the gaps between them start to fill. When the first soliton becomes apparent (light blue, z = 9.2), the gaps are completely filled. Then the slope of the decay of the spectrum reduces, as seen in the green curve when the second soliton is just formed (z = 11.16), to reach eventually an ultrabroad spectrum (magenta, z = 20).

The mKdV model

We solve in a same way the mKdV equation ( 4) with c 1 = 0, starting from the same initial data with FWHM = 100 and ν = 0.3. At high amplitude level (A ∼ 2), an initial 100 fs pulse splits into a set of FCP solitons, which are mKdV breathers (Fig. 17). The FCPs form in a first stage until z ≃ 12 then slowly go away one from the other. A huge spectral broadening occurs when they form. Comparison between temporal profile and spectrum shows indeed that the spectrum extends at the very point when the solitons separate (Fig. 18). At a lower amplitude (A ∼ 0.1, 0.2), it may happen that only two solitons form. They interact and form a kind of superbreather (Fig. 19). In the meantime, the spectrum oscillates, showing a recurrence to the initial state that recalls the Fermi-Pasta-Ulam recurrent dynamics (Fig. 20), see Ref. [27].

Conclusion

We have presented three generic non-SVEA models: modified Korteweg-de Vries, sine-Gordon, and mKdV-sG models. Initially developed in the frame of a scalar model, of a two-component medium in which each component obeys a two-level model, these models have been generalized to more complicated Hamitonians, and to vector fields. We showed numerically that circularly polarized few-cycle solitons are particularly robust, but when sub-cycle duration is considered, for a spectral bandwidth spanning over several octaves, these solitons become unstable and decay into one or two linearly polarized single-humped solitons.

The non-SVEA models numerically predict supercontinuum generation over several octaves, and we identified several mechanisms in the process. First, the sG term accounts for Raman broadening. Second, a modulation instability mechanism entails spectral broadening at the beginning of the process. Generation of high-harmonics has also been evidenced, with its relation to the formation of the supercontinuum. Finally, we observed that FCP solitons are formed in the process, and interact nonlinearly together: their formation is closely related to that of the ultrabroad spectrum, in correspondence with the known impact of picosecond/subpicosecond soliton for-mation for supercontinuum generation from GNLSE modeling. At present time, a definite advantage of non-SVEA models is that they do not require a high number of terms and coefficients, while they procure the possibility to explore few-cycle to sub-cycle ultrashort pulse dynamics on safer theoretical grounds. However, the tuning of non-SVEA models and parameters to realistic material properties and actual experiments calls for further studies.

Fig. 1 .

 1 Fig. 1. Schematic representation of (a) envelope vs. (b) solitary wave types of solitons.

Fig. 2 .

 2 Fig. 2. General scheme of the material spectral model.

Fig. 3 .

 3 Fig. 3. The typical analytical breather solution.

Fig. 4 .

 4 Fig.[START_REF] Russell | Meeting of the British Association for the Advancement of Science[END_REF]. Breathers of the generalized double sG equation. Left: ν = 2, q = 0.2; right: ν = √ 3, q = 0.4.

Fig. 5 .Fig. 6 .

 56 Fig. 5. The stable circularly polarized FCP: amplitude |u|.

Fig. 7 .

 7 Fig. 7. The lifetime of the circularly polarized FCP, measured in units of the dispersion length, against the number of optical cycles in the pulse (Red line: numerical data. Light blue line: exponential fit).

Fig. 8 .Fig. 9 .

 89 Fig. 8. Trajectory of the electric field vector (u, v) = (Re(f ), Im(f )) in the transverse plane of the wave. Blue: initial, one circularly polarized pulse, red: final, two linearly polarized pulses.

Fig. 10 .

 10 Fig. 10. a: The mKdV model, input with FWHM = 100, ν = 0.3, A = 0.7, for z = 0, 3.1, 6.2, 9.3, 12.4, and 15.5. b: The sG model, input with FWHM = 100, ν = 0.3, A = 2.5, for z = 0, 0.8, 1.6, 2.4, 3.2, and 4.

Fig. 11 .

 11 Fig. 11. Analytical self-phase modulation effect (red dotted lines), vs the mKdV model (blue solid lines); input with FWHM = 100, ν = 0.3, A = 0.7, and for z = 0, 2, 4, 10, and 20.

Fig. 12 .

 12 Fig. 12. Initial pulse with FWHM = 80, ν = 0.375, and A = 5, according to the sG model. Notice the high number of harmonics, and the oscillations of their amplitudes.

Fig. 13 .

 13 Fig. 13. Blue: mKdV equation, Red solid lines: Gaussian with width increasing as √ n, Green dotted lines: Gaussian with constant width. Here FWHM = 100, ν = 0.3, and A = 0.7, for z = 2.

Fig. 14 .

 14 Fig. 14. Evolution of the profile of an initial pulse with FWHM = 100, ν = 0.3, and A = 1, according to the sG model at the end of our computation (z = 20).

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Initial pulse with FWHM = 100, ν = 0.3, and A = 2.5, according to the sG model. a) FCP solitons form, don't separate, but interact; b) generation of an ultrabroad supercontinuum.

Fig. 17 .

 17 Fig. 17. Evolution of an initial pulse with FWHM = 100, ν = 0.3, and A = 0.7, according to the mKdV model (same as in Fig. 18). a) initial pulse, z = 0, b) set of emerging FCP solitons, z = 999, c) zoom on some FCP solitons in b) (the soliton speed is compensated for clarity).

Fig. 18 .

 18 Fig. 18. Initial pulse with FWHM = 100, ν = 0.3, and A = 0.7, according to the mKdV model. left: temporal evolution, right: spectrum evolution.

Fig. 19 .Fig. 20 .

 1920 Fig.19. Initial pulse with FWHM = 100, ν = 0.3, and A = 0.17. Two FCP solitons are formed, and propagate together as a higher order NLS 'super-breather'.