Approximate Solutions to Nonlinear Optimal Control Problems in Astrodynamics

Abstract

A method to solve nonlinear optimal control problems is proposed in this work. The method implements an approximating sequence of time-varying linear quadratic regulators that converge to the solution of the original, nonlinear problem. Each subproblem is solved by manipulating the state transition matrix of the state-costate dynamics. Hard, soft, and mixed boundary conditions are handled. The presented method is a modified version of an algorithm known as “approximating sequence of Riccati equations.” Sample problems in astrodynamics are treated to show the effectiveness of the method, whose limitations are also discussed.

1. Introduction

Optimal control problems are solved with indirect or direct methods. Indirect methods stem from the calculus of variations [1, 2]; direct methods use a nonlinear programming optimization [3, 4]. Both methods require the solution of a complex set of equations (Euler-Lagrange differential equations or Karush-Kuhn-Tucker algebraic equations) for which iterative numerical methods are used. These iterative procedures implement some form of Newton’s method to find the zeros of a nonlinear function. They are initiated by providing an initial guess solution. Guessing an appropriate initial solution is not trivial and requires a deep knowledge of the problem at hand. In indirect methods, the initial value of the Lagrange multiplier has to be provided, whose lack of physical meaning makes it difficult to formulate a good guess. In direct methods, the initial trajectory and control have to be guessed at discrete points over the whole time interval.

This paper presents an approximate method to solve nonlinear optimal control problems. This is a modification of the method known as “approximating sequence of Riccati equations” (ASRE) [5, 6]. It transforms the nonlinear dynamics and objective function into a pseudolinear and quadratic-like structure, respectively, by using state- and control-dependent functions. At each iteration, these functions are evaluated by using the solutions at the previous iteration, and therefore, a series of time-varying linear quadratic regulators is treated. This sequence is solved with a state transition matrix approach, where three different final conditions are handled: final state fully specified, final state not specified, and final state not completely specified. These define hard, soft, and mixed constrained problems, respectively.

The main feature of the presented method is that it does not require guessing any initial solution or Lagrange multiplier. In fact, iterations start by evaluating the state- and control-dependent functions using the initial condition and zero control, respectively. The way the dynamics and objective function are factorized recalls the state-dependent Riccati equations (SDRE) method [7–9]. These two methods possess some similarities, although the way they solve the optimal control problem is different. As the method is approximated, suboptimal solutions are derived. These could be used as first guess solutions for either indirect or direct methods.

2. The Nonlinear Optimal Control Problem

The optimal control problem requires that, given a set of 𝑛 first-order differential equations 𝐱̇=𝐟(𝐱,𝐮,𝑡), (1) the 𝑚 control functions 𝐮(𝑡) must be determined within initial, final time 𝑡_𝑖, 𝑡_𝑓, such that the performance index 𝐽=𝜑𝐱𝑡_𝑓,𝑡_𝑓+^𝑡𝑓_𝑡𝑖𝐿(𝐱,𝐮,𝑡)d𝑡 (2) is minimized while satisfying 𝑛+𝑞 two-point conditions 𝐱𝑡_𝑖=𝐱_𝑖,𝝍𝐱𝑡_𝑓,𝑡_𝑓=0. (3)

