# A High-Precision Single Shooting Method for Solving Hypersensitive Optimal Control Problems

Authors: Binfeng Pan; Yang Wang; Shaohua Tian
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Solving hypersensitive optimal control problems is a long-standing challenge for decades in optimization engineering, mainly due to the possible nonexistence of the optimal solution to meet the required error tolerance under double-precision arithmetic and the hypersensitivity of the optimal solution with respect to the initial conditions. In this paper, a new high-precision single shooting method is presented to address the above two difficulties. Multiple-precision arithmetic and Taylor series method are introduced to provide the accurate optimal solution with arbitrary higher significant digits and arbitrary higher integral accuracy, respectively. Besides, a new modified bidirectional single shooting method is developed, which fully utilizes the three-segment structure of the hypersensitive optimal control problems and provides appropriate initial guess that is close to the optimal solutions. Numerical demonstrations in a typical hypersensitive optimal control problem are presented to illustrate the effectiveness of this new method, which indicates that the accurate optimal solution of this challenging problem can be easily solved by this simple single shooting method within several iterations.


---

## Body

## 1. Introduction

An optimal control problem and its associated Hamiltonian boundary value problem (HBVP) are called hypersensitive if the time interval of interest is long relative to the rates of expansion and contraction of the linearized Hamiltonian dynamics in certain directions in a neighborhood of the optimal solution. Solving hypersensitive optimal control problems (HOCPs) is a challenging task, which has been studied extensively for decades [1–6].Typically, there are two types of methods, usually categorized as direct and indirect methods, for solving HOCPs. The direct method converts the optimal control problem into a nonlinear programming problem (NLP) by appropriate discretization [7]. Various direct methods have been developed during the past few decades, such as direct multiple shooting method [8], gradient-based optimization techniques [9, 10], differential inclusion method [11], interior point method [7, 12], collocation method [13, 14], pseudospectral method [4, 15–17], and control parameterization method [18–20]. Direct methods are straightforward and robust to accommodate complex conditions, such as path constraints and intermediate constraints. When applying direct methods to solve HOCPs, the key is to use a higher density of nodes in the first and third segment, according to the three-segment structure characteristic of HOCPs, which is described qualitatively as “take-off,” “cruise,” and “landing” by analogy to an airport-to-airport trajectory for an aircraft [1]. Reference [3] introduced a simple method for mesh point distribution for solving HOCPs based on density functions, and the problem of mesh refinement is subsequently converted to a problem of finding an appropriate density function. Reference [4] presented a hp-adaptive pseudospectral method to solve HOCPs, which iteratively determines the number of segments, the width of each segment, and the polynomial degree required in each segment. Reference [5] proposed a symplectic algorithm with nonuniform grids combining with density functions to solve HOCPs. Later, [21] applied GPOPS-II, a commercial software program that employs a Legendre-Gauss-Radau quadrature orthogonal collocation method to solve HOCPs. In their method, an adaptive mesh refinement method is implemented to determine the number of mesh intervals and the degree of the approximating polynomial within each mesh interval to achieve a specified accuracy. However, there are several disadvantages for solving HOCPs by direct methods. The first one is that the obtained optimal solution by direct methods is merely approximate to the accurate solution. The error between these two solutions is usually small and acceptable for most optimal control problems; however, this error may be extremely large and unacceptable for HOCPs, as demonstrated in this paper later. Besides, direct methods cannot provide insight regarding the optimal solutions [6].Unlike direct methods, indirect methods transform the original problems to HBVPs by first-order optimality conditions [22], which are helpful to provide the dynamical information of the systems. Besides, the optimal solution is guaranteed to be a least extremal. However, it is extremely difficult to solve HOCPs by indirect methods, mainly due to the high sensitivity to the initial unknowns. Based on the three-segment structure, [1, 2, 23, 24] developed a dichotomic basis method to solve HOCPs. The core idea is the decomposition of contracting and expanding behaviors of the system by coordinate transformation, and the iterative process is required due to the approximation of dichotomic basis. Recently, [6] presented a finite-time Lyapunov analysis based method to construct an approximate dichotomic basis and developed a corresponding manifold-following solution approximation method. However, the optimal solutions obtained by the above methods are only near-optimal, and, as pointed in [6], the sensitivity of the optimal solution to the errors in the unknown boundary conditions may cause ill-conditioning, which potentially exceeds the available precision. Thus, to the best of the authors’ knowledge, the above difficulties are not well solved in the existing literature, and obtaining the accurate optimal solutions of HOCPs by indirect methods is still extremely daunting.In this paper, a new high-precision single shooting method is presented to address the above two main difficulties for solving HOCPs by indirect methods:(1) the possible nonexistence of the optimal solution to meet the required error tolerance under double-precision arithmetic; (2) the hypersensitivity of the optimal solution with respect to the initial conditions. Multiple-precision arithmetic and Taylor series method, which provide higher significant digits and higher integral accuracy, respectively, are introduced to address the first difficulty. A modified bidirectional single shooting method (MBISSM) is developed to settle the second one, which introduces an intermediate point along the equilibrium segment as the starting point and shoots to both boundary points simultaneously. Finally, numerical simulations are provided to demonstrate the effectiveness and efficiency of the proposed method, which reveals that the challenging HOCPs can be easily solved within several iterations by this proposed method.


## 2. Hypersensitive Optimal Control Problem

