Optimization along Families of Periodic and Quasiperiodic Orbits in Dynamical Systems with Delay

TL;DR

This paper extends a continuation-based optimization method to find and analyze periodic and quasiperiodic solutions in delay-differential equations, including complex cases like quasiperiodic tori, using boundary-value problem formulations.

Contribution

It introduces a novel approach for optimizing solutions along families of periodic and quasiperiodic orbits in delay systems, including the quasiperiodic case, which is a new contribution.

Findings

Successfully applied to delay-differential equations using coco software

Demonstrated optimization on quasiperiodic invariant tori

Handled complex boundary-value problems with changing segment structures

Abstract

This paper generalizes a previously-conceived, continuation-based optimization technique for scalar objective functions on constraint manifolds to cases of periodic and quasiperiodic solutions of delay-differential equations. A Lagrange formalism is used to construct adjoint conditions that are linear and homogenous in the unknown Lagrange multipliers. As a consequence, it is shown how critical points on the constraint manifold can be found through several stages of continuation along a sequence of connected one-dimensional manifolds of solutions to increasing subsets of the necessary optimality conditions. Due to the presence of delayed and advanced arguments in the original and adjoint differential equations, care must be taken to determine the degree of smoothness of the Lagrange multipliers with respect to time. Such considerations naturally lead to a formulation in terms of multi-segment boundary-value problems (BVPs), including the possibility that the number of segments may change, or that their order may permute, during continuation. The methodology is illustrated using the software package coco on periodic orbits of both linear and nonlinear delay-differential equations, keeping in mind that closed-form solutions are not typically available even in the linear case. Finally, we demonstrate optimization on a family of quasiperiodic invariant tori in an example unfolding of a Hopf bifurcation with delay and parametric forcing. The quasiperiodic case is a further original contribution to the literature on optimization constrained by partial differential BVPs.

Figures

Equations

$$\dot{z}=-z-z\left(t-1\right)+\text{cos}\omega t,$$

$$z(t)=r(\omega )\cos (\omega t-\theta (\omega )),$$

$$r(\omega )={\left[2+{\omega }^2-2\omega \sin \omega +2\cos \omega \right]}^{-1/2}$$

$$\cos \theta (\omega )=-\frac{1+\cos \omega }{r^3(\omega )},\sin \theta (\omega )=\frac{\sin \omega -\omega }{r^3(\omega )}.$$

$$\text{mbox}maximize{\mu }_A=x(0)$$

$$x^{\prime }$$

$$\tau \in (0,1/T)$$

$$x^{\prime }$$

$$\tau \in (1/T,1)$$

$$0$$

$$L\left(x(\cdot ),\varphi ,T,{\mu }_A,{\lambda }_1(\cdot ),{\lambda }_2,{\lambda }_3,{\eta }_A\right)={\mu }_A+{\eta }_A\left(x\left(0\right)-{\mu }_A\right)$$

$$+\underset{0}{\overset{1/T}{\int }}{\lambda }_1\left[x^{\prime }+T\left[x+x\left(\tau +1-1/T\right)-\cos \left(2\pi \tau +\varphi \right)\right]\right]\mathrm{d}\tau$$

$$+\underset{1/T}{\overset{1-2/T}{\int }}{\lambda }_1\left[x^{\prime }+T\left[x+x\left(\tau -1/T\right)-\cos \left(2\pi \tau +\varphi \right)\right]\right]\mathrm{d}\tau$$

$$+\underset{1-2/T}{\overset{1-1/T}{\int }}{\lambda }_1\left[x^{\prime }+T\left[x+x\left(\tau -1/T\right)-\cos \left(2\pi \tau +\varphi \right)\right]\right]\mathrm{d}\tau$$

$$+\underset{1-1/T}{\overset{1}{\int }}{\lambda }_1\left[x^{\prime }+T\left[x+x\left(\tau -1/T\right)-\cos \left(2\pi \tau +\varphi \right)\right]\right]\mathrm{d}\tau$$

$$+{\lambda }_2\left(x\left(0\right)-x\left(1\right)\right)+{\lambda }_3\left(x\left(0\right)+x\left(1-1/T\right)-\text{cos}\left(\varphi \right)\right),$$

$$-{\lambda }_1^{\prime }+T{\lambda }_1+T{\lambda }_1\left(\tau +1/T\right)=0$$

$$-{\lambda }_1^{\prime }+T{\lambda }_1+T{\lambda }_1\left(\tau -1+1/T\right)=0$$

$$0$$

$$0=\underset{0}{\overset{1}{\int }}T{\lambda }_1\sin \left(2\pi \tau +\varphi \right)\mathrm{d}\tau +{\lambda }_3\sin \left(\varphi \right)$$

$$0$$

$$+\underset{1/T}{\overset{1}{\int }}{\lambda }_1\left(x(\tau -1/T)+x^{\prime }(\tau -1/T)/T\right)\mathrm{d}\tau$$

$$+\underset{0}{\overset{1}{\int }}{\lambda }_1\left(x-\cos (2\pi \tau +\varphi )\right)\mathrm{d}\tau +{\lambda }_3{T}^{-2}{x}^{\prime }(1-1/T),$$

$$1-{\eta }_A=0.$$

$$\delta =x\left(0\right)+x\left(1-1/T\right)-\cos \varphi$$

$${\mu }_A=\left(\cos \varphi +\cos (\omega +\varphi )+\omega \sin \varphi \right)/r(\omega ),$$

$$\delta =\omega \left(\sin \varphi +\sin (\omega +\varphi )-\omega \cos \varphi \right)/r(\omega )$$

$$\dot{z}=f\left(t,z,z\left(t-\alpha \right),p\right),$$

$$x^{\prime }$$

$$x\left(0\right)-x\left(1\right)=0.$$

$$2\alpha /T<\beta <1-\alpha /T.$$

$$L_{\mathrm{BVP}}\left(x(\cdot ),\alpha ,T,p,{\lambda }_f(\cdot ),{\lambda }_{\mathrm{bc}}\right)=$$

$${\lambda }_{\mathrm{bc}}^{\text{T}}\left(x(0)-x(1)\right)+\underset{0}{\overset{\alpha /T}{\int }}{\lambda }_f^{\text{T}}\left(x^{\prime }-T{f}_1\right)\mathrm{d}\tau$$

Full text

1 Introduction

