ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia Detailed study of low-thrust trajectory optimization problem in Russia begins in the end of 1950 th. The monographs [1,2] summarize results obtained to the early 1970 th. The wide range of optimization techniques was studied, including parametric optimization, direct, indirect and combined optimization algorithms. Most popular techniques during these years were gradient descent versions (both in parametric space and space of initial co-states), modified Newton algorithms, linearization and asymptotic technique [1-4]. First studies in this field were carried out under strong impact of the proposed manned mission to Mars. Feasibility study of this mission was carried at the Rocket Space Corporation Energia since The nuclear electric propulsion was considered as main propulsion candidate for this mission since 1960 ( see the page background ). Later scientific planetary missions and problem of spacecraft delivery into a high-altitude orbits become to be considered as a main motive of the optimization technique development. Somewhat separately the solar sail trajectory optimization developed [1, 2, 8]. At present, there are at least three research teams in Russia studying low-thrust trajectory optimization problems: Keldysh Institute of Applied Mathematics (KIAM), Russian Academy of Sciences Team leader is Academician T.M. Eneev. The members of team are G.B. Efimov, R.Z. Akhmetshin, G.S. Zaslavsky, V.V. Ivashkin, V. Zharov, and others. Moscow Aviation Institute (MAI) and Khrunichev State Research and Production Space Center (KSR&PSC) Team leader is Prof. M.S. Konstantinov, the members of team are G.G. Fedotov, V.G. Petukhov, and others. Samara State Aerospace University (SSAU) Team leader is Prof. V.V. Salmin, the members of team are V.V. Vasilyev, V.V. Yurin, S.A. Ishkov. In addition to above mentioned teams, there are researchers from the Moscow State University, Space Research Institute, and others institutes and companies who are working actively in the low-thrust optimization. The theme Low-thrust trajectory optimization in Russia has extremely wide range. It should be referred to the half-century history. Therefore this poster presents only some recent studies which are known to author. LOW-THRUST TRAJECTORY OPTIMIZATION IN RUSSIA Vyacheslav PETUKHOV * * Khrunichev State Research and Production Space Center Novozavodskaja st., 18, Moscow, , Russia.

1. CONSIDERED OPTIMAL CONTROL PROBLEMS (OCP) Two conventional mathematical models of low-thrust engine are under consideration. There are variable and constant specific impulse models. The jet power is fixed or it is given function of spacecraft position and time within frame of variable specific impulse problem. There are not any other restrictions on the thrust and specific impulse values. These assumptions lead to the power-limited problem (LP-problem). The exhaust velocity is fixed or it is given function of spacecraft position and time within frame of constant specific impulse problem. This assumption leads to the constant exhaust velocity problem (CEV-problem). There are two principally different CEV-problems, namely the minimum-time and minimum-propellant problems. In case of absence an external restrictions the minimum-time problem leads to the continuous thrusting during whole transfer. On the contrary, the minimum-propellant problem leads to the optimal trajectories having thrust switching. Conventional indirect approach to solve these problems supposes using calculus of variations or maximum principle. All these problems are reduced to the boundary value problems using maximum principle. In the simplest case there are two-point boundary value problems (TPBVP). Typical LP-problem and minimum-time CEV-problem lead to the smooth dependency of TPBVP residual vector versus unknown TPBVP parameters. On the contrary, this dependency is essential non-smooth within the minimum-propellant CEV-problem. Last fact is main reason of numerical instability of most numerical algorithm using indirect approach, derivatives of residuals with respect to unknown TPBVP parameters or any its approximations. Therefore, LP- and minimum-propellant CEV-problems should be solved using essential different numerical solvers. The expertise of low-thrust optimization in Russia confirms the robustness and efficiency of the continuation (and, in some cases, modified-Newton and shooting) techniques for LP-problem and minimum-time CEV-problem solving. The minimum-propellant problems can be solved using these algorithms also. But their convergence domain usually becomes extremely small in this case. More robust techniques are required in this case. As a rule, there are algorithms using direct or combined indirect-direct approaches. 2 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Purpose: To minimize for a dynamical system obeying the differential equations and having following boundary conditions: t = 0: x(0) =x 0, dx(0)/dt = v 0, m(0) = m 0, t = T: x(T) = x f, dx(T)/dt = v f.. Here T is final time, a is thrust acceleration, x and v are vectors of spacecraft position and velocity respectively, is gravity potential. Let apply Pontryagins maximum principle to reduce this OCP to the TPBVP. Hamiltonian of this dynamical system is H = -a T a/2 + p x T v + p v T x + p v T a. So, optimal control is a = p v and equations of optimal motion become following: The boundary conditions for rendezvous mission has form: t = 0: x(0) =x 0, dx(0)/dt = v 0, t = T: x(T) = x f, dx(T)/dt = v f. In fact, it is necessary to solve equation f(z) = 0, where is vector of residuals, is vector of unknown TPBVP parameters, p x = -dp v /dt Typical LP-Problem 2. MATHEMATICAL STATEMENT OF LOW-THRUST OPTIMIZATION PROBLEM: INDIRECT APPROACH Here conventional mathematical formulation of main low-thrust trajectory optimization problems is presented. Pontryagins maximum principle is used to reduce OCP to the boundary value problem. For the sake of simplicity, the minimal set of constraints is considered. It is considered the simplest case of boundary conditions corresponding to rendezvous problem. 3 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