Consider the optimal control problem with the dynamical equations as follows:(1)x˙t=fxt,ut,where x(t)∈Rn and u(t)∈Rm are the state and control vector at time t, respectively, and f:Rn×Rm→Rn is the smooth function. The performance index is(2)J=∫t0tfLxt,utdt,where t0 and tf are the prescribed initial time and terminal time, respectively. The boundary conditions are defined as(3)xt0=x0,(4)xtf=xf,where x0 and xf are known.According to the optimal control theory [22], the corresponding Hamiltonian function is implicitly dependent on time, which is expressed as(5)Hxt,ut,λt=Lxt,ut+λtTfxt,ut,where λ(t)∈Rn is the costate vector. The Hamiltonian differential equations are as follows:(6)x˙t=∂Hxt,ut,λt∂λt,λ˙t=-∂Hxt,ut,λt∂xtand the optimal control vector u∗(t) satisfies (7)∂Hxt,u∗t,λt∂u∗t=0.Denotep=(x,λ), p˙=G(p) as an alternate expression for the Hamiltonian system in (6). For sufficiently large value of tf, that is, tf is long relative to the rates of both expansion and contraction in certain directions in the neighborhood of the optimal solution; this problem becomes completely hypersensitive. It is a degenerate class of two time-scale HBVPs that is composed of fast and slow subsystems. Otherwise, the problem is partial hypersensitive [1, 2]. In this paper, the completely hypersensitive case is considered, and the saddle type equilibrium point peq=(xeq,λeq) is supposed to exist. As illustrated in [1, 2], the optimal trajectory shows the three-segment structure geometrically, which is described by (8)p∗t=pi∗t,0≤t≤tibl,pe∗t,tibl≤t≤tf-tfbl,pt∗t,tf-tfbl≤t≤tf,where tibl and tfbl denote the durations of the initial and terminal segment, respectively. pi∗(t) is the optimal solution of the initial segment with initial conditions pi(0)=(x0,λ0), pt∗(t) is the optimal solution of the terminal segment with terminal conditions pt(tf)=(xf,λf), and pe∗(t) is the optimal solution of equilibrium segment with pe∗(t)≈peq. Please notice that tibl and tfbl should be large enough to allow the boundary-layer segments to reach peq with sufficient accuracy [2].Based on the observation thatpe∗(t)≈peq, [1, 23] presented a dichotomic basis method for HOCPs. In their method, the original optimal control problem is spitted into three subproblems.(a) Subproblem A. Solve Hamiltonian system p˙=G(p) with the optimal control determined by (7), initial conditions defined in (3), and terminal conditions given as (9)pt¯ibl=peq,where t¯ibl≥tibl is the guessing time duration of the initial segment.(b) Subproblem B. Solve Hamiltonian system p˙=G(p) without control effort, that is,(10)u=0together with the initial conditions and terminal conditions as follows:(11)pt¯ibl=peq,ptf-t¯fbl=peq,where t¯fbl≥tfbl is the time duration of the terminal segment. Since the control variables are totally determined, this subproblem is actually an initial value problem solved by a numerical integration method.(c) Subproblem C. Solve Hamiltonian system p˙=G(p) with the optimal control determined by (7), the initial conditions defined in (4), and terminal conditions given as (12)ptf-t¯fbl=peq.Please notice that subproblem C has the same structure as subproblem A, except that the time direction is opposite. Thus it requires a backward integration to complete numerical computation.Dichotomic basis method is elegant and efficient; however, the obtained optimal solution by this method is only near-optimal. Thus, obtaining the accurate optimal solution of HOCPs by indirect methods is still an open challenge.


## 3. High-Precision Single Shooting Method

In this section, a new high-precision single shooting method is developed to obtain the accurate optimal solution for HOCPs, which combines with three important themes: multiple-precision arithmetic, Taylor series method, and a MBISSM. In this proposed method, multiple-precision arithmetic and Taylor series method provide the ability of high-precision computation with higher significant digits and higher integral accuracy, respectively, and the MBISSM provides appropriate initial guess that is in the close neighborhood of the optimal solution.


### 3.1. Multiple-Precision Arithmetic

Currently, most optimization applications and algorithms utilize double-precision (64-bit) accuracy, which allows numerical computations to be accurate to ~15 significant digits and is more than sufficient for most optimal control problems. However, for a rapidly growing body of applications, such as HOCPs considered in this paper, double-precision arithmetic may not be sufficient any more. Therefore, a higher level of numeric precision is required, and multiple-precision arithmetic is an effective way.Typically, the ability of multiple-precision arithmetic is commonly provided in software, as hardware support for higher precision is severely lacking. It is highly likely that a commodity processor supporting higher floating-point precision will be developed in the near future [25]. Quadruple (128-bit or ~30 significant digits) or octal (256-bit or ~60 significant digits) precision arithmetic is a typical choice. For some specific problems, hundreds or even more significant digits are required to obtain numerically meaningful results.Nowadays there are several commercial or freely available multiple-precision software packages, which usually provide high-level language interfaces, custom data types, and operator overload features. In this paper, a Multiprecision Computing Toolbox developed by Advanpix [26] is applied for high-precision computation. This software is a MATLAB extension for computing with arbitrary precision; that is, the toolbox equips MATLAB with a new multiple-precision numeric type and extensive set of mathematical functions that are capable of computing with arbitrary precision. Thus many existing MATLAB programs can be directly converted to achieve arbitrary precision with small efforts.


### 3.2. Taylor Series Method

Numerical methods for integrating ordinary differential equations have been studied since the end of the last century, and a large number of integration formulas have been developed. Typically, most numerical integration methods are designed to approximate the exact solution of the ordinary differential equations, which gives rise to truncation errors. Generally, the truncation errors remain to be sufficiently small for most optimal control problems with common numerical integration methods, such as the well-known Runge-Kutta methods. However, for HOCPs considered in this paper, the truncation errors would be amplified to be incredibly large with these methods, even when multiple-precision arithmetic is applied.The truncation errors can be reduced by two different ways: by reducing the step size and by using the higher-order integration formula. If the step size between two adjacent values becomes smaller, the truncation errors of the numerical integrations would decay. The disadvantage is that the step size may decrease to be so small that it may take extremely long time to complete the integration, and the round-off error would be increased [27]. Thus, a higher-order integration formula, called Taylor series method [28], is preferred to provide accurate numerical integration solution. Taylor series method is one of the oldest methods, which traces back to Newton and Euler. It has an advantage that its formula at an arbitrarily high order can be easily expressed in the united form with large step sizes and thus small values of computing time is required [29]. Therefore, from the viewpoint of numerical simulations, it is rather easy and attractive to use this method at a very high order so as to reduce the truncation errors to a required level.Consider the initial value problem as follows:(13)y˙t=gt,yt,yt0=y0,where y(t)∈Rn, t0 is the specified initial time, y0 is the prescribed initial conditions, and g[t,y(t)] is the smooth function. According to the rules of Taylor series method, the value of the solution y(ti+1) at ti+1=ti+hi+1 is approximated by the nth degree Taylor series of y(t) developed at ti and evaluated at t=ti+1, that is,(14)yti+1≃yti+dytidthi+1+⋯+1n!dnytidtnhi+1n≃yi+gti,yihi+1+⋯+1n!dn-1gti,yidtn-1hi+1n=defyi+1,where hi+1 is the step size and y(t0) is defined as(15)yt0=y0.In this paper, a variable step size and variable order (VSVO) scheme, which aims to control the integration error automatically by adjusting the integration step size through the criteria related to the prescribed tolerance and variable order formulation, is implemented during the integration by Taylor series method. The details of VSVO scheme can be found in [30].


