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

For nonlinear time delay systems, domains of attraction are rarely studied despite their importance for technological applications. The present paper introduces a novel method to determine an upper bound on the radius of attraction by a numerical approach. 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 is undertaken with the time delay as bifurcation parameter. This analysis results in a possible improvement of the previous estimation. An example of a time-delayed swing equation demonstrates the overall methodology.


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]. Indeed, for a nonlinear system the Principle of Linearized Stability [17] 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 [6,[18][19][20][21][22][23][24]. Time delay systems require an initial function instead of the initial value which would be sufficient in ODEs. Hence, the state space, as well as the domain of attraction as a part of it, is infinite dimensional. The present paper addresses such domains of attraction for equilibria in autonomous retarded functional differential equations (RFDE).
To this end, consider the spectral abscissa, which by the Principle of Linearized Stability [17] decides about asymptotic stability and instability. It should be noted that this decisive measure 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. there is no relation to the size of the domain of attraction [25]. 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 normbased 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 [26,27].
(ii) Numeric Basin stability [28] in its original form for 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 [29] gives a generalization to time delay systems. Obviously, the scalar measure can contribute to a comparison of different system configurations as required in (R4).
(iii) The level set of a Lyapunov-Krasovskii functional [6,18] 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 [30]. 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 term which is suitable for the linearized system. The latter approach is wellknown for ordinary differential equations (ODE) [31]. However, is has to be reckoned with conservative estimations. Lyapunov-Krasovskii functionals are in general com- plicated expressions, such that the raw level 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 [32] can be used as well, since there are also Razumikhin based generalizations of LaSalle's invariance principle [33]. However, the same discussion as in (iii) holds.
(v) Bounds in a finite dimensional parameter space [34], 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 which 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 means a norm criterion, which ensures convergence to the attractor (interior of the green dashed circle with radiusŘ 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. Lyapunov-Krasovskii results can be translated to such easier to handle norm criteria, which goes along with further conservativeness [18,19].
(vii) An overestimation of the radius of attraction (interior of the red dashed circle with radiusR 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). Albeit 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). As soon as an initial state, which 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.

Notation and Preliminaries
Consider the mathematical nomenclature: open ball in a normed space X with center c ∈ X and radius r ∈ R >0 B X (r) In the following, we revisit some required preliminaries about autonomous retarded functional differential equations (RFDE). Subsequently, X denotes a normed space on the delay interval [−τ, 0], τ ∈ R >0 , e.g., the Banach space of continuous functions The notation of a state as x t ∈ X, is common. It represents a restriction of the solution function x(t) to the time horizoñ t ∈ [t − τ,t] ( Figure 2). Hence, a forward orbit of an initial function describes the set of all states which arise over the whole forward time in which the solution exists t ∈ [0,t max ). Thereby, x t (·; φ ) refers to solutions of the Cauchy probleṁ with an initial state φ ∈ X and x(t) ∈ R n ,() = d(·) dt , F : X → R n continuously differentiable, F(0 [−τ,0] ) = 0 n . For theorems about well-posedness in certain Banach spaces see, e.g., Hale and Verduyn Lunel [6] or Kharitonov [35]. We are mainly interested in autonomous differential difference equations [36] with a single discrete delay τ > 0, as a special type of RFDE (3) where f : R n × R n → R n . An equilibrium state is defined as φ e (θ ) ≡ x e : x t (·; φ e ) = φ e (·), ∀t ≥ 0 with a so-called equilibrium solution x(t; φ e ) ≡ 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 [37], 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. [17], Chapter VII, Theorem 6.8).

Definition 1.1 (Local Asymptotic Stability (LAS) in X)
The zero equilibrium of the autonomous system (3) is locally asymptotically stable and locally attractive in X, i.e.
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 ∈ [0, ∞). This is ensured by continuous dependence on initial data.

