Mesh Refinement Method for Solving Bang-Bang Optimal Control Problems Using Direct Collocation

TL;DR

This paper introduces a mesh refinement technique for solving bang-bang optimal control problems with direct collocation, accurately identifying switch times and improving computational efficiency.

Contribution

It presents a novel mesh refinement approach that automatically estimates control discontinuities and switch times, enhancing solution accuracy and efficiency.

Findings

Accurately determines control switch times using roots of switching functions.

Improves computational efficiency over previous methods.

Demonstrates effectiveness on three example problems.

Abstract

A mesh refinement method is developed for solving bang-bang optimal control problems using direct collocation. The method starts by finding a solution on a coarse mesh. Using this initial solution, the method then determines automatically if the Hamiltonian is linear with respect to the control, and, if so, estimates the locations of the discontinuities in the control. The switch times are estimated by determining the roots of the switching functions, where the switching functions are determined using estimates of the state and costate obtained from the collocation method. The accuracy of the switch times is then improved on subsequent meshes by dividing the original optimal control problem into multiple domains and including variables that define the locations of the switch times. While in principle any collocation method can be used, in this research the previously developed Legendre-Gauss-Radau collocation method is employed because it provides an accurate approximation of the costate which in turn improves the approximation of the switching functions. The method of this paper is designed to be used with a previously developed mesh refinement method in order to accurately approximate the solution in segments where the solution is smooth. The method is demonstrated on three examples where it is shown to accurately determine the switching structure of a bang-bang optimal control problem. When compared with previously developed mesh refinement methods, the results demonstrate that the method developed in this paper improves computational efficiency when solving bang-bang optimal control problems.

Tables

$𝐚$	=	vector field for right-hand side of dynamics
$𝐀$	=	matrix defining vector field for right-hand side of dynamics at collocation points
$𝐛$	=	vector field for boundary conditions
$𝐜$	=	vector field for path constraints
$𝐂$	=	matrix defining vector field for path constraints at collocation points
$𝐃$	=	Legendre-Gauss-Radau differentiation matrix
${𝐃}_{N+1}$	=	last column of Legendre-Gauss-Radau differentiation matrix
$f$	=	$𝐯$ independent contribution of $\mathscr{H}$
$\mathscr{H}$	=	Hamiltonian of optimal control problem
$𝒥$	=	objective functional
$K$	=	number of mesh intervals used
${\ell }_j^{\left(k\right)}$	=	Lagrange polynomial $j$ of mesh interval $k$
$ℒ$	=	vector field for integrand appearing in Lagrange cost
$𝐋$	=	matrix defining vector field for integrand at collocation points
$ℳ$	=	Mayer cost
$M$	=	number of mesh iterations
$n_y$	=	number of state components
$n_u$	=	number of control components
$n_c$	=	number of path constraints
$n_s$	=	number of switching time parameters
$N_k$	=	number of collocation points used in mesh interval $k$
$N$	=	total number of collocation points used
$N_f$	=	total number of collocation points used for problem on final mesh
${𝒫}_d$	=	$d^{th}$ time domain
$Q$	=	number of constant control domains used to approximate bang-bang control profile
${𝒮}_k$	=	mesh interval $k$
$t_0$	=	initial time
$t_f$	=	final time
$t_s^{\left[S\right]}$	=	$S^{th}$ switch time parameter
$t_d^i$	=	estimated discontinuity time for $i^{th}$ control component
$t_{\sigma }^i$	=	midpoint time of two adjacent points containing a change in sign of ${\sigma }_i$
$t_u^i$	=	midpoint time of two adjacent points containing largest absolute difference
		in value of $i^{th}$ control component
$𝒯$	=	computation time
$𝐮\left(t\right)$	=	control in time horizon
$𝐔$	=	matrix of control parameterization at collocation points
$𝐯$	=	linearly dependent components of control in Hamiltonian
$𝐰$	=	vector of corresponding LGR weights at collocation points
$𝐖$	=	diagonal matrix of corresponding LGR weights at collocation points
$𝐲\left(t\right)$	=	state in time horizon
$𝐘\left(\tau \right)$	=	state approximation in $\tau$ domain
$𝐘$	=	matrix of state approximation at discretized points
$𝐳$	=	nonlinearly dependent components of control in Hamiltonian
$𝚫$	=	defect constraints matrix
$𝝀\left(t\right)$	=	costate in time horizon
$𝝀$	=	matrix of costate estimates at collocation points
${𝝀}_{N+1}$	=	vector of costate estimates at non-collocated end point
$𝚲$	=	matrix of NLP multipliers corresponding to defect constraints at collocation points
$𝝈$	=	vector field for switching functions of $\mathscr{H}$
$\tau$	=	domain used for Legendre-Gauss-Radau collocation
${\tau }_j^{\left(k\right)}$	=	support point $j$ of mesh interval $k$

Table 2. Table 2: Mesh refinement performance results for Example 1 using h ​ p ℎ 𝑝 hp -BB, h ​ p ℎ 𝑝 hp -I, h ​ p ℎ 𝑝 hp -II, h ​ p ℎ 𝑝 hp -III, and h ​ p ℎ 𝑝 hp -IV mesh refinement methods.

	$hp$-BB	$hp$-I	$hp$-II	$hp$-III	$hp$-IV
$M$	$2$	$10$	$4$	$5$	$5$
$N_f$	$40$	$85$	$86$	$115$	$95$
$𝒯$ (s)	$0.1234$	$0.4859$	$0.2826$	$0.4132$	$0.3054$

Table 3. Table 3: Mesh refinement performance results for Example 2 using h ​ p ℎ 𝑝 hp -BB, h ​ p ℎ 𝑝 hp -I, h ​ p ℎ 𝑝 hp -II, h ​ p ℎ 𝑝 hp -III, and h ​ p ℎ 𝑝 hp -IV mesh refinement methods.

	$hp$-BB	$hp$-I	$hp$-II	$hp$-III	$hp$-IV
$M$	$2$	$7$	$5$	$6$	$4$
$N_f$	$60$	$80$	$78$	$114$	$85$
$𝒯$ (s)	$0.2183$	$0.6234$	$0.5000$	$0.8461$	$0.4156$

Table 4. Table 4: Mesh refinement performance results for Example 3 using h ​ p ℎ 𝑝 hp -BB, h ​ p ℎ 𝑝 hp -I, h ​ p ℎ 𝑝 hp -II, h ​ p ℎ 𝑝 hp -III, and h ​ p ℎ 𝑝 hp -IV mesh refinement methods.

	$hp$-BB	$hp$-I	$hp$-II	$hp$-III	$hp$-IV
$M$	$2$	$16$	$9$	$8$	$8$
$N_f$	$90$	$204$	$235$	$283$	$194$
$𝒯$ (s)	$0.7019$	$8.4331$	$4.0833$	$4.7948$	$3.6157$

Equations

$$\mathcal{J}=\mathcal{M}(\mathbf{y}(t_0),t_0,\mathbf{y}(t_f),t_f)+{\int }_{t_0}^{t_f}\mathcal{L}(\mathbf{y}(t),\mathbf{u}(t),t)dt,$$

$$\frac{d\mathbf{y}}{dt}=\mathbf{a}(\mathbf{y}(t),\mathbf{u}(t),t),$$

$${\mathbf{c}}_{\min }\le \mathbf{c}(\mathbf{y}(t),\mathbf{u}(t),t)\le {\mathbf{c}}_{\max },$$

$${\mathbf{b}}_{\min }\le \mathbf{b}(\mathbf{y}(t_0),t_0,\mathbf{y}(t_f),t_f)\le {\mathbf{b}}_{\max }.$$

$$\bigcup _{d=1}^Q{\mathcal{P}}_d=[t_0,t_f],\bigcap _{d=1}^Q{\mathcal{P}}_d=\{t_s^{[1]},\dots ,t_s^{[Q-1]}\},$$

