Helmer Keator B David

G Karl

Satrajit S Ghosh

Auer Tibor

Camille Maumet

Das Samir

Flandin Guillaume

Nichols Thomas

J Krzysztof

Jessica Gorgolewski

Turner

Turner J 10

Kennedy D 11

Poline J B 12

Describing Assessments and Experiment Metadata with the Neuroimaging Data Model (NIDM)

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A fundamental goal of organizing and annotating scientific data is to provide investigators with an effective mechanism to access and share data with the wider scientific community. To attain these goals, the data needs to have both adequate descriptions to make it useable [START_REF] Poline | Data sharing in neuroimaging research[END_REF] and be represented in a structured form, accessible to computational tools. The Neuroimaging Data Model (NIDM; http://nidm.nidash.org/ ) is an ongoing effort to represent, in a single data model, the different components of a research activity (e.g., participants, project information, derived data), their relations, and provenance [START_REF] Keator | Towards structured sharing of raw and derived neuroimaging data across existing resources[END_REF]. NIDM has a modular design currently consisting of NIDM-Results, NIDM-Workflows, and NIDM-Experiment. NIDM-Workflows is focused on the description of data analysis pipelines and detailed software-specific variations. NIDM-Results is focused on the representation of mass-univariate neuroimaging results using a common descriptive standard across neuroimaging software packages. NIDM-Experiment is focused on the representation of the experiment design, the source data collected during the experiment and information on the participants, including generic assessments, demographics, and visit information. This abstract presents our prototyping work on representing assessments and metadata collected during the course of a typical neuroimaging experiment.

Methods

NIDM is developed using a community-driven process that engages stakeholders to participate in the identification of use-cases that drive development. In-person workshops and weekly video conferences are used to maintain communication while example NIDM documents, specifications, and software are developed by the INCF-NIDASH task force members. NIDM is an extension of the W3C recommended Provenance Data Model (PROV-DM ; www.w3.org/TR/prov-dm/ ) for derived data and experiment descriptions in neuroimaging related fields. NIDM Experiment captures details about an investigation (Figure 1). This pattern models metadata from neuroimaging data management systems at three levels: Investigation, Session, and Series. The Investigation level captures administrative information, while the Session and Series levels model assessments and MRI sessions. The NIDM-Assessment object model (Figure 2) provides a generic model for describing the data dictionary and acquired assessment data. Using the NIDM-Assessment model, an assessment "DataStructure" is composed of individual "DataElements" corresponding to the questions in the assessment. When an assessment's question contains multiple response choices, it is coded as a "ValueSet" in NIDM-Assessment and documents the answer choices and their coded values (Figure 3). Data collected for a particular assessment are stored as acquisition objects of the assessment type described by the corresponding "DataStructure". These objects are associated with an acquisition activity found at the Session or Series levels of the NIDM-Experiment model, depending on when the assessment was administered. The NIDM-Assessment representations are serialized using the Terse RDF Triple Language ( http://www.w3.org/TeamSubmission/turtle/ ) and queried using SPAQRL ( https://www.w3.org/TR/rdf-sparql-query/ ).

Results

The NIDM-Experiment and NIDM-Assessment models are actively under development and currently being evaluated by mapping multiple datasets (e.g. FBIRN [START_REF] Keator | The Function Biomedical Informatics Research Network Data Repository[END_REF], OpenfMRI [START_REF] Poldrack | Toward open sharing of task-based fMRI data: the OpenfMRI project[END_REF], NCANDA [START_REF] Rohlfing | N-CANDA data integration: anatomy of an asynchronous infrastructure for multi-site, multi-instrument longitudinal data capture[END_REF], etc.) into the structures. Our examples are made available through the NIDM GitHub respository: https://github.com/incf-nidash/nidm . The Conte Center on Brain Programming in Adolescent Vulnerabilities at the University of California, Irvine ( http://contecenter.uci.edu ) has made use of the early NIDM-Experiment related models in building its production informatics resource, which provides center investigators with inter-project data, including models for maternal and fetal heart rate, fetal movement, maternal blood oxygenations levels, derived DTI and fMRI measurements, and brain connectivity graphs. Our initial experiments have shown the model to be flexible and applicable to a wide range of heterogenous data types.

Conclusions

NIDM Experiment supports a variety of use cases focused on representing primary and derived data organization with the intent of simplifying data exchange, integration, and sharing. Our preliminary results indicate a modeling approach based on the W3C PROV-DM and NIDM is appropriate for modeling the heterogeneous data often collected in the course of a neuroimaging-related investigation and demonstrates the benefit of using semantic-web methods for data annotation. Continued work will include finalizing the NIDM-Experiment model and creating applications to aid investigators in converting their existing data to the NIDM-Experiment representation, and using the representations in analysis tools and information management applications.

Figure 1 :

 1 Figure 1: The NIDM-Experiment core structure. The Investigation level describes the overall project including investigators and project personnel, a description of the project, consent forms, etc. The Session level describes the data acquisition activities for a subject visits. The Series level describes the imaging series and behavorial/clinical assessments administered during a session.

Figure 2 :

 2 Figure 2: The NIDM-Assessment object model. An assessment consists of a DataStructure entity describing the overall assessment and one or more DataElement and/or CodedDataElements, describing the assessment questions. CodedDataElements are described using ValueSets.

Figure 3 :

 3 Figure 3: Partial NIDM representation of demographics assessment data dictionary.

Acknowledgments

We acknowledge the work of all INCF task force members and many other colleagues who have helped in this effort. Further, we acknowledge support by the UCI Conte Center (1P50MH096889-01A1), the NCANDA Consortium (U01AA021697-01), and the ReproNim Center for Reproducible Neuroimaging Computation (1P41EB019936-01A1).


Lotfi Braham

Ayoub Mediouni

Marius Nouri

Mahieddine Bozga

Axel Dellabani

Saddek Legay

Bensalem

Bozga

Braham Lotfi Mediouni

Ayoub Nouri

Marius Bozga

Mahieddine Dellabani

Axel Legay

Saddek Bensalem

Statistical Model Checking Stochastic Real-time Systems

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Statistical Model Checking (SMC) is a powerful alternative to classical numerical probabilistic model-checking that generally fail to handle large state-space systems. SMC was successfully applied in the assessment of different real-life systems in various application domains. Classical model checkers [START_REF] Kwiatkowska | Prism 4.0: verification of probabilistic real-time systems[END_REF][START_REF] David | Uppaal smc tutorial[END_REF] now include SMC as part of their analysis engines and have been recently joined by a variety of specialized ones [START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF][START_REF] Jegourel | A platform for high performance statistical model checking -plasma[END_REF][START_REF] Ballarini | Hasl: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF][START_REF] Nouri | Statistical model checking QoS properties of systems with SBIP[END_REF]. All these tools mainly differ in their modeling and properties specification formalisms. Uppaal-smc [START_REF] David | Uppaal smc tutorial[END_REF] considers Networks of Priced Timed Automata, which are high-level representations of D/CTMCs for system modeling, and weighted MTL for properties specification. Prism [START_REF] Kwiatkowska | Prism 4.0: verification of probabilistic real-time systems[END_REF] treats in addition Markov Decision Processes and Probabilistic Timed Automata for modeling, and Probabilistic Computation Tree, Continuous Stochastic Logic (CSL), and LTL for specification. Plasma Lab [START_REF] Jegourel | A platform for high performance statistical model checking -plasma[END_REF] is a modular statistical model checker that allows to use external simulators and checkers. Its default configuration supports DTMCs specified in a Prism dialect and bounded LTL. Ymer [START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF] is one of the rare tools to implement SMC (Hypothesis testing) for GSMPs and CSL, however it is no more maintained. Finally, COSMOS [START_REF] Ballarini | Hasl: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF] relies on Generalized Stochastic Petri Nets as input models and Hybrid Automata Stochastic Logic, a more expressive formalism, for properties specifications.

In this paper, we present the newest release of SBIP, a statistical model checker that enriches the existing BIP tool-set [START_REF] Basu | Rigorous component-based system design using the bip framework[END_REF] with statistical analyses. BIP provides a general framework to support design activities ranging from specification and validation to implementation and deployment in a rigorous way. To implement this vision, a rich tool-set was built for modeling, languages embedding, functional validation, models transformation and distributed code generation.

In its previous version [START_REF] Nouri | Statistical model checking QoS properties of systems with SBIP[END_REF], SBIP was limited to the analysis of DTMCs with respect to bounded LTL properties. In this release, it was redesigned and extended to support GSMPs, CTMCs, MTL, parametric exploration of LTL and MTL properties and analysis of rare events. The tool has also benefited from a major revision of its workflows and GUI. It now provides an Integrated Development Environment (IDE) where one can edit, compile, simulate models, and perform analyses. Additionally, SBIP is now organized around wellstructured projects that enclose models, properties and traces. It also includes support for graphical visualization of analysis results.

SBIP Design and Functionalities

SBIP is fully developed in Java and runs on GNU/Linux. It is freely available at http://www-verimag.imag.fr/Statistical-Model-Checking.html. The tool is distributed with a large set of case studies and a detailed documentation (e.g., user manual, installation details, video tutorials). For the sake of simplicity, we also provide a virtual machine with a pre-installed version of the tool.

This new release was designed in a modular fashion to allow more flexibility and extensibility. As depicted in Fig. 1, SBIP consists of three generic functional modules: Stochastic Simulation Engine, Monitoring, and Statistical Analyses that currently include Hypothesis Testing (HT), Probability Estimation (PE), Parametric Exploration (PX) and Importance Splitting (IP) for rare events analysis. All these modules are fully independent and interact through well-defined Java interfaces. The latter also define a clean and easy way to extend the tool with further modules (simulators, monitors and analyzers). In practice, statistical analysis algorithms trigger the stochastic simulation engine to produce a new execution trace which is monitored against an input property to produce a local verdict. Depending on the used analysis method, several iterations are generally required, to produce the final verdict. The proposed design allows to perform different analyses in separate workflows, namely simple simulation, standard SMC analyses, parametric SMC exploration and analysis of rare events. These workflows rely on common features such as models and properties edition, compilation and generated traces inspection.

