Pierre-Loïc Garoche

Xavier Thirioux

Loïc

Je Remercie

Dtis

Guillaume

Hugo

Kosta, Nour, Jeanne, César, Nathanaël, Pauline, Sébastien, Valentin, Sovanna, Oktay, Thomas, Quentin David

Félix

Margot Arthur

Clément Matthieu Mathias

Et Emmanuelle

Éliane, Rémy, Nicolas Juliette Anaïs

Paul Rousse

Tout d'abord, je souhaite remercier mes encadrants Pierre-Loïc Garoche et Didier Henrion. Leur expertise académique, leur soutien constant et leur optimisme sans faille durant mes travaux de thèse m'ont été précieux. Vous m'avez donné beaucoup de liberté dans le choix des sujets que j'ai abordés et je vous en suis infiniment reconnaissant. Enfin, vous m'avez permis de rencontrer différentes communautés de chercheurs (lors de collaborations et lors de conférences), c'était pour moi un sacré coup de pouce.

Cette thèse a été possible grâce au soutien financier de la Direction Générale de l'Armement (DGA) et celui de l'Office National d' Études et de Recherches Aérospatiales (ONERA). Je remercie également l' École Doctorale Aéronautique Astronautique (EDAA) et l' École Supérieure d'Aéronautique(ISAE-Supaéro).

I would like to thanks Thao Dang and Anders Rantzer for their review of this manuscript and their interesting feedback on my work. I also thank each member of the jury for attending my thesis defense which was a wonderful moment for me despite the complex sanitary condition.

Je remercie grandement Alexandre Chapoutot et Julien Alexandre Dit Sandretto de m'avoir accueilli dans leur équipe à l'ENSTA ParisTech. Travailler à vos côtés est toujours un réel plaisir, et nos discussions sont souvent fructueuses.

Notations

N natural numbers, R reals, R + positive reals, R * + strictly positive reals, R n real-valued vectors of dimension n, R n×m real-valued matrix of dimension n × m, S n×n real-valued symmetric matrix of dimension n × n, ℘(A) power set of A,

P n conic sets in R n , IR interval sets in R, x norm of IR, x = max x∈[x] x , for [x] ∈ IR n , [x] for [x] ∈ IR, let [x] = sup x∈[x] x, x
Euclidean norm of x ∈ R n , A 0 semidefinite positive matrix A ∈ S n×n , A 0 iff ∀x ∈ R n , x ⊤ Ax 0, A B semidefinite positive partial order A, B ∈ S n×n , A B iff A -B 0, T r trace, det determinant, A ⊤ transpose, A Frobenius norm, A = tr(A ⊤ A), L(I; X) continuous-time functions from I = [0, T ] ⊆ R + , T ∈ R + ∪ {∞} to X, l(I; X) discrete-time functions from I = {1, . . . , T } ⊆ N, T ∈ N ∪ {∞} to X, C 1 (I; X) continuously differentiable functions from I to X, f ∞ the ∞-norm, f ∞ = max t∈R + f (t) , L ∞ (I; X) functions from I to X, bounded by the

• ∞ norm, f the 2-norm, f = T r ∞ 0 f ⊤ (t)f (t)dt 1/2 ,
L 2 (I; X) functions from I to X, square integrable over I, L 2,loc (I; X) functions from I to X, locally square integrable over I, f, g scalar product f, g = T r

Introduction

In a project life cycle, there is a high interest to find errors at the early stages of the development. Take the example of software development. When the error is detected during the test phase, every new fix in the code requires reallocating human resources. Subparts of the faulty element might as well require new development to correct the error. After the error is fixed, the project should go back through the entire manufacturing, integration, and testing phases cycle. Thus, the cost-tofix is consequent. The cost-to-fix during the project has been analyzed within the aerospace industry in [START_REF] Haskins | Error Cost Through the Project Life Cycle[END_REF]. Table 1 reports the average cost-to-fix of an error during the project phase in software development for aerospace companies. This "exponential growth" of the cost justifies investing in solutions that can identify these errors at an early stage in the development. In the case of software development for the automotive industry, these solutions should analyze a complex class of systems: embedded systems.

An embedded system generally represents a mechanical system controlled by a computer program, e.g. an automatic car drove by an embedded computer. At a given time, such a system can be represented by its state (e.g. the position of a car, its velocity, and the state of the on-board computer). The evolution of the state Project Phase Cost-to-fix Factor Requirements 1× Design 4× Build 16× Test 61× Operations 157×

Table 1: Average cost-to-fix of errors during the project phase in software development (Table 8 in [START_REF] Haskins | Error Cost Through the Project Life Cycle[END_REF]).

1 Contents 2 through time is described by the system dynamic (e.g. a part of the dynamic of the car might be represented with the equation of motion described with the theory of rigid body, the other part might be represented as a state machine modeling the computer program). Any trajectory of the system is then a time-dependent function that associates to a time instant, the system state evolving according to the system dynamic.

Embedded systems are usually subject to rigorous safety requirements. Often, these requirements can be expressed as a property over the state space of the system. The system is safe if its trajectories avoid the set of states, the unsafe set, where this property no longer holds. If there are only a few system trajectories, it is sufficient to compute them by simulating the system in order to prove or disprove that the system is safe. When such a solution is not possible (because it is too computationally demanding, or because there is an infinite number of trajectories), verifying the safety requirement in every possible scenario becomes more challenging.

A possible approach is to verify only a subset of the trajectories. This approach is formalized using Monte-Carlo methods in [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF]. In this framework, the probability distribution of the system trajectories is estimated by randomly choosing a set of trajectories with respect to probability distribution of uncertain system parameters or inputs. Then the probability of a safety property can be numerically estimated. The Monte-Carlo method is versatile and can be applied to a large variety of systems to verify complex safety properties (e.g. temporal properties for hybrid systems in [START_REF] Sankaranarayanan | Falsification of temporal properties of hybrid systems using the crossentropy method[END_REF]). These simulationbased methods require a sample of input sets. Usually, these input sets are large. For example, they can be uncountable sets (e.g. the initial position of a car) and/or they can have an infinite number (e.g. a continuous-time real signal might model the longitudinal effect of the wind over the car). It is therefore not possible to cover all the behaviors of the system.

Another method is to manipulate infinite sets of trajectories instead of reasoning over a finite set of trajectories. In such an approach, trajectories are represented as a subset of the time-and state-space: the reachable tube. Then the system satisfies the safety properties if the reachable tube does not intersect with the unsafe set (i.e. the set of states violating the safety property).

Computing this reachable tube is usually difficult. For most systems, there exist no "out-of-the-box" methods to compute the reachable set. The first step is to abstract the dynamical systems with a class of uncertain systems that encapsulate the set of behaviors. However, even for simple systems, the reachable set is a complex geometrical object that does not admit a simple geometrical representation (e.g. the reachable set of a linear time-invariant system is not a semi-algebraic). Most of the time, the reachable tube is therefore only overapproximated, and if the safety properties are satisfied by the overapproximation, then the system satisfies these safety

Paul Rousse

Set-based Cosimulation 1. Related works properties.

In practice, to compute these overapproximations, we use computer-representable geometrical sets (e.g. intersection of hyperplanes, ellipsoidal sets, superlevel set of polynomial functions). The most precise is the overapproximation, the more properties can be proved over the reachable set. Therefore, there is a high interest to overapproximate the reachable set with a geometrical template that correctly fits the reachable set. Finding the "good" geometrical template is complex. Since these templates do not exactly represent the set of reachable states, they introduce some non-existing trajectories into the overapproximation. If these overapproximations are too pessimistic, it might be impossible to prove that the system satisfies some safety property. At the same time, these geometrical templates should not be too costly to compute. The "good" template is highly dependent on the system (its dynamic, external noises sets, the initial set of states).

This thesis proposes methods to overapproximate the reachable tube for specific classes of embedded systems.

Related works

In the case of dynamical systems where trajectories are solutions to an ordinary differential equation with bounded unknown input disturbances, the set of reachable states can be computed by solving an optimal control problem [Lee andMarkus, 1969, Gusev andZykov, 2018]. For a given state and a given cost function that associates to each initial state a positive cost (and a negative cost if the state is outside the set of initial states), if the maximal cost leading to a given state is positive, then this state is reachable. When such optimal control problem is solved (using Hamilton-Jacobi-Bellman -HJB-viscosity subsolutions, see [Soravia, 2000]), the set of states associated with a positive cost corresponds to the reachable set of the system. However, HJB solutions are difficult to compute. They rely on numerical integration of (partial) differential equations and these solutions are usually not regular.

HJB based methods propagate a cost function along with the flow of the dynamical system. Occupation measures and barrier certificates methods aim at finding constraints over the reachable tube of a dynamical system: [START_REF] Wang | Safety verification of state/time-driven hybrid systems using barrier certificates[END_REF] uses Integral Quadratic Constraints (IQC) for verification purposes using barrier certificates where the positivity of the energetic state is ensured by using a nonnegative constant multiplier: [Henrion andKorda, 2014, Korda, 2016] use an occupation measure approach where the IQC can potentially be incorporated as a constraint over the moment of the trajectories (note however that these references do not deal explicitly with IQCs). A hierarchy of semi-definite conditions is derived for polynomial dynamics. Then, off-the-shelf Semi-Definite Programming (SDP) solvers are used to solve

Set-based Cosimulation

Paul Rousse

Contents the feasibility problem. Optimization-based methods do not usually take advantage of the model structure as they consider a large class of systems (convex, Lipschitz, or polynomial dynamics for example). Similarly than for HJB methods, moment-based methods can be used for a large class of systems, but they do not scale well, i.e. they are limited to systems with a small number of states. When the dynamical system has a convex reachable set, the set can be described by an intersection of hyperplanes. These hyperplanes are obtained by finding trajectories maximizing a linear cost. This optimization problem can be solved in practice using the Pontryagin Maximum Principle -PMP-(see [START_REF] Lee | Foundations of Optimal Control Theory[END_REF], Graettinger and Krogh, 1991, Varaiya, 2000, Gusev and Zykov, 2018]). By solving this problem for different cost direction, it is possible to describe the reachable set as an intersection of hyperplanes.

The case of Linear Time-Varying (LTV) systems with ellipsoidal bounded inputs is studied in [Chernous'ko, 1999, Kurzhanski and Varaiya, 2002, Kurzhanskiy and Varaiya, 2007]. Such systems can model infinity norm bounded input-output Linear Time-Invariant (LTI) systems. The reachable set (which is convex and bounded; see [START_REF] Kurzhanski | On ellipsoidal techniques for reachability analysis. Part II: Internal approximations boxvalued constraints[END_REF]) can be overapproximated with time-varying ellipsoidal sets. Each ellipsoid is described by its parameters (center and radius) that are solutions to an Initial Value Problem (IVP). These parameters produce tight ellipsoids (i.e., ellipsoids touching the reachable set) which are external approximations of the reachable set. When multiple ellipsoids with different touching trajectories are considered, their intersection is a strictly smaller overapproximation of the reachable set. The accuracy of the overapproximation can be made arbitrarily small by adding more well-chosen ellipsoids. The exact representation of the reachable set is possible by using an uncountable set of ellipsoids.

The study of LTI systems with IQC constraint is closely related to the Linear Quadratic Regulator (LQR) problem. In the LQR problem, a quadratic integral is minimized at the terminal time. Optimal trajectories belong to a time-varying parabolic surface, whose quadratic coefficients are a solution to a Differential Riccati Equation (DRE). [START_REF] Savkin And Petersen ; Savkin | Robust state estimation for uncertain systems with averaged integral quadratic constraints[END_REF], Guseinov and Nazlipinar, 2011, Gusev and Zykov, 2018] describe the reachable set of LTI systems with terminal IQC. [Jönsson, 2002] formalizes the problem with a game theory approach. [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF] solves the differential Riccati inequality over a finite horizon using a basis of polynomial functions, then an SDP solver (such as Sedumi in [Sturm, 1999]) searches for a solution that minimizes the final volume of the overapproximation. This algorithm has been implemented in available tools (see LTVTools toolbox, [START_REF] Seiler | LTVTools (Beta), A MATLAB Toolbox for Linear Time-Varying System[END_REF]). In all these works, the overapproximation of the reachable set is conditioned by the existence of a solution to the DRE over the interval of integration. In the case of unstable systems, there exists no stable solution to the continuous algebraic Riccati equation. Any reachable set overapproximation is then defined only over a

Paul Rousse

Set-based Cosimulation

Outline and Contribution finite interval of time.

A second approach to reachability analysis is set-based methods. Here, we try to find a set invariant to the trajectories of the system. Usually, these methods are expressed as a greatest fixed point equation and are iteratively solved by removing infeasible trajectories from the set. In these frameworks, the sets are not necessarily defined as level sets and operations are focused over set operations (e.g. intersection, union, Minkowski sum, and inclusion). Contrary to previous level set methods, the family of systems that can be analyzed is usually larger since no system structure is assumed. [START_REF] Moore | Introduction to Interval Analysis[END_REF] introduced the validated numerical integration framework where numerical integration schemes (such as the Runge-Kutta one) are redefined over the set of real intervals. The system trajectories are overapproximated as a union of intervals in the state-and time-space. Other shapes (e.g. zonotopes -the projection of a hypercube-, polytopes) are used for their ability to compute the Minkowski (see [Girard, 2005]).

These fix-point algorithms are also used in compositional methods where the system is decomposed into a closed-loop composition of systems as in [START_REF] Chen | What is ellipsoidal modelling and how to use it for control and state estimation?[END_REF]. In such cases, a prior overapproximation of each internal signal is iteratively refined through computation. In such an approach, interacting dynamics between subsystems are neglected, leading to more conservative overapproximations, however, systems of higher dimensions can be treated. It has been used for stability analysis of continuous-time systems [START_REF] Platzer | [END_REF]Clarke, 2009, Eqtami andGirard, 2019] in level set-based approaches for reachability analysis.

Outline and Contribution

The thesis is organized in three parts (see Table 2), Part I and Part II derive methods to overapproximate the reachable tube for linear systems subject to bounded disturbances (bounded by a quadratic inequality over the signal space). Part I overapproximates by mean of quadratic superlevel sets, Part II overapproximates with interval sets. Part III describes a general framework to reason about the interconnection of systems.

In Part I, we study linear systems subject to bounded disturbances. The disturbance set is defined with a quadratic inequality between the disturbance signal, the state signal, and the input signal. Chapter 1 derives the general framework of overapproximation of the reachable set using the level set method with quadratic forms. Chapter 2 and Chapter 3 are applying previously derived results for two specific disturbance sets. Chapter 2 details the case where the set of disturbances is bounded at any time. Since there an infinite number of overapproximations, we show how the minimal volume overapproximation can be computed using a contin-

Set-based Cosimulation

Paul Rousse uation method. Chapter 3 studies the case where the disturbance is constrained by an Integral Quadratic Constraint (IQC). We show that the reachable set is exactly characterized as an intersection of previously defined overapproximations. The work presented in Chapter 3 has been published in the paper "Rousse, P., Garoche, P.-L., and [START_REF] Rousse | Parabolic set simulation for reachability analysis of linear time invariant systems with integral quadratic constraint[END_REF]. Parabolic set simulation for reachability analysis of linear time invariant systems with integral quadratic constraint. In 2019 18th European Control Conference, ECC 2019", and in the paper "Rousse, P., Garoche, P.-L., and [START_REF] Rousse | Parabolic Set Simulation for Reachability Analysis of Linear Time-Invariant Systems with Integral Quadratic Constraint[END_REF]. Parabolic Set Simulation for Reachability Analysis of Linear Time-Invariant Systems with Integral Quadratic Constraint. European Journal of Control".

In Part II, we study how the validated numerical integration method builds upon interval arithmetic can be used to overapproximate the reachable set of an IQC system. Chapter 4 presents the classical interval arithmetic framework and the guaranteed integration framework, and extend it to study a nonlinear IQC system. This work has been published in the paper "Rousse, P., Alexandre dit Sandretto, J., Chapoutot, A., and Garoche, P.-L. (2020a). Guaranteed Simulation of Dynamical Systems with Integral Constraints and Application on Delayed Dynamical Systems. In Lecture Notes in Computer Science, volume 11971 LNCS".

In Part III, we more specifically study the interconnection of systems In Chapter 5, we present the classical abstract interpretation framework and introduce the concrete semantic of the interconnection of systems, then, we describe abstract domains for signal spaces that will be used to represent the trajectories of the interconnection of the system. Finally, we treat several examples to compute the reachable set of the interconnection of systems. This chapter proposes a common framework for ellipsoidal methods applied for the reachability analysis of linear time-varying systems subject to bounded disturbances. Ellipsoidal methods have been studied in the '80s mainly for dynamical systems subject to ∞-norm bounded disturbances (at any time the disturbance belongs to a bounded ellipsoidal set). Since then, it has been extended disturbances subjects to state-input-disturbance inequality, namely to QC disturbances (Quadratic Constraint disturbances) and IQC disturbances (Integral Quadratic Constraint disturbances). These models have been widely used in the robust control community for stability analysis of nonlinear systems (among other applications). We show that such sets of disturbances can be described by a set of quadratic constraints in the signal space. It results in an elegant approach to present ellipsoidal methods for a wide family of models.

Part I Ellipsoidal Methods

Section 1.1 defines the system of interest, the disturbance set, and the set of reachable states. Section 1.2 introduces the level set approach in order to overapproximate the set of reachable states. Section 1.3 introduces conic templates and there corresponding quadratic value functions. Section 1.4 applies the results of Section 1.2 to define overapproximating time-varying conic sets. The time-varying coefficients of their corresponding value function are the solution of a DRE (Differential Riccati Equation) parametrized by positive multipliers. This section provides two results about overapproximating the reachable set (Theorem 1.1 and Theorem 1.2). The domain of definition of the DRE is then analyzed in Section 1.5. This chapter ends with Section 1.6 with a discussion about the conservatism of the approach (the sources of pessimism) and a comparison between the proposed approach with state of the art found in the literature.

Theorem 1.1 and Theorem 1.2 are then applied within Chapter 2 for QC systems (i.e. with ∞-norm constrained disturbances) and in Chapter 3 for IQC systems (i.e. with 2-norm constrained disturbances).

Dynamical systems

In this part, the system of interest is a linear time-varying system perturbed by an unknown disturbance that satisfies a set of constraints with the state and input trajectories. The constraints are expressed as a quadratic form in the signal space for the Definition 1.1. Dynamical system For a given input signal u

∈ L ∞ (R + ; R p ), given time-varying matrices A ∈ L ∞ (R + ; R n×n ), B ∈ L ∞ (R + ; R n×m ), C ∈ L ∞ (R + ; R n×p
), let the system S be defined as

S : ẋ(t) = A(t)x(t) + B(t)w(t) + C(t)u(t) y = (x, u, w) ∈ D (1.1)
where w ∈ L 2 (R + ; R m ) is an unknown disturbance described by D ⊆ L 2 (R + ; R n+m+p ), the set of disturbances D = {y ∈ L 2 (R + ; R n+m+p )|∀µ ∈ D * , ∀t > 0, y, y M,µ|t ≥ 0} where y, y M,µ|t = t 0 y(s) ⊤ M (s)y(s)µ(s)ds
where M is a quadratic form negative in the disturbance dimension and

D * ⊆ L 2,loc (R + ; R * + ) is a subset of positive functions from R + → R * + locally integrable over R + .

Paul Rousse

Set-based Cosimulation

Dynamical systems

The system S is a linear system since for any (x, u, w) ∈ S and α > 0, α(x, u, w) ∈ S. The set of disturbance is chosen causal since the constraint satisfied by any signal y in D, i.e. y, y M,µ|t ≥ 0, for all t ≥ 0, only depends on past values of y. One should notice that any signal y of D satisfies ∀t > 0, y, y M,η|t ≥ 0 where η = aµ + bν with µ, ν ∈ D * and a, b ≥ 0. Without loss of generality, we can therefore assume that D * is a convex cone.

Let a block decomposition of M be

M (t) =   M x (t) M x,u (t) M x,w (t) M ⊤ x,u (t) M u (t) M u,w (t) M ⊤ x,w (t) M ⊤ w,u (t) M w (t)   . (1.2)
Let the decomposition of M in the basis [x, 1, w] be such that

M x1 (t) = π ⊤ x1 M (t)π x1 (t), M w,x1 (t) = π ⊤ w M (t)π x1 (t), M w (t) = π ⊤ w M (t)π ⊤ w , M x1w (t) = π ⊤ x1w (t)M (t)π x1w (t) (1.3)
with the projections

π w =   0 0 I w   , π x1 (t) =   I x 0 0 u(t) 0 0   and π x1w (t) = π x1 (t) π w (the input u(•) is contained in the definition of M x1 (•), M w,x1 (•), and M x1w (•)).
The set of reachable states is defined as the time-varying set that associates to a time-instant t ≥ 0 the set x of states that are reachable starting from a given set of initial states X 0 subset of R n : Definition 1.

Reachable set

The set of reachable states is

R(t; X 0 ) = {x(t) | (x, u, w) ∈ S, x(0) ∈ X 0 }
where X 0 ⊂ R n is the set of initial states.

This chapter and the two following one propose a method to overapproximate the set of reachable states R(t; X 0 ).

Problem 1.1. Reachability problem

Find a P time-dependent set t ∈ R + → P(t) ⊆ R n that overapproximates R(t; X 0 ) at any time t ≥ 0, i.e. R(t; X 0 ) ⊆ P(t).

Set-based Cosimulation

Paul Rousse

Set-based simulation

In this part, the reachable set R(t; X 0 ), for t ∈ R + , is overapproximated by the superlevel set

P(t) = {x ∈ R n | p(t, x) ≥ 0} of a value function p : R + × R n → R
expressed in the time and state space domain. P(t) overapproximates R(t; X 0 ) for any t ≥ 0 if for any system trajectory x of S and any time instant t ≥ 0, it holds p(t, x(t)) ≥ 0, i.e.

x ∈ S ⇒ ∀t ∈ R + , p(t, x(t)) ≥ 0.

(1.4)

When the set of disturbances is reduced to a singleton w = 0, the conditions 1. p(0, x 0 ) is positive for any initial state x 0 ∈ X 0 ;

2. t → p(t, x(t)
) is an increasing function over R + for any system trajectory x ∈ S are sufficient to enforce (1.4), and therefore to have R(t; X 0 ) ⊆ P(t) (i.e. P(t) is an overapproximation of the set of reachable states of S, see Figure 1.1). Condition 2)

Figure 1.1: The value function p is positive (increasing and starting from zero at the initial time) along each trajectory x of the system S. The superlevel set P(t) of x → p(t, x) overapproximates the reachable set R(t; X 0 ) of the system S at any given time t > 0.

is expressed over the system trajectories. To avoid enforcing condition 2) over the set of trajectories, it can be replaced with a stronger one: for any t ∈ R + and any x ∈ R n , the value function of p increases along the flow of S. When p is continuous and differentiable in time and space, this condition is written

∂p ∂t + ∂p ∂x ⊤ f (t, x) ≥ 0, for all x ∈ R n and all t ∈ R + (1.5)

Paul Rousse

Set-based Cosimulation

Set-based simulation

where f (t, x) = A(t)x + C(t)u(t) is the vector flow of S (for the case where w = 0). When the set D is not trivial, it is still possible to use the previous approach. A sufficient condition such that (1.4) holds can be derived using a S-procedure approach (see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], Section 2.6.3). If there is a µ ∈ D * such that for any t ∈ R + p(t, x(t)) ≥ λ(t) y, y M,µ|t (1.6) where y = (x, u, w) is a system trajectory and λ : R + → R * + is a strictly positive time-dependent function, then (1.4) holds and R(t; X 0 ) ⊆ P(t). Similarly than in the trivial case, the following conditions are sufficient for (1.6) to hold:

1. p(0, x 0 ) is positive for any initial state x 0 ∈ X 0 2. t → p(t, x(t))λ(t) y, y M,µ|t is an increasing function over R + for any system trajectory x ∈ S.

When p is continuous and differentiable in time and space, and when λ is continuous and differentiable, a sufficient condition for 2) to hold is for any w t ∈ R m and any x t ∈ R n , it implies that (1.7) holds. Contrary to (1.7), (1.8) is expressed over the state space and the disturbance space. Finally, the problem of overapproximation is reduced to proving the positivity of some function over a space of finite dimension.

∂p ∂t + ∂p ∂x ⊤ f (t, x t , w t ) ≥ λ(t) y, y M,µ|t + λ(t)µ(t)y(t) ⊤ M (t)y(t) (1.7) for all x t ∈ R n and all t ∈ R + , where f (t, x t , w t ) = A(t)x t + B(t)w t + C(t)u(t)
In the general case, proving the positivity of some function is difficult. In the next part, this problem will be solved for a specific family of p, namely the time-varying quadratic forms for which proving the positivity of a function can be equivalently solved in its dual form by proving that the minimum is positive.

Remark 1.1. Relationship with Liouville equations

Equation (1.7) can be interpreted in many ways. It corresponds to the Koopman equation (dual of the Liouville equation). It is as well the value function of the minimization problem min x 0 ,y=(x,u,w) p(t, x) such that y is a system trajectory. We do not investigate further these connections in this manuscript.

Set-based Cosimulation Paul Rousse

We summarize the above discussion in Proposition 1.1.

Proposition 1.1. Level-set overapproximation For a given µ ∈ D * , and a given λ : R + → R + * continuous, differentiable, and increasing over R + , let p : R + × R n → R be a continuous and differentiable function over

R + × R n satisfying ∂p ∂t + ∂p ∂x ⊤ f (t, x, w) ≥ λ(t)λ(t) -1 p(t, x(t)) + λ(t)µ(t)y(t) ⊤ M (t)y(t) for all t ∈ R + , x ∈ R n and w ∈ R m where f (t, x, w) = A(t)x + B(t)w + C(t)u(t), y(t) = x ⊤ (t) u ⊤ (t) w ⊤ (t) ⊤ ,
Then, the 0-superlevel set P(t) of x → p(t, x), for a given t ≥ 0, defined by

P(t) = {x ∈ R n | p(t, x) ≥ 0},
overapproximates the set of reachable states of S, i.e. R(t; X 0 ) ⊆ P(t), for any t ≥ 0.

Among the trajectories of the system, some might belongs to the boundary ∂P(t) of an overapproximation P at any time t ≥ 0. Such trajectories are called touching trajectories.

Definition 1.3. Touching trajectory

A trajectory x of the system S is a touching trajectory of P if x(t) belongs to the surface of P(t) at every time t ∈ I, i.e. x(t) ∈ ∂P(t).

Since a touching trajectory is a trajectory of the system, it belongs as well to the set of reachable states. Therefore, the overapproximation P locally touches the set of reachable states.

Conic sets

The next sections apply Proposition 1.1 to time-varying conic sets.

Definition 1.4. Conic set

Let

P = P ⊤ ∈ R (n+1)×(n+1) be the coefficient of the quadratic form over R n p : R n →R x → x 1 ⊤ P x 1
.

Paul Rousse

Set-based Cosimulation

Conic sets

We define the conic set as

P = {x ∈ R n | p(x) ≥ 0}.
Let P be the set of conic sets of R n and let Conic : S (n+1)×(n+1) → P be the function that associates to any quadratic coefficient P ∈ S (n+1)×(n+1) the conic set P = Conic(P ).

Let a block decomposition of P be

P = E f f ⊤ g . (1.9) P is a conic subset of R n centered around x c = -E -1
f . E corresponds to the (signed) curvature of P. When E is not negative definite, i.e. E ≺ 0, P is unbounded, otherwise P is bounded (see Figure 1. 2). When E ≺ 0, P is a (bounded) ellipsoidal subset of R n . In this case,

P = ∅ iff g -f ⊤ E -1 f ≥ 0 (i.e.
x c belongs to P). In such a case, the ellipsoidal set can be equivalently described by the relation

x ∈ P ⇔ (x -x c ) ⊤ Q -1 (x -x c ) ≤ 1 where Q = (g -f ⊤ Ef )(-E) -1 .
The volume of P is then equal to

Vol(P ) = π n/2 (g -f ⊤ Ef )det(-E) -1 Γ(n/2 + 1)
where Γ is the gamma function.

In this report, the volume is only used as a minimizing criterion within sets of ellipsoids of fixed dimensions. Most of the time the constant in n is being neglected and since we only want to minimize the volume, we use the pseudo volume that we define by Vol(P ) = (gf ⊤ Ef )det(-E) -1 .

(1.10)

Let TrSq : S (n+1)×(n+1) → R be the map that associates to an ellipsoid P parametrized by P ∈ S (n+1)×(n+1) the squared sum of its semi-axes:

TrSq(P ) = (g -f ⊤ Ef )trace(-E -1 ).
(1.11) Remark 1.2. Representation of the conic set Conic sets could be represented differently, i.e. their center and the curvature arguments. However, we observed that alternative representations are less convenient. The time-varying conic sets are not bounded in the general case and their center is not a continuous function of time. Associated ordinary equations are most of the time difficult to manipulate and analyze.

Set-based Cosimulation

Paul Rousse When E (see the block decomposition of P in (1.9)) is negative definite (subfigure 1.2a) and P = ∅, P is an ellipsoidal set, to this respect P is a bounded and convex set. When E is not sign-definite (subfigure 1.2b), P and P = ∅, the P is unbounded.

Overapproximation with time-varying conics

This section applies Proposition 1.1 to overapproximate R(t; X 0 ) with time-varying conic sets (defined in Section 1.3). The coefficients of the time-varying conics are expressed as the solution of a DRE (Differential Riccati Equation) that is parametrized by the positive multipliers λ and µ introduced in Section 1.2. When generating all the time-varying conics for all the possible parameters (λ, µ), one obtains a family of time-dependent overapproximations of R(t; X 0 ), we define their intersection that is a tighter overapproximation. The next section will further analyze the DRE to determine the domain of definitions of these overapproximations.

Let the time-dependent conic P : R + → P associated with its value function p : R + × R n → R defined by

p(t, x) = x 1 ⊤ P (t) x 1
where P : R + → S n×n the time-varying coefficient associated with the quadratic form p.

We will use Proposition 1.1 to find sufficient conditions over Ṗ such that P overapproximates the set of reachable states R(t; P 0 ). Let q be the overapproximation of the variation of p(t, x(t))λ(t) y, y M,µt along a system trajectory y = (x, u, w) ∈ S as defined in Property 1.1

q(t, x, w) = ∂p ∂t + ∂p ∂x ⊤ f (t, x, w) -λ(t)λ(t) -1 p(t, x(t)) -λ(t)µ(t)y ⊤ M (t)y. (1.12)

Paul Rousse

Set-based Cosimulation

Overapproximation with time-varying conics

Contrary to p(t, x(t))λ(t) y, y M,µt which depends on trajectories, q only depends on states at the given time t. Moreover, q is a quadratic function of

  x 1 w   q(t, x, w) = x 1 ⊤ Ṗ (t) x 1 + 2 x 1 ⊤ P (t) A(t)x + B(t)w + C(t)u(t) 0 -λ(t)λ(t) -1 x 1 ⊤ P (t) x 1 -λ(t)µ(t)   x u(t) w   ⊤ M (t)   x u(t) w   . (1.13) Let us fix a pair (t, x) ∈ R + × R n .
To prove the positivity of h : w → q(t, x, w), we prove that the minimum of h is positive. For any t ≥ 0, M w (t) ≤ 0, therefore the minimum of the quadratic function w → q(t, x, w) exists. The minimum is reached at w * , w * (t, x) = arg min

w∈R p q(t, x, w) = λ(t) -1 µ(t) -1 M w (t) -1 B ⊤ 1 (t)P (t) -λ(t)µ(t)M w,x1 (t) x 1 .
(1.14)

Let q * (t, x) = q(t, x, w * (t, x)), we can rewrite q * as q * (t,

x) = x 1 ⊤ Q(t) x 1 (1.15)
where

Q(t) = Ṗ (t) + P (t)A 1 (t) + A ⊤ 1 (t)P (t) -λ(t) λ(t) -1 P (t) -λ(t) µ(t) M x1 (t) + λ(t) -1 µ(t) -1 B ⊤ 1 (t)P (t) -λ(t)µ(t)M w,x1 (t) ⊤ M -1 w (t) B ⊤ 1 (t)P (t) -λ(t)µ(t)M w,x1 (t)
with the following matrices

A 1 (t) = A(t) C(t)u(t) 0 0 and B 1 (t) = B(t) 0 (1.16)
and the projection and block decomposition of M into the basis x ⊤ 1 w ⊤ ⊤ is defined in (1.3). When Q(t) = 0 for all t ≥ 0, Proposition 1.1 holds.

Definition 1.5. Time-varying conic For a given µ ∈ D * , λ a continuous, differentiable function strictly positive and

Set-based Cosimulation

Paul Rousse increasing over R + , a given initial condition P 0 ∈ S (n+1)×(n+1) , let P : [0, T ] → S (n+1)×(n+1) be the solution of the initial value problem (time dependences are omitted for readability)

0 = Ṗ + P A 1 + A ⊤ 1 P -λ λ -1 P -ν M x1 + ν -1 B ⊤ 1 P -νM w,x1 ⊤ M -1 w B ⊤ 1 P -νM w,x1
(1.17)

where ν(•) = λ(•) µ(•) and γ(•) = λ(•) λ -1 (•), with initial condition P (0) = P 0 .
The ordinary differential equation (1.17) is a Differential Riccati Equation (DRE). Since all the parameters of (1.17) are measurable, the solution P to (1.17) exists locally and is unique. The convergence properties and continuity of the solution to the DRE (1.17) are studied in [Kučera, 1973]. Depending on the initial condition P 0 and on the parameters of the system, the solution P might diverge in finite-time. In such a case, the corresponding overapproximation P(•) is only defined over a finite time-horizon as well. In Section 1.5, we will show that an appropriate choice of λ and ν allows to have a solution P to DRE (1.17) defined over any time-horizon.

We can then define the time-varying conic set as

P(t) = Conic(P (t)). (1.18)
By Proposition 1.1, since along any trajectory x ∈ S, t → p(t, x) is an increasing function of time, the following property holds Theorem 1.1. Overapproximation with a time-varying conic The set of reachable states R(t; P 0 ) of S is overapproximated at any time instant t ∈ R + and for any conic set of initial of states P 0 ∈ P R(t; P 0 ) ⊆ P(t)

where P is the time-varying paraboloid defined in (1.18) with the time-varying coefficient P solution of the DRE (1.17) of Definition 1.5, A 1 (•) and B 1 (•) are defined in (1.16).

The ODE (1.17) depends on the multipliers λ ∈ C 1 (R + ; R * + ) and µ ∈ D * . For each multiplier, Theorem 1.1 ensures that the corresponding P(t; λ, µ) is a valid overapproximation of the reachable set R(t; P 0 ). We define the time-varying set Π that associates to any time

-instant t > 0 a subset Π(t) ⊆ R n Π(t) = λ ∈ C 1 (R + ; R * + ) µ ∈ D * P λ,µ (t) (1.19)

Paul Rousse

Set-based Cosimulation 1.5. Domain of definition of the overapproximations as well as the set time-varying paraboloids:

