Nelson F Barroso

Rosane Ushirobira

Denis Efimov

Universal formula for robust stabilization of affine nonlinear multistable systems

In this paper, the problem of robust stabilization of affine nonlinear multistable systems with respect to disturbance inputs is studied. The results are obtained using the framework of input-to-state stability (ISS) and integral input-to-state stability (iISS) for systems with multiple invariant sets. The notions of ISS and iISS control Lyapunov functions as well as the small control property are extended within the multistability framework. It is verified that the universal control formula can be applied to yield the ISS (iISS) property to the closedloop system. The efficiency of the proposed control Lyapunov function in the multistable sense is illustrated in two academic examples.

I. INTRODUCTION

Due to its importance for several scientific disciplines ranging from mechanics and electronics [19], [START_REF] Angeli | Characterization of input-to-state stability for systems with multiple invariant sets[END_REF] to biology [26],

[20], [START_REF] Pchelkina | Control of oscillatory behavior of multispecies populations[END_REF] and neurosciences [25], the analysis of stability and robustness properties of multistable systems have become increasingly attractive from the perspective of systems and control theory. Multistable systems include bistable systems (with at least two stable equilibria) [START_REF] Yakubovich | Stability of Stationary Sets in Control Systems with Discontinuous Nonlinearities[END_REF], [START_REF] Chaves | Bistable biological systems: A characterization through local compact input-to-state stability[END_REF], almost globally stable systems (with only one attracting invariant set) [START_REF] Angeli | An almost global notion of input-to-state stability[END_REF], and nonlinear systems with generic invariant sets [START_REF] Angeli | Detection of multistability, bifurcations and hysteresis in a large class of biological positivefeedback systems[END_REF], [START_REF] Dudkowski | Hidden attractors in dynamical systems[END_REF], [14], [START_REF] Forni | Differential analysis of nonlinear systems: Revisiting the pendulum example[END_REF], [18], [START_REF] Rumyantsev | Stability and stabilization of motion with respect to part of variables[END_REF], [START_REF] Stan | Analysis of interconnected oscillators by dissipativity theory[END_REF], [33]. In [12], it was proposed a global asymptotic stability notion as well as the necessary and sufficient Lyapunov characterization for multistable systems, having as the object of investigation all compact invariant solutions of the system (including locally stable and unstable ones).

By virtue of nontrivial relationships between different regions that compose its state space, and the complex intertwined boundaries between then, multistable systems are extremely sensitive to initial conditions and perturbations. The characterization introduced in [12] has been also proved to be useful in robustness analysis with respect to external disturbances. In fact, it was made clear in [3] that the most natural way to solve this problem is to relax the Lyapunov stability requirement on relatively mild additional assumptions on the decomposition of invariant sets. This intuitive path has led to a new line of research which starts from the characterizations of input-to-state stability (ISS) for this class of systems in terms of usual Lyapunov dissipation N. F. Barroso, R. Ushirobira and D. Efimov are with Inria, University of Lille, CNRS, UMR 9189 -CRIStAL, F-59000 Lille, France. D. Efimov is also with ITMO University, 197101 Saint Petersburg, Russia. This work was supported by the ANR project WaQMoS (ANR 15 CE 04 0002). It was also supported by the Ministry of Science and Higher Education of Russian Federation, passport of goszadanie no. 2019-0898. N. F. Barroso is partly supported by the regional council of Hauts-de-France.

inequalities, generalizing the classical ISS theory [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF], [START_REF] Sontag | New characterizations of the input-to-state stability[END_REF],

[8], [3], [START_REF] Angeli | A characterization of integral input-to-state stability[END_REF]. In its turn, integral input-to-state stability (iISS) characterization, which is weaker than the classical one given in [START_REF] Sontag | Coments on integral variants of iss[END_REF], [21], [5], was extended in [16] for systems with multiple invariant sets.

Once ISS and iISS characterizations in the multistable sense are already available, the research interest on the problem of designing robust stabilizing control laws in this framework naturally increases. In this setting, such a problem consists in finding state feedback control laws that make the closed loop system ISS or iISS stable with respect to a family of finite disjoint compact invariant sets against external disturbances effects. In the classical approach, most of this activity is centered around the control Lyapunov function (CLF) theory [START_REF] Artstein | Stabilization with relaxed controls[END_REF], [START_REF] Sontag | A "universal" construction of artstein's theorem on nonlinear stabilization[END_REF], [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF], [11], [10]. In [START_REF] Artstein | Stabilization with relaxed controls[END_REF], it was shown that the existence of a CLF leads to an explicit formulation for stabilizing control laws. Similar results were proven in [START_REF] Liberzon | Universal contruction of feedback laws achieving iss and integral-iss disturbance attenuation[END_REF] for the ISS (iISS) case resulting in appropriated universal formulation for the assignment of a ISS (iISS) CLF rendering ISS (iISS) properties to the closedloop system. In the present work, we are interested in the robust stabilization of multistable affine nonlinear systems with respect to disturbance inputs. Our approach is based on the theory developed in [12], [3], [16] and it aims to find conditions of CLF existence in the context of systems with multiple invariant sets (compact and maybe disconnected) and show how such a control can be explicitly designed for robust stabilization in ISS (iISS) sense.

The outline of this work is as follows. The main definitions and problem statement are given in Section II, while the obtained results are presented in Section III. In section IV are presented two examples of application of the proposed CLF approach. Final remarks and discussion are summarized in Section V.

II. DEFINITIONS AND PROBLEM STATEMENT

Let M be a n-dimensional smooth manifold without boundary, equipped with a metric δ : M × M → R + , R + = {s ∈ R : s ≥ 0}. Consider a nonlinear model of dynamical systems evolving on this manifold:

ẋ(t) = f (x(t), u(t)), ∀ t ∈ R + , (1)
y(t) = h(x(t)), (2)
where x(t) ∈ M is the state vector, u(t) ∈ U ⊆ R m is the input vector, u is an element of U, the set of admissible controls R + → U (locally essentially bounded and measurable signals), and y(t) ∈ Y ⊆ R p is the output vector. Let f : M×U → T x M be a locally Lipschitz continuous function on M (here T x M denotes the tangent space of M at x), and assume that h : M → Y is continuously differentiable, h(0) = 0 and f (0, 0) = 0. Denote by x(t, x 0 ; u) the uniquely defined solution of (1) at time t ≥ 0 such that x(0) = x 0 under the input u ∈ U. For the unperturbed system, i.e. the system (1) with u ≡ 0, we have:

ẋ(t) = f (x(t), 0), t ≥ 0, (3)
and we say that S ⊂ M is invariant if for all x 0 ∈ S, x(t, x 0 ; 0) ∈ S for all t ∈ R.

For a set S ⊂ M and x ∈ M define the corresponding distance as

|x| S = inf a∈S δ(x, a).
For a measurable function g : R

+ → R m , define its L ∞ - norm as ||g|| ∞ = ess sup t≥0 |g(t)|.

A. Decomposition of a compact invariant set

Let Λ ⊂ M be a compact invariant set for the unperturbed system (3). To characterize the evolution of this system along M, it is useful to decompose Λ and explicitly determine the existence of solutions traveling between different components of its decomposition.

Definition 1. [23] A decomposition of Λ is a finite, disjoint family of compact invariant sets Λ 1 , • • • , Λ k such that Λ = k i=1 Λ i .
For an invariant set Λ, its attracting and repulsing subsets can be defined, respectively, as follows:

A(Λ) = {x 0 ∈ M : |x(t, x 0 ; 0)| Λ → 0 as t → +∞}, R(Λ) = {x 0 ∈ M : |x(t, x 0 ; 0)| Λ → 0 as t → -∞}.
Based on these definitions, we can define a relation between two invariant sets W ⊂ M and D ⊂ M by W ≺ D if A(W) ∩ R(D) = ∅. This relation implies that there is a solution connecting the set D with the set W. A collection of r disjoint sets that can be reached from one to another in a loop by a suitable concatenation of systems solutions is called r-cycle.

Definition 2. [23] Let Λ 1 , • • • , Λ k be a decomposition of Λ. 1) An r-cycle (r ≥ 2) is an ordered r-tuple of distinct indices i 1 , • • • , i r such that Λ i1 ≺ • • • ≺ Λ ir ≺ Λ i1 . 2) A 1-cycle is an index i such that (R(Λ i ) ∩ A(Λ i )) \ Λ i = ∅.
3) A filtration ordering is a numbering of the Λ i so that