Stochastic Simulation Engine. Currently, SBIP allows to use two different stochastic simulators, namely, for classical stochastic BIP [START_REF] Nouri | Statistical model checking QoS properties of systems with SBIP[END_REF] that enables to model discrete-time systems (DTMCs) and for the newly implemented Stochastic Real-Time BIP [START_REF] Nouri | Performance evaluation of stochastic real-time systems with the sbip framework[END_REF] for continuous-time systems with arbitrary distributions (GSMPs and CTMCs) 3 . The former produces untimed traces needed to verify bounded LTL properties (and to guarantee backward compatibility), whereas We implemented simulators to produce traces in different modes, i.e., symbol-wise, piece-wise and trace-wise. We use the first mode for online monitoring and to be able to interrupt simulations as soon as a verdict is obtained. The second is primordial for rare events analysis and allows to generate traces as a concatenation of trace-fragments. Finally, we use the third mode for offline monitoring.

Monitor. The new release of the tool implements monitoring capabilities for MTL and bounded LTL formulas. Our monitoring algorithms are inspired from the rewrite-based procedures introduced in [START_REF] Bulychev | Rewrite-based statistical model checking of wmtl[END_REF][START_REF] Roşu | Rewriting-based techniques for runtime verification[END_REF]. Given a formula and a trace, the monitor alternates rewriting and simplification phases. Rewriting consumes a symbol of the trace and partially evaluates the formula by unfolding temporal operators and evaluating atomic propositions to their truth value. Simplification applies Boolean reduction rules to the formula in order to conclude or to simplify it. The implemented MTL/LTL grammars and monitors allow for expressing properties with nested operators and having parameters, i.e., variables used to represent a range of properties in a compact way.

Statistical Analyses. In addition to classical SMC algorithms, i.e., HT [START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF] and PE [START_REF] Hérault | Approximate Probabilistic Model Checking[END_REF], we propose in this release two additional analyses (exploitable via independent workflows) for the exploration of properties parameters, Parametric Exploration (PX), and for rare events analysis, Importance Splitting (IP) [START_REF] Jegourel | Importance splitting for statistical model checking rare properties[END_REF]. To recall, HT allows to answer qualitative queries, i.e., given a stochastic system S and a property φ, it enables to assess whether the probability for S to satisfy φ is greater or equal to a given threshold θ. PE addresses quantitative queries, that is to compute a probability estimate p for S to satisfy φ.

Parametric Exploration (PX) is an automated way to perform statistical model checking on a family of properties, in a batch mode. A family of properties is specified in a compact way as a parametric property φ(x), where x is an integer parameter ranging over a finite instantiation domain Π. Similarly to Prism, our implemented algorithm returns a set of SMC verdicts corresponding to the verification of the parametric property instances φ(v x ) with respect to v x ∈ Π.

This can be very useful when exploring unknown system parameters such as, buffers sizes guaranteeing no overflow, or the amount of consumed energy. It automates the exploration for large parameters domains as opposed to tedious and time consuming manual procedures. This exploration differs from Uppaalsmc parametric SMC which explores the parameters of the input model.

Importance Splitting (IP) overcomes the problem of estimating the probability P (S |= φ) of a system S to satisfy a property φ representing a rare event. This is done by considering a set of intermediate levels l i that corresponds to less rare properties φ i , s.t., φ n ⇒ φ n-1 ⇒ . . . ⇒ φ 1 , where φ n = φ. P (S |= φ) is therefore computed as the product of the conditional probabilities to reach l i from l i-1 , i.e., Π n i=1 P (S |= φ i | S |= φ i-1 ). In our implementation, the intermediate levels l i and associated φ i are defined via a score function given as input. To evaluate a system trace with respect to φ, we implemented a procedure that tells the level reached by the trace, i.e., the intermediate property it satisfies. Our algorithm is similar to the analysis procedure proposed in Plasma Lab. It iterates over levels, and for each one, it simulates m trace prefixes among which m s reach the next level and m f do not. The conditional probability to reach the next level is thus estimated as the ratio ms /m. In the next iteration, the simulation of successful prefixes is resumed, while the rest (m f ) are replaced by successful ones sampled uniformly. We note that IP is currently limited to the analysis of DTMCs.

Case Studies

In this section, we briefly present experiments performed using SBIP4 . Different case studies covering various application domains were considered to validate the new release of the tool. We implemented models for communication protocols, namely Firewire, Bluetooth, and the Precision Time Protocol (PTP), for a vehicle gear controller, a Pacemaker and a mutual exclusion scenario. All the experiments were performed on a Dell Latitude 5480 with an i7-7820HQ processor and 32 GB of RAM, running Ubuntu 16.04.

On these models, we tackled different types of requirements. For the Firewire case study, we focused on analyzing its leader election protocol in different topologies (2, 3 and 5 nodes) with respect to convergence time, by considering the impact of contention (φ 1,2,3 ) and regarding the impact of a node position on its probability to become the leader (φ 4 ). In this study, except φ 3 performed using PE, the other properties were performed using PX. We also built a parametric model of the Bluetooth device discovery mechanism with one sender and one receiver that can be either in an active (v1) or a sniff mode (v2). For this model, we were interested in studying the energy consumption of the receiver in both modes (φ 6 ) in addition to the convergence time (φ 5 ). The PTP protocol was subject to the analysis of the maximal drift between the master and the slave clocks (φ 7 ). We also verified requirements regarding the time relationships between atrial and ventricular events in the pacemaker model (φ 9,10 ). Analyses of the Bluetooth, PTP and the gearbox models were performed using PX, while we used PE for the Pacemaker. We also considered a model of three concurrent processes arbitrarily requesting access to a shared resource. In this case study, the goal was to estimate the probability that each process is able to access the resource 10 times within 30 system steps (rare property φ 11 ). Using our IP implementation, we obtained 2.35 × 10 -7 in less than 13s, while it was not possible to observe the rare event using P E upon 3 minutes of execution.

Case study Model

In addition to these experiments summarized in Table 1, we report in the last two columns some performance measures of the tool, namely, the number of SMC loops performed for parametric exploration, and the average SMC time for a single loop. We observed that depending on the model size and the property complexity, the time varies from some seconds to a dozen of minutes, except for the pacemaker model where it took more than an hour. In this particular case, PE required 4883 long execution traces, representing approximately 8 minutes of real system execution.

Discussion

Most SMC tools [START_REF] Kwiatkowska | Prism 4.0: verification of probabilistic real-time systems[END_REF][START_REF] David | Uppaal smc tutorial[END_REF][START_REF] Younes | Verification and Planning for Stochastic Processes with Asynchronous Events[END_REF][START_REF] Jegourel | A platform for high performance statistical model checking -plasma[END_REF], 1] use dedicated abstract models as input for verification. In contrast, SBIP uses BIP, a full-fledged expressive component-based framework developed to support system design from specification to analysis and implementation. It allows for incrementally building complex systems from elementary components and offers real-time capabilities, in addition to highlevel coordination and synchronization primitives e.g. multi-party interactions and priorities. Furthermore, it enables including external C++ code, e.g. for modeling complex data structures and integrating legacy code.

We briefly discuss SBIP capabilities with respect to major SMC tools. Regarding the analyses, SBIP implements the HT and P E algorithms similarly to Uppaal-smc [START_REF] David | Uppaal smc tutorial[END_REF], Prism [START_REF] Kwiatkowska | Prism 4.0: verification of probabilistic real-time systems[END_REF] and Plasma Lab [START_REF] Jegourel | A platform for high performance statistical model checking -plasma[END_REF]. Besides, only Prism offers a parametric functionality similar to P X. Furthermore, to the best of our knowledge only Plasma Lab and COSMOS [START_REF] Ballarini | Hasl: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF] support rare events analysis. The former is the only one implementing IP as in our tool, while the latter rather relies on importance sampling. Our underlying modeling formalism allows for expressing arbitrary probability distributions over time. It offers built-in standard distributions, e.g. Normal, and a simple mechanism for specifying custom distributions. In contrast, Prism is restricted to uniform and exponential distributions, whereas in Uppaal-smc one need to define such distributions manually by using a subset of the C language. The expressiveness of BIP together with the reliance on concrete executions result in lower runtime performance compared to Uppaal-smc and Prism. Comparatively, the authors of Plasma Lab chose to focus on modularity at the expense of performance. In the future, we plan to optimize our simulation engine to improve the overall performance.

Fig. 1 :

 1 Fig. 1: SBIP architecture

Table 1 :

 1 Summary of performance

	For the gearbox system, we in-vestigated the minimum and maximum time required to		φ Analysis	#smc loops	avg smc time
	complete a gear change (φ 8 ).		φ 1	P X	11	1m 21s
	Firewire(2) CTMC	φ 2 φ 3	P X P E	9 -	1m 59s 2m 28s
			φ 4	P X	2	3m 27s
			φ 1	P X	17	1m 53s
	Firewire(3) CTMC	φ 2 φ 3	P X P E	11 -	3m 34s 3m 38s
			φ 4	P X	3	4m 43s
			φ 1	P X	18	3m 54s
	Firewire(5) CTMC	φ 2 φ 3	P X P E	17 12m 36s -7m 23s
			φ 4	P X	5	10m 16s
	Bluetooth v1 CTMC	φ 5 φ 6	P X P X	9 16	2m 27s 3m 11s
	Bluetooth v2 CTMC	φ 5 φ 6	P X P X	11 14 13m 05s 3m 0s
	PTP	GSMP φ 7	P X	15	8m 42s
	Gear Control CTMC φ 8	P X	11	54s
	Pacemaker CTMC	φ 9 φ 10	P E P E	--	1h 28m 1h 30m
	Mutual Exclusion	DTMC φ 11	IP P E	--	13s 3m 37s

