Johannes Schiffer
email: j.schiffer@leeds.ac.uk

Denis Efimov

Romeo Ortega

Nikita Barabanov

An input-to-state stability approach to verify almost global stability of a synchronous-machineinfinite-bus system

Keywords: Power systems, control systems, nonlinear dynamics Power system dynamics, power system stability, nonlinear control systems, input-to-state stability

Conditions for almost global stability of an operating point of a realistic model of a synchronous generator with constant field current connected to an infinite bus are derived. The analysis is conducted by employing the recently proposed concept of input-to-state stability (ISS)-Leonov functions, which is an extension of the powerful cell structure principle developed by Leonov and Noldus to the ISS framework. Compared to the original ideas of Leonov and Noldus, the ISS-Leonov approach has the advantage of providing additional robustness guarantees. The efficiency of the derived sufficient conditions is illustrated via numerical experiments.

Introduction (a) Motivating challenges

Modern industrialized societies heavily rely on the use and steady supply of electric energy. Therefore, it is of paramount importance to guarantee the stable, reliable and efficient operation of our power systems. However, the latter are bulk, complex and highly nonlinear systems that are continuously subjected to a large variety of disturbances and contingencies [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF]. Hence, ensuring the stability of a power system is a daunting task, which has been at the core of power system operation since its early days in the 1920s [START_REF] Rüdenberg | Elektrische Schaltvorgänge und verwandte Störungserscheinungen in Starkstromanlagen[END_REF][START_REF]Transient Performance of Electric Power Systems[END_REF]-see [21, Chapter 1] for a review of the research history on power system stability analysis.

Over the past decades, the complexity of this task has continuously exacerbated due to the penetration of large shares of volatile renewable energy sources. While this is a very desirable development from environmental and societal factors, the reliable and efficient integration of high amounts of renewable energy sources represents a major technical challenge for the operation of power systems [START_REF] Farhangi | The path of the smart grid[END_REF][START_REF] Winter | Pushing the limits: Europe's new grid: Innovative tools to combat transmission bottlenecks and reduced inertia[END_REF]. In particular, the increasing presence of such units results in power systems more frequently operating closer to their security margins [START_REF] Farhangi | The path of the smart grid[END_REF][START_REF] Winter | Pushing the limits: Europe's new grid: Innovative tools to combat transmission bottlenecks and reduced inertia[END_REF]. As a consequence, the problem of power system stability can be expected to become even more relevant in the coming years. This makes it important to derive easily and quickly verifiable analytic conditions for stability, which is the topic addressed in the present work.

A further motivation for our work is the fact that-despite the rich and extensive literature on power system stability-even today many basic questions remain open. This applies to the analysis of vast networks, but also to individual network components and is explained by the complexity of their dynamics. Thus, researchers and practitioners alike are forced to invoke several assumptions simplifying the mathematical tasks of stability analysis and control design. Amongst the most prevalent assumptions are neglecting fast dynamics [START_REF] Schiffer | A survey on modeling of microgrids-from fundamental physics to phasors and voltage sources[END_REF][START_REF] Venkatasubramanian | Fast time-varying phasor analysis in the balanced three-phase large electric power system[END_REF][START_REF] Wang | Analysis and Damping Control of Power System Low-frequency Oscillations[END_REF], constant voltage amplitudes and small frequency variations [START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF]Chapter 11], [START_REF] Wang | Analysis and Damping Control of Power System Low-frequency Oscillations[END_REF]. Invoking such assumptions allows to derive reduced-order synchronous generator (SG) models [START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF]Chapter 11] and employ algebraic line models [START_REF] Schiffer | A survey on modeling of microgrids-from fundamental physics to phasors and voltage sources[END_REF][START_REF] Venkatasubramanian | Fast time-varying phasor analysis in the balanced three-phase large electric power system[END_REF][START_REF] Wang | Analysis and Damping Control of Power System Low-frequency Oscillations[END_REF], which significantly simplifies the analysis. Unfortunately, most of the employed assumptions are not physically justifiable in generic operation scenarios. A paramount example for this is the common approximation of the motion of the machine rotor, i.e., the swing equation, in terms of mechanical and electrical power instead of their corresponding torques. However, as this approximation is only valid for small frequency variations around the nominal frequency, it is not an adequate representation of the true physical system when the latter is subjected to large disturbances [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Kundur | Power System Stability and Control[END_REF][START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF]; see also the discussions in [START_REF] Caliskan | Uses and abuses of the swing equation model[END_REF][START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF][START_REF] Fiaz | A port-hamiltonian approach to power network modeling and analysis[END_REF][START_REF] Venezian | A warning about the use of reduced models of synchronous generators[END_REF]. In summary, a main limitation of current power system stability analysis methods is that they are based on reduced-order models only valid within a limited range of operating conditions. A direct consequence of this is that any stability assessment relying on such reducedorder models may lead to erroneous predictions [START_REF] Caliskan | Uses and abuses of the swing equation model[END_REF][START_REF] Venezian | A warning about the use of reduced models of synchronous generators[END_REF], which-in the worst case-could lead to a blackout.

(b) Existing literature and contributions

In the present work, we perform a global stability analysis of a classical and very well known system called the single-machine-infinite-bus (SMIB) model [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Kundur | Power System Stability and Control[END_REF]. Thereby, the main distinction to existing results [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF] is that, motivated by the above discussion, we consider a fourth-order nonlinear SG model derived from first principles-instead of the usual simplified SG models. A thorough analysis and understanding of the SMIB scenario is a fundamental step in the development of a systematic mathematical framework for the analysis of more complex power system dynamics under less restrictive and more realistic assumptions as in [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF][START_REF] Fiaz | A port-hamiltonian approach to power network modeling and analysis[END_REF] or in the classical literature [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Kundur | Power System Stability and Control[END_REF][START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF]. Furthermore, the insights that can be gained from our analysis can provide valuable practical conclusions for the design, control and operation of SGs as well as for the implementation of an emerging strand of control concepts that aim at operating gridconnected inverters such that they mimic the dynamics of conventional SGs [START_REF] Beck | Virtual synchronous machine[END_REF][START_REF] Visscher | Virtual synchronous machines for frequency stabilisation in future grids with a significant share of decentralized generation[END_REF][START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF]. The latter objective also has motivated the related works [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF]. The analysis in [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF] proceeds along the classical lines of constructing an integro-differential equation resembling the forced pendulum equation and, subsequently, showing that the SGIB system is almost globally asymptotically stable if and only if an equilibrium of that forced pendulum equation is almost globally asymptotically stable. In [START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF] the same authors provide slightly simpler conditions for stability resulting from verifying if a real-valued nonlinear map defined on a finite interval is a contraction. But, as stated in [START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF], these conditions are hard to verify analytically. Furthermore, and perhaps more importantly, the geometric tools employed to establish the results in [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF], don't seem to be applicable to a multi-machine power system. In [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF] a scenario similar to that of the SMIB system is analyzed. However, the analysis in [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF] is conducted under very stringent assumptions on the specific form of the infinite bus voltage, as well as the steady-state values of the mechanical torque. A nonlinear analysis of the SMIB system is also conducted in [START_REF] Wang | Nonlinear dynamic analysis of a single-machine infinite-bus power system[END_REF], yet using a power rather than a torque balance in the swing equation and neglecting the stator dynamics. Then the SMIB model reduces to a nonlinear forced pendulum, as also studied in [START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF][START_REF] Efimov | Robustness of delayed multistable systems with application to droop-controlled inverter-based microgrids[END_REF].

To achieve our main objective, i.e. the provision of sufficient conditions for almost global stability of the attractive equilibrium set of the SMIB system, we employ the recently developed framework of input-to-state (ISS) stability for periodic systems with multiple invariant sets [START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF]. Therein, the ISS approach from [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF][START_REF]Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF] is combined with the cell structure principle developed in [START_REF] Gelig | Stability of nonlinear systems with nonunique equilibrium position[END_REF][START_REF] Leonov | Phase synchronisation. theory and applications[END_REF][START_REF]On the boundedness of the trajectories of phase systems[END_REF][START_REF] Noldus | New direct Lyapunov-type method for studying synchronization problems[END_REF][START_REF] Yakubovich | Stability of stationary sets in control systems with discontinuous nonlinearities[END_REF] to derive necessary and sufficient conditions for ISS of periodic systems. This permits to relax the usual sign definiteness requirements on the Lyapunov function and its timederivative by exploiting the periodicity of the system. This relaxation is essential to establish the main result of the present paper. More precisely, the main contributions of the paper are three-fold.

• Provide a sufficient condition for ISS of the equilibrium set of the SMIB system.

• By using this result, establish almost global asymptotic stability of the attractive equilibrium set of the SMIB system, i.e., we show that for all initial conditions, except a set of measure zero, the solutions of the SMIB system asymptotically converge to an asymptotically stable equilibrium point. • Illustrate the efficiency of the proposed conditions numerically via a benchmark example taken from [START_REF] Anderson | Power System Control and Stability[END_REF]Table D.2].

Differently from our previous related work [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF][START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF], in which LaSalle's invariance principle is employed, the present ISS-based approach has the advantage of providing additional robustness guarantees with respect to exogeneous inputs. Both properties-stability and robustness-are highly desirable from a practical point of view to characterize a power system's performance under persistent disturbances, such as load variations [START_REF] Kundur | Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions[END_REF]Section D]. Yet, they are also deemed to be hard to establish, in particular in the global setting considered in the present paper [START_REF] Kundur | Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions[END_REF]Section D].

The remainder of the paper is structured as follows. The SMIB model is introduced in Section 2. The stability analysis is conducted in Section 3. Section 4 presents a benchmark numerical example. The paper is concluded in Section 5 with a summary and an outlook on future work.

Notions of ISS of systems with multiple invariant sets.

For the notions of ISS of systems with multiple invariant sets, we follow the conventions and definitions in [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF][START_REF]Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF][START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF] for dynamical systems evolving on a Riemannian manifold M of dimension n. Hence, the distance of a point x ∈ M from the set S ⊂ M is denoted by |x| S = min a∈S ∆(x, a), where the symbol ∆(x 1 , x 2 ) denotes the Riemannian distance between x 1 and x 2 . Furthermore, |x| = |x| {0} denotes the usual Euclidean norm of a vector x ∈ R n . For further details on properties of dynamical systems and, in particular, systems with multiple invariant sets the reader is referred to [START_REF]Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF][START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF][START_REF] Nitecki | Filtrations, decompositions, and explosions[END_REF] and references therein.