The optimization of time-delay systems has been the subject of intensive research for many years. Such systems arise naturally in control applications where unmodeled actuator dynamics results in delays between input signals and actuator responses sieber2008tracking , car following models that account for driver reaction times orosz2005bifurcations , and machine tool dynamics due to the regenerative effect tlusty1963stability . The wide range of applications has motivated the development of novel techniques for their optimization. For example, Göllmann et al. gollmann2009optimal used a formulation based on the Pontryagin minimum principle to derive necessary optimality conditions for optimal control problems with delays in state and control variables. The obtained equations were discretized and transformed into a large-scale nonlinear programming model, which was then solved using off-the-shelf solvers. In another investigation, Yusoff and Sims yusoff2011optimisation combined the semi-discretization method insperger2002semi for time-periodic delay equations with differential evolution to optimize a variable helix/pitch tool geometry for regenerative chatter mitigation. Their results were also validated experimentally, confirming the predicted significant improvements in chatter stability. This problem of optimal selection of parameters for subtractive manufacturing was also reported in iglesias2019optimum ; iglesias2014optimisation ; wojciechowski2017optimisation . Liao et al. liao2014nonlinear developed an optimization technique for periodic solutions of delay differential equations using the harmonic balance method and continuation techniques. They posed an amplitude optimization problem subject to the algebraic constraints obtained by substitution of a truncated Fourier representation in the governing equation along with the stability conditions. The sensitivity expressions were analytically derived, and the optimization problem was then solved for the unknown Fourier coefficients and the unknown parameters. The delayed Duffing oscillator was used to validate the methodology.

The calculus of variations serves as a useful tool for constrained optimization problems. Here, a Lagrangian functional is constructed by combining the objective function with the imposed constraints using Lagrange multipliers (adjoint variables) as coefficients. The vanishing of the variations of the Lagrangian with respect to the design variables and the Lagrange multipliers then yields the necessary optimality conditions for a stationary point. In general, these equations cannot be solved directly. Instead, nonlinear solvers may be applied to various finite-dimensional discretizations. A major challenge with this approach is the selection of a good initial guess which converges to the desired solution. A resolution built on principles of parameter continuation was originally proposed in the work of Kernévez and Doedel kernevez1987optimization . There, a sequence of properly initialized stages of continuation along one-dimensional manifolds of solutions to a subset of the necessary optimality conditions was used to connect the local extremum to an initial solution guess with vanishing Lagrange multipliers. This methodology was recently revisited by Li and Dankowicz li2017staged and there cast in terms of partial Lagrangians relevant to the general context of constrained optimization of integro-differential boundary-value problems without delay. Importantly, this work showed how the Lagrangian structure was consistent with a staged construction paradigm implemented in the software package coco.

In this work, we generalize the successive continuation approach of Kernévez and Doedel to optimization along families of periodic and quasiperiodic orbits in dynamical systems with delay. We derive the necessary optimality conditions from a suitably constructed Lagrangian without first discretizing the governing equations and unknowns. This approach is in contrast to other studies rubino2018adjoint , in which the discretization of the governing equations is first carried out and then the Lagrangian is constructed based on the discretized equations. In our formulation, the Lagrange multipliers satisfy coupled, piecewise-defined, boundary-value problems with both delayed and advanced arguments. Depending on the imposed constraints, the Lagrange multipliers may be discontinuous or nonsmooth at the interval boundary points, naturally resulting in a multi-segment problem calver2017numerical .

We first motivate our interest and approach with the problem of optimization of the response amplitude of a harmonically-forced, scalar, linear, delay-differential equation in Sect. 2. The general framework for problems with single delays is then considered in Sect. 3, first for periodic orbits and subsequently for families of two-dimensional quasiperiodic invariant tori. As discussed in detail, the latter optimization problem falls into the category of constrained optimization for partial differential equations (PDEs) matthias ; hinze2008optimization ; olikara2016computation , for which the necessary optimality conditions take the form of coupled, piecewise-defined PDEs with non-local coupling, as well as associated boundary and interval conditions representing periodicity in one dimension and rotation in the other. Subsections of Sect. 3 consider example applications to the search for a saddle of the response amplitude of a harmonically-forced Duffing oscillator subject to delayed feedback control and a geometric fold along a family of quasiperiodic trajectories for constant rotation number. Analysis using the coco software package validates the successive continuation approach, as well as the simultaneous discretization of the dynamic constraints and adjoint equations. A number of additional considerations and opportunities for future work are considered in the concluding section.

2 Motivating Example

We illustrate the general framework for optimization along families of solutions to delay-differential equations (DDEs) by first considering periodic responses $z(t)$ of frequency $\omega$ for a harmonically-forced, scalar, linear, delay-differential equation

[TABLE]

where we omit (here, and throughout the paper) functional arguments when they are obvious from the context. It follows from the method of undetermined coefficients that such responses are of the harmonic form

[TABLE]

where

[TABLE]

and

[TABLE]

Let us consider the optimization problem of finding the forcing frequency $\omega$ for which such a periodic response has maximum amplitude. It follows from (3) that the maximum amplitude $r_{\mathrm{crit}}\doteq r(\omega_{\mathrm{crit}})\approx 0.89$ is achieved for $\omega_{\mathrm{crit}}\approx 1.72$ (cf. Fig. 1), and that $z(t_{\mathrm{crit}})=r_{\mathrm{crit}}$ at time $t_{\mathrm{crit}}\doteq\theta(\omega_{\mathrm{crit}})/\omega_{\mathrm{crit}}\approx 2.24$ (up to multiples of the period $T_{\mathrm{crit}}\doteq 2\pi/\omega_{\mathrm{crit}}\approx 3.65$). Hence, for this simple optimization problem all components of the solution are known exactly, enabling a comparison with the results of numerical algorithms.

2.1 Formulation as a constrained optimization problem

We transform the above optimization problem into a format suitable for a general numerical solver by introducing the excitation period $T=2\pi/\omega$ as an unknown ($T$ replaces $\omega$) and rescaling time (calling the new time $\tau$) such that $x(\tau)\doteq z(T\tau+T\phi/2\pi)$. Here, the free phase $\phi$ is to be chosen so as to shift the time on the interval $[0,1]$ when the periodic solution $x$ has a critical point to $\tau=0$. Thus, we are seeking a solution to the constrained optimization problem

[TABLE]

with respect to a continuous function $x$ on $[0,1]$, as well as the variables $T$ and $\phi$, subject to the equality constraints

[TABLE]

Here, the constraints (6) and (7) impose the original delay-differential equation on the interval $(0,1)$. They rely on periodicity to wrap the delayed argument back into this interval assuming that $T>1$. The constraints (8) and (9) are boundary conditions. Constraint (8) imposes periodicity also on the interval boundary, while (9) is a phase condition that ensures that $x^{\prime}(0)=0$, consistent with $x$ having a critical point at $\tau=0$ and justifying the maximization of $x(0)$ as a substitute for the amplitude. By continuity of $x$ on $[0,1]$ and (8) it follows that $x$ is, in fact, a smooth function on $[0,1]$. Indeed, from the explicit solution in the previous section, it follows that $x(\tau)=r(2\pi/T)\cos 2\pi\tau$ and $\phi=\theta(2\pi/T)$ and, in particular, that optimality is obtained for $x(\tau)=x_{\mathrm{crit}}(\tau)\doteq r_{\mathrm{crit}}\cos 2\pi\tau$ and $\phi=\phi_{\mathrm{crit}}\doteq\theta(2\pi/T_{\mathrm{crit}})$ for $T=T_{\mathrm{crit}}$.