2.2. Typical CEV-Problem with Thrust Switching 2.3. Typical Minimum-Time CEV-Problem (continuous thrust) Purpose: To minimize. As result, 1,, so o.d.e. for m and p m can be excluded from consideration. The vectors of residuals and TPBVP parameters become following, where is Hamiltonian. Purpose: To minimize for a dynamical system obeying the differential equations and having following boundary conditions: t = 0: x(0) =x 0, dx(0)/dt = v 0, m(0) = m 0, t = T: x(T) = x f, dx(T)/dt = v f. Here is step-like thrusting function, P – thrust, w – exhaust velocity, m – spacecraft mass. t = 0: x(0) =x 0, dx(0)/dt = v 0, m(0) = m 0, t = T: x(T) = x f, dx(T)/dt = v f, p m (T) = 0, where step-like thrusting function and switching function In fact, it is necessary to solve equation f(z) = 0, where is vector of residuals, and Pontryagins maximum principle reduces this OCP to the following TPBVP: 4 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

3. KELDYSH INSTITUTE OF APPLIED MATHEMATICS KIAM team developed algorithm of transporting trajectories in the early 1960 th [1, 2, 9, 10]. The essence of this algorithm is linearization of differential equations of optimal motion in the vicinity of a conic-patched trajectory (transporting trajectories). Typically, the conic- patched problem was calculated as solution of the Lambert problem. Within frame of this technique, the two-burn transporting trajectory only was considered. The linear TPBVP was solved in closed form. This method brings to the satisfactory results if trajectory angular distance not exceeds ~220°. Otherwise, the relative error becomes too big. More power technique was developed later [11]. The new approach uses following 3-phases algorithm: 1 st phase. The multi-burn conic-patched trajectory is optimized. 2 nd phase. The TPBVP is linearized in the vicinity of obtained conic-patched trajectory. This linear TPBVP is solved using conventional numerical algorithms, e.g. sweep technique. 3 rd phase. Solution of the linearized TPBVP is used as initial approximation for complete LP-problem. The modified Newton technique is used to solve the TPBVP of LP-problem Algorithms for LP-Problem KIAM team begins to study low-thrust trajectory optimization problem in early 1960 th. This team carried out feasibility studies of numerous missions. These studies required developing of mathematical tools and optimization technique. Many studies cannot be presented here. There are development and applying the computer algebra tools for trajectory optimization problem, analytical trajectory optimization, parametric optimization using local-optimal thrust steering. KIAM team is carrying out studies of geocentric and interplanetary low-thrust trajectories, including multi-target trajectories to asteroids and comets. Since 1994 many feasibility studies are carried out by joint team of KIAM and MAI. 5 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

KIAM team uses combined indirect-direct approach to solve CEV-problem. The program of low-thrust direction, which is not restricted, was defined using the necessary optimal conditions of maximum principle. The program of thrust value was determined directly as an ordinary function of few parameters. The optimal conditions for these parameters combined with terminal conditions forms the extended TPBVP which is easier than the traditional one. Solution of corresponding LP-problem provides initial conditions for co-states variables and additional parameters, which fix the program of thrust value. The modified Newton algorithm is used for the TPBVP solving. These ideas applying leads to excluding co-state variable p m and switching function s from the consideration. Instead, the step-like thrusting function is defined explicitly as function of trajectory arc number i, e.g. = [1 + (-1) i ] / 2. Time points of thrust switching are considered as an external parameters which should be optimized. Therefore, following constrained minimization problem is under consideration: Here. Of course, corresponding differential equations for state and co-state variables (except for p m ) should hold true. This NLP-problem reduces to the system of nonlinear equations using Lagrange multipliers. Last system is solved with respect to z and t i, i = 1,…, n-1, using modified Newton algorithm. R.Z. Akhmetshin developed in analogous algorithm for optimizing multirevolutional low-thrust transfer from an elliptical orbit into GEO. His paper will be published in Cosmic Research in P t LP-problem CEV-problem t0t0 t1t1 t2t2 t3t3 t4t4 t5t5 Arc 1 Arc 2 Arc 3 Arc 4 Arc Algorithms for CEV-Problem with Thrust Switching 6 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