Λ i ≺ Λ j ⇒ i ≤ j.
However, what qualifies a decomposition for its treatment by means of Lyapunov-like analytical tools is the absence of these cycles. Therefore, we will next consider the following assumption: Assumption 1. [3] A compact invariant set W containing all αand ω-limit sets of the unperturbed system (3), admits a finite decomposition without cycles: W = k i=1 W i for some non-empty disjoint compact sets W i which forms a filtration ordering of W, as detailed in Definition 2.

B. Robust stability notions for a decomposable compact invariant set W

In this subsection, we list several ISS and iISS stability properties for (1) with respect to W satisfying Assumption 1. Most of these properties are direct extensions of the classical ISS and iISS notions introduced in [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF], [START_REF] Sontag | New characterizations of the input-to-state stability[END_REF], [21], [5]. The definitions of function classes K and K ∞ can be found in [8]. A function V : M → R + is called positive definite if it vanishes only at the origin, and proper unbounded if V (x) → +∞ for |x| W → +∞. The Lie derivative of a continuously differentiable function V along a vector field f : M → R n is denoted by:

DV (x)f (x, u) = ∂V (x) ∂x f (x, u). Definition 3. [3],
[16] The system (1) has the practical asymptotic gain (pAG) property if there exist η ∈ K ∞ and q ≥ 0 such that for all x 0 ∈ M and all u ∈ U, the solutions are defined for all t ≥ 0 and the following holds:

lim sup t→+∞ |x(t, x 0 ; u)| W ≤ η(||u|| ∞ ) + q. (4)
If q = 0, then we say that the asymptotic gain (AG) property holds. Moreover, if (4) is satisfied for q = 0 for the system (3) only, then we will say that (1) has the zero-global attraction (0-GATT) property with respect to a compact invariant set W.

Definition 4.

[3] The system (1) has the limit property (LIM) with respect to W if there exists µ ∈ K ∞ such that for all x 0 ∈ M and all u ∈ U the solutions are defined for all t ≥ 0 and the following holds:

inf t≥0 |x(t, x 0 ; u)| W ≤ µ(||u|| ∞ ).
Definition 5.

[3] The system (1) has the practical global stability (pGS) property with respect to W if there exists β ∈ K ∞ and c ≥ 0 such that for all x 0 ∈ M and all u ∈ U, the following holds for all t ≥ 0:

|x(t, x 0 ; u)| W ≤ β (max{|x 0 | W + c, ||u|| ∞ }) . Definition 6. [3], [16] A C 1 function V : M → R is a practical ISS Lyapunov function for (1) if there exist K ∞ functions α 1 , α 2
, α 3 and γ, and q ≥ 0, c ≥ 0 such that

α 1 (|x| W ) ≤ V (x) ≤ α 2 (|x| W ) + c (5)
and the following dissipation inequality holds for all (x, u) ∈ M × U:

DV (x)f (x, u) ≤ -α 3 (|x| W ) + γ(|u|) + q. ( 6
)
If ( 6) holds for q = 0, then V is said to be an ISS Lyapunov function. If (6) holds for q = 0 and a positive definite function α 3 : R + → R + , then V is said to be an iISS Lyapunov function.

Definition 7.

[16] The system (1) has the uniform boundedenergy bounded-state (UBEBS) property if for some α, γ, σ ∈ K ∞ and some positive constant c, the following estimate holds for all t ≥ 0, all x 0 ∈ M and all u ∈ U:

α(|x(t, x 0 ; u)| W ) ≤ γ(|x 0 | W ) + t 0 σ(|u(τ )|)dτ + c.
Definition 8.

[16] The system (1), ( 2) has the smooth dissipativity property if there exist a

C 1 function V : M → R + , α 1 , α 2 , σ ∈ K ∞ ,
|x| W = 0 ⇒ h(x) = 0, ∀x ∈ M
such that (5) is satisfied for all x ∈ M and the following dissipation inequality holds for all (x, u) ∈ M × U:

DV (x)f (x, u) ≤ -α 4 (|h(x)|) + σ(|u|). (7)
Definition 9.

[16] The system (1), ( 2) has the weak zerodetectability property if the following relation holds:

h(x(t, x 0 ; 0)) ≡ 0, ∀t ≥ 0 ⇒ |x(t, x 0 ; 0)| W → 0 as t → ∞.
The principal results connecting these properties are as following: Theorem 1.

[3] Consider a nonlinear system [START_REF] Ahmed | Robust global synchronization of Brockett oscillators[END_REF] and let W be as in Assumption 1. Then the following are equivalent:

1) The system enjoys the pAG or AG property.

2) The system admits an ISS Lyapunov function.

3) The system admits an ISS Lyapunov function constant on invariant sets. 4) The system admits a practical ISS Lyapunov function. 5) The system enjoys the LIM property and the pGS. The system as in [START_REF] Ahmed | Robust global synchronization of Brockett oscillators[END_REF] that satisfies these properties will be called ISS in the multistable sense with respect to the set W and input u.

Theorem 2.

[16] Consider a nonlinear system (1) and let W be as in Assumption 1. Then the following facts are equivalent:

1) 0-GATT and UBEBS properties.

2) Existence of an iISS Lyapunov function V such that DV (x) = 0 for all x ∈ W.

3) Existence of an iISS Lyapunov function V . 4) Existence of an output function that makes the system

smoothly dissipative and weakly zero-detectable. The system as in [START_REF] Ahmed | Robust global synchronization of Brockett oscillators[END_REF] that satisfies these properties will be called iISS in the multistable sense with respect to the set W and the input u.

C. Problem statement

In this paper, we deal with a subclass of nonlinear dynamical systems (1), (2) affine in the input of the following form:

ẋ(t) = f (x(t), v(t)) + G(x(t))u(t), t ∈ R + , (8) y(t) = h(x(t)), (9)
where x(t) ∈ M is the state vector, v(t) ∈ R k is a disturbance, u(t) ∈ U ⊆ R m is the input vector, u is an element of U, the set of admissible controls R + → U (locally essentially bounded and measurable signals), and y(t) ∈ Y ⊆ R p is the output vector. For this system f : M × R k → R n and the columns of the matrix function G : M → R n×m are assumed to be locally Lipschitz continuous on M, h : M → R p is continuously differentiable, h(0) = 0, and f (0, 0) = 0. Systems in the form (8), ( 9) are said to be ISS (iISS) stabilizable if there exists a control law u = K(x), with K : M → U ⊆ R m , so that the closed-loop system has the ISS (iISS) property in v ∈ R k with respect to W. Therefore, the problem studied in this work can be formally determined as follows:

Problem. Consider the systems described by the affine nonlinear model (8), [START_REF] Dudkowski | Hidden attractors in dynamical systems[END_REF]. Under Assumption 1, find conditions for these systems to be ISS (iISS) stabilizable.

III. MAIN RESULTS

As for ISS (iISS) stabilization with respect to a compact set [21], [START_REF] Sontag | A "universal" construction of artstein's theorem on nonlinear stabilization[END_REF], we will look for existence conditions of a stabilizing feedback using the CLF framework. To present our main result, we introduce a suitable notion of ISS (iISS) CLF from the point of view of multistability framework. Definition 10. A practical ISS CLF for the system (8), [START_REF] Dudkowski | Hidden attractors in dynamical systems[END_REF] and control u(t) ∈ U ⊆ R m is a differentiable function V : M → R + satisfying:

1) There exist class K ∞ functions α 1 and α 2 such that (5) holds for all x ∈ M. 2) There exist class K ∞ functions χ and α 3 and a constant q ≥ 0 such that for all x ∈ M and all

v ∈ R k inf u∈U {a(x, v) + b(x)u} ≤ -α 3 (|x| W ) + χ(|v|) + q, ( 10
)
where a(x, v) = DV (x)f (x, v), b(x) = DV (x)G(x). If (10) holds for q = 0, then V is an ISS CLF. In addition, such V is an iISS CLF if the both items are satisfied for q = 0 and a continuous positive definite function α 3 .

The inequality (10) can be rewritten as follows:

inf u∈U {a(x, v) -χ(|v|) -q + b(x)u} ≤ -α 3 (|x| W ),
where χ ∈ K ∞ is such that χ(s) ≥ χ(s) for all s ∈ R + and there exists a function Finally, define

