A method for the optimization of nonlinear systems with delays that guarantees stability and robustness

TL;DR

This paper introduces a novel method for optimizing nonlinear delay differential systems that guarantees stability and robustness against parameter uncertainties, enabling practical and reliable system design.

Contribution

The method guarantees stability and robustness through constraints rather than cost functions, allowing simultaneous optimization of economic goals and system stability.

Findings

Optimized a laser diode to maximize intensity while ensuring stability.

Guarantees stability despite uncertain parameters within specified intervals.

Enables practical robust optimization for nonlinear delay systems.

Abstract

We present a method for the steady state optimization of nonlinear delay differential equations. The method ensures stability and robustness, where a system is called robust if it remains stable despite uncertain parameters. Essentially, we ensure stability of all steady states of the nonlinear system on the steady state manifold that results from the variation of the uncertain parameters. The uncertain parameters are characterized by finite intervals, which may be interpreted as error bars and therefore are of immediate practical relevance. Stability despite uncertain parameters can be guaranteed by enforcing a lower bound on the distance of the optimal steady state to submanifolds of saddle-node and Hopf bifurcations on the steady state manifold. We derive constraints that ensure this distance. The proposed method differs from previous ones in that stability and robustness are guaranteed with constraints instead of with the cost function. Because the cost function is not required to enforce stability and robustness, it can be used to state economic or similar goals, which is natural in applications. We illustrate the proposed method by optimizing a laser diode. The optimization finds a steady state of maximum intensity while guaranteeing asymptotic or exponential stability despite uncertain model parameters.

Figures

Tables

Table 1. Table 1. System parameters. Nominal values are taken from Verheyden et al.  [ 29 ] . The abbreviations ’s.t.’ and ’opt.’ are short for ’subject to’ and ’optimization’, respectively.

symbol	parameter	uncertainty	nom. value
$j$	pump current	$\pm 0.01$	$j^{\left(0\right)}$, s.t. opt.
$\eta$	feedback strength	$\pm 0.0001$	${\eta }^{\left(0\right)}$, s.t. opt.
${\alpha }_{LW}$	linewidth enhancement factor	$\pm 1$	as in [29]
$\varphi$	feedback phase	$\pm 0.1\pi$	as in [29]
$\tau$	feedback delay	$\pm 10$	${\tau }^{\left(0\right)}$, s.t. opt.
$T$	carrier to photon lifetime	$\pm 10$	as in [29]
$d$	diffusion constant	$\pm 0.001$	as in [29]

Equations

$$\dot{x}(t)=f(x(t),x(t-{\tau }_1),\dots ,x(t-{\tau }_m),\alpha )$$

$${\tau }_i={\tau }_i(x(t),\alpha ),i=1,\dots ,m.$$

$$\det \left(\lambda I-A_0-\sum _{i=1}^m{A}_i\exp (-\lambda {\tau }_i)\right)=0,$$

$$\mathrm{s}(\sigma ,\omega ,\tau )$$

$$\mathrm{c}(\sigma ,\omega ,\tau )$$

$$E(t)\in \mathbb{C}$$

$$\frac{\mathrm{d}A(t)}{\mathrm{d}t}=$$

$$+\eta {\mathrm{e}}^{-\mathrm{i}(\varphi -\Omega \tau )}A(t-\tau ),$$

$$T\frac{\mathrm{d}{N}_k(t)}{\mathrm{d}t}=$$

$$+P_k-F_k[1+2N_k(t)]|A(t){|}^2$$

$$T\frac{\mathrm{d}{N}_1(t)}{\mathrm{d}t}=$$

$$+P_1-F_1[1+2N_1(t)]|A(t){|}^2$$

$$T\frac{\mathrm{d}{N}_K(t)}{\mathrm{d}t}=$$

$$+P_K-F_K[1+2N_K(t)]|A(t){|}^2.$$

$$P_{k}=\left\{\begin{array}[]{ll}1.075j&\text{ for }|k\cdot h|<0.2\\ -0.8&\text{ else }\end{array}\right\},~{}k=1,\dots,K\,.$$

$$\zeta (t)=2h\sqrt{\frac{400}{\pi }}\sum _{k=1}^K{F}_k{N}_k(t)$$

$$F_k=\exp (-400(k\cdot h{)}^2),k=1,\dots ,K.$$

