#  Polyhedral Estimation of $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$
Incremental Gains of Nonlinear Systems 111This work was supported by the
Engineering and Physical Sciences Research Council (EPSRC) of the United
Kingdom - Project Reference 1950140.

Dimitris Kousoulidis<EMAIL_ADDRESS>Fulvio Forni<EMAIL_ADDRESS>University of Cambridge, Department of Engineering, Trumpington Street, CB2
1PZ Cambridge, United Kingdom

###### Abstract

We provide novel dissipativity conditions for bounding the incremental
$\mathcal{L}_{1}$ gain of systems. Moreover, we adapt existing results on the
$\mathcal{L}_{\infty}$ gain to the incremental setting and relate the
incremental $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gain bounds through
system adjoints. Building on work on optimization based approaches to
constructing polyhedral Lyapunov functions, we make use of these conditions to
obtain a Linear Programming based algorithm that can provide increasingly
sharp bounds on the gains as a function of a given candidate polyhedral
storage function or polyhedral set. The algorithm is also extended to allow
for the design of linear feedback controllers for performance, as measured by
the bounds on the incremental gains. We apply the algorithm to a couple of
numerical examples to illustrate the power, as well as some limitations, of
this approach.

###### keywords:

Incremental gains , Polyhedra , Linear programming , Nonlinear control

##  1 Introduction

When we treat open dynamical systems as operators $\mathcal{G}$ that map input
signals $w$ onto output signals $z=\mathcal{G}w$, the induced operator norms
of $\mathcal{G}$ quantitatively measure the amplification factor of the input-
output relationship, which we call the system gain. Gains provide a useful
tool for system specifications, allowing us to quantitatively characterize
performance and to pose control design problems to optimize them in closed
loop [28]. Gains are also used for studying the stability of interconnections
through the small gain theorem, which states that the feedback interconnection
of two systems $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ will be stable as long
as $\|\mathcal{G}_{1}\|\|\mathcal{G}_{2}\|<1$ [12, Chapter 3]. For these
purposes, gains based on various norms have been extensively studied and used
in systems and control theory. The selection of the norm needs to balance the
needs of capturing performance versus our ability to compute with it [8, 28].
In this paper we will be using incremental gains. Incremental gains are a
classic concept in control literature [12]. They bound the incremental
amplification between any pair of input signals. Incremental gains maintain
the interpretations and properties of gains discussed above, including a form
of the small gain theorem [12, Chapter 3]. For linear systems, gains and
incremental gains are equivalent.

We will be computing bounds on the incremental $\mathcal{L}_{1}$ and the
$\mathcal{L}_{\infty}$ gains:

$\|\mathcal{G}\|_{1}=\sup_{w_{1}-w_{2}\in\mathcal{L}_{1}\setminus\\{0\\}}\frac{\|z_{1}-z_{2}\|_{1}}{\|w_{1}-w_{2}\|_{1}},$
(1)

and

$\|\mathcal{G}\|_{\infty}=\sup_{w_{1}-w_{2}\in\mathcal{L}_{\infty}\setminus\\{0\\}}\frac{\|z_{1}-z_{2}\|_{\infty}}{\|w_{1}-w_{2}\|_{\infty}},$
(2)

respectively, where

$\displaystyle\|w\|_{1}$
$\displaystyle=\int_{0}^{\infty}|w|_{1}dt=\int_{0}^{\infty}\sum_{i}|w_{i}|dt,$
$\displaystyle\|w\|_{\infty}$
$\displaystyle=\sup_{t}|w|_{\infty}=\sup_{i,\,t}|w_{i}|.$

$w\in\mathcal{L}_{1}$ if $\|w\|_{1}<\infty$, and $w\in\mathcal{L}_{\infty}$ if
$\|w\|_{\infty}<\infty$.

The $\mathcal{L}_{1}$ norm assigns lower magnitude to short, high amplitude
signals over long, low amplitude signals; while the $\mathcal{L}_{\infty}$
norm assigns lower magnitude to long, low amplitude signals over short, high
amplitude signals. The $\mathcal{L}_{1}$ gain is most suitable in quantifying
specifications in which the signals are each directly related to the
consumption of a resource, for example fuel use in a rocket or transaction
fees in a trading system, while the $\mathcal{L}_{\infty}$ gain is most
suitable in quantifying specifications about the worse-case amplitudes of
signals. Moreover, for linear systems, upper bounds on the $\mathcal{L}_{1}$
and $\mathcal{L}_{\infty}$ gains, $\gamma_{1}$ and $\gamma_{\infty}$, provide
upper bounds for any $\mathcal{L}_{p}$ gain with $1\leq p\leq\infty$,
$\gamma_{p}$ [12]

$\gamma_{p}\leq\gamma_{1}^{1/p}\gamma_{\infty}^{1/q},\ \text{where
}1/p+1/q=1,$ (3)

offering another reason to study this combination of gains.

$\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains have been mostly studied in
the context of positive systems. This is because of the ease with which they
can be expressed and their efficient calculation [10, 23]. This has also led
to conditions for robust stability and stabilization via integral linear
constraints, a linear analogue to integral quadratic constraints [10].
Moreover, weighted $\mathcal{L}_{1}$ gains can be used to obtain a complete
characterization of stability for positive interconnections of positive
systems [14].

Beyond positive systems, $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains of
linear systems are computed through numerical integration of the impulse
response [1, 22, 26]. Current results for $\mathcal{L}_{1}$ and
$\mathcal{L}_{\infty}$ bounds on gains for nonlinear systems include [20,
Theorem 5.1], where a uniform bound for all $\mathcal{L}_{p}$ gains in terms
of a provided Lyapunov function is derived. The work of [4] and [15] presents
algorithms for computing the $\mathcal{L}_{\infty}$ gain of linear
time/parameter varying (LTV/LPV) systems via invariant polyhedral sets.

The main contributions of this paper are novel conditions for computing
incremental $\mathcal{L}_{1}$ gain bounds based on dissipation inequalities
with respect to polyhedral storage functions. Conditions for incremental
$\mathcal{L}_{\infty}$ gain bounds using polyhedral set invariance are also
presented. The $\mathcal{L}_{\infty}$ conditions are closely related to prior
work [4, 15], but tweaked for the incremental gain setting. Both sets of
conditions are framed as optimization problems and solutions to either one of
the gains can be mapped to solutions of the other for an adjoint system,
allowing for any algorithm designed to compute incremental $\mathcal{L}_{1}$
gains to be used for computing incremental $\mathcal{L}_{\infty}$ gain bounds
and vice versa. This generalizes the known relationship between
$\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains for linear systems [9].
When the polyhedral function/set is fixed the conditions are reduced to linear
programming problems. This allows us to adapt our previous work on finding
polyhedral Lyapunov functions of fixed complexity [21] to this problem for
both analysis and synthesis. We provide definitions and the conditions at an
abstract level in Section 2, specialize them to polyhedral functions and
provide our novel conditions in Section 3, provide an overview of our adapted
analysis and synthesis algorithms in Section 4, and apply the algorithms on
some numerical examples in Section 5.

Notation: We use calligraphic letters to denote sets and operators, and
capital letters to denote matrices. We use $\bm{1}$ to denote vectors of all
ones and $[k]$ to denote set $\\{1,\dots,k\\}$. Operations applied to a set
are meant to denote operations applied to each element of the set; for
example, we use $\mathcal{A}^{T}$ to denote $\\{A^{T}:A\in\mathcal{A}\\}$.
$A^{i}$ denotes the $i^{th}$ row of matrix $A$, while $A_{j}$ denotes the
$j^{th}$ column. $[A,B]$ denotes the matrix obtained by stacking matrices $A$
and $B$ _horizontally_ , while $[A;B]$ denotes the matrix obtained by stacking
matrices $A$ and $B$ _vertically_. The vectorization operator
$\operatorname{vec}(\cdot)$ converts a $n\times m$ matrix into a column vector
with $nm$ entries by stacking up its columns. We use $\|x\|$ to denote the
norm of a _signal_ $x$ and $|x|$ to denote the norm of a _vector_ $x$. The
sign function $\operatorname{sgn}(x)$ is set to $-1\text{ if }x<0,0\text{ if
}x=0,\text{and }1\text{ if }x>0$. For any given matrix $A$, $A\geq 0$ means
that all the elements of $A$ are non-negative.

##  2 Gains and Incremental Gains

###  2.1 Gains for Linear Systems

For linear systems, the $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains can
be characterized in terms of the matrix of impulse responses $H$. We denote
the entry in the $i^{th}$ row and $j^{th}$ column of $H$ as $h_{ij}$. This
corresponds to the impulse response of the $i^{th}$ output to the $j^{th}$
input.

The $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains of a system with $m$
outputs and $n$ inputs are then [10]

$\|\mathcal{G}\|_{1}=\max_{1\leq j\leq
n}\left(\sum_{i=1}^{m}\left(\int_{0}^{\infty}|h_{ij}(t)|dt\right)\right)$ (4)

and

$\|\mathcal{G}\|_{\infty}=\max_{1\leq i\leq
m}\left(\sum_{j=1}^{n}\left(\int_{0}^{\infty}|h_{ij}(t)|dt\right)\right).$ (5)