Π * = {P λ,µ (t)|λ ∈ C 1 (R + ; R * + ), µ ∈ D * }. (1.20)
The following theorem is a direct consequence of Theorem 1.1.

Theorem 1.2. Intersection of overapproximations

The set of reachable states R(t; P 0 ) of S is overapproximated by the intersection of time-varying paraboloids Π(t) generated by the sets of multipliers

(λ, µ) ∈ C 1 (R + ; R * + ) × D * , i.e. R(t; P 0 ) ⊆ Π(t),
for every t ≥ 0.

Proof. Direct consequence of Theorem 1.1 and the definition (1.19) of Π. ♦

Domain of definition of the overapproximations

In the case of reachability analysis, the existence and boundedness of the overapproximation over a given time-horizon is of great interest. This part first expresses the block decomposition of P (t). We explicit the ODEs satisfied by each block. This formulation shows the structure of the DRE (1.17) and is then used to study the domain of definition of the overapproximations P and therefore of their point-wise intersection Π.

Coefficient expansion of the DRE

The expression of (1.17) highlights the "Riccati formulation" of the ODE satisfied by P , it conveniently formulates the ODE in a one-line equation. In this part, we give an alternative useful form of (1.17) that will help to characterize the domain of definition of the overapproximations P.

Let P (•) be a solution of (1.17) for a given initial condition P 0 ∈ P. The associated

value function p(t, x) = x 1 ⊤ P (t) x 1 is a quadratic function of R n . Let E : R + →
S n×n be the time-varying quadratic coefficient of p, let f : R + → R n be the affine coefficient and let g : R + → R be the constant coefficient. Namely,

P (t) = E(t) f (t) f ⊤ (t) g(t)
.

Set-based Cosimulation Paul Rousse

Chapter 1. Ellipsoidal Methods

20

In this section, we expand (1.17) to express the ODE satisfied by E, f , and g. Expressions (1.21), (1.22), and (1.23) are more verbose than (1.17) but exhibit a more meaningful structure. Using (1.17), E satisfies DRE

0 = Ė(t) + E(t)A(t) + A ⊤ (t)E(t) -λ(t)λ(t) -1 E(t) -λ(t)µ(t)M x (t) + λ(t) -1 µ(t) -1 (B ⊤ (t)E(t) + M x,w (t)) ⊤ M w (t) -1 (B ⊤ (t)E(t) + M x,w (t)).
(1.21)

E is independent of f and g, therefore the quadratic coefficient of p is independent of the center and the offset of P 0 . f satisfies the linear time-varying differential equation

0 = ḟ (t) + A ⊤ (t)f (t) -λ(t)λ(t) -1 f (t) -(M x,u (t) + E(t)C(t))u(t) + (E(t)B(t) + M x,w )M w (t) -1 (B ⊤ (t)f (t) -M ⊤ u,w u(t)).
(1.22)

When E is measurable over R + , there exists a solution to (1.22) over R + which is bounded over any interval [0, T ], T > 0. g satisfies the linear time-varying differential equation

0 = ġ(t) -λ(t)λ(t) -1 g(t) + f (t) u(t) ⊤ G(t) f (t) u(t) (1.23)
where

G(t) = B(t)M w (t) -1 B ⊤ (t) C(t) -B(t)M w (t) -1 M u,w (t) (C(t) -B(t)M -1 w M ⊤ u,w (t)) ⊤ -M u (t) + M u,w (t)M -1 w (t)M ⊤ u,w (t)
.

Similarly than for f , when E is measurable over R + , there exists a solution f to (1.23) over R + , this solution is bounded over any interval [0, T ], T > 0.

The domain of definition of P is therefore only dependent over the domain of definition of (1.21). The differential equation (1.21) has been well studied in control. In particular, it is known that solutions of (1.21) might diverge in finite-time (independently of the regularity of the coefficient). The next section shows that under some hypothesis over D * , we can show that Π is bounded at any time t ∈ R + .

Domain of definition of the time-varying ellipsoids

Since the domain of definition of P is the domain of definition of its quadratic parameter E solution of (1.21), we study the domain of definition of E. Let the dynamical function of E in (1.21) be Ricc :

S n×n × R + × R + such that (1.21) is equivalently formulated by Ė = Ricc(E(t), λ(t), µ(t)).
A sufficient condition for E to exist over R + is that E t → Ricc(E t , λ(t), µ(t)) is Lipschitz over S n×n for every t ≥ 0 (see the Cauchy-Lipschitz Theorem in Proposition 1 of

Paul Rousse

Set-based Cosimulation 1.5. Domain of definition of the overapproximations [Zeidler, 1995a]). This is not the case since Ricc is a quadratic function in E t and thus there exists no constant that can bound the slope of this operator. However, if the solution E to (1.21) is bounded over R + , then E t → Ricc(E t , λ(t), µ(t)) is Lipschitz over the bounded set of possible solutions for every t ≥ 0.

Assumption 1.1. Well-posed disturbance set There is a K > 0 such that for any κ ≥ K, µ κ ∈ D * where µ κ (t) = exp(-κt).

Such an assumption about D * will be satisfied in the cases presented in Chapter 2 and Chapter 3. When Assumption 1.1 holds, for a κ ≥ K such that for the pair

(λ κ , µ κ ) ∈ C 1 (R + ; R * + ) × D * , with λ κ (t) = exp(κt), µ κ (t) = exp(-κt), it holds λκ (t) λ κ (t) -1 = κ and λ κ (t)µ κ (t) = 1 for any t ≥ 0.
In this section, we show that, when Assumption 1.1 holds and when E 0 ≻ 0, then there is a solution E to (1.21) that is upper-bounded E(t) E(t) (see Proposition 1.2) and lower-bounded E(t) E(t) (see Proposition 1.3) at every t ≥ 0. These two bounds are sufficient to prove that E is bounded by the matrix norm • (see Proposition 1.4). We conclude that the solution E exists over R + (Proposition 1.5).

To prove that E is upper-bounded, we use (1.21) and the fact that M w ≺ 0. By integration of (1.21) over the interval [0, t], t > 0, it holds

ψ(t, 0) ⊤ E(t)ψ(t, 0) -E 0 = t 0 ψ(t, τ ) ⊤ -M x + (B ⊤ E(τ ) + M ⊤ xw ) ⊤ M -1 w (B ⊤ E(τ ) + M ⊤ xw ) ψ(t, τ )dτ (1.24)
where ψ is the transition matrix of t → A(t) -κ 2 I. Since M w ≺ 0,

ψ(t, 0) ⊤ E(t)ψ(t, 0) E 0 - t 0 ψ(t, τ )M x ψ(t, τ )dτ.
The transition matrix ψ is invertible over [0, T ], therefore the following property holds Proposition 1.2. E's upper-bound For any t ≥ 0, E(t) E(t) where

E(t) = (ψ(t, 0) -1 ) ⊤ E 0 - t 0 ψ(t, τ )M x ψ(t, τ )dτ ψ(t, 0) -1 .
To prove that E is lower-bounded, we show that for some κ ≥ K large enough, E 0 is a lower-bound to E.

For any E 0 ≻ 0, since A and B are bounded at any time (see Definition 1.1), there is a κ ≥ K large enough s.t.

-E 0 A(t)-A(t) ⊤ E 0 -M x +(B(t) ⊤ E 0 +M ⊤ xw ) ⊤ M -1 w (B(t) ⊤ E 0 +M ⊤ xw )+κE 0 0 (1.25)
Set-based Cosimulation Paul Rousse for every t ≥ 0 Let ∆(t) = E(t) -E 0 . By definition of E, ∆(0) = 0. Using 1.25, ∆ satisfies the following ordinary differential inequality

-E 0 A(t)-A(t) ⊤ E 0 -M x +(B(t) ⊤ E 0 +M ⊤ xw ) ⊤ M -1 w (B(t) ⊤ E 0 +M ⊤ xw )+κE 0 Ė(t)-∆(t).
Using the definition of E and the one of ∆, it holds

∆(t) -∆(t)A(t) -A(t) ⊤ ∆(t) + (B(t) ⊤ E(t) + M ⊤ xw ) ⊤ M -1 w (B(t) ⊤ E(t) + M ⊤ xw ) -(B(t) ⊤ E 0 + M ⊤ xw ) ⊤ M -1 w (B(t) ⊤ E 0 + M ⊤ xw ) (1.26) We use the following expansion in (1.26) Y ⊤ RY -X ⊤ RX = -Γ ⊤ RY -Y ⊤ RΓ -Γ ⊤ RΓ with Γ = X -Y , X = B ⊤ E 0 + M ⊤ xw , Y = B ⊤ E + M ⊤ xw , and R = M -1 w . Then, (1.26) gives ∆(t) -∆(t) Ã(t) -Ã(t) ⊤ ∆(t) -∆(t) ⊤ B(t)M -1 w B(t) ⊤ ∆(t) where Ã(t) = A(t) -(B(t) ⊤ E(t) + M ⊤ xw )M -1 w for every t ≥ 0. Since M w ≺ 0, it holds ∆(t) + ∆(t) Ã(t) + Ã(t) ⊤ ∆(t) 0
By integration over [0, t], it holds ψ(t, 0) ⊤ ∆(t) ψ(t, 0) 0 where ψ(t, 0) is the transition matrix of Ã(•) from 0 to t. Since ψ(t, 0) is invertible, ∆(t) 0, i.e. E(t) E 0 .

Proposition 1.

E's lower-bound

There exists a κ ≥ K such that the solution E to (1.21) with λ(t) = exp(κt) and µ(t) = exp(-κt) is lower-bounded by E 0 , i.e. for every t ≥ 0, E(t) E(t) where

E(t) = E 0 .
Properties 1.2 and 1.3 provide an upper and lower bound for the positive semidefinite matrix order. Proposition 1.4 shows that it is a sufficient condition for E to be bounded by the • norm. 1.6. Discussion

Proof. B -B = B-A 2 + B-C 2 .
The inequality satisfied by A, B, and C gives 0 2, 1.3, and 1.4, when condition (1.25)

B -A C -A and 0 C -B C -A. Since for any X, Y ∈ S n×n , 0 X Y implies that X ≤ Y , we have B -B = B -A 2 + B -C 2 ≤ B -A 2 + B -C 2 ≤ C -A 2 + A -C 2 ≤ C -A ♦ Using Properties 1.
holds, E is bounded over its domain of definition [0, T ]. It implies that E t → Ricc(E t , λ κ (t), µ κ (t)) is Lipschitz over the set E = t∈[0,T ] E t ∈ S n×n | E t -Ẽ(t) ≤ E(t) -E(t) , where Ẽ(t) = E(t)+E(t) 2 for time instants in [0, T ]. By using a contradiction argument, if E diverges at a time instant t d ∈ R + , then E is not continuous at t d ∈ [0, T ].
Since Ricc is Lipschitz over E , this contradicts the Cauchy-Lipschitz Theorem (see Proposition 1 in [Zeidler, 1995a]) and thus E is defined over R + .

Proposition 1.5. Domain of definition of the solution to the DRE When E 0 ≻ 0 and when Assumption 1.1 holds, there is a pair of multipliers (λ, µ), λ ≥ 0 and µ ∈ D * s.t. E is defined over R + .

Discussion

Related works

Linear systems subject to disturbances bounded by quadratic constraints have been studied in the verification of dynamical systems (as in [Chaudenson, 2013]), in guaranteed state estimation (as in [Bertsekas andRhodes, 1971, Savkin andPetersen, 1995]) and in stability analysis (as in [Jönsson, 1996]). Reachability analysis of such systems has been derived within three different approaches: a set-based approach, an optimal control approach, and a level-set approach. Each method leads to fundamentally the same result which is a time-varying ellipsoid overapproximating the reachable tube and whose parameter is the solution to a Differential Riccati Equation.

In the set-based approach, the reachable set is expressed as operations (Minkowski sum and affine transformation) over the set of initial states and the disturbance set.

Set-based Cosimulation

Paul Rousse [Chernous'ko, 1999] studies the case of linear time-varying system subjects to a disturbance bounded by an ∞-norm constraint (i.e. w ⊤ (t)R(t)w(t) ≤ 1 with R(t) 0). The set of trajectories is overapproximated by a time-varying ellipsoid. Its center and radius satisfy an ordinary differential equation similar to (1.17). These differential equations are obtained using a set-based reasoning. The reachable set is described with set operations involving the initial set and the set of disturbance. More precisely the reachable set is expressed as a Minkowski sum between the flow of the initial set and the flow of the disturbance set. The first set corresponds to the image of the initial state through the autonomous dynamic and is simply an affine transformation of the set of initial states. The second set describes the influence of the exogenous disturbance w. For small time-step increase, the reachable set can be soundly approximated with an ellipsoidal set whose center and radius evolve according to an ordinary differential equation. This differential equation is parametrized by a free positive time-varying signal.

The original ellipsoidal method was firstly derived in [Schweppe, 1973]. Since then, this set-based approach has been extended to different geometrical shapes (such as zonotopes [Girard, 2005]) where the set operations could be overapproximated. The differential equation satisfied by the radius of the time-varying ellipsoid is a Differential Riccati Equation (DRE). This DRE and its associated Continuous Algebraic Riccati Equation (CARE, the equilibrium solutions of the DRE) are crucial for the control community and have been extensively studied in many works. The reader can refer to [START_REF] Bittanti | The Riccati Equation[END_REF] for an exhaustive survey.

Such an approach is difficult to reproduce when it comes to more complex disturbances as the one studied in this thesis.

The optimal control approach uses the following observation: let a real-valued function V defined over the state space be such that V (x 0 ) ≥ 0 for any x 0 in the initial state, for any reachable state x t ∈ R n , there exists a system trajectory x(•) such that x(t) = x t and V (x(0)) ≥ 0. Then, (t, x t ) belongs to the reachable set iff max V (x(0)) for x(•) a system trajectory s.t. x(t) = x t When the system is subject to bounded disturbances, the optimization problem is constrained.

Such a constrained optimization problem can be studied using the Hamilton-Jacobi-Bellman (HJB) equation. The HJB equation is a nonlinear partial differential equation (PDE). This PDE models the propagation of the cost function along the flow of the system. Solutions to the PDE are difficult to approximate in practice (its solution is often not smooth).

When the cost is quadratic, the constrained optimization problem is known as the constrained Linear Quadratic Regulator problem (LQR). [Matveev and Yakubovich,

Paul Rousse

Set-based Cosimulation 1.6. Discussion 1997] studied it using a Lagrangian relaxation. When the Lagrangian multipliers are independent of the state, i.e. are time-dependent functions, the relaxed optimization problem is the well-known LQR problem that has been extensively studied in the literature (see [START_REF] Lee | Foundations of Optimal Control Theory[END_REF], Chapter 5). [Jönsson, 2002] applied this approach to study Integral Quadratic Constraint systems (that fall into the class of system of interest of this thesis). Lagrangian multipliers are here constant weights.

The level-based approach represents an overapproximation of the reachable tube as the superlevel-set of a real-valued function defined over the time and space state. The overapproximation relationship is obtained by choosing a level-set function positive over the initial set and of increasing value along each system trajectory (see Section 1.2). Contrary to the optimal control approach, the level-set function does not solve any optimization problem.

In [START_REF] Yin | Reachability analysis using dissipation inequalities for uncertain nonlinear systems[END_REF], the authors overapproximate the reachable set of a nonlinear system subject to IQC constraints using polynomial level sets. [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF] shows that the solution to the DRE is overapproximated by solutions to a Differential Riccati Inequality (DRI). This DRI can be equivalently expressed as Differential Linear Matrix Inequality (DLMI) by using the Schur complement. Then, the authors find a solution that minimizes the input-output gain. Since the solution to the DLMI is a time-dependent matrix, the optimization problem has decision variables and constraints that belong to a signal space (of infinite dimension). To solve this optimization problem in practice, the authors express the signals in a finite-dimensional signal basis (with splines) and defined a time-sampled version of the DLMI constraint. The resulting optimization problem can be solved with semidefinite programming solvers (such as [Sturm, 1999]). The solution to the DRE can then be chosen by successively finding an optimal multiplier value with the DLMI and the DRE to minimize the input-output gain.

In previously cited works, two representation of the ellipsoidal sets are used:

• the centered representation where the ellipsoidal set E is characterized by their radius (or their curvature, i.e. the inverse of the radius), i.e.

E = {x ∈ R n | (x - c) ⊤ R -1 (x -c) ≤ 1};
• the homogeneous coordinates representation where the ellipsoidal set is expressed as the superlevel set of a quadratic form y → y ⊤ P y, P ∈ S n+1 in

homogeneous coordinates y = x 1 , E = {x ∈ R n | x 1 ⊤ P x 1 ≥ 0 in coordi- nates y = x 1 .

Set-based Cosimulation Paul Rousse

The center-radius representation is frequently used in a set-based approach (as in [Chernous'ko, 1999, Kurzhanski andVaraiya, 2014]) whereas the center-curvature is more frequent in optimal control and level set approaches (as in [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF]).

The homogeneous coordinate representation is less used in the literature (see [Savkin andPetersen, 1995, Savkin andPetersen, 1996a]).

The ODE satisfied by the coefficient of the time-varying ellipsoid depends on the chosen representation. When the disturbance set is centered, the center and radius of the ellipsoid are independent from each other (i.e. the center x c satisfies an ODE ẋc = f x (t, x c ) independent of the radius, and the radius Q satisfies another ODE Q = f Q (t, Q) independent from the center x c ). Such representation is therefore convenient for centered disturbances. However, in the general case (where disturbances are not necessarily centered), the center and radius ODEs are coupled. In our work, we found the homogeneous coordinate representation to lead to well-defined ODEs compared to the centered representation. Also, in this representation, the ODE satisfied by the coefficient of time-varying ellipsoid has an elegant form which is a DRE.

Conservatism

In this chapter, the conservatism of the ellipsoidal method is not addressed. It is studied in the two following chapters. Here, we give few hints about the different sources of conservatism and their impact over the overapproximation Π.

Incomplete dual description Few sources of uncertainties due to an incomplete dual space characterization can be identified. These incomplete dual only have a consequence of more conservatism. Since we are only interested in overapproximations, it does not compromise our analysis.

The disturbance constraint is taken into account in (1.6) by using a positive multiplier. A more general approach would choose the multiplier in the dual set of functions. And therefore, λ would be a strictly positive function of (t, x, w). Instead, (1.6) uses a multiplier λ chosen in a set of time-dependent functions. This approach allows us to derive the equations in the case of conic overapproximation case but induces some conservatism. Typically, the impact of disturbances correlated to the state is neglected.

The characterization of D * is provided to the model. In practice, disturbances are rarely described by their dual. Different disturbances will be addressed in the case of ∞-norm constraints (Chapter 2) and 2-norm constraints (Chapter 3). In practice, an underapproximation D * of D * produces an overapproximation D of D. Therefore, the reachable set R(t; X 0 ) of S with the set of disturbances D is overapproximated by the reachable set R(t; X 0 ) of S with the set of disturbances D. Overapproximation

Paul Rousse

Set-based Cosimulation 1.6. Discussion P(t) computed with Theorem 1.1 or Π(t) of R(t; X 0 ) are valid overapproximation of R(t; X 0 ).

Conic overapproximation

For some specific D * , the reachable set R(t; P 0 ) is not a pure conic set when P 0 ∈ P (see Chapter 2). Several facts motivate the choice of conic sets to describe the set of reachable states:

• when the set of disturbances is trivial (singleton set, as at the beginning of Section 1.2), the set of reachable states is exactly described by a time-varying ellipsoidal set;

• any set can be described as a (possibly non-finite, uncountable) intersection of conic sets;

• computing the minimum (and therefore proving the positivity) of a quadratic function can be done analytically;

Contrary to many works in reachability analysis, we use the general set of conics (instead of the subset of ellipsoidal sets). It allows to describe non-convex reachable sets and hopefully reduce the conservatism of the global approach.

Conclusion

This chapter proposed to overapproximate the reachable set of a linear time-varying system with time-varying conic sets. The overapproximation relationship is obtained with a "Lyapunov" approach: the value function associated with the time-varying conic is chosen to be increasing along every trajectory and positive over the set of initial states. The coefficient (in homogeneous coordinates) of the time-varying conic satisfies a parametrized Differential Riccati Equation (DRE). When some assumption about the system holds, these parameters can be chosen such that the DRE has a solution (i.e. the overapproximation exists) over the entire time-domain. Such a result is obtained by finding a lower and an upper bound to the DRE.

The next two chapters apply Theorem 1.1 and Theorem 1.2 to two subclasses of the system defined in (1.1). Chapter 2 studies linear systems subject to Quadratic Constraint (QC systems), i.e. systems where the state-input-disturbance signal y(t) = (x(t), u(t), w(t)) satisfies the quadratic constraint y ⊤ (t)M (t)y(t) ≥ 0 at any time t ≥ 0. Chapter 3 studies linear systems subject to Integral Quadratic Constraint (IQC), i.e. In Chapter 1, we formalized the general framework of reachability analysis by means of the ellipsoidal method. This chapter applies this framework to a specific subclass of disturbances: point-wise constrained disturbances. Such a model frequently appears for verification purposes since local linearization and overapproximation of the system dynamic provides such an abstraction.

This chapter is organized as follows. The QC system is presented in Section 2.1, Section 2.2 provides motivating examples and computes their reachable set. Section 2.3 applies Theorem 1.1, Chapter 1, to the QC system. Section 2.4 identifies touching trajectories and their associated support ellipsoid. Section 2.5 studies the problem of finding optimal time-varying ellipsoids (e.g. that minimize the end volume of the overapproximation), Section 2.5.2 provides a continuation method to find this optimal overapproximation. Section 2.6 presents several study cases. This chapter ends with a discussion in Section 2.7 and a conclusion in Section 2.8.

Definition of the system

In this chapter, we study an LTV system subject to a disturbance that satisfies a point-wise quadratic constraint with the state and the input. Definition 2.1. Quadratic constraint system

Let the system S be defined by

S :          ẋ = A(t)x + B(t)w + C(t)u(t)   x(t) u(t) w(t)   ⊤ M (t)   x(t) u(t) w(t)   ≥ 0 for any t ≥ 0 (2.1)
where M is a quadratic form that is negative in the disturbance dimension and s.t.

M x1\w (t) = M x1 (t) -M x1,w (t)M w (t) -1 M ⊤ x1,w (t) ≻ 0 (2.2)
for any t ≥ 0 (with the representation of M and u defined in (1.3), M x1\w is the Schur complement of the block decomposition of M ).

The assumption (2.2) ensures that the disturbance set is not empty at any time.

Abstraction

At a given time, admissible states and disturbances belong to a conic set D 1 (t) ⊂ R n×m described by

(x t , w t ) ∈ D 1 (t) ⇔   x t 1 w t   ⊤ M x1w (t)   x t 1 w t   ≥ 0.
For a given x t , t ≥ 0, let W (t, x t ) ⊂ R m be the set of admissible disturbances w t at time instant t and state x t . Using the QC constraint in Definition 2.1, W (t, x t ) is equivalently described by

w t ∈ W (t, x t ) ⇔ (w t -r t ) ⊤ R t (w t -r t ) ≤ 1 where R t = - x t 1 ⊤ M x1\w (t) x t 1 -1 M w and r t = M -1 w M ⊤ x1,w
x t 1 . Since M w ≺ 0 and M x1\w ≻ 0, it holds R t ≻ 0. To this respect, W (t, x t ) is an ellipsoidal set, bounded (since R t ≻ 0) and not empty (since R t ≺ ∞).

The following paragraphs show how the QC systems can be used as an abstraction to nonlinear systems.

Let the following system

ẋ = Ax + Bsat(c ⊤ x)
where sat is the saturation operator defined by sat : sat(y) = y when -1 ≤ y ≤ 1 sat(y) = sign(y) otherwise

Set-based Cosimulation

Paul Rousse

Let w = sat(y) and y = c ⊤ x, for any x ∈ R n , it holds

w ≥ 0 ⇒ w ≤ y w ≤ 0 ⇒ w ≥ y It implies (y -w)w ≥ 0. Which corresponds to M p =   0 0 c 2 0 0 0 c ⊤ 2 0 -1   .
The associated set of disturbance D For any y ∈ R, it also holds that w ≤ 1 and w ≥ -1. Therefore, it holds 1-w 2 ≥ 0 and

M c =   0 0 0 0 1 0 0 0 -1   .
is a valid QC constraint for the system G sat . The associated set of disturbance D 1,c (t) is a cylinder (represented in Figure 2.1b).

Any convex combinations of the QC constraint is a valid constraint as well, and therefore M λ = (1λ)M p + λM c for any λ ∈ [0, 1] is a valid QC constraint for G sat as well. For any given λ ∈ [0, 1], the associated set D 1,λ is a hyperbolic set (represented in Figure 2.1c).

A simple example

In the case of a 1-dimensional system, the reachable set of S of Definition 2.1 can be easily characterized. The reachable set corresponds to a time-varying interval whose bounds are trajectories of S. These extremal trajectories x and x are generated by disturbances that steer the state away from a center trajectory. Example 2.1 and Example 2.2 compute the reachable set for two different disturbance sets (i.e. two different QC constraints). Example 2.1 compute the reachable set for a sector inequality, Example 2.2 for an exogenous disturbance.

Example 2.1.

Let the following 1-dimensional LTI system parametrized by a k > 0

S 1 k :      ẋ = -x + w x(0) ∈ [-1, 1] w(t) 2 ≤ kx(t)
= ( √ k -1)x x(0) = 1, ẋ = ( √ k -1)x x(0) = -1.
Every other disturbance w ∈ [w, w] steers the state in-between these two trajectories x and x. The reachable set can be exactly derived

R 1 k (t) = [x, x].
We can then identify the following cases:

• when √ k -1 > 0, trajectories x and x are exponentially unstable, the reachable set is unbounded when t → ∞ and the system S 1 k is said unstable; • when √ k -1 > 0, trajectories x and x are exponentially unstable, the reachable set is unbounded when t → ∞ and the system S 1 k is said stable;

• when √ k -1 = 0, extremal trajectories satisfies x = -1 and x = 1, the reachable set is a constant tube R(t) = [-1, 1]. Example 2.2.
Let the following 1-dimensional LTI system, parametrized by a k > 0, be defined by

S 2 k :      ẋ = -x + w x(0) ∈ [-1, 1] w(t) 2 ≤ k, for all t ≥ 0 (2.4)
As in Example 2.1, extremal trajectories of S 2 k corresponds to disturbances that steer the state away of 0. I.e. the disturbance is equal to

w(t) = √ k, when x ≥ 0 and w(t) = - √ k, when x ≤ 0 Set-based Cosimulation Paul Rousse (a) S 1 0.25 (b) S 1 1 (c) S 1 1.44 Figure 2.2: Reachable set of dynamical system S 1 k , defined in (2.4), for different values of k. For k = 0.5, S 1 0.5 is a stable system, the reachable set R(t) is a tube converging to {0} when t → ∞. For k = 1, S 1 1 the reachable set R(t) is a constant tube equals to [-1, 1] for all t ≥ 0. For k = 1.2, S 1 0.
5 is an unstable system, the reachable set R(t) is a tube diverging when t → ∞.

Figure 2.3:

Reachable set of dynamical system S 2 0.04 defined in (2.4). The reachable tube is not contracting when t → ∞.

Therefore, extremal trajectories x and x satisfies

ẋ = -x + √ k, x(0) = 1, ẋ = -x - √ k, x(0) = -1.
Every other disturbance w ∈ [w, w] steers the state in-between these two trajectories x and x. The reachable set can be exactly derived

R 2 k (t) = [x, x].

Paul Rousse

Set-based Cosimulation 2.3. Ellipsoidal method

Ellipsoidal method

In order to apply Theorem 1.1 of Chapter 1 for the QC system of Definition 2.1, we first identify the dual space D * of the set of signals.

The set positive functions can be fully characterized as the subset of squareintegrable signals f ∈ L 2 (R + ; R) such that the scalar product with any positive measurable square-integrable signal µ ∈

L 2 (R + ; R) is positive. Proposition 2.1. Positive signal duality Any square-integrable measurable f ∈ L 2 (R + ; R) signal is positive over R + iff f, µ ≥ 0 for any square-integrable measurable µ positive over R + . Proof. (⇒) if f (•), µ(•) ≥ 0 over R + , then f, µ ≥ 0 (⇐) if there is a t ≥ 0 s.t. f (t) ≤ 0, then since f is measurable, there is an interval I ⊂ R + centered over t
where f, 1 I ≤ 0, therefore the property holds for µ = δ I . ♦ Therefore, the system S of Definition 2.1 is strictly equivalent to the system of Definition 1.1 where D * is the set of square-integrable functions, measurable and positive over R + . Therefore, the following Theorem 1.1 holds.

Corollary 2.1. Overapproximation of the reachable set, application of Theorem 1.1

The set of reachable states R(t; P 0 ) of S of Definition 2.1 is overapproximated at any time instant t ∈ R + and for any conic set of initial of states P 0 ∈ P with the coefficient P 0 ∈ S (n+1)×(n+1) R(t; P 0 ) ⊆ P(t)

where P is the time-varying paraboloid defined by its time-varying coefficient P solution of the DRE (time dependence is omitted for readability)

0 = Ṗ + P A 1 + A ⊤ 1 P -µ M x1 + µ -1 B ⊤ 1 P -µM w,x1 ⊤ M -1 w B ⊤ 1 P -µM w,x1 (2.5)
with initial condition P (0) = P 0 where µ is any measurable function strictly positive over R + , A 1 (•), and B 1 (•) are defined in (1.16).

Corollary 2.1 is a weaker form of Theorem 1.1 since the multiplier λ is taken as the constant function λ(t) = 1, t ≥ 0. The following parts establish that this subset

Set-based Cosimulation

Paul Rousse of overapproximations is satisfactory enough to overapproximate the reachable set R(t; P 0 ). Next sections will manipulate the dynamic function of P in (2.5). We name it below. For given t ≥ 0, P t ∈ S (n+1)×(n+1) , µ t > 0, let the operator Ricc :

R + × S (n+1)×(n+1) × R * + → S (n+1)×(n+1) be defined by Ricc(t, P t , µ t ) =P t A 1 (t) + A ⊤ 1 (t)P t -µ t M x1 (t) + µ -1 t B ⊤ 1 (t)P t -µ t M w,x1 (t) ⊤ M -1 w (t) B ⊤ 1 (t)P t -µ t M w,x1 (t) .
(2.6) Then, for a given µ(t) ∈ D * , solutions of (2.5) satisfies Ṗ (t) + Ricc(t, P (t), µ(t)) = 0 at any t ≥ 0.

Remark 2.2. Comparison with the ellipsoidal method

In the case where

M x = 0, M x,u = 0, M u = diag([0, . . . , 0, 1]), M w,u = M -1
w w c . The system S of Definition 2.1 falls into the scope of ellipsoidal methods developed in [Chernous'ko, 1999]. In these works, the reachable sets are overapproximated with time-varying ellipsoids defined by their time-dependent center x c (•) and timedependent radius Q(•). (x c , Q) and (E, f, g) are linked by the following equations

Q = (g -f ⊤ Ef ) -1 E -1 x c = -E -1 f
By deriving (x c , Q) according to time and using (2.5) (with the block decomposition (1.9)), we obtain similar differential equations than the one presented in [Chernous'ko, 1999]. Such a remark gives further insight into (2.5). The µ corresponds to a weight that either drives the system toward the disturbance direction or toward the direction of the system dynamic.

Support ellipsoids

This section defines touching trajectories and support conics. Touching trajectories are trajectories of S staying in contact with the boundary of the reachable set R(t; P 0 ). Such trajectories are associated with a time-dependent conic overapproximation P (defined in Corollary 2.1) defined such that the state x(t) of the touching trajectory stays on the boundary of P(t) at any time t ≥ 0. Since the reachable set R(t; P 0 ) and its overapproximation P(t) touches at x(t), their normals coincide at x(t) and the conic P(t) is said to be supported by the reachable set R(t; P 0 ) (see Figure 2.4).

The optimal disturbance defined in (1.14) when λ = 1 is

w(t) = µ(t) -1 M w (t) -1 B ⊤ 1 (t)P (t) -µ(t)M w,x1 (t) x 1 .
(2.7)

Paul Rousse

Set-based Cosimulation

Support ellipsoids

Figure 2.4: The touching trajectory x touches the boundary of the reachable set R(t; P 0 ) and the boundary of the overapproximation P(t) at any t ≥ 0.

Let x be the trajectory generated by the feedback w defined in (2.7). The integration of (1.15) gives p(t, x(t)) = y, y M,µt + p(0, x(0)).

(2.8)

Any point x t ∈ R n such that p(t, x t ) = 0 belongs to the boundary of the conic set P(t). Therefore, using (2.8), sufficient conditions for x to be a touching trajectory are a) y ⊤ (t)M (t)y(t) = 0 for all t ≥ 0, b) p(0, x(0)) = 0, i.e. x(0) belongs to the boundary of P(0).

Since w is chosen such that (2.7), the condition a) corresponds to an equation that depends over µ(t) -1 , with µ(t) > 0. The expansion of a) using (2.7) gives

y ⊤ (t)M (t)y(t) = z ⊤ t M x1 (t)z t + 2z ⊤ t M x1,w (t)M w (t) -1 µ(t) -1 B ⊤ 1 (t)P (t) -M w,x1 (t) z t +z ⊤ t µ(t) -1 B ⊤ 1 (t)P (t) -M w,x1 (t) M w (t) -1 µ(t) -1 B ⊤ 1 (t)P (t) -M w,x1 (t) z t (2.9)
where

z t = x(t) 1 .
If, at any time t ≥ 0, there is a strictly positive root µ * (t) to the quadratic equation (2.9), then a) is satisfied. We now show that such a µ * (t) ∈ (0, ∞) exists. Let h : R → R be the quadratic function defined by h(µ(t) -1 ) = y ⊤ (t)M (t)y(t).

When µ(t) -1 = 0, expression (2.9) becomes

h(0) = z ⊤ t (M x1 (t) -M x1,w (t)M w (t) -1 M w,x1 (t)) z t Set-based Cosimulation Paul Rousse Since M x1 (t) -M x1,w (t)M w (t) -1 M ⊤ w,x1 (t) ≻ 0, it holds h(0) > 0. Since M w ≺ 0, when µ(t) → 0, h(µ(t) -1 ) → -∞. h is a continuous function, h(0) > 0 and h(∞) = -∞,
by the intermediate value theorem, there is always a strictly positive solution µ * (t) to the equation y ⊤ (t)M (t)y(t) = 0 and therefore condition a) can be satisfied.

The above discussion is summarized in the following property.

Proposition 2.2. Touching trajectories and support conic When w, µ are chosen such that conditions a) and b) holds, then x is a touching trajectory and P is a support conic.

Proposition 2.2 defines some touching trajectory of the reachable set. Touching trajectories are defined as long as their associated support conic is defined.

Optimal ellipsoids

In this section, we address the problem of finding overapproximations minimizing some given criteria. Necessary optimality conditions are derived out of the Pontryagin's Maximum Principle (PMP) in Section 2.5.1. These conditions lead to locally optimal solutions of the initial optimization problem. Section 2.5.2 provides a numerical method to compute sub-optimal overapproximations.