### 3.3. Modified Bidirectional Single Shooting Method

Single shooting methods are one of the most commonly used numerical methods to solve HBVPs, which attempt to identify appropriate initial conditions for a related initial value problem to provide the solution of the original boundary value problem. If the boundary conditions are not satisfied to the desired accuracy, the selecting process of initial conditions is repeated with a new set of initial conditions until the desired accuracy is achieved. Single shooting methods are typically classified as forward single shooting method (FSSM), backward single shooting method (BSSM), and bidirectional single shooting method (BISSM), according to their different integration directions. As illustrated in Figure1, FSSM and BSSM, which start from the initial and terminal boundary points, respectively, are unidirectional. Their solving procedures are all the same except for the opposite integral direction. Other than the above unidirectional shooting methods, BISSM starts from both boundary points and shoots to an arbitrary intermediate point [31]. However, these methods are not appropriate to solve HOCPs, mainly due to the extreme difficulty of selecting proper initial guess, even if multiple-precision arithmetic and Taylor series method are incorporated.Figure 1
Typical single shooting methods.
(a)
FSSM (b)
BSSM (c)
BISSMIn this paper, a new MBISSM is presented to address the above difficulty. As illustrated in Figure2, this MBISSM starts from an arbitrary intermediate point inside the time interval and shoots to both boundary points simultaneously. Denote an intermediate point p(ti) at ti∈(t0,tf) along the equilibrium segment as the starting point, which is the unknown to be determined numerically by the shooting method. According to the observation that this optimal solution is near the equilibrium value peq during equilibrium segment as suggested by three-segment structure, it is reasonable to set peq as the initial guess, which is expressed as follows:(16)pguessti=peq.Figure 2
Modified bidirectional single shooting method (MBISSM).In this MBISSM, the boundary conditions are given as(17)xt0=x0,xtf=xfwhich are both required to be satisfied by backward and forward integrations from the starting point, respectively. Denote p(t0∣p(ti))={x(t0∣p(ti)),λ(t0∣p(ti))} and p(tf∣p(ti))={x(tf∣p(ti)),λ(tf∣p(ti))} to be the state and costate vectors at t=t0 and t=tf by integrating p(ti) backward and forward, respectively. Thus, the optimization problem is stated as follows:(18)Findptisuchthatxt0∣pti=x0,xtf∣pti=xf.Once the optimal solution at the starting point is obtained by the MBISSM, denoted asp∗(ti), the states at the initial and terminal time are completely determined by integrating from p∗(ti), denoted as x∗(t0∣p∗(ti)) and x∗(tf∣p∗(ti)), respectively. The initial state errors, which are the differences between x∗(t0∣p∗(ti)) and the given initial value x0 defined in (3), can be expressed as(19)Δxt0∣p∗ti=x∗t0∣p∗ti-x0.Denotex(tf∣x0,λ0∗) as the state that is integrated from x0 and λ0∗=λ∗(t0∣p∗(ti)). The terminal state errors are denoted as(20)Δxtf∣x0,λ0∗=xtf∣x0,λ0∗-xf.Due to the existence of the initial state errors defined in (19), the Euclidean norm of the terminal state errors may be extremely larger for the HOCPs, which may violate the error tolerance. In this case, a lower error tolerance for solving (18) is required with the help of multiple-precision arithmetic, until the required accuracy is satisfied.Although multiple shooting methods are wildly demonstrated to be superior to the single shooting methods in terms of the convergence domain [32], they do not provide any information regarding the flow structure in the phase space, which aids the development of a simpler approximate solution method for solving HOCPs [6]. Thus, multiple shooting methods are not considered in this paper.


## 3.1. Multiple-Precision Arithmetic

Currently, most optimization applications and algorithms utilize double-precision (64-bit) accuracy, which allows numerical computations to be accurate to ~15 significant digits and is more than sufficient for most optimal control problems. However, for a rapidly growing body of applications, such as HOCPs considered in this paper, double-precision arithmetic may not be sufficient any more. Therefore, a higher level of numeric precision is required, and multiple-precision arithmetic is an effective way.Typically, the ability of multiple-precision arithmetic is commonly provided in software, as hardware support for higher precision is severely lacking. It is highly likely that a commodity processor supporting higher floating-point precision will be developed in the near future [25]. Quadruple (128-bit or ~30 significant digits) or octal (256-bit or ~60 significant digits) precision arithmetic is a typical choice. For some specific problems, hundreds or even more significant digits are required to obtain numerically meaningful results.Nowadays there are several commercial or freely available multiple-precision software packages, which usually provide high-level language interfaces, custom data types, and operator overload features. In this paper, a Multiprecision Computing Toolbox developed by Advanpix [26] is applied for high-precision computation. This software is a MATLAB extension for computing with arbitrary precision; that is, the toolbox equips MATLAB with a new multiple-precision numeric type and extensive set of mathematical functions that are capable of computing with arbitrary precision. Thus many existing MATLAB programs can be directly converted to achieve arbitrary precision with small efforts.


## 3.2. Taylor Series Method