Note that, other than for the computation of the $m\times n$ integrals
$\int_{0}^{\infty}|h_{ij}(t)|dt$, the remaining computations for (4) and (5)
amount to summations and maximizations over a finite number of scalars. As
such evaluating $\int_{0}^{\infty}|h_{ij}(t)|dt$ numerically is the main
computational challenge for evaluating the two gains and various methods have
been proposed for computing it [1, 22, 26].

It follows from (4) and (5) that the $\mathcal{L}_{1}$ gain of a system with
impulse responses $H$ is equivalent to the $\mathcal{L}_{\infty}$ gain of a
system with impulse responses $H^{T}$. This also means that the
$\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains of single-input single-
output or self-adjoint linear systems are equivalent.

###### Remark 1

For the special case of positive linear systems, all $h_{ij}(t)\geq 0$ for all
$t\geq 0$ and so
$\int_{0}^{\infty}|h_{ij}(t)|dt=\int_{0}^{\infty}h_{ij}(t)dt$, which is also
equivalent to the frequency response of $\mathcal{G}$ evaluated at $0$. This
enables efficient calculation and design for $\mathcal{L}_{1}$ and
$\mathcal{L}_{\infty}$ gains [10, 23].

###  2.2 Incremental Gains for Nonlinear Systems

The approach in Section 2.1 relies on concepts from linear time invariant
systems theory that are hard to generalize to other classes of systems. A more
general approach to computing incremental gains (or, at least, upper bounds on
incremental gains) is through dissipativity theory [27], which we apply to
more general, minimal realization state-space systems of the form

$\dot{x}=f(x)+Bw\qquad z=Cx,$ (6)

with $n$ states, $n_{w}$ disturbance inputs, and $n_{z}$ outputs. To obtain a
well defined operator we also set $x(0)=0$.

Under this framework we replace conditions over all signals, such as those in
(1) and (2), with point-wise in time conditions over all time, inputs, and
state variables – called _dissipation inequalities_. This framework then
requires finding a positive function of the state, the _storage function_ ,
which plays a central role in the dissipation inequalities. For simplicity, in
this paper we only consider linear input/output system matrices $B$ and $C$.
However, the approach can be in principle be further generalized to nonlinear
dependencies as well.

For incremental $\mathcal{L}_{1}$ gains the suitable dissipation inequality
for all $w_{1},\,w_{2}$ and $x_{1},\,x_{2}$ is

$\dot{\psi}\leq-|z_{1}-z_{2}|_{1}+\gamma_{1}|w_{1}-w_{2}|_{1},$ (7)

where $\psi(x_{1},x_{2})$ is a storage function. To make it simpler to verify
(7), we limit ourselves to considering storage functions of the form
$\psi(x_{1}-x_{2})$ that satisfy $\psi(x_{1}-x_{2})\geq 0$ for all
$x_{1},x_{2}\in\mathbb{R}^{n}$, and $\psi(0)=0$. Integrating (7) in time from
$0$ to $T$ gives

$\|z_{1}-z_{2}\|_{1}\leq-\psi(x_{1}(T)-x_{2}(T))+\gamma_{1}\|w_{1}-w_{2}\|_{1}$

and so, since $\psi(x_{1}(T)-x_{2}(T))\geq 0$,

$\frac{\|z_{1}-z_{2}\|_{1}}{\|w_{1}-w_{2}\|_{1}}\leq\gamma_{1},$

giving us $\gamma_{1}$ as an upper bound for the incremental $\mathcal{L}_{1}$
gain of the system.

For $\mathcal{L}_{\infty}$ there is no need to integrate the signals involved
because they directly correspond to point-wise in time conditions which are
conveniently expressed in terms of set invariance [4, Section 6.4]. What we
need is a limit on the maximum incremental amplification of the output
relative to the input. We capture this by confirming that inputs with bounded
incremental norms lead to states with bounded difference, i.e. confined to a
bounded set $\mathcal{X}$, for all $t\geq 0$,

$\displaystyle\|w_{1}-w_{2}\|_{\infty}\leq 1,\,x_{1}(0)$
$\displaystyle=0,\,x_{2}(0)=0$
$\displaystyle\implies(x_{1}(t)-x_{2}(t))\in\mathcal{X}$ (8a)
$\displaystyle(x_{1}(t)-x_{2}(t))\in\mathcal{X}$
$\displaystyle\implies\|z_{1}-z_{2}\|_{\infty}\leq\gamma_{\infty},$ (8b)

giving $\gamma_{\infty}$ as an upper bound for the incremental
$\mathcal{L}_{\infty}$ gain of the system. Here the set $\mathcal{X}$ plays a
role similar to that of a storage function.

We next specialize (7) and (8) to the family of polyhedral functions and sets
respectively.

##  3 Polyhedral Approximation of Gains

###  3.1 Polyhedral Preliminaries

To compute gains we use C-sets and gauge functions. C-sets are compact convex
sets that contain $0$ in their interior, and gauge functions are subadditive
positively homogeneous functions (of degree $1$). We can extract a gauge
function from a C-set through its Minkowski functional and a C-set from a
gauge function through its unit ball, so each C-set induces a gauge function
and vice versa [4].

Polyhedra have two distinct representations, which we call V-representation
and H-representation, V-rep and H-rep for short. A V-rep polyhedral set is
represented by its extremal vertices, while a H-rep polyhedral set is
represented by intersections of affine half-spaces. We adopt the convention of
using $\mathcal{V}$ to denote V-rep sets constructed from numerical matrix
$V$. Likewise, we use $\mathcal{H}$ to denote H-rep sets constructed from
numerical matrix $H$.

V-rep

given $n\times m$ matrix $V$,

$\mathcal{V}=\\{x\in\mathbb{R}^{n}:x=Vp,\,p\geq 0,\,\bm{1}^{T}p=1\\}$ (9)

with interior

$\operatorname{int}\mathcal{V}=\\{x:x=Vp,\,p>0,\,\bm{1}^{T}p=1\\}$ (10)

and corresponding gauge function

$\displaystyle\psi_{\mathcal{V}}(x)=$
$\displaystyle\min_{p\in\mathbb{R}^{m}}\\{\bm{1}^{T}p:x=Vp,\,p\geq 0\\}$ (11a)
$\displaystyle=$
$\displaystyle\max_{h\in\mathbb{R}^{n}}\\{h^{T}x:h^{T}V\leq\bm{1}^{T}\\}\,$
(11b)

where (11b) is the dual problem to (11a) [7].

H-rep

given $m\times n$ matrix $H$,

$\mathcal{H}=\\{x\in\mathbb{R}^{n}:Hx\leq\bm{1}\\}$ (12)

with interior

$\operatorname{int}\mathcal{H}=\\{x:Hx<\bm{1}\\}.$ (13)

$\mathcal{V}$ is the convex hull of the columns of $V$, which are called the
vertices of $\mathcal{V}$. On the other hand, $\mathcal{H}$ is the
intersection of half-spaces based on the rows of $H$. We use gauge functions
of V-rep sets for computing $\mathcal{L}_{1}$ gains and H-rep sets for
computing $\mathcal{L}_{\infty}$ gains.

We take a set-theoretic approach to computing $\mathcal{L}_{\infty}$ gains so
we omit the functional perspective for brevity. On the other hand, because we
compute the $\mathcal{L}_{1}$ gain through dissipation inequalities, we need
to characterize the rate of change of V-rep gauge functions, which will act as
the candidate storage functions in (7). We do this through the use of
subdifferentials [24]:

$\partial\psi(x)=\\{h:\forall v\in\mathbb{R}^{n},\,\psi(v)-\psi(x)\geq
h^{T}(v-x)\\}.$ (14)

To provide some intuition about subdifferentials, the subdifferential of an
everywhere differentiable convex function consists at every point of a unique
element that is equivalent to its gradient. Unfortunately,
$\psi_{\mathcal{V}}(x)$ is not everywhere differentiable but we can still
obtain a convenient numerical representation of
$\partial\psi_{\mathcal{V}}(x)$.

######  Proposition 1 (Subdifferential of $\psi_{\mathcal{V}}$)

$\partial\psi_{\mathcal{V}}(x)=\\{h:h^{T}x=\psi_{\mathcal{V}}(x),\,h^{T}V\leq\bm{1}^{T}\\}.$
(15)

###### Proof 1

We first rewrite $\psi_{\mathcal{V}}$ in the following equivalent form

$\psi_{\mathcal{V}}(x)=\max_{h\in\mathcal{Z}}h^{T}x$

where $\mathcal{Z}=\\{h:V^{T}h\leq\bm{1}\\}$. We then apply Danskin’s Theorem
[3, Proposition 4.5.1] to derive the expression of the subdifferential. To
apply the theorem we must show that $\mathcal{Z}$ is compact. $\mathcal{Z}$ is
the intersection of closed half-spaces so is closed and, using [24, Corollary
14.5.1], it must also be bounded because $0\in\operatorname{int}\mathcal{V}$
and $\mathcal{V}$ is a C-set. From [3, Proposition 4.5.1], the subdifferential
is then given by the convex hull of all $h$ such that $h\in\mathcal{Z}$ and
maximize $h^{T}x$. The latter condition corresponds to
$h^{T}x=\psi_{\mathcal{V}}(x)$. Combining these constraints together we get
the set $\\{h:h^{T}x=\psi_{\mathcal{V}}(x),\,h^{T}V\leq\bm{1}^{T}\\}$, which
is already convex. This leads to (15). $\Box$