In this section, we are interested in solving the following optimization problem Minimize J(P (•)) such that P (•) is a solution of (2.5) with µ(t) > 0 over [0, T ]

(2.10)

where

J : L 2 ([0, T ]; S (n+1)×(n+1) ) → R associates to a time-dependent conic P(•) a cost in R.
The cost is composed of an integral cost and a final cost as follow

J(P (•)) = Ψ(P (T )) + T 0 L(P (t))dt
The following specializations of (2.10) will be detailed within this section: a) the minimal pseudo-volume (value proportional to the squared volume, defined in (1.10)) of the terminal conic

Ψ(P T ) = Vol(P T ) L(P t ) = 0 (2.11)

Paul Rousse

Set-based Cosimulation b) the minimal pseudo-volume of the tube t → P(t)

Ψ(P T ) = 0 L(P t ) = Vol(P t ) (2.12)
c) the sum of squared semi-axis (defined in (1.11))

Ψ(P T ) = TrSq(P T ) L(P t ) = 0 (2.13) Figure 2.5 illustrates the difference between the different objectives of cases a) and c).

Application of the Pontryagin's Minimum Principle

Necessary conditions at the optimal solution of Problem (2.10) can be derived using the Minimum Pontryagin's Principle for matrices as presented in [Athans, 1967].

The PMP is an approach to solve optimal control problems, i.e. to find a control signal for system trajectory such that some cost function is minimized. To find the necessary conditions of optimality, the Hamiltonian of the dynamical system together with a co-state variable are introduced. The Hamiltonian is a storage function that preserves the final cost of a trajectory. The co-state variable measures the sensitivity

Set-based Cosimulation

Paul Rousse of the final cost to variation in the control signal. The co-state is the solution of a linear time-varying system with final state conditions. Necessary conditions of optimality imply constraints between the Hamiltonian, the co-state, and the control signal.

In the present case, the dynamical system is the time-dependent parameter of the conic and is described by differential equations (2.5), the control input is the µ parameter. Contrary to the classical application of the PMP for dynamical systems where the state is represented with vectored variables, P is a symmetric matrix realvalued variable. The co-state can be expressed as a solution to a matrix ordinary differential equation. The scalar product corresponds to the inner product of the space, i.e. trace(A ⊤ B) where A, B ∈ S n×n . Note that the matrix representation is equivalent to a vectored representation of each variable. However, the formulation is more convenient in the matrix form.

In what follows, we define the Hamiltonian H, the co-state Q, the optimal control µ * . This section ends with the necessary conditions satisfied by an optimal solution P * to (2.10). Those necessary conditions are formulated as a Boundary Value Problem (BVP) in Theorem 2.1. The next subsection provides numerical methods to compute sub-optimal solutions to (2.10) using the BVP.

Let the Hamiltonian be defined by

H(P (t), Q(t), µ(t), t) = trace(Q(t)Ricc(t, P (t), µ(t))) + L(P (t)) (2.14)
where Ricc is defined in (2.6). Let Q be the co-state solution of

Q -A Q (t)Q -QA ⊤ Q (t) -U P (t) = 0 (2.15)
where

A Q (t) = A 1 (t) + µ(t) -1 B ⊤ 1 (t)P (t) -µ(t)M w,x1 (t) ⊤ M -1 w (t)B ⊤ 1 (t)
and U P (t) = L P (P (t))

with final conditions Q(T ) = Ψ P (P (T )).

Functions L P and Ψ P can be easily computed in cases a), b) and c) ( Vol P and T P are respectively defined in (2.22) and ( 2.23), their differential are defined in Appendix 2.A):

a) L P (P t ) = 0 Ψ P (P T ) = Vol P (P T ) Paul Rousse Set-based Cosimulation 2.5. Optimal ellipsoids b) L P (P t ) = Vol P (P t ) Ψ P (P T ) = 0 c) L P (P t ) = 0 Ψ P (P T ) = TrSq P (P T )
The optimal control µ * is a solution of the following minimization problem

µ * (t) = arg min µ(t)>0 H(P (t), Q(t), µ(t), t) (2.16) Let the function h t : µ t → H(P (t), Q(t), µ t , t), it holds h t =h 0 t -µ t trace Q(t)M x1\w (t) + µ -1 t trace (Q(t)P (t)B 1 (t)M -1 w (t)B 1 (t)P (t))
where h 0 t is a term independent of µ t . Deriving h t gives the following minimizer µ * (t) for (2.16)

µ * (t) = -trace Q(t)P (t)B 1 (t)M -1 w (t)B ⊤ 1 (t)P (t) trace Q(t)M x1\w (t) 1/2 (2.17) Since, M w (t) ≺ 0 and M x1\w (t) ≻ 0, when Q(t) ≻ 0, µ * (t)
is correctly defined (i.e, the square root exists) and the minimization problem (2.16) does have a solution.

Theorem 2.1. Necessary conditions for optimality

Let P * be an optimal trajectory associated with the optimal control µ * as defined in (2.17) of (2.10).

Necessary conditions for P * are expressed as the existence of a solution to a boundary value problem. The co-state is a symmetric matrix function.

A continuation method to solve the PMP

The BVP defined in Theorem 2.1 is not easy to solve. Solutions of (2.5) might have a finite escape time. Even when the solution E(•) is correctly defined for a given µ(•), there are no guarantees about the sign of E(•). Considering the different optimality criterion presented in Section 2.5.1, when E(•) is not invertible along the trajectory, the integration of the co-state variable is compromised.

In this part, we investigate the use of a continuation method to solve the BVP in Theorem 2.1.

The BVP can be equivalently defined as the following root-finding problem parametrized by the final time of integration T > 0

F (Q 0 , T ) = Q(T ) -Ψ P (P (T )).
(2.18)

Set-based Cosimulation Paul Rousse where P is a solution of (2.5) for a given initial state Q 0 , Q is a solution of (2.15) with the initial state

Q(0) = Q 0 .
The continuation method proposes to find the solution Q 0,T f of (2.18) where T = T f , T f > 0 given, by following the curve of solutions

Q : [0, T f ] → S (n+1)×(n+1) T → Q 0,T (2.19)
from a trivial point (that will be at T = 0 in our case) to the point of interest T = T f . When the curve Q is regular enough, a prediction correction algorithm can be used to follow the curve over [0, T f ] (see Figure 2.6).

When T = 0, the root Q 0 of the equation F (Q 0 , 0) = 0 is indeed trivial (using (2.18)).

Proposition 2.3. Initial point the co-state

The solution of (2.18) for T = 0 is Q(0) = Ψ P (P 0 ).

The regularity of the curve Q is difficult to assess globally since the curve is an implicit function; it would require to study the regularity of the inverse function of

Q 0 → F (Q 0 , T ) for T ∈ [0, T f ]. However, for a given point (T, Q 0,T ) on the curve Q, if F is linearizable at (T, Q(T )) and if Q 0 → F (Q 0 , T
) can be inverted, then the curve Q exists in the neighborhood of T . This is stated in the implicit function theorem as described in [Zeidler, 1995b, Chapter 4.8].

The linearization F of F can be derived using variational calculus (see Fréchet derivatives [Zeidler, 1995b, Chapter 2.1]). Let R : [0, T ] → S (n+1)×(n+1) and S :

Paul Rousse

Set-based Cosimulation

Example

[0, T ] → S (n+1)×(n+1) (resp.) be variations of P and Q (resp.) for a given variation S 0 ∈ S (n+1)×(n+1) to the initial co-state Q 0 . R and S are solutions of an LTV system, details of how (R, S) are derived are given in Annexe 2.B. The linearization F along (S 0 , dt) is

F = ∆F Q 0 • S 0 + ∆F T dt where ∆F Q 0 • S 0 = R(T ) -∆Ψ P (P (T )) • S(T ) ∆F T = Q(T ) -∆Ψ P (P (T )) • Ṗ (T )
where the Fréchet derivatives ∆Ψ P (P ) • S are defined in Appendix 2.A.

Proposition 2.4. Tangent of the co-state When ∆F Q 0 is invertible, the tangent S 0 of Q satisfies 0 = ∆F Q 0 • S 0 + ∆F T . (2.20)
The BVP is solved for an increasing sequence of time-horizons {T k } until the time-horizon T f is reached. For each time -horizon T k , we compute the solution of the BVP problem, i.e. the initial state of the co-state Q 0,T k . We use the continuity of the curve (guaranteed by using the inverse function theorem) to predict the next Q 0,T k+1 .

Remark 2.3. Corrector step

The corrector step numerically solving (2.18) uses a BVP solver such as bvp5c [START_REF] Kierzenka | A BVP solver that controls residual and error[END_REF].

Remark 2.4. Complexity of the BVP algorithm P , Q, R and S are symmetric matrices. Therefore, the differential equations have a 2n(n + 1) dimensional state.

Example

We study a linear time-invariant system of two dimensions defined by (2.1) with the following parameters The gray area corresponds to the minimum volume ellipsoid. Since the Pontryagin's Maximum Principle only gives necessary conditions for optimality, and since the ellipsoidal method only provides overapproximations, the E vol is not the actual minimum volume ellipsoid overapproximating the reachable tube R(T ; E 0 ). The dotted line represents an overapproximating hyperbola touching the reachable set.

A = -2 1 -3 -3 B = 1 1 ⊤ C = 1 1 ⊤ u = 0 M = diag 1 2 0 -1 (2.21) Set-based Cosimulation Paul Rousse
The reachable set at T = 1 and the minimal volume overapproximating ellipsoid are drawn in Figure 2.7. The minimal volume overapproximating ellipsoid is drawn at different times in Figure 2.8.

Discussion

The QC constraint is a ∞-norm constraint over the state-input-disturbance signal. Such QC system (or subfamilies of QC systems) has been studied many times in the literature. First studies goes back to the 60's with the unknown-but-bounded disturbance model (in [Schweppe, 1973], Section 7.5), in this model the disturbance belongs to an ellipsoidal set at any time, i.e. (w(t) -

w c (t)) ⊤ R(t)(w(t) -w c (t)) ≤ 1.
In other works, the inequality

αv 2 ≤ vw ≤ βv 2 which is equivalent to the QC constraint w -α+β 2 ≤ β-α 2 v
, is referred as the sector inequality (as in [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF]). [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF] (in Section 4.2.3) studies the stability of the norm-bound linear time-invariant system, a system subject to disturbance w = ∆x where ∆

(t) ≤ 1. Such constraint is equivalent to the QC constraint w(t) ≤ x(t) . Any QC constraint y(•) ⊤ M (•)y(•) ≥ 0 is equivalent to the set of IQC constraints ∀µ(•) ≥ 0, T 0 y(t) ⊤ M (t)y(t)µdt ≥ 0 for any T ∈ R + ∪ {+∞}. For this reason,

Paul Rousse

Set-based Cosimulation

2.7. Discussion (a) t = 0 (b) t = 0.2 (c) t = 0.4 (d) t = 0.6 (e) t = 0.8 (f ) t = 1
Figure 2.8: Reachable set R (in red) of the dynamical system S, defined in Section 2.6, at different time instant t. V opt is the level-set associated with the minimal volume overapproximating ellipsoids and x * is its associated touching trajectory.

Set-based Cosimulation

Paul Rousse the study of QC system is closely related to the study of IQC systems. Moreover, the integration of the QC constraint is a necessary step that happens during the Lagrangian relaxation of the LQR problem. [START_REF] Bertsekas | Recursive State Estimation for a Set-Membership Description of Uncertainty[END_REF]] identified this link in the case of guaranteed state estimation where the system is subject to hard bound constraint or to energy constraint. In the robust control community, a lot of effort in the study of QC constraints has been dedicated to the search of a class of multiplier µ that can be expressed as state-space model of finite dimension (see [START_REF] Veenman | Robust stability and performance analysis based on integral quadratic constraints[END_REF]] Section 6).

The most common way to handle QC constraints is the use of constant multipliers, i.e. µ(•) = µ 0 > 0 as in [Jönsson, 1996]. Such an approach is commonly referred in the literature as the S-procedure (see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], Section 2.6.3).

More recently, [START_REF] Fetzer | Invariance with dynamic multipliers[END_REF], Veenman et al., 2016] used so-called dynamic multipliers that allowed to reduce the conservatism of robust stability analysis and reachability analysis of IQC systems. The inequality w(t) 2 ≤ v(t) 2 is equivalently replaced by w, Hw ≤ v, Hv where H is a positive operator over the set of signal (i.e. a positive transfer function).

For QC systems, multipliers belong to the set of positive functions over the time interval. Each multiplier is associated with a valid overapproximation of the reachable set. Since the set of multipliers is infinite, it is interesting to find one "good" overapproximation.

In [Chernous'ko, 1999], the author tries to search for minimal volume overapproximation (for the case of exogenous disturbances, i.e. w ⊤ (t)R(t)w(t) ≤ 1 with R(t) 0). This optimization problem is an optimal control problem for which necessary conditions can be derived using the Pontryagin's Maximum Principle (PMP, as in Section 2.5). In our work, we consider a wider family of systems, namely the LTV systems with a disturbance that satisfies a quadratic constraint at any time. The ellipsoidal method is then a sub-case of our approach.

In the context of discrete-time IQC models, [START_REF] Fry | IQC-based robustness analysis of discrete-time linear time-varying systems[END_REF] chooses the overapproximation that maximizes the input-to-output robust L 2 gain (where the L 2 norm of the signal is computed over a finite time interval). A similar approach method is applied in [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF] for the case of continuous-time system with IQC systems. In both of these works, since the Riccati equation (DRE in the continuoustime case, difference Riccati equation for the discrete-time case) is expressed as an LMI constraint, SDP solvers can then be used to find a solution minimizing a given criterion by using an SDP solver.

Conclusion

Conclusion

This chapter applied the framework introduced in Chapter 1 for the specific case of disturbances with a point-wise quadratic constraint. Since the set of time-varying conic overapproximations is infinite, we proposed two methods to find one "tight" overapproximation of the reachable set. In the first method, we compute a touching time-varying conic overapproximation, this overapproximation touches the reachable set along system trajectories (so-called touching trajectories). In the second method, we find a time-varying paraboloid that minimizes a given criterion (e.g. its volume at a given time). We derive necessary conditions by using the Pontryagin's Maximum Principle (PMP) and proposed a continuation method to solve its associated boundary value problem.

Future Works In this chapter, we only computed overapproximations minimizing their volume at a given time. An interesting extension of this chapter would be to compute be the time-varying conic minimizing its volume at each time over the interval of integration. To address this problem, it is possible to parametrize optimal overapproximations P(t, t * ) with two time indexes t and t * . t * corresponds to the time where the time-varying conic t → P(t, t * ) minimizes its volume, t corresponds to the regular time index. In such a situation, the overapproximation t * → P(t * , t * ) is a time-varying conic with a minimal volume at each time instant t * . Another problem of interest would be to adapt the work of [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF], where optimal multipliers are searched by solving a Differential Linear Matrix Inequality (DLMI), to our optimization problem in order to find the optimal positive multiplier µ.

Set-based Cosimulation

Paul Rousse

Appendices of Chapter 2

2.A Derivatives of the cost functions

We hereby detail the computations of the first derivative Vol P and TrSq P (resp.) of Vol, defined in (1.10), and TrSq, defined in (1.11) (resp.). Then, the Fréchet derivatives of Vol P and TrSq P at P ∈ S (n+1)×(n+1) in a given direction S ∈ S (n+1)×(n+1) .

Derivative of Vol The derivative of Vol at P in the direction δP ∈ S (n+1)×(n+1) is Vol P (P )

• δ P = det(-E) -1 (δ g -2δ ⊤ f Ef -f ⊤ δ E f ) -det(-E) -1 (g -f ⊤ Ef )trace(-E -1 δ E ).
Vol P (P ) • δ P can be expressed with the inner product Vol P (P )

• δ P = trace det(-E) -1 -f f ⊤ + (g -f ⊤ Ef )E -1 -Ef -(Ef ) ⊤ 1 δ P
Therefore, we can express the differential Vol P (P ) of Vol at P in its matrix form (we hereby use an abuse of notation)

Vol P (P ) = det(-E) -1 -f f ⊤ + (g -f ⊤ Ef )E -1 -Ef -(Ef ) ⊤ 1 (2.22)
We compute the Fréchet derivative ∆ P Vol P • S of Vol P in the direction S ∈ R (n+1)×(n+1) . Let the block decomposition of

S S = T u u ⊤ v
To compute ∆ P Vol P • S, we first the following functions

d(P ) = det(-E) -1 r(P ) = g -f ⊤ Ef V P (P ) = -f f ⊤ + r(P )E -1 -Ef -(Ef ) ⊤ 1
then Vol P = d(P )V P (P ). Fréchet derivatives of d, r and V P in the direction S are Derivative of TrSq similarly than for V , we can derive the differential TrSq P of TrSq at P (in its matrix representation)

∆ P d(P ) • S = -det(-E) -1 trace(E -1 T ) ∆ P r(P ) • S = v -u ⊤ Ef -f ⊤ T f -f ⊤ Eu ∆ P V P (P ) • S = -uf ⊤ -f u ⊤ + ∆r(P ) • S E -1 -r(P )E -1 T E -1 -T f -Eu -(T f + Eu) ⊤ 0 Paul Rousse Set-based Cosimulation
TrSq P (P ) = trace(E -1 ) f f ⊤ Ef (Ef ) ⊤ -1 + (g -f ⊤ Ef ) E -2 0 0 0 (2.23)
Similarly, we can compute the Fréchet derivative of TrSq P at P in the direction S.

TrSq =trace(-E -1 T E -1 ) f f ⊤ Ef (Ef ) ⊤ -1 + trace(E -1 ) uf ⊤ + f u ⊤ , T f + Eu (T f + Eu) ⊤ 0 -∆r(P ) • S E -2 0 0 0 -(g -f ⊤ Ef ) -E -1 T E -2 -E -2 T E -1 0 0 0 .

2.B Variations of the state and co-state

Hereby, we explicit the computation of the variation (R, S) of the state and co-state (P, Q) with respect to the initial state (0, S 0 ). R and S are solutions of

Ṙ + A ⊤ Q (t)R + RA Q (t) + V R (t) + trace(ST S + RT R )H R = 0 Ṡ -A Q (t)S -SA ⊤ Q (t) + V S (t) + trace(ST S + RT R )H S = 0 (2.24) with initial conditions R(0) = 0 S(0) = S 0
For the following computations, all signals are time-dependent. We omit its notation to simplify the following computations.

Let k = a b , µ = √ k, then dµ = da a µ-db b µ, where a = -trace(QM x1\w ) b = trace(QP B 1 M -1 w B ⊤ 1 P ) da = ∆a • (R, S) = -trace(SM x1\w ) db = ∆b • (R, S) = trace(SP B 1 M -1 w B ⊤ 1 P + R B 1 M -1 w B ⊤ 1 P Q + QP B 1 M -1 w B ⊤ 1 R) T S = µ M x1\w trace(QM x1\w ) - µ P B 1 M -1 w B ⊤ 1 P trace(QP B 1 M -1 w B ⊤ 1 P ) T R = - µ (B 1 M -1 w B ⊤ 1 P Q + QP B 1 M -1 w B ⊤ 1 ) trace(QP B 1 M -1 w B ⊤ 1 P ) Set-based Cosimulation Paul Rousse Finally, dµ = ∆b • (R, S) = trace(S T S + R T R ). H R = -M x1\w -µ -2 P B 1 M -1 w B ⊤ 1 P H S = -µ -2 (P B 1 M -1 w B ⊤ 1 Q + QB 1 M -1 w B ⊤ 1 P ) V S = ∆ P L P • R + µ -1 RB 1 M w -1B ⊤ 1 Q + µ -1 QB 1 M w -1B ⊤ 1 R

2.C Contracting property in the centered case

We show that in the specific case of centered stable LTI system (see the definition on the paragraph below), we can find an overapproximation that is contracting as t goes to ∞, i.e. P(t) → {0}. The following chapter shows that this result is not available in the case of IQC systems. Let S a centered stable LTI QC system be a QC system as in Definition 2.1 with:

1. time-invariant A, B, C and M matrices (the linear time-invariant system part);

2. with a centered set of initial states (i.e. with f 0 = 0 in the block decomposition of P 0 ); 3. with a centered set of disturbances (i.e. M x1,w (•) = 0); 4. with a null input signal, i.e. u(•) = 0;

5. such that S is exponentially stable.

These assumptions are restrictive but will simplify the following proofs. When this specific case, DRE (1.21) is

0 = Ė + EA + A ⊤ E -µ(t)M x + µ(t) -1 EBM -1 w B ⊤ E. (2.25)
Assumptions 1 and 5 let us express a negative definite solution Ē(•) to the DRE (2.25) with initial condition Ē(0) = Ē0 ≺ 0 such that Ē-1 converges to 0 when t → ∞. Assumptions 2, 3 and 4 are chosen such that any overapproximation P is centered at any time, when they holds, f (•) = 0 and g(•) = g 0 . Since Ē(•) ≤ 0, P(•) is an ellipsoidal set at any time. The radius of P is √ Ē-1 , and if Ē(t) -1 → 0 when t → ∞, then P → {0}. The following paragraph will demonstrate this result. First, we express such a Ē and show that Ē(t) -1 → 0 when t → ∞. Then we show that Ē is an upperbound for solutions of (2.25) with any initial conditions E 0 ≺ 0.

Paul Rousse

Set-based Cosimulation

2.C. Contracting property in the centered case

Let Ē(t) = exp λt Ē0 and µ(t) = μ(t) = κ exp λt with λ, κ > 0. Ē satisfies (2.25) if and only if Ē0 is solution to the following Continuous Algebraic Riccati Equation 0

= λ Ē0 + Ē0 A + A ⊤ E -κM x + κ -1 Ē0 BM -1 w B ⊤ Ē0 .
(2.26)

Since the system is exponentially stable, there is a Ē0 ≺ 0 and κ > 0 such that (2.26) is solved for λ = 0. By continuity, there is a Ē0 ≺ 0, a κ > 0 and a λ > 0 such that (2.26) holds.

It will appear that Ē(t) is an upperbound for any solution to (2.25). Let ∆(t) = Ē(t) -E(t), we aim at proving that ∆(t) ≻ 0. Using (2.25) and the definition of Ē, it holds

∆ = -A ⊤ ∆ -∆A -μ(t) -1 Ē(t)BM -1 w B ⊤ Ē(t) + μ(t) -1 E(t)BM -1 w B ⊤ E(t).
Using the following identity

ĒR Ē -ERE = ∆R∆ + ER∆ + ∆RE where R = -μ(t) -1 BM -1 w B ⊤ , R(t) ≻ 0, we have ∆ + A R (t) ⊤ ∆ + ∆A R (t) = ∆(t)R(t)∆(t)
By integration, it holds

Φ ⊤ R (T )∆(T )Φ R (T ) = ∆(0) + T 0 Φ ⊤ (t)∆(t)R(t)∆(t)Φ ⊤ (t)dt
where Φ R is the transition matrix associated with A R . Therefore, when ∆(0

) ≻ 0, since R ≻ 0, it holds ∆ ≻ 0. By definition, Ē-1 → 0 when t → ∞. Since Ē(t) ≻ E(t), it holds 0 ≺ E(t) -1 ≺ Ē(t) -1
therefore, E(t) -1 → 0. By consequence, P(t) → {0} when t → ∞.

Set-based Cosimulation Paul Rousse

Chapter 3

Integral Quadratic Constraints Chapter 1 introduced the general framework of a set-based simulation method for linear systems subject to bounded disturbances. This chapter applies it for a specific subclass of systems: LTV system subjects to a disturbance that satisfies an Integral Quadratic Constraint (IQC) on the state and input signal. IQC models are a classical tool of robust control theory (see e.g. [Megretski, 2010, Megretski and[START_REF] Megretski | [END_REF]). They can model infinite-dimensional states, nonlinear dynamics, delays, rate limiters, uncertain systems (see [Helmersson, 1999], [START_REF] Rantzer | Analysis of Rate Limiters Using Integral Quadratic Constraints[END_REF]], [START_REF] Peaucelle | Integral Quadratic Separators for performance analysis[END_REF] and [START_REF] Ariba | Stability analysis of time-delay systems via Bessel inequality: A quadratic separation approach[END_REF]).

In Section 3.1, we first define IQC systems. In Section 3.2, we motivate the use of IQC systems as a modeling tool. In Section 3.3, we apply Theorem 1.1 of Chapter 1 to define overapproximations of the reachable set of an IQC system. Contrary to overapproximations defined in Chapter 2 that are supported by the reachable set at any time, it is not possible to satisfy such property for an IQC system. The next sections are then dedicated to the investigation of another approach in order to find tight overapproximations to the reachable set. In Section 3.4, we define an extended system of the original IQC system. In Section 3.5, we show that the reachable set of this extended system can be overapproximated with time-varying paraboloids and we identify touching trajectories and their associated overapproximations supported by the reachable set. In Section 3.6, we define the intersection of all the supporting timevarying paraboloids, Theorem 3.1 states that the reachable set of the extended system is exactly characterized by this intersection of overapproximations. A projection of this overapproximation is the exact reachable set of the initial IQC system. In Section 3.7, we detail the practical implementation of the overapproximation of the reachable set of the IQC system. In Section 3.8, we present some examples. In Section 3.9 and Section 3.10, we discuss and conclude this chapter.

Definition of the system

In this chapter, the system of interest is an LTV system subjects to a disturbance that satisfies an Integral Quadratic Constraint (IQC) with the state and input signal. This IQC constraint can be seen as an energetic relationship on the disturbance and the state/input signals (see Remark 3.1). This system is a subclass of the system described in Definition 1.1.

Definition 3.1. IQC system

Let the system S be defined by

S :      ẋ(t) = A(t)x(t) + B(t)w(t) + C(t)u(t) t 0 x(τ ) u(τ ) w(τ ) ⊤ M (τ ) x(τ ) u(τ ) w(τ ) dτ ≥ 0 for every t ≥ 0 (3.1)
where M is a quadratic form negative in the disturbance dimension, i.e. such

Paul Rousse

Set-based Cosimulation n+m+p) .

(t) =   M x (t) M x,u (t) M x,w (t) M ⊤ x,u (t) M u (t) M u,w (t) M ⊤ x,w (t) M ⊤ w,u (t) M w (t)   ∈ S (n+m+p)×(
Since for two system trajectories (x 1 , w 1 , u 1 ) and (x 2 , w 2 , u 2 ) of S, the trajectory (x 1 , w 1 , u 1 ) + (x 2 , w 2 , u 2 ) does not necessarily satisfy the IQC constraint, i.e. the system S is not linear. However, (α x 1 , α w 1 , α u 1 ), for α ∈ R, is a system trajectory, and therefore the set of trajectories of S is conic.

Remark 3.1. Energetic constraint

The IQC constraint in (3.1) can be reformulated as follows

t 0 (w(τ ) -w c (τ )) ⊤ R (w(τ ) -w c (τ )) dτ ≤ t 0 x(τ ) u(τ ) ⊤ Q x(τ ) u(τ ) dτ (3.2)
where R = -M w ≻ 0,

Q = M x M xu M ⊤ xu M u - M xw M uw ⊤ M -1 w M xw M uw
and

w c = M -1 w M xw M uw .
Then, (3.2) can be reformulated in terms of norm 2 constraints

w -w c , R(w -w c ) t ≤ x q0 + [ x u ] , Q [ x u ] t where y, y t = t 0 y(τ ) ⊤ y(τ )dτ. Let v = √ R(w -w c ), the term v, v 1/2 t
corresponds to the energy of v over [0, t]. The IQC constraint (3.1) can therefore be understood as the energetic constraint

v, v 1/2 t ≤ (x q0 + [ x u ] , Q [ x u ] t ) 1/2 .

Set-based Cosimulation

Paul Rousse The second trajectory (dashed red line) is delayed, with a delay of τ = 4, compared to the first trajectory (plain blue line). This is true for any delay τ . Therefore, the state of such a system can be steered away from its equilibrium position at any time in the future.

Motivation and examples of IQC systems

IQC models are widely used in robust control theory to assert the stability of dynamical systems. Usually, they model nonlinearities or uncertainty within the model. They have strong connections with fundamental system theory and physical concepts (energy and passivity, see Remark 3.1). IQC models can be used to model physical quantities such as the error of an observation (as in [START_REF] Savkin And Petersen ; Savkin | Robust state estimation for uncertain systems with averaged integral quadratic constraints[END_REF]). They also arise when only the H ∞ gain of a feedback loop is available (see [Megretski, 2010]).

Contrary to examples of QC systems in Section 2.2.2, we could not find a simple expression of the reachable set. Hereby, we represent a few trajectories for two systems of one state dimension. The first system corresponds to a linear stable LTI system disturbed by a unit energy disturbance. The second system corresponds to the same stable LTI system with a feedback loop of unit H ∞ gain.

Unit energy exogenous input: the system described by (3.3) corresponds to a linear stable system disturbed by a unit energy noise.

         ẋ = -2x + w t 0 w(τ ) 2 dτ ≤ 1 x(0) = 0 (3.3)

Paul Rousse

Set-based Cosimulation

Overapproximation of the reachable set with time-varying conics

The reachable set of such a system is known to be bounded (see [Boyd et al., 1994, Chapter 6.1.1]). For a trajectory (x, w) of the system (3.3), let (x τ , w τ ) be the translated signal, i.e. (x τ (t + τ ), w τ (t + τ )) = (x(t), w(t)), for any t > 0, (x τ (t), w τ (t)) = (0, 0) elsewhere (see Figure 3.1a). Since the system is time-invariant, the translated signal (x τ , w τ ) is as well a trajectory of (3.3). The state can be steered away from its stable equilibrium at every time τ > 0. Therefore, we expect the reachable set of such an IQC system to be not contractive when t → ∞.

Unit H ∞ gain feedback loop: the system described by (3.4) corresponds to a linear stable system with a unit H ∞ feedback loop.

         ẋ = -2x + w t 0 w(τ ) 2 ≤ t 0 x(τ ) 2 x(0) = 1 (3.4)
Contrary to the system described by (3.3) any translated signal of a system trajectory (x, w) is not necessarily a trajectory of (3.4). However, the same observation applies to the reachable set. For any trajectory, since x(0) = 0, for a null w over a given interval, the system accumulates energy along time. This energy can then be used at any time in the future (see Figure 3.1b). As for the Unit energy case, we expect the reachable set of such an IQC system to be not contractive when t → ∞.

Overapproximation of the reachable set with time-varying conics

Theorem 1.1 of Chapter 1 derives a time-varying conic overapproximation to the reachable set of a system subject to a bounded disturbance. The coefficient of the time-varying conic is the solution of a differential equation parametrized by an element of the dual space D * of the disturbance set D. In this section, we apply Theorem 1.1 for the IQC system of Definition 3.1. To do so, we first identify the dual space D * .

The dual space D * is characterized by for all y ∈ D. Hereby, a given signal y = (x, u, w) belongs to the disturbance set D iff F (t) = t 0 f (τ )dτ is positive for every t ≥ 0, where f (τ ) = y(τ ) ⊤ M (τ )y(τ ). This condition is strictly equivalent to f, µ t ≥ 0 for every t ≥ 0 where

µ ∈ D * iff ∞ 0 y(τ ) ⊤ M (τ )y(τ )µ(τ )dτ ≥ 0 Set-based Cosimulation Paul Rousse
µ t (τ ) = 1 if τ ≤ t 0 otherwise.
Therefore, µ t ∈ D * for any t ≥ 0. Since the weighted sum of µ t and µ t ′ belongs to D * , every function ν defined by ν(t) = ∞ 0 κ(τ )µ τ (t)dτ , where κ is a positive function over R + , belongs to D * as well (see Figure 3.2). To this respect, we can choose D * as the set of functions from R + to R that are positive and decreasing over R + .

Proposition 3.1. Positive integral duality

Any square-integrable measurable signal

f ∈ L 2 (R + ; R) is of positive integral over any interval [0, t], t > 0, iff f, µ ≥ 0
for any square-integrable measurable µ positive and decreasing over R + .

Proof. Depending on the values of λ and µ, the overapproximation is bounded over the interval of integration (P s in blue) or is diverges in finite-time (P u in red).

(⇒) if f (•), µ(•) ≥ 0 over R + , then f, µ ≥ 0. (⇐) if there is a t ≥ 0 s.t. f (t) ≤ 0, then since f is measurable, there is an interval I ⊂ R + centered
The system S of Definition 3.1 is strictly equivalent to the system of Definition 1.1 where D * is the set of square-integrable functions, measurable and positive over R + . Therefore, the following Theorem 1.1 holds.

Corollary 3.1. Application of Theorem 1.1

The set of reachable states R(t; P 0 ) of S of Definition 3.1 is overapproximated at any time instant t ∈ R + and for any conic set of initial of states P 0 ∈ P with coefficient P 0 ∈ S (n+1)×(n+1) , i.e.

R(t; P

0 ) ⊆ P(t) Example 3.1. Let A = -1, B = 1, M = 1 0 0 0 1 0 0 0 -2 , C = 0 and u : [0, ∞[ → 0.
Two overapproximations, as defined in Definition 1.5, are represented in Figure 3.3. The overapproximation P s corresponds to the case where λ(t) = exp(0.3t) and µ(t) = 1, for all t ≥ 0. The overapproximation P u corresponds to the case where λ(t) = 0 and µ(t) = 1, for all t ≥ 0. P s is defined over R + and P u diverges in finite-time.

To assert the quality of overapproximations defined in Corollary 3.1, as for QC systems in Chapter 1, we would like to identify touching trajectories and their associated support conic. To find such trajectories, we should find a valid couple (λ, µ) and an optimal disturbance w such that the level set function associated with the support conic is equal to zero. According to (1.12), it implies that λ is a positive increasing function of R + and where µ is a positive and decreasing function over R + (i.e. µ ∈ D * ), such that y(t) ⊤ M (t)y(t) = 0. It results that µ = µ * where µ * solves (2.9). However, since µ * might be increasing, µ * does not necessarily belong to D * . Therefore, we are only able to identify a subset of touching trajectories. It is not possible to find all conics supported by the reachable set in a similar way than for the QC system in Chapter 2. The next sections show that it is possible to define a time-varying paraboloid overapproximating and supported by the reachable set R of an extended system S. We show that the intersection of these overapproximations leads to the exact characterization of the reachable set of S.

Extended system

In the following sections, we study the system S which extends S, composed of the state signal x solution of

ẋ(τ ) = Ax(τ ) + Bw(τ ) + Cu(τ ) with τ ∈ [0, t] x(0) = x 0 (3.5) with a signal x q ∈ L 2 ([0, t]; R) (corresponding to the IQC constraint) defined for τ ∈ [0, t] by x q (τ ) = x q0 + τ 0 x(s) u(s) w(s) ⊤ M x(s) u(s) w(s) ds, (3.6)
and that satisfies the state constraint

x q (τ ) ≥ 0 for all τ ∈ [0, t].

(3.7)

The constrained dynamical system S(Z 0 , t) is then defined for a given set of initial states Z 0 ⊂ R n × R and a terminal time t > 0

z = (x, x q ) ∈ S(Z 0 , t) ⇔       
x solves (3.5) and x q solves (3.6) with (x 0 , x q0 ) ∈ Z 0 , and x q satisfies (3.7)

(3.8) Let the reachable set of S(Z 0 , t) be R(Z 0 , t) = z(t) z ∈ S(Z 0 , t) . (3.9) Then, R(Z 0 , t) ⊆ Z + where Z + = R n × R + .

Paul Rousse

Set-based Cosimulation 3.4. Extended system

Paraboloids

We overapproximate the reachable set R(Z 0 , t) of S(Z 0 , t) with paraboloids.

Definition 3.2. Paraboloid Given (E, f, g) ∈ S n×n × R n × R, we define the value function p by p : R n × R → R x = (x, x q ) → p(x) -x q ,
and the paraboloid:

P = x = (x, x q ) ∈ R n+1 p( x) ≥ 0 . Definition 3.3.

Scaled paraboloid

For P = Parab(E, f, g) ∈ P and a scaling factor λ > 0, let λP ∈ P be the scaled paraboloid defined by λP = Parab(λE, λf, λg).

Scaled paraboloids satisfy the following:

Proposition 3.2. Overapproximation relationship of a scaled paraboloid Given P ∈ P and λ ≥ 1, it holds P ∩ Z + ⊆ λP ∩ Z + .

Proof. Let h and h ′ (resp.) be the value functions of P = Parab(E, f, g) and λP (resp.) evaluated at (x, x q ) ∈ P. Since (x, x q ) ∈ P, h ≤ 0, i.e. x ⊤ Ex -2f ⊤ x + g ≤ -x q . Then, h ′ = λ(x ⊤ Ex -2f ⊤ x + g) + x q ≤ -(λ -1)x q . Since (x, x q ) ∈ Z + and since λ -1 ≥ 0, we have (λ -1)x q ≥ 0 i.e. h ′ ≤ 0 meaning that (x, x q ) ∈ λP ∩ Z + . ♦

Remark 3.2. On the paraboloid sets

In what follow, a paraboloid P t = Parab(E t , f t , g t ) (see Definition 3.2) is used to overapproximate the reachable set R( P 0 , t) (for a given P 0 ∈ P) at a given t > 0. P t is a paraboloid centered around the ray x q → (x t , x q ) with x t = E -1 t f t with a summit at (x t , x qs ) where x qs = f ⊤ t E -1 t f tg t . Since R( P 0 , t) ⊂ Z + , depending on the parameters (E t , f t , g t ) of P t = Parab(E t , f t , g t ), the overapproximation P t ∩ Z + might describe an empty, unbounded, convex or not convex set. Studying the subset P t ∩ Z * of R n where Z * = R n × {0} (i.e. the section of the paraboloid in the cone R * + = R n × {0}) gives more insight into the shapes of these overapproximations. When x qs > 0 and

E t ≻ 0, P t ∩ Z * is an ellipsoid of quadratic coefficient (f ⊤ t E -1 t f t -g t ) -1 E t . When E t is not sign-definite, P t ∩ Z * is unbounded and not convex. When E t ≺ 0 and x qs ≥ 0, P t ∩ Z * = R n .
Contrary to ellipsoidal sets used in e.g. [START_REF] Savkin And Petersen ; Savkin | Robust state estimation for uncertain systems with averaged integral quadratic constraints[END_REF], Savkin and Petersen, 1996a, Scherer and Veenman, 2018], overapproximations used

Set-based Cosimulation

Paul Rousse in this work are not necessarily bounded and convex (e.g. when E t is signundefined). Since the reachable set is not always convex (see Example 3.3 and Section 3.8.2), it allows us to define tight overapproximations that are in contact with the reachable set even where the surface of the reachable set is not locally convex.

Overapproximation with time-varying paraboloids

In this section, we define time-varying paraboloids that overapproximate the reachable set S. This overapproximation relationship arises from Theorem 1.1 of Chapter 1.

These overapproximations are then used in Section 3.6 to define tighter overapproximations of the reachable set. We prove that the overapproximations P are tight since there are touching trajectories of R( P 0 , t) that both belong to the surface of P(t) and to the surface of R( P 0 , t) for t ∈ I. Finally, the method is presented for a simple example.

To overapproximate the reachable set R, we choose λ and µ such that λ(t) µ(t) = 1, for every t ≥ 0. This constraint is compatible with the necessity for λ to be positive and increasing and µ to be strictly positive and decreasing. Let γ = λ λ -1 , by hypothesis over λ, γ can be any positive function from R + to R + .

Definition 3.4. Time-varying paraboloid

For an initial paraboloid P 0 ∈ P, let the time-varying paraboloid P be defined as

P : I → P t → Parab(P (t))
where the time-varying coefficient P solves

0 = Ṗ + P A 1 + A ⊤ 1 P -γ P -M x1 + B ⊤ 1 P -M w,x1 ⊤ M -1 w B ⊤ 1 P -M w,x1 (3.10)
with the initial condition P (0) = P 0 , where Parab(P 0 ) = λ 0 P 0 . Let T be the function that associates to the initial paraboloid P 0 ∈ P the time-varying paraboloid P. Let T P (P) = T E (E 0 ) and I(P) = [0, T P (P)[ be the interval of definition of P.

Then, by integration of (1.7) along a system trajectory, for the case described in Definition 3.4, it holds d p(t, x) ≥ γ ( p(t, x) + x q ) where x = (x, x q ). Since γ(•) ≥ 0, it holds that p(0, x(0)) ≥ 0 implies that for any t ≥ 0, then p(t, x(t)) ≥ 0.

Paul Rousse

Set-based Cosimulation

Overapproximation with time-varying paraboloids

For γ ∈ L 2,loc (R + ; R + ), λ 0 ≥ 1, let P = T ( P 0 , λ 0 , γ) be the time-varying paraboloid with time-varying parameters defined by (3.10) for initial conditions defined by P 0 . Equation (3.28) in Appendix 3.A proves that when the quadratic coefficient E diverges in finite-time, at t * > 0, the interval of definition of the time-varying paraboloids, which is [0, t * [ can be as well prolongated to the closed interval [0, t * ].

The worst disturbance, given by (1.14), is associated with the maximal variation of the value function.

Proposition 3.3. Variation of the value function

For γ ∈ L 2,loc (R + ; R + ), λ 0 ≥ 1, and P 0 ∈ P, let P = T ( P 0 , λ 0 , γ). For an optimal trajectory z * generated by the disturbance w defined in (1.14) s.t. z * (0) ∈ ∂ P 0 , for any t ≥ 0 it holds

˙ p z * (t) = γ(t)( p z * (t) -x * q (t)). (3.11)
Proof. Direct derivation from (3.10). ♦

For any t ∈ I(P), the solution to the ODE (3.11) is

p z * (t) = (1 -λ 0 )(x * q (0) -p 0 z * )e t 0 γ (r)dr - t 0 γ(s)x * q (s)e t s γ (r)dr ds (3.12)
where p 0 z * is the evaluation at z * (0) of the value function of P 0 . Equations (3.11) and (3.12) provide a convenient way to a) prove the overapproximation relationship; b) identify touching trajectories; c) reject invalid trajectories. a) p z (t) ≤ 0: for a valid system trajectory z, since z(0) ∈ P 0 (i.e. p 0 z ≤ 0) and (3.7) holds, (3.12) ensures that z(t) ∈ P(t) for any t (stated in Corollary 3.2); b) p z (t) = 0: optimal trajectories z * as defined in Proposition 3.3 are touching trajectories when z(0) belongs to the boundary of P 0 and γ(t)x * q (t) = 0, for all t ≥ 0, and (1λ 0 )x * q (t) (stated in Proposition 3.4); c) p z (t) > 0: when z violates the constraint (3.7), there is a t ≥ 0 s.t. x * q (t) < 0, one can choose a γ(t) > 0 such that ˙ p z (t) ≥ 0 eventually leading to p z (t ′ ) > 0, t ′ > t, proving that z(t ′ ) / ∈ R(t; P 0 ) (this will be used in Proposition 3.7).