Model of a synchronous generator connected to an infinite bus

The SMIB system is introduced in this section based on [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF][START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF][START_REF] Grainger | Power System Analysis[END_REF]. We employ a generator reference direction, i.e, current flowing out of the SG terminals is counted positively. Following [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF][START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF], we make some standard assumptions on the SG [START_REF] Grainger | Power System Analysis[END_REF][START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF], which are also used in [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF][START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF]: the rotor is round; the machine has one pole pair per phase; there are no damper windings and no saturation effects; there are no Eddy currents; the rotor current i f is a real constant. The latter can be achieved by choosing the excitation voltage such that i f is kept constant, see [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF]. Furthermore, we assume that all three-phase signals are balanced [START_REF] Grainger | Power System Analysis[END_REF][START_REF] Schiffer | A survey on modeling of microgrids-from fundamental physics to phasors and voltage sources[END_REF]. With regards to the SG this is equivalent to assuming a "perfectly built" machine connected in star with no neutral line, see also [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF]. We denote the three-phase voltage at the infinite bus by

v abc := √ 2V    sin(δg) sin(δg -2π 3 ) sin(δg + 2π 3 )    , (2.1)
where V ∈ R >0 is the root-mean-square (RMS) value of the constant voltage amplitude (line-toneutral) and

δg = δg(0) + ω s t ∈ R, (2.2)
with the grid frequency ω s being a positive real constant. Furthermore, the electrical rotor angle of the SG is denoted by δ : R ≥0 → R and the electrical frequency by ω = δ. As illustrated in Fig. 1, δ is the angle between the axis of coil a of the SG and the d-axis. Recall that the rotor current i f is constant by assumption and denote the peak mutual inductance by M f ∈ R >0 . Then the three-phase electromotive force (EMF) e abc : R ≥0 → R 3 induced in the stator is given by [START_REF] Grainger | Power System Analysis[END_REF][START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF]

e abc = M f i f ω    sin(δ) sin(δ -2π 3 ) sin(δ + 2π 3 )    . (2.3)
We denote the stator resistance by R ∈ R >0 and the stator inductance by L = Ls -Ms, where Ls ∈ R >0 is the self-inductance and Ms ∈ R >0 the mutual inductance. In the SMIB scenario, the inductance and resistance of the transmission line connecting the SG to the infinite bus are often included in the parameters L and R, see [START_REF] Anderson | Power System Control and Stability[END_REF]Section 4.13.1]. Also note that L > 0, since in practice Ls > Ms. Then the electrical equations describing the dynamics of the three-phase stator current i abc : R ≥0 → R 3 are given by

L di abc dt = -Ri abc + e abc -v abc . (2.4)
The SMIB model is completed with the electromechanical equations describing the rotor dynamics, given by

δ = ω, J ω = -Dω + Tm -Te, (2.5)
where the electrical torque Te is computed as [START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF] Te = ω -1 i abc e abc .

(2.6) Furthermore, J ∈ R >0 is the total moment of inertia of the rotor masses, D ∈ R >0 is the damping coefficient and Tm ∈ R ≥0 is the mechanical torque provided by the prime mover, which we assume is constant.

For our stability analysis, it is convenient to perform a coordinate transformation that maps sinusoidal three-phase waveforms with constant amplitude and frequency-such as (2.1)-to constant signals. This can be achieved by representing the system (2.4), (2.5) and (2.6) in dqcoordinates via the dq-transformation given in Appendix (a). In the present case, we choose the transformation angle ϕ := ω s t (2.7)

and denote the angle difference between the rotor angle δ and the dq-transformation angle ϕ by

θ := δ -ϕ.
In dq-coordinates, the grid voltage (2.1) is given by the following constant vector (see [START_REF] Schiffer | A survey on modeling of microgrids-from fundamental physics to phasors and voltage sources[END_REF]),

v dq = v d vq = √ 3V sin(δg -ϕ) cos(δg -ϕ) = √ 3V sin(δg(0)) cos(δg(0)) ,
where the second equality follows from (2.2). Likewise, the EMF in dq-coordinates is given by

e dq = e d eq = bω sin(θ) bω cos(θ) , (2.8)
where we have defined

b := 3/2M f i f . (2.9)
Hence, the electrical torque Te in (2.6) is given in dq-coordinates by

Te = ω -1 i abc e abc = ω -1 i dq e dq = b(iq cos(θ) + i d sin(θ)).
(2.10)

Furthermore, with ϕ given in (2.7) and T dq (•) given in (5.1), we have that

dT dq (ϕ) dt i abc = ω s -iq i d ,
see [29, equation (4.8)]. Consequently, by replacing the rotor angle dynamics, i.e., δ, with the relative rotor angle dynamics, i.e., θ, the SMIB model given by (2.4), (2.5) and (2.6) becomes in dq-coordinates

θ = ω -ω s , J ω = -Dω + Tm -b (iq cos(θ) + i d sin(θ)) + d 1 (t), L id = -Ri d -Lω s iq + bω sin(θ) -v d + d 2 (t), L iq = -Riq + Lω s i d + bω cos(θ) -vq + d 3 (t), (2.11)
where Remark 1. In the related works [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF][START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF] the dq-transformation angle ϕ is chosen as the rotor angle, i.e., ϕ = δ. This yields e d = 0 in (2.8). Furthermore, in [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF] a dq-transformation is chosen in which the d-axis lags the q-axis by π/2. In our notation, this coordinate frame can be obtained by setting the transformation angle to ϕ = δ + π. Hence, eq = bω in [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF] and eq = -bω in [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF][START_REF]A method for proving the global stability of a synchronous generator connected to an infinite bus[END_REF]. The analysis reported in the paper can be conducted in any coordinate frame. Yet, the one employed here is independent of the rotor angle of a specific machine and, hence, seems to be more suitable to extend the results to the multi-machine case, which motivates our choice of ϕ.

d(t) = d 1 (t), d 2 (t), d 3 (t) ∈ R 3 is

Almost global stability

This section is dedicated to the analysis of almost global stability of an equilibrium point of the system (2.11). Our analysis is conducted by employing the recently proposed framework of ISS-Leonov functions [START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF]. To this end, we perform the following three steps. First, we construct a suitable error system in Section 3(a). Then, we show in Section 3(b) that the set of equilibria of the For the analysis, we employ the following definition.

Definition 1. An equilibrium point (modulo 2π) of the system (2.11) is almost globally asymptotically stable if it is asymptotically stable and for all initial conditions, except those contained in a set of zero Lebesgue measure, the solutions of the system converge to that equilibrium point (modulo 2π).

(a) Error coordinates

It is convenient to introduce the two constants

c : = b (v 2 d + v 2 q )((Lω s ) 2 + R 2 ), P := 1 c -b 2 ω s R + (Tm -Dω s )((Lω s ) 2 + R 2 ) . (3.1)
Note that (2.9) implies that c is nonzero if the rotor current i f is nonzero, which is satisfied in any practical scenario. The existence of isolated equilibria is a natural prerequisite for any stability analysis. Therefore, we make the following natural assumption that by [4, Proposition 1] implies the existence of two unique steady-state solutions (modulo 2π) of the system (2.11).

Assumption 1. The parameters of the system (2.11) are such that i f > 0 and that |P| < 1 if d = 0 3 for all t ∈ R.

With Assumption 1, we denote one of the two equilibrium points of the system (2.11) by (θ s 1 , ω s , i s d,1 , i s q,1 ) and introduce the error states

x := θ, ω, ĩd , ĩq , θ = θ + θ s 1 , ω = ω + ω s , i d = ĩd + i s d,1 , iq = ĩq + i s q,1 . (3.2)
In the error coordinates ( θ, ω, ĩd , ĩq) the dynamics (2.11) take the form θ =ω,

J ω = -D(ω + ω s ) -b(i s d,1 + ĩd ) sin(θ s 1 + θ) -b(i s q,1 + ĩq) cos(θ s 1 + θ) + Tm + d 1 (t), L id = -R(i s d,1 + ĩd ) -Lω s (i s q,1 + ĩq) + b(ω + ω s ) sin(θ s 1 + θ) -v d + d 2 (t), L iq = -R(i s q,1 + ĩq) + Lω s (i s d,1 + ĩd ) + b(ω + ω s ) cos(θ s 1 + θ) -vq + d 3 (t). (3.3)
Moreover, in the equilibrium point (θ s 1 , ω s , i s d,1 , i s q,1 ) with d(t) = 0 3 , we have from (2.11) that

0 = -Dω s + Tm -b i s q,1 cos(θ s 1 ) + i s d,1 sin(θ s 1 ) , 0 = -Ri s d,1 -Lω s i s q,1 + bω s sin(θ s 1 ) -v d , 0 = -Ri s q,1 + Lω s i s d,1 + bω s cos(θ s 1 ) -vq. (3.4)
Hence, by using (3.4) the system (3.3) becomes in error coordinates θ =ω,

J ω = -D ω -b ĩd sin(θ s 1 + θ) -bi s d,1 (sin(θ s 1 + θ) -sin(θ s 1 )) -b ĩq cos(θ s 1 + θ) -bi s q,1 (cos(θ s 1 + θ) -cos(θ s 1 )) + d 1 (t), L id = -R ĩd -Lω s ĩq + bω s (sin(θ s 1 + θ) -sin(θ s 1 )) + bω sin(θ s 1 + θ) + d 2 (t), L iq = -R ĩq + Lω s ĩd + bω s (cos(θ s 1 + θ) -cos(θ s 1 )) + bω cos(θ s 1 + θ) + d 3 (t), (3.5)
the equilibrium of which is now shifted to the origin. The remainder of this section is devoted to the analysis of this system.

Remark 2.

As can be seen from (2.9), the signs of c and, hence, P defined in (3.1) depend on the sign of the constant rotor current i f . For the subsequent analysis, it is important to know the sign of c, as it determines which of the two equilibria of the system (2.11) is stable. Assumption 1 is made without loss of generality because the analysis for i f < 0 follows verbatim.

(b) Decomposable invariant set containing all equilibria

With Assumption 1, the system (3.5) has a finite and disjoint set of equilibria (modulo 2π) in the sense of [9, Definition 3]. Denote this set by

W = 0 4 ∪ [θ s 2 , ω s , i s d,2 , i s q,2 ] -[θ s 1 , ω s , i s d,1 , i s q,1 ] . ( 3.6)
The lemma below shows that-for bounded solutions-W contains all αand ω-limit sets of the system (3.5)-see [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF][START_REF]Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF][START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF] for details. This feature is instrumental to apply the ISS-Leonov framework [START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF] and establish our main result.