Numerical methods for integrating ordinary differential equations have been studied since the end of the last century, and a large number of integration formulas have been developed. Typically, most numerical integration methods are designed to approximate the exact solution of the ordinary differential equations, which gives rise to truncation errors. Generally, the truncation errors remain to be sufficiently small for most optimal control problems with common numerical integration methods, such as the well-known Runge-Kutta methods. However, for HOCPs considered in this paper, the truncation errors would be amplified to be incredibly large with these methods, even when multiple-precision arithmetic is applied.The truncation errors can be reduced by two different ways: by reducing the step size and by using the higher-order integration formula. If the step size between two adjacent values becomes smaller, the truncation errors of the numerical integrations would decay. The disadvantage is that the step size may decrease to be so small that it may take extremely long time to complete the integration, and the round-off error would be increased [27]. Thus, a higher-order integration formula, called Taylor series method [28], is preferred to provide accurate numerical integration solution. Taylor series method is one of the oldest methods, which traces back to Newton and Euler. It has an advantage that its formula at an arbitrarily high order can be easily expressed in the united form with large step sizes and thus small values of computing time is required [29]. Therefore, from the viewpoint of numerical simulations, it is rather easy and attractive to use this method at a very high order so as to reduce the truncation errors to a required level.Consider the initial value problem as follows:(13)y˙t=gt,yt,yt0=y0,where y(t)∈Rn, t0 is the specified initial time, y0 is the prescribed initial conditions, and g[t,y(t)] is the smooth function. According to the rules of Taylor series method, the value of the solution y(ti+1) at ti+1=ti+hi+1 is approximated by the nth degree Taylor series of y(t) developed at ti and evaluated at t=ti+1, that is,(14)yti+1≃yti+dytidthi+1+⋯+1n!dnytidtnhi+1n≃yi+gti,yihi+1+⋯+1n!dn-1gti,yidtn-1hi+1n=defyi+1,where hi+1 is the step size and y(t0) is defined as(15)yt0=y0.In this paper, a variable step size and variable order (VSVO) scheme, which aims to control the integration error automatically by adjusting the integration step size through the criteria related to the prescribed tolerance and variable order formulation, is implemented during the integration by Taylor series method. The details of VSVO scheme can be found in [30].


## 3.3. Modified Bidirectional Single Shooting Method

Single shooting methods are one of the most commonly used numerical methods to solve HBVPs, which attempt to identify appropriate initial conditions for a related initial value problem to provide the solution of the original boundary value problem. If the boundary conditions are not satisfied to the desired accuracy, the selecting process of initial conditions is repeated with a new set of initial conditions until the desired accuracy is achieved. Single shooting methods are typically classified as forward single shooting method (FSSM), backward single shooting method (BSSM), and bidirectional single shooting method (BISSM), according to their different integration directions. As illustrated in Figure1, FSSM and BSSM, which start from the initial and terminal boundary points, respectively, are unidirectional. Their solving procedures are all the same except for the opposite integral direction. Other than the above unidirectional shooting methods, BISSM starts from both boundary points and shoots to an arbitrary intermediate point [31]. However, these methods are not appropriate to solve HOCPs, mainly due to the extreme difficulty of selecting proper initial guess, even if multiple-precision arithmetic and Taylor series method are incorporated.Figure 1
Typical single shooting methods.
(a)
FSSM (b)
BSSM (c)
BISSMIn this paper, a new MBISSM is presented to address the above difficulty. As illustrated in Figure2, this MBISSM starts from an arbitrary intermediate point inside the time interval and shoots to both boundary points simultaneously. Denote an intermediate point p(ti) at ti∈(t0,tf) along the equilibrium segment as the starting point, which is the unknown to be determined numerically by the shooting method. According to the observation that this optimal solution is near the equilibrium value peq during equilibrium segment as suggested by three-segment structure, it is reasonable to set peq as the initial guess, which is expressed as follows:(16)pguessti=peq.Figure 2
Modified bidirectional single shooting method (MBISSM).In this MBISSM, the boundary conditions are given as(17)xt0=x0,xtf=xfwhich are both required to be satisfied by backward and forward integrations from the starting point, respectively. Denote p(t0∣p(ti))={x(t0∣p(ti)),λ(t0∣p(ti))} and p(tf∣p(ti))={x(tf∣p(ti)),λ(tf∣p(ti))} to be the state and costate vectors at t=t0 and t=tf by integrating p(ti) backward and forward, respectively. Thus, the optimization problem is stated as follows:(18)Findptisuchthatxt0∣pti=x0,xtf∣pti=xf.Once the optimal solution at the starting point is obtained by the MBISSM, denoted asp∗(ti), the states at the initial and terminal time are completely determined by integrating from p∗(ti), denoted as x∗(t0∣p∗(ti)) and x∗(tf∣p∗(ti)), respectively. The initial state errors, which are the differences between x∗(t0∣p∗(ti)) and the given initial value x0 defined in (3), can be expressed as(19)Δxt0∣p∗ti=x∗t0∣p∗ti-x0.Denotex(tf∣x0,λ0∗) as the state that is integrated from x0 and λ0∗=λ∗(t0∣p∗(ti)). The terminal state errors are denoted as(20)Δxtf∣x0,λ0∗=xtf∣x0,λ0∗-xf.Due to the existence of the initial state errors defined in (19), the Euclidean norm of the terminal state errors may be extremely larger for the HOCPs, which may violate the error tolerance. In this case, a lower error tolerance for solving (18) is required with the help of multiple-precision arithmetic, until the required accuracy is satisfied.Although multiple shooting methods are wildly demonstrated to be superior to the single shooting methods in terms of the convergence domain [32], they do not provide any information regarding the flow structure in the phase space, which aids the development of a simpler approximate solution method for solving HOCPs [6]. Thus, multiple shooting methods are not considered in this paper.


## 4. Numerical Demonstration

Consider the following hypersensitive optimal control problem adopted from [1], which is to determine the control, u(t), on the time interval t∈[0,tf] to minimize the performance index as follows:(21)J=12∫0tfx12t+x22t+u2tdt.The dynamic constraints are expressed as(22)x˙1t=x2t,x˙2t=-x1t-x13t+utwith the boundary conditions defined as(23)x10=1,(24)x20=0,(25)x1tf=0.75,(26)x2tf=0,where tf is fixed. It is known that, for sufficiently large value of tf, this optimal control problem becomes hypersensitive [1]. Thus, the Hamiltonian function is formulated as(27)Hxt,ut,λt=12x12t+x22t+u2t+λ1tx2t+λ2t-x1t-x13t+ut,where λ(t)=[λ1(t),λ2(t)]T is the costate vector associated with x(t), the governing costate differential equations of which are given as(28)λ˙1t=-∂Hxt,ut,λt∂x1t=-x1t+λ2t+3λ2tx12t,λ˙2t=-∂Hxt,ut,λt∂x2t=-x2t-λ1t.According to optimal control theory [22], the optimal control is determined by(29)∂Hxt,u∗t,λt∂u∗t=u∗t+λ2t=0,that is,(30)u∗t=-λ2t.Substituting (30) into the state equations in (22) yields(31)x˙1t=x2t,x˙2t=-x1t-x13t-λ2twhich constitute the Hamiltonian differential equations together with the costate equations in (28), and the corresponding equilibrium value is peq=[0,0,0,0]T.In this paper, all computations are executed in a desktop computer with a 2.20 GHz CPU. All of the codes are implemented under MATLAB, and the required accuracy is set as10-12.