The problem consists in finding a solution that represents a stationary point of the augmented performance index 𝐽=𝜑𝐱𝑡_𝑓,𝑡_𝑓+𝝂^𝑇𝝍𝐱𝑡_𝑓,𝐮𝑡_𝑓,𝑡_𝑓+^𝑡𝑓_𝑡𝑖𝐿(𝐱,𝐮,𝑡)+𝝀^𝑇(𝐟(𝐱,𝐮,𝑡)−𝐱̇)d𝑡, (4) where 𝝀 is the vector of costate and 𝝂 is the multiplier of the boundary condition. The necessary conditions for optimality, also referred to as Euler-Lagrange equations, are 𝐱̇=𝜕𝐻,̇𝜕𝝀𝝀=−𝜕𝐻,𝜕𝐱𝜕𝐻𝜕𝐮=0, (5) where 𝐻, the Hamiltonian, is 𝐻(𝐱,𝝀,𝐮,𝑡)=𝐿(𝐱,𝐮,𝑡)+𝝀^𝑇𝐟(𝐱,𝐮,𝑡). (6) The differential-algebraic system (5) must be solved together with the final boundary conditions (3) and the transversality conditions 𝝀𝑡_𝑓=𝜕𝜑+𝜕𝐱𝜕𝝍𝜕𝐱^𝑇𝝂_{𝑡=𝑡𝑓}, (7) which define a differential-algebraic parametric two-point boundary value problem whose solution supplies 𝝂 and the functions 𝐱(𝑡), 𝝀(𝑡), 𝐮(𝑡), 𝑡∈[𝑡_𝑖,𝑡_𝑓].

3. The Approximating Sequence of Riccati Equations

Let the controlled dynamics (1) be rewritten in the form 𝐱̇=𝐴(𝐱,𝑡)𝐱+𝐵(𝐱,𝐮,𝑡)𝐮, (8) and let the objective function (2) be rearranged as 1𝐽=2𝐱^𝑇𝑡_𝑓𝑆𝐱𝑡_𝑓,𝑡_𝑓𝐱𝑡_𝑓+12^𝑡𝑓_𝑡𝑖𝐱^𝑇𝑄(𝐱,𝑡)𝐱+𝐮^𝑇𝑅(𝐱,𝐮,𝑡)𝐮d𝑡, (9) where the operators 𝐴, 𝐵, 𝑆, 𝑄, and 𝑅 have appropriate dimensions. The nonlinear dynamics (8) and the performance index (9) define an optimal control problem. The initial state, 𝐱_𝑖, is assumed to be given, while the final condition (𝝍 in (3)) can assume three different forms (see Section 4). The problem is formulated as an approximating sequence of Riccati equations. This method reduces problem (8)-(9) to a series of time-varying linear quadratic regulators that are defined by evaluating the state- and control-dependent matrices using the solution of the previous iteration (the first iteration considers the initial condition and zero control).

The initial step consists in solving problem 0, which is defined as follows: 𝐱̇⁽⁰⁾=𝐴𝐱_𝑖,𝑡𝐱⁽⁰⁾+𝐵𝐱_𝑖,𝟎,𝑡𝐮⁽⁰⁾,1𝐽=2𝐱^(0)𝑇𝑡_𝑓𝑆𝐱_𝑖,𝑡_𝑓𝐱⁽⁰⁾𝑡_𝑓+12^𝑡𝑓_𝑡𝑖𝐱^(0)𝑇𝑄𝐱_𝑖,𝑡𝐱⁽⁰⁾+𝐮^(0)𝑇𝑅𝐱_𝑖,𝟎,𝑡𝐮⁽⁰⁾d𝑡. (10) Problem 0 is a standard time-varying linear quadratic regulator (TVLQR), as the arguments of 𝐴, 𝐵, 𝑆, 𝑄, and 𝑅 are all given except for the time. This problem is solved to yield 𝐱⁽⁰⁾(𝑡) and 𝐮⁽⁰⁾(𝑡), 𝑡∈[𝑡_𝑖,𝑡_𝑓], where the superscript denotes the problem that the solution refers to.