Lemma 1. Consider the system (3.5) with Assumption 1. Suppose that all its solutions are bounded and that

4RD[(Lω s ) 2 + R 2 ] > (Lbω s ) 2 . (3.7)
Then the set W is decomposable in the sense of [9, Definition 3] and contains all α-and ω-limit sets of the system (3.5).

Proof.

In order to show that the set W contains all α-and ω-limit sets of (3.5) and that it is decomposable, note that the electrical dynamics in (3.5) takes the form

L idq = u w + bω sin(θ s 1 + θ) bω cos(θ s 1 + θ) , (3.8)
where we defined

u w := -R -Lω s Lω s -R ĩdq +
bω s (sin( θ + θ s 1 ) -sin(θ s 1 )) bω s (cos( θ + θ s 1 ) -cos(θ s 1 ))

.

Following [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF], consider the function

W (x) = L 2 (u 2 + w 2 ) + J ω2 2 (Lω s ) 2 + R 2 + q( θ), (3.9)
where x is defined in (3.2),

q( θ) = c[-cos(θ s 1 -φ + θ) -sin(θ s 1 -φ) θ + cos(θ s 1 -φ)] (3.10) and 1 φ = arctan Lω s v d + Rvq Lω s vq -Rv d . (3.11)
Clearly, W is continuously differentiable as a function W : R 4 → R (but not as a function W :

S × R 3 → R,
where S denotes the unit circle). With (3.5), a lengthy but straightforward calculation yields

Ẇ = -R[u 2 + w 2 ] -D((Lω s ) 2 + R 2 )ω 2 + ωu Lbω s cos( θ + θ s 1 ) -ωw Lbω s sin( θ + θ s 1 ) + d 1 ω((Lω s ) 2 + R 2 ) + d 2 (Lω s w -Ru) -d 3 (Lω s u + Rw), = -ω u w M ω u w -ω u w M d d,
where

M =     D((Lω s ) 2 + R 2 ) - Lbω s cos( θ+θ s 1 ) 2 Lbω s sin( θ+θ s 1 ) 2 - Lbω s cos( θ+θ s 1 ) 2 R 0 Lbω s sin( θ+θ s 1 ) 2 0 R     , M d =    -(Lω s ) 2 + R 2 0 0 0 R Lω s 0 -Lω s R    .
(3.12) In order to show that W indeed contains all αand ω-limit sets of (3.5), we consider the nominal version of the system (3.5), i.e., with d(t) = 0 3 for all t ∈ R. Since R > 0, the Schur-complement implies that M > 0 if and only if (3.7) holds. Then, we have immediately that (3.7) implies that

Ẇ ≤ 0.
Consequently, by invoking LaSalle's invariance principle [START_REF] Van Der Schaft | L2-Gain and Passivity Techniques in Nonlinear Control[END_REF] we conclude that all bounded solutions of the nominal system (3.5) converge to the set where ω = w = u = 0 ∀t ∈ R.

(3.13) From (3.8) we have that ω = w = u = 0 implies ĩdq constant. Since ω = 0 also implies θ constant, W is the only invariant set under the restriction (3.13) and it contains all αand ω-limit sets of (3.5). Consequently, if d(t) = 0 3 for all t ∈ R, then Ẇ ≤ 0 implies that W is decomposable in the sense of [9, Definition 3], completing the proof.

(c) A condition for ISS

Recall the function q( θ) defined in (3.10) and let

h( θ) = q( θ) -2 θ2 , (3.14)
where is a positive real parameter. Denote

min := inf ∈ R >0 | h( θ) ≤ 0 ∀ θ ∈ R . (3.15)
It is shown in [4, Lemma 1] that c > 0 defined in (3.1) is an upper bound for min , i.e., the infimum in (3.15) is indeed always achievable. The lemma below is useful to establish our main ISS result.

Lemma 2. Let > min . Then for all θ ∈ R and any ε ∈ [0,min ),

-2 θ2 ≤ h( θ) ≤ - ε 2 θ2 .
Proof. By the definition of min we have that

h( θ) := q( θ) -min 2 θ2 ≤ 0.
Hence,

h( θ) = h( θ) - 1 2 ( -min ) θ2 ≤ - 1 2 ( -min ) θ2 ≤ - ε 2 θ2
for any ε ∈ [0,min ). To derive the lower bound for h, we note that since > min ,

-2 θ2 ≤ h( θ) ⇔ 0 ≤ h( θ) + 2 θ2 = h( θ) - 1 2 ( -min ) θ2 + 2 θ2 ≤ - 1 2 ( -min ) θ2 + 2 θ2 ≤ 2 θ2 ,
completing the proof.

To streamline the presentation of our result the following assumption is needed. A physical interpretation of this key assumption is given in the discussion in Remark 4 below. Assumption 2. There exist constants χ > 0 and ψ > 0 satisfying

Q =     χ 2 min 1 0 1 (Lω s ) 2 + R 2 D -χJ 2 -ψ 4 Lbω s 0 Lbω s 4 R -χL 2 -ψ     > 0.
(3.16)

In addition, we borrow the following result from [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF].

Lemma 3.

[4] Consider the system (3.5) verifying Assumptions 1 and 2. Then, all solutions x = θ, ω, ĩd , ĩq of the system (3.5) are bounded.

We are now in a position to present our ISS result.

Theorem 1. Consider the system (3.5) with Assumptions 1 and 2. The system (3.5) is ISS with respect to the set W defined in (3.6).

Proof. Note that Lemma 3 implies that all solutions of the system (3.5) are bounded under Assumption 2. Moreover, we see that Q > 0 in (3.16) reduces to (3.7) for χ = ψ = 0. As Q > 0 for some χ > 0, ψ > 0 by assumption, the condition of Lemma 1 is met. Consequently, the set W is decomposable and contains all αand ω-limit sets of the system (3.5). To this end, we need to construct a function V, that-along all solutions of the system (3.5)-satisfies the following two properties

α 1 (|x| W ) -σ(| θ|) ≤ V (x) ≤ α 2 (|x| W + g) (3.17) and V (x) ≤ -λ(V (x)) + γ(|d|), (3.18)
where g ≥ 0 is a scalar, α 1 ∈ K∞, α 2 ∈ K∞, σ ∈ K∞, γ ∈ K∞, and λ : R → R is a continuous function satisfying λ ∈ K∞ for nonnegative arguments. Following [9, Definition 7], such function V is called an ISS-Leonov function 2 .

Recall W defined in (3.9), h defined in (3.14) and consider the following ISS-Leonov function candidate V : R 4 → R,

V (x) = L 2 (u 2 + w 2 ) + J ω2 2 (Lω s ) 2 + R 2 + h( θ). (3.19)
In the first step of the proof, we show that V satisfies condition (3.17). Note that Lemma 2 implies that for all

x ∈ R 4 , L 2 (u 2 + w 2 ) + J ω2 2 (Lω s ) 2 + R 2 -2 | θ| 2 ≤ V (x) ≤ L 2 (u 2 + w 2 ) + J ω2 2 (Lω s ) 2 + R 2 - ε 2 | θ| 2 ≤ κ|( θ, ω, ĩd , ĩq)|W
for some real κ > 0. Thus, property (3.17) is satisfied.

The second and last step of the proof consists in verifying condition (3.18). Recall M and M d given in (3.12). Differentiating V given in (3.19) along solutions of the system (3.5) yields

V = -    ω u w    M    ω u w    -    ω u w    M d d -θ ω = - x d Qx Q xd Q xd γI 3 :=Q x d -χ L 2 (u 2 + w 2 ) + J ω2 2 (Lω s ) 2 + R 2 - ε 2 | θ| 2 + γ|d| 2 , (3.20)
where

Qx =       1 2 χε 1 2 0 0 1 2 ((Lω s ) 2 + R 2 )(D -χJ 2 ) - Lbω s cos( θ+θ s 1 ) 2 Lbω s sin( θ+θ s 1 ) 2 0 - Lbω s cos( θ+θ s 1 ) 2 R -χL 2 0 0 Lbω s sin( θ+θ s 1 ) 2 0 R -χL 2       , Q xd = 0 1×3 1 2 M d .
We need to verify that Q ≥ 0. Since γ > 0, the Schur complement implies that Q ≥ 0 if and only if

Qx - 1 γ Q xd Q xd = Qx - 1 4γ      0 0 0 0 0 ((Lω s ) 2 + R 2 ) 2 0 0 0 0 (Lω s ) 2 + R 2 0 0 0 0 (Lω s ) 2 + R 2      ≥ 0. (3.21)
2 From (3.17), we see that the sign definiteness requirements of an ISS-Leonov function are relaxed compared to the standard Lyapunov function, because V does not have to be positive definite with respect to the variable θ, i.e., the variable with respect to which the dynamics (3.5) is periodic. Furthermore, the time-derivative of V only needs to be negative definite for positive values of V. See [START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF] for further details.

Furthermore, since χ > 0, ε > 0 and γ > 0, applying the Schur complement also to the matrix in (3.21) implies that Q ≥ 0 if and only if

R - χL 2 - 1 4γ (Lω s ) 2 + R 2 > 0 (3.22)
and

((Lω s ) 2 + R 2 ) D - χJ 2 - 1 4γ ((Lω s ) 2 + R 2 ) 2 - 2 2χε - (Lbω s ) 2 4(R -χ L 2 -1 4γ ((Lω s ) 2 + R 2 )) ≥ 0.
(3.23) Recall from Lemma 2 that 0 < ε <min . Hence, the left-hand side of (3.23) attains its maximum with respect to ε for ε =min . We select this value for ε and proceed by requiring (3.23) to be satisfied with strict inequality. This ensures that there exists a 0 < ε <min for which Q ≥ 0. It is then straightforward to verify that this implies that the expression on the left-hand side of (3.23) attains its maximum with respect to for = 2 min > min . For this choice of (3.23) becomes

((Lω s ) 2 + R 2 ) D - χJ 2 - 1 4γ ((Lω s ) 2 + R 2 ) 2 - 2 min χ - (Lbω s ) 2 4(R -χL 2 ) -1 γ ((Lω s ) 2 + R 2 ) > 0.
(3.24) Hence, by defining the auxiliary variable ψ = ((Lω s ) 2 + R 2 )γ -1 > 0, (3.23) being satisfied with strict inequality is equivalent to positive definiteness of the matrix Q in (3.16). Clearly, Q > 0 also implies that (3.22) is satisfied. Recall that Q > 0 by Assumption 2. Hence, Q ≥ and