### 4.1. Optimal Solutions with Double-Precision Arithmetic

In this section, several numerical tests are implemented to demonstrate the performance of the proposed MBISSM with double-precision arithmetic. The middle point withti=tf/2 is chosen as the starting point for the MBISSM. The FFSM and a commercial software GPOPS-II [21], which is a Legendre-Gauss-Radau quadrature orthogonal collocation based direct method, are implemented for the comparison. The simulation results of the BSSM and the BISSM are not provided in this section, as they are quite close with the FFSM’s. Although GPOPS-II adopts a direct method to solve the original optimal control problem described in ((21)–(26)), it provides highly accurate costates estimation, which permits the ability to compare its optimal solution directly with other indirect methods. MATLAB function ode45, which is based on an explicit Runge-Kutta (4,5) formula with adaptive step control, is used as numerical integration algorithm and the error tolerance for all of these three methods is set at 10-12.Eight test cases with varying terminal timetf in the range of [10,200] are implemented, and the numerical solutions are provided in Table 1. In Table 1, comparisons of MBISSM, FFSM, and GPOPS-II in the Euclidean norm of terminal state errors and convergent CPU time are illustrated. The character “X” is used to denote the divergent cases. The initial guesses for these methods are all the same, that is, the equilibrium value peq. As illustrated in Table 1, all the test cases by the FFSM are divergent, and 6 cases are successfully convergent by the MBISSM, which reveals that the MBISSM has much better convergence property when compared with the FSSM. GPOPS-II has the best convergence property without a failure. However, the corresponding Euclidean norm of terminal state errors is much larger than the MBISSM’s. As described in the previous section, the terminal state errors by MBISSM are mainly caused by the tiny initial state errors defined in (19). Unlike the MBISSM, the terminal state errors by the GPOPS-II are totally determined by the number of mesh intervals and the degree of the approximating polynomial within each mesh interval.Table 1
Simulation results by MBISSM, FFSM, and GPOPS-II with double-precision arithmetic.
t f (s) Euclidean norm of terminal state errors Convergent CPU time (s) MBISSM FFSM GPOPS-II MBISSM FFSM GPOPS-II 10 2.819×10-10 X 1.658×10-7 116.6 X 3.4 20 2.237×10-7 X 1.847×10-4 31.2 X 4.2 30 3.632×10-5 X 1.092×10-1 23.1 X 4.3 40 1.45×10-3 X 9.661×101 20.3 X 4.4 50 1.499×101 X 9.829×104 17.1 X 4.3 60 6.472×103 X 5.618×107 17.8 X 4.5 110 X X ∞ X X 5.1 200 X X ∞ X X 5.2Takingtf=50s case, for example, the optimal solution at the middle point is obtained as follows:(32)x1∗tf2=6.45536833498639×10-8,x2∗tf2=-6.61537533239143×10-8,λ1∗tf2=2.99999567350630×10-8,λ2∗tf2=6.45965578930295×10-8.Obviously, the above solution is approximate zero, which indicates that the optimal solution at the middle point is indeed in the close neighborhood of the equilibrium value with peq=[0,0,0,0]T and guarantees the success of the proposed MBISSM with the equilibrium values as the initial guess. Due to the fact that the initial guess is extremely close to the accurate optimal solution, this HOCP can be easily solved in only 6 iterations by the MBISSM. With the above optimal solution at the middle point as the starter, the optimal solution obtained at t0 by the MBISSM can be calculated by backward integration, given as(33)x1∗0=0.999999999999180,x2∗0=-0.000000000000271,λ1∗0=2.603163099258960,λ2∗0=0.236067977499456.Thus the initial state errors defined in (19) are(34)Δx10=x1∗0-x10=-8.2×10-13,Δx20=x2∗0-x20=-2.71×10-13,where the initial conditions x1(0) and x2(0) are defined in ((23)-(24)). The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 3, which illustrates that the terminal state errors defined in (20) reach two orders of magnitude and the maximum state errors are three orders of magnitude for the whole trajectory. Fortunately, with the help of multiple-precision arithmetic, these errors can be significantly reduced to meet the required error tolerance, which will be discussed later.Figure 3
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by MBISSM for terminal timetf=50s with double-precision arithmetic.For comparison, the corresponding initial optimal costates withtf=50s by GPOPS-II are given as follows:(35)λ1∗0=2.603163100374669,λ2∗0=0.236067977366574which are provided by GPOPS-II via a highly accurate costates estimation method. The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 4, which reveals that the terminal state errors reach as high as four orders of magnitude and the maximum state errors are five orders of magnitude for the whole trajectory. Because the large terminal state errors by GPOPS-II are directly determined by the discrete characteristics of direct methods, they cannot be reduced effectively even with multiple-precision arithmetic.Figure 4
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by GPOPS-II for terminal timetf=50s with double-precision arithmetic.


### 4.2. Optimal Solutions with Multiple-Precision Arithmetic