Domain of Attraction
The present paper deals with the domain of attraction D X ⊆ 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 φ e = 0 [−τ,0] . Thus, the domain of attraction collects all initial functions φ ∈ X, which 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 x t exists on t ∈ [0, ∞) and lim t→∞ x t X = 0}.

Radius of Attraction
Based on the notion of a norm ball around 0 with a strong simplification can be derived:

Definition 2.2 (Radius of Attraction, Ball of Attraction)
The radius of attraction in X of an asymptotically stable equilibrium is defined as and the ball of attraction B X (R X ) denotes the largest norm ball inside the domain of attraction D X .
The radius of attraction is schematically represented by the black circle in Figure 1. Thereby, the definition of attractivity in (6) ensures B(R X ) = / 0 by R X ≥ δ a > 0. Of course, the exact value of 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 [38]. 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 [39,40] 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 D c X to the origin. This leads to the equivalent definition R X = inf φ ∈X { φ X : x(t; φ ) → 0}, which is helpful for numeric estimations.

Definition 2.3 (Smallest overestimation of R X )
Assume Φ ⊂ X is a set of tested initial functions. The smallest overestimation of the radius of attractionR X ≥ R X is given bŷ If no divergent solution can be found, the above definition results in The following simple scalar example demonstrates how sensitive such estimations are to the choice of initial functions. Figure 3 shows some solutions, where initial functions in the domain of attraction φ ∈ D PC are printed in black. Four families of initial functions in PC([−τ, 0], R) are considered. We are interested in the radius of attraction R PC , i.e. the maximum uniform norm φ C below which all solutions still converge to the zero equilibrium. Clearly, the tested set of initial functions yields In the case of a delay value τ = 1 the family of step functions (Fig. 3c) shows the smallest attracted set and the family of decreasing linear functions with φ (0) = 0 ( Fig. 3d) is converging for all tested functions. However, in the case of τ = 5 the step functions ( Fig. 3g) even show the largest attracted set, whereas the decreasing linear functions (Fig. 3h) already diverge for smaller φ C values. Hence, the class of functions which 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 .

Generalizations of Domain and Radius of Attraction
The radius of attraction R X is dictated by the smallest initial function which 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 D c X is possibly only a small exception within the closure of the domain of attraction. On the one hand, R X would be unnecessarily restrictive in such applications, where a statement about generic convergence to the attractor suffices. On the other hand, such exceptions will typically also not be captured by arbitrarily chosen initial functions which are used in simulations. For ODEs, generalizations of the domain of attraction based on topological concepts like quasi-stability [26,41,42] or based on measure theoretic concepts like almost global stability [43] are well known (see Appendix). In this light, we propose definitions of generalized radii of attraction R gen X in time delay systems in the appendix. It can be expected, that upper boundsR X , which are based on arbitrarily selected initial functions, give only an estimation of the such introduced R gen X ( Figure 1).

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 analogue of the Lebesgue measure on infinite dimensional spaces (Sullivan [44], Theorem 2.38). As a consequence, there is no universally accepted notion of "almost every function" [45].

On the Selection of the State Space
There are various normed spaces, predominantly Banach spaces which come into question as state space for RFDEs. In literature, the space of continuous functions C([−τ, 0], R n ) endowed with the uniform norm is mostly used [6,10,17,36,46]. Also frequent [47][48][49][50] are spaces A L p space alone would not be sufficient, since a RFDE requires the point value x(0) = φ (0) ∈ R n at t = 0 to start with, which is a set of Lebesgue measure zero. In M p both, the initial function φ ∈ L p and the initial value φ (0) ∈ 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 choice of the normed space has far reaching consequences and should be well considered.