4. MOSCOW AVIATION INSTITUTE AND KHRUNICHEV STATE RESEARCH AND PRODUCTION SPACE CENTER Low-thrust trajectory optimization problems were studied in MAI since 1960 th. Some early results are presented in [4]. Both direct and indirect approaches were considered. Later direct approach was rejected in fact, because complexity of solving large NLP-problems. Using indirect approach, there were derived robust algorithms and software for solving both LP- and CEV-problems. Special techniques and software were developed for very multirevolutional transfers. At present, in the connection with computational productivity increasing and new numerical techniques developing, direct and, especially, combined direct-indirect approaches are attracted attention again. MAI team is carrying out studies of geocentric and interplanetary low-thrust trajectories. Carried out studies includes optimization of flyby, rendezvous, gravity-assisted trajectories, multi-target missions, trajectories to Moon and libration points, escape trajectories, transfers between orbits, nuclear and solar electric propulsion missions. Recently, there were carried out the optimization of Phobos-Soil mission, optimization of insertion into GEO using electric propulsion for advanced Russian and foreign spacecrafts. Since 1994 many feasibility studies are carried out by joint team of KIAM and MAI Algorithms for LP-Problem At present, the continuation (homotopy) technique is used for LP-problem solving. Corresponding algorithms and software were developed by author [12-14]. Purpose of new continuation method Regularization of numerical trajectory optimization, i.e. elimination (if possible) the methodical deffects of numerical optimization. Particularly, the was stated and solved problem of trajectory optimization using trivial initial approximation (the coasting along the initial orbit for example). Applied problems of trajectory optimization under consideration 1. Planetary low thrust trajectory optimization; 2. Lunar low thrust trajectory optimization within the frame of restricted problem of three bodies. 7 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Problem: to solve non-linear system (1) with respect to vector z Let z 0 - initial approximation of solution. Then (2) where b - residuals when z = z 0. Let consider immersion z( ), where is a scalar parameter, and equation (3) with respect to z( ). Obviously, z(1) is solution of eq. (1). Let differentiate eq. (2) on and resolve it with respect to dz/d : (4) Just after integrating eq. (4) from 0 to 1 we have solution of eq. (1). Equation (4) is the differential equation of continuation algorithm (the formal reduction of non-linear system (1) into initial value problem (4)). Differential equations (4) are integrated using a high-order Runge-Kutta method with adaptive step size control Continuation (Homotopy) Technique Application to Optimal Control Problem (OCP) Optimal motion equations (after applying the maximum principle): Boundary conditions (an example): Boundary value problem parameters and residuals: Sensitivity matrix: Associated system of optimal motion o.d.e. and perturbation equations for residuals and sensitivity matrix calculation: Extended initial conditions: 8 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

OCP-Solver Based on Continuation Algorithm Reduction of optimal control problem to the boundary value problem by maximum principle Initial approximation z 0 Initial residuals b calculation by optimal motion o.d.e. integrating for given guess value z 0 of boundary value problem parameters Associated integrating of optimal motion equations and perturbation equations for current z( ) to calculate current residuals f(z, ) and sensitivity matrix f z (z, ) Continuation methods o.d.e. integrating with respect to from 0 to 1 Integrating of optimal motion equations for current z( ) to calculate current residuals f(z, ) and for pertubed z( ) to calculate f z (z, ) by finite- difference Solution z(1) CONTINUATION METHOD 1st version of o.d.e. right- hand side calculation 2nd version of o.d.e. right-hand side calculation 9 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Continuation with Respect to Gravity Parameter (1/3) Occasional reasons of continuation algorithm failure: sensitivity matrix degeneration (bifurcation of optimal solutions) Mostly bifurcations of optimal planetary trajectories are connected with different number of complete orbits. If angular distance will remain constant during continuation, the continuation way in the parametric space will not cross boundaries of different kinds of optimal trajectories. So, the sensitivity matrix will not degenerate. The purpose of the technique modification - to fix angular distance of transfer during continuation Earth-to-Mars, rendezvous, launch date June 1, 2000, V = 0 m/s, T=300 days 1 - coast trajectory ( 1 = 0) intermediate trajectories (0 < 2 < 3 < 4 < 1) 5 - final (optimal) trajectory ( 5 = 1) Sequence of trajectory calculations using basic continuation method Sequence of trajectory calculations using continuation with respect to gravity parameter ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Let x 0 (0), x 0 (T) - departure planet position when t=0 and t=T; x k - target planet position when t=T. Let suppose primary gravity parameter to be linear function of, and let choose initial value of this gravity parameter 0 using following condition: 1) angular distances of transfer are equal when =0 and =1; 2) When =1 primary gravity parameter equals to its real value (1 for dimensionless equations) The initial approximation is SC coast motion along departure planet orbit. Let the initial true anomaly equals to 0 at the start point S, and the final one equals to k = 0 + at the final point K ( is angle between x 0 and projection of x k into the initial orbit plane). The solution of Kepler equation gives corresponding values of mean anomalies M 0 and M k (M=E-e sinE, where E=2 arctg{[(1-e)/(1+e)] 0.5 tg( /2)} is eccentric anomaly). Mean anomaly is linear function of time at the keplerian orbit: M=M 0 +n (t-t 0 ), where n=( 0 /a 3 ) 0.5 is mean motion. Therefore, the condition of angular distance invariance is M k +2 N rev =nT+M 0, where N rev is number of complete orbits. So initial value of the primary gravity parameter is 0 =[( M k +2 N rev - M 0 )/T] 2 a 3, and current one is ( )= 0 +(1- 0 ). The shape and size of orbits should be invariance with respect to, therefore v(t, )= ( ) 0.5 v(t, 1) Continuation with Respect to Gravity Parameter (2/3) 11 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Continuation with Respect to Gravity Parameter (3/3) z = (p v (0), dp v (0)/dt) T = b = f(z 0 ) Equations of motion: Boundary conditions: Residuals: Boundary value problem parameters: Equation of continuation method: where 12 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