$$\begin{array}[]{lcl}t&=&\displaystyle\frac{t_{s}^{[d]}-t_{s}^{[d-1]}}{2}\tau+\frac{t_{s}^{[d]}+t_{s}^{[d-1]}}{2}~{},\vspace{3pt}\\ \tau&=&\displaystyle 2\frac{t-t_{s}^{[d-1]}}{t_{s}^{[d]}-t_{s}^{[d-1]}}-1~{}.\end{array}$$

$$\bigcup _{k=1}^K{\mathcal{S}}_k=[-1,+1],\bigcap _{k=1}^K{\mathcal{S}}_k=\{T_1,\dots ,T_{K-1}\},$$

$${\mathbf{y}}^{(k)}(\tau )\approx {\mathbf{Y}}^{(k)}(\tau )=\sum _{j=1}^{N_k+1}{\mathbf{Y}}_j^{(k)}{\ell }_j^{(k)}(\tau ),{\ell }_j^{(k)}(\tau )=\prod _{\stackrel{l=1}{l\ne j}}^{N_k+1}\frac{\tau -{\tau }_l^{(k)}}{{\tau }_j^{(k)}-{\tau }_l^{(k)}},$$

$$\frac{d{\mathbf{Y}}^{(k)}(\tau )}{d\tau }=\sum _{j=1}^{N_k+1}{\mathbf{Y}}_j^{(k)}\frac{d{\ell }_j^{(k)}(\tau )}{d\tau }.$$

$$\sum _{j=1}^{N_k+1}{D}_{lj}^{(k)}{\mathbf{Y}}_j^{(k)}-\frac{t_f-t_0}{2}\mathbf{a}\left({\mathbf{Y}}_l^{(k)},{\mathbf{U}}_l^{(k)},t({\tau }_l^{(k)},t_0,t_f)\right)=0,l\in \{1,\dots ,N_k\},$$

$$D_{lj}^{(k)}=\frac{d{\ell }_j^{(k)}({\tau }_l^{(k)})}{d\tau },l\in \{1,\dots ,N_k\},j\in \{1,\dots ,N_k+1\},$$

$${\mathbf{Y}}_{N^{[d-1]}+1}^{[d-1]}={\mathbf{Y}}_1^{[d]},$$

$$N^{[d]}=\sum _{k=1}^{K^{[d]}}{N}_k^{[d]}.$$

$$\mathcal{J}=\mathcal{M}({\mathbf{Y}}_1^{[1]},t_0,{\mathbf{Y}}_{N^{[Q]}+1}^{[Q]},t_f)+\sum _{d=1}^Q\frac{t_s^{[d]}-t_s^{[d-1]}}{2}{\left[{\mathbf{w}}^{[d]}\right]}^{\mathsf{T}}{\mathbf{L}}^{[d]},$$

$${\mathbf{\Delta }}^{[d]}={\mathbf{D}}^{[d]}{\mathbf{Y}}^{[d]}-\frac{t_s^{[d]}-t_s^{[d-1]}}{2}{\mathbf{A}}^{[d]}=0,d\in \{1,\dots ,Q\},$$

$${\mathbf{c}}_{\min }\le {\mathbf{C}}_j^{[d]}\le {\mathbf{c}}_{\max },j\in \{1,\dots ,N^{[d]}\},d\in \{1,\dots ,Q\},$$

$${\mathbf{b}}_{\min }\le \mathbf{b}({\mathbf{Y}}_1^{[1]},t_0,{\mathbf{Y}}_{N^{[Q]}+1}^{[Q]},t_f)\le {\mathbf{b}}_{\max }.$$

$${\mathbf{Y}}_{N^{[d-1]}+1}^{[d-1]}={\mathbf{Y}}_1^{[d]},d\in \{2,\dots ,Q\},$$

$${\mathbf{A}}^{[d]}=\left[\begin{array}{c}\mathbf{a}\left({\mathbf{Y}}_1^{[d]},{\mathbf{U}}_1^{[d]},t_1^{[d]}\right)\\ ⋮\\ \mathbf{a}\left({\mathbf{Y}}_{N^{[d]}}^{[d]},{\mathbf{U}}_{N^{[d]}}^{[d]},t_{N^{[d]}}^{[d]}\right)\end{array}\right]\in {\mathbb{R}}^{N^{[d]}\times n_y},$$

$${\mathbf{C}}^{[d]}=\left[\begin{array}{c}\mathbf{c}\left({\mathbf{Y}}_1^{[d]},{\mathbf{U}}_1^{[d]},t_1^{[d]}\right)\\ ⋮\\ \mathbf{c}\left({\mathbf{Y}}_{N^{[d]}}^{[d]},{\mathbf{U}}_{N^{[d]}}^{[d]},t_{N^{[d]}}^{[d]}\right)\end{array}\right]\in {\mathbb{R}}^{N^{[d]}\times n_c},$$

$${\mathbf{L}}^{[d]}=\left[\begin{array}{c}\mathcal{L}\left({\mathbf{Y}}_1^{[d]},{\mathbf{U}}_1^{[d]},t_1^{[d]}\right)\\ ⋮\\ \mathcal{L}\left({\mathbf{Y}}_{N^{[d]}}^{[d]},{\mathbf{U}}_{N^{[d]}}^{[d]},t_{N^{[d]}}^{[d]}\right)\end{array}\right]\in {\mathbb{R}}^{N^{[d]}\times 1},$$

$${\mathbf{Y}}^{[d]}=\left[\begin{array}{c}{\mathbf{Y}}_1^{[d]}\\ ⋮\\ {\mathbf{Y}}_{N^{[d]}+1}^{[d]}\end{array}\right]\text{ and }{\mathbf{U}}^{[d]}=\left[\begin{array}{c}{\mathbf{U}}_1^{[d]}\\ ⋮\\ {\mathbf{U}}_{N^{[d]}}^{[d]}\end{array}\right],$$

$${\mathbf{\lambda }}^{[d]}=({\mathbf{W}}^{[d]}{)}^{-1}{\mathbf{\Lambda }}^{[d]},{\mathbf{\lambda }}_{N^{[d]}+1}^{[d]}=({\mathbf{D}}_{N^{[d]}+1}^{[d]}{)}^{\mathsf{T}}{\mathbf{\Lambda }}^{[d]},$$

$$\mathcal{H}(\mathbf{y}(t),\mathbf{\lambda }(t),\mathbf{u}(t),t)=\mathcal{L}(\mathbf{y}(t),\mathbf{u}(t),t)+\mathbf{\lambda }(t){\mathbf{a}}^{\mathsf{T}}(\mathbf{y}(t),\mathbf{u}(t),t).$$

$$\mathcal{H}(\mathbf{y}(t),\mathbf{\lambda }(t),\mathbf{u}(t),t)=f(\mathbf{y}(t),\mathbf{\lambda }(t),\mathbf{z}(t),t)+\mathbf{\sigma }(\mathbf{y}(t),\mathbf{\lambda }(t),t){\mathbf{v}}^{\mathsf{T}}(t),$$

