The Adjoint Method for Time-Optimal Control Problems

Among all optimal control problems, the solution of time-optimal problems is probably the most challenging ones. In this paper, we consider a dynamic system described by a set of minimal coordinates and given initial values. We are looking for the optimal control of the system such that the time for a prescribed process or maneuver becomes a minimum.

In order to solve such problems, fundamental books on optimal control as, e.g., Refs. [1] and [2] discuss direct and indirect methods. In summary, direct methods, also addressed as nonlinear programing, transform the dynamic problem to a static problem by a parameterization of the control variables, see Ref. [3], for example. Numerous methods follow a direct approach, e.g., interior point algorithms, active set, sequential quadratic programing, pseudo-spectral method, to name but a few. Indirect methods are based on the solution of the optimality conditions stated by Pontryagin resulting in a two-point boundary value problem. Solution strategies for the boundary value problem are discussed by various authors [4–7] and are already well established in a wide range in engineering sciences. The boundary value problem can, e.g., be solved with shooting methods or with discretization techniques in which both states and controls are discretized. For more details to direct and indirect methods, a good overview of both is given by Bianco et al. [8]. In a previous work, Bottasso et al. [9] presented a combined indirect approach for solving inverse and trajectory optimization problems in multibody dynamics, which also corresponds to the ideas presented in Ref. [5]. The latter mentioned work deals with the time-optimal control of a vehicle using an indirect approach to solve the boundary value problem.

A first significant formulation of the time-optimal control problem of a general vehicle maneuver has been stated by Hendrikx et al. [10]. This was followed by further articles concerning motorcycle lap time simulations in Ref. [4] and car lap time simulations discussed in Refs. [6] and [7]. Basically, all the mentioned authors have solved the boundary value problems associated with the optimal control problems. Time-optimal control problems play a significant role as well in robotic applications, see, e.g., Ref. [11] for the time-optimal control problem with specified trajectory. More recently, in Ref. [12] time-optimal motions following a predefined end-effector path are considered for kinematically redundant robots as well as for nonredundant robots. Here, a multiple shooting approach is used to solve the boundary value problem.

An indirect alternative to solving the underlying boundary value problem is given by gradient methods, independently developed by Kelley [13] and by Bryson and Denham [14] already in the sixties of the past century.

Here, the basic idea is to find a control history leading to a minimum of a cost functional. Numerous strategies are available to compute the control for which a functional attains a minimum, as, e.g., the method of the steepest descent, the conjugate gradient method, the gradient projection method, the Gauss-Newton method, or quasi-Newton methods like the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm when estimating the Hessian, see, e.g., Ref. [15] for an overview and comparison of performance behavior. In any cases, the gradient of the cost functional with respect to the control history has to be computed. A control signal is often time discretized in order to obtain a finite dimensional problem. In case of a time discretized control variable, the computational efforts are extremely high due to the numerical evaluation of the gradient at every time-step. However, the effort of the gradient computation can be significantly reduced by using the adjoint method, see, e.g., the pioneering work by Bryson and Ho [16]. The basic idea of this approach is the introduction of additional adjoint variables determined by a set of adjoint differential equations from which the gradient can be computed more efficiently than, e.g., with a numerical (direct) differentiation. To determine the partial derivatives with respect to each time-discrete control parameter by means of a numerical differentiation, the entire forward system must be solved once for a variation of every discretization parameter. With the adjoint method, however, the gradient can be computed very efficiently by solving two systems of ordinary differential equations only, i.e., by the solutions of the system equations and of the so-called adjoint equations. With the gradient information, the control input can be improved iteratively converging to an optimal control. Even by choosing a suboptimal initial guess for the control inputs, in general, the proposed method converges to the (local) optimal solution.

Various authors have already utilized this method in the field of inverse dynamics of multibody systems more than 20 years ago, as, e.g., Refs. [17] and [18], but nevertheless, the topic still provides open research potential.