EPOCH Application for Planetary Low-Thrust Missions Optimization (LP-problem) V. Petukhov, Moscow Aviation Institute, Khrunichev Space Research and Production Space Center Application EPOCH (Electric Propulsion trajectories - Optimal CHoice) optimizes interplanetary low-thrust trajectories using continuation algorithm. The variable specific impulse power-limited problem is considered. Guess values for initial co-states are not required. EPOCH allows to optimize rendezvous and flyby trajectories to planets, asteroids, comets, and user- defined objects. Both basic continuation algorithm and continuation w.r.t. gravity parameter are implemented. Both NEP and SEP missions are supported. Polynomial approximation of the solar array efficiency can be used. Asteroid and comet database is supported. EPOCH uses standard MS Windows user interface. 3D-graphics is supported. Mouse-based rotation facility is implemented for ease of handling the 3D-views Software 13 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Examples of EPOCH Numerical Examples Optimal Trajectories to Mercury and Near-Earth Asteroids (CPU time for AMD Athlon XP 1600+) 30+ orbits, CPU time 60 s 3+ orbits, CPU time 1.2 s 3+ orbits, CPU time 1.1 s 3+ orbits, CPU time 1.2 s 14 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Optimal Control in the Restricted Problem of Three Bodies: Earth-Moon L 2 Rendezvous Using Moon Gravity Assisted Maneuver Thrust acceleration, mm/s Т, days 95 Moon orbit Final Moon position Initial Moon position Initial L 2 position Final L 2 position Initial orbit Gravity assisted maneuver Thrust acceleration, mm/s Т, days orbits7.5 orbits 15 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

4.2. Algorithms and Software for CEV-Problem Two essentially different CEV-problems are supported by algorithms and software developed by the team members. There are conventional finite-thrust problem and averaged low-thrust problem for very multirevolutional trajectories. All algorithms uses indirect approach. G.G. Fedotov uses modified Newton technique in combination with line search. He developed a quite robust software allowing to optimize both interplanetary and interorbital trajectories, including gravity-assisted trajectories, missions using combined chemical/electric propulsion or multi-mode electric propulsion, etc. [14, 19-21, 23, 24]. His software used, for example, for Phobos-Soil trajectory optimization, feasibility study of numerous planetary missions, solar thermal propulsion mission optimization. Yu.A. Zakharov ( 1996) developed an efficiency technique to continue a conic-patched trajectory into optimal finite-thrust one [5, 6]. This algorithm was used both for interplanetary and interorbital trajectory optimization also, including multi-target trajectory to asteroids and comets. M.S. Konstantinov developed non-iterative procedure for approximate optimization of multirevolutional minimum-time problem [18, 22]. The essence of this method is using a simplified model problem which has unique solution. For the model problem the optimal feedback control is found. This feedback control is applied to the real CEV-problem. The model problem becomes equivalent to real CEV-problem on the near-circular orbit. Therefore the feedback control provides more and more precise guidance when spacecraft approaches to the final circular orbit. As result, the boundary conditions are satisfied exactly and optimality conditions are satisfied approximately. The difference between exact and approximate optimal solution typically does not exceed few percents. This technique was used to study low-thrust missions to GEO. V.G. Petukhov developed robust continuation algorithm for multirevolutional minimum-time problem [13, 17]. 3D transfers between arbitrary elliptical orbits are supported. The factored secant update algorithm is used to solve minimum-propellant problem. Developed software was used for many missions analysis, e.g. [15, 16]. This software allowed to design interpolation algorithm for quick and precise mission performance estimation. There were analyzed: the structure of optimal thrust steering for different boundary conditions; coasting arcs evolution when transfer duration is increased; the bifurcation of optimal solutions. 16 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