Intuitively, time-varying paraboloid P is contracting for valid trajectories (i.e. when z(0) ∈ P 0 and (3.7) holds) and expanding for invalid trajectories (i.e. when either z(0) ∈ P 0 or either (3.7) is violated).

Set-based Cosimulation Paul Rousse

Since t 0 γ(s)x * q (s)e t s γ (r)dr might not be equal to 0, trajectories generated by the worst-case disturbance w * do not necessarily stay in contact with the time-varying paraboloid and therefore are not touching trajectories (see Definition 1.3). For this reason, we call optimal trajectories the trajectories generated by w * given by (1.14).

When some conditions hold (see below), an optimal trajectory might be a touching trajectory.

Proposition 3.4. Touching optimal trajectories Let z * be an optimal trajectory of P s.t. z * (0) ∈ ∂ P 0 , if (1λ 0 )x * q (0) = 0, and γ(t)x * q (t) = 0 for any t ≥ 0, then z * is a touching trajectory of P .

Proof. Since, z * (0) ∈ ∂ P 0 , it holds p 0 z * = 0. Using (3.12), we get that p z * (t) = 0 for all t ≥ 0. ♦ Theorem 1.1 can be rewritten in the specific IQC case for the extended system S Corollary 3.2. Overapproximation of the set of reachable states For a set of initial states P 0 , a time-varying multiplier γ ∈ L 2,loc (R + ; R + ) and an initial multiplier λ 0 ≥ 1, let P = T ( P 0 , λ 0 , γ). The reachable set R( P 0 , t) of S( P 0 , t), t > 0, is overapproximated by P(t), i.e.

∀t ∈ I(P), R( P 0 , t) ⊆ P(t) ∩ Z + .

Proof. For any trajectory z of S( P 0 , t), z(0) ∈ P 0 implies that p z * (0) ≤ 0, since (3.7) holds, using (3.12), p z * (t) ≤ 0, i.e. z(t) ∈ P(t). Any trajectory z satisfies the state constraint (3.7) over R + , so z(t) ∈ Z + , for any t ∈ I(P) and therefore

R( P 0 , t) ⊆ Z + . ♦ Example 3.2. Let A = -1, B = 1, M = 1 0 0 0 1 0 0 0 -2 , C = 0 and u : [0, ∞[ → 0. Solutions to IVP (1.21) (that is Ė = -1 2 E 2 + 2E -1)
diverge when E 0 ≺ E -and converge to E + when E 0 ≻ E -(see Figure 3.4) where E -≺ E + are the roots of the equation

-1 2 E 2 + 2E -1 = 0 for E ∈ R, E -= 2 - √ 2 and E + = 2 + √ 2.
Since time-varying paraboloids of Definition 3.4 are defined over the domain of definition of E, for a time-varying paraboloid P with the initial value (E 0 , f 0 , g 0 ), depending on whether E 0 is in the stable region (i.e. E 0 ≻ E -) or the unstable region (i.e. E 0 ≺ E -), the time-varying paraboloid might be defined over R + or a finite horizon only. Figure 3.5 shows the trajectory of the time-varying paraboloid P for E 0 in the stable region. Figure 3.6 shows the trajectory of the paraboloid for E 0 in the unstable region. .5: Time-varying paraboloid (its boundary is the green line) overapproximating the reachable set (the gray shaded area with the red boundary) at different time instants t in {0.00, 0.91, 1.62, 10.00} for an initial maximum energetic level of x q,0 = 0.06. The solution to (1.21) converges to a constant value when t → +∞. The shaded regions are the reachable set R( P 0 , t) (with P 0 = (E 0 , f 0 , g 0 ), E 0 = 0.6, f = 0 and g 0 = -0.06), the thin lines are the boundary of the overapproximation P(t) of Corollary 3.2.

Overapproximation with time-varying paraboloids

Figure 3.6: Time-varying paraboloid (its boundary is the green line) overapproximating the reachable set (the grey shaded area with the red boundary) at different time instants t in {0.00, 0.91, 1.62, 10.00} for an initial maximum energetic level of x q,0 = 0.03. The solution to (1.21) has a finite escape time and diverges at t = 1.68. The shaded regions are the reachable set R( P 0 , t) (with P 0 = (E 0 , f 0 , g 0 ), E 0 = 0.3, f = 0 and g 0 = -0.03), the thin lines are the boundary of the overapproximation P(t) of Corollary 3.2.

Set-based Cosimulation

Paul Rousse

Remark 3.3. Domain of definition

Previous works in reachability analysis of IQC systems (such as [START_REF] Savkin | Recursive State Estimation for Uncertain Systems with an Integral Quadratic Constraint[END_REF][START_REF] Savkin And Petersen ; Savkin | Robust state estimation for uncertain systems with averaged integral quadratic constraints[END_REF], Jönsson, 2002, Seiler et al., 2019]) derive overapproximations as in Section 3.5 but with a scaling function γ(t) = 0, for all t ≥ 0. In such a case, the DRE (3.10) might have a finite escape time depending on the initial set (more precisely, depending on the initial condition E 0 of E) and on the system's parameters. When the initial set belongs to the unstable region, no time-varying paraboloid is defined over R + (even if the reachable set is bounded at any time). In Example 3.2, the overapproximations are derived using a null scaling function γ(t) = 0, their domain of definition that depends on the initial condition E 0 of the associated coefficient of E (see the convergence analysis in Figure 3.4). When E 0 is greater than the unstable equilibrium of DRE (1.21), E is defined over R + , otherwise, E has a finite escape time. Such an unstable equilibrium of DRE (1.21) exists for stable systems, but it does not exist for unstable systems: the DRE (1.21) does not have any equilibrium at all. In this case, solutions E diverge in finite-time for every E 0 ∈ R. Example 3.3 shows that by using a nonzero scaling function γ in the DRE (1.21), it is possible to find a γ large enough such that the overapproximation is defined over the entire time-interval R + .

Thus, our method to overapproximate the reachable set of an IQC system extends the scope of systems that can be studied to unstable systems.

Example 3.3.

Let the system S( P 0 ) defined by parameters:

A = -1, B = 1, B u = 0, M = 1 0 0 0 1 0 0 0 -0.9
and a zero input signal u, and let the set of initial states be the paraboloid P 0 = Parab(E 0 , f 0 , g 0 ) with parameters:

E 0 = 1, f 0 = 0 and g 0 = -0.015.
Let the disturbance signal w = x √ 0.9 , w satisfy the IQC (3.7) since t 0 x 2 (τ ) -0.9w 2 (τ )dτ = 0 for every t ≥ 0. Trajectories associated with such w satisfy ẋ = αx with α = 1 √ 0.9 -1. Since α > 0, every trajectory starting from a non zero initial condition x(0) = x 0 = 0 diverges when t → ∞, and the system is said unstable.

Figure 3.7 shows plots of the reachable set of this system and several overapproximations at different time instants. Overapproximations are derived using Definition 3.4 for different scaling functions and initial scaling factors. These

Paul Rousse

Set-based Cosimulation The time-varying paraboloid P ns is defined over [0, 3] whereas P s is defined over R + . overapproximations have different domains of definitions depending on these scalings. The solution to DRE (3.10) for γ(t) = 0, for all t ≥ 0, and λ 0 = 0 has a finite escape time and diverges at T P (P ns ) = 1.7. The solution to DRE (3.10) for γ(t) = 1, for all t ≥ 0, and λ 0 = 0 is defined over R + .

Their domains of definitions can be studied by conducting a stability analysis of DRE (1.21). Figure 3.9 plots the phase portrait of DRE (1.21) for different constant scaling functions. Figure 3.8 plots the domain of convergence of DRE (1.21) depending on the scaling function γ and the initial condition E 0 . For γ = 0, solutions to the DRE escape in finite-time for every initial value E 0 . When γ(t) = κ, for all t ≥ 0, with κ ≥ 0, for every E 0 , there is a value of κ such that the associated solution E converges to a stable equilibrium of the DRE.

Set-based Cosimulation

Paul Rousse

Exact reachable set

The previous section introduced time-varying paraboloids that overapproximate the reachable set R of S. In this section, for the set of time-varying paraboloid defined in Definition 3.4, at each time instant, we show that the intersection of these paraboloids is an overapproximation of the reachable set (in Section 3.6.2). Then, we prove that when some topological assumption holds about the reachable set, our overapproximation is equal to the reachable set (in Sections 3.6.3, 3.6.4, and 3.6.5).

Overapproximation with an intersection of time-varying paraboloids

In this section, a set of time-varying paraboloids is defined. At a given time, the intersection of the paraboloids gives better overapproximations of the reachable set.

With additional assumptions about the topology of the reachable set, the reachable set is exactly characterized. This approach relies on the use of Corollary 3.2 and preliminary results showing that for any state of the overapproximation, there exists a trajectory in S( P 0 , t), t > 0, leading to this state. Let Π * be defined as follows

Π * = {T ( P 0 , λ 0 , γ) | γ ∈ L 2,loc (R + ; R + ), γ ≥ 0, λ 0 ∈ R, λ 0 ≥ 1}. ( 3
P(t) (3.15)
the intersection of all the defined time-varying paraboloids P of Π * at time t (see Figure 3.10). Since Π * (•) is defined over R + , Π(•) is defined over R + .

Figure 3.10: For a given t ≥ 0, consider three time-varying paraboloids P i ∈ Π * , i = 1, 2, 3. Light color shaded areas are their corresponding parabolic sets P i (t) at t, i = 1, 2, 3. The grey color shaded area is their intersection. By (3.15), Π(t) is a subset of P 1 (t) ∩ P 2 (t) ∩ P 3 (t).

Since Assumption 1.1 is satisfied in the case of IQC systems (indeed, t → exp(-kt) is a decreasing and positive function), therefore Proposition 1.5 holds, and the following holds:

Corollary 3.3.

Domain of definition of the intersection of time-varying paraboloid

When E 0 ≻ 0 (i.e. the set of initial states is bounded), Π is defined over R + .

Proof. By Proposition 1.5. ♦

We now prove that, when some assumptions about the topology of Π hold (Assumption 3.1 and 3.2), we have R( P 0 , t) = Π(t) ∩ Z + , for any t ≥ 0 (Theorem 3.1, Section 3.6.5). To achieve that:

Set-based Cosimulation

Paul Rousse

Figure 3.11: Past trajectory z * constructed to prove that any state within Π(t) ∩ Z + belongs to R( P 0 , t) (in grey color shaded area with red boundary). The trajectory z * connects the endpoint z(t ′ ) to trajectory z(t) which is a touching trajectory of some time-varying paraboloid P ∈ Π * (its boundary is represented with the plain black line).

• we prove the overapproximation relationship R( P 0 , t) ⊆ Π(t) (Section 3.6.2);

• we prove that any state (x, x q ) ∈ Π(t) is reachable from a state (x, x ′ q ) ∈ ∂ Π(t) with x q ≤ x ′ q (Section 3.6.3);

• for a state z t ∈ ∂ Π(t), we find a touching trajectory z * = (x * , x * q ) of Π such that z * (t) = z t . This touching trajectory (x * , x * q ) of Π satisfies the state constraint x q (•) ≥ 0 over [0, t] (Section 3.6.4);

• finally, we conclude that any z t ∈ Π(t) is reachable from P 0 , thus R( P 0 , t) = Π(t) ∩ Z + (Section 3.6.5).

The three last bullet-points characterize a system trajectory z that associates a given state (x, x q ) ∈ Π(t) to an initial state that belongs to the set of initial states. We illustrate this in Figure 3.11.

Overapproximation relationship

Corollary 3.2 states that each time-varying paraboloid of Definition 3.4 is an overapproximation of the reachable set. An intersection of many time-varying paraboloids is as well an overapproximation of the reachable set.

Proposition 3.5. Overapproximation relationship R( P 0 , t) ⊆ Π(t) ∩ Z + for any t ≥ 0.

Paul Rousse

Set-based Cosimulation 3.6. Exact reachable set Figure 3.12: Time-varying paraboloids (its boundary is the green line) overapproximating the reachable set (the grey shaded area with the red boundary) at different time instants t in {0.00, 0.91, 1.62, 10.00} for different initial multipliers. Time-varying multipliers (see Definition 3.4) are equal to the null signal, γ i = 0, and initial multiplier λ 0 i ≥ 1 are respectively equal to 1.0, 1.6, 2.2, 2.7 and 3.3 for i = 0, . . . , 4. The shaded regions are the reachable set R( P 0 , t) (with P 0 = (E 0 , f 0 , g 0 ), E 0 = 0.3, f = 0 and g 0 = -0.03), the thin lines are the boundary of the overapproximation P(t) of Theorem 1.1.

Proof. This is a direct consequence of Corollary 3.2 and (3.15).

♦ Example 3.4.
Example continued from Example 3.2. In the case where the solution to (1.21) does not converge (i.e. E 0 < E -), Figure 3.12 shows several paraboloid trajectories with different initial multipliers. Time-varying multipliers are equal to the 0 function and initial multiplier λ i are greater than 1, P 0 ∩ Z + ⊂ λ i P 0 ∩ Z + . Therefore, each time-varying paraboloid is a valid constraint that bounds R( P 0 , t), t ∈ I( Π * ) (Theorem 1.1). Therefore, R( P 0 , t) ⊆ P * (t) = P 0 (t)∩P

1 (t)∩• • •∩P 4 (t)
where P i = T ( P 0 , λ i , 0), and λ i are resp. equal to 1, 1.6, 2.2, 2.7 and 3.3 for i = 0, . . . , 4. In this case, the overapproximation P * (t) is strictly included in P 0 (t).

Observations in Example 3.4 motivate the use of multiple time-varying paraboloids to get better overapproximations of the reachable set R( P 0 , t), t > 0.

Past trajectory for states in the interior of the overapproximation

Proposition 3.6 shows that the state (x, αx q ) is reachable from the given state (x, x q ) for any given α ∈ [0, 1].

Proposition 3.6. Reachability of states in the interior of the overapproximation For t ≥ 0, if (x, x q ) ∈ R( P 0 , t) then (x, αx q ) ∈ R( P 0 , t) for all α ∈ [0, 1].

Set-based Cosimulation

Paul Rousse

Proof. Let f : t, x → Ax + Bw(t) + Cu(t). Since for any (t,

x) ∈ R + × R n , f (., x) is locally measurable over R + , f (t, .
) is Lipschitz over R n , (3.5) has a unique solution x (see [Schuricht and von der Mosel, 2000], Theorem 1.1) that is time-continuous. Therefore, for a trajectory (x, x q ) ∈ S( P 0 , T ), T > 0, x is time-continuous.

For ǫ > 0, let w ∈ L 2 ([0, t + ǫ]; R m ), s.t. w ⊤ (s)M w w(s) = -(1 -α)x q (t) 1 ǫ when s ∈ [t, t + ǫ]. Then t+ǫ t w ⊤ (s)M w w(s)ds → -(1 -α)x q (t) when ǫ → 0. Using the Cauchy-Schwartz inequality t+ǫ t (-M w ) 1 2 w(s)ds ≤ √ ǫ t+ǫ t -w T (s)M w w(s)ds
and the time-continuity of x, the quantity

t+ǫ t x(s) u(s) 0 ⊤ M x(s) u(s) w(s)
ds → 0 when ǫ → 0. By integration, x q (t + ǫ) → αx q (t) when ǫ → 0. Since x is timecontinuous, x(t + ǫ) → x(t) when ǫ → 0. Since u is bounded at any time (u ∈ L ∞ (R + ; R p )) and since x is continuous, w is bounded over [t, t + ǫ]. Therefore, u, x and w are bounded over [t, t + ǫ] x q is continuous over [t, t + ǫ]. Then, there exists a t ′ ∈ [t, t + ǫ] such that x q (τ ) ≥ αx q (t) ≥ 0 for all τ ∈ [t, t ′ ] and x q (t ′ ) → αx q (t) when ǫ → 0. Therefore, the constraint x q (•) ≥ 0 is satisfied over [t, t ′ ] and the trajectory (x, x q ) is a valid trajectory of S( P 0 , t ′ ) for all t ≤ t ′ . ♦

Past trajectory for states on the boundary of the overapproximation

In this section, touching trajectories of Π are identified. We show that all these touching trajectories satisfy the state constraint (3.7).

The value function h of a time-varying paraboloid P ∈ Π * can be approximated at the first-order along a touching trajectory z * of another time-varying paraboloid P ∈ Π * when their time-varying multiplier γ and γ and initial multiplier λ0 and λ 0 are close to each other. In this part, we compute this first-order approximation when γ = γ + δ and λ0 = λ 0 + δ 0 for small variations δ ∈ L 2 (R + ; R) and δ 0 ∈ R (i.e. when δ + |δ 0 | tends to 0).

Paul Rousse

Set-based Cosimulation

Exact reachable set

To show that the reachable set R( P 0 , t), t > 0, is exactly described by Π(t), we use a proof by contradiction. For any optimal trajectory z * of P ∈ Π s.t. z * is violating the constraint (3.7), there is a P ∈ Π * such that the endpoint z * (t) does not belong to P (t) and therefore to Π(t). To do so, we will study the value function of a time-varying paraboloid P for touching trajectories of P . where

β(t) = δ 0 λ 0 x * q (0) + t 0 δ(s)ψ(s)e t s (γ (r)+δ(r))dr ds and ψ(s) = x * q (s) - t s γ(τ )x * q (τ )e s τ γ (r)dr dτ.
Proof. Let (E , f , g) and ( Ẽ , f , g) (resp.) be parameters of P and P (resp.), and

ν = (E -Ẽ )x * -(f -f ).
Using (3.5), (1.14), and (3.10), ν satisfies the linear time-varying differential equation

ν(τ ) = A ν (τ )ν(τ ) -δ(τ )n(τ ). (3.16) with n = E x * -f and A ν (τ ) = -A ⊤ + M xw M -1 w B ⊤ + Ẽ (τ )BM -1 w B ⊤ + γ(τ )I
. By Definition 3.4, initial values of P and P satisfy

1 λ 0 P (0) = 1 λ0 P (0) = P 0 , therefore ν(0) = δ 0 (E 0 x * (0) -f 0 ) where (E 0 , f 0 , g 0 ) = P 0 .
ν(t) -ν(0) ≤ t 0 L(τ ) ν(τ ) dτ + t 0 |δ(τ )| n(τ ) dτ. Since l → l 0 |δ(τ )| n(τ )
dτ is a non-decreasing function over [0, t], by applying the Grönwall inequality, we get

ν(t) ≤   ν(0) + t 0 |δ(τ )| n(τ ) dτ   e t 0 L(τ )dτ . (3.18) Let q(t) = h t (z * (t)) -ht (z * (t)). (3.19)
z * is an optimal trajectory of P s.t. z * (t) ∈ ∂P (t), therefore, h t (z * (t)) = 0, therefore, using (3.11)

h τ (z * (τ )) = t τ γ(s)x * q (s)e
(τ ) = -ν ⊤ (τ )BM -1 w B ⊤ ν(τ ) + γ(τ )h τ (z * (τ )) -γ(τ ) hτ (z * (τ )) -γ(τ )x * q (τ ) + γ(τ )x * q (τ ).
Using (3.19) and (3.20):

q(τ ) = -ν ⊤ (τ )BM -1 w B ⊤ ν(τ ) + δ(τ )ψ(τ ) + (γ(τ ) + δ(τ ))q(τ ) (3.21)
where

ψ(τ ) = x * q (τ ) - t τ γ(s)x * q (s)e s τ γ (r)dr ds

Paul Rousse

Set-based Cosimulation 3.6. Exact reachable set with the initial condition q(0) = h 0 (z * (0))-h0 (z * (0)). Since z * is a touching trajectory of P , it holds h 0 (z * (0)) = 0, therefore

x * q (0) = -λ 0 (x * (0) ⊤ E 0 x * (0) -2f ⊤ 0 x * (0) + g 0 ).
Therefore, h0 (z * (0)) satisfies

h0 (z * (0)) = - δ 0 λ 0 x * q (0)
and q(0) = δ 0 λ 0

x * q (0). Since t belongs to the open set of I( P ), the optimal trajectory z * and ν are defined and continuous over [0, t]. Moreover, since γ and δ are measurable over [0, t], the solution to the linear time-varying equation (3.21) exists over [0, t] and is

q(τ ) = -h0 (z * (0)) + τ 0 -ν(s) ⊤ BM -1 w B ⊤ ν(s) + δ(s)ψ(s) e t s (γ (r)+δ(r))dr ds.
Then, using (3.18)

q(t) - δ 0 λ 0 x * q (0) - t 0 δ(s)ψ(s)e
ψ(s) = x * q (s) - t s γ(τ )x * q (τ )e s τ γ (r)dr dτ.
Proof. This is a direct consequence of Proposition 3.7 and of the property:

(|a -b| ≤ c) ∧ (c < |b|) ⇒ sign(a) = sign(b) for a, b, c ∈ R. ♦
Provided the existence of a (δ, δ 0 ) ∈ L 2 (R + ; R) × R such that γ + δ ≥ 0 and λ 0 + δ 0 ≥ 1, the first-order approximation of the value function of P = T ( P 0 , λ 0 +δ 0 , γ +δ) gives a way to identify time-varying paraboloids P that belongs to Π * such that an invalid trajectory with an end state z t ∈ ∂P(t) (meaning with an initial state outside of the initial set P 0 or a trajectory violating the constraint) does not belong to P (t) and therefore, does not belong to Π(t). Proposition 3.9 states that the touching trajectories of Π satisfy the state constraint (3.7). Proposition 3.9 is proven by choosing a valid trajectory candidate. If this trajectory violates the state constraint (3.7), then Proposition 3.8 provides a proof that this trajectory does not belong to the overapproximation Π.

Proposition 3.9. Valid touching trajectory

For P ∈ Π, if z t ∈ ∂ Π(t) and z t ∈ ∂P (t) for t in the open set of I(P ), then the optimal trajectory z * of P such that z * (t) = z t is a valid touching trajectory of P and Π.

Proof. Let ψ : R + → R be defined for s ≥ 0 by

ψ(s) = x * q (s) - t s γ(τ )x * q (τ )e s τ γ (r)dr dτ.
Let τ ∈ [0, t] and I = [τ, t].

• Case 1, x * q (0) < 0: with δ 0 > 0, using Proposition 3.8, z * (t) / ∈ P (t) where P ∈ Π * since δ 0 + λ 0 ≥ 1, so z * (t) / ∈ Π(t).

Paul Rousse

Set-based Cosimulation

Exact reachable set

• Case 2, ψ(•) < 0 over I: any δ(•) ≥ 0 over I and δ(•) = 0 elsewhere such that t 0 δ(s)ψ(s)ds = 0 and for δ 0 = 0, using Proposition 3.8, z * (t) / ∈ P (t) where P ∈ Π * since γ + δ ≥ 0, so z * (t) / ∈ Π(t). • Case 3, ψ(•) > 0 over the open of I and there is a l ∈ I, s.t. s∈[τ,l] γ(s)x * q (s)ds = 0: since γ ≥ 0, there exists a δ ≤ 0 such that γ + δ ≥ 0 and for δ 0 = 0, using Proposition 3.8, z * (t) / ∈ Π(t). • Case 4, ψ(•) = 0 over I: since x * q is continuous over R + , and since γ is locally square-integrable, ψ(t) = 0 ⇒ x * q (t) = 0, therefore x * q (•) = 0 over I. Consequently h τ (z * (τ )) = 0 for τ ∈ I. Cases 1 to 4 show that for z * (t) ∈ ∂ Π(t):

• either ∀l ∈ I, l τ γ(τ )x * q (τ )dτ = 0 and ψ(•) = x * q (s) > 0; • nor x * q (l) = 0 for l ∈ I. Let a partition [0, t] = i∈N I i be such that over each open interval I i , ψ(•) ⊲⊳ i 0 with ⊲⊳ i ∈ {<, >, =}. We deduce that for every s ∈ [0, t], x * q (s) ≥ 0 and I γ(τ )x * q (τ )dτ = 0. z * is a valid trajectory, i.e. the constraint (3.7) is satisfied. By (3.11), z * is a touching trajectory of P. Moreover, since z * (0) ∈ P 0 , z * is as well a touching trajectory of Π.

♦

Since Π(t) is an intersection of closed sets, Π(t) is closed as well. In the general case, for an infinite intersection Y * = i∈N Y i of closed sets Y i , i ∈ N, any boundary point y ∈ ∂Y * does not necessarily belong to the boundary of any Y i , i ∈ N (e.g.

ǫ∈]1,2] [-ǫ, ǫ] = [-1, 1], but there is no ǫ ∈]1, 2] such that 1 ∈ ∂[-ǫ, ǫ]).
The following assumption states that for every state on the boundary of the overapproximation Π(t), t > 0, there exists a time-varying paraboloid P such that this state belongs as well to the boundary of the P (t).

Assumption 3.1. Closedness of Π(t).

For every z t ∈ ∂ Π(t), there is a P ∈ Π * such that z t ∈ ∂P(t).

This assumption is not a strong one and it is satisfied for simpler cases (see [Rousse et al., 2019, Property 11]).

In Proposition 3.9, the existence of γ and λ 0 is conditioned by t belonging to the open domain I( P ); to ensure this, E(•) is assumed to be bounded over [0, T ] (by considering the case where t is in the open set of I(P )). In the general case, the boundedness of E(•) is not granted (see the unstable case in Example 3.2 and Figure 3.6). Assumption 3.2 states that for any state on the boundary of the overapproximation Π(t), t > 0, there is a neighbor state on the boundary of P (t) where P is a time-varying paraboloid of Π * not diverging at t (i.e. t belongs to the interior of T P (P )).

Set-based Cosimulation

Paul Rousse Assumption 3.2. Unbounded time-varying paraboloids For t > 0, for all ǫ > 0, for all z t ∈ ∂ Π(t) such that z t ∈ ∂P (t), P ∈ Π * with P (t) unbounded, there is a zt that belongs to the boundary of P (t), zt ∈ ∂ P (t), such that z tzt < ǫ.

Lemma 1 shows that all state z t ∈ ∂ Π(t) (with t ∈ I( Π * ) given) is the terminal state of a touching trajectory z * of Π with an initial state z * (0) ∈ ∂ Π(0). Proof. Let z t / ∈ R( P 0 , t) such that for any ǫ > 0, there is a touching trajectory z of P , P finite, with z(t) ∈ R( P 0 , t) and z(t)z t < ǫ. For x q (t) > 0, we can define the optimal trajectory z * with z * (t) = z t . For all τ ∈ [0, t], P (τ ) is finite. Then, Proposition 3.9 can be used over [0, τ ]. Therefore, if z t ∈ ∂P such that P diverges at t, it holds

z t ∈ ∂ Π(t) ⇔ z t ∈ R( P 0 , t)
For states not belonging to a diverging time-varying paraboloid, the property is a direct consequence of Assumption 3.1, Proposition 3.9. ♦ Lemma 1 shows that any point on the boundary belongs to the reachable set since, for any given terminal state, we found a past trajectory (the touching trajectory) that satisfies the constraint (3.7) and with an initial condition in the set of initial states.

Exact reachable set

We now state the main result of the chapter.

Theorem 3.1. Exact reachability