Given any direction $z\in\mathbb{R}^{n}$, $\psi_{\mathcal{V}}(x)$ also has a
well defined directional derivative [19, Theorem D1.2.2]

$\displaystyle\psi_{\mathcal{V}}^{\prime}(x;z)\\!$
$\displaystyle=\max_{h\in\partial\psi_{\mathcal{V}}(x)}\ h^{T}z$
$\displaystyle=\max_{h}\\{h^{T}z:h^{T}x=\psi_{\mathcal{V}}(x),\,h^{T}V\leq\bm{1}^{T}\\}$
(16a)
$\displaystyle=\min_{p,k}\\{\bm{1}^{T}p\\!+\\!\psi_{\mathcal{V}}(x)k:z\\!=\\!Vp\\!+\\!kx,p\\!\geq\\!0\\}$
(16b)

where (16b) is the dual problem to (16a).

We can integrate (16) to obtain the total change in $\psi_{\mathcal{V}}$ under
some dynamics $\dot{x}=f(x)$

$\psi_{\mathcal{V}}(x(t))-\psi_{\mathcal{V}}(x(0))=\int_{0}^{t}\psi_{\mathcal{V}}^{\prime}(x(\tau);f(x(\tau)))d\tau.$
(17)

This follows from [4, Theorem 2.10 and 2.11]; see also [25, Appendix 1 Theorem
4.3] and [24, Corollary 24.2.1].

Eq. 16a shows that the rate of change of V-rep polyhedral functions can be
characterized in terms of linear programming.

###  3.2 Polyhedral Gain Conditions

To compute bounds on incremental $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$
gains we will use the variational system of (6), along generic solution pairs
$(x,w):\mathbb{R}\to\mathbb{R}^{n}\times\mathbb{R}^{n_{w}}$ of (6):

$\dot{\delta x}=\partial f(x)\delta x+B\delta w,\,\delta z=C\delta x.$ (18)

Along a generic solution $(x,w)$, (18) corresponds to a time-varying linear
system with $n$ (variational) states, $n_{w}$ inputs, and $n_{z}$ outputs,
which captures the effect of (infinitesimal) perturbations on $(x,w)$ [11, 16,
17].

In what follows, we will take advantage of the fact that a uniform bound
$\gamma$ that satisfies $\|\delta z\|\leq\gamma\|\delta w\|$ for all
input/output pairs of the variational system and uniformly with respect to all
solution pairs $(x,w)$ of (6), is also a bound on the incremental gain of (6).
For instance, consider a generic pair of inputs signal $w_{0}$ and $w_{1}$,
and define the parameterized signal $w_{s}=sw_{1}+(1-s)w_{0}$, for $0\leq
s\leq 1$. Using (6), compute the state trajectory $x_{s}$ and the output
$z_{s}$ related to the input $w_{s}$ and to the initial condition
$x_{s}(0)=0$. For any $0\leq s\leq 1$, the quantities $\frac{d}{ds}w_{s}$ and
$\frac{d}{ds}z_{s}$ correspond to the input and output signals of (18)
computed along the solution pair $(x_{s},w_{s})$. Thus,

$\begin{split}\|z_{1}-z_{0}\|&=\left\|\int_{0}^{1}\frac{d}{ds}z_{s}ds\right\|\
\leq\int_{0}^{1}\left\|\frac{d}{ds}z_{s}\right\|ds\\\
&\leq\gamma\int_{0}^{1}\left\|\frac{d}{ds}w_{s}\right\|ds\
=\gamma\int_{0}^{1}\|w_{1}-w_{0}\|ds\\\ &=\gamma\|w_{1}-w_{0}\|,\end{split}$
(19)

where the second inequalities follows from the uniform bound $\|\delta
z\|\leq\gamma\|\delta w\|$.

Together with the use of the variational system, we assume the existence of a
V-rep polyhedral relaxation of its differential dynamics

$\mathcal{A}=\\{A_{1},\dots,A_{k}\\}\ \text{such that, for all $x$,}\ \partial
f(x)\in\operatorname{\mathbf{Co}}\mathcal{A},$ (20)

where

$\operatorname{\mathbf{Co}}\mathcal{A}=\left\\{\sum_{i=1}^{k}\lambda_{i}A_{i}:A_{i}\in\mathcal{A},\,\lambda_{i}\geq
0,\,\sum_{i=1}^{k}\lambda_{i}=1\right\\}.$

Namely, this means that, for any $x$, $\partial
f(x)=\sum_{i=1}^{k}\lambda_{i}(x)A_{i}$ where $A_{i}\in\mathcal{A}$,
$\lambda_{i}(x)\geq 0$, and $\sum_{i=1}^{k}\lambda_{i}(x)=1$. This relaxation
is inspired by a standard step in the literature for computation of nonlinear
gains by linear matrix inequalities (LMIs) [6]. We don’t present a general
algorithm for computing such relaxations since it is a non-trivial problem.
This is related to the trade-off between the quality of the relaxation and the
number of elements in it. Instead, in our examples, we work directly with sets
$\mathcal{A}$, with the results then applying to any nonlinear system that
satisfies (20).

For a linear system, (18) has the same dynamics as (6) with $f(x)=\partial
f(x)=A$ and so $\mathcal{A}=\\{A\\}$.

######  Theorem 1 (Bound on incremental $\mathcal{L}_{1}$ gain)

Consider variables
$\eta_{w}\in\mathbb{R},\,\eta_{z}\in\mathbb{R},\,V\in\mathbb{R}^{n\times
m},\,P\in\mathbb{R}^{m\times
2n_{w}},\,\mathcal{M}=\\{M_{1}\in\mathbb{R}^{m\times
m},\dots,M_{k}\in\mathbb{R}^{m\times m}\\}$ (where $k$ is the size of
$\mathcal{A}$) and vector valued function of $V$,

$\hat{z}_{V}=[|CV_{1}|_{1},\dots,|CV_{m}|_{1}]$

(where $V_{j}$ is the $j^{th}$ column of $V$).

An upper bound for the incremental $\mathcal{L}_{1}$ gain of (8) is given by

$\displaystyle\min_{\eta_{w},\eta_{z},V,P,\mathcal{M}}\quad$
$\displaystyle\>\eta_{w}/\eta_{z}$ $\displaystyle\begin{split}P\geq 0,\
[B,-B]=VP,&\ \eta_{w}\bm{1}^{T}=\bm{1}^{T}P\end{split}$ (21a)
$\displaystyle\eta_{z}>0,\ \operatorname{off-diag}(M_{i})\geq 0,\ A_{i}$
$\displaystyle V=VM_{i},\ -\eta_{z}\hat{z}_{V}=\bm{1}^{T}M_{i},$ for all
$\displaystyle i\in[k].$ (21b)

Moreover, the gauge function $\psi=\psi_{\mathcal{V}}/\eta_{z}$ satisfies (7)
for $\gamma_{1}=\eta_{w}/\eta_{z}$. $\lrcorner$

######  Theorem 2 (Bound on incremental $\mathcal{L}_{\infty}$ gain)

Consider variables
$\eta_{z}\in\mathbb{R},\,\eta_{w}\in\mathbb{R},\,H\in\mathbb{R}^{m\times
n},\,P\in\mathbb{R}^{2n_{z}\times
m},\,\mathcal{M}=\\{M_{1}\in\mathbb{R}^{m\times
m},\dots,M_{k}\in\mathbb{R}^{m\times m}\\}$ and vector valued function of $H$,

$\hat{w}_{H}=[|B^{T}(H^{1})^{T}|_{1};\dots;|B^{T}(H^{m})^{T}|_{1}]$

(where $H^{j}$ is the $j^{th}$ row of $H$).

An upper bound for the incremental $\mathcal{L}_{\infty}$ gain of (8) is given
by

$\displaystyle\min_{\eta_{z},\eta_{w},H,P,\mathcal{M}}\quad$
$\displaystyle\>\eta_{z}/\eta_{w}$ $\displaystyle\eta_{w}>0,\
\operatorname{off-diag}(M_{i})\geq 0,\ H$ $\displaystyle A_{i}=M_{i}H,\
-\\!\eta_{w}\hat{w}_{H}=M_{i}\bm{1},$ for all $\displaystyle i\in[k]$ (22a)
$\displaystyle\begin{split}P\geq 0,\ [C;-C]=PH,&\
\eta_{z}\bm{1}=P\bm{1}.\end{split}$ (22b)

Moreover, the polyhedral set $\mathcal{X}=\mathcal{H}/\eta_{w}$ satisfies (8)
for $\gamma_{\infty}=\eta_{z}/\eta_{w}$. $\lrcorner$