$$\mathrm{Re}\{A(t)\}-\mathrm{Im}\{A(t)\}=0.$$

$$f(x(t),x(t-\tau ),\alpha )=\left[\begin{array}{c}\text{eqref}eq:discElectricalField\\ ⋮\\ \text{eqref}eq:phaseCondLaser\end{array}\right]$$

$$\underset{x,j,\eta ,\tau }{\min }-|A{|}^2$$

$$\text{s.t.}0$$

$$0.02\le$$

$$0.0001\le$$

$$220\le$$

$$\mathrm{Re}\{{\lambda }_i\}<0\text{ for all }{\lambda }_i.$$

$$\mathrm{Re}\{{\lambda }_i\}<\sigma <0\text{ for all }{\lambda }_i.$$

$${\alpha }^{(0)}-{\alpha }^{(c)}+d\frac{r}{||r||}=0.$$

$${\alpha }_i\in [{\alpha }_i^{(0)}-\Delta {\alpha }_i,{\alpha }_i^{(0)}+\Delta {\alpha }_i],$$

$$\frac{{\alpha }_i}{\Delta {\alpha }_i}\in \left[\frac{{\alpha }_i^{(0)}}{\Delta {\alpha }_i}-1,\frac{{\alpha }_i^{(0)}}{\Delta {\alpha }_i}+1\right].$$

$$0=G(x^{(c)},{\alpha }^{(c)},u^{(c)}),$$

$$0=H(x^{(c)},{\alpha }^{(c)},u^{(c)},\kappa ,r),$$

$$f(x^{(c)},...,x^{(c)},{\alpha }^{(c)})$$

$$\sigma w-A_{0}w-\sum _{i=1}^m{A}_i\exp (-\sigma {\tau }_i)w$$

$$w^{\prime }w-1$$

$$\left[\begin{array}{ccc}{\nabla }_{x^{(c)}}{f}^{\prime }& B_{12}& 0\\ 0& B_{22}& 2w\end{array}\right]\kappa$$

$$\left[\begin{array}{ccc}{\nabla }_{\alpha }{f}^{\prime }& B_{32}& 0\end{array}\right]\kappa -r$$

$$r^{\prime }r-1$$

$$B_{12}=$$

Full text

1. Introduction

Delays can significantly influence the behavior of a dynamical system. In fact, delays can have both a stabilizing and a destabilizing effect [25]. Technically, it is more difficult to analyze a delayed systems, because even a single delay leads to an infinite dimensional dynamical system [7].

Delays are present in various kinds of dynamical systems. They are, for example, part of traffic models where they model reaction times of human drivers [9], they influence population dynamics in multi-stage populations [3, 1], they cause limit cycles in milling and turning processes [23, 10], and they can destabilize supply chains [24]. The delay present in external cavity lasers causes very interesting dynamical behaviour [16, 19, 8]. External cavity lasers can exhibit even quasi-chaotic dynamics. The relevant time constants are typically dominated by photon lifetimes of around 1ps [2, p. 232]. As a consequence, it is challenging to implement feedback control, and thus designs and operating modes that are intrinsically stable and robust are of interest.

Parameter variations that lead to a significant change in the behavior of semiconductor laser dynamics can be associated with the occurrence of bifurcations. While a Hopf bifurcation is desirable at the initial lasing threshold, further bifurcations might render the single-wavelength operation unstable. Since a loss of stability involves a real or conjugate pair of eigenvalues crossing the imaginary axis into the right half complex plane, it is obvious to search for laser parameters that result in modes of operations with eigenvalues with sufficiently negative real parts. Vanbiervliet et al. [28] use nonlinear programming methods to find laser diode parameters that result in the minimal (i.e., maximally negative) leading eigenvalue.

The method proposed in the present paper is also based on nonlinear programming and on characterizing stability and robustness with the distance of the leading eigenvalue to the imaginary axis. Our method is fundamentally different from the one by Vanbiervliet et al. [28], however. We do not use the cost function to enforce the desired eigenvalues, but we impose additional constraints that establish robust stability (see Section 3). Moreover, we use the distance to the closest bifurcation in the parameter space as a measure for robustness (see the next paragraph for a brief explanation). Since the cost function is not required to enforce stability, other optimization goals can be stated; in our case this is the laser intensity. Essentially, we can systematically determine the mode of the highest laser intensity that is sufficiently robust by optimizing with respect to the laser intensity and enforcing a lower bound on the distance to all stability boundaries. Note that this is different from determining the most robust mode by pushing all eigenvalues as far into the left half of the complex plane as possible.

