We want to stabilize this system using a dynamic output feedback controller of
order two (i.e., $n_{c}=2$) but with the restriction that the matrices $B_{c}$
and $D_{c}$ are diagonal. We again start with creating a representation for
the open-loop system using `tds_create`.

⬇

A0=[-3.5 -2 3 -1.5;3 -6.5 3 -1;4 -2 -7 2.5;-2.5 1.5 -3.5 1];

A1=[-0.25 0.5 -0.25 0.25;0.5 -0.5 1 -0.5;

1 -1 0.5 -0.15;-0.15 -0.15 0.1 0.2];

B = [-1 2;-2 1;3 -2.5;-4 0];

C = [0 1 1 0;1 0 0 1];

sys = tds_create({A0,A1},[0 1],{B},5,{C},0);

Next, we create the `mask` for our controller, i.e., an `ss`-object with
entries `1` and `0` to indicate which entries are allowed to be optimized and
which should remain fixed.

⬇

mask = ss(ones(2),eye(2),ones(2),eye(2));

To specify the basis values for the fixed entries, we create a `ss`-object
`basis`, whose elements dictate the values of the fixed entries.

⬇

basis = ss(zeros(2),zeros(2),zeros(2),zeros(2));

Next we run the optimization procedure starting from the initial controller

$A_{c}=-I_{2}\text{ and }B_{c}=C_{c}=D_{c}=I_{2}.$

⬇

initial = ss(-eye(2),eye(2),eye(2),eye(2));

options_stabopt = tds_stabopt_options('nstart',1);

[K,cl]␣=␣tds_stabopt_dynamic(sys,2,'options',options_stabopt,...

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣'initial',initial,'basis',basis,'mask',mask);

Inspecting the resulting controller using the matlab Command Window, we see
that the resulting controller indeed satisfies the specified pattern.

$\displaystyle A_{c}=\begin{bmatrix}-9.2431&11.7967\\\
0.8572&-3.1225\end{bmatrix}\text{, }B_{c}=\begin{bmatrix}-2.5544&0\\\
0&0.2507\end{bmatrix}\text{, }$ $\displaystyle
C_{c}=\begin{bmatrix}0.8863&-2.3617\\\ -15.8203&12.1480\end{bmatrix}\text{ and
}D_{c}=\begin{bmatrix}0.8420&0\\\ 0&0.0265\end{bmatrix}.$

Analyzing the resulting closed-loop system, we see that it is (strongly)
exponential stable as its (strong) spectral abscissa is equal to -0.1786. ∎

####  3.3.5 PID controller

Proportional-Integral-Derivative (PID) controllers are frequently used in
industry to stabilize a set-point. The feedback signal of a PID controller
consists of an affine combination of a term proportional to the difference
between the output and the desired set point, a term proportional to the
accumulation of the error over time (i.e., its integral), and a term
proportional to the rate of change of this error (i.e., its time-derivative):

$u(t)=K_{p}\,y(t)+K_{i}\int_{0}^{t}y(s)\,\mathop{}\\!\mathrm{d}s+K_{d}\,\dot{y}(t).$
(3.22)

The derivative component improves the responsiveness of the controller as it
allows to anticipate the future value of the error. The integrator on the
other hand allows to eliminate a steady-state error for a non-zero reference
signal. In the following example we will demonstrate how the design of a PID
controller can be transformed into the design of a dynamic output feedback
controller of the form (3.2) by augmenting the state and output variables and
by imposing a certain structure on some of the controller matrices.

######  Example 3.8.

Consider the following time-delay system:

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&\underbrace{\begin{bmatrix}-3&\phantom{-}2&\phantom{-}6\\\
-3&-4&\phantom{-}1\\\
\phantom{-}1&-1&-2\end{bmatrix}}_{A_{0}}x(t)+\underbrace{\begin{bmatrix}\phantom{-}0.5&-0.7&-0.8\\\
-0.3&\phantom{-}0.6&\phantom{-}0.2\\\
-0.2&\phantom{-}0.6&-0.4\end{bmatrix}}_{A_{1}}x(t-1)+\underbrace{\begin{bmatrix}1&\phantom{-}0\\\
1&\phantom{-}1\\\ 1&-1\end{bmatrix}}_{B}u(t-0.2)\\\\[28.45274pt]
y(t)&=&\underbrace{\begin{bmatrix}1&1&1\end{bmatrix}}_{C}x(t).\end{array}\right.$
(3.23)

To express a controller of the form (3.22) as a dynamic output feedback
controller, we first need to extend the state-space model (3.23). As in
Example 3.6, we will use a DDAE reformulation of the plant. More specifically,
we will introduce the auxiliary variable $\zeta$, representing the output
derivative, and the novel output $y_{e}$