The columns of the matrix $V$ in (21) are the vertices of the level set of the
storage function. The intuition is that (21a) quantifies the impact that the
(incremental) input has on the growth of the storage function. In the proof of
Theorem 1 this is well captured by (24). Likewise, (21b) characterizes how the
internal dynamics of the system and its output features affect the storage, as
formalized below by (25). This splitting is made possible by the use of
storages with polyhedral representation. In a similar way, the rows of the
matrix $H$ in (22) characterize a polyhedral set on the system state space as
an intersection of half-spaces. The goal of (22a) is to bound the expansion of
this polyhedral set under the action of the forced system dynamics, in the
sense of (8a). Later, (22b) maps this bound on the system state into a bound
on the output, as clarified by (8b). We refer the reader to the proofs of the
theorems for a more detailed discussion on the connection between the
conditions of the theorems and the incremental gains of the system.

The optimization problems in (21) and (22) are adjoint. We mean this in the
sense that each solution of (21) for system matrices $\\{\mathcal{A},B,C\\}$
with variables $(\eta_{w},\eta_{z},V,P,\mathcal{M})$ is also a solution of
(22) for system matrices $\\{\mathcal{A}^{T},C^{T},B^{T}\\}$ (the adjoint
system) with variables $(\eta_{z},\eta_{w},V^{T},P^{T},\mathcal{M}^{T})$ and
vice versa. This matches the expected relationship between the two gains for
linear systems ((4) and (5)). It also means that it is sufficient to implement
only one of (21) or (22).

###### Remark 2

Because the approach relies on the set inclusion of the variational dynamics
(18) into a polytopic set $\operatorname{\mathbf{Co}}\mathcal{A}$, the same
methodology can also be directly applied to any system whose dynamics can be
characterized by polytopic linear differential inclusions (LDIs), including
linear time and parameter varying (LTV/LPV) systems.

###### Remark 3

Both (21) and (22) make use of two positive constants instead of one, which
might appear unexpected. This is because, in each case, one of the two
constants ($\eta_{w}$ in (21a) and $\eta_{z}$ in (22b)) act as normalizers for
the scale of the storage functions/sets. The normalizing constants can be
fixed to $1$ with no loss of generality but we keep the additional degrees of
freedom to better match the presentation in Section 4.

###  3.3 Proofs

_Proof of Theorem 1:_ Let $w=w_{1}-w_{2}$, $x=x_{1}-x_{2}$, and
$z=z_{1}-z_{2}$. Our aim is to show that every solution
$(\eta_{w},\eta_{z},V,P,\mathcal{M})$ of (21) satisfies, for all $w$ and $x$

$\dot{\psi}_{\mathcal{V}}\leq-\eta_{z}|z|_{1}+\eta_{w}|w|_{1}$ (23)

for system (6) and subject to (20), where

$\dot{\psi}_{\mathcal{V}}=\psi_{\mathcal{V}}^{\prime}(x;f(x_{1})-f(x_{2})+Bw).$

But, from (16a),

$\displaystyle\psi_{\mathcal{V}}^{\prime}(x;f(x_{1})-$ $\displaystyle
f(x_{2})+Bw)$ $\displaystyle\leq\psi_{\mathcal{V}}^{\prime}(x;$ $\displaystyle
f(x_{1})-f(x_{2}))+\psi_{\mathcal{V}}^{\prime}(x;Bw),$

so (23) is implied by

$(\forall w,x)\ \psi_{\mathcal{V}}^{\prime}(x;Bw)-\eta_{w}|w|_{1}\leq 0$ (24)

and

$(\forall x_{1},x_{2})\
\psi_{\mathcal{V}}^{\prime}(x;f(x_{1})-f(x_{2}))+\eta_{z}|Cx|_{1}\leq 0;$ (25)

which we address separately.

For (24), we have that

$\displaystyle\psi_{\mathcal{V}}^{\prime}(x;Bw)$
$\displaystyle=\max_{h}\\{h^{T}Bw:h^{T}x=\psi_{\mathcal{V}}(x),h^{T}V\leq\bm{1}^{T}\\}$
$\displaystyle\leq\max_{h}\\{h^{T}Bw:h^{T}V\leq\bm{1}^{T}\\}=\psi_{\mathcal{V}}(Bw).$

Treating the 1-norm as the polyhedral gauge function with V-rep matrix
$[I,-I]$, for any $w$ there exists a $p^{*}\geq 0$ such that $w=[I,-I]p^{*}$
and $\bm{1}^{T}p^{*}=|w|_{1}$. Then, from (21a), $Bw=[B,-B]p^{*}=VPp^{*}$,
$(Pp^{*})\geq 0$. Finally, from (11a) and (21a)

$\psi_{\mathcal{V}}(Bw)\leq\bm{1}^{T}Pp^{*}=\eta_{w}|w|_{1}.$

To show (25) we first bound $\psi_{\mathcal{V}}^{\prime}(x;f(x_{1})-f(x_{2}))$
in terms of $\mathcal{A}$. Recall that
$\partial\psi_{\mathcal{V}}(x)=\\{h:h^{T}x=\psi_{\mathcal{V}}(x),h^{T}V\leq\bm{1}^{T}\\}$.
Then

$\displaystyle\psi_{\mathcal{V}}^{\prime}(x;f(x_{1})-f(x_{2}))=\max_{h\in\partial\psi_{\mathcal{V}}(x)}\
h^{T}\big{(}f(x_{1})-f(x_{2})\big{)}$
$\displaystyle=\max_{h\in\partial\psi_{\mathcal{V}}(x)}\
h^{T}\left(\int_{0}^{1}\partial f\big{(}sx_{1}+(1-s)x_{2}\big{)}xds\right)$
$\displaystyle\leq\int_{0}^{1}\max_{h_{s}\in\partial\psi_{\mathcal{V}}(x)}\
h_{s}^{T}\big{(}\partial f(sx_{1}+(1-s)x_{2})x\big{)}ds$
$\displaystyle=\int_{0}^{1}\max_{h_{s}\in\partial\psi_{\mathcal{V}}(x)}\
h_{s}^{T}\left(\sum_{i=1}^{k}\lambda_{i}(x_{1},x_{2},s)A_{i}x\right)ds$
$\displaystyle\leq\int_{0}^{1}\sum_{i=1}^{k}\lambda_{i}(x_{1},x_{2},s)\max_{h_{i}\in\partial\psi_{\mathcal{V}}(x)}\
h_{i}^{T}\left(A_{i}x\right)ds$
$\displaystyle\leq\int_{0}^{1}\left(\sum_{i=1}^{k}\lambda_{i}(x_{1},x_{2},s)\right)\max_{i}\left(\max_{h_{i}\in\partial\psi_{\mathcal{V}}(x)}\
h_{i}^{T}A_{i}x\right)ds$
$\displaystyle=\max_{i}\left(\max_{h_{i}\in\partial\psi_{\mathcal{V}}(x)}\
h_{i}^{T}A_{i}x\right)=\max_{i}\left(\psi_{\mathcal{V}}^{\prime}(x;A_{i}x)\right).$

From (11a), for any $x$ there exists a $p^{*}\geq 0$ such that $x=Vp^{*}$ and
$\bm{1}^{T}p^{*}=\psi_{\mathcal{V}}(x)$. Combining with (21b) and rewriting
$M_{i}=Q_{i}+k^{*}I$ for some $Q_{i}\geq 0$ and some $k^{*}\in\mathbb{R}$, we
have $A_{i}x=A_{i}Vp^{*}=VM_{i}p^{*}=V(Q_{i}+k^{*}I)p^{*}=VQ_{i}p^{*}+k^{*}x$,
where $(Q_{i}p^{*})\geq 0$. Finally, from (16b), (21b), and the subadditivity
of the 1-norm

$\displaystyle\psi_{\mathcal{V}}^{\prime}(x;A_{i}x)$
$\displaystyle=\min_{p,k}\\{\bm{1}^{T}p\\!+\\!\psi_{\mathcal{V}}(x)k:A_{i}x\\!=\\!Vp\\!+\\!kx,p\\!\geq\\!0\\}$
$\displaystyle\leq\bm{1}^{T}Q_{i}p^{*}\\!+\\!\bm{1}^{T}p^{*}k^{*}=\bm{1}^{T}M_{i}p^{*}$
$\displaystyle=-\eta_{z}\sum_{j}|CV_{j}|_{1}p_{j}^{*}\leq-\eta_{z}|Cx|_{1}.$

$\Box$

_Proof of Theorem 2:_ This is a straightforward variant of a known result in
the closely related LTV/LPV setting [4], which we prove here for completeness.

Let $w=w_{1}-w_{2}$, $x=x_{1}-x_{2}$, and $z=z_{1}-z_{2}$. Our aim is to show
that every solution $(\eta_{z},\eta_{w},H,P,\mathcal{M})$ of (22) satisfies,
for all $t\geq 0$,

$\displaystyle\|w\|_{\infty}\leq\eta_{w},x(0)=0$ $\displaystyle\implies
x(t)\in\mathcal{H}$ (26a) $\displaystyle x(t)\in\mathcal{H}$
$\displaystyle\implies\|z\|_{\infty}\leq\eta_{z}.$ (26b)