When Assumption 3.1 and Assumption 3.2 hold, the reachable set R( P 0 , t) of system S( P 0 , t) (defined in Section 3.4) is equal to the set Π defined in (3.15), namely Π(t) = R( P 0 , t) for all t ≥ 0.

Proof. Corollary 3.2 states that R( P 0 , t) ⊆ Π(t). By Proposition 3.6, for z t ∈ Π(t), we can construct a trajectory z such that z

(t) = z t , z(t -) = z * t ∈ ∂ Π(t) (Proposition 3.6). Since z * t ∈ ∂ Π(t)
, using Lemma 1, there exists a trajectory z such that z(t -) = z * t and z is a touching trajectory of Π on [0, t[. Since z is a touching trajectory of Π, z(0) ∈ ∂ Π(0) with Π(0) = P 0 = R(0). By Proposition 3.8, the trajectory z is valid (i.e. satisfies the energy constraint (3.7)) z t ∈ R( P 0 , t). ♦

Paul Rousse

Set-based Cosimulation 3.7. Implementation

Implementation

In this part, we discuss the practical implementation of the reachable sets overapproximation using Theorem 3.1. To do so, we compute a subset Π * of Π * , Π * ⊆ Π * .

(3.24) Π * corresponds to the time-varying paraboloid set generated by a finite subset of time-varying multiplier and initial multiplier. Then, the intersection of each timevarying paraboloid evaluated at a given t > 0 is an overapproximation of the reachable set R( P 0 , t). Finally, the DRE numerical integration is detailed for the case of non-negative solutions to the DRE. We propose an algorithm (Algorithm 3.1) that computes Π * . Its implementation in Matlab is available online [Rousse, 2019].

Subset of time-varying multipliers and initial multipliers In this work, the time-horizon is partitioned into intervals of width T c > 0 and the construction of the over-approximation is restarted for each interval with an appropriate initial multiplier while the time-varying multiplier is set to the zero function γ(•) = 0. A multiplier is then described by a sequence of initial multipliers

{λ k } k∈N , λ k ≥ 1, k ∈ N.
In the ideal case, the multipliers (time-varying and initial) would be chosen such that the following property is satisfied

∃ǫ > 0, ∀τ ∈ [t, t + ǫ], x * q (τ ) ≥ 0 (3.25)
where (x * , x * q ) corresponds to the touching trajectory associated with the time-varying multiplier γ and initial multiplier λ 0 and such that (x * (t), x * q (t)) = (x, x q ). In practice, since there might be an infinite number of states (x, x q ) satisfying ẋq ≥ 0, only a finite number of states are checked. These states are chosen as projections of a given point in given directions over Z * ∩ ∂Π where Z * = R n × {0}. These points are then used to evaluate a range of initial multipliers λ 0 to enforce ẋ * q (kT c ) ≥ 0. γ is not used.

Paraboloid numerical integration Consider two paraboloids P = T ( P 0 , 1, γ), P = T ( P 0 , 1, γ). If γ(.) = γ(.) over an interval [0, t i ], t i > 0, then P (.) = P (.) over [0, t i ]. Let t i ≥ 0 correspond to the maximal time instant where there is P ∈ Π * such that P | [0,t i ] = P | [0,t i ] (i.e. such that the restriction of P on [0, t i ] is equal to the one of P on the same interval). And let t f ≥ 0 correspond either to the integration horizon T > 0, or to the maximal of the interval of definition of P . For implementation purposes, each time-varying paraboloid is defined over the interval

[t i , t f ] ⊆ [0, T ].
The solution to P = T ( P 0 , 1, γ) is then described by parameters (E, f, r) with

(E(t), f (t), g(t)) = (E k (t), f k (t), g k (t)) for each t ∈ [kT c , (k + 1)T c ], k ∈ N.
Set-based Cosimulation Paul Rousse The increasing and decreasing parts of E are respectively represented by the terms L and K. Our Chandrasekhar inspired method performs better since the square root term L and K are much smaller than E, and they produce less numerical errors.

For f and g, integration of the ODE as given in (1.22) and ( 1

Examples

Algorithm 3.1 deduced from Corollary 3.2 and Theorem 3.1 is used to compute the overapproximation Π defined in (3.24) (subset of Π * defined in (3.14)) of the reachable set R( P 0 , t) of the system S( P 0 , t) (described in Section 3.4), t ≥ 0. Several examples are treated. With these examples, we provide some performance evaluations of our approach.

Examples from COMPleib

To evaluate the performance of our approach, we compute an overapproximation of the reachable set for several real-life systems from the COMPleib library [Leibfritz, 2006]. For each system, a stabilizing controller is generated for the generalized plant (by using the h2syn function of Matlab), then the system is reduced using a balanced

Paul Rousse

Set-based Cosimulation

Examples

truncation method to a given state space size (by using the balred function of Matlab). The set of initial states is chosen such that the quadratic term belongs to the set of stable solutions to the associated Continuous Algebraic Riccati Equation. The simulation are ran for an input u(t) = [ 1 ... 1 ] ⊤ exp(-t) for t ∈ [0, 2]. Each ODE is numerically integrated using the ode113 solver in Matlab. Finally, we run the simulation with one time-varying paraboloid and then multiple time-varying paraboloids. CPU time performances for a computer with an Intel i5 2.5GHz are presented in Table 3.1 and Table 3.2.

In Figure 3.13, we show several runs for the examples. Each paraboloid is overapproximated with a box, we show the intersection of these intervals.

Performance is mainly dependent on the number of paraboloids that we consider, and our ability to efficiently solve the DRE. Table 3.1: Computation times (in seconds) of the overapproximation for different systems sizes, using a unique time-varying paraboloid. (When the size of the original system is smaller than the required reduced system size, then the model reduction is not applicable -n.a.-.)

System verification

We study the stable IQC system S( P 0 , t), defined in (3.8), at a given time t in [0, 1], for a parabolic set of initial states P 0 = Parab(E 0 , f 0 , g 0 ), with E 0 = a+b a a a+b , f 0 = [ 0 0 ] , g 0 = -0.015, a = 10 -2 and b = 10 -6 , and for the following parameters

A = -I, B = I, C = 0, M = I 0 0 0 1 0 0 0 -2I
where I = [ 1 0 0 1 ], and with a zero input signal u. The reachable set R( P 0 , t) of S( P 0 , t), defined in (3.9), is computed using (3.13) and Theorem 3. 4) n.a.

Table 3.2: Computation times (in seconds) and number of paraboloids (in parenthesis) of the overapproximation for different systems sizes. (When the size of the original system is smaller than the required reduced system size, then the model reduction is not applicable -n.a.-.) R( P 0 , t) set at time t = 0.794 and its projection R( P 0 , t)| x over the LTI state space (i.e. projection over (x 1 , x 2 ) states). In Figure 3.14b, the constraints boundaries ∂P(t) (for P ∈ Π * , Π * defined in Section 3.6) are touching the reachable set R( P 0 , t). The non-convexity of R( P 0 , t) arises from the non-positive solutions to the DRE (1.21). Figure 3.14c represents the projection of the reachable tube t → R( P 0 , t) projected over the LTI dimension (x 1 , x 2 ).

Delayed system

IQC relationships can be derived for delayed systems (see [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF]). The delay operator is a linear time-invariant system that has a state of infinite dimension. The states dynamic can be described by a wave partial differential equation. Projections of the state over a base of Legendre polynomials (of maximal degree r ∈ N) have a linear dynamic that only depends on smaller degree projection. Moreover, the error between the true state and the projections satisfies energetic constraints (that is derived from Jensen inequality). By increasing the maximal degree r of the considered polynomials, for similar input, the set of reached output is strictly reduced. For each degree r, the IQC falls into the context of this work since M r w ≺ 0. In the sequel, the reachable set of this overapproximating model is computed and plots of the reachable outputs are given. Consider the following delayed system where s ∈ C is the Laplace variable, τ > 0 and a delay h > 0 and the input signal is

Set-based Cosimulation

Paul Rousse defined as follows

u(t) = 0 if t < 0 or t > 5 0.3(1 -t 5 )t + 0.1sin(2πt) otherwise.
Let S r ( P 0 , T ) for g ∈ N be the relationship of maximal polynomial degree r, and let R r (Z 0 , t) be the associated reachable set at time t > 0 for a given set of initial states. The reachable set R r for a set of initial states P 0 is computed using (3.13) and Theorem 3.1 for different orders r of the hierarchy. Figure 3.15 is a plot of the projections of the reachable set over the output map (where r = 1, 2, 3) together with the trajectory of the delayed system. We have the following relationship: R 3 ( P 0 , t) ⊂ R 2 ( P 0 , t) ⊂ R 1 ( P 0 , t).

Discussion

IQC theory was introduced in [Rantzer and [START_REF] Rantzer | [END_REF]] as a method to assert robust stability and performance for uncertain LTI systems. It merges the results from [Yakubovich, 1967] and [START_REF] Zames | On the Input-Output Stability of Time-Varying Nonlinear Feedback Systems-Part II: Conditions Involving Circles in the Frequency Plane and Sector Nonlinearities[END_REF], Zames, 1966a]. The uncertain system is represented as an interconnection between an LTI system with an unknown block that is described by a set of to input-output energy relationships. The Kalman-Yakubowich-Popov lemma gives equivalence between the stability of the interconnection and a Linear Matrix Inequality. These IQCs are represented in a frequency domain. However, some of them can be equivalently expressed over a finite-horizon in the temporal domain. Such IQCs

Paul Rousse

Set-based Cosimulation 3.9. Discussion are called hard IQCs (or complete IQCs) and have been studied in [Megretski, 2010]. These IQCs allows to derive bounds over the state at any given time and are therefore of high interest in the verification literature. More recently, [START_REF] Scherer | Stability analysis by dynamic dissipation inequalities: On merging frequency-domain techniques with time-domain conditions[END_REF], showed that any IQCs can always be expressed as a finite-horizon IQC plus a state-dependent quantity. Such an inequality allows to derive bounds over the reachable tube of the uncertain system.

When computing bounds over the reachable set of an IQC system, the IQC is integrated into a storage function. The IQC is scaled with multiplier and then integrated into the storage function. The class of multipliers used to handle IQC constraint is closely related to the method used to solve the DRE. [Jönsson, 2002, Savkin and Petersen, 1995, Savkin and Petersen, 1996a[START_REF] Savkin And Petersen ; Savkin | Robust state estimation for uncertain systems with averaged integral quadratic constraints[END_REF] use constant multipliers, i.e. the IQC constraint is integrated using the S-procedure (see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], Section 2.6.3). In the study of LTI IQC systems, multipliers can be chosen as positive operators over the signal space (i.e. positive LTI systems, see [START_REF] Fetzer | Invariance with dynamic multipliers[END_REF], Scherer and Veenman, 2018, Yin et al., 2020]). Recursive methods based over integration of the DRE allow to study LTV systems and thus time-varying multipliers.

The existence of an overapproximation over a given time-horizon is dependent on the existence of a solution to a DRE or a DLMI.

The existence of a solution to the DRE has been extensively studied in several works (see a review of the Riccati equation in [Kučera, 1973] and in [START_REF] Bittanti | The Riccati Equation[END_REF]). The DRE equation is known to diverge in finite-time depending on its initial condition and on the parameters. [START_REF] Savkin | Recursive State Estimation for Uncertain Systems with an Integral Quadratic Constraint[END_REF] provides the existence of a solution to the DRE over any time-horizon for a subclass of IQC. [Jönsson, 2002] gives equivalent conditions between the existence of a solution to the DRE and a full-rank condition over the Hamiltonian. [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF] shows (in Theorem 1) that the DRE is equivalent to a DRI (Differential Riccati Inequality). Contrary to the DRE, the DRI can be expressed as a DLMI by using the Schur complement, thus the existence of a solution to the DRE is then expressed as a feasibility problem of an infinite-dimensional linear problem. The problem can then be solved over a finite basis of functions. An SDP solver can then be used to find the optimal weights and solution of the DRE.

In our work, we provide (for the LTI case) a theorem stating that overapproximations of the reachable always exist over any time-horizon. Such a result is obtained by using time-varying multipliers for the IQC constraint.

Set-based Cosimulation

Paul Rousse

Chapter 3. Integral Quadratic Constraints 88

Conclusion

This chapter applied the framework described in Chaper 1 to overapproximate the reachable sets of a linear time-varying systems with an unknown input bounded by integral quadratic constraints, modeling e.g. delay, rate limiter, or energy bounds. We define a family of paraboloidal overapproximations. These paraboloids are supported by the reachable tube on touching trajectories. Parameters of each paraboloid are expressed as a solution to an initial value problem. Our integration scheme is not guaranteed and the paraboloids we compute are subject to the error of the differential equation numerical integration. The next chapter, Chapter 4, investigates the use of interval arithmetic and validated integration scheme to derive a guaranteed overapproximation to the reachable set of a system subject to an IQC disturbance.

Future work

We showed that the reachable set can be described as an intersection of uncountably many paraboloids. In our implementation, a subset of these time-varying paraboloids is computed to overapproximate the reachable set. Then, we compute a minimal volume paraboloid that contains the intersection of all the paraboloids. The computation time of our method is directly dependent on the number of time-varying paraboloids. Finding only one time-varying paraboloid which minimizes its end volume would avoid integrating multiple time-varying paraboloids. Solutions exist for this optimization problem.

The differential Riccati equation can be weakly solved using a basis of polynomial solutions (as in [START_REF] Seiler | Finite horizon robustness analysis of LTV systems using integral quadratic constraints[END_REF]). Then Sum-Of-Square relaxation provides a suboptimal overapproximating paraboloid. Previous works implementing this approach use conservative overapproximations that do not fully incorporate the state constraint. In future works, we will develop such an approach with the results presented in this chapter. A locally optimal solution of the optimization problem can be derived using the Pontryagin's Maximum Principle. We provided such a solution for in the previous chapter, Chapter 2. This could be adapted as well for the IQC case.

In our implementation, the time-varying multipliers and initial multiplier are chosen such that some touching trajectory validate the constraint in the future. Other criteria could be derived such as studying the average behaviors of the trajectories. Since the computational complexity is linear in the number of time-varying paraboloids that need to be simulated, an efficient choice of the multipliers can lead to algorithms that demand less computational resources.

Paul Rousse

Set-based Cosimulation (3.10) can be derived solving the following optimal control problem (for t > 0) max

w∈L 2 ([0,t];R m ) t 0 x(τ ) u(τ ) w(τ ) M x(τ ) u(τ ) w(τ ) dτ -x q,t s.t. ẋ = Ax + Bw + Cu x(t) = x t
for given (x t , x q,t ) ∈ Z + . This is a special instance of the LQR problem (see e.g. [START_REF] Savkin And Petersen ; Savkin | Robust state estimation for uncertain systems with averaged integral quadratic constraints[END_REF]). For x ∈ L 2 ([0, T ]; R n ) a touching trajectory, let

n = Ex -f. (3.26)
Using (3.10), n satisfies the following differential equation

ẋ ṅ = L x n + N u where L = A -BM -1 w M ⊤ xw -BM -1 w B ⊤ -(M x -M xw M -1 w M ⊤ xw ) -A ⊤ + M xw M -1 w B ⊤ and N = C -BM -1 w M ⊤ uw -(M xu -M xw M -1 w M ⊤ uw )
.

The value function evaluated along the touching trajectory x is then obtained by introducing the parameter

r = g -f ⊤ x (3.27) Set-based Cosimulation Paul Rousse which satisfies ṙ = u H R x n u with H = M uw M -1 w M xw -M ⊤ xu -(C -M uw M -1 w B ⊤ ) and R = M u -M uw M -1 w M ⊤ uw .
Using (3.26) and (3.27), the value function satisfies:

p(t, z(t)) = x(t) ⊤ n(t) + r(t) + x q (t).
Let the time-varying paraboloid P = T ( P 0 ) such that P diverges in finite-time, i.e. T P ( P 0 ) < ∞. Since all the touching trajectories are continuous in time, each touching trajectory is defined over [0, T P ( P 0 )]. Their corresponding value function h evaluated along the touching trajectory is as well continuous over [0, T P ( P 0 )]. Therefore, one can extend the definition of P until T P ( P 0 ) using the continuity of the value function

P(T ) = {z ∈ R n+1 | lim t→T t<T p(t, z) ≤ 0}.
(3.28)

where T = T P ( P 0 ). We state this result in the following property Proposition 3.10. Continuous extension For any P = T ( P 0 ), if the quadratic coefficient of the time-varying paraboloid set P diverges in finite-time, then the extension to the right of P exists and is defined by (3.28).

Part II Interval Analysis Methods

91

In Part I, we applied the levelset method to overapproximate the reachable tube of a linear time-varying system subject to bounded disturbances. In this part, we propose another approach based on the framework of Interval Arithmetic and validated integration to overapproximate the reachable tube of a nonlinear system subject to a disturbance bounded by an Integral Quadratic Constraint (IQC, see Chapter 3).

The Interval Arithmetic has been introduced in the '50s and '60s as a numerical method to evaluate mathematical expressions while embedding rounding errors. These rounding errors were represented as intervals of values and propagated by structural induction within the mathematical expression. Thus, the exact evaluation of this expression is guaranteed to belongs to this interval despite the use of approximated numerical real arithmetic (e.g. the floating-point arithmetic). The initial Interval Arithmetic was designed to evaluate expressions with additions, subtractions, and multiplications, and was later on extended to divisions, Taylor series and integrals. More complicated problems have been studied, and in particular, root equations and fixed-point equations can as well be addressed within this framework. For such equations, the exact solution is often not computer-representable (it can be an irrational number) and Interval Arithmetic provides a practical tool to find bounds for which the exact solution is guaranteed to lie in.

In dynamical systems analysis, continuous trajectories are solutions to an Initial Value Problem (IVP)

ẋ(t) = f (t, x) x(0) = x 0
where x ∈ L 2,loc (I; R n ), with I is the domain definition of the solution x. Such an IVP can be highly sensitive on the numerical error and providing bounds containing the exact solutions is often necessary. This IVP can, in fact, be conveniently and equivalently expressed by the following fixed-point equation:

x(t) = x 0 + t 0 f (τ, x(τ ))dτ.
In this form, it is possible to use the Interval Arithmetic framework to compute a guaranteed solution to the IVP. [START_REF] Moore | Introduction to Interval Analysis[END_REF] described a method where the trajectory x is represented as a union of interval in the time and state space. Since then, more problems in the field of analysis of dynamical systems have been addressed such as Differential Algebraic Equations (DAE), guaranteed estimations (where the measurement noise is chosen in an interval of value), and, in particular, for system verifications. The interval evaluation allows to model unknown inputs such as unknown initial value (when x 0 belongs to a set of initial states) and input disturbances (such as a disturbance w(•) where w(t) ∈ [-1, 1] m at every t ∈ I).

In such a case, the same validated numerical integration method can be applied to compute an overapproximation of the reachable tube of a dynamical system.

In this part, we overapproximate the reachable set of a dynamical system subject to a disturbance bounded by an integral constraint. Contrary to Chapter 3, the system of interest has a non-linear dynamic, and the set of disturbances is described by a non-linear integral constraint.

Non-linear systems can already be studied with IQC models. To do so, the image of the non-linear block is modeled as an unknown but bounded disturbance. Then, the methods presented in Part I can be used for the class of systems studied in this chapter. However, there is a practical difficulty. Getting this IQC model cannot be done automatically out of the dynamical function. Most of the time, the IQC model is obtained manually by identifying the non-linearities in the model and by overapproximating them with known IQCs.

To automatize the analysis of such systems, a possible approach is to study the syntax of the dynamical system models. Usually, the non-linear system is given as a formula involving only a few elementary operations. When these operations can be overapproximated, the resulting expression can be as well overapproximated. The interval arithmetic and its associated validated numerical integration framework use such an approach to overapproximate the reachable tube of a non-linear system. The dynamical function (whether it is linear or not) can be automatically overapproximated by syntactic decomposition over its expression.

In the classical validated numerical integration framework introduced in [START_REF] Moore | Introduction to Interval Analysis[END_REF], models of interest are usually dynamical systems with an unknown disturbance bounded by an ∞-norm constraint. These frameworks do not take into account disturbances defined by an integral constraint (as presented in Chapter 3). However such models are interesting as they can model complex systems (such as systems with internal delays).

In Chapter 4, we present the classical interval arithmetic framework and the validated numerical integration framework, then, we extend the last to overapproximate the reachable set of non-linear system subject to a disturbance described by a nonlinear integral constraint.

Paul Rousse

Set-based Cosimulation Chapter 4

Validated Integration of Nonlinear Systems subject to Integral Constraint In this chapter, we present a method to compute the reachable set of a nonlinear dynamical system subject to an unknown disturbance described by an integral constraint between the disturbance and the state trajectory. The Interval Arithmetic and validated numerical integration frameworks are used. With additional assumptions about the dynamic of the disturbance, the integral constraint gives bounds over the set of disturbances. A contractor over the set of reachable states is defined out of these bounds. This contractor is then used in a fixed point algorithm with a propagation step (as described in [Alexandre dit Sandretto and Chapoutot, 2016]). Our algorithm is implemented using the DynIbex library [START_REF] Dit | [END_REF] and applied to overapproximate the reachable tube of a dynamical system with an internal delay.

In Section 4.1, we describe the basics of interval analysis. In Section 4.2, we give a short introduction of validated numerical integration. In Section 4.3, we present the nonlinear system with disturbances subject to an integral constraint. In Section 4.4, we extend the validated numerical integration presented in Section 4.2. In Section 4.5, we present an example and compare it to the results obtained in Chapter 3. In Section 5.10 we conclude this chapter.

Interval arithmetic

A simple and common way to represent and manipulate sets of values is interval arithmetic (see [START_REF] Moore | Introduction to Interval Analysis[END_REF]). An interval [x i ] = [x i , x i ] defines the set of reals x i such that x i ≤ x i ≤ x i . IR denotes the set of all intervals over reals. The size (or width) of [x i ] is denoted by w([

x i ]) = x i -x i .
Interval arithmetic extends to IR elementary functions over R. For instance, the interval sum, i.e., [

x 1 ]+[x 2 ] = [x 1 +x 2 , x 1 +x 2 ]
, encloses the image of the sum function over its arguments.

An interval vector or a box [x] ∈ IR n , is a Cartesian product of n intervals. The enclosing property defines what is called an interval extension or an inclusion function.

Definition 4.1. Inclusion function

Consider a function f : R n → R m , then [f ] : IR n → IR m is said to be an inclusion function of f to intervals if

∀[x] ∈ IR n , [f ]([x]) ⊇ {f (x), x ∈ [x]} .
It is possible to define inclusion functions for all elementary functions such as ×, ÷, sin, cos, exp, etc. The natural inclusion function is the simplest to obtain: all occurrences of the real variables are replaced by their interval counterpart and all arithmetic operations are evaluated using interval arithmetic. More sophisticated inclusion functions such as the centered form, or the Taylor inclusion function may also be used (see [START_REF] Jaulin | Applied Interval Analysis[END_REF] for more details).

Paul Rousse

Set-based Cosimulation In the first example of division, the result is the interval containing all the real numbers because the denominator contains 0.

As an example of inclusion function, we consider a function p defined by p(x, y) = xy + x .

The associated natural inclusion function is

[p]([x], [y]) = [x][y] + [x],
in which variables, constants and arithmetic operations have been replaced by its interval counterpart. And so

p([0, 1], [0, 1]) = [0, 2] ⊆ {p(x, y) | x, y ∈ [0, 1]} = [0, 2].

System of equations

The Interval Arithmetic framework can be used to find bounds containing the solutions of a set of equations. In Section 4.1.1, we present a method to solve a fixed-point equation. In Section 4.1.1, we present a method to solve a system of equations that is not given in the form of a fixed point equation.

Fixed-point equation

To solve the fixed-point equation

x = f (x) (4.1)
where f : R n → R n , it is possible to compute the fixed-point iterations

x k+1 = f (x k ) Set-based Cosimulation Paul Rousse for x 0 ∈ R n . If the self-map f is contractive, i.e. if f (x) -f (y) ≤ k x -y
where 0 ≤ k < 1, then

• the equation ( 4.1) has a unique fixed point x * , and

• the sequence {x k } k∈N converges to this fixed point x * .
This result is known as the Banach Theorem (Section 1.6 in [Zeidler, 1995a]). It gives a practical algorithm to compute an approximate solution to (4.1). This iteration method has been adapted to the framework of interval arithmetic (e.g. in Chapter 6 of [START_REF] Moore | Introduction to Interval Analysis[END_REF]). Let [f ] be the interval evaluation of f . If [f ] is contractive (with respect to the norm dist1 ), i.e. for every

[x] ∈ IR n , it holds dist ([f ]([x]), [f ]([y])) ⊆ k dist ([x], [y]), where 0 ≤ k < 1, then the sequence of [x k+1 ] defined by [x k+1 ] = [f ]([x k ]) [x 0 ] ∈ IR n (4.2)
converges to the singleton {x * }, i.e.

[x k ] → {x * }. (4.
3)

The convergence property (4.3) is not satisfying since the sequence of [x k ] are not guaranteed to contain the fixed point solution x * . Such a result can be obtained by observing that

x * ∈ [x] ⇒ x * ∈ [f ]([x])
for every [x] ∈ IR n . Then, the interval version of the iteration method can be applied for the case where

• [f ] is contractive (with respect to the set inclusion), i.e. [f ]([x]) ⊆ [x]
, and

• [x 0 ] contains the solution x * to the fixed-point equation (4.1).

in such a case, the sequence of [x k ] defined by

[x k+1 ] = [f ]([x k ]) (4.4) satisfies x * ∈ [x k+1 ] ⊆ [x k ] ⊆ [x k-1 ] ⊆ • • • ⊆ [x 0 ].

Interval arithmetic

Thus, each iterate [x k ] is a sound approximation of x * . When the contractive property is not available, the operator [f ] can be enforced to be contractive by using the operator [f ] defined by

[f ]([x]) = [f ]([x]) ∩ [x].
This iterative approach to solve the fixed point equation ( 4.1) has been improved by using other existing fixed-point algorithm. Chapter 8 of [START_REF] Moore | Introduction to Interval Analysis[END_REF] The iteration sequence converges to [x ∞ ] = 0.666666 ± 4e -6. The error is due to numerical imprecision introduced by the floating-point arithmetic.

Contractors When the equation to solve can be represented as a fixed point equation, the previous section proposes iterative methods that can refine an a priori bound containing a solution of this equation. When the equation is not a fixed point equation, it is possible to use similar concepts: the contractor. A contractor is an operator that associates to every given set, a subset that contains all the points where the constraint is verified (see [START_REF] Chabert | Contractor programming[END_REF]).

Definition 4.2. Contractor

For a constraint f that maps R n to a truth value, a contractor Ctc of f associates to a subset of R n to a subset of R n . For any [b], [b ′ ] ∈ IR n , Ctc must verify the following properties:

• the contraction: Ctc ([b]) ⊆ [b],
• the conservativeness:

∀x ∈ [b] \ Ctc ([b]) , f (x) is not satisfied, • the monotonicity: [b ′ ] ⊆ [b] ⇒ Ctc ([b ′ ]) ⊆ Ctc ([b])
Set-based Cosimulation Paul Rousse As for the fixed point equation and thanks to the contraction and monotonicity property, if a solution x * belongs to the truth value of f , then every sequence of [x k ] defined by

[x k+1 ] = Ctc ([x k+1 ]) , such that x * ∈ [x 0 ], satisfies x * ∈ [x k+1 ] ⊆ [x k ] ⊆ [x k-1 ] ⊆ • • • ⊆ [x 0 ].
The precision and the speed of convergence of the sequence of [x k ] depends on the contraction of Ctc.

Example 4.3.

To find the set of x such that f (x) ≥ 0 where

f (x) = 2 -x 4 (4.6)
it is possible to define a contractor that exploit the concavity of f (see Figure 4.2).

The resulting contractor produces a decreasing (in the set inclusion sens) sequence of iterates 1.189207, 1.189207] (see Figure 4.3).

