L2 Norm-Based Control Regularization for Solving Optimal Control Problems

Solutions to practical optimal control problems (OCPs) may consist of control profiles that switch between control limits or assume values interior to their admissible set, either due to activation of inequality state path constraints or existence of singular control arcs. Abrupt switches in the control (i.e., bang-bang control) jeopardizes the numerical solution of OCPs unless care is taken to isolate precise time transition points where sharp switches occur (excluding the chattering phenomenon). We propose a novel control regularization method, called Bang-Bang Singular Regularization (BBSR), based on L2 norm-based regularization. We present an analysis on the L2 norm-based regularization at two levels: 1) its connection to trigonometric regularization and 2) its ability to approximate regular and singular control arcs. The utility of the method is demonstrated in solving three classes of trajectory optimization problems: 1) space minimum-fuel low-thrust trajectories with bang-bang thrust profiles, 2) the Goddard rocket problem with its known bang-singular-bang control structure, and 3) minimum-time spacecraft reorientation with both bang-bang and second-order singular arcs. The results demonstrate the utility of the BBSR method in approximating extremal control profiles that may consist of pure regular and/or mixed regular and singular control arcs.


I. INTRODUCTION
Solutions to trajectory design and optimal control problems (OCPs) can be achieved using either direct or indirect optimization methods [1], [2].The solution procedure consists of a number of well-known computational hurdles, which when overcome, can have a significant positive impact on mission operation and vehicle design.Different OCPs can be solved using a number of software packages that have been developed to automate the solution process and enable researchers to solve practical OCPs [3], [4], [5], [6], [7].
Indirect methods have been shown to be of great importance as a fundamental formulation of optimal control necessary conditions, for example, to optimize spacecraft low-thrust trajectories [8].Specifically, utilization of calculus of variations combined with Pontryagin's Minimum Principle (PMP) leads to two-or multi-point boundary-value problems (TPBVPs/MPBVPs).Analytic solutions to the resulting The associate editor coordinating the review of this manuscript and approving it for publication was Xiwang Dong.
boundary-value problems exist only for a limited number of applications and under simplifying assumptions [9].In practice, TPBVPs are solved through a variety of numerical methods [2].For nonlinear OCPs, all known methods must come to terms, in a given application, with the issue of whether a particular extremal satisfying the local necessary conditions is, in fact, the global extremal [10], [11].Another major difficulty frequently encountered during the solution procedure is associated with the structure of the extremal controls that introduces non-smoothness into the underlying dynamics.Therefore, numerical propagation of the dynamics is plagued by the presence of control non-smoothness, which reduces the domain of convergence of the Hamiltonian boundary-value problems [12].
Despite the fact that physical systems are inherently nonlinear, in practice, a significant number of systems have control-affine dynamics.Depending on the form of the performance index, their Hamiltonian, H , may turn out to be affine in control.The mere existence of an affine control structure in the Hamiltonian does not mean that the optimal control definitely contains singular arcs [41].The strong optimality conditions, ∂H /∂u = 0 and ∂ 2 H /∂u 2 > 0, give no information on the control for the cases that the Hamiltonian is linear in the controls, whereas the weak form of the optimality conditions (i.e., PMP), in most applications, determines the control without ambiguity.These typical and ideal cases correspond to problems where ∂H /∂u = 0 never occurs during one or multiple finite time intervals.
A particular difficulty arises if the extremal trajectory entails singular arcs on finite time intervals where the so-called control switch function vanishes and when the PMP is not sufficient for characterizing an extremal control.These cases correspond to situations where the PMP necessary conditions for optimality fail to define the ''optimal'' control.The presence of singular arcs creates numerical issues for indirect solvers because the PMP, without auxiliary logic, does not provide a unique solution [42], [43].For these cases, a necessary, but not sufficient condition for existence of a singular arc is the vanishing of the switching function over a finite time interval.Without further utilization of higher-order optimality conditions a chattering phenomenon in some circumstances, is observed in the control profile [44].OCPs with mixed regular and singular control arcs are typically treated by dividing the entire trajectory into distinct phases of regular (i.e., non-singular) and singular arcs.Several approaches have been developed for solving singular OCPs, using both indirect [45], [46], [47], [48], [49], [50] and direct methods [51], [52], [53].
In general, there appears to be no justification for ignoring the possibility of the existence of singular control arcs.It cannot be stated a priori that a singular control sub-arc (even when a feasible singular control exists) will always be part of the optimal solution [54].In fact, there are examples in which a singular control is not optimal even though it is a feasible control for some choices of system parameters and boundary conditions.However, there are many examples in which singular control does appear in the optimal solution, and for this reason, the possibility of singular arcs should be considered along with the Necessary Conditions (NCs) for optimality, when we deal with a control-affine Hamiltonian.Modern optimization problems are often of the control-affine Hamiltonian type and when singular control arcs exist, they can be found for many of these problems using existing approaches as well as the methods we discuss herein.The Goddard rocket problem [55], periodic optimal flight [56], [57], and energy extraction of wave-energy converters [58], [59] are a few example OCPs to name.
Aside from resolving the ambiguity in the optimality of singular arcs and their appearance in the optimal trajectory, it is important to develop a strategy to detect the presence of singular arcs.A common practice is to perform continuation over a strictly convex (modified) Hamiltonian.A majority of these methods typically employ various regularization approaches to solve for the singular control [60], [61].A regularization method transforms the singular control problem into a series of nonsingular problems by minimizing the sum of the original objective and a regularization term, where the regularization term is typically a quadratic function of the control.This artistic approach has demonstrated capability of revealing the presence of singular arcs [62].Construction of a convex Hamiltonian can be achieved, for instance, by modifying the Lagrangian term (running cost) of standard cost functionals via the introduction of a quadratic term.This quadratic term is multiplied by a scalar continuation parameter, through which the contribution of the convex quadratic penalty term will be reduced to approach zero.This is usually achieved through a sequence of neighboring discrete choices over the admissible values of the continuation parameter.This continuation is performed specifically to gain sufficient information regarding the structure of the control, and possible existence of singular arcs, which will be used to initialize a shooting method tailored for handling singular arcs.Hybrid direct/indirect methods are also developed [63].Recently, inspired by the concept of error controls introduced in [13], the Epsilon-Trig regularization method is devised that modifies the dynamics and uses trigonometry to approximate controls; this method has been shown capable of capturing (to a high accuracy) optimal ''bang-singular-bang'' discontinuous control structures [14], [64].A recent improvement over the Epsilon-Trig method is the unified trigonometric method (UTM) [65] proposed by Mall and Taheri [65].The UTM alleviates some of the main implementation difficulties associated with the Epsilon-Trig method and makes it possible to solve a variety of OCPs with regular and singular control arcs, state-only inequality and mixed state-control inequality path constraints [66], [67].
The main contributions of the paper are as follows.First, we present an easy-to-implement control regularization to facilitate the numerical solution of OCPs with extremal 125960 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.controls that may consist of both regular (i.e., non-singular) and singular arcs or any combination of them.A prominent advantage of the proposed control regularization method is its ability to approximate both regular and singular control arcs using one unified algebraic form.The OCP NCs are fundamental as usual and the proposed method embeds the optimal control NCs into a one-parameter family that can be swept arbitrarily close to the OCP NCs.Second, the method is easy to implement in the sense that one does not need to modify the cost functional or system dynamics and the method acts as a filter that is applied directly to the control.This simplicity makes the proposed method to be on par with the hyperbolic tangent smoothing (HTS) method, which is an easy-to-implement control regularization method [39].However, the proposed BBSR method can accurately approximate a larger class of mixed regular/singular control arcs.We have shown that with our proposed control regularization method, standard boundary-value solvers (e.g., MATLAB's bpv4c) can be used for solving challenging OCPs.This makes the BBSR method accessible to practitioners from academia and industry.Third, the proposed method embeds complex MPBVPs into a neighboring family of smooth TPBVPs.This embedding is achieved through the so-called control switching function and introduction of a smoothing parameter.In essence, the complex control, which may consist of regular and singular control arcs, is embedded into a one-parameter family of smooth curves.By performing a standard continuation over the smoothing parameter, structures of control are approximated and matched to within a high accuracy, which represents an important advantage.Fourth, three important classes of trajectory optimization problems are solved from space flight to atmospheric flight domains to demonstrate the broad applicability of the method in accurate approximation of complex extremal control structures.The numerical results are compared and validated against the solutions in the literature.
The paper is organized as follows.First, a brief review of formulation of OCPs with control-affine Hamiltonian is given in Section II.This is followed by a review of the UTM method along with the PMP that characterizes the optimal control expression.A comparison between the control expressions of the UTM and the L2 norm-based regularization is presented.Section III presents an analysis on what L2 norm-based regularization achieves and its ability in approximating both bang-bang and singular control arcs.A comparison between the HTS-based and the proposed L2 norm-based regularization methods is also presented.Section IV presents the details of forming one-parameter TPBVPs associated with the three OCPs.Numerical results of the considered test problems are presented in Section V. Finally, concluding remarks are given in Section VI.

II. PROBLEM FORMULATION
The BBSR method is motivated by a key observation on the control expressions obtained with the UTM method.In this section, we first consider a simple OCP and briefly demonstrate the application of the UTM for solving OCPs with a control-affine Hamiltonian.Then, we establish a connection between the proposed L2 norm-based regularization and its practical implications.
Consider a general Bolza-form OCP written as where φ and L denote terminal and running costs, respectively, and f 0 and f 1 denote the drift and control-influence parts of a control-affine dynamical system, respectively.We assume that the running cost L depends affinely on u, i.e., L can be written as We also assume that the control u is bounded (i.e., without loss of generality, we can write |u(t)| ≤ 1).We proceed to analyze the Hamiltonian and extremal control expressions and drop the arguments of L 0 , L 1 , f 0 and f 1 for making the resulting algebraic equations more concise.Let λ(t) denote the costate associated with the state x(t).Since the Hamiltonian has a control-affine form, we have where H 0 = L 0 + λf 0 denotes the collection of terms that do not depend on control and S = L 1 +λf 1 is the so-called control switching function.The extremal (superscript '*') control can be derived according to PMP as,

A. REVIEW OF THE UTM
The key step, in the UTM [65], is to use trigonometry and to reparameterize the control, u(t) in Eq. ( 1), in terms of two new orthogonal control inputs collectively denoted as Then, the idea is to set u(t) = u 1 (t) = sin(θ(t)) while minimizing the second component of the new control, cos(θ(t)), which serves as an ''error-control'' term along the solution.Thus, the OCP defined in Eq. ( 1) can be written as, where . Note that a penalization factor, ρ, is multiplied to the orthogonal control term, cos(θ(t)).

B. APPLICATION OF THE INDIRECT METHOD
We proceed by applying the standard indirect method to the OCP in Eq. ( 4).Let λ(t) denote the costate associated with the state x(t).The explicit dependence of θ on t is dropped.
Since the Hamiltonian has a control-affine form, we have where the new control is θ.Since control is smooth and unbounded, we can use Weierstrass optimality condition and determine the extremal control [9].An alternative approach is to use the idea of Hamiltonian vectorization proposed in [68].
Let S = L 1 + λf 1 denote the coefficient of sin(θ), the Hamiltonian in Eq. ( 5) can be written as, Minimization of the Hamiltonian is achieved automatically upon using the following expressions [68], Remark: While the original control, u, is parameterized in terms of the two new orthogonal controls, after the derivation of the new extremal controls is completed and the set of NCs of optimality are established, we only use the regularized control, i.e., u = sin(θ), which means that u * (S) is replaced with u(S; ρ) = −S/ S 2 + ρ 2 .Thus, one could have directly used u(S; ρ) = −S/ S 2 + ρ 2 for the control in the formulation of the OCP.However, the steps outlined under the UTM section clarify that the L2 norm-based regularization actually tries to minimize the orthogonal (error) component of the control when control trigonometrization is used.However, the equivalence of the extremal control is established between the two approaches in this section.In other words, one could simply consider the L2 norm-based regularization to be used as a filter that directly applies to the control of the original indirect method, given in Eq. (3).Direct application of the L2 norm-based regularization to the control makes it similar to the HTS method [39] in that it only regularizes the control without affecting the standard steps taken to derive the NCs.However, as will be shown, the L2 norm-based regularization is capable of approximating both bang-bang and singular control arcs, which offers a significant advantage over the HTS method while retaining its appealing ease-of-implementation feature.

III. BANG-BANG AND SINGULAR CONTROL ARCS
Analyzing the Hamiltonian for OCPs associated with the standard indirect method (Eq.( 2)) and UTM (Eq.( 6)), we are able to perform additional analysis on the control expression, u 1 = sin(θ) that approximates the optimal control as, where ρ ∈ (0, ∞) is a smoothing parameter.Equation ( 8) approximates the optimal control in Eq. (3) because The Taylor series expansion of u(S; ρ) about ρ = 0 yields where (•) ′ = d dρ , and O(ρ 3 ) denotes the collection of higherorder terms.From Eq. ( 10) it is clear that when ρ is small and |S| is significantly greater than ρ, then, u(S; ρ) approximately equals − S |S| , which corresponds to the regular arcs of the optimal control u * (S) in Eq. ( 3).When ρ is small and |S| has the same order of magnitude as ρ and satisfies |S| ≥ ρ 2 , then, u(S; ρ) takes a value in [−1, 1], which corresponds to the singular arc of the optimal control u * in Eq. ( 3).
Therefore, it can be seen that on regular (bang-bang) arcs, u(S; ρ) converges to u * (S) as ρ decreases to 0. On singular arcs, if they exist, we expect that u(S; ρ) also converges to an optimal control, namely, to have the following condition in which u * = u * (x, λ, t) denotes the singular control.We have observed this convergence on singular arcs in our extensive numerical experiments, as is reported in Section V. Suppose after repeatedly differentiating the control switching function, ∂H ∂u = S, with respect to time t for k times, the control u explicitly appears (i.e., a finite-order singular arc), On a singular arc, S and all its time derivatives must be 0, i.e., Now suppose we are able to solve S (k) (x, λ, t, u) = 0 and obtain an expression for u on a singular arc as a function of x, λ, and t, denoted as, Under the assumption that the generalized Legendre-Clebsch second-order condition (a.k.a. the Kelley condition) is satisfied, we can substitute Eq. ( 14) into Eq.( 11) and have Based on Eq. ( 15), in a numerical smoothing procedure, we may treat as a residual and monitor ε(ρ) as ρ → 0. In Eq. ( 16), T sing denotes the time interval corresponding to a singular arc, i.e., the residual ε(ρ) measures the maximum difference between u(S; ρ) = − S √ S 2 +ρ 2 and u * (x, λ, t) on a singular arc.If the residual ε(ρ) converges to 0 as ρ → 0, then the converged control should indeed be an optimal control on both regular and singular arcs.Because T sing is typically not known a priori, in implementation of Eq. ( 16), T sing may be replaced with a conservative estimate of the singular arc time interval, Tsing , such that Tsing ⊆ T sing .Such an estimate Tsing may be obtained numerically based on S ≈ 0 (because on an exact singular arc S = 0) or based on |u(S; ρ)| ≪ 1 (because on exact regular arcs, optimal control u * is bang-bang, i.e., u * = ±1).

IV. FORMULATION OF BOUNDARY-VALUE PROBLEMS
In this section, formulations of the boundary-value problems (BVPs) associated with three OCPs are presented when the proposed BBSR method is used for solving them.

A. LOW-THRUST TRAJECTORY OPTIMIZATION
Formulation of a minimum-fuel low-thrust trajectory optimization problem is given.The spacecraft motion dynamics are expressed in terms of modified equinoctial orbital elements (MEEs) [69] defined as x = [p, f , g, h, k, L] ⊤ .MEEs are chosen due to their superior numerical stability and enhanced convergence performance in solving TPBVPs [17], [24].The state-space representation of the system is expressed as, where m is the spacecraft mass, α is the thrust steering unit vector (i.e., || α|| = 1), δ T ∈ [0, 1] is the engine throttle input, T max is the maximum available thrust, and c = I sp g 0 is the engine exhaust velocity with I sp denoting the constant specific impulse value.The matrix mapping the control input to the states A(x MEE ) ∈ R 6×3 and the vector representing the unperturbed motion b(x MEE ) ∈ R 6×1 (each nonlinearly dependent on the MEEs) are given as, where s L = sin(L), c L = cos(L), s 2 = 1 + h 2 + k 2 , and w = 1 + fc L + gs L .Also, arguments of A and b are dropped for brevity in the remainder of the paper.
We consider fixed-time, rendezvous-type, minimum-fuel trajectories.The cost functional is written in Mayer form as minimize The (optimal control) Hamiltonian can be formed as where ⊤ is the costate vector.To solve for the optimal controls α * and δ * T , one has to use Pontryagin's minimum principle.In addition, following Lawden's primer vector theory, the optimal thrust steering unit vector, α * , and optimal thrust throttle input, δ * T , can be derived [17] as, where S T denotes the thrust switching function that governs the thrust throttle value δ * T .Note that there exists a possibility for singular control arcs (i.e., when S T = 0 for one or multiple non-zero finite time intervals), however, singular control arcs are rare in space flights [65].Because of the bounded discontinuous behavior of 0 ≤ δ * T ≤ 1, we will employ the BBSR method (Eq.( 17)) to smoothly embed the non-smooth throttle control into the continuous-time dynamics as, where the smoothing parameter, ρ ∈ (0, ∞), is introduced whose value determines the ''smoothness'' of the throttle profile.Equation ( 23) can simply be substituted into the Hamiltonian.The costate dynamics can be derived using the Euler-Lagrange equation as, with H MF defined in Eq. ( 21).The algebraic relations of the costates differential equations are not presented for conciseness, but we used MATLAB's symbolic tool to derive costate equations.The initial and final position vectors, initial and final velocity vectors, and initial mass are all known, but final mass is unknown.Thus, a condition on the final mass costate can be derived through the transversality condition.The seven initial costates, λ(t 0 ), are still unknown and must be determined to satisfy the final constraints, where x MEE T denotes the target MEEs.In summary, the TPBVP consists of state dynamics given in Eq. ( 18), costate equations, Eq. ( 24), the steering α * given in Eq. ( 22), smooth throttle in Eq. ( 23), and final constraints given in Eq. ( 25).There are various numerical methods (e.g., single-or multiple-shooting schemes) to solve the resulting one-parameter family of TPBVPs.In a single-shooting scheme and for a relatively large initial value of the smoothing parameter, ρ, the unknown vector, λ(t 0 ), has to be determined that obeys the NCs.Then, the converged solution is used in a recursive manner to solve another TPBVP with a smaller value of ρ.A numerical continuation on the value of ρ can be performed to reduce its value to ρ ≈ 1.0 × 10 −5 .It is shown that difference in the value of the cost functional of the exact solution and the regularized solutions with small values of ρ is insignificant for engineering applications [39].

B. GODDARD ROCKET PROBLEM
Consider the vertical ascent of a variable-mass rocket under the inverse-square gravitational acceleration.Let x(t) = [h, v, m] ⊤ , where h, v and m denote the altitude, velocity, and mass, respectively.The non-dimensional equations of motion can be expressed as [70]: where V e is a constant exhaust velocity, and T (u) is thrust of the rocket, D(x) is drag, and q(x) is dynamic pressure with their relations given as, where T max denotes the maximum thrust value and u ∈ [−1, 1], which corresponds to T ∈ [0, T max ].Drag force, D(x), is a function of the drag coefficient, C D , reference area, A, and dynamic pressure, q(x); all normalized by the maximum gravitational force expressed in terms of the initial mass of the rocket, m 0 , and acceleration due to sea-level gravity constant, g 0 .An exponential atmospheric model is used to model the dynamic pressure, q(x), with the change in h and v, with respect to the air density at sea level, ρ 0 , and initial height, h 0 , from the Earth's surface.The air density drop with altitude is given as a decay rate, β.The dynamical model and values were adopted from the work of Graichen and Petit [71], originally given by Seywald and Cliff [70].
The scaling of parameters can be performed according to Ref. [72].For the sake of comparison, the following values have been used in this work: β = 500, V e = 0.5, The Goddard rocket problem is a free-final-time OCP, with a fixed initial state vector, x(0) = [h 0 , v 0 , m 0 ] ⊤ , fixed initial time, t 0 = 0, and fixed final mass, m f = 0.6 (all parameters, including time, are non-dimensional).The OCP cost functional can be stated in Mayer form as, minimize where The state and time boundary conditions are given as: denote the costate vector.The (optimal control) Hamiltonian can be formed as, The costate dynamics can be obtained by using the Euler-Lagrange equation, λ = −[∂H G /∂x] ⊤ .The Hamiltonian and the control switching function can be written as, 125964 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Thus, the extremal control is characterized as, The extremal control can be directly approximated as, Since the Hamiltonian does not explicitly depend on time and the final time, t f , is free, the Hamiltonian value has to remain zero along any extremal solution.Since h(t f ) and v(t f ) are free, we obtain that λ h (t f ) = −1 and λ v (t f ) = 0 due to transversality conditions.The TPBVP consists of the state dynamics, Eq. ( 26), costate dynamics, λ = −[∂H G /∂x] ⊤ , the smooth approximation of the extremal control, Eq. ( 32), the condition on the final value of the Hamiltonian, H G (t f ) = 0, and the set of final conditions written as, which can be solved using a single-shooting or collocationbased methods.For a single-shooting scheme, a ''good'' estimate for the missing initial costates, λ(t 0 ), and final time, t f , have to be provided to the solver.These solvers iterate over the unknown values, using Newton-or Quasi-Newton based methods, to satisfy the residual vector in Eq. ( 33) to within a specified user-defined tolerance.

C. MINIMUM-TIME REORIENTATION OF SPACECRAFT
Motivated by the ongoing interest in the optimal control of rigid bodies [73], [74], we revisit the minimum-time reorientation of an axisymmetric rigid spacecraft under the influence of three control torques.We consider an axisymmetric rigid-body nonrest-to-rest maneuver, which can have complex regular and singular control profiles.
Let u I = [u , u 2 u 3 ] ⊤ denote the control vector and let x I = [x 1 , x 2 , ω 1 , ω 2 , ω 3 ] ⊤ denote the state vector.Time rate of change of the states can be written as, where ⊤ denotes the angular velocity vector.Here, x 1 and x 2 are used for determining the relative position of the inertially fixed ''3''-axis as viewed by an observer fixed in the body frame of the spacecraft [75], [76].In Eq. ( 34), a = (I 2 − I 3 )/I 1 where (I 1 , I 2 , I 3 ) are the principal-axis moments of inertia (with a = 0.5).The derivation details are given in [74].
The OCP that corresponds to the minimum-time reorientation of an axisymmetric rigid spacecraft can be written as, minimize where U = {u I | |u i | ≤ 1, for i ∈ {1, 2, 3}} denotes the admissible control set.The boundary conditions are summarized as: and t f is free.We will consider two cases.In Case I, we include ω 3 and u 3 , but don't enforce ω 3 (t) = 0 constraint for t ∈ [t 0 , t f ].In Case II, we enforce ω 3 (t) = 0 (corresponding to u 3 (t) = 0).This corresponds to removing ω3 from the set of differential equations and its boundary conditions.The TPBVPs associated with both cases are presented and numerical results for both cases are given in Section V.
Case I (ω ⊤ denote the costate vector.Upon forming the Hamiltonian, H I = λ ⊤ I f I , the costate dynamics can be obtained by using the Euler-Lagrange equation, λI = −[∂H I /∂x I ] ⊤ .The extremal controls, u i for i ∈ {1, 2, 3}, are characterized as, The non-smooth controls are approximated using the BBSR method as, Since the Hamiltonian does not explicitly depend on time and the final time, t f , is free, we can obtain a condition on the final value of the Hamiltonian as H I (t f ) = −1.The TPBVP consists of the state dynamics, Eq. ( 34), costate dynamics, λI = −[∂H I /∂x I ] ⊤ , the extremal control, Eq. ( 37), the condition on the final value of the Hamiltonian, H I (t f ) = −1, and the set of final conditions written as, Case II (ω The time rate of change of the states can be written as, where ω = [ω 1 , ω 2 ] ⊤ denotes the angular velocity vector. Let The Hamiltonian can be formed as The boundary conditions are similar to Case I except for the removal of the conditions on ω 3 .The non-smooth controls are approximated using the BBSR method (with S i defined in Eq. ( 36)) as, The TPBVP consists of the state dynamics, Eq. ( 39), costate dynamics, λII = −[∂H II /∂x II ] ⊤ , the extremal control, Eq. ( 40), the condition on the final value of the Hamiltonian, H II (t f ) = −1, and the set of final conditions written as,

V. NUMERICAL RESULTS
The numerical results for the three test problems are presented in this section.All problems are solved in MATLAB running on a Windows 10 Enterprise OS with 32 GB of RAM and Interl(R) Xeon(R) Gold CPU@ 2.3GHz.

A. MINIMUM-FUEL LOW-THRUST TRAJECTORY
We consider a benchmark minimum-fuel low-thrust trajectory optimization from the Earth to asteroid Dionysus (referred to as the E2D problem) [17].This asteroid has high eccentricity and inclination values of 0.542 and 13.54 deg, respectively.The boundary conditions and parameters are chosen to compare the results with those reported in [17].The Figure 2 shows the three-dimensional minimum-fuel trajectory.Figure 3 shows the time histories of the engine throttle and thrust switching function with ρ = 1.0 × 10 −5 .The optimal throttle profile consists of 12 switches, with late-departure and early-arrival coast arcs.The final mass is m(t f ) = 2718.32kg, which is the global optimal solution for this benchmark problem [39].We performed 50  Here, r(t f ) and v(t f ) denote the position and velocity at the end of the low-thrust trajectory and r D and v D denote the position and velocity vectors of Dionysus [24].The results are consistent with the solutions obtained and reported in the literature using the HTS method, in which the only difference is in the problem formulation for throttle regularization, which can be written as .

B. THE GODDARD ROCKET PROBLEM
The TPBVPs associated with the Goddard rocket problem are solved in MATLAB using its built-in boundary-value problem solver, bvp4c.The RelTol and AbsTol fields in options are set to 1.0 × 10 −4 , respectively.The bvp4c is initialized with all costate values set to 0 and the guess for t f is 0.5.The initial value of the smoothing parameter is 1.0 and its value is decreased by multiplying it with a 0.85 factor.The continuation is performed until ρ ≤ 1.0 × 10 −4 .Thus, a set of 57 TPBVPs are solved.To improve the accuracy of the solution, the last solution of the bvp4c is resolved using bvp5c with RelTol and AbsTol fields in options are set to 1.0 × 10 −7 .It takes nearly 3 seconds to obtain the solutions.To demonstrate the ability of the BBSR method in approximating regular and singular control arcs, a finite set of T max ∈ {1.5, 2.5, 3.5, 4.5} is considered.Figure 4 depicts the time histories of the states and control.Figure 5 shows the time histories of the switching function  16)) versus ρ.
for each T max value.As it is known, the control profile will change and for different values of T max , it is possible that a singular control arc becomes part of the extremal control.For the chosen problem parameters, and for T max ≥ 2.2, singular control becomes part of the extremal control.Also, as the value of T max increases beyond the critical T max = 2.2, the duration of the singular control arc increases.The results illustrate the ability of the proposed BBSR method in approximating regular and singular control arcs using the L2 norm-based control regularization given in Eq. (32).Table 1 summarizes the values of t f and h(t f ) for different T max values.For the case of T max = 3.5, the results are consistent with those reported in [71] for the same thrust value.As the value of T max increases, the flight time gets shorter while the maximum altitude value increases.
In order to verify Eq. ( 16), we derived the expression for the singular control arc for the Goddard rocket problem.It is known that Goddard rocket problem has a finite-order singular arc.The order of the singular control is q = 1, i.e., d 2q dt 2q S G (t) = 0 is an explicit function of control, which we denote as u sing .The algebraic expression of the singular control is lengthy, but we have evaluated it along the solution that is obtained using the BBSR method.Figure 6 shows the time histories of T , S G and T max u sing for t ∈ [0, 0.1].It is clear that during the time that S G (t) = 0, the control lies on the singular control arc. Figure 7 shows the profile of the control error (defined in Eq. ( 16)) versus ρ values.As the value of ρ decreases, the error becomes smaller.

C. MINIMUM-TIME REORIENTATION OF SPACECRAFT
Case I: The resulting TPBVPs are solved in MATLAB using its built-in boundary-value problem solver, bvp4c.The RelTol and AbsTol fields in the options are set to 1.0 × 10 −2 , respectively.The bvp4c is initialized with random values for all costate in range of [0,1] and the guess for t f is 5.The initial value of the smoothing parameter is 1.0 and its value is decreased by multiplying it with a 0.95 factor.The continuation is performed until ρ ≤ 1.0 × 10 −6 .Thus, a set of 270 TPBVPs are solved.To improve the accuracy of the solution, the last solution of the bvp4c is resolved using bvp5c with RelTol and AbsTol fields in the options set to 1.0 × 10 −12 .It takes approximately 6 seconds to obtain the solutions.The shortest maneuver time is t f = 2.479.Figure 8 shows the time histories of the states.Figure 9 shows the time histories of the three controls and their corresponding switching functions.The extremal solution consists of regular control arcs with two control switches in u 1 , one control switch in u 2 and two control switches in u 3 .Figure 10 shows the time history of the Hamiltonian, which verifies the that the solution is an extremal solution.The spikes in the plot correspond to the control switches that affect the value of the Hamiltonian in the sixth significant digit.
Case II: We adopted a solution strategy similar to Case I.However, it is theoretically proven that the solution to Case II (i.e., when ω 3 (t) = 0) will consist of a singular control arc for u * 1 while u * 2 will have a regular control profile.We will present the results for two small values of the smoothing parameter.Figure 11 shows the time histories of the states.Figure 12 shows the time histories of the two controls and their corresponding switching functions.The extremal solution consists of mixed regular and singular control arcs in u 1 and regular control arcs in u 2 with one switch.The minimum maneuver time is t f = 2.883892.
Figures 13 and 14 show the time histories of u 1 and S 1 for ρ = 9.56 × 10 −9 and ρ = 9.73 × 10 −15 , respectively.The numerical value of the maneuver time is still the same for both of the smoothing values, t f = 2.883892.However, the transition from a regular control into a second-order singular arc is known to lead to control chattering [2].In particular, compared to the u 1 profile in Figure 12, the control has two more control switches when u 2 becomes -1 prior to the singular control.However, Figure 14 and its enlarged subplot exhibit the control chattering when the value of ρ is reduced even further.

VI. CONCLUSION
A method is proposed for a novel unified control regularization of regular (i.e., bang-bang) and singular control arcs that arise in solving optimal control problems.The proposed method called, Bang-Bang Singular Regularization (BBSR), is based on a key observation in the structure of the control when the original control is parameterized in terms of two orthogonal trigonometric-based components.An equivalence is established between the trigonometric-based regularization and the L2 norm-based regularization.Leveraging the existing equivalence, the easy-to-implement BBSR control regularization is applied directly to the control of the original non-smooth optimal control problem.In other words, using the algebraic expression for the control switching function, an L2 norm-based control regularization is used to construct a one-parameter family of smooth controls that approximate both regular and singular control arcs in a compact unified manner.
Application of the BBSR method is demonstrated successfully by solving three important classes of optimal control problems: 1) low-thrust trajectory design, 2) the Goddard rocket problem with a known ''bang-singularbang'' control structure, and 3) minimum-time nonrest-torest attitude reorientation.Numerical results demonstrated the flexibility of the method in generating extremal control structures that consist of mixed regular (i.e., bang-bang) and singular control arcs.

TABLE 1 .
Summary of results for the Goddard rocket problem.FIGURE 8. Case I: spacecraft reorientation problem: time histories of the states with ρ ≈ 9.66 × 10 −7 .