We first show that (22a) implies (26a). Then we show that (22b) implies (26b).

To show (26a) we write the conditions for the forward invariance of
$\mathcal{H}$ subject to the incremental dynamics
$\dot{x}=f(x_{1})-f(x_{2})+Bw$ and $|w|_{\infty}\leq\eta_{w}$. Let
$\mathcal{C}_{j}=\\{x:H^{j}x=1,Hx\leq\bm{1}\\}$. Via Nagumo’s Theorem [4,
Section 4.2], (26a) holds if

$\max_{x_{1}-x_{2}\in\mathcal{C}_{j},|w|_{\infty}\leq\eta_{w}}\
H^{j}(f(x_{1})-f(x_{2}))+H^{j}Bw\leq 0$ (27)

for all $j\in[m]$. In fact, the meaning of the above inequality is that at
each point $x=x_{1}-x_{2}$ on the boundary of $\mathcal{H}$, the vector field
$\dot{x}=f(x_{1})-f(x_{2})+Bw$ is directed towards the inside of
$\mathcal{H}$. This guarantees that $\mathcal{H}$ is a forward invariant set
for the incremental dynamics for any $\|w\|_{\infty}\leq\eta_{w}$.

Eq. 27 holds if the following inequality holds

$\max_{x_{1}-x_{2}\in\mathcal{C}_{j}}\
H^{j}(f(x_{1})-f(x_{2}))\leq-\eta_{w}|B^{T}(H^{j})^{T}|_{1},$ (28)

since, from Hölder’s inequality,

$|H^{j}\\!Bw|\leq\sum|(H^{j}\\!B)^{T}_{i}w_{i}|\leq\eta_{w}|B^{T}\\!(H^{j})^{T}|_{1}\
.$

We now prove that (22a) implies (28) (and thus (26a)). First,

$\displaystyle\max_{x_{1}-x_{2}\in\mathcal{C}_{j}}\ H^{j}(f(x_{1})-f(x_{2}))$
$\displaystyle=\max_{x_{1}-x_{2}\in\mathcal{C}_{j}}\
H^{j}\left(\int_{0}^{1}\partial f(sx_{1}+(1-s)x_{2})xds\right)$
$\displaystyle=\max_{x_{1}-x_{2}\in\mathcal{C}_{j}}\
H^{j}\left(\int_{0}^{1}\sum_{i=1}^{k}\lambda_{i}(x_{1},x_{2},s)A_{i}xds\right)$
$\displaystyle=\max_{x_{1}-x_{2}\in\mathcal{C}_{j}}\
\left(\int_{0}^{1}\sum_{i=1}^{k}\lambda_{i}(x_{1},x_{2},s)H^{j}A_{i}xds\right)$
$\displaystyle\leq\max_{x_{1}-x_{2}\in\mathcal{C}_{j}}\
\int_{0}^{1}\left(\sum_{i=1}^{k}\lambda_{i}(x_{1},x_{2},s)\right)\max_{i}\left(H^{j}A_{i}x\right)ds$
$\displaystyle=\max_{i}\left(\max_{x\in\mathcal{C}_{j}}\ H^{j}A_{i}x\right).$

We then have that, through (22a) and by rewriting $M_{i}=Q_{i}+k^{*}I$ for
some $Q_{i}\geq 0$ and some $k^{*}\in\mathbb{R}$, for each $i$

$\displaystyle\max_{x\in\mathcal{C}_{j}}\ H^{j}A_{i}x$
$\displaystyle=\max_{x\in\mathcal{C}_{j}}\
(M_{i}Hx)_{j}=\max_{x\in\mathcal{C}_{j}}\ (Q_{i}Hx)_{j}+k^{*}H^{j}x$
$\displaystyle\leq(Q_{i}\bm{1})_{j}+k^{*}=(M_{i}\bm{1})_{j}=-\eta_{w}|B^{T}(H^{j})^{T}|_{1}.$

For (26b) we need to show that $\|z\|_{\infty}=\|Cx\|_{\infty}\leq\eta_{z}$
for all trajectories $x$ such that $x(t)\in\mathcal{H}$. Since $\mathcal{H}$
is forward invariant, this is verified if $[C;-C]x\leq\eta_{z}\bm{1}$ for all
$x$ such that $Hx\leq\bm{1}$. Indeed, from (22b),

$[C;-C]x=PHx\leq P\bm{1}=\eta_{z}\bm{1}.\vspace*{-5mm}$

$\Box$

###### Remark 4

Theorems 1 and 2 can be written in terms of the Jacobian $\partial f(x)$,
without explicitly taking into account any convex relaxation $\mathcal{A}$.
This would make both theorems more general but not tractable, since we would
have an infinite set of constraints to satisfy, one for each $x$. In practice,
to compute the gains, we would still need to rely on convex relaxations.

##  4 Computing Gain Bounds

###  4.1 Analysis

Because of the connection between (21) and (22), we only need one algorithm
for computing both gain bounds. We choose to design an algorithm for the
$\mathcal{L}_{1}$ problem and apply the algorithm to the adjoint system if we
are after $\mathcal{L}_{\infty}$ gains bounds. A comparison between our
algorithm and other approaches is provided in Section 4.3.

Because having variables multiplying each other makes a problem nonconvex and
nonlinear, directly solving (21) with the polyhedral candidate storage
function represented by $V$ as a variable is not computationally tractable.
However, for a fixed $V$, (21) can be efficiently solved using linear
programming (LP). This suggests a solution strategy based on an alternating
iteration to produce increasingly tight bounds on the gains by locally
optimizing a candidate $V$ (building on [21] that uses the same approach to
find polyhedral Lyapunov functions for autonomous systems). We alternate
between two steps, which we call the _gain estimation_ step and the
_polyhedral modification_ step. Both steps are LP problems. During the gain
estimation step we consider $V$ fixed and compute a bound on the gain. During
the polyhedral modification step, we introduce a small variation in our
candidate storage function representation $\delta V$ with the goal of enabling
a smaller gain bound to be verified.

Gain estimation step: after fixing $V$, variables $(\eta_{w},P)$ and
$(\eta_{z},\mathcal{M})$ of (21) become independent of each other. We can then
solve (21) by solving two LP sub-problems: minimize $\eta_{w}$ subject to
(21a) (with variables $\eta_{w},P$), and maximize $\eta_{z}$ subject to (21b)
(with variables $\eta_{z},\mathcal{M}$). The optimal $\eta_{w}/\eta_{z}$ is
then equal to the solution to (21) for the given $V$. Moreover, if we fix $V$
such that $\psi_{\mathcal{V}}$ is a decaying incremental Lyapunov function for
the internal dynamics, (21b) is guaranteed to be feasible (since there must
then exist $M_{i}$ such that $\operatorname{off-diag}(M_{i})\geq
0,A_{i}V=VM_{i}$ and $1^{T}M_{i}<0$) and so the gain estimation step is
guaranteed to obtain a finite gain bound (since (21a) is always feasible and
implies $\eta_{w}>0$).

Polyhedral modification step: we replace our objective $\min\eta_{w}/\eta_{z}$
with $\min\log(\eta_{w}/\eta_{z})=\min\\{\log(\eta_{w})-\log(\eta_{z})\\}$ to
obtain a simpler problem. We then search for a small perturbation around
solution $(\eta_{w},\eta_{z},V,P,\mathcal{M})$ in the direction that most
improves the linear approximation of our objective. Working with small
perturbations allows us to express the modification step as a LP in the
(variational) variables
$\delta\eta_{w}\in\mathbb{R},\,\delta\eta_{z}\in\mathbb{R},\,\delta
V\in\mathbb{R}^{n\times m},\,\delta P\in\mathbb{R}^{m\times
2n_{w}},\,\delta\mathcal{M}=\\{\delta M_{1}\in\mathbb{R}^{m\times
m},\dots,\delta M_{k}\in\mathbb{R}^{m\times m}\\}$ (where $k$ is the size of
$\mathcal{A}$) and vector valued function of $V$,

$\delta\hat{z}_{V}=[\operatorname{sgn}(CV_{1})^{T}C\delta
V_{1},\dots,\operatorname{sgn}(CV_{m})^{T}C\delta V_{m}]$

(where $V_{j}$ and $\delta V_{j}$ are the $j^{th}$ columns of $V$ and $\delta
V$ respectively):

$\displaystyle\min_{\delta\eta_{w},\delta\eta_{z},\delta V,\delta
P,\delta\mathcal{M}}\quad\>\delta\eta_{w}/\eta_{w}$
$\displaystyle-\delta\eta_{z}/\eta_{z}$ (29a)
$\displaystyle\begin{split}P+\delta P\geq 0,\ 0=\delta VP+V\delta P,\
\delta\eta_{w}&\bm{1}^{T}=\bm{1}^{T}\delta P\end{split}$ (29b)
$\displaystyle\begin{split}\operatorname{off-diag}(M_{i}+\delta M_{i})\geq
0,\,A_{i}\delta V=\delta V&M_{i}+V\delta M_{i},\\\ \eta_{z}+\delta\eta_{z}>0,\
-(\delta\eta_{z}\hat{z}+\eta_{z}\delta\hat{z})&=\bm{1}^{T}\delta M_{i},\\\
\text{for all }&i\in[k]\end{split}$ (29c)
$\displaystyle|\operatorname{vec}(\delta V)|_{\infty}$
$\displaystyle\leq\varepsilon.$ (29d)