The constrained optimization problem (5)–(9) gives rise to the Lagrangian

[TABLE]

where the Lagrange multipliers are $\lambda_{1}(\tau)$ (a function on $[0,1]$) for the DDE constraints (6) and (7), $\lambda_{2}$ and $\lambda_{3}$ for the boundary conditions (8) and (9), and $\eta_{A}$ for the relationship between the fitness $\mu_{A}$ and $x(0)$ in (5).

In (2.1), the integral for the pairing between $\lambda_{1}$ and the DDE constraints has been split into $4$ parts, one for each of the intervals $(0,1/T)$, $(1/T,1-2/T)$, $(1-2/T,1-1/T)$, and $(1-1/T,1)$, reflecting different functional forms of the differential equations (6) and (7) for $x$ on $(0,1/T)$ and $(1/T,1)$, respectively, and anticipating possible discontinuities in $\lambda_{1}$ and $\lambda_{1}^{\prime}$. For example, the split at $\tau=1-1/T$ is in anticipation of a potential discontinuity of the Lagrange multiplier $\lambda_{1}$ at this instant caused by the imposition of a constraint on $x$ evaluated at this time in (9). This discontinuity implies a potential discontinuity of $\lambda_{1}^{\prime}$ at $\tau=1-2/T$. For the same reason, the appearances of $x(0)$ in (5) and (9) suggest that $\lambda_{1}(0)\neq\lambda_{1}(1)$ resulting in a potential discontinuity of $\lambda_{1}^{\prime}$ at $\tau=1-1/T$. All functions are assumed to be continuously differentiable on the partition implied by the integrals in (2.1). The ordering of the discontinuity points assumes that $T>3$ (Fig. 1 shows that the optimal $T$ is in this range).

Imposing vanishing variations of the Lagrangian $L$ with respect to variations in all its arguments recovers the original constraints (5)–(9) and the following adjoint system determining the Lagrange multipliers. Specfically, vanishing variations with respect to $x$ imply

[TABLE]

for $\tau\in\left(0,1/T\right)\cup\left(1/T,1-2/T\right)\cup\left(1-2/T,1-1/T\right)$ and

[TABLE]

for $\tau\in\left(1-1/T,1\right)$. Boundary and interface conditions for these equations are obtained by considering variations with respect to $x(0)$, $x(1/T)$, $x(1-2/T)$, $x(1-1/T)$, and $x(1)$, corresponding in that order to

[TABLE]

using the convention that $\lambda_{1}(\tau^{*})_{\pm}\doteq\lim_{{\tau\to\tau^{*}}\pm}\lambda_{1}(\tau)$ and recalling that $x(\tau)$ is continuous on $[0,1]$. Vanishing variations with respect to $\phi$ and $T$ imply the integral constraints

[TABLE]

and

[TABLE]

respectively. Finally, vanishing variation with respect to $\mu_{A}$ implies that

[TABLE]

In summary, the system of original contraints (5)–(9) and adjoint equations (11)–(20) is a nonlinear integro-differential boundary-value problem (BVP) defining the critical points of the Lagrangian $L$ and the constrained optimization problem (5)–(9).

In this example, the dimension of the manifold on which the constrained optimization problem is posed equals $1$, corresponding to the numbers of degrees of freedom of the nonlinear subsystem (5)–(9) (with variables $x$, $T$ and $\phi$). In contrast, the full system (5)–(9), (11)–(20) has no such degrees of freedom and, consequently, generically has only isolated solutions. Several properties put it beyond the reach of “off-the-shelf” BVP solvers:

It consists of differential equations on multiple intervals (thus, the problem is called a multi-segment BVP) with differential functional forms and continuous “right-hand sides”. The number and length of these intervals is strongly problem dependent, and may even change during the optimization process. 2. 2.

The differential equations evaluate their right-hand sides at times deviating from $\tau$ (delayed or advanced arguments). 3. 3.

The second point leads to nonlocal coupling across segments that is not restricted to coupling at the boundaries of the intervals. For example, (11) couples the values of $\lambda_{1}$ in $(1-2/T,1-1/T)$ to values of $\lambda_{1}$ in $(1-1/T,1)$.

On the other hand, the system (5)–(9), (11)–(20) has some additional structure that aids both in its construction and solution:

The equations are only forward coupled in that a solution to the original constraints (5)–(9) can be obtained independently of the values of the Lagrange multipliers, while a solution to the adjoint equations (11)–(20) requires knowledge of $x$, $T$, and $\phi$, and generically exists, at best, only for isolated choices of $x$, $T$, and $\phi$. 2. 2.

The adjoint equations (11)–(2.1) (thus excluding (20)) are linear and homogeneous in the Lagrange multipliers $\lambda_{j}$ ($j=1,2,3$) and $\eta_{A}$. A trivial solution of this subset of the adjoint system is therefore given by vanishing Lagrange multipliers for any $x$, $T$, and $\phi$. 3. 3.

The adjoint equation (20) is trivial both in construction and solution. Imposing its solution ($\eta_{A}=1$) on the remaining adjoint system, however, renders the latter nonhomogeneous.

This structure will also be present for more general cases than the example and can be exploited in the search for solutions, as well as to generate the adjoint equations (11)–(2.1) automatically during a staged construction of the optimization problem similar to li2017staged .

In this example, a few facts about the Lagrange multipliers may be deduced directly from the adjoint equations. It follows immediately from (14) and (15) that $\lambda_{1}$ is continuous at $\tau=1/T$ and $\tau=1-2/T$, and from (11) that $\lambda_{1}^{\prime}$ is continuous at $\tau=1/T$. Moreover, using the explicitly known solution for $x$, it follows that the Lagrange multiplier $\lambda_{3}$ must equal [math] at a local extremum. Indeed, substitution of the modified phase condition

[TABLE]

in lieu of (9) implies that

[TABLE]

where $\phi$ is implicitly determined by

[TABLE]

for $\delta\approx 0$. Implicit differentiation of both conditions with respect to the residual $\delta$ shows that the rate of change of $\mu_{A}$ with respect to $\delta$ equals [math] at $\delta=0$. This, in turn, implies that that $\lambda_{3}=0$ at an extremum, i.e., that $\lambda_{1}$ is, in fact, continuous also at $\tau=1-1/T$ and, consequently, continuously differentiable also at $\tau=1-2/T$. In contrast, $\lambda^{\prime}_{1}$ experiences a discontinuity at $\tau=1-1/T$ for nonzero $\eta_{A}=\lambda_{1}(0)-\lambda_{1}(1)$.

2.2 Simple continuation