ψ(x) = sup v∈R k {a(x, v) -χ(|v|) -q},
Ψ(x) = ψ(x) + 1 2 α 3 (|x| W ),
then for all x ∈ M and v ∈ R k :

Ψ(x) ≥ 1 2 α 3 (|x| W ) + a(x, v) -χ(|v|) -q, inf u∈U {Ψ(x) + b(x)u} ≤ - 1 2 α 3 (|x| W ).
Definition 11. A differentiable function V : M → R possesses a variant of the classical small control property (SCP) with respect to x ∈ M if it takes constant values on any W i with i = 1, . . . , k, from the decomposition of W, and for each ε > 0 there is a δ > 0 such that whenever 0 < |x| W < δ there exists some control u ∈ U ⊆ R m with |u| < ε such that

Ψ(x) + b(x)u < 0. (11)
Note that for the both above cases (ISS/iISS), the condition (10) can be replaced with another one:

b(x) ≡ 0 ⇒ Ψ(x) ≤ - 1 2 α 3 (|x| W ), (12)
for all x ∈ M, and α 3 is a class K ∞ function or a positive definite function for the cases of ISS or iISS, respectively. Therefore, the condition (12) formulates the main restriction to check for ISS (iISS) CLF in the system (8), [START_REF] Dudkowski | Hidden attractors in dynamical systems[END_REF]. Following the CLF framework from [START_REF] Artstein | Stabilization with relaxed controls[END_REF], [START_REF] Sontag | A "universal" construction of artstein's theorem on nonlinear stabilization[END_REF], we will use the "universal" control formula, which in our case takes the form:

u = K(x) = κ(Ψ(x), |b(x)| 2 )b(x) T , (13)
where

κ(s, r) = -s+ √ s 2 +r 2 r if x / ∈ W 0 if x ∈ W .
Then, according with the above developments we propose the following theorem and lemma. We omit the proofs due to space limitations.

Theorem 3. Let an affine system (8) admit Assumption 1.

If the function V : M → R + is an ISS (iISS) CLF, then the feedback law [START_REF] Angeli | Characterization of input-to-state stability for systems with multiple invariant sets[END_REF] is continuous in M \ W and it provides the ISS (iISS) property in v with respect to W for (8). If such a CLF satisfies the SCP given in Definition 11, then the control ( 13) is continuous on M.

Note that V (x) = 0 for x ∈ A (it reaches its minimum), when A ⊂ W, then the control (13) will provide almost global or local attraction of the set A in the case v = 0.

Our further directions of research will apply the smooth dissipativity property for iISS stabilization of (8), ( 9), but here we will formulate the following consequence of our main result: Lemma 1. Consider an affine system (8) that satisfies Assumption 1. Then there exists a continuous control that provides the ISS (iISS) property in v with respect to W for (8) if and only if the system admits an ISS (iISS) CLF with SCP.

IV. EXAMPLES

In this section we will illustrate peculiarities and features of the proposed CLF in the multistable framework on two academic examples.

A. Duffing oscillator

Consider the following variant of Duffing oscillator:

ẋ1 (t) = x 2 (t) + u(t) + v 1 (t), ẋ2 (t) = -x 2 (t) + x 1 (t)[x 2 1 (t) -1] + v 2 (t), where x(t) = (x 1 (t) x 2 (t)) ⊤ ∈ M = R 2 , u(t) ∈ R and v(t) = (v 1 (t) v 2 (t)) ⊤ ∈ R 2 have
the same sense as before. It is straightforward to check that for u = 0 and v = 0 the system has three equilibria, thus, W = {(-1, 0), (0, 0), (1, 0)}. Linearization shows that the origin is a locally asymptotically stable, and the two equilibria, (-1, 0) and (1, 0) are unstable foci (see Fig. 1).

Select A = {(-1, 0), (1, 0)} as the set that will be almost globally attractive in the closed-loop system. To this end, as a CLF candidate, choose

V (x) = 1 4 (x 2 1 -1) 2 + 1 2 x 2 2 ,
whose derivative for the oscillator takes the usual form V = a(x, v) + b(x)u, where

a(x, v) = 2x 1 (x 2 1 -1)x 2 -x 2 2 + x 1 (x 2 1 -1)v 1 + x 2 v 2 , b(x) = x 1 (x 2 1 -1). Select χ(s) = s 2 , then a(x, v) -χ(|v|) =     x 1 (x 2 1 -1) x 2 v 1 v 2     ⊤ Q     x 1 (x 2 1 -1) x 2 v 1 v 2     - 1 4 x 2 2 + 2.3x 2 1 (x 2 1 -1) 2 ,
where

Q =     -2.3 1 0.5 0 1 -0.75 0 0.5 0.5 0 -1 0 0 0.5 0 -1    
is a negative definite matrix. Hence, we can take

ψ(x) = - 1 4 x 2 2 + 2.3x 2 1 (x 2 1 -1) 2 ≥ a(x, v) -χ(|v|). Next, let |x| W = x 2 2 + x 2 1 (x 2 1 -1) 2
and α 3 (s) = 1 4 s 2 , then we obtain

Ψ(x) = ψ(x) + 1 2 α 3 (|x| W ) ,
and it is easy to check that b

(x) = 0 ⇒ Ψ(x) ≤ - 1 2 α 3 (|x| W ) .
Consequently, V is an ISS CLF for the original system and for any ε > 0 there is a δ > 0 such that whenever 0 < |x| W < δ there exists some control |u| < ε such that

Ψ(x) + b(x)u < 0,
then SCP holds, and the control ( 13) should be continuous in R 2 . The behavior of the controlled system is shown in Fig. 1. Note that with the applied control the two unstable foci become almost attractive in the closed-loop system. Observe that by starting the system from the same initial conditions the control action forces the trajectories to different attractors. To generate the plots the noise signals were chosen v 1 (t) = v 2 (t) = 0.2 sin(10t).

B. Brockett oscillator

This system has the following model:

ẋ1 (t) = x 2 (t) + v 1 (t), ẋ2 (t) = -x 1 (t) -βx 2 (t) |x(t)| 2 -1 +αu(t) + v 2 (t),
where again

x(t) = (x 1 (t) x 2 (t)) ⊤ ∈ M = R 2 , u(t) ∈ R and v(t) = (v 1 (t) v 2 (t)) ⊤ ∈ R 2 ;
α > 0 and β > 0 are constant parameters. For u = 0 and v = 0 the system has the equilibrium at the origin and an attracting from almost all initial conditions limit cycle on the unit sphere S = {x ∈ R 2 : |x| = 1}, thus, W = {(0, 0), S} [START_REF] Ahmed | Robust global synchronization of Brockett oscillators[END_REF].

Assume that A = {(0, 0)}, which in our example will be locally attracting, and choose in this case as a CLF:

V (x) = 1 2 (x 1 + x 2 ) 2 + 1 2 x 2 1 ,
whose derivative for the Brockett oscillator admits

a(x, v) = -x 2 1 + (x 1 + x 2 )x 2 1 + β -β|x| 2 +(2x 1 + x 2 )v 1 + (x 1 + x 2 )v 2 , b(x) = α(x 1 + x 2 ).
Select χ(s) = s 2 , and similarly we get

a(x, v) -χ(|v|) =     x 1 x 1 + x 2 v 1 v 2     ⊤ Q     x 1 x 1 + x 2 v 1 v 2     - 1 2 x 2 1 + (x 1 + x 2 ) 2 + (x 1 + x 2 )x 2 [1 + β -β|x| 2 ],
where

Q =     -0.5 0 0.5 0 0 -1 0.5 0.5 0.5 0.5 -1 0 0 0.5 0 -1    
is a negative definite matrix. Therefore, we can accept

ψ(x) = - 1 2 x 2 1 + (x 1 + x 2 ) 2 + (x 1 + x 2 )x 2 1 + β -β|x| 2 ,
and, finally take is satisfied and in fact, V is an ISS CLF for the original system. However, it is impossible to demonstrate the SCP property in this case due to the term (x 1 + x 2 )x 2 1 + β -β|x| 2 in Ψ(x) which is not approaching to zero while x becomes close to S. Therefore, in this example the control (13) will be discontinuous on S.