The distance of a candidate steady state of operation to the closest bifucation manifold can serve as a measure for robustness. Dobson showed this distance occurs along a particular normal vector to the bifurcation manifold [4]. Mönnigmann and Marquardt incorporated these normal vectors to state robustness in economic optimization problems [20]. Their normal vector approach has been extended in different contexts, for example, to achieve specified transient behavior [11], to the stability of periodically operated systems [12] and, more recently, to delayed systems [13, 21, 22]. The Lang-Kobayashi type laser diode models, which serve as examples in the present paper, belong to the latter problem class.

After introducing the system class in the remainder of Section 1, we present a sample system, a laser model, and discuss the optimization task in Section 2. In Section 3, we propose normal vector constraints as an instrument for achieving robust stability and derive constraints for robust asymptotic and robust exponential stability. We apply these normal vector constraints to the laser optimization in Section 4. Conclusions and an outlook are given in Section 5.

1.1. System Class and notation

We consider the class of delay differential equations (DDEs) with multiple uncertain and state dependent delays

[TABLE]

with the state vector $x\in\mathbb{R}^{n_{x}}$, uncertain parameters $\alpha\in\mathbb{R}^{n_{\alpha}}$, where $f$ is smooth and maps from $\mathbb{R}^{n_{x}(m+1)}\times\mathbb{R}^{n_{\alpha}}$, or an open subset thereof, into $\mathbb{R}^{n_{x}}$. We assume there exist $m$ state and parameter dependent delays

[TABLE]

We refer to $\lambda\in\mathbb{C}$ as an eigenvalue at a steady state $x$ of (1), if

[TABLE]

where $A_{0}$ and $A_{i}$ are the Jacobians of the right hand side of (1) with respect to $x(t)$ and $x(t-\tau_{i})$ respectively (see, e.g., Engelborghs and Roose [6]). The delay $\tau_{0}=0$ is merely introduced to simplify the notation. This allows us to replace $x(t)$ by $x(t-\tau_{0})$ and refer to $x(t-\tau_{i})$, $i=0,\dots,m$ instead of referring to $x(t)$ and $x(t-t_{i})$, $i=1,\dots,m$ separately. Furthermore, we introduce the abbreviations

[TABLE]

A solution of a set of nonlinear equations is called regular if the Jacobian of the nonlinear equations evaluated at this solution has full rank.

2. Laser diode model and problem outline

We use a laser diode model to motivate the use of the normal vector method to be introduced in Section 3. The model was proposed by Verheyden et al. [29]. It is of the Lang-Kobayashi type and models the laser dynamics in the presence of an external cavity. The external cavity feeds a small fraction of the emitted light back into the semiconductor with a delay, thus necessitating a delay differential equation (DDE). Since the semiconductor is represented by a diffusion equation, the overall system model requires coupling a delay differential to a partial differential equation. We use a method-of-lines approximation of the partial differential equation as proposed in [29]. A technical complication arises because of the rotational symmetry of the steady state solutions of the Lang-Kobayashi model [8, 14].

2.1. Dynamical laser diode model

[TABLE]

The energy powering the laser is supplied through the semiconductor. The partial differential equation modeling the carrier diffusion can be spatially discretized with a method-of-lines approach. According to Verwer and Sanz-Serna [30], the discretization converges to the underlying PDE solution for a sufficient number of discrete points. Let $k=1,2,...,K$ refer to the points resulting from the spatial discretization. Then

[TABLE]

result for the discretized carrier density $N_{k}(t)$ in the interior ($k=2,3,\dots,K-1$), where the constant $h=\frac{0.5}{K-1}$ describes the distance of the discrete spatial points. Assuming Neumann boundary conditions, the remaining two equations for $k=1$ and $k=K$ are

[TABLE]

and

[TABLE]

The carrier diffusion is driven by an anisotropic pump current. Its discretized spatial distribution is

[TABLE]

The modal gain

[TABLE]

depends on the spatial distribution of the electrical field

[TABLE]

