**5. Numerical Example**

In a slight refashioning of the notation in the Section 2.2, Equation (12), let the parameter vector *θ* be defined by all the *unknown* parameters defining the interaction functions. Assuming prior distribution *φ*(*θ*) over these unknowns and parameter bounds Θ, we construct the following optimal control problem for robustness against the unknown parameters.

**Problem SD** (Swarm Defense)**.** For *K* defenders and *N* attackers, determine the defender controls *uk*(*t*) that minimize:

$$J = \int\_{\theta} \left[ 1 - P\_0(t\_{f'}, \theta) \right] \phi(\theta) d\theta \tag{27}$$

subject to:

$$\begin{cases} \dot{y}\_k(t) = f(y\_k(t), u\_k(t)), & y\_k(0) = y\_{k0} \\ \ddot{x}\_j(t, \theta) = \\ F\_S(t, \theta) + F\_{HVU}(t, \theta) + F\_D(t, \theta), & x\_j(0, \theta) = x\_{j0}(\theta) \\ Q\_j(t, \theta) = \\ -Q\_j(t, \theta) \sum\_{k=1}^K P\_k(t, \theta) d\_y^{j,k}(x\_j(t, \theta), y\_k(t)), & Q\_j(0, \theta) = 1 \\ P\_k(t) = \\ -P\_k(t, \theta) \sum\_{j=1}^N Q\_j(t, \theta) d\_x^{k,j}(y\_k(t), x\_j(t, \theta)), & P\_k(0, \theta) = 1 \end{cases}$$

for swarm attackers *j* = 1, . . . , *N* and controlled defenders *k* = 1..., *K*.

We implement Problem **SD** for both swarm models in Section 2.1, for a swarm of *N* = 100 attackers and *K* = 10 defenders.

#### *5.1. Example Model 1: Virtual Body Artificial Potential*

The cooperative swarm forces *FS* are defined with the Virtual Body Artificial Potential of Section 2.1 with parameters *α*, *d*<sup>0</sup> and *d*1. In lieu of a potential for the virtual leaders, we assign the HVU tracking function:

$$f\_{HVII} = -\frac{K\_1(x\_i - y\_0)}{||x\_i - y\_0||}\tag{28}$$

where *<sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>R</sup><sup>3</sup> is the position of the HVU. The dissipative force *fvi* <sup>=</sup> <sup>−</sup>*K*2*x*˙*<sup>i</sup>* is employed to guarantee stability of the swarm system. *K*<sup>1</sup> and *K*<sup>2</sup> are positive constants. The swarm's collision avoidance response to the defenders is defined by Equation (4) with parameters *αh*, *h*<sup>0</sup> and *h*1. Since there is only a repulsive force between swarm members and defenders, not an attractive force, we set *h*<sup>1</sup> = *h*0. For attrition, we use the the damage function defined in Equation (21) of [20]:

$$F\_D = \lambda \Phi \left(\frac{F - ar^2}{\sigma}\right), \qquad r = ||x\_i - y\_j||\_2 \tag{29}$$

where Φ is the cumulative normal distribution and ·<sup>2</sup> is the vector 2-norm. This function smoothly penalizes proximity, with the impact decreasing with distance. The parameters *λ*, *F*, *a*, and *σ* shape the steepness of this function and the decline of damage over distance. For the damage rate of defenders inflicted on attackers, we calibrate by the parameters *λD*, *σD*. For the damage rate of attackers inflicted on defenders, we calibrate by the parameters *λA*, *σA*. In both cases, the parameters *F* and *a* in [20] are set to *F* = 0, *a* = 1. Table 1 provides the parameter values that remain fixed in each simulation, and and Table 2 provides the parameters we consider as uncertain.




**Table 2.** Model 1 Varied Parameter Values.

We first use the nominal parameter values provided in Tables 1 and 2 to find a nominal solution defender trajectories that result in the minimum probability of HVU destruction. With the results of these simulations as a reference point, we consider as uncertain each of the parameters that define attacker swarm model and weapon capabilities. In this simulation, these parameters are considered individually. The number of discretization nodes for parameter space was chosen by examination of the Hamiltonian. To illustrate this method and the results obtained in Section 4 we compute Hamiltonians for the Problem **SD** and Model 1 with *θ* = *d*0, *d*<sup>0</sup> ∈ [0.5, 1.5] and *M* = [5, 8, 11]. As *M* increases the sequence of Hamiltonians should converge to the optimal Hamiltonian for the Problem **SD**. For Problem **SD** that should result in a constant, zero-valued Hamiltonian. Figure 3 shows the respective Hamiltonians for *M* = [5, 8, 11]. The value *M* = 11 was chosen for simulations, based on the approximately zero-valued Hamiltonian it generates.

**Figure 3.** Convergence of Hamiltonion as number of parameter nodes *M* increases.

We compare the performance of the solution generated using uncertain parameter optimal control Problem **SD** versus a solution obtained with the nominal values. Figure 4 shows the nominal solution trajectories. The comparitive results of the nominal solutions vs the uncertain parameter control solutions are shown in Figure 5, where the performance of each is shown for different parameters values.

**Figure 4.** Shown are four snapshots during a simulations at *t* = 0, 15, 30, and 45 (time units are arbitrary). Defenders are represented by blue spheres and attackers by red spheres. The HVU is the yellow sphere. Below these snapshots, we show full trajectories for the entire simulation, which is the result of an optimization protocol using the parameters shown in Table 1.

**Figure 5.** Performance of Solutions of Swarm Model 1 as parameter values are varied. Each panel illustrates a different varied parameter, stated on the *x*-axis.

As seen in Figure 5 the trajectories generated by optimization using the nominal values perform poorly over a range of *α*, *d*0, *σA*, *α<sup>k</sup>* and *h*0. In the case of *h*0, for example, this is because the attackers are less repelled by the defenders when *h*<sup>0</sup> is decreased, and they are more able to destroy the HVU from a longer distance as *σ<sup>A</sup>* is increased. The parameter uncertainty solution, however, demonstrates that using the uncertain parameter optimal control framework a solution can be provided which is robust over a range of parameter values. We contrast these results with the case of uncertain parameters *d*<sup>1</sup> and *λA*, also shown in Figure 5. It can be seen that robustness improvements are modest to non-existent for these parameters. This suggests an insensitivity of the problem *d*<sup>1</sup> and *λ<sup>A</sup>* parameters. This kind of analysis can be used to guide inference and observability priorities.

#### *5.2. Example Model 2: Reynolds Boid Model*

To demonstrate flexibility of the proposed framework to include diverse swarm models we have applied the same analysis as was done in Section 5.1 to the Reynolds Boid Model introduced in Section 2.1. We apply the same HVU tracking function as Equation (28). The herding force *FD* of the defenders repelling attackers is applied as a separation force in the form of Equation (10). The fixed parameter values are the same as those in Table 1; the uncertain parameters and ranges are given in Table 3. The results are shown in Figure 6. Again, we see that the tools developed in this paper can be used to gain an insight into the robustness properties of the nominal versus uncertain parameter solutions. For example, we can see that the uncertain parameter solutions perform much better than the nominal ones for the cases where *λ*, *σ* and *wI* are uncertain.

**Figure 6.** Performance of Solutions of Swarm Model 2 as parameter values are varied.



#### **6. Conclusions**

In this paper, we have built on our previous work on developing an efficient numerical framework for solving uncertain parameter optimal control problems. Unlike uncertainties introduced into systems due to stochastic "noise", parameter uncertainties do not average or cancel out in regard to their effects. Instead, each possible parameter value creates a specific profile of possibility and risk. The uncertain optimal control framework which has been developed for these problems exploits this inherent structure by producing answers which have been optimized over all parameter profiles. This approach takes into account the possible performance ranges due to uncertainty, while also utilizing what information is known about the uncertain features, such as parameter domains and prior probability distributions over the parameters. Thus, we are able to contain risk, while providing plans which have been optimized for performance under all known conditions. The results reported in this paper include analysis of the consistency of the adjoint variables of the numerical solution. In addition, the paper includes a numerical analysis of a large scale adversarial swarm engagement that clearly demonstrates the benefits of using the proposed framework.

There are many directions for future work for the topics of this paper. The numerical simulations in this paper consider the parameters individually, as one-dimensional parameter spaces. However, Problem **P** allows for multi-dimensional parameter spaces. A more dedicated implementation, taking advantage of the parallelizable form of Equation (16), for example, could certainly manage several simultaneous parameters. Exponential growth as parameter space dimension increases is an issue for both the quadrature format of Equation (15) and handling of the state space size for Equation (16). This can be somewhat mitigated by using sparse grid methods for high-dimensional integration to define the nodes in Equation (15). For large enough sizes, Monte Carlo sampling, rather than quadrature might be more appropriate for designating parameter nodes.

A further direction for future work would be to incorporate these methods into the design of more responsive closed-loop control solutions. The optimization methods in this paper provide open-loop controls. While useful, closed-loop controls would be more ideal for dynamic situations with uncertainty. There are many ways, however, that open-loop solutions can provide stepping stones to developing closed-loop solutions. For instance, Ref. [26] utilizes closed-loop solutions to train a neural network to learn an optimal closed-loop control strategy. Open-loop solutions can also be used to provide initial guesses to discretized closed-loop optimizations, seeding the optimization algorithm.

Another direction for future work is in the greater application of the duality results of Section 4. The numerical results in this paper simply utilize the Hamiltonian consistency. The proof of Theorem 1, however, additionally demonstrates the consistency of the adjoint variables for the problem. As the results demonstrate, parameter sensitivity for these swarm models is highly non-linear. The numerical solutions of Section 5 are able to demonstrate this sensitivity by applying the solution to varied parameter values. However, this is actually a fairly expensive method for a large swarm, as it involves re-evaluation of the

swarm ODE for each parameter value. More importantly, it would not be scalable to highdimensional parameter spaces, as the exponential growth of that approach to sensitivity analysis would be unavoidable. The development of an analytical adjoint sensitivity method for this problem could be of great utility for paring down numerical simulations to only focus on the parameters most relevant to success.

**Author Contributions:** Formal analysis, T.T.; Writing—original draft, C.W. and I.K.; Writing review & editing, Q.G. and A.H.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Office of Naval Research under Grant No. N0001421WX01974. Additionally, this work was supported by the Department of the Navy, Office of Naval Research, Consortium for Robotics Unmanned Systems Education and Research at the Naval Postgraduate School.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We would like to thank the University of Texas at San Antonio, the University of California Santa Cruz, and the Naval Postgraduate School for their administrative and research support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Assumptions and Definitions**

We we impose the assumptions in Section 2 of [4]. The definition of accumulation point used in the following proof can be found in Definition 3.2 of [4]. The following assumption is placed on the choice of numerical integration scheme to be utilized in approximating Problem **P**:

**Assumption A1.** *For each <sup>M</sup>* <sup>∈</sup> <sup>N</sup>*, there is a set of nodes* {*θ<sup>M</sup> <sup>i</sup>* }*<sup>M</sup> <sup>i</sup>*=<sup>1</sup> ⊂ Θ *and an associated set of weights* {*α<sup>M</sup> <sup>i</sup>* }*<sup>M</sup> <sup>i</sup>*=<sup>1</sup> <sup>⊂</sup> <sup>R</sup>*, such that for any continuous function h* : <sup>Θ</sup> <sup>→</sup> <sup>R</sup>*,*

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

This is the same as Assumption 3.1 of [4]; we include it for reference. In additions to the assumptions of [4], we also impose the following:

**Assumption A2.** *The functions <sup>f</sup> and <sup>r</sup> are <sup>C</sup>*1*. The set* <sup>Θ</sup> *is compact and <sup>x</sup>*<sup>0</sup> : <sup>Θ</sup> #→ <sup>R</sup>*nx is continuous. Moreover, for the compact sets X and U defined in Assumptions 2.3 and 2.4 of [4], and for each t* ∈ [0, *T*]*, θ* ∈ Θ*, the Jacobians rx and fx are Lipschitz on the set X* × *U, and the corresponding Lipschitz constants Lr and Lf are uniformly bounded in θ and t. The function F is <sup>C</sup>*<sup>1</sup> *on X for all <sup>θ</sup>* <sup>∈</sup> <sup>Θ</sup>*; in addition, F and Fx are continuous with respect to <sup>θ</sup>.*

### **Appendix B. Main Theorem Proof**

The theorem relies on the following lemma:

**Lemma A1.** *Let* {*uM*} *be a sequence of optimal controls for Problem* **PM** *with an accumulation point <sup>u</sup>*<sup>∞</sup> *for the infinite set <sup>V</sup>* <sup>⊂</sup> <sup>N</sup>*. Let* (*x*∞(*t*, *<sup>θ</sup>*), *<sup>λ</sup>*∞(*t*, *<sup>θ</sup>*)) *be the solution to the dynamical system:*

$$\begin{cases} \dot{\mathbf{x}}^{\infty}(t,\theta) = f(x^{\infty}(t,\theta), u^{\infty}(t), \theta) \\ \dot{\lambda}^{\infty}(t,\theta) = -\frac{\partial H(x^{\infty}(t,\theta), \lambda^{\infty}(t,\theta), u^{\infty}(t), t, \theta)}{\partial x} \end{cases} \tag{A1}$$

$$\begin{cases} \mathbf{x}^{\infty}(0,\theta) = \mathbf{x}\_0(\theta) \\ \lambda^{\infty}(T,\theta) = \frac{\partial F(\mathbf{x}^{\infty}(T,\theta),\theta)}{\partial \mathbf{x}} \end{cases} \tag{A2}$$

*where H is defined as per Equation (19), and let* {(*xM*(*t*, *θ*), *λM*(*t*, *θ*))} *for M* ∈ *V be the sequence of solutions to the dynamical systems:*

$$\begin{cases}
\dot{\mathbf{x}}\_M(t,\theta) = f(\mathbf{x}\_M(t,\theta), \boldsymbol{\mu}\_M(t), \theta) \\
\dot{\lambda}\_M(t,\theta) = -\frac{\partial H(\mathbf{x}\_M(t,\theta), \boldsymbol{\lambda}\_M(t,\theta), \boldsymbol{\mu}\_M(t), t, \theta)}{\partial \boldsymbol{x}}
\end{cases} \tag{A.3}$$

$$\begin{cases} \mathbf{x}\_M(0, \theta) = \mathbf{x}\_0(\theta) \\ \lambda\_M(T, \theta) = \frac{\partial F(\mathbf{x}\_M(T, \theta), \theta)}{\partial \mathbf{x}} \end{cases} \tag{A4}$$

*Then, the sequence* {(*xM*(*t*, *<sup>θ</sup>*), *<sup>λ</sup>M*(*t*, *<sup>θ</sup>*))} *converges pointwise to* (*x*∞(*t*, *<sup>θ</sup>*), *<sup>λ</sup>*∞(*t*, *<sup>θ</sup>*)) *and this convergence is uniform in θ.*

**Proof.** The convergence of {*xM*(*t*, *θ*)} is given by Lemmas 3.4 and 3.5 of [4]. The convergence of the sequence of solutions {*λM*(*t*, *θ*)} is guaranteed by the optimality of {*uM*}. The convergence of {*λM*(*t*, *θ*)} then follows the same arguments given the convergence of {*xM*(*t*, *θ*)}, utilizing the regularity assumptions placed on the derivatives of *F*, *r*, and *f* with respect to *x* to enable the use of Lipschitz conditions on the costate dynamics and transversality conditions.

**Remark A1.** *Note that λM*(*t*, *θ*) *is not a costate of Problem* **P***λM, since it is a function of θ. However, when θ* = *θ<sup>M</sup> <sup>i</sup> , then <sup>λ</sup>M*(*t*, *<sup>θ</sup><sup>M</sup> <sup>i</sup>* ) = *<sup>λ</sup>*˜ *<sup>M</sup> <sup>i</sup>* (*t*)*, where <sup>λ</sup>*˜ *<sup>M</sup> <sup>i</sup> is the costate of Problem* **<sup>P</sup>***λ<sup>M</sup> generated by the pair of solutions to Problem* **PM***,* (*x*˜*<sup>M</sup> <sup>i</sup>* , *u*<sup>∗</sup> *<sup>M</sup>*) *. In other words, the function λM*(*t*, *θ*) *matches the costate values at all collocation nodes. Since these values satisfy the dynamics equations of Problem* **P***λM, a further implication of this is that the values of λM*(*t*, *θ<sup>M</sup> <sup>i</sup>* ) *produce feasible solutions to Problem* **P***λM.*

**Remark A2.** *Since the functions* {(*xM*(*t*, *<sup>θ</sup>*), *<sup>λ</sup>M*(*t*, *<sup>θ</sup>*))} *obey the respective identities xM*(*t*, *<sup>θ</sup><sup>M</sup> <sup>i</sup>* ) = *x*˜*<sup>M</sup> <sup>i</sup>* (*t*) *and <sup>λ</sup>M*(*t*, *<sup>θ</sup><sup>M</sup> <sup>i</sup>* ) = *<sup>λ</sup>*˜ *<sup>M</sup> <sup>i</sup>* (*t*)*, their convergence to* (*x*∞(*t*, *<sup>θ</sup>*), *<sup>λ</sup>*∞(*t*, *<sup>θ</sup>*)) *also implies the convergence of the sequence of discretized primals and duals,* {*X*˜ *<sup>M</sup>*} *and* {Λ˜ *<sup>M</sup>*}*, to accumulation points given by the relations*

$$\lim\_{M \in \tilde{V}} \tilde{\boldsymbol{x}}\_l^M(t) = \boldsymbol{x}^{\infty}(t, \boldsymbol{\theta}\_l^M), \quad \lim\_{M \in \tilde{V}} \tilde{\boldsymbol{\lambda}}\_l^M(t) = \boldsymbol{\lambda}^{\infty}(t, \boldsymbol{\theta}\_l^M).$$

We now prove Theorem 1. Let {(*xM*(*t*, *θ*), *λM*(*t*, *θ*))} for *M* ∈ *V* be the sequence of solutions defined by Equation (A3) and let (*x*∞(*t*, *θ*), *λ*∞(*t*, *θ*)) be the accumulation functions defined by Equation (A1). Incorporating Remarks A1 and A2, we have:

$$\begin{aligned} &\lim\_{M\in V} H^M(\mathcal{X}\_{M\_I}\bar{\Lambda}\_{M\_I}u\_{M\_I}, t) = \\ &\lim\_{M\in V} \sum\_{i=1}^M a\_i^M \left[ \bar{\lambda}\_i^M(t) f(\bar{\mathfrak{x}}\_i^M(t), u(t), \theta\_i^M) + r(\bar{\mathfrak{x}}\_i^M(t), u(t), t, \theta\_i^M) \right] \\ &= \lim\_{M\in V} \sum\_{i=1}^M a\_i^M \left[ \lambda\_M(t, \theta\_i^M) f(\ge\_M (t, \theta\_i^M), u(t), \theta\_i^M) + r(\ge\_M (t, \theta\_i^M), u(t), t, \theta\_i^M) \right] \end{aligned}$$

Due to the results of Lemma A1, and applying Remark 1 of [4] on the convergence of the quadrature scheme for uniformly convergent sequences of continuous functions, we find that:

$$\lim\_{M \in V} H^M(\mathcal{R}\_{M\prime} \bar{\mathcal{A}}\_{M\prime} u\_{M\prime} t) = 0$$

$$\int\_{\Theta} [\lambda^{\infty}(t,\theta)f(\mathbf{x}^{\infty}(t,\theta),\mu^{\infty}(t),\theta) + r(\mathbf{x}^{\infty}(t,\theta),\mu^{\infty}(t),t,\theta)]d\theta = \mathbf{H}(\mathbf{x}^{\infty},\lambda^{\infty},\mu^{\infty},t)$$

thus proving the theorem.