According to the properties enumerated above, a solution to (5)–(9), (11)–(20) may be sought using a method of successive continuation kernevez1987optimization ; li2017staged with an embedded multi-segment boundary-value problem implementation that permits evaluation of the right-hand side at arguments shifted by arbitrary times. Specifically, this method overcomes the problem of initializing a nonlinear solver for the full system by defining a sequence of continuation problems with one-dimensional solution manifolds that connect an initial solution guess with Lagrange multipliers all equal to [math] with the sought critical point for which $\eta_{A}$ must equal $1$.

To this end, we consider the system given by the relationship between $\mu_{A}$ and $x(0)$ in (5), the boundary-value problem constraints (6)–(9), and the adjoint integral-differential boundary-value problem (11)–(2.1), but purposely omit the algebraic constraint (20). Although we anticipate that $\lambda_{3}$ will equal [math] throughout the analysis, we keep $\lambda_{3}$ as an unknown and monitor its value during continuation. By linearity and homogeneity of the adjoint subsystem in the Lagrange multipliers $\lambda_{j}$ and $\eta_{A}$, it follows that solutions to the full system lie on either of two one-dimensional manifolds. The first of these consists of functions $x(\tau)=r(2\pi/T)\cos 2\pi\tau$ with corresponding $T$, $\phi$, and $\mu_{A}=x(0)=r(2\pi/T)$, and with vanishing Lagrange multipliers. The second manifold consists of the periodic solution $x_{\mathrm{crit}}(\tau)=r_{\mathrm{crit}}\cos 2\pi\tau$ with corresponding $T=T_{\mathrm{crit}}$, $\phi=\phi_{\mathrm{crit}}$, and $\mu_{A}=r_{\mathrm{crit}}$, and with varying Lagrange multipliers proportional to $\eta_{A}$. The two manifolds clearly intersect at the local extremum of $\mu_{A}$ along the first manifold. The sought solution to the complete set of equations (5)–(9), (11)–(20) corresponds to the point along the second manifold where $\eta_{A}=1$.

In this example, the solutions along the first manifold are known explicitly. In other cases, an initial periodic response may be approximately obtained from the dynamically stable solution by direct simulation. Given such an initial solution guess for $x$, $T$, and $\phi$, a nonlinear solver may be employed to converge to a point on the manifold. A numerical continuation algorithm (e.g., pseudo-arclength continuation) may then be used to generate a sequence of points along the manifold, meanwhile monitoring for local extrema of $\mu_{A}$ and singular points for the system Jacobian (corresponding to branch points on the manifold). As shown above, and true also in the general case, these coincide. Using standard techniques, numerical continuation may proceed from such a branch point along the secondary manifold with the help of a candidate direction of continuation, for example, one that is i) transversal to the tangent direction to the original solution manifold and ii) in the plane spanned by the tangent directions to the two manifolds at the branch point.

Continuation using such an implementation in the coco software package COCO approximately locates an extremum (in the form of a fold point in $\mu_{A}$ along the solution manifold) at $T\approx 3.6515$ as shown in Fig. 2. Branch switching from the nearby branch point (exact coincidence is lost due to discretization) and continuation until $\eta_{A}=1$ yields the graphs of $x(\tau)$ and $\lambda_{1}(\tau)$ shown in Fig. 3. As seen in the bottom panel, $\lambda_{1}(\tau)$ is approximately continuous at $1-1/T\approx 0.73$, albeit with discontinuous derivative at this point, since $\lambda_{1}(0)-\lambda_{1}(1)=1$.

3 General Optimization Framework

In this section, we discuss the general methodology for optimization on periodic and quasiperiodic solutions $z(t)\in\mathbb{R}^{n}$ of delay-differential equations with a single delay of the form

[TABLE]

where $f:\mathbb{R}^{1}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{q}\rightarrow\mathbb{R}^{n}$ is periodic in its first argument with period $T$. Here, $\alpha$ and $p$ denote the time delay and the problem parameters (excluding $T$), respectively.

As the motivating example in the previous section illustrates, the problem Lagrangian and, by implication, the adjoint equations are linear in the Lagrange multipliers. The adjoint equations may therefore be constructed term-by-term by successively deriving the contributions from disjoint collections of constraints from the corresponding partial Lagrangians associated with a subset of the Lagrange multipliers. Until the full set of constraints has been considered, the adjoint equations are not completely known. The following subsections discuss the partial Lagrangians and the implied contributions to the adjoint equations resulting from the DDE constraints and boundary conditions associated with periodic and quasiperiodic orbits. For particular examples, we indicate the additional contributions associated with problem-specific constraints that complete the construction of the adjoint equations. In all cases, the contribution from the objective function to the Lagrangian implies the algebraic adjoint condition that the corresponding Lagrange multiplier ($\eta_{A}$ in the previous section) must equal $1$ at a stationary point.

3.1 Periodic orbits

Suppose first that $T>\alpha$ and consider the problem of optimizing a scalar-valued objective functional on a family of continuous solutions $x(\tau)$ to the differential equations

[TABLE]

and the boundary conditions

[TABLE]

By a rescaling of the independent variable by $T$, such solutions correspond to periodic solutions of (24) with period $T$. By continuity and periodicity, such solutions must be continuously differentiable to all orders.

Suppose, in fact, that $T>3\alpha$ and that the objective functional and any additional constraints depend on pointwise values of $x(\tau)$ only at $\tau=0$, $\tau=1$, and $\tau=\beta$ for some $\beta=\beta(\alpha,T)$ such that

[TABLE]

As we show below, such dependence results in an additional adjoint equation associated with variations with respect to $x(\beta)$. Other pointwise dependencies of the objective functional would be treated similarly, while dependence on an integral over the entire interval $[0,1]$ of a function of $x$ would not result in additional adjoint equations. We may formulate a corresponding partial Lagrangian

[TABLE]

where $f_{j}(\tau)=f\left(T\tau,x\left(\tau\right),x\left(\tau+j-\alpha/T\right),p\right)$. Here, $\lambda_{f}(\tau)$ and $\lambda_{\mathrm{bc}}$ are the Lagrange multipliers associated with the imposition of the differential equations and boundary conditions, respectively, and each integrand is assumed to be continuously differentiable on the corresponding interval. The splitting of the integral is here motivated by an anticipated discontinuity of $\lambda_{f}$ at $\tau=\beta$ and, consequently, of $\lambda_{f}^{\prime}$ at $\tau=\beta-\alpha/T$, the different functional forms of the original DDEs on the intervals $(0,\alpha/T)$ and $(\alpha/T,1)$, and an anticipated discontinuity in $\lambda_{f}^{\prime}$ also at $\tau=1-\alpha/T$.

By the stated assumptions on the objective function and any additional constraints, it is easy to show that, at a stationary point of the total Lagrangian, $\lambda_{f}(\tau)$ must be continuous at $\tau=\alpha/T$, $\tau=\beta-\alpha/T$, and $\tau=1-\alpha/T$. Using the notation