When stated in rotating coordinates as in (5a), the rotational symmetry of the Lang-Kobayashi model (see [14]) is evident from the fact that $A(t)\exp(i\theta)$ is a steady state solution to (5a) for any $\theta\in\mathbb{R}$, if $A(t)\in\mathbb{C}$ is a steady state solution. We follow Verheyden et al. [29] in removing this indeterminacy with the additional condition

[TABLE]

It is known that the stability properties of the steady states are not affected by this additional algebraic equation (see [8], Section II). Consequently, steady states of (1) correspond to solutions $0=f(x,x,\alpha)$ of the $n_{x}=K+3$ equations

[TABLE]

in the $n_{x}$ variables $x=[\mathrm{Re}\{A\},\mathrm{Im}\{A\},N_{1},\dots,N_{K},\Omega]^{\prime}$. The model has $n_{\alpha}=7$ parameters $\alpha=[j,\eta,\alpha_{LW},\phi,\tau,T,d]^{\prime}$. We consider all parameters to be uncertain in the sense that they are known up to an error bar only. The uncertainties or error bars are given in Table 1.

2.2. Optimization problem outline

It is the objective of the following optimization to find the parameter configuration for which the laser intensity becomes maximal. A high laser intensity is important for applications relying on the high power density of lasers such as laser cutting [17] or laser cladding [27]. We can state the task of maximizing the stationary laser intensity as

[TABLE]

The bounds on the feedback strength $\eta$ model limitations in the coating process for the semi-reflective mirror at the end of the external cavity. Bounds on $\tau$ are necessary to account for geometric limitations of the external cavity length. The pump current bounds represent a limited range of a current source.

2.3. Limitations of naive steady state optimizations of dynamical systems

2.3.1. Optimization without stability constraints result in unstable operation

The numerical solution of the the optimization problem (6) for $K=31$ results in the optimal parameters $j=1.0265$, $\eta=0.02$ and $\tau=220$. These parameters correspond to a stationary intensity of $|A|^{2}=6.6622$. While optimal, this steady state is unstable, which is evident from Figure 1.

This result illustrates an optimization such as (6) that ignores the stability properties of the dynamical system is likely to fail, because an optimal but unstable, or more generally, an optimal but not sufficiently robust mode of operation is useless. Arguably, the intensity could first be optimized, and the resulting steady state could be stabilized by adding feedback control a posteriori. For systems with very fast intrinsic time scales this is not a viable option, however. Moreover, even if such an a posteriori stabilization was possible, it would raise the question whether there exists a competitive intrinsically stable and robust mode of operation that does not require adding a controller device to the laser diode system. The remainder of the paper is devoted to an optimization method that systematically takes stability and robustness into account.

3. Normal vector method

We give an informal introduction to the simple geometric idea of the proposed method first. Subsequent sections then use more precise notions from bifurcation theory.

3.1. Normal vectors to enforce stability

A steady state is asymptotically stable, if the eigenvalues defined by (3) are all in the left half-plane, i.e.,

[TABLE]

The system loses stability if a variation of the parameters $\alpha_{i}$ causes one or more eigenvalues to cross from the left into the right half plane. Assuming no higher codimension bifurcation points appear, the critical points constitute a manifold that locally separates the parameter space into a stable and an unstable part [15] (see Figure 2). We sometimes prefer to enforce exponential stability with a prescribed decay rate $\sigma<0$ and replace (7) by

[TABLE]

In this case the existence of real eigenvalue or complex conjugate pair of eigenvalues with real part $\sigma$ defines the separating manifold. We call this manifold critical manifold in both cases, (7) and (8), by a slight abuse of terminology.

The closest distance of a point $\alpha^{(0)}$ in the parameter space to the critical manifold can serve as a measure for robustness (see Figure 2). The closest distance occurs along a direction that is normal to the critical manifold [4]. This direction and the closest critical point are labeled $r$ and $\alpha^{(c)}$, respectively, in Figure 2. It is evident from the figure that the distance $d$ between $\alpha^{(0)}$ and $\alpha^{(c)}$ respects

[TABLE]

Stability can be enforced with the constraint $d>0$ in an optimization.

Robustness can also be enforced easily as follows. We assume the uncertainty is quantified by error bars around the nominal values $\alpha^{(0)}$

[TABLE]

where $\Delta\alpha_{i}$ measures the uncertainty of the respective $\alpha_{i}$. Rescaling the parameter space by these uncertainty intervals111Alternatively, the $\Delta\alpha_{i}$ can be used to define a metric. transforms the hyperrectangular uncertainty region (10) into a hypercube,

