On Norm-Based Estimations for Domains of Attraction in Nonlinear Time-Delay Systems

TL;DR

This paper introduces a numerical method to estimate the domain of attraction in nonlinear time-delay systems, combining initial function selection, simulations, and bifurcation analysis for improved bounds.

Contribution

It presents a novel methodology for estimating attraction domains in nonlinear time-delay systems using a combination of numerical simulations and bifurcation analysis.

Findings

Numerical upper bounds on attraction domains can be refined iteratively.

Bifurcation analysis improves the estimation of the domain of attraction.

Application to a time-delayed swing equation demonstrates the method's effectiveness.

Abstract

For nonlinear time-delay systems, domains of attraction are rarely studied despite their importance for technological applications. The present paper provides methodological hints for the determination of an upper bound on the radius of attraction by numerical means. Thereby, the respective Banach space for initial functions has to be selected and primary initial functions have to be chosen. The latter are used in time-forward simulations to determine a first upper bound on the radius of attraction. Thereafter, this upper bound is refined by secondary initial functions, which result a posteriori from the preceding simulations. Additionally, a bifurcation analysis should be undertaken. This analysis results in a possible improvement of the previous estimation. An example of a time-delayed swing equation demonstrates the various aspects.

Equations

$$x_t(\theta )\stackrel{\text{def}}{=}x(t+\theta ),\theta \in [-\tau ,0],$$

$${\gamma }^{+}(\varphi )\stackrel{\text{def}}{=}\{x_t(\cdot ;\varphi )\in X:t\in [0,t_{max})\}$$

$$\dot{x}(t)$$

$$x(\theta )$$

$$F(x_t)=f(x(t),x(t-\tau )),$$

$$\forall \epsilon >0,\exists \delta (\epsilon )>0:\Vert \varphi {\Vert }_X\le \delta \Rightarrow \Vert x_t{\Vert }_X\le \epsilon ,\forall t\ge 0,$$

$$\exists {\delta }_a>0:\Vert \varphi {\Vert }_X\le {\delta }_a\Rightarrow \underset{t\to \infty }{\lim }\Vert x_t{\Vert }_X=0.$$

$${\mathcal{D}}_X=\{\varphi \in X:x_t\text{ exists on }t\in [0,\infty )\text{ and }\underset{t\to \infty }{\lim }\Vert x_t{\Vert }_X=0\}.$$

$${\mathcal{B}}_X(r)$$

$$\text{with }{\mathcal{B}}_X({\varphi }_0;r)$$

$$R_X=\sup \{r>0:{\mathcal{B}}_X(r)\subseteq {\mathcal{D}}_X\}$$

$${\hat{R}}_X(\Phi )=\underset{\varphi \in \Phi }{\inf }\{\Vert \varphi {\Vert }_X:x(t;\varphi )\nrightarrow 0\}.$$

$${\hat{R}}_X(\Phi )=\inf \varnothing =\infty .$$

$$\dot{x}(t)$$

$$x(\theta )$$

$$\frac{\text{d}\hat{x}}{\text{d}\hat{t}}(\hat{t})$$

$$\hat{x}(\hat{\theta })$$

$$\Vert \hat{\varphi }{\Vert }_{C[-1,0]}=\underset{\hat{\theta }\in [-1,0]}{\sup }\varphi (\tau \hat{\theta })=\Vert \varphi {\Vert }_{C[-\tau ,0]}$$

$$\Vert \hat{\varphi }{\Vert }_{M^2[-1,0]}$$

$$=\Vert \varphi (0){\Vert }_2+\sqrt{\frac{1}{\tau }}{\left({\int }_{-\tau }^0\Vert \varphi (\theta ){\Vert }_2^2\text{d}\theta \right)}^{\frac{1}{2}}$$

$$\ne \Vert \varphi {\Vert }_{M^2[-\tau ,0]}.$$

$$\left[\begin{array}{c}{\dot{x}}^{\text{I}}(t)\\ {\dot{x}}^{\text{II}}(t)\end{array}\right]$$

$$\left[\begin{array}{c}{x}^{\text{I}}(\theta )\\ x^{\text{II}}(\theta )\end{array}\right]$$

$$\displaystyle\|\phi\|_{Q}:=\left\|\begin{array}[]{c}\|\phi^{\textrm{\scriptsize I}}\|_{X}\\ \|\phi^{\textrm{\scriptsize II}}(0)\|\end{array}\right\|$$

$$Q:=X/\{\varphi \in X:{\varphi }^{\text{I}}(\theta )\equiv 0_{n-q},{\varphi }^{\text{II}}(0)=0_q\}$$

$${\hat{\mathcal{D}}}_X\stackrel{\text{def}}{=}\{\varphi \in X:r(t_{\text{num}};\varphi )<{\delta }_{\text{num}}\},$$

$${\Phi }^k\stackrel{\text{def}}{=}\left\{{\varphi }^k(\cdot ;p)\in X:p\in {\mathbb{R}}^m\right\}$$

$$\Phi \stackrel{\text{def}}{=}\bigcup _k{\Phi }^k.$$

$$\displaystyle y^{(n)}(t)+g\big{(}$$

$$\displaystyle y^{(n-1)}(t-\tau),\ldots,\dot{y}(t-\tau),y(t-\tau)\big{)}=0.$$

$$\varphi (\theta ;p)=[y_0(\theta ;p),\dots ,y_0^{(n-1)}(\theta ;p){]}^{\top }.$$

$$\underset{t\to \infty }{\lim }\Vert x_t(\cdot ;\varphi (\theta ;p)){\Vert }_X=$$

$${\hat{R}}_X^{\mathcal{T}}(\Phi )\stackrel{\text{def}}{=}\underset{t\ge 0}{\min }\{\Vert x_t(\cdot ;\varphi ){\Vert }_X:x(t;\varphi )\nrightarrow 0,\varphi \in \Phi \},$$

$$\varphi (\theta ;p)=\sum _{i=1}^d{p}_i{b}_i(\theta )$$

Full text

1 Introduction

Time delays due to communication, measurement, data processing, delayed actuator reactions, or transport processes are omnipresent in technical systems. That is why stability analysis of time-delay systems has gathered increasing interest in the last decades. However, the focus was mainly on stability criteria for linear systems [1, 2, 3, 4, 5], be it Lyapunov-Krasovskii functionals [6, 7, 8], Lyapunov-Razumikhin functions [6], comparison principles [9, 10, 11], input-output approaches [12], or eigenvalue calculations [13, 14, 15, 16]. By contrast, in technical applications, we encounter nonlinear systems [17, 18, 19, 20, 21]. Indeed, for a nonlinear system the Principle of Linearized Stability [22] allows to deduce stability or instability of equilibria, but the question arises: what is the practical relevance of knowledge about asymptotic stability if there is no knowledge about the domain of attraction? In this light, small initial disturbances can still result in departing trajectories. In relation to the results listed above, estimations on domains of attraction in time-delay systems are rare in literature [23, 24, 25, 26, 27]. Time-delay systems require an initial function instead of the initial value, which would be sufficient in ODEs. Hence, the state space, which contains the domain of attraction, is infinite dimensional. The present paper addresses such domains of attraction for equilibria in autonomous retarded functional differential equations (RFDE) with a focus on systems with a discrete time delay.

To this end, consider the spectral abscissa, which in an eigenvalue-based stability analysis of the linearized system decides about asymptotic stability and instability (Principle of Linearized Stability, Diekmann [22], Ch. VII, Theorem 6.8). It should be noted that this decisive number has neither an implication to allowed perturbations of the RFDE right-hand side or system parameters, i.e. robustness, nor to allowed initial perturbations, i.e. the domain of attraction [28]. However, in practical applications various requirements concerning the domain of attraction occur and necessitate adequate approaches.

(R1) (test) Let a certain initial function be given. Is this initial function an element of the domain of attraction? (e.g. contingency tests in power systems)

