**3. Problem Formulation**

The above models depend on a large number of parameters. The dynamic swarming model coupled with attrition functions results in over a dozen key parameters, and many more would result from a non-homogeneous swarm. A concern would be that this adds too much model specificity, making optimal defense strategies lack robustness due to sensitivity to the specific set of model parameters. This concern turns out to be justified. When defense strategies are optimized for fixed, nominal parameter values, they display catastrophic failure for small perturbations of certain parameters, as can be seen in Figure 1. In fact, the plots included in Figure 1 clearly demonstrate that the sensitivity of the cost with respect to the uncertain parameters is highly non-linear. Thus, generating robust defense strategies requires a more sophisticated formalism introduced in the next Section 3.1.

**Figure 1.** Example performance of solutions calculated using nominal values when parameter value is varied. Calculated using values in Section 5.1. Magenta dot marks the nominal value used in the optimization problem. In the left panel, *d*<sup>0</sup> is varied as the parameter; in the right panel *h*<sup>0</sup> is varied.

#### *3.1. Uncertain Parameter Optimal Control*

The class of problems addressed by the computational algorithm is defined as follows:

**Problem P.** Determine the function pair (*x*, *<sup>u</sup>*) with *<sup>x</sup>* <sup>∈</sup> *<sup>W</sup>*1,∞([0, *<sup>T</sup>*] <sup>×</sup> <sup>Θ</sup>; <sup>R</sup>*nx* ), *<sup>u</sup>* <sup>∈</sup> *L*∞([0, *T*]; R*nu* ) that minimizes the cost

$$J[\mathbf{x}, \boldsymbol{u}] = \int\_{\Theta} \left[ F(\mathbf{x}(T, \theta), \theta) + \int\_{0}^{T} r(\mathbf{x}(t, \theta), \boldsymbol{u}(t), t, \theta) dt \right] d\theta \tag{13}$$

subject to the dynamics

$$\frac{\partial \mathbf{x}}{\partial t}(t, \theta) = f(\mathbf{x}(t, \theta), \mathbf{u}(t), \theta), \tag{14}$$

initial condition *x*(0, *θ*) = *x*0(*θ*), and the control constraint *g*(*u*(*t*)) ≤ 0 for all *t* ∈ [0, *T*]. The set *<sup>L</sup>*∞([0, *<sup>T</sup>*]; <sup>R</sup>*nu* ) is the set of all essentially bounded functions, *<sup>W</sup>*1,∞([0, *<sup>T</sup>*] <sup>×</sup> <sup>Θ</sup>; <sup>R</sup>*nx* ) the Sobolev space of all essentially bounded functions with essentially bounded distributional derivatives, and *<sup>F</sup>* : <sup>R</sup>*nx* <sup>×</sup> <sup>R</sup>*n<sup>θ</sup>* #→ <sup>R</sup>, *<sup>r</sup>* : <sup>R</sup>*nx* <sup>×</sup> <sup>R</sup>*nu* <sup>×</sup> <sup>R</sup> <sup>×</sup> <sup>R</sup>*n<sup>θ</sup>* #→ <sup>R</sup>, *<sup>g</sup>* : <sup>R</sup>*nu* #→ <sup>R</sup>*ng* . Additional conditions imposed on the state and control space and component functions are specified in Appendix A.

In Problem **<sup>P</sup>**, the set <sup>Θ</sup> is the domain of a parameter *<sup>θ</sup>* <sup>∈</sup> <sup>R</sup>*n<sup>θ</sup>* . The format of the cost functional is that of the integral over Θ of a Mayer–Bolza type cost with parameter *θ*. This parameter can represent a range of values for a feature of the system, such as in ensemble control [21], or a stochastic parameter with a known probability density function.