[TABLE]

Enforcing a parametric distance of $d>\sqrt{n_{\alpha}}$ leads to a robustly stable steady state, because the uncertainty region (11) is enclosed in a hypersphere of radius $\sqrt{n_{\alpha}}$ that cannot cross the stability boundary [20]. The geometry of the hypercubic uncertainty region, the enclosing hypersphere and the stability boundary are illustrated in Figure 2 for $n_{\alpha}=2$.

Applied bifurcation theory states systems of equations that characterize critical manifolds for the case (7) in the form of augmented systems (see, e.g., [15]). They have the form

[TABLE]

where $x^{(c)}$ is a steady state on the stability boundary for the critical parameters $\alpha^{(c)}$. The vector $u^{(c)}$ collects various auxiliary variables such as eigenvectors that belong to the critical eigenvalues. We also use the term augmented systems for systems of equations used to characterize critical manifolds for (8). Section 3.2 explains how to extend these systems for the calculation of the normal vectors $r$ required for the robustness constraints introduced with Figure 2. These extensions, which can be found with a procedure stated in [20], result in systems of nonlinear equations of the form

[TABLE]

where $r\in\mathbb{R}^{n_{\alpha}}$ is the normal vector at the critical point $x^{(c)}$, $\alpha^{(c)}$ and $u^{(c)}$ from (12). The vector $\kappa$ collects additional auxiliary variables. Just as in (12), the number of equations in $0=G$ and $0=H$ depends on the type of critical point. This particular types of interest here are treated in the next section.

3.2. Derivation of normal vector systems

If a steady state belongs to a critical manifold for stability, there exists a zero eigenvalue or a complex conjugate pair on the imaginary axis. Critical manifolds for stability can therefore be characterized with the augmented systems for fold and Hopf bifurcations. Critical manifolds of steady states with a prescribed decay rate $\sigma<0$ can be characterized correspondingly, i.e., by augmented systems that state the existence of a real eigenvalue $\sigma$ or a complex conjugate pair with real part $\sigma$. Because these systems obviously resemble those for fold and Hopf bifurcations closely, we refer to them as augmented systems for modified fold and modified Hopf points, respectively.

3.2.1. Normal vectors system for fold and modified fold manifolds

Assume $x^{(c)}$ is a steady state for parameter values $\alpha^{(c)}$. If $\sigma$ is real and the leading eigenvalue at the steady state $x^{(c)}$, then there exists a $w\in\mathbb{R}^{n}$ such that $x^{(c)}$, $\alpha^{(c)}$, $a$ and $b$ obey the following equations, which are a simple extension of the standard augmented system for saddle-node bifurcations of delay differential equations (see, e.g., [5]):

[TABLE]

Equations (14a) holds, because $x^{(c)}$ is a steady state by assumption. Since $\sigma$ is a real eigenvalue by assumption, (3) holds. Consequently, there exists an eigenvector $w\in\mathbb{R}^{n}$ that satisfies (14b). Equation (14c) normalizes this eigenvector to unit length.

The following proposition extends a result stated in [22] to the case with delays that depend on uncertain parameters $\alpha$.

Proposition 1** (normal vector to manifold of modified fold points).**

Let $\sigma\leq 0$ be arbitrary but fixed. Assume that $(x^{(c)},\alpha^{(c)},w)$ is a regular solution to (14) in the $2n_{x}+1$ variables $x^{(c)},w$ and one of the elements of $\alpha^{(c)}$. Then $r$ that obeys the following equations is normal to the manifold of modified fold points at this solution:

[TABLE]

where

[TABLE]

Proof.

Consider (14) as $2n+1$ equations in the $2n+n_{\alpha}$ variables $x^{(c)}$, $w$ and $\alpha^{(c)}$. These equations define an $(n_{\alpha}-1)$-dimensional manifold of modified fold points in the neighborhood of the known regular solution. Evaluated at any point on this manifold, the rows of the Jacobian of (14) with respect to $x^{(c)}$, $w$ and $\alpha^{(c)}$ span the normal space to the manifold [20][26, Chapter 3]. Since we would like to span the normal space with column vectors, it is more convenient to work with the transposed Jacobian, which we denote $B$. Formally, $B$ can be stated as the outer product

