Parabolic Set Simulation for Reachability Analysis of Linear Time   Invariant Systems with Integral Quadratic Constraint

Paul Rousse; Pierre-Lo\"ic Garoche; Didier Henrion (LAAS-MAC)

arXiv:1902.10982·math.OC·February 14, 2020

Parabolic Set Simulation for Reachability Analysis of Linear Time Invariant Systems with Integral Quadratic Constraint

Paul Rousse, Pierre-Lo\"ic Garoche, Didier Henrion (LAAS-MAC)

PDF

TL;DR

This paper introduces a novel parabolic set simulation method for reachability analysis of LTI systems with integral quadratic constraints, providing a practical way to overapproximate reachable sets in energy-constrained systems.

Contribution

It extends ellipsoidal techniques to parabolic sets for LTI systems with IQCs, including a method to generate supporting paraboloids for exact reachable set characterization.

Findings

01

Overapproximates reachable sets with time-varying parabolic sets.

02

Proves the relationship between paraboloids and reachable sets.

03

Applicable to systems with delays, rate limiters, and energy constraints.

Abstract

This work extends reachability analyses based on ellipsoidal techniques to Linear Time Invariant (LTI) systems subject to an integral quadratic constraint (IQC) between the past state and disturbance signals , interpreted as an input-output energetic constraint. To compute the reachable set, the LTI system is augmented with a state corresponding to the amount of energy still available before the constraint is violated. For a given parabolic set of initial states, the reachable set of the augmented system is overapproximated with a time-varying parabolic set. Parameters of this paraboloid are expressed as the solution of an Initial Value Problem (IVP) and the overapproximation relationship with the reachable set is proved. This paraboloid is actually supported by the reachable set on so-called touching trajectories. Finally , we describe a method to generate all the supporting…

Tables2

Table 1. Table 1 : Computation times (in seconds) of the overapproximation for different systems sizes, using a unique time-varying paraboloid. (When the original system’s size is smaller than the required reduced system size, then the model reduction is not applicable - n.a. -.)

System size	Helicopter (HE7)	Aircraft (AC10)	Coupled Spring (CSE1)
5	4.32	4.64	3.65
10	5.12	5.96	3.86
19	7.42	10.62	7.92
30	n.a.	28.85	n.a.
40	n.a.	50.66	n.a.
49	n.a.	88.00	n.a.

Table 2. Table 2 : Computation times (in seconds) and number of paraboloids (in parenthesis) of the overapproximation for different systems sizes. (When the original system’s size is smaller than the required reduced system size, then the model reduction is not applicable - n.a. -.)

System size	Helicopter (HE7)	Aircraft (AC10)	Coupled Spring (CSE1)
5	83.63 (66)	36.64 (13)	213.88 (232)
10	89.55 (57)	25.77 (9)	261.32 (197)
19	167.53 (52)	27.67 (4)	21.97 (4)
30	n.a.	113.72 (7)	n.a.
40	n.a.	117.60 (4)	n.a.

Equations184

{\dot{\mathbold x} (τ) \mathbold x (0) = A \mathbold x (τ) + B \mathbold w (τ) + B_{u} \mathbold u (τ) = x_{0} with τ \in [0, t]

{\dot{\mathbold x} (τ) \mathbold x (0) = A \mathbold x (τ) + B \mathbold w (τ) + B_{u} \mathbold u (τ) = x_{0} with τ \in [0, t]

x_{q 0} + 0 \int τ [\mathbold x (s) \mathbold u (s) \mathbold w (s)]^{⊤} M [\mathbold x (s) \mathbold u (s) \mathbold w (s)] d s \geq 0 for all τ \in [0, t]

x_{q 0} + 0 \int τ [\mathbold x (s) \mathbold u (s) \mathbold w (s)]^{⊤} M [\mathbold x (s) \mathbold u (s) \mathbold w (s)] d s \geq 0 for all τ \in [0, t]

M = M_{x} M_{xu}^{⊤} M_{x w}^{⊤} M_{xu} M_{u} M_{u w}^{⊤} M_{x w} M_{u w} M_{w} \in S^{n + m + p}

M = M_{x} M_{xu}^{⊤} M_{x w}^{⊤} M_{xu} M_{u} M_{u w}^{⊤} M_{x w} M_{u w} M_{w} \in S^{n + m + p}

\mathbold x_{q} (τ) = x_{q 0} + 0 \int τ [\mathbold x (s) \mathbold u (s) \mathbold w (s)]^{⊤} M [\mathbold x (s) \mathbold u (s) \mathbold w (s)] d s,

\mathbold x_{q} (τ) = x_{q 0} + 0 \int τ [\mathbold x (s) \mathbold u (s) \mathbold w (s)]^{⊤} M [\mathbold x (s) \mathbold u (s) \mathbold w (s)] d s,

\mathbold x_{q} (τ) \geq 0 for all τ \in [0, t] .

\mathbold x_{q} (τ) \geq 0 for all τ \in [0, t] .

\mathbold{z}=(\mathbold{x},\mathbold{x_{q}})\in\mathscr{S}(\mathcal{Z}_{0},t)\Leftrightarrow\left\{\begin{array}[]{l}\textrm{$\mathbold{x}$ solves \eqref{eq:ode}}\\ \textrm{and $\mathbold{x_{q}}$ solves \eqref{eq:xq_int}}\\ \textrm{with $(x_{0},x_{q0})\in\mathcal{Z}_{0}$}\\ \textrm{$\mathbold{x_{q}}$ satisfies \eqref{eq:state_cons}}\end{array}\right.

\mathbold{z}=(\mathbold{x},\mathbold{x_{q}})\in\mathscr{S}(\mathcal{Z}_{0},t)\Leftrightarrow\left\{\begin{array}[]{l}\textrm{$\mathbold{x}$ solves \eqref{eq:ode}}\\ \textrm{and $\mathbold{x_{q}}$ solves \eqref{eq:xq_int}}\\ \textrm{with $(x_{0},x_{q0})\in\mathcal{Z}_{0}$}\\ \textrm{$\mathbold{x_{q}}$ satisfies \eqref{eq:state_cons}}\end{array}\right.

R (Z_{0}, t) = {\mathbold z (t) ∣ \mathbold z \in S (Z_{0}, t)} .

R (Z_{0}, t) = {\mathbold z (t) ∣ \mathbold z \in S (Z_{0}, t)} .

\begin{array}[]{rcl}h:&\mathbb{R}^{n}\times\mathbb{R}&\to\mathbb{R}\\ &(x,x_{q})&\mapsto x^{\top}Ex-2f^{\top}x+g+x_{q},\end{array}

\begin{array}[]{rcl}h:&\mathbb{R}^{n}\times\mathbb{R}&\to\mathbb{R}\\ &(x,x_{q})&\mapsto x^{\top}Ex-2f^{\top}x+g+x_{q},\end{array}

P (E, f, g) = {(x, x_{q}) \in R^{n + 1} h (x, x_{q}) \leq 0} .

P (E, f, g) = {(x, x_{q}) \in R^{n + 1} h (x, x_{q}) \leq 0} .

\dot{\mathbold E} (t) = - \mathbold E (t) A - A^{⊤} \mathbold E (t) - M_{x} + (B^{⊤} \mathbold E (t) + M_{x w}^{⊤})^{⊤} M_{w}^{\scalebox 0.6 - 1} (B^{⊤} \mathbold E (t) + M_{x w}^{⊤}) .

\dot{\mathbold E} (t) = - \mathbold E (t) A - A^{⊤} \mathbold E (t) - M_{x} + (B^{⊤} \mathbold E (t) + M_{x w}^{⊤})^{⊤} M_{w}^{\scalebox 0.6 - 1} (B^{⊤} \mathbold E (t) + M_{x w}^{⊤}) .

\dot{\mathbold f} (t) = - A^{⊤} \mathbold f (t) + (M_{xu} + \mathbold E (t) B_{u}) \mathbold u (t) + (\mathbold E (t) B + M_{x w}) M_{w}^{\scalebox 0.6 - 1} (B^{⊤} \mathbold f (t) - M_{u w}^{⊤} \mathbold u (t)) .

\dot{\mathbold f} (t) = - A^{⊤} \mathbold f (t) + (M_{xu} + \mathbold E (t) B_{u}) \mathbold u (t) + (\mathbold E (t) B + M_{x w}) M_{w}^{\scalebox 0.6 - 1} (B^{⊤} \mathbold f (t) - M_{u w}^{⊤} \mathbold u (t)) .

\mathbold g (t) = g_{0} + 0 \int t [\mathbold f (τ) \mathbold u (τ)]^{⊤} G [\mathbold f (τ) \mathbold u (τ)] d τ

\mathbold g (t) = g_{0} + 0 \int t [\mathbold f (τ) \mathbold u (τ)]^{⊤} G [\mathbold f (τ) \mathbold u (τ)] d τ

G = [B M_{w}^{\scalebox 0.6 - 1} B^{⊤} (B_{u} - B M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤})^{⊤} B_{u} - B M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤} - M_{u} + M_{u w} M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤}] .

G = [B M_{w}^{\scalebox 0.6 - 1} B^{⊤} (B_{u} - B M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤})^{⊤} B_{u} - B M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤} - M_{u} + M_{u w} M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤}] .

\begin{array}[]{rl}P\colon I&\to\mathbb{P}\\ t&\mapsto\mathcal{P}(\mathbold{E}(t),\mathbold{f}(t),\mathbold{g}(t))\end{array}

\begin{array}[]{rl}P\colon I&\to\mathbb{P}\\ t&\mapsto\mathcal{P}(\mathbold{E}(t),\mathbold{f}(t),\mathbold{g}(t))\end{array}

h_{z} (t) = [\mathbold x (t) 1]^{⊤} [\mathbold E (t) - \mathbold f^{⊤} (t) - \mathbold f (t) \mathbold g (t)] [\mathbold x (t) 1]^{⊤} + \mathbold x_{q} (t)

h_{z} (t) = [\mathbold x (t) 1]^{⊤} [\mathbold E (t) - \mathbold f^{⊤} (t) - \mathbold f (t) \mathbold g (t)] [\mathbold x (t) 1]^{⊤} + \mathbold x_{q} (t)

\dot{h}_{z} (t)

\dot{h}_{z} (t)

+ 2 [\mathbold x (t) 1]^{⊤} [\mathbold E (t) - \mathbold f^{⊤} (t) - \mathbold f (t) \mathbold g (t)] [\dot{\mathbold x} (t) 0]

+ [\mathbold x (t) \mathbold u (t) w_{t}]^{⊤} M [\mathbold x (t) \mathbold u (t) w_{t}] .

k (w_{s}) =

k (w_{s}) =

+ w_{t}^{⊤} M_{w} w_{t}

\mathbold w^{*} = - M_{w}^{\scalebox 0.6 - 1} (B^{⊤} (\mathbold E \mathbold x - \mathbold f) + [M_{x w} M_{u w}]^{⊤} [\mathbold x \mathbold u]) .

\mathbold w^{*} = - M_{w}^{\scalebox 0.6 - 1} (B^{⊤} (\mathbold E \mathbold x - \mathbold f) + [M_{x w} M_{u w}]^{⊤} [\mathbold x \mathbold u]) .

\forall t \in I (P), R (P_{0}, t) \subseteq P (t) \cap Z_{+} .

\forall t \in I (P), R (P_{0}, t) \subseteq P (t) \cap Z_{+} .

\mathbold z (0) \in P (0) \Rightarrow \mathbold z (t) \in P (t) for all t \in I (P) .

\mathbold z (0) \in P (0) \Rightarrow \mathbold z (t) \in P (t) for all t \in I (P) .