For computation of numerical solutions, we introduce an approximation of Problem **P**, referred to as Problem **PM**. Problem **PM** is created by approximating the parameter space, Θ, with a numerical integration scheme. This numerical integration scheme is defined in terms of a finite set of *<sup>M</sup>* nodes {*θ<sup>M</sup> <sup>i</sup>* }*<sup>M</sup> <sup>i</sup>*=<sup>1</sup> and an associated set of *<sup>M</sup>* weights {*α<sup>M</sup> <sup>i</sup>* }*<sup>M</sup> <sup>i</sup>*=<sup>1</sup> <sup>⊂</sup> <sup>R</sup> such that

$$\int\_{\Theta} h(\theta)d\theta = \lim\_{M \to \infty} \sum\_{i=1}^{M} h(\theta\_i^M) \alpha\_i^M. \tag{15}$$

given certain function smoothness assumptions. See Appendix A Assumption A1 for formal assumptions. Throughout the paper, *M* is used to denote the number of nodes used in this approximation of parameter space.

For a given set of nodes {*θ<sup>M</sup> <sup>i</sup>* }*<sup>M</sup> <sup>i</sup>*=1, and control *<sup>u</sup>*(*t*), let *<sup>x</sup>*¯*<sup>M</sup> <sup>i</sup>* (*t*), *i* = 1, ... , *M*, be defined as the solution to the ODE created by the state dynamics of Problem **P** evaluated at *θ<sup>M</sup> i* :

$$\begin{cases} \frac{d\mathfrak{x}\_i^M}{dt}(t) = f(\mathfrak{x}\_i^M(t), u(t), \theta\_i^M) \\ \mathfrak{x}\_i^M(0) = \mathfrak{x}\_0(\theta\_i^M), \end{cases} \quad i = 1, \dots, M. \tag{16}$$

Let *X*¯ *<sup>M</sup>*(*t*)=[*x*¯*<sup>M</sup>* <sup>1</sup> (*t*), ... , *<sup>x</sup>*¯*<sup>M</sup> <sup>M</sup>*(*t*)]. The system of ODEs defining *<sup>X</sup>*¯ *<sup>M</sup>* has dimension *nx* × *M*, where *nx* is the dimension of the original state space and *M* is the number of nodes. The numerical integration scheme for parameter space creates an approximate objective functional, defined by:

$$\int \mathbf{J}^{M}[\bar{X}^{M}, \mathbf{u}] = \sum\_{i=1}^{M} \left[ F\left(\mathfrak{X}\_{i}^{M}(T), \theta\_{i}^{M}\right) + \int\_{0}^{T} r\left(\mathfrak{X}\_{i}^{M}(t), \mathbf{u}(t), t, \theta\_{i}^{M}\right) dt \right] \mathbf{a}\_{i}^{M}. \tag{17}$$

In [4], the consistency of **PM** is proved. This is the property that, if optimal solutions to Problem **PM** converge as the number of nodes *<sup>M</sup>* <sup>→</sup> <sup>∞</sup>, they converge to feasible, optimal solutions of Problem **P**. See [4] for detailed proof and assumptions.

#### *3.2. Computational Efficiency*

The computation time of the numerical solution to the discretized problem defined in Equations (16) and (17) will depend on the value of *M*. Ideally, it should be sufficiently small so as to allow for a fast solution. On the other hand, a value of *M* that is too small will result in a solution that is not particularly useful, i.e., too far from the optimal. Naturally, the question arises: how far is a particular solution from the optimal? One tool for assessing this lies in computing the Hamiltonian and is addressed in Section 4.

#### **4. Consistency of Dual Variables**

The dual variables provide a method to determine the solution of an optimal control problem or a tool to validate a numerically computed solution. For numerical schemes based on direct discretization of the control problem, analyzing the properties of the dual variables and their resultant Hamiltonian may also lead to insight into the validity of an approximation scheme [22,23]. This could be especially helpful in highdimensional problems, such as swarming, where parsimonious discretization is crucial to computational tractability.

Previous work shows the consistency of the primal variables in approximate Problem **PM** to the original parameter uncertainty framework of Problem **P**. Here, we build on that and prove the consistency of the dual problem of Problem **P** as well. This theoretical contribution is diagrammed in Figure 2. The consistency of the dual problem in parameter space enables approximate computation of the Hamiltonian from numerical solutions.