PL Application for Planetary Low-Thrust Missions Optimization (CEV-problem) G. Fedotov, Moscow Aviation Institute G.G. Fedotov (MAI) developed an robust software to solve CEV-problem. This application evolves for a twenty years. It uses an indirect approach and modified Newton technique in combination with line search to solve equations for residuals. The equations of motion are represented in the spherical coordinates. Both minimum-time and minimum-propellant problems are supported. The power and thrust throttling performance can be taken into account, so NEP, SEP, and conventional propulsion missions can be considered. Various kinds of boundary conditions are supported allowing to consider flyby, rendezvous, gravity-assisted trajectories, and transfer between orbits. At present, this application runs under MS-DOS and has a good user interface and graphical environment. This application can be considered as Russian analogue of JPLs VARITOP and SEPTOP. Of course, this application inherits common shortcomings of Newton-like algorithms, namely a relatively small convergence region, instability if trajectory structure (i.e. sequence of burn and coast arcs) is changed during solving process Software Based on Modified Newton algorithm 17 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Examples of PL Numerical Results 18 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia Optimal low-thrust transfer from inclined LEO into GEO: 53 orbits, 54 thrust arcs Optimal low-thrust trajectory to Mercury using Venus swingby

Yu.A. Zakharov (MAI) used a robust discrete continuation algorithm to solve CEV-problem. The essence of this technique is continuation from optimal conic-patched solution into the finite-thrust CEV-problem. Corresponding TPBVP is solved using discrete version of predictor-corrector continuation algorithm. Continuation parameter is inverse initial thrust acceleration. This continuation parameter equals zero on the conic-patched solution. At the first phase of TPBVP solving the optimal conic-patched solution and Lawdens primer-vector time history is calculated. Then the step on continuation parameter is carried out (predictor). TPBVP having new value of continuation parameter is solved using modified Newton (Broyden-update) technique (corrector). The guess values for the Broyden-update algorithm is TPBVP parameters getting by linear or quadratic extrapolation of TPBVP parameters from previous values of continuation parameter. The manual support only of stepsize control wrt. continuation parameter was realized. This technique was applied to the both planetary and interorbital missions optimization. It was extended to the multi-target trajectory to the asteroids and comets, including intermediate flybys and rendezvous. Yu.A. Zakharov developed an effective DISIPO application based on this technique. This algorithm is presented in [5, 6] Discrete Continuation from Optimal Conic-Patched Solution for CEV-Problem Algorithm Description [6] 1 st phase: Conic-patched trajectory optimization. NLP-problem: V min, x(t 0 ) = x 0 (t 0, p 0 ), v(t 0 ) = v 0 (t 0, p 0 ), x(t 0 +T) = x f (t 0 +T, p 0 ), v(t 0 +T) = v f (t 0 +T, p 0 ). Numerical technique: Quasi-Newton algorithms (Davidon-Fletcher-Powell update) or specialized algorithm for Lambert problem if two-burn transfer is considered and boundary points and transfer duration are fixed. 19 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

A, B, C, D, E, F – constants of integrations, U, V – radial and circumferential components of velocity, p – latus rectum, e – eccentricity, - true anomaly, t p – time of pericenter crossing. Combining (1) and (2) we get linear system with respect to constants of integration on each conic arc. Values p v (t 0 ), p v (t 0 ) are calculated using (2) after computation of integration constants on the first conic arc. The next step is transformation these vectors to the geocentric fixed reference frame using corresponding transformation matrix. Primer-vector p v and its negative derivative p x has following solution on the Keplerian trajectory arc (D.F. Lawden): where p v and p x are presented in the orbital reference frame (radial, circumference, and normal orths), 2 nd phase: Calculation of the primer-vector initial value on the conic-patched trajectory. At each point of velocity increment p v (t i ) = V i /| V i |. (1) (2) 3 rd phase: Conic-patched trajectory optimality checking. The condition |p v (t)| 1 is checked. If this inequality is violated, the conic-patched solution is not optimal. This means the optimal conic- patched trajectory has boundary coast arcs or it should have more burns. In this case it should be changed corresponding parameters of transfer, e.g. transfer duration, boundary conditions or number of burns. 4 th phase: Continuation from the conic-patched trajectory to the finite-thrust trajectory. Let z( ) = (p v (t 0 ) T, p x (t 0 ) T ) T – unknown TPBVP parameters, = 1/a 0 – continuation parameter (inverse value of initial thrust acceleration), f(z, ) - TPBVP residuals vector. Initial value of TPBVP parameters z 0 = z(0) was computed in the 2 nd phase. Let suppose i = i-1 +, where is small enough. Let consider TPBVP for thrust acceleration value a 0 = 1/ i. Let put to z( i-1 ) the guess values of TPBVP parameters. Then Broyden update algorithm is applied to solve this TPBVP. Last step should be repeated until required thrust acceleration value will be achieved. 20 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Smooth Continuation from LP-problem into CEV-problem 21 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia Author developed technique of smooth continuation from the LP-problem into the CEV-problem. The equations of optimal motion is modified such a way that they become depending on continuation parameter, these equation is coincided with equations of LP-problem at start value = 0 and they are coincided with equations of CEV-problem at final value = 1. Unfortunately, this technique works proper only if trajectory structure (i.e. sequence of thrust and coast arcs) does not changed during optimization. Earth Jupiter Thrust acceleration, mm/s 2 t, days LP-problem Earth Jupiter Thrust Coast t, days CEV-problem In-plane angle Out-of-plane angle t, days