Breaking (29) into four parts, (29a) is the linear approximation of the
variation in $\delta\log(\eta_{w}/\eta_{z})$, while (29b) and (29c)
approximate (21a) and (21b) respectively. Eq. 29d ensures that the variations
considered are small because it also bounds all other variables.

To account for the effect of the approximations and ensure that the gain bound
remains sound we only apply the variation to the candidate storage function to
obtain $V_{\textrm{new}}:=V+\delta V$ and then solve the gain estimation
problem again for $V_{\textrm{new}}$. For small enough $\varepsilon$, this is
guaranteed to lead to an improved gain bound. We set $\varepsilon$ at the
$N^{th}$ iteration of our algorithm to $\varepsilon_{0}/N$, where
$\varepsilon_{0}$ is the initial step-size. Additionally, if, after a
polyhedral modification step, a worse gain bound is computed during the gain
estimation step or there exist redundant vertices in the polyhedral set we
halve $\varepsilon$ and repeat the polyhedral modification step (a vertex
$V_{j}$ is redundant if $\psi_{\mathcal{V}}(V_{j})<1$). We stop our algorithm
when $\varepsilon$ becomes smaller than a threshold.

Finally, to bootstrap the gain estimation step to a feasible solution we first
use the algorithm from [21] to obtain a candidate $V$ that represents a
polyhedral incremental Lyapunov function for the internal dynamics. This stage
also controls the complexity of Lyapunov functions, given by the number of
vertices of $\mathcal{V}$ (equivalently, by the number of columns of $V$).
This number, $m$, is an input to the algorithm. In general, internal dynamics
that are more oscillatory will require functions with higher complexity [2]
(so if $m$ is set too low the algorithm in [21] might not be able to find a
polyhedral Lyapunov function). We also observe that higher complexity
functions are better able to approximate arbitrary gauge functions and so can
provide tighter bounds, but this comes at an increased computational cost. We
briefly explore this trade-off in Section 5.1.

The full algorithms for the $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$
incremental gains are summarized below.

Data: Matrices $(\mathcal{A},B,C)$,
polyhedral function complexity $m$,
and initial and minimum step-size parameters
$(\varepsilon_{0},\varepsilon_{min})$
Result: Upper bound on $\mathcal{L}_{1}$ incremental gain,
and corresponding storage function representation $V$
Procedure:
$(\varepsilon,N,\gamma^{\prime},\delta V):=(\varepsilon_{0},1,\infty,0)$
$V:=V$ from [21] LDI algorithm with inputs $(\mathcal{A},m)$
while True:
_Gain Estimation Step:_
$\eta_{w},P:=\operatorname{argmin}_{\eta_{w},P}\,\\{\eta_{w}:\eqref{eq:incr-
lone-hard-input}\\}$ using current $V$
$\eta_{z},\mathcal{M}:=\operatorname{argmax}_{\eta_{z},\mathcal{M}}\,\\{\eta_{z}:\eqref{eq:incr-
lone-hard-state}\\}$ using current $V$
if $\eta_{w}/\eta_{z}>\gamma^{\prime}$ or $\texttt{containsRedundant}(V)$:
$\varepsilon,V:=(\varepsilon/2,V-\delta V)$
else:
$(\varepsilon,N,\gamma^{\prime}):=(\varepsilon_{0}/N,N+1,\eta_{w}/\eta_{z})$
if $\varepsilon<\varepsilon_{min}$:
return $(\gamma^{\prime},V)$
_Polyhedral Modification Step:_
$V:=V+\delta V$ from (29) using
$(\eta_{w},\eta_{z},\varepsilon,V,P,\mathcal{M})$

Algorithm 1 Incremental $\mathcal{L}_{1}$ gain bound computation

Data: Matrices $(\mathcal{A},B,C)$,
polyhedral function complexity $m$,
and initial and minimum step-size parameters
$(\varepsilon_{0},\varepsilon_{min})$
Result: Upper bound on $\mathcal{L}_{\infty}$ gain,
and corresponding invariant set representation $H$
Procedure:
$\gamma_{\infty},H^{T}:=$ Algorithm 1 with matrix inputs
$(\mathcal{A}^{T},C^{T},B^{T})$
return $(\gamma_{\infty},H)$

Algorithm 2 Incremental $\mathcal{L}_{\infty}$ gain bound computation

###  4.2 Synthesis

After fixing $V$ the constraints in (21b) are convex with respect to the
relaxed dynamics $\mathcal{A}$. This makes it easy to adapt the algorithms of
Section 4.1 to linear feedback control synthesis by adding extra free
variables corresponding to the control parameters.

We consider (6) with added controlled inputs $u$ and measured outputs $y$

$\displaystyle\dot{x}$ $\displaystyle=f(x)+Bw+B_{u}u,\,$ $\displaystyle z$
$\displaystyle=Cx,$ $\displaystyle y$ $\displaystyle=C_{y}x.$

State feedback design is obtained by taking $C_{y}=I$ but our procedure works
with any output matrix, allowing for the design of static output feedback
controllers of the form $u=Ky$.

The procedure is based on the extension of the gain estimation step of Section
4.1, by replacing each constraint $A_{i}V=VM_{i}$ with

$(A_{i}+B_{u}KC_{y})V=VM_{i}.$ (30)

Since $K$ enters (30) linearly (for fixed $V$), it is also straightforward to
consider additional linear constraints on $K$ with the goal of enforcing
bounds on the control gains or sparsity patterns. By adding sums of
constraints of the form $\Gamma_{l}^{T}K\Gamma_{r}\leq g$, where $\Gamma_{l}$
and $\Gamma_{r}$ are matrices of suitable size, we can limit the range of
(linear combinations of) the coefficients of $K$. Likewise, constraints of the
form $\Gamma_{l}^{T}K\Gamma_{r}=0$ can be used to enforce specific sparsity
patterns on $K$.

The introduction of the additional variable $K$ also imposes minor changes to
the polyhedral modification step. The small perturbations in (29) must take
into account the extended solution manifold

$(\eta_{w},\eta_{z},V,P,\mathcal{M},K).$

This means that we need to add $\delta K$ as a new variable to (29) and
replace the constraints $A_{i}\delta V=\delta VM_{i}+V\delta M_{i}$ with

$(A_{i}+B_{u}KC_{y})\delta V+B_{u}\delta KC_{y}V=\delta VM_{i}+V\delta M_{i}.$
(31)

We also add constraint $\|\operatorname{vec}(\delta
K)\|_{\infty}\leq\varepsilon$ to guarantee that only small perturbations are
considered. Finally, we add constraints $\Gamma_{l}^{T}(K+\delta
K)\Gamma_{r}\leq g$ and $\Gamma_{l}^{T}(K+\delta K)\Gamma_{r}=0$ for each
constraint $\Gamma_{l}^{T}K\Gamma_{r}\leq g$ and
$\Gamma_{l}^{T}K\Gamma_{r}=0$, respectively.

###  4.3 Comparison with Other Approaches

The use of polyhedral sets for the computation of gains is not new. They have
been used in [4, Section 6.4] and [15] for the computation of the
$\mathcal{L}_{\infty}$ gain of discrete time LTV/LPV systems.

In [4], the algorithm computes the largest set included in
$\|Cx\|_{\infty}\leq\mu$ that is invariant under all possible
$\|u\|_{\infty}\leq 1$. This is done by starting with polyhedral set
$\|Cx\|_{\infty}\leq\mu$ and propagating the set backwards in time until the
set becomes invariant, in which case $\mu$ is an upper bound for the
$\mathcal{L}_{\infty}$ gain, or no longer contains $0$, in which case $\mu$ is
a lower bound for the $\mathcal{L}_{\infty}$ gain. This is repeated for
different values of $\mu$ to obtain a sharp bound.

The approach of [15] computes the reachable set from the origin by forward
propagation in time. After each step, the smallest $\mu$ for which set
$\|Cx\|_{\infty}\leq\mu$ contains the reachable set is a lower bound on the
gain. [15] also provide bounds on the error between the computed lower bound
and the actual gain of the system, as a function of the number of forward
steps taken and of the system spectral radius. This allows the algorithm to
terminate in a pre-determined finite number of steps as a function of the
required accuracy. A drawback of the approach is that it scales worse for
inputs of large dimension. This is because V-rep polyhedra must be used
throughout, including to represent the $\mathcal{L}_{\infty}$ unit ball of the
image of the input $Bw$, which by itself requires $2^{n_{w}}$ vertices.