[x k ] for an initial [x 0 ] = [-2, 2]. The sequence converges to [x ∞ ] = [-

Paul Rousse

Set-based Cosimulation where

Φ(y)(t) = y 0 + t 0 F (s, y(s))ds with F : R + × R n → R n , y 0 ∈ R n
, and with the unknown y : R + → R n . The function y : R + → R n can be represented with a union of intervals in the time and state space, i.e.

{(t, y(t)) | t > 0} ∈ k [ y k ]
where [ y k ] ∈ IR n+1 (see Figure 4.4). With such a representation, one can solve the fixed point equation (4.8) using the fixed point iteration proposed in Section 4.1.1. This section details this approach for a more general IVP than (4.7).

Mathematically, differential equations have no explicit solutions, except for a few particular cases. Nevertheless, the solution can be numerically approximated with the help of integration schemes such as Taylor series [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] or Runge-Kutta methods [Sandretto andChapoutot, 2016, Dit Sandretto and[START_REF] Sandretto | [END_REF].

In the following, we consider a generic parametric differential equation as an interval initial value problem (IIVP) defined by

     ẏ = F (t, y, x, p, u) 0 = G(t, y, x, p, u) y(0) ∈ Y 0 , x(0) ∈ X 0 , p ∈ P, u ∈ U, t ∈ [0, t end ] , (4.9)
with

F : R × R n × R m × R r × R s → R n and G : R × R n × R m × R r × R s → R m .
The vector variable y of dimension n is the differential variable while the variable x is an algebraic variable of dimension m with an initial condition y(0) ∈ Y 0 ⊆ R n and x(0) ∈ X 0 ⊆ R m . In other words, differential-algebraic equations (DAE) are considered, and in the case of m = 0, this differential equation simplifies to an ordinary differential equation (ODE). Note that usually, the initial values of algebraic variable

Paul Rousse

Set-based Cosimulation

Validated numerical integration methods

x are computed by numerical algorithms used to solve DAE but we consider it fixed here for simplicity. Variable p ∈ P ⊆ R r stands for parameters of dimension r and variable u ∈ U ⊆ R s stands for a control vector of dimension s. We assume standard hypotheses on F and G to guarantee the existence and uniqueness of the solution to such a problem.

A validated simulation of a differential equation consists of a discretization of time, such that t 0 • • • t end , and a computation of enclosures of the set of states of the system y 0 , . . . , y end , by the help of a validated integration scheme. In details, a validated integration scheme is made of • an integration method Φ(F, G, y j , t j , h), starting from an initial value y j at time t j and a finite time horizon h (the step-size), producing an approximation y j+1 at time t j+1 = t j + h, of the exact solution y(t j+1 ; y j ), i.e., y(t j+1 ; y j ) ≈ Φ(F, G, y j , t j , h);

• a truncation error function lte Φ (F, G, y j , t j , h), such that

y(t j+1 ; y j ) = Φ(F, G, y j , t j , h) + lte Φ (F, G, y j , t j , h).
Basically, a validated numerical integration method is based on a numerical integration scheme such as Taylor series [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] or Runge-Kutta methods [Sandretto andChapoutot, 2016, Dit Sandretto and[START_REF] Sandretto | [END_REF] which is extended with interval analysis tools to bound the local truncation error, i.e., the distance between the exact and the numerical solutions. Such methods work in two stages at each integration step, starting from an enclosure [y j ] ∋ y(t j ; y 0 ) at time t j of the exact solution, we proceed by: i. a computation of an a priori enclosure [ y j+1 ] of the solution y(t; y 0 ) for all t in the time interval [t j , t j+1 ]. This stage allows one to prove the existence and the uniqueness of the solution.

ii. a computation of a tightening of state variable [y j+1 ] ∋ y(t j+1 ; y 0 ) at time t j+1 using [ y j+1 ] to bound the local truncation error term lte Φ (F, G, y j , t j , h).

A validated simulation starts with the interval enclosures [y(0)], [x(0)], [p] and [u] of respectively, Y 0 , X 0 , P, and U . It produces two lists of boxes:

• the list of discretization time steps: {t 0 , . . . , t end };

• the list of state enclosures at the discretization time steps: {[y 0 ], . . . , [y end ]};

• the list of a priori enclosures: {[ y 0 ], . . . , [ y end ]}.

Set-based Cosimulation

Paul Rousse

System with integral constraint over the state

We presented the classical framework of interval arithmetic and the framework of validated numerical integration. Such a framework handle systems that are subject to disturbances bounded by the ∞-norm (i.e. nonlinear systems subject to a disturbance w such that w(t) ≤ 1 for example). However, the case where the disturbance is subject to an integral constraint has never been addressed until now.

The next sections propose a method to overapproximate the set of reachable states of a nonlinear system subject to a disturbance bounded by an integral constraint. We assume the disturbance and its time derivative to be bounded in ∞-norm (i.e. to be bounded at any time). These bounds are used to get a first overapproximation of the reachable tube. This coarse overapproximation might contain a set of trajectories not satisfying the integral constraint. We use a contractor operator (as introduced in Section 4.1.1) in order to reduce the initial prior overapproximation of the disturbance set. This new disturbance is then reused to get a new (smaller) overapproximation of the reachable tube. These two operations are iterated until a fixed point is reached.

Set-based Cosimulation

Paul Rousse with parameters (4.11). The system has a unique trajectory (in red) that is overapproximated by the union of intervals [ y k ] (in grey).

Consider the system:

ẋ = f (t, x, w) x(0) ∈ x 0 (4.13)
where w is an unknown disturbance in L 2,loc (R + ; R m ) that satisfies the integral constraint, for any τ ≥ 0:

τ 0 w(s) 2 ds ≤ τ 0 g(s, x(s))ds (4.14)
where g : R + × R n is a given function.

Remark 4.1. IQC bounds

The integral constraint does not give any bound on the disturbance as it can be easily understood from the unit energy disturbed system ) is satisfied for every ǫ > 0, however, no bounds can be determined for w since w (0) → ∞ when ǫ → 0. Please note that the system defined in (4.15) has a bounded reachable set even if the disturbance cannot be bounded at any given time (see [Boyd et al., 1994, Chapter 8.1.2]).

         ẋ = -x + w x(0) = 0 1 ≥ 1 0 w 2 (τ )dτ . ( 4

Validated numerical integration for dynamical systems subject to integral constraints

This section presents the main contribution of this chapter. For a system described by ( 4.13) and subjects to the integral constraint defined by (4.14), we compute an overapproximation of its reachable tube over the time domain [0, T ], where the time horizon T > 0 is given. A first overapproximation of the reachable tube is computed using pessimistic bounds over the disturbances. The integral constraint in (4.14) is used to derive a contractor. This contractor and a propagation step are applied in a fixed point algorithm until a contraction factor is reached. We run the algorithm over a simple example.

Extended system

We extend the state of the system with the integral value corresponding to the integral constraint in (4.14):

ż(t) = g(t, x(t)) -w(t) 2 z(0) = 0 (4.16)
Then, (4.14) can be equivalently expressed for z:

∀t ∈ R + , z(t) ≥ 0. (4.17)
As mentioned in Remark 4.1, no L ∞ bounds can be derived for L 2 bounded signals. To study such systems, we make further assumptions about the disturbance: Assumption 4.1. Continuous disturbance signal w is continuous, differentiable, and of continuous derivative over R + .

Set-based Cosimulation Paul Rousse

This assumption seems reasonable in the case of real systems modeling since disturbances modeled by integral constraints correspond to physical quantities. Since the continuity of a function over a closed interval implies its boundedness, Assumption 4.1 implies that the signal w is bounded and of bounded variation over [0, T ]. Therefore, there exists [w] ∈ IR m and [w ′ ] ∈ IR m such that for all t ∈ [0, T ]:

w(t) ∈ [w] ẇ(t) ∈ [w ′ ] (4.18)
Using Assumption 4.1 and (4.16), the following system will be studied:

S :                          ẋ(t) = f (t, x(t), w(t)) ż(t) = g(t, x(t), w(t)) -w(t) 2 ẇ(t) ∈ [w ′ ] x(0) ∈ [x 0 ] z(0) = 0 0 ≤ z(t) w(t) ∈ [w] (4.19)
where [x 0 ] ∈ IR n is the set of initial states. We use the following notation (x, z, w)

∈ S iff (x, w) ∈ L 2,loc ([0, T ]; R n ) × L 2,loc ([0, T ]; R) × L 2,loc ([0, T ]; R m ) is a trajectory of S.
(4.18) gives prior bounds on the disturbance w. They can be used to propagate the trajectories using standard validated numerical integration frameworks. Thanks to this, we get a first a priori overapproximation of the reachable set. In the next section, we use this first overapproximation and a contractor (defined out of the integral inequality) in a fixed point algorithm in order to get a tighter overapproximation of the reachable set.

Bounds over the disturbance

In this section, (4.18) and the integral constraint in (4.17) are used to derive bounds over the disturbance w. These bounds are then used to define a contractor over the a priori enclosure of the trajectories.

We present a preliminary result before deriving bounds over the disturbance w:

Proposition 4.1. Overapproximate intersection For [v] ∈ IR p , p ∈ N and r > 0, if v ≤ r then [v] ⊂ [-r, r] p .
Proof. In an Euclidean space, the norm 1 and norm 2 satisfy

v 2 1 + • • • + v 2 p ≤ |v 1 | + • • • + |v p | for any (v 1 , . . . , v p ) ∈ R p . ♦ Paul Rousse Set-based Cosimulation
h 0 w(τ ) 2 dτ ≤ r ⇒ ∀τ ∈ [0, h], w(τ ) ∈ [W r ],
where

[W r ] = [-k, k] n with k = r h + h 2 w ′ (where w ′ is the maximum Euclidean norm over the elements of [w ′ ]).
Proof. By applying the Cauchy-Schwartz inequality between the signal w and t → 1 for the inner product of square-integrable function, we have:

h 0 w(τ )dτ 2 ≤ h h 0 w(τ ) 2 dτ ≤ hr. By (4.18), w(τ ) = w 0 + τ 0 w 1 (κ)dκ with w 0 ∈ [w] and w 1 (•) ∈ [w ′ ].
Using the reverse triangular inequality, we have:

h 0 w 0 dτ ≤ √ rh + h 0 τ 0 w 1 (κ)dκ .
Then, we get:

hw 0 ≤ √ hr + h 2 2 w ′ . (4.20)
This relationship is derived over [0, h] but is also valid for any time interval [t, t + h] of width h, t > 0. Therefore, by using Proposition 4.1 and (4.20), we have:

∀τ ∈ [0, h], w(τ ) ∈ [W r ]. ♦
We then use Proposition 4.2 to derive bounds in the specific case of (4.16). Consider a system trajectory (x, z, w) ∈ S, such that at a given t ∈ [0, T ] and h > 0 s.t. t + h ∈ [0, T ], and for all τ ∈ [t, t + h]:

(x(t), z(t), w(t)) ∈ [y t ] (x(τ ), z(τ ), w(τ )) ∈ [ y t ]
where

[y t ] = [x t ] × [z t ] × [w t ] [ y t ] = [ x t ] × [ z t ] × [ w t ] . (4.21)
The trajectories belong to [y t ] at t and are in [ y t ] between [t, t + h]. At t + h, for a given t ≥ 0 and a given h ≥ 0, (4.19) implies that z satisfies:

z(t + h) = z(t) + t+h t g(t, x(t))dτ - t+h t w(τ ) 2 dτ.

Set-based Cosimulation Paul Rousse

By applying (4.17) at t+h implies that z(t+h) ≥ 0, we have the following relationship:

t+h t w(τ ) 2 dτ ≤ z(t) + t+h t g(τ, x(τ ))dτ. (4.22) Let a function q(z, x) = z + t+h t g(τ, x(τ ))dτ. (4.23)
By using an interval evaluation [q] of q, the upperbound of q(z, x) can be evaluated

for z ∈ [z t ] and x ∈ [ x t ]. We denote by [q]([z t ], [ x t ]) this upperbound. For any w ∈ L 2,loc ([t, t + h], [ w t ]), 4.22 implies: t+h t w(τ ) 2 dτ ≤ [q]([z t ], [ x t ]).
Then, Proposition 4.2 can be used to derive bounds over the disturbance w:

Proposition 4.3. Disturbance overapproximation For a w ∈ L 2,loc ([t, t + h]; R m ) defined over an interval of length h > 0, t > 0. If w satisfies (4.18) (with given bounds [w], [w ′ ] ∈ IR m ), then for any τ ∈ [t, t + h]: w(τ ) ∈ [W q ], (4.24)
where

[W q ]([ x t ], [z t ]) = [-r, r] m with r = [q]([z t ],[ xt])
h and q defined in (4.23).

Proof. This is a direct application of Proposition 4.

♦

We then define the operator over [y t ] and [

y t ] C([y t ], [ y t ]) = ([y t ] ∩ [Y g ]([ x t ], [z t ]), [ y t ] ∩ [Y g ]([ x t ], [z t ])) (4.25)
where [y t ] and [ y t ] are defined in (4.21),

[Y g ] = [-∞, ∞] n × [0, ∞] × [W q ],
with [W q ] defined in Proposition 4.3.

Proposition 1. C defined in (4.25) is a contractor. Proof. By Proposition 4.3, we have, for τ ∈ [t, t + h], w(τ ) ∈ [W q ],

Paul Rousse

Set-based Cosimulation 4.4. Dynamical systems subject to integral constraints i.e., all the disturbance signals of S belong to [W q ], so the contractor is conservative. Since the contractor is defined as an intersection with [y t ] and [ y t ] respectively, we have

([y t ], [ y t ]) ⊆ C([y t ], [ y t ]), C is contractive. For any ([y ′ t ], [ y ′ t ]) such that [y ′ t ] ⊆ [y t ] and [ y ′ t ] ⊆ [ y t ], C([y ′ t ], [ y ′ t ]) ⊆ C([y t ], [ y t ]),
i.e. C is monotone. ♦

Integral constraint propagation

The contractor defined by (4.25) is used in a fixed point algorithm as in [Alexandre dit Sandretto and Chapoutot, 2016]. A priori enclosure of the trajectory is computed using bounds (4.18) over w. The integration algorithm gives

• the discretization time steps: {t 0 , . . . , t end };

• the state enclosure at the discretization time steps:

Y 0 = {[y 0 0 ], . . . , [y 0 end ]};
• the a priori enclosures:

Y 0 = {[ y 0 0 ], . . . , [ y 0 end ]}.
We then apply the contractor over each couple of discretized time-step boxes [y 0 i ] ∈ Y 0 and their associated a priori enclosures [ y 0 i ] ∈ Y 0 . These two steps are repeated in a fixed point algorithm until the contraction factor is lower than a given value. In this approach, time steps are computed at the first iteration of the algorithm and are not updated.

Example 4.5.

We study the following linear time-invariant system disturbed by an unknown signal w constrained by a 2-norm inequality:

         ẋ(t) = -x(t) + w(t) t 0 w(τ ) 2 dτ ≤ t 0 0.01x(τ ) 2 dτ x(0) ∈ [-1, 1] (4.26) with [w] = [-1, 1] and [w ′ ] = [-1, 1] in (4.18) for t ∈ [0, 2.5]. Figure 4
.7 shows the reachable set of this dynamical system computed with the method described in this section.

Set-based Cosimulation

Example

In this part, we present an application of the method described in Section 4.4 to a system with a time delay.

For u, v ∈ L 2,loc (R + ; R), the delay operator D h over an input signal u is defined by the following relationship:

v = D h (u) ⇔ v(t) = u(t -h) for all t ≥ h v(t) = 0 otherwise. (4.27)
Validated numerical integration of differential equations with delays is challenging. Since they act as a memory of the past input signal over an interval of width h, the state of the delay belongs to L 2,loc ([0, h], R). The dimension of the system state space is therefore non-finite.

The stability of linear time-invariant (LTI) systems with internal delays is studied in [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF]. The state of the delay is projected over a finite Legendre polynomial basis. These projections are time-dependent values since the state of the delay is also time-varying. The time derivative of these projections only depends on the input of the delay operator. Then the norm of the state is overapproximated using a Bessel inequality. By integrating this inequality, we get an Integral Quadratic Constraint (IQC) between the output of the delay operator, its

Paul Rousse

Set-based Cosimulation 4.5. Example input, the derivative of its inputs, the projections over the truncated basis of Legendre polynomial and an error signal. The IQC models the energy of the remaining of the Legendre expansion (i.e. the error signal). In [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF], the stability of the delayed LTI system is assessed for all possible error signals satisfying the derived IQC. We use this IQC to overapproximate the reachable set of such systems.

In what follows, we use the first order of the IQC relationship described in [Seuret and Gouaisbaut, 2015, Theorem 5]. The state ξ corresponds to the average value of the state of the delay. The remaining energy of the state is bounded by an integral quadratic constraint.

     ξ(t) = -15ξ(t) + 1.5v(t) -w(t) with ξ(0) = 0 under the IQC t 0 w(s) 2 ds ≤ t 0 0.0025 v(s) 2 -0.75 (v (s) -ξ (s)) 2 ds (4.
28)

The IQC system (4.28) is used to overapproximate the delay in the following system:

ẋ = -x -k c D h (x) + u x(0) = 0 (4.29)
where k c = 4, h = 0.01 and u(t) = 1t. (4.27,4.28,4.29) are then combined in a unique linear time-invariant system with an integral quadratic constraint.

             Ẋ(t) = AX + B w w(t) + B u u(t) X(0) = 0 0 t 0 w(τ ) 2 dτ ≤ t 0 X(τ ) u(τ ) ⊤ M X(τ ) u(τ ) (4.30)
where the matrices are defined by The bounds in Eq.( 4.24) are [w] = [-10, 10] and

A = 1.
[w ′ ] = [-1, 1]. The initial distur- bance set is defined such that [w 0 ] = [w].
Figure 4.8 corresponds to the reachable tube of the delayed system modeled with the integral quadratic constraint. Y IQC is the reachable tube of the corresponding system.

Set-based Cosimulation

Paul Rousse .8: Computation of the reachable tube of the system (4.30) using the validated numerical integration framework described in the previous section and the contractor C introduced in Section 4.4. y (the blue line) corresponds to the response of the delayed system. Y IQC is the exact reachable tube of system computed using the paraboloid method presented in Chapter 3.

Discussion One motivation of this work is to use Integral quadratic constraint (IQC) models in a validated numerical integration framework. IQC models are widely used in the robust control community for stability analysis of dynamical systems. When the IQC system is stable, there exists an invariant over the set of states (x, z), and the maximal reachable z value (i.e. the maximal integral value reachable) is bounded for any trajectory.

In our approach, such an invariant does not exist. The overapproximation of the maximal reachable z is constantly increasing in size. Consequently, bounds provided by the fixed point algorithm are also strictly increasing in size. When these bounds reach the prior bounds given by (4.18) over the disturbance, the reachable set tends to the reachable set computed without the integral constraint. Figure 4.9 corresponds to the reachable set of Example 4.5 for a larger horizon of integration. The integral constraints provide bounds over w. However, when the energy level is too high, these bounds are strictly included in the bounds given by (4.18). At t = 15s, the reachable set converges to the reachable set of the system with no integral constraint between the disturbance and the state.

The bounds of the input disturbance depend on the result of the used guaranteed set integration method. Therefore, if the later is too pessimistic, the proposed contraction method will only rely on the bounds [w] and [w ′ ] of Eq.( 4.18).

In our approach, a larger class of systems is considered compared to the linear case treated in Chapter 3. Contrary to IQC models, only the dependence in the disturbance needs to be quadratic for the integral constraint.

In terms of scalability, our approach needs the state of the original dynamical

Paul Rousse

Set-based Cosimulation

Conclusion

Figure 4.9: Computation of the reachable tube of the system (4.15) in Example 4.5 over [0, 20] using the validated numerical integration framework described in Chapter 4 and the contractor C introduced in Section 4.4. In blue, the reachable set when only (4.18) is used (i.e., when the integral constraint is not used). In green, the reachable set of the system when the integral constraint is taken into account.

system to be extended from n variables to p = n + m + 1 variables (m states for w, 1 state for z). Since the disturbance signal spans in a subspace of R n , m is always smaller than n. Since m is often close to 1 (the delay modeled as an integral quadratic constraint introduces a 1-dimensional disturbance signal), p is close to n (or 2n in the worst case). However, only the integration part can suffer from the dimension of the system. Based on the advantage of our approach, a less expansive integration method can be used for large systems for a similar result.

Conclusion

We presented a method to compute an overapproximation of the reachable tube for dynamical systems with integral constraints over the input set. To overapproximate the reachable set, we use a Runge-Kutta validated numerical integration scheme with pessimistic bounds over the input. It provides a first conservative bound over the reachable tube. Then, the integral constraint is used to define a contractor over the reachable tube. This contractor and a propagation step are successively applied to the overapproximation until a fixed point is reached. We evaluated our algorithm with DynIbex library to simulate a delayed system, i.e., an infinite-dimensional system that can be modeled as a linear time-invariant system subject to an integral quadratic constraint.

Set-based Cosimulation Paul Rousse

Future works The method developed in this chapter is guaranteed (we compute an overapproximation of the reachable tube). However, our overapproximations tend to constantly grow in size, even when the reachable set is known to be bounded (see Section 4.5). Such an issue originates from the use of intervals to overapproximate the reachable set. Intervals are too conservative and the use of affine forms to overapproximate the reachable set should provide better-conditioned overapproximations. One purpose of this work was to investigate how models introduced by the robust control community can be used in the field of validated numerical integration. We proposed the use of an IQC model that bounds the input to output L 2 gain of a system with an internal constant-delay. Many complex systems can be modeled in a similar approach. More specifically, the error of approximation in a reduced system can be expressed with an input to output L 2 gain constraint. Simplification of models is very appealing for validated numerical integration since the computational time is mainly dependent on the system dimension. Being able to reduce the order of the system and to bound the error with a 2-norm gain would lead to a more efficient algorithm.

Part I presents a set-based simulation method that uses time-varying ellipsoids to compute overapproximations of the reachable set of a system. The systems of interest are continuous-time linear time-varying systems subject to disturbances bounded by quadratic constraints in a Hilbert space. Two subclasses of models have been more specifically developed, the (point-wise) quadratic constraint model and the Integral Quadratic Constraint (IQC) case. These classes of systems are of particular interest since they are widely studied in the field of robust control theory. In particular, these models can be used as abstractions to complex models that do not fall into the initial scope of considered systems (e.g. nonlinear systems, systems with internal delays).

Many other reachability analysis frameworks are available. They differ by the family of systems they can analyze, the geometrical sets used to overapproximate the reachable set, and the available tools to compute these overapproximations. Complex systems might involve different family of dynamical systems (discrete-time and continuous-time, nonlinear and linear behaviors for example). In this case, one might need to use several reachability analysis frameworks to compute the overapproximation of the reachable set.

In this part, we propose to analyze an interconnection of systems. Each subsystem in the interconnection corresponds to a system that can be analyzed with its associated reachability analysis framework (such as the ones developed in Part I and Part II). This interconnection of systems can be described with two basic operations, a composition of systems and a feedback operator. One theoretical and practical difficulty in the analysis of such an interconnection is to "close the loop", i.e. to study the following system for every input u in a set of inputs U . Chapter 5 describes the classical framework of abstract interpretation and its application to study an interconnection of systems.

Set-based Cosimulation

Paul Rousse

This chapter presents the abstract interpretation framework applied to the analysis of an interconnection of systems. Such interconnection frequently appears in system modeling and usually involves systems of heterogeneous nature. For example, a continuous plant controlled by a discrete-time controller is an interconnection of two systems. Hereby, the behavior of the continuous plant might be described as solutions to an ordinary differential equation and the behavior controller by a computer program. Even if there exists tools to individually analyze each system, analyzing their interconnection is still challenging. In this chapter, we show that the abstract interpretation framework is suited to address such issues.

Abstract interpretation is a common tool in the verification of computer programs. To prove that a program satisfies a set of properties, abstract interpretation starts from these two following observations: a computer program can be exactly described by its concrete semantic, i.e. the set of its executions; and the specification to verify is, as well, a subset in the domain of the concrete semantic. If the concrete semantic is a subset of this property, then the program is valid. However, this concrete semantic is usually not computable. Abstract interpretation aims at computing an abstract semantic, a superset of this concrete semantic. The abstract semantic is also called sound approximation (or over-approximation) of the concrete semantic. Contrary to the concrete semantic, the abstract semantic is designed to be a simpler and hopefully a computable mathematical object. If the abstract semantic is a subset of the property, then the concrete semantic, i.e. the program, verifies the property.

Abstract interpretation has been mainly used in program analysis where the concrete semantic domain corresponds to the set of finite or infinite sequences of symbols. However, this notion of semantic is versatile enough and can be used in other fields. In the case of dynamical systems, the semantic domain corresponds to a signal space (such as the set of finite or infinite horizon, discrete or continuous-time, real-valued vectored signals). This chapter is organized as follows. In Section 5.1, we define a language for an interconnection of systems. In Section 5.2, we define the concrete semantic of an interconnection of systems. Section 5.2, Section 5.3 and Section 5.4 present the usual framework and classical results of Abstract Interpretation. Most of their content is available in [START_REF] Cousot | Systematic design of program analysis frameworks[END_REF] and is restated hereby for completeness. In Section 5.2, we define the concrete semantic that is equivalent to the semantic of the interconnection of systems. In Section 5.3, we define abstract domains. In Section 5.4, we define the abstract semantic. In Section 5.5, Section 5.6, and Section 5.7, we present abstracts domains that we use to describe signals. In Section 5.8, we apply the framework to find a sound approximation of the reachable tube of several interconnections of systems. In Section 5.9 and Section 5.10, we present the works related to the presented framework and conclude.

Notations are not classical and are chosen to be consistent with the rest of the document. Instead of the classical • # for abstract elements and • ♭ for concrete elements, we use a starred notation for abstract elements (x * , X * , . . . ) and no superscript notation (x, X, . . . ) for the concrete elements.

Paul Rousse

Set-based Cosimulation

We assume that any ISys expressed in the syntax described by Table 5.1 identify each connection (a serial connection or a feedback connection) with a unique label in the set of labels V (in other words, a label does not appear two times in the formula). An interconnection of systems ISys is then described with a set of systems, a set of sources, and connections in-between each of them.

Signal The syntax describing an interconnection of systems is used to define the connections in-between each system. A connection is labeled with a u ∈ V and associated with a signal u that is an element of a domain Dom u . These signals are chosen as a concatenation of constant (e.g. a parameter, an initial state), continuous-time (e.g. the output or input of a continuous-time dynamical system), and/or discrete-time signals (e.g. a signal generated by a computer program). Signals might be defined over a finite or infinite time-horizon. Then, Dom u corresponds to the Cartesian product of each domain of each type that composes u. All the signals of the interconnection of systems are then described with a single variable ρ, called an environment, ρ = (u, x, y, . . . ), V = {u, x, y, . . . }, that belongs to the Cartesian product of all the domains Dom = v∈V Dom v .

For an environment ρ ∈ Dom, a set of labels V, and a label u ∈ V, let ρ u be the projection of ρ over Dom u .

Source A source U ⊆ Dom u is a subset associated with a signal u ∈ Dom u , and a label u ∈ V. A source describes an exogenous input such as an initial state of a dynamical system (or of a program), or an unknown bounded input signal (e.g. U might be the set of signals from R + to R bounded at any time by [-1, 1]).

Sink The sink • is a terminal point used such that every signal in the interconnection is connected to another element.

System The systems are defined as operators between an input domain and an output domain. For example, the system S in " . . . u S y . . . ", where u, y ∈ V, defines a mapping S : Dom u → Dom y from Dom u (the domain of the signal u associated with the label u ∈ V) and Dom y (the domain of the signal y associated with the label y ∈ V). An interconnection of systems is then described with serial connections S 1 v S 2 , feedback connections µ x {S} and exogenous inputs in U . The connections between systems are assumed to be correctly defined such that every types are compatible. A serial connection

Paul Rousse

Set-based Cosimulation The resulting closed-loop system " . . . u µ x {S} y . . . " defines a system mapping Dom u to Dom y , Interconnection of systems An interconnection of systems S of type ISys is then described as a source Src connected to a sink Snk.

Semantic of the interconnection of systems For an interconnection of systems ISys denoted by S and expressed in the syntax defined in Table 5.1, the semantic of S corresponds to the set of environments [S ] ⊆ Dom defined by

[S ] = {ρ ∈ Dom | [S ](ρ; ∅, ∅) is True}
where ρ → [S ](ρ; ∅, ∅) is defined by Table 5.2. For each expression S of the syntax, Table 5.2 defines a function

[S] : Dom × V × V → {True, False}.

by structural induction over the syntax of S. The set V is defined by

V = {∅} ∪ V ∪ V 2 , (5.1)
V contains a symbol ∅, the labels in V, and the pairs of labels in V 2 . For a ρ ∈ Dom, each operator [S](ρ; u, y) defined in Table 5.2 evaluates the truth value that the signal ρ y is an output of the system S for the input ρ u . The symbol ∅ is associated with expression in the syntax which ends with a sink symbol or begin with a source. A Src denoted by U has only an output connection and is therefore associated with the function [U ](ρ; ∅, u). A Snk denoted by Y has only an input connection and is therefore associated with the function [Y](ρ; y, ∅). A system Sys denoted by S has an input and an output connection and is therefore associated with the function [S](ρ; u, y). The pair of labels are used in the concatenation " (U , S)" of a source U and a system S.

Set-based Cosimulation

Paul Rousse where S ol,x and S ol,y are projections of the system S ol over x and y (resp.). In such a case, the semantic of the system is defined as the set of environments ρ = (u, x, y) ∈ Dom such that (5.2a,5.2b,5.2c) are satisfied. The goal of this chapter is therefore to find the set of solutions to the fixed-point equation (5.2b) for all the inputs u in U . However, the fixed-point equation (5.2b) is difficult to solve in practice. This is especially true when Dom x is a combination of signals of heterogeneous types (e.g.

Paul Rousse

Set-based Cosimulation 5.1. Syntax and semantic of an interconnection of systems when x is a combination of discrete-time and continuous-time signals). The following section expresses the system semantic with operators lifted to sets. This semantic results as well in a fixed-point equation in a particular structured set: a partially ordered set. In such a structure, there exist practical ways to solve this fixed-point equation. The system can be equivalently described by " U u µ r {S ol } z • " where the block diagram of the open-loop system S ol is represented in Figure 5.2.

Remark 5.1. Construction of the semantic

For interconnections of systems as well as for computer programs, the formalization of the semantic of a syntax is not unique. [START_REF] Bouissou | A hybrid denotational semantics for hybrid systems[END_REF] considers an interconnection of a computer program and a dynamical system. The semantic of the computer program part is constructed by extending each trace with the next possible reachable state (by using the transition function). In between two transitions of the program, the dynamical system is an autonomous system and can just run toward the future until a new event happens. Then, the dynamical system can just be modeled as a discrete-time system, and the semantic of the interconnection is expressed in a similar way than for the computer program.

Set

Such an approach is possible because the continuous-time system (the dynamical system) and the discrete-time system (the computer program) defines equations that can be "decoupled" in time. In our case, we do not only consider interconnections of one discrete-time system with one continuous-time system but also interconnections of multiple instances of each type. Therefore, signals can interact with each other, and it is not possible to define a semantic that uses a forward propagation scheme for each systems. We define the semantic of the interconnection of systems as the set of signals that satisfies a set of equations. x(t) = x 0 + t 0 f (x(s), u(s))ds for all t ≥ 0 for a given initial state x 0 ∈ R nx , a given observation function h : R nu × R nx → R ny , and a given dynamical function f : R nx × R nu → R nx . Then, the signal x of the feedback connection of " . . . u µ x {S ol } y . . . " satisfies the fixed-point equation

x(t) = x 0 + t 0 f (x(s), u(s))ds
for any t ∈ I. When the solution x is differentiable, this fixed-point equation is equivalent to the initial value problem ẋ(t) = f (x(t), u(t)), for all t ∈ I x(0) = x 0 .

Paul Rousse

Set-based Cosimulation 5.2. Concrete domain and semantic

Concrete domain and semantic

Now that we have defined the syntax for an interconnection of systems and its associated semantic, we introduce the concrete domain and concrete semantic. The concrete semantic is defined to be equivalent to the system semantic except that its construction involves operators over sets of signals. The next section defines sufficient conditions to derive a sound approximation of this concrete semantic: the abstract semantic.

Extended interconnection of systems In order to solve the fixed-point equation with a set-based approach, we study an extended version of the interconnection of systems (see Remark 5.3). To do so, for each signal u in Dom u associated with the label u ∈ V, we introduce a new signal u ′ in Dom u ′ = Dom u associated with the new label

u ′ ∈ V ′ = {v ′ | v ∈ V}.
This extended version of the interconnection of systems redirects any output y associated with y ∈ V of each system to the new signal y ′ associated with y ′ ∈ V. Then, the solutions to the fixed-point equation are searched in the set of environments ρ satisfying ρ y = ρ y ′ for any label y ∈ V. This extended system has a set of labels V = V ∪ V ′ , and a domain Dom = Dom × Dom.

In the rest of this document, the interconnection of systems implicitly refers to this extended interconnection of systems. We keep the notation V and Dom to respectively refer to the set of labels (i.e. to V) and to the domain (i.e. to Dom). and its set of solutions X ⊆ Dom, where f : Dom → Dom. Such a fixed-point equation can be studied in a set-based approach. Let the lift to sets F of f be defined by F (X) = {f (x) | x ∈ X}.

(5.5)

The solutions to the fixed-point equation X = F (X) (5.6) are closely related to the fixed-points of f . Every subset X ⊆ X of fixed-points to (5.4) is a fixed-point of F . For every pair of fixed-points X, Y of F , X ∪ Y is as well a fixed-point of F . And therefore, if Z is a fixed-point of F , then Y = Z ∪ X is as well a fixed-point of F , i.e. there is a fixed-point Y to F larger (or equal) than X. Thus, (5.6) where λ = (x, x ′ ) ∈ Dom 2 and g is an extended version of f defined by

g(λ) = (x ′ , f (x)).
In fact, λ = (x, x ′ ) is a fixed-point of g iff x is a fixed-point of f (the proof is direct from (5.7) as x = x ′ = f (x)). Let Λ be the set of fixed-points of g.

Similarly than for f , the lift to sets G of g can be used to characterize Λ. Let

Σ x = {(x, x ′ ) ∈ Dom 2 |x = x ′ } and Λ be a fixed-point of G such that Λ ⊆ Σ x . For λ 1 = (x 1 , x ′ 1 ) ∈ Λ, since Λ is a fixed-point of G, there is a λ 2 ∈ Λ such that λ 1 = g(λ 2 ), i.e. x 1 = x ′ 2 x ′ 1 = f (x 2 ) Since λ 1 and λ 2 belong to Λ ⊆ Σ x , it holds x 1 = x ′
1 and x 2 = x ′ 2 and therefore,

x ′ 1 = f (x 1 ).
To summarize, when G(Λ) = Λ and Λ ⊆ Σ x , every x, such that λ = (x, x ′ ) ∈ Λ, is a fixed-point of f . The projection Λ| x of Λ is a subset of X, i.e. Λ| x ⊆ X. For this reason, Λ (and Λ) can be characterized by searching for the largest fixed-point of G lower than Σ x .

Example 5.2.

Let f be the polynomial defined by f (x) = x 2 -1. The fixed-points of f are the roots of the second degree equation f (x)x = 0. Thus, f has two fixedpoints X = {x 1 , x 2 }, where

x 1 = -1+ √ 5 2 and x 2 = -1- √ 5 2
. Similarly, the fixedpoints of f • f are the roots of the fourth degree equation f • f (x)x = 0. f • f has four fixed-points which are X = {-1, 0, x 1 , x 2 } (see Figure 5.3). Let the lift to sets of f be F (as defined in Remark 5.3). X is a fixed-point of F since F (X) = {0, -1, x 1 , x 2 } = X. In a more general way, all the sets defined by X = {x 1 , x 2 , . . . , x m } that are composed of elements of a periodic sequence {x n } n∈N (of period m ∈ N) defined by

x n+1 = f (x n ) are fixed-points of F (indeed F (X) = {f (x 1 ), f (x 2 ), . . . , f (x m )} = {x 2 , x 3 , . . . , x 1 } = X).

Concrete domain

The concrete domain A is defined as the powerset of environments, i.e. A = ℘(Dom). We equip this concrete domain with a complete lattice structure.

Paul Rousse

Set-based Cosimulation

∀x ∈ A, x ⊑ x (reflectivity) ∀x, y ∈ A, x ⊑ y ∧ y ⊑ x ⇒ x = y (antisymmetry) ∀x, y, z ∈ A, x ⊑ y ∧ y ⊑ z ⇒ x ⊑ z (transitivity) Definition 5.

Lattice

A lattice is a tuple (A, ⊑, ⊔, ⊓) with (A, ⊑) a poset, ⊔ a binary upperbound, i.e. ∀x,

y ∈ A, x ⊑ x ⊔ y ∧ y ⊑ x ⊔ y ∀x, y, z ∈ A, x ⊑ z ∧ y ⊑ z ⇒ x ⊔ y ⊑ z ⊓ a binary lowerbound, i.e. ∀x, y ∈ A, x ⊓ y ⊑ x ∧ x ⊓ y ⊑ y ∀x, y, z ∈ A, z ⊑ x ∧ z ⊑ y ⇒ z ⊑ x ⊔ y Definition 5.3. Complete lattice A complete lattice is a tuple (A, ⊑, ⊔)
where A a set, ⊑ a partial order relationship (reflective, antisymmetric and transitive), ⊔ a lowerbound, i.e., for all subset S of A, ∀a ∈ S, a ⊆ ⊔S ∀b ∈ S, (∀a ∈ S, a ⊑ b) ⇒ ⊔S ⊑ b

The concrete domain (P(Dom), ⊆, ∪, ∩, Dom, ∅) is a complete lattice.

Set-based Cosimulation Paul Rousse

Concrete semantic

Each construction in the syntax (defined in Section 5.1) can be associated with an operator over the concrete domain. The concrete semantic is then defined by structural induction over the syntax. We will see that this concrete semantic can be expressed as the greatest fixed-point of a monotonic operator in the concrete domain. In the first paragraph, we present known results for the computation of such a greatest fixedpoint. In the second paragraph, we define the concrete semantic of the interconnection of systems.

Operators and fixed-points in a complete lattice Let a complete lattice (A, ⊑ , ⊓, ⊤, ⊥) and an operator F : A → A. F is monotonic when F preserves the ordering relationship

, i.e. X ⊑ Y ⇒ F (X) ⊑ F (Y ) for every X, Y ∈ A. A fixpoint Y of F is an element that satisfies Y = F (Y ). Let gfp X {F } be the greatest fixpoint of F that is lower than X. Y is a pre-fixpoint of F if Y ⊑ F (Y ), and Y is a post-fixpoint of F if F (Y ) ⊑ Y .
The following theorem provides the existence of such a fixed-point.

Theorem 5.1. Knaster-Tarski

The set of fixed-points of a monotonic operator F in a complete lattice is a complete lattice. Moreover, the greatest fixed-point of F that is lower than X is

gfp X {F } = {Y | Y ⊑ X and Y ⊑ F (Y )}.
Proof. See [Tarski, 1955]. ♦ By Theorem 5.1, the set of fixed-points of F lower than X ∈ A is a complete lattice, there is a lowest and a greatest fixed-point. The monotonicity of the operator F can be used to define a decreasing sequence of Y k in A starting from a post-fixpoint of F . When the initial element Y 0 is greater than X, then each element is greater than gfp X {F }.

Proposition 5.1. Descending chains

If F is a monotonic operator and X is a post-fixpoint of F , then

gfp X {F } ⊑ F k+1 (X) ⊑ F k (X) ⊑ • • • ⊑ F (X) ⊑ X .
Proof. It is a direct consequence of the monotonicity of F applied to the post-fixpoint X of F . ♦

Paul Rousse

Set-based Cosimulation

Concrete domain and semantic

Concrete semantic

The concrete semantic is defined in Definition 5.4 as the greatest fixed-point of the concrete operator X → S (X; ∅, ∅) defined by Table 5.3 by structural induction over the syntax of an interconnection of systems S (as introduced in Table 5.1). Similarly than for the semantic of the interconnection of systems, for every expression S of the syntax of Table 5.1, we define an operator S (R; u, y) in the concrete semantic that associates to a set of environments R ⊆ Dom, an input label u, and an output label y, another set of environments in Dom:

S : ℘(Dom) × V × V → ℘(Dom)
where V is defined by

V = {∅} ∪ V ∪ V 2 . • (R; y, ∅) := R U b (R; ∅, u) := {ρ ∈ R | ρ u ′ ∈ U} S b (R; v, w) := {ρ[w ′ ← S b (ρ v ), v ← ρ v ′ ] | ρ ∈ R} S 1 v S 2 (R; u, w) := S 2 S 1 (R; u, v); v, w µ x {S} (R; u, y) := S (R; (u, x), (y, x)) (U , S) (R; u, (v, y)) := S (R; v, y) ∩ U (R; ∅, u)
Table 5.3: Concrete operator for an interconnection of systems. U b ⊆ Dom u , S b : Dom v → Dom w , S 1 and S 2 are either a Src and a Sys, two Sys, or a Sys and a Snk, S is a Sys, and U is a Src. u, v, w, x, and y are labels in V ∪ {∅}.

We now study the fixed-points lower that Σ ⊆ Dom of the concrete operator X → S (X) where Σ is defined as the set of environments where each signal is equal to its prime signal (see Remark 5.3)

, i.e. Σ = {ρ ∈ Dom | ∀v ∈ V, ρ v = ρ v ′ } .
(5.8)

Each environment ρ of the semantic is a fixed-point of any concrete operator since, by construction, ρ ∈ Σ, and since each proposition satisfied in Table 5.1 implies that {ρ} is a fixed-point of its corresponding operator in Table 5.2. More generally, any subset of the semantic of the interconnection is a fixed-point of the concrete semantic. Set-based Cosimulation Paul Rousse

Abstract domains

Expressing the semantic as the fixed-point has been initially introduced to study computer program [START_REF] Cousot | Systematic design of program analysis frameworks[END_REF]. In these works, the concrete domain is represented with the sets of finite and/or infinite traces, the fixed-point is iteratively constructed by extending each trace through time. Then, the fixed-point of interest is the lowest fixed-point that contains these trajectories. [START_REF] Bouissou | A hybrid denotational semantics for hybrid systems[END_REF] has extended this approach to the case of an interconnection of a controller (modeled as a computer program) and dynamical systems. The trajectories of the dynamical system are simulated over a small time-horizon.

As we pointed out in Remark 5.1, such an approach is possible because systems are be decoupled in time. In our case, the system trajectories are expressed as the solution to a fixed-point equation that depends over an input. For each input, there might be one (or more) solution to this fixed-point equation. The semantic of the system is then defined as the union of all these fixed-point solutions for every possible input. The system semantic can be equivalently defined with operators over sets. For these monotonic operators, the union of fixed-points corresponds to the greatest fixed-point (by the Theorem 5.1).

Abstract domains

The previous section defines the concrete semantic as the greatest fixed-point of an operator in the concrete domain. It is then possible to use an iterative method to overapproximate this fixed-point (see Proposition 5.1). However, elements of the concrete domain are too complex to be computer-represented and therefore, the concrete semantic cannot be calculated. In this section, we soundly approximate it with elements chosen in a subset of the concrete domain. This sound approximation of the concrete domain is called the abstract domain. The abstract domain is supplied with a complete lattice structure where elements and operators (resp.) can be represented and calculated (resp.) on a computer program. The abstract semantic is then deduced from the concrete semantic. Each operator of the abstract semantic mimics its associated operator in the concrete semantic.

The concrete domain and abstract domain are linked with a so-called Galois connection. This Galois connection enforces the soundness property. Each element of the concrete domain is associated with an abstract element that is a sound approximation (i.e. an overapproximation). The abstract semantic is derived such that the Galois connection is preserved through each evaluation in the concrete semantic.

To simplify notations, the complete lattice structure associated with the concrete domain (℘(Dom), ⊆, ∪, ∩) is denoted with (A, ⊑, ⊔, ⊓). Each element in A corre-

Set-based Cosimulation

Paul Rousse sponds to a set of signals of the interconnection. Abstracting the concrete semantic corresponds to describing this set of signals with elements in a A ⊆ A. In practice, elements of A are chosen to be computer representable (such as time-varying intervals or time-varying ellipsoids for example). By using an isomorphism from Ā to a set A * , each element P ∈ Ā is associated with an abstract element P * of an abstract domain A * .

Definition 5.5. Galois connection Let (A, ⊑, ⊔, ⊓) and (A * , ⊑ * , ⊔ * , ⊓ * ) be two complete lattices. A pair of functions (α, γ), with α : A → A * and γ :

A * → A, is a Galois connection if it holds ∀x ∈ A, ∀x * ∈ A * , α(x) ⊑ * x * ⇔ x ⊑ γ(x * )
We can enumerate a few properties of the Galois connection.

Theorem 5.2. Galois connection Let (A, ⊑, ⊔, ⊓) and (A * , ⊑ * , ⊔ * , ⊓ * ) be two complete lattices and (α, γ) be a Galois connection between A and A

* (A, ⊑, ⊔, ⊓) γ ← ----→ α (A * , ⊑ * , ⊔ * , ⊓ * ).
Then, the following properties are satisfied α is a monotonic function, γ is a monotonic function, ∀x * ∈ A * , α • γ(x * ) ⊑ * x * , and ∀x ∈ A, x ⊑ γ • α(x).

The abstraction function α associates to any concrete element x ∈ A an abstract element x * = α(x). The concretisation function γ associates to any abstract element x * ∈ A * a concrete element x = γ(x * ).

Sometimes the existence of a Galois connection between two complete lattices is too strong of a requirement as the abstraction function α might not exist. For such cases, [START_REF] Cousot | Abstract interpretation frameworks[END_REF] proposes to relax the Galois connection framework to work only with the concretisation function γ.

Definition 5.6. Concretisation function

Let (A, ⊑, ⊔, ⊓) and (A * , ⊑ * , ⊔ * , ⊓ * ) be two complete lattices. A concretisation function is monotonic function γ : A * → A. X * is an abstraction of X when X ⊆ γ(X * ).

Paul Rousse

Set-based Cosimulation

Abstract semantic

Partial narrowing operator In practice, when Y * 0 is an abstraction of gfp X {F }, it is always possible to refine Y * 0 with a sound approximation of Y → F (Y ) ⊓ Y . The produced sequence decreases more rapidly than the sequence derived by (5.9) and is still a sound approximation of gfp X {F }. An abstract operator △ * is a narrowing operator if it is a sound approximation of the operator ⊓ * and such that for any X * , Y * ∈ A * it holds X * △ * Y * ⊑ * Y * . Classical definitions of the narrowing operator usually require that the sequence {Y * k } generated by

Y * k+1 = F * (Y * k ) △ * Y * k Y * 0 = D *
is ultimately stationary in finite-time (see Definition 2.2.4 in [Miné, 2004]), such property ensure the termination of iterates.

(Y * ⊓ * X * ) ⊑ * (Y * △ * X * ), and (Y * △ * X * ) ⊏ ∼ * X * .

Abstract semantic

In Section 5.4.1 and Section 5.4.2, we defined the abstract counterparts of each operator of the concrete semantic (as defined in Section 5.2). The abstract semantic is then defined by structural induction overt the syntax (defined in Section 5.1). When there is a Galois connection between the concrete domain and when each operator used within the abstract semantic is sound, then the abstract semantic is sound with respect to the concrete semantic.

Abstractions for vectorial space of finite dimension

The previous sections detailed the classical framework of abstract interpretation and the concrete semantic of an interconnection of dynamical systems. The following two sections describe abstract domains used to abstract elements of the concrete domains. This section details abstractions for constant signals such as system parameters or an initial state of a system. The next section uses a point-wise lift of these abstract domains to represents time-varying signals.

Ellipsoidal domain

In this section, we define the domain of convex sets, this domain can be supplied with a complete lattice structure. Since every convex set can be described as an intersection of ellipsoids, we then define the abstract domain of ellipsoidal sets. Finally, we define a narrowing operator that allows to describe a sound approximation of convex sets with fewer ellipsoidal sets.

Let Conv n be the set of convex subsets of R n . (Conv n , ⊆) is a partially ordered set. Since an intersection of convex sets is convex, c 1 , c 2 ∈ Conv n implies that c 1 ∩ c 2 ∈ Conv n . A union c 1 ∪ c 2 of convex sets c 1 and c 2 is not convex in the general case. However, its convex hull hull(c 1 ∪ c 2 ) (where hull(X) is defined as the intersection of all the convex sets greater than X, hull(X) = X⊆C C) is convex and it soundly approximates the union ∪. Let ∪ c be the convex hull-union operator, i.e. c 1 ∪ c c 2 = hull(c 1 ∪ c 2 ). To this respect Proposition 5.6. Convex set domain

(Conv n , ⊆, ∪ c , ∩, ∅, R n ) is a complete lattice.
Proposition 5.7. Galois connection of the convex set domain The pair of functions α C = hull and γ C = id is a Galois connection between the two complete lattices

(℘(R n ), ⊆, ∪, ∩, ∅, R n ) and (Conv n , ⊆, ∪ c , ∩, ∅, R n ).
Proof. Since hull(X) is defined as the intersection of all the convex sets greater than X, the equivalence hull(X) ⊆ C ⇔ X ⊆ C holds for all X ⊆ R n and C ∈ Conv n . ♦ Any convex set C ∈ Conv n can be described as an intersection of ellipsoids, i.e. C = E∈E E for E ⊆ Elli n . But contrary to Conv n , ellipsoids in Elli n can be represented with a symmetric matrix in S (n+1)×(n+1) . Since matrices can be represented and manipulated in a computer program, the manipulation of convex sets as an intersection of ellipsoidal sets is easier.

Paul Rousse

Set-based Cosimulation The power set of ellipsoids ℘(Elli n ) can be provided with a complete lattice structure. For E 1 , E ∈ ∈ ℘(Elli n ), we define the following relationships

E 1 ⊑ E E 2 iff E 2 ⊆ E 1 , E 1 ⊔ E E 2 = E 1 ∩ E 2 , and E 1 ⊓ E E 2 = E 1 ∪ E 2 . Proposition 5.8. Ellipsoidal domain (℘(Elli n ), ⊑ E , ⊔ E , ⊓ E , ∅, R n ) is a complete lattice.
Proposition 5.9. Galois connection of ellipsoidal domain

The pair of functions α

E (C) = {E ∈ Elli n | C ⊆ E} and γ E (E) = E∈E E is a Galois connection between the two complete lattices (Conv n , ⊆, ∪ c , ∩, ∅, R n ) and (℘(Elli n ), ⊑ E , ⊔ E , ⊓ E , ∅, R n ).
The Galois connections can be composed, and therefore, there is a Galois connection between the complete lattice (℘ Figure 5.4).

(R n ), ⊆, ∪, ∩, ∅, R n ) and (℘(Elli n ), ⊑ E , ⊔ E , ⊓ E , ∅, R n ) (see
Narrowing operator The intersection ⊓ E introduces many new terms. Each time, we compute the intersection between two elements E 1 and E 2 of ℘(Elli n ), E 1 ∩ E 2 requires as many ellipsoids as the one describing E 1 plus the ones that describe E 2 . To avoid having abstract elements that accumulate terms because of the intersection ⊓ E , we introduce a narrowing operator that soundly approximates the intersection ∩. This narrowing operator associates an intersection of ellipsoids with a sound approximating

Set-based Cosimulation

Paul Rousse ellipsoid. Since there is an infinite number of ellipsoids that overapproximates an intersection of ellipsoids, and since the set of ellipsoids greater than this intersection (with respect to ⊆) is not a pointed ordered set (i.e. there is no least element), we choose the least element with respect to another partial order (such as the volume of the ellipsoid for example).

We will talk more specifically about the case where ellipsoids are ordered by a cost function J : Elli n → R. For example, J might be the volume of the ellipsoids. Then, we define the narrowing operator △ * J by

Y * △ * J X * = arg inf J(Z * ) s.t. Y * ∩ X * ⊆ Z * for Y * , X * , Z * ∈ ℘(Elli n ). By construction, Y * ∩ X * ⊆ Y * △ * J X * , moreover, since X * belongs to the feasible set of this optimization problem, J(Y * △ * J X * ) ≤ J(X * ). △ *
J is a narrowing operator with respect to the order ⊏ ∼ * J defined by Z * ⊏ ∼ * J X * iff J(Z * ) ≤ J(X * ).

Abstractions for time-varying signals

In this section, we detail an abstract domain used to describe continuous-time or discrete-time signals.

Let Dom x be the concrete domain of a time-varying signal x (continuous-or discretetime) of the interconnection of systems. The signal x is a function that associates to a time-domain T (that is a subset of R + ) a value in R n . We can abstract such signals with the domain of time-varying sets which corresponds to a point-wise lift of the complete lattice (℘(R n ), ⊆, ∪, ∩, ∅, R n ) over the time-domain T . Definition 5.9. Point-wise lifting If (A, ⊑, ⊔, ⊓, ⊥, ⊤) is a complete lattice and T is a set, then (T → A, ⊑, ⊔, ⊓, ⊥, ⊤) is also a complete lattice if we define the dotted operator by point-wise lifting:

(T → A, ⊑, ⊔, ⊓, ⊥, ⊤), X ⊑Y is defined by ∀t ∈ T , X(t) ⊑ Y (t), ( ˙ X )(t) is defined by {X(t) | X ∈ X }, ⊤(t) = ⊤, ( ˙ X )(t) is defined by {X(t) | X ∈ X }, ⊥(t) = ⊥.
Then, we define the Time-varying set domain as the complete lattice (T → A, ⊑, ⊔, ⊓, ⊥, ⊤). that corresponds to the point-wise lift of the complete lattice (A, ⊑, ⊔, ⊓, ⊥, ⊤). Time-varying set domains are an abstraction of time-varying signals domains.

Paul Rousse

Set-based Cosimulation where X k ∈ A x . By Definition 5.9, the PWL domain inherits from the complete lattice structure of (A, ⊆, ∪, ∩, ⊥, ⊤). We define a concretization function γ PWL . This concretization function associates to Z PWL ∈ A PWL , an element of A, and is defined .11) where

γ PWL (Z PWL ) = j∈I γ PWL (Z PWL , j) (5
γ PWL (Z PWL , j) =    k∈K j λ j,k • Z PWL (k) ∀k ∈ K j , λ j,k ≥ 0 and k∈K j λ j,k = 1    where k∈K j λ j,k • Z PWL (k) is a weighted Minkowski sum
Z = {(p, x) ∈ R 2 | x = sin(1.5p) 2 }, ( 5
.12)

Figure 5.7: The PWL abstraction Z PWL in (5.13) defines a sound approximation γ PWL (Z PWL ) of the set Z in (5.12).

with p ∈ P = [-1, 1]. P is partitioned with a regular mesh {p k } k∈K = {-1.0, -0.6, -0.2, 0.2, 0.6, 1.0}, where K = 1, . . . , 6. The set Z can be soundly approximated by γ PWL (Z PWL ) where Z PWL is defined by

Z PWL (k) → {p k } × X k (5.13) with X 1 = X 6 = [1.09, 0.96] X 2 = X 5 = [0.71, 0.55] X 3 = X 4 = [0.12, -0.03]
where X k ⊆ R (see Figure 5.7).

Operator evaluation Let a function

F : A → A. When F satisfies F k λ k Y k ⊆ k λ k F (Y k ) ( 5

Examples

In this section, we present two examples that explains the concepts introduced in the previous sections. In Section 5.8.1, we present a simple closed-loop LTI discretetime system analyzed through Interval Arithmetic (presented in Chapter 4), In Section 5.8.2, we present an application of the piecewise linear abstract domain.

Closed-loop of a discrete-time system

We study the interconnection of systems S (see its block diagram in Figure 5.8) over the discrete-time interval 5.16) with a = 0.9 and b = 0.02. The set of initial states belongs to

T d = {1, . . . , T } ⊆ N for T > 0 G X i y k x i w ′ k w k S ol S
G :          x t+1 = ax t + bw t x 1 = x i y t = x t w ′ t = x t (
X 1 = [-1, 1].
Thanks to its simple form (G is 1 state dimension discrete-time linear system), its associated lifted operator G can be exactly computed using the interval arithmetic abstract domain as introduced in Chapter 4. Its exact computation is given in Remark 5.5.

Remark 5.5. Exact fixed-point expression

For this specific system, the S represents a discrete-time system

     x t+1 = (a + b)x t x 1 = x i ∈ X i y t = x t
(5.17) Dom can be described as a 1 + 2T vectorial space where each element can be represented by a vector (x i , w 1 , w 2 , . . . , w T , y 1 , y 2 , . . . , y T ). Each trajectory {x t } of S is a sequence that linearly depends on the initial value x 1 = x i ∈ X i . For every k = 1, . . . , T , it holds x t = x i c k-1 , where c = a + b, and thus, the semantic of the interconnection of systems [S ] ⊆ Dom is then exactly equal to

[S ] = {x i • (1, 1, c, . . . , c T -1 , 1, c, . . . , c T -1 ) | x i ∈ X i }.
In what follows, we detail the computation of the concrete semantic S of the interconnection of systems S . The set of labels of S is V = {x i , w, y, x ′ i , w ′ , y ′ }, each label is associated with its respective domain Dom xi = Dom x ′ i = R, Dom l = l(T d ; R) for l ∈ {w, y, w ′ , y ′ }. The concrete domain A = ℘(Dom) corresponds to the power set of the domain of all internal signals. The concrete semantic S ⊆ Dom is then deduced by syntactic decomposition of the expression (5.15). It corresponds to the greatest

Paul Rousse

Set-based Cosimulation As the definition of the concretisation function γ suggests, for ρ = (x i , w, y, x ′ i , w ′ , y ′ ) ∈ Dom, the signals x i and x ′ i are represented with the interval γ Int (X * (1)) ⊆ R, and the signals w,w ′ , y, and y ′ are represented with the interval set γ (T d →Int) (X * ).

Let S be an expression appearing in the syntactic decomposition of S (S is expressed in the syntax given in Table 5.1). Since we are searching for a fixed-point lower than Σ, each operator X → S (X; ∅, ∅) in the concrete semantic can be soundly approximated by S Σ (X; u, y) = S (X; u, y) ∩ Σ. (5.18) where u, y ∈ V (where V is defined by (5.1)). A sound approximation S * Σ of S Σ can be evaluated within the abstract domain of interval arithmetic as defined in Section 4.1 of Chapter 4. The greatest fixed-point is computed using descending chains as defined in Proposition 5.1. Starting with a post-fixpoint X * 0 (k) = [-10, 10] for all t ∈ T d , we compute the sequence X * k defined by

X * k+1 = S * Σ (X * k ; ∅, ∅)
for few iterations. due to the form of G), these iterates are Y k = γ T d →Int (X * k ), for k = 0, 1, 2, . . . and we compare them to Y cl that corresponds to the exact semantic of the interconnection of systems (projected over w = w ′ = y = y ′ ). Since every function evaluation is sound, by Proposition 5.5, and if X * 0 is a post-fixpoint of gfp Σ * {X → S * Σ (X; ∅, ∅)}, then γ T d →Int (X * k ) are a sound approximation of [S ]. Also, the time-varying interval Y k is a sound approximation of y = y ′ and w = w ′ where ρ = (x i , w, y, x i , w, y) ∈ [S ]. Figure 5.10 compares the volume of the overapproximation with the volume of the exact reachable set of S (see Remark 5.5) for the 6 first iterates. After three iterations, the volume error is less than 1%. The error exponentially decreases with iterations and converges to a value close to the numerical precision of the floatingpoint arithmetic.

Paul Rousse

Set-based Cosimulation

Examples

The set X init of initial states is defined by (5.21) Let the open-loop system S ol be S ol : " (X init , id u ) (x0,u) G y f nl ".

X init = {(x 1,init , 0) | x 1,init ∈ [-1, 1]}.
(5.22)

The system S can equivalently be expressed with " µ u {S ol } ". Its domain is Dom = Dom x 0 ∈ X init y = G(x 0 , u) u = f nl (y)

Concrete semantic

The concrete domain A is the power set of Dom, i.e. A = ℘(Dom). Each element ρ ∈ X in a concrete element X ∈ A is an environment ρ = (x 0 , u, y, x ′ 0 , u ′ , y ′ ) (i.e. an environment of the system extended with the "primes variables"). The concrete semantic is then defined by induction over the syntax.

Abstract semantic

We study two possible abstractions of the interconnection of systems S .

The first abstraction (Y 0 , U 0 ) makes use of the time-varying ellipsoidal domain. Since y and u are 1-dimensional signals, the time-varying ellipsoidal domain coincides with the interval domain. The operator G is a linear time-invariant system and can be abstracted using the ellipsoidal method presented in Chapter 1. The operator f nl is a nonlinear system, we abstract it with U 0 as the set of signals with values belonging to [-1, 1].

The second abstraction (Y i , U i ) (for i in {5, 10, 20}) makes use of the piecewise linear template allowing to describe the relationship between the initial condition x 0 and the signals (u, y). The set of initial condition is partitioned using a uniform grid {x k 0 } k∈K i where I have i elements. To each initial condition, we associate a trajectory (y k , u k ). Then, the signals u and y are abstracted using the ellipsoidal domain, and, as above, the operator G is abstracted using the ellipsoidal method as presented in Chapter 1. The operator f nl is abstracted using a linear piecewise model between each centered trajectories (y k , u k ).

Figure 5.13 and Figure 5.14 show the overapproximation of signals y and u (resp.) with the ellipsoidal domain (resp. Y 0 and U 0 ) and with the piecewise linear abstraction (resp. Y i and U i for i in {5, 10, 20}). Figure 5.12 compares the volume of the overapproximations with the volume of the exact reachable set of Y.

Set-based Cosimulation

Paul Rousse

Figure 5.12: Comparison of the precision when using a PWL abstraction. Y i for i in {5, 10, 20} corresponds to overapproximations of Y using the PWL abstraction (i corresponds to the number of points in the mesh of the initial state), Y 0 corresponds to the ellipsoidal domain.

Thanks to the Galois connection between the concrete and abstract domain, we know that the abstract semantic is a sound approximation of the concrete semantic. Since Y spans over [-0.7, 0.7] (see Figure 5.12), using (5.20), U spans [-1, 1]. Therefore, the abstraction U 0 of U with the ellipsoidal domain cannot be smaller than U . And to this respect, the smallest fixed-point using an ellipsoidal domain is often strictly larger than the actual set of trajectories. The ellipsoidal fail to capture the dependencies between the signals u and y.

When using the PWL, abstractions Y i and U i better represent the set of trajectories of the system S (see Figure 5.13,Figure 5.14 and Figure 5.12).

Paul Rousse

Set-based Cosimulation

Paul Rousse

Set-based Cosimulation 5.9. Related works 5.9 Related works [START_REF] Oulamara | Abstract interpretation with higher-dimensional ellipsoids and conic extrapolation[END_REF] introduces the gauge domain to overapproximate the reachable tube of discrete-time systems (that are computer programs in this case).

The reachable tube is abstracted with ellipsoids and ellipsoidal cones oriented in the time axis direction. In this work, we consider a larger class of abstractions: time-varying subsets of the state space, examples treated make use of time-varying ellipsoids as well as intervals.

An efficient reachability analysis framework is often specialized for a very narrow family of models. However, modern automotive systems are a mixture of many systems that are modeled in different frameworks. Trying to find one common model is often cumbersome. In this work, we propose an approach combining different analysis frameworks to overapproximate the reachable set of an interconnection of complex systems. These frameworks are articulated around the abstract interpretation framework.

In our work, the interconnection of systems is expressed with two connections: the serial connection of systems and the feedback connection. Such constructions have been studied in the field of reachability analysis under the name of compositional methods. [START_REF] Eqtami | A quantitative approach on assume-guarantee contracts for safety of interconnected systems[END_REF] computes invariants for an interconnection of two systems. This invariant is expressed as a level set in the state space. Then each system of the interconnection is used to refine a prior overapproximation (i.e. postfixpoint) until a fixed-point is found. [START_REF] Chen | What is ellipsoidal modelling and how to use it for control and state estimation?[END_REF] applies the same reasoning to reachability analysis in the case of multi-robot applications. The reachable set of each robot is overapproximated using a prior overapproximation of the reachable set of other robots. The reachable tubes are computed using intervals propagated with Taylor models (using the tool Flow* described in [START_REF] Chen | Flow*: An analyzer for non-linear hybrid systems[END_REF]). Previously cited works use compositional methods in order to analyze higher dimensional systems. Our motivation is different, we use compositional methods to couple different reachability analysis frameworks. [START_REF] Bogomolov | Reachability analysis of linear hybrid systems via block decomposition[END_REF] expresses the reachable set of an hybrid system as a fixedpoint equation over sets of system trajectories. The hybrid system is expressed as an interconnection of a continuous system and a discrete system. Then, trajectories are abstracted as a union of convex sets and the reachable tube is expressed as the least fixed-point of a discrete-time operator and the continuous-time operator.

In our approach, we produce a decreasing sequence of sets that overapproximates the reachable tube of the interconnection of systems. It is possible to think of these iterates as a set-based version of the Picard iterates that arise in ordinary differential equation integration (see Remark 5.2). In the field of validated numerical integration, [START_REF] Moore | Introduction to Interval Analysis[END_REF] has introduced a set-based Picard operator evaluated by using

Set-based Cosimulation

Paul Rousse the interval arithmetic framework (as presented in Chapter 4).

Conclusion

In this chapter, we study an interconnection of systems formulated within an abstract interpretation approach. The interconnection of systems is described with a syntax that defines a set of systems (a mapping from an input signal to an output signal), a set of sources (i.e. a set of signal), and a set of connections between a pair of these two (a serial connection and a feedback connection). Each connection is associated with: a signal, that is a constant, or a time-varying, continuous-or discrete-time variable; and a relationship that this signal satisfies with other signals in the interconnection of systems. We show that the set of signals of the interconnection of systems is the greatest fixed-point of a monotonic operator in a complete lattice. Such a problem is suited to the abstract interpretation framework. We introduce the point-wise lift abstract domains in order to represent time-varying signals. The method is applied to two toy examples.

Future works This work can be extended in several ways. We have abstracted signals with time-varying sets. When studying an interconnection of systems, this approach neglects the correlation between two systems. The piecewise linear (PWL) abstraction allows to represent this correlation between two variables. However, we only applied the PWL abstraction between a constant signal and a time-varying signal. This use is easier since the constant variable belongs to a space of low dimensions for which it was straightforward to define a partition. The use of the PWL abstraction between two time-varying signals is more complex since time-varying signals belong to spaces of higher dimensions (of infinite dimension continuous-time signals and of large dimension for discrete-time signals). Future works could focus on an extension of the PWL abstraction for time-varying signals.

The precision of our method has been evaluated for specific systems where the actual reachable tube was computable. For more complex systems, the precision can be evaluated between each system by computing the distance between some existing output signals and the actual overapproximation of the trajectory transformer. Future work would focus on a more systematic approach to evaluate the precision of produced overapproximations.

Paul Rousse

Set-based Cosimulation

Perspectives Conclusion

During this thesis, we developed three methods to overapproximate the reachable set of an embedded system.

Our first contribution concerns results around overapproximations with timevarying conics of the reachable set of a Linear Time-Varying (LTV) system subject to bounded disturbances. The time-varying coefficient of this conic is the solution to an initial value problem. Contrary to previous works for such a class of systems, we do not restrict the overapproximations to time-varying ellipsoids (i.e. overapproximations might be unbounded), and we work with homogeneous coordinates. In this coordinate system, the initial value problem has an elegant expression which is a Differential Riccati Equation (DRE). Such expression is interesting as the DRE has been heavily studied. We more specifically studied two subclasses of systems: LTV system subjects to Quadratic Constrained (QC) disturbances, and LTV systems subjects to Integral Quadratic Constrained (IQC) disturbances. For these subclasses, we proved that it exists an overapproximation over any time-horizon. Such a guaranty of existence has been rarely addressed. Usually, the existence of an overapproximation is either conditioned by the non-emptiness of a feasible set of a linear optimization problem or is either conditioned by the existence of a non-diverging solution to a DRE.

For QC systems, we proposed the use of positive time-varying multipliers. Each multiplier is associated with a time-varying conic overapproximation. Since the set of multipliers is infinite, the set of time-varying conics is as well infinite. We provide two methods to choose a multiplier that produce a tight time-varying conic overapproximation. In the first method, the multiplier is chosen such that the boundary of the overapproximation touches the reachable set along a so-called touching trajectory. In the second method, the multiplier is chosen such that the overapproximation has a minimal volume at a given time. Suboptimal solutions of this minimization problem are expressed with an application of the Pontryagin's Maximum Principle and solved with a continuation algorithm. These two methods were already developed for the case of systems with disturbances bounded by an exogenous signal, we extended it for the more general case of QC disturbances.

For the IQC systems, we introduced the extended system. This extended system embeds the integral quadratic into a new state of the system and the IQC is then expressed as a constraint over this new state. We then compute overapproximations of the reachable set with time-varying paraboloids (which corresponds to an extension of the previously introduced time-varying conics). For this extended system, it is possible to define touching time-varying paraboloid overapproximations and their associated touching trajectories. As for the QC system, the boundary of the touching time-varying paraboloids stays in contact with the reachable set of the extended system on the touching trajectory. We prove that the intersection of all the time-varying paraboloid overapproximations (generated by all the multipliers) exactly describes the reachable set of the extended system. This result is interesting as it shows that the set of multipliers is sufficient to exactly describe this reachable set. The existing literature in IQC systems never addressed the precision of overapproximations in such a way. Usually, since the IQC framework originated from stability analysis, multipliers where chosen in frequency domains, and touching trajectories and touching time-varying paraboloids could not be defined. The precision of the overapproximations was only appreciated on a practical example by comparing the different obtained overapproximations. By choosing multipliers in a temporal domain, we can define these touching trajectories and touching time-varying paraboloids. We developed and implemented an algorithm to compute a tight overapproximation of this reachable set. The set of multipliers is chosen according to the behavior of the touching trajectories that violate the state constraint in the future.

Possibles extensions of these works include the use of semidefinite solvers to derive overapproximations (optimal in volume for example). Such methods have been already developed for multipliers chosen in a frequency basis and can be adapted to the temporal domain multipliers used in this thesis. Optimal time-varying conics can as well be improved. We only considered optimal time-varying conic that minimizes the volume at a given time. A natural extension of our method would be to compute time-varying overapproximation minimizing the volume at any time. For such a case, the coefficient of the overapproximation satisfies a partial differential equation that involves two time-indexes, the regular time-index and the time-index at which the time-varying conic is of minimal volume.

Our second contribution is an extension of the validated numerical integration framework based upon interval arithmetic to the analysis of nonlinear systems subject to a disturbance bounded by an integral constraint. Validated numerical integration frameworks provide methods to compute a valid tube that contains the solution to an

Paul Rousse

Set-based Cosimulation

Conclusion

initial value problem. These methods can be used as well to compute overapproximations of the reachable set of a dynamical system subject to an unknown but bounded input. Previous works in this field studied the case of input with a point-wise bound (i.e. bounded at any time), we studied the case of an input bounded by an integral constraint. Such integral models are interesting as they can model complex systems (such as systems with inner delay). To do so, we use the integral constraint to define a contractor over the set of system trajectories. The contractor associates with a set of trajectory, a subset of these trajectories that satisfies the integral constraint. We then use this contractor in a fixed point algorithm to refine a prior overapproximation of the reachable tube.

A natural extension would explore the use of other geometrical shapes to overapproximate the reachable set. In particular, affine forms can be appropriate to express the relationship between the integral value and the state of the system (this relationship is not handled by the arithmetic interval framework and introduces some conservatism). Furthermore, in our work, we assume that the unknown disturbance evolves according to a model, such a hypothesis allows us to define the contractor out of the integral constraint. Future works should consider to remove such an assumption.

Our last contribution is a framework to study interconnections of systems. This framework is formulated within an abstract interpretation approach.

The interconnection of systems is described by a set of systems (which are mappings from a signal space to another signal space), a set of sources (a set of input signals), and a set of connections between these two. A connection is associated with a signal that corresponds to a constant, or time-varying (discrete-time or continuoustime) variable. Each connection defines a relationship between signals of the interconnection. The problem of interest is then to overapproximate the set of signals satisfying these equations. We show that these equations can be expressed as a fixed point equation over the signal space and that the set of signals of the interconnection of systems is the greatest fixed point of the lift to sets of this fixed point equation. Contrary to other approaches in this field, this set of signals is not defined as a sequence that evolves according to a transition function (as it is usually done in analysis of computer programs), a dynamic function (as for dynamical systems), or a combination of these two (as for hybrid systems). The semantic cannot be iteratively constructed, and furthermore, the concrete semantic is not defined as the least fixed point of a monotonic operator.

To overapproximate the reachable set of signals, we compute the greatest fixed point in an abstract domain. We propose to abstract time-varying signals with a point-wise lifting of an interval domain, or of an ellipsoidal domain, and to abstract each system with an overapproximating method (such as the two first methods pro-

Set-based Cosimulation

Paul Rousse posed in this thesis). Then, we overapproximate the greatest fixed point in the abstract domain with a refinement approach, i.e. starting from a prior overapproximation of the greatest fixed point and by refining it until a fixed point is reached.

To enhance the precision of our overapproximations, we introduce a so-called piecewise linear (PWL) abstract domain which corresponds to abstractions defined by a weighted sum (in the sense of the Minkowski sum) at chosen interpolating points. This PWL abstraction improved the analysis of systems where signals were correlated. Extensions of this work could investigate new abstract domains for time-varying signals. Since abstract interpretation has been introduced to study computer programs, there are not so many options to abstract continuous-time signals. We used a point-wise lift of the subsets of the state space. One could think to provide a frequency-based abstraction where signals are characterized by there frequency spectrum. Also, we used the PWL abstraction to represents the relationship between a scalar and a time-varying signal. A similar approach to model the relationship between two time-varying signals is necessary to improve the conservatism of overapproximations.

Introduction 1 1

 11 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Outline and Contribution . . . . . . . . . . . . . . . . . . . . . . . . Dynamical systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Set-based simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Conic sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Overapproximation with time-varying conics . . . . . . . . . . . . . . 1.5 Domain of definition of the overapproximations . . . . . . . . . . . . 1.5.1 Coefficient expansion of the DRE . . . . . . . . . . . . . . . . 1.5.2 Domain of definition of the time-varying ellipsoids . . . . . . . 1.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.

3 . 4

 34 Extended system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Paraboloids . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Overapproximation with time-varying paraboloids . . . . . . . . . . . 3.6 Exact reachable set . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.1 Overapproximation with an intersection of time-varying paraboloids 68 3.6.2 Overapproximation relationship . . . . . . . . . . . . . . . . . 3.6.3 Past trajectory for states in the interior of the overapproximation 3.6.4 Past trajectory for states on the boundary of the overapproximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.5 Exact reachable set . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.1 Examples from COMPleib . . . . . . . . . . . . . . . . . . . . 3.8.2 System verification . . . . . . . . . . . . . . . . . . . . . . . . 3.8.3 Delayed system . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices 3.A Continuous extension of the domain of definition of the time-varying paraboloids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  systems . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Set-based simulation . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 Conic sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 Overapproximation with time-varying conics . . . . . . . . . . . 16 1.5 Domain of definition of the overapproximations . . . . . . . . . . 19 1.5.1 Coefficient expansion of the DRE . . . . . . . . . . . . . . 19 1.5.2 Domain of definition of the time-varying ellipsoids . . . . . 20 1.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.6.1 Related works . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.6.2 Conservatism . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

  Figure 1.2: When E (see the block decomposition of P in (1.9)) is negative definite (subfigure 1.2a) and P = ∅, P is an ellipsoidal set, to this respect P is a bounded and convex set. When E is not sign-definite (subfigure 1.2b), P and P = ∅, the P is unbounded.

Proposition 1 . 4 .

 14 Boundedness of an upper and lower bounded matrixLet A, B, C ∈ S n×n , if A B C then B -B ≤ A -C where B = A+C 2 .

Figure 2 . 5 :

 25 Figure 2.5: Ellipsoids overapproximating a set of points: minimum volume ellipsoid (in red, corresponds to the objective a)), minimum sum of square semi-axis (in green, corresponds to the objective c)).

Figure 2 . 6 :

 26 Figure 2.6: In the continuation method, we follow the curve Q over [0, T f ] by using a prediction-correction algorithm.

Figure 2 . 7 :

 27 Figure 2.7: Reachable set of the system of Definition 2.1 with parameters defined in(2.21). The gray area corresponds to the minimum volume ellipsoid. Since the Pontryagin's Maximum Principle only gives necessary conditions for optimality, and since the ellipsoidal method only provides overapproximations, the E vol is not the actual minimum volume ellipsoid overapproximating the reachable tube R(T ; E 0 ). The dotted line represents an overapproximating hyperbola touching the reachable set.

3. 1 .

 1 Definition of the system that M w (•) ≺ 0 with the block decomposition of M M

( a )

 a Unit energy exogenous input. (b) Unit H ∞ norm feedback.

Figure 3 . 1 :

 31 Figure 3.1: Two trajectories of the unit energy response system and the unit H ∞ feedback loop. The second trajectory (dashed red line) is delayed, with a delay of τ = 4, compared to the first trajectory (plain blue line). This is true for any delay τ . Therefore, the state of such a system can be steered away from its equilibrium position at any time in the future.

Figure 3 . 2 :

 32 Figure 3.2: An infinite sum of functions corresponds to the set of positive and decreasing functions D * .

Figure 3 . 3 :

 33 Figure 3.3: Two overapproximations of the reachable set for the system of Example 3.1. Depending on the values of λ and µ, the overapproximation is bounded over the interval of integration (P s in blue) or is diverges in finite-time (P u in red).

Figure 3 . 4 : 2 Figure 3

 3423 Figure 3.4: Convergence analysis of the DRE for Example 3.2

3. 5 .Figure 3 . 7 :

 537 Figure 3.7: Stabilized constraint P s (t) versus finite escape time constraints P ns (t).The time-varying paraboloid P ns is defined over [0, 3] whereas P s is defined over R + .

Figure 3 . 8 :

 38 Figure 3.8: Convergence analysis of the DRE (1.21) for the unstable system defined in Example 3.3 for different scaling functions γ(•) = κ, with κ ∈ {0, 0.5, 1}. Stable equilibriums are marked with circles, unstable equilibrium are marked with stars. For κ = 0, the DRE does not have any stable equilibrium, and solutions have a finite escape time. When κ ∈ {0.5, 1}, the DRE has a stable equilibrium, the solution E with initial condition E 0 converges to this equilibrium and it is defined over R + .

Figure 3 . 9 :

 39 Figure 3.9: Domain of convergence of the solutions to the DRE (1.21) for the unstable system defined in Example 3.3 for different scaling functions γ(•) = κ, and different initial conditions E 0 . The red area corresponds to solutions with finite escape time, the yellow area corresponds to solutions converging to a stable equilibrium of the DRE.

Proposition 3 .

 3 7. First-order expansion of the value function For γ ∈ L 2,loc (R + , R + ) and λ 0 ≥ 1, let the corresponding time-varying paraboloid P = T ( P 0 , λ 0 , γ). For any t in the open set of I(P ), it exists ǫ > 0 and H > 0 s.t. for any δ ∈ L 2 (R + ; R), δ ≤ ǫ, for any δ 0 ∈ R, |δ 0 | ≤ ǫ, P = T ( P 0 , λ0 , γ) where γ = γ + δ and λ0 = λ 0 + δ 0 s.t. t belongs to the open set of I( P ), let ht be the value function of P (t) and z * = (x * , x * q ) an optimal trajectory of P , it holds ht (z * (t))β(t) ≤ Hǫ 2

ts

  (γ (r)+δ(r))dr ds ≤ Hǫ 2 with H = tRK 2 N ( n(0) 2 + n 2 ) (3.22) a finite constant where R = BM -1 w B ⊤ , K = exp t 0 L(τ )dτ and N = t 0 e t s (γ (r)+δ(r))dr ds. This ends the proof. ♦ Proposition 3.8 gives conditions where the sign of ht (z * (t)) is only determined by its first-order approximation defined in Proposition 3.7. Proposition 3.8. Sign of the value function of perturbed trajectories Let z * be a touching trajectory of P = T ( P 0 , λ 0 , γ) for γ ∈ L 2,loc (R + ; R + ), λ 0 ≥ 1 given and t ∈ I(P ) given. If there is a δ ∈ L 2,loc (R + ; R) and a δ 0 ∈ R, s.t. δ ≤ ǫ and |δ 0 | ≤ ǫ and t ∈ I( P ) (where P = T ( P 0 , λ 0 + δ 0 , γ + δ)) and r)+δ(r))dr ds is equal to the sign of ht (z * (t)) where ht is the value function of P (t) and H > 0 defined in (3.22) and