V ≤ -χV + γ|d| 2 , (3.25)
which is property (3.18). Consequently, with Assumption 2, V satisfies conditions (3.17) and (3.18), which implies that V is an ISS-Leonov function for the system (3.5). Hence, by invoking [

) it follows that the system (3.5) is ISS with respect to the set W. This completes the proof.

Remark 3. The ISS-Leonov function V defined in (3. [START_REF]On the boundedness of the trajectories of phase systems[END_REF]) is not periodic in θ. Thus, it does not qualify as a usual ISS Lyapunov function [START_REF]Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF]. Yet, this non-periodicity of V in θ is essential for being able to treat the constant terms in the dynamics (3.5) and to establish the ISS claim in Proposition 1. Furthermore, in order to ensure continuity of the ISS-Leonov function in (3.19) it is necessary to consider the angle θ evolving in R rather than, as usual, on the torus S. The same applies to the function W defined in (3.9). In R, θ is not bounded a priori. While it is shown in [9, Theorem 3] that existence of an ISS-Leonov function with θ ∈ R ensures boundedness of the distance |x| W with x ∈ S × R 3 , this does not imply boundedness of |x| W for x ∈ R 4 . Since boundedness of solutions is a fundamental prerequisite to invoke Lemma 1 (that ensures decomposability of W), we need the auxiliary Lemma 3 to establish boundedness of solutions of the system (3.5) evolving in R 4 .

Remark 4. By using the standard definitions for the reactance X, the conductance G and the susceptance B [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Kundur | Power System Stability and Control[END_REF], i.e.,

X = ω s L, G = R R 2 + X 2 , B = - X R 2 + X 2 , condition (3.7)-which is implied by Assumption 2-can be rewritten as 4DG(R 2 + X 2 ) 2 > X 2 b 2 ⇔ 4D R X > |B|b 2 .
This shows that a high damping factor D and a high R/X ratio, i.e., a high electrical dissipation, are beneficial to ensure ISS of the system (3.5) with respect to the equilibrium set W. On the contrary, a high value of |b|, i.e., a high excitation and consequently large EMF amplitude, deteriorate the likelihood of verifying the ISS property for the system (3.5). Both observations are sensible from a physical perspective and consistent with practical experience.

(d) Main stability result

We give now conditions for almost global asymptotic stability of one of the two equilibria (modulo 2π) of the system (2.11). This result follows as a corollary of Theorem 1 after the relaxation of Assumption 2-setting ψ = 0-as indicated below.

Corollary 1. Consider the system (3.5) with Assumption 1. Suppose that d(t) = 0 3 for all t ∈ R and that Assumption 2 is satisfied for ψ = 0. Then, the equilibrium point (θ s 1 , ω s , i s d , i s q ) satisfying |θ s 1 -φ| < π 2 (modulo 2π) with φ defined in (3.11) is almost globally asymptotically stable, i.e., for all initial conditions, except a set of measure zero, the solutions of the system (2.11) asymptotically converge to that equilibrium point.

Proof. For d(t) = 0 3 for all t ∈ R, we have that γ = 0 in Q in (3.20). Hence, Q ≥ 0 is equivalent to Qx ≥ 0, see (3.20). By following the proof of Theorem 1, with = 2 min , we see that Qx ≥ 0 is implied by Q > 0 for ψ = 0. As the latter holds by assumption, it follows that, see (3.25),

V ≤ -χV. (3.26) This shows that for all initial conditions, the solutions of the system (3.5) asymptotically converge to the set W. Recall from the proof of Theorem 1 that, under the standing assumptions, Ẇ ≤ 0. Furthermore, by the definition of q( θ) in (3.10), we have that q(0) = 0, q (0) = 0, q (0) = c cos(θ s 1 -φ).

Hence, the function W defined in (3.9) has a local minimum at the origin if |θ s 1 -φ| < π 2 , while it has a local maximum at |θ s 1 -φ| > π 2 . Therefore, if |θ s 1 -φ| < π 2 , the origin is locally asymptotically stable, while it is unstable if |θ s 1 -φ| > π 2 . To establish almost global asymptotic stability of the equilibrium point satisfying |θ s 1 -φ| < π 2 , we follow [START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF] and invoke [START_REF] Monzón | Local and global aspects of almost global stability[END_REF]Proposition 11]. To this end, consider the Jacobian of the system (3.5) evaluated at an equilibrium point, i.e.,

J =      0 1 0 0 b J i s q sin(θ s ) -i s d cos(θ s ) -D J -b J sin(θ s ) -D J cos(θ s ) b L ω s cos(θ s ) b L sin(θ s ) -R L -ω s -b L ω s sin(θ s ) b L cos(θ s ) ω s -R L     
.

By using (3.4) and (3.11), we obtain

det(J ) = b cos(θ s -φ) J R 2 + (Lω s ) 2 .
Hence, det(J ) < 0 if |θ s -φ| > π 2 , i.e., if the considered equilibrium point is unstable. Furthermore, J is a real-valued matrix of dimension 4. Consequently, if det(J ) < 0, then J has at least one positive real eigenvalue and from [START_REF] Monzón | Local and global aspects of almost global stability[END_REF]Proposition 11] it follows that the region of attraction of the unstable equilibrium has zero Lebesgue measure. Hence, for all initial conditions, except a set of measure zero, the solutions of the system (3.5) asymptotically converge to the equilibrium point (θ s 1 , ω s , i s d , i s q ) satisfying |θ s 1 -φ| < π 2 (modulo 2π). This completes the proof.

Numerical example

The analysis is illustrated on a numerical benchmark example with data taken from [1, ) is a linear matrix inequality (LMI) with decision variables χ > 0 and ψ > 0. Hence, its feasibility can be verified efficiently using standard software tools, such as Yalmip [START_REF] Löfberg | YALMIP : a toolbox for modeling and optimization in MATLAB[END_REF]. For the subsequent experiments we take D as additional decision variable and seek to determine the minimal damping coefficient D min required for our condition (3.16) to be satisfied. Thereby, following our previous work [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF][START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF], we consider a large variety of operating scenarios. More precisely, we modify the machine loading by varying Tm, as well as the excitation by varying i f . All other parameters are kept constant. Furthermore, we define the error coordinates x in (3.2) with respect to the unstable equilibrium point, i.e., satisfying |θ s 1 -φ| > π 2 . The resulting values for D min are shown in Figure 2a. Since the machine torque Te is normalized with respect to the single-phase power base, the rated torque corresponds to T s e = 3.0 pu. Also, usually it is only possible to load large SGs at their rated power if the amplitude of the EMF e abc is higher than the rated voltage [START_REF] Umans | Fitzgerald and Kingsley's Electric Machinery[END_REF]Chapter 5]. This explains the chosen ranges of values for E s and T s e in Figure 2a. By Theorem 1 the SMIB system (3.5) is ISS with respect to the corresponding equilibrium set W (if it exists) for any damping coefficient D ≥ D min . Furthermore, by Corollary 1 we conclude that-for d(t) = 0 3 for all t ≥ 0-the equilibrium point with |θ s 1 -φ| < It can be seen from the simulation results that the trajectories of the perturbed system (i.e., d(t) = 0 3 ) remain close to those of the unperturbed system (i.e., d(t) = 0 3 ), as predicted by the ISS property.

Moreover, in line with the observations made in [START_REF] Barabanov | Almost global attractivity of a synchronous generator connected to an infinite bus[END_REF][START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF], Figure 2a shows that the required damping increases with the machine loading and the excitation magnitude. This is consistent with the well-known power-angle characteristic that can be employed for reduced-order SG models and states that the likelihood of an SG to undergo instability after a change in load increases with the SG being operated closer to its generation limit [START_REF] Grainger | Power System Analysis[END_REF][START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF].

Conclusions

A global stability analysis of the equilibria of a realistic SMIB model has been presented. The main result is established by using the powerful property of ISS, which is verified using an ISS-Leonov function for periodic systems with multiple invariant sets-an approach first advocated in [START_REF] Efimov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF]. Compared to the usual Lyapunov functions, ISS-Leonov functions have the advantage that the the conditions of sign definiteness imposed are obviated. This key feature turns out to be 3 The scaling factor 3 in J and D originates from the following fact. In [START_REF] Anderson | Power System Control and Stability[END_REF], the mechanical equation (2.5) is expressed in pu with respect to the 3-phase base power S 3φ . Hence, the pu values of J, D, Tm and Te are also expressed with respect to S 3φ .

In the model (2.5), the electrical torque Te in (2.10) is expressed with respect to the single-phase power S 3φ 3 . Consequently, one way to match our model (2.5) with the parameters in [START_REF] Anderson | Power System Control and Stability[END_REF] is to scale J and D by a factor 3, i.e., to represent the mechanical equation (2.5) with respect to the power base essential for the analysis of the SMIB system carried out in this paper. The conservativeness of the derived conditions has been assessed via extensive numerical evaluations on a benchmark problem.

Our current research is directed towards the development of transient stability analysis methods for power system models with time-dependent transmission lines and physically more realistic SG models (e.g., considering torque and voltage control loops as well as saturation effects), which is a natural extension of the presented results. Given the "scalable" nature of the analysis tools employed here this seems a feasible-albeit difficult-task. Another objective is to further reduce the conservativeness of the conditions by improving the construction of the employed Lyapunov function.

  rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Schematic representation of a two-pole round-rotor SG based on [14, Figure 3.2]. The abc-axes correspond to the axes of the stator coils denoted by a 1,2 , b 1,2 and c 1,2 . The axis of coil a is chosen at θ = δ -ϕ = 0. The rotor current i f flows through the rotor windings. The dq-axes denote the rotating axes of the dq-frame corresponding to the mapping T dq (ϕ) with ϕ given in (2.7) and T dq (•) given in (5.1) of Appendix (a).

  rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The remainder of the claim is established by invoking [9, Theorem 3].

  rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

π 2 (

 2 modulo 2π) with φ defined in (3.11) is almost globally asymptotically stable. Both properties are further illustrated in Figure 2b for D = 62.4Dnom, T s e = 3pu and |E s | = 1.566pu. The bold lines correspond to the case d(t) = 0 3 , while the dashed lines correspond to the case d 1 = 0.25 sin(3t), d 2 = 0.5 sin(5t), d 3 = 0.3 sin(t).