The effect of the noise for both the open-loop and closedloop systems with α = β = 1 is shown in Fig. 2. Again, observe that by initiating the system from the same initial conditions the control action forces the trajectories to different attractors. In simulation the noise signals also were chosen as v 1 (t) = v 2 (t) = 0.2 sin(10t). In this work we dealt with a subclass of nonlinear multistable dynamical systems affine in the input. The main problem addressed here was to establish conditions for the existence of a feedback control that renders the ISS (iISS) stability property with respect to disturbances to the closedloop system. In this direction, we have properly extended the notions of ISS (iISS) control Lyapunov functions and small control property within the multistability framework. In the multistable sense, the existence of a ISS (iISS) CLF satisfying the SCP implies the existence of a feedback control law that can be explicitly designed by the universal formula strategy. To exemplify the ISS (iISS) CLF application such a feedback control was designed to robustly stabilize two different nonlinear multistable systems.

  a continuous positive definite function α 4 , and a continuous output map h : M → R p with

  e.g., take χ(s) = max{χ(s), sup |x|,|v|≤s |a(x, v)|} [21], then with such a choice ψ(x) = sup |v|≤|x| {a(x, v)χ(|v|) -q}. By construction ψ(x) ≤ 0 for all x ∈ W, and since the function a is at least Lipschitz continuous in x, then the function ψ admits the same property. In addition, for all x ∈ W inf u∈U {ψ(x) + b(x)u} ≤ -α 3 (|x| W ).

Figure 1 :

 1 Figure 1: A trajectory convergence for the system with disturbance: without control (dashed line), with control (solid line)

|x| W = |x| 1 +

 1 |x| 2 |1 -|x| 2 | with α 3 (s) = 1 4 s 2 .Then, as beforeΨ(x) = ψ(x) + 1 2 α 3 (|x| W ) ,and the main restriction for V to be a CLF, i.e.b(x) = 0 ⇒ Ψ(x) ≤ -1 2 α 3 (|x| W )

Figure 2 :

 2 Figure 2: A trajectory convergence for the system with disturbance: without control (dashed line), with control (solid line)


A Anantharajah

F Kwabia Tchana

L Manceron

M Ndao

J.-M Flaud

X Landsheere

Integrated band intensities and absorption cross sections of phosgene (Cl 2 CO) in the mid-infrared at 199, 250 and 300 K

Keywords: Phosgene, Cl 2 CO, Fourier transform infrared (FTIR) spectroscopy, Integrated band intensities, absorption cross sections, Temperature dependence, Earth's atmosphere

) in the mid-infrared at 199, 250 and 300 K