Techniques and Software for Multirevolutional Transfers Optimization [13, 17] Author uses equinoctial orbital elements: Here p, e,,, i, are Keplerian elements, - primary gravity parameter. It is considered conventional CEV-problem without any constraints on thrust direction. Maximum principle reduces the problem into TPBVP. The numerical averaging over the orbital period is used for computational consumption reducing and numerical stability increasing. A numerous versions of boundary conditions were considered. The continuation (homotopy) procedure was used to solve minimum-time problem (see 4.1.1). The simple typical guess values: p h0 = 1, p ex0 = p ey0 = p ix0 = p iy0 = 0 (initial values of co-state variables), T = 1 (dimensionless orbital period referred to the initial orbit) as a rule provides stable convergence of optimization. Of course, initial values of co-states from the OCP solution having close boundary conditions provides improved convergence. Solver of minimum-propellant problem uses minimum-time solution as initial approximation. The factored secant update algorithm is used for minimum-propellant problem. Both techniques demonstrates theirs robustness and efficiency and there were used for a numerous applied problems. E-ProTO software was developed based on these algorithms.,,,,, 22 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

E-ProTO Application for Multirevolutional Geocentric Low-Thrust Trajectories Optimization (CEV-problem) and ASTROLABE Visualization Tool V. Petukhov, Moscow Aviation Institute, Khrunichev State Research and Production Space Center Application E-ProTO (Electric Propulsion Trajectory Optimization) optimizes and simulates multirevolutional geocentric trajectories. The constant ejection velocity problem is considered. Two techniques of TPBVP solving are implemented: the continuation algorithm is used for the minimum-time problem and the factored secant update algorithm is used for minimum-propellant problem. Averaged equations of motion are used during optimization. E-ProTo allows optimize both minimum-time and minimum-propellant problems. J 2 can be taken into account during optimization. E-ProTO implements trajectory simulation using optimal control program and various perturbations including primarys gravity terms of arbitrary degree and order, moon/sun and planets gravity, air drag, earth precession and nutation, coasting during eclipses. 23 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

E-ProTO Numerical Results (1/4) Evolution of Optimal Trajectory Structure (transfer from inclined HEO into GEO) / /0-2/2-1/ /0-2/2-1a/1-2/2-1/ /0-2/2-1a/1-2/ /2-1a/1-2/ p/1-2/2-1a/1-2/2-1p/ T, days Trajectory Diagram Arcs SequenceV req, m/s (N thrust arcs / N coast arcs) ° t, days Coast Number of thrust arcs per orbit (p if unique thrust arc is placed at perigee and a if unique thrust arc is placed at apogee) Number of coast arcs per orbit 24 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

The trajectory optimization of SMART-1 spacecraft is considered. The initial orbit has following parameters: perigee radius km, apogee radius km, inclination 6.655°, right ascension of ascending node °, and argument of perigee °. The parameters of the final orbit are following: apogee radius km, inclination 5.49°, right ascension of ascending node 0.25°, argument of perigee 79.66°. The final perigee radius is free. Initial spacecraft mass equals to kg, transfer duration is fixed and equals to 284 days. The electric propulsion provides thrust mN and exhaust velocity m/s. Optimization E-ProTO led to the following results: final spacecraft mass equals to kg and final perigee radius equals to ~24500 km. All orbits but one have single thrust arc. SMART-1 Trajectory Optimization E-ProTO Numerical Results (2/4) 25 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

Minimum-time transfer between non-coplanar circular orbits E-ProTO Numerical Results (3/4) Two kind of optimal trajectories were found. Eccentricity remains close to zero on the C-trajectories (optimal trajectory of 1 st kind ) and averaged eccentricity has a finite maximum on the E-trajectories (optimal trajectories of 2 nd kind). The critical inclination existence was found. E-trajectories exists if angle between initial and final orbits exceeds ~47.3. Analogous kinds of trajectories were found in case of transfer from an inclined elliptical orbit into a circular one. 26 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia h 0 = km h 0 = km h 0 = km E-solution C-solution Bifurcation point i ° Minimum-time transfer from inclined circular orbit into GEO (thrust N, specific impulse 1500 s. initial spacecraft mass 1320 kg)

