Finite Horizon Robustness Analysis of LTV Systems Using Integral Quadratic Constraints

# Finite Horizon Robustness Analysis of LTV Systems Using Integral Quadratic Constraints

Pete Seiler111P. Seiler is with the Department of AEM at the University of Minnesota, MN, USA. seile017@umn.edu    Robert M. Moore222R. M. Moore, C. Meissen, and A. Packard are with the Department of Mechanical Engineering at the University of California, Berkeley, CA 94720, USA.
{max.moore,cmeissen,apackard}@berkeley.edu
Chris Meissen222R. M. Moore, C. Meissen, and A. Packard are with the Department of Mechanical Engineering at the University of California, Berkeley, CA 94720, USA.
{max.moore,cmeissen,apackard}@berkeley.edu
Murat Arcak 333M. Arcak is with the Department of Electrical Engineering and Computer Science at the University of California, Berkeley, CA 94720, USA. arcak@berkeley.edu    Andrew Packard222R. M. Moore, C. Meissen, and A. Packard are with the Department of Mechanical Engineering at the University of California, Berkeley, CA 94720, USA.
{max.moore,cmeissen,apackard}@berkeley.edu
###### Abstract

The goal of this paper is to assess the robustness of an uncertain linear time-varying (LTV) system on a finite time horizon. The uncertain system is modeled as a connection of a known LTV system and a perturbation. The input/output behavior of the perturbation is described by time-domain, integral quadratic constraints (IQCs). Typical notions of robustness, e.g. nominal stability and gain/phase margins, can be insufficient for finite-horizon analysis. Instead, this paper focuses on robust induced gains and bounds on the reachable set of states. Sufficient conditions to compute robust performance bounds are formulated using dissipation inequalities and IQCs. The analysis conditions are provided in two equivalent forms as Riccati differential equations and differential linear matrix inequalities. A computational approach is provided that leverages both forms of the analysis conditions. The approach is demonstrated with two examples.

Finite Horizon Robustness Analysis of LTV Systems Using Integral Quadratic Constraints
 Pete Seiler111P. Seiler is with the Department of AEM at the University of Minnesota, MN, USA. seile017@umn.edu and Robert M. Moore222R. M. Moore, C. Meissen, and A. Packard are with the Department of Mechanical Engineering at the University of California, Berkeley, CA 94720, USA. {max.moore,cmeissen,apackard}@berkeley.edu and Chris Meissen222R. M. Moore, C. Meissen, and A. Packard are with the Department of Mechanical Engineering at the University of California, Berkeley, CA 94720, USA. {max.moore,cmeissen,apackard}@berkeley.edu and Murat Arcak 333M. Arcak is with the Department of Electrical Engineering and Computer Science at the University of California, Berkeley, CA 94720, USA. arcak@berkeley.edu and Andrew Packard222R. M. Moore, C. Meissen, and A. Packard are with the Department of Mechanical Engineering at the University of California, Berkeley, CA 94720, USA. {max.moore,cmeissen,apackard}@berkeley.edu

## 1 Introduction

This paper develops theoretical and computational methods to analyze the robustness of linear time-varying (LTV) systems over finite time horizons. Motivating applications for this work include robotic systems [19, 28] and space launch / re-entry vehicles [16, 32] both of which undergo finite-time trajectories. Typical notions of robustness, e.g. nominal stability and gain/phase margins, can be insufficient for such systems. For example, one approach is to evaluate the stability and performance of the LTV system at “frozen” time instances along the trajectory. However, there are LTV systems for which is stable for each frozen time but with trajectories that grow unbounded [12].

This paper moves beyond the frozen analysis technique and instead evaluates time-domain metrics over the finite horizon. The analysis is performed on an uncertain LTV system modeled, as shown in Figure 1, by an interconnection of a known, nominal LTV system and a perturbation . The perturbation is used to model difficult to analyze elements including nonlinearities and dynamic or parametric uncertainty. The input-output properties of are characterized by integral quadratic constraints (IQCs) [17]. An extensive list of IQCs for various classes of perturbations is given in [17, 33]. The main result in [17] is an (infinite-horizon), input-output stability theorem using frequency domain IQCs. The proof relies on an operator theoretic, homotopy method.