Integrated band intensities and absorption cross sections of phosgene (

Introduction

Phosgene (Cl 2 CO) is used as an industrial chemical in the manufacture of plastics, insecticides, pharmaceuticals, and herbicides. Its high toxicity and adverse effects required reducing its use over the years by working with safer substitutes [START_REF] Fu | Global phosgene observations from the Atmospheric Chemistry Experiment (ACE) mission[END_REF]. In the Earth's atmosphere, phosgene is formed in the troposphere by OH-initiated oxidation of chlorinated hydrocarbons, namely : chloroform (CHCl 3 ), methylchloroform (CH [START_REF] Helas | On sources and sinks of phosgene in the troposphere[END_REF] For infrared remote sensing, the 11.8 µm (830-860 cm -1 ) spectral region is mainly used, but its strong absorption, corresponding to the ν 5 band, is masked by the ν 4 band of the Freon-11 (CFC-11 or CCl 3 F), a greenhouse gas that absorbs strongly in the same spectral region. This interference leads to an overestimation of the concentration of CFC-11 in the atmosphere, if the absorption of Cl 2 CO is not taken into account [START_REF] Toon | Atmospheric COCl 2 measured by solar occultation spectrometry[END_REF]. So, modelling the Cl 2 CO contribution in atmospheric spectra is also essential to correct errors in CFC-11 retrievals. Spectroscopically, phosgene is an asymmetric top molecule which, as listed in Table 1, has six infrared fundamental vibrational modes: two in the mid-infrared ν 5 (851.013 cm -1 ) and ν 1 (1828.202 cm -1 ) and four low energy modes at 301.546 cm -1 (ν 3 ), 443.172 cm -1 (ν 6 ), 572.526 cm -1 (ν 2 ) and 582.089 cm -1 (ν 4 ). Its rotational constants have been determined for the ground and first excited states by various microwaves studies [15-21]. Among the fundamental vibrational modes, ν 1 and ν 5 bands are of major interest since they are used for the remote sensing of Cl 2 CO, thanks to their strong absorption in the atmospheric "windows" located at 5.5 and 11.8 µm, respectively. However, in addition to the infrared active fundamentals, the spectrum of Cl 2 CO contains various overtones, combinations and hot bands that originate from the states below 600 cm -1 (ν 2 , ν 3 , ν 4 and ν 6 ). Due to these low energy states, the Cl 2 CO cold bands are expected to be overlapped by many hot bands which contribute significantly to the absorption. Indeed, only ~58% of the population is in the ground state at 300 K (see Table 1). As a consequence the infrared spectrum of Cl 2 CO is very crowded at room temperature.

At 199 K, the ground state population is 82% and thus hot bands are strongly reduced, compared to the situation at room temperature. This is an important point when estimating the contribution of hot bands to the band intensity. This is why in order to reduce the spectral density and then facilitate the spectroscopic analysis, high-resolution studies of the fundamentals of phosgene, considering both the 35 Cl 2 CO and 35 Cl 37 ClCO isotopologues, have been carried out at low temperatures [START_REF] Tchana | High-resolution analysis of the ν 1 and ν 5 bands of phosgene 35 Cl 2 CO and 35 Cl 37 ClCO[END_REF][START_REF] Flaud | First farinfrared high-resolution analysis of the ν 2 and ν 4 bands of phosgene 35 Cl 2 CO and 35 Cl 37 ClCO[END_REF][START_REF] Ndao | First far-infrared high resolution analysis of the ν 3 band of phosgene 35 Cl 2 CO and 35 Cl 37 ClCO[END_REF][START_REF] Flaud | First far-infrared high resolution analysis of the ν 6 band of phosgene 35 Cl 2 CO and 35 Cl 37 ClCO[END_REF][START_REF] Flaud | First infrared high resolution analysis of the ν 5 band of phosgene 37 Cl 2 CO[END_REF]. From these studies, very accurate band centers, rotational, centrifugal distortion and interaction constants have been determined for the six vibrational states of the two isotopologues of Cl 2 CO.

Besides these spectroscopic studies dealing with line positions, the quantitative intensity aspect needs to be investigated for this molecule. However measuring individual Cl 2 CO transitions line-by-line is extremely challenging, even at low sample pressures and with Doppler-limited resolution because, as already said the phosgene spectrum is very dense: Three isotopologues, strong hot bands, etc. A more realistic approach for quantitative measurements is then to determine integrated band intensities and/or cross-sections which are of interest for atmospheric remote sounding applications. Up to now and to the best of our

Experimental details

Absorption spectra of phosgene have been recorded in the range 525 to 2400 cm -1 using the high-resolution Bruker IFS125HR Fourier transform spectrometer (FTS) located at the LISA facility in Créteil (France). The instrument was equipped with a silicon carbide Globar source, a KBr/Ge beamsplitter and a liquid nitrogen cooled HgCdTe (MCT) detector. The diameter of the entrance aperture of the spectrometer was set to 3.15 mm to maximize the intensity of infrared radiation falling onto the MCT detector without saturation or loss of spectral resolution. Interferograms were recorded with a 40 kHz scanner frequency and a maximum optical path difference (MOPD) of 9 cm. According to the Bruker definition (resolution = 0.9/MOPD), this corresponds to a resolution of 0.1 cm -1 . The spectra were obtained by Fourier transformation of the averaged interferograms (100 scans) using a Mertz phase correction with a 1 cm -1 resolution, a zero-filling factor of 2 and no apodization (boxcar option).

For the measurements, a small cryogenic cell (5.10 ± 0.01 cm path length) was used, a schematic viewcut of which is given in Figure 1. The cell is made of stainless steel, to minimize adsorption and corrosion issues, and fitted with 9.5 mm diameter, 0.4 mm thick, wedged diamond windows (E6, Netherlands). The cell is surrounded by an outer stainless steel jacket, in which cold nitrogen (N 2 ) gas is flowed in and out (Figure 1). The temperature control is obtained mainly by varying the cold gas flow rate, but a minor correction can be obtained using an additional heating element (Thermocoax, France) wrapped around the gas cell itself. The temperature measurement is made with four separate Pt-100 sensors: one is inserted in the gas to be measured, a few millimeters above the infrared beam, two others are placed on different sides of the cell body, and a fourth one on the gas tube inlet used for filling the cell. The whole assembly is placed directly inside the sample compartment of the spectrometer, and the vacuum level inside the spectrometer is thus critical to avoid spurious ice absorptions. The FTS was thus evacuated below 3×10 -4 hPa by an all-magnetic bearing turbomolecular pump in order to minimize absorption by atmospheric gases while avoiding excessive mechanical vibrations.

Triphosgene was purchased from Sigma Aldrich and heated to generate phosgene.

Phosgene production begins at a temperature of 80°C, below which the gas is not produced.

The reaction proceeds cleanly up to 110°C at a steady rate. Before use, the sample was purified at -130°C by freeze-and-thaw cycles to remove any non-condensable impurities.

Samples were prepared afterwards by applying the following procedures: first, a mixture of 1.7% Cl 2 CO in nitrogen (Air Liquide, France, 99.9998% purity) was prepared; next, the sample was stabilized for two hours to obtain a dilute and homogeneous sample; and from this, an amount, corresponding to the desired Cl 2 CO amount, was extracted and diluted with about 500 hPa of nitrogen as a broadening gas. Six samples were introduced in the cell at 199 ± 1 K, nine samples at 250 ± 1 K, and six samples at 300 ± 1 K. The sample pressure in the cell was measured using a precision capacitive pressure gauge (MKS Baratron, USA, 10 and 1100 hPa ranges) with a reading accuracy of 0.12% according to the manufacturer's specification. Considering the uncertainty arising from small variations of the pressure during the recording (~ 0.8%), we estimated the measurement uncertainty on the pressure to be equal to 1%. A summary of the experimental conditions is provided in Table 2. Each spectrum was ratioed against a background spectrum of the empty cell, which was recorded at a resolution of 0.1 cm -1 with a large number of scans, in order to ensure the best possible signal-to-noise in the ratioed spectra.

Figure 2 gives an overview of the N 2 -broadened Cl 2 CO spectrum recorded at 300 K (spectrum S1 in Table 2). As displayed in Figure 2, the observed bands are ν 2 /ν 4 at 577.31 cm -1 , ν 5 at 851.01 cm -1 , ν 2 +ν 6 at 1014.78 cm -1 , ν 2 +ν 5 at 1414.92 cm -1 , 2ν 5 at 1681.38 cm -1 , and ν 1 at 1828.20 cm -1 . Note a small impurity of CO 2 observed in the spectrum with the spectral feature at 2350 cm -1 corresponding to the ν 3 band of CO 2 . As the total sample pressure P t includes CO 2 impurity, it is thus required to assess the CO 2 partial pressure contribution in order to determine the true Cl 2 CO partial pressure. To calculate the CO 2 partial pressure, the ν 3 band intensity 9.0710 -17 cm -1 /molecule.cm -2 reported in Ref.

Uncertainty assessment

Examination of the standard deviations of the fits shows that they are generally less than A realistic assessment of systematic photometric errors is not an easy task and probably strongly depends on the setup used. We base our present estimate on the arguments developed in Ref.

[25] using instruments comparable to ours. In particular, HgCdTe detectors are used for their high sensitivity, but are also known for their deviation to a linear behavior. When using an HgCdTe for mid-infrared measurements, deviations up to 2.3% on retrieved band intensities were found when compared to that derived using more linear detectors such as InSb or DTGS ones. We therefore retained a 3% value for the contribution of systematic errors. For each band, we have then estimated the accuracy of the integrated band intensities from the uncertainties on the individual experimental parameters, i.e. ε si (sample purity), ε t (temperature), ε p (pressure), ε pl (path length), ε fit (standard deviation from fit) and ε sys (systematic errors), assuming that these uncertainties are uncorrelated:

Using this equation and the experimental relative uncertainties involved, the estimated overall uncertainty (estimated accuracy) for each of the integrated band intensity measured was calculated. They vary from 3.4 to 6.2% but, on average, the estimated accuracy on the reported integrated band intensity is equal to 4%, and this for the three temperatures considered in this work.

Variations of the absorption spectra and of integrated band intensities with respect to temperature

Examples of the absorption cross sections measured in this work, for all observed bands are presented in Figures 5 and6 for T = 199 and 300 K. As can be seen, when decreasing the temperature, the absorption peaks increase and sharpen leading to a reduction of the absorption widths of the bands. One can see this effect very well in the spectral region of the ν 5 band. It is mainly due to the variations in the vibrational and rotational populations with temperature. It is also worth noticing that hot bands persist in the spectra even for the coldest temperature since their contributions are still amounting to 18% at 199 K (see Table 1).

Figure 7 presents the temperature dependence of the measured integrated band intensities for the various phosgene bands. Changing the temperature from 300 to 199 K leads to non-negligible variations in the integrated band intensities with maximum differences reaching 31% for the weak ν 2 +ν 5 band. Even if these variations are higher than our expected measurement uncertainties, it is difficult to assess if the integrated band intensities are temperature-dependent or not. Indeed, for the ν 1 , ν 2 /ν 4 and ν 5 bands, the variations are less than 10%. Such variations have been observed in the literature for similar heavy molecules. procedures can be used: if only the pure absorber gas is introduced in the cell, given the needed pressure range, thermomolecular effects can distort the pressure reading, especially for relatively low pressures (a few hPa). Instead, if, as we chose to do in this study, the absorbing gas is introduced from a previously done mixture with an inert buffer gas, the error on the total gas pressure will be relatively small, as the total pressure is high (hundreds of hPa range), however the partial pressure of the absorbing gas can no longer be directly measured.

This may be a problem if the molecule can adsorb on the cell walls, an effect which depends on the cell wall nature and probably, to some extent, on the temperature. This probably affects our measurements, especially the first two data points at low partial pressures (1.894 and 2.439 hPa) but, above these pressures, data points present a consistent linear dependence. We have retained all data to illustrate this difficulty, but use only the higher pressure data points for the determination of band region integrated intensities. The second comment relates to the apparent temperature dependence of the band region integrated intensities. It is expected that the fundamental band region integrated intensities should not present such dependences, as it is expected that fundamental bands (ν i ) and their associated hot bands (2ν i -ν i and so on)

should have the same transition moment (neglecting anharmonicity). In practice, it seems not always to be the case, as other transitions can fall by chance in the vicinity of the fundamental transitions and have different temperature dependences than that of the fundamentals. These can be difference bands of the type (ν i + ν j ) -ν k . For instance, if the weak features observed on the high frequency side of the ν 1 band is obviously not a cold band (perhaps ν 5 +ν 2 -ν 6 ), other small features show a marked, although different temperature dependences. The feature near 977 cm -1 is likely to be the ν 1 -ν 5 difference band and its linear decrease with temperature is consistent with the calculated decrease in v 5 = 1 state population. The weak feature near 827 cm -1 is likely to be the ν 1 -(ν 2 +ν 6 ) as it clearly presents a higher order temperature dependence. This is consistent with the expected decrease in v 2 = v 6 = 1 vibrational state population and the experimentally observed temperature decrease for this spectral region (see

  photochemical decay of carbon tetrachloride (CCl 4 ), as well as by oxidation of its tropospheric parent compounds [5,6]. With a lifetime of about 70 days, tropospheric phosgene is mostly removed by reaction in water droplets or by ocean deposition, whereas stratospheric phosgene, due to its weak absorption in the UV, is slowly oxidized to form ClO x and has a lifetime of several years [1,3]. It is then one of the chlorinated compounds involved in stratospheric ozone depletion. It has been observed by stratospheric balloons using the solar occultation technique [6-8] and by satellite instruments such as the Atmospheric Chemistry Experiment -Fourier Transform Spectrometer (ACE-FTS) and the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) for Cl 2 CO vertical distribution retrievals [4,9,10]. For infrared remote sensing, the 11.8 µm (830-860 cm -1 ) spectral region is mainly

  knowledge, three studies have published integrated band intensities of Cl 2 CO at room temperature [23-25]. However, according to Toon et al. [8], further laboratory measurements of the integrated band intensities of the mid-infrared bands of Cl 2 CO as a function of temperature are needed, on the one hand to resolve the large discrepancies between the existing room temperature measurements [23-25], and on the other hand to understand the contribution of hot bands. Indeed, according to Toon et al. [8], not taking into account properly the hot bands contributions leads to a systematic error of ~20% on the concentrations of phosgene retrieved from atmospheric spectra. In order to help solving the questions raised previously the present work reports at temperatures of 199, 250 and 300 K accurate integrated band intensities derived from laboratory spectra covering the spectral range 525-2400 cm -1 . In the next chapters we describe the measurements, estimate their accuracy, and compare the results with the previous sets of measurements at room temperature [23-25].

[ 26 ] 3 . Integrated band intensities and uncertainty assessment 3 . 1 5 and ν 1 .

 2633151 was used. The calculated CO 2 partial pressure was then used to determine the partial pressure of phosgene in the sample. Pressure values for all the samples are reported in Table 2. Determination of the integrated band intensities Absorption measurements were performed for various pressures to determine the integrated intensities for all bands of phosgene indicated above i.e. ν 2 /ν 4 , ν 5 , ν 2 +ν 6 , ν 2 +ν 5 , 2ν The derived integrated absorbances (cm -1 ) are plotted as a function of the Cl 2 CO partial pressure and show a clear linear dependence to pressure (Figure 3). For all selected spectral regions, high correlation coefficients of ~0.99 are obtained. The slope of the linear fit (cm -1 /hPa) is used to determine the integrated band intensity (cm -1 /molecule.cm -2 ), according to the following equation: Where = 2.68676×10 19 /cm 3 is the Loschmidt number at standard conditions, = 1013.25 hPa, = 273.15 K, l = 5.10 cm , is the measurement temperature in K and (cm - 1 /hPa) the spectral integrated intensity measured in decadic absorbance and normalized to the partial pressure. Wavenumber integration limits, retrieved values of the integrated band intensities for the three measurement temperatures, 199, 250 and 300 K and corresponding uncertainties are listed in Table 3. Also for comparison, values derived at 300 K, from the PNNL database [25], Hopper et al. [23] and Lovell and Jones [24] are given (See section 4 for a detailed discussion).

3.5% (see Figure 4 )

 4 for all temperatures. However, estimating the accuracy of the measured integrated band intensities requires considering the uncertainties effecting all physical parameters and contributions from possible systematic errors [25, 27] together with the standard deviation of the fits. These various sources of error and their associated uncertainties, expressed relative to the integrated band intensities, are given in Figure 4, for all the bands studied in this work at 300K. It is clear that the systematic errors and the standard deviations of the fits are the main sources of errors. These contributions are also predominant for the remaining two low temperatures which are not illustrated here. The dominant contributions to systematic errors come from the location of the 100 and 0% transmission photometric levels, as well as from the electronic and detector nonlinearities [25, 27].

For

  example, for the N 2 O 5 molecule, deviations of 6 to 10% are observed in the values of the integrated band intensities when the temperature varies from 233 to 296 K [28]. For ClONO 2 , deviations of 6 to 31% are observed in the values of the integrated band intensities when the temperature varies from 213 to 296 K [29]. The largest difference (31%) was observed for a very weak band, as in our study. Finally, for the BrONO 2 molecule, a deviation of 27% was observed in the value of the relative absorption cross sections when the temperature varies from 213.9 to 293.3 K [30]. More generally, the precise determination of integrated band intensities for reactive or adsorbing volatile organic compounds at variable temperatures poses specific experimental challenges and our results calls for a few remarks. Two

4 .

 4 Comparison with previous works at room temperature Table 3 reports the phosgene integrated band region intensities derived in this work for the three measurement temperatures, 199, 250 and 300 K and values derived at 300 K from the PNNL database [25], Hopper et al. [23] and Lovell and Jones [24]. Our results agree reasonably well (better than 9%) with the PNNL data [25], and are consistent within the combined measurement uncertainties (Figure 8), except for the weak bands ν 2 /ν 4 , ν 2 +ν 5 and ν 2 +ν 6 where the PNNL results are higher by 14.1%, 16.4% and 14.5%, respectively. On the other hand, an excellent agreement (better than 4%) is observed with the values from Hopper et al. [23]. Finally a strong disagreement is observed when comparing with the values from Lovell and Jones [24] which exhibit 20 to 43% discrepancies with our results and also with those of Hopper et al. [23] and the PNNL data [25]. As a conclusion one can say that three measurements are consistent and therefore that the results of Lovell and Jones [24] can likely be discarded.5. SummaryIntegrated band intensities of phosgene have been measured in the 525-2400 cm -1 spectral region at 199, 250 and 300 K. The measurements were performed by recording N 2 -broadened spectra at a resolution of 0.1 cm -1 , over a wide range of pressures, using Fourier transform infrared spectroscopy. Our results at room temperature were compared with previously published measurements. The results from this work show an excellent agreement (better than 4%) with the results of Hopper et al.[START_REF] Hopper | Vibrational intensities. XVI. COF 2 , COCl 2 , and COBr 2[END_REF] and are consistent to within their combined measurement uncertainties with the PNNL data[START_REF] Sharpe | Gas-phase databases for quantitative infrared spectroscopy[END_REF]. On the contrary they disagree with the measurements of Lovell and Jones[START_REF] Lovell | Infrared intensities and bond moments in phosgene and thiophosgene[END_REF] with a maximum difference of ~43%. With respect to the temperature, the integrated band intensities show a non-negligible dependence (7.5 to ~31%), when decreasing temperature from 300 to 199 K, probably due to the presence of overlapping weak features arising from neighboring bands such as overtones, combination/difference and hot bands. It would be important to check this assumption by recording high resolution pure samples of Cl 2 CO over the range of temperatures used for the present measurements.

  line-by-line analysis from high resolution spectra at variable temperature, cold enough to eliminate hot bands and preferably with monoisotopic molecules. This goes, however, much beyond the scope of our present study. To summarize, the integrated band intensity of each band system should in principle be independent of temperature. However, overlapping weak features arising from neighboring bands such as combination/difference and hot bands can be responsible for small variations of the integrated intensities with temperature. In fact, each spectral region presents a specific temperature dependence, which cannot be clearly understood in absence of a complete spectral assignment of all unresolved weaker features. To investigate this dependence, it would be necessary to record pure Cl 2 CO high resolution spectra at low and room temperatures.

	Figure 5, middle trace). This was discussed, for instance, for measurements of acetone at
	variable temperatures [31]. A later study alleviated this problem by simply renormalizing all