Minimum-time transfer from Inclined Elliptical Orbit into GEO (thrust N, specific impulse 1500 s. initial spacecraft mass 1320 kg) Isosurfaces of transfer duration vs. initial perigee/apogee altitude and inclination T = 120 days T = 150 days T = 180 days h 0 i0i0 T = 90 days E-ProTO Numerical Results (4/4) 27 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

5. SAMARA STATE AEROSPACE UNIVERSITY Monographs [7] presents main results of low-thrust trajectory optimization carried out by SSAU team. Typical approach of their study can be described by following pattern [7]: 1) Mathematical statement of problem is carried out in the most general form; 2) Equations of motion and constraints are simplified using physical features of problem; 3) Simplified problem is optimized; 4) Proximity of computed solution to optimal one is estimated. In practice, 2 nd stage means using, for example, equations of near-circular motion or averaged equation of motion using pre-defined (in parametric form) thrust steering. Parametric representation of thrust steering can be fulfilled using local optimization (maximization of derivatives in the local point) or other physical reasons. During 3 rd stage, OCP is interpreted as optimization of slow controls, i.e. parameters of control law. Maximum principle is applied for slow controls optimization. Such simplified approach allows to use various numerical method. As a rule modified Newton method is used. Main feature of this studies is using of OCP sufficient conditions [25] to estimate proximity of computed solution to the optimal one. SSAU team studied problems of optimal associated attitude and trajectory control, optimization of multirevolutional transfers between non-coplanar near-circular or coplanar elliptical orbits, including problems with gravity and drag perturbation forces. It was proposed discrete mathematical representation of orbital low-thrust motion using, in fact, averaging technique. This representation was used to solve OCP using dynamic programming algorithm. 28 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

6. FUTURE VISION Russian expertise shows viability and efficiency of indirect approach to the low-thrust trajectory optimization. Novel indirect optimization techniques (most of all, homotopy techniques) demonstrates robustness at least comparable with direct approaches. These new methods keep main common advantages of indirect approach, namely they require relatively low computational resources and they satisfy exactly to necessary optimal conditions. Unfortunately, these methods are robust for smooth problems only, namely for LP-problem and minimum-time CEV-problem. It remains actual the development of robust indirect technique for minimum-propellant CEV-problem. In authors opinion, main obstacle to development such technique is fluctuating of trajectory structure (i.e. sequence of thrust and coast arcs) during optimization. This fluctuating leads to the discontinues of sensitivity matrix f z and, therefore, to the numerical failure. Therefore, novel technique should provide freezing of trajectory structure. It can be achieved, for example, by using combined indirect-direct method similar to presented in section 3.2 but applying continuation technique instead modified Newton one. Of course, necessary optimal conditions should be checked after TPBVP solving. If they are not satisfied, the pre-defined trajectory structure should be changed. Another way, may be, is regularization of minimum-propellant CEV-problem, i.e. problem smoothing. For example, change of independent variable t s obeying to differential expression dt = s 2 ds, where s is switching function, removes discontinuities from the equations of optimal motion. Another very promising way of indirect optimization is connected with multiply shooting method. The interval optimization methods are very interesting because they can support non-smooth problems. In authors opinion, direct and indirect approaches are supplemented each other because they typically have opposite shortcomings and advantages. Optimization of low-thrust trajectories to Moon using SOCS is brilliant example of direct approach successful applying [26]. Without doubt, it would be a very hard work to repeat these results using indirect approach. But the computation of unique optimal low- thrust trajectory to Moon using SOCS consumes ~1.5 days of CPU time [26]! In the other hand, calculating of optimal 30-orbits trajectory to Mercury using EPOCH requires only 60 s of CPU time, GTO-to-GEO low-thrust multirevolutional trajectory optimization using E-ProTO requires only a minutes for very rough guess values and a seconds otherwise. Feasibility study of advanced missions, as a rule, requires a lot of trajectories calculation. For example, there were performed optimization of ~1000 multirevolutional trajectories using E-ProTO to analyze continuous low-thrust mission from an inclined elliptical orbits into GEO (see isosurfaces in section ). All calculations consume only ~6 hours of CPU time. In case of direct approach about 3 years would be required to do the same work. We rarely have such time resources. 29 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