At a generic, subsequent iteration, problem 𝑘 has to be solved. This is defined as follows: 𝐱̇^(𝑘)=𝐴𝐱^(𝑘−1)(𝑡),𝑡𝐱^(𝑘)+𝐵𝐱^(𝑘−1)(𝑡),𝐮^(𝑘−1)(𝑡),𝑡𝐮^(𝑘),1𝐽=2𝐱^(𝑘)𝑇𝑡_𝑓𝑆𝐱^(𝑘−1)𝑡_𝑓,𝑡_𝑓𝐱^(𝑘)𝑡_𝑓+12^𝑡𝑓_𝑡𝑖𝐱^(𝑘)𝑇𝑄𝐱^(𝑘−1)(𝑡),𝑡𝐱^(𝑘)+𝐮^(𝑘)𝑇𝑅𝐱^(𝑘−1)(𝑡),𝐮^(𝑘−1)(𝑡),𝑡𝐮^(𝑘)d𝑡. (11) Problem 𝑘 is again a TVLQR; note that 𝐱^(𝑘−1) and 𝐮^(𝑘−1) are the solutions of problem 𝑘−1 achieved at previous iteration. Solving problem 𝑘 yields 𝐱^(𝑘)(𝑡) and 𝐮^(𝑘)(𝑡), 𝑡∈[𝑡_𝑖,𝑡_𝑓].

Iterations continue until a certain convergence criterion is satisfied. In the present implementation of the algorithm, the convergence is reached when 𝐱^(𝑘)−𝐱^(𝑘−1)_∞=max_{𝑡∈[𝑡𝑖,𝑡𝑓]}𝑥_𝑗^(𝑘)(𝑡)−𝑥_𝑗^(𝑘−1)(𝑡),𝑗=1,…,𝑛≤𝜀, (12) where 𝜀 is a prescribed tolerance. That is, iterations terminate when the difference between each component of the state, evaluated for all times, changes by less than 𝜀 between two consecutive iterations.

4. Solution of the Time-Varying Linear Quadratic Regulator by the State Transition Matrix

With the approach sketched in Section 3, a fully nonlinear optimal control problem is reduced to a sequence of time-varying linear quadratic regulators. These can be solved a number of times to achieve an approximate solution of the original, nonlinear problem. This is done by exploiting the structure of the problem as well as its state transition matrix. This scheme differs from that implemented in [5, 6], and, in part, is described in [1].

Suppose that the following dynamics are given: 𝐱̇=𝐴(𝑡)𝐱+𝐵(𝑡)𝐮, (13) together with the quadratic objective function 1𝐽=2𝐱^𝑇𝑡_𝑓𝑆𝑡_𝑓𝐱𝑡_𝑓+12^𝑡𝑓_𝑡𝑖𝐱^𝑇𝑄(𝑡)𝐱+𝐮^𝑇𝑅(𝑡)𝐮d𝑡, (14) where 𝑄, 𝑅, and 𝑆 are positive semidefinite and positive definite time-varying matrices with appropriate dimensions, respectively. The Hamiltonian of this problem is 1𝐻=2𝐱^𝑇𝑄(𝑡)𝐱+𝐮^𝑇𝑅(𝑡)𝐮+𝝀^𝑇[𝐴(𝑡)𝐱+𝐵(𝑡)𝐮], (15) and the optimality conditions (5) read ̇(𝑡)(𝑡)̇(𝑡)𝐱=𝐴𝐱+𝐵𝐮,𝝀=−𝑄𝐱−𝐴^𝑇(𝑡)(𝑡)𝝀,0=𝑅𝐮+𝐵^𝑇(𝑡)𝝀. (16)(17)(18) From (18), it is possible to get 𝐮=−𝑅⁻¹(𝑡)𝐵^𝑇(𝑡)𝝀, (19) which can be substituted into (16)-(17) to yield 𝐱̇=𝐴(𝑡)𝐱−𝐵(𝑡)𝑅⁻¹(𝑡)𝐵^𝑇(𝑡)̇(𝑡)𝝀,𝝀=−𝑄𝐱−𝐴^𝑇(𝑡)𝝀. (20) In a compact form, (20) can be arranged as 𝐱̇𝝀̇=𝐴(𝑡)−𝐵(𝑡)𝑅⁻¹(𝑡)𝐵^𝑇(𝑡)−𝑄(𝑡)−𝐴^𝑇(𝑡)𝐱𝝀. (21) Since (21) is a system of linear differential equations, the exact solution can be written as 𝐱(𝑡)=𝜙_𝑥𝑥𝑡_𝑖,𝑡𝐱_𝑖+𝜙_𝑥𝜆𝑡_𝑖,𝑡𝝀_𝑖,𝝀(𝑡)=𝜙_𝜆𝑥𝑡_𝑖,𝑡𝐱_𝑖+𝜙_𝜆𝜆𝑡_𝑖,𝑡𝝀_𝑖, (22)(23) where the functions 𝜙_𝑥𝑥, 𝜙_𝑥𝜆, 𝜙_𝜆𝑥, and 𝜙_𝜆𝜆 are the components of the state transition matrix, which can be found by integrating the following dynamics: 𝜙̇_𝑥𝑥𝜙̇_𝑥𝜆𝜙̇_𝜆𝑥𝜙̇_𝜆𝜆=𝐴(𝑡)−𝐵(𝑡)𝑅⁻¹(𝑡)𝐵^𝑇(𝑡)−𝑄(𝑡)−𝐴^𝑇(𝑡)𝜙_𝑥𝑥𝜙_𝑥𝜆𝜙_𝜆𝑥𝜙_𝜆𝜆, (24) with the initial conditions 𝜙_𝑥𝑥𝑡_𝑖,𝑡_𝑖=𝜙_𝜆𝜆𝑡_𝑖,𝑡_𝑖=𝐼_𝑛×𝑛,𝜙_𝑥𝜆𝑡_𝑖,𝑡_𝑖=𝜙_𝜆𝑥𝑡_𝑖,𝑡_𝑖=0_𝑛×𝑛. (25)