SRT-BIP sources are available at https://gricad-gitlab.univ-grenoble-alpes. fr/verimag/bip/compiler/tree/stochastic-real-time

See details in http://www-verimag.imag.fr/TR/TR-2018-5.pdf

The research leading to these results has received funding from the EU's H2020 programme under grant agreements no. 700665 (CITADEL), 7300080 (ESROCOS)


Interval estimation for second-order delay differential equations with delayed measurements and uncertainties

The interval estimation design is studied for a secondorder delay differential equation with position delayed measurements, uncertain input and initial conditions. The proposed method contains two consecutive interval observers. The first one estimates the interval of admissible values for the position without delay for each instant of time using new delay-dependent conditions on positivity. Then derived interval estimates of the position are used to design the second observer estimating an interval of admissible values for the velocity of the considered dynamical system. The results are illustrated by numerical experiments for an example.

I. INTRODUCTION

Delays appear in many control systems at different levels of conception: at the state dynamics, in the control channel, or in the measurement signals. Sensors, actuators, and communication networks that are involved in the feedback loops usually introduce delays. Arising in differential equations, they may cause instability and oscillation of the solutions of the considered system [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF], [START_REF] Sipahi | Stability and stabilization of systems with time delay limitations and opportunities[END_REF]. For such models, which are infinite-dimensional in contrast with ordinary differential equations (ODEs), the analysis and design are much more complicated and require specially developed concepts and algorithms [START_REF] Richard | Time delay systems: an overview of some recent advances and open problems[END_REF]. For instance, the observability and methods for estimation for delayed systems with unknown inputs and nonlinearities are considered in [START_REF] Zheng | On observation of time-delay systems with unknown inputs[END_REF], [START_REF] Califano | On the observer canonical form for nonlinear time-delay systems[END_REF], application of an algebraic approach for observer design in LPV time-delay systems is presented in [START_REF] Briat | Design of LPV observers for LPV time-delay systems: an algebraic approach[END_REF], an estimation problem for positive systems with time-varying unknown delays is studied in [START_REF] Rami | Estimation of linear positive systems with unknown time-varying delays[END_REF]. Due to presence of uncertain initial conditions and disturbances, the exact estimation may be impossible, and that is why there exists another popular solution to use interval observers for uncertain systems, which provides at each instant of time a set of admissible values (an interval) for the state whose size is consistent with the system uncertainty and obtained measurements [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF], [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous time systems[END_REF], [START_REF] Kieffer | Guaranteed nonlinear state estimator for cooperative systems[END_REF], [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF], [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF], [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF].

In this work a simple benchmark problem is investigated of an unstable second order delayed system with delayed measurements (e.g. a delayed model for motion of a single mass point)

ẍ(t) = -ax(t -τ ) + f (t), t ≥ 0, y(t) = x(t -θ) + ν(t), τ = θ,
where x(t) ∈ R is the position, y(t) ∈ R is the measured output signal, f : R+ → R and ν : R+ → R are state perturbation and measurement noise (unknown bounded signals), τ > 0 and θ > 0 are the known delays. Our goal is to design an interval observer for this system, but the main issue is that beside stability conditions, to construct an interval observer it is necessary to check that the estimation error dynamics possess the positivity property. The existing solutions in the field [START_REF] Mazenc | Exponentially stable interval observers for linear systems with delay[END_REF], [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF], [START_REF] Efimov | Interval observer design for estimation and control of time-delay descriptor systems[END_REF], [START_REF] Polyakov | Output stabilization of time-varying input delay systems using interval observation technique[END_REF] are based on the delay-independent conditions of positivity from [START_REF] Haddad | Stability theory for nonnegative and compartmental dynamical systems with time delay[END_REF], [START_REF] Rami | Stability analysis and synthesis for linear positive systems with time-varying delays[END_REF]. Some results on interval observer design for uncertain timevarying delay can be found in [START_REF] Rami | Estimation of linear positive systems with unknown time-varying delays[END_REF], [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF]. In [START_REF] Efimov | Linear interval observers under delayed measurements and delaydependent positivity[END_REF] the delay-dependent conditions on positivity are introduced for the case with equal delay in the state and in the output, that may correspond to a delayfree system with delayed measurements. Therefore, the obstacle for an interval observer design is that for the considered system, by its form, the existing conditions on stability and positivity are not met. To this end, in this paper the measurement delays are supposed to be different from state delays. Using the theory of non-oscillatory solutions for second-order functional differential equations [START_REF] Agarwal | Nonoscillation theory of functional differential equations with applications[END_REF], [START_REF] Domoshnitsky | Nonoscillation, maximum principles, and exponential stability of second order delay differential equations without damping term[END_REF], the new conditions on positivity of estimation error dynamics are developed and, furthermore, a full-state interval estimation technique for considered class of systems is designed, which consists in consequent two interval observers. The first one generates the interval estimates for the position, and the second interval observer estimates the system velocity. The efficiency of proposed technique is shown by a numerical example.

The outline of this paper is as follows. After preliminaries in Section II, and an introduction of the considered time delay system properties in Section III, the interval observer design is given in Section IV. The results of the motivation example simulation are presented in Section V.

II. PRELIMINARIES

• R is the Euclidean space (R+ = {τ ∈ R : τ ≥ 0}), C n τ = C([-τ, 0], R n ) is the set of continuous maps from [-τ, 0] into R n for n ≥ 1; C n τ + = {y ∈ C n τ : y(s) ∈ R n + , s ∈ [-τ, 0]}; • xt is an element of C n τ defined as xt(s) = x(t + s) for all s ∈ [-τ, 0]; • |x| denotes the absolute value of x ∈ R, ||x||2

A. Interval relations

Given a matrix A ∈ R m×n , define A + = max{0, A}, A -= A + -A (similarly for vectors) and denote the matrix of absolute values of all elements by |A| = A + + A -.

Lemma 1. [23]

Let A ∈ R m×n be a constant matrix and x ∈ R n be a vector variable and x ≤ x ≤ x for some x, x ∈ R n , then

A + x -A -x ≤ Ax ≤ A + x -A -x. (1)

B. Delay-independent positivity

Consider a time-invariant linear system with time-varying delay:

ẋ(t) = A0x(t) + A1x(t -τ (t)) + Bf (t), t ∈ [0, +∞), (2)
x(h) = φ(h) for -τ ≤ h ≤ 0, φ ∈ C n τ , (3)
where

x(t) ∈ R n , xt ∈ C n τ is the state function; τ : R+ → [-τ , 0] is the time-varying delay, a Lebesgue measurable function of time, τ ∈ R+ is the maximum delay; f ∈ L m ∞ is the input; the constant matrices A0, A1 and B have appropriate dimensions; φ : (-τ , 0] → R n is a Borel measurable bounded function of initial conditions. Definition 1. [21] Function x : R → R n , which is locally absolutely continuous on [0, ∞), is called a solution of problem (2), (3) if it satisfies (2) for almost all t ∈ [0, ∞) and equality (3) for t ≤ 0.
The matrix A0 is called Metzler if all its off-diagonal elements are nonnegative. The system (2) is called positive if for φ ∈ C n τ + it has the corresponding solution x(t) ≥ 0 for all t ≥ 0. Lemma 2. [START_REF] Haddad | Stability theory for nonnegative and compartmental dynamical systems with time delay[END_REF], [START_REF] Rami | Stability analysis and synthesis for linear positive systems with time-varying delays[END_REF] The system (2) is positive iff A0 is Metzler, A1 ≥ 0 and Bf (t) ≥ 0 for all t ≥ 0. A positive system (2) is asymptotically stable with f (t) ≡ 0 for all t ∈ R+ iff there are p, q ∈ R n + (p > 0 and q > 0) such that

p T [A0 + A1] + q T = 0.
Under conditions of the above lemma the system has bounded solutions for f ∈ L m ∞ [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF]. Note that for linear time-invariant systems the conditions of positive invariance of polyhedral sets have been similarly given in [START_REF] Dambrine | Feedback control of timedelay systems with bounded control and state[END_REF], as well as conditions of asymptotic stability in the nonlinear case have been considered in [START_REF] Dambrine | Stability analysis of time-delay systems[END_REF], [START_REF] Dambrine | Stability and stability domains analysis for nonlinear differential-difference equations[END_REF], [START_REF] Borne | Stability Theory at the end of the XXth Century, ch. Vector Lyapunov Function: Nonlinear, Time-Varying[END_REF].

C. Representation of the solution for delay differential equations

Since τ is the maximum delay for τ (t) define a bounded set T = {s ∈ (0, τ ] : s ≤ τ (s)}, then the non-zero initial value problem (2), (3) can be rewritten to have zero initial conditions for t < 0 and with the same solution x(t) for all t ≥ 0:

ẋ(t) = A0x(t) + A1x(t -τ (t)) + Bf (t) + f * (t), t ∈ [0, +∞), (4)
x(h) = 0 for h < 0, where

f * (t) = A1φ(t -τ (t)) t ∈ T 0 otherwise .
In addition to the problem (2), (3), where x, f and φ are vector signals or functions, we will consider a problem where the solution is a n × n matrix function. For example, the n × n matrix function C(t, s) = X(t)X -1 (s), where X(t) satisfies a homogeneous initial value problem

Ẋ(t) = A0X(t) + A1X(t -τ (t)), t ≥ s, (5)
X(θ) = 0 for θ < s, X(s) = In for each s ≥ 0, is called the Cauchy matrix of [START_REF] Califano | On the observer canonical form for nonlinear time-delay systems[END_REF]. By construction C(t, s) = 0 for 0 ≤ t < s. Using this Cauchy matrix C(t, 0), a unique solution of non-homogeneous system (4) will take the form

x(t) = C(t, 0)x(0) + t 0 C(t, s) (Bf (s) + f * (s)) ds, (6)
which is also the solution for the representation (2) with the initial conditions (3).

Based on this idea, for non-zero initial function in [START_REF] Efimov | Linear interval observers under delayed measurements and delaydependent positivity[END_REF] the conditions on delay-dependent positivity are introduced by verifying an additional constrain on the first interval t ∈ [0, τ ]: Lemma 3. [START_REF] Efimov | Linear interval observers under delayed measurements and delaydependent positivity[END_REF] The system [START_REF] Sipahi | Stability and stabilization of systems with time delay limitations and opportunities[END_REF] with Bf (t) ≥ 0 for all t ≥ 0, x(0) ∈ R n + , with a Metzler matrix A0, A1 ≥ 0 and 0 ≤ (A0)i,i ≤ e(A1)i,i < (A0)i,i + τ -1 for all i = 1, . . . , n, has the corresponding solution x(t) ≥ 0 for all t ≥ 0 provided that

Bf (t) ≥ -f * (t) ∀t ∈ [0, τ ].

D. Conditions on positivity of a second order system

Following the result of [START_REF] Domoshnitsky | Nonoscillation, maximum principles, and exponential stability of second order delay differential equations without damping term[END_REF], consider the second-order delay differential equation with an input signal:

ẍ(t) + a(t)x(t -τ ) -b(t)x(t -θ) = f (t) t ∈ [0, +∞), (7)
x(h) = φ(h) for h ≤ 0, ẋ(0) ∈ R with constant delays τ, θ ≥ 0 and nonnegative functions a, b, f ∈ L∞ and φ ∈ C max{τ,θ} . The corresponding homogeneous equation is considered as

ẍ(t) + a(t)x(t -τ ) -b(t)x(t -θ) = 0 t ∈ [0, +∞). ( 8
)
For a signal q ∈ L∞ denote further the following short hands q * = essinf t≥0 q(t), q * = esssup t≥0 q(t).

Theorem 1. [START_REF] Domoshnitsky | Nonoscillation, maximum principles, and exponential stability of second order delay differential equations without damping term[END_REF] Assume that 0 ≤ τ < θ and there exists ε > 0 such that the inequalities

ε ≤ {a(t) -b(t)} ≤ 1 4 b 2 * (θ -τ ) 2 , ∀t ∈ [0, +∞) (9)
and

1 √ a * exp b * (θ-τ ) 2 4 arctan b * (θ -τ ) 2 √ a * exp b * (θ-τ ) 2 4 > θ -τ (10)
are fulfilled. Then (1) the elements C11(t, s) and C12(t, s) of the Cauchy matrix

C(t, s) = C11(t, s) C12(t, s) C21(t, s) C22(t, s) ∈ R 2×2
of (8) are nonnegative for 0 ≤ s < t < +∞;

(2) the Cauchy function C12(t, s) of ( 8) satisfies the exponential estimate

|C12(t, s)| ≤ N e -α(t-s) ∀0 ≤ s ≤ t < +∞
for some N > 0, α > 0 and the integral estimate

sup t≥0 t 0 C12(t, s)ds ≤ 1 ε ;
(3) if there exists limt→∞{a(t) -b(t)} = k, then equality

lim t→∞ t 0 C12(t, s)ds = 1 k is fulfilled.
According to this theorem, the solution x(t) of the equation ( 7) being written in the form (6) has nonnegative elements C11(t, s) and C12(t, s) of the Cauchy matrix function, and with nonnegative f (t) and the proper choice of parameters a(t), b(t) and delays τ , θ the system (7) will have a bounded solution. The Theorem 1 also establishes the exponential estimate for the Cauchy function C12(t, s) and its convergence rate. This result concerns the exponential stability of [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF], which are based on the maximum principles for the second order delay differential equation ( 8) [START_REF] Domoshnitsky | Nonoscillation, maximum principles, and exponential stability of second order delay differential equations without damping term[END_REF]. Using the representation of solutions [START_REF] Briat | Design of LPV observers for LPV time-delay systems: an algebraic approach[END_REF] (see also [START_REF] Mazenc | Exponentially stable interval observers for linear systems with delay[END_REF] below), we can see that this principle is reduced to positivity of the Cauchy matrix element C11(t, s) and C12(t, s) for nonnegative initial conditions. Proposition 1. Let 0 < τ < θ in the system [START_REF] Rami | Estimation of linear positive systems with unknown time-varying delays[END_REF] and the conditions ( 9), (10) of the Theorem 1 hold. If f (t) + f * (t) ≥ 0 for all t ≥ 0, ẋ(0) ∈ R+ and φ(h) ≥ 0 for all h ∈ [-θ, 0], where

f * (t) =      -a(t)φ(t -τ ) + b(t)φ(t -θ) t ∈ [0, τ ] b(t)φ(t -θ) t ∈ (τ, θ] 0 t > θ ,
then the corresponding solution satisfies x(t) ≥ 0 for all t ≥ 0.

III. PROBLEM STATEMENT

The main object of study in this note is the delay differential equation of the second order, which represents an unstable delayed model of motion of a single mass point:

ẍ(t) = -a(t)x(t -τ ) + f (t), t ∈ [0, +∞) (11)
where x(t) ∈ R is the position of the point, ẋ(t) ∈ R and ẍ(t) ∈ R are the velocity and the acceleration of the point motion, respectively; τ ≥ 0 is a known constant state delay, a(t) ∈ R+ with a ∈ L∞ is the parameter function, f ∈ L∞ is the external input. The equation ( 11) can be presented in a time-varying version of (2) with the corresponding matrices:

A0 = 0 1 0 0 , A1(t) = 0 0 -a(t) 0 , B = 0 1 . (12)
The initial conditions for [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF] are considered as a scalar function only for position:

x(h) = φ(h) for -τ ≤ h ≤ 0, φ ∈ Cτ ; ẋ(0) ∈ R, (13)
then the solution (6) can be rewritten to describe only the position of (11) as follows:

x(t) = C11(t, 0)x(0) + C12(t, 0) ẋ(0) (14)
+ t 0 C12(t, s) (f (s) + f * (s)) ds,
where

C(t, s) = C11(t, s) C12(t, s) C21(t, s) C22(t, s) ∈ R 2×2
with the initial condition C(0, 0) = I2, and

f * (t) = -a(t)φ(t -τ ) t ∈ [0, τ ] 0 t > τ .
Furthermore, we consider that the position of ( 11) is available for measurements with some known constant delay θ, θ ≥ τ :

y(t) = x(t -θ) + ν(t), (15)
where y(t) ∈ R is an output with measurement noise ν ∈ L∞. It is worth stressing that since θ ≥ τ , the initial conditions for (11), x(h) = φ(h), should be defined for h ∈ [-θ, 0]. To continue the analysis with this data, we also need to introduce the following hypothesis:

Assumption 1. There exist known functions φ, φ ∈ C θ such that φ(h) ≤ φ(h) ≤ φ(h) for all h ∈ [-θ, 0], and ẋ0 ≤ ẋ(0) ≤ ẋ0 for some known ẋ0 , ẋ0 ∈ R.

The assumption about a known set [φ, φ] for the initial conditions φ is standard for the interval or set-membership estimation theory [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF], [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF], [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous time systems[END_REF], [START_REF] Kieffer | Guaranteed nonlinear state estimator for cooperative systems[END_REF], [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF]. We will assume that the parameter a is known and the instant values of the signals f (t) and ν(t) are unavailable: Assumption 2. There exist known signals f , f ∈ L∞ and a constant ν0 > 0 such that f (t) ≤ f (t) ≤ f (t) and |ν(t)| ≤ ν0 for all t ≥ 0.

Therefore, the uncertain inputs f (t) and ν(t) in ( 11) and ( 15) belong to the known intervals [f (t), f (t)] and [-ν0, ν0] respectively for all t ≥ 0.

The goal is to design an interval observer for ( 11), ( 15)

ξ(t) = F [ξt, f (t), f (t), ν0, y(t)], ξt ∈ C s θ , s > 0, x(t) ẋ(t) = G[ξt, f (t), f (t), ν0, y(t)], x(t) ẋ(t) = G[ξt, f (t), f (t), ν0, y(t)]
such that for all t ≥ 0

x(t) ≤ x(t) ≤ x(t), ẋ(t) ≤ ẋ(t) ≤ ẋ(t)
provided that assumptions 1 and 2 are satisfied and x -x, xx, ẋ -ẋ, ẋ -ẋ ∈ L∞. A similar problem has been studied in [START_REF] Mazenc | Exponentially stable interval observers for linear systems with delay[END_REF], [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF], [START_REF] Efimov | Delay-dependent positivity: Application to interval observers[END_REF].

IV. INTERVAL OBSERVER DESIGN

In this section we will present two steps to design interval observers for the system [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF] and [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF]. First, using the result of Theorem 1 an observer will be given for interval estimation of the position x(t). Second, another interval observer will be designed to obtain the interval inclusion for ẋ(t).

A. The first observer for the position

The corresponding homogeneous equation

ẍ(t) = -a(t)x(t -τ ), t ∈ [0, +∞)
for [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF] has unbounded solutions in case of a constant parameter a [START_REF] Domoshnitsky | Nonoscillation, maximum principles, and exponential stability of second order delay differential equations without damping term[END_REF], which means that (11) is unstable. For the output y(t) given in [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF] the observer for [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF] can be constructed as follows: [START_REF] Efimov | Interval observer design for estimation and control of time-delay descriptor systems[END_REF] or in the state space form:

ẍ = -a(t)x(t-τ )+ (t)x(t-θ)+f (t)-(t)y(t)+ (t)ν(t),
ẋ1 = x 2 , ẋ2 = -a(t)x 1 (t -τ ) + (t)x 1 (t -θ) + f (t) -(t)y(t) + (t)ν(t),
where ∈ L∞ is an observer gain to be designed.

Remark 1. Note that, as it has been explained in Section III, the equation ( 11) can be presented in the form (2) with time-varying matrices A0, A1 and B as in [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF], and with the output in form [START_REF] Efimov | Interval estimation for uncertain systems with time-varying delays[END_REF] with different measurement delay θ ≥ τ , then the matrix A1 can be nonnegative only when θ = τ , in other cases it is always A1 < 0, then the system (11) and ( 16) do not posses the delay-independent or delay-dependent positivity properties according to lemmas 2 and 3, respectively.

Let us consider the first observer for the position x(t) of the system [START_REF] Efimov | Interval observer design for estimation and control of time-delay descriptor systems[END_REF] in the form:

ẍ-(t) = -a(t)x -(t -τ ) + (t)x -(t -θ) +f (t) -(t)y(t) -(t)ν0 -(t), (17)
ẍ+ (t) = -a(t)x + (t -τ ) + (t)x + (t -θ) +f (t) -(t)y(t) + (t)ν0 + (t),
where x -(t), x + (t) ∈ R are the estimates for the position of motion [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF] for t ∈ [0, +∞) with initial conditions

x -(h) = φ(h), x + (h) = φ(h) ∀h ∈ [-θ, 0], ẋ-(0) = ẋ0 , ẋ+ (0) = ẋ0
from Assumption 1 and

(t) = a(t)[φ(t -τ ) -φ(t -τ )] t ≤ τ 0 t > τ .
Proposition 2. Let the measurement delays satisfy the relation θ ≥ τ > 0, and assumptions 1, 2 be satisfied. For the system [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF] with initial conditions [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF] and the observer (17) select the observer gain (t) ≥ 0 to satisfy the conditions ( 9), (10) of Theorem 1 with b(t) = (t) for all t ≥ 0. Then its position satisfies the interval inclusion

x -(t) ≤ x(t) ≤ x + (t) ∀t ∈ [0, +∞) (18)
and x + -x, x -x -∈ L∞.

Using this observer it is possible to derive the interval estimates for the position x(t) without delay, but the velocity is not yet estimated since the matrix A1 < 0.

B. The second observer for the velocity

As mentioned above, the second order delay differential equation ( 11) can be presented in form of (2) with matrices [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF]. Let us consider the delayed term -a(t)x(t -τ ) as a disturbance and rewrite (2) for this case:

ẋ(t) = Ax(t) + ρ (t, x(t), x(t -τ ), f (t)) , (19)
where x(t) = [x1(t) x2(t)] T is the state vector, x1(t) is the estimated position by [START_REF] Haddad | Stability theory for nonnegative and compartmental dynamical systems with time delay[END_REF], x2(t) is a velocity of motion [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF]; According to this assumption there is a nonsingular matrix S ∈ R 2×2 such that D = S(A0 -KC)S -1 [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF], and the new coordinates z = Sx can be introduced transforming the system [START_REF] Rami | Stability analysis and synthesis for linear positive systems with time-varying delays[END_REF] to the form:

A = A0 -KC = -k1 1 -k2 0 , ρ (t, x(t), x(t -τ ), f (t)) = k1x1(t) k2x1(t) + f (t) -a(t)x1(t -τ ) , K = [k1 k2] T is
ż(t) = Dz(t) + Sρ (t, x(t), x(t -τ ), f (t)) . (20)
Using Lemma 1 we obtain that

ρ(t) ≤ ρ (t, x(t), x(t -τ ), f (t)) ≤ ρ(t),
where the functions ρ(t) and ρ(t) depend only on available information (the variables x - 1 (t), x + 1 (t), x - 1 (t -τ ) and x + 1 (t -τ ) are given by the first observer [START_REF] Polyakov | Output stabilization of time-varying input delay systems using interval observation technique[END_REF] for all t ∈ [0, +∞)):

ρ(t) = k1x - 1 (t) k2x - 1 (t) + f (t) -a(t)x + 1 (t -τ ) , ρ(t) = k1x + 1 (t) k2x + 1 (t) + f (t) -a(t)x - 1 (t -τ )
. Now, applying the results of [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF], [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with l 1 /l 2 performance[END_REF] two estimates z, z ∈ R 2 can be calculated based on the available information on these intervals (the interval inclusion (18) for x1(t) without delay), such that

z(t) ≤ z(t) ≤ z(t) ∀t ∈ [0, +∞). (21)
In other words, an interval observer can be designed for the transformed dynamics [START_REF] Efimov | Linear interval observers under delayed measurements and delaydependent positivity[END_REF]:

ż(t) = Dz(t) + S + ρ(t) -S -ρ(t), ż(t) = Dz(t) + S + ρ(t) -S -ρ(t); z(h) = S + Φ(h) -S -Φ(h), (22) z
(h) = S + Φ(h) -S -Φ(h), ∀h ∈ [-τ, 0]; x(t) = (S -1 ) + z(t) -(S -1 ) -z(t), x(t) = (S -1 ) + z(t) -(S -1 ) -z(t),
where

Φ(h) =    φ(h) 0 h < 0 ẋ0 h = 0    , Φ(h) =    φ(h) 0 h < 0 ẋ0 h = 0   
and the relations (1) are used to calculate the initial conditions for z(h), z(h) at h ∈ [-θ, 0] and the estimates x, x. Proposition 3. Let assumptions 1, 2 and 3 be satisfied. Then for the second order delay equation [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF], presented in form (2), with initial conditions 3, and with the interval observer [START_REF] Domoshnitsky | Nonoscillation, maximum principles, and exponential stability of second order delay differential equations without damping term[END_REF] the relations for the velocity

x 2 (t) ≤ x2(t) ≤ x2(t), x(t) ∈ R 2 , ∀t ∈ [0, +∞) (23)
are fulfilled provided that the conditions of Proposition 2 are verified.

V. EXAMPLE

To show the efficiency of the proposed observers we consider the motivation example [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF]:

ẍ(t) = -ax(t -τ ) + f (t), t ∈ [0, +∞)
with the values of parameter a = 2, the internal delay τ = 0.2, the perturbation f (t) = 0.5(cos(2t) + 0.3cos(10t)); the measurement delay θ = 0.5 and the noise ν(t) = ν0sin(60t) for ν0 = 0.07. The initial conditions (3) x(h) = φ(h) for h ∈ [-θ, 0], where φ(h) = 0.1 sin( 1 4πθ t), ẋ(0) = 0. The Assumption 2 is satisfied

VI. CONCLUSION

The problem of interval estimation for a second-order delay differential equation with position delayed measurements, uncertain input and initial conditions is studied in this work. The proposed approach consists in two consecutively connected interval observers. The first one estimates the set of admissible values for the position without delay using new delay-dependent conditions on positivity of a second order system. Then derived interval estimates of the position are used to design the second observer evaluating an interval of admissible values for the velocity of the considered dynamical system. The results are illustrated by numerical experiments for an example. The future directions of research will include generalization of this approach to a generic model of linear systems and control design procedure development.

Assumption 3 .

 3 a new observer gain, C = 1 0 . There are K ∈ R 2 and a Metzler matrix D ∈ R 2×2 such that the matrices A0 -KC and D have the same eigenvalues and the pairs (A0 -KC, χ1) and (D, χ2) are observable for some χ1 ∈ R 1×2 , χ2 ∈ R 1×2 .

Figure 1 .

 1 Figure 1. The result of simulation of observer (17) for the position x(t) of (11).

Figure 2 .

 2 Figure 2. Errors of estimation e 1 (t), e 2 (t) of the first observer (17).

Figure 3 .

 3 Figure 3. The result of simulation of observer (22) for the velocity ẋ(t) of (11).

Figure 4 .

 4 Figure 4. Errors of estimation e 1 (t), e 2 (t) of the second observer (22).

  Kharkovskaia is with Department of Control Systems and Informatics, Saint Petersburg State University of Information Technologies Mechanics and Optics (ITMO), 49 Kronverkskiy av., 197101 Saint Petersburg, Russia.

T. D. Efimov, A. Polyakov and J.-P. Richard are at Inria, Univ. Lille, Centrale Lille, CNRS, UMR 9189 -CRIStAL, F-59000 Lille, France E. Fridman is with School of Electrical Engineering, Tel-Aviv University, Tel-Aviv 69978, Israel. This work was partially supported by the Government of Russian Federation (Grant 08-08), the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031) and by Israel Science Foundation (grant No 1128/14).


Johannes Schiffer
email: schiffer@b-tu.de

Denis Efimov
email: denis.efimov@inria.fr

Romeo Ortega
email: ortega@lss.supelec.fr

Almost Global Synchronization in Radial Multi-Machine Power Systems

Sufficient conditions for almost global synchronization of second-order multi-machine power systems with radial topology are provided. The analysis is based on the recently proposed multivariable cell structure approach using Leonov functions-an extension of the powerful cell structure principle developed by Leonov and Noldus to nonlinear systems, the dynamics of which are periodic with respect to several state variables and possess multiple invariant solutions. The efficiency of the derived conditions is illustrated via a numerical example.

I. INTRODUCTION

A. Motivation

Undoubtly, power systems are one of the most complex man-made engineering networks [START_REF] Machowski | Power system dynamics: stability and control[END_REF]. At the same time, our modern societies crucially rely on a secure and reliable electricity supply, i.e., on a well-functioning power system. Yet, in addition to their mere extension in terms of size and number of components, power systems exhibit a rich nonlinear behavior and are permanently subjected to disturbances [START_REF] Machowski | Power system dynamics: stability and control[END_REF]- [START_REF] Kundur | Power system stability and control[END_REF]. Therefore, the task of instantaneously supplying the customer's load demand, while ensuring a safe and stable system operation is very challenging. This challenge has constantly been increasing in the past years as more and more (deterministic) conventional generation units are replaced by (volatile) renewable energy sources [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 number of transient disturbances acting on the power network has increased drastically and the system is operated under higher stress [START_REF] Winter | Pushing the limits: Europe's new grid: Innovative tools to combat transmission bottlenecks and reduced inertia[END_REF].

One of the most important problems in power system operation is to-following a disturbance-maintain all rotational generation units in synchrony, i.e., spinning at the same speed [START_REF] Kundur | Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions[END_REF]. If large disturbances (e.g., a short circuit) are considered, this is called the transient stability problem [START_REF] Kundur | Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions[END_REF]. As a consequence of a large disturbance, the resulting system's trajectories are highly influenced by the nonlinear system behavior and the generator rotor angles may exhibit severe changes [START_REF] Kundur | Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions[END_REF]. This directly implies that a local analysis is, in general, insufficient to fully reveal the systems' stability properties and dynamic behavior. Motivated by this, the present paper is devoted to the global transient stability analysis of a power system, in which synchronous generators (SGs) are represented by the standard swing equation.

B. Existing Literature

The available approaches in the literature usually depart from the following setup [START_REF] Kundur | Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions[END_REF], [START_REF] Arapostathis | Global analysis of swing dynamics[END_REF]. Consider a power system and denote its state vector by x. Suppose the system is operating in a synchronized state x pre . Then a fault occurs and is cleared after an interval of time. The post-fault power system might be different from the pre-fault system, since the faulted component may have been disconnected. Assume the post-fault system possesses a synchronized state x post . The key questions addressed in transient stability studies are if x post is attractive and if the initial condition x 0 of the postfault system is within the attractor of x post .

There are two main approaches for transient stability studies: time domain simulations and analytic Lyapunovbased approaches (also called direct energy methods) [START_REF] Pavella | Transient stability of power systems: a unified approach to assessment and control[END_REF]. Compared to a simulation-based stability assessment, Lyapunov methods have the advantages that system stability can be verified in a rigorous manner without explicitly solving the system dynamics and without screening over a large range of, possibly harmless, contingencies [START_REF] Pavella | Transient stability of power systems: a unified approach to assessment and control[END_REF]. To this end, there exists a variety of techniques, such as the identification of the (closest) unstable equilibrium point (UEP) [START_REF] Chang | Direct stability analysis of electric power systems using energy functions: theory, applications, and perspective[END_REF] or algorithms based on convex optimization [START_REF] Vu | Lyapunov functions family approach to transient stability assessment[END_REF] as well as consensus-inspired approaches [START_REF] Dörfler | Synchronization and transient stability in power networks and non-uniform Kuramoto oscillators[END_REF].

All the results discussed above pursue a local analysis of the system based on the following conjecture. Although a power system usually admits many synchronized states [START_REF] Baillieul | Geometric critical point analysis of lossless power system models[END_REF], the standard implicit assumption adopted in the literature is that the state x post of interest is the one "closest" to the prefault one x pre . Thus, also x 0 will presumably be "close" to 1x post , hence justifying the focus on local analyses [START_REF] Arapostathis | Global analysis of swing dynamics[END_REF].

Despite this (well-known) observation, there are only few works addressing a global analysis of the power system dynamics. Some global properties of power systems with generators modeled by the swing equation have been reported in [START_REF] Arapostathis | Global analysis of swing dynamics[END_REF]. By using more detailed SG models, a global stability analysis is performed in [START_REF] Caliskan | Compositional transient stability analysis of multimachine power networks[END_REF]. Yet, the result critically relies on the construction of a very specific synchronized state and a particular value for the mechanical torque of each SG in the system, both of which limit the practicality of the approach. Related results have been obtained under less stringent assumptions in [START_REF] Natarajan | Almost global asymptotic stability of a constant field current synchronous machine connected to an infinite bus[END_REF]- [START_REF] Schiffer | An input-to-state stability approach to verify almost global stability of a synchronousmachine-infinite-bus system[END_REF], but these are limited to the single-machine-infinite-bus scenario.

C. Contributions

The main contribution of the present paper is to provide sufficient conditions for almost global synchronization of a multi-machine power system with radial topology and in which SGs are represented by the standard swing equation. This is done by using the recently developed framework of multivariable cell structures and the concept of a Leonov function [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. The work in [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF] is inspired by the cell structure framework introduced for systems with a scalar periodic variable by Leonov and co-workers [START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF], [START_REF]Phase synchronisation. theory and applications[END_REF] as well as in [START_REF] Noldus | New direct Lyapunov-type method for studying synchronization problems[END_REF]. Compared to the original cell structure framework, the approach in [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF] is applicable to nonlinear systems, the dynamics of which are periodic with respect to several state variables and which possess multiple invariant solutions. Clearly, both latter properties are inherent features of multi-machine power systems.

Compared to standard Lyapunov theory, the use of Leonov functions permits to relax the usual sign definiteness requirements on the Lyapunov function and its time-derivative by exploiting the periodicity of the system. This relaxation is essential to provide conditions for global boundedness of trajectories in the present paper. By using this result together with LaSalle's invariance principle [START_REF] Van Der Schaft | L2-gain and passivity techniques in nonlinear control[END_REF], we establish almost global asymptotic stability of the attractive equilibrium set of the power system, i.e., we show that for all initial conditions, except a set of measure zero, the solutions of the power system asymptotically converge to a stable equilibrium.

At this point, it seems convenient to discuss an important technical issue with regard to our main result. As any stability analysis (numerical or analytical), the present analysis is model-based. We are aware that the employed power system model is simplified and does not necessarily capture the global behavior of a true physical power system. Nonetheless, we believe the present analysis makes a significant contribution towards a more complete understanding of the synchronization problem in large-scale power systems by providing a very different (global) perspective on the problem compared to the available literature. In addition, the results are directly applicable to the synchronization problem of second-order Kuramoto oscillators [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF], [START_REF] Rodrigues | The Kuramoto model in complex networks[END_REF] and our employed approach has also the potential to be used-with appropriate modifications-in the global analysis of more generic complex oscillator networks [START_REF] Tang | Synchronization in complex networks and its application-a survey of recent advances and challenges[END_REF]. Notation. We define the sets R ≥0 := {x ∈ R|x ≥ 0}, R >0 := {x ∈ R|x > 0} and S := [0, 2π). The set of nonnegative integers is denoted by Z ≥0 . For a set V, |V| denotes its cardinality and [V] k denotes the set of all subsets of V that contain k elements. For a set of, possibly unordered, positive natural numbers V = {l, k, . . . , n}, the short-hand i ∼ V denotes i = l, k, . . . , n. Given a positive integer n, we use 0 n to denote the vector of all zeros, 1 n the vector with all ones and I n the n × n identity matrix. Let x = col(x 1 , . . . , x n ) ∈ R n denote a column vector with entries x i ∈ R. Let diag(a i ) ∈ R n×n denote a diagonal matrix with entries a i ∈ R. For a matrix A ∈ R n×n , λ min (A) denotes its minimum eigenvalue. For a function f : R n → R, ∇f denotes its gradient and ∇ 2 f its Hessian. We employ the notation ∇f

(x * ) := ∇f (x)| x=x * and ∇ 2 f (x * ) := ∇ 2 f (x)| x=x * .
Let j denote the imaginary unit. For x ∈ R n , the element-wise sine, cosine, and tanh functions are denoted by sin

(x) ∈ [-1, 1] n , cos(x) ∈ [-1, 1] n and tanh(x) ∈ [-1, 1] n ,
respectively. Also, tanh 2 (x) denotes the square function applied element-wise to tanh(x). Furthermore, |x| = √ x x denotes the usual Euclidean norm of a vector x ∈ R n and |x| ∞ = max i |x i | its infinity norm.

II. PRELIMINARIES

A. Power System Model

We consider a radial multi-machine power system model with N > 1 nodes. The topology of the electrical network is described by an undirected and connected graph G = (N , E), where the set of nodes is denoted by N = {1, 2, . . . , N } and the set of edges (representing power lines) by E = {e 1 , . . . , e N -1 }, see [START_REF] Schiffer | A survey on modeling of microgrids -from fundamental physics to phasors and voltage sources[END_REF]. Furthermore, by associating an arbitrary ordering to the edges, we introduce the node-edge incidence matrix B ∈ R N ×N -1 , the entries of which are defined as b il = 1 if node i is the source of the l-th edge e l , b il = -1 if i is the sink of e l and b il = 0 otherwise.

To each node i ∈ N , we associate a phase angle θ i : R ≥0 → R and the corresponding electrical frequency ω i = θi . In order to establish a fundamental claim in this paper, namely convergence of bounded solutions, it is necessary to define the angles θ i on the real line, rather than (as usual) on the torus 2 . Yet, this implies that the angles are not bounded a priori, making it necessary to establish boundedness of trajectories separately, which is done in Section III-C.

Following the standard approach in transient stability studies [START_REF] Machowski | Power system dynamics: stability and control[END_REF], [START_REF] Kundur | Power system stability and control[END_REF], [START_REF] Anderson | Power System Control and Stability[END_REF], we assume that the voltage amplitudes at all nodes are positive real constants, i.e., V i ∈ R >0 and that the line admittances are purely inductive. Hence, if node i ∈ N is connected to node k ∈ N , this is represented by the nonzero susceptance B ik = B ki ∈ R <0 . If these nodes are not connected via a power line, then B ik = 0. The set of neighbors of a node i ∈ N is denoted by

N i := {k k ∈ N , k = i, B ik = 0}
. With these considerations, the active power flow P i : R |Ni| → R at node i is given by 3 [START_REF] Kundur | Power system stability and control[END_REF]

P i = G ii V 2 i + k∼Ni V i V k |B ik | sin(θ ik ), (II.1)
where we have introduced the short-hand θ ik = θ iθ k and G ii ∈ R ≥0 represents the shunt conductance at node i. Furthermore, the dynamics of the unit at the i-th node, i ∈ N , are given by the standard swing equation [START_REF] Machowski | Power system dynamics: stability and control[END_REF], [START_REF] Kundur | Power system stability and control[END_REF]

θi = ω i , M i ωi = -D i ω i -P i + P d i , (II.2)
where M i ∈ R >0 is the inertia constant, D i ∈ R >0 is the damping and droop coefficient, P d i ∈ R is the active power setpoint and the active power flow P i is given by (II.1).

In order to write the system (II.1), (II.2), i ∼ N , compactly, we introduce the matrices 2 Defining the phase angles in Euclidean space is also a common step in local Lyapunov-based stability analysis of power systems [START_REF] Vu | Lyapunov functions family approach to transient stability assessment[END_REF], [START_REF] Chiang | A BCU method for direct analysis of power system transient stability[END_REF]. 3 To simplify notation the time argument of all signals is omitted in the sequel.

M := diag(M i ) ∈ R N ×N >0 , D := diag(D i ) ∈ R N ×N >0 ,
and the vectors

θ := col(θ i ) ∈ R N , ω := col(ω i ) ∈ R N , P net := col(P d i -G ii V 2 i ) ∈ R N . With a m = V i V k |B ik |, m = 1, . . . , N -1,
we define the diagonal matrix of power line weights

A := diag(a m ) ∈ R N -1 >0 , (II.3)
and the potential function U : R N → R,

U (θ) := - {i,k}∈N ×N V i V k |B ik | cos(θ ik ) = -1 N -1 Acos(B θ).
(II.4)

The potential U (θ) and its gradient,

∇U (θ) = BAsin(B θ), (II.5)
possess certain symmetry properties, i.e., for any α ∈ R,

U (θ + α1 N ) = U (θ), ∇U (θ + α1 N ) = ∇U (θ). (II.6)
These symmetry properties arise from the fact that the power flows (II.1) only depend upon angle differences. Furthermore, since

1 N B = 0 N , 1 N ∇U (θ) = 0. (II.7)
Then, the system (II.1), (II.2), i ∼ N , can be written as

θ = ω, M ω = -Dω -∇U (θ) + P net . (II.8)
Clearly, the dynamics (II.8) are 2π-periodic in the variables θ i , i = 1, . . . , N , and hence satisfy [18, Assumption 1]. Remark 2.1: In addition to a bulk power system, the model (II.8) can also represent a microgrid with constant voltage amplitudes and power setpoints P d , see [START_REF] Schiffer | A survey on modeling of microgrids -from fundamental physics to phasors and voltage sources[END_REF], [START_REF] Schiffer | Conditions for stability of droop-controlled inverter-based microgrids[END_REF], as well as a second-order Kuramoto model [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF]. For N = 1 the model (II.8) reduces to a nonlinear pendulum.

B. Synchronized Motions and Equilibria

We employ the definition below to characterize desired invariant solutions of the system (II.8), see also [START_REF] Schiffer | Conditions for stability of droop-controlled inverter-based microgrids[END_REF].

Definition 2.2: The system (II.8) admits a synchronized motion if it has an invariant solution of the form

θ s (t) = ω s t + θ s 0 , ω s = ω * 1 N , ∀t ≥ 0, where ω * ∈ R and θ s 0 ∈ R n such that |θ s 0,i -θ s 0,k | < π 2 ∀i ∈ N , ∀k ∈ N i .
Note that the properties (II.6) have the following implication for synchronized motions of the system (II.8). If the system (II.8) possesses a synchronized motion θ s (t) = ω s t+ θ s 0 + α1 N , ω s = ω * 1 N with α = 0, then it always admits an infinite number of synchronized solutions for α ∈ R. From (II.7) and (II.8), it follows that for ωs = 0 N [START_REF] Schiffer | Conditions for stability of droop-controlled inverter-based microgrids[END_REF],

ω * = 1 N P net 1 N D1 N . (II.9)
Therefore, ω * is uniquely defined by P net and D.

Motivated by these observations and inspired by [START_REF] Arcak | Passivity as a design tool for group coordination[END_REF]- [START_REF] Trip | An internal model approach to (optimal) frequency regulation in power grids with time-varying voltages[END_REF], we introduce the new variable

η = B θ ∈ R N -1 , (II.10)
where we recall that B is the network incidence matrix and E the set of power lines. Thus, η is a projection of θ on the subspace orthogonal to 1 N and defines the phase angle differences between the nodes. Furthermore, with the change of variables (II.10) it also follows from (II.4) that

U (η) = -1 N -1 Acos(η), ∇U (η) = Asin(η), ∇ 2 U (η) = Acos(η), (II.11)
and thus from (II.5) that

∇U (θ) = BAsin(B θ) = B∇U (η).
In the new coordinates, the dynamics (II.8) are given by

η =B ω, M ω = -Dω -B∇U (η) + P net .
(II.12)

Clearly, appearance of a synchronized motion of the system (II.8) corresponds to an equilibrium of (II.12), i.e.,

η * = B (θ s 0 +ω s t+α1 N ) = B (θ s 0 +(ω * t+α)1 N ) = B θ s 0
and asymptotic stability of col(η * , ω s ) implies asymptotic convergence of the solutions col(θ, ω) to col(θ s , ω s ) up to a constant uniform shift in all angles θ av 1 N .

III. ALMOST GLOBAL SYNCHRONIZATION OF MULTI-MACHINE POWER SYSTEMS -A LEONOV

FUNCTION APPROACH This section is dedicated to the analysis of global synchronization in multi-machine power systems. Our analysis is conducted by employing the recently proposed framework of Leonov functions [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF].

A. Error Coordinates

Recall from (II.9) that the synchronization frequency ω * is uniquely defined. However, the question whether the system (II.12) possesses one or more equilibrium solutions (modulo 2π) is very hard to answer in a general setting. In fact, that problem is identical to the problem of existence of solutions to the nonlinear active power flow problem [START_REF] Baillieul | Geometric critical point analysis of lossless power system models[END_REF], [START_REF] Araposthatis | Analysis of power-flow equation[END_REF] and to that of determining the equilibria of a network of Kuramoto oscillators [START_REF] Mehta | Algebraic geometrization of the kuramoto model: Equilibria and stability analysis[END_REF]-both of which are long-standing active research areas on their own, see [START_REF] Mehta | Algebraic geometrization of the kuramoto model: Equilibria and stability analysis[END_REF]. Furthermore, it follows from [31, Theorem 2] that any equilibrium satisfying |η * 0i | < π 2 , i = 1, . . . , N -1, is locally asymptotically stable. As a consequence of the abovementioned facts and as existence of isolated equilibria is a natural prerequisite for any stability analysis, we make the following assumption.

Assumption 3.1: The system (II.12) only possesses isolated equilibria col(η * , ω * 1 N ) ∈ R 2N -1 , at least one of which is locally asymptotically stable. The Jacobian matrix of the dynamics (III.3) evaluated at any unstable equilibrium point has at least one eigenvalue with positive real part.

With Assumption 3.1, we denote an asymptotically stable equilibrium point of the system (II.12) by col(η * , 1 N ω * ) and introduce the error states

η(t) := η(t) -η * , ω(t) := ω(t) -1 N ω * . (III.1)
Furthermore, by introducing the short-hand

ζ(η) := ∇U (η + η * ) -∇U (η * ), (III.2)
the system (II.12) becomes in error coordinates

η = B ω, M ω = -D ω -Bζ(η), (III.3)
the nominal equilibrium of which is now shifted to the origin and, because of Assumption 3.1, is isolated and asymptotically stable. The remainder of this section is devoted to the stability analysis of equilibria of this system.

B. Leonov Function Candidate

The notion of a Leonov function is introduced following [17], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. For its presentation, we define two auxiliary sets:

W := {col(η, ω) ∈ R 2N -1 : |η| ∞ = c}, π ≤ c < 2π, U := ∪ r∈Z ≥0 U r , U r := {col(η, ω) ∈ R 2N -1 : ω = 0 N , |η| ∞ = 2rπ}. (III.4)
Then, a Leonov function for the system (III.3) is defined as follows [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. Definition 3.2 ( [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]): A C 1 function V : R 2N -1 → R is a Leonov function for the system (III.3) if there exist a constant g ≥ 0, functions α ∈ K ∞ , ψ ∈ K and a continuous function λ : R → R, satisfying λ(0) = 0 and λ(s)s > 0 for all s = 0, such that

α(|ω|)-ψ(|η|)-g ≤ V (η, ω) ∀col(η, ω) ∈ R 2N -1 , inf col(η,ω)∈W V (η, ω) > 0, sup col(η,ω)∈U V (η, ω) ≤ 0 (III.5)
and the following dissipation inequality holds:

V + λ(V ) ≤ 0 ∀col(η, ω) ∈ R 2N -1 . (III.6)
From (III.5) and (III.6), we see that the sign definiteness requirements of a Leonov function are relaxed compared to the standard Lyapunov function [START_REF] Van Der Schaft | L2-gain and passivity techniques in nonlinear control[END_REF], because the function V in Definition 3.2 does not have to be positive definite with respect to the variable η, i.e., the variable with respect to which the dynamics (III.3) are periodic. Furthermore, the time-derivative of V only needs to be negative definite for positive values of V. See [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF] for further details.

In addition, the requirements on (III.6) have been further relaxed in [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]Corollary 3]. This result is used to establish the main result of the present paper and hence recalled here. For this purpose, we introduce the following sets:

Ω = {col(η, ω) ∈ R 2N -1 : V ≤ 0}, Ω ε,c = {col(η, ω) ∈ R 2N -1 : V ≤ ε, |η| ∞ < c}, Z = {col(η, ω) ∈ R 2N -1 : |ω| > ξ}, (III.7)
with c defined in (III.4) and for some ε ∈ R >0 and ξ ∈ R >0 . Corollary 3.3 ( [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]): Suppose that there exists a Leonov function V : R 2N -1 → R for the system (III.3), such that sup η∈R N -1 ψ(|η|) < +∞ and the inequality (III.6) is verified only for col(η, ω)

∈ (R 2N -1 \ Ω) ∩ (Z ∪ Ω ε,c ).
Then for all initial conditions col(η(0), ω(0)) ∈ R 2N -1 the corresponding trajectories col(η, ω) are bounded ∀t ≥ 0.

d dt h = -αDM -1 h + (α(α -1)D 2 M -1 + BS(η + η * ))ω, S(η + η * ) := ∇ 2 U (η + η * )B = Acos(η + η * )B . (III.9)
Let κ ∈ R ≥0 be a parameter and Φ ∈ R N ×N , Φ = Φ > 0, be a design matrix. Then our proposed Leonov function candidate for the system (III.3) is given by

V (η, ω) :=ω M ω + h (η, ω)Φh(η, ω) -κ + 2 U (η + η * ) -U (η * ) -∇U (η * )tanh(η) = Bζ(η) ω Φ αΦD αDΦ M + α 2 DΦD :=Ψ Bζ(η) ω -κ + 2[U (η + η * ) -U (η * ) -∇U (η * )tanh(η)]. (III.10)
We also make use of this lower bound for V in the sequel:

V (η) =2[U (η + η * ) -U (η * ) -∇U (η * )tanh(η)] + λ min (Ψ)|Bζ| 2 -κ ≤ V, ∀(η, ω) ∈ R 2N -1 , (III.11)
and, in particular, of its behavior in the set

Ω = {col(η, ω) ∈ R 2N -1 : V ≤ 0, 0 < |η| ∞ ≤ c}. (III.12)

C. A Condition for Global Boundedness of Solutions

A sufficient condition for global boundedness of trajectories is presented by deriving conditions under which V in (III.10) is a Leonov function for the dynamics (III.3). To streamline the presentation of our result the following assumption is needed. Assumption 3.4: Recall the sets W and Ω and the function V defined in (III.4), (III.12) and (III.11). Consider the matrix Q(η) given in (III.13). There exist parameters

c ∈ [π, 2π), Φ > 0, α ∈ [0, 1], ν > 0, β > 0 and µ > 0, such that Q(η) >0, inf col(η,ω)∈W V (η) > 0, sup col(η,0 N )∈Ω -β|Bζ(η+η * )| 2 + 1 ν |Bdiag(tanh 2 (η))∇U (η * )| 2 < 0.
(III.14) Proposition 3.5: Consider the system (III.3) with Assumptions 3.1 and 3.4. Select

κ = 2 N -1 i=1 |a i sin(η * i )|. (III.15)
Then, the function V in (III.10) is a Leonov function for the system (III.3). Furthermore, all solutions of the system (III.3) are bounded.

Q(η) = α(ΦDM -1 + DM -1 Φ) -Φ(α(α -1)D 2 M -1 + BS(η + η * )) (-Φ(α(α -1)D 2 M -1 + BS(η + η * ))) 2D -(ν + µ)I N -β I N -αD I N -αD .
(III.13)

Proof: The claim is established by invoking Corollary 3.3. Recall the sets W and U defined in (III.4). Since tanh(η) ∈ [-1, 1] N -1 , by choosing κ as specified in (III.15), we ensure that sup col(η,ω)∈U V ≤ 0. By assumption, inf col(η,ω)∈W V > 0, which due to (III.11) implies that inf col(η,ω)∈W V > 0. Furthermore, it can be seen from (III.10) in a straightforward manner that the required functions α ∈ K ∞ and ψ ∈ K as well as the constant g ≥ 0 exist for the proposed V . Hence, the conditions in (III.5) are satisfied. Moreover, since U (η +η * ), tanh(η) and ζ(η + η * ) are bounded functions, the additional requirement sup η∈R N -1 ψ(|η|) < +∞ of Corollary 3.3 is satisfied.

Next, with (III.9) and

d dt tanh(η) = diag(1 N -1 -tanh 2 (η)) η
we have that

V ≤ - h ω Q(η) h ω -µ|ω 2 | -β|Bζ| 2 + 1 ν Bdiag(tanh 2 (η))∇U (η * ) 2 ≤ -µ|ω 2 | -β|Bζ| 2 + 1 ν Bdiag(tanh 2 (η))∇U (η * ) 2 ,
where Q(η) is given in (III.13), ν, µ as well as β are positive parameters and the last inequality follows since Q(η) > 0 by assumption.

Recall that both ζ and tanh are bounded functions, µ > 0 and |B∇U (η * )| is a constant. Hence, it is evident that there exist ξ > 0 and χ > 0, such that V + χV ≤ 0 for all col(η, ω) ∈ Z defined in (III.7). Assumption 3.4 and the facts that V (η, 0 N ) ≥ V (η, ω) for all col(η, ω) ∈ R 2N -1 and V ≤ 0 for all col(0 N -1 , ω) ∈ R 2N -1 imply that there exist (sufficiently small) parameters ε > 0 and χ > 0, such that V ≤ -χV for all

(η, ω) ∈ {col(η, ω) ∈ R 2N -1 : V ≤ V ≤ ε, |η| ∞ < c},
with V defined in (III.11). This, in turn, implies the existence of Ω ε,c with the desired properties.

Thus, all conditions of Corollary 3.3 are satisfied. Hence, V is a Leonov function for the system (III.3) and all solutions col(η, ω) ∈ R 2N -1 are bounded for all t ≥ 0.

Remark 3.6: Physically, the term B∇U (η * ) in Assumption 3.4 corresponds to the stationary network power flows. Thus, the conditions for global boundedness of trajectories in Proposition 3.5 are more likely to be satisfied in lightly loaded operating conditions. This seems reasonable from a practical point of view.

D. Main Synchronization Result

Recall that with Assumption 3.1, all equilibria of the system (III.3) are isolated and at least one of these equilibria is asymptotically stable. As the system (III.3) is continuous, this implies necessarily that some of the remaining equilibria are unstable. For an illustration see, e.g., the related numerical experiments in [START_REF] Mehta | Algebraic geometrization of the kuramoto model: Equilibria and stability analysis[END_REF]. Denote by X the set of asymptotically stable equilibria of the system (III.3).

The result below shows that, in addition to boundedness of solutions, the conditions of Proposition 3.5 also imply almost global asymptotic stability of the set X .

Theorem 3.7: Consider the system (III.3) with Assumptions 3.1 and 3.4. The set X is almost globally asymptotically stable, i.e., for all initial conditions, except a set of measure zero, the solutions of the system (III.3) asymptotically converge to a point in X .

Proof: The proof is established by using LaSalle's invariance principle [START_REF] Van Der Schaft | L2-gain and passivity techniques in nonlinear control[END_REF] together with the function

W (η, ω) = 1 2 ω M ω + U (η + η * ) -∇U (η * )η
and following the steps of [START_REF]Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF]Theorem 1]. The details are omitted for space reasons.

IV. NUMERICAL EXAMPLE The analysis is illustrated via an example derived from the IEEE 9 bus test system with three generators [28, Chapters 2 and 9]. We consider the network corresponding to the pre-fault configuration as described in [START_REF] Anderson | Power System Control and Stability[END_REF]Chapter 9]. In accordance with our assumptions, we neglect the transfer conductances in the off-diagonal entries of the reduced admittance matrix. Furthermore, the network topology in [START_REF] Anderson | Power System Control and Stability[END_REF]Chapter 9] is meshed. To obtain a radial topology, we eliminate the (smallest) off-diagonal entry of the admittance matrix, i.e., the one between generators 2 and 3. Following standard practice, we set the droop gains to D i = 1 0.05 [pu] (with respect to the rated machine powers). All other system parameters are as given in [START_REF] Anderson | Power System Control and Stability[END_REF]Chapter 9].