The main contribution of this paper is an algorithm to compute robustness metrics on a finite horizon. First, nominal finite-horizon LTV performance is reviewed (Section 2) focusing on induced and -to-Euclidean gains. The -to-Euclidean gain is useful for computing bounds on the reachable set of states. Next, sufficient conditions are given (Section 3) to bound the induced and -to-Euclidean gains for uncertain LTV systems. The analysis is formulated with dissipation inequalities and IQCs. This yields performance conditions in the form of infinite-dimensional differential linear matrix inequalities (DLMIs). These conditions can be equivalently re-written with a Riccati Differential Equation (RDE). This equivalence is based on a variation of the strict bounded real lemma (Theorem 1 in Section 2) which generalizes existing results in [29, 24, 10, 4]. The algorithm to compute finite horizon robustness metrics (Section 4) leverages both forms of these conditions. The proposed approach is demonstrated with two examples (Section 5).

This paper adds to existing results to assess robustness of time-varying systems. The most closely related work is [11] which also considers finite horizon robustness with IQCs. However, the work in [11] does not consider disturbances and the theoretical approach / resulting numerical algorithm are different than given here. The work in [23] and [9] is also relevant. Uncertain linear parameter varying systems are considered in [23]. The analysis is formulated using dissipation inequalities and IQCs. A similar approach is used in [9] to provide robustness conditions for discrete-time LTV systems. The theoretical and computational approaches provided here differ from these previous works in order to handle continuous-time LTV systems on finite horizons. Other related work includes robustness analysis for finite horizon batch processes [14], nonlinear systems [30], and periodic LTV systems via time-domain lifting [6, 15, 13] or harmonic balance / frequency domain lifting [34, 7, 27].

Notation: Let and denote the sets of -by- real matrices and -by- real, symmetric matrices. The finite-horizon norm of a signal is defined as . If is finite then . denotes the set of rational functions with real coefficients that have no poles on the imaginary axis. is the subset of functions in that are analytic in the closed right-half of the complex plane. Finally, let and be given with and . Define an ellipsoid in terms of as

## 2 Nominal Performance

### 2.1 Finite Horizon LTV Systems

Consider an LTV system defined on :

 ˙x(t) =A(t)x(t)+B(t)d(t) (1) e(t) =C(t)x(t)+D(t)d(t) (2)

is the state, is the input, and is the output. The state matrices , , , and are piecewise-continuous (bounded) functions of time. It is assumed throughout that . Thus implies and are in for any [3].

Many different performance metrics can be defined for this (nominal) finite-horizon LTV system. This paper mainly focuses on two specific metrics. First, the finite-horizon induced -gain of is

 ∥G∥2,[0,T]:=sup{∥e∥2,[0,T]∥d∥2,[0,T] ∣∣∣ x(0)=0,0≠d∈L2,[0,T]}.

As noted above, implies . Thus the gain is finite for any fixed horizon .

Next, assume . Then the finite-horizon -to-Euclidean gain of is

 ∥G∥E,[0,T]:=sup{∥e(T)∥2∥d∥2,[0,T] ∣∣∣ x(0)=0,0≠d∈L2,[0,T]}.

The -to-Euclidean gain depends on the system output only at the final time . The assumption that ensures this gain is well-defined.

The -to-Euclidean gain can be used to bound the set of states reachable by disturbances of a given norm. This reachable set is formally defined as follows:

 Rβ:={x(T) ∣∣ x(0)=0,∥d∥2,[0,T]≤β}.

If and then . In this special case, if then . This implies the reachable set is contained in a sphere of radius . More general ellipsoidal bounds on the reachable set can be obtained by proper selection of the output matrices. For example, select and for some given with . With these choices implies an ellipsoidal bound on the reachable set: . The size of the reachable set scales with the norm of the disturbance input. The state at intermediate times can similarly be bounded using the -to-Euclidean gain .

### 2.2 Generic Quadratic Cost