integrated band intensities on those deduced from PNNL spectra measured at room temperature

[START_REF] Harrison | Mid infrared absorption cross sections for acetone[END_REF]

. We chose to avoid this procedure to propose a new, independent measurement. It is however clear that a better, more precise intensity determination should proceed through a

Acknowledgments

The authors are grateful to Dr. A. Jolly (LISA) for helpful discussions.


François Boige
email: boige@laplace.univ-tlse.fr

Frédéric Richardeau

Stéphane Lefebvre

GLOBAL ELECTRO-THERMAL MODELLING AND CIRCUIT-TYPE SIMULATION OF SIC MOSFET POWER DEVICES IN SHORT-CIRCUIT OPERATION FOR CRITICAL SYSTEM ANALYSIS

Keywords: electrothermal simulation, model validation, MOSFET, short-circuit, silicon carbide

come

Global electro-thermal modelling and circuit-type simulation of SiC Mosfet power devices in short-circuit operation for critical system analysis

INTRODUCTİON

Nowadays, due to mass production from major manufacturers, SiC power MOSFET is gradually used in several power applications such as inverter for PV panels and grid, UPS and railway traction drive. In that kind of application, short-circuit (SC) faults can appear leading to large transient power dissipation and a risk of failure. Examples of these faults-types are available on Fig. 1. In order to develop a fast and robust protection against these faults, studies on the ruggedness of SiC MOSFET under short-circuit have been carried out [START_REF] Chen | Study of short-circuit robustness of SiC MOSFETs, analysis of the failure modes and comparison with BJTs[END_REF]- [START_REF] Nguyen | Gate Oxide Reliability Issues of SiC MOSFETs Under Short-Circuit Operation[END_REF]. Nevertheless, none electro-thermal models describe globally and accurately SiC MOSFET behaviour under short-circuit faults. Particularly, the gate-leakage current increases significantly during short-circuit, as depicted in Fig. 2, and the saturation current behaviour is unusual compared to silicon (Si) devices. The significant gate-leakage current observed during SC can directly affect the gatevoltage control of the device through the voltage drop across the gate driver resistor. Besides, the temperature-junction estimation in short-circuit mode appears as an important topic to determine the maximum delay time of the gate drive protection circuit. The goal of the proposed study is to develop a new global and macro-model for SiC MOSFET and particularly applied, in this paper, for the SCT2080KE device from ROHM™. Emphasis has been put on finding and fitting accurate models for transient gate leakage current during SC phase (IG), drain-source saturation current (IDsat) and threshold voltage (VGS(th)) according to complex and several parameters, particularly using estimated maximum junction temperature (Tj) at the top of the chip. The original contribution of the paper to the classical MOSFET model is highlighted in Fig. 14 at the end of the paper. Combining the gate-current and drain-source current behaviour, the aim of this paper is to present how the electro-thermal modelling methodology of a SiC MOSFET under extreme stress has been built and evaluated for a practical use. In Section 2, after an introduction, a first dipole-type gate-leakage current model is presented function of four parameters (dissipated energy, gate driver voltage, gate resistance and drain to source voltage). In Section 3, a 1D estimation model of the temperature distribution inside the chip is described. In Section 4, with the results from this estimation, a physical model is fitted in order to simulate the drain-source saturation current evolution during the SC phase. Finally, in Section 5, the of the global behaviour, using PLECS software, are presented and compared with experimentations.