[TABLE]

for $j=0,1$ and $q=\alpha,T$ ($\partial_{k}f$ is the partial derivative of $f$ with respect to its $k$th argument, $\mathrm{d}/\mathrm{d}q$ is the total derivative of an expression with respect to $q$), the contributions to the necessary adjoint conditions for a stationary point of the total Lagrangian are given by

[TABLE]

for variations with respect to $x(\tau)$ on $\tau\in\left(0,\alpha/T\right)$;

[TABLE]

for variations with respect to $x(\tau)$ on $\tau\in\left(\alpha/T,\beta-\alpha/T\right)\cup\left(\beta-\alpha/T,\beta\right)\cup\left(\beta,1-\alpha/T\right)$;

[TABLE]

for variations with respect to $x(\tau)$ on $\tau\in\left(1-\alpha/T,1\right)$;

[TABLE]

for variations with respect to $x(0)$, $x(\beta)$, and $x(1)$, respectively;

[TABLE]

for variations with respect to $\alpha$;

[TABLE]

for variations with respect to $T$; and

[TABLE]

for variations with respect to $p$. The terms $f_{j,T}$ and $f_{j,\alpha}$ in (36) and (37) both contain time derivatives $x^{\prime}$ with delayed or advanced arguments, since $T$ and $\alpha$ both appear in the evaluation of $x$ in the third arguments of $f_{0}$ and $f_{1}$.

As previously anticipated, the explicit dependence of the Lagrangian on the internal state point $x(\beta)$ results in a potential discontinuity of the Lagrange multiplier $\lambda_{f}(\tau)$ at $\tau=\beta$. Continuous differentiability of $x(\tau)$ on $[0,1]$ and of $\lambda_{f}(\tau)$ on $(0,\beta-\alpha/T)$, $(\beta-\alpha/T,\beta)$, $(\beta,1-\alpha/T)$, and $(1-\alpha/T,1)$ implies that the necessary conditions for an extremum are in the form of a multi-segment boundary-value problem in a single trajectory segment for $x(\tau)$ and four coupled trajectory segments for $\lambda_{f}(\tau)$. A similar result is obtained, for example, in the limiting case when $\beta=1-\alpha/T$. This case specializes to the example discussed in the previous section, since there $\alpha=1$, $\beta=1-1/T$, and $T>3$. In contrast, when $\beta$ is either [math] or $1$, i.e., when there is no dependence of the objective function or any additional constraints on an internal point, then we obtain a single trajectory segment for $x(\tau)$ and three coupled trajectory segments for $\lambda_{f}(\tau)$ with both variables continuous throughout the interval $[0,1]$.

3.2 A Duffing oscillator with delayed PD control

As an application of the general methodology when $\beta=0$, consider the harmonically-forced Duffing oscillator with delayed state (proportional and derivative; PD) feedback given by the DDE

[TABLE]

Inspired by hu1998resonances , for fixed $\zeta$, $\mu$, $a$, $b$, and $\gamma$, we seek a delay $\alpha$ that minimizes the maximum amplitude of oscillation along a family of periodic responses of this system under variations in the excitation period $T$. Since the optimization problem involves minimizing a maximum, it corresponds to the search for a saddle point in the value of the oscillation amplitude on the two-dimensional constraint manifold.

Following Section 2.1, let $x_{1}(\tau)\doteq z(T\tau+T\phi/2\pi)$ and $x_{2}(\tau)\doteq\dot{z}(T\tau+T\phi/2\pi)$ represent the displacement and velocity, respectively, on the rescaled time interval $[0,1]$. The phase $\phi$ is again to be chosen so as to shift the time on this interval when the oscillator reaches its maximum displacement to $\tau=0$. It follows that the objective functional is given by

[TABLE]

for solutions of (25)–(27) subject to the phase condition

[TABLE]

and corresponding to the vector field

[TABLE]

where $p=\phi$.

The partial Lagrangian for the objective functional and phase condition is

[TABLE]

where $\lambda_{\mathrm{ph}}$ and $\eta_{A}$ are additional Lagrange multipliers. This partial Lagrangian adds the term $(\eta_{A},\lambda_{\mathrm{ph}})^{\text{T}}$ to the variation with respect to $x(0)$ in (35) (first term) and results in the algebraic adjoint constraint

[TABLE]

assuming no additional dependence of the problem Lagrangian on $\mu_{A}$.

Since neither the objective functional nor the additional phase condition depend on $x$ evaluated at an interior point of the interval $[0,1]$, it follows that $\lambda_{f}$ is continuous on the entire interval. This simplifies the partial Lagrangian $L_{\mathrm{BVP}}$ in (29) as the two integrals with boundary $\beta$ can be combined, and the resulting adjoint DDE contribution (33) can be applied on the combined interval $(\alpha/T,1-\alpha/T)$, such that $\lambda_{f}(\tau)$ is in fact continuously differentiable on $(0,1-\alpha/T)$. Correspondingly, the middle adjoint condition in (35) can be omitted. Moreover, like in the motivating example in Section 2.1, it is easy to see that the rate of change of $\mu_{A}=x_{1}(0)$ with respect to $\delta=x_{2}(0)$ vanishes at $\delta=0$. We conclude that $\lambda_{\mathrm{ph}}=0$ at a stationary point of the Lagrangian. This implies that $\lambda_{f,2}(1)=\lambda_{f,2}(0)$ and, by inspection of (33) and (34), that both components of $\lambda_{f}$ are actually continuously differentiable throughout the interval $[0,1]$.

Since the dimension $n_{\mathrm{opt}}$ of the optimization manifold equals $2$, the successive continuation approach proposed by Kernévez and Doedel kernevez1987optimization requires multiple stages (in contrast to the motivating example in Section 2.1, where $n_{\mathrm{opt}}=1$): one initially optimizes only with respect to one variable, following a curve in the optimization manifold, keeping $n_{\mathrm{opt}}-1$ variables fixed. At each successive stage of continuation one releases one further optimization variable, until all variables are free. In this analysis, we propose to keep $\alpha$ fixed during the initial stage of continuation, corresponding to the imposition of a constraint on the set of unknowns. To this end, we consider the additional partial Lagrangian

[TABLE]

This partial Lagrangian adds the constraint

[TABLE]

and the algebraic adjoint equation (for vanishing variation with respect to $\mu_{\alpha}$)

[TABLE]

and adds $\eta_{\alpha}$ to the adjoint variations with respect to $\alpha$ in (36). The total problem Lagrangian is now given by

[TABLE]

The necessary conditions for an extremum of the total Lagrangian are then given by (i) the original differential equations and boundary conditions, (25)–(27), (40), (41), and (46), and (ii) the various adjoint equations, including (44) and (47), assembled in stages as constraints and variables are added, setting the sums of all resulting contributions equal to [math]. Although we anticipate that $\lambda_{\mathrm{ph}}$ will equal [math] throughout the analysis, we keep $\lambda_{\mathrm{ph}}$ as an unknown and monitor its value during continuation.