S 3φ 3 .

 3 rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Minimum required damping coefficient D min compared to nominal damping coefficient Dnom for condition(3.16) to be satisfied in a broad range of operating conditions. Simulation results for the system (3.5).The solid lines show the state trajectories of ω, ĩd , ĩq for the case d(t) = 0 3 . The dashed lines correspond to the case d 1 = 0.25 sin(3t), d 2 = 0.5 sin(5t), d 3 = 0.3 sin(t).

Figure 2 :

 2 Figure 2: Numerical evaluation of the conditions (3.16) and illustration of the ISS property.

Table 1 :

 1 Table D.2, H16]. All system parameters are given in per unit (pu) including normalized time (for ω Base = 2π60rad/sec), see Table 1. Furthermore, based on [1, Example 5.1], we assume a line resistance of Parameters of the system (3.5) taken from [1, Table D.2, H16]

	Parameter Numerical value Parameter Numerical value
	ω s	in pu 1	J	in pu 3 • 2558.9
	Dnom	3 • 2	V	1
	R	0.0221	L	1.21
	R = 0.02 pu and a line inductance of L = 0.4 pu. The nominal damping coefficient is given by
	Dnom = 3 • 2 pu (see [1, Table D.2]) 3 .
	Condition (3.16

Here arctan(•) denotes the standard arctangent function defined on the interval [0, π). Since the case v d = vq = 0 is excluded by definition of v d and vq, φ is well-defined.

Appendix (a) The dq-transformation

The dq-transformation employed in the model derivation of the SMIB model in Section 2 is stated, see [START_REF] Anderson | Power System Control and Stability[END_REF][START_REF] Machowski | Power System Dynamics: Stability and Control[END_REF][START_REF] Schiffer | A survey on modeling of microgrids-from fundamental physics to phasors and voltage sources[END_REF]. Let x : R ≥0 → R 3 be a balanced three-phase signal, cf. (2.3), (2.1), and : R ≥0 → R.

Consider the mapping T

is called dq-transformation. Competing Interests. The author(s) declare that they have no competing interests.

Funding. This work was partially supported by the Government of Russian Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031).


Pierre Gloaguen

Marie-Pierre Etienne

Sylvain Le Corff

Online Sequential Monte Carlo smoother for partially observed stochastic differential equations

Keywords: Stochastic differential equations, Smoothing, Sequential Monte Carlo Methods

This paper introduces a new algorithm to approximate smoothed additive functionals for partially observed stochastic differential equations. This method relies on the recent procedure introduced in [24] which allows to compute such approximations online, i.e. as the observations are received, and with a computational complexity growing linearly with the number of Monte Carlo samples. The algorithm of [24] cannot be used in the case of partially observed stochastic differential equations since the transition density of the latent data is usually unknown. We prove that a similar algorithm may still be defined for partially observed continuous processes by replacing this unknown quantity by an unbiased estimator obtained for instance using general Poisson estimators. We prove that this estimator is consistent and its performance are illustrated using data from two models.

Introduction

This paper introduces a new algorithm to solve the smoothing problem for partially observed continuous time stochastic processes. In this setting, the hidden state process (X t ) t≥0 is assumed to be a solution to a stochastic differential equation (SDE) and the only information available is given by noisy observations (Y k ) 0≤k≤n of the states (X k ) 0≤k≤n at some discrete time points (t k ) 0≤k≤n . The bivariate stochastic process {(X k , Y k )} 0≤k≤n is a state space model such that conditional on the state sequence (X k ) 0≤k≤n the observations (Y k ) 0≤k≤n are independent and for all 0 ≤ ≤ n the conditional distribution of Y given {X k } 0≤k≤n depends on X only.

Statistical inference for partially observed state sequences often requires to solve bayesian filtering and smoothing problems, i.e. the computation of the posterior distributions of sequences of hidden states given observations. The filtering problem refers to the estimation, for each 0 ≤ k ≤ n, of the distributions of the hidden state X k given the observations (Y 0 , . . . , Y k ). Smoothing stands for the estimation of the distributions of the sequence of states (X k , . . . , X p ) given observations (Y 0 , . . . , Y ) with 0 ≤ k ≤ p ≤ ≤ n. These posterior distributions are crucial to compute maximum likelihood estimators of unknown parameters using the observations (Y 0 , . . . , Y n ) only. For instance, the E-step of the EM algorithm introduced in [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] boils down to the computation of a conditional expectation of an additive functional of the hidden states given all the observations up to time n. Similarly, by Fisher's identity, recursive maximum likelihood estimates may be computed using the gradient of the loglikelihood which can be written as the conditional expectation of an additive functional of the hidden states. See [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]Chapter 10 and 11], [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF][START_REF] Corff | Convergence of a particle-based approximation of the block online Expectation Maximization algorithm[END_REF][START_REF] Corff | Online Expectation Maximization based algorithms for inference in hidden Markov models[END_REF][START_REF] Poyiadjis | Particle approximations of the score and observed information matrix in state space models with application to parameter estimation[END_REF] for further references on the use of these smoothed expectations of additive functionals applied to maximum likelihood parameter inference in latent data models.

The exact computation of these expectations is usually not possible in the case of partially observed diffusions. In this paper, we propose to use Sequential Monte Carlo (SMC) methods to approximate smoothing distributions with random particles associated with importance weights. [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian bayesian state estimation[END_REF][START_REF] Kitagawa | Monte-Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF] introduced the first particle filters and smoothers for state space models by combining importance sampling steps to propagate particles with resampling steps to duplicate or discard particles according to their importance weights. Unfortunately, these methods cannot be applied directly to partially observed stochastic differential equations since some elementary quantities, such as transition densities of the hidden states, are not available explicitly. Discretization procedures may be used to approximate transition densities, for instance the Euler-Maruyama method, the Ozaki discretization which proposes a linear approximation of the drift coefficient between two observations [START_REF] Ozaki | A bridge between nonlinear time series models and nonlinear stochastic dynamical systems: a local linearization approach[END_REF][START_REF] Shoji | Estimation for nonlinear stochastic differential equations by a local linearization method 1[END_REF], or Gaussian based approximations using Taylor expansions of the posterior mean and variance of an observation given the observation at the previous time step, [START_REF] Kessler | Estimation of an ergodic diffusion from discrete observations[END_REF][START_REF] Kessler | Statistical methods for stochastic differential equations[END_REF][START_REF] Uchida | Adaptive estimation of an ergodic diffusion process based on sampled data[END_REF]. Other approaches based on Hermite polynomials expansion were also introduced by [START_REF] Ait-Sahalia | Transition densities for interest rate and other nonlinear diffusions[END_REF][START_REF] Ait-Sahalia | Maximum likelihood estimation of discretely sampled diffusions: a closed-form approximation approach[END_REF][START_REF] Ait-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF] and extended in several directions recently, see [START_REF] Li | Maximum-likelihood estimation for diffusion processes via closed-form density expansions[END_REF] and all the references on the approximation of transition densities therein. However, even the most recent discretization based approximations of the transition densities induce a systematic bias of particle based approximations of posterior distributions, see for instance [START_REF] Del Moral | The Monte Carlo method for filtering with discrete-time observations[END_REF]. To overcome this difficulty, [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF] proposed to solve the filtering problem by combining SMC methods with an unbiased estimate of the transition densities based on the generalized Poisson estimator (GPE). In this case, only the Monte Carlo error has to be controlled as there is no Taylor expansion to approximate unknown transition densities.

The only solution to solve the smoothing problem for partially observed SDE using SMC methods has been proposed in [START_REF] Olsson | Particle-based likelihood inference in partially observed diffusion processes using generalised Poisson estimators[END_REF] and extends the fixed-lag smoother of [START_REF] Olsson | Sequential monte carlo smoothing with application to parameter estimation in nonlinear state space models[END_REF]. Using forgetting properties of the hidden chain, the algorithm improves the performance of [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF] to approximate smoothing distributions but at the cost of a bias that does not vanish as the number of particles grows to infinity. In the case of discrete time state space models, approximations of the smoothing distributions may also be obtained using the Forward Filtering Backward Smoothing algorithm (FFBS) and the Forward Filtering Backward Simulation algorithm (FFBSi) developed respectively in [START_REF] Kitagawa | Monte-Carlo filter and smoother for non-Gaussian nonlinear state space models[END_REF][START_REF] Hürzeler | Monte Carlo approximations for general state-space models[END_REF][START_REF] Doucet | On sequential Monte-Carlo sampling methods for Bayesian filtering[END_REF] and [START_REF] Godsill | Monte Carlo smoothing for non-linear time series[END_REF]. Both algorithms require first a forward pass which produces a set of particles and weights approximating the sequence of filtering distributions up to time n. Then, a backward pass is performed to compute new weights (FFBS) or sample trajectories (FFBSi) in order to approximate the smoothing distributions. Recently, [START_REF] Olsson | Efficient particle-based online smoothing in general hidden Markov models: the PaRIS algorithm[END_REF] proposed a new SMC algorithm, the particle-based rapid incremental smoother (PaRIS), to approximate on-the-fly (i.e. using the observations as they are received) smoothed expectations of additive functionals. Unlike the FFBS algorithm, the complexity of this algorithm grows only linearly with the number of particles N and contrary to the FFBSi algorithm, no backward pass is required.

In this paper, we extend the use of PaRIS algorithm to partially observed SDE. The proposed algorithm allows to approximate smoothed expectations of additive functionals online and with a complexity growing only linearly with the number of particles. The crucial and simple result (Lemma 1) of the application of PaRIS algorithm to SDE is that the accept reject mechanism introduced in [START_REF] Douc | Sequential Monte Carlo smoothing for general state space hidden Markov models[END_REF] ensuring the linear complexity of the procedure is still correct when the transition densities are replaced by unbiased estimates. The usual FFBS and FFBSi algorithms may not be extended this easily since they both require the computation of weights defined as ratios involving the transition densities, thus replacing these unknown quantities by unbiased estimates does not lead to unbiased estimators of the weights. The proposed Generalized Random version of PaRIS algorithm, hereafter named GRand PaRIS algorithm, may be applied to general hidden Markov models whose Markovian dynamics is ruled by a stochastic differential equation (one of the first two domains defined in [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF]) but also to any general state space model where the transition density of the hidden chain may be estimated unbiasedly.

Section 2 describes the proposed algorithm to approximate smoothed additive functionals using unbiased estimates of the transition density of the hidden states and details the application of this algorithm when the transition density may be approximated using a GPE. In Section 3, classical convergence results for SMC smoothers are extended to the setting of this paper and illustrated with numerical experiments in Section 4. All proofs are postponed to Appendix A.