More recently, Refs. [19] and [20] show how the adjoint method can be applied efficiently to a multibody system described by differential-algebraic equations of index three. A piecewise adjoint method is presented in Ref. [21], in which the coordinate partitioning underlying ordinary differential equations are formulated as a boundary value problem, and are solved by multiple shooting methods. It has to be emphasized that the adjoint method can be utilized to solve even highly nonlinear systems with a vast number of degrees-of-freedom in a straight forward way, as, e.g., in the application of second-order adjoint sensitivity analysis of large multibody systems in Ref. [22], in which a comparison between the adjoint method and direct differentiation shows a significant reduction in computational costs.

In this paper, a gradient-based solution approach using the adjoint method for time-optimal control problems for dynamic systems is presented, which has not yet been published in the open literature so far. The main novelty lies in the adaption of the adjoint gradient computation to a functional with variable integration limits. For that purpose, the adjoint method is formulated for two different domain approaches concerning a transformation between the time coordinate and a coordinate within a fixed interval essential for free final time problems. On the one hand, the adjoint method in a hybrid domain gives the advantage of solving the system and adjoint equations in time domain using classical solver and step size strategies, but defining the control and the gradient in the space domain. On the other hand, the adjoint method is presented in complete state space domain, in which the state and adjoint equations, cost functional, and control undergo a total transformation to a coordinate with fixed boundaries. Both methods arise as innovative tools for the efficient computation of the gradient of a cost functional for all kinds of optimization methods. Due to the novel consistent formulation, both presented methods even allow the application of preferable quasi-Newton methods using the BFGS-update for the Hessian. Therefore, time-optimal control problems can now be considered in an efficient way concerning convergence properties.

where $P(x(t),u(t))$ is a proper penalty function which is zero if the states and control variables satisfy the constraint conditions, and increases rapidly if they are not satisfied. Moreover, P should be continuously differentiable at least once. Unless violating the constraints, the cost functional J is simply the length of the time interval $tf−t0$⁠. The time-optimal control problem is then formulated as the problem of finding controls $u(t)$ in the time interval $t∈[t0,tf]$ which minimizes the functional J.

So far, we have not made any restrictions to the adjoint variables, now we shall find it convenient to choose $p(t)$ such that all terms multiplied with $δx$ vanish. Thus, the adjoint variables are defined by the linear final value problem

where ε must be a sufficiently small positive number.

A problem with the presented adjoint method arises because the modified control $u*(t)=u(t)+δu(t)$ also causes a change of the time interval $[t0,tf]$ to $[t0,tf*]$⁠. If $tf*<tf$⁠, which should be the case for a control variation of the form like in Eq. (14), the improved control $u*(t)$ can be truncated at $t=tf*$ although the information about $u*(t)$ is available also for $tf*<t<tf$⁠. However, if ε was selected too large, the terminal condition may be eventually met at a later time instance $tf*>tf$⁠. In this case, we have no information about $u*(t)$ for $tf<t<tf*$⁠.

which ranges between the fixed boundaries $s0=Φ(x(t0))$ and $sf=Φ(x(tf))$⁠. Based on this transformation, two solution strategies concerning the adjoint gradient method are pursued in this paper: (a) a hybrid formulation, where the state and adjoint equations are solved in the original time domain and only the controls and the gradient are defined in the fixed interval, and (b) the complete transformation of the time coordinate to a fixed interval. Both methods are discussed in detail below.

To eliminate the influence of all terms multiplied with $δx(t)$⁠, we now define the adjoint variables in the following way:

where ε is a sufficiently small positive number.

With this result, we are ready to apply a gradient-based optimization procedure for finding a minimum of J, which can be summarized as follows:

1. Choose an initial guess for the control history $u(s)$ which causes the states to satisfy the terminal condition in Eq. (2).

2. Solve Eq. (1) for the state variables $x(t)$ until Eq. (2) is satisfied. During the numerical integration of the state equations, the argument s in $u(s)$ must be replaced by the original time t with Eq. (15).

3. Solve Eq. (21) for the adjoint variables $p(t)$ by a backward integration starting at t = t_f.

4. Select a value for ε and compute the control variation $δu(s)$ for every s from Eq. (25) by involving the inverse function of s(t). Alternatively, a quasi-Newton method can be applied to the update formula in Eq. (25), e.g., the BFGS-method.