$\left\\{\begin{array}[]{rcl}\underbrace{\begin{bmatrix}I&0\\\
C&0\end{bmatrix}}_{\texttt{Et}}\begin{bmatrix}\dot{x}\\\
\dot{\zeta}\end{bmatrix}(t)&=&\underbrace{\begin{bmatrix}A_{0}&0\\\
0&I\end{bmatrix}}_{\texttt{A0t}}\begin{bmatrix}x\\\
\zeta\end{bmatrix}(t)+\underbrace{\begin{bmatrix}A_{1}&0\\\
0&0\end{bmatrix}}_{\texttt{A1t}}\begin{bmatrix}x\\\
\zeta\end{bmatrix}(t-\tau_{1})+\underbrace{\begin{bmatrix}B\\\
0\end{bmatrix}}_{\texttt{Bt}}u(t-\tau_{2})\\\\[28.45274pt]
y_{e}(t)&=&\underbrace{\begin{bmatrix}C&0\\\
0&I\end{bmatrix}}_{\texttt{Ct}}\begin{bmatrix}x\\\
\zeta\end{bmatrix}(t),\end{array}\right.$ (3.24)

in which `Et`, `A0t`, `A1t`, `Bt`, and `Ct` correspond to the matrices in the
code below. A PID controller now corresponds to a structured dynamic output
feedback controller:

$\left\\{\begin{array}[]{rcl}\dot{x}_{c}(t)&=&0\,x_{c}(t)+\begin{bmatrix}I&0\end{bmatrix}y_{e}(t)\\\\[6.0pt]
u(t)&=&K_{I}x_{c}(t)+\begin{bmatrix}K_{P}&K_{D}\end{bmatrix}y_{e}(t),\end{array}\right.$

i.e., the matrix $A_{c}$ should remain fixed to the zero matrix and the matrix
$B_{c}$ should always equal $\begin{bmatrix}I&0\end{bmatrix}$. The entries of
the matrices $C_{c}$ and $D_{c}$ on the other hand can be freely tuned. To
design such a controller in we can use the following code.

⬇

% System matrices

A0 = [-3 2 6;-3 -4 1;1 -1 -2];

A1 = [0.5 -0.7 -0.8;-0.3 0.6 0.2;-0.2 0.6 -0.4];

B = [1 0;1 1;1 -1]; C = [1 1 1];

tau1 = 1; tau2 = 0.2;

%Define transformed system matrices

n = size(A0,1); p = size(B,2); q = size(C,1);

Et = [eye(n) zeros(n,q);C zeros(q)];

A0t = [A0 zeros(n,q);zeros(q,n) eye(q)];

A1t = [A1 zeros(n,q); zeros(q,n+q)];

Bt = [B; zeros(q,p)]; Ct = [C zeros(q);zeros(q,n) eye(q)];

% Use tds_create_ddae to create a representation for (3.24)

tds = tds_create_ddae(Et,{A0t A1t},[0 tau1],{Bt},tau2,{Ct});

% Define the mask and basis arguments

mask=ss(zeros(q),zeros(q,2*q),ones(p,q),ones(p,2*q));

basis=ss(zeros(q),[eye(q) zeros(q)],zeros(p,q),zeros(p,2*q));

% Design a structured controller

initial = ss(zeros(q),[eye(q) zeros(q)],ones(p,q),ones(p,2*q));

options_stabopt = tds_stabopt_options('nstart',1);

nc␣=␣q;

[K,cl]␣=␣tds_stabopt_dynamic(tds,nc,'mask',mask,'basis',basis,...

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣'initial',initial,'options',options_stabopt);

The resulting controller `K` is an `ss`-object with

⬇

K.A =

0

K.B =

1 0

K.C =

-5.1159

-2.5347

K.D =

-2.3480 -0.2611

1.6036 -0.2695

meaning that

$K_{P}=\begin{bmatrix}-2.3480\\\
\phantom{-}1.6036\end{bmatrix},\,K_{I}=\begin{bmatrix}-5.1159\\\
-2.5347\end{bmatrix},\text{ and }K_{D}=\begin{bmatrix}-0.2611\\\
-0.2695\end{bmatrix}.$ (3.25)

Using `tds_strong_sa`, we find that the strong spectral abscissa of the
corresponding closed loop is equal to -1.1563 meaning that the closed-loop
system is strongly exponentially stable.

∎

####  3.3.6 Networked systems and decentralized controllers

In the following example we will consider a controller design problem for a
networked system, i.e., a system that consists of multiple interacting
subsystem. In such a context, centralized control, in which one global
controller steers all subsystems — see Figure 3.3(a), is often impractical or
even infeasible, due to the high communication and processing requirements.
Therefore, decentralized control, in which each subsystem is controlled by
using only local information — see Figure 3.3(b), is often preferred as it
lowers the communication requirements. Furthermore, the local controllers can
process the available information in parallel. However, as each controller can
only use local information, this comes at the cost of a decrease in
performance. Distributed and overlapping controllers form a middle ground
between centralized and decentralized control. Each subsystem is still
controlled by a local controller, but these local controllers can now exchange
information with neighboring controllers or neighboring subsystems, see Figure
3.3(c). This allows for a trade-off between communication and processing
requirements on the one hand, and performance on the other hand.

$P_{1}$$P_{2}$$P_{3}$$P_{4}$$P_{5}$$P_{6}$$K$ (a) Centralized control

$P_{1}$$P_{2}$$P_{3}$$P_{4}$$P_{5}$$P_{6}$$K_{1}$$K_{2}$$K_{3}$$K_{4}$$K_{5}$$K_{6}$
(b) Decentralized control

$P_{1}$$P_{2}$$P_{3}$$P_{4}$$P_{5}$$P_{6}$$K_{1}$$K_{2}$$K_{3}$$K_{4}$$K_{5}$$K_{6}$
(c) Distributed/overlapping control

Figure 3.3: Different approaches for controlling networked systems.

######  Example 3.9 ([9, Section 5.1]).

Consider a networked system that consists of a ring of $N$ interconnected
subsystems (similar as in Figure 3.3). The dynamics for the $i$th subsystem
are given by

$\left\\{\begin{array}[]{rcl}\dot{x}_{i}(t)&=&\begin{bmatrix}0&0.5\\\
0.5&-3\end{bmatrix}x_{i}(t)+\begin{bmatrix}0\\\
5\end{bmatrix}u_{i}(t-0.1)+0.5x_{l,i}(t-0.3)+0.5x_{r,i}(t-0.3)\\\\[15.0pt]
y_{i}(t)&=&\begin{bmatrix}1&1\\\ 0&1\end{bmatrix}x_{i}(t),\end{array}\right.$
(3.26)

with $x_{l,i}(t)$ the state of its left neighbor and $x_{r,i}(t)$ the state of
its right neighbor. To create a representation for the overall system we use
the `kron` function (which returns the Kronecker tensor product of two
matrices).

⬇

% Number of subsystems

N = 5;

% Delays and system matrices of the individual subsystems

tau1 = 0.1; tau2 = 0.3;

Ai = [0 0.5;0.5 -3]; Adi = 0.5*eye(2);

Bi = [0;5]; Ci = [1 1;0 1];

% Adjacency matrix of the network

PN = diag(ones(1,N-1),-1) + diag(ones(1,N-1),1);

PN(1,N) = 1; PN(N,1) = 1;

% System matrices of the overall system

A0 = kron(eye(N),Ai); A1 = kron(PN,Adi);

B = kron(eye(N),Bi); C = kron(eye(N),Ci);

% Create open-loop system

P = tds_create({A0,A1},[0 tau2],{B},tau1,{C},0);

Next we will compare different controller structures for this system. We start
with designing a centralized static output feedback controller. As initial
controller, we use the values given in [9, Section 5.1].

⬇

%% STEP 2: Centralized controller

o = tds_stabopt_options('nstart',1);

Kinit␣=␣blkdiag([-0.8045␣0.6966],[0.8351␣-0.2437],...

␣␣␣␣␣␣␣␣␣␣[-1.1480␣0.1049],[0.7223␣2.5855],[0.2157␣-1.1658]);

[K1,cl1]␣=␣tds_stabopt_static(P,'options',o,'initial',Kinit);

We find the densely filled feedback matrix

$K=\begin{bmatrix}-4.010&\phantom{-}4.770&-6.227&\phantom{-}5.330&-1.081&0.276&2.486&-0.993&-11.033&\phantom{-}8.613\\\
-0.001&\phantom{-}1.657&-11.586&\phantom{-}9.225&-6.228&\phantom{-}6.175&\phantom{-}1.599&-0.682&-5.337&\phantom{-}3.453\\\
\phantom{-}5.321&-3.217&-6.622&\phantom{-}5.387&-7.855&\phantom{-}6.30&\phantom{-}0.323&\phantom{-}0.974&-7.508&\phantom{-}5.130\\\
\phantom{-}0.166&\phantom{-}0.532&-0.160&\phantom{-}0.085&-0.478&-0.055&-8.242&\phantom{-}7.473&-5.090&\phantom{-}4.458\\\
-1.129&\phantom{-}1.674&-0.242&\phantom{-}0.391&\phantom{-}3.320&-4.119&-2.007&\phantom{-}2.357&-9.736&\phantom{-}8.144\par\LTX@newpage\end{bmatrix}$

and the corresponding (strong) spectral abscissa is equal to -5.3382. Next, we
design a set of decentralized controllers by imposing a block diagonal
structure on the feedback matrix.

⬇

mask = kron(eye(N),[1 1]);

basis = zeros(N,N*2);

[K2,cl2] = tds_stabopt_static(P,'options',o,'initial',Kinit,...

␣␣␣␣␣␣␣␣␣'mask',mask,'basis',basis);

We obtain the feedback matrices

$\begin{array}[]{lll}K_{1}=\begin{bmatrix}-9.7969&9.0638\end{bmatrix},&K_{2}=\begin{bmatrix}-17.0255&15.4809\end{bmatrix},&K_{3}=\begin{bmatrix}-7.1284&6.0666\end{bmatrix},\\\\[7.02625pt]
K_{4}=\begin{bmatrix}-8.7850&7.7052\end{bmatrix},&K_{5}=\begin{bmatrix}-15.2369&13.5828\end{bmatrix},\end{array}$

and the (strong) spectral abscissa of the corresponding closed-loop system is
equal to -4.1872. As we have less free parameters in our feedback matrix, it
is not surprising that the damping rate decreases. The resulting controller is
however easier to implement. Finally, we design an overlapping controller in
which each local controller also receives the measurements of the neighboring
systems.

⬇

mask = kron(eye(N),[1 1]) + kron(PN,[1 1]);

basis = zeros(N,N*2);

[K3,cl3] = tds_stabopt_static(P,'options',o,'initial',Kinit,...

␣␣␣␣␣␣␣␣'mask',mask,'basis',basis);

We obtain the feedback matrix

$\begin{array}[]{lll}K=\begin{bmatrix}-6.511&5.532&-2.394&1.993&0&0&0&0&-7.581&7.459\\\
-1.459&1.377&-7.959&6.919&-0.892&0.899&0&0&0&0\\\
0&0&-10.15&10.12&-11.63&10.28&-4.757&4.991&0&0\\\
0&0&0&0&-4.103&3.760&-8.242&7.408&-5.409&5.102\\\
-2.404&2.791&0&0&0&0&1.354&-1.041&-9.455&7.691\end{bmatrix}\end{array}$

and the corresponding strong spectral abscissa is equal to -4.7371. ∎

###  3.4 Obtaining a state-space representation from a SISO transfer function

In this subsection we will briefly discuss how one can create a tds_ss-
representation for a time-delay system described by a SISO transfer function.
More specifically, the function tds_create_tf allows to represent transfer
functions of the following form

$T(s)=\frac{P(s)}{Q(s)}+D(s)$

with

$\textstyle
P(s)=\sum\limits_{k=1}^{m}\left(\sum\limits_{i=1}^{n_{P}}p_{k,i}\lambda^{n_{P}-(i-1)}\right)e^{-s\tau_{k}}\text{,
}Q(s)=\sum\limits_{k=1}^{m}\left(\sum\limits_{i=1}^{n_{Q}}q_{k,i}\lambda^{n_{Q}-(i-1)}\right)e^{-s\tau_{k}},$

and

$\textstyle D(s)=\sum\limits_{k=1}^{m}d_{k}e^{-\lambda\tau_{k}}$

in which $p_{k,i}$, $q_{k,i}$, and $d_{k}$ are real-valued coefficients and
$\tau_{1},\dots,\tau_{m}$ are (non-negative) delay values. To avoid advanced
systems, at least one delay must be equal to zero and the corresponding
leading coefficient of $Q(s)$ must be non-zero, i.e., let $\tau_{k}=0$ then
$q_{k,1}\neq 0$. Furthermore, the transfer function must be proper, i.e.,
$n_{P}\leq n_{Q}$. The quasi-polynomial $D(s)$ allows to specify direct
feedthrough terms if $\frac{P(s)}{Q(s)}$ is strictly proper, i.e.,
$n_{P}<n_{D}$.

The function tds_create_tf either takes three or four arguments:

  * •

`tds_create_tf(P,Q,TAU)` with `P` an $m\times(n_{P}+1)$-dimensional real-
valued matrix containing the coefficients $p_{k,i}$, `Q` an
$m\times(n_{Q}+1)$-dimensional real-valued matrix containing the coefficients
$q_{k,i}$, and `TAU` an $m$-dimensional vector containing the delays;

  * •

`tds_create_tf(P,Q,D,TAU)` with `P` an $m\times(n_{P}+1)$-dimensional real-
valued matrix containing the coefficients $p_{k,i}$, `Q` an
$m\times(n_{Q}+1)$-dimensional real-valued matrix containing the coefficients
$q_{k,i}$, `D` an $m$-dimensional real-valued vector containing the
coefficients $d_{k}$ and `TAU` an $m$-dimensional vector containing the
delays.

The return argument is a `tds_ss`-object that represents, depending on the
transfer function, a retarded, neutral or delay-descriptor system. This object
can subsequently be used as an input for the controller design functions
mentioned above. Let us now illustrate this function using two examples.

######  Example 3.10.

As a first example we will consider the following strictly proper transfer
function, corresponding to a retarded time-delay system

$T(s)=\frac{(2s-1)+3se^{-s}}{(4s^{2}-2s+1)+(3s-1)e^{-s}}.$

Extracting the coefficients of the quasi-polynomials and the delay values, we
get

⬇

P = [2 -1;3 0]; Q = [4 -2 1;0 3 -1]; TAU = [0 1];

sys = tds_create_tf(P,Q,TAU)

with sys a tds_ss_retarded-object with inputs and outputs:

⬇

sys =

LTI Retarded Time-Delay State-Space System with properties:

A: {[2x2 double] [2x2 double]}

hA: [0 1]

B1: {[2x1 double] [2x1 double]}

hB1: [0 1]

C1: {[0 1]}

hC1: 0

in which some parts of the output were omitted for notational convenience.

Similiarly we can create a representation for the following neutral system
with explicit direct feedthrough terms

$T(s)=\frac{(6s-2)e^{-s}+se^{-2s}+(s-1)}{(5s^{2}-3s-2)e^{-2s}+(-2s^{2}+3s-2)}-2e^{-s}+3e^{-2s}+1$

using the code

⬇

P = [6 -2;1 0;1 -1]; Q = [0 0 0;5 -3 -2;-2 3 -2];

D = [-2; 3; 1]; TAU = [1 2 0];

sys = tds_create_tf(P,Q,D,TAU);

which returns a tds_ss_neutral-object

⬇

sys =

LTI Neutral Time-Delay State-Space System with properties:

H: {[2x2 double]}

hH: 2

A: {[2x2 double] [2x2 double]}

hA: [0 2]

B1: {[2x1 double] [2x1 double] [2x1 double]}

hB1: [0 1 2]

C1: {[0 1]}

hC1: 0

D11: {[1] [-2] [3]}

hD11: [0 1 2]

∎

##  Chapter 4 Performance criteria and robust control

In this chapter, we will discuss how the $\mathcal{H}_{\infty}$-norm, the
pseudospectral abscissa and the distance to instability of a LTI time-delay
system can be computed within . Besides, we will show how a static or dynamic
output feedback controller that (locally) minimizes the closed-loop
$\mathcal{H}_{\infty}$-norm can be designed. First, Section 4.1 will introduce
the $\mathcal{H}_{\infty}$-norm for time-delay systems. As was the case for
the spectral abscissa of a neutral DDE, the $\mathcal{H}_{\infty}$-norm of a
time-delay system that is not strictly proper, might suffer from a fragility
problem, in the sense that the $\mathcal{H}_{\infty}$-norm might be
discontinuous with respect to the delay values. will therefore consider the
strong $\mathcal{H}_{\infty}$-norm (introduced in [16, Definition 4.4]), which
is a continuous function of both the entries of the system matrices and the
delays. Besides we also will discuss how the $\mathcal{H}_{\infty}$-norm can
be used for examining the robust stability of feedback systems in the presence
of uncertainties and for characterizing input-to-output performance. Next,
Section 4.2 will consider several $\mathcal{H}_{\infty}{}$-optimization
problems arising in model order reduction and the robust control framework.
For more details on robust control theory we refer the interested reader to
[43] or the derived textbook [44], which discuss various aspects of robust
control theory focusing on multiple-input-multiple-output systems without
delays. For specific results related to delay systems, we refer to [29] and
[42]. Finally, Section 4.3 introduces the pseudospectral abscissa and the
distance to instability which allow to examine the robust stability of a
dynamical system in the presence of structured uncertainties [19].

###  4.1 Computing the H-infinity norm

In this subsection we will discuss how the $\mathcal{H}_{\infty}$-norm of a
(stable) LTI time-delay system can be computed in . More specifically, we will
consider systems of the following form

$\left\\{\begin{array}[]{rcl}\displaystyle
E\dot{x}(t)&=&\displaystyle\sum\limits_{k=1}^{m_{A}}A_{k}\,x(t-h_{A,k})+\sum_{k=1}^{m_{B_{2}}}B_{2,k}\,w(t-h_{B_{2},k})-\sum_{k=1}^{m_{H}}H_{k}\dot{x}(t-h_{H,k})\\\\[13.0pt]
z(t)&=&\displaystyle\sum\limits_{k=1}^{m_{C_{2}}}C_{2,k}\,x(t-h_{C_{2},k})+\sum_{k=1}^{m_{D_{22}}}D_{22,k}\,w(t-h_{D,k})\end{array}\right.$
(4.1)

with $x\in\mathbb{R}^{n}$ the state variable, $w\in\mathbb{R}^{p_{2}}$ the
performance input (representing e.g., external disturbances) and
$z\in\mathbb{R}^{q_{2}}$ the performance output (representing e.g., error
signals). As in (3.1), all matrices are real-valued and of appropriate
dimensions; all delays are non-negative (except for $h_{H,1}$,…,$h_{H,H}$
which should be positive); and the matrix $E$ can be singular.

Note: In contrast to (3.1) the input and output channels are now denoted by
$w$ and $z$ instead of $u$ and $y$. In (3.1) the input and output channels
($u$ and $y$) were used for feedback control, while here the input and output
channels ($w$ and $z$) are used to quantify the input-to-output behavior of
the system. As these input and output channel serve a difference purpose, we
will use a different symbol to denote them. In the next section (Section 4.2),
where we will design output feedback controllers that (locally) minimize the
$\mathcal{H}_{\infty}$-norm, we will need both feedback channels, denoted by
$u$ and $y$, and performance channels, denoted $w$ and $z$. Accordingly, we
use $B_{2,k}$, $C_{2,k}$, and $D_{22,k}$ instead of $B_{1,k}$, $C_{1,k}$, and
$D_{11,k}$ to denote the input, output and direct feedthrough matrices. This
change in notation has also consequences on the syntax that needs to be used
when constructing a `tds_ss`-object for the system, as spelled out below.

We can again use the functions tds_create, tds_create_neutral and
tds_create_ddae to represent (4.1) in . However, we now need to explicitly
indicate that we want to specify the performance input, output and direct
feedthrough terms.

  * •

tds_create(A,hA,’B2’,B2,hB2,’C2’,C2,hC2,’D22’,D22,hD22) allows to create
retarded state-space models, i.e., $E=I_{n}$ and $m_{H}=0$.

  * •

tds_create_neutral(H,hH,A,hA,’B2’,B2,hB2,’C2’,C2,hC2,’D22’,D22,hD22) allows to
create neutral state-space models, i.e.,  $E=I_{n}$ and $m_{H}>0$.

  * •

tds_create_ddae(E,A,hA,’B2’,B2,hB2,’C2’,C2,hC2,’D22’,D22,hD22) allows to
create delay descriptor systems, i.e., $E$ not necessarily identity and
potentially singular, and $m_{H}=0$. Important: Recall that all DDAEs and
delay descriptor systems considered in should satisfy 2.1.

Notice that we now have to use the flags `’B2’`, `’C2’`, and `’D22’` to
indicate that we want to set the fields related to the inputs $w$ and outputs
$z$. In all function calls above, the arguments `D22` and `hD22` are again
optional. Further, the arguments `H`, `A`, `B2`, `C2`, and `D22` should be
cell arrays containing the system matrices. The arguments `hH`, `hA`, `hB2`,
`hC2`, and `hD22` should be arrays containing the delay values. Finally, the
argument `E` of tds_create_ddae should be a (real-valued) matrix.

For the $\mathcal{H}_{\infty}$-norm computation, will internally transform all
representable systems to the following standard form

$\left\\{\begin{array}[]{rcl}\displaystyle
E_{s}\dot{x}_{s}(t)&=&A_{s,0}\,x_{s}(t)+\displaystyle\sum\limits_{k=1}^{m}A_{s,k}\,x_{s}(t-\tau_{k})+B_{s}w(t)\\\\[13.0pt]
z(t)&=&\displaystyle C_{s}x_{s}(t).\end{array}\right.$ (4.2)

with $\vec{\tau}=[\tau_{1},\dots,\tau_{m}]$ the union of all delays and
$x_{s}$ the augmented state variable to which additional states are added to
eliminate the input and output delays, the direct feedthrough terms and the
delayed state-derivatives (i.e., the “neutral” terms). For ease of notation,
we will use this standard form in the following discussion, but note that the
software can deal with all systems that can be created using tds_create,
tds_create_neutral and tds_create_ddae.

Let us start with defining the $\mathcal{H}_{\infty}$-norm in terms of the
standard form (4.2). To this end, let us first introduce the corresponding
transfer function from input $w$ to output $z$, which is given by

$T_{zw}(s;\vec{\tau})=C_{s}\left(sE_{s}-A_{s,0}-\sum\limits_{k=1}^{m}A_{s,k}\,e^{-s\tau_{k}}\right)^{-1}B_{s}.$
(4.3)

This transfer function describes the input-to-output behavior of the system
for a zero initial state. Furthermore, if the null solution of the following
DDAE

$\displaystyle
E_{s}\dot{x}_{s}(t)=A_{s,0}\,x_{s}(t)+\displaystyle\sum\limits_{k=1}^{m}A_{s,k}\,x_{s}(t-\tau_{k})$

is exponentially stable, the $\mathcal{H}_{\infty}$-norm of (4.2) is defined
as:

$\|T_{zw}(s;\vec{\tau})\|_{\mathcal{H}_{\infty}}:=\sup_{\begin{subarray}{c}s\in\mathbb{C}\\\
\Re(s)>0\end{subarray}}\|T_{zw}(s;\vec{\tau})\|_{2}=\sup_{\omega\in\mathbb{R}^{+}}\|T_{zw}(\jmath\omega;\vec{\tau})\|_{2},$
(4.4)

in which the last equality follows from an extension of the maximum modulus
principle. The $\mathcal{H}_{\infty}$-norm of a stable system is thus equal to
the maximal value of the frequency response (i.e., the transfer function
evaluated on the imaginary axis) measured in spectral norm. On the other hand,
we will say that the $\mathcal{H}_{\infty}$-norm of an unstable system is
equal to $\infty$.

The $\mathcal{H}_{\infty}$-norm has several applications in (robust) control
theory.

  * •

Firstly, the $\mathcal{H}_{\infty}$-norm can be used to measure the
performance of a dynamical system. More precisely, it quantifies the ‘worst-
case’ suppression (or amplification) from the input signal $w$ to the output
$z$ (assuming zero initial conditions), in the sense that

$\|T_{zw}(s;\vec{\tau})\|_{\mathcal{H}_{\infty}}=\limsup_{w\in
L_{2}}\frac{\|z\|_{L_{2}}}{\|w\|_{L_{2}}},$

with $L_{2}$ the (Hilbert) space of square-integratable functions on the
interval $[0,+\infty)$ and $\|\cdot\|_{L_{2}}$ the associated norm, i.e.,
$f\in L_{2}$ if
$\|f\|_{L_{2}}:=\left(\int_{0}^{+\infty}\|f(t)\|^{2}_{2}\,\mathop{}\\!\mathrm{d}t\right)^{1/2}<\infty$.
For more details see [43, Chapter 4]. By appropriately choosing the
performance inputs and outputs, the $\mathcal{H}_{\infty}$-norm can thus be
used to quantify the reference tracking, the disturbance rejection, and the
noise suppression capabilities of the system. We will further examine this
interpretation in Examples 4.5 and 4.6.

Note. By viewing the system (4.2) as a (linear) operator mapping input signals
in $L_{2}$ to output signals in $L_{2}$, the $\mathcal{H}_{\infty}$-norm can
thus be interpreted as an induced norm.

  * •

Secondly, the $\mathcal{H}_{\infty}$-norm allows to quantify the effect of
uncertainties on the stability of a dynamic system. More specifically, the
$\mathcal{H}_{\infty}$-norm is equal to the reciprocal of the structured
complex distance to instability of an exponentially stable system, i.e.,

$\|T_{zw}\|_{\mathcal{H}_{\infty}}=\inf\\{\|\Delta\|_{2}:\Delta\in\mathbb{C}^{p_{2}\times
q_{2}}\text{ such that }\\\
E_{s}\dot{x}_{s}(t)=\sum\limits_{k=0}^{m}A_{s,k}\,x_{s}(t-\tau_{k})+B_{s}\Delta
C_{s}\,x_{s}(t)\text{ is not exponentially stable }\\}^{-1},$

with $1/0=+\infty$, to be consistent with the definition of the
$\mathcal{H}_{\infty}$-norm given above.

The $\mathcal{H}_{\infty}$-norm is thus a measure for the robust stability of
the system. Note however that the $\mathcal{H}_{\infty}$-norm is only a rough
estimate for the distance to instability of the system in a practical sense as
it only considers one complex (matrix-valued) uncertainty. In Section 4.3, we
will see how the structured distance to instability with respect to real-
valued uncertainties acting on both the system matrices and the delays can be
computed in , which is a better measure for robust stability in practice.

  * •

Thirdly, the $\mathcal{H}_{\infty}$-norm can be used as an error measure for a
reduced-order approximation of a dynamical system [43, Chapter 7]. For
example, consider the following exponentially stable time-delay system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&A_{0}\,x(t)+\sum\limits_{k=1}^{m}A_{k}\,x(t-\tau_{k})+B_{2}\,w(t)\\\
z(t)&=&C_{2}\,x(t)\end{array}\right.$

and the following delay-free approximation

$\left\\{\begin{array}[]{rcl}\dot{x}_{a}(t)&=&A_{a}\,x_{a}(t)+B_{a}\,w(t)\\\
z(t)&=&C_{a}\,x_{a}(t)\end{array}\right.$ (4.5)

The closeness of these two models can then be quantified by

$\|T(s)-T_{\mathrm{approx.}}(s)\|_{\mathcal{H}_{\infty}},$ (4.6)

with $T(s)$ and $T_{\mathrm{approx.}}(s)$ the transfer function of the
original system and the approximation, respectively. Note that (4.6)
corresponds to the maximal mismatch of the original frequency response and its
approximation in spectral norm. We will consider this application in more
detail in Example 4.7. Furthermore, in the same example we will also see that
finding the state-space system of the form (4.5) that minimizes the mismatch
in terms of (4.6) can be reformulated as a $\mathcal{H}_{\infty}{}$-controller
design problem (which will be considered in Section 4.2). As such, can also be
used to derive a delay-free approximation for a time-delay system.

  * •

Finally, in robust control theory the small gain theorem is frequently used to
guarantee stability in the presence of uncertainties [43, Chapter 9].
Application of this small gain theorem typically gives rise to constraints in
terms of the $\mathcal{H}_{\infty}$-norm of an appropriately defined transfer
function.

To compute the $\mathcal{H}_{\infty}$-norm of a time-delay system, provides
the function tds_hinfnorm. This function takes one mandatory argument, the
tds_ss object of which we want to compute the $\mathcal{H}_{\infty}$-norm.
Additional options can be specified using the optional argument `options`. For
more details on the available options, see help tds_hinfnorm or Section B.8.

######  Remark 4.1.

The function tds_hinfnorm employs a combination of the BBBS-algorithm [6, 4]
and Newton’s method to compute the $\mathcal{H}_{\infty}$-norm of a time-delay
system. More specifically, in a first step it uses the BBBS-algorithm to
compute the $\mathcal{H}_{\infty}$-norm of a finite dimensional approximation
of the system (obtained using a similar spectral discretisation approach as
was used in (2.9) to approximate the characteristic roots of a DDE). This
result is subsequently refined using Newton corrections. The BBBS-algorithm
itself consists of two step that are alternated iteratively. First, it
computes the frequency intervals at which the spectral norm of the frequency
response is larger than a given level. Next, it updates this level with the
maximum value of the spectral norm of the frequency response at the geometric
mean of these intervals. This two-step process is repeated until the change in
the level is sufficiently small.

######  Example 4.1.

Consider the following single-input single-output (SISO) time-delay system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&\begin{bmatrix}-4&2\\\
-3&-3\end{bmatrix}x(t)+\begin{bmatrix}-2&-1\\\
3&-2\end{bmatrix}x(t-1)+\begin{bmatrix}1\\\
-1\end{bmatrix}w(t-2),\\\\[10.03749pt]
z(t)&=&\begin{bmatrix}-2&1\end{bmatrix}x(t).\end{array}\right.$ (4.7)

As the underlying RDDE is exponentially stable, the asymptotic input-to-output
behavior of this system can be characterized in terms of the following
transfer function:

$T_{zw}(s)=-\frac{9+3s+5e^{-s\tau_{1}}}{s^{2}+7s+18+(4s+5)e^{-s\tau_{1}}+7e^{-s2\tau_{1}}}e^{-s\tau_{2}}.$
(4.8)

The corresponding $\mathcal{H}_{\infty}$-norm is found using the following
code.

⬇

% Create a representation for the state-space system

A0 = [-4 2;-3 -3]; A1 = [-2 -1;3 -2];

B = [1;-1]; C = [-2 1];

tds = tds_create({A0 A1},[0 1],'B2',{B},[2],'C2',{C},[0]);

%␣Compute␣the␣H-infinity␣norm

[hinf,wPeak]␣=␣tds_hinfnorm(tds);

The function `tds_hinfnorm` returns both the $\mathcal{H}_{\infty}$-norm
(`hinf`) and the angular frequency at which this maximal gain is achieved
(`wPeak`), i.e., the $\omega$ that maximizes
$\|T_{zw}(\jmath\omega;\vec{\tau})\|_{2}$ (note that `wPeak` might be equal to
$\infty$ for transfer functions that are not strictly proper, meaning that the
supremum in (4.4) is only attained in the limit $\omega\rightarrow\infty$). In
this case we find that the $\mathcal{H}_{\infty}$-norm is equal to 1.5388 and
is attained at $\omega=3.5571$.

We can use the function tds_sigma(tds,omega) to verify the obtained result.
This function computes the singular values of the transfer function associated
with tds at the given angular frequencies omega. For a SISO-system, this plot
corresponds to the magnitude of the frequency response.

⬇

% Plot the magnitude of the frequency response

omega = logspace(-2,2,10001);

[sv] = tds_sigma(tds,omega);

semilogx(omega,sv)

hold on

plot(xlim,[hinf hinf])

plot(wPeak,hinf,'x')

Figure 4.1 shows the result and indicates the values `hinf` and `wPeak`
obtained above. (Luckily,) we see that the obtained result is indeed correct.

Figure 4.1: The singular value plot of (4.8) for
$s\in\jmath\,[10^{-2},\,10^{2}]$ (black), the value of the
$\mathcal{H}_{\infty}$-norm (red) and the angular frequency at which the
maximal gain is attained (cross).

∎

As mentioned before, the $\mathcal{H}_{\infty}$-norm might be discontinuous
with respect to the delays. This observation implies that an infinitesimal
change on the delays might cause a jump in the value of the
$\mathcal{H}_{\infty}$-norm, as illustrated in the example below.

######  Example 4.2.

Consider the following single-input single-output (SISO) system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&\begin{bmatrix}-4&2\\\
-3&-3\end{bmatrix}x(t)+\begin{bmatrix}-2&-1\\\
3&-2\end{bmatrix}x(t-\tau_{1})+\begin{bmatrix}1\\\
-1\end{bmatrix}w(t-\tau_{2})\\\\[13.04874pt]
z(t)&=&\begin{bmatrix}-2&1\end{bmatrix}x(t)+w(t)+w(t-\tau_{1})-2w(t-\tau_{2})\end{array}\right.$
(4.9)

with $\tau_{1}=1$ and $\tau_{2}=2$. Notice that the considered system is quite
similar to the one in Example 4.1. The only difference is that we have added
three direct feedthrough terms. The associated transfer function is now given
by:

$T_{zw}(s)=-\frac{9+3s+5e^{-s\tau_{1}}}{s^{2}+7s+18+(4s+5)e^{-s\tau_{1}}+7e^{-s2\tau_{1}}}e^{-s\tau_{2}}+1+e^{-s\tau_{1}}-2e^{-s\tau_{2}},$
(4.10)

which, in contrast to (4.8), is not strictly proper. Furthermore, it can be
shown that the $\mathcal{H}_{\infty}$-norm of (4.9) is equal to 3.7019
(rounded to 4 decimal places). However, if we make a small perturbation to
$\tau_{1}$, e.g., $\tau_{1}=0.995$, the $\mathcal{H}_{\infty}$-norm suddenly
jumps to approximately 4. Furthermore, even if we reduce the size of the
perturbation, e.g., $\tau_{1}=0.999$, we still observe a jump in the
$\mathcal{H}_{\infty}$-norm to approximately 4. To understand what happens,
let us use the tds_sigma-function to plot the magnitude of the frequency
response.

⬇

% Create a representation for the time-delay system

A0 = [-4 2;-3 -3]; A1 = [-2 1;3 -2];

B = [1;-1]; C = [-2 1];

D0 = 1; D1 = 1; D2 = -2;

tau1 = 1; tau2 = 2;

tds = tds_create({A0 A1},[0 tau1],'B2',{B},tau2,'C2',{C},0,...

␣␣␣␣␣␣␣'D22',{D0,D1,D2},[0␣tau1␣tau2]);

%␣Plot␣the␣magnitude␣of␣the␣frequency␣response

omega␣=␣logspace(-2,3.5,10001);

[sv]␣=␣tds_sigma(tds,omega);

figure

semilogx(omega,sv)

From Figure 4.2(a), we observe that the frequency responce indeed attains a
maximum of approximately $3.7019$ (at $\omega=2.7509$). Next, we repeat the
same procedure for $\tau_{1}=0.995$.

⬇

% Create a representation for the time-delay system

tau1 = 0.995;

tds2 = tds_create({A0 A1},[0 tau1],'B2',{B},tau2,'C2',{C},0,...

␣␣␣␣␣␣␣␣'D22',{D0,D1,D2},[0␣tau1␣tau2]);

%␣Plot␣the␣magnitude␣of␣the␣frequency␣responce

[sv2]␣=␣tds_sigma(tds2,omega);

figure

semilogx(omega,sv2)

and $\tau_{1}=0.999$

⬇

tau1 = 0.999;

tds3 = tds_create({A0 A1},[0 tau1],'B2',{B},tau2,'C2',{C},0,...

␣␣␣␣␣␣␣␣'D22',{D0,D1,D2},[0␣tau1␣tau2]);

%␣Plot␣the␣magnitude␣of␣the␣frequency␣responce

[sv3]␣=␣tds_sigma(tds3,omega);

figure

semilogx(omega,sv3)

Figures 4.2(b) and 4.2(c) show the result. We observe that in both cases the
frequency response at high frequencies significantly differs from the one in
Figure 4.2(a) and that the $\mathcal{H}_{\infty}$-norm in both cases jumps to
approximately 4.

(a) $\tau_{1}=1$ and $\tau_{2}=2$.

(b) $\tau_{1}=0.995$ and $\tau_{2}=2$.

(c) $\tau_{1}=0.999$ and $\tau_{2}=2$.

Figure 4.2: The (maximum) singular value plot of (4.10) for several values of
$\tau_{1}$ and $\tau_{2}$.

∎

We can draw several parallels with Example 2.6 and Remark 2.2. Firstly, as the
spectral abscissa of a NDDE, the $\mathcal{H}_{\infty}$-norm of a time-delay
system that is not strictly proper might be a discontinuous function of the
delay values, even if the system is strongly exponentially stable. However, if
we look at a fixed angular frequency, the singular values of the frequency
response change in a continuous way when the delay parameters are varied.
(Recall that we made a similar observation for the individual characteristic
roots of a NDDE.) Nonetheless, we observe from Figures 4.2(b) and 4.2(c) that
the larger the angular frequency, the larger the sensitivity (i.e.,
derivative) of the singular values of the frequency response to delay
perturbations. In the example above, this derivative grows unbounded as
$\omega$ goes to infinity. Furthermore, it can be shown that for any
$\epsilon>0$ and $\delta>0$ (no matter how small), there exists a delay
perturbation smaller than $\epsilon$ for which the $\mathcal{H}_{\infty}$-norm
of (4.9) jumps to at least $4-\delta$. Furthermore as this perturbation bound
is decreased, the frequencies at which such a peak in the singular value plot
can be attained, moves to infinity (as indicated with the arrow in Figure
4.2(c)). Note again the similarity with Remark 2.2. Below, we will therefore
introduce the strong $\mathcal{H}_{\infty}$-norm and the asymptotic transfer
function, which have similar roles as the strong spectral abscissa and the
delay difference equation underlying the NDDE in Section 2.3.

The asymptotic transfer function of (4.2) is given by:

$T_{zw,a}(s):=-C_{s}^{(2)}\left(A_{0}^{(22)}+\sum_{k=1}^{m}A_{k}^{(22)}e^{-s\tau_{k}}\right)^{-1}B_{s}^{(2)},$
(4.11)

with

$A_{k}^{(22)}=U^{T}A_{k}V\text{, }C_{s}^{(2)}=C_{s}V\text{ and
}B_{s}^{(2)}=U^{T}B_{s}$

in which the columns of the matrices $U$ and $V$ span an orthonormal basis for
the left and right null space of $E$, respectively, see Section 2.4.

This asymptotic transfer function owes its name to the fact that it converges
to (4.3) at high frequencies (see [16, Equations (3.2)-(3.3) and Proposition
3.3].

######  Remark 4.2.

For time-delay systems of the form

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&\displaystyle\sum_{k=1}^{m_{A}}A_{k}\,x(t-h_{A,k})+\sum_{k=1}^{m_{B_{2}}}B_{2,k}\,w(t-h_{B_{2},k})-\sum_{k=1}^{m_{H}}H_{k}\,\dot{x}(t-h_{H,k})\\\\[12.045pt]
z(t)&=&\displaystyle\sum_{k=1}^{m_{C_{2}}}C_{2,k}\,x(t-h_{C_{2},k})+\sum_{k=1}^{m_{D_{22}}}D_{22,k}\,w(t-h_{D_{22},k}),\end{array}\right.$
(4.12)

the asymptotic transfer function is given by

$T_{zw,a}(s)=\sum_{k=1}^{m_{D_{22}}}D_{22,k}e^{-sh_{D_{22},k}}.$

For example, the asymptotic transfer function associated with (4.9) is equal
to

$T_{zw,a}(s)=1+e^{-s\tau_{1}}-2e^{-s\tau_{2}}.$ (4.13)

Next, we define the strong $\mathcal{H}_{\infty}$-norm of (4.3) as the
smallest upper bound for the $\mathcal{H}_{\infty}$-norm which insensitive to
(sufficiently) small delay changes.

######  Definition 4.1 ([16, Definitions 4.2 and 4.4]).

The strong $\mathcal{H}_{\infty}$-norm of (4.3) is given by

$|||T_{zw}(s;\vec{\tau})|||_{\mathcal{H}_{\infty}}:=\lim_{\epsilon\searrow
0+}\sup\\{\|T_{zw}(\cdot;\vec{\tau}_{\epsilon})\|_{\mathcal{H}_{\infty}}:\vec{\tau}_{\epsilon}\in\mathcal{B}(\vec{\tau},\epsilon)\cap\mathbb{R}_{+}^{m}\\}$

with $\mathcal{B}(\vec{\tau},\epsilon)$ as defined in Definition 2.7 and
$T_{zw}(\cdot;\vec{\tau}_{\epsilon})$ the transfer function in (4.3) where the
delays are replaced by $\vec{\tau}_{\epsilon}$. Similarly, the strong
$\mathcal{H}_{\infty}$-norm of (4.11) is given by

$|||T_{zw,a}(s;\vec{\tau})|||_{\mathcal{H}_{\infty}}:=\lim_{\epsilon\searrow
0+}\sup\\{\|T_{zw,a}(\cdot;\vec{\tau}_{\epsilon})\|_{\mathcal{H}_{\infty}}:\vec{\tau}_{\epsilon}\in\mathcal{B}(\vec{\tau},\epsilon)\cap\mathbb{R}_{+}^{m}\\}.$

In contrast to the nominal $\mathcal{H}_{\infty}$-norm, the strong
$\mathcal{H}_{\infty}$-norm continuously depends on both the entries of the
matrices in and the delay values [16, Thoerem 4.5].

Next, we will give a characterization for the strong
$\mathcal{H}_{\infty}$-norm of both (4.3) and (4.11) that is useful from a
computational point of view.

######  Proposition 4.1.

Under the assumption of strong exponential stability, the strong
$\mathcal{H}_{\infty}$-norm of the asymptotic tranfer function (4.11) is given
by

$|||T_{zw,a}|||_{\mathcal{H}_{\infty}}=\max_{\vec{\theta}\in[0,2\pi)^{m}}\left\|C_{s}^{(2)}\left(A_{0}^{(22)}+\sum_{k=1}^{m}A_{k}^{(22)}e^{-\jmath\theta_{k}}\right)^{-1}B_{s}^{(2)}\right\|_{2}.$
(4.14)

For state-space systems of the form (4.12) this expression reduces to,

$|||T_{zw,a}(s)|||_{\mathcal{H}_{\infty}}=\max_{\vec{\theta}\in[0,2\pi)^{m_{D_{22}}}}\left\|\sum_{k=1}^{m_{D_{22}}}D_{22,k}e^{-\jmath\theta_{k}}\right\|_{2}.$

The strong $\mathcal{H}_{\infty}$-norm of (4.3) is given by

$|||T_{zw}||_{\mathcal{H}_{\infty}}=\max\left\\{\|T_{zw}\|_{\mathcal{H}_{\infty}},|||T_{zw,a}|||_{\mathcal{H}_{\infty}}\right\\}.$

The strong $\mathcal{H}_{\infty}$-norm of (4.3) can thus be computed using a
two-step approach. First the strong $\mathcal{H}_{\infty}$-norm of the
associated asymptotic transfer function is computed by solving the
optimization problem in (4.14). Next, it is verified whether the spectral norm
of the frequency response attains a higher value at a finite frequency using
the approach mentioned above Example 4.1. This procedure is described in more
detail in [16].

As small delay perturbations are inevitable in practice, the function
`tds_hinfnorm` will always return the strong $\mathcal{H}_{\infty}$-norm of a
time-delay system. In the following example, we apply this function to the
system in (4.9).

######  Example 4.3.

Let us reconsider (4.9) from Example 4.2 with $\tau_{1}=1$ and $\tau_{2}=2$.
The corresponding asymptotic transfer function is given by (4.13). It follows
from Proposition 4.1 that

$|||T_{zw,a}|||_{\mathcal{H}_{\infty}}=\max_{\vec{\theta}\in[0,2\pi)^{2}}\|1+e^{-\jmath\theta_{1}}-2e^{-\jmath\theta_{2}}\|_{2}=4.$

Figure 4.2(a) shows that the spectral norm of the frequency response does not
attain a higher value, meaning that the strong $\mathcal{H}_{\infty}$-norm of
(4.9) is equal to 4. Let us verify this result using the
`tds_hinfnorm`-function.

⬇

>> [hinf,wPeak] = tds_hinfnorm(tds)

hinf =

4

wPeak =

Inf

Notice that wPeak now equals $\infty$, indicating that strong
$\mathcal{H}_{\infty}$-norm of (4.10) is equal to the strong
$\mathcal{H}_{\infty}$-norm of the associated asymptotic transfer function in
(4.13). (Note however that it does not imply that
$\lim_{\omega\rightarrow\infty}\sigma_{1}\left(T_{zw}(\jmath\omega\right)$
exists.) ∎

######  Remark 4.3.

Similar to the strong spectral abscissa of the associated delay-difference
equation ($C_{D}$), the strong $\mathcal{H}_{\infty}$-norm of the asymptotic
transfer function can be computed separately in . To this end, we first obtain
a representation for the asymptotic transfer function using the function
get_asymptotic_transfer_function. Next, we can apply tds_hinfnorm on the
resulting object to compute its strong $\mathcal{H}_{\infty}$-norm.

⬇

Ta = get_asymptotic_transfer_function(tds);

[hinfa,wPeaka] = tds_hinfnorm(Ta)

hinfa =

4

wPeaka =

Inf

In this section, we introduced the $\mathcal{H}_{\infty}$-norm and the strong
$\mathcal{H}_{\infty}$-norm and briefly discussed their applicability. Next,
we will discuss how can be used for $\mathcal{H}_{\infty}$\- synthesis
problems.

###  4.2 Designing controllers that minimize the strong
$\mathcal{H}_{\infty}$-norm

This subsection describes how can be used to design output feedback
controllers that (locally) minimize the strong $\mathcal{H}_{\infty}$-norm of
the corresponding closed-loop system. To this end, we will consider the
following state-space representation

$\left\\{\begin{array}[]{rcc}\displaystyle
E\dot{x}(t)&=&\displaystyle\sum\limits_{k=1}^{m_{A}}A_{k}\,x(t-h_{A,k})+\sum_{k=1}^{m_{B_{1}}}B_{1,k}\,u(t-h_{B_{1},k})+\sum_{k=1}^{m_{B_{2}}}B_{2,k}\,w(t-h_{B_{2},k})\\\\[13.0pt]
&&\displaystyle-\sum_{k=1}^{m_{H}}H_{k}\dot{x}(t-h_{H,k})\\\\[15.0pt]
y(t)&=&\displaystyle\sum\limits_{k=1}^{m_{C_{1}}}C_{1,k}\,x(t-h_{C_{1},k})+\sum_{k=1}^{m_{D_{11}}}D_{11,k}\,u(t-h_{D_{11},k})+\sum_{k=1}^{m_{D_{12}}}D_{12,k}\,w(t-h_{D_{12},k})\\\\[15.0pt]
z(t)&=&\displaystyle\sum\limits_{k=1}^{m_{C_{2}}}C_{2,k}\,x(t-h_{C_{2},k})+\sum_{k=1}^{m_{D_{21}}}D_{21,k}\,u(t-h_{D_{21},k})+\sum_{k=1}^{m_{D_{22}}}D_{22,k}\,w(t-h_{D_{22},k})\end{array}\right.$
(4.15)

with $x\in\mathbb{R}^{n}$ the state variable, $u\in\mathbb{R}^{p_{1}}$ the
control input (related to the “actuators” that steer the system),
$y\in\mathbb{R}^{q_{1}}$ the measured output (representing the available
measurements), $w\in\mathbb{R}^{p_{2}}$ the performance input (representing
e.g., noise signals), and $z\in\mathbb{R}^{q_{2}}$ the performance output
(representing e.g., error signals). As always, the system matrices are real-
valued and have appropriate dimensions, and the delays are non-negative
(positive). To create a representation for a time-delay system of the form
(4.15), we again use the functions `tds_create`, `tds_create_neutral` and
`tds_create_ddae`.

  * •

tds_create(A,hA,B1,hB1,C1,hC1,D11,hD11,B2,hB2,C2,hC2,D12,hD12,D21,...
hD21,D22,hD22) allows to create retarded state-space models, i.e., $E=I_{n}$
and $m_{H}=0$.

  * •

tds_create_neutral(H,hH,A,hA,B1,hB1,C1,hC1,D11,hD11,B2,hB2,C2,hC2,...D12,hD12,D21,hD21,D22,hD22)
allows to create neutral state-space models, i.e., $E=I_{n}$ and $m_{H}>0$.

  * •

tds_create_ddae(E,A,hA,B1,hB1,C1,hC1,D11,hD11,B2,hB2,C2,hC2,D12,...
hD12,D21,hD21,D22,hD22) allows to create delay descriptor systems, i.e., $E$
is not necessarily the identity matrix can can be singular and $m_{H}=0$.

The arguments H, A, B1, C1, D11, B2, C2, D12, D21, and D22 should be cell
arrays containing the system matrices. The arguments hH, hA, hB1, hC1, hD11,
hB2, hC2, hD12,hD21, and hD22 should be arrays containing the delay values.
The argument E of tds_create_ddae should be a (real-valued) matrix. As often
only a limited number of fields are present, it is also possible to specify
the necessary fields using key-value pairs. Furthermore, when using the key-
value based syntax, zero delays do not need to be specified. For example, to
create a representation for the retarded state-space system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&A_{1}x(t-h_{A,1})+A_{2}x(t-h_{A,2})+B_{1}u(t)+B_{2}w(t)\\\
y(t)&=&C_{1}x(t)+D_{12,1}w(t-h_{D_{12,1}})+D_{12,2}w(t-h_{D_{12,2}})\\\
z(t)&=&C_{2}x(t)\end{array}\right.$

the following syntax can be used

⬇

tds_create(A,hA,B1,[0],C1,[0],'B2',B2,'C2',C2,'D12',D12,hD12)

with A, B1, C1, B2, C2, and D12 cell arrays and hA and hD12 arrays.

The considered $\mathcal{H}_{\infty}$-synthesis problem consists of designing
a feedback control law of the form (3.2) or (3.3) (which generates inputs $u$
based on the measurements $y$) such that the (strong)
$\mathcal{H}_{\infty}$-norm of the resulting closed-loop system from $w$ to
$z$ is minimized, as depicted in Figure 4.3.

Plant
---
Dynamics: (4.15)
Controller
---
Dynamics: (3.2) or (3.3)
$w$$z$$y$$u$ Figure 4.3: Set-up for the design of
$\mathcal{H}_{\infty}$-controllers.

As for stabilization in Chapter 3, we will employ an optimization approach
that consists of directly minimizing the (strong) $\mathcal{H}_{\infty}$-norm
with respect to the controller parameters. To this end, the functions
tds_hiopt_static and tds_hiopt_dynamic can be used. These functions take the
same arguments as respectively tds_stabopt_static and tds_stabopt_dynamic.
More specifically, the first argument of both functions is a `tds_ss`
representation for the open-loop plant, while tds_hiopt_dynamic takes the
desired controller order as a second mandatory argument. The remaining
arguments should be passed as key-value pairs.

  * •

`’options’` allows to specify additional options. The `’options’`-argument
should be a `struct` with the necessary fields. A valid `’options’`-structure
can be created using the tds_hiopt_options-function (for more details see
Section B.10).

  * •

`’initial’` allows to specify the initial values for the optimization
variables. The `’initial’`-argument should be a cell array containing the
initial controllers.

  * •

`’mask’` allows to specify which entries of the matrices in (3.2) or (3.3)
that should remain fixed to their basis value (specified using the `’basis’`
argument), i.e., the entries that should not be touched by the optimization
algorithm.

  * •

`’basis’` allows to specify the value of the entries of (3.2) or (3.3) that
are not optimized (indicated using the `’mask’` argument).

######  Remark 4.4.

The functions tds_hiopt_static and tds_hiopt_dynamic will test whether the
provided initial controllers are strongly stabilizing. If an initial
controller is not strongly stabilizing, the functions tds_stabopt_static and
tds_stabopt_dynamic are used to design a strongly stabilizing controller
starting from the initial value. If this does not result in a strongly stable
controller, this initial controller is skipped in the minimization of the
$\mathcal{H}_{\infty}$-norm. If no strongly stabilizing controller can be
found, an error is thrown.

######  Example 4.4 ([16, Section 7.1]).

Consider the following delay-descriptor system.

$\left\\{\begin{array}[]{rcl}\begin{bmatrix}1&0\\\
0&0\end{bmatrix}\dot{x}(t)&=&\begin{bmatrix}-0.1&-1\\\
1&-1\end{bmatrix}x(t)+\begin{bmatrix}0\\\
1\end{bmatrix}u(t)+\begin{bmatrix}0\\\ 1\end{bmatrix}w(t)\\\\[8.53581pt]
y(t)&=&\begin{bmatrix}0&1\\\
0&0\end{bmatrix}x(t-\tau_{1})+\begin{bmatrix}0&0\\\
0&1\end{bmatrix}x(t-\tau_{2})\\\\[11.38109pt]
z(t)&=&\begin{bmatrix}2&-1\end{bmatrix}x(t)\\\ \end{array}\right.$ (4.16)

with $\tau_{1}=1$ and $\tau_{2}=2$. As always, we start with creating a
`tds_ss` representation.

⬇

E = [1 0;0 0];

A = {[-0.1 -1;1 -1]}; B1 = {[0;1]}; B2 = {[0;1]};

C1 = {[0 1; 0 0],[0 0;0 1]};

C2 = {[2 -1]};

sys = tds_create_ddae(E,A,0,B1,0,C1,[1 2],'B2',B2,'C2',C2);

Next we run the optimization algorithm to design a static feedback controller.

⬇

options = tds_hiopt_options('nstart',1);

K_initial␣=␣[0.25␣-0.5];

[K,cl]␣=␣tds_hiopt_static(sys,'options',options,'initial',K_initial);

We obtain the following feedback law

$u(t)=[-0.3533-0.1012]\ y(t),$ (4.17)

and the strong $\mathcal{H}_{\infty}$-norm of the corresponding closed-loop
system is equal to 1.8333. Figure 4.4(a) shows the spectral norm of the
closed-loop frequency response in function of $\omega$. We observe that it
attains a maximum of 1.833 at $\omega=0$. However, Figure 4.4(b) shows the
same plot after we made a small perturbation to the first delay. At high
frequencies, we now see a jump to almost the $\mathcal{H}_{\infty}$-norm. Let
us therefore compute the strong H-infinity norm of the asymptotic transfer
function using the following code.

⬇

>> Ta = get_asymptotic_transfer_function(cl);

>> hinfa = tds_hinfnorm(Ta)

hinfa =

1.8331

We conclude that (for this example) the standard $\mathcal{H}_{\infty}$-norm
and the strong $\mathcal{H}_{\infty}$-norm of the asymptotic transfer function
(almost) coincide.

(a) $\tau_{1}=1$ and $\tau_{2}=2$

(b) $\tau_{1}=0.999$ and $\tau_{2}=2$

Figure 4.4: Singular value plot of the frequency response of the feedback
interconnection of (4.16) and (4.17) for different values of $\tau_{1}$.

∎

Next, we illustrate how to design a structured controller using
`tds_hiopt_dynamic`.

######  Example 4.5.

In this academic example we will revisit the set-up from Example 3.7, in which
we designed a stabilizing dynamic output feedback controller of order two
under the constraint that the matrices $B_{c}$ and $D_{c}$ are diagonal.
Although the obtained controller was stabilizing, we have no indication of the
performance of the closed-loop system, e.g., “How well are external
disturbances attenuated?”. For the sake of this example, imagine therefore
that the output measurements of the plant are disturbed by a noise signal $w$
(see Figure 4.5).

Plant$+$ $w$$y$Structured controller$u$$z$ Figure 4.5: Situation sketch for
Example 4.5: the control output is disturbed by measurement noise.

We now are interested in quantifying the effect of this measurement noise on
the “real” output of the plant ($z$). To this end, let us extend the state-
space model in (3.21) with the performance input and output channels:

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&A_{0}x(t)+A_{1}x(t-1)+Bu(t-5)\\\
y(t)&=&Cx(t)+w(t)\\\ z(t)&=&Cx(t),\end{array}\right.$ (4.18)

with the matrices $A_{0}$, $A_{1}$, $B$, and $C$ as defined in (3.21). In this
case, $|||T_{zw}(s;\vec{\tau})|||_{\mathcal{H}_{\infty}}$ thus quantifies the
worst-case amplification of measurement noise at the output of the plant. For
the controller parameters obtained in Example 3.7 the strong
$\mathcal{H}_{\infty}$-norm is equal to 187.6219. Next we will design a
dynamic output feedback controller (with the same structure) that locally
minimizes the closed-loop $\mathcal{H}_{\infty}$-norm. To this end, we can use
the following code. (For illustrative purposes, the initial controller is
chosen such that the optimization procedure converges quickly.)

⬇

% Create a representation for (4.18)

A0 = [-3.5 -2 3 -1.5;3 -6.5 3 -1;4 -2 -7 2.5;-2.5 1.5 -3.5 1];

A1 = [-0.25 0.5 -0.25 0.25;0.5 -0.5 1 -0.5;...

1 -1 0.5 -0.15;-0.15 -0.15 0.1 0.2];

B = [-1 2;-2 1;3 -2.5;-4 0];

C = [0 1 1 0;1 0 0 1];

D12 = eye(2);

sys=tds_create({A0,A1},[0 1],{B},5,{C},0,'D12',{D12},'C2',{C});

%␣Define␣the␣mask␣and␣basis␣arguments

mask␣=␣ss(ones(2),eye(2),ones(2),eye(2));

basis␣=␣ss(zeros(2),zeros(2),zeros(2),zeros(2));

%␣Design␣a␣structured␣controller

nc␣=␣2;

initial␣=␣ss([-0.25␣0.75;-1.5␣0.5],diag([0␣-0.5]),...

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣[-0.1␣0.1;-1␣3.5],diag([0.01␣-1]));

options_stabopt␣=␣tds_stabopt_options('nstart',1);

[K,cl]␣=␣tds_hiopt_dynamic(sys,nc,'options',options_stabopt,...

␣␣␣␣␣␣␣␣␣␣␣␣␣'initial',initial,'basis',basis,'mask',mask);

The obtained controller matrices are

$\begin{array}[]{ll}A_{c}=\begin{bmatrix}-0.3029&0.7356\\\
-1.502&0.4897\end{bmatrix}\text{, }&B_{c}=\begin{bmatrix}-0.03063&0\\\
0&-0.271\end{bmatrix}\text{, }\\\\[12.0pt]
C_{c}=\begin{bmatrix}-0.06593&0.04455\\\ -0.9919&3.475\end{bmatrix}\text{, and
}&D_{c}=\begin{bmatrix}0.01272&0\\\ 0&-0.944\end{bmatrix}.\end{array}$

The corresponding closed-loop system is strongly exponentially stable (the
strong spectral abscissa is equal to -0.1046) and the strong
$\mathcal{H}_{\infty}$-norm is now equal to 3.2424. The controller designed
above thus (significantly) reduces the strong $\mathcal{H}_{\infty}$-norm, at
the cost of an increase in the strong spectral abscissa. ∎

Below we retake the previous example in a more practical setting. In a typical
$\mathcal{H}_{\infty}$-controller design problem we not only want to suppress
measurement noise, but also have good disturbance rejection and good reference
tracking. Furthermore, it is often desirable to limit the control effort to
avoid saturating the actuators and to satisfy the physical constraints of the
system. For good practical performance, one thus has to take multiple
objectives into account at the same. However, considering these different
objectives separately can lead to a waterbed effect in which reducing one
objective, increases another one. For example, to have good disturbance
rejection, a large loop gain is required. But, a large loop gain implies that
sensor noise is barely suppressed and that the system is more sensitive to
unmodelled dynamics. Therefore, a trade-off between the different control
goals needs to be made. To achieve this, we will optimize a weighted
combination of these performance objectives. Furthermore, this weighting needs
to be frequency-dependent, as different objectives are relevant in different
frequency ranges. For example, disturbance signals are typically significant
for low frequencies, while measurement noise and modelling uncertainties are
typically significant for higher frequencies.

The design approach sketched above gives rise to a so called mixed-sensitivity
$\mathcal{H}_{\infty}$-synthesis problem, i.e., find a controller that
minimizes the following objective function

$\left\|\begin{bmatrix}W_{1}S\\\ W_{2}KS\\\
W_{3}T\end{bmatrix}\right\|_{\mathcal{H}_{\infty}}$ (4.19)

with $S:=(I-PK)^{-1}$ the sensitivity function, $T:=(I-PK)^{-1}PK$ the
complementary sensitivity function, $P$ the transfer function of the plant,
$K$ the transfer function of the controller and $W_{1}$, $W_{2}$, and $W_{3}$
appropriately chosen weighting functions. The three components in this
objective function, each play a different role.

  * •

By chosing $W_{1}$ large for low frequencies we obtain a small $S$, which
assures good reference tracking and a good disturbance rejection.

  * •

By chosing $W_{2}$ large we obtain a small $KS$, which limits the control
effort.

  * •

By chosing $W_{3}$ large for high frequencies we obtain a small $T$, which
implies good robustness against model uncertainties and good measurement noise
attenuation.

Below we briefly illustrate how such a design problem can be tackled in . For
an extensive overview of how such mixed-sensitivity synthesis problems pop up
in (robust and optimal) controller design and their interpretation in terms of
performance and robust stability, we refer the interested reader to [43] and
[29].

Note. In the current implementation of the employed weighting functions need
to be stable. (A similar constraint applies to the mixsyn function in matlab.)

######  Example 4.6.

In this example we again consider the system from Example 3.7. We now want to
find controller parameters that (locally) minimize (4.19) for

$W_{1}(s)=0.25\frac{s+30}{s+2.5}I_{2}\text{, }W_{2}(s)=0\text{, and
}W_{3}(s)=2\frac{s+125}{s+1000}.$

To tackle this design problem in , we consider the set-up depicted in Figure
4.6. The objective function now corresponds to the $\mathcal{H}_{\infty}$-norm
of the closed-loop transfer function from $w$ to
$\begin{bmatrix}z_{1}&z_{2}&z_{3}\end{bmatrix}^{T}$.

$P$$+$$W_{2}$$W_{1}$$W_{3}$$w$$z_{1}$$z_{2}$$z_{3}$$K$$y$$u$ Figure 4.6:
Situation sketch for Example 4.6.

To incroporate the performance outputs $z_{1}$ and $z_{3}$111In this example,
we ignore $z_{2}$ as $W_{2}(s)=0$, but a similar approach can be used to
incorporate this output., we need to augment the model with two auxiliary
state variables $\xi_{1}$ and $\xi_{3}$:

$\left\\{\begin{array}[]{rcl}\dot{\xi}_{1}(t)&=&-2.5\xi_{1}(t)+(Cx(t)+w(t))\\\
z_{1}(t)&=&6.875\xi_{1}(t)+0.25(Cx(t)+w(t))\end{array}\right.$

and

$\left\\{\begin{array}[]{rcl}\dot{\xi}_{3}(t)&=&-1000\xi_{3}(t)+Cx(t)\\\
z_{3}(t)&=&-1750\xi_{3}(t)+2Cx(t).\end{array}\right.$

Adding these dynamics to (3.21), we obtain the following state-space model

$\left\\{\begin{array}[]{rcl}\begin{bmatrix}\dot{x}(t)\\\ \dot{\xi}_{1}(t)\\\
\dot{\xi}_{3}(t)\end{bmatrix}&=&\begin{bmatrix}A_{0}&0&0\\\ C&-2.5I&0\\\
C&0&-1000I\end{bmatrix}\begin{bmatrix}x(t)\\\ \xi_{1}(t)\\\
\xi_{3}(t)\end{bmatrix}+\begin{bmatrix}A_{1}&0&0\\\ 0&0&0\\\
0&0&0\end{bmatrix}\begin{bmatrix}x(t-1)\\\ \xi_{1}(t-1)\\\
\xi_{3}(t-1)\end{bmatrix}+\begin{bmatrix}B\\\ 0\\\
0\end{bmatrix}u(t-5)+\begin{bmatrix}0\\\ I\\\ 0\end{bmatrix}w(t)\\\\[15.0pt]
y(t)&=&\begin{bmatrix}C&0&0\end{bmatrix}\begin{bmatrix}x(t)\\\ \xi_{1}(t)\\\
\xi_{3}(t)\end{bmatrix}+w(t)\\\\[15.0pt] \begin{bmatrix}z_{1}(t)\\\
z_{3}(t)\end{bmatrix}&=&\begin{bmatrix}0.25C&6.875I&0\\\
2C&0&-1750I\end{bmatrix}\begin{bmatrix}x(t)\\\ \xi_{1}(t)\\\
\xi_{3}(t)\end{bmatrix}+\begin{bmatrix}0.25I\\\
0\end{bmatrix}w(t).\end{array}\right.$

Now, we can apply the same procedure as before.

⬇

% Open-loop system matrices and delays

A0 = [-3.5 -2 3 -1.5;3 -6.5 3 -1;4 -2 -7 2.5;-2.5 1.5 -3.5 1];

A1 = [-0.25 0.5 -0.25 0.25;0.5 -0.5 1 -0.5;...

1 -1 0.5 -0.15;-0.15 -0.15 0.1 0.2];

B = [-1 2;-2 1;3 -2.5;-4 0]; C = [0 1 1 0;1 0 0 1];

tau1 = 1; tau2 = 5; n = 4; p = 2; q = 2;

% weighting functions

K1 = 0.25; a1 = 30; b1 = 2.5; K3 = 2; a3 = 125; b3 = 1000;

% Augmented system matrices

A0e=[A0 zeros(n,2*q);C -b1*eye(q) zeros(q);...

C zeros(q) -b3*eye(q)];

A1e = zeros(n+2*q); A1e(1:n,1:n) = A1;

B1e = [B;zeros(2*q,p)]; B2e = [zeros(n,q);eye(q);zeros(q)];

C1e = [C zeros(q,2*q)];

C2e = [K1*C K1*(a1-b1)*eye(q) zeros(q);...

K3*C zeros(q) K3*(a3-b3)*eye(q)];

D12 = eye(q); D22 = [K1*eye(q);zeros(q)];

sys = tds_create({A0e,A1e},[0 tau1],{B1e},tau2,{C1e},0,...

'D12',{D12},'B2',{B2e},'C2',{C2e},'D22',{D22});

%%␣Define␣the␣mask␣and␣basis␣arguments

mask␣=␣ss(ones(2),eye(2),ones(2),eye(2));

basis␣=␣ss(zeros(2),zeros(2),zeros(2),zeros(2));

%%␣Design␣a␣structured␣controller

nc␣=␣2;

initial␣=␣ss([-0.2␣0.5;-5␣-8],diag([0␣-2]),...

␣␣␣␣␣␣␣␣␣␣␣␣␣[-0.2␣-0.5;-6␣-5],diag([0.01␣-3]));

options_stabopt␣=␣tds_stabopt_options('nstart',1);

[K,cl]␣=␣tds_hiopt_dynamic(sys,nc,'options',options_stabopt,...

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣'initial',initial,'basis',basis,'mask',mask);

We obtain the following controller parameters

$\begin{array}[]{ll}A_{c}=\begin{bmatrix}-0.2215&0.6798\\\
-4.992&-7.975\end{bmatrix}\text{, }&B_{c}=\begin{bmatrix}0.07326&0\\\
0&-2.139\end{bmatrix}\text{, }\\\\[12.0pt]
C_{c}=\begin{bmatrix}-0.04723&-0.07836\\\ -6.029&-5.074\end{bmatrix}\text{,
and }&D_{c}=\begin{bmatrix}5.226\times 10^{-4}&0\\\
0&-2.628\end{bmatrix}\end{array}$

and the corresponding value of the objective function is equal to 10.7073.

∎

To conclude this section, we will discuss how the function tds_hiopt_dynamic
can be used to obtain a delay-free approximation for a stable and strictly
proper time-delay system. To this end we will reformulate the minimization of
the mismatch error in (4.6) as a controller design problem by defining
appropriate control inputs and outputs.

######  Example 4.7.

In this example we look for a delay-free approximation of the form (4.5) for
the following time-delay system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&A_{0}x(t)+A_{1}x(t-1)+B_{2}w(t)\\\
z(t)&=&C_{2}x(t)\end{array}\right.$ (4.20)

with

$A_{0}=\begin{bmatrix}-3&1\\\ 2&-3\end{bmatrix},\
A_{1}=\begin{bmatrix}-2&0.5\\\ 0.7&-1\end{bmatrix},\ B_{2}=\begin{bmatrix}0\\\
1\end{bmatrix}\text{ and }C_{2}=\begin{bmatrix}1&0\end{bmatrix},$

such that the approximation error in (4.6) is minimized. To be able to tackle
this problem using `tds_hiopt_dynamic`, we need to reformulate the
approximation problem as a controller design problem. Let us therefore append
(4.20) with the following control input and output channels:

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&A_{0}x(t)+A_{1}x(t-1)+B_{2}w(t)\\\
y(t)&=&w(t)\\\ z(t)&=&C_{2}x(t)-u(t).\end{array}\right.$ (4.21)

Eliminating the control signals $u$ and $y$, the dynamics of the closed loop
of (4.21) and (3.2) (with $D_{c}=0$) can be described by the

$\left\\{\begin{array}[]{rcl}\begin{bmatrix}\dot{x}(t)\\\
\dot{x}_{c}(t)\end{bmatrix}&=&\begin{bmatrix}A_{0}&0\\\
0&A_{c}\end{bmatrix}\begin{bmatrix}x(t)\\\
x_{c}(t)\end{bmatrix}+\begin{bmatrix}A_{1}&0\\\
0&0\end{bmatrix}\begin{bmatrix}x(t-1)\\\
x_{c}(t-1)\end{bmatrix}+\begin{bmatrix}B_{2}\\\
B_{c}\end{bmatrix}w(t)\\\\[10.03749pt]
z(t)&=&\begin{bmatrix}C_{2}&-C_{c}\end{bmatrix}\begin{bmatrix}x(t)\\\
x_{c}(t)\end{bmatrix}.\end{array}\right.$

The corresponding transfer function from $w$ to $z$ is given by

$\begin{bmatrix}C_{2}&-C_{c}\end{bmatrix}\left(\begin{bmatrix}sI-
A_{0}-A_{1}e^{-s}&0\\\ 0&sI-
A_{c}\end{bmatrix}\right)^{-1}\begin{bmatrix}B_{2}\\\ B_{c}\end{bmatrix}\\\
=C_{2}\left(sI-A_{0}-A_{1}e^{-s}\right)^{-1}B_{2}-C_{c}\left(sI-
A_{c}\right)^{-1}B_{c}.$

The $\mathcal{H}_{\infty}$-norm of this closed-loop transfer function thus
corresponds to the approximation error defined in (4.6). The output feedback
controller with $D_{c}=0$ that minimizes the closed-loop
$\mathcal{H}_{\infty}$-norm of (4.21) thus corresponds to the dynamic system
of the form (4.5) that best approximates the delay system in the sense of
(4.6). Below we compute a delay-free approximation for (4.20) with $n_{a}=5$.

⬇

% Create a representation for the augmented system

A0 = [-3 1;2 -3]; A1 = [-2 0.5;0.7 -1];

B = [0;1]; C = [1 0];

sys = tds_create({A0 A1},[0 1],'D12',1,'B2',B,'C2',C,'D21',-1);

%␣Compute␣the␣approximationn

na␣=␣5;␣%␣order␣of␣the␣delay-free␣approximatio

p␣=␣1;␣q␣=␣1;␣%␣dimensions␣of␣the␣input␣and␣output␣channels

%␣Enforce␣D_c␣=␣0

mask␣=␣ss(ones(na),ones(na,p),ones(q,na),zeros(q,p));

basis␣=␣ss(zeros(na),zeros(na,p),zeros(q,na),zeros(q,p));

%␣Set␣initial␣controller

initial=ss([-0.8␣0␣-0.35␣0.14␣-0.1;-0.3␣-0.5␣-0.5␣0.14␣0.7;...

␣␣␣␣␣␣␣␣0.18␣0.21␣-0.5␣-0.3␣-1.2;-0.04␣-0.2␣-0.14␣-0.9␣0.4;...

␣␣␣␣␣␣␣␣-0.2␣-0.86␣1␣-0.33␣-0.4],[0.66;-0.5;0.04;0.2;-0.3],...

␣␣␣␣␣␣␣␣[2␣0␣1␣2␣0.4],1);

options␣=␣tds_hiopt_options('nstart',1);

[approx,CL]␣=␣tds_hiopt_dynamic(sys,na,'options',options,...

␣␣␣␣␣␣␣␣␣␣'initial',initial,'mask',mask,'basis',basis);

Figure 4.7 compares the original time-delay system and its delay-free
approximation for $n_{a}=5$. We see that the approximation matches the
original system quite well and the approximation error
$\|T(s)-T_{\mathrm{approx.}}(s)\|_{\mathcal{H}_{\infty}}$ is equal to
7.6316e-03.

Note.  Due to the employed definition of the $\mathcal{H}_{\infty}$-norm for
unstable systems i.e., the $\mathcal{H}_{\infty}$-norm of an unstable system
is equal to infinity, the obtained approximation is always stable.

Figure 4.7: The magnitude frequency response of (4.20) and its delay-free
approximation with $n_{a}=5$.

∎

####  4.2.1 Mixed performance criteria

It is also possible to design controllers that minimize a combination of the
strong spectral abscissa and the strong $\mathcal{H}_{\infty}$-norm. More
specifically, we will consider the objective function

$\alpha C+(1-\alpha)|||T_{zw}(s;\vec{\tau})|||_{\mathcal{H}_{\infty}}$ (4.22)

with $\alpha\in[0,1]$ a weighting coefficient that determines the trade-off
between minimizing the strong spectral abscissa and the strong
$\mathcal{H}_{\infty}$-norm. We will demonstrate how we can optimize this
alternative objective function using in the following example.

######  Example 4.8.

Consider the following time-delay system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&\begin{bmatrix}1&-2&4\\\ 3&0.5&-1\\\
2&0.4&-2\end{bmatrix}x(t)+\begin{bmatrix}1.5&0.3&2\\\ 0.7&-0.8&0.4\\\
0.5&0.4&-0.9\end{bmatrix}x(t-1)\\\\[17.07164pt]
&&\qquad\qquad\qquad\qquad+\begin{bmatrix}0.3&0.4\\\ -0.7&-0.5\\\
0.7&-0.1\end{bmatrix}u(t-0.1)+\begin{bmatrix}-0.7\\\ -0.5\\\
-0.3\end{bmatrix}w(t)\\\\[19.91692pt] y(t)&=&\begin{bmatrix}-1&0.3&0\\\
0.4&0.9&1\end{bmatrix}x(t)\\\\[11.38109pt]
z(t)&=&\begin{bmatrix}3&-5&-4\end{bmatrix}x(t).\end{array}\right.$ (4.23)

The code below designs three static output feedback controllers for (4.23) by
minimizing the objective function in (4.22) for $\alpha=1$ (i.e., minimizing
the spectral abscissa), for $\alpha=0$ (i.e., minimizing the
$\mathcal{H}_{\infty}$-norm), and for $\alpha=0.5$ (i.e., a trade-off between
minimizing the spectral abscissa and the $\mathcal{H}_{\infty}$-norm).

⬇

% Create a representation for the time-delay system

A0 = [1 -2 4;3 0.5 -1;-2 0.4 -2];

A1 = [1.5 0.3 2;0.7 -0.8 0.4;0.5 0.4 -0.9];

B1 = [0.3 0.4;-0.7 -0.5;0.7 -0.1]; B2 = [-0.7;-0.5;-0.3];

C1 = [-1 0.3 0;0.4 0.9 1]; C2 = [3 -5 -4];

sys=tds_create({A0,A1},[0 1],{B1},0.1,{C1},0,'B2',{B2},'C2',{C2});

%␣alpha␣=␣1

options␣=␣tds_hiopt_options('nstart',1,'alpha',1);

[Dc1,cl1]=tds_hiopt_static(sys,'options',options,'initial',zeros(2));

%␣alpha␣=␣0

options␣=␣tds_hiopt_options('nstart',1,'alpha',0);

[Dc2,cl2]=tds_hiopt_static(sys,'options',options,'initial',Dc1);

%␣alpha␣=␣0.5

options␣=␣tds_hiopt_options('nstart',1,'alpha',0.5);

[Dc3,cl3]=tds_hiopt_static(sys,'options',options,'initial',Dc1);

The obtained controllers are

$\displaystyle D_{c,1}$ $\displaystyle\,=\begin{bmatrix}8.4197&-0.4036\\\
4.3451&10.6842\end{bmatrix}\text{,}$ (4.24) $\displaystyle D_{c,2}$
$\displaystyle\,=\begin{bmatrix}7.0877&\phantom{-}0.0571\\\
6.5345&13.3927\end{bmatrix},\text{ and }$ (4.25) $\displaystyle D_{c,3}$
$\displaystyle\,=\begin{bmatrix}8.1324&-0.7980\\\
4.7536&11.2210\end{bmatrix}.$ (4.26)

Table 4.1 gives the correspoding closed-loop spectral abscissa and
$\mathcal{H}_{\infty}$-norm. As expected, when decreasing $\alpha$ more
emphasis is put on minimizing the $\mathcal{H}_{\infty}$-norm at the cost of
an increase in the spectral abscissa.

Table 4.1: The spectral abscissa and the $\mathcal{H}_{\infty}$-norm of the feedback interconnection of (4.23) and the controllers in (4.24), (4.25), and (4.25). | spectral abscissa |  $\mathcal{H}_{\infty}$-norm
---|---|---
$\alpha=1$ | -0.8751 | 1.8061
$\alpha=0$ | -0.6153 | 1.2908
$\alpha=0.5$ | -0.8285 | 1.4869

∎

As mentioned before, the $\mathcal{H}_{\infty}$-norm corresponds to the
reciprocal of the distance to instability with respect to a complex-valued
perturbation acting on the zero-delay state matrix. From a practical point of
view, the $\mathcal{H}_{\infty}$-norm thus only provides a crude approximation
of the perturbation that is necessary to destabilize the system, as the
perturbations on the model are typically real-valued and may act on all system
matrices and even the delays. In the next section we will therefore consider a
more realistic distance to instability.

###  4.3 Pseudospectrum and the distance to instability

In this section we will consider a class of uncertain retarded time-delay
systems with bounded perturbations on both the state matrices and the delays.
These perturbations allow to account for mismatches between model and reality.
For example, to describe a model parameter that is only known up to a certain
precision (e.g., due to an imprecise measuring device), we add a real-valued
scalar perturbation to the model. Now instead of a single model, we obtain a
collection of possible models. To examine the stability of this uncertain
system for all possible parameter values, this section introduces the notions
of the pseudospectral abscissa and the distance to stability.

######  Example 4.9 (Inspired by [3]).

Consider the following simple time-delay model for a turning process

$\ddot{d}(t)+2\xi\omega\,\dot{d}(t)+\omega^{2}\,d(t)=\frac{k}{m}\,\big{(}d(t-\tau)-d(t)\big{)}$

with $\omega$ the natural angular frequency, $\xi$ the damping ratio, $m$ the
modal mass, $k$ the cutting force coefficient, and $\tau$ the time of one
revolution of the work-piece. This model can be transformed into a standard
first-order delay differential equation of retarded type by considering the
state vector $x(t)=\begin{bmatrix}d(t)&\dot{d}(t)\end{bmatrix}^{\top}$

$\dot{x}(t)=\begin{bmatrix}0&1\\\
-\omega^{2}-\frac{k}{m}&-2\xi\omega\end{bmatrix}x(t)+\begin{bmatrix}0&0\\\
\frac{k}{m}&0\end{bmatrix}x(t-\tau).$ (4.27)

In this section we will examine the stability of (4.27) when the parameters
$k$, $\xi$, and $\tau$ are uncertain. More specifically, we will assume that
the parameter $k$ lies in the interval $[k-\Delta_{k},k+\Delta_{k}]$, $\xi$ in
$[\xi-\Delta_{\xi},\xi+\Delta_{\xi}]$, and $\tau$ in
$[\tau-\Delta_{\tau},\tau+\Delta_{\tau}]$. To model this situation we need two
perturbations ($\delta_{1}$ and $\delta_{2}$) acting on the state matrices and
one perturbation acting on the delay value $\delta\tau_{1}$:

$\small\dot{x}(t)=\left(\begin{bmatrix}0&1\\\
{\small-\omega^{2}-\frac{k}{m}}&{\small-2\xi\omega}\end{bmatrix}+\begin{bmatrix}0\\\
\frac{-1}{m}\end{bmatrix}\delta_{1}\begin{bmatrix}\Delta_{k}&0\end{bmatrix}+\begin{bmatrix}0\\\
-2\omega\end{bmatrix}\delta_{2}\begin{bmatrix}0&\Delta_{\xi}\end{bmatrix}\right)x(t)+\\\
\left(\begin{bmatrix}0&0\\\ \frac{k}{m}&0\end{bmatrix}+\begin{bmatrix}0\\\
\frac{1}{m}\end{bmatrix}\delta_{1}\begin{bmatrix}\Delta_{k}&0\end{bmatrix}\right)x\big{(}t-(\tau+\Delta_{\tau}\delta\tau_{1})\big{)}$
(4.28)

with $|\delta_{1}|\leq 1$, $|\delta_{2}|\leq 1$, and $|\delta\tau_{1}|\leq 1$.
Notice that in this case the state matrix $A_{0}$ is affectected by both
perturbations and that the perturbation $\delta_{1}$ affects both system
matrices. ∎

Let us start by introducing the proper notation. Firstly, we will consider the
set

$\delta:=(\delta_{1},\dots,\delta_{L_{M}})$ (4.29)

which consists of $L_{M}$ matrix perturbations that can have aribrary
dimensions and can either be real or complex-valued, i.e.,

$\delta_{l}\in\mathbb{R}^{\rho_{l}\times\varsigma_{l}}\text{ or
}\delta_{l}\in\mathbb{C}^{\rho_{l}\times\varsigma_{l}}\text{ for
}l=1,\dots,L_{M}.$

Similarly, we consider the set

$\delta\tau:=(\delta\tau_{1},\dots,\delta\tau_{L_{D}})$

which consists of $L_{D}$ real-valued delay perturbations. Further, we assume
that these perturbations are bounded, i.e., there exists a bound $\epsilon$
such that

$\|\delta_{l}\|_{2/\mathrm{fro}}\leq\epsilon\text{ for }l=1,\dots,L_{M}$

and

$|\delta\tau_{l}|\leq\epsilon\text{ for }l=1,\dots,L_{D},$

in which $\|\delta_{l}\|_{2/\mathrm{fro}}$ denotes that we can either use the
spectral norm (the induced 2-norm) or the Frobenius norm (the element-wise two
norm). (Furthermore, for different matrix perturbations, different norms may
be used.) Now that we have introduced these perturbations, we can consider the
following uncertain dealy differential equation:

$\dot{x}(t)=\big{(}A_{0}+\tilde{A}_{0}(\delta)\big{)}\,x(t)+\sum_{k=1}^{m}\big{(}A_{k}+\tilde{A}_{k}(\delta)\big{)}\,x\big{(}t-(h_{A,k}+\tilde{h}_{A,k}(\delta\tau))\big{)}$
(4.30)

in which the matrices $A_{0},\dots,A_{m}$ all belong to $\mathbb{R}^{n\times
n}$ and represent the nominal state matrices. Similarly, the scalars
$h_{1,m}\dots,h_{A,m}$ represent the nominal delay values. The uncertain
matrices $\tilde{A}_{0}(\delta),\dots,\tilde{A}_{m}(\delta)$ are defined as

$\tilde{A}_{k}(\delta)=\textstyle\sum\limits_{j=1}^{N_{A_{k}}}G_{k,j}\
\delta_{\texttt{ind}_{A_{k}}[j]}\ H_{k,j}$

with $N_{A_{k}}$ the number of uncertain terms, $\texttt{ind}_{A_{k}}[j]$ the
$j$th element of the array $\texttt{ind}_{A_{k}}$ that indicates the index of
the uncertainty in $\delta$ and $G_{k,j}$ and $H_{k,j}$ real-valued shape
matrices of appropriate dimensions. (Note that this seemingly cumbersome
notation is necessary to account for the fact that one matrix may be affected
by multiple uncertainties and that one uncertainty may affect multiple
matrices.) Similarly, the uncertain delays $\tilde{h}_{A,k}(\delta\tau)$ are
defined as

$\tilde{h}_{A,k}(\delta\tau)=\max\left\\{-h_{A,k},\textstyle\sum\limits_{j=1}^{N_{\tau_{k}}}w_{k,j}\,\delta\tau_{\texttt{ind}_{h_{A,k}}[j]}\right\\}$

with $N_{h_{A,k}}$ the number of uncertain terms, $\texttt{ind}_{h_{A,k}}[j]$
the $j$th element of the array $\texttt{ind}_{\tau_{k}}$ that indicates the
index of the uncertainty in $\delta\tau$ and $w_{k,j}$ a real-valued weighting
coefficient.

Given the definitions for the perturbations above and the uncertain DDE
(4.30), the first question that we want to address is: “Given a bound
$\epsilon$, is the uncertain system (4.30) exponentially stable for all
instances of the uncertainties smaller than $\epsilon$?”. To answer this
question, we will introduce the $\epsilon$-pseudospectral abscissa as the
worst-case value of the spectral abscissa of (4.30) over all uncertainty
values.

######  Definition 4.2.

The $\epsilon$-pseudospectral abscissa of (4.30) will be denoted by
$c^{\mathrm{ps}}(\epsilon)$ and is defined as follows:

$c^{\mathrm{ps}}(\epsilon):=\max\\{c(\delta,\delta\tau):\|\delta_{l}\|\leq\epsilon\text{
for }l=1,\dots,L_{M}\ \&\ |\delta\tau_{l}|\leq\epsilon\text{ for
}l=1,\dots,L_{D}\\}$

with $c(\delta,\delta\tau)$ the spectral abscissa of (4.30) for a given
instance of the uncertainties.

The uncertain system (4.30) is thus exponentially stable for all instances of
the uncertainties bounded by $\epsilon$ if (and only if) the
$\epsilon$-pseudospectral abscissa is strictly negative.

Secondly, we are interested in answering the following related question:
“Given a set of uncertainties and an uncertain system of the form (4.30), what
is the largest $\epsilon$ such that (4.30) is exponentially stable for all
instances of the uncertainties smaller than $\epsilon$?” or equivalently,
“What is the smallest $\epsilon$ such that there exists an instance of (4.30)
with the norm uncertainties smaller than or equal to $\epsilon$ that is not
exponentially stable?”. This leads us to the notion of distance to
instability.

######  Definition 4.3.

The distance to instability of (4.30) is defined as the smallest $\epsilon$
for which the spectral abscissa is positive, or in other words

$d_{\mathrm{INS}}=\min\\{\epsilon\geq 0:c^{ps}(\epsilon)\geq 0\\}.$

Important: Currently only supports the computation of the pseudospectral
abscissa and the distance to instability for uncertain systems with real-
valued (`’r’`) uncertainties bounded in Frobenius norm (`’f’`).

######  Example 4.10.

In this example, we will retake the set-up in Example 4.9 using the following
parameter values: $k=8\times 10^{6}\,\mathrm{N}\mathrm{m}^{-1}$,
$\omega=775\,\mathrm{s}^{-1}$ , $\tau=0.012\,\mathrm{s}$, $\xi=0.05$,
$m=50\,\textrm{kg}$, $\Delta_{k}=1\times 10^{6}\,\mathrm{N}\mathrm{m}^{-1}$,
$\Delta_{\xi}=0.005$, and $\Delta_{\tau}=0.001\,\mathrm{s}$.

Let us start by representing the nominal system (4.27) in

⬇

k=8e6; omega=775; xi=0.05; tau=0.008; m=50;

A0 = [0 1;-omega^2-k/m -2*xi*omega]; A1 = [0 0;k/m 0];

tds = tds_create({A0,A1},[0 tau]);

and computing its nominal spectral abscissa

⬇

>> tds_strong_sa(tds,-30)

ans =

-16.3646

We conclude that the system is exponentially stable.

Next we will see how the uncertain RDDE (4.28) can be represented in . We
start by creating a representation for the uncertainties $\delta_{1}$ and
$\delta_{2}$ using the function tds_create_delta

⬇

delta=tds_create_delta({1,1},{1,1},{'r','r'},{'f','f'});

This function takes four cell arrays of length $L_{M}$ (here equal to 2) as
input. The first cell indicates the row dimension of the uncertainties, the
second the column dimension, the third whether the uncertainty is real (`’r’`)
or complex-valued (`’c’`), and the fourth whether its size should be measured
using the spectral (`’s’`) or Frobenius norm (`’f’`).
Next, we create objects that represents the uncertain matrices
$\tilde{A}_{0}(\delta)$ and $\tilde{A}_{1}(\delta)$ using the function
tds_uncertain_matrix:

⬇

Delta_k = 1e6; Delta_xi = 0.005;

G01 = [0;-1/m]; H01 = [Delta_k 0];

G02 = [0;-2*omega]; H02 = [0 Delta_xi];

uA0 = tds_uncertain_matrix({1,2},{G01,G02},{H01,H02});

G11 = [0;1/m]; H11 = [Delta_k 0];

uA1 = tds_uncertain_matrix({1},{G11},{H11});

The function `tds_uncertain_matrix` takes three cell arrays as input argument:
the first argument is the array $\texttt{ind}_{A_{k}}$, the second one is a
cell array that stores the left shape matrices, and the third one is a cell
array storing the right shape matrices.
Subsequently we create the uncertain delay term
$\tilde{\tau}_{1}(\delta\tau)=\Delta_{\tau}\delta\tau_{1}$:

⬇

Delta_tau = 0.001;

uhA1 = tds_uncertain_delay({1},{Delta_tau});

in which the function `tds_uncertain_delay` takes two cell arrays as input:
the first is the array $\texttt{ind}_{\tau_{k}}$ and the second is a cell
array with the weights. Finally we use the method `add_uncertainty` to add
these uncertainties to the nominal time delay system `tds`.

⬇

utds = tds.add_uncertainty(delta,1,{uA0,uA1},{[],uhA1});

In this case, the method `add_uncertainty` takes four input arguments: the set
of matrix uncertainties $\delta$, the number of delay uncertainties (in this
case 1) and two cell arrays of length $m_{A}$, with respectively the uncertain
matrices and delays222If a term is not affected by the uncertainties, then one
can just use an empty array instead of the uncertain matrix or delay.. The
function returns an `utds_ss`-object that represents the uncertain delay
differential equation of which we can compute the $\epsilon$-pseudospectral
abscissa using the function `tds_psa`.

⬇

>> epsilon = 1;

>> tds_psa(utds,epsilon)

ans =

20.7022

We thus conclude that the uncertain delay differential equation (4.28) is not
stable for all admissible uncertainty values. Let us therefore compute the
corresponding distance to instability.

⬇

tds_dins(utds,[0 1]);

We find that $r_{\mathrm{ins}}=0.4125$, meaning that the uncertain system is
stable for all $k$ in
$(k-r_{\mathrm{ins}}\Delta_{k},k+r_{\mathrm{ins}}\Delta_{k})$, $\xi$ in $(\xi-
r_{\mathrm{ins}}\Delta_{\xi},\xi+r_{\mathrm{ins}}\Delta_{\xi})$, and $\tau$ in
$(\tau-r_{\mathrm{ins}}\Delta_{\tau},\tau+r_{\mathrm{ins}}\Delta_{\tau})$.

∎

## Acknowledgments

The software package builds on twenty years of research in the group of Wim
Michiels from the Numerical Analysis and Applied Mathematics (NUMA) section of
KU Leuven, Belgium. The package integrates several methods and algorithms,
some of which gave rise to specific software tools available from
https://twr.cs.kuleuven.be/research/software/delay-control/ Writing the
package would not have been possible without the contribution of several
(former) group members and collaborators. Therefore, in the first place, the
authors want to thank Suat Gumussoy, who strongly contributed to the
methodology and software tools for $\mathcal{H}_{\infty}$ analysis and design.
Also the adopted specification of systems and controllers in the package
carries the stamp of Suat’s work. Secondly, the authors want to acknowledge
the contributions of Haik Silm from which this software package greatly
benefited. Special thanks further goes to Zhen Wu (algorithm for automatically
determining the resolution of the discretization such that all characteristic
roots in a given half plane are captured), Francesco Borgioli (computation of
real structured pseudospectral abscissa and stability radii), Joris
Vanbiervliet (the first work on stabilization of time-delay systems via the
minimization of the spectral abscissa) and Deesh Dileep (design of
structurally constrained controllers). Dan Pilbauer, Adrian Saldanha and Tomáš
Vyhlídal contributed to the validation and first applications of the design
tools to DDAE models, in the context of vibration control applications.

uses `HANSO` (version 3.0) to solve the non-smooth, non-convex optimization
problems arising in the context of controller-design. `HANSO` is available for
download from https://cs.nyu.edu/~overton/software/hanso/. `HANSO` is based on
the methods described in [23, 7] and is licensed under the GNU General Public
License version 3.0.

## Bibliography

  * [1] P. Apkarian and D. Noll  “Nonsmooth H$\infty$ Synthesis”  In _IEEE Transactions on Automatic Control_ 51.1, 2006, pp. 71–86
  * [2] “Differential-Difference Equations” 6, Mathematics in Science and Engineering  New York (NY): Academic Press, 1963
  * [3] Francesco Borgioli and Wim Michiels  “A novel method to compute the structured distance to instability for combined uncertainties on delays and system matrices”  In _IEEE Transactions on Automatic Control_ 65.4 IEEE, 2019, pp. 1747–1754
  * [4] Stephen Boyd and Venkataramanan Balakrishnan  “A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L$\infty$-norm”  In _Systems & Control Letters_ 15.1 Elsevier, 1990, pp. 1–7
  * [5] Dimitri Breda, Stefano Maset and Rossana Vermiglio  “Stability of linear delay differential equations: A numerical approach with MATLAB”  New York (NY): Springer, 2014
  * [6] NA Bruinsma and M Steinbuch  “A fast algorithm to compute the H$\infty$-norm of a transfer function matrix”  In _Systems & Control Letters_ 14.4 Elsevier, 1990, pp. 287–293
  * [7] James V Burke et al.  “Gradient sampling methods for nonsmooth optimization”  In _Numerical Nonsmooth Optimization_ Springer, 2020, pp. 201–225
  * [8] J.V. Burke, D. Henrion, A.S. Lewis and M.L. Overton  “HIFOO - a MATLAB package for fixed-order controller design and H$\infty$ optimization” 5th IFAC Symposium on Robust Control Design  In _IFAC Proceedings Volumes_ 39.9, 2006, pp. 339–344
  * [9] Deesh Dileep et al.  “Design of robust decentralised controllers for MIMO plants with delays through network structure exploitation”  In _International Journal of Control_ 93.10 Taylor & Francis, 2020, pp. 2275–2289
  * [10] Thomas Erneux  “Applied delay differential equations”
  * [11] Emilia Fridman  “Stability of linear descriptor systems with delay: a Lyapunov-based approach”  In _Journal of Mathematical Analysis and Applications_ 273.1 Elsevier, 2002, pp. 24–44
  * [12] Emilia Fridman  “Introduction to time-delay systems: Analysis and control”  Cham: Birkhäuser, 2014
  * [13] Emilia Fridman and Uri Shaked  “H$\infty$-control of linear state-delay descriptor systems: an LMI approach”  In _Linear algebra and its applications_ 351 Elsevier, 2002, pp. 271–302
  * [14] Keqin Gu, Jie Chen and Vladimir L Kharitonov  “Stability of time-delay systems”  Boston (MA): Birkhäuser, 2003
  * [15] Suat Gümüşsoy, Didier Henrion, Marc Millstone and Michael L Overton  “Multiobjective robust control with HIFOO 2.0”  In _IFAC Proceedings Volumes_ 42.6 Elsevier, 2009, pp. 144–149
  * [16] Suat Gümüşsoy and Wim Michiels  “Fixed-order H-infinity control for interconnected systems using delay differential algebraic equations”  In _SIAM Journal on Control and Optimization_ 49.5 SIAM, 2011, pp. 2212–2238
  * [17] Jack K Hale and Sjoerd M Verduyn Lunel  “Introduction to functional differential equations” 99, Applied Mathematical Sciences  New York (NY): Springer, 1993
  * [18] Jack K Hale and Sjoerd M Verduyn Lunel  “Strong stabilization of neutral functional differential equations”  In _IMA Journal of Mathematical Control and Information_ 19.1,2 OUP, 2002, pp. 5–23
  * [19] Diederich Hinrichsen and Anthony J Pritchard  “Mathematical Systems Theory I: Modelling, State Space Analysis, Stability and Robustness”  Heidelberg (Germany): Springer Berlin, 2005
  * [20] Tamás Insperger and Gábor Stépán  “Stability of the milling process”  In _Periodica Polytechnica Mechanical Engineering_ 44.1, 2000, pp. 47–57
  * [21] In _Stability of Functional Differential Equations_ 180, Mathematics in Science and Engineering  Elsevier, 1986
  * [22] Nikolai Nikolaevich Krasovskii  “Stability of motion: Applications of Lyapunov’s Second Method to Differential Systems and Equations with Delay” (translated from Russian)  In _Stanford University Press_ , 1963
  * [23] Adrian S Lewis and Michael L Overton  “Nonsmooth optimization via quasi-Newton methods” Code available from https://cs.nyu.edu/~overton/software/hanso/.  In _Mathematical Programming_ 141.1 Springer, 2013, pp. 135–163
  * [24] Wim Michiels  “Spectrum-based stability analysis and stabilisation of systems described by delay differential algebraic equations”  In _IET control theory & applications_ 5.16, 2011, pp. 1829–1842
  * [25] Wim Michiels and Silviu-Iulian Niculescu  “On the delay sensitivity of Smith Predictors”  In _International Journal of Systems Science_ 34.8-9 Taylor & Francis, 2003, pp. 543–551  DOI: 10.1080/00207720310001609057
  * [26] Wim Michiels and Silviu-Iulian Niculescu  “Stability, control, and computation for time-delay systems. An eigenvalue based approach”  Philadelphia (PA): SIAM, 2014
  * [27] Wim Michiels and Bin Zhou  “On the strong relative degree of time-delay systems with non-commensurate delays”  In _SIAM Journal on Control and Optimization_ 60.6, 2022, pp. 3384–3408  DOI: 10.1137/21M1442656
  * [28] N. Olgac and B.T. Holm-Hansen  “A Novel Active Vibration Absorption Technique: Delayed Resonator”  In _Journal of Sound and Vibration_ 176.1, 1994, pp. 93–104
  * [29] Hitay Özbay, Suat Gümüşsoy, Kenji Kashima and Yutaka Yamamoto  “Frequency Domain Techniques for H$\infty$ Control of Distributed Parameter Systems”  Philadelphia (PA): SIAM, 2018  DOI: 10.5555/3317188
  * [30] K. Pyragas  “Continuous control of chaos by self-controlling feedback”  In _Physics Letters A_ 170.6, 1992, pp. 421–428  DOI: 10.1016/0375-9601(92)90745-8
  * [31] Adrián Ramírez, Sabine Mondié, Rubén Garrido and Rifat Sipahi  “Design of proportional-integral-retarded (PIR) controllers for second-order LTI systems”  In _IEEE Transactions on Automatic Control_ 61.6 IEEE, 2015, pp. 1688–1693
  * [32] B.S. Razumikhin  “On the Stability of Systems with Delay”  In _Prikladnaya Matematika i Mekhanika_ 20, 1956, pp. 500–512
  * [33] Rifat Sipahi et al.  “Stability and stabilization of systems with time delay”  In _IEEE Control Systems Magazine_ 31.1 IEEE, 2011, pp. 38–65
  * [34] O. J. Smith  “Closer Control of Loops with Dead Time”  In _Chemistry Engineering Progress_ 53.5 American Institute of Chemical Engineers (AIChE), 1957, pp. 217–219
  * [35] Joris Vanbiervliet, Koen Verheyden, Wim Michiels and Stefan Vandewalle  “A nonsmooth optimisation approach for the stabilisation of time-delay systems”  In _ESAIM: Control, Optimisation and Calculus of Variations_ 14.3 EDP Sciences, 2008, pp. 478–493
  * [36] Koen Verheyden, Tatyana Luzyanina and Dirk Roose  “Efficient computation of characteristic roots of delay differential equations using LMS methods”  In _Journal of Computational and Applied Mathematics_ 214.1, 2008, pp. 209–226
  * [37] Erik I. Verriest and Wim Michiels  “Inverse Routh table construction and stability of delay equations”  In _Systems & Control Letters_ 55.9, 2006, pp. 711–718  DOI: https://doi.org/10.1016/j.sysconle.2006.02.002
  * [38] Tomáš Vyhlídal  “Analysis and synthesis of time delay system spectrum”, 2003
  * [39] Tomáš Vyhlídal, Wim Michiels and Paul McGahan  “Synthesis of strongly stable state-derivative controllers for a time-delay system using constrained non-smooth optimization”  In _IMA Journal of Mathematical Control and Information_ 27.4, 2010, pp. 437–455  DOI: 10.1093/imamci/dnq025
  * [40] Tomáš Vyhlídal, Nejat Olgac and Vladimír Kučera  “Delayed resonator with acceleration feedback - Complete stability analysis by spectral methods and vibration absorber design”  In _Journal of Sound and Vibration_ 333.25, 2014, pp. 6781–6795  DOI: https://doi.org/10.1016/j.jsv.2014.08.002
  * [41] Zhen Wu and Wim Michiels  “Reliably computing all characteristic roots of delay differential equations in a given right half plane using a spectral method”  In _Journal of Computational and Applied Mathematics_ 236.9 Elsevier, 2012, pp. 2499–2514
  * [42] Qing-Chang Zhong  “Robust Control of Time-delay Systems”  Londen (UK): Springer-Verlag, 2006
  * [43] Kemin Zhou, John C. Doyle and Keith Glover  “Robust and Optimal Control”  Upper Saddle River (NJ): Prentice Hall, 1996
  * [44] Kemin Zhou and John Comstock Doyle  “Essentials of robust control”  Upper Saddle River (NJ): Prentice hall, 1998

##  Appendix A Installation

The package is available for downloaded from
https://gitlab.kuleuven.be/u0011378/tds-control/-/archive/main/tds-control-
main.zip?path=tds-control. To install the package, simply add the directory
code| and all its subdirectories to your matlab search path:

⬇

>> addpath(genpath('code'));

Important. To be able to use the controller design algorithms, such as
tds_stabopt_static and tds_stabopt_dynamic,

⬇

>> addpath(genpath('hanso3_0'));

To add and

⬇

>> savepath

which saves the current search path.

Important. utilizes functionality from the MathWorks® Control System Toolbox™,
such as the ss-class and the sigma-function. If these functions are not
present in your matlab search path, will revert to a fallback implementation.
Note however, both performance and accuracy might be reduced.

##  Appendix B Documentation

###  B.1 tds_roots

The function tds_roots computes the characteristic roots of a given delay
differential equation or time-delay system in a specified right half-plane or
rectangular region. To this end, it implements the method from [41]. First, a
generalized eigenvalue problem of the form

$\Sigma_{N}\,v=\lambda\Pi_{N}\,v$ (B.1)

with $\Sigma_{N}$ and $\Pi_{N}$ belonging to
$\mathbb{R}^{(N+1)n\times(N+1)n}$, is constructed such that the eigenvalues of
(B.1) approximate the desired characteristic roots. More specifically, the
generalized eigenvalue problem (B.1) corresponds to a spectral discretization
of the infinitesimal generator of the solution operator underlying the delay
differential equation. (For more details on spectral discretisation methods
for computing the characteristic roots of time-delay systems, we refer the
interested reader to [5].) The approximations for the characteristic roots of
the time-delay system obtained by computing the eigenvalues of (B.1) are
subsequently improved by applying Newton’s method. Note that the dimensions of
the matrices $\Sigma_{N}$ and $\Pi_{N}$ depend on a parameter $N$, the degree
of the spectral discretisation. This $N$ is related to the quality of the
spectral approximation: the larger $N$, the better individual characteristic
roots are approximated and the larger the region in the complex plane in which
the characteristic roots are well approximated. However, the larger $N$, the
larger the size of the generalized eigenvalue problem and hence the longer the
computation time. An interesting feature of the method from [41] is that it
includes an empirical strategy to find a suitable $N$ such that the desired
characteristic roots are well-approximated, while keeping the generalized
eigenvalue problem as small as possible.

Note: As solving large eigenvalue problems is computationally costly, the
maximal size of the generalized eigenvalue problem (B.1) is limited. More
specifically, if the inequality $N(n+1)>{\tt options.max_size_evp}$ holds,
then $N$ is lowered to the largest $N$ for which $N(n+1)\leq{\tt
options.max_size_evp}$ holds. However, lowering $N$ means that certain
characteristic roots in the desired region are not sufficiently well
approximated and might be missed by tds_roots. See also Example 2.2.

####  B.1.1 Documentation

  * •

L = tds_roots(tds,region) returns a vector L containing the characteristic
roots of tds in the right half-plane $\\{z\in\mathbb{C}:\Re(z)\geq\text{
{region} }\\}$ if region is a scalar and in the rectangular region
$[\texttt{region}(1),\,\texttt{region}(2)]\times\jmath\,[\texttt{region}(3),\,\texttt{region}(4)]$
if region is a vector (of length 4).

  * •

[L,LD] = tds_roots(tds,region) also returns a vector LD containing the
eigenvalues of (B.1).

  * •

[L,LD,info] = tds_roots(tds,region) also returns a structured array info
containing the following information:

    * –

N \- the employed degree of the spectral discretisation.

    * –

N_original \- the degree for the spectral discretization deemed necessary to
capture all characteristic roots in the desired region (N is smaller than
N_original if the maximal size of the generalized eigenvalue problem would be
exceeded for N_original).

    * –

max_size_evp_enforced \- flag indicating that for N_original, the maximal size
for the eigenvalue problem (B.1) would be exceeded and that the degree of the
spectral discretisation is therefore lowered to N.

    * –

gamma_r_exceeds_one \- flag indicating that $\gamma(\text{ {region} })\geq 1$,
meaning that the considered right half-plane contains infinitely many
characteristic roots or that there exist infinitesimal delay perturbations
such that it does. The default value $N=15$ is therefore used. See also the
warning Case: gamma(region) >= 1.

    * –

newton_initial_guesses \- array with the initial guesses for Newton’s method.

    * –

newton_final_values \- array of the same length with the values after Newton’s
method.

    * –

newton_residuals \- array of the same length with the norm of the residuals
after the Newton corrections.

    * –

newton_unconverged_initial_guesses \- array containing the indices of the
initial guesses in newton_initial_guesses for which Newton’s method failed to
converge.

    * –

newton_large_corrections \- array containing the indices for which the result
after Newton’s method lies far from the initial guess.

  * •

[..] = tds_roots(tds,region,options) allows to specify additional options.

    * –

fix_N: fixes the parameter $N$ of the generalized eigenvalue problem (B.1).
Note that the constaint on the maximal size of the generalized eigenvalue
problem is now ignored.

    * –

max_size_evp: gives the maximal size of (B.1) (if fix_N is not specified).
Default value -  600.

    * –

newton_tol tolerance on the norm of the residuals for Newton’s method.
Default value -  1e-10.

    * –

newton_max_iter: maximum number of Newton iterations for correcting a given
root.
Default value -  20.

    * –

commensurate_basic_delay: value of the basic delay $h$ in case of commensurate
delays. May lead to a speed-up for computing N if the ratio $\tau_{m}/h$ is
small.

    * –

quiet: suppress all output (including warnings).
Default value -  false.

####  B.1.2 Warnings

The function tds_roots may give the following warnings.

  * •

The provided time-delay system is an ODE/a DAE. The characteristic function of
the provided system does not contain delays. falls back to eig for computing
the characteristic roots. Note that only the characteristic roots in the
desired region are returned in L.

  * •

Size of the discretized EVP would exceed its maximum value. The degree of the
spectral discretisation needed to accurately capture all characteristic roots
in the specified region (N_original) would result in an eigenvalue problem
(B.1) that is larger than options.max_size_evp. The degree of the spectral
discretisation is therefore lowered to N.

  * •

Case: gamma(r) >= 1 (i.e., CD>r). This warning is thrown when region is a
scalar (one is interested in the characteristic roots in a given right half-
plane) and $\gamma(\text{ {region} })\geq 1$. This means that the desired
right half-plane contains infinitely many characteristic roots or that there
exist infinitesimal perturbations on the delays such that it does. As a
consequence, the default value $N=15$ is used (lowered if $N(n+1)$ exceeds
options.max_size_evp, see warning above). This means that not all desired
characteristic roots may be captured. The user can manually increase $N$ using
the option

###  B.2 tds_gamma_r

The function tds_gamma_r computes the quantity $\gamma(r)$ (as defined in
(2.23) for neutral DDEs and in (2.30) for DDAEs). As mentioned in Proposition
2.4, the quantity $\gamma(0)$ gives a necessary condition for strong stability
(i.e., stability is preserved for sufficiently small delay perturbations). On
the other hand, the inequality $\gamma(r)<1$ guarantees that the right half-
plane $\\{\lambda\in\mathbb{C}:\Re(\lambda)>r\\}$ only contains finitely many
characteristic roots (even when infinitesimal delay perturbations are
considered). To compute the quantity $\gamma(r)$ for a NDDE, the code will
first obtain a prediction for the solution of the optimization problem in
(2.23) by restricting the search space for $\vec{\theta}$ to a grid in
$\\{0\\}\times[0,2\pi)^{m_{H}-1}$ (note that $\theta_{1}$ can be fixed to
zero). The resulting value is subsequently improved by solving the following
overdetermined system of nonlinear equations in the unknowns
$v\in\mathbb{C}^{n}$, $u\in\mathbb{C}^{n}$, $\lambda\in\mathbb{C}$, and
$\theta\in\\{0\\}\times[0,2\pi)^{m_{H}-1}$

$\left\\{\begin{array}[]{l}\left(H_{1}e^{-rh_{H,1}}+\sum_{k=2}^{m_{H}}H_{k}e^{-rh_{H,k}}e^{\jmath\theta_{k}}\right)v-\lambda
v=0\\\
u^{H}\left(H_{1}e^{-rh_{H,1}}+\sum_{k=2}^{m_{H}}H_{k}e^{-rh_{H,k}}e^{\jmath\theta_{k}}\right)-u^{H}\lambda=0\\\
n(v)-1=0\\\ u^{H}v-1=0\\\
\Im\left(\bar{\lambda}(u^{H}H_{k}v)e^{-rh_{H,k}}e^{\jmath\theta_{k}}\right)=0,\text{
for }k=2,\dots,m_{H}\end{array}\right.$ (B.2)

with $n(v)-1=0$ a normalization constraint, using the fsolve function in
Matlab. A solution $(v,u,\lambda,\theta)$ of (B.2) corresponds to an
eigenvalue $\lambda$ of the matrix

$H_{1}e^{-rh_{H,1}}+\sum_{k=2}^{m_{H}}H_{k}e^{-rh_{H,k}}e^{\jmath\theta_{k}}$

with $u$ and $v$, respectively, the corresponding left and right eigenvectors
such that the derivatives of the modulus of this eigenvalue with respect to
$\theta_{2},\dots,\theta_{m_{H}}$ equals zero. The approach for computing the
quantity $\gamma(r)$ of a DDAE is similar.

####  B.2.1 Documentation

  * •

gammar = tds_gamma_r(tds,r) computes the quantity $\gamma(\texttt{r})$ of tds.

  * •

gammar = tds_gamma_r(tds,r,options) allows to specify additional options.
options must be a structured array created using the function
tds_gamma_r_options. The available options are:

    * –

quiet: suppress all output (including warnings).
Default value -  false.

    * –

Ntheta: the number of grid points in each direction.
Default value -  10.

###  B.3 tds_CD

The function tds_CD computes the strong spectral abscissa of the delay
difference equation underlying a delay differential equation. To compute this
quantity we use the expression for $C_{D}$ in Proposition 2.5. More
specifically, we use fsolve to find a zero crossing of the function (2.22).

####  B.3.1 Documentation

  * •

CD = tds_CD(tds) computes the strong spectral abscissa of the delay difference
equation associated with tds.

  * •

CD = tds_CD(tds,options) allows to specify additional options. options must be
a structured array created using the function tds_CD_options. The available
options are:

    * –

tol: tolerance
Default value -  1e-10

    * –

quiet: suppress all output (including warnings)
Default value -  false

    * –

Ntheta: the number of grid points in each direction
Default value -  10

    * –

CD0: initial estimate for the strong spectral abscissa of the underlying delay
difference equation
Default value -  0

###  B.4 tds_sa

The function tds_sa computes the spectral abscissa (see Definitions 2.2 and
2.5) of a time-delay system. To this ends it computes all characteristic roots
in a given right half-plane using tds_roots and then returns the real part of
the rightmost root.

####  B.4.1 Documentation

  * •

C = tds_sa(tds,r) computes the spectral abscissa of tds by computing its
characteristic roots in the right half plane $\\{z\in\mathbb{C}:\Re(z)\geq
R\\}$.

  * •

C = tds_sa(tds,r,options) computes the spectral abscissa with the default
options for the tds_roots replaced by options.

###  B.5 tds_strong_sa

The function tds_strong_sa computes the strong spectral abscissa of a time-
delay system. Recall from Definition 2.7 that the strong spectral abscissa is
defined as the smallest upper bound for the spectral abscissa which is
insensitive to infinitesimal delay changes. Furthermore, remember that for
a(n) (essentially) retarded time-delay system the strong spectral abscissa
corresponds to the spectral abscissa and that for a(n) (essentially) neutral
time-delay system the strong spectral abscissa is equal to the maximum of the
spectral abscissa and the strong spectral abscissa of the associated
difference equation (which can be computed using the function tds_CD). To
compute the strong spectral abscissa, the code will therefore first determine
the strong spectral abscissa of the associated delay difference equation (CD).
Subsequently, it checks for characteristic roots in the right half-plane
$\\{z\in\mathbb{C}:\Re(z)>\texttt{CD}\\}$.

####  B.5.1 Documentation

  * •

C = tds_strong_sa(tds,r) computes the strong spectral abscissa of tds by first
computing CD = tds_CD(tds) and subsequently computing the characteristic roots
of tds in the right half-plane
$\\{z\in\mathbb{C}:\Re(z)>\max(\texttt{CD},\texttt{r})\\}$.

  * •

C = tds_strong_sa(tds,r,options) allows to specify additional options. The
argument options should be a structure with the necessary fields. The
available options are:

    * –

roots: a structure created using tds_roots_options containing the options for
tds_roots.

    * –

CD - a structure created using tds_CD_options containing the options for
tds_CD.

###  B.6 tds_tzeros

The function tds_tzeros allows to compute the transmission zeros of a SISO
time-delay system, i.e., a point $z_{0}$ in the complex plane for which the
transfer function becomes zero. To this end, the function exploits the fact
that computing the transmission zeros of a SISO transfer function can be
reformulated as computing the characteristic roots of a DDAE (see Section
2.6.1).

####  B.6.1 Documentation

  * •

tzeros = tds_tzeros(tds,rect) returns a vector tzeros containing the
transmission zeros of the time-delay system tds in the rectangular region
[rect(1) rect(2)] $\times\jmath$ [rect(3) rect(4)].

  * •

tzeros = tds_tzeros(tds,rect,options) allows to specify additional options
created using the function tds_roots_options.

###  B.7 tds_stabopt_static and tds_stabopt_dynamic

The functions tds_stabopt_static and tds_stabopt_dynamic look for (strongly)
stabilizing output feedback controller for a given time-delay system by
minimizing a certain objective function (see below) in the entries of the
controller matrices, i.e., the matrices in (3.3) and (3.2). Below, we will
briefly discuss the core working principles of these functions, for more
details see [35, 39, 24]. For ease of notation, we will restrict our attention
to the dynamic feedback control of the following retarded time-delay system

$\left\\{\begin{array}[]{rcl}\dot{x}(t)&=&A_{0}\,x(t)+\sum_{k=1}^{m}A_{k}x(t-\tau_{k})+B_{1,0}\,u(t)+\sum_{k=1}^{m}B_{1,k}\,u(t-\tau_{k}),\\\\[9.0pt]
y(t)&=&C_{1,0}\,x(t)+\sum_{k=1}^{m}C_{1,k}\,x(t-\tau_{k})+D_{11,0}\,u(t)+\sum_{k=1}^{m}D_{11,k}\,u(t-\tau_{k}).\end{array}\right.$

The discussion for static control or neutral and delay descriptor systems is
however similar. The dynamics of the feedback interconnection of the system
above with the (3.2) can be described by

$\begin{bmatrix}I&0&0&0\\\ 0&0&0&0\\\ 0&0&I&0\\\
0&0&0&0\end{bmatrix}\begin{bmatrix}\dot{x}(t)\\\ \dot{y}(t)\\\
\dot{x}_{c}(t)\\\
\dot{u}_{u}(t)\end{bmatrix}=\begin{bmatrix}A_{0}&0&0&B_{1,0}\\\
C_{1,0}&-I&0&D_{11,0}\\\ 0&0&0&0\\\
0&0&0&-I\end{bmatrix}\begin{bmatrix}x(t)\\\ y(t)\\\ x_{c}(t)\\\
u(t)\end{bmatrix}\ +\\\\[5.01874pt]
\sum_{k=1}^{m}\begin{bmatrix}A_{k}&0&0&B_{1,k}\\\ C_{1,k}&0&0&D_{11,k}\\\
0&0&0&0\\\ 0&0&0&0\end{bmatrix}\begin{bmatrix}x(t-\tau_{k})\\\
y(t-\tau_{k})\\\ x_{c}(t-\tau_{k})\\\
u(t-\tau_{k})\end{bmatrix}+\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\
0&B_{c}&A_{c}&0\\\ 0&D_{c}&C_{c}&0\\\ \end{bmatrix}\begin{bmatrix}x(t)\\\
y(t)\\\ x_{c}(t)\\\ u(t)\end{bmatrix}.$ (B.3)

In other words, the closed-loop system corresponds to a delay differential
algebraic equation of the form

$E_{cl}\,\dot{x}_{cl}(t)=\sum_{k=1}^{m_{cl}}A_{cl,k}(\mathbf{p})\
x_{cl}(t-\tau_{cl,k})$ (B.4)

in which $\mathbf{p}$ denotes the vector containing the free control
parameters (i.e., the optimization variables). In the case of an unstructured
controller design problem, $\mathbf{p}$ consist of the entries of the matrices
$A_{c}$, $B_{c}$, $C_{c}$ and $D_{c}$. For the design of a structured
controller (see Section 3.3) only certain elements of these matrices are
contained in $\mathbf{p}$. Note that the state matrices
$A_{cl,1},\dots,A_{cl,m}$ depend affinely on the controller parameters in
$\mathbf{p}$.

Next, we will give more information on the objective functions that are
minimized. If (B.4) is essentially retarded (see Section 2.4),
tds_stabopt_dynamic attempts to find a stabilizing controller by minimizing
the spectral abscissa, or in other words, the objective function is given by

$f(\mathbf{p}):=c\left(E_{cl},A_{1,cl}(\mathbf{p}),\dots,A_{m_{cl},cl}(\mathbf{p})\right).$

If (B.4) is essentially neutral, closed-loop stability might be sensitive to
infinitesimal delay perturbations (see Section 2.3.2). To guarantee strong
exponential stability (as defined in Definition 2.6), we can consider either
of the following two approaches (see also Section 3.2).

  * •

Approach 1: minimize the strong spectral abscissa of (B.4), or in other words,
the objective function is given by

$f(\mathbf{p}):=C\left(E_{cl},A_{0,cl}(\mathbf{p}),\dots,A_{k,cl}(\mathbf{p})\right).$

  * •

Approach 2: minimize the spectral abscissa of (B.4) under the constraint that
the inequality

$\gamma_{0}(A_{0,cl}^{(22)}(\mathbf{p}),\dots,A_{m_{cl},cl}^{(22)}(\mathbf{p}))<1$

with $A_{0,cl}^{(22)},\dots,A_{m_{cl},cl}^{(22)}$ defined in a similar vain as
in Section 2.4, holds. To solve this constraint optimization problem, uses a
log-barrier method, i.e., a barrier function which grows to infinity as
$\gamma(0)$ goes 1, is added to the cost function. More specifically,
considers the following barrier function

$-w_{2}\log\big{(}1-w_{1}-\gamma_{0}(A_{0,cl}^{(22)}(\mathbf{p}),\dots,A_{m_{cl},cl}^{(22)}(\mathbf{p}))\big{)}$
(B.5)

in which $0\leq w_{1}<1$ dictates the acceptable region for $\gamma_{0}$ and
$w_{2}>0$ determines the weight of the barrier function in the overall cost
function, i.e., the larger $w_{2}$, the more emphasis is put on keeping
$\gamma_{0}$ small. Note that in order to use this approach, an initial
controller for which the inequality $\gamma_{0}<1-w_{1}$ holds, is required.
implements therefore an initial optimization step in which
$\gamma_{0}(E_{cl},A_{0,cl}(\mathbf{p}),\dots,A_{k,cl}(\mathbf{p}))$ is
minimized with respect to the control parameters. Once a controller found for
which $\gamma_{0}<1-w_{1}$ the following extended cost function is minimized
in $\mathbf{p}$:

$f(\mathbf{p}):=c\left(E_{cl},A_{0,cl}(\mathbf{p}),\dots,A_{k,cl}(\mathbf{p})\right)\\\
-w_{2}\log\big{(}1-w_{1}-\gamma_{0}(A_{0,cl}^{(22)}(\mathbf{p}),\dots,A_{m_{cl},cl}^{(22)}(\mathbf{p}))\big{)}.$

In each case thus solves an unconstrained optimization problem of the form

$\min_{\mathbf{p}}f(\mathbf{p}).$ (B.6)

The objective function $f$ is typically neither convex nor smooth. (However
the object function is typically almost everywhere differentiable.) The
therefore relies on HANSO v3.0 [23, 7] for solving these nonsmooth, nonconvex
optimization problems.

Important. To avoid ending up in a local optimum of (B.6), the optimization
routine is restarted from different initial values for $\mathbf{p}$.

####  B.7.1 Documentation

  * •

K = tds_stabopt_static(tds) attempts to design a strongly stablizing static
output feedback gain matrix K for tds.

  * •

K = tds_stabopt_static(tds,’key1’,value1,’key2’,value2,...) allows to pass the
following optional arguments using key-value pairs:

    * –

options structured array created using tds_stabopt_options containing
additional options

    * –

initial cell array containing the initial optimization values

    * –

mask logical array indicating which entries of the feedback gain matrix are
allowed to change (true) and which should remain fixed to their original value
(false)

    * –

basis array specifying the values of the fixed entries of the feedback gain
matrix

  * •

[K,CL] = tds_stabopt_static(tds,...) also returns a tds_ss_ddae-object that
represents the closed-loop system.

  * •

K = tds_stabopt_dynamic(tds,nc) attempts to design a strongly stablizing
dynamic output feedback controller K of order nc for tds.

  * •

K = tds_stabopt_dynamic(tds,nc,’key1’,value1,’key2’,value2,...) allows to pass
the following optional arguments using key-value pairs:

    * –

options structured array created using tds_stabopt_options containing
additional options

    * –

initial cell array containing the initial optimization values

    * –

mask an ss-object indicating which entries of the feedback gain matrix are
allowed to change (1) and which should remain fixed to their original value
(0)

    * –

basis an ss-object specifying the fixed entries of the feedback gain matrix

  * •

[K,CL] = tds_stabopt_dynamic(tds,nc,...) also returns a tds_ss_ddae-object
that represents the closed-loop system.

####  B.7.2 Options

The function tds_stabopt_options can be used to create a structured array
specifying additional options for tds_stabopt_static and tds_stabopt_dynamic.

  * •

nstart: number of starting points for the optimization procedure
Default value:  5.

  * •

roots: options for tds_roots (see Section B.1)

  * •

Ntheta: number of grid points in each direction for computing $\gamma_{0}$
(Approach 2, see Section B.2) or $C_{D}$ (Approach 1, see Section B.3)
Default value -  10.

  * •

cpumax: quit the BFGS-procedure if the CPU time (in seconds) exceeds the given
value
Default value:  inf.

  * •

fvalquit: quit the BFGS-procedure if the object values drops below the given
value
Default value:  -inf.

  * •

gradient_sampling: apply the gradient sampling algorithm from [23] after the
BFGS-procedure. By default, gradient sampling is disabled as it is typically
quite slow.
Default value:  false.

  * •

method: for essentially neutral closed-loop systems select which cost function
is used:

    * –
    * –
    * –

Default value:

  * •

Default value:  1e-3.

  * •

Default value:  1e-3.

  * •

print_level: 0 (no printing), 1 (minimal printing), 2 (verbose), 3 (debug)
Default value:  2.

####  B.7.3 Warnings

The functions tds_stabopt_static and tds_stabopt_dynamic may give the
following warnings.

  * •

Warning: Resulting controller is not stabilizing. This warning indicates that
tds_stabopt_static or tds_stabopt_dynamic failed to find a (strongly)
stabilizing controller. Note that this warning does not imply that there does
not exist a (strongly) stabilizing controller. It merely indicates that
starting from the chosen initial values, the optimization procedure ended up
in local optima that do not correspond with a strictly negative (strong)
spectral abscissa. A possible solution is to sample more initial values using
the option

###  B.8 tds_hinfnorm

The function tds_hinfnorm computes the (strong) $\mathcal{H}_{\infty}$-norm of
a time-delay system. To this end, it implements the method from [16]. This
method can be summerized as follows:

  1. 1.

Compute the strong $\mathcal{H}_{\infty}$-norm of the associated asymptotic
transfer function by solving the optimization problem in (4.14).

     * •

Prediction step. As (4.14) is generally not convex and thus may have multiple
local optima, the object function is first sampled on a grid in the
$\theta$-space with

     * •

Correction step. Improve the obtained result using Newton’s method.

  2. 2.

Verify whether the maximal singular value function attains a higher value at a
finite frequency.

     * •

Prediction step. Form a descriptor system that approximates the time-delay
system using a spectral discretisation of degree Note. In contrast to

     * •

Correction step. Improve the obtained result using Newton’s method.

Important: If the provided time-delay system is not strongly exponentially
stable, the function tds_hinfnorm will return infinity.

####  B.8.1 Documentation

  * •

[HINF] = TDS_HINFNORM(TDS) returns the strong $\mathcal{H}_{\infty}$-norm of
TDS.
Note. HINF=Inf when TDS is not strongly exponentially stable.

  * •

[HINF,WPEAK] = tds_hinfnorm(TDS) also returns the angular frequency at which
the H-infinity norm is attained. WPEAK=Inf when the strong H-infinity norm of
TDS is equal to the strong H-infinity norm of the underlying asymptotic
transfer function. WPEAK=NaN when TDS is not strongly exponentially stable.

  * •

[..] = tds_hinfnorm(TDS,OPTIONS) allows to specify additional options (see
below).

####  B.8.2 Options

The function

  * •

Default value:  20.

  * •

Default value:  20

  * •

Default value:

  * •

Default value:  1e-3

  * •

Default value:  1e-6

####  B.8.3 Warnings

  * •
  * •

###  B.9 tds_sigma

The function tds_sigma computes the singular values of the frequency response
matrix of a time-delay system at the given angular frequencies.

####  B.9.1 Documentation

  * •

SV = tds_sigma(TDS,W) returns the singular values SV of the frequency response
of TDS at the provided angular frequency vector W. The matrix SV has length(W)
columns and the column SV(:,k) gives the singular values (in descending order)
of the frequency response of TDS at the frequency W(k).

###  B.10 tds_hiopt_static and tds_hiopt_dynamic

The functions tds_hiopt_static and tds_hiopt_dynamic tune static and dynamic
output feedback controllers for time-delay systems by minimizing the closed-
loop $\mathcal{H}_{\infty}$-norm with respect to the controller parameters.
Its working principles are quite similar to those of the functions
tds_stabopt_static and tds_stabopt_dynamic (see Section B.7). More
specifically, its core working principles can be summarized as follows:

  1. 1.

Create a delay descriptor representation for the closed-loop system (depicted
in Figure 4.3).

  2. 2.

Take the user-provided initial values for the controller parameters or
randomly generate them.

  3. 3.

Check whether the closed-loop system is strongly exponentially stable for the
initial controller values. If not run tds_stabopt_static or
tds_stabopt_dynamic with the entries of this controller as initial values to
obtain a stabilizing starting point. If the result is not stabilizing, discard
this initial values.

  4. 4.

For all stable starting points, minimize the closed-loop
$\mathcal{H}_{\infty}$-norm with respect to the controller parameters using a
BFGS algorithm with weak line search.

  5. 5.

Optional: improve the obtained optimum using gradient sampling.

For more details we refer the interested reader to [16].

Note. The relies on HANSO v3.0 for the last two steps in the description
above.

####  B.10.1 Documentation

  * •

K = tds_hiopt_static(tds) designs a strongly stablizing static output feedback
gain matrix K for tds by minimizing the closed-loop
$\mathcal{H}_{\infty}$-norm with respect to the entries of K.

  * •

K = tds_hiopt_static(tds,’key1’,value1,’key2’,value2,...) allows to pass the
following optional arguments using key-value pairs:

    * –

options structured array created using tds_hiopt_options containing additional
options

    * –

initial cell array containing the initial optimization values

    * –

mask logical array indicating which entries of the feedback gain matrix are
allowed to change (true) and which should remain fixed to their original value
(false)

    * –

basis array specifying the values of the fixed entries of the feedback gain
matrix

  * •

[K,CL] = tds_hiopt_static(tds,...) also returns a tds_ss_ddae-object that
corresponds to the resulting closed-loop system.

  * •

K = tds_hiopt_dynamic(tds,nc) designs a strongly stablizing dynamic output
feedback controller K of order nc for tds by minimizing the closed-loop
$\mathcal{H}_{\infty}$-norm with respect to the controller parameters.

  * •

K = tds_hiopt_dynamic(tds,nc,’key1’,value1,’key2’,value2,...) allows to pass
the following optional arguments using key-value pairs:

    * –

options structured array created using tds_hiopt_options containing additional
options

    * –

initial cell array containing the initial optimization values

    * –

mask an ss-object indicating which entries of the feedback gain matrix are
allowed to change (1) and which should remain fixed to their original value
(0)

    * –

basis an ss-object specifying the fixed entries of the feedback gain matrix

  * •

[K,CL] = tds_hiopt_dynamic(tds,nc,...) also returns a tds_ss_ddae-object that
corresponds to the resulting closed-loop system.

####  B.10.2 Options

The function tds_hiopt_options can be used to create a structured array
specifying additional options for tds_hiopt_static and tds_hiopt_dynamic.

  * •

nstart: number of starting points for the optimization procedure
Default value:  5.

  * •

cpumax: quit the BFGS-procedure if the CPU time (in seconds) exceeds the given
value
Default value:  inf.

  * •

fvalquit: quit the BFGS-procedure if the object values drops below the given
value
Default value:  -inf.

  * •

gradient_sampling: apply the gradient sampling algorithm from [23] after the
BFGS-procedure. By default, gradient sampling is disabled as it is typically
quite slow.
Default value:  false.

  * •

print_level: 0 (no printing), 1 (minimal printing), 2 (verbose), 3 (debug)
Default value:  2.

  * •

alpha: the value of $\alpha$ in (4.22)
Default value:  0.

  * •

roots: options for tds_roots (see Section B.1)

####  B.10.3 Warnings and errors

The functions tds_hiopt_static and tds_hiopt_dynamic may give the following
warnings and errors.

  * •

Failed to find a strongly stabilizing initial controller. As the strong
$\mathcal{H}_{\infty}$-norm of a system that is not strongly exponentially
stable is equal to $\infty$, these two functions check for each starting point
whether the resulting closed-loop system satisfies this stability requirement.
If this is not the case, they will run tds_stabopt_static or
tds_stabopt_dynamic starting the given starting point, to find a stabilizing
initial controller. If for none of the starting points, a strongly stabilizing
initial controller can be found, this error is thrown.

###  B.11 tds_create_delta

The function tds_create_delta can be used to create a structured array that
stores the relevant information of a set of matrix perturbations $\delta$, as
defined in (4.29), such as the dimensions, whether the perturbations are real
or complex-valued and the norm in which the size of the perturbations should
be measured.

####  B.11.1 Documentation

  * •

delta = tds_create_delta(dim1,dim2,type,norm) with dim1, dim2, type and norm
cell arrays of equal length. The elements of dim1 and dim2 give respectively
the first and second dimensions of the different perturbations. The elements
of type indicate which the perturbations are real (delta is a structured array
storing the relevant information of the set of matrix perturbations.

###  B.12 tds_uncertain_matrix and tds_uncertain_delay

The function tds_uncertain_matrix creates a structured array that stores all
relevant information related to an uncertain matrix of the form:

$\tilde{A}_{k}(\delta)=\textstyle\sum\limits_{j=1}^{N_{A_{k}}}G_{k,j}\
\delta_{\texttt{ind}_{A_{k}}[j]}\ H_{k,j}.$ (B.7)

The function tds_uncertain_delay creates a structured array that stores all
relevant information related to an uncertain delay of the form:

$\tilde{h}_{A,k}(\delta\tau)=\max\left\\{-h_{A,k},\textstyle\sum\limits_{j=1}^{N_{\tau_{k}}}w_{k,j}\,\delta\tau_{\texttt{ind}_{h_{A,k}}[j]}\right\\}.$
(B.8)

####  B.12.1 Documentation

  * •

umat = tds_uncertain_matrix(ind,G,H) creates a structured array umat that
stores the relevant information of an uncertain matrix of the form (B.7) with
ind a cell array storing the positions of the perturbations in the
corresponding set of matrix perturbations (created using tds_create_delta),
and G and H cell arrays containing the shape matrices.

  * •

udelay = tds_uncertain_delay(ind,w) creates a structured array udelay that
stores the relevant information of an uncertain delay of the form (B.8) with
ind a cell array storing the positions of the perturbations in the
corresponding set of delay perturbations and w a cell array containing the
weights.

###  B.13 tds_psa

The function tds_psa computes the $\epsilon$-pseudospectral abscissa of the
uncertain retarded delay differential equation:

$\dot{x}(t)=\big{(}A_{0}+\tilde{A}_{0}(\delta)\big{)}\,x(t)+\sum_{k=1}^{m}\big{(}A_{k}+\tilde{A}_{k}(\delta)\big{)}\,x\big{(}t-(h_{A,k}+\tilde{h}_{A,k}(\delta\tau))\big{)}.$
(B.9)

Recall that the $\epsilon$-pseudospectral abscissa was defined as the worst-
case value of the spectral abscissa over all allowable uncertainty values,
i.e.,

$c^{\mathrm{ps}}(\epsilon):=\max\\{c(\delta,\delta\tau):\|\delta_{l}\|\leq\epsilon\text{
for }l=1,\dots,L_{M}\ \&\ |\delta\tau_{l}|\leq\epsilon\text{ for
}l=1,\dots,L_{D}\\}$ (B.10)

with $c(\delta,\delta\tau)$ the spectral abscissa of (B.9) for a given
instance of the uncertainties.

####  B.13.1 Documentation

  * •

psa = tds_psa(utds,epsilon) returns the epsilon-pseudospectral abscissa, as
defined in (B.10), of the uncertain retarded DDE represented by utds.

###  B.14 tds_dins

The function tds_dins computes the distance to instability of the uncertain
RDDE (B.9). Recall that the distance to instability was defined as the
smallest $\epsilon$ for which the $\epsilon$-pseudospectral abscissa is
positive, i.e.,

$d_{\mathrm{INS}}=\min\\{\epsilon\geq 0:c^{ps}(\epsilon)\geq 0\\}.$ (B.11)

####  B.14.1 Documentation

  * •

dins = tds_dins(utds,interval) computes the distance to instability, as
defined in (B.11), of utds by finding the zero-crossing of the function
$\epsilon\mapsto c^{\mathrm{ps}}(\epsilon)$ in the search interval specified
by interval.