The two nominal performance metrics introduced above are generalized in Section 3 to assess robustness of uncertain systems. A generic quadratic cost function is defined next in order to unify these various nominal and robust performance metrics. Specifically, let , , , and be given. are assumed to be piecewise continuous (bounded) functions. A quadratic cost function is defined by as follows:

 J(d) :=x(T)TFx(T)+∫T0[x(t)d(t)]T[Q(t)S(t)S(t)TR(t)][x(t)d(t)]dt subject to: Eq.~{}??? with x(0)=0 (3)

The finite-horizon induced gain of can be related to the quadratic cost by proper choice of . In particular, let be given and select , , , and . This yields the following cost function

 J(d)=∥e∥22,[0,T]−γ2∥d∥22,[0,T] (4)

Thus if and only if .

The finite-horizon -to-Euclidean gain of can also be related to the quadratic cost but with different choices for . Let be given and select , , , and . This yields the following cost function

 J(d)=∥e(T)∥22−γ2∥d∥22,[0,T] (5)

Thus if and only if .

### 2.3 Strict Bounded Real Lemma

The next theorem states an equivalence between a bound on the quadratic cost and the existence of a solution to a Riccati Differential Equation (RDE) or Riccati Differential Inequality (RDI). The theorem is expressed in terms of strict inequalities and generalizes existing results for the induced gain of LTV systems [29, 24, 10, 4].

###### Theorem 1.

Let be given with for all . The following statements are equivalent:

1. such that .

2. There exists a differentiable function such that and

 ˙Y+ATY+YA+Q−(YB+S)R−1(YB+S)T=0

This is a Riccati Differential Equation (RDE).

3. There exists and a differentiable function such that and

 ˙P +ATP+PA+Q −(PB+S)R−1(PB+S)T⪯−ϵI

This is a strict Riccati Differential Inequality (RDI).

###### Proof.

The proof of () is given below as it highlights the dissipation inequality framework. The remainder of the proof is in Appendix A for completeness.

By the Schur complement lemma [2], the RDI and imply such that

 [˙P+ATP+PAPBBTP0]+[QSSTR]⪯−~ϵI (6)

Next define a a quadratic storage function . Let be a solution of the LTV system (Equation 1) starting from and forced by some input . Multiply Equation 6 on the left and right by and its transpose to obtain the following dissipation inequality:

 ˙V+[xd]T[QSSTR][xd]≤−~ϵ[xd]T[xd] (7)

Integrate the dissipation inequality from to :

 V(x(T),T)−V(x(0),0)

Apply the boundary condition to obtain:

 J(d)≤V(x(0),0)−~ϵ∥d∥22,[0,T] (8)

This bound is valid for any . Hence, applying yields .

This theorem assumes zero initial conditions . This implies that the initial stored energy is zero and hence is dropped from Equation 8 in the proof. Non-zero initial conditions can be incorporated by retaining this initial stored energy in the performance bound.

Nominal performance is most easily assessed using the RDE. The performance is achieved if the associated RDE exists on when integrated backward from . The assumption ensures is invertible and hence the RDE is well-defined for all . Thus the solution of the RDE exists on unless it grows unbounded. As a concrete example, it was noted in Section 2.2 that for specific choices of . The matrix , and hence the RDE, depends on the choice of . For a fixed , the performance is achieved if the associated RDE exists on when integrated backward from . The smallest bound on the induced gain can be found via bisection. The RDI will be used later to assess the robustness of uncertain LTV systems.

## 3 Robust Performance

### 3.1 Uncertain LTV Systems

An uncertain, finite-horizon LTV system is given by the interconnection of a nominal LTV system and a perturbation as shown in Figure 1. The LTV system is described by the following state-space model:

 ˙xG(t)=AG(t)xG(t)+BG1(t)w(t)+BG2(t)d(t)v(t)=CG1(t)xG(t)+DG11(t)w(t)+DG12(t)d(t)e(t)=CG2(t)xG(t)+DG21(t)w(t)+DG22(t)d(t) (9)