Lemma 1 .

 1 If Assumption 3.1 and Assumption 3.2 hold, every state z t ∈ ∂ Π(t) has a past touching trajectory z * of Π s.t. z * (t) = z t .

3 5 Sim 8 AddAlgorithm 3 . 2 :

 35832 .23) is used. Algorithm 3.1 summarizes the computation of Π* . An implementation on Matlab is available online[Rousse, 2019]. Algorithm 3.1: Computation of Π * defined by (3.24), in Section 3.7, as the subset of Π * defined by (3.14), in Section 3.6. Inputs : -A paraboloid P 0 ∈ P of initial of states -A time-horizon of simulation T > 0 -Sample time T c > 0 of constraint addition -The maximum N P cardinal of Π * Sim Parab = {( P 0 , 0)} 4 while t < T do /* Find the new time-varying paraboloids to consider and add them to Sim Parab (see Algorithm 3.2) */ Parab = U pdate Sim P arab(Sim Parab, Π) /* Simulate the paraboloid for T c */ 6 for (P τ , τ ) ∈ Sim Parab do 7 Simulate P(•) over [t, t + T c ] with P(τ ) = P τ Computation of the new paraboloids New Parab and update the set of paraboloids Sim Parab. Parameters: -Searching directions Search Dir ⊂ R n to add constraints -The maximum N new paraboloids to add at each step Inputs : -The current overapproximation Π * t -The updated set of paraboloids Sim Parab Output : -The updated set of paraboloids to simulateSim Parab Algorithm : 1 New Parab = {} 2 for n ∈ Search Dir do 3 project x c on ∂Π t in the direction n 4 let x * be this projection and P * ∈ Π * its corresponding touching paraboloid 5 compute λ given (x * ,P * ) 6 for λ = 1 + dλ, 1 + 2dλ, . . . , λ do 7 add (P * , λ) to New Parab 8 Sort New Parab according to the values of λ 9 Keep N new elements of New Parab with the largest value of λ 10 for (P * , λ) ∈ New Parab do 11 add (λP * (t), t) to Sim Parab