Consequences of the Selection
According to (6) there is a non-zero radius of attraction R X = 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 point-wise perturbations ( · C ) or address derivations in energy ( · L 2 as part of · 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 space of continuous functions C([−τ, 0], R n ) is frequently not sufficient. However, the space of piecewise continuous functions PC([−τ, 0], R n ) equipped with the uniform norm is a valid alternative for RFDEs [35,51] 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.
The radius of attraction with respect to the transformed system refers toφ (θ ). Obviously, by 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 As a consequence, the radius of attraction in M 2 of the transformed system has no meaning for the original RFDE.

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 Example 6.1 below). Hence, the history of some state variables ∈ R n does not influence the system dynamics at all, i.e.
Usually, there is no distinction in the domains of the initial function components, such Obviously, the radius of attraction as a norm criterion on φ X would be unnecessarily restrictive.
This situation motivates the use of with a certain choice of the outer norm (which is not that important by equivalence of norms in finite dimensions). However, · Q is only a seminorm in X. Therefore, we propose to introduce the quotient space endowed with the norm · Q . It regards functions to be in the same equivalence class, if they differ in values φ II (θ ) for θ < 0. Hence, these values are fully ignored in the sufficient criterion for zero-convergence φ Q < R Q .

Remark 3.2 (Differently Delayed States in Other Contexts)
Coexisting state variables with and without delay or different maximum delays are rarely addressed in literature. Gu [52] discusses multiple delay channels with focus on stability. In the context of controller synthesis, Delfour, Lee and Manitius introduce the so-called F-reduction [53,54] to overcome the problem.

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 ofR 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 parametrized 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 [55].
If application-adapted functions occur exclusively, only the corresponding subspace of X is of interest. Consequently, the intersection of the domain of attraction with this subspace 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, which does not belong to the domain of attraction, is of interest. According to Definition 2.3, it gives the desired overestimationR X ≥ R X of the radius of attraction.

Primary Initial Functions
We use the term primary initial functions to describe simple functions φ k ( · ; p), k = 1, . . . , n k like constant functions, jump functions, polynomial or trigonometric functions which are parametrized in p ∈ R n p . Denote the set of all elements in a family {φ k } p of initial functions of type k shortly as and the set of all considered functions as 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 [56] or two [55,57] parameters (p 1 , p 2 ) ∈ R 2 . In the latter case a graphical representation of the (p 1 , p 2 ) plane is common. By marking those parameter combinations which result in an attractor-convergent solution, the finite dimensional intersection Φ k ∩ 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 zeroconvergence. However, the states x t in a delay system are whole solution segments which are probably not in the chosen subset of initial functions x t ∈ Φ k . In view of avoiding unnecessary computational burden, the examination of the whole parameter space is not undertaken: to derive estimationsR X for the radius of attraction, it is sufficient to increase the parameters until zero-convergence fails. This can be done in a star-like scheme in which the parameter modulus along the rays 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. Figure 4: Intersection of the domain of attraction with certain families of primary initial functions Φ k ∩ D X . The attractor is the trivial equilibrium in (28), (29). Black pixels indicate parameter combinations which lead to zero-convergent solutions, i.e.

Remark 4.1 (Fractal Domain of Attraction Boundaries)
The boundary of domains of attraction in RFDEs is possibly fractal [55][56][57][58]. 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.2 (Higher Order Derivatives)
Assume the state space representationẋ(t) = F(x t ) with x(t) ∈ R n stems from a system description in a scalar variable y(t) ∈ R with higher order derivatives y (n) (t) + g y (n−1) (t), . . . ,ẏ(t), y(t), Then, the transformation to x, e.g. x = [y,ẏ, . . . , y (n−1) ] , goes along with restrictions to the space of initial functions. Hence, the n-dimensional vector function φ ( · ; p) = x 0 ( · ; p) cannot freely be chosen. Instead, only a scalar initial function y 0 ( · ; p) for y(t) must be specified. For instance, with the transformation from above, this results in