where is the state. The inputs are and while and are outputs. The state matrices are piecewise continuous (bounded) functions of time with appropriate dimensions, e.g. . The perturbation is a bounded, causal operator. Well-posedness of the interconnection is defined as follows.

###### Definition 2.

is well-posed if for all and there exists unique solutions , , , and satisfying Equation (9) and with a causal dependence on .

The perturbation can have block-structure as is standard in robust control modeling [37]. It can include blocks that are hard nonlinearities (e.g. saturations) and infinite dimensional operators (e.g. time delays) in addition to true system uncertainties. The term “uncertainty” is used for simplicity when referring to .

### 3.2 Integral Quadratic Constraints (IQCs)

IQCs [17] are used to describe the input/output behavior of . They can be formulated in either the frequency or time domain. The time domain formulation is more useful for analysis of uncertain time-varying systems. This formulation is based on the graphical interpretation in Figure 2. The inputs and outputs of are filtered through an LTI system with zero initial condition . The dynamics of are given as follows:

 ˙xψ(t)=Aψxψ(t)+Bψ1v(t)+Bψ2w(t)z(t)=Cψxψ(t)+Dψ1v(t)+Dψ2w(t) (10)

where is the state. A time domain IQC is an inequality enforced on the output over a finite horizon. The formal definition is given next.

###### Definition 3.

Let and be given with piecewise continuous. A bounded, causal operator satisfies the time domain IQC defined by if the following inequality holds for all and :

 ∫T0z(t)TM(t)z(t)dt≥0 (11)

where is the output of driven by inputs with zero initial conditions .

The notation is used when satisfies the corresponding IQC. Time domain IQCs, as defined above, are specified as finite-horizon constraints on the outputs of . These are often referred to as hard IQCs [17]. The definition given here only requires the IQC to hold over the analysis horizon . This is in contrast to hard IQCs used for infinite horizon analysis which require the constraint to hold over all finite time horizons. Two examples of time domain IQCs are provided below.

###### Example 4.

Consider an LTI uncertainty with . Let be given with for all . Then the following frequency domain IQC holds and

 (12)

where and are Fourier transforms of and . This IQC corresponds to the use of -scales used in structured singular value analysis [25, 5, 20, 37]. A factorized representation for yields a valid time domain IQC. Specifically, let where

 (13)

It is shown in [1] that is a valid time domain IQC for over any finite horizon .

###### Example 5.

Time domain IQCs are often specified with as an LTI system and as a constant matrix. Definition 3 above allows to be time-varying. This generalization is useful for time-varying real parameters. Let where and for all . Define and where is piecewise continuous and satisfies . Then satisfies the time domain IQC defined by . Time-varying IQCs can be defined for other uncertainties, e.g. see related work in [22].

An extensive library of IQCs is provided in [17] for various types of perturbations. Most IQCs are specified in the frequency domain using a multiplier . Under some mild assumptions, a valid time-domain IQC can be constructed from via a -spectral factorization [26]. This allows the library of known (frequency domain) IQCs to be used for time-domain, finite-horizon analysis. More general IQC parameterizations are not necessarily “hard” but can be handled with the method in [8].

### 3.3 Robust Induced L2 Gain

The robustness of is analyzed using the interconnection shown in Figure 3. The extended system of (Equation 9) and the IQC filter (Equation 10) is governed by the following state space model:

 ˙x(t)=A(t)x(t)+B(t)[w(t)d(t)]z(t)=C1(t)x(t)+D1(t)[w(t)d(t)]e(t)=C2(t)x(t)+D2(t)[w(t)d(t)] (14)

The extended state vector is where . The state-space matrices are given by (dropping the dependence on time ):

 A:=[AG0Bψ1CG1Aψ],B:=[BG1BG2Bψ1DG11+Bψ2Bψ1DG12] D1:=[Dψ1DG11+Dψ2Dψ1DG12] D2:=[DG21DG22]

The actual system to be analyzed is with input and initial condition . The analysis is instead performed with the extended LTV system (Equation 14) and the constraint . The constrained extended system has inputs and initial condition . The IQC implicitly constrains the input . The IQC covers such that the constrained extended system without includes all behaviors of the original system .