**Figure 2.** Diagram of primal and dual relations for parameter uncertainty control. Red lines designate the contribution of this paper.

In [24], necessary conditions for Problem **P** were established. These conditions are as follows:

**Problem P***<sup>λ</sup>* [([24], pp. 80–82)]**.** If (*x*∗, *u*∗) is an optimal solution to Problem **P**, then there exists an absolutely continuous costate vector *λ*∗(*t*, *θ*), such that for *θ* ∈ Θ:

$$\frac{\partial \lambda^\*}{\partial t}(t,\theta) = -\frac{\partial H(\mathbf{x}^\*, \lambda^\*, \mathbf{u}^\*, t, \theta)}{\partial \mathbf{x}},$$

$$\lambda^\*(T, \theta) = \frac{\partial F(\mathbf{x}^\*(T, \theta), \theta)}{\partial \mathbf{x}}\tag{18}$$

where *H* is defined as:

$$H(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}, \mathbf{t}, \boldsymbol{\theta}) = $$

$$\lambda f(\mathbf{x}(t, \boldsymbol{\theta}), \boldsymbol{\mu}(t), \boldsymbol{\theta}) + r(\mathbf{x}(t, \boldsymbol{\theta}), \boldsymbol{\mu}(t), \mathbf{t}, \boldsymbol{\theta}). \tag{19}$$

Furthermore, the optimal control *u*∗ satisfies

$$\mu^\*(t) = \underset{\mathfrak{u} \in U}{\text{arg min}} \mathbf{H}(\mathfrak{x}^\*, \lambda^\*, \mathfrak{u}, t),$$

where **H** is given by

$$\mathbf{H}(\mathbf{x}, \lambda, \mu, t) = \int\_{\Theta} H(\mathbf{x}, \lambda, \mu, t, \theta) d\theta. \tag{20}$$

Because Problem **PM** is a standard non-linear optimal control problem, it admits a dual problem as well. Problem **P***Mλ*, provided by the Pontryagin minimum principle (a survey of minimum principle conditions is given by [25]). Applied to **PM** this generates:

**Problem P***Mλ***.** For feasible solution (*X*¯ *<sup>M</sup>*, *u*) to Problem **PM**, find Λ¯ (*t*)=[*λ*¯ *<sup>M</sup>* <sup>1</sup> (*t*), ... *<sup>λ</sup>*¯ *<sup>M</sup> <sup>M</sup>*(*t*)], *λ*¯ *<sup>M</sup> <sup>i</sup>* : [0, *<sup>T</sup>*] <sup>→</sup> <sup>R</sup>*Nx* , that satisfies the following conditions:

$$\frac{d\check{\lambda}\_i^M}{dt}(t) = -\frac{\partial H^M(\vec{x}\_i^M, \vec{\lambda}\_i^M, u, t)}{\partial \mathfrak{x}\_i^M},$$

$$\lambda\_i^M(T) = \mathfrak{a}\_i^M \frac{\partial F(\vec{x}\_i^M, \theta\_i^M)}{\partial \vec{x}\_i^M},\tag{21}$$

where *H*¯ *<sup>M</sup>* is defined as:

$$\bar{H}^{M}(\bar{X}^{M}, \bar{\Lambda}\_{M\prime} u, t) = \sum\_{i=1}^{M} \left[ \bar{\lambda}\_{i}^{M} f(\bar{x}\_{i}^{M}(t), u(t), \theta\_{i}^{M}) + a\_{i}^{M} r(\bar{x}\_{i}^{M}(t), u(t), t, \theta\_{i}^{M}) \right]. \tag{22}$$

An alternate direction from which to approach solving Problem **P** overall is to approximate the necessary conditions of Problem **P** , i.e., Problem **P***λ*, directly rather than to approximate Problem **P**. This creates the system of equations:

$$\frac{d\lambda}{dt}(t,\theta\_i^M) = -\frac{\partial H(\mathbf{x}, \mathbf{u}, t, \theta\_i^M)}{\partial \mathbf{x}}$$