[TABLE]

which results in the block matrix

[TABLE]

where the block columns contain $n$, $n$ and $1$ columns and the block rows contain $n$, $n$ and $n_{\alpha}$ rows, respectively. The elements in the second column of (17) remain to be determined. The block $B_{22}$ is the Jacobian of $\sigma w-A_{0}^{\prime}\,w-\sum_{i=1}^{m}A_{i}^{\prime}\exp(-\sigma\tau_{i})\,w$ with respect to $w$, which results in (16b). The block matrices $B_{12}$ and $B_{32}$ in (17) require the application of the chain and product rules, since $\tau_{k}$ is, as described in (2), a function of $x^{(c)}$ and $\alpha^{(c)}$. Therefore,

[TABLE]

which yields (16a) and $B_{32}$ as stated in (16c) results when $\nabla_{x^{(c)}}$ is replaced by $\nabla_{\alpha^{(c)}}$. The columns of $B$ span the normal space to the critical manifold in the space with the components $(x,w,\alpha)$. Because we need to measure the distance to the critical manifold in the space of the parameters $\alpha$, we seek the particular normal direction that only has nonzero elements in the directions of the parameters $\alpha$. This corresponds to the linear combination $\kappa$ with

[TABLE]

which is equivalent to (15b) and (15c). These equations determine the normal direction $r$ up to its length. The last condition (15d) fixes the length of $r$. ∎

The augmented system for a fold bifurcation manifold is the special case of (14) with $\sigma=0$, which is the well-known augmented system (see, e.g., [15])

[TABLE]

We can derive the normal vector system for the manifold defined by (18) as a corollary to Prop. 1.

Corollary 1** (normal vector to manifold of fold bifurcations).**

If $(x^{(c)}$, $\alpha^{(c)}$, $w)$ is a regular solution to (14) for $\sigma=0$ in the $2n_{x}+1$ variables $x^{(x)}$, $w$ and one of the elements of $\alpha^{(c)}$, then $r$ that obeys the following equations is normal to the manifold of fold points at this solution:

[TABLE]

The corollary follows directly from Prop. 1 by substituting $\sigma=0$.

3.2.2. Normal vector systems for Hopf and modified Hopf manifolds

This section first briefly summarizes the augmented systems for Hopf and modified Hopf points as needed for the paper. The normal vector systems are derived subsequently.

Assume $x^{(c)}$ is a steady state for parameter values $\alpha^{(c)}$. If $\lambda=\sigma\pm\mathrm{i}\omega$ are the leading eigenvalues at this steady state, then there exists a $w=a+\mathrm{i}b\in\mathbb{C}^{n}$ such that $x^{(c)}$, $\alpha^{(c)}$ and $w$ obey the equations

[TABLE]

which are a simple extension of the augmented system for Hopf bifurcations of delay differential equations (see, e.g., [5]). Equations (20a) hold, because $x^{(c)}$ is a steady state by assumption. Since $\lambda=\sigma\pm\mathrm{i}\omega$ are eigenvalues at this steady state by assumption, (3) holds. Consequently, there exists a $w=a+\mathrm{i}b\in\mathbb{C}^{n}$ such that

[TABLE]

and separating this equation into its real and imaginary part for $\lambda=\sigma+\mathrm{i}\omega$ yields (20b) and (20c), respectively. The eigenvector is normalized to unit length by (20d) and its complex phase is fixed by (20e).

The normal vector system can now be derived based on (20).

Proposition 2** (normal vector to manifold of modified Hopf points).**

Let $\sigma\leq 0$ be arbitrary but fixed. Let $\omega$, $a$ and $b$ be as in (20). Assume that $(x^{(c)},\alpha^{(c)},\omega,a,b)$ is a regular solution to (20) in the $3n_{x}+2$ variables $x^{(c)},\omega,a,b$ and one of the elements of $\alpha^{(c)}$. Then $r$ that obeys the following equations is normal to the manifold of modified Hopf points at this solution:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Proof.

Equations (20) comprise $3n+2$ equations that depend on $3n+n_{\alpha}+1$ variables and parameters $x^{(c)}$, $a$, $b$, $\omega$ and $\alpha^{(c)}$. In the neighborhood of the given regular solution, these equations define an $(n_{\alpha}-1)$-dimensional manifold of modified Hopf points. By the same arguments as in the proof to Prop. 1, the columns of the matrix $B$ that results from the outer product