If both 𝐱_𝑖 and 𝝀_𝑖 were given, it would be possible to compute 𝐱(𝑡) and 𝝀(𝑡) through (22)-(23), and therefore the optimal control function 𝐮(𝑡) with (19). The initial condition is assumed to be given, whereas the computation of 𝝀_𝑖 depends on the final condition, which, in the present algorithm, can be defined in three different ways.

4.1. Hard Constrained Problem

In a hard constrained problem (HCP), the value of the final state is fully specified, 𝐱(𝑡_𝑓)=𝐱_𝑓, and therefore, (14) does not account for 𝑆. The value of 𝝀_𝑖 can be found by writing (22) at final time 𝐱_𝑓=𝜙_𝑥𝑥𝑡_𝑖,𝑡_𝑓𝐱_𝑖+𝜙_𝑥𝜆𝑡_𝑖,𝑡_𝑓𝝀_𝑖 (26) and by solving for 𝝀_𝑖; that is, 𝝀_𝑖𝐱_𝑖,𝐱_𝑓,𝑡_𝑖,𝑡_𝑓=𝜙⁻¹_𝑥𝜆𝑡_𝑖,𝑡_𝑓𝐱_𝑓−𝜙_𝑥𝑥𝑡_𝑖,𝑡_𝑓𝐱_𝑖. (27)

4.2. Soft Constrained Problem

In a soft constrained problem (SCP), the final state is not specified, and thus 𝑆 in (14) is an 𝑛×𝑛 positive definite matrix. The transversality condition (7) sets a relation between the state and costate at final time 𝝀𝑡_𝑓=𝑆𝑡_𝑓𝐱𝑡_𝑓, (28) which can be used to find 𝝀_𝑖. This is done by writing (22)-(23) at final time and using (28) 𝐱𝑡_𝑓=𝜙_𝑥𝑥𝑡_𝑖,𝑡_𝑓𝐱_𝑖+𝜙_𝑥𝜆𝑡_𝑖,𝑡_𝑓𝝀_𝑖,𝑆𝑡_𝑓𝐱𝑡_𝑓=𝜙_𝜆𝑥𝑡_𝑖,𝑡_𝑓𝐱_𝑖+𝜙_𝜆𝜆𝑡_𝑖,𝑡_𝑓𝝀_𝑖. (29) Equations (29) represent a linear algebraic system of 2𝑛 equations in the 2𝑛 unknowns {𝐱(𝑡_𝑓),𝝀_𝑖}. The system can be solved by substitution to yield 𝝀_𝑖𝐱_𝑖,𝑡_𝑖,𝑡_𝑓=𝜙_𝜆𝜆𝑡_𝑖,𝑡_𝑓−𝑆𝑡_𝑓𝜙_𝑥𝜆𝑡_𝑖,𝑡_𝑓⁻¹×𝑆𝑡_𝑓𝜙_𝑥𝑥𝑡_𝑖,𝑡_𝑓−𝜙_𝜆𝑥𝑡_𝑖,𝑡_𝑓𝐱_𝑖. (30)