GATE-LEAKAGE CURRENT MODELLING

MODEL DETERMİNATİON

After few micro-seconds of short-circuit, a great transient current surge appears across the thin MOSFET gate oxide (50nm). This current is most likely caused by a highly temperature dependant tunnelling considering the important energy dissipated inside the chip (up to 10J/cm² in 12µs). However, the SC duration is much lower than the thermal time constant of the component (𝜏 𝑑𝑖𝑒 ≅ 350𝜇𝑠 ≫ 𝜏 𝑆𝐶 = 12𝜇𝑠). In these conditions, dissipated energy can be considered instead of the junction temperature in order to explain the time dependence of the gate leakage current during SC. Furthermore, as depicted in Fig. 2, the gate-leakage current dynamic also depends on the voltage bias imposed by the gate-driver (V buffer ), drain-source voltage (V DS ) (not displayed) and the external gate resistor (R G ) (not displayed). The aim of this section is to explain how the influence of each parameter has been identified and modelled. In order to reduce the complexity of the study, only gate current increase, along the sequence n°3 in Fig. 3, is modelled. As a result, the surge occurs from a fixed energy (0.8J or 5.9J/cm² of active chip area for all measurements) so the model will be valid from this value only. For the other parameters, the model validity has to be large: A "dipole-type" model has been chosen (fig. 15) as an equivalent non-linear controlled current-load FIGleak between gate and source electrodes supplied by the driver voltage. This voltage source is made of a voltage source (Vbuffer) and the series gate resistor (RG). Only the current-load equation has to be determined. Of course, all experimental data follow source characteristic by a single crossing point (circles in Fig. 3). However, in Fig. 3., it is remarkable that for a given energy and different buffer bias a load straight line is identifiable and is significantly energy (temperature) dependant. In order to analyse the energy influence on the currentload equation, load straight lines are displayed in 3D representation (IG, VGS, En) in Fig. 4. This representation is basically the same as depicted in Fig. 3. But the load characteristic of the gate driver can be visualised as a surface. Therefore, a polynomial surface (1) has been chosen and fitted with MATLAB™ to model the current-load variation function of VGS and dissipated energy (En) on model validity area.