The approaches of both [4] and [15] can be applied to continuous time systems
through a forward Euler discretization scheme (from $\dot{x}=Ax$ to
$x^{+}=[I+\tau A]x$). Discretizing the systems in this way guarantees that the
set invariance properties necessary to bound the $\mathcal{L}_{\infty}$ gain
will hold for the continuous system if they hold for the discretized system
[4]. By using smaller discretization time constants the bound on the gain can
then be made increasingly tighter, but at an increasing computational cost.
This is a sharp difference with our approach, which is based on conditions
that apply directly to continuous time systems, without the need of any
discretization step.

The approaches of [4] and [15] iteratively shrink or grow polyhedra at every
step of computation. This may lead to increasingly large matrices to represent
the polyhedra, with the matrix size growing exponentially with the number of
steps taken. In contrast, our approach does not shrink or grow polyhedra but
reshapes them at each step, keeping constant the size of the representation
matrix. This is important since the computational cost at each iteration
scales with the size of the polyhedral representation. As such, while for the
approaches of [4] and [15] there is guaranteed convergence to the actual
LTV/LPV gain in a finite number of steps for discrete or discretized systems,
each step could become computationally intense, with no upper bound obtained
until the process terminates. In contrast, there is no guaranteed convergence
to the actual LTV/LPV gain for our approach but each step always has a fixed
cost and produces a tighter upper bound. The complexity of the polyhedra,
characterized by the number of vertices $m$, is governed by the user.

Another important unique strength of our optimization-based approach is that
it allows for synthesis problems to be tackled directly.

##  5 Examples

###  5.1 DC Motor

We revisit the DC motor example of [21] for the computation of incremental
system gains. In the first part, we estimate a bound on the incremental
$\mathcal{L}_{\infty}$ gain between a disturbance acting on an open loop
(linear) DC motor and its output speed. Later, we design a feedback controller
that stabilizes the motor position while minimizing the effect of disturbances
in closed loop.

The nominal parameters are: inertia $J_{0}=.01$, viscous friction constant
$b_{0}=.1$, electromagnetic field constant $E_{0}=.01$, resistance $R_{0}=1$,
and inductance $L_{0}=.5$. In each case, we run the algorithm 10 times with
$\epsilon_{0}=1/5$ and different random seeds for $V$ and report the result of
the best run.

For the analysis scenario, the variational dynamics of the DC motor are

$\displaystyle\dot{\delta x}$ $\displaystyle=\begin{bmatrix}-b/J&E/J\\\
-E/L&-R/L\end{bmatrix}\delta x+\begin{bmatrix}0\\\ 1\end{bmatrix}\delta w$
(32) $\displaystyle\delta z$
$\displaystyle=\begin{bmatrix}1&0\end{bmatrix}\delta x,$

where the first state corresponds to the speed of the motor and the second to
the current. We first consider the case where each parameter in (32) is set to
its nominal value, $b=b_{0}$, $E=E_{0}$, $J=J_{0}$, and $R=R_{0}$. We apply
Algorithm 2 to find H-rep sets with three and four half-spaces (shown in Fig.
1) for which the bounds obtained from Theorem 2 are around $0.083$ and $0.050$
respectively. In this case, because the impulse response is positive, the
actual $\mathcal{L}_{\infty}$ gain can be analytically computed to be
$1/20.02$ using Remark 1 and matches the bound for the four half-space set to
4 significant figures.

Figure 1:  Three and four sided H-rep sets found using Algorithm 2 bounding
the $\mathcal{L}_{\infty}$ gain of (32) (with nominal parameter values) to
$\leq 0.083$ and $\leq 0.050$ respectively.

We then consider (32) with uncertain parameters

$J\\!\in\\![J_{0}/\delta_{s},\delta_{s}J_{0}],\,b\\!\in\\![b_{0}/\delta_{s},\delta_{s}b_{0}],\,E\\!\in\\![E_{0}/\delta_{s},\delta_{s}E_{0}]$
(33)

for $\delta_{s}=8$, leading to a relaxation set $\mathcal{A}$ composed of $8$
vertices. Because this is not a linear system, the incremental
$\mathcal{L}_{\infty}$ gain can’t be computed analytically in this case. A
lower bound for it can be obtained by computing the gain for linear systems
with fixed parameters picked from the vertices of the sets in (33). Out of
these, the highest occurs when
$(J,b,E)=(J_{0}/\delta_{s},b_{0}/\delta_{s},\delta_{s}E_{0})$ and can be
analytically computed to be around $2.12$. We use Algorithm 2 to compute an
upper bound. We work with H-rep sets with four, six, eight, and ten sides,
leading to upper bounds on the incremental gain of $6.6$, $5.2$, $4.8$, and
$4.4$ respectively. This shows how increasing the complexity of the polyhedral
representation leads to less conservative results. The polyhedra found are
shown in Fig. 2.

Figure 2:  Four, six, eight, and ten sided H-rep polyhedral sets found using
Algorithm 2 bounding the incremental $\mathcal{L}_{\infty}$ gain of (32)-(33)
(with $\delta_{s}=8$) to $6.6$, $5.2$, $4.8$, and $4.4$ respectively.

For the synthesis scenario we design a linear state feedback controller and a
linear output feedback controller, the latter using motor position and motor
current. The goal is to minimize the worse case effect that bounded
disturbance forces have on the motor position.

The variational dynamics reads

$\displaystyle\dot{\delta x}$ $\displaystyle=\begin{bmatrix}0&1&0\\\
0&-b\\!/J&E\\!/J\\\ 0&-E\\!/L&-R\\!/L\end{bmatrix}\delta
x+\\!\begin{bmatrix}0\\\ 0\\\ 1\\!/L\end{bmatrix}\delta
u+\\!\begin{bmatrix}0\\\ 1\\\ 0\end{bmatrix}\delta w$ (34)
$\displaystyle\delta z$ $\displaystyle=\begin{bmatrix}1&0&0\end{bmatrix}\delta
x,$ (35)

where the first and second states are related to the motor position and
velocities, and the last state is related to the motor current. We also
consider the following parametric uncertainties

$J\in[J_{0}/\delta_{p},\delta_{p}J_{0}],\,E\in[E_{0}/\delta_{p},\delta_{p}E_{0}]$
(36)

where $\delta_{p}=1.4$. We first look into state feedback, using the output

$y_{s}=-Ix.$ (37)

We then move to output feedback design, from

$y_{o}=-\begin{bmatrix}1&0&0\\\ 0&0&1\end{bmatrix}x.$ (38)

In both case, we adopt a polyhedral representation given by 12 half-spaces.
With the modified conditions of Section 4.2, our algorithm proves a closed-
loop $\mathcal{L}_{\infty}$ bound of $0.17$, for the state feedback controller

$u=\begin{bmatrix}27.5&3.04&1.41\end{bmatrix}y_{s}.$ (39)

Likewise, the algorithm derives a closed-loop $\mathcal{L}_{\infty}$ gain of
$0.26$ for the output feedback controller

$u=\begin{bmatrix}21.8&1.74\end{bmatrix}y_{o}.$ (40)

The associated polyhedra are shown in Fig. 3.

Figure 3:  12 sided H-rep sets found simultaneously to feedback schemes
(37),(39) (left) and (38),(40) (right) using Algorithm 2 adapted to synthesis.
They bound the incremental $\mathcal{L}_{\infty}$ gain for (34)-(36) (with
$\delta_{p}=1.4$) to $\leq 0.17$ and $\leq 0.26$ respectively.

###  5.2 Compartmental System

We next consider a small compartmental system example. Compartmental systems
are naturally suited for pharmacological and epidemiological models, water
reservoirs and queuing systems, and as discretizations of partial differential
equations (PDEs) [18]. Each compartment can have inflows and outflows, both
internal and external. In our example we consider a basic compartmental system
where matter moves among compartments and its quantity is conserved through
internal flows but can be added or removed from the overall system through
external inflows or outflows. We will once again be running the algorithm 10
times with $\varepsilon_{0}=1/5$ and different random seeds for $V$ and
reporting the result of the best run.

Specifically, we consider the design of output feedback controllers with
integral action for a compartmental model that approximates diffusion through
a pipe. The pipe is modelled by 5 compartments. The size of the state space of
the closed-loop system is $6$, with the sixth state needed for the integral
action of the controller. The objective is to minimize the effects of external
disturbances on the measured output associated to the last compartment. We
express this as incremental $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gain
synthesis problems for system

$\begin{split}\dot{x}_{1}&=-\sigma x_{1}+x_{2}+w_{1}+u\\\ \dot{x}_{2}&=\sigma
x_{1}-2x_{2}+x_{3}+w_{2}\\\ \dot{x}_{3}&=x_{2}-2x_{3}+x_{4}+w_{3}\\\
\dot{x}_{4}&=x_{3}-2x_{4}+x_{5}+w_{4}\\\ \dot{x}_{5}&=x_{4}-2x_{5}+w_{5}\\\
\dot{q}&=r-x_{5}\\\ z&=x_{5}-r,\end{split}$ (41)

where each state $x_{i}$ corresponds to the matter in the $i$th compartment,
and $q$ corresponds to the integral error $r-x_{5}$, for the desired reference
value $r$. The input $u$ represents the rate at which matter is being injected
at the beginning of the pipe. We also set $\sigma\in[0.5,1]$, encompassing
dynamics where the flow from the first compartment to the second becomes
partially obstructed. Matter is conserved everywhere except in the last
compartment, where it is escaping with a rate of $1$. Variables $w_{i}$
represent independent external disturbances that are acting on each
compartment. The proportional-integral controllers $u=Ky$ use the measured
output