$$\lambda(T, \theta\_i^M) = \frac{\partial F(\mathbf{x}(T, \theta\_i^M), \theta\_i^M)}{\partial \mathbf{x}}\tag{23}$$

for *i* = 1, . . . , *M*, where *H* is defined as:

$$H(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{u}, t, \boldsymbol{\theta}) = \lambda f(\mathbf{x}(t, \boldsymbol{\theta}), \boldsymbol{u}(t), \boldsymbol{\theta}) + r(\mathbf{x}(t, \boldsymbol{\theta}), \boldsymbol{u}(t), t, \boldsymbol{\theta}).$$

This system of equations can be re-written in terms of the quadrature approximation of the stationary Hamiltonian defined in Equation (20). Define

$$H^M(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}, t) := \sum\_{i=1}^M \boldsymbol{\alpha}\_i^M H(\mathbf{x}(t, \theta\_i^M), \boldsymbol{\lambda}(t, \theta\_i^M), \boldsymbol{\mu}(t), t, \theta\_i^M).$$

Let

$$\tilde{\Lambda}(t) = [\tilde{\lambda}\_1^M(t), \dots, \tilde{\lambda}\_M^M(t)] = [\lambda(t, \theta\_1^M), \dots, \lambda(t, \theta\_M^M)]^T$$

and let

$$\mathfrak{X}\_{\mathcal{M}} = [\mathfrak{x}\_1^{\mathcal{M}}(t), \dots, \mathfrak{x}\_{\mathcal{M}}^{\mathcal{M}}(t)]$$

denote the semi-discretized states from Equation (16). Equation (23) can then be written as:

$$\frac{d\bar{\lambda}\_i^M}{dt}(t) = -\frac{1}{\mathfrak{a}\_i^M} \cdot \frac{\partial H^M(\bar{\mathcal{X}}\_M, \bar{\mathcal{A}}, \boldsymbol{\mu}, t)}{\partial \mathfrak{x}\_i^M}$$

$$\bar{\lambda}\_i^M(T) = \frac{\partial F(\mathfrak{x}\_i^M(T), \theta\_i^M)}{\partial \mathfrak{x}\_i^M} \tag{24}$$

for *i* = 1, . . . , *M*. Thus, we reach the following discretized dual problem:

**Problem P***λM***.** For feasible controls *u* and solutions *X*˜ *<sup>M</sup>* to Equation (16), find Λ˜ (*t*)=[*λ*˜ *<sup>M</sup>* <sup>1</sup> (*t*),... *<sup>λ</sup>*˜ *<sup>M</sup> <sup>M</sup>*(*t*)], *<sup>λ</sup>*˜ *<sup>M</sup> <sup>i</sup>* : [0, *<sup>T</sup>*] <sup>→</sup> <sup>R</sup>*nx* , that satisfies the following conditions:

$$\frac{d\bar{\lambda}\_i^M}{dt}(t) = -\frac{1}{a\_i^M} \cdot \frac{\partial H^M(\bar{\mathcal{X}}\_M, \bar{\mathcal{A}}, \mu, t)}{\partial \bar{x}\_i^M},$$

$$\bar{\lambda}\_i^M(T) = \frac{\partial F(\mathfrak{x}\_i^M, \theta\_i^M)}{\partial \mathfrak{x}\_i^M},\tag{25}$$

where *H*˜ *<sup>M</sup>* is defined as:

$$\tilde{H}^{M}(\bar{\mathbf{X}}\_{M}, \bar{\mathbf{A}}\_{M}, \boldsymbol{\mu}, \boldsymbol{t}) = \sum\_{i=1}^{M} \left[ a\_{i}^{M} \tilde{\lambda}\_{i}^{M} f(\bar{\mathbf{x}}\_{i}^{M}(t), \boldsymbol{\mu}(t), \boldsymbol{\theta}\_{i}^{M}) + a\_{i}^{M} r(\bar{\mathbf{x}}\_{i}^{M}(t), \boldsymbol{\mu}(t), \boldsymbol{t}, \boldsymbol{\theta}\_{i}^{M}) \right]. \tag{26}$$