$$u_{i}^{*}(t)=\left\{\begin{array}[]{lcl}u_{i}^{\min}&,&\sigma_{i}({\bf{y}}(t),\boldsymbol{\lambda}(t),t)>0,\\ u_{i}^{\max}&,&\sigma_{i}({\bf{y}}(t),\boldsymbol{\lambda}(t),t)<0,\end{array}\right.\quad(i=1,\ldots,I)~{},$$

$$(t_d^i{)}^{(k)}=\frac{(t_{\sigma }^i{)}^{(k)}+(t_u^i{)}^{(k)}}{2},$$

$$\mathcal{J}=r_1{N}_1(t)+r_2{N}_2(t)+r_3{N}_3(t)+{\int }_0^T{u}_1(t)dt,$$

$$\begin{array}[]{l}\begin{array}[]{lcl}\dot{N}_{1}(t)&=&-a_{1}N_{1}(t)+2a_{3}N_{3}(t)(1-u_{1}(t))~{},\\ \dot{N}_{2}(t)&=&-a_{2}N_{2}(t)u_{2}(t)+a_{1}N_{1}(t)~{},\\ \dot{N}_{3}(t)&=&a_{3}N_{3}(t)+a_{2}N_{2}(t)u_{2}(t)~{},\end{array}\end{array}$$

$$\begin{array}[]{rclcl}0&\leq&u_{1}(t)&\leq&1~{},\\ u_{2}^{\min}&\leq&u_{2}(t)&\leq&1~{},\end{array}$$

Full text

Nomenclature

1 Introduction

Optimal control problems arise frequently in many engineering applications due to the need to optimize performance of a controlled dynamical system. In general, optimal control problems do not have analytic solutions and, thus, must be solved numerically. Numerical methods for optimal control fall into two broad categories: indirect methods and direct methods. In an indirect method, the first-order variational optimality conditions are derived, and the optimal control problem is converted to a Hamiltonian boundary-value problem (HBVP). The HBVP is then solved numerically using a differential-algebraic equation solver. In a direct method, the state and control are approximated, and the optimal control problem is transcribed into a finite-dimensional nonlinear programming problem (NLP) [1]. The NLP is then solved numerically using well-developed software [2, 3, 4].

Over the past two decades, a particular class of direct methods, called direct collocation methods, has been used extensively for solving continuous optimal control problems. A direct collocation method is an implicit simulation method where the state and control are parameterized, and the constraints in the continuous optimal control problem are enforced at a specially chosen set of collocation points. Traditional direct collocation methods take the form of an $h$ method (for example, Euler or Runge-Kutta methods) where the domain of interest is divided into a mesh, the state is approximated using the same fixed-degree polynomial in each mesh interval, and convergence is achieved by increasing the number and placement of the mesh points [1]. In contrast to an $h$ method, in recent years so-called $p$ methods have been developed. In a $p$ method, the number of intervals is fixed, and convergence is achieved by increasing the degree of the approximation in each interval. To achieve maximum effectiveness, $p$ methods have been developed using orthogonal collocation at Gaussian quadrature points [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. For problems whose solutions are smooth and well-behaved, a Gaussian quadrature orthogonal collocation method converges at an exponential rate [20, 21, 22, 23, 24]. Gauss quadrature collocation methods use either Legendre-Gauss (LG) points [8, 12, 13, 20, 21, 22, 23], Legendre-Gauss-Radau (LGR) points [11, 12, 13, 14, 20, 23, 24], or Legendre-Gauss-Lobatto (LGL) points [7].

Various $h$ or $p$ direct collocation methods have been developed previously. Reference [10] describes what is essentially a $p$ method where a differentiation matrix is used to identify switches, kinks, corners, and other discontinuities in the solution. Reference [25] develops a fixed-order method that uses a density function to generate a sequence of non-decreasing size meshes on which to solve the optimal control problem. Finally, in Ref. [1] an error estimate is developed by integrating the difference between an interpolation of the time derivative of the state and the right-hand side of the dynamics. The error estimate developed in Ref. [1] is predicated on the use of a fixed-order method (for example, trapezoid, Hermite-Simpson, Runge-Kutta) and computes a low-order approximation of the integral of the aforementioned difference.

Although $h$ methods have been used extensively and $p$ methods are useful on certain types of problems, both the $h$ and $p$ approaches have limitations. In the case of an $h$ method, it may be required to use an extremely fine mesh to improve accuracy. In the case of a $p$ method, it may be required to use an unreasonably large degree polynomial to improve accuracy. In order to reduce significantly the size of the finite-dimensional approximation, and thus improve computational efficiency of solving the NLP, in recent years the new class of $hp$ collocation methods has been developed for solving an optimal control problem. In an $hp$ method, both the number of mesh intervals and the degree of the approximating polynomial within each mesh interval are allowed to vary. While $hp$ methods were originally developed as finite-element methods for solving partial differential equations [26, 27, 28, 29, 30], over the past several years $hp$ methods have been developed for solving optimal control problems [15, 16, 19, 31, 32]. The methods described in Refs. [15] and [16] describe methods where the error estimate is based on the difference between an approximation of the time derivative of the state and the right-hand side of the dynamics midway between the collocation points. Next, Ref. [19] develops an error estimate based on the difference between the state interpolated using an increased number of Legendre-Gauss-Radau points in each mesh interval and the state obtained by integrating the dynamics on the solution using the interpolated state and control. Similar to the methods of Refs. [15] and [16], however, the method of Ref. [19] can only increase the size of the mesh. Additionally, Ref. [31] develops a method that adjusts the number of mesh intervals and the degree of the approximating polynomial within a mesh interval based on a proven convergence rate of the $hp$-adaptive LGR collocation method. Moreover, Ref. [31] describes an approach for reducing the size of the mesh. Finally, Ref. [32] describes an approach for adjusting the mesh based on the decay rate of the coefficients of a Legendre polynomial expansion of the state.

While the aforementioned $hp$ mesh refinement methods as described in Refs. [15, 16, 19, 31, 32] can improve accuracy and computational efficiency when compared with traditional $h$ or $p$ methods, a key missing aspect of these methods is that they do not exploit the structure in the optimal solution. In particular, in the case where the optimal solution is nonsmooth, rather than refining the mesh based on knowledge of the structure of the solution, these previously developed methods improve accuracy by increasing the number of collocation points in the vicinity of a discontinuity. As a result, these methods often place an undesirably large number of collocation points in the neighborhood of a discontinuity. Consequently, optimal control problems whose controls are discontinuous can require a large number of mesh refinement iterations to converge using the previously developed $hp$-adaptive methods.

A particular class of optimal control problems whose optimal solutions are discontinuous is the class of bang-bang optimal control problems. Bang-bang optimal control problems arise in a wide variety of well-known application areas [33]. A key feature of an optimal control problem with a bang-bang optimal control is that the Hamiltonian is linear with respect to one or more components of the control. Due to the linear dependence of the Hamiltonian on the control and under the assumption that the solution does not contain any singular arcs, Pontryagin’s minimum principle applies such that the optimal control lies at either its minimum or maximum limit. In the context of a mesh refinement method using collocation, if the method could algorithmically detect the bang-bang structure of the optimal solution, it may be possible to obtain an accurate solution in a more computationally efficient manner than would be possible using a standard mesh refinement method that does not exploit the structure of the solution.

The objective of this research is to develop a mesh refinement method for solving bang-bang optimal control problems by algorithmically exploiting the structure of the optimal solution. In particular, this research focuses on the development of a method that significantly improves computational efficiency while simultaneously reducing the mesh size and the number of mesh refinement iterations required in order to obtain a solution to a bang-bang optimal control problem. Previous research for solving bang-bang optimal control problems has been conducted using indirect methods. Refs. [34, 35, 36, 37, 38, 39, 40, 41] employ indirect shooting while Refs. [40, 41] include the second-order optimality conditions for bang-bang optimal control problems. Furthermore, Refs. [42, 43, 44] employ a direct shooting method for solving bang-bang optimal control problems by parameterizing the control using piecewise constants where the durations of the intervals are added as optimization parameters in order to solve for the bang-bang control profile. In particular, Ref. [44] combines the switch time computation and the time-optimal switching method developed in Refs. [42] and [43], respectively. In the context of $hp$-adaptive mesh refinement methods, knotting methods have been developed in Refs. [45, 46] which allow discontinuities in optimal control profiles to be taken into account by explicitly introducing a switch time variable to the problem definition. In particular, Ref. [45] utilizes the concept of super-elements which introduce switch time variables that are taken into account as parameters in the resulting optimization process. The modified Legendre pseudospectral scheme described in Ref. [33] also uses a similar approach by handling bang-bang optimal control problems using a knotting method that solves for the optimal control using an assumed number of switch times and constant control arcs. The assumed number of switch times is then increased or decreased for each iteration based on the approximated solution, with the scheme converging upon the solution for the optimal number of switch times. Additionally, Ref. [47] describes a slight modification of the $hp$-adaptive mesh refinement derived in Ref. [32], where new mesh points are placed at discontinuity locations that are estimated based on the switching functions of the Hamiltonian of the optimal control problem. Finally, Ref. [48] describes a mesh refinement method that is used in conjunction with the aforementioned $hp$-adaptive mesh refinement methods and detects discontinuities via jump function approximations.

In this paper a new direct collocation mesh refinement method is developed for solving optimal control problems whose solutions have a bang-bang structure. While in principle the method of this paper can be used with any collocation method, in this paper the Legendre-Gauss-Radau collocation method [11, 12, 13, 14, 20, 23, 24] is employed because it produces accurate state, control, and costate approximations [11, 12, 13, 14, 24]. Moreover, the approach of this paper is designed to be used in conjunction with a previously developed $hp$ mesh refinement method such as those described in Refs. [15, 16, 19, 31, 32]. First, a solution is obtained on a coarse mesh. Next, using the solution on this coarse mesh, the costate is estimated at the collocation points using the methods developed in Refs. [12, 13, 14]. Then, the state and costate approximations on the coarse mesh are used to determine algorithmically if the Hamiltonian is linear with respect to one or more components of the control. Using the state and costate approximations, the switching functions are estimated at the collocation points for those components of the control for which the Hamiltonian depends upon linearly. The estimates of the switching functions are then used to estimate a discontinuity in the control between any two collocation points where a switching function changes sign. Using these estimates of the control discontinuities, the locations of the switch times in the control are then introduced as optimization variables and the optimal control problem is divided into multiple domains. Within each domain, those components of the control that have a bang-bang solution structure are fixed at either their lower or upper limits depending upon the sign of the switching function. It is important to note that those control components that do not have a bang-bang structure remain free to vary within their defined bounds. Furthermore, it is noted that the multiple-domain formulation is analogous to using super-elements as described in Ref. [45], with a key difference being that the structure of the bang-bang optimal control profile has been algorithmically detected using the estimates of the switching functions found on the initial mesh. The multiple-domain optimal control problem is then solved using LGR collocation where it is noted again that the switch times are determined as part of the optimization.

This paper is organized as follows. Section 2 defines the general single-phase optimal control problem in Bolza form. Section 3 describes the rationale for using Legend-Gauss-Radau collocation points as the set of nodes to discretize the continuous optimal control problem and the Legendre-Gauss-Radau collocation method. Section 4 briefly overviews the form of bang-bang optimal control problems. Section 5 describes the bang-bang mesh refinement method of this paper. Section 6 demonstrates the performance of the mesh refinement method of this paper when solving bang-bang optimal control problems as compared to the four previously developed $hp$-adaptive mesh refinement methods of Refs. [19, 15, 31, 32]. Section 7 provides a discussion of both the approach and the results. Finally, Section 8 provides conclusions on this research.

2 Single-Phase Optimal Control Problem

Without loss of generality, consider the following general single-phase optimal control problem in Bolza form defined on the time horizon $t\in[t_{0},t_{f}]$. Determine the state ${\bf{y}}(t)\in\mathbb{R}^{1~{}\times~{}n_{y}}$, the control ${\bf{u}}(t)\in\mathbb{R}^{1~{}\times~{}n_{u}}$, the start time $t_{0}\in\mathbb{R}$, and the terminus time $t_{f}\in\mathbb{R}$ that minimize the objective functional

[TABLE]

subject to the dynamic constraints

[TABLE]

the inequality path constraints

[TABLE]

and the boundary conditions

[TABLE]

3 Legendre-Gauss-Radau Collocation

In order to develop the method described in this paper, a direct collocation method must be chosen. While in principle any direct collocation method can be used to approximate the optimal control problem given in Section 2, in this research the previously developed Legendre-Gauss-Radau (LGR) collocation method [11, 12, 13, 14, 19, 24] will be used because it has been shown that the LGR collocation method produces an accurate state, control, and costate [11, 12, 13, 14, 24]. It is noted that the accuracy of the costate estimate of the LGR collocation method plays an important role in the method of this paper because the accuracy of the costate directly influences the accuracy of the estimates of the switching functions of the Hamiltonian which are in turn used to estimate the switch times in the optimal control (forming the basis of the method of this paper).

In the context of this research, a multiple-interval form of the LGR collocation method is chosen. The time horizon $t\in[t_{0},t_{f}]$ may be divided into $Q$ time domains, ${\cal{P}}_{d}=[t_{s}^{[d-1]},t_{s}^{[d]}]\subseteq[t_{0},t_{f}],~{}d\in\{1,\ldots,Q\}$, such that

[TABLE]

where $t_{s}^{[d]},~{}d\in\{1,\ldots,Q-1\}$ are the switch time variables of the problem, $t_{s}^{[0]}=t_{0}$, and $t_{s}^{[Q]}=t_{f}$. Thus in the case where $Q=1$ the domain consists of only a single domain ${\cal{P}}_{1}=[t_{0},t_{f}]$ and $\{t_{s}^{[1]},\ldots,t_{s}^{[Q-1]}\}=\emptyset$. In the multiple-interval LGR collocation method, each of the time domains ${\cal{P}}_{d}=[t_{s}^{[d-1]},t_{s}^{[d]}],~{}d\in\{1,\ldots,Q\}$, is converted into the domain $\tau\in[-1,+1]$ using the affine transformation,

[TABLE]

The interval $\tau\in[-1,+1]$ for each domain ${\cal{P}}_{d}$ is then divided into $K$ mesh intervals, ${\cal{S}}_{k}=[T_{k-1},T_{k}]\subseteq[-1,+1],\;k\in\{1,\ldots,K\}$ such that

[TABLE]

and $-1=T_{0}<T_{1}<\ldots<T_{K-1}<T_{K}=+1$. For each mesh interval, the LGR points used for collocation are defined in the domain of $[T_{k-1},T_{k}]$ for $k\in\{1,\ldots,K\}$. The state of the continuous optimal control problem is then approximated in mesh interval ${\cal{S}}_{k},\;k\in\{1,\ldots,K\}$, as

[TABLE]

where $\ell_{j}^{(k)}(\tau)$ for $j\in\{1,\ldots,N_{k}+1\}$ is a basis of Lagrange polynomials on ${\cal{S}}_{k}$, $\left(\tau_{1}^{(k)},\ldots,\tau_{N_{k}}^{(k)}\right)$ are the set of $N_{k}$ Legendre-Gauss-Radau (LGR) [49] collocation points in the interval $[T_{k-1},T_{k})$, $\tau_{N_{k}+1}^{(k)}=T_{k}$ is a non-collocated support point, and ${\bf{Y}}_{j}^{(k)}\equiv{\bf{Y}}^{(k)}(\tau_{j}^{(k)})$. Differentiating ${\bf{Y}}^{(k)}(\tau)$ in Eq. (8) with respect to $\tau$ gives

[TABLE]

The dynamics are then approximated at the $N_{k}$ LGR points in mesh interval $k\in\{1,\ldots,K\}$ as

[TABLE]

where

[TABLE]

are the elements of the $N_{k}\times(N_{k}+1)$ Legendre-Gauss-Radau differentiation matrix [12] in mesh interval ${\cal{S}}_{k}$, $\;k\in\{1,\ldots,K\}$, and ${\bf{U}}_{l}^{(k)}$ is the approximation of the control at the $l^{th}$ collocation point in mesh interval ${\cal{S}}_{k}$. It is noted that continuity in the state and time between mesh intervals ${\cal{S}}_{k-1}$ and ${\cal{S}}_{k}$, $k\in\{1,\ldots,K\}$, is enforced by using the same variables to represent ${\bf{Y}}_{N_{k-1}+1}^{(k-1)}$ and ${\bf{Y}}_{1}^{(k)}$, while continuity in the state between the domains ${\cal{P}}_{d-1}$ and ${\cal{P}}_{d}$, $d\in\{2,\ldots,Q\}$, is achieved using the additional continuity constraint

[TABLE]

where the superscript $[d]$ is used to denote the $d^{th}$ time domain, ${\bf{Y}}_{j}^{[d]}$ denotes the value of the state approximation at the $j^{th}$ discretization point in the time domain ${\cal{P}}_{d}$, and $N^{[d]}$ is the total number of collocation points used in time domain ${\cal{P}}_{d}$ computed by

[TABLE]

The Legendre-Gauss-Radau approximation of the multiple-domain optimal control problem then leads to the following nonlinear programming problem (NLP). Minimize the objective function

[TABLE]

subject to the defect constraints

[TABLE]

the path constraints

[TABLE]

the boundary conditions

[TABLE]

and the continuity constraints

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

${\bf{D}}^{[d]}\in\mathbb{R}^{N^{[d]}~{}\times~{}[N^{[d]}+1]}$ is the LGR differentiation matrix in time domain ${\cal{P}}_{d},~{}d\in\{1,\ldots,Q\}$, and ${\bf{w}}^{[d]}\in\mathbb{R}^{N^{[d]}~{}\times~{}1}$ are the LGR weights at each node in time domain ${\cal{P}}_{d},~{}d\in\{1,\ldots,Q\}$. It is noted that ${\bf{a}}\in\mathbb{R}^{1~{}\times~{}n_{y}}$, ${\bf{c}}\in\mathbb{R}^{1~{}\times~{}n_{c}}$, and ${\cal{L}}\in\mathbb{R}^{1~{}\times~{}1}$ correspond, respectively, to the vector fields that define the right-hand side of the dynamics, the path constraints, and the integrand of the optimal control problem, where $n_{y}$ and $n_{c}$ are, respectively, the number of state components and path constraints in the problem. Additionally, the state matrix, ${\bf{Y}}^{[d]}\in\mathbb{R}^{[N^{[d]}+1]~{}\times~{}n_{y}}$, and the control matrix, ${\bf{U}}^{[d]}\in\mathbb{R}^{N^{[d]}~{}\times~{}n_{u}}$, in time domain ${\cal{P}}_{d},~{}d\in\{1,\ldots,Q\}$, are formed as

[TABLE]

respectively, where $n_{u}$ is the number of control components in the problem. Finally, as described in Ref. [12], estimates of the costate may be obtained at each of the discretization points in the time domain ${\cal{P}}_{d},d\in\{1,\ldots,Q\}$, using the transformation

[TABLE]

where $\boldsymbol{\lambda}^{[d]}\in\mathbb{R}^{N^{[d]}~{}\times~{}n_{y}}$ is a matrix of the costate estimates at the collocation points in time domain ${\cal{P}}_{d}$, ${\bf{W}}^{[d]}=\text{diag}({\bf{w}}^{[d]})$ is a diagonal matrix of the LGR weights at the collocation points in time domain ${\cal{P}}_{d}$, $\boldsymbol{\Lambda}^{[d]}\in\mathbb{R}^{N^{[d]}~{}\times~{}n_{y}}$ is a matrix of the NLP multipliers obtained from the NLP solver corresponding to the defect constraints at the collocation points in time domain ${\cal{P}}_{d}$ , $\boldsymbol{\lambda}_{N^{[d]}+1}^{[d]}\in\mathbb{R}^{1~{}\times~{}n_{y}}$ is a row vector of the costate estimates at the non-collocated end point in time domain ${\cal{P}}_{d}$, and ${\bf{D}}_{N^{[d]}+1}^{[d]}\in\mathbb{R}^{N^{[d]}~{}\times~{}1}$ is the last column of the LGR differentiation matrix in time domain ${\cal{P}}_{d}$.

4 Control-Linear Hamiltonian

The Hamiltonian of the Bolza optimal control problem defined in Eqs. (1)–(4) is given as

[TABLE]

Assume now that the Hamiltonian given in Eq. (23) has the following form:

[TABLE]

where $f\in\mathbb{R}$, $\boldsymbol{\sigma}\in\mathbb{R}^{1~{}\times~{}I}$, ${\bf{v}}(t)=(u_{1}(t),\ldots,u_{I}(t))\in\mathbb{R}^{1\times I}$ and ${\bf{z}}(t)=(u_{I+1}(t),\ldots,u_{n_{u}}(t))\in\mathbb{R}^{1\times(n_{u}-I)}$. It is noted that ${\bf{v}}(t)$ and ${\bf{z}}(t)$ are vectors that correspond, respectively, to those components of the control upon which the Hamiltonian depends linearly and nonlinearly, while $\boldsymbol{\sigma}=(\sigma_{1},\ldots,\sigma_{I})$ is a vector that defines the $I$ switching functions. In other words, the components $(u_{1}(t),\ldots,u_{n_{u}}(t))$ of the control are ordered such that any control-linear component has an index that is lower than any control-nonlinear component and $(\sigma_{1},\ldots,\sigma_{I})$ are functions of the state and the costate that are used to determine locations where a switch in the control may occur. Then, assuming that the optimal solution contains no singular arcs, the optimal control for any control-linear component, $u_{i}^{*}(t),~{}i\in\{1,\ldots,I\}$, is obtained from Pontryagin’s minimum principle as

[TABLE]

where the signs of the switching functions, $\sigma_{i},~{}i\in\{1,\ldots,I\}$, are determined by the state ${\bf{y}}(t)$, the costate $\boldsymbol{\lambda}(t)$, and the time $t$. Thus, the optimal control will have a bang-bang structure where discontinuities in the solution occur whenever a switching function, $\sigma_{i},~{}i\in\{1,\ldots,I\}$, changes sign (again, assuming no singular arcs).

Assume now that the optimal control problem has been approximated using the LGR collocation method as described in Section 3. Furthermore, assume that the NLP resulting from LGR collocation has been solved to obtain estimates of the state, control, and costate as given in Eqs. (21) and (22). Then the approximation can be utilized to detect possible switch times in the optimal control. These possible discontinuity locations can be detected by evaluating the switching functions $\sigma_{i},~{}i\in\{1,\ldots,I\}$, at each of the collocation points and checking for sign changes in $\sigma_{i}$ between two adjacent collocation points. In particular, a sign change in $\sigma_{i}$ between any two adjacent collocation points $t_{k}$ and $t_{k+1}$ indicates that $\sigma_{i}$ has a root in the time interval $t\in[t_{k},t_{k+1}]$, thus indicating a switch in the control at a time $t\in[t_{k},t_{k+1}]$. As discussed in Ref. [50], the rationale for using changes in the sign of $\sigma_{i}$ to detect switches in the control is due to the fact that the gradient of the objective with respect to a switching point in a bang-bang control is proportional to the coefficient $\sigma_{i}$ of the control in the Hamiltonian. Consequently, the gradient of the objective with respect to a switching point only vanishes at values of time where the coefficient of the control vanishes, and the change in sign of $\sigma_{i}$ between a pair of points provides an estimate of a switch time in the component of the control corresponding to the switching function $\sigma_{i}$.

5 Bang-Bang Mesh Control Refinement Method

In this section, the bang-bang mesh refinement method is developed. The method consists of several steps. The first step, described in Section 5.1 determines if the Hamiltonian is linear in one or more components of the control. The second step, described in Section 5.2 computes estimates of the locations of discontinuities in each component of the control upon which the Hamiltonian depends upon linearly. The third step, described in Section 5.3, reformulates the optimal control problem into a multiple-domain optimal control problem where the domains are determined based on the estimates of the discontinuity locations obtained in the second step of the method. The interior endpoints of these domains are then treated as optimization variables by introducing an appropriate number of switch time parameters which use the estimates of the discontinuity locations as initial guesses for these new variables. In addition, any component of the control for which the Hamiltonian has been determined to depend upon linearly is then fixed at either its minimum or maximum value in each domain depending upon the direction of the switch in sign of the corresponding switching function component identified in the second step. The multiple-domain optimal control problem is then solved iteratively using the LGR collocation method described in Section 3 together with a previously developed $hp$-adaptive mesh refinement method (to refine the mesh based on the error in the solution for those control components for which the Hamiltonian does not depend upon linearly and provide sufficient resolution in the state approximation). Finally, Section 5.4 provides a summary of the bang-bang mesh refinement method.

5.1 Method for Identifying Bang-Bang Optimal Control Problems

Assume now that the optimal control problem as formulated in Section 2 has been transcribed into a nonlinear programming problem using the LGR collocation method as developed in Section 3 assuming that $Q=1$ (that is, a single domain is used). Furthermore, assume that the resulting NLP has been solved to obtain estimates of the state, control, and costate as given in Eqs. (21) and (22), and that the mesh refinement accuracy tolerance has not been satisfied and, thus, mesh refinement is required. As a result, it is possible that the optimal solution may possess a bang-bang optimal control. The first step in determining if a bang-bang optimal control is a possibility is to determine if the Hamiltonian is linear in one or more of the components of the control. Moreover, the method should be able to detect a control-linear Hamiltonian without requiring any external intervention.

In this research, the determination of a control-linear Hamiltonian is made using hyper-dual derivative approximations [51]. In particular, using the state and costate approximations obtained from the LGR collocation method, the hyper-dual derivative approximation of Ref. [51] is used to compute the first- and second-derivatives of the Hamiltonian with respect to each component of the control. In order to identify linearity in the Hamiltonian with respect to a given control component, sample values of the control component are taken between the bounds of the control and partial derivatives of the Hamiltonian are then taken with respect to that control component while holding all other variables constant. If the partial derivatives of the Hamiltonian obtained using the various sample values of the control component are found to be constant, then the Hamiltonian is identified as possibly being linear with respect to that control component. Moreover, if any control component has been identified to have a zero second derivative and has zero second-order partial derivatives with respect to any other control components [that is, all cross partial derivatives $\frac{\partial^{2}H}{\partial u_{i}\partial u_{j}}=0,~{}(i\neq j)$], then the Hamiltonian is linear with respect to that component of the control and is thus a candidate for bang-bang control mesh refinement. Estimates of the switching functions of the Hamiltonian are then obtained by computing the partial of the Hamiltonian with respect to each control-linear component using hyper-dual derivative approximations. Finally, it is noted that the hyper-dual derivative approximation is not subject to truncation error up to second derivatives and, thus, provides exact (that is, to machine precision) first- and second-derivatives of the Hamiltonian with respect to the control.

5.2 Estimating Locations of Switches in Control

Once it has been determined that the Hamiltonian is linear in at least one control component, the next step is to estimate times at which a control component may switch between its lower and upper limit, thus leading to a discontinuity in the control solution profile. Assuming that the optimal solution contains no singular arcs, a discontinuity in the control will occur when a switching function $\sigma_{i},~{}i\in\{1,\ldots,I\}$ changes sign. Given that the solution of the optimal control problem has been approximated using the LGR collocation as given in Section 3, an estimate of a discontinuity in a control component will be when a switching function $\sigma_{i}$ changes sign between two adjacent collocation points. Furthermore, because any switching function, $\sigma_{i}$ may change sign one or more times, it is possible that any or all components of the control ${\bf{v}}(t)$ may have one or more discontinuities and that the discontinuity time for a particular component, $v_{i}(t)$, may differ from the discontinuity time of another control-linear component, $v_{j}(t),~{}(i\neq j)$.

In order to accommodate the possibility of multiple discontinuities within a given mesh interval, ${\cal{S}}_{k},~{}k\in\{1,\ldots,K\}$, each discontinuity estimate time $(t_{d}^{i})^{(k)}$ in ${\cal{S}}_{k}$ for a given control component $u_{i},~{}i\in\{1,\ldots,I\}$ is computed using

[TABLE]

where $(t_{\sigma}^{i})^{(k)}$ is the midpoint time between the two adjacent discretization points showing a change in sign of the corresponding switching function $\sigma_{i}$ within ${\cal{S}}_{k}$, and $(t_{u}^{i})^{(k)}$ is the midpoint time between the two adjacent discretization points displaying the largest absolute difference in the value of the $i^{th}$ control-linear component $u_{i}$ within ${\cal{S}}_{k}$. An example schematic of estimating the discontinuity time, $(t_{d}^{i})^{(k)}$, for control-linear component, $u_{i}$, for mesh interval ${\cal{S}}_{k}$ is shown in Fig. 1.

By taking into account both of the midpoint times $(t_{\sigma}^{i})^{(k)}$ and $(t_{u}^{i})^{(k)}$ in estimating the discontinuity time $(t_{d}^{i})^{(k)}$, the bang-bang control profile can be properly maintained relative to all control-linear components within a mesh interval ${\cal{S}}_{k}$ containing multiple discontinuities. Furthermore, in the event that the switching function $\sigma_{i}$ changes sign across two adjacent mesh intervals ${\cal{S}}_{k}$ and ${\cal{S}}_{k+1}$, the mesh point $T_{k}$ that lies at the interface between the mesh intervals ${\cal{S}}_{k}$ and ${\cal{S}}_{k+1}$ is used as the estimate for the discontinuity time $(t_{d}^{i})^{(k)}$ (that is $(t_{d}^{i})^{(k)}=T_{k}$) as shown in Fig. 2.

After checking all mesh intervals ${\cal{S}}_{k},~{}k\in\{1,\ldots,K\}$, for possible discontinuities, the computed discontinuity estimates, $(t_{d}^{i})^{(k)},~{}i\in\{1,\ldots,I\},~{}k\in\{1,\ldots,K\}$, are arranged in ascending order and used as initial guesses for the switch time parameters $t_{s}^{[S]},~{}S\in\{1,\ldots,n_{s}\}$, that are to be introduced in the subsequent mesh iterations, where $n_{s}$ is equal to the total number of discontinuities detected in the solution obtained on the initial mesh. Finally, the limit of each control-linear component $u_{i},~{}i\in\{1,\ldots,I\},$ to be used on the subsequent mesh iterations in each of the newly created time domains ${\cal{P}}_{d},~{}d\in\{1,\ldots,Q\}$, is identified by checking the sign of the corresponding switching function $\sigma_{i}$ within ${\cal{P}}_{d}$, while each control-nonlinear component $u_{i},~{}i\in\{I+1,\ldots,n_{u}\},$ is left free to vary between its defined bounds, where $Q=n_{s}+1$ on the subsequent mesh iterations. Thus, by using the estimates of the switching functions, the structure of the bang-bang optimal control profile has been automatically detected and used to set up the subsequent mesh iterations to include $n_{s}$ switch time parameters.

5.3 Reformulation of Optimal Control Problem Into Multiple Domains

Assuming the optimal control problem has been identified suitable for bang-bang mesh refinement as described in Section 5.1, the method for automatically detecting the structure of the bang-bang control as described in Section 5.2 may be employed. Once acquired, the detected structure of the bang-bang control may be used to introduce the appropriate number of switch time parameters, $t_{s}^{[S]},~{}S\in\{1,\ldots,n_{s}\}$, to be solved for on subsequent mesh iterations, where the initial guess for each switch time parameter $t_{s}^{[S]}$ is the estimated discontinuity time that was found using the process described in Section 5.2. The switch time parameters are included by using variable mesh points between the newly created time domains, ${\cal{P}}_{d},~{}d\in\{1,\ldots,Q\}$, where $Q=n_{s}+1$ . Specifically, the variable mesh points are employed by dividing the time horizon $t=[t_{0},t_{f}]$ of the optimal control problem identified as a bang-bang control into $Q=n_{s}+1$ time domains as described in Section 3, where $n_{s}$ is the number of switch time parameters introduced based on the detected structure of the bang-bang control profile for the solution on the initial mesh. Each of the time domains ${\cal{P}}_{d}=[t_{s}^{[d-1]},t_{s}^{[d]}],~{}d\in\{1,\ldots,Q\}$, has bounds enforced on $t_{s}^{[d-1]}$ and $t_{s}^{[d]}$ which appropriately bracket the estimated discontinuity times detected from the structure. Additionally, the bounds of the control-linear components $u_{i},~{}i\in\{1,\ldots,I\},$ within each time domain ${\cal{P}}_{d},~{}d\in\{1,\ldots,Q\}$, are set to the identified constant bang control limit for each control-linear component during the corresponding spans of time based on the previously detected structure, while the bounds of the control-nonlinear components $u_{i},~{}i\in\{I+1,\ldots,n_{u}\},$ are left unchanged. Thus the bang-bang optimal control problem is effectively transcribed into a multiple-domain optimal control problem employing constant values for the control-linear components in each time domain ${\cal{P}}_{d}~{},d\in\{1,\ldots,Q\},$ such that the optimal switch times are solved for on the subsequent mesh iterations while the control-nonlinear components are left free to vary between their respective bounds. A schematic for the process of dividing the identified bang-bang optimal control problem into a multiple-domain optimal control problem employing constant values for the control-linear components $u_{i},~{}i\in\{1,\ldots,I\},$ in each time domain is shown in Fig. 3.

5.4 Summary of Bang-Bang Control Mesh Refinement Method

A summary of the bang-bang control mesh refinement method developed in this paper appears below. Here $M$ denotes the mesh number, and in each loop of the method, the mesh number increases by $1$. The method terminates in Step 4 when the error tolerance is satisfied or when $M$ reaches a prescribed maximum $M_{\max}$.

\MakeFramed\FrameRestore

Bang-Bang Control Mesh Refinement Method

1. Step 1:

Specify initial mesh. 2. Step 2:

Solve LGR collocation NLP of Eqs. (13)–(17) on current mesh. 3. Step 3:

Compute relative error $e$ on current mesh. 4. Step 4:

If $e<\epsilon$ or $M>M_{\max}$, then quit. Otherwise, go to Step 5. 5. Step 5:

If $M=1$, determine the number $I$ of control components for which Hamiltonian is linear in the control in manner described in Section 5.1. Otherwise, go to Step 6. 6. Step 6:

If $I=0$ or $M>1$, employ standard mesh refinement and return to Step 2. 7. Step 7:

If $I>0$ and $M=1$, employ bang-bang mesh refinement using the following steps:

1. (a)

Estimate $n_{s}$ discontinuities in control components using method described in Section 5.2. 2. (b)

Partition time domain into multiple ($Q=n_{s}+1$) domains using method of Section 5.3. 3. (c)

Solve multiple-domain optimal control problem that includes following features:

1. (i)

Include $n_{s}$ variable mesh points in multiple-domain formulation. 2. (ii)

Fix $I$ bang-bang control components at either lower/upper limit in each of the $Q$ domains. 4. (d)

Increment $M\longrightarrow M+1$ and return to Step 2.

\endMakeFramed

6 Examples

In this section, three nontrivial bang-bang optimal control problems are solved using the bang-bang mesh refinement method described in Section 5. The first example is the three compartment model problem taken from Ref. [40]. The second example is the robot arm problem taken from Ref. [52]. The third example is the free-flying robot problem taken from Ref. [53]. The efficiency of the bang-bang mesh refinement method developed in this paper is evaluated and compared against four previously developed $hp$-adaptive mesh refinement methods described in Refs. [19, 15, 31, 32]. For problems requiring more than a single mesh refinement to meet accuracy tolerance, the bang-bang mesh refinement method will utilize the $hp$-adaptive method described in Ref. [32] to further refine the phases not meeting the specified mesh accuracy tolerance. It is noted that any of the four previously developed $hp$-adaptive mesh refinement methods may be paired with the bang-bang mesh refinement method. An initial coarse mesh of ten intervals with five collocation points in each interval is used for each problem. Furthermore, upon identification of the bang-bang control solution profile, the bang-bang mesh refinement method initially employs two mesh intervals to discretize each of the newly created time domains employing constant values for the control-linear components, with five collocation points in each mesh interval for the first and second examples and six collocation points in each mesh interval for the third example. The performance of the mesh refinement methods are evaluated based on the number of mesh iterations, total number of collocation points used, and the total computation time for the problem to be solved satisfactorily for the specified mesh accuracy tolerance. All plots were created using MATLAB Version R2016a (build 9.0.0.341360).

The following conventions are used for all of the examples. First, $M$ denotes the number of mesh refinement iterations required to meet the mesh refinement accuracy tolerance of $\epsilon=10^{-6}$, where $M=1$ corresponds to the initial mesh. Second, $N_{f}$ denotes the total number of collocation points on the final mesh. Third, ${\cal{T}}$ denotes the total computation time required to solve the optimal control problem. Furthermore, the $hp$-adaptive mesh refinement methods described in Refs. [19, 15, 31, 32] are referred to, respectively, as the $hp$-I, $hp$-II, $hp$-III, and $hp$-IV mesh refinement methods. Additionally, the bang-bang mesh refinement method developed in this paper is referred to as the $hp$-BB mesh refinement method. Finally, for all mesh refinement methods, a minimum and maximum number of three and ten collocation points in each mesh interval is enforced, respectively. All results shown in this paper were obtained using the C++ optimal control software $\mathbb{CGPOPS}$ [54] using the NLP solver IPOPT [4] in full Newton (second-derivative) mode with an NLP solver optimality tolerance of $10^{-9}$. All first- and second-derivatives required by the NLP solver were computed using the hyper-dual derivative approximation method described in Ref. [51]. All computations were performed on a 2.9 GHz Intel Core i7 MacBook Pro running MAC OS-X version 10.13.6 (High Sierra) with 16GB 2133MHz LPDDR3 of RAM. C++ files were compiled using Apple LLVM version 9.1.0 (clang-1000.10.44.2).

6.1 Example 1: Three Compartment Model Problem

Consider the following optimal control problem taken from Ref. [40]. Minimize the objective functional

[TABLE]

subject to the dynamic constraints

[TABLE]

the control inequality constraints

[TABLE]

and the boundary conditions

[TABLE]

where $a_{1}=0.197$, $a_{2}=0.395$, $a_{3}=0.107$, $r_{1}=1$, $r_{2}=0.5$, $r_{3}=1$, $T=7$, and $u_{2}^{\min}=0.70$. The optimal control problem given in Eqs. (27) – (30) was solved using each of the five mesh refinement methods $hp$-BB, $hp$-I, $hp$-II, $hp$-III, and $hp$-IV. The solutions obtained using any of these five mesh refinement methods are in close agreement and match the solution given in Ref. [40] (see pages $200$ – $203$ of Ref. [40]), and a summary of the performance of each method is shown in Table 2. In particular, it is seen in Table 2 that the $hp$-BB method is more computationally efficient, requires fewer mesh iterations, and results in a smaller final mesh to meet the mesh refinement accuracy tolerance of $10^{-6}$ when compared with any of the other mesh refinement methods. Furthermore, as shown in Fig. 4, the solution obtained using the $hp$-BB mesh refinement method accurately captures the bang-bang control profile of the optimal solution. Specifically, it is seen that the optimal switch times are obtained to nearly machine precision. Furthermore, estimates of the switching functions obtained using the solution on the initial mesh are shown in Fig. 5, where it is seen that the roots of the switching functions are in close proximity to the locations of the optimal switch times. Thus, while the previously developed $hp$-adaptive mesh refinement methods are able to satisfy the mesh refinement accuracy tolerance in a reasonable number of iterations, the $hp$-BB method outperforms all of these methods. Moreover, unlike the other methods where a large number of collocation points are placed in the vicinity of a discontinuity (as seen in Fig. 4), the $hp$-BB method places no unnecessary collocation points at the discontinuity due to the fact that optimal control problem has been divided into multiple domains and variable mesh points are included that define the locations of the discontinuities in the control (again, see Fig. 4).

6.2 Example 2: Robot Arm Problem

Consider the following optimal control problem taken from Ref. [52]. Minimize the objective functional

[TABLE]

subject to the dynamic constraints

[TABLE]

the control inequality constraints

[TABLE]

and the boundary conditions

[TABLE]

where

[TABLE]

The optimal control problem given in Eqs. (31)–(35) was solved using each of the five mesh refinement methods $hp$-BB, $hp$-I, $hp$-II, $hp$-III, and $hp$-IV. The solutions obtained using any of these five mesh refinement methods are in close agreement and match the solution given in Ref. [52] (see page $20$ of Ref. [52]), and a summary of the performance of each method is shown in Table 3. In particular, it is seen in Table 3 that the $hp$-BB method is more computationally efficient, requires fewer mesh iterations, and results in a smaller final mesh to meet the mesh refinement accuracy tolerance of $10^{-6}$ when compared with any of the other mesh refinement methods. Furthermore, as shown in Fig. 6, the solution obtained using the $hp$-BB mesh refinement method accurately captures the bang-bang control profile of the optimal solution. Specifically, it is seen that the optimal switch times are obtained to nearly machine precision. Furthermore, estimates of the switching functions obtained using the solution on the initial mesh are shown in Fig. 7, where it is seen that the roots of the switching functions are in close proximity to the locations of the optimal switch times. Thus, while the previously developed $hp$-adaptive mesh refinement methods are able to satisfy the mesh refinement accuracy tolerance in a reasonable number of iterations, the $hp$-BB method outperforms all of these methods. Moreover, unlike the other methods, where a large number of collocation points are placed in the vicinity of a discontinuity (as seen in Fig. 6), the $hp$-BB method places no unnecessary collocation points at the discontinuity due to the fact that optimal control problem has been divided into multiple domains and variable mesh points are included that define the locations of the discontinuities in the control (again, see Fig. 6).

6.3 Example 3: Free-Flying Robot Problem

Consider the following optimal control problem taken from Ref. [53]. Minimize the cost functional

[TABLE]

subject to the dynamic constraints

[TABLE]

the control inequality constraints

[TABLE]

and the boundary conditions

[TABLE]

where $F_{1}(t)=u_{1}(t)-u_{2}(t)$ and $F_{2}(t)=u_{3}(t)-u_{4}(t)$ are the real control, $\alpha=0.2$, and $\beta=0.2$.

The optimal control problem given in Eqs. (36) – (39) was solved using each of the five mesh refinement methods $hp$-BB, $hp$-I, $hp$-II, $hp$-III, and $hp$-IV. The solutions obtained using any of these five mesh refinement methods are in close agreement and match the solution given in Ref. [52] (see pages $328$ – $329$ of Ref. [1]), and a summary of the performance of each method is shown in Table 4. In particular, it is seen in Table 4 that the $hp$-BB method is more computationally efficient, requires fewer mesh iterations, and results in a smaller final mesh to meet the mesh refinement accuracy tolerance of $10^{-6}$ when compared with any of the other mesh refinement methods. Furthermore, as shown in Fig. 8, the solution obtained using the $hp$-BB mesh refinement method accurately captures the bang-bang control profile of the optimal solution. Specifically, it is seen that the optimal switch times are obtained to nearly machine precision. Furthermore, estimates of the switching functions obtained using the solution on the initial mesh are shown in Fig. 9, where it is seen that the roots of the switching functions are in close proximity to the locations of the optimal switch times. Thus, while the previously developed $hp$-adaptive mesh refinement methods are able to satisfy the mesh refinement accuracy tolerance in a reasonable number of iterations, the $hp$-BB method outperforms all of these methods. Moreover, unlike the other methods, where a large number of collocation points are placed in the vicinity of a discontinuity (as seen in Fig. 8), the $hp$-BB method places no unnecessary collocation points at the discontinuity due to the fact that optimal control problem has been divided into multiple domains and variable mesh points are included that define the locations of the discontinuities in the control (again, see Fig. 8).

7 Discussion

The results of Section 6 demonstrate the effectiveness of the mesh refinement method developed in this paper for problems whose optimal control has a bang-bang structure. In particular, the results of Section 6 show that, while the previously developed mesh refinement methods are able to find a solution that meets a specified mesh refinement accuracy tolerance, these methods place an unnecessarily large number of collocation points in the vicinity of a discontinuity in the control. In addition, these methods often require a large amount of mesh refinement to meet a desired accuracy tolerance. On the other hand, for problems whose optimal control has a bang-bang structure, the mesh refinement method developed in this paper locates discontinuities accurately. This improved accuracy (over a standard mesh refinement method) is due to the fact that accurate estimates are obtained of the switching functions associated with those components of the control upon which the Hamiltonian depends linearly (where an accurate estimate of the switching functions is obtained because the costate of the optimal control problem is approximated accurately using the LGR collocation method). Then, by partitioning the horizon into multiple domains, introducing variables that define the locations of the switch times in the control-linear components, and fixing the control-linear components to lie at either its lower or upper limit in each domain, the method developed in this paper accurately identifies the switch times. Moreover, solutions that meet the specified accuracy tolerance are obtained in fewer mesh refinement iterations when compared with using one of the previously developed mesh refinement methods.

8 Conclusions

A mesh refinement method has been described for solving bang-bang optimal control problems using direct collocation. First, the solution of the optimal control problem is approximated on a coarse mesh. If the approximation on the coarse mesh does not satisfy the specified accuracy tolerance, the method determines automatically if the Hamiltonian of the optimal control problem is linear with respect to any components of the control. Then, for any control component upon which the Hamiltonian depends linearly, the locations of the discontinuities in the control are obtained by estimating the roots of the switching functions associated with any component of the control that appears linearly in the Hamiltonian using estimates of the switching functions obtained using the state and costate obtained from the solution on the initial mesh. The estimates of the switching functions are then used to determine the bang-bang structure of the optimal solution. Using estimates of the locations of discontinuities in the control obtained from the detected structure, the horizon is partitioned into multiple domains and parameters corresponding to the the locations of the switch times are introduced as variables in the optimization. Then, by fixing in each domain any component of the control that has a bang-bang structure to lie at either its lower or upper limit, the multiple-domain optimal control problem is solved to accurately determine the switch times. The method has been demonstrated on three examples where it has been shown to efficiently obtain accurate approximations to solutions of bang-bang optimal control problems.

Acknowledgments

The authors gratefully acknowledge support for this research from the U.S. Office of Naval Research under grant N00014-15-1-2048 and from the U.S. National Science Foundation under grants DMS-1522629, DMS-1924762, and CMMI-1563225.