$y=\begin{bmatrix}-x_{1}&r-x_{5}&q\end{bmatrix}^{T}.$ (42)

The variational dynamics associated to the system reads

$\begin{split}\dot{\delta x}_{1}&=-\sigma\delta x_{1}+\delta x_{2}+\delta
w_{1}+\delta u\\\ \dot{\delta x}_{2}&=\sigma\delta x_{1}-2\delta x_{2}+\delta
x_{3}+\delta w_{2}\\\ \dot{\delta x}_{3}&=\delta x_{2}-2\delta x_{3}+\delta
x_{4}+\delta w_{3}\\\ \dot{\delta x}_{4}&=\delta x_{3}-2\delta x_{4}+\delta
x_{5}+\delta w_{4}\\\ \dot{\delta x}_{5}&=\delta x_{4}-2\delta x_{5}+\delta
w_{5}\\\ \dot{\delta q}&=-\delta x_{5}\\\ \delta z&=\delta x_{5},\end{split}$
(43)

which, for $\sigma\in[0.5,1]$, leads to a relaxation set $\mathcal{A}$
composed of 2 vertices.

We use Algorithms 1 and 2 to design controllers that minimize the incremental
$\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains from the disturbance input
$w$ to the performance output $z$. The results are summarized in Table 1.

| $K$ | Upper Bound | Lower Bound
---|---|---|---
$\mathcal{L}_{1}$ | $\begin{bmatrix}2.51&1.28&0.27\end{bmatrix}$ | 2.47 | 1.47
$\mathcal{L}_{\infty}$ | $\begin{bmatrix}9.10&6.45&0.23\end{bmatrix}$ | 5.06 | 3.54
Table 1: Table of control parameters and incremental gain bounds found using
Algorithms 1 and 2 (with $m=18$) for minimizing each target gain bound for the
linear case. The lower bounds are computed from the vertices of the closed-
loop linear systems satisfying (43) with $\sigma\in[.5,1]$ using (4) and (5)
combined with [26].

Both performance specifications have interesting applications. For example,
the $\mathcal{L}_{\infty}$ gain bound bounds the maximum amount by which the
regulated output will vary about its reference value when subjected to
independent bounded disturbances on each of the compartments.

Beyond their direct interpretations in terms of performances, the gains also
quantify the robustness of the closed-loop system to unmodeled dynamics, via
the small gain theorem. For example, the gain in Table 1 guarantees closed-
loop incremental stability even if the last compartment is leaking at a rate
$\in[0.6,1.4]$. This is obtained by setting $\delta w_{5}$ to
$\pm(1/\gamma)\delta x_{5}$, where $\gamma$ is the smallest upper bound in
Table 1.

##  6 Conclusions

We have discussed theoretical conditions for the estimation of bounds on the
incremental $\mathcal{L}_{1}$ and $\mathcal{L}_{\infty}$ gains of nonlinear
systems. Their computation is performed via a two-step procedure based on
linear programming problems. To the best of our knowledge, the
characterization of $\mathcal{L}_{1}$ gain bounds is a novel contribution
while the $\mathcal{L}_{\infty}$ bounds are closer to results available in the
literature, based on set invariance. The algorithms outlined in the paper are
sound but susceptible to local optimal. They can be used both for system
analysis and for closed-loop control synthesis. The gains computed can
quantify performance specifications or be used in conjunction with the small
gain theorem to prove stability of interconnections. The practical use of the
algorithms was illustrated by two numerical examples.

More work is needed on optimizing, tweaking, and formally evaluating the
properties of the optimization procedure used in the algorithm. From the modes
of failure observed in sub-optimal solutions, it seems that explicitly
enforcing symmetry and adding the dual constraints of the solutions from the
gain estimation step to the polyhedral modification step helps avoid a big
class of local optima. Other extensions include weighted norms and synthesis
of nonlinear controllers [5]. This is particularly relevant because it is well
known that linear state feedback controllers are sub-optimal for the
$\mathcal{L}_{\infty}$ problem [13].

## References

  * [1] V. Balakrishnan and S. Boyd.  On computing the worst-case peak gain of linear systems.  Systems & Control Letters, 19(4):265–269, 1992.
  * [2] L. Benvenuti and L. Farina.  Eigenvalue regions for positive systems.  Systems & Control Letters, 51(3-4):325–330, 2004.
  * [3] D. Bertsekas, A. Nedic, and A. Ozdaglar.  Convex Analysis and Optimization.  Athena Scientific, 2003.
  * [4] F. Blanchini and S. Miani.  Set-Theoretic Methods in Control.  Systems & Control: Foundations & Applications. Birkhäuser, Cham, second edition, 2015.
  * [5] F. Blanchini and M. Sznaier.  Persistent disturbance rejection via static-state feedback.  IEEE Transactions on Automatic Control, 40(6):1127–1131, 1995.
  * [6] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan.  Linear Matrix Inequalities in System and Control Theory.  Studies in Applied and Numerical Mathematics. Society for Industrial and Applied Mathematics, 1994.
  * [7] S. Boyd and L. Vandenberghe.  Convex Optimization.  Cambridge University Press, Cambridge, UK ; New York, 2004.
  * [8] S. P. Boyd and C. H. Barratt.  Linear Controller Design: Limits of Performance.  Prentice Hall, Englewood Cliffs, N.J, 1991.
  * [9] C. Briat.  Robust stability analysis of uncertain linear positive systems via integral linear constraints: L1- and L$\infty$-gain characterizations.  In 2011 50th IEEE Conference on Decision and Control and European Control Conference, pages 6337–6342, 2011.
  * [10] C. Briat.  Robust stability and stabilization of uncertain linear positive systems via integral linear constraints: L1-gain and L$\infty$-gain characterization.  International Journal of Robust and Nonlinear Control, 23(17):1932–1954, 2013.
  * [11] P. E. Crouch and A. J. van der Schaft.  Variational and Hamiltonian Control Systems.  Number 101 in Lecture Notes in Control and Information Sciences. Springer, Berlin, 1987.
  * [12] C. A. Desoer and M. Vidyasagar.  Feedback Systems: Input–Output Properties.  Classics in Applied Mathematics. Elsevier, 1975.
  * [13] I. J. Diaz-Bobillo and M. A. Dahleh.  State feedback l1-optimal controllers can be dynamic.  Systems & Control Letters, 19(2):87–93, 1992.
  * [14] Y. Ebihara, D. Peaucelle, and D. Arzelier.  L1 gain analysis of linear positive systems and its application.  In 2011 50th IEEE Conference on Decision and Control and European Control Conference, pages 4029–4034, 2011.
  * [15] I. Fialho and T. Georgiou.  On the L1 norm of uncertain linear systems.  IEEE Transactions on Automatic Control, 40(6):1142–1147, 1995.
  * [16] F. Forni and R. Sepulchre.  A Differential Lyapunov Framework for Contraction Analysis.  IEEE Transactions on Automatic Control, 59(3):614–628, 2014.
  * [17] F. Forni, R. Sepulchre, and A. van der Schaft.  On differential passivity of physical systems.  In 52nd IEEE Conference on Decision and Control, pages 6580–6585, 2013.
  * [18] W. M. Haddad, V. Chellaboina, and Q. Hui.  Nonnegative and Compartmental Dynamical Systems.  Princeton University Press, 2010.
  * [19] J.-B. Hiriart-Urruty and C. Lemaréchal.  Fundamentals of Convex Analysis.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2001.
  * [20] H. K. Khalil.  Nonlinear Systems.  Prentice Hall, 2002.
  * [21] D. Kousoulidis and F. Forni.  Polyhedral Lyapunov Functions with Fixed Complexity.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 3293–3298, 2021.
  * [22] A. Linnemann.  Computing the L1-Norm of Continuous-Time Linear Systems.  SIAM Journal on Control and Optimization, 46(6):2052–2070, 2007\.
  * [23] A. Rantzer and M. E. Valcher.  Scalable Control of Positive Systems.  Annual Review of Control, Robotics, and Autonomous Systems, 4(1):319–341, 2021.
  * [24] R. T. Rockafellar.  Convex Analysis.  Princeton University Press, 1970.
  * [25] N. Rouche, P. Habets, and M. Laloy.  Stability Theory by Liapunov’s Direct Method, volume 22 of Applied Mathematical Sciences.  Springer New York, New York, NY, 1977.
  * [26] N. K. Rutland and P. G. Lane.  Computing the 1-norm of the impulse response of linear time-invariant systems.  page 11, 1995.
  * [27] J. C. Willems.  Dissipative dynamical systems part I: General theory.  Archive for Rational Mechanics and Analysis, 45(5):321–351, 1972\.
  * [28] K. Zhou, J. C. Doyle, and K. Glover.  Robust and Optimal Control.  Pearson, Upper Saddle River, N.J, 1st edition edition, 1995.