Secondary Initial Functions
We use the term secondary initial functions for functions which result a posteriori from previous simulations. Let T (t) : X → X, φ → T (t)φ := x t (·; φ ) denote the solution operator. As described above, the semigroup property T (s + t) = T (t)T (s), t, s ≥ 0, does not help much for considerations of Φ k ∩ D X , which are addressed in graphical representations like Figure 4. Nevertheless, it does help for considerations of the whole domain of attraction D X , which indeed is addressed by the radius of attraction. Actually, all segments of simulated solutions x(t) over an interval t ∈ [s − τ, s], i.e. states x s , can be interpreted as further initial functions φ = x s . All these intermediate results have automatically been tested for zero convergence as well, since lim t→∞ x t (θ ; φ (θ ; p)) = lim t→∞ x t (θ , x s (θ ; φ (θ ; p))).
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 for the second order Example 6.1 below, how a oscillating solution segment with a smaller norm arises from a constant primary initial function. We denote the minimum norm value of diverging trajectories (red dashed line in Figure 5) aŝ where Φ is the set of primary initial functions (20). Obviously, without a noteworthy need of additional computational effort,R T X can only be equal or smaller than the estimation gained by primary initial functionsR X . Indeed, in Example 6.1 a tighter upper  (28), (29). The constant primary initial function φ = x 0 corresponds to a red point in Figure 7 and 4a. Not φ = x 0 but a state x t , t > 0 takes the minimum norm value and thus is called secondary initial function. The components φ = [φ 1 , φ 2 ] of both are shown in subfigures, where the point φ 1 (0) fully represents the equivalence class of φ 1 (θ ) in X = Q as defined in (17).
bound on the radius of attraction than the minimum norm of all considered primary initial functions will be achieved.

Extended Primary Initial Functions
A basis extension can give further improvements of the estimation for the radius of attraction. Thereby we mean to construct with basis functions b i ( · ) ∈ X up to order d > 0. The difference to primary initial functions (Section 4.1) lies in the higher number of parameters [p 1 , . . . , p d ] =: p ∈ R n×d and thus in the computational effort. This ansatz is also taken by Leng et al. [29] 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 ∈ B d ⊂ R n×d , which 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.3 (Legitimation of polynomial approximations in the space of continuous functions)
As mentioned above, a restriction of the tested initial functionsφ to a certain subset of X like the n + 1-dimensional vector space P n ([−τ, 0], R) ⊂ C ∞ ⊂ C of polynomials with degree deg(φ ) ≤ n is convenient. Is such a procedure even theoretically capable to yield an arbitrarily good approximation of the radius of attraction?
The question implies two aspects: Firstly, assume that a set of continuous functions D C ⊆ C is given. Is it possible to derive an arbitrarily good estimation of the measure R C based on an approximation of D C ⊆ C by a set of polynomialsD n C ⊆ P n ? Polynomial approximations and in an analogous manner trigonometric approximations, are legitimated by the density of P n ([−τ, 0], R) in C([−τ, 0], R) (Stone-Weierstrass theorem). Thus, D C can be approximated in the space of polynomials byD n C ⊆ P n such that ∀ε > 0, ∃N : ∀n ≥ N : infφ Secondly, is a polynomial approximation of a zero-convergent function as well zeroconvergent? In other words, the defining property of D C shall be preserved forφ ∈ P n under the given uniform convergence ofφ → φ with increasing polynomial degree n. Indeed, this is the case if continuous dependence on the initial conditions of the RFDE is ensured (Hale and Verduyn Lunel [6], Ch.2, Theorem 2.2).