2 The Generalized Random PaRIS algorithm (X t ) t≥0 is defined as a weak solution to the following SDE in R d :

X 0 = x 0 and dX t = α(X t )dt + dW t , (1)
where (W t ) t≥0 is a standard Brownian motion. It is assumed that α is of the form α(x) = ∇ x A(x) where A : R d → R is a twice continuously differentiable function.

The solution to (1) is supposed to be partially observed at times t 0 = 0, . . . , t n through an observation process (Y k ) 0≤k≤n in (R m ) n+1 . For all 0 ≤ k ≤ n, the distribution of Y k given X k := X t k has a density with respect to a reference measure λ on R m given by g(X k , •) = g k (X k ). The distribution of X 0 has a density with respect to a reference measure µ on R d given by χ. For all 0 ≤ k ≤ n -1, the conditional distribution of X k+1 given X k has a density q k (X k , •) with respect to µ. Let 0 ≤ k ≤ k ≤ n, the joint smoothing distributions of the hidden states are defined, for all measurable function h on (R d ) k -k+1 , by:

φ k:k |n [h] = E [h(X k , . . . , X k )|Y 0:n ] .
For all 0 ≤ k ≤ n, φ k = φ k:k|k denote the filtering distributions. The aim of this section is to detail the extension of PaRIS algorithm to approximate expectations of the form

φ 0:n|n [H n ] = E [H n (X 0:n )|Y 0:n ] where H n = n-1 k=0 h k (X k , X k+1 ) , (2)
when the transition density of the hidden states is not available explicitly and where {h k } n-1 k=0 are given functions on R d × R d . The algorithm is based on the following link between the filtering and smoothing distributions for additive functionals, see [START_REF] Olsson | Efficient particle-based online smoothing in general hidden Markov models: the PaRIS algorithm[END_REF]:

φ 0:n|n [h] = φ n [T n [h]] , where T n [h](X n ) = E [h(X 0:n )|X n , Y 0:n ] . (3)
The approximation of (3) requires first to approximate the sequence of filtering distributions. Sequential Monte Carlo methods provide an efficient and simple solution to obtain these approximations using sets of particles

{ξ k } N =1 associated with weights {ω k } N =1 , 0 ≤ k ≤ n. At time k = 0, N particles {ξ 0 } N
=1 are sampled independently according to ξ 0 ∼ η 0 , where η 0 is a probability density with respect to µ. Then, ξ 0 is associated with the importance weights ω 0 = χ(ξ 0 )g 0 (ξ 0 )/η 0 (ξ 0 ). For any bounded and measurable function h defined on R d , the expectation φ 0 [h] is approximated by

φ N 0 [h] = 1 Ω N 0 N =1 ω 0 h ξ 0 , Ω N 0 := N =1 ω 0 . Then, for 1 ≤ k ≤ n, using {(ξ k-1 , ω k-1 )} N =1 , the auxiliary particle filter of [26] samples pairs {(I k , ξ k )} N
=1 of indices and particles using an instrumental transition density p k on R d × R d and an adjustment multiplier function ϑ k on R d . Each new particle ξ k and weight ω k at time k are computing following these steps:

-choose a particle index I k at time k -1 in {1, . . . , N } with probabilities proportional to ω j k-1 ϑ k (ξ j k-1 ), for j in {1, . . . , N } ;

-sample ξ k using this chosen particle according to ξ k ∼ p k (ξ

I k k-1 ,
•) ; -associate the particle ξ k with the importance weight:

ω k := q k (ξ I k k-1 , ξ k )g k (ξ k ) ϑ k (ξ I k k-1 )p k (ξ I k k-1 , ξ k ) . ( 4
)
The expectation φ k [h] is approximated by

φ N k [h] := 1 Ω N k N =1 ω k h ξ k , Ω N k := N =1 ω k .
PaRIS algorithm uses the same decomposition as the FFBS algorithm introduced in [START_REF] Doucet | On sequential monte-carlo sampling methods for bayesian filtering[END_REF] and the FFBSi algorithm proposed by [START_REF] Godsill | Monte Carlo smoothing for non-linear time series[END_REF] to approximate smoothing distributions. It combines both the forward only version of the FFBS algorithm with the sampling mechanism of the FFBSi algorithm. It does not produce an approximation of the smoothing distributions but of the smoothed expectation of a fixed additive functional and thus may be used to approximate [START_REF] Ait-Sahalia | Maximum likelihood estimation of discretely sampled diffusions: a closed-form approximation approach[END_REF]. Its crucial property is that it does not require a backward pass, the smoothed expectation is computed on-the-fly with the particle filter and no storage of the particles or weights is needed. PaRIS algorithm relies on the following fundamental property of

T k [H k ] when H k is as in (2): T k [H k ](X k ) = E [T k-1 [H k-1 ](X k-1 ) + h k-1 (X k-1 , X k )|X k , Y 0:k-1 ] , = φ k-1 (dx k-1 )q k-1 (x k-1 , X k ) {T k-1 [H k-1 ](x k-1 ) + h k-1 (x k-1 , X k )} φ k-1 (dx k-1 )q k-1 (x k-1 , X k ) .
Therefore, [START_REF] Olsson | Efficient particle-based online smoothing in general hidden Markov models: the PaRIS algorithm[END_REF] introduces sufficient statistics τ i k (starting with

τ i 0 = 0, 1 ≤ i ≤ N ), approximating T k [H k ](ξ i k ), for 1 ≤ i ≤ N and 0 ≤ k ≤ n. First, replacing φ k-1 by φ N k-1 in the last equation leads to the following approximation of T k [H k ](ξ i k ): T N k [H k ](ξ i k ) = N j=1 Λ N k-1 (i, j) T k-1 [H k-1 ](ξ j k-1 ) + h k-1 (ξ j k-1 , ξ i k ) , (5)
where

Λ N k (i, ) = ω k q k (ξ k , ξ i k+1 ) N =1 ω k q k (ξ k , ξ i k+1 ) , 1 ≤ ≤ N . (6)
Computing exactly these approximations would lead to a complexity growing quadratically with N because of the normalizing constant in [START_REF] Del Moral | The Monte Carlo method for filtering with discrete-time observations[END_REF]. Therefore, PaRIS algorithm samples particles in the set {ξ j k-1 } N j=1 with probabilities Λ N k (i, •) to approximate the expectation (5) and produce τ i k . Choosing Ñ ≥ 1, at each time step 0 ≤ k ≤ n -1 these statistics are updated according to the following steps.

(i) Run one step of a particle filter to produce

{(ξ k , ω k )} for 1 ≤ ≤ N . (ii) For all 1 ≤ i ≤ N , sample independently J i, k in {1, . . . , N } for 1 ≤ ≤ N with probabilities Λ N k (i, •),
given by ( 6).

(iii) Set

τ i k+1 := 1 N N =1 τ J i, k k + h k ξ J i, k k , ξ i k+1 .
Then, ( 2) is approximated by

φ N 0:n|n [τ n ] = 1 Ω N n N i=1 ω i n τ i n .
As proved in [START_REF] Olsson | Efficient particle-based online smoothing in general hidden Markov models: the PaRIS algorithm[END_REF], the algorithm is asymptotically consistent (as N goes to infinity) for any precision parameter Ñ . However, there is a significant qualitative difference between the cases Ñ = 1 and Ñ ≥ 2. As for the FFBSi algorithm, when there exists σ + such that 0 < q k < σ + , PaRIS algorithm may be implemented with O(N ) complexity using the accept-reject mechanism of [START_REF] Douc | Sequential Monte Carlo smoothing for general state space hidden Markov models[END_REF].

In general situations, PaRIS algorithm cannot be used for stochastic differential equations as q k is unknown. Therefore, the computation of the importance weights ω k and of the acceptance ratio of [START_REF] Douc | Sequential Monte Carlo smoothing for general state space hidden Markov models[END_REF] is not tractable. Following [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF][START_REF] Olsson | Particle-based likelihood inference in partially observed diffusion processes using generalised Poisson estimators[END_REF], filtering weights can be approximated by replacing q k (ξ k , ξ i k+1 ) by an unbiased estimator q k (ξ k , ξ i k+1 ; ζ k ), where ζ k is a random variable in R q such that:

q k (ξ k , ξ i k+1 ; ζ k ) > 0 a.s and E q k (ξ k , ξ i k+1 ; ζ k ) G N k+1 = q k (ξ k , ξ i k+1 ) ,
where, for all 0 ≤ k ≤ n,

F N k = σ Y 0:k ; (ξ u , ω u , τ u ); J ,j v ; 1 ≤ ≤ N, 0 ≤ u ≤ k, 1 ≤ j ≤ N , 0 ≤ v < k , G N k+1 = F N k ∨ σ Y k+1 ; (ξ k+1 , ω k+1 ); 1 ≤ ≤ N .
Practical choices for ζ k are discussed below, see for instance [START_REF] Doucet | On sequential Monte-Carlo sampling methods for Bayesian filtering[END_REF] which presents the choice made for the implementation of such estimators in our context. In the case where q k is unknown, the filtering weights in (4) then become:

ω k := q k (ξ I k k-1 , ξ k ; ζ k )g k (ξ k ) ϑ k (ξ I k k-1 )p k (ξ I k k-1 , ξ k ) . (7)
Therefore, to obtain a generalized random version of PaRIS algorithm, we only need to be able to sample from the discrete probability distribution Λ N k (i, •) in the case of SDE based HMM. Consider the following assumption: for all 0 ≤ k ≤ n, there exists a random variable σk + measurable with respect to G N k+1 such that,

sup x,y,ζ q k (x, y; ζ) ≤ σk + . (A1)
Lemma 1. Assume that A 1 holds for some 0 ≤ k ≤ n -1. For all 1 ≤ i ≤ N , define the random variable J i k as follows:

repeat Sample independently ζ, U ∼ U[0, 1] and J ∈ {1, . . . , N } with probabilities proportional to { ω 1 k , . . . , ω N k }. until U ≤ q k (ξ J k , ξ i k+1 , ζ)/σ k + . Set J i k = J. Then, the conditional probability distribution given G N k+1 of J i k is Λ N k (i, •). Proof. See Appendix A.
Note that Lemma 1 still holds if assumption (A1) is relaxed and replaced by one of the two following assumptions:

sup j,y,ζ q k (ξ j k , y, ζ) ≤ σk + . (A2) sup i,j,ζ q k (ξ j k , ξ j k+1 , ζ) ≤ σk + . (A3)
It is worth noting that under assumptions (A1) or (A2), the linear complexity property of PaRIS algorithm still holds, whereas if only assumption (A3) holds, the algorithm has a quadratic complexity.