[TABLE]

span the normal space to this manifold. Evaluating (22) results in the matrix shown in (21b), which can be seen as follows. Let $B$ be the block matrix with blocks $B_{jk}$, $j=1,\dots,4$, $k=1,\dots,5$ as in (21b) and blocks $B_{jk}$, $j=5$, $k=1,\dots,5$ as in (21c). Block matrices $B_{14}$, $B_{15}$, $B_{21}$, $B_{31}$, $B_{41}$, $B_{44}$ and $B_{45}$ are zero, because the respective functions in (22) do not depend on the variables with respect to which the derivatives are calculated. The expressions for the block matrices $B_{11}$, $B_{24}$, $B_{34}$, $B_{25}$, $B_{35}$ are also obvious. $B_{22}$ results from the outer product

[TABLE]

and $B_{23}$, $B_{32}$, $B_{33}$, $B_{42}$ and $B_{43}$ can be obtained from similar calculations. The remaining blocks, which contain the derivatives of (20b) and of (20c) with respect to $x^{(c)}$ and $\alpha^{(c)}$, require application of the chain and product rule, because the delays $\tau_{k}$ and the Jacobians $A_{k}$ are functions of $x^{(c)}$ and $\alpha^{(c)}$. This yields

[TABLE]

Substituting

[TABLE]

results in the expression for $B_{12}$ stated in the proposition. The derivatives of (20b) with respect to $\alpha^{(c)}$, and of (20c) with respect to $x^{(c)}$ and $\alpha^{(c)}$, can be calculated accordingly and result in the expressions stated in the proposition for $B_{13}$ and, respectively, $B_{52}$ and $B_{53}$. Now recall the columns of $B$ span the normal space to the critical manifold in the space with the components $(x,a,b,\omega,\alpha)$. For the same reasons as stated in the proof of Prop. 1, we seek the particular linear combination $\kappa$ of the columns of $B$ that only has nonzero contributions in the subspace of the parameters $\alpha$. This corresponds to the linear combination $\kappa$ that solves (21b) and (21c). If $\kappa$ solves these equations $\rho\kappa$ solves the equations for any $\rho\in\mathbb{R}$, $\rho\neq 0$. Because unit length of $r$ is enforced with (21d), $\kappa$ and $r$ are uniquely defined. ∎

The augmented system (23) and Prop. 2 treat the case of modified Hopf points with leading eigenvalues $\lambda=\sigma\pm\mathrm{i}\omega$, $\sigma\leq 0$. The case of Hopf bifurcations, i.e., $\sigma=0$, is obviously included. The augmented system for Hopf bifurcations of the treated class of delay differential equations reads (see, e.g., [6, 5])

[TABLE]

where $a$, $b$ and $\omega$ are as in (20). We state the normal vector system for completeness in the following corollary, which immediately follows from Prop. 2 with $\sigma=0$.

Corollary 2** (normal vector to manifold of Hopf bifurcations).**

Assume $(x^{(c)}$, $\alpha^{(c)}$, $\omega$, $a$, $b)$ is a regular solution to (23) in the variables $x^{(c)}$, $\omega$, $a$, $b$ and one of the elements of $\alpha^{(c)}$. Then $r$ that obeys (23), (21b), (21c) and (21d) is normal to the Hopf bifurcation manifold at this solution, where the block matrices $B_{jk}$ from Prop. 2 have to be replaced by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

4. Laser diode optimization with normal vector constraints

With the optimization goal stated in Section 2.1 and the constraints for robust stability described in the previous section, we are able to state the full optimization problem with constraints for robustness. It reads

[TABLE]

where $x=[\mathrm{Re}\{A\},\mathrm{Im}\{A\},N_{1},\dots,N_{K},\Omega]^{\prime}$, $\alpha=[j,\eta,\alpha_{LW},\phi,\tau,T,d]^{\prime}$, and the superscripts $(c)$ and $(0)$ refer to the critical and nominal points, respectively, Equation (24b) enforces steady state operation. Constraints (24c)–(24e) are the robust counterparts to (6c)–(6e). Constraints (24f) and (24g) are the normal vector constraints. Since a robust distance to a manifold of generalized Hopf points has to be enforced in this particular example, they are of the type (20) and (21b)–(21d), respectively. Just as in the motivating example in Section 2, there exist $n_{\alpha}=7$ uncertain parameters (see Table 1). Four of these uncertain parameters are assumed to be uncertain but fixed, whereas the remaining three, $j$, $\eta$ and $\tau$, can be altered to maximize the laser intensity.