**Lemma 1.** *The mapping:*

$$(\mathfrak{x}\_i^{\mathcal{M}}, \vec{u}) \mapsto (\mathfrak{x}\_i^{\mathcal{M}}, \vec{u})\_{\prime} \quad \frac{\bar{\lambda}\_i^{\mathcal{M}}}{\mathfrak{a}\_i^{\mathcal{M}}} \mapsto \bar{\lambda}\_i^{\mathcal{M}}.$$

*for i* = 1, . . . , *M is a bijective mapping from solutions of Problem* **P***M<sup>λ</sup> to Problem* **P***λM.*

**Proof.** In the case of this particular problem, unlike standard control, the collocation of the relevant dynamics involves no approximation of differentiation (since the discretization is in the parameter domain rather than the time domain), and, thus, the mapping of covectors between Problem **P***M<sup>λ</sup>* and Problem **P***H*˜ *<sup>M</sup>*(*X*˜ *<sup>M</sup>*,Λ˜ *<sup>M</sup>*,*u*,*t*)=*λ<sup>M</sup>* is straightforward and simply constructively provided by the lemma itself. The two mappings of the lemma, (*x*¯*<sup>M</sup> <sup>i</sup>* , *u*¯) #→ (*x*˜*<sup>M</sup> <sup>i</sup>* , *<sup>u</sup>*˜) (identity mapping) and *<sup>λ</sup>*¯ *<sup>M</sup> i α<sup>M</sup> i* #→ *<sup>λ</sup>*˜ *<sup>M</sup> <sup>i</sup>* (scaling by <sup>1</sup> *α<sup>M</sup> i* ) are both bijections.

**Theorem 1.** *Let* {*X*˜ *<sup>M</sup>*, <sup>Λ</sup>˜ *<sup>M</sup>*, *uM*}*M*∈*<sup>V</sup> be a sequence of solutions for Problem* **<sup>P</sup>***λ<sup>M</sup> with an accumulation point* {*X*˜ <sup>∞</sup>, <sup>Λ</sup>˜ <sup>∞</sup>, *<sup>u</sup>*∞}*. Let* (*x*∞, *<sup>λ</sup>*∞, *<sup>u</sup>*∞) *be the solutions to Problem* **<sup>P</sup>***<sup>λ</sup> for the control u*∞*. Then*

$$\lim\_{M \in V} \check{H}^M(\check{X}\_{M\prime}\check{\Lambda}\_{M\prime}u\_{M\prime}t) = \mathbf{H}(\boldsymbol{\pi}^{\infty}, \boldsymbol{\lambda}^{\infty}, \boldsymbol{\mu}^{\infty}, t)$$

*where H*˜ *<sup>M</sup> is the Hamiltonian of Problem* **P***λ<sup>M</sup> as defined by Equation (26) and* **H** *is the Hamiltonian of Problem* **P** *as defined by Equation (20). The proof of this theorem can be found in the Appendix B.*

The convergence of the Hamiltonians of the approximate, standard control problems to the Hamiltonian of the general problem, **H**(*x*∞, *λ*∞, *u*∞, *t*), means that many of the useful features of the Hamiltonians of standard optimal control problems are preserved. For instance, it is straightforward to show that the satisfaction of Pontryagin's minimum principle by the approximate Hamiltonians implies minimization of **H**(*x*∞, *λ*∞, *u*∞, *t*) as well. That is, that

$$\mathbf{H}(\mathfrak{x}^{\infty}, \lambda^{\infty}, \mathfrak{u}^{\infty}, t) \le \mathbf{H}(\mathfrak{x}^{\infty}, \lambda^{\infty}, \mathfrak{u}, t)$$

for all feasible *u*. Furthermore, when applicable, the stationarity properties of the standard control Hamiltonian, such as a constant-valued Hamiltonian in time-invariant problems, or stationarity with respect to *u*(*t*) in problems with open control regions, are also preserved.