As demonstrated previously, although the proposed MBISSM has much higher convergence property than the FSSM, it still fails to be convergent for larger terminal time and the terminal state errors of convergent case may also be extremely large. Fortunately, both of the above shortcomings can be fully addressed with multiple-precision arithmetic and Taylor series method.Still taking thetf=50s case, for example, the same problem is resolved by the MBISSM with multiple-precision arithmetic and Taylor series method as integration algorithm. The corresponding iterative formulation of Taylor coefficients is(36)x1k+1=1k+1x2k,x2k+1=1k+1-x1k-∑n=0kx1k-n∑p=0nx1n-px1p-λ2k,λ1k+1=1k+1-x1k+λ2k+3∑n=0kλ2k-n∑p=0nx1n-px1p,λ2k+1=1k+1-x2k-λ1k,where k=0,1,…,m-1, m is the order of Taylor series. Multiple-precision arithmetic with 30 significant digits is adapted for high-precision optimization. The integral accuracy for the Taylor series method and the error tolerance for the MBISSM are set at 10-30 and 10-27, respectively. This HOCP can be easily solved in only 10 iterations by the MBISSM, and the convergent CPU time is about 980 s. The optimal solution at the middle point is obtained as follows:(37)x1∗tf2=6.45536835897444817728091151098×10-8,x2∗tf2=-6.61537535985934926089400692506×10-9,λ1∗tf2=2.99999567382193482349955449625×10-8,λ2∗tf2=6.4596557861702446143763973621×10-8.With the above optimal solution at the middle point as the starter, the optimal solution at initial time can be calculated by backward integration, that is,(38)x1∗0=1.00000000000000000000000000000,x2∗0=-2.49674476502450236412812090912×10-28,λ1∗0=2.60316309926254841367437467726,λ2∗0=2.360679774997916414724683759212×10-1.Thus the initial state errors defined in (19) are(39)Δx10=x1∗0-x10=0,Δx20=x2∗0-x20=-2.49674476502450236412812090912×10-28.The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 5, and the Euclidean norm of terminal state errors defined in (20) is about 8.7455×10-14. As the error tolerance is set at a higher accuracy of 10-27, the obtained optimal solution accuracy 8.7455×10-14 still satisfies the required accuracy 10-12. If the accuracy cannot meet the required accuracy, a higher significant digits and a higher integral accuracy are required.Figure 5
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by MBISSM for terminal timetf=50s with 30 significant digits.Another case with terminal timetf=200s is also implemented, which fails with double-precision arithmetic as illustrated in Table 1. For this case, 90 significant digits are used for high-precision computation. The integral accuracy for the Taylor series method and the error tolerance for this single shooting method are set at 10-90 and 10-87, respectively. This problem is easily solved by the proposed high-precision single shooting method within 10 iterations, and the convergent CPU time is about 35 hours. The obtained optimal solution at t=0 is given as(40)x1∗0=1,(41)x2∗0=-5.57488926652873071685938602638128927705372911638218741237980296620560640489334856949249866×10-88,(42)λ1∗0=2.60316309926255757333545904831747724935957041065781988207219177679102777660491239842757615,(43)λ2∗0=2.36067977499789696409173668731276235440618359611525724270862087044025986692985578709359839×10-1.Thus the corresponding initial state errors defined in (19) are(44)Δx10=x1∗0-x10=0,Δx20=x2∗0-x20=-5.57488926652873071685938602638128927705372911638218741237980296620560640489334856949249866×10-88.The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 6, and the Euclidean norm of terminal state errors defined in (20) is about 1.5841×10-29, which meets the required accuracy 10-12.Figure 6
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by MBISSM for terminal timetf=200s with 90 significant digits.Lower significant digits are used, that is, the first 55 significant digits ofλ1∗(0) and λ2∗(0), which are (45)λ1∗0=2.603163099262557573335459048317477249359570410657819882,λ2∗0=2.360679774997896964091736687312762354406183596115257242×10-1.With (23)-(24) and (45) as the initial starter of trajectory integration, the terminal state errors defined in (20) are provided in Table 2. The integral accuracy for the Taylor series method remains unchanged, which is 10-87. As illustrated in Table 2, the Euclidean norm of terminal state errors increases rapidly along with the deceasing of significant digits, which reaches as high as about 2.10452×103 for the 55 significant digits case.Table 2
Comparison of Euclidean norm of terminal state errors by MBISSM with different significant digits.
Significant digits Terminal state errors 85 5.866875762628193513683714399836181990749066000576664570626536775036996028775307020565×10-27 80 1.2729135942414578912275331034096639697677188650860207411613500771425023492253212×10-21 75 8.32800893193697128535236564651012081921145896143028040583575297035124161334×10-17 70 1.100962781462392971848264179350886082349518980048139710140583435933765×10-11 65 1.2582068831722865274799518796876135992863954934328056170904360060×10-6 60 2.67148789475288413499359451519048554798453717463573873843464×10-2 55 2.104524873696096706792124497334559383710609981183820732×103If lower integral accuracies of Taylor series method are used with 90 significant digits, the terminal state errors defined in (20) are provided in Table 3, taking (23)-(24) and (42)-(43) as the initial starter of trajectory integration. As illustrated in Table 3, the terminal state errors increase rapidly along with the raising of integral accuracy. It is amazingly found that the terminal state errors are four orders of magnitude large if the integral accuracy is set at 10-50, which obviously cannot meet the required error tolerance. Based on the above discussions, it can be concluded that no optimal solution may exist to meet the required error tolerance for lower significant digits and lower integral accuracy, which explains why the MBISSM fails to solve the same problem with double-precision arithmetic in section.Table 3
Comparison of the Euclidean norm of terminal state errors by MBISSM with different integral accuracies of Taylor series method.
Integral accuracy Terminal state errors 10-85 1.58319633679772217944824131350505236807202739692078827429906812432096942480802199853304275×10-29 10-80 3.72525147090499732555084530502769934173983207752374155382811547722333645443103106733571409×10-27 10-75 3.81437611392372690108266166748050457344969141070869052138948837041124626661690251076409706×10-22 10-70 1.13268310209557943195606354923552665677244000325670771361005795912304125278779540527586414×10-16 10-65 5.03071114527927667540803559488399586213189718894019354724386973494151029974516002963841624×10-12 10-60 4.71217759584233098925431899981529901744063117797769340938845177106494597110360867737808901×10-7 10-55 5.88575074985461247363291315566785488575603113859669264078215794228281756486137147031727330×10-2 10-50 8.64647234216576719278782329532914219591521188043541077338325766415336171636888697563052874×103


## 4.1. Optimal Solutions with Double-Precision Arithmetic