While the optimization (6) resulted in an unstable laser operation, for which small disturbances triggered quasi-chaotic behavior, robust stability is guaranteed in (24). The robustness constraints (24f)-(24j) prevent crossing the Hopf bifurcation manifolds shown in Figure 3. The parameter space shown in this figure is spanned by the optimization variables. It is evident from the figure that Hopf bifurcations always appear as the feedback strength $\eta$ is increased. The 3-dimensional plot in Figure 3 also underlines the necessity to interpret 2-dimensional visualizations with care, especially as the the space of uncertain variables is 7-dimensional in this example.

The normal vector constraints (24f)-(24j) take only one critical manifold into account, while multiple such manifolds obviously exist (see Figure 3). Multiple manifolds can be treated rigorously with an approach that resembles active set optimization methods [18]. We do not discuss this aspect here in detail, however, since only one Hopf manifold plays a role.

Solving (24) results in the optimal parameters $j^{(0)}=1$, $\eta^{(0)}=1.911\cdot 10^{-3}$ and $\tau^{(0)}=246.46$. The laser diode operates at an intensity of $|A|^{2}=6.268$ at the optimal robust steady state. We checked the stability of the steady state solutions at the $2^{n_{\alpha}}=128$ vertices of the robustness hypercube (10) to corroborate the optimization result. If compared to the previous unstable result, the feedback strength $\eta$ is smaller. Its upper bound was active in the previous optimization without robustness constraints. In the current robustly stable optimization, it is bounded by the stability requirement. In both cases, the optimal pump current $j$ and the feedback delay $\tau$ are bounded by their upper and lower bounds, respectively. However, the bounds on $j$ and $\tau$ have to be tightened in the current optimization to guarantee that the bounds are met even in the presence of uncertainty (cf. Figure 4b-5b). These parameter differences result in different intensities $|A|^{2}$. The restrictions leading to a robustly stable laser operation result in an intensity reduction of about $6\%$. This reduction can be interpreted to be the price for stability and robustness.

Figures 4–6 visualize the robustly stable optimal point and its robustness regions in the parameter space. All figures contain cuts of the manifolds shown in Figure 3.

Figure 4 illustrates the restricting nature of the stability constraints in the $\eta$-$\tau$-plane. The value of $\eta$ cannot be increased further without violating the robustness constraint. Parts of the hyperspherical outer approximation of the uncertainty region appear to lie in the unstable region. We stress, however, the hypersphere remains completely in the stable region. The impression conveyed by Figure 4 is merely caused by the 2-dimensional cut of the 7-dimensional parameter space.

Figure 5 shows the optimal point in the $j$-$\tau$-plane. It is evident from Figure 5b that, in addition to the lower bound from (24e) on $\tau$, the upper bound from (24c) on $j$ is active at the optimal point.

Figure 6 illustrates the optimal point in the $\phi$-$\alpha_{LW}$-plane. This plane is chosen, because the components of the normal vector along the $\phi$- and $\alpha_{LW}$-axes are the largest ones. The absence of contour lines in this figure accounts for the fact, that $\alpha_{LW}$ and $\phi$ are assumed to be uncertain, but cannot be modified for optimization.

The dynamical behavior of the system at the optimal point is illustrated with Figure 7 for $K=31$. The laser diode converges to its steady state, despite the initial disturbance.

5. Conclusion

We introduced a method for the robust optimization of parametrically uncertain finite-dimensional nonlinear dynamical systems with state- and parameter-dependent delays. Stability and robustness are guaranteed with constraints that enforce a lower bound on the distance of the optimal point to submanifolds of saddle-node and Hopf bifurcations on the steady state manifold. The lower bound on the distance can be interpreted as a finite variation of the model parameters under which stability must be guaranteed; it therefore is a useful robustness measure in applications. We illustrated the proposed method with the optimization of a laser diode system with a delay due to an external cavity.

Acknowledgements

Funding by Deutsche Forschungsgemeinschaft (grant number MO 1086/13-1) is gratefully acknowledged.