(R2) (prove) Let a bound for the initial functions in a certain norm be given. Do all initial perturbations within this bound lead to attractor-convergent solutions? (e.g. requirements in controller design)

(R3) (disprove) Proving non-fulfillment of (R2).

(R4) (compare) Let different systems or controller configurations be given. Which of the systems is the better one w.r.t. the domain of attraction in the sense of some order relation?

(R5) (parametrize) Let a system or controller with free parameters be given. How is the domain of attraction influenced by the parameters? (e.g. system analysis and controller synthesis)

Different approaches are possible in order to meet the requirements listed above.

(i) A single time-domain simulation for each specific initial function of interest is able to fulfill requirement (R1) in a hardly conservative way. However, in contrast to norm-based criteria, which - once formulated - allow to conclude convergence directly, the computational effort of such individual simulations should be kept in mind. This is for instance the reason why in (delay-free) power system stability analysis, methods have been introduced, which do no longer rely only on on-demand time-domain simulations [29, 30].

(ii) Numerical Basin stability [31] in its original form for time-delay-free systems approximates the Lebesgue measure of the finite dimensional domain of attraction (blue shaded area in Figure 1).

Hence, there is a relation to the probability of randomly chosen initial values to be within the domain of attraction. Leng et al. [32] give a generalization to time-delay systems. The derived number can contribute to a comparison of different system configurations as required in (R4).

(iii) The level set of a Lyapunov-Krasovskii functional [6, 23] can be used as estimation for the domain of attraction (represented by the inner of $V=V_{crit}$ in Figure 1). Consequently, it yields a scalar criterion in order to test whether an initial function surely belongs to the domain of attraction (R1). The background is a generalization of LaSalle’s invariance principle to RFDEs [33]. However, already for linear systems Lyapunov-Krasovskii functionals are challenging. Vast efforts have been made in view of the possible delay dependence of stability in such linear RFDEs. The challenge is then to derive functionals that conclude stability even for time delay values near the bound at which stability is lost. Nonlinear systems necessitate individual approaches. Alternatively, instead of searching for such special functionals, there is also the ansatz to use a functional that is suitable for the linearized system. The latter approach is well-known for ordinary differential equations (ODE) [34]. However, it has to be reckoned with conservative estimations. Lyapunov-Krasovskii functionals are in general complicated expressions, such that the raw sublevel set criterion itself is not very intuitive. Therefore, it is usually translated to a norm criterion, see (vi).

(iv) The level set of a Lyapunov-Razumikhin function [35] can be used as well, since there are also Razumikhin-based generalizations of LaSalle’s invariance principle [36].

(v) Bounds in a finite dimensional parameter space [37], by analytical or numerical methods, come into play if only a certain class of parametrizable initial functions is relevant. For instance, a jump from a predefined initialization to a variable value is a class of initial functions that is fully determined by a single parameter. Then, the space of initial functions is indeed finite dimensional and the finite number of parameters is decisive.

(vi) An underestimation of the radius of attraction yields a norm criterion that ensures convergence to the attractor (interior of the green dashed circle with radius $\check{R}_{X}$ in Figure 1). It is essential in (R2) and yields a very simple criterion for testing in (R1). Of course, usually nothing is known about the conservativeness of such criteria. Results based on Lyapunov-Krasovskii functionals can be translated to such easier to handle norm criteria, which goes along with further conservativeness [23, 24].

(vii) An overestimation of the radius of attraction (interior of the red dashed circle with radius $\hat{R}_{X}$ in Figure 1) means an upper bound on possible norm criteria described above. Such an overestimation is needed in order to prove that any larger norm requirement certainly cannot be fulfilled, i.e. (R3). Although no real relations between domains of attraction can be deduced, the knowledge that one system surely possesses a small domain of attraction might also give hints for comparison purposes (R4). Furthermore, there is usually no means of getting an impression how close a lower bound might be to the real radius of attraction, unless an upper bound is provided additionally. As soon as an initial state that does not belong to the domain of attraction is found, its norm is known to be such an upper bound on the radius of attraction. This is what simulative approaches are predestined for and what the present paper addresses.

The present paper is organized as follows. The introduction closes with used notations and needed preliminaries in the context of autonomous retarded functional differential equations. Section 2 addresses domain and radius of attraction as well as generalized concepts thereof. Section 3 shows why the state space selection should be well considered and proposes how the frequent situation of differently delayed state variables can be appropriately treated. Section 4 provides the reader with instructions on how a simulation-based approach should be tackled. In section 5, a bifurcation analysis-based approach is introduced, which results in a possible refinement of previous estimations. Finally, the overall methodology is demonstrated with the delayed swing equation as example system.

1.1 Notation and Preliminaries

Consider the mathematical nomenclature:

$C([a,b],\mathbb{R}^{n})$ space of continuous functions with

$\|\phi\|_{C}\stackrel{{\scriptstyle\textrm{def}}}{{=}}\max\limits_{\theta\in[a,b]}(\|\phi(\theta)\|_{2})$

$PC([a,b],\mathbb{R}^{n})$ space of piecewise continuous functions,

$\|\cdot\|_{PC}=\|\cdot\|_{C}$

$L^{p}([a,b],\mathbb{R}^{n})$ $L^{p}$-space with

$\|x\|_{L^{p}}=(\int_{a}^{b}\|x(\theta)\|^{p}\textrm{d}\theta)^{1/p}$

$\mathbb{R}$, $\mathbb{R}_{>0}$ real (positive) numbers

$\mathcal{B}_{X}(c;r)$ open ball in a normed space $X$

with center $c\in X$ and radius $r\in\mathbb{R}_{>0}$

$\mathcal{B}_{X}(r)$ $=\mathcal{B}_{X}(0_{X};r)$

$0_{[a,b]}$; $0_{n}$ zero function $[a,b]\to\mathbb{R}^{n}$; zero in $\mathbb{R}^{n}$

$\chi_{M}(t)$ indicator function w.r.t. a set $M$