4.3. Mixed Constrained Problem

In a mixed constrained problem (MCP), some components of the final state are specified and some are not. Without any loss of generality, let the state be decomposed as 𝐱=(𝐲,𝐳), where 𝐲 are the 𝑝 known components at final time, 𝐲(𝑡_𝑓)=𝐲_𝑓, and 𝐳 are remaining 𝑛−𝑝 elements. The costate is decomposed accordingly as 𝝀=(𝝃,𝜼). With this formalism, 𝑆 in (14) is (𝑛−𝑝)×(𝑛−𝑝), and it is pre- and postmultiplied by 𝐳(𝑡_𝑓). The transversality condition (7) is 𝜼(𝑡_𝑓)=𝑆(𝑡_𝑓)𝐳(𝑡_𝑓).

The MCP is solved by partitioning the state transition matrix in a suitable form such that, at final time, (22)-(23) read 𝐲𝑡_𝑓𝐳𝑡_𝑓=𝜙_𝑦𝑦𝜙_𝑦𝑧𝜙_𝑧𝑦𝜙_𝑧𝑧𝐲_𝑖𝐳_𝑖+𝜙_𝑦𝜉𝜙_𝑦𝜂𝜙_𝑧𝜉𝜙_𝑧𝜂𝝃_𝑖𝜼_𝑖𝝃𝑡_𝑓,𝜼𝑡_𝑓=𝜙_𝜉𝑦𝜙_𝜉𝑧𝜙_𝜂𝑦𝜙_𝜂𝑧𝐲_𝑖𝐳_𝑖+𝜙_𝜉𝜉𝜙_𝜉𝜂𝜙_𝜂𝜉𝜙_𝜂𝜂𝝃_𝑖𝜼_𝑖, (31)(32) where the dependence of the state transition matrix components on 𝑡_𝑖, 𝑡_𝑓 is omitted for brevity. From the first row of (31), it is possible to get 𝝃_𝑖=𝜙⁻¹_𝑦𝜉𝐲_𝑓−𝜙_𝑦𝑦𝐲_𝑖−𝜙_𝑦𝑧𝐳_𝑖−𝜙⁻¹_𝑦𝜉𝜙_𝑦𝜂𝜼_𝑖, (33) which can be substituted in the second row of (31) to yield 𝐳𝑡_𝑓=𝜙_𝑧𝑦−𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝑦𝐲_𝑖+𝜙_𝑧𝑧−𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝑧𝐳_𝑖+𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝐲_𝑓+𝜙_𝑧𝜂−𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝜂𝜼_𝑖. (34) Equations (33)-(34), together with the transversality condition 𝜼(𝑡_𝑓)=𝑆(𝑡_𝑓)𝐳(𝑡_𝑓), can be substituted in the second row of (32) to compute the component of the initial costate 𝜼_𝑖𝐱_𝑖,𝐲_𝑓,𝑡_𝑖,𝑡_𝑓=𝜙_𝜂𝜂⁻¹𝐰𝐱_𝑖,𝐲_𝑓,𝑡_𝑖,𝑡_𝑓, (35) where 𝜙_𝜂𝜂=𝜙_𝜂𝜂−𝜙_𝜂𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝜂−𝑆𝜙_𝑧𝜂−𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝜂,𝐰𝐱_𝑖,𝐲_𝑓,𝑡_𝑖,𝑡_𝑓=𝑆𝜙_𝑧𝑦−𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝑦−𝜙_𝜂𝑦+𝜙_𝜂𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝑦𝐲_𝑖+𝑆𝜙_𝑧𝑧−𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝑧+𝜙_𝜂𝑧+𝜙_𝜂𝜉𝜙⁻¹_𝑦𝜉𝜙_𝑦𝑧𝐳_𝑖+𝑆𝜙_𝑧𝜉𝜙⁻¹_𝑦𝜉−𝜙_𝜂𝜉𝜙⁻¹_𝑦𝜉𝐲_𝑓. (36) Once 𝜼_𝑖 is know, the remaining part of the initial costate, 𝝃_𝑖, is computed through (33), and therefore, the full initial costate is obtained as a function of the initial condition, given final condition, initial and final time; that is,  𝝀_𝑖(𝐱_𝑖,𝐲_𝑓,𝑡_𝑖,𝑡_𝑓)=(𝝃_𝑖(𝐱_𝑖,𝐲_𝑓,𝑡_𝑖,𝑡_𝑓),𝜼_𝑖(𝐱_𝑖,𝐲_𝑓,𝑡_𝑖,𝑡_𝑓)).