\begin{array}[]{rl}\displaystyle\max_{\mathbold{w}\in\mathbf{L}_{2}([0,t];\mathbb{R}^{m})}{}&\displaystyle\int\limits_{0}^{t}\left[\begin{smallmatrix}\mathbold{x}(\tau)\\ \mathbold{u}(\tau)\\ \mathbold{w}(\tau)\end{smallmatrix}\right]M\left[\begin{smallmatrix}\mathbold{x}(\tau)\\ \mathbold{u}(\tau)\\ \mathbold{w}(\tau)\end{smallmatrix}\right]d\tau-x_{q,t}\\ \textrm{s.t.}&\dot{\mathbold{x}}=A\mathbold{x}+B\mathbold{w}+B_{u}\mathbold{u}\\ &\mathbold{x}(t)=x_{t}\\ \end{array}

\begin{array}[]{rl}\displaystyle\max_{\mathbold{w}\in\mathbf{L}_{2}([0,t];\mathbb{R}^{m})}{}&\displaystyle\int\limits_{0}^{t}\left[\begin{smallmatrix}\mathbold{x}(\tau)\\ \mathbold{u}(\tau)\\ \mathbold{w}(\tau)\end{smallmatrix}\right]M\left[\begin{smallmatrix}\mathbold{x}(\tau)\\ \mathbold{u}(\tau)\\ \mathbold{w}(\tau)\end{smallmatrix}\right]d\tau-x_{q,t}\\ \textrm{s.t.}&\dot{\mathbold{x}}=A\mathbold{x}+B\mathbold{w}+B_{u}\mathbold{u}\\ &\mathbold{x}(t)=x_{t}\\ \end{array}

\mathbold n = \mathbold E \mathbold x - \mathbold f .

\mathbold n = \mathbold E \mathbold x - \mathbold f .

[\dot{\mathbold x} \dot{\mathbold n}] = L [\mathbold x \mathbold n] + N \mathbold u

[\dot{\mathbold x} \dot{\mathbold n}] = L [\mathbold x \mathbold n] + N \mathbold u

L = [A - B M_{w}^{\scalebox 0.6 - 1} M_{x w}^{⊤} - (M_{x} - M_{x w} M_{w}^{\scalebox 0.6 - 1} M_{x w}^{⊤}) - B M_{w}^{\scalebox 0.6 - 1} B^{⊤} - A^{⊤} + M_{x w} M_{w}^{\scalebox 0.6 - 1} B^{⊤}]

L = [A - B M_{w}^{\scalebox 0.6 - 1} M_{x w}^{⊤} - (M_{x} - M_{x w} M_{w}^{\scalebox 0.6 - 1} M_{x w}^{⊤}) - B M_{w}^{\scalebox 0.6 - 1} B^{⊤} - A^{⊤} + M_{x w} M_{w}^{\scalebox 0.6 - 1} B^{⊤}]

N = [B_{u} - B M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤} - (M_{xu} - M_{x w} M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤})] .

N = [B_{u} - B M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤} - (M_{xu} - M_{x w} M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤})] .

\dot{\mathbold r} = \mathbold u (H R) [\mathbold x \mathbold n \mathbold u]

\dot{\mathbold r} = \mathbold u (H R) [\mathbold x \mathbold n \mathbold u]

H = (M_{u w} M_{w}^{\scalebox 0.6 - 1} M_{x w} - M_{xu}^{⊤} - (B_{u} - M_{u w} M_{w}^{\scalebox 0.6 - 1} B^{⊤}))

H = (M_{u w} M_{w}^{\scalebox 0.6 - 1} M_{x w} - M_{xu}^{⊤} - (B_{u} - M_{u w} M_{w}^{\scalebox 0.6 - 1} B^{⊤}))

R = M_{u} - M_{u w} M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤} .

R = M_{u} - M_{u w} M_{w}^{\scalebox 0.6 - 1} M_{u w}^{⊤} .

h_{t} (\mathbold x (t)) = \mathbold x (t)^{⊤} \mathbold n (t)^{⊤} + \mathbold r (t) .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Parabolic Set Simulation for Reachability Analysis of Linear Time-Invariant Systems with Integral Quadratic Constraint

Paul Rousse

[email protected]

Pierre-Loïc Garoche

[email protected]

Didier Henrion

[email protected]

ONERA, 31400 Toulouse, France

LAAS-CNRS, Université Toulouse, CNRS, Toulouse, France and Faculty of Electrical Engineering, Czech Technical University in Prague, Czechia,

Abstract

This paper describes the computation of reachable sets and tubes for linear time-invariant 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. Compared to previous methods based on the classical linear quadratic regulator, our approach can be applied to unstable systems as well. We tested our approach on large scale systems.

keywords:

Reachability analysis, Set-based simulation, Integral Quadratic Constraints, Uncertain systems

1 Introduction

We consider the reachability problem for Linear Time-Invariant (LTI) with Integral Quadratic Constraints (IQC). Reachable set computation is an active field of research in control theory (see [1]). It has many applications such as state estimation (see [2]) or verification (see [3]) of dynamical systems. IQC is a classical tool of robust control theory (see e.g. [4, 5]). It can model infinite dimensional states, non-linear dynamics, delays, rate limiters, uncertain systems (see [6, 7, 8] and [9]). Up to now, IQCs have mainly been used to evaluate the stability of systems. Despite their ability to model complex systems, we are still lacking results: we do not have a proper characterization of the reachable set.

In this paper, we extend reachability analysis based on ellipsoidal techniques (see e.g. [10, 11, 12]) for LTI systems subject to an IQC. This IQC is a trajectory constraint (i.e. valid at any time) between past state-trajectory, input signals, and unknown disturbance signals. To override dealing with constraints over the state-trajectories, we study the LTI system augmented with a state corresponding to the integral term in the IQC. For a given parabolic set of initial states, the reachable set of the augmented system is overapproximated by a time-varying parabolic set. Parameters of this paraboloid are expressed as the solution to an Initial Value Problem (IVP) partially described by a Differential Riccati Equation (DRE). The paraboloid is a tight overapproximation as it stays in contact with the boundary of the reachable set on so-called touching trajectories. By studying touching trajectories that are close to violating the constraint, we find conditions to generate all the supporting time-varying parabolic sets. At a given time, the intersection of these supporting parabolic sets is an exact representation of the reachable set. A wise choice of these paraboloids can be made to always have a bounded overapproximation of the reachable set. An algorithm to overapproximate the reachable set is provided (with an adaptation from [13] of the Chandrasekhar method for numerical integration of the DRE in the case where the solution is not sign-definite). An open source version of our implementation is available on-line [14].

Related work

Reachability analysis of LTI systems with ellipsoidal bounded inputs is studied in [10, 11, 12]. Such systems can model infinity norm bounded input-output LTI systems. The reachable set (which is convex and bounded; see [11]) can be overapproximated with time-varying ellipsoidal sets. Each ellipsoid is described by its parameters (center and radius) that are solutions to an 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. Our work presents a similar approach for a different class of systems, namely LTI with IQC constraint; and different sets, instead of ellipsoids, we use paraboloids.

An optimal control formulation of the reachable set problem is also possible [15, 16]. For a given state, if the maximal integral cost leading to this state violates the constraint, then this state is unreachable. It can be solved (using Hamilton-Jacobi-Bellman -HJB- viscosity subsolutions, see [17]) leading to global constraints over the reachable set. If the reachable set can be expressed as the intersection (possibly uncountable) of elements of the chosen function family, then the intersection of the resulting constraints gives an exact representation of the reachable set. However, HJB solutions are known to not scale well with the system dimension. Our specialized solution showed good results for large systems.

HJB based methods propagate the constraints along the flow of the dynamical system. Occupation measures and barrier certificates methods aim at finding constraints over the reachable tube of a dynamical system: [18] uses IQCs for verification purposes using barrier certificates where the positivity of the energetic state is ensured by using a nonnegative constant multiplier: [19, 20] use an occupation measure approach where the IQC can potentially be incorporated as a constraint over the moment of the trajectories (note however these references do not deal explicitly with IQCs). A hierarchy of semi-definite programs are derived for polynomial dynamics. Then, off-the-shelf Semi-Definite Programming (SDP) solvers are used to solve 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. However, few properties about the reachable set can be formalized such as its domain of existence, its boundedness. We hereby study a narrow class of systems that both have practical and theoretical importance (LTI+IQC system). Our methods scale favorably with the number of states, and the existence and boundedness of the overapproximation is always granted.

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 DRE. References [21, 22, 16] describe the reachable set of LTI systems with terminal IQC. Reference [23] formalizes the problem with a game theory approach. Reference [24] solves the differential Riccati inequality over a finite horizon using a basis of polynomial functions, then an SDP solver search for a solution that minimizes the final volume of the overapproximation. This algorithm has been implemented in available tools (see LTVTools toolbox, [25]). In all these works, the overapproximation of the reachable set is conditioned by the existence of a solution to the DRA 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 finite interval of time. By taking into account the constraint over the full integration range, we can always find an overapproximation valid over any interval of integration. Therefore, our method is less conservative than the standard LQR approach.

Contributions

We study the reachable set computation of an LTI system with IQC. To the knowledge of the authors, this is the first paper to provide a set-based solution for reachable set computation for LTI systems with IQC that make use of the constraint over the full integration time (not only at the terminal time). In the conference paper [26], we presented our preliminary results. Under strong assumptions on the system, we proved that our overapproximation is the exact reachable set computation. In the current paper:

•

We have an exact characterization of the reachable set for a larger class of systems compared to our previous results presented in [26].

•

We extend the existing ellipsoidal method [presented in 10, 11, 12] for the reachability analysis of bounded-input LTI systems to the reachability analysis of LTI systems with IQC. These parabolic constraints are defined by time-varying parameters which are the solution to an IVP. Part of this IVP (the quadratic coefficient of the parabolic constraint) is a DRE. The IVP convergence property is obtained thanks to the convergence property of the DRE.

•

The constraint in the IQC system is modeled as a constraint over the state of an augmented system. Unlike in [27, 21, 28, 23], where the constraint is only used at the terminal time, we rather use the constraint on the entire time-domain. The differential equations of the paraboloid’s parameters are then differently defined than in previous works. The differential equations depend on a so-called scaling function: a time-varying parameter defined by the user. This scaling function can always be chosen such that the overapproximation is defined over any time domain.

•

Our method can efficiently be used for large scale systems for which moment based methods [19, 20] could not be applied.

Outline

The LTI system with IQC and the reachability analysis problem are introduced (Section 1). Parabolic constraints and their associated parameter IVP are defined, their domain of definition is analyzed, the overapproximation property is formulated, as well as the touching trajectories (Section 2). A method to generate a set of time-varying parabolic constraints is described. The intersection of these paraboloids exactly describes the reachable set of the system (Section 3). An algorithm to overapproximate the reachable set is described (Section 4). This algorithm is evaluated for numerical examples and large-scale LTI systems (Section 5).

1.1 Notation