The following differential matrix inequality is used to assess the robust performance of  111The notation in (15) corresponds to an omitted factor required to make the corresponding term symmetric.:

 [˙P+ATP+PAPBBTP0]+[QSSTR]+(⋅)TM[C1D1]⪯−ϵI (15)

This inequality depends on the extended system, IQC, and quadratic cost . It is compactly denoted as . This notation emphasizes that the constraint is a differential linear matrix inequality (DLMI) in for fixed and . The dependence on and is not explicitly denoted but will be clear from context.

The next theorem formulates a sufficient condition to bound the (robust) induced gain of . The proof uses IQCs and a standard dissipation argument [31, 35, 36, 12]. For induced gains the quadratic cost matrices are chosen as and

 Q(t):=C2(t)TC2(t),S(t):=C2(t)TD2(t)R(t):=D2(t)TD2(t)−γ2[0nw00Ind] (16)
###### Theorem 6.

Let be an LTV system defined by (9) and be a bounded, causal operator. Assume is well-posed and . If there exists , and a differentiable function such that and

 DLMIRob(P,M,γ2,t)⪯−ϵI∀t∈[0,T] (17)

then .

###### Proof.

Let and be given. By well-posedness, has a unique solution . As noted above, the extended system and IQC “cover” this system. In particular, forcing with from yields . Define . Then are a solution of the extended system (14) with inputs and initial condition . Moreover, satisfies the the IQC defined by .

Define a storage function by . Left and right multiply the DLMI (15) by and its transpose to show that satisfies the following dissipation inequality for all :

 ˙V+[x[wd]]T[QSSTR][x[wd]]+zTMz≤−ϵdTd (18)

Use the choices for in (16) to rewrite the second term as . Integrate over to obtain:

 x(T)TP(T)x(T)−xTG,0P11(0)xG,0+∫T0zT(t)M(t)z(t)dt −(γ2−ϵ)∥d∥22,[0,T]+∥e∥22,[0,T]≤0.

Apply and to conclude:

 ∥e∥22,[0,T]≤xTG,0P11(0)xG,0+(γ2−ϵ)∥d∥22,[0,T] (19)

Finally, if then .

The effect of non-zero initial conditions is captured in Equation 19. If , this simplifies to implying . Furthermore, on this implies for all . The IQC is valid only for and hence need not hold in general.

### 3.4 Robust L2-to-Euclidean Gain

A similar theorem provides a bound on the -to-Euclidean gain of . This requires the additional assumptions that and so that . Hence and the gain from to is well-defined. To assess the robust -to-Euclidean gain define as:

 (20)

With these choices for the next theorem is a minor adaptation of Theorem 6 and the proof is omitted.

###### Theorem 7.

Let be an LTV system defined by (9) with and . Let be a bounded, causal operator. Assume is well-posed and . If there exists , , and a differentiable function and such that and

 DLMIRob(P,M,γ2,t)⪯−ϵI∀t∈[0,T] (21)

then .

The condition in Theorem 7 robustly bounds the states reachable by disturbances for any uncertainty . Note so that only depends on . The IQC filter is used only for analysis and is neglected in the bound.

Robust reachable sets with non-zero initial conditions can be computed with minor modifications. For example, assume the initial condition of lies within the ellipsoid for some . The IQC still requires . Next, enforce for some (in addition to the conditions in Theorem 7). It follows from the dissipation inequality proof that the terminal state of is bounded by . Additional variations on robust reachable sets with non-zero initial conditions can be found in Chapter 2 of [18].

### 3.5 RDE Condition for Robust Performance

Theorems 6 and 7 provide a DLMI (15) to bound the induced and -to-Euclidean gain of . More general robust performance conditions can be formulated by proper choice of . The numerical algorithm proposed in Section 4 relies on an equivalence between the DLMI (15) and a related RDE condition. This equivalence is demonstrated with an extended quadratic cost function that combines the performance specification and the IQC . Specifically, define with the extended dynamics in (14): . The cost matrices are chosen as:

 [QSSTR]:=(⋅)TM[C1D1]+[QSSTR]F:=F (22)