5. Update the control history by setting $u(s)→u(s)+δu(s)$⁠. With the updated control $u(s)$⁠, the cost functional should be reduced, if ε was selected sufficiently small.

6. Vary ε, resolve the state equations and evaluate the cost functional until an optimal value for ε is found, which minimizes the cost functional along the search direction. For that purpose, a line search algorithm as described in Ref. [24] can be used.

7. Replace the control history $u(s)$ by the updated control and repeat steps 2–6. Terminate the optimization procedure when Eq. (24) is sufficiently small (zero gradient).

In summary, the advantage of the presented method is that the equations of motion and also the adjoint system can be solved in the time domain and only the control variables and the gradient formula are transformed into the state space. Unfortunately, this requires the derivative of the controls with respect to s which can become infinite for bang–bang controls and could cause an undesirable excitation term in the adjoint equations in Eq. (21a). In the next paragraph, we show that this problem can be overcome by an elimination of the time coordinate also in the state equations in Eq. (1a). However, the price to pay is that a standard multibody simulation software is probably no longer applicable.

The adjoint variables $p(s)$ and $ξ(s)$ can now be again defined such that all terms with $δx(s)$ and $δt(s)$ disappear. This is the case if $ξ′(s)=0$ and $ξ(sf)=1$⁠, i.e., $ξ(s)=1$ for $s∈[s0,sf]$⁠, and if

in which ε is again a small positive number determining the size of the control update.

Finally, we summarize the main steps of the algorithm for the s-domain:

1. Select a control history $u(s)$ in the s-interval $[s0,sf]$⁠.

2. Solve the state equations in Eq. (27) for $x(s)$ and the differential equation in Eq. (28) for t(s), simultaneously.

3. Solve the adjoint differential equations in Eq. (34) for the adjoint variables $p(s)$ backward in s.

4. Compute the control update $δu(s)$ for every s by using Eq. (36). Note, instead of using the update $δu(s)$ directly, Eq. (36) can be used to apply a quasi-Newton method.

5. Improve the control history by setting $u(s)→u(s)+δu(s)$⁠. If ε was selected sufficiently small, the updated control $u(s)$ reduces the cost functional in Eq. (29), in which the final time is given by $tf=t(sf)$⁠.

6. Find a suitable factor ε, e.g., by using a line search algorithm.

7. Repeat steps 2–6 and stop the algorithm when $δJ¯$ in Eq. (35) is sufficiently small.

Now we consider two numerical examples in order to apply the proposed methods for computing the time-optimal control for dynamic systems. For the first example, the adjoint method in a hybrid domain is used for trajectory planning of a space flight, while in the second example the adjoint method in the s-domain is utilized for computing time-optimal vehicle maneuvers.

As a first example, we consider the time-optimal flight to the sphere of influence. The sphere of influence is a domain around a celestial body that indicates where the gravity of the body has a dominating effect on the dynamics of other space objects [25]. We are interested in the optimal control of a space vehicle so that we reach the sphere of influence in minimal time. We consider only planar motions described in polar coordinates r(t) and $φ(t)$⁠, see Fig. 3.

i.e., the length of the time interval $tf−t0$⁠.

For the numerical example, the following parameters are used: The specific impulse $4.21×104 s$⁠, the thrust force $F=1033 kN$⁠, the radius of the Earth $ρ=6375 km$⁠, the gravitational acceleration $g=9.8106×10−3 km/s2$⁠, and the radius of the sphere of influence $rfinal=9.25×105 km$⁠. Initially, the space vehicle is in a circular Earth orbit with altitude $h=300 km$⁠. Hence, the initial conditions at $t=t0=0$ are given by: $r(0)=ρ+h, φ(0)=0, vr(0)=0, vt(0)=gρ2/(ρ+h)$ and $m(0)=115 t$⁠. We consider two different initial controls, which are both far away from the optimal solution: (a) u(r) = 0 and (b) $u(r)=r π/(2rfinal)$⁠. In the first case, the thrust always acts in tangential direction. In the latter case, the thrust acts tangentially at the beginning and radially at the end, i.e., the sphere of influence is crossed in normal direction. The controls u(r) are discretized by 500 data points in time. Notice that a traditional gradient computation by direct differentiation would require 500 additional forward solutions of the state equations, whereas with our method the gradient is obtained already from one solution of the adjoint system equations. Figure 4 depicts the decrease of the final time with an increasing number of quasi-Newton iterations for both initial controls. The final time $1.8106×104 s$ for the initial control (a) and $1.5704×104 s$ for the initial control (b) are reduced to $tf=1.3187×104 s$ after the optimization process.