V DS ∈ [400V, 600V],V buffer ∈ [18V, 21.5V], Energy ∈ [0.8J, 1.5J].
I g = a 1 (𝐸 𝑛 ). V GS + a 2 (𝐸 𝑛 ) (

RESULTS

In

TEMPERATURE ESTİMATİON

MODEL CONSTRUCTİON

During a short-circuit event, strong electrical field is distributed along the low doped epi-layer and a strong current density flows through the chip. This particular configuration induces a high and fast transient temperature gradient inside the chip within few micro seconds. With encapsulated components, the temperature gradient is not measurable in real time because the external case temperature do not represent the junction temperature of the chip [START_REF] Castellazzi | Short-circuit robustness of SiC Power MOSFETs: Experimental analysis[END_REF]. In order to fit physical models, such as temperature dependant electron mobility or threshold voltage, an accurate temperature estimation is needed. All these physical models are required for the drainsource current proposed model. Moreover, the thermal model must be accurate and low simulation time consuming. To do that, the best compromise is a one dimensional transient thermal model across the chip depth, as shown in Fig. 6. During short-circuit transient, the dc bus voltage, VDC, is applied across the power device and leads to a depletion layers with a width xp in the P-well (highly doped Na) and xn in the N-drift region (low doped Nd). The resulting electric-field distribution in the space charge region is given by [START_REF] Wang | Temperature-Dependent Short-Circuit Capability of Silicon Carbide Power MOSFETs[END_REF], neglecting the effect of the current on the depleted charges.

for: x p  𝑥  x j E(𝑥) = qN a ϵ s (𝑥 + Δx p -x j ) (4)

for :

x j  𝑥  x n E(𝑥) = - qN d ϵ s (𝑥 -Δx n -x j ) (5)
for : 𝑥  x p and for 𝑥  x n E(𝑥) = 0 (6) Where ϵ s is the dielectric constant for 4H-SiC material and q is the electron charge. VDS is the electric field integral along the space dimension.

V ds = E(x j ) • Δx n 2 + E(x j ) • Δx p 2 ( 7
)
Doping values and dimensions in Fig 6 . are chosen according to the literature data [START_REF] Romano | A Comprehensive Study of the Short-Circuit Ruggedness of Silicon Carbide Power MOSFETs[END_REF]. It results a Δx n = 7.6µm and Δx p = 8nm at VDS=600V.

The 1D approximation of the heat flux profile in the volume Q(x, t) [W/m 3 ] is depending on the electric field [V/m] function of the variable "x" and the drainsaturation current function of the variable "t" with the following formulation:

Q(𝑥, 𝑡) = E(𝑥) • I DS (𝑡) S chip (8)
Where I ds (𝑡) is the experimental short-circuit current and S chip (13.6mm 2 ) the active chip area.

Fig. 6 depicts the boundary conditions. The top surface of the chip is assumed adiabatic due to the thick plastic coating around the die and the bottom is considered at a constant temperature (ambient temperature) due to the fact that the heat diffusion depth is much lower than the thickness of the die at the end of the short-circuit stress and has been confirmed by changing the boundary conditions (not presented here). This 1D model is implemented with COMSOL® software.

RESULTS

In fig. 7 and8, Tj(x,t) the temperature distribution in the die and the maximum temperature are presented function of the short circuit duration. Estimated temperature increases up to 1300K for 12µs of SC duration. Estimated temperature are representative to results already observed in the literature [START_REF] Wang | Temperature-Dependent Short-Circuit Capability of Silicon Carbide Power MOSFETs[END_REF], [START_REF] Duong | Electro-thermal simulation of 1200 V 4H-SiC MOSFET short-circuit SOA[END_REF] but obtained with other estimation methods which tend to confirm the results. Moreover, temperature is maximal for x=xj where the electrical field is also maximum. For the following thermal modelling, only the maximum estimated temperature will be used.

DRAİN-SOURCE CURRENT MODELLİNG

Fig.

Saturation current determination methodology

To accurately model the drain-source saturation current as it is shown in Fig. 2b, a representative physical model is needed. Available models depend of a lot of physical parameters that are not known over a wide temperature range. In order to estimate them, data from experiment are needed. Methodology is summarized in Fig. 9. Different models have been tested but finally, the following model has been chosen due to its classical form [START_REF] Baliga | Modern power devices[END_REF]:

I Dsat (𝑇 𝑗 ) = K gm (𝑇 𝑗 ) 2 • (𝑉 𝑔𝑠 -V gs(th) (𝑇 𝑗 )) 2 (9)
With

K gm (𝑇 𝑗 ) = μ n (𝑇 𝑗 ) • C ox • Z L (10)
Where, 𝐶 𝑜𝑥 is the oxyde capacitance, Z and L are the width and length of the channel respectively. The carrier mobility in the channel, μ n , is at least the parallel combinations of three different mobility functions with different sensitivity toward temperature [START_REF] Potbhare | A Physical Model of High Temperature 4H-SiC MOSFETs[END_REF]. The Matheissen's rule gives:

μ n (𝑇 𝑗 ) = μ a0 ( 𝑇 𝑗 T 0 ) a 1 + μ a0 μ b0 ( 𝑇 𝑗 T 0 ) b-a + μ a0 μ c0 ( 𝑇 𝑗 T 0 ) c-a (11)
Where μ a0 , μ b0 , μ c0 are mobility parameters at ambient temperature and a,b,c are sensitive thermal coefficients that have to be estimated to obtain a wide range temperature model. Otherwise, literature [START_REF] Baliga | Modern power devices[END_REF], gives a linear threshold voltage dependence with temperature:

V gs(th) (𝑇 𝑗 ) = x 1 • 𝑇 𝑗 + x 2 (12) Where x1 and x2 are unknown parameters. All the previous parameters are determined in two steps. At first, K gm and V gs(th) are determined function of the temperature. To do that, an algorithm has been developed to find the best pair of parameters (x1, x2) which minimize the difference between calculated K gm with the equation ( 9) for different experimental data. Results are displayed in Fig. 10, with the couple (x1=-0.025, x2=15.4) and Kgm remain the same for different buffer bias (19V, 20V, 21V). In a second time, the mobility model ( 11) is fitted on K gm profile previously determined using MATLAB's function fit. Results are evaluated in Fig. 11. Unfortunately, the model is not robust to Vds variation. Other models have been tested to take this parameter into account but without success. The next step is to take into account a more complex model and use a metaheuristic method, such as genetic algorithm, to find the optimum and to precisely take into account variation of the DC bus voltage.

IMPLEMENTATION IN PLECS

The models have been implemented in PLECS™, a software for simulations of electrical circuits. PLECS presents the advantages of being easily pluggable with MATLAB and of quickly simulating complex topologies. The aim of this section is to present all the previous results implemented into a single circuit-type model as summarized in Fig. 14. The classical MOSFET model (in Fig. 14 Where 𝜌 is the density of the SiC material (kg.m -3 ), S is the chip area (m 2 ), P(t) is the power dissipated into the chip (W), c the specific heat of SiC (J.Kg - 1 .K -1 ), α the thermal diffusivity of SiC (m².s -1 ) and K is a dimensionless corrective factor. By fitting K factor, we show that this simple model can give results close to those obtained by the COMSOL model. The results of the whole simulation are shown in Fig. 13. Our objective is fulfil: the simulation is well working and will be able to be used in more complex circuit, for example, short-circuit detection and more generally concerning the safety analysis of critical system. Moreover, the model clearly represent the effect of the buffer voltage bias variation. All the modelling procedure is automatized in MATLAB and this model could be easily obtained for other SiC MOSFET components. However, strong coupling with PLECS circuit and COMSOL 1D model should be implemented and drain saturation current robustness to VDS voltage variation has to be also investigated in future works.

CONCLUSION

The aim of this paper was presenting a complete and new modelling methodology for SiC power MOSFET device in high-current and hightemperature conditions. Firstly, the gate transient leakage current is modelled. In a second step the temperature is estimated inside the chip and finally a physical model of the drain saturation current is given. After validating these different models all the previous models are embedded into a circuit-type MOSFET model. This "block" could be used, as it is, to simulate fault into complex circuits using SiC MOSFET components like a motor drive, for example. Moreover, the proposed model enhancement does not add lot of calculation time and is fully compatible with classical numerical or circuit-type MOSFET models. Results are presented for the ROHM component SCT2080KE but the method is clearly practical for other manufacturers such as CREE or STMicroletronics.

Fig. 1 .

 1 Fig. 1. (a) inverter-leg cross-conduction (Type I SC). (b) Internal short-circuit by electrical breakdown of one of the two leg-switches (Type II SC).

Fig. 2 .

 2 Fig. 2. Experimental waveforms, from a lab. set-up, during SiC MOSFET short-circuit. (a) gate-leakage current and dissipated energy inside the chip. (b) drain current during SC (saturation current) and drain-source voltage.

Fig. 3 .

 3 Fig. 3. Measured data points at the intersection between the driver characteristic and the equivalent current-load characteristic FIGleak in the plane (IG,VGS) for RG=47Ω, VDS=600V and different dissipated energy values.

Fig. 4 .

 4 Fig. 4. Measured data points at the intersection between the voltage driver and the current-load characteristic in the space (IG,VGS, En) for RG=47Ω and VDS =600V.

  order to determine the model robustness, two different tests were made. Firstly, the proposed model was compared with experimental results which are different from those used for the model fitting (Fig.5.(a)). Secondly, the proposed model was compared with data obtained from another component SCH2080KE, i.e. for new experimental results (Fig.5.(b)). The simulation results are, in each case, very close to experimental data. In conclusion, it can be considered that the numerical proposed model is clearly validated and really robust in its full domain of validity.

Fig. 5 .

 5 Fig. 5. Model (red lines) and experimental (blue lines) results. (a) robustness results obtained from test data which are not used for the model fitting but with the same device. (b) robustness Results obtained from test data not used for the model fitting and with another device.

Fig. 6 .

 6 Fig. 6. 1D thermal model of the SiC MOSFET.

Fig. 7 .

 7 Fig. 7. Estimated T(x, t) temperature distribution vertically along the component for a short-circuit at Vds=600V, Vbuffer=21V, Tinit (t=0)=298,15K.

Fig. 8 .

 8 Fig. 8. 1D Estimation of the maximum junction temperature Tj (xj, t) for different buffer-voltages bias at Vds=600V, Tjinit (t=0s)=298,15K.

Fig. 11 .

 11 Fig. 11. Saturation current results obtained after fitting for different gate bias (experimental results in blue lines and simulation results in red lines).

Fig. 10 .

 10 Fig. 10. Kgm and Vgs(th) fitting based on the 1D thermal model and experimental data acquired with different gate bias and VDS=600V.

Fig. 12 .

 12 Fig. 12. Proposed SiC MOSFET model in PLECS software.

  black line) is completed (in red line) adding the previous gate leakage current modelling by "FIGleak", the drainsource current modelling by "FIdsat" and finally the gate threshold voltage by "FVGth". The different models are performed by PLECS function block or by code-C script. The PLECS schematic is displayed in Fig.12. Today, the co-simulation between COMSOL and PLECS has not been implemented. In order to estimate the junction temperature inside the chip without COMSOL the following simplified model is proposed based from [9]: Δ𝑇 𝑗 (t) 𝐾 𝑆. 𝜌. 𝑐. √𝜋. 𝛼 1 √𝑡 ∫ 𝑃(𝜏) 𝑑𝜏 (13) 𝑡 0

Fig. 13 .

 13 Fig. 13. Simulated waveforms from the proposed PLECS circuit-type model compared with measured ones for different buffer voltages bias (@Vds=600V, Tjinit=298,15K (25°C) and Rg = 47Ω)

1 )

 1 With a 1 (𝐸 𝑛 ) = p 01 + p 11 𝐸 𝑛 + p 21 𝐸 𝑛

	𝑝00	174.2161 -0.5505. 𝑉 𝑑𝑠
	𝑝10	192.8417 -0.5720. 𝑉 𝑑𝑠
	𝑝01	-11.8436 + 0.0356. 𝑉 𝑑𝑠
	𝑝20	52.8609 -0.1482. 𝑉 𝑑𝑠
	𝑝11	-13.7143 + 0.0386. 𝑉 𝑑𝑠
	𝑝30	-0.9317 + 0.0026. 𝑉 𝑑𝑠
	𝑝21	-4.1480 + 0.0111. 𝑉 𝑑𝑠
		2	(2)
		a 2 (𝐸 𝑛 ) = p 00 + p 10 𝐸 𝑛 + p 20 𝐸 𝑛 2 + p 30 𝐸 𝑛 3 (3)
		(p 00 , p 01 , … ) are the estimated factors.
		This "response surface" approach has two main
		advantages. On one hand, MATLAB™ polynomial
		surface fitting algorithm from the fit function is
		powerful and robust. On an other hand, the

polynomial function and dipole type model makes it easily computable and "pluggable" into a commercial-type software. At this point, the model represents effects of Vbuffer, RG and En on gate leakage current. However, VDS variation is not yet added in the model. To do that, factors (p 00 , p 01 , … ) are, themselves, estimated for different VDS bias and precalculated in a table. It appears that these factors are linearly varying with VDS. Then, the variation has been estimated and added to the model.

I. Estimated factors function of Vds

ACKNOWLEDGEMENTS

This research work received financial support from the French National Research Agency (ANR).