\chi_{M}(t)=\left\{\begin{array}[]{ll}1,&\textrm{if }t\in M\\ 0,&\textrm{otherwise }\end{array}\right.

dom domain of a function

$\partial M$, $\overline{M}$ boundary, closure of a set $M$

$x(t;\phi)$ solution at time $t$ for initial function $\phi$

$x_{t}=x_{t}(\cdot;\phi)$ state, i.e. delay-width segment of $x$ at $t$

$x_{t}\colon[-\tau,0]\to\mathbb{R}^{n}$, $\theta\mapsto x_{t}(\theta):=x(t+\theta)$

In the following, we revisit some required preliminaries concerning autonomous retarded functional differential equations (RFDE). Subsequently, $X$ denotes a normed space on the delay interval $[-\tau,0]$, $\tau\in\mathbb{R}_{>0}$, e.g., the Banach space of continuous functions $X=(C([-\tau,0],\mathbb{R}^{n}),\|\cdot\|_{C})$ endowed with the uniform norm $\|\phi\|_{C}\stackrel{{\scriptstyle\textrm{def}}}{{=}}\max\limits_{\theta\in[-\tau,0]}(\|\phi(\theta)\|_{2})$. The notation of a state as $x_{t}\in X$,

[TABLE]

is common. It represents a restriction of the solution function $x(\tilde{t})$ to the time horizon $\tilde{t}\in[t-\tau,t]$ (Figure 2).

Hence, a forward trajectory of an initial function $x_{0}=\phi$

[TABLE]

describes the set of all states that arise over the whole forward time in which the solution exists $t\in[0,t_{max})$. Thereby, $x_{t}(\cdot;\phi)$ refers to solutions of the Cauchy problem

[TABLE]

with an initial state $\phi\in X$ and $x(t)\in\mathbb{R}^{n}$, $\dot{()}=\frac{\textrm{d}(\cdot)}{\textrm{d}t}$, $F\colon X\to\mathbb{R}^{n}$ continuously differentiable, $F(0_{[-\tau,0]})=0_{n}$. For theorems about well-posedness in certain Banach spaces, see, e.g., Hale and Verduyn Lunel [6] or Kharitonov [38]. We are mainly interested in autonomous differential difference equations [39] with a single discrete delay $\tau>0$, as a special type of RFDE (3), where

[TABLE]

$f\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$. An equilibrium state of (3) is defined as $\phi_{e}(\theta)\equiv x_{e}:x_{t}(\cdot;\phi_{e})=\phi_{e}(\cdot)$, $\forall t\geq 0$ with an equilibrium solution $x(t;\phi_{e})\equiv x_{e}$. The value $x_{e}$ is simply named equilibrium. W.l.o.g. we assume $x_{e}=0_{n}$ for the equilibrium of interest. The focus of the present paper lies on the domain of attraction of this equilibrium. Hence, the trivial equilibrium is required to be attractive. Although it is well known that attractivity and stability according to Lyapunov are actually independent properties in nonlinear systems [40], we are only concerned with their combination, i.e. asymptotic stability. The latter has the advantage that it can easily be proven by the Principle of Linearized Stability (Diekmann et al. [22], Chapter VII, Theorem 6.8) if the equilibrium is hyperbolic.

Definition 1.1 (Local Asymptotic Stability (LAS) in $X$)

The zero equilibrium of the autonomous system (3) is locally asymptotically stable (LAS) in $X$ if

it is stable in $X$, i.e.

[TABLE]

and locally attractive in $X$, i.e.

[TABLE]

Definition 1.1 implicitly assumes forward completeness of solutions in an environment around the zero solution, i.e. the solutions under consideration exist for all $t\in[0,\infty)$.

2 Domain and Radius of Attraction

2.1 Domain of Attraction

The present paper deals with the domain of attraction $\mathcal{D}_{X}\subseteq X$ (also referred to as region of attraction, basin of attraction, or region of asymptotic stability in literature) of an asymptotically stable equilibrium state in $\phi_{e}=0_{[-\tau,0]}$. Thus, the domain of attraction collects all initial functions $\phi\in X$ that lead to a zero-convergent solution.

Definition 2.1 (Domain of Attraction)

The domain of attraction of an asymptotically stable zero equilibrium is defined as

[TABLE]

2.2 Radius of Attraction

Based on the notion of a norm ball around $0_{[-\tau,0]}=\phi_{0}\in X$,

[TABLE]

a strong simplification can be derived:

Definition 2.2 (Radius of Attraction and Ball of Attraction)

The radius of attraction in $X$ of an asymptotically stable equilibrium is defined as

[TABLE]

and the ball of attraction $\mathcal{B}_{X}(R_{X})$ denotes the largest norm ball inside the domain of attraction $\mathcal{D}_{X}$.

The radius of attraction is schematically represented by the black circle in Figure 1. Thereby, the definition of attractivity in (6) ensures $\mathcal{B}(R_{X})\neq\emptyset$ by $R_{X}\geq\delta_{a}>0$. Of course, the exact value of $\mathcal{B}(R_{X})$ remains usually unknown.

Remark 2.1 (Notion of the Radius of Attraction)

The term radius of attraction is analogously used in the context of forward / pullback attractors in nonautonomous delay-free systems [41]. The definition is in accordance with the well-known universal definition of a stability radius as largest ball, whose elements satisfy a certain stability condition. However, it should be carefully distinguished between the question of robustness, i.e. parameter perturbations, which is often addressed by this concept [42, 43], and the domain of attraction, i.e. initial value perturbations.

While the radius of attraction (9) is defined as largest ball in the interior of the domain of attraction, it can also be considered as the minimum distance of the complement $\mathcal{D}_{X}^{c}$ to the origin. This leads to the equivalent definition $R_{X}=\inf_{\phi\in X}\{\|\phi\|_{X}:x(t;\phi)\not\to 0\}$, which is helpful for numerical estimations.

Definition 2.3 (Upper Bound on the Radius of Attraction Based on a Set of Initial Functions)

Assume $\Phi\subset X$ is a set of tested initial functions. The resulting upper bound on the radius of attraction ${\hat{R}_{X}(\Phi)}\geq R_{X}$ is given by

[TABLE]

If no divergent solution can be found, the above definition yields

[TABLE]

The following simple scalar example demonstrates how sensitive such estimations are to the choice of initial functions.

Example 2.1 (Scalar Nonlinear RFDE)

Consider

[TABLE]

with $x(t)\in\mathbb{R}$ and $\phi\in PC([-\tau,0],\mathbb{R})$. Figure 3 shows some solutions, where initial functions in the domain of attraction $\phi\in\mathcal{D}_{PC}$ are printed in black.

Four families of initial functions in $PC([-\tau,0],\mathbb{R})$ are considered. We are interested in the radius of attraction $R_{PC}$, i.e. the maximum uniform norm $\|\phi\|_{C}$ below which all solutions still converge to the zero equilibrium. Clearly, the tested set of initial functions yields $R_{PC}\leq{\hat{R}_{PC}(\Phi)}=\min_{i}\hat{R}_{PC}^{i}=1.1$ if $\tau=1$ and ${\hat{R}_{PC}(\Phi)}=0.4$ if $\tau=5$. In the case of a delay value $\tau=1$ the family of step functions (Fig. 3c) shows the smallest attracted set and the family of decreasing linear functions with $\phi(0)=0$ (Fig. 3d) is converging for all tested functions. However, in the case of $\tau=5$ the step functions (Fig. 3g) even show the largest attracted set, whereas the decreasing linear functions (Fig. 3h) already diverge for smaller $\|\phi\|_{C}$ values. Hence, the class of functions that leads to the best estimation is not even constant for one type of RFDE under parameter variations. As a consequence, there is no obvious relation between the class of initial functions and the conservativeness of estimations on $R_{X}$.

2.3 Generalizations of Domain and Radius of Attraction

The radius of attraction $R_{X}$ is dictated by the smallest initial function that leads to a non-zero-convergent solution - no matter how relevant this initial function might be. As schematically shown by the intruding blue line in Figure 1, the most critical element in the complement $\mathscr{D}_{X}^{c}$ is possibly only a small exception within the closure of the domain of attraction.

For ODEs, there are generalizations of the domain of attraction. Such generalizations are either based on topological concepts, like quasi-stability domains [44, 45, 29], or based on measure theoretic concepts, like almost global stability [46] (see Appendix). These generalizations form supersets of the exact domain of attraction. They ignore that some small sets in the interior actually do not lead to the desired zero-convergence, and, hence, do not belong to the exact domain of attraction. As a consequence, inner estimations of the generalized domain of attraction might be considerably larger than inner estimations of the exact domain of attraction - for the price of marginal exceptions that might be overseen.

In this light, we propose definitions of generalized radii of attraction $R_{X}^{gen}$ in time-delay systems in the Appendix. It can be expected that upper bounds $\hat{R}_{X}$ that are based on arbitrarily selected initial functions, give only an estimation of the such introduced $R_{X}^{gen}$ (Figure 1). However, $R_{X}^{gen}$ might even be the actual number of interest: if a statement about generic convergence to the attractor suffices, $R_{X}$ would be unnecessarily restrictive.

Remark 2.2 (No Lebesgue Measure in Infinite Dimensional Spaces)

The notion of a property to be valid for almost every (a.e.) point in the state space or a subdomain is convenient in the context of ordinary differential equations. If not otherwise stated, exclusions of Lebesgue measure zero are meant. However, there is no analog of the Lebesgue measure on infinite dimensional spaces (Sullivan [47], Theorem 2.38). As a consequence, there is no universally accepted notion of “almost every function” [48].

3 On the Selection of the State Space

There are various normed spaces - predominantly Banach spaces - that come into question as state space for RFDEs. In literature, the space of continuous functions $C([-\tau,0],$ $\mathbb{R}^{n})$ endowed with the uniform norm is mostly used [6, 22, 39, 10, 49]. Also frequent [50, 51, 52, 53] are spaces $M^{p}=\mathbb{R}^{n}\times L^{p}([-\tau,0],\mathbb{R}^{n})$ with $1\leq p\leq\infty$, especially $M^{2}$. A $L^{p}$ space alone would not be sufficient, since a RFDE requires the point value $x(0)=\phi(0)\in\mathbb{R}^{n}$ at $t=0$ to start with, which is a set of Lebesgue measure zero. In $M^{p}$ both, the initial function $\phi\in L^{p}$ and the initial value $\phi(0)\in\mathbb{R}^{n}$, are given separately. Space selections in literature seem at most to be aligned to the mathematical theory. However, due to non-equivalence of norms in infinite dimensional spaces, the selection of the normed space has far reaching consequences and should be well considered.

3.1 Consequences of the Selection

According to (6), there is a non-zero radius of attraction $R_{X}\neq 0$ in any normed space, in which the equilibrium under consideration can be proven to be attractive. However, the interpretation as well as the size of the radius of attraction depend heavily on the chosen norm. For instance, the radius of attraction might be a measure for maximum allowed pointwise perturbations ($\|\cdot\|_{C}$) or address derivations in energy ($\|\cdot\|_{L^{2}}$ as part of $\|\cdot\|_{M^{2}}$). Thus, the technical or physical application is of interest.

The space also determines which initial functions come into question. Consider delayed input in a control law. It requires a certain initialization until first data is received. As soon as the signal becomes available, this initialization is replaced by the actual value, such that a jump discontinuity occurs unless the values are equal by chance. If the start-up process or a communication interruption is only one of many possible scenarios, not only the finite dimensional case of initial jumps, but an infinite dimensional function space, which allows such discontinuous functions, has to be taken into account. This example shows that the Banach space of continuous functions $C([-\tau,0],\mathbb{R}^{n})$ is frequently not sufficient. However, the Banach space of piecewise continuous functions $PC([-\tau,0],\mathbb{R}^{n})$ equipped with the uniform norm is a valid alternative for RFDEs [54, 38] and does not alter the meaning of the radius of attraction. Obviously, other spaces like $M^{2}$ also allow discontinuous functions, but the interpretation of the radius of attraction is completely different.

Remark 3.1 (Invariance under Time Scale Transformations)

Time scale transformations are convenient to get dimensionless parameter values in appropriate scales and to lower the number of parameters. However, by affecting the time delay value, time scale transformations do also change the domain of the initial function. What are the consequences for the radius of attraction? For instance, a delay normalizing transformation is achieved by $t=\tau\hat{t}$, $x(\tau\hat{t})=\hat{x}(\hat{t})$ and $\theta=\tau\hat{\theta}$, $\phi(\tau\hat{\theta})=\hat{\phi}(\hat{\theta})$, respectively. Hence, $\dot{x}=f(x(t),x(t-\tau)),x_{0}=\phi$ becomes

[TABLE]

The radius of attraction with respect to the transformed system refers to $\hat{\phi}(\hat{\theta})$. Obviously, by

[TABLE]

the uniform norm of the initial function and thus the radius of attraction in $C$ or $PC$ is invariant under such transformations. However, for instance the norm in $M^{2}=\mathbb{R}^{n}\times L^{2}$ is not since

[TABLE]

As a consequence, the radius of attraction in $M^{2}$ of the transformed system has no meaning for the original RFDE.

3.2 Quotient Space to Incorporate Differently Delayed States

In physical or technical problems it is common that some state components occur in a delayed form, but others contribute only by their instantaneous values (see e.g. [18, Scholl.2019b, 55] or Example 6.1 below). Hence, the history of some state variables $x^{\textrm{\scriptsize II}}(t)\in\mathbb{R}^{q}$, $q<n$ in $x(t)=[x^{\textrm{\scriptsize I}}(t)^{\top},x^{\textrm{\scriptsize II}}(t)^{\top}]^{\top}\in\mathbb{R}^{n}$ does not influence the system dynamics at all, i.e.

[TABLE]

Usually, there is no distinction in the domains of the initial function components, such that $\textrm{dom}(\phi^{\textrm{\scriptsize II}})=\textrm{dom}(\phi^{\textrm{\scriptsize I}})=[-\tau,0]$. Thereby, it is ignored that some of the required initial data $[\phi^{\textrm{\scriptsize I}{\top}},\phi^{\textrm{\scriptsize II}{\top}}]^{\top}=\phi\in X([-\tau,0],\mathbb{R}^{n})$ remains unused (cmp. [56]). Thus, the Banach space norm $\|\phi\|_{X}$ depends on irrelevant values of $\phi^{\textrm{\scriptsize II}}$ for $\theta<0$. Consider $\|\phi^{\textrm{\scriptsize II}}\|_{X}>\|\phi^{\textrm{\scriptsize II}}(0)\|$. Obviously, the radius of attraction as a norm criterion on $\|\phi\|_{X}$ would be unnecessarily restrictive.

This situation motivates the use of

[TABLE]

with a certain choice of the outer norm (which is not that important by equivalence of norms in finite dimensions). However, $\|\cdot\|_{Q}$ is only a seminorm in $X$. Therefore, we propose to introduce the quotient space

[TABLE]

endowed with the norm $\|\cdot\|_{Q}$. It regards functions to be in the same equivalence class if they differ in values $\phi^{\textrm{\scriptsize II}}(\theta)$ for $\theta<0$. Hence, these values are fully ignored in the sufficient criterion for zero-convergence $\|\phi\|_{Q}<R_{Q}$.

4 Estimations by Time-Forward Simulations

The state space $X$ is infinite dimensional. However, a simulative approach means to select certain initial functions. Thereby, Example 2.1 demonstrates that an obvious relation between the class of initial functions and the conservativeness of the overestimation of $\hat{R}_{X}$ does not exist. Hence, the question arises, which initial functions should be taken into account.

The physical or technical context might already motivate certain initial functions.

(i) Predefined parameterized initial functions might be of particular importance, e.g., jump functions for delayed controller input (cf. Section 3.1).

*(ii) *In the case of a switched system, the initial function should correspond to the dynamics of the previous system definition.

*(iii) *There is also the proposal to use solution segments built from the solutions of the system without delayed terms [57].

If application-adapted functions occur exclusively, only the corresponding subset of $X$ is of interest. Consequently, the intersection of the domain of attraction with this subset suffices. Otherwise, a large span of different families of initial functions should be tested. In the following section various possibilities of initial functions are classified. The norm of the smallest found initial function that does not belong to the domain of attraction is of interest. According to Definition 2.3, it gives the desired overestimation $\hat{R}_{X}\geq R_{X}$ of the radius of attraction.

Results in the present paper are generated by the Matlab solver dde23, which is based on a Runge-Kutta triple [58]. Stiff problems might require alternative approaches [59]. For further details on numerical methods we refer to [60]. No matter which solver is chosen, a tuning of the discretization step size or error bounds is inevitable. Furthermore, the criterion according to which a numerical solution is assumed to approximate a zero-convergent solution has to be selected. Consider a numerical solution over the finite time span $t\in[-\tau,t_{\textrm{num}}]$ for a given initial function $\phi$. A convergence decision must be based on the approximation of the last state $x_{t_{\textrm{num}}}(\cdot;\phi)$, i.e. on the numerical data for $t\in[t_{\textrm{num}}-\tau,t_{\textrm{num}}]$. Let $r(t_{\textrm{num}};\phi):=\|\hat{x}_{t_{\textrm{num}}}(\cdot;\phi)\|_{X}$ be the respective norm of the interpolated numerical values. The numerical domain of attraction is

[TABLE]

cmp. (6), with $\delta_{\textrm{num}}>0$ suitably chosen. If a lower bound on the radius of attraction $\check{R}_{X}\leq R_{X}$ is available $\delta_{\textrm{num}}$ should be based on this result. Furthermore, $\delta_{\textrm{num}}$ can be used as a termination criterion in the simulation to lower the numerical effort. In the following, with a slight abuse of notation, we identify the numerical domain of attraction $\hat{\mathcal{D}}_{X}$ with $\mathcal{D}_{X}$.

Remark 4.1 (Numerical Errors)

All results on the domain of attraction that are based on time simulations [32, 19, 21, 61, 62, 57, 63, 64] rely on the distinction between convergence and divergence from numerical results. However, this asymptotic behavior $t\to\infty$ is usually not covered by any numerical bounds, and thus the classification cannot be considered as proven. Concerning the radius of attraction, trajectories wrongly classified as convergent will not contribute to the upper bound estimation, but trajectories wrongly classified as non-zero-convergent might give a too small upper bound. While from a mathematical point of view this is not satisfactory, it should be noted that if the numerical errors (e.g. numerical damping) are able to generate non-convergence for a given initial function, then comparable errors in system parameters (e.g. physical damping), equation structure, or external perturbations will do so a fortiori.

4.1 Primary Initial Functions

We use the term primary initial functions to describe simple functions $\phi^{k}(\,\cdot\,;p)$, $k\in\{1,\ldots,K\}$ like constant, jump, polynomial, or trigonometric functions, which are parameterized in $p\in\mathbb{R}^{m}$. Denote the set of all elements in a family $\{\phi^{k}\}_{p}$ of initial functions of type $k$ shortly as

[TABLE]

and the set of all considered functions as

[TABLE]

There are few considerations of domains of attraction in literature. However, the existing ones are mostly based on scalar linear or sinusoidal initial functions dependent on one [63] or two [57, 64] parameters $(p_{1},p_{2})\in\mathbb{R}^{2}$. In the latter case a graphical representation of the $(p_{1},p_{2})$ plane is common. By marking those parameter combinations that result in an attractor-convergent solution, the finite dimensional intersection $\Phi^{k}\cap\mathcal{D}_{X}$ of the domain of attraction with the chosen family of primary initial functions becomes visible (black pixels in Figure 4).

These plots are by far more computationally expensive than phase plots for ordinary differential equations. By the semigroup property, each point $x(t)$ of an ODE solution can be seen as new initial point and does not need to be checked again for zero-convergence. However, the states $x_{t}$ in a delay system are whole solution segments that are probably not in the chosen subset of initial functions $x_{t}\not\in\Phi^{k}$. In view of avoiding unnecessary computational burden, the examination of the whole parameter space is not undertaken: to derive estimations $\hat{R}_{X}(\Phi)$ for the radius of attraction, it is sufficient to increase the parameters until zero-convergence fails. This can be done in a circular or star-like scheme, in which the parameter modulus is increased (Figure 7 below). Of course, symmetries should be taken into account to reduce the computational effort further. In addition, the problem is well suited for parallel computing.

Remark 4.2 (Fractal Domain of Attraction Boundaries)

The boundary of a domain of attraction is possibly fractal [62, 57, 63, 64]. Fractality results in a high sensitivity of simulative results on parameter variations in the initial function. However, it does not influence the validity of gained upper bounds on the radius of attraction.

It should be noted that physics might give restrictions on how the state components are related to each other.

Remark 4.3 (Higher Order Derivatives)

Assume that the state space representation $\dot{x}(t)=F(x_{t})$ with $x(t)\in\mathbb{R}^{n}$ stems from a system description in a scalar variable $y(t)\in\mathbb{R}$ with higher order derivatives

[TABLE]

Then, the transformation to $x$, e.g., $x=[y,\dot{y},\ldots,y^{(n-1)}]^{\top}$, goes along with restrictions to the space of initial functions. Hence, the $n$-dimensional vector function $\phi(\,\cdot\,;p)=x_{0}(\,\cdot\,;p)$ cannot freely be chosen. Instead, only a scalar initial function $y_{0}(\,\cdot\,;p)$ for $y(t)$ must be specified. For instance, with the transformation from above, this results in

[TABLE]

4.2 Secondary Initial Functions

We use the term secondary initial functions for functions that arise a posteriori from previous simulation results. As described in Section 4.1, the semigroup property $x_{s+t}(\cdot;\phi)=x_{t}(\cdot;x_{s}(\cdot;\phi))$ does not help much with respect to considerations of $\Phi^{k}\cap\mathcal{D}_{X}$, which are addressed in graphical representations like Figure 4. Nevertheless, it does help for considerations of the whole domain of attraction $\mathcal{D}_{X}$, which indeed is addressed by the radius of attraction. Actually, all segments of simulated solutions $x(t)$ over an interval $t\in[s-\tau,s]$, i.e. states $x_{s}$, can be interpreted as further initial functions $\phi=x_{s}$. All these intermediate results have automatically been tested for zero convergence as well since

[TABLE]

As a consequence, it is worthwhile to further examine divergent solutions. Already Figure 3g indicates that a diverging trajectory might temporarily get closer to the equilibrium than the initial function itself. Figure 5 shows how an oscillating solution segment with a smaller norm can arise from a constant primary initial function.

We denote the minimum norm value of diverging trajectories (red dashed line in Figure 5) as

[TABLE]

where $\Phi$ is the set of primary initial functions (23). Obviously, without a noteworthy need of additional computational effort, $\hat{R}_{X}^{\mathcal{T}}(\Phi)$ can only be equal or smaller than the estimation gained by primary initial functions ${\hat{R}_{X}(\Phi)}$. Indeed, in Example 6.1 below, a tighter upper bound on the radius of attraction than the minimum norm of all considered primary initial functions will be achieved.

4.3 Extended Primary Initial Functions

A basis extension can give further improvements of the estimation for the radius of attraction. Thereby, we mean to construct

[TABLE]

with basis functions $b_{i}\in X([-\tau,0],\mathbb{R})$ up to order $d>0$. The difference to primary initial functions (Section 4.1) lies in the higher number of parameters $[p_{1},\ldots,p_{d}]=:p\in\mathbb{R}^{n\times d}$ and thus in the computational effort. This ansatz is also taken by Leng et al. [32] who are not interested in an estimation for $R_{X}$ but in basin stability. The latter is up to normalization approximated by the percentage of those randomly chosen parameter values $p$ in a predefined set $p\in\mathcal{B}_{d}\subset\mathbb{R}^{n\times d}$ that lead to a zero-convergent solution (Section 1 (ii)). For instance, trigonometric, Legendre or Bernstein polynomials can be used as basis functions. Of course, the results allow again to derive secondary initial functions.

Remark 4.4 (Approximation properties)

Similarily to most inner approximations of domains of attraction, no statement about the tightness of the derived upper bound can be given. Neither a statement about convergence for extended primary initial functions is possible. Indeed, arbitrarily close approximations of initial functions in $X=C([-\tau,0],\mathbb{R}^{n})$ are achievable if polynomials with sufficiently high degree are considered (Stone-Weierstrass theorem). Moreover, if continuous dependence on the initial conditions is ensured (Hale and Verduyn Lunel [6], Ch.2, Theorem 2.2), then a sufficiently close approximation of an element in $\mathcal{D}_{X}$ will also reach the ball $\mathcal{B}_{X}(0;\delta_{a})$ of (6) in finite time. Consequently, zero-convergence follows and thus the approximation also belongs to $\mathcal{D}_{X}$. However, this argument only holds for elements in $\mathcal{D}_{X}$, while the radius of attraction is determined by those initial functions that lie on the boundary $\partial\mathcal{D}_{X}$. Asymptotic behavior of an element in this boundary set might not be reflected by any approximation. Hence, in view of the arbitrary difference between the radius of attraction $R_{X}$ and the potentially larger radius of attraction $R_{X}^{qs}$ w.r.t. the quasi-stability domain $\textrm{int}(\overline{\mathcal{D}_{X}})$ (Definition 8.1), even a large set of tested initial functions might not result in tight estimations of the radius of attraction. In the end, the combination of an analytically derived lower bound $\check{R}_{X}$ and a numerically derived upper bound $\hat{R}_{X}$, such that  $R_{X}\in[\check{R}_{X},\hat{R}_{X}]$, is most meaningful. The gap $\hat{R}_{X}-\check{R}_{X}$ provides a tightness bound for both estimations.

5 Estimations Based on Bifurcation Analysis

A positively invariant or even invariant set with some distance to the attractor cannot be contained in the domain of attraction and thus yields an upper bound on the radius of attraction.

Proposition 5.1 (Upper Bound on the Radius of Attraction Based on a Positively Invariant Set)

Assume $S$ is a positively invariant set and $\overline{S}\cap\{0_{[-\tau,0]}\}=\emptyset$. Then $S\cap\mathcal{D}_{X}=\emptyset$ and

[TABLE]

Proof.

If $\phi\in S$ is an element of a positively invariant set $S$ with some distance to the attractor in $\phi_{e}=0_{[-\tau,0]}$, the defining property of the domain of attraction $\lim_{t\to\infty}\|x_{t}(\theta;\phi)\|_{X}=0$ cannot be fulfilled. ∎

In the approaches described above, non-zero-convergent forward trajectories from time-domain simulations represent this set $S$. Furthermore, these forward trajectories approach $\omega$-limit sets, which are themselves invariant sets (complete continuity of $F$ in (3) and boundedeness of trajectories presumed; Smith [65], Corollary 5.6). However, there might also be equilibria, limit cycles, invariant tori, or chaotic sets, which are hidden [66] in a certain sense or unstable and thus will not be approached in time-domain simulations. However, it is challenging to locate these limit sets and even to prove their existence. In the following, we try to get a priori knowledge of some of these sets.

Remark 5.1 (Number of Limit Cycles)

It should be noted that the procedure described below will not necessarily capture all limit cycles. Even in the supposedly simple case of planar ODEs with polynomial vector fields of degree $d>1$ no finite upper bound for the number of limit cycles has yet been found. This is known as Hilbert’s 16th problem [67].

We propose to start with the delay-free, thus finite dimensional and manageable system that results in setting the delay value or the coefficient of the delayed terms to zero. All additional equilibria, limit cycles, invariant tori, or chaotic sets have to emerge with an increase of these parameters. Hence, we consider the delay parameter or the respective coefficient as bifurcation parameter in a bifurcation analysis (for bifurcation theory of RFDEs see [6, 22, 68, 69, 70]). It is well known, that the loss and regain of stability by parameter variation is generically accompanied by a Hopf bifurcation. A subcritical Hopf bifurcation means that an unstable limit cycle occurs for parameter values at which the equilibrium is stable (Figure 6). In a finite dimensional system in the plane, the situation would be clear: the unstable limit cycle, which surrounds the asymptotically stable equilibrium, represents the boundary of the domain of attraction. The infinite dimensionality of time-delay systems makes conclusions harder. However, by Proposition 5.1 the knowledge about such limit sets can be used for estimations of the radius of attraction.

As in the ODE case, it is possible to formulate a boundary value problem (BVP) that is solved by the periodic solutions of a limit cycle. Periodicity of the whole solution $x(t)=x(T+t),t\in[-\tau,\infty)$ is by semigroup property already ensured, if $x(\theta)=x(T+\theta),\theta\in[-\tau,0]$, i.e. by equality of the two states $x_{0}=x_{T}$. Hence, instead of an initial condition $x_{0}=\phi$ the boundary condition $x_{0}=x_{T}$ is required, such that a two-point BVP over one period $t\in(0,T]$ has to be solved. The dependence of the boundary position on the beforehand unknown parameter $T$ can be removed by a time scale transformation. Consider $\dot{x}=f(x(t),x(t-\tau))$. With $t={\hat{t}}T$ and $x(\hat{t}T)=:z(\hat{t})$ the two-point BVP over the normalized period $\hat{t}\in(0,1]$ becomes

[TABLE]

Still, it is ambiguous which delay-spanning segment of the periodic solution is addressed. This translational symmetry has to be resolved by an additional phase shift condition leading to a unique solution. However, the above boundary value problem can only be solved numerically if an appropriate initial estimation of the limit cycle is given. To this end, we consider a numerical bifurcation analysis with time delay as bifurcation parameter. Usually, the latter is used if parameter variations in a given range are of interest. However, as a byproduct, effects of parameter variation on domains of attraction might become visible [71, 72, 73]. We are mainly interested in a fixed parametrization, which at a first glance is not related to bifurcation analysis. However, since in a bifurcation analysis limit cycles are tracked from their emergence, (30) is well initializable in an iterative manner based on successive parameter steps until the parametrization of interest is reached.

We make use of an implementation in the Matlab toolbox DDE-BIFTOOL [74], which is based on a collocation method with piecewise polynomials [75]. The latter provides a numerical approximation of the periodic solution.

Already for simple lower order systems, the infinite dimensionality of time-delay systems frequently results in complex phenomena and a rich variety of bifurcations [76, 77, 72, 73, Scholl.2019b]. Of course, bifurcations of the above described limit cycles might occur with further delay parameter variation. For instance, limit cycles vanish if they collide in a fold bifurcation, a new limit cycle might branch off in a period doubling bifurcation, or an invariant torus, which surrounds a limit cycle, might stem from a Neimark-Sacker bifurcation [78]. A bifurcation analysis aims at detecting such phenomena. Approximations of Floquet multipliers allow to conclude stability. After the numerical analysis, if a certain fixed delay value is given, some occurring invariant sets can be concluded. The periodic solution leads to the limit cycle trajectory (set of all delay-width solution segments) that we use as invariant set in Proposition 5.1. Even more, the bifurcation analysis results in valuable insights how the system behavior changes with increasing delay. Thereby, also the variation of the derived upper bound for the radius of attraction with increasing delay values can be observed.

Remark 5.2 (Stable Manifolds)

Generically, unstable equilibria or unstable limit cycles in RFDEs are saddle-type sets with a finite number of unstable eigenvalues or Floquet multipliers, respectively, but an infinite number of stable ones. Hence, there is also a stable manifold corresponding to the stable eigenfunctions of the linearized system [22]. Indeed, the stable manifold of a saddle-type set or the domain of attraction of a competing attractor is also a possible invariant set in Proposition 5.1. However, these sets remain usually unknown. Nevertheless, knowledge from ODEs suggests that solutions nearby the stable manifold approach the saddle-type set temporarily in time-forward simulations. As a result, the saddle-type set can to some extent already be reflected in secondary initial functions. This conjecture is supported by comparing the secondary initial function in Figure 5 and the limit cycle solution in Figure 6a.

Remark 5.3 (Relation to Bi- or Multistability)

A coexisting second attractor is an example of an invariant set in Proposition 5.1 and thus implies that the first attractor has a finite radius of attraction $R_{X}<\infty$ and cannot be globally asymptotically stable. The coexistence of two or more attractors is commonly denoted as bistability or multistability, respectively [79, 63, 80, 81, 82]. Thereby, it is a matter of definition whether - based on a space compactification - infinity can also be considered as an attractor in a non-dissipative system. All unbounded trajectories form its domain of attraction (without a meaningful radius of attraction). Conversely, non-global asymptotic stability of an equilibrium (i.e. $R_{X}<\infty$) does not imply the existence of a second attractor (including an attractor at infinity). To this end, consider the scalar ODE $\{\dot{x}=-x$, if $x\leq 1$ and $\dot{x}=0$, if $x>1\}$, which has only one attractor at $x_{e}=0$ with domain of attraction $(-\infty,1)$ and radius of attraction $R_{\mathbb{R}}=1$ (the other equilibria are non-isolated and thus non-attractive).

6 Demonstrative Example

The following example demonstrates the proposed approach for estimating the radius of attraction. Second order ordinary differential equations with periodic nonlinearity are relevant for various non-related applications like mechanical pendulum mechanisms, Josephson junction circuits, phase-locked loops, or synchronous generators [34]. At least in the latter case, the impact of delayed damping on the domain of attraction is of high technical interest [21]. Schaefer et al. [21] give a simulative analysis of basin stability for such a delayed swing equation.

Example 6.1 (Swing Equation with Time Delay)

Consider the swing equation or constantly driven pendulum equation with additional delayed damping described by $\ddot{y}(t)+a\dot{y}(t)+\tilde{a}\dot{y}(t-\tau)+\sin(y(t))=w$. A state space representation in $x=[y-y_{e},\dot{y}]$, such that the equilibrium $y_{e}:=\arcsin(w)$ is shifted to the origin $x_{e}=[0,0]^{\top}$, reads

[TABLE]

with $x(t)\in\mathbb{R}^{2}$, $a,\tilde{a},\tau\in\mathbb{R}_{>0}$, $w\in(0,1)$, and $\phi\in X$ in a Banach space $X$ to be selected. Unless otherwise stated, results refer to the parameter set

[TABLE]

The aim is to find an upper bound $\hat{R}_{X}>R_{X}$ for the radius of attraction $R_{X}$ in the chosen norm.

(i) Selection of the State Space $X$

The history of the first component $x_{1}$ does not contribute to the system dynamics at all. According to (20), this motivates the use of a quotient space $X=Q$. We define

[TABLE]

Remark 6.1 (Definition on a Manifold)

There is a rotational symmetry in $x_{1}$, such that - instead of $x(t)\in\mathbb{R}^{2}$ - the instantaneous values can also be defined on a cylinder, i.e. $[x_{1}(t),x_{2}(t)]^{\top}\in S^{1}\times\mathbb{R}$ and consequently

[TABLE]

Thereby, the number of stable equilibria reduces from infinite to one in Example 6.1 since all belong to the same equivalence class. Stability theory of RFDEs can be transferred to such systems on manifolds [6, 83]. (The respective far reaching consequences are not in the scope of this paper.)

(ii) Selection of Primary Initial Functions

Actually, a physical meaning of $x_{1}$ and $x_{2}$ as angle and angular velocity in (31) implies $\frac{\textrm{d}}{\textrm{d}t}\phi_{1}=\phi_{2}$ (Remark 4.3). However, because of the irrelevance of the $x_{1}$-history in the above introduced quotient space (33), this relation has no consequence for the selection of initial functions. As representative of the first component $\phi_{1}\colon[-\tau,0]\to\mathbb{R}$ we use

[TABLE]

where $p_{1}\in\mathbb{R}$ parametrizes this initial function. For the second component $\phi_{2}\colon[-\tau,0]\to\mathbb{R}$, we select a finite number of primary initial functions parameterized by $p_{2}\in\mathbb{R}$ (Section 4.1):

[TABLE]

In order to simplify the evaluation, the primary initial functions above are chosen such that $\|\phi_{2}^{k}(\,\cdot\,;p_{2})\|_{C}=|p_{2}|$, $k\in\{1,\ldots,4\}$. This parametrization is advantageous since $\|\cdot\|_{Q}$ norms of the initial functions $\phi\in\Phi^{k}=\{[\phi_{1},\phi_{2}^{k}]^{\top}\}_{k}$ equal Euclidean norms in the $(p_{1},p_{2})$ plane, i.e.

[TABLE]

Thus, a norm ball $\mathcal{B}_{Q}(r){=}\{\phi\in Q:\|\phi\|_{Q}<r\}$ corresponds to the inner of a circle in the parameter plane.

Remark 6.2 (Case $\tau\to 0$)

If the delay $\tau\to 0$ vanishes, the domain of $\phi$ reduces to a singleton $[-\tau,0]\to\{0\}$ and the state space $Q$ is replaced by $\mathbb{R}^{2}$, which is fully covered by the two values $[p_{1},p_{2}]^{\top}=[\phi_{1}(0),\phi_{2}(0)]^{\top}=[x_{1}(0),x_{2}(0)]^{\top}$ describing $\Phi$. Then $\|\cdot\|_{Q}$ is replaced by $\|\cdot\|_{2}$, such that $\hat{R}_{2}=R_{2}$ describes exactly the maximum inner circle $\mathcal{B}_{2}(R_{2})$ within the domain of attraction of the ODE in the phase plane.

(iii) Simulative Results for Primary Initial Functions

The domain of attraction for the $k$-th family of initial functions (22) is given by

[TABLE]

Figure 4 shows results in the $(p_{1},p_{2})$ parameter plane. To get an approximation of the relevant parts of the boundary $\partial(\mathcal{D}_{Q}\cap\Phi^{k})$, parameters are increased in a star-like scheme with a finite number of equidistant angles, until convergence to the origin fails. For visualization purposes, Figure 7 contains the (possibly non-convex) hull of the outer parameter values for each family of initial functions.

Obviously, not the inner circle of the visualized hull but the minimum norm of the single vertices results in the upper bound

[TABLE]

All four families of initial functions $\Phi^{k}$, $k\in\{1,\ldots,4\}$ taken in Figure 7 lead to $\min_{k}\hat{R}^{k}_{Q}=\hat{R}^{3}_{Q}=0.65$. The latter corresponds to the most central blue point.

(iv) Secondary Initial Functions

By consideration of secondary initial functions, as described in Section 4.2, an improvement of this estimation is possible without the need for further simulations. Instead of only considering the initial states $x_{0}=\phi\in\Phi^{k}$, the estimation $\hat{R}_{Q}^{k,\mathscr{T}}:=\hat{R}_{Q}^{\mathscr{T}}(\Phi^{k})$ uses the minimum norm of all solution states $x_{t}(\cdot\,;\phi)$, $t\geq 0$ that result from non-zero-convergent primary initial functions $\phi\in\Phi^{k}$. The final estimation $\hat{R}_{Q}^{\mathscr{T}}=\min_{k}\hat{R}_{Q}^{k,\mathscr{T}}$ merges the results of all families of initial functions. The dashed line in Figure 7 indicates that there are indeed diverging solutions that come even closer to the attractor than all the considered primary initial functions. Secondary initial functions of all chosen families $\Phi^{k}$, $k\in\{1,\ldots,4\}$ yield similar minimum norm values $\hat{R}_{Q}^{k,\mathscr{T}}$. The best result is gained by a constant initial function, which leads to a diverging solution with minimum Q-norm $\hat{R}^{\mathscr{T}}_{Q}=\hat{R}^{1,\mathscr{T}}_{Q}=0.503$, as shown in Figure 5.

(v) Bifurcation Analysis

The trivial equilibrium is alternately stable and unstable with increasing delay. A bifurcation analysis of (31) reveals that each regain of stability is accompanied by a subcritical Hopf bifurcation (cmp. [84]). Vertical slices in Figure 6 show the emerging unstable limit cycle in the $(x_{1},x_{2})$ plane for different delay values. For the chosen delay value $\tau=20$ the corresponding periodic solution in Figure 6a is remarkably similar to the state that was gained as secondary initial function in Figure 5. This can be explained by the saddle character of the limit cycle (Remark 5.2). A comparison of the norm values is given in Figure 8. Thereby, the role of the the quotient space norm becomes visible. The radius of attraction $R_{PC}$ in the space $PC([-\tau,0],\mathbb{R}^{2})$ is identical to $R_{Q}$ in the quotient space. However, the uniform norm of limit cycle states would lead to unnecessarily conservative estimations $\hat{R}^{LC}_{PC}$ (Figure 8). The minimum $Q$-norm value of the limit cycle states $x_{t}^{LC}$ is $\hat{R}^{LC}_{Q}:=\min_{t\in[0,T)}\|x_{t}^{LC}\|_{Q}=0.507$.

Hence, for the delay value under consideration, secondary initial functions and the bifurcation analysis based approach yield almost equal upper bounds $\hat{R}^{LC}_{Q}\approx\hat{R}^{\mathscr{T}}_{Q}$ on the radius of attraction $R_{Q}\leq\hat{R}_{Q}=:\hat{R}^{\mathscr{T}}_{Q}$ in Example 6.1. Obviously, this does not always hold true. For small delay values in Figure 8, already the chosen simple primary initial functions lead to better estimations. This means that primary initial functions with smaller norm values than the limit cycle lead to diverging solutions. Hence, these must cross the limit cycle in the $(x_{1},x_{2})$ plane. Time delay enables such crossings since forward uniqueness in RFDEs does not refer to instantaneous values $x(t)\in\mathbb{R}^{2}$ but to solution segments, i.e. to states $x_{t}\in Q$ in the infinite dimensional space $X=Q$. Additionally, Figure 8 illustrates that secondary initial functions can only lead to improvements. However, Figure 8 also demonstrates that the bifurcation analysis based approach is a worthwhile supplement. Compared to the search for diverging primary initial functions in time-domain simulations, the numerical solution of boundary value problems for limit cycles is less computationally expensive. In addition, the bifurcation analysis gives further insights. The growing limit cycles, which are revealed by a bifurcation analysis, result in an enlargement of the upper bound on the radius of attraction with increasing delay.

7 Conclusion

Though as important as in delay-free systems, domains of attraction in time-delay systems still find relatively little consideration in literature. The previous sections contribute the following:

• •

*Estimating Domains of Attraction:

Approaches and Quantification*

Various requirements concerning the domain of attraction are important in practical applications. An overview of possible approaches for estimations of the domain of attraction in time-delay systems is given. The radius of attraction in a certain Banach space represents an easy to handle number. In many cases a lower bound on the radius of attraction is of interest. These estimations are mostly conservative and will be subject of further considerations. Nevertheless, such results leave unavoidably open whether a small radius means a conservative estimation or indeed a small domain of attraction. In order to prove the latter, upper bounds on the radius of attraction become important. These can be obtained by simulative means.

• •

*Generalized Concepts for the Domain of Attraction:

No Lebesgue Measure in Infinite Dimensional Spaces*

It should be kept in mind that simulation usually does only cover generic cases. In the context of ordinary differential equations, genericity is addressed by concepts like quasi-stability domains or almost global stability. The paper proposes to generalize these definitions to the infinite dimensional case, thereby the fact that there is no Lebesgue measure in infinite dimensional spaces must be taken into account.

• •

*Selection of the Banach Space:

Consequences for the Radius of Attraction*

The radius of attraction and its physical interpretation depend heavily on the chosen norm. It is pointed out that the radius of attraction in some spaces is not invariant under time scale transformations.

• •

*Quotient Spaces:

Incorporate that States are Differently Delayed*

We propose to choose a quotient space if some of the state variables do not occur delayed in the system equations. Otherwise unused initial data leads to unnecessary restrictions.

• •

*Secondary Initial Functions:

Improvement of the Examined Set of Initial Functions *

An important question is how to choose the initial functions that are numerically tested for belonging to the domain of attraction. Besides physically or technically motivated functions, primary functions in the sense of simple function families and the construction of initial functions with more degrees of freedom by certain basis functions are convenient. The paper proposes to take into account all available solution segments from simulations as well. These secondary initial functions are able to improve estimations for the radius of attraction without notable computational effort.

• •

*Bifurcation Analysis:

Getting A Priori Knowledge of Other Invariant Sets*

Competing equilibria, limit cylces, invariant tori, chaotic sets, or more general invariant sets with some distance to the attractor can also contribute to upper bounds on the radius of attraction. While the delay-free system might be well understood, time delay is able to induce complex dynamics. In order to detect some further limit sets, we propose a numerical bifurcation analysis. By starting from a delay-free system, a bifurcation analysis allows firstly to detect their emergence and then to track them numerically.

• •

*Demonstrative Example:

Swing Equation with Additional Delayed Damping*

An example demonstrates significant improvements in the estimation of the radius of attraction by secondary initial functions. It also confirms that the proposed bifurcation analysis based approach is a worthwhile additional element for upper bound estimations of the radius of attraction.

Acknowledgement

The authors acknowledge support by the Helmholtz Association under the Joint Initiative "Energy Systems Integration" (ZT-0002).

8 Appendix

8.1 Definitions of Generalized Domains and Radii of Attraction

For some purposes it is sufficient to describe a generalized domain of attraction in which "generically" convergence to the attractor occurs. There might be not attracted sets $M\not\subset\mathcal{D}_{X}$ of first (Baire) category in the interior of $\overline{\mathcal{D}_{X}}$ (see Chiang [45], Figure 2 for an ODE example). A numerical testing for convergence will typically not capture such sets. At the same time they form essential restrictions to analytical estimations that do not allow any exceptions (Figure 1). Nevertheless, meager sets are also of less relevance in practical applications. As a consequence, the discussion in Section 1 must be extended by requirements like

(R1’)

testing whether a certain initial function is within a generalized domain of attraction,

(R2’)

proving generic convergence to the attractor,

(R3’)

disproving generic convergence to the attractor,

(R4’)

comparing systems w.r.t. generalized domains of attraction.

Such generalizations of the domain of attraction are well known from ODEs. The above described topological point of view is only one possibility. It is taken by Zaborszky et al. [44], who introduce the topological concept of a quasi-stability boundary, which is defined as boundary of $\textrm{int}(\overline{\mathcal{D}_{X}})$. The main application is (delay-free) power system stability [45, 29]. In contrast to this topological concept, there is also the measure theoretical one, which is based on convergence for almost every (a.e.) initial value in a domain. Rantzer et al. [46] introduce the approach of Lyapunov densities for almost global stability of equilibria in ODEs and thereby ignore non-attracted sets of Lebesgue measure zero. Measure theoretic considerations can further be transferred to a probability space such that (R1’) means convergence to the attractor with probability one / almost sure convergence (equivalently). It should be noted that the topological definition of a set of first category / a meager set on the one hand and the measure theoretical definition of a a null set / a set of measure zero on the other hand are not equivalent. Neither includes the other [85]. The authors of the present paper assume that it is an open question whether for subsets of the boundary $\partial\mathcal{D}_{X}$ equivalence can be assumed since it would require more knowledge about possible fractal structures. Below, we propose definitions of generalized domains of attraction in time-delay systems. Thereby, we use the concept of prevalence by Hunt [86] to address the lack of a Lebesgue measure in infinite dimensional spaces (Remark 2.2).

Definition 8.1 (Generalized Concepts for Domains of Attraction)

i) The quasi-stability domain is defined as

[TABLE]

where $\overline{\mathcal{D}_{X}}$ denotes the closure of the domain of attraction.

ii) An almost-everywhere stability domain is any set $\mathcal{D}_{X}^{ae}\supseteq\mathcal{D}_{X}$ with

[TABLE]

where a.e. is for infinite dimensional spaces understood in the sense of prevalence [86].

Obviously, these generalizations are not smaller than the actual domain of attraction, i.e. $\mathcal{D}_{X}\subseteq\mathcal{D}_{X}^{qs}$ as well as $\mathcal{D}_{X}\subseteq\mathcal{D}_{X}^{ae}$. We denote the corresponding radius of attraction by

[TABLE]

with $gen\in\{qs,ae\}$.