REFERENCES Monographs on Low-Thrust Trajectory Optimization (in Russian) 1. Grodzovsky, G.L., Yu.N. Ivanov, V.V. Tokarev. Space Flight Mechanics Using Low-Thrust. Moscow, Nauka, 1966, 680 pp. 2. Grodzovsky, G.L., Yu.N. Ivanov, V.V. Tokarev. Space Flight Mechanics. Optimization problems. Moscow, Nauka, 1975, 703 pp. 3. Lebedev, V.N. Computation of Low-Thrust Spacecraft Motion. Moscow, Comput. Center of USSR Academy of Sciences, 1968, 108 pp. 4. Konstantinov, M.S. Mathematical Programming Techniques for Flying Vehicles Design. Moscow, Mashinostroeniye, 1975, 164 pp. 5. Zakharov, Yu.A. Interorbital Spacecrafts Design. Moscow, Mashinostroeniye, 1984, 174 pp. 6. Grishin, S.D., Yu.A. Zakharov, V.K. Odelevsky. Design of Spacecrafts with Low-Thrust Engines. Moscow, Mashinostroeniye, 1990, 224 pp. 7. Salmin, V.V. Low-Thrust Transfers Optimization. Moscow, Mashinostroeniye, 1987, 208 pp. 8. Polyakhova, E.N. Space Flight with Solar Sail: Problems and Prospect. Moscow, Nauka, 1986, 304 pp. Other Publications 9. Beletsky, V.V., V.A. Egorov. Interplanetary Flights Using Fixed-Power Engines. Cosmic Research, 2, No. 3, Beletsky, V.V., V.A. Egorov, V.G. Ershov. Analysis of Interplanetary Trajectories Using Fixed-Power Engines. Cosmic Research, 3, No. 3, Glazkov, A.I. Solving of Two-Point Boundary Value Problem for Low-Thrust Trajectories Optimization. Cosmic Research, 22, No. 6, Petukhov, V. One Numerical Method to Calculate Optimal Power-Limited Trajectories. Moscow, IEPC , Petukhov, V. Low-Thrust Trajectory Optimization. Presentation at the seminar on Space Flight Mechanics, Control, and Information Science of Space Research Institute (IKI), Moscow, June 14, 2000 ( 30 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia

14. Eneev T.M., M.S. Konstantinov, V.A. Egorov, R.Z. Akhmetshin, G.B. Efimov, G.G. Fedotov, V.G. Petukhov. Some Methodical Problems of Low-Thrust Trajectory Optimization. Preprint of KIAM No. 110, Moscow, Medvedev A., V. Khatulev, V. Yuriev, V. Petukhov, M. Konstantinov. Combined flight profile to insert telecommunication satellite into geostationary orbit using Rockot light-weight class launch vehicle. 51 st International Astronautical Congress. IAF-00-V.2.09, Rio de Janeiro, Brasilia, October 2-6, Konstantinov M.S., G.G. Fedotov, V.G. Petukhov, et al. Electric Propulsion Mission to GEO Using Soyuz/Fregat Launch Vehicle. 52 nd International Astronautical Congress. IAF 01 V.3.02, Toulouse, France, October 1-5, Petukhov V.G. Optimal Multirevolutional Transfers between Non-Coplanar Elliptical Orbits. Proceedings of 2 nd International Symposium Low Thrust Trajectories (LOTUS-2). IAS, Toulouse, France, June 18-20, Konstantinov M.S. Optimization of Low Thrust Transfer from Elliptical Orbit into Non-Coplanar Circular Orbit. Proceedings of 2 nd International Symposium Low Thrust Trajectories (LOTUS-2). IAS, Toulouse, France, June 18-20, Fedotov, G.G., M.S. Konstantinov, V.G. Petukhov, G.A. Popov. Low-cost electric propulsion mission to the near-Earth asteroids. IAF-95-A.6.07, Oslo, Norway, Fedotov G.G., M.S. Konstantinov, V.G. Petukhov. Electric propulsion mission to Jupiter. IAF-96-A.1.10, China, October Fedotov G.G., M.S. Konstantinov, V.G. Petukhov, G.A. Popov. Electric propulsion mission to the main belt asteroid: optimal parameters of the ion engine. IAF-97-S.3.03, Turin, Italy, October Konstantinov M.S. Optimization of low thrust transfer between noncoplanar elliptic orbits. IAF-97-A.6.06, Turin, Italy, October Konstantinov M., G. Fedotov. Electric propulsion mission to Mercury. Second European Spacecraft Propulsion Conference, May, 1997 (ESA SP-398, Aug.1997). 24. Konstantinov M.S., G.G. Fedotov. Asteroid mission at use of two-mode electric thrusters. IAF-99-Q.5.04, Amsterdam, The Netherlands, October 4-8, Krotov, V.F., V.I. Gurman. Optimal Control Methods and Problems. Moscow, Nauka, 1973, 446 pp. (in Russian). 26. Betts, J.T., S.O. Erb. Optimal Low Thrust Trajectories to Moon ( ) 31 ESA/ESTEC Workshop on Trajectory Design and Optimization. Noordwijk, October 24-25, 2002 V. Petukhov. Low-Thrust Trajectory Optimization in Russia