Bounded estimator of q k

For x, y ∈ R d , by Girsanov and Ito's formulas, the transition density q k (x, y) of (1) satisfies, with ∆ k = t k+1 -t k ,

q k (x, y) = ϕ ∆ k (x, y) exp {A(y) -A(x)} E W x,y,∆ k exp - ∆ k 0 φ(w s )ds ,
where W x,y,∆ k is the law of Brownian bridge starting at x at 0 and hitting y at ∆ k , (w t ) 0≤t≤∆ k is such a Brownian bridge, ϕ ∆ k (x, y) is the p.d.f. of a normal distribution with mean x and variance ∆ k , evaluated at y and φ : R d → R is defined as:

φ(x) = α(x) 2 + A(x) /2 ,
with the Laplace operator. Assume that there exist random variables L w and U w such that for all 0 ≤ s ≤ ∆ k , L w ≤ φ(w s ) ≤ U w . Let κ be a random variable taking values in N with distribution µ and (U j ) 1≤j≤κ be independent uniform random variables on [0, ∆ k ], and ζ k = {κ, w, U 1 , . . . , U κ } . As shown in [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF], a positive unbiased estimator is given by

q k (x, y; ζ k ) = ϕ ∆ k (x, y) exp {A(y) -A(x)} × exp {-U w ∆} ∆ κ k µ(κ)κ! κ j=1 U w -φ(w U j ) . ( 8
)
Interesting choices of µ are discussed in [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF] and we focus here on the so called GPE-1, where µ is a Poisson distribution with intensity (U w -L w )∆ k . In that case, the estimator (8) becomes:

q k (x, y; ζ k ) = ϕ ∆ k (x, y) exp {A(y) -A(x) -L w ∆ k } κ j=1 U w -φ(w U j ) U w -L w . (9)
On the r.h.s. of (9), the product over κ elements is bounded by 1, therefore, a sufficient condition to satisfy of the assumptions (A1)-(A3) is that the function:

ρ ∆ k : R d × R d → R (x, y) → ϕ ∆ k (x, y) exp {A(y) -A(x) -L w ∆ k } (10)
is upper bounded almost surely by σk + . In particular, if L w is bounded almost surely, (10) always satisfies assumption (A3) and Algorithm 1 can be used. This condition is always satisfied for models in the domains D 1 and D 2 defined in [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF], i.e. domains for which the exact algorithms EA1 and EA2 can be used.

When (A1) or (A2) holds, it can be nonetheless of practical interest to choose the bound σk + corresponding to (A3). Indeed, this might increase significantly the acceptance rate of the algorithm, and therefore reduce the number of drawings of the random variable ζ, which has a much higher cost than the computation of ρ, as it requires simulations of Brownian Bridges. Moreover, this latter option can also avoid numerical optimization if no analytical expression of σk + is available. In practice, we found this option more efficient in terms of computation time when N has moderate values.

Algorithm 1 GRand PaRIS algorithm

for all i ∈ 1, . . . , N do Sample ξ i 0 ∼ η 0 , τ i 0 = 0 and ω i 0 = g 0 (ξ i 0 )χ 0 (ξ i 0 )/η 0 (ξ i 0 ). end for for k ∈ 0, . . . , n -1 do for all i ∈ 1, . . . , N do Set τ i k+1 = 0; Sample I i k+1 in {1, . . . , N } with probabilities proportional to

{ ω 1 k ϑ k+1 (ξ 1 k ), . . . , ω N k ϑ k+1 (ξ N k )}. Sample ξ i k+1 ∼ p k (ξ I i k+1 k , •). For all 1 ≤ m ≤ M , sample independently ζ m k = (κ m , w m , (U m j ) 1≤j≤κm ) with κ m ∼ µ, w m ∼ W X k+1 k and (U m j ) 1≤j≤κm ∼ U[0, ∆ k ] ⊗κm . Compute ω i
k+1 using equation [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. for all ∈ 1, . . . , N do Sample J i, k as in Lemma 1. Update τ i k+1 = τ i k+1 + (τ

J i, k k + h k (ξ J i, k k , ξ i k+1
))/ Ñ . end for end for end for [START_REF] Ait-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF]

Convergence results

Consider the following assumptions.

H1 (i) For all k ≥ 0 and all x ∈ R d , g k (x) > 0. (ii) sup k≥0 |g k | ∞ < ∞. H2 sup k≥1 |ϑ k | ∞ < ∞, sup k≥1 |p k | ∞ < ∞ and sup k≥1 | ω k | ∞ < ∞, where ω 0 (x) = χ(x)g 0 (x) η 0 (x) and for k ≥ 1 ω k (x, x ; z) = q k (x, x ; z)g k+1 (x ) ϑ k+1 (x)p k (x, x ) . Lemma 2. For all 0 ≤ k ≤ n -1, the random variables { ω i k+1 τ i k+1 } N i=1 are independent conditionally on F N k and E ω 1 k+1 τ 1 k+1 F N k = φ N k [ϑ k+1 ] -1 φ N k q k (•, x)g k+1 (x) {τ k (•) + h k+1 (•, x)} dx .
Proof. See appendix A Proposition 1. Assume that H1 and H2 hold and that for all 1 ≤ k ≤ n, osc(h k ) < +∞. For all 0 ≤ k ≤ n and all N ≥ 1, there exist b k , c k > 0 such that for all N ≥ 1 and all ε ∈ R + ,

P φ N k [τ k ] -φ k [T k h k ] ≥ ε ≤ b k exp -c k N ε 2 .
Proof. See appendix A

Numerical experiments

This section investigates the performance of the proposed algorithm with the sine and log-growth models. In both cases, the proposal distribution p k is chosen as the following approximation of the optimal filter (or the fully adapted particle filter in the terminology of [START_REF] Pitt | Filtering via simulation: Auxiliary particle filters[END_REF]):

p k (x k-1 , x k ) ∝ qk (x k-1 , x k )g k (x k ) ,
where qk (x k-1 , x k ) is the p.d.f. of Gaussian distibution with mean α(x k-1 )∆ k and variance ∆ k I d , i.e. the Euler approximation of equation [START_REF] Ait-Sahalia | Transition densities for interest rate and other nonlinear diffusions[END_REF]. As the observation model is linear and Gaussian, the proposal distribution is therefore Gaussian with explicit mean and variance. In order to evaluate the performance of the proposed algorithm, the following strategy has been chosen. We compare the estimation of the EM intermediate quantity with the one obtained by the fixed lag method of [START_REF] Olsson | Particle-based likelihood inference in partially observed diffusion processes using generalised Poisson estimators[END_REF], for different values of the lag (namely, 1,2,5,10,50). The particle approximation of Q(θ, θ) for each model is computed using each algorithm, see Figure 1 for the SINE model (and respectively Figure 3 for the log-growth model). This estimation is performed 200 times to obtain the estimates Q 1 , . . . , Q 200 , using Ñ = 2 particles for PaRIS algorithm, and M = 30 replications for the Monte Carlo approximation q k of each q k . Moreover, the E step requires the computation of a quantity such as (2) with h k = log g k +log q k . log q k is not available explicitly and is approximated using the unbiased estimator proposed in [23, Appendix B] based on 30 independent Monte Carlo simulations. The intermediate quantity of the EM algorithm is also estimated with our algorithm 30 times using N = 5000 particles, the reference value is then computed as the arithmetic mean of these 30 estimations, and denoted by Q . Figure 1 (resp. 3 ) shows this estimate for an example on one simulated data set. The GRand Paris algorithm is performed using N = 400 particles in both cases, the fixed lag technique using N = 1600 so that both estimations require similar computational times.

The SINE model

The performance of the GRand PaRIS algorithm are first highlighted using the SINE model, where (X t ) t≥0 is supposed to be the solution to:

dX t = sin (X t -θ) dt + dW t , X 0 = x 0 . (11)
This simple model has no explicit transition density, however GPE estimators may be computed by simulating Brownian bridges. The process solution to [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF] is observed regularly at times t 0 = 0, . . . , t 100 = 50 through the observation process (Y k ) 0≤k≤100 :

Y k = X k + ε k , (12)
where the (ε) 0≤k≤100 are i.i.d. N (0, 1). In the example displayed on Figure 1, we set θ = 0. In that case, the function ρ ∆ k defined in [START_REF] Doucet | On sequential monte-carlo sampling methods for bayesian filtering[END_REF] can be upper bounded either on (x, y) or only on y, the GRand PaRIS algorithm has therefore a linear complexity. This same experiment was reproduced on 100 different simulated data sets. For each simulation s, the empirical absolute relative bias arb s and the empirical absolute coefficient of variation acv s are computed as

arb s = |m( Q s ) -Q s | | Q s | (13)
acv s = σ( Q s ) |m( Q s )| (14)
where m( Q s ) and σ( Q s ) are the empirical mean and standard deviation of the sample Q s 1 , . . . , Q s 200 . For each estimation method, the resulting distributions of arb 1 , . . . , arb 100 and acv 1 , . . . , acv 100 are shown on Figure 2.

The GRand PaRIS algorithm outperforms the fixed lag methods for any value of the lag as the bias is the lowest (it is already negligible for N = 400) and with a lower variance than fixed lag estimates with negligible bias (i.e., in this case, lags larger than 10). Small lags lead to strongly biased estimates for the fixed lag method, and unbiased estimates are at the cost of a large variance. It is worth noting here that the lag for which the bias is small is model dependent.

Log-growth model

Following [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF] and [START_REF] Olsson | Efficient particle-based online smoothing in general hidden Markov models: the PaRIS algorithm[END_REF], the performance of the proposed algorithm are also illustrated with the log-growth model defined by: Particle smoother for SDE

dZ t = κZ t 1 - Z t γ dt + σZ t dW t , Z 0 = z 0 . (15)
In order to use the exact algorithms of [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF] and the GPE of [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF], we consider (15) after the Lamperti transform, i.e., the process defined by X t = η(Z t ), with η(z) := -log(z)/σ, which satisfies the following SDE:

dX t = :=α(Xt) σ 2 - κ σ + κ γσ exp (-σX t ) dt + dW t , X 0 = x 0 = η(z 0 ). ( 16
)
In this case, the conditions of the Exact Algorithm 2 defined in [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF] are satisfied, as for any m ∈ R there exists U m such that for all x ≥ m, ψ(x) := α 2 (x)+α (x) ≤ U m . Moreover, ψ is lower bounded uniformly by L. Then, GPE estimators may be computed by simulating the minimum of a Brownian bridge, and simulating Bessel bridges conditionally to this minimum, as proposed by [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF].