The time-optimal control angles u(r) computed from the initial controls (a) and (b) are shown in Fig. 5. The state histories r(t) and $φ(t)$ for the space vehicle are summarized in Fig. 6 and the flight trajectory of the space vehicle is plotted in Fig. 7. They are almost identical for both initial controls. The spacecraft moves nearly on a straight flight path. Notice that for the initial control (a) some small oscillations remain in the control u(r) at the end of the maneuver. The reason for this numerical artifact is that the final value $u(rfinal)$ is never updated by the gradient formula in Eq. (41) due to the zero values for p₃ and p₄ at t = t_f and, hence, the control keeps its original value u = 0. Clearly, the final control value has no influence on the flight time. For the initial control (a), the algorithm converges to a step function at $r=rfinal$ which requires a higher number of iterations. However, as it can be seen from the convergence plot in Fig. 4, the cost functional no longer decreases already after 20 iterations. This can be explained by the fact that the control history close to t = t_f has only small influence on the cost functional. Hence, the oscillations in u(r) do not change the flight time and the flight path very much. However, Fig. 4 shows that choosing the initial control (b), for which the final value at t = t_f coincides with the expected value $π/2$⁠, increases the convergence rate and the oscillations at the end of the optimal control history in Fig. 5 vanish.

As a second example, we consider a mechanical model of a car for which accelerating, braking, and steering should be determined such that the time for a vehicle maneuver becomes minimal.

To describe the physical behavior of the vehicle, we use a single-track model and a Pacejka tire model for the road to wheel contact, see Fig. 8. The position of the vehicle is described by the absolute coordinates x and y of its center of gravity (CoG) and by the orientation angle ψ with respect to the x-axis. Moreover, m is the total mass of the vehicle and J denotes its moment of inertia. Both wheels have the effective radius R and the moment of inertia J_w about their rotation axis. The lengths a and b are the distances measured from the center of gravity to the rear axle and to the front axle. The control inputs of the model are the steering angle u₁ of the front axis and the driving/braking torque u₂ on the rear axis.

Since, for the Pacejka tire model, the tire forces $F1,X, F1,Y, F2,X$⁠, and $F2,Y$ must be described in the body-fixed frame of the chassis, it is convenient to formulate the equations of motion with respect to body-fixed coordinates, as proposed in Refs. [6] and [26].

where $ωψ$ is the angular velocity $ψ̇$ of the chassis about the vertical and ω₁ and ω₂ denote the angular velocities of the wheels.

with $i=X,Y$⁠. The numerical values for the tire parameters a_i, b_i, c_i, μ_i, k_i, d_i, and e_i are shown in Table 1.

where $xc(s)$ and $yc(s)$ are the coordinates of the centerline in arc length description and primes denote derivatives with respect to s. Using, for example, a Newton scheme, s and r can be expressed by x and y from Eq. (46).

Here, $r¯=w−h$⁠, where w is the half road width and h is a small thickness parameter to activate the penalty function before the track constraint is violated. The penalty function is zero if $|r|<r¯$ and increases rapidly if $|r|>r¯$⁠. Note, that $P1(r)$ is quadratic in r if $r¯<|r|<r*$⁠, whereas it increases linearly if $|r|>r*$⁠. Hence, $P1(r)$ is a once continuously differentiable function which is necessary to avoid jumps in the adjoint equations.

Here, $u¯2$ is the maximum driving torque and $u2*$ separates the domains of quadratic and linear increase.