The quadratic cost associated with these choices is:

 J([wd]) :=x(T)TFx(T)+∫T0zT(t)M(t)z(t)dt +∫T0⎡⎣x(t)[w(t)d(t)]⎤⎦T[Q(t)S(t)S(t)TR(t)]⎡⎣x(t)[w(t)d(t)]⎤⎦dt

The next corollary states the equivalence between the DLMI and RDE conditions. The DLMI can be rewritten as an RDI by the Schur complement lemma [2]. Hence the corollary follows directly from Theorem 1.

###### Corollary 8.

Let be given by (22). The following are equivalent for any and :

1. There exists a differentiable function such that and .

2. for all . In addition, there exists a differentiable function such that and

 ˙Y+ATY+YA+Q−(YB+S)R−1(YB+S)T=0

## 4 Computational Approach

This section describes computational details and presents an algorithm that combines complementary aspects of the DLMI and RDE robust performance conditions.

### 4.1 IQC Parameterization

There is typically an infinite set of valid IQCs for a given uncertainty . Numerical implementations using IQCs often involve a fixed choice for and optimization subject to convex constraints on [17, 33, 21]. The algorithms given in the following sections will use such parameterizations. Two examples are given below.

###### Example 9.

Consider an LTI uncertainty with . By Example 4, satisfies any IQC with , , and . A typical choice for is [33]:

 Ψv11:=[1,1(s+p),…1(s+p)v]T with p>0 (23)

The robustness analysis is performed by selecting to obtain (fixed) and optimizing over . The results depend on the choice of . Larger values of represent a richer class of IQCs and hence yield less conservative results but with increasing computational cost. Further details on this parameterization are given in [33].

###### Example 10.

The analysis can incorporate conic combinations of multiple IQCs. Let and define valid IQCs for . Hence where is the output driven by and . For any the two constraints can be combined to yield:

 ∫T0λ1zT1M1z1+λ2zT2M2z2dt≥0 (24)

Thus a valid time-domain IQC for is given by

 (25)

The analysis is performed by selecting and optimizing over .

### 4.2 Analysis with the DLMI Condition

Assume the IQC is with fixed and constrained to lie within a feasible set described by LMIs. The DLMI (15) has the same form for induced and -to-Euclidean gains but with different choices of . In both cases the DLMI is linear in for fixed . The dependence on enters via . The best (smallest) bound on the robust gain can be computed from a convex semidefinite program (SDP):

 minγ2 subject to: M∈M,P(T)⪰F DLMIRob(P,M,γ2,t)⪯−ϵI∀t∈[0,T]

There are two main issues with solving this SDP. First, the DLMI corresponds to an infinite number of constraints since it must hold for all . This can be approximated by enforcing the DLMI on a finite time grid .

Second, the optimization requires a search over the space of functions . This issue is addressed by restricting to be a linear combination of differentiable basis functions. Specifically, let and be given scalar and matrix differentiable basis functions. The storage function and its derivative are given by:

 P(t) =Ns∑j=1hj(t)Xj+Nm∑k=1Hk(t)xk (26) ˙P(t) =Ns∑j=1˙fj(t)Xj+Nm∑k=1˙Fk(t)xk (27)

Here and are optimization variables. Many choices are possible for the basis functions. Initial work in [18] used scalar basis functions generated with a cubic spline and no matrix basis functions. The spline is constructed by selecting an interpolation time grid where . Note, the spline grid is distinct from the DLMI grid . The spline consists of cubic functions defined on the intervals . It interpolates the decision variables , i.e. . The cubic functions satisfy boundary conditions to ensure continuity of the spline and its first/second derivatives at the interval endpoints. The corresponding spline basis functions are not easy to express in explicit form but they can be evaluated numerically at any . Additional details are given in [18]. The algorithm proposed below also uses a matrix basis function generated by the RDE condition.

The approximations for the DLMI and lead to a finite dimensional SDP in variables , , , and