The conditions in Assumption 3.4 are evaluated for a range of operating points. At first we determine the parameters α, β, ν and Φ, such that Q(η) ≤ 0. From Assumption 3.4 it is evident that we seek to maximize β and ν. This is done via a polytopic approach and by implementing the corresponding matrix inequalities in Yalmip [START_REF] Löfberg | YALMIP : a toolbox for modeling and optimization in MATLAB[END_REF] with µ = 0. With the given system data, we obtain β = 2.6 • 10 -2 and ν = 33.76

Next we evaluate the feasibility of the remaining two conditions in Assumption 3.4 for a wide range of different operating points. We find that the conditions are feasible for values up to |η * 1 | = 28 °and |η * 2 | = 25 °. This shows that a reasonable range of operating points can be guaranteed to be almost globally asymptotically stable with the conditions of Proposition 3.5. In Fig. 1 almost global synchronization in multi-machine power systems with radial topology. The result has been established by combining LaSalle's invariance principle with the recently developed concept of Leonov functions [START_REF] Efimov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF], [START_REF]On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. The analysis has been illustrated via numerical experiments and it has been shown that the proposed conditions can be verified in a reasonably broad range of operating scenarios.

Usually, LaSalle-based convergence claims do not inherently provide some kind of robustness guarantees. In addition, the model employed in the analysis is a simplified representation of a true power system. Hence, two natural extensions of the presented results are to provide additional robustness measures to account for model uncertainties and extend the analysis to more detailed models with meshed topologies. Both directions are currently under investigation.

Fig. 1 .

 1 Fig. 1. Contour plot of -β|Bζ(η+η * )| 2 + 1 ν |Bdiag(tanh 2 (η))∇U (η * )| 2 for |η|∞ ≤ 2π and η * = col(22 °, -25 °). The dashed curves represent the level set V = 0. The plots show that with c = π the conditions on Ω and W in Assumption 3.4 are satisfied.

  The Leonov function candidate employed in our subsequent global analysis of invariant solutions of the dynamics (III.3) is introduced below. To this end, let

	h(η, ω) := αD ω + Bζ(η),	(III.8)
	with α ∈ [0, 1], ζ(η) defined in (III.2) and

Note that, in general, x 0 = xpre.

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

J. Schiffer acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 734832.