Figure 3 .

 3 Figure 3.13: Overapproximation of the output reachable set (projection of the reachable set R(t) through the observation map; the red area) of the AC10 example from the COMPleib library. The plain black line corresponds to the unperturbed trajectory of the system.

  Figure 3.14: The green surface in (a) is the reachable set R( P 0 , t) at t = 0.794 of S( P 0 , t) computed using Theorem 3.1. Its projection over the LTI state space (x 1 , x 2 ) (in solid red line) is shown in (b), each green line corresponds to one constraint P ∈ Π * computed with Theorem 1.1. (c) is the reachable tube t → R( P 0 , t) of S( P 0 , t) projected over the LTI state space (x 1 , x 2 ) for t ∈ [0, 1]. The red section corresponds to the time t = 0.794.

Figure 3 . 15 :

 315 Figure 3.15: Reachable set computation of the delayed system

3 .

 3 A. Continuous extension of the domain of definition of the time-varying paraboloids Appendices of Chapter 3 3.A Continuous extension of the domain of definition of the time-varying paraboloids This part introduces an intermediate result that is used for the proof of Proposition 3.7 in Section 3.6. By Definition 3.4, the domain of a time-varying paraboloid P is the domain of its quadratic time-varying coefficient E. Since the solution of the DRE (1.21) might diverge in a finite-time T E (E 0 ) < ∞ (where E 0 is the initial condition of (1.21)), P is defined only in the right-open interval [0, T E (E 0 )[. In this part, we show that since the touching trajectories of P are defined over the closed interval [0, T E (E 0 )], i.e. the definition of P can be prolonged to the same closed interval.

Contents 4 . 1

 41 Interval arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.1.1 System of equations . . . . . . . . . . . . . . . . . . . . . . 97 4.2 Validated numerical integration methods . . . . . . . . . . . . . 102 4.3 System with integral constraint over the state . . . . . . . . . . 105 4.4 Validated numerical integration for dynamical systems subject to integral constraints . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.4.1 Extended system . . . . . . . . . . . . . . . . . . . . . . . 107 4.4.2 Bounds over the disturbance . . . . . . . . . . . . . . . . . 108 4.4.3 Integral constraint propagation . . . . . . . . . . . . . . . 111 4.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

4. 1 .

 1 Interval arithmetic A few examples of arithmetic operations between interval values are given [-2, 5] + [-8, 12] = [-10, 17] [-10, 17] -[-8, 12] = [-10, 17] + [-12, 8] = [-22

5 )

 5 uses a Newton gradient descent to improve the contractive properties of [f ]. Example 4.2.Let f be a univariate function f : R → R defined byf (x) = 0.5(2x).(4.The solution x * to the fixed point equation (4.1) is x * = 2 3 . We compute the fixed-point iterates [x k ] defined by (4.4) starting from [x 0 ] = [0.5, 1]. The first iterates are (see Figure 4.1) [x 0 ] = [0.500, 1.000] = 0.750 ± 5.00e -01 [x 1 ] = [0.500, 0.750] = 0.625 ± 2.50e -01 [x 2 ] = [0.625, 0.750] = 0.688 ± 1.25e -01 [x 3 ] = [0.625, 0.688] = 0.656 ± 6.25e -02 [x 4 ] = [0.656, 0.688] = 0.672 ± 3.12e -02

Figure 4 . 1 :

 41 Figure 4.1: Sequence of [x k ] defined by (4.4) overapproximating the fixed point solution x * of (4.1) with f defined in (4.5) for [x 0 ] = [0.5, 1].

Figure 4 . 2 :

 42 Figure 4.2: The red area corresponds to the set of points (x, f (x)) where f (x) ≥ 0. The blue area corresponds to the set of points under the linear approximation of the concave function f at the boundaries of [x].

Figure 4 . 3 :

 43 Figure 4.3: The red area corresponds to the set of points x where f (x) ≥ 0..

4. 3 .Figure 4 . 5 :

 345 Figure 4.5: Results of the validated simulation for the system described in (4.10) with parameters (4.11). The system trajectories are overapproximated by the union of intervals [ y k ] (in grey). Two trajectories are represented (in red): one with y 0 = 0.4 and w(t) = cos(t); and another one with y 0 = 2 and w(t) = 1.

Figure 4 . 6 :

 46 Figure 4.6: Results of the validated simulation for the system described in (4.10) with parameters (4.11). The system has a unique trajectory (in red) that is overapproximated by the union of intervals [ y k ] (in grey).

4 . 2

 42 Validated numerical integration for dynamical systems subject to integral constraintsLet w be defined for any ǫ > 0 by  (τ )dτ = 1, the inequality in(4.15

4. 4 .

 4 Dynamical systems subject to integral constraints When w satisfies (4.18) and a given integral constraint, hard bounds (i.e. in ∞-norm) can be derived over w: Proposition 4.2. Disturbance bounds For a w ∈ L 2,loc ([0, h]; R m ) defined over an interval of length h > 0. If w satisfies (4.18) (with given bounds [w], [w ′ ] ∈ IR m ), then for any r > 0:

Figure 4 . 7 :

 47 Figure 4.7: Computation of the overapproximation of the reachable set of Example 4.5 using the algorithm presented in Section 4.4. Blue boxes correspond to the a priori enclosures at the first iteration of the algorithm Y 0 , green boxes are the a priori enclosure at the 3 rd iteration Y 3 of the algorithm.

Figure 4

 4 Figure 4.8: Computation of the reachable tube of the system (4.30) using the validated numerical integration framework described in the previous section and the contractor C introduced in Section 4.4. y (the blue line) corresponds to the response of the delayed system. Y IQC is the exact reachable tube of system computed using the paraboloid method presented in Chapter 3.

Figure 4 . 10 :

 410 Figure 4.10: Block diagram of a closed-loop system

125 5 . 1 .

 51 Syntax and semantic of an interconnection of systems " . . . u S 1 v S 2 w . . . " involves two systems S 1 : Dom u → Dom v and S 2 : Dom v → Dom w (i.e. the output set Dom v of S 1 corresponds to the input set of S 2 ). A feedback connection " . . . u µ x {S} y . . . " involves an open-loop system S : Dom u ×Dom x → Dom y × Dom x where the input and output signals both contain the same state signal x associated with the label x ∈ V.

[

  •](ρ; y, ∅) := True [U b ](ρ; ∅, u) := ρ u ∈ U b [S b ](ρ; v, w) := ρ w = S b (ρ v ) [S 1 v S 2 ](ρ; u, w) := [S 1 ](ρ; u, v) ∧ [S 2 ](ρ; v, w) [µ x {S}](ρ; u, y) := [S](ρ; (u, x), (y, x)) [(U , S)](ρ; u, (v, y)) := [S](ρ; v, y) ∧ [U ](ρ; ∅, u) Table 5.2: An environment ρ belongs to the semantic [S ] of an interconnection of systems S whenever [S ](ρ; ∅, ∅) is true. [S ](ρ; ∅, ∅) is computed by structural induction over the syntax of S . U b ⊆ Dom u , S b : Dom v → Dom w , S 1 and S 2 are either a Src and a Sys, two Sys, or a Sys and a Snk, S is a Sys, and u, v, w, x, and y are labels in V ∪ {∅}.In practice, the semantic on the interconnection of systems defines a set of equations over the signals. These equations are of three kinds: the inclusion (defined by the sources " U b u . . . ", i.e. u ∈ U b ), the mapping (defined by the system operator, " . . . u S b y . . . ", i.e. y = S ol (u)), and a fixed-point equation (defined by the feedback connection, " . . . u µ x {S} y . . . " i.e. (y, x) = S(u, x)). For any given expression, it is possible to compose the systems, and concatenate the signals such that the semantic is defined by " U u µ x {S ol } y • " (seeExample 5.1), which is equivalent to the following set of equations x = S ol,x (u, x) y = S ol,y (u, x)

Figure 5 . 1 :

 51 Figure 5.1: Block diagram of the system described by the formula (5.3).

Figure 5 . 2 :

 52 Figure 5.2: Open-loop system S ol of (5.3).

Remark 5 . 2 .

 52 Initial value problem as fixed-point equation For V = {x, u, y}, Dom x = L 2 (I; R nx ), Dom u = L 2 (I; R nu ), Dom y = L 2 (I; R ny ), and I = [0, T ] with T > 0. Let S ol be the system that associates to a time-varying input (u, x) ∈ Dom u × Dom y , the output (y, x) ∈ Dom y × Dom x where    y(t) = h(x(t), u(t))

Remark 5 . 3 .

 53 Lift to sets of a fixed-point equation Consider the fixed-point equation x = f (x)(5.4)

5. 2 .Figure 5 . 3 :

 253 Figure 5.3: Fixed-points of f and f • f in Example 5.2.

Proposition 5 . 2 .

 52 Sets of trajectories are fixed-points of the semantic Every R ⊆ [S] is a fixed-point of X → S (X). Proposition 5.3. Fixed-point of the concrete semantic Every fixed-point R of X → S (X) is a subset of [S ].

Proposition 5 . 5 .

 55 Soundness of the abstract semanticProvided the abstract operator are sound, the abstract semantic is sound with respect to the concrete one: S ⊑ γ S * .

5. 5 .Figure 5 . 4 :

 554 Figure 5.4: The set X ⊆ R 2 is soundly approximated by its convex envelop C = α C (X) = hull(X). The set C is equivalently described with C = γ E (E) the intersection of all ellipsoids overapproximating C (such as E 1 , E 2 ∈ E).

5. 7 .

 7 Figure 5.6: A partition of P = [-1, 1] 2 . P i is a cell, Pj is the neighbor cells of the vertex p k .

  .14) then the function F PWL defined by W PWL = F PWL (Z PWL ) Paul Rousse Set-based Cosimulation 5.8. Examples with W PWL (k) = F • Z PWL (k) is a sound approximation of F . The condition (5.14) is held for a particular class of functions. Let f : Dom → Dom be an affine function over the vectorial space Dom = Dom p × Dom x . Let F the lift to sets of f , i.e. F (A) = {f (a) | a ∈ A}. Since f is affine, f (λa + µb) = λf (a) + µf (b). For λ, µ ≥ 0 s.t. λ + µ = 1, the Minkowski sum λA + µB is well defined. For x ∈ F (λA + µB), there is a a ∈ A and b ∈ B such that x = f (λa + µb) = λf (a) + µf (b) ∈ λF (A) + µF (B). Therefore, F (λA + µB) = λF (A) + µF (B).

Figure 5 . 8 :

 58 Figure 5.8: Block diagram of the closed-loop discrete-time LTI system.

  of a monotonic operator X → S (X; ∅, ∅). Its detailed computation isS = gfp Σ {X → S (X; ∅, ∅)} S (X; ∅, ∅) = • ( * ; y, ∅) • µ w {G} ( * ; x i , y) • X i (X; ∅, x i ) X i (X; ∅, x i ) = {ρ ∈ X | ρ x ′ i ∈ X i } G (X; (x i , w), (y, w)) = {ρ[(x ′ i , w ′ ) ← G(ρ y , ρ w ), (x i , w) ← (ρ x ′ i , ρ w ′ )] | ρ ∈ X} • (X; y, ∅) = X µ w {G} (X; x i , y) = G ol (X; (x i , w), (y, w))whereΣ = {ρ | σ l = σ l ′ , l ∈ {x i , w, y}}.The concrete domain A is abstracted with the time-varying abstraction (as introduced in Section 5.6) where sets are intervals and the time domain isT d A * = (T d → Int).Since the interval domain is a complete lattice structure, this instance of the timevarying set abstraction has a complete lattice structure as well. The two domains A and A * are linked by the concretisation function γ(X * ) = (γ Int (X * (1)) × (γ (T d →Int) (X * )) 2 ) 2

Figure 5 .Figure 5 . 9 :

 559 Figure 5.9: Decreasing sequence {Y k } of overapproximations of the reachable tube Y cl of the signal y of the interconnection of systems S cl .

  x0 × Dom u × Dom y , where Dom x0 = R 2 , and Dom u = Dom y = L(T ; R). An environment ρ = (x 0 , u, y) belongs to the semantic of the system [S ] (i.e. [S ](ρ; ∅, ∅) is true) whenever the following relations are satisfied 

Figure 5 . 13 :Figure 5 . 14 :

 513514 Figure 5.13: Sound approximation of the relationship y(t) → u(t).

  Definition of the system . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Motivation of the QC model . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Abstraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 A simple example . . . . . . . . . . . . . . . . . . . . . . . . . Ellipsoidal method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Support ellipsoids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Optimal ellipsoids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Application of the Pontryagin's Minimum Principle . . . . . . 2.5.2 A continuation method to solve the PMP . . . . . . . . . . . . 2.6 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices 2.A Derivatives of the cost functions . . . . . . . . . . . . . . . . . . . . . 2.B Variations of the state and co-state . . . . . . . . . . . . . . . . . . . 2.C Contracting property in the centered case . . . . . . . . . . . . . . . . Definition of the system . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Motivation and examples of IQC systems . . . . . . . . . . . . . . . .

	Contents	vi
	2.3 3 Integral Quadratic Constraints
	3.1
	2 Quadratic Constraints
	2.1 v

1 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.2 Conservatism . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Overapproximation of the reachable set with time-varying conics . . .

Table 2 :

 2 Summary of the systems studied in this thesis. Abbreviations: CT (Continuous Time), LTV (Linear Time-Varying), NL (Non Linear), QC (Quadratic Constraint), IQC (Integral Quadratic Constraint), and NLIC (Non Linear Integral Constraint).

	Contents

  Definition of the system . . . . . . . . . . . . . . . . . . . . . . . 2.2 Motivation of the QC model . . . . . . . . . . . . . . . . . . . . 2.2.1 Abstraction . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 A simple example . . . . . . . . . . . . . . . . . . . . . . . 2.3 Ellipsoidal method . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Support ellipsoids . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Optimal ellipsoids . . . . . . . . . . . . . . . . . . . . . . . . . .

	Future works The DRE of interest satisfies many interesting properties one of
	them being that when t → E(t) solves the DRE (1.21), P ′ : t → k(t)P (t) solves the same DRE (1.21) provided some conditions over k and the parameters of the
	system hold, in other words, the DRE is separable. Such a result provides a lot of
	information on how the reachable tube of such systems behaves depending on the
	initial set of states. On another note, the DRE 1.21 is closely related to its associated Chapter 2
	Continuous Algebraic Riccati Equation (CARE; the CARE corresponds to constant
	solutions to DRE 1.21 with a free initial condition). The CARE appears in the
	Kalman-Yakubovitch-Popov Lemma when studying the robust stability of uncertain systems. In this case, the quantity λλ -1 is related to the input-to-output L 2 gain. A Quadratic Constraints
	clear understanding of this connection would help to choose an appropriate value for
	λ.
	Contents
	2.1 2.5.1 Application of the Pontryagin's Minimum Principle . . . .
	2.5.2 A continuation method to solve the PMP . . . . . . . . .
	2.6 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.A Derivatives of the cost functions . . . . . . . . . . . . . . . . . .
	2.B Variations of the state and co-state . . . . . . . . . . . . . . . .
	2.C Contracting property in the centered case . . . . . . . . . . . . .
	t 0 y ⊤ (s)M (s)y(s)ds ≥ 0 at each t ≥ 0.
	Set-based Cosimulation Paul Rousse	Paul Rousse Set-based Cosimulation

  Overapproximation with time-varying paraboloids . . . . . . . . 3.6 Exact reachable set . . . . . . . . . . . . . . . . . . . . . . . . . Exact reachable set . . . . . . . . . . . . . . . . . . . . . . 3.7 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.1 Examples from COMPleib . . . . . . . . . . . . . . . . . . 3.8.2 System verification . . . . . . . . . . . . . . . . . . . . . . 3.8.3 Delayed system . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.A Continuous extension of the domain of definition of the timevarying paraboloids . . . . . . . . . . . . . . . . . . . . . . . . .

Contents 3.1 Definition of the system . . . . . . . . . . . . . . . . . . . . . . . 3.2 Motivation and examples of IQC systems . . . . . . . . . . . . . 3.3 Overapproximation of the reachable set with time-varying conics 3.4 Extended system . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Paraboloids . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 3.6.1 Overapproximation with an intersection of time-varying paraboloids . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Overapproximation relationship . . . . . . . . . . . . . . . 3.6.3 Past trajectory for states in the interior of the overapproximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.4 Past trajectory for states on the boundary of the overapproximation . . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.5

  .13) corresponds to the set of all time-varying paraboloids with initial conditions P 0 , generated by the set of positive time-varying multipliers γ ∈ L 2,loc (R + ; R + ), and the set of initial multipliers λ 0 ≥ 1. Let

		3.6. Exact reachable set
	Π Π * (t) = P ∈ Π * t ∈ I(P)	(3.14)
	be the set of all the defined time-varying paraboloids at time t ≥ 0. For t ≥ 0, let Π(t) =
	Paul Rousse	Set-based Cosimulation

* P∈ Π * (t)

  might introduce elements of Dom that are not solutions to

	Chapter 5. Set-based Cosimulation	130
	the fixed-point equation (see Example 5.2). To avoid this, we reformulate (5.4)
	with
	λ = g(λ)	(5.7)
	Set-based Cosimulation	Paul Rousse

  In this work, we use a relaxed form of the narrowing operator where the ultimately stationary property is not required and where the property X

* △ * Y * ⊑ * Y * for every X * and Y * in A * holds for a different partial order ⊏ ∼ . Then, iterates Y * k decreases with respect to the partial ordering ⊏ ∼ . Definition 5.8. Partial narrowing operator An abstract binary operator △ * is a narrowing with respect to the partial order ⊏ ∼ * if and only if, for all X * , Y * ∈ A * , it holds

  PWL ) whenever (p, x) can be expressed as a convex combination of the points (p k , x k ), with k ∈ K j , where p k is a vertex of P j and where (p k , x k ) ∈ Z PWL (k) (seeExample 5.4).

	Example 5.4.
	Let the set Z be

1

. In other words, let (p, x) ∈ Dom such that p belongs to the cell P j , j ∈ I. The point (p, x) belongs to γ PWL (Z

Paul RousseSet-based Cosimulation

Set-based CosimulationPaul Rousse

where dist([x], [y]) = sup x∈[x],y∈[y] xy .Paul RousseSet-based Cosimulation

For two subsets A and B of a vectorial space, x ∈ λ • A + µ • B whenever x = λa + µb for some a ∈ A and b ∈ B, λ, µ ≥ 0 and λ + µ = 1.Set-based CosimulationPaul Rousse

Remerciements

Cardinal limitation of Π *

In order to have a tractable integration of the reachable set computation, we limit the cardinality of Π * in the following way

• at each time step kT c , we consider only the N new scaled paraboloids of the largest initial multiplier;

• Π * below N P , oldest time-varying paraboloids are dismissed in benefit of more recent ones;

N new and N P are user-defined parameters. Choosing the paraboloids with this heuristic showed good results in practice. These rules try to consider only elements of Π * that are more stable. Since for two solutions E and Ẽ of (1.21) respectively defined over [0, T ] and [0, T ] where T , T ∈ R ∪ {∞}, if E (0) Ẽ (0), then E (t) Ẽ (t) for t in the interval of definition of E and Ẽ , we have T ≥ T (this property follows directly by writing the corresponding value function of the basic LQR optimization problem). Therefore, for a time-varying paraboloid that is positive definite at t > 0, its scaled time-varying paraboloid at t will be defined for a longer time horizon. DRE numerical integration DRE integration is subject to numerical instability. Numerical integration of the DRE (1.21) does not produce good results in practice (see [START_REF] Kenney | Numerical Integration of the Differential Matrix Riccati Equation[END_REF]). Experiments presented in this works make use of the Chandrasekhar method [Lainiotis, 1976]. This method integrates the Ordinary Differential Equation (ODE) (1.21) E using an intermediate ODE over the timedependent matrix L in L 2 (R + , R n×n ) as follow

where Ė0 = Ė(0) given by (1.21). Then E is a solution to (1.21). Since L(t)L(t) ⊤ 0, this method is only applicable to strictly increasing solutions of the DRE. As seen in Example 3.2, the solutions to ODE (1.21) are not strictly increasing over the time horizon, even for a positive definite initial condition. Therefore, the Chandrasekhar method cannot be used directly. We instead use the following approach, let L, K ∈ L 2 (R + ; R n×n ) be such as:

Syntax and semantic of an interconnection of systems

We study an interconnection of systems that can be expressed with two constructions between subsystems: a serial connection and a feedback connection. Each connection is associated with a signal and an equation satisfied by this signal. The problem of interest is then to identify the set of signals of the interconnection that satisfies this set of equations.

In this part, we present the syntax, or language, used to describe our interconnection of systems. This syntax is introduced for two reasons. First, it identifies the class of systems we can analyze. Second, it gives a structure to the set of equations described by the interconnection of systems. This structure is, later on, reused to analyze the interconnection of systems. Once the syntax is defined, the actual mathematical meaning of the syntax is expressed by its semantic. We associate to each construct in our syntax an equation that should be satisfied by the signals of the interconnection.

An interconnection of systems is described with an ISys expressed within the syntax detailed in Table 5.1.

V a set of labels.

(Labels)

Table 5.1: Syntax of an interconnection of systems ISys.

Set-based Cosimulation Paul Rousse

Proof. Let ρ ∈ R = S (R), we want to show that ρ is a system trajectory, i.e.

[S ](ρ; ∅, ∅) is true. Since R is a fixed-point of X → S (X) and each operator only modifies a projection of R, and R is as well a fixed-point for each operator. Moreover, for every label u ∈ V, R ⊆ Σ, i.e. ρ u = ρ u ′ . We prove that ρ is a trajectory by induction. For each operator, we prove that R = S (R; u, y) implies that [S](ρ; u, y) is true:

• the last operators have a similar proof than the one for the serial connection. ♦ Proposition 5.2 and Proposition 5.

. Therefore, the semantic of the interconnection of systems is equal to the greatest fixed-point of X → S (X). Definition 5.4. Concrete semantic as the greatest fixed-point of a monotonic operator Let Σ ⊆ Dom defined as in (5.8). The concrete semantic S is defined by

where the monotonic operator X → S (X; ∅, ∅) is defined in Table 5.3 by structural induction over the syntax of the interconnection of systems S .

Proposition 5.4. Concrete semantic

The semantic of the interconnection of systems S is equal to the concrete semantic, i.e.

[S ] = S .

Remark 5.4. Concrete semantic as the greatest fixed-point of a monotonic operator

The approach to express the system trajectories (or the computer program traces) as a fixed-point of a monotonic operator in a complete lattice structure is classical. However, the fixed-point of interest is usually the lowest fixed-point of a monotonic operator.

Paul Rousse

Set-based Cosimulation 5.4. Abstract semantic

Abstract semantic

The previous section showed that the concrete domain and abstract domain can be both provided with a mathematical structure of complete lattices. In what follows, we assume to have a Galois connection between the concrete domain (A, ⊑, ⊔, ⊓) and an abstract domain

Now that abstract domains are defined, we should define the abstract semantic. We will see that the abstract semantic mimics the concrete semantic: each function evaluation in the concrete domain gives rise to a function evaluation in the abstract domain, each fixed-point computation in the concrete domain gives rise to a fixedpoint computation in the abstract domain. For a given concrete semantic S , we can compute an abstract semantic S * . This abstract semantic is built by structural induction over the syntax of the system interconnection. This semantic is sound when α( S ) ⊑ * S * , or equivalently when S ⊑ γ( S * ). The soundness property is obtained by using sound operator evaluation and fixed-point computation over abstract elements.

In Section 5.4.1, we detail the abstract counterpart of the concrete functions. In Section 5.4.2, we detail the abstract fixed-point computation.

Abstract evaluation

For a function F : A → A in the concrete domain, it is sometimes possible to define F * : A * → A * an "abstract evaluation" of F . To ensure the soundness of the abstract semantic, the soundness of F * with respect to the Galois connection and F should be ensured.

Definition 5.7. Sound approximation of a function

For a function

Abstract fixed-points

The previous section details how the abstract counterpart of concrete functions can be defined. We now detail the computation of the greatest fixed-point that appears in the concrete semantic. Since the greatest fixed-point in the concrete domain exists but is not necessarily computable, we aim at computing a sound approximation of the greatest fixed-point in the abstract domain by using a fixed-point transfer theorem:

Set-based Cosimulation Paul Rousse

for F * : A * → A * a monotonic function and F : A → A a monotonic function.

When F * is a sound approximation of F then

Proof. See Theorem 2 in [Cousot, 2002] ♦

When each operator is sound, the greatest fixed-point in the concrete domain is soundly approximated by the abstract greatest fixed-point in the abstract domain.

Practical computation of fixed-point

Theorem 5.3 gives a way to soundly approximate the fixed-point of the concrete domain by computing a greatest fixed-point in the abstract domain. Contrary to the concrete domain, the abstract domain is computer-representable and abstract evaluations are computable. This greatest fixed-point can be computed using an iterative approach as in Property 5.1: let assume that for a D * ∈ A * , γ(D * ) overapproximates gfp X {F }, i.e. gfp X {F } ⊑ γ(D * ), that D * is a post-fixpoint of F * , i.e.

and that F * is a monotonic operator, then the sequence {Y * k } of iterates, defined by

is monotonically decreasing, i.e. Y * k+1 ⊑ * Y * k for every k ≥ 0, and each iterate is a sound approximation of gfp X {F }, i.e.

they correspond to a refinement of γ(D * ).

Paul Rousse

Set-based Cosimulation

Piecewise linear abstraction

Proposition 5.10. Galois connection of the domain of time-varying sets If the pair (α, γ) defines a Galois connection between the two complete lattices (℘(R n ), ⊆, ∪, ∩, ⊥, ⊤) and (A * , ⊑ * , ⊔ * , ⊓ * , ⊥ * , ⊤ * ), then (α TV , γ TV ), defined by

defines a Galois connection between the domain of time-varying signals

Proof. ♦

Time-varying signals can be approximated with time-varying ellipsoids which is a point-wise lift of the ellipsoidal domain introduced in Section 5.5.1.

Example 5.3.

The set of signals

can be soundly approximated with the time-varying ellipsoid X * = α TV (X) (see Figure 5.5) defined by

where E is the time-varying ellipsoid defined by its time-varying center c(t) = 0 and its time-varying radius r(t) = e t 10 , for ever t ≥ 0. The set γ TV (X * ) contains all the signals x : R + → R belonging to (x(t)c(t)) 2 ≤ r(t) 2 at every t ∈ R + .

Piecewise linear abstraction

Let us assume that the interconnection of systems involves a source signal " P p " where p is of finite dimension (e.g. p is a real-valued vector) and where P is bounded (e.g. p could be an unknown initial state or an uncertain parameter of the system). The concrete domain Dom corresponds to the Cartesian product of Dom p and of the domain Dom x of other signals x labeled by x in the interconnection. Each value p ∈ P is associated with a subset of the semantic, i.e. a subset of Dom x . In this section, we propose to define a piecewise linear abstraction that expresses this relationship between p (over a partition of P) and x in the concrete semantic. Section 5.8.2 uses such abstraction to derive tight approximation.

We first define the partition of P, then we introduce the piecewise linear abstraction, and finally, we detail the Galois connection of this abstraction with the concrete domain.

Set-based Cosimulation

Paul Rousse Partition Let a partition of P be described by a finite set of cells P i ⊆ P, i ∈ I, with I a set of labels associated with each cell, such that:

• cells are convex polytopes, to this respect, each cell P i can be described as the convex hull of a finite set of points; • cells are disjoint, i.e. i = j ⇒ P o i ∩ P o j = ∅ (where P o is the open set associated with the set P ); • cells cover the entire set P, i.e. P = i∈I P i . Let {p k } k∈K be the vertices of the partition with K a set of labels identifying the vertices. For a cell P i , let K i ⊆ K the labels of the vertices of P i . For each vertex p k , let the cell neighborhood of p k , denoted by Pk , be the union of cells touching p k (see Figure 5.6).

Piecewise linear abstraction Let the three complete lattices (A p , ⊆, ∪, ∩, ⊥, ⊤), (A x , ⊆, ∪, ∩, ⊥, ⊤), and (A, ⊆, ∪, ∩, ⊥, ⊤) where A p = ℘(Dom p ), A x = ℘(Dom x ), and A = ℘(Dom p × Dom x ) such that Dom p and Dom x are vector spaces. Let the piecewise linear (PWL) domain be defined as the point-wise lift (introduced in Definition 5.9) of A x over the set of vertices K, i.e.

(A PWL , ⊆, ∪, ∩, ⊥, ⊤).

where A PWL is a subset of K → A, the set of functions that associates to each K an element of A. Elements Z PWL ∈ A PWL associates to a k ∈ K an element {p k } × X k

Paul Rousse

Set-based Cosimulation

Piecewise linear models

In this section, we study a nonlinear system described as the interconnection of a linear time-invariant system with a nonlinear feedback. The response of the linear system is overapproximated with results of Part I, and the one of the nonlinear feedback is overapproximated with the interval arithmetic framework (as presented in Section 4.1 of Chapter 4). The concrete domain of this system is abstracted with a PWL abstraction as presented in Section 5.7.

Description of the interconnection of systems

We study the interconnection of systems S described by the expression (5.19) (see