As in the previous section, we may locate an extremum of $L$ by several successive stages of continuation, in each stage omitting one or both of the adjoint conditions (44) and (47). In particular, by holding $\mu_{\alpha}$ fixed and letting $\eta_{A}$ vary freely, we may arrive at a solution with $\eta_{A}=1$ in two stages of continuation: first, by continuing along a one-dimensional manifold with vanishing Lagrange multipliers, and next by branch-switching at a local extremum of $\mu_{A}$ to a secondary branch along which only the Lagrange multipliers vary and, in fact, do so proportionally to $\eta_{A}$. A final stage of continuation then proceeds from the point on this second manifold where $\eta_{A}=1$, but this time with $\eta_{A}$ fixed at $1$ and $\mu_{\alpha}$ free to vary. A sought extremum is obtained when $\eta_{\alpha}=0$.

An example of such an analysis for the case when $\zeta=0.05$, $\mu=0.05$, $a=0.05$, $b=-0.05$, and $\gamma=0.5$ is shown in Fig. 4 (projected into the $(\alpha,2\pi/T,\mu_{A})$ space). Here, the full integro-differential boundary-value problem is discretized and analyzed using the coco COCO package following the methodology discussed in dankowicz2013recipes in terms of continuous, piecewise-polynomial approximants on a uniform partition of every solution segment into $N=10$ mesh intervals, resulting in a large system of nonlinear algebraic equations. The successive continuation approach then proceeds along the following stages:

• •

Initial guess. An initial solution guess for $x(\tau)$ near the first manifold is first constructed using direct simulation with $\alpha=0.1$ and $T=2\pi$, after which $\phi$ is adjusted such that the maximum of $x_{1}$ occurs at $\tau\approx 0$. We finally let $\mu_{A}=x_{1}(0)$ and $\mu_{\alpha}=\alpha$.

• •

Stage 1: Continuation along manifold with vanishing Lagrange multipliers. The delay $\alpha$ is held constant by fixing $\mu_{\alpha}$ at its initial value. Continuation proceeds along the blue curve in Fig. 4, monitoring for branch points (coincident with extrema in $\mu_{A}$ up to discretization errors).

• •

Stage 2: Continuation along manifold with varying Lagrange multipliers. Branch off at the discovered branch point (labeled BP in Fig. 4) with $\mu_{\alpha}$ still fixed, stopping when $\eta_{A}$ reaches $1$. During this continuation all primary variables $x(\cdot)$, $\phi$, $T$, $\alpha$ stay constant. Only Lagrange multipliers change their values. This continuation does not change any coordinates in Fig. 4 (we remain at the point BP).

• •

Stage 3: Continuation with varying $\mu_{\alpha}$. Fix $\eta_{A}$ at $1$ and allow $\mu_{\alpha}$ (and, consequently, $\alpha$) to vary. Continue while monitoring $\eta_{\alpha}$ for zero crossings (red curve in Fig. 4). The point where $\eta_{\alpha}=0$ along the red curve is labelled “Local Optimum”. At this point all necessary conditions for a stationary point of $L$ are satisfied, including $\eta_{A}=1$ and $\eta_{\alpha}=0$.

The end point of stage 3 corresponds to a critical point at $\alpha\approx 0.7824$, $\phi\approx 1.488$, and $T\approx 5.88$ (which Fig. 4 confirms to be a saddle point). We may compare the resulting optimal delay with the prediction from a first-order multiple-scales perturbation analysis for the weakly nonlinear (small $\mu$), weakly damped (small $\zeta$), and weakly forced (small $\gamma$) case, which predicts a maximal (with respect to $T$) response amplitude

[TABLE]

(independent of $\mu$, see the appendix for intermediate steps and hu1998resonances ; rusinek2014dynamics ). The computed optimal delay $\alpha\approx 0.7824$ is in close agreement with the predicted optimal delay $\pi/4\approx 0.7854$ obtained from (49) for the case that $b=-a$. The optimal displacement profile $x_{1}(\tau)$ and the components of $\lambda_{f}(\tau)$ are shown in Fig. 5. The top panel shows close agreement between the results obtained using continuation and the harmonic response obtained from the perturbation analysis, at the computed optimal values of $\alpha$, $T$, and $\phi$. Panel (b) of Fig. 5 shows the functional Lagrange multipliers $\lambda_{f}$, confirming that they are approximately smooth in this example (since the objective does not depend on $\beta\in(0,1)$) but with $\lambda_{f,1}(1)\neq\lambda_{f,1}(0)$ and $\lambda_{f,2}(1)\approx\lambda_{f,2}(0)$ (since the objective functional and the phase constraint depend on $x(0)$ and $\lambda_{\mathrm{ph}}\approx 0$).

Further comparisons between the results obtained using the successive continuation approach and those predicted by the perturbation analysis are shown in Figs. 6 and 7 for the case when the oscillator is only subjected to displacement feedback, i.e., when $b=0$, with weak ($\mu=0.05$) and strong ($\mu=1$) nonlinearity, respectively. In each case, the perturbation analysis predicts a saddle in the response amplitude for $\alpha=\pi/2\approx 1.5708$, while the computational results are $\alpha\approx 1.4712$ and $0.8712$, respectively. For the case of weak nonlinearity depicted in Fig. 6, there is still close agreement between the optimal time histories for $x_{1}(\tau)$, while this is no longer true for the case of strong nonlinearity shown in Fig. 7. The frequency-response curves shown in the lower panels of Figs. 6 and 7 were obtained using numerical continuation for the computed and predicted critical values of $\alpha$. In the case of the weak nonlinearity, we note a weak dependence on the location and magnitude of the peak on the value of the delay, while the differences are stark in the case of the strong nonlinearity. In the latter case, the optimal delay predicted by the perturbation analysis produces a peak amplitude more than $50\%$ larger than that obtained using the numerical method.

3.3 Quasiperiodic Orbits

We proceed to consider the problem of optimizing a scalar-valued objective functional on a family of quasiperiodic solutions of (24), for which there exists an irrational rotation number $\varrho$ and a smooth function $Z:\mathbb{S}\times\mathbb{S}\rightarrow\mathbb{R}^{n}$ (here, $\mathbb{S}$ denotes the unit circle) such that

[TABLE]

in terms of the period $T$ of the vector field $f$ in its first argument. Let subscripts ${}_{\theta_{1}}$ and ${}_{\theta_{2}}$ denote partial derivatives with respect to the corresponding arguments. Substitution into the governing equation then yields the partial differential equation (PDE)

[TABLE]

on the two-dimensional torus $\mathbb{S}\times\mathbb{S}$.

We decompose this PDE along its characteristics. To this end, consider the continuous function $V:\mathbb{S}\times[0,1]\rightarrow\mathbb{R}^{n}$ given by