Estimations Based on Bifurcation Analysis
An invariant set with some distance to the attractor cannot be contained in the domain of attraction and thus yields a natural upper bound on the radius of attraction.
Proof. Obviously, if φ ∈ S is an element of an invariant set S with some distance to the attractor in 0 [−τ,0] , the defining property of the domain of attraction lim t→∞ x t (θ ; φ ) = 0 [−τ,0] cannot be fulfilled.
In the previously described approaches, non-zero-convergent forward orbits from time domain simulations represent this invariant set. Furthermore, these forward orbits approach ω-limit sets which are themselves invariant sets (complete continuity of F in (3) and boundedeness of trajectories presumed; Smith [59], Corollary 5.6). However, there might also be equilibria, limit cycles, invariant tori or chaotic sets, which are hidden [60] 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 [61].
(a) Initial state φ ∈ Q and limit cycle solution x(t) for τ = 20 with period T = 7.4.
(b) At delay values which mark a regain of stability, subcritical Hopf bifurcations occur. For a certain delay value, the resulting limit cycle becomes visible in the (x 1 , x 2 ) plane. Coloring of limit cycles indicates time t modulo period T . Emerging limit cycles are only shown as long as the equilibrium is attractive. Figure 6: Unstable limit cycle in (28), (29).
We propose to start with the delay-free, hence finite dimensional and manageable system, which results in setting the delay parameter to zero. All additional equilibria, limit cycles, invariant tori or chaotic sets have to emerge with an increase of the delay parameter. Hence, we consider the delay parameter as bifurcation parameter in a bifurcation analysis (for bifurcation theory of RFDEs see [6,17,[62][63][64]). 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 delay 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) which is solved by the periodic solutions of a limit cycle. Periodicity of the whole solution x(t) = x(T + t),t ∈ [−τ, ∞) is by semigroup property already ensured, if x(θ ) = x(T + θ ), θ ∈ [−τ, 0], i.e. by equality of the two states x 0 = x T . Hence, instead of an initial condition x 0 = φ the boundary condition x 0 = x T is required, such that a twopoint BVP over one period t ∈ (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ẋ = f (x(t), x(t − τ)). With t =tT and x(tT ) =: z(t) the two-point BVP over the normalized periodt ∈ (0, 1] becomes Still it is ambiguous which delay-spanning segment of the periodic orbit is addressed. This translational symmetry has to be resolved by an additional phase shift condition leading to a unique solution. Since in a bifurcation analysis limit cycles are tracked from their emergence with varying delay parameter, the problem is well initializable and numerically solvable. We make use of an implementation in the Matlab toolbox DDE-BIFTOOL [65] which is based on a collocation method with piecewise polynomials [66]. A numeric approximation of the periodic solution results. Of course, limit cycle bifurcations 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 [67]. A bifurcation analysis aims at detecting such phenomena. Approximations of Floquet multipliers allow to conclude stability. After the numeric analysis, if a certain fixed delay value is given, some occurring invariant sets can be concluded. The periodic solution x(t) yields the limit cycle orbit (set of all delay-width solution segments) which 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 (36), corresponding to the stable eigenfunctions of the linearized system [17]. 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 remain usually unknown. Nevertheless, knowledge from ODEs suggests that solutions nearby the stable manifold approach temporarily the saddle-type set in time forward simulations. Thus, the saddle-type set can to some extend 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.

Demonstrative Example
The following example demonstrates the proposed approach for estimating the radius of attraction. Second order ordinary differential equations with periodic nonlinear-ity are relevant for various non-related applications like mechanical pendulum mechanisms, Josephson junction circuits, phase-locked loops or synchronous generators [31]. At least in the latter case, the impact of delayed damping on the domain of attraction is of high technical interest [68]. Schaefer et al. [68] give a simulative analysis of basin stability for such a delayed swing equation.
The aim is to find an upper boundR X > R X for the radius of attraction R X in the chosen norm.
(i) Selection of the State Space X The history of x 1 does not contribute to the system dynamics at all. According to (18), this motivates the use of a quotient space X = Q. We define

Remark 6.1 (Definition on a Manifold)
There is a rotational symmetry in x 1 , such that instead of x(t) ∈ R 2 the instantaneous values can also be defined on a cylinder, i.e. [x 1 (t), x 2 (t)] ∈ S 1 × R and consequently 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,69]. (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 (28) implies d dt φ 1 = φ 2 (Remark 4.2). However, because of the irrelevance of the x 1 -history in the above introduced quotient space (30), this relation has no consequence for the selection of initial functions. As representative of the first component φ 1 : where p 1 ∈ R parametrizes this initial function. For the second component φ 2 : [−τ, 0] → R, we select a finite number of primary initial functions parametrized by p 2 ∈ R (Section 4.1): To simplify the evaluation, the primary initial functions above are chosen such that φ k 2 ( · ; p 2 ) C = |p 2 |, k ∈ {1, . . . , 4}. This parametrization is advantageous since · Q norms of the initial function equal Euclidean norms in the (p 1 , p 2 ) plane, i.e.
Thus, a norm ball B Q (r)={φ ∈ Q : φ Q < r} corresponds to the inner of a circle in the parameter plane.
Then · Q is replaced by · 2 , such thatR 2 = R 2 describes exactly the maximum inner circle 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 (19) is given by Figure 4 shows results in the (p 1 , p 2 ) parameter plane. To get an approximation of the relevant parts of the boundary ∂ (Φ k ∩ D Q ), 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 shows 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 boundR Thereby, the superscript k indicates thatR k Q refers to the k-th family of initial functions Φ k , k ∈ {1, . . . , 4}. All four families of initial functions taken in Figure 7 Figure 7 indicates that there are indeed diverging solutions which come even closer to the attractor than all the considered primary initial functions. Secondary initial functions of all chosen families Φ k , k ∈ {1, . . . , 4} yield similar minimum norm valuesR k,T Q . The best result is gained by a constant initial function which leads to a diverging solution with minimum Q-norm R T Q =R 1,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 (28) [70] reveals that each regain of stability is accompanied by a subcritical Hopf bifurcation. 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 τ = 20 the corresponding periodic solution in Figure 6a is remarkably similar to the state which 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([−τ, 0], R 2 ) is identical to R Q in the quotient space. However, the uniform norm of limit cycle states would lead to unnecessarily conservative estimationsR LC PC (Figure 8). The minimum Q-norm value for the limit cycle solution isR LC Q = 0.507.
Hence, for the delay value under consideration, secondary initial functions and the bifurcation analysis based approach yield almost equal upper boundsR LC Q ≈R T Q on the radius of attraction R Q ≤R Q =:R 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 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) ∈ R 2 , but to solution segments, i.e. to states x t ∈ 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 numeric 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.

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 Scalar Measures 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 measure. 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 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, which are simulatively 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 also all available solution segments from simulations into account. 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 or chaotic sets 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 of these sets, we propose a numeric bifurcation analysis. By starting from the delay-free system, a bifurcation analysis allows first to detect their emergence and then to track them numerically for increasing time delay values.
• 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 bifurcation analysis based approach is a worthwhile additional element of the new method. 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. [41], who introduce the topological concept of a quasi-stability boundary, which is defined as boundary of int(D X ). The main application is (delay-free) power system stability [26,42]. 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. [43] 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 meagre 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 [71]. The authors of the present paper assume that it is an open question whether for subsets of the boundary ∂ 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 [72] to address the lack of a Lebesgue measure in infinite dimensional spaces (Remark 2.2).

Definition 8.1 (Generalized Domain of Attraction Concepts)
i) The quasi-stability domain is defined as where D X denotes the closure of the domain of attraction. ii) An almost-everywhere stability domain is any set D ae X ⊇ D X with φ ∈ D X for a.e. φ ∈ D ae , where a.e. is for infinite dimensional spaces understood in the sense of prevalence [72].
Obviously, these generalizations are not smaller than the actual domain of attraction, i.e. D X ⊆ D qs X as well as D X ⊆ D ae X . We denote the corresponding radius of attraction by R gen X = sup{r > 0 : B X (r) ⊆ D gen X }.