Let $\mathbb{S}^{n}\subset\mathbb{R}^{n\times n}$ the set of real valued symmetric square matrices of size $n$ . For $A\in\mathbb{S}^{n}$ , we write $A\succ 0$ (resp. $A\prec 0$ ) when $A$ is positive definite (resp. negative definite). We define the matrix norm $\left\lVert A\right\rVert=\sqrt{\mathrm{tr}(A^{\top}A)}$ for $A\in\mathbb{R}^{n\times m}$ , where $\mathrm{tr}(B)$ is the trace of $B\in\mathbb{R}^{n\times n}$ . Let a signal be a function that associates to a time instant in $\left[0,+\infty\right[$ a vector from $\mathbb{R}^{n}$ . For a given interval $I\subseteq\mathbb{R}$ , let $\mathbf{L}_{2}(I;\mathbb{R}^{n})$ denote the Hilbert space of signals equipped with the norm:

$\left\lVert\mathbold{u}\right\rVert=\sqrt{\int_{t\in I}\mathbold{u}^{T}(t)\mathbold{u}(t)dt}<\infty.$ Let $\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ be the set of locally square integrable signals from $\mathbb{R}^{+}$ to $\mathbb{R}^{+}$ . Let $|X|$ the cardinality of a countable set $X$ . For a set $\Omega\subset\mathbb{R}^{n}$ , let $\partial\Omega$ denote its boundary. Let $\mathscr{C}^{1}(I;\mathbb{R}^{n})$ the set of functions from $I$ to $\mathbb{R}^{n}$ which are continuous and differentiable with continuous derivative.

Let $\mathrm{sign}:\mathbb{R}\mapsto\{-1,0,1\}$ such that $\mathrm{sign}(s)=1$ is $s>0$ , $\mathrm{sign}(s)=-1$ if $s<0$ , $\mathrm{sign}(s)=0$ otherwise.

1.2 System

For a given input signal $\mathbold{u}\in\mathscr{C}^{1}(\mathbb{R}^{+};\mathbb{R}^{p})\cap\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{p})$ , given matrices $A\in\mathbb{R}^{n\times n}$ , $B\in\mathbb{R}^{n\times m}$ , $B_{u}\in\mathbb{R}^{n\times p}$ , and a given terminal time $t>0$ , we study the trajectories $\mathbold{x}\in\mathbf{L}_{2}([0,t];\mathbb{R}^{n})$ of the LTI system:

[TABLE]

where $\mathbold{w}\in\mathbf{L}_{2,loc}([0,t];\mathbb{R}^{m})$ is an unknown disturbance that satisfies:

[TABLE]

for given initial conditions $(x_{0},x_{q0})\in\mathbb{R}^{n}\times\mathbb{R}^{+}$ , and given symmetric matrix

[TABLE]

with $M_{w}\prec 0$ . Many IQC models verify such condition, as the one described in [4].

In this work, the constraint (2) is expressed as a constraint over a state $\mathbold{x_{q}}\in\mathbf{L}_{2}([0,t];\mathbb{R})$ defined for $s\in[0,t]$ by:

[TABLE]

then

[TABLE]

The constrained dynamical system $\mathscr{S}(\mathcal{Z}_{0},t)$ is then defined for a given set of initial states $\mathcal{Z}_{0}\subset\mathbb{R}^{n}\times\mathbb{R}$ and a terminal time $t>0$ :

[TABLE]

Define the reachable set:

[TABLE]

Then, $\mathcal{R}(\mathcal{Z}_{0},t)\subseteq\mathcal{Z}_{+}$ where $\mathcal{Z}_{+}=\mathbb{R}^{n}\times\mathbb{R}^{+}$ , and let $\mathcal{Z}_{*}=\mathbb{R}^{n}\times\{0\}$ .

1.3 Paraboloids

We overapproximate the reachable set $\mathcal{R}(\mathcal{Z}_{0},t)$ of $\mathscr{S}(\mathcal{Z}_{0},t)$ with paraboloids:

Definition 1 (Paraboloid).

Given $(E,f,g)\in\mathbb{S}^{n}\times\mathbb{R}^{n}\times\mathbb{R}$ , define the value function:

[TABLE]

and the paraboloid:

[TABLE]

Let $\mathbb{P}=\left\{\mathcal{P}(E,f,g)\middle|E\in\mathbb{S}^{n},f\in\mathbb{R}^{n},g\in\mathbb{R}\right\}$ be the set of paraboloids. Elements of $\mathbb{P}$ are not strictly speaking paraboloids since $E\succeq 0$ is not assumed in the definition.

Definition 2 (Scaled Paraboloid).

For $\mathcal{P}\in\mathbb{P}$ with parameters $(E,f,g)$ and a scaling factor $\gamma>0$ , let $\gamma\mathcal{P}\in\mathbb{P}$ be the scaled paraboloid defined by parameters $(\gamma E,\gamma f,\gamma g)$ .

Scaled paraboloids satisfy the following:

Property 1.

Given $\mathcal{P}\in\mathbb{P}$ and $\gamma\geq 1$ , it holds $\mathcal{P}\cap\mathcal{Z}_{+}\subseteq\gamma\mathcal{P}\cap\mathcal{Z}_{+}$ .

Proof 1.

Let $h$ and $h^{\prime}$ (resp.) the value functions of $(E,f,g)=\mathcal{P}$ and $\gamma\mathcal{P}$ (resp.) evaluated at $(x,x_{q})\in\mathcal{P}$ . Since $(x,x_{q})\in\mathcal{P}$ , $h\leq 0$ , i.e. $x^{\top}Ex-2f^{\top}x+g\leq-x_{q}$ . Then, $h^{\prime}=\gamma(x^{\top}Ex-2f^{\top}x+g)+x_{q}\leq-(\gamma-1)x_{q}$ . Since $(x,x_{q})\in\mathcal{Z}_{+}$ and since $\gamma-1\geq 0$ , we have $(\gamma-1)x_{q}\geq 0$ i.e. $h^{\prime}\leq 0$ meaning that $(x,x_{q})\in\gamma\mathcal{P}\cap\mathcal{Z}_{+}$ .

For $P$ a function that associates to a time $t$ of a time-interval $I\subset\mathbb{R}^{+}$ a set of states $P(t)\subset\mathbb{R}^{n+1}$ . We define a touching trajectory:

Definition 3 (Touching Trajectory).

A trajectory $\mathbold{z}^{*}$ solution to (1, 4) is a touching trajectory of $P$ when $\mathbold{z}^{*}(t)$ belongs to the surface of $P(t)$ at any time $t\in I$ , i.e. $\mathbold{z}^{*}(t)\in\partial P(t)$ .

1.4 Problem Statement

We are now ready to state the two problems studied in this work.

Problem 1.

Find an overapproximation of the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ at any $t>0$ for a given paraboloid of initial conditions $\mathcal{P}_{0}\in\mathbb{P}$ .

Theorem 1 in Section 2 solves Problem 1. It restates classical results about LQR systems applied to reachability analysis of IQC systems when the constraint (2) is a terminal time constraint. Theorem 2 in Section 3 is another solution to Problem 1. It extends the result in Theorem 1 taking into account the constraint over the entire interval of integration.

Problem 2.

Find a sequence of overapproximations that converges to the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ at any $t>0$ for a given paraboloid of initial conditions $\mathcal{P}_{0}\in\mathbb{P}$ .

Theorem 3 in Section 3 solves Problem 2. We prove for any given state on the boundary of $\mathcal{R}(\mathcal{P}_{0},t)$ that there exists a tight overapproximation touching $\mathcal{R}(\mathcal{P}_{0},t)$ on this state. Then, we can describe the reachable set as the intersection of all the possible overapproximations.

2 Overapproximation with Paraboloids

In this section, Problem 1 is solved using time-varying paraboloids $P:I\rightarrow\mathbb{P}$ where $I$ is the interval of definition of $P$ . Time-varying parameters $(\mathbold{E},\mathbold{f},\mathbold{g})$ of $P$ satisfies a differential equation that guarantees an overapproximation relationship with the reachable set, i.e. $\mathcal{R}(\mathcal{P}_{0},t)\subseteq P(t)$ for any $t\in I$ . We express existence and domain of definition $I$ of the time-varying paraboloid $P$ . We prove that the overapproximations $P$ are tight since there are touching trajectories of $\mathcal{R}(\mathcal{P}_{0},t)$ that both belong to the surface of $P(t)$ and to the surface of $\mathcal{R}(\mathcal{P}_{0},t)$ for $t\in I$ . Finally, the method is presented for a simple toy example.

Parameters of $P$ are expressed as solutions to an initial value problem. For given $E_{0}\in\mathbb{S}^{n}$ , let $\mathbold{E}$ be the solution to the following Differential Riccati Equation (DRE) with initial condition $\mathbold{E}(0)=E_{0}$ :

[TABLE]

Let $T_{E}(E_{0})\in\mathbb{R}^{+}\cup\{+\infty\}$ be defined for the initial condition $E_{0}$ s.t. $\left[0,T_{E}(E_{0})\right[$ is the interval of definition of the solution of (8) (existence, uniqueness, convergence properties and continuity of the solution are studied in [29]). Let $\mathbold{f}$ denote the solution to the following IVP with initial condition $\mathbold{f}(0)=f_{0}$ :

[TABLE]

$\mathbold{f}$ satisfies a linear varying parametric differential equation with a continuous input signal. On $\left[0,T_{E}(E_{0})\right[$ , solution $\mathbold{f}$ to (9) exists, is unique and continuous. By continuity of $\mathbold{f}$ and $\mathbold{u}$ over $\left[0,T_{E}(E_{0})\right[$ , $\mathbold{g}$ is defined on $\left[0,T_{E}(E_{0})\right[$ . For $t\in\left[0,T_{E}(E_{0})\right[$ , let:

[TABLE]

where

[TABLE]

Definition 4 (Time-Varying Paraboloid).

For an initial paraboloid $\mathcal{P}_{0}\in\mathbb{P}$ , let the time-varying paraboloid $P$ be defined as:

[TABLE]

where the time-varying coefficients $(\mathbold{E},\mathbold{f},\mathbold{g})$ are solutions to (8,9,10) with initial condition $\mathcal{P}(E_{0},f_{0},g_{0})=\mathcal{P}_{0}$ . Let $\mathcal{T}$ be the function that associates to the initial paraboloid $\mathcal{P}_{0}\in\mathbb{P}$ the time-varying paraboloid $P$ . Let $T_{P}(P)=T_{E}(E_{0})$ and $\mathcal{I}(P)=\left[0,T_{P}(P)\right[$ be the interval of definition of $P$ .

For $P=\mathcal{T}(\mathcal{P}_{0})$ , let $h(t,\cdot)$ be the value function of $P(t)$ at $t\in\mathcal{I}(P)$ . For $z_{t}=(x_{t},x_{q,t})\in\mathbb{R}^{n+1}$ , $w_{t}\in\mathbb{R}^{m}$ , $h_{z}(t)=h(t,\mathbold{z}(t))$ is the value function along the trajectory $\mathbold{z}=(\mathbold{x},\mathbold{x_{q}})$ solution to (1, 4) generated by $\mathbold{w}$ such that $\mathbold{w}(t)=w_{t}$ and $\mathbold{z}(t)=z_{t}$ .

Property 2.

The maximum time derivative of the value function $h(t,\mathbold{z}(t))$ along the trajectories $\mathbold{z}$ for a disturbance $w_{t}$ exists for all $t\in\mathcal{I}(P)$ and it is equal to zero.

Proof 2.

At a time $t\in\mathcal{I}(P)$ , it holds

[TABLE]

and the time derivative of $h_{z}$ at $t$ is

[TABLE]

Therefore, $k(w_{t})=\dot{h}_{z}(t)$ is a quadratic function of $w_{t}$ :

[TABLE]

where $\mathbold{d}$ is function of $\mathbold{E}(t)$ , $\mathbold{f}(t)$ , $\mathbold{g}(t)$ , $\mathbold{x}(t)$ and $\mathbold{u}(t)$ and system parameters. Since $M_{w}\prec 0$ , the supremum of $w_{t}\mapsto k(w_{t})$ exists and is attained for $w_{t}=\mathbold{{w^{*}}}(t)=\operatorname*{arg\,max}_{w_{t}\in\mathbb{R}^{m}}k(w_{t})$ with:

[TABLE]

Since $M_{w}\prec 0$ , $M_{w}^{\scalebox{0.6}{-1}}$ is well defined. Using (8,9,10) in (11), we get $\max_{w_{t}\in\mathbb{R}^{m}}\dot{h}_{z}=0$ . Therefore, $\dot{h}_{z}\leq 0$ for any $w_{t}\in\mathbb{R}^{m}$ .

We can now state one of our main results:

Theorem 1 (Solution to Problem 1)

Let $P=\mathcal{T}(\mathcal{P}_{0})$ for a set of initial states $\mathcal{P}_{0}$ . For all $t\in\mathcal{I}(P)$ , the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ of $\mathscr{S}(\mathcal{P}_{0},t)$ , is overapproximated by $P(t)$ , i.e.:

[TABLE]

Proof 3.

Using Property 2, by integration of $\dot{h}_{z}$ , if $h_{z}(0)\leq 0$ then $\forall t\in\mathcal{I}(P),h_{z}(t)\leq 0$ , i.e.:

[TABLE]

The constraint (5) ensures that $\mathbold{z}(t)\in\mathcal{Z}_{+}$ .

Property 3.

Let $\mathbold{z^{*}}$ be a trajectory generated by $\mathbold{{w^{*}}}$ defined in (12) such that initial condition satisfies $\mathbold{z^{*}}(0)\in\partial\mathcal{P}_{0}$ . At any time $t\in\mathcal{I}(P)$ , it holds $\mathbold{z^{*}}(t)\in\partial P(t)$ .

Proof 4.

$\mathbold{z^{*}}$ * is the trajectory generated by the optimal disturbance $\mathbold{{w^{*}}}$ . Using Property 2, $\dot{h}_{z^{*}}(t)=0$ for any $t\geq 0$ . Since $h_{z^{*}}(0)=0$ , by integration, $h_{z^{*}}(t)=0$ .*

Trajectories generated by $\mathbold{{w^{*}}}$ defined in (12) stay in contact with the surface of their time-varying paraboloids. Touching trajectories of $P$ do not necessarily belong to $\mathscr{S}(\mathcal{P}_{0},t)$ , $t\in\mathcal{I}(P)$ , as the energetic constraint might be locally violated.

Remark 1 (Representation of paraboloids).

In [21], the time-varying value function is a quadratic function defined by its quadratic coefficient $\mathbold{S}$ , its center $\mathbold{x_{c}}$ and its value at the center $\mathbold{\rho}$ , all satisfying an IVP. In [21], the center $\mathbold{x_{c}}$ can diverge when the determinant of $\mathbold{S}$ vanishes. However, the corresponding time-varying value function is time-continuous and can be extended continuously. In this paper, we choose to work with variables $\mathbold{E}$ , $\mathbold{f}$ and $\mathbold{g}$ (see Definition 1) to avoid this issue.

Example 1.

Let $A=-1$ , $B=1$ , $M=\left[\begin{smallmatrix}1&0&0\\ 0&1&0\\ 0&0&-2\end{smallmatrix}\right]$ , $B_{u}=0$ and $\mathbold{u}:\left[0,\infty\right[\mapsto 0$ . Solutions to IVP (8) (that is $\dot{\mathbold{E}}=-\frac{1}{2}\mathbold{E}^{2}+2\mathbold{E}-1$ ) diverge for $E_{0}\prec E^{-}$ (see Figure 1) where $E^{-} \prec E^{+}$ are the roots of the equation $-\frac{1}{2}E^{2}+2E-1=0$ for $E\in\mathbb{R}$ , $E^{-}=2-\sqrt{2}$ and $E^{+}=2+\sqrt{2}$ .

Figure 2 shows the trajectory of the paraboloid for $E_{0}$ in the stable region $E_{0}\succ E^{-}$ while Figure 3 shows the trajectory of the paraboloid for $E_{0}$ in the unstable region $E_{0}\prec E^{-}$ .

2.1 Domain of definition of $P$

By Definition 4, the domain of a time-varying paraboloid $P$ is the domain of its quadratic time-varying coefficient $\mathbold{E}$ . Since the solution of the DRE (8) might diverge in a finite time $T_{E}(E_{0})<\infty$ (where $E_{0}$ is the initial condition of (8)), $P$ is defined only in the right-open interval $\left[0,T_{E}(E_{0})\right[$ . In this part, we show that since the touching trajectories of $P$ are defined over the closed interval $\left[0,T_{E}(E_{0})\right]$ ,the definition of $P$ can be prolongated to the same closed interval.

Property 2 and (8,9,10) can be derived solving the following optimal control problem (for $t>0$ ):

[TABLE]

for given $(x_{t},x_{q,t})\in\mathcal{Z}_{+}$ . This is a special instance of the LQR problem (see e.g. [21]). For $x\in\mathbf{L}_{2}(T;\mathbb{R}^{n})$ a touching trajectory, let

[TABLE]

be the normal to the paraboloid surface. Using (8,9,10), $\mathbold{n}$ satisfies the following differential equation:

[TABLE]

where

[TABLE]

and

[TABLE]

The value function evaluated along the touching trajectory $\mathbold{x}$ is then obtained by introducing the parameter $\mathbold{r}=\mathbold{g}-\mathbold{f}^{\top}\mathbold{x}$ which satisfies:

[TABLE]

with

[TABLE]

and

[TABLE]

The value function is then equal to:

[TABLE]

Let the time-varying paraboloid $P=\mathcal{T}(\mathcal{P}_{0})$ diverge in finite time, i.e. $T_{P}(\mathcal{P}_{0})<\infty$ . Since all the touching trajectories are continuous in time, each touching trajectory is defined over $[0,T_{P}(\mathcal{P}_{0})]$ . Their corresponding value function $h$ evaluated along the touching trajectory is as well continuous over $[0,T_{P}(\mathcal{P}_{0})]$ . Therefore, one can prolongates the definition of $P$ until $T_{P}(\mathcal{P}_{0})$ by continuity of the value function:

[TABLE]

where $T=T_{P}(\mathcal{P}_{0})$ . We state this result in the following property:

Property 4.

For any $P=\mathcal{T}(\mathcal{P}_{0})$ , if the quadratic coefficient of the time-varying paraboloid set $P$ diverges in finite time, then the prolongation to the right of $P$ is defined:

[TABLE]

3 Exact Reachable Set

In this section, we first show that the state constraint (2) can be used to redefine the time-varying paraboloids (in Sections 3.1 and 3.2). Then we define a set of time-varying paraboloid (in Section 3.3). At each time instant, the intersection of these paraboloids is an overapproximation of the reachable set (in Section 3.4). Finally, we prove that when some topological assumption holds about the reachable set, our overapproximation is equal to the reachable set (in Sections 3.5, 3.6 and 3.7).

3.1 Scaled Paraboloids

Property 1 ensures that the overapproximation relationship is still valid when we scale the paraboloid. In this section, a definition of time-varying paraboloid with continuous time scaling is given.

For a given scaling function $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ (non-negative and locally square integrable function), an initial scaling factor ${\gamma}_{0}\geq 1$ and given initial conditions $\mathcal{P}_{0}=(E_{0},f_{0},g_{0})\in\mathbb{P}$ , as in differential equations (8,9,10) in Section 2, we can similarly define the initial value problem:

[TABLE]

with

[TABLE]

The differential equation (13) is similar to (8,9,10) except that a multiplicator to the constraint is added to the value function and the initial condition of the time-varying paraboloid is scaled by ${\gamma}_{0}$ .

For $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ , ${\gamma}_{0}\geq 1$ , let $P=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0},\mathbold{{\gamma}})$ be the time-varying paraboloid with time-varying parameters defined by (13) for initial conditions defined by $\mathcal{P}_{0}$ . The Hamiltonian form of the equation can be as well defined and when the quadratic coefficient $\mathbold{E}$ diverges in finite time, the interval of definition of the time-varying paraboloids can be as well prolongated to the closed interval.

The worst disturbance is still expressed by (12) and Property 2 can be restated as:

Property 5.

For $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ , ${\gamma}_{0}\geq 1$ , and $\mathcal{P}_{0}\in\mathbb{P}$ , let $P=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0},\mathbold{{\gamma}})$ . For an optimal trajectory $\mathbold{{z}^{*}}$ generated by the disturbance $\mathbold{w}$ defined in (12) s.t. $\mathbold{{z}^{*}}(0)\in\partial{P}(0)$ , for any $t\geq 0$ it holds:

[TABLE]

Proof 5.

Direct derivation from (13).

When $\mathbold{{z}^{*}}(0)\in\partial{P}(0)$ , the solution to the ODE (15) is, for $t\in \mathcal{I}(P)$ :

[TABLE]

Since $\int_{0}^{T}\gamma(t)\mathbold{{x}_{q}^{*}}(t)dt$ might not be equal to [math], trajectories generated by the worst case disturbance $\mathbold{{w^{*}}}$ do not necessarily stay in contact with the time-varying paraboloid and therefore are not touching trajectories. For this reason, we call optimal trajectories the trajectories generated by $\mathbold{{w^{*}}}$ given in (12) .

Property 6.

Let $\mathbold{{z}^{*}}$ an optimal trajectory of ${P}$ s.t. $\mathbold{{z}^{*}}(0)\in\partial{P}(0)$ , if:

[TABLE]

for any $t\geq 0$ and $h_{z^{*}}(0)=0$ , then $\mathbold{{z}^{*}}$ is a touching trajectory of ${P}$ .

Proof 6.

Using (16).

Therefore, for any other trajectory of the constrained system $\mathscr{S}$ , $h$ is a decreasing function of time along the trajectory. In this case, Theorem 1 can be rewritten for continuously scaled time-varying paraboloid:

Theorem 2

For a set of initial states $\mathcal{P}_{0}$ , a scaling function $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ and an initial scaling factor ${\gamma}_{0}\geq 1$ , let $P=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0},\mathbold{{\gamma}})$ . The reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ of $\mathscr{S}(\mathcal{P}_{0},t)$ , $t>0$ , is overapproximated by $P(t)$ , i.e.:

[TABLE]

Proof 7.

By integration of (16) over the interval $[0,t]$ .

3.2 Definition domain of time-varying paraboloid

In the case where the DRE in (8) does not have any convergent solution for any positive definite initial condition (i.e. when system defined in Section 1.1 is unstable), works based on the LQR formulation of the IQC cannot overapproximate the reachable set for any $t\geq 0$ (for any initial condition of the DRE, the solution to the DRE is defined over a finite escape time). In this part, we show that for any system (whether it is stable or unstable) there is always a $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ such that the corresponding time-varying paraboloid is defined over $\mathbb{R}^{+}$ . This is one of the key advantages of our approach.

For a given positive definite initial condition $E_{0}\succeq 0$ , for the scaling factor

[TABLE]

over $\mathbb{R}^{+}$ and a scaling factor ${\gamma}_{0}=1$ , the corresponding solution $\mathbold{E}$ to the DRE (13a) does not diverge over $\mathbb{R}^{+}$ if:

[TABLE]

Since $-M_{w}\succ 0$ , the Schur complement of (17) leads to the equivalent non-negativity condition:

[TABLE]

where

[TABLE]

By choosing $\kappa$ such that

[TABLE]

then there is a convergent solution $\mathbold{E}$ to the DRE (13a). Therefore, for any given $E_{0}\succ 0$ , there exists a $\kappa>0$ such that (17) is satisfied.

Property 7.

There is a $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ such that $P=\mathcal{T}(\mathcal{P}_{0},1,\mathbold{{\gamma}})$ is defined over $\mathbb{R}^{+}$ .

By Property 1, ${P}(t)\cap\mathcal{Z}_{+}$ is an overapproximation of the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ for any $t\geq 0$ .

Remark 2.

At a given time, the reachable set of an IQC system is always bounded. Whether the DRE (8) has a solution or not over a given interval $[0,T]$ , $T>0$ , we can always bound the set of reachable states. This result is the main difference with other works in reachable set overapproximation for IQC systems (see [21, 22, 16, 23, 24]). Up to the knowledge of the authors, other works reachable set overapproximation for IQC systems only use results from Theorem 1. These results are dependent on the existence of a solution to the DRE. By taking into account the constraint over the interval of integration and not only at the terminal time, we have less conservative results.

Example 2.

Figure 4 shows plots of the reachable set of the unstable system $\mathscr{S}(\mathcal{P}_{0})$ defined by parameters:

[TABLE]

and a zero input signal $u$ . The set of initial states $\mathcal{Z}_{0}$ is a paraboloid $\mathcal{Z}_{0}=\mathcal{P}_{0}=\mathcal{P}(\mathbold{E}_{0},\mathbold{f}_{0},\mathbold{g}_{0})$ with:

[TABLE]

The solution to DRE (13a) for $\mathbold{{\gamma}}=0$ and ${\gamma}_{0}=0$ has a finite escape time and diverges at $T_{P}(P_{ns})=1.7$ . The solution to DRE (13a) for $\mathbold{{\gamma}}=0$ and ${\gamma}_{0}=1$ is defined over $\mathbb{R}^{+}$ .

3.3 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 Property 1 and preliminary results showing that for any state of the overapproximation, there exists a trajectory in $\mathscr{S}(\mathcal{P}_{0},t)$ , $t>0$ , leading to this state.

Let $\Pi$ be defined as follows:

[TABLE]

$\Pi$ corresponds to the set of all time-varying paraboloids with initial conditions $\mathcal{P}_{0}$ and generated by the set of non-negative scalings $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ and the set of initial scaling factors ${\gamma}_{0}\geq 1$ . Let

[TABLE]

the set of all the defined time-varying paraboloids at time $t\geq 0$ . By Property 7, $\Pi(t)\neq\emptyset$ for any $t\geq 0$ . In other words, $\Pi$ is defined over $\mathbb{R}^{+}$ as well.

For $t\geq 0$ , let

[TABLE]

the intersection of all the defined time-varying paraboloids $P$ of $\Pi$ at time $t$ (see Figure 5). Since $\Pi^{\cap}(\cdot)$ is defined over $\mathbb{R}^{+}$ , $\Pi^{\cap}(\cdot)$ is defined over $\mathbb{R}^{+}$ .

We now prove that, when some assumptions about the topology of $\Pi^{\cap}$ hold (Assumption 1 and 2), we have $\mathcal{R}(\mathcal{P}_{0},t)=\Pi^{\cap}(t)\cap\mathcal{Z}_{+}$ , for any $t\geq 0$ (Theorem 3, Section 3.7). To achieve that:

•

we prove the overapproximation relationship $\mathcal{R}(\mathcal{P}_{0},t)\subseteq\Pi^{\cap}(t)$ (Section 3.4);

•

we prove that any state $(x,x_{q})\in\Pi^{\cap}(t)$ is reachable from a state $(x,x_{q}^{\prime})\in\partial\Pi^{\cap}(t)$ with $x_{q}\leq x_{q}^{\prime}$ (Section 3.5);

•

for a state $z_{t}\in\partial\Pi^{\cap}(t)$ , we find a touching trajectory $\mathbold{z^{*}}=(\mathbold{x^{*}},\mathbold{x_{q}^{*}})$ of $\Pi^{\cap}$ such that $\mathbold{z^{*}}(t)=z_{t}$ . This touching trajectory $(\mathbold{x^{*}},\mathbold{x_{q}^{*}})$ of $\Pi^{\cap}$ satisfies the state constraint $\mathbold{x_{q}}(\cdot)\geq 0$ over $[0,t]$ (Section 3.6);

•

finally, we conclude that any $z_{t}\in\Pi^{\cap}(t)$ is reachable from $\mathcal{P}_{0}$ , thus $\mathcal{R}(\mathcal{P}_{0},t)=\Pi^{\cap}(t)\cap\mathcal{Z}_{+}$ (Section 3.7).

3.4 Overapproximation Relationship

Theorem 2 states that each time-varying paraboloid defined in Section 3.1 is an overapproximation of the reachable set. An intersection of many time-varying paraboloids is as well an overapproximation of the reachable set.

Property 8.

$\mathcal{R}(\mathcal{P}_{0},t)\subseteq\Pi^{\cap}(t)\cap\mathcal{Z}_{+}$ * for any $t\geq 0$ .*

Proof 8.

This is a direct consequence of Theorem 2 and (20).

Example 3 (Continued from Example 1).

In the case where the solution to (8) does not converge (i.e. $E_{0}<\mathbold{E}^{-}$ ), Figure 6 shows several paraboloid trajectories with different initial scaling factors. Scaling functions are equal to 0 and initial scaling factors $\gamma_{i}$ are greater than $1$ , $\mathcal{P}_{0}\cap\mathcal{Z}_{+}\subset\gamma_{i}\mathcal{P}_{0}\cap\mathcal{Z}_{+}$ . Therefore, each time-varying paraboloid is a valid constraint that bounds $\mathcal{R}(\mathcal{P}_{0},t)$ , $t\in\mathcal{I}(\Pi)$ (Theorem 1). Therefore, $\mathcal{R}(\mathcal{P}_{0},t)\subseteq P^{\cap}(t)=P_{0}(t)\cap P_{1}(t)\cap\dots\cap P_{4}(t)$ where $P_{i}=\mathcal{T}(\mathcal{P}_{0},\gamma_{i},\mathbold{0})$ , and $\gamma_{i}$ are resp. equal to $1$ , $1.6$ , $2.2$ , $2.7$ and $3.3$ for $i=0,\dots,4$ . In this case, the overapproximation $P^{\cap}(t)$ is strictly included in $P_{0}(t)$ .

Observations in Example 3 motivate the use of multiple time-varying paraboloids to get better overapproximations of the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ , $t>0$ .

3.5 Past trajectory for states in the overapproximation’s interior

Property 9 is already presented in [26, Property 7], we restate it hereby for completeness.

Property 9 shows that the state $(x,\alpha x_{q})$ is reachable from the given state $(x,x_{q})$ for any given $\alpha\in[0,1]$ .

Property 9.

For $t\geq 0$ , if $(x,x_{q})\in\mathcal{R}(\mathcal{P}_{0},t)$ then $(x,\alpha x_{q})\in\mathcal{R}(\mathcal{P}_{0},t)$ for all $\alpha\in[0,1]$ .

Proof 9.

Let $f:t,x\mapsto Ax+B\mathbold{w}(t)+B_{u}\mathbold{u}(t)$ . Since for any $(t,x)\in\mathbb{R}^{+}\times\mathbb{R}^{n}$ , $f(.,x)$ is locally measurable over $\mathbb{R}^{+}$ , $f(t,.)$ is Lipschitz over $\mathbb{R}^{n}$ , (1) has a unique solution $\mathbold{x}$ (see [30], Theorem 1.1) that is time-continuous. Therefore, for a trajectory $(\mathbold{x},\mathbold{x_{q}})\in\mathscr{S}(\mathcal{P}_{0},T)$ , $T>0$ , $\mathbold{x}$ is time-continuous.

For $\epsilon>0$ , let $\mathbold{w}\in\mathbf{L}_{2}([0,t+\epsilon];\mathbb{R}^{m})$ , s.t. $\mathbold{w}^{\top}(s)M_{w}\mathbold{w}(s)=-(1-\alpha)\mathbold{x_{q}}(t)\frac{1}{\epsilon}$ when $s\in[t,t+\epsilon]$ . Then

[TABLE]

when $\epsilon\rightarrow 0$ . Using Cauchy-Schwartz inequality:

[TABLE]

and the time-continuity of $\mathbold{x}$ , the quantity

[TABLE]

when $\epsilon\rightarrow 0$ . By integration, $\mathbold{x}_{q}(t+\epsilon)\rightarrow\alpha\mathbold{x_{q}}(t)$ when $\epsilon\rightarrow 0$ . Since $\mathbold{x}$ is time-continuous, $\mathbold{x}(t+\epsilon)\rightarrow\mathbold{x}(t)$ when $\epsilon\rightarrow 0$ . By continuity of $\mathbold{u}$ , $\mathbold{x}$ and $\mathbold{w}$ over $[t,t+\epsilon]$ , $\mathbold{x_{q}}$ is continuous over $[t,t+\epsilon]$ . Then, there exists a $t^{\prime}\in[t,t+\epsilon]$ such that $\mathbold{x_{q}}(\tau)\geq\alpha\mathbold{x_{q}}(t)\geq 0$ for all $\tau\in[t,t^{\prime}]$ and $\mathbold{x_{q}}(t^{\prime})\rightarrow\alpha\mathbold{x_{q}}(t)$ when $\epsilon\rightarrow 0$ . Therefore, the constraint $\mathbold{x_{q}}(\cdot)\geq 0$ is satisfied over $[t,t^{\prime}]$ and the trajectory $(\mathbold{x},\mathbold{x_{q}})$ is a valid trajectory of $\mathscr{S}(\mathcal{P}_{0},t^{\prime})$ for all $t\leq t^{\prime}$ .

3.6 Past trajectory for states in the overapproximation’s boundary

In this section, touching trajectories of $\Pi^{\cap}$ are identified. We show that all these touching trajectories satisfies the state constraint (5).

The value function $\tilde{h}$ of a time-varying paraboloid $\tilde{P}\in\Pi$ can be approximated at the first order along a touching trajectory $\mathbold{z^{*}}$ of another time-varying paraboloid ${P}\in\Pi$ when their scaling functions $\mathbold{\tilde{\gamma}}$ and $\mathbold{{\gamma}}$ and initial scaling factor $\tilde{\gamma}_{0}$ and ${\gamma}_{0}$ are close. In this part, we compute this first order approximation when $\mathbold{\tilde{\gamma}}=\mathbold{{\gamma}}+\mathbold{\delta}$ and $\tilde{\gamma}_{0}={\gamma}_{0}+\delta_{0}$ for small variations $\mathbold{\delta}\in\mathbf{L}_{2}(\mathbb{R}^{+};\mathbb{R})$ and $\delta_{0}\in\mathbb{R}$ (i.e. when $\left\lVert\mathbold{\delta}\right\rVert+\left|\delta_{0}\right|$ tends to [math]).

To prove that the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ , $t>0$ , is exactly described by $\Pi^{\cap}(t)$ , we show that for any optimal trajectory $\mathbold{{z}^{*}}$ of ${P}\in\Pi^{\cap}$ s.t. $\mathbold{{z}^{*}}$ is violating the constraint (5), there is a $\tilde{P}\in\Pi$ such that the end point $\mathbold{{z}^{*}}(t)$ does not belong to $\tilde{P}(t)$ and therefore to $\Pi^{\cap}(t)$ . To do so, we will study the value function of a time-varying paraboloid $\tilde{P}$ for touching trajectories of ${P}$ .

Property 10.

For $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+},\mathbb{R}^{+})$ and ${\gamma}_{0}\geq 1$ , let the corresponding time varying paraboloid ${P}=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0}, \mathbold{{\gamma}})$ . For any $t$ in the open set of $\mathcal{I}({P})$ , it exists $\epsilon>0$ and $H>0$ s.t. for any $\mathbold{\delta}\in\mathbf{L}_{2}(\mathbb{R}^{+};\mathbb{R})$ , $\left\lVert\mathbold{\delta}\right\rVert\leq\epsilon$ , for any $\delta_{0}\in\mathbb{R}$ , $\left|\delta_{0}\right|\leq\epsilon$ , $\tilde{P}=\mathcal{T}(\mathcal{P}_{0},\tilde{\gamma}_{0},\mathbold{\tilde{\gamma}})$ where $\mathbold{\tilde{\gamma}}=\mathbold{{\gamma}}+\mathbold{\delta}$ and $\tilde{\gamma}_{0}={\gamma}_{0}+\delta_{0}$ s.t. $t$ belongs to the open set of $\mathcal{I}(\tilde{P})$ , let $\tilde{h}_{t}$ the value function of $\tilde{P}(t)$ and $\mathbold{{z}^{*}}=(\mathbold{{x}^{*}},\mathbold{{x}_{q}^{*}})$ an optimal trajectory of ${P}$ , it holds:

[TABLE]

where

[TABLE]

and

[TABLE]

Proof 10.

Let $({\mathbold{E}},{\mathbold{f}},{\mathbold{g}})$ and $(\tilde{\mathbold{E}},\tilde{\mathbold{f}},\tilde{\mathbold{g}})$ (resp.) be parameters of ${P}$ and $\tilde{P}$ (resp.), and

[TABLE]

Using (1, 12, 13), $\mathbold{\nu}$ satisfies the linear time varying differential equation:

[TABLE]

with $\mathbold{n}={\mathbold{E}}\mathbold{{x}^{*}}-\mathbold{f}$ and $A_{\nu}(\tau)=-A^{\top}+M_{xw}M_{w}^{\scalebox{0.6}{-1}}B^{\top}+\tilde{\mathbold{E}}(\tau)BM_{w}^{\scalebox{0.6}{-1}}B^{\top}+\mathbold{{\gamma}}(\tau)I$ . By (14), initial values of ${P}$ and $\tilde{P}$ satisfies

[TABLE]

therefore $\mathbold{\nu}(0)=\delta_{0}(E_{0}\mathbold{{x}^{*}}(0)-f_{0})$ where $(E_{0},f_{0},g_{0})=\mathcal{P}_{0}$ . Since $t$ belongs to the open set of $\mathcal{I}(\tilde{P})$ , $\tilde{\mathbold{E}}(\cdot)$ is bounded over $[0,t]$ (the discontinuity of $\tilde{\mathbold{E}}$ can only occur at the final integration time). By time-continuity of $\tilde{\mathbold{E}}(\cdot)$ over $[0,t]$ , there is a scalar $K>0$ that bounds $\left\lVert\tilde{\mathbold{E}}(\cdot)\right\rVert$ over $[0,t]$ . Then, since $\mathbold{{\gamma}}$ is measurable, there exists a measurable function $L\in\mathbf{L}_{2}([0,t];\mathbb{R}^{+})$ such that:

[TABLE]

over $\tau\in[0,t]$ . We can integrate (22):

[TABLE]

Since $l\mapsto\int_{0}^{l}\left|\mathbold{\delta}(\tau)\right|\left\lVert\mathbold{n}(\tau)\right\rVert d\tau$ is a non decreasing function over $[0,t]$ , by applying the Grönwall inequality, we get:

[TABLE]

Let

[TABLE]

$\mathbold{{z}^{*}}$ * is an optimal trajectory of ${P}$ s.t. $\mathbold{{z}^{*}}(t)\in\partial{P}(t)$ , therefore, ${h}_{t}(\mathbold{{z}^{*}}(t))=0$ , therefore, using (15):*

[TABLE]

Using (1, 12, 13, 25), $q$ satisfies:

[TABLE]

Using (24) and (25):

[TABLE]

where

[TABLE]

with initial condition $q(0)={h}_{0}(\mathbold{{z}^{*}}(0))-\tilde{h}_{0}(\mathbold{{z}^{*}}(0))$ . Since $\mathbold{{z}^{*}}$ is a touching trajectory of ${P}$ , it holds ${h}_{0}(\mathbold{{z}^{*}}(0))=0$ , therefore:

[TABLE]

Therefore, $\tilde{h}_{0}(\mathbold{{z}^{*}}(0))$ satisfies:

[TABLE]

and $q(0)=\frac{\delta_{0}}{{\gamma}_{0}}\mathbold{{x}_{q}^{*}}(0)$ .

Since $t$ belongs to the open set of $\mathcal{I}(\tilde{P})$ , the optimal trajectory $\mathbold{{z}^{*}}$ and $\mathbold{\nu}$ are defined and continuous over $[0,t]$ . Moreover, since $\mathbold{{\gamma}}$ and $\mathbold{\delta}$ are measurable over $[0,t]$ , the solution to the linear time-varying equation (26) exists over $[0,t]$ and is:

[TABLE]

Then, using (23):

[TABLE]

with

[TABLE]

a finite constant where $R=\left\lVert BM_{w}^{\scalebox{0.6}{-1}}B^{\top}\right\rVert$ , $K=\exp{\int_{0}^{t}L(\tau)d\tau}$ and $N=\int_{0}^{t}e^{\int_{s}^{t}(\mathbold{{\gamma}}(r)+\mathbold{\delta}(r))dr}ds$ . This ends the proof.

Remark 3.

When $\mathbold{{\gamma}}=\mathbold{0}$ , Property 10 matches with [26, Property 9] where their initial scaling correspond to our initial scaling factor ${\gamma}_{0}$ .

Property 11 gives conditions where the sign of $\tilde{h}_{t}(\mathbold{{z}^{*}}(t))$ is only determined by its first order approximation defined in Property 10.

Property 11.

Let $\mathbold{z^{*}}$ a touching trajectory of ${P}=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0},\mathbold{{\gamma}})$ for $\mathbold{{\gamma}}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R}^{+})$ , ${\gamma}_{0}\geq 1$ given and $t\in\mathcal{I}({P})$ given. If there is a $\mathbold{\delta}\in\mathbf{L}_{2,loc}(\mathbb{R}^{+};\mathbb{R})$ and a $\delta_{0}\in\mathbb{R}$ , s.t. $\left\lVert\mathbold{\delta}\right\rVert\leq\epsilon$ and $\left|\delta_{0}\right|\leq\epsilon$ and $t\in\mathcal{I}(\tilde{P})$ (where $\tilde{P}=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0}+\delta_{0},\mathbold{{\gamma}}+\mathbold{\delta})$ ) and

[TABLE]

then the sign of

[TABLE]

is equal to the sign of $\tilde{h}_{t}(\mathbold{{z}^{*}}(t))$ where $\tilde{h}_{t}$ is the value function of $\tilde{P}(t)$ and $H>0$ defined in (27) and

[TABLE]

Proof 11.

This is a direct consequence of Property 10 and of the property: $(\left|a-b\right|\leq c)\land(c<\left|b\right|)\Rightarrow\mathrm{sign}(a)=\mathrm{sign}(b)$ for $a,b,c\in\mathbb{R}$ .

Provided the existence of a $(\mathbold{\delta},\delta_{0})\in\mathbf{L}_{2}(\mathbb{R}^{+};\mathbb{R})\times\mathbb{R}$ such that $\mathbold{{\gamma}}+\mathbold{\delta}\geq 0$ and ${\gamma}_{0}+\delta_{0}\geq 1$ , the first order approximation of the value function of $\tilde{P}=\mathcal{T}(\mathcal{P}_{0},{\gamma}_{0}+\delta_{0},\mathbold{{\gamma}}+\mathbold{\delta})$ gives a way to identify time varying paraboloids $\tilde{P}$ that belongs to $\Pi$ such that an invalid trajectory with an end state $z_{t}\in\partial P(t)$ (meaning with initial state outside of the initial set $\mathcal{P}_{0}$ or a trajectory violating the constraint) does not belongs to $\tilde{P}(t)$ and therefore, does not belong to $\Pi^{\cap}(t)$ .

Property 12 states that the touching trajectories of $\Pi^{\cap}$ satisfy the state constraint (5). Property 12 is proven by choosing a valid trajectory candidate. If this trajectory violates the state constraint (5), then Property 11 provides a proof that this trajectory does not belongs to the overapproximation $\Pi^{\cap}$ .

Property 12.

For ${P}\in\Pi^{\cap}$ , if $z_{t}\in\partial\Pi^{\cap}(t)$ and $z_{t}\in\partial{P}(t)$ for $t$ in the open set of $\mathcal{I}({P})$ , then the optimal trajectory $\mathbold{{z}^{*}}$ of ${P}$ such that $\mathbold{{z}^{*}}(t)=z_{t}$ is a valid touching trajectory of $P$ and $\Pi^{\cap}$ .

Proof 12.

Let $\psi:\mathbb{R}^{+}\mapsto\mathbb{R}$ defined for $s\geq 0$ by:

[TABLE]

Let $\tau\in[0,t]$ and $I=[\tau,t]$ .

•

Case 1, $\mathbold{{x}_{q}^{*}}(0)<0$ **: with $\delta_{0}>0$ , using Property 11, $\mathbold{{z}^{*}}(t)\notin\tilde{P}(t)$ where $\tilde{P}\in\Pi$ since $\delta_{0}+{\gamma}_{0}\geq 1$ , so $\mathbold{{z}^{*}}(t)\notin\Pi^{\cap}(t)$ .

•

Case 2, $\psi(\cdot)<0$ over $I$ **: any $\mathbold{\delta}(\cdot)\geq 0$ over $I$ and $\mathbold{\delta}(\cdot)=0$ elsewhere such that $\int_{0}^{t}\mathbold{\delta}(s)\psi(s)ds\neq 0$ and for $\delta_{0}=0$ , using Property 11, $\mathbold{{z}^{*}}(t)\notin\tilde{P}(t)$ where $\tilde{P}\in\Pi$ since $\mathbold{{\gamma}}+\mathbold{\delta}\geq 0$ , so $\mathbold{{z}^{*}}(t)\notin\Pi^{\cap}(t)$ .

•

Case 3, $\psi(\cdot)>0$ over the open of $I$ and there is a $l\in I$ , s.t. $\int_{s\in[\tau,l]}\mathbold{{\gamma}}(s)\mathbold{{x}_{q}^{*}}(s)ds\neq 0$ **: since $\mathbold{{\gamma}}\geq 0$ , there exists a $\mathbold{\delta}\leq 0$ such that $\mathbold{{\gamma}}+\mathbold{\delta}\geq 0$ and for $\delta_{0}=0$ , using Property 11, $\mathbold{{z}^{*}}(t)\notin\Pi^{\cap}(t)$ .

•

Case 4, $\psi(\cdot)=0$ over $I$ **: since $\mathbold{{x}_{q}^{*}}$ is continuous over $\mathbb{R}^{+}$ , and since $\mathbold{{\gamma}}$ is locally square integrable, $\psi(t)=0\Rightarrow\mathbold{{x}_{q}^{*}}(t)=0$ , therefore $\mathbold{{x}_{q}^{*}}(\cdot)=0$ over $I$ . Consequently ${h}_{\tau}(\mathbold{{z}^{*}}(\tau))=0$ for $\tau\in I$ .

Cases 1 to 4 show that for $\mathbold{{z}^{*}}(t)\in\partial\Pi^{\cap}(t)$ :

•

either $\forall l\in I,\int_{\tau}^{l}\mathbold{{\gamma}}(\tau)\mathbold{{x}_{q}^{*}}(\tau)d\tau=0$ and $\psi(\cdot)=\mathbold{{x}_{q}^{*}}(s)>0$ ;

•

nor $\mathbold{{x}_{q}^{*}}(l)=0$ for $l\in I$ .

*Let a partition $[0,t]=\bigcup_{i\in\mathbb{N}}I_{i}$ such that over each open interval $I_{i}$ , $\psi(\cdot)\bowtie_{i}0$ with $\bowtie_{i}\in\{<,>,=\}$ . We deduce that for any $s\in[0,t]$ : $\mathbold{{x}_{q}^{*}}(s)\geq 0$ and $\int_{I}\mathbold{{\gamma}}(\tau)\mathbold{{x}_{q}^{*}}(\tau)d\tau=0$ . $\mathbold{{z}^{*}}$ is a valid trajectory, i.e. the constraint (5) is satisfied. By 15, $\mathbold{{z}^{*}}$ is a touching trajectory of $P$ . Moreover, since $\mathbold{{z}^{*}}(0)\in\mathcal{P}_{0}$ , $\mathbold{{z}^{*}}$ is as well a touching trajectory of $\Pi^{\cap}$ . *

Since $\Pi^{\cap}(t)$ is an intersection of closed sets, $\Pi^{\cap}(t)$ is closed as well. In the general case, for an infinite intersection $\mathcal{Y}^{\cap}=\bigcap_{i\in\mathbb{N}}Y_{i}$ of closed sets $Y_{i}$ , $i\in\mathbb{N}$ , any boundary point $y\in\partial\mathcal{Y}^{\cap}$ does not necessarily belongs to the boundary of any $Y_{i}$ , $i\in\mathbb{N}$ (e.g. $\bigcap_{\epsilon\in]1,2]}[-\epsilon,\epsilon]=[-1,1]$ , but there is no $\epsilon\in]1,2]$ such that $1\in\partial[-\epsilon,\epsilon]$ ). The following assumption states that for every state on the boundary of the overapproximation $\Pi^{\cap}(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 1.

For any $z_{t}\in\partial\Pi^{\cap}(t)$ , there is a ${P}\in\Pi$ such that $z_{t}\in\partial P(t)$ .

This assumption is not a strong one and has been proved for simpler cases (see [26, Property 11]).

In Property 12, the existence of $\mathbold{\tilde{\gamma}}$ and ${\gamma}_{0}$ is conditioned by $t$ belonging to the open domain $\mathcal{I}(\tilde{P})$ ; to ensure this, $\left\lVert\mathbold{E}(\cdot)\right\rVert$ is assumed to be bounded over $[0,T]$ (by considering the case where $t$ is in the open set of $\mathcal{I}({P})$ ). In the general case, the boundedness of $\left\lVert\mathbold{E}(\cdot)\right\rVert$ is not granted (see the unstable case in Example 1 and Figure 3). Assumption 2 states that for any state on the boundary of the overapproximation $\Pi^{\cap}(t)$ , $t>0$ , there is neighbor state on the boundary of $\tilde{P}(t)$ where $\tilde{P}$ is a time-varying paraboloid of $\Pi$ not diverging at $t$ (i.e. $t$ belongs to the interior of $T_{P}(P)$ ).

Assumption 2.

For $t>0$ , for all $\epsilon>0$ , for any $z_{t}\in\partial\Pi^{\cap}(t)$ such that $z_{t}\in\partial{P}(t)$ , ${P}\in\Pi$ with ${P}$ unbounded, there is a $\tilde{z}_{t}$ that belongs to the boundary of $\tilde{P}(t)$ , $\tilde{z}_{t}\in\partial\tilde{P}(t)$ , such that $\left\lVert z_{t}-\tilde{z}_{t}\right\rVert<\epsilon$ .

Lemma 1 shows that any state $z_{t}\in\partial\Pi^{\cap}(t)$ (with $t\in\mathcal{I}(\Pi)$ given) is the terminal state of a touching trajectory $\mathbold{{z}^{*}}$ of $\Pi^{\cap}$ with initial state $\mathbold{{z}^{*}}(0)\in\partial\Pi^{\cap}(0)$ .

Lemma 1.

If Assumptions 1 and 2 hold, any state $z_{t}\in\partial\Pi^{\cap}(t)$ has a past touching trajectory $\mathbold{{z}^{*}}$ of $\Pi^{\cap}$ s.t. $\mathbold{{z}^{*}}(t)=z_{t}$ .

Proof 13.

Let $z_{t}\notin\mathcal{R}(\mathcal{P}_{0},t)$ such that for any $\epsilon>0$ , there is a touching trajectory $\mathbold{\tilde{z}}$ of $\tilde{P}$ , $\tilde{P}$ finite, with $\mathbold{\tilde{z}}(t)\in\mathcal{R}(\mathcal{P}_{0},t)$ and $\left\lVert\mathbold{\tilde{z}}(t)-z_{t}\right\rVert<\epsilon$ . For $\mathbold{{x}_{q}}(t)>0$ , we can define the optimal trajectory $\mathbold{{z}^{*}}$ with $\mathbold{{z}^{*}}(t)=\mathbold{{z}}_{t}$ . For any $\tau\in[0,t]$ , ${P}(\tau)$ is finite. Then, Property 12 can be used over $[0,\tau]$ . Therefore, if $z_{t}\in\partial{P}$ such that ${P}$ diverges at $t$ , it holds:

[TABLE]

For states not belonging to a diverging time-varying paraboloid, the property is a direct consequence of Assumption 1, Property 12.

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 (2) and with initial condition in the set of initial states.

3.7 Exact Reachable Set

We now state the main result of the paper:

Theorem 3 (Exact reachability, solution of Problem 2)

When Assumptions 1 and 2 hold, the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ of system $\mathscr{S}(\mathcal{P}_{0},t)$ (defined in Section 1.1) is equal to the set $\Pi^{\cap}$ defined in (20), namely

[TABLE]

for all $t\geq 0$ .

Proof 14.

Theorem 1 states that $\mathcal{R}(\mathcal{P}_{0},t)\subseteq\Pi^{\cap}(t)$ . By Property 9, for $z_{t}\in\Pi^{\cap}(t)$ , we can construct a trajectory $\mathbold{{z}}$ such that $\mathbold{{z}}(t)=z_{t}$ , $\mathbold{{z}}(t^{-})=z_{t}^{*}\in\partial\Pi^{\cap}(t)$ (Property 9). Since $z_{t}^{*}\in\partial\Pi^{\cap}(t)$ , using Lemma 1, there exists a trajectory $\mathbold{{z}}$ such that $\mathbold{{z}}(t^{-})=z_{t}^{*}$ and $\mathbold{{z}}$ is a touching trajectory of $\Pi^{\cap}$ on $\left[0,t\right[$ . Since $\mathbold{{z}}$ is a touching trajectory of $\Pi^{\cap}$ , $\mathbold{{z}}(0)\in\partial\Pi^{\cap}(0)$ with $\Pi^{\cap}(0)=\mathcal{P}_{0}=\mathcal{R}(0)$ . By Property 11, the trajectory $\mathbold{{z}}$ is valid (i.e. satisfies the energy constraint (5)) $z_{t}\in\mathcal{R}(\mathcal{P}_{0},t)$ .

4 Implementation

In this part, we discuss the practical implementation of the reachable sets overapproximation using Theorem 3. To do so, we compute a subset $\widetilde{\Pi}$ of $\Pi$ :

[TABLE]

$\widetilde{\Pi}$ corresponds to the time-varying paraboloid set generated by a finite subset of scaling functions and initial scaling. Then, the intersection of each time-varying paraboloid evaluated at a given $t>0$ is an overapproximation of the reachable set $\mathcal{R}(\mathcal{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 LABEL:algo:approx) that computes $\widetilde{\Pi}$ , its implementation in Matlab is available online [14].

Subset of scaling functions and initial scaling factors:

In this work, we choose to consider discrete scalings for the time-varying paraboloids. The scalings are applied at each $kT_{c}$ , for $T_{c}>0$ given and $k$ in $\mathbb{N}$ . A scaling is then described by a sequence of scaling factors $\{\lambda_{k}\}_{k\in\mathbb{N}}$ , $\lambda_{k}\geq 1$ , $k\in\mathbb{N}$ . The scaling functions are not used: $\mathbold{{\gamma}}(\cdot)=0$ .

In the ideal case, the scaling (function and factors) would be chosen such that the following property is verified:

[TABLE]

where $(\mathbold{{x}^{*}},\mathbold{{x}_{q}^{*}})$ corresponds to the touching trajectory associated with the scaling function $\mathbold{{\gamma}}$ and scaling factor ${\gamma}_{0}$ and such that $(\mathbold{{x}^{*}}(t),\mathbold{{x}_{q}^{*}}(t))=(x,x_{q})$ . In practice, since there might be an infinite number of states $(x,x_{q})$ verifying $\dot{x}_{q}\geq 0$ , only a finite number of states are checked. These states are chosen as projections of a given point in given directions over $\mathcal{Z}_{*}\cap\partial\widetilde{\Pi}^{\cap}$ . These points are then used to evaluate a range of scaling factors ${\gamma}_{0}$ to enforce $\dot{\mathbold{{x}}_{q}^{*}}(kT_{c})\geq 0$ . $\mathbold{{\gamma}}$ is not used.

Paraboloid numerical integration:

Let two paraboloids ${P}=\mathcal{T}(\mathcal{P}_{0},1,\mathbold{{\gamma}})$ , $\tilde{P}=\mathcal{T}(\mathcal{P}_{0},1,\mathbold{\tilde{\gamma}})$ . If $\mathbold{{\gamma}}(.)=\mathbold{\tilde{\gamma}}(.)$ over an interval $[0,t_{i}]$ , $t_{i}>0$ , then ${P}(.)=\tilde{P}(.)$ over $[0,t_{i}]$ . Let $t_{i}\geq 0$ corresponds to the maximal time instant where there is $\tilde{P}\in\Pi$ such that $\tilde{P}|_{[0,t_{i}]}={P}|_{[0,t_{i}]}$ (i.e. such that the restriction of $\tilde{P}$ on $[0,t_{i}]$ is equal to the one of ${P}$ on the same interval). And let $t_{f}\geq 0$ corresponds 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}]\subseteq[0,T]$ .

Since $\mathbold{{\gamma}}(.)=0$ over $\left]kT_{c},(k+1)T_{c}\right[$ , for any $\tau\in[0,T_{c}]$ , $k \in\mathbb{N}$ . Over $\left]kT_{c},(k+1)T_{c}\right[$ , the IVP (13) is then equivalent to the ”unscaled” IVP (8,9,10) for $(\mathbold{E}_{k},\mathbold{f}_{k},\mathbold{g}_{k})$ over each time interval $\left[kT_{c},(k+1)T_{c}\right]$ The solution to ${P}=\mathcal{T}(\mathcal{P}_{0},1,\mathbold{{\gamma}})$ is then described by parameters $(\mathbold{E},\mathbold{f},r)$ with

[TABLE]

for each $t\in[kT_{c},(k+1)T_{c}]$ , $k\in\mathbb{N}$ .

Cardinal limitation of $\widetilde{\Pi}$ :

In order to have a tractable integration of the reachable set computation, we limit the cardinality of $\widetilde{\Pi}$ in the following way:

•

at each time step $kT_{c}$ , we consider only $N_{new}$ scaled paraboloids of highest scaling factor;

•

$\widetilde{\Pi}$ 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 only consider elements of $\Pi$ that are more stable. Since for 2 solutions ${\mathbold{E}}$ and $\tilde{\mathbold{E}}$ of (8) respectively defined over $[0,{{T}}]$ and $[0,{\tilde{T}}]$ where ${{T}},{\tilde{T}}\in\mathbb{R}\cup\{\infty\}$ , if ${\mathbold{E}}(0)\preceq\tilde{\mathbold{E}}(0)$ , then ${\mathbold{E}}(t)\preceq\tilde{\mathbold{E}}(t)$ for $t$ in the interval of definition of ${\mathbold{E}}$ and $\tilde{\mathbold{E}}$ , we have ${{T}}\geq{\tilde{T}}$ (these property follow 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. A direct integration of the DRE (8) does not produce good results in practice (see [31]). Experiments presented in this works make use of the Chandrasekhar method [13]. This method integrates the Ordinary Differential Equation (ODE) (8) $\mathbold{E}$ using an intermediate ODE over the time-dependent matrix $\mathbold{L}$ in $\mathbf{L}_{2}(\mathbb{R}^{+},\mathbb{R}^{n\times n})$ :

[TABLE]

with

[TABLE]

where $\dot{\mathbold{E}}_{0}=\dot{\mathbold{E}}(0)$ given by (1). Then $\mathbold{E}$ is a solution to (8).

Since $\mathbold{L}(t)\mathbold{L}(t)^{\top}\succeq 0$ , this method is only applicable to strictly increasing solutions of the DRE. As seen in Example 1, the solutions to ODE (8) 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 $\mathbold{L},\mathbold{K}\in\mathbf{L}_{2}(\mathbb{R}^{+},\mathbb{R}^{n\times n})$ such as:

[TABLE]

with

[TABLE]

where $\dot{\mathbold{E}}_{0}=\dot{\mathbold{E}}(0)=\dot{\mathbold{E}}_{0}^{+}+\dot{\mathbold{E}}_{0}^{-}$ given by (8), with $\dot{\mathbold{E}}_{0}^{+}\succeq 0$ and $\dot{\mathbold{E}}_{0}^{-}\preceq 0$ The increasing and decreasing parts of $\mathbold{E}$ are respectively represented by the terms $\mathbold{L}$ and $\mathbold{K}$ . Our Chandrasekhar inspired method performs better since the square root term $\mathbold{L}$ and $\mathbold{K}$ are much smaller than $\mathbold{E}$ and produces less numerical errors.

For $\mathbold{f}$ and $\mathbold{g}$ , integration of the ODE as given in (9) and (10) is used.

Algorithm LABEL:algo:approx summarizes the computation of $\widetilde{\Pi}$ . An implementation on Matlab is available online [14].

5 Examples

Algorithm LABEL:algo:approx deduced from Theorem 2 and 3 is used to compute the overapproximation $\Pi^{\cap}$ defined in (29) (subset of $\Pi$ defined in (19)) of reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ of the system $\mathscr{S}(\mathcal{P}_{0},t)$ (described in Section 1.1), $t\geq 0$ . Several examples are treated. With these examples, we provide some performance evaluations of our approach.

5.1 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 [32]. For each system, a stabilizing controller is generated for the generalized plant using the h2syn function of Matlab, then the system is reduced using a balanced truncation method to a given state space size. 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)=\left[\begin{smallmatrix}1&\dots&1\end{smallmatrix}\right]^{\top}\exp(-t)$ for $t\in[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 2.

In Figure 7, 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.

5.2 System Verification

We study the stable IQC system $\mathscr{S}(\mathcal{P}_{0},t)$ , defined in (6), at a given time $t$ in $[0,1]$ , for a parabolic set of initial states $\mathcal{P}_{0}=\mathcal{P}(E_{0},f_{0},g_{0})$ , with $E_{0}=\left[\begin{smallmatrix}a+b&a\\ a&a+b\end{smallmatrix}\right],\,f_{0}=\left[\begin{smallmatrix}0\\ 0\end{smallmatrix}\right],\,g_{0}=0.015,\,a=10^{-2}\textrm{ and }b=10^{-6}$ , and for the following parameters

[TABLE]

where $I=\left[\begin{smallmatrix}1&0\\ 0&1\end{smallmatrix}\right]$ , and with a zero input signal $u$ .

The reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ of $\mathscr{S}(\mathcal{P}_{0},t)$ , defined in (7), is computed using (18) and Theorem 3, for $t\in[0,1]$ . Figures 8(a) and 8(b) show the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ set at time $t=0.794$ and its projection $\mathcal{R}(\mathcal{P}_{0},t)|_{x}$ over the LTI state space (i.e. projection over $(x_{1},x_{2})$ states). In Figure 8(b), the constraints boundaries $\partial P(t)$ (for $P\in\Pi$ , $\Pi$ defined in Section 3) are touching the reachable set $\mathcal{R}(\mathcal{P}_{0},t)$ . The non-convexity of $\mathcal{R}(\mathcal{P}_{0},t)$ arises from the non-positive solutions to the DRE (8). Figure 8(c) represents the projection of the reachable tube $t\mapsto\mathcal{R}(\mathcal{P}_{0},t)$ projected over the LTI dimension $(x_{1},x_{2})$ .

6 Conclusion

IQCs are 2-norm constraints (i.e. energetic constraints) between signals. Classical models in the robust control community involve 2-norm constraint or $\infty$ -norm constraint (i.e. hard bounds between signal). $\infty$ -norm bounds over the signals have been treated (e.g. the ellipsoidal method). However, complex systems can be usually described with many relationships. In future works, the computation of the reachable set for systems with multiple 2-norm and $\infty$ -norm relationships will be investigated. Also, such models will be used to define abstractions for nonlinear systems. The dynamic of the system will be linearized and the non-linearity modeled as a bounded (2-norm and/or $\infty$ -norm) disturbance.

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 of 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 exists for this optimization problem.

The differential Riccati equation can be weakly solved using a basis of polynomial solutions. 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 approach with the results presented in this paper.

Locally optimal solution of the optimization problem can be derived using the maximum of Pontryagin principle. Such solution are already available for the ellipsoidal method. The adaptation to the paraboloidal method will be the topic of a future work.

Our integration scheme is not guaranteed and the paraboloids we compute are subject to the error of the differential equation numerical integration. Guaranteed integration scheme exists to overapproximate the reachable set of linear time-invariant systems. The image of an ellipsoidal set through the matrix exponential is overapproximated. This result is then used to overapproximate the reachable set of a linear system. In future works, we will develop a similar approach for the paraboloidal method.

In our implementation, the scaling functions and initial scaling factor 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 most of the computational effort are linear in the number of time-varying paraboloids that needs to be simulated, an efficient choice of the scaling factor can lead to algorithms that demand less computational resources.

In control applications where the system is described by a partial differential equation, a linear approximation of the model can be derived by projecting the state over a finite basis of function. The approximation is then described by an ordinary differential equation that usually have a high number of states (several order of magnitudes). In its current implementation, our method proved to be efficient for systems of less than a hundred states. For higher system dimension, the numerical integration of the differential Riccati equation might cumbersome. In future works, numerical integration of sparse differential Riccati equation could be used to treat such examples.

Bibliography32

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] F. Blanchini, S. Miani, Set-theoretic methods in control, Springer, Boston: Birkhäuser, 2008.
2[2] L. Jaulin, M. Kieffer, O. Didrit, E. Walter, Applied interval analysis: with examples in parameter and state estimation, robust control and robotics, Vol. 1, Springer-Verlag London, London, 2001.
3[3] A. M. Bayen, I. M. Mitchell, M. K. Osihi, C. J. Tomlin, Aircraft autolander safety analysis through optimal control-based reach set computation, Journal of Guidance, Control, and Dynamics 30 (1) (2007) 68–77.
4[4] A. Megretski, A. Rantzer, System analysis via integral quadratic constraints, IEEE Transactions on Automatic Control 42 (6) (1997) 819–830.
5[5] A. Megretski, KYP lemma for non-strict inequalities and the associated minimax theorem, ar Xiv preprint ar Xiv:1008.2552.
6[6] A. Helmersson, An IQC-based stability criterion for systems with slowly varying parameters, in: International Federation of Automatic Control, Vol. 32: 14th World Congress, Elsevier, 1999, pp. 3183–3188.
7[7] A. Megretski, Integral quadratic constraints for systems with rate limiters, Tech. Rep. LIDS-P-2407, Massachusetts Institute of Technology, Laboratory for Information and Decision Systems, Cambridge, Massachusetts (1997).
8[8] D. Peaucelle, L. Baudouin, F. Gouaisbaut, Integral quadratic separators for performance analysis, in: European Control Conference, Budapest, 2009.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Parabolic Set Simulation for Reachability Analysis of Linear Time-Invariant Systems with Integral Quadratic Constraint

Abstract

keywords:

1 Introduction

Related work

Contributions

Outline

1.1 Notation

1.2 System

1.3 Paraboloids

Definition 1** (Paraboloid).**

Definition 2** (Scaled Paraboloid).**

Property 1**.**

Proof 1**.**

Definition 3** (Touching Trajectory).**

1.4 Problem Statement

Problem 1**.**

Problem 2**.**

2 Overapproximation with Paraboloids

Definition 4** (Time-Varying Paraboloid).**

Property 2**.**

Proof 2**.**

Theorem 1** **(Solution to Problem 1)

Proof 3**.**

Property 3**.**

Proof 4**.**

Remark 1** (Representation of paraboloids).**

Example 1**.**

2.1 Domain of definition of PPP

Property 4**.**

3 Exact Reachable Set

3.1 Scaled Paraboloids

Property 5**.**

Proof 5**.**

Property 6**.**

Proof 6**.**

Theorem 2

Proof 7**.**

3.2 Definition domain of time-varying paraboloid

Property 7**.**

Remark 2**.**

Example 2**.**

3.3 Overapproximation with an intersection of time-varying paraboloids

3.4 Overapproximation Relationship

Property 8**.**

Proof 8**.**

Example 3** (Continued from Example 1).**

3.5 Past trajectory for states in the overapproximation’s interior

Property 9**.**

Proof 9**.**

3.6 Past trajectory for states in the overapproximation’s boundary

Property 10**.**

Proof 10**.**

Remark 3**.**

Property 11**.**

Proof 11**.**

Property 12**.**

Proof 12**.**

Assumption 1**.**

Assumption 2**.**

Lemma 1**.**

Proof 13**.**

3.7 Exact Reachable Set

Theorem 3** **(Exact reachability, solution of Problem 2)

Proof 14**.**

4 Implementation

Subset of scaling functions and initial scaling factors:

Paraboloid numerical integration:

Cardinal limitation of Π~\widetilde{\Pi}Π:

DRE numerical integration:

5 Examples

5.1 Examples from COMPleib

5.2 System Verification

6 Conclusion

Definition 1 (Paraboloid).

Definition 2 (Scaled Paraboloid).

Property 1.

Proof 1.

Definition 3 (Touching Trajectory).

Problem 1.

Problem 2.

Definition 4 (Time-Varying Paraboloid).

Property 2.

Proof 2.

Theorem 1 (Solution to Problem 1)

Proof 3.

Property 3.

Proof 4.

Remark 1 (Representation of paraboloids).

Example 1.

2.1 Domain of definition of $P$

Property 4.

Property 5.

Proof 5.

Property 6.

Proof 6.

Proof 7.

Property 7.

Remark 2.

Example 2.

Property 8.

Proof 8.

Example 3 (Continued from Example 1).

Property 9.

Proof 9.

Property 10.

Proof 10.

Remark 3.

Property 11.

Proof 11.

Property 12.

Proof 12.

Assumption 1.

Assumption 2.

Lemma 1.

Proof 13.

Theorem 3 (Exact reachability, solution of Problem 2)

Proof 14.

Cardinal limitation of $\widetilde{\Pi}$ :