5. Numerical Examples

Two simple problems with nonlinear dynamics are considered to apply the developed algorithm. These correspond to the controlled relative spacecraft motion and to the controlled two-body dynamics for low-thrust transfers.

5.1. Low-Thrust Rendezvous

This problem is taken from the literature where a solution is available, for comparison’s sake [10, 11]. Consider the planar, relative motion of two particles in a central gravity field expressed in a rotating frame with normalized units: the length unit is equal to the orbital radius, the time unit is such that the orbital period is 2𝜋, and the gravitational parameter is equal to 1. In these dynamics, the state, 𝐱=(𝑥₁,𝑥₂,𝑥₃,𝑥₄), represents the radial, tangential displacements (𝑥₁,𝑥₂) and the radial, tangential velocity deviations (𝑥₃,𝑥₄), respectively. The control, 𝐮=(𝑢₁,𝑢₂), is made up by the radial and tangential accelerations, respectively.

The equations of motion are 𝑥̇₁=𝑥₃,𝑥̇₂=𝑥₄,𝑥̇₃=2𝑥₄−1+𝑥₁1𝑟³−1+𝑢₁,𝑥̇₄=−2𝑥₃−𝑥₂1𝑟³−1+𝑢₂, (37) with 𝑟=(𝑥₁+1)²+𝑥²₂. The initial condition is 𝐱_𝑖=(0.2,0.2,0.1,0.1). Two different problems are solved to test the algorithm in both hard and soft constrained conditions.

Hard Constrained Rendezvous. The HCP consists in minimizing 1𝐽=2^𝑡𝑓_𝑡𝑖𝐮^𝑇𝐮d𝑡 (38) with the final, given condition 𝐱_𝑓=(0,0,0,0) and 𝑡_𝑖=0, 𝑡_𝑓=1.

Soft Constrained Rendezvous. The SCP considers the following objective function: 1𝐽=2𝐱^𝑇𝑡_𝑓𝑆𝐱𝑡_𝑓+12^𝑡𝑓_𝑡𝑖𝐮^𝑇𝐮d𝑡, (39) with 𝑆=diag(25,15,10,10), 𝑡_𝑖=0 and 𝑡_𝑓=1 (𝐱_𝑓 is free).

The differential equations (37) are factorized into the form of (8) as 𝑥̇₁𝑥̇₂𝑥̇₃𝑥̇₄=𝑓𝑥00100001₁,𝑥₂11+𝑥₁𝑥0020𝑓₁,𝑥₂−20_𝐴(𝐱)×𝑥₁𝑥₂𝑥₃𝑥₄+00001001_𝐵𝑢₁𝑢₂, (40) with 𝑓(𝑥₁,𝑥₂)=−1/[(𝑥₁+1)²+𝑥²₂]^3/2+1. Thus, the problem is put into the pseudo-LQR form (8)-(9) by defining 𝐴(𝐱) and 𝐵 as in (40) and by setting 𝑄=0_4×4 and 𝑅=𝐼_2×2.