[TABLE]

such that $\tau=t/T$, $\theta_{1}(0)=\varphi$, and without loss of generality $\theta_{2}(0)=0$. Shifting and wrapping of arguments between and along characteristics will occur several times below. To simplify notation, suppose that $T>\alpha$ and introduce the wrapping operation $W$ for a function $V$ on $\mathbb{S}\times[0,1]$ as

[TABLE]

for $\tau-a+j\in[0,1]$ and all $\varphi\in\mathbb{S}$. It follows by periodicity that

[TABLE]

for $\tau-a+j\in[0,1]$ and all $\varphi\in\mathbb{S}$. Differentiation and use of (51) then implies that

[TABLE]

along with the boundary conditions

[TABLE]

Equations (55)–(57) are a family of coupled DDE BVPs in time $\tau$, parametrized by the continuous periodic angle $\varphi$. A family of orbit segments $\mathbb{S}\times[0,1]\ni(\varphi,\tau)\mapsto V(\varphi,\tau)\in\mathbb{R}^{n}$ solving this family of BVPs then spans the sought quasiperiodic invariant torus. Such a family is unique only up to a shift of its argument $\varphi\in\mathbb{S}$. We isolate a locally unique solution by introducing the integral phase condition

[TABLE]

in terms of a given continuously-differentiable reference function $V^{*}:\mathbb{S}\rightarrow\mathbb{R}^{n}$ that is either fixed throughout the analysis or updated as appropriate. For fixed values of the problem delay $\alpha$, excitation period $T$, and problem parameters $p$, the resultant integro-differential BVP (55)–(58) defining the quasiperiodic response is over-determined (recall that the rotation number $\varrho$ is fixed) such that one has to leave at least one system parameter free to vary to obtain isolated solutions. For example, for fixed $\alpha$, we thus expect to obtain a one-dimensional manifold of quasiperiodic invariant tori under simultaneous variations in $T$ and a single element of $p$.

We now apply the construction of the Lagrangian and adjoint equations to this family of DDE BVPs to formulate optimization problems with constraints of the form (55)–(58), following the procedure from section 3.1. We assume that neither the objective functional nor any additional constraints depend on $V$ evaluated for $\tau$ on the interior of the interval $[0,1]$, and that they only depend on $V$ on the boundaries $\tau=0$ and $\tau=1$ through integrals over $\varphi$. In this case, the Lagrange multipliers $\lambda_{f}$ for the DDE constraint (55) will be continuous on the domain $\mathbb{S}\times[0,1]$ (including periodicity in their first argument $\varphi$). The partial Lagrangian for the constraints (55)–(58) is then given by

[TABLE]

where we abbreviate $\int_{\varphi}=\int_{0}^{2\pi}\mathrm{d}\phi$ and, similarly to section 3.1, let $f_{j}(\varphi,\tau)=f(T\tau,V,W^{j}_{\alpha/T}V,p)$. The vector-valued functions $\lambda_{f}(\varphi,\tau)$ and $\lambda_{\mathrm{rot}}(\varphi)$, and the scalar $\lambda_{\mathrm{ph}}$ are the Lagrange multipliers associated with the imposition of the differential equations (55) and (56), boundary conditions (57), and the integral phase condition (58), respectively. Each integrand is assumed to be continuously differentiable on the corresponding interval, and $\lambda_{f}$ and $\lambda_{\mathrm{rot}}$ are assumed to be continuous and, hence, periodic in $\varphi$ for all $\tau$. It is again straightforward to show that $\lambda_{f}$ must be continuous in $\tau$ on $\tau=\alpha/T$ and $\tau=1-\alpha/T$ at a stationary point of the total Lagrangian. In this case, $\lambda_{f}$ is continuously differentiable in $\tau$ everywhere except at $\tau=1-\alpha/T$, where a slope discontinuity is anticipated from the boundary conditions (57).

Analogously to Section 3.1, consider the notation

[TABLE]

for $j=0,1$ and $q=\alpha,T$. Then, the contributions to the necessary adjoint conditions for a stationary point of the total Lagrangian are given by

[TABLE]

for variations with respect to $V(\varphi,\tau)$ on $(\varphi,\tau)\in\mathbb{S}\times\left(0,\alpha/T\right)$;

[TABLE]

for variations with respect to $V(\varphi,\tau)$ on $(\varphi,\tau)\in\mathbb{S}\times\left(\alpha/T,1-\alpha/T\right)$;

[TABLE]

for variations with respect to $V(\varphi,\tau)$ on $(\varphi,\tau)\in\mathbb{S}\times\left(1-\alpha/T,1\right)$;

[TABLE]

for variations with respect to $V(\varphi,0)$ on $\varphi\in\mathbb{S}$;

[TABLE]

for variations with respect to $V(\varphi,1)$ on $\varphi\in\mathbb{S}$;

[TABLE]

for variations with respect to $\alpha$;

[TABLE]

for variations with respect to $T$; and

[TABLE]

for variations with respect to $p$.

3.4 A Hopf unfolding with delay and forcing

Consider, for example, the problem of finding a local maximum in $\omega$ along a family of quasiperiodic invariant tori of the vector field

[TABLE]

where $r=\sqrt{u_{1}^{2}+u_{2}^{2}}$, $\alpha=1$, and $p=\omega$. Notably, an example in the tutorial for the coco trajectory segment toolbox COCO shows that no such local maximum exists when $\alpha=0$, since then $\omega T=2\pi\varrho$. In the present case, we consider the optimization problem

[TABLE]

subject to the constraints (55)–(58) (the coupled DDEs with boundary conditions and phase condition, depending on $\varphi$). The problem Lagrangian is then given by

[TABLE]

where $L_{\mathrm{BVP}}$ is given in (59) and $\eta_{\omega}$ is the additional Lagrange multiplier. The necessary conditions for an extremum of the total Lagrangian are then given by (i) the original differential equations and boundary conditions (55)–(58); (ii) the adjoint conditions (excluding (67)) obtained by appending $\eta_{\omega}$ to the variation with respect to $p$ (69) and setting all the resulting contributions equal to [math]; and (iii) the condition that $\eta_{\omega}=1$.

As in previous examples, we immediately note that $\lambda_{\mathrm{ph}}$ must equal [math] at a stationary point of the Lagrangian, since the objective function is clearly independent of the particular choice of family $(\varphi,\tau)\mapsto V(\varphi,\tau)$ selected by the phase condition. The adjoint boundary conditions (65) and (66) then imply that

[TABLE]

Moreover, direct computation using (72) and the boundary condition (57) shows that

[TABLE]

It follows from (63) and (64) that

[TABLE]

i.e., that $\lambda_{f}$ is continuously differentiable in $\tau$ on the entire interval $[0,1]$.