In this section, several numerical tests are implemented to demonstrate the performance of the proposed MBISSM with double-precision arithmetic. The middle point withti=tf/2 is chosen as the starting point for the MBISSM. The FFSM and a commercial software GPOPS-II [21], which is a Legendre-Gauss-Radau quadrature orthogonal collocation based direct method, are implemented for the comparison. The simulation results of the BSSM and the BISSM are not provided in this section, as they are quite close with the FFSM’s. Although GPOPS-II adopts a direct method to solve the original optimal control problem described in ((21)–(26)), it provides highly accurate costates estimation, which permits the ability to compare its optimal solution directly with other indirect methods. MATLAB function ode45, which is based on an explicit Runge-Kutta (4,5) formula with adaptive step control, is used as numerical integration algorithm and the error tolerance for all of these three methods is set at 10-12.Eight test cases with varying terminal timetf in the range of [10,200] are implemented, and the numerical solutions are provided in Table 1. In Table 1, comparisons of MBISSM, FFSM, and GPOPS-II in the Euclidean norm of terminal state errors and convergent CPU time are illustrated. The character “X” is used to denote the divergent cases. The initial guesses for these methods are all the same, that is, the equilibrium value peq. As illustrated in Table 1, all the test cases by the FFSM are divergent, and 6 cases are successfully convergent by the MBISSM, which reveals that the MBISSM has much better convergence property when compared with the FSSM. GPOPS-II has the best convergence property without a failure. However, the corresponding Euclidean norm of terminal state errors is much larger than the MBISSM’s. As described in the previous section, the terminal state errors by MBISSM are mainly caused by the tiny initial state errors defined in (19). Unlike the MBISSM, the terminal state errors by the GPOPS-II are totally determined by the number of mesh intervals and the degree of the approximating polynomial within each mesh interval.Table 1
Simulation results by MBISSM, FFSM, and GPOPS-II with double-precision arithmetic.
t f (s) Euclidean norm of terminal state errors Convergent CPU time (s) MBISSM FFSM GPOPS-II MBISSM FFSM GPOPS-II 10 2.819×10-10 X 1.658×10-7 116.6 X 3.4 20 2.237×10-7 X 1.847×10-4 31.2 X 4.2 30 3.632×10-5 X 1.092×10-1 23.1 X 4.3 40 1.45×10-3 X 9.661×101 20.3 X 4.4 50 1.499×101 X 9.829×104 17.1 X 4.3 60 6.472×103 X 5.618×107 17.8 X 4.5 110 X X ∞ X X 5.1 200 X X ∞ X X 5.2Takingtf=50s case, for example, the optimal solution at the middle point is obtained as follows:(32)x1∗tf2=6.45536833498639×10-8,x2∗tf2=-6.61537533239143×10-8,λ1∗tf2=2.99999567350630×10-8,λ2∗tf2=6.45965578930295×10-8.Obviously, the above solution is approximate zero, which indicates that the optimal solution at the middle point is indeed in the close neighborhood of the equilibrium value with peq=[0,0,0,0]T and guarantees the success of the proposed MBISSM with the equilibrium values as the initial guess. Due to the fact that the initial guess is extremely close to the accurate optimal solution, this HOCP can be easily solved in only 6 iterations by the MBISSM. With the above optimal solution at the middle point as the starter, the optimal solution obtained at t0 by the MBISSM can be calculated by backward integration, given as(33)x1∗0=0.999999999999180,x2∗0=-0.000000000000271,λ1∗0=2.603163099258960,λ2∗0=0.236067977499456.Thus the initial state errors defined in (19) are(34)Δx10=x1∗0-x10=-8.2×10-13,Δx20=x2∗0-x20=-2.71×10-13,where the initial conditions x1(0) and x2(0) are defined in ((23)-(24)). The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 3, which illustrates that the terminal state errors defined in (20) reach two orders of magnitude and the maximum state errors are three orders of magnitude for the whole trajectory. Fortunately, with the help of multiple-precision arithmetic, these errors can be significantly reduced to meet the required error tolerance, which will be discussed later.Figure 3
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by MBISSM for terminal timetf=50s with double-precision arithmetic.For comparison, the corresponding initial optimal costates withtf=50s by GPOPS-II are given as follows:(35)λ1∗0=2.603163100374669,λ2∗0=0.236067977366574which are provided by GPOPS-II via a highly accurate costates estimation method. The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 4, which reveals that the terminal state errors reach as high as four orders of magnitude and the maximum state errors are five orders of magnitude for the whole trajectory. Because the large terminal state errors by GPOPS-II are directly determined by the discrete characteristics of direct methods, they cannot be reduced effectively even with multiple-precision arithmetic.Figure 4
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by GPOPS-II for terminal timetf=50s with double-precision arithmetic.


## 4.2. Optimal Solutions with Multiple-Precision Arithmetic