The two problems have been solved with the developed method. Table 1 reports the details of the HCP and SCP, whose solutions are shown in Figures 1 and 2, respectively. In Table 1, 𝐽 is the objective function at the final iteration, “Iter” is the number of iterations, and the “CPU time” is the computational time (this refers to an Intel Core 2 Duo 2 GHz with 4 GB RAM running Mac OS X 10.6). The termination tolerance 𝜀 in (12) is 10⁻⁹. The optimal solutions found replicate those already known in the literature [10, 11], indicating the effectiveness of the developed method.

Table 1 

Rendezvous solutions details.

(a) 𝑥1 versus 𝑥2

(b) 𝑥3 versus 𝑥4

(c) 𝑢1 versus 𝑢2

Figure 1 

Hard constrained rendezvous.

(a) 𝑥1 versus 𝑥2

(b) 𝑥3 versus 𝑥4

(c) 𝑢1 versus 𝑢2

Figure 2 

Soft constrained rendezvous.

5.2. Low-Thrust Orbital Transfer

In this problem, the controlled, planar Keplerian motion of a spacecraft in polar coordinates is studied. The dynamics are written in scaled coordinates, where the length unit corresponds to the radius of the initial orbit, the time unit is such that its period is 2𝜋, and the gravitational parameter is 1. The state, 𝐱=(𝑥₁,𝑥₂,𝑥₃,𝑥₄), is made up by the radial distance from the attractor (𝑥₁), the phase angle (𝑥₂), the radial velocity (𝑥₃), and the transversal velocity (𝑥₄), whereas the control, 𝐮=(𝑢₁,𝑢₂), corresponds to the radial and transversal accelerations, respectively [12, 13]. The equations of motion are 𝑥̇₁=𝑥₃,𝑥̇₂=𝑥₄,𝑥̇₃=𝑥₁𝑥²₄−1𝑥²₁+𝑢₁,𝑥̇₄=−2𝑥₃𝑥₄𝑥₁+𝑢₂𝑥₁, (41) and the objective function is 1𝐽=2^𝑡𝑓_𝑡𝑖𝐮^𝑇𝐮d𝑡, (42) with 𝑡_𝑖=0 and 𝑡_𝑓=𝜋. The initial state corresponds to the conditions at the initial orbit; that is, 𝐱_𝑖=(1,0,0,1). Two different HCPs are solved, which correspond to the final states 𝐱_𝑓=(1.52,𝜋,0,1.52^−3/2) and 𝐱_𝑓=(1.52,1.5𝜋,0,1.52^−3/2), respectively. This setup mimics an Earth-Mars low-thrust transfer. The dynamics (41) and the objective function (42) are put in the form (8)-(9) by defining 𝑄=0_4×4, 𝑅=𝐼_2×2, and 𝐴(𝐱)=100100001𝑥−𝑥³₁00𝑥₁𝑥₄00−2𝑥₄𝑥₁0,𝐵(𝐱)=01000010𝑥₁. (43)

The two HCPs have been solved with the developed method. The solutions’ details are reported in Table 2, whose columns have the same meaning as in Table 1. It can be seen that more iterations and an increased computational burden are required to solve this problem. The solution with 𝑥_2,𝑓=1.5𝜋 is reported in Figure 3.

Table 2 

Earth-Mars transfer details.

(a) Transfer trajectory

(b) Control profile

Figure 3 

Orbital transfer with 𝑥2,𝑓=1.5𝜋.

6. Conclusion

In this paper, an approximated method to solve nonlinear optimal control problems has been presented, with applications to sample cases in astrodynamics. With this method, the nonlinear dynamics and objective function are factorized in a pseudolinear and quadratic-like forms, which are similar to those used in the state-dependent Riccati equation approach. Once in this form, a number of time-varying linear quadratic regulator problems are solved. A state transition matrix approach is used to deal with the time-varying linear quadratic regulators. The results show the effectiveness of the method, which can be used to either have suboptimal solutions or to provide initial solutions to more accurate optimizers.