We proceed to locate an extremum by applying the successive continuation technique to the set of equations obtained by omitting the trivial adjoint condition that $\eta_{\omega}=1$. To this end, we approximate $V(\varphi,\tau)$, $\lambda_{f}(\varphi,\tau)$, and $\lambda_{\mathrm{rot}}(\varphi)$ by truncated Fourier series in $\varphi$ with $\tau$-dependent Fourier coefficient functions, as appropriate, approximated by continuous piecewise-polynomial interpolants on the interval $[0,1]$. Although we anticipate that $\lambda_{\mathrm{ph}}$ will equal [math] throughout continuation, we keep $\lambda_{\mathrm{ph}}$ as an unknown and monitor its value during continuation. We first continue along a one-dimensional manifold along which the Lagrange multipliers always equal [math], and then branch switch at a local maximum of $\mu_{\omega}$ to a secondary branch along which the family $V$ remains unchanged, while the Lagrange multipliers vary linearly in $\eta_{\omega}$. The solution to the necessary conditions for a local stationary point is then obtained once $\eta_{\omega}=1$ along the secondary branch.

The results of such an analysis using coco is shown in Figs. 8 and 9 for the case that $\varrho\approx 0.6618$. Here, dependence on $\varphi$ is approximated using a Fourier series truncated at the fifth harmonic corresponding to $11$ trajectory segments on the torus based at $\varphi=(i-1)/11$, $i=1,\ldots,11$. Each $\tau$ dependent Fourier coefficient is discretized using polynomials of degree $4$ on a uniform mesh with $10$ intervals. The one-dimensional family of quasiperiodic orbits in Fig. 8 along the first manifold with vanishing Lagrange multipliers indicates the existence of a local maximum in $\mu_{\omega}\approx 0.43685$ for $T\approx 5.3153$. Branch switching from the nearby branch point (as before, exact coincidence is lost due to discretization) and continuing until $\eta_{\omega}=1$ yields the approximate torus and the corresponding Lagrange multipliers $\lambda_{f}$ and $\lambda_{\mathrm{rot}}$ shown in Fig. 9.

As an aside, direct numerical simulation using initial conditions predicted by the continuation analysis suggest that quasiperiodic tori found on the lower half of the one-dimensional family shown in Fig. 8 are stable to sufficiently small perturbations, while the tori found on the upper half are unstable, with a critical loss of stability coincident with the peak value of $\mu_{\omega}$.

4 Conclusions

The various examples in previous sections illustrate the successful application to the case with single time delays of the general methodology to optimization along implicitly defined solutions to integro-differential boundary-value problems first proposed by Kernevez and Doedel kernevez1987optimization for ordinary differential equations. Here, the partial Lagrangian approach introduced in li2017staged was used to derive adjoint conditions that were linear and homogeneous in the unknown Lagrange multipliers. This allowed a search for local extrema to proceed along a connected sequence of one-dimensional manifolds of solutions to the necessary conditions for such extrema minus the trivial algebraic adjoint conditions on a subset of the Lagrange multipliers: first, along a branch with vanishing Lagrange multipliers, then switching to a branch with linearly varying Lagrange multipliers, and then along additional branches until all the previously omitted trivial algebraic adjoint conditions were satisfied.

In contrast to the case of ordinary differential equations, the presence of time delays introduces potential discontinuities that must be accounted for in any numerical solution strategy. By the properties of differential equations with time-shifted arguments, such discontinuities propagate across time, gaining an order of continuity for each iteration. Here, we have only accounted for zeroth- or first-order discontinuities in the formulation of the governing boundary-value problems. On each segment along which a function was shown to be continuously differentiable, we have approximated such a function by a continuous piecewise-polynomial function of degree $4$ in each mesh interval, ignoring continuity in the first derivative across mesh boundaries or discontinuities of order two or higher within each mesh interval. The piecewise-polynomial approximants have been used to impose a discretization of the governing differential equations at a set of collocation nodes within each interval and to evaluate functions with time-shifted arguments on the same or other intervals. Such a collocation strategy is consistent with the approach in engelborghs2002stability , and there compared to an alternative mesh strategy that depends on the delay. We have not undertaken a detailed analysis of the sensitivity of the results to the numerical mesh or polynomial degree. Notably, while we rely in this paper invariably on uniform meshes, it is common to consider adaptive meshes for which the number of intervals and their relative size may change during continuation. We leave such an implementation for future work.

In all the examples, a Lagrange multiplier associated with a phase condition was found to equal [math] on a local extremum of the corresponding Lagrangian. As stated previously, we nevertheless retained this Lagrange multiplier as an unknown and monitored its value during continuation. Experiments with the number of mesh intervals were used to determine whether this value was effectively [math] also in the computational analysis. An alternative would have been to eliminate this variable from the set of adjoint equations while simultaneously eliminating one of the adjoint conditions. In a single instance, this may indeed be useful, but when relying on a general-purpose implementation as envisioned in a planned future implementation of coco, it is better to retain the variable and use its numerical value as an indicator of the accuracy of the solution.

There are a number of directions to go in future work. These include consideration of circumstances in which the ratio $\alpha/T$ violates one or several of the inequalities assumed in the previous sections during continuation. Such violations may necessitate a piecewise definition of the Lagrangian across parameter space with different segmentations of the governing differential equations in each region. Problems with multiple delays, as well as problems with state- or time-dependent delays could also be explored as motivated by particular applications.

Acknowledgements.

We would like to thank Mingwu Li for insightful discussions during the course of this study. J.S. gratefully acknowledges the financial support of the EPSRC via grants EP/N023544/1 and EP/N014391/1, and from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement No 643073.

Appendix

We review the application of the method of multiple scales to the optimal selection of a time delay that results in a minimal peak displacement amplitude in the harmonically-forced response of a Duffing oscillator under delayed displacement and velocity feedback, as discussed in Sect. 3.1.

Consider the delay-differential equation

[TABLE]

for $0<\epsilon\ll 1$. We seek an approximate solution of the form

[TABLE]

where $T_{i}=\epsilon^{i}t$. To leading order in $\epsilon$,

[TABLE]

where $cc$ denotes complex conjugate terms. Elimination of secular terms at higher orders in $\epsilon$ then yields a set of conditions on the derivatives of the complex amplitude $A$ with respect to its arguments. In particular, if we let

[TABLE]

it follows from the first-order analysis that steady-state oscillations with angular frequency $1+\epsilon\sigma$ result provided that

[TABLE]

Elimination of $\varphi$ yields the desired, implicit, frequency-amplitude relationship

[TABLE]

from which we deduce the maximum value of $\rho$ given by

[TABLE]

obtained when

[TABLE]

In the special case that $b=-a$, the maximum value of $\rho$ achieves the local minimum $\gamma/2(\zeta+\sqrt{2}a)^{2}$ for $\alpha=\pi/4$, while for $b=0$, the local minimum $\gamma/2(\zeta+a)^{2}$ is obtained when $\alpha=\pi/2$.