The process solution to ( 16) is observed regularly at times t 0 = 0, . . . , t 50 = 100 through the observation process (Y k ) 0≤k≤50 defined as:

Y k = X k + ε k , (17)
where the (ε k ) 0≤k≤50 are i.i.d. N (0, σ 2 obs ). The parameters are given by θ = (κ = 0.1, σ = 0.1, γ = 1000, σ 2 obs = 4) .

In that case, the ρ ∆ k function defined in [START_REF] Doucet | On sequential monte-carlo sampling methods for bayesian filtering[END_REF] can be upper bounded as a function of y when x ∈ {ξ 1 k , . . . , ξ N k }, the GRand PaRIS algorithm has therefore a linear complexity. The intermediate quantity of the EM algorithm is evaluated as for the SINE model, see Figures 3 and4.

The results for the fixed lag technique are similar to the ones presented in [23, Figure 1] on the same model. For small lags, the variance of the estimates is small, but the estimation is highly biased. The bias rapidly decreases as the lag increases, together with a great increase of variance. Again, the GRand PaRIS algorithm outperforms the fixed lag smoother as it shows a similar (vanishing) bias as the fixed lag for the largest lag and a smaller variance than the fixed lags estimates with negligible bias.

Conclusions

This paper presents a new online SMC smoother for partially observed differential equations. This algorithm relies on an accept-reject procedure inspired from the recent PaRIS algorithm. The main result of the article for practical applications is that the mechanism of this procedure remains valid when the transition density is approximated by a an unbiased positive estimator. The proposed procedure outperforms the existing fixed lag smoother for SDE of [START_REF] Olsson | Particle-based likelihood inference in partially observed diffusion processes using generalised Poisson estimators[END_REF], as it does not introduce an intrinsic and non vanishing bias. In addition, numerical simulations highlight a better variance using data from two different models. It can be implemented for the class of models D 1 and D 2 defined in [START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discusion)[END_REF] with a linear complexity in N .

A Proofs

Proof of Lemma 1. Let τ be the first time draws are accepted in the accept-reject mechanism. For all ≥ 1, write

A k = U < q k (ξ J k , ξ i k+1 , ζ k )/σ k + .
Let h be a function defined on {1, . . . , N },

E h(J i,j k ) G N k+1 = m≥1 E h(J m )1 τ =m G N k+1 , = m≥1 h(m) m-1 =1 E 1 (A k ) c G N k+1 E h(J m )1 A k m G N k+1 , = m≥1 m-1 =1 E 1 - q k (ξ J k , ξ i k+1 ; ζ k ) σk + G N k+1 × E h(J m ) q k (ξ Jm k , ξ i k+1 ; ζ m k ) σk + G N k+1 , = m≥1 E 1 - q k (ξ J 1 k , ξ 1 k+1 ) σk + G N k+1 m-1 E h(J 1 ) q k (ξ J 1 k , ξ 1 k+1 ) σk + G N k+1 , = E h(J 1 )q k (ξ J 1 k , ξ i k+1 ) G N k+1 /E q k (ξ J 1 k , ξ i k+1 ) G N k+1 , = N =1 h( )ω k-1 q k (ξ k , ξ i k+1 ) N m=1 ω m k-1 q k (ξ m k , ξ i k+1 ) , = N =1 Λ N k-1 (i, )h( ) ,
which concludes the proof.

Proof of Lemma 2. The independence is ensured by the mechanism of SMC methods. By [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF],

E ω i k+1 τ i k+1 F N k = E   q k (ξ I i k+1 k , ξ i k+1 ; ζ k )g k+1 (ξ i k+1 ) ϑ k+1 (ξ I i k+1 k )p k+1 (ξ I i k+1 k , ξ i k+1 ) τ i k+1 F N k   .
Particle smoother for SDE Note that by Lemma 1,

E τ i k+1 G N k+1 = N =1 ω k q k (ξ k , ξ i k+1 ) τ k + h k (ξ k , ξ i k+1 ) N =1 ω k q k (ξ k , ξ i k+1 ) , E q k (ξ I i k+1 k , ξ i k+1 ; ζ k ) G N k+1 = q k (ξ I i k+1 k , ξ i k+1 ) .
Since τ i k+1 and ζ k are independent conditionally to G N k+1 :

E τ i k+1 q k (ξ I i k+1 k , ξ i k+1 ; ζ k ) G N k+1 = q k (ξ I i k+1 k , ξ i k+1 ) N =1 ω k q k (ξ k , ξ i k+1 ) τ k + h k (ξ k , ξ i k+1 ) N =1 ω k q k (ξ k , ξ i k+1 )
.

Moreover, conditionally to F N k , the probability density function of (ξ i k+1 , I i k+1 ) is given by (x, j) → ω j k ϑ k+1 (ξ j k )p k (ξ j k , x)

Ω k φ N k [ϑ k+1 ]
.

Therefore, this yields:

E ω i k+1 τ i k+1 F N k = φ N k [ϑ k+1 ] -1 N j=1 ω j k Ω k ϑ k+1 (ξ j k )
q k (ξ j k , x)g k+1 (x) ϑ k+1 (ξ j k )p k (ξ j k , x)

× N =1 ω k q k (ξ k , x) τ k + h k (ξ k , x) N =1 ω k q k (ξ k , x) p k (ξ j k , x)dx , = φ N k [ϑ k+1 ] -1 × N =1 ω k Ω k N j=1 ω j k q k (ξ j k , x) N =1 ω k q k (ξ k , x) g k+1 (x)q k (ξ k , x) τ k + h k (ξ k , x) dx = φ N k [ϑ k+1 ] -1 φ N k q k (•, x)g k+1 (x) {τ k (•) + h k (•, x)} dx ,
which concludes the proof.

Proof of Proposition 1. The results is proved by induction. At time k = 0, the result holds using that for all 1 ≤ i ≤ N , ρ i 0 = 0 and the convention T 0 [h 0 ] = 0. In addition, φ N 0 is a standard importance sampler estimator of φ 0 with ω i 0 ≤ | ω 0 | ∞ so that for any bounded function h on X,

P φ N 0 [h] -φ 0 [h] ≥ ε ≤ b 0 exp -c 0 N ε 2 .
Assume the results holds for k ≥ 1 and that ϑ k+1 = 1 for simplicity. Write

φ N k+1 [τ k+1 ] -φ k+1 [T k+1 [h k+1 ]] = a N /b N ,
where a N = N -1 N i=1 ω i k+1 τ i k+1 -φ k+1 [T k+1 [h k+1 ]] and b N = N -1 N i=1 ω i k+1 . By Lemma 2, the random variables { ω i k+1 τ i k+1 } N i=1 are independent conditionally on F N k and by H2,

ω i k+1 τ i k+1 -φ k+1 [T k+1 [h k+1 ]] ≤ 2| ω k+1 | ∞ |H k+1 | ∞ .
Therefore, by Hoeffding inequality,

P a N -E a N F N k ≥ ε = E P a N -E a N F N k ≥ ε F N k ≤ 2 exp -c k N ε 2 .
On the other hand,

E a N F N k = φ N k [Υ k ] , where Υ k (x k ) = q k (•, x)g k+1 (x) (τ k (x k ) + h k+1 (x k , x) -φ k+1 [T k+1 [h k+1 ]]) dx .
P b N -E b N F N k ≥ ε = E P b N -E b N F N k ≥ ε F N k ≤ 2 exp -c k N ε 2 .
Note that

E b N F N k = φ N k q k (•, x)g k+1 (x)dx .
By the induction assumption,

P E b N F N k -φ k q k (•, x)g k+1 (x)dx ≥ ε ≤ b k exp -c k N ε 2 .
The proof is completed using Lemma 3. Then,

P a N b N > ≤ B exp -CN β 2M 2 .
Proof. See [START_REF] Douc | Sequential Monte Carlo smoothing for general state space hidden Markov models[END_REF].

Figure 1 :

 1 Figure 1: SINE model. Process X solution to the SDE (balls) and observations Y (circles) at times t 0 = 0, . . . , t 100 = 50 [left]. Estimation of the EM intermediate quantity Q(θ, θ) using the fixed lag (FL) technique for 5 different lags, and the GRand PaRIS algorithm using 200 replicates [right]. The whiskers represent the extent of the 95% central values. The dot represents the empirical mean over the 200 replicates. The dotted line shows the reference value, computed using the GRand PaRIS algorithm with N = 5000 particles.

Figure 2 :

 2 Figure 2: SINE model. Distribution of the empirical absolute relative bias [left] and of the empirical absolute coefficient of variation [right] for each method.

Figure 3 :

 3 Figure 3: Log-growth model. Process X solution to the SDE (balls) and observations Y (circles) at times t 0 = 0, . . . , t 100 = 50 [left]. Estimation of the EM intermediate quantity Q(θ, θ) using the fixed lag (FL) technique for 5 different lags, and the GRand PaRIS algorithm using 200 replicates [right]. The whiskers represent the extent of the 95% central values. The dot represents the empirical mean over the 200 replicates. The dotted line shows the reference value, computed using the GRand PaRIS algorithm with N = 5000 particles.

Figure 4 :

 4 Figure 4: Log-growth model. Distribution of the empirical absolute relative bias [left] and of the empirical absolute coefficient of variation [right] for each method.

By [ 24 ,

 24 Lemma 11], φ k [Υ k ] = 0 which implies by the induction assumption thatP E a N F N k ≥ ε ≤ b k exp -c k N ε 2 .Then,P (|a N | ≥ ε) ≤ b k exp -c k N ε 2 . Similarly, as b N ≤ | ω k | ∞ , by Hoeffding inequality,

Lemma 3 .

 3 Assume that a N , b N , and b are random variables defined on the same probability space such that there exist positive constants β, B, C, and M satisfying Particle smoother for SDE (i) |a N /b N | ≤ M , P-a.s. and b ≥ β, P-a.s., (ii) For all > 0 and all N ≥ 1, P [|b N -b| > ] ≤ B exp (-CN 2 ), (iii) For all > 0 and all N ≥ 1, P [|a N | > ] ≤ B exp -CN ( /M ) 2 .

AgroParistech, UMR MIA 518, F-75231 Paris, France.

Laboratoire de Mathématiques d'Orsay, Univ. Paris-Sud, CNRS, Université Paris-Saclay.