As demonstrated previously, although the proposed MBISSM has much higher convergence property than the FSSM, it still fails to be convergent for larger terminal time and the terminal state errors of convergent case may also be extremely large. Fortunately, both of the above shortcomings can be fully addressed with multiple-precision arithmetic and Taylor series method.Still taking thetf=50s case, for example, the same problem is resolved by the MBISSM with multiple-precision arithmetic and Taylor series method as integration algorithm. The corresponding iterative formulation of Taylor coefficients is(36)x1k+1=1k+1x2k,x2k+1=1k+1-x1k-∑n=0kx1k-n∑p=0nx1n-px1p-λ2k,λ1k+1=1k+1-x1k+λ2k+3∑n=0kλ2k-n∑p=0nx1n-px1p,λ2k+1=1k+1-x2k-λ1k,where k=0,1,…,m-1, m is the order of Taylor series. Multiple-precision arithmetic with 30 significant digits is adapted for high-precision optimization. The integral accuracy for the Taylor series method and the error tolerance for the MBISSM are set at 10-30 and 10-27, respectively. This HOCP can be easily solved in only 10 iterations by the MBISSM, and the convergent CPU time is about 980 s. The optimal solution at the middle point is obtained as follows:(37)x1∗tf2=6.45536835897444817728091151098×10-8,x2∗tf2=-6.61537535985934926089400692506×10-9,λ1∗tf2=2.99999567382193482349955449625×10-8,λ2∗tf2=6.4596557861702446143763973621×10-8.With the above optimal solution at the middle point as the starter, the optimal solution at initial time can be calculated by backward integration, that is,(38)x1∗0=1.00000000000000000000000000000,x2∗0=-2.49674476502450236412812090912×10-28,λ1∗0=2.60316309926254841367437467726,λ2∗0=2.360679774997916414724683759212×10-1.Thus the initial state errors defined in (19) are(39)Δx10=x1∗0-x10=0,Δx20=x2∗0-x20=-2.49674476502450236412812090912×10-28.The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 5, and the Euclidean norm of terminal state errors defined in (20) is about 8.7455×10-14. As the error tolerance is set at a higher accuracy of 10-27, the obtained optimal solution accuracy 8.7455×10-14 still satisfies the required accuracy 10-12. If the accuracy cannot meet the required accuracy, a higher significant digits and a higher integral accuracy are required.Figure 5
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by MBISSM for terminal timetf=50s with 30 significant digits.Another case with terminal timetf=200s is also implemented, which fails with double-precision arithmetic as illustrated in Table 1. For this case, 90 significant digits are used for high-precision computation. The integral accuracy for the Taylor series method and the error tolerance for this single shooting method are set at 10-90 and 10-87, respectively. This problem is easily solved by the proposed high-precision single shooting method within 10 iterations, and the convergent CPU time is about 35 hours. The obtained optimal solution at t=0 is given as(40)x1∗0=1,(41)x2∗0=-5.57488926652873071685938602638128927705372911638218741237980296620560640489334856949249866×10-88,(42)λ1∗0=2.60316309926255757333545904831747724935957041065781988207219177679102777660491239842757615,(43)λ2∗0=2.36067977499789696409173668731276235440618359611525724270862087044025986692985578709359839×10-1.Thus the corresponding initial state errors defined in (19) are(44)Δx10=x1∗0-x10=0,Δx20=x2∗0-x20=-5.57488926652873071685938602638128927705372911638218741237980296620560640489334856949249866×10-88.The integrated trajectory with x1(0), x2(0), λ1∗(0), and λ2∗(0) as the initial starter is shown in Figure 6, and the Euclidean norm of terminal state errors defined in (20) is about 1.5841×10-29, which meets the required accuracy 10-12.Figure 6
States profiles along the trajectory integrated from the fixed initial states and the corresponding optimal initial costates obtained by MBISSM for terminal timetf=200s with 90 significant digits.Lower significant digits are used, that is, the first 55 significant digits ofλ1∗(0) and λ2∗(0), which are (45)λ1∗0=2.603163099262557573335459048317477249359570410657819882,λ2∗0=2.360679774997896964091736687312762354406183596115257242×10-1.With (23)-(24) and (45) as the initial starter of trajectory integration, the terminal state errors defined in (20) are provided in Table 2. The integral accuracy for the Taylor series method remains unchanged, which is 10-87. As illustrated in Table 2, the Euclidean norm of terminal state errors increases rapidly along with the deceasing of significant digits, which reaches as high as about 2.10452×103 for the 55 significant digits case.Table 2
Comparison of Euclidean norm of terminal state errors by MBISSM with different significant digits.
Significant digits Terminal state errors 85 5.866875762628193513683714399836181990749066000576664570626536775036996028775307020565×10-27 80 1.2729135942414578912275331034096639697677188650860207411613500771425023492253212×10-21 75 8.32800893193697128535236564651012081921145896143028040583575297035124161334×10-17 70 1.100962781462392971848264179350886082349518980048139710140583435933765×10-11 65 1.2582068831722865274799518796876135992863954934328056170904360060×10-6 60 2.67148789475288413499359451519048554798453717463573873843464×10-2 55 2.104524873696096706792124497334559383710609981183820732×103If lower integral accuracies of Taylor series method are used with 90 significant digits, the terminal state errors defined in (20) are provided in Table 3, taking (23)-(24) and (42)-(43) as the initial starter of trajectory integration. As illustrated in Table 3, the terminal state errors increase rapidly along with the raising of integral accuracy. It is amazingly found that the terminal state errors are four orders of magnitude large if the integral accuracy is set at 10-50, which obviously cannot meet the required error tolerance. Based on the above discussions, it can be concluded that no optimal solution may exist to meet the required error tolerance for lower significant digits and lower integral accuracy, which explains why the MBISSM fails to solve the same problem with double-precision arithmetic in section.Table 3
Comparison of the Euclidean norm of terminal state errors by MBISSM with different integral accuracies of Taylor series method.
Integral accuracy Terminal state errors 10-85 1.58319633679772217944824131350505236807202739692078827429906812432096942480802199853304275×10-29 10-80 3.72525147090499732555084530502769934173983207752374155382811547722333645443103106733571409×10-27 10-75 3.81437611392372690108266166748050457344969141070869052138948837041124626661690251076409706×10-22 10-70 1.13268310209557943195606354923552665677244000325670771361005795912304125278779540527586414×10-16 10-65 5.03071114527927667540803559488399586213189718894019354724386973494151029974516002963841624×10-12 10-60 4.71217759584233098925431899981529901744063117797769340938845177106494597110360867737808901×10-7 10-55 5.88575074985461247363291315566785488575603113859669264078215794228281756486137147031727330×10-2 10-50 8.64647234216576719278782329532914219591521188043541077338325766415336171636888697563052874×103


## 5. Conclusion

In this paper, a new high-precision single shooting method is presented to solve the challenging hypersensitive optimal control problem, which consists of multiple-precision arithmetic, Taylor series method, and a new MBISSM. Multiple-precision arithmetic and Taylor series method are introduced to provide higher significant digits and integral accuracy, and the MBISSM is developed to generate appropriate initial guess, which is extremely close to the optimal solution by utilizing three-segment structure characteristic of HOCPs. Thus, the two main difficulties, the possible nonexistence of the optimal solution to meet the required tolerance under double-precision arithmetic and the hypersensitivity of the optimal solution with respect to the initial conditions by indirect methods, are satisfactorily solved. The numerical simulations demonstrate the efficiency of this proposed method, the optimal solution of which can be easily solved within several iterations. The simulation results also reveal the fact that no optimal solution may exist to meet the required error tolerance for lower significant digits and lower integral accuracy, which is still not widely recognized yet in optimization engineering.The main disadvantage of the proposed high-precision single shooting method is that it may require much larger CPU time. Taking thetf=200s case with 90 significant digits, for example, it may take several hours to complete a single iteration. Thus, a faster computation ability is required for future researches.



---

*Source: 7908378-2018-04-15.xml*