Trajectory Planning for Spacecraft Rendezvous with On / Off Thrusters

Trajectory Planning for Spacecraft Rendezvous with On/Off Thrusters ? Rafael Vazquez ∗Francisco Gavilan ∗Eduardo F. Camacho ∗∗ ∗Departamento de Ingenieria Aerospacial ∗∗ Departamento de Ingenieria de Sistemas y Automatica Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Sevilla, Spain (e-mails: [email protected], [email protected], eduar[email protected]) Abstract: The objective of this work is to present a trajectory planning algorithm for spacecraft rendezvous that is able to incorporate Pulse-Width Modulated (PWM) control signals. The algorithm is based on linearization around a previously computed solution. To initialize the algorithm, a first solution needs to be obtained. To do so, the trajectory planning problem is solved using Pulse-Amplitude Modulated (PAM) control signals; these are then converted to PWM signals, which are used as an initial guess. Iterating, the solution is refined until an optimal value is reached. Simulations show that this method converges after a few iterations. The algorithm is simple and fast, hence it could be implemented online or used together with a Model Predictive Controller. Keywords: Spacecraft autonomy, Space robotics, Pulse-width modulation, Trajectory planning, Optimal trajectory, Linearization. 1. INTRODUCTION Technology enabling simple autonomous spacecraft rendezvous and docking is becoming a growing necessity as access to space continues increasing. After decades of development, many approaches have been proposed; see Woffinden and Geller (2008) for an historical account or Fehse (2003) for the basics. Classicaly (Wie (1998)), the problem is modeled using impulsive maneuvers, computing an initial and final impulse (∆V) to achieve rendezvous. Impulsive changes of speed are also used in recent works (Tong et al. (2007); Geller (2006)). Lately, approaches based on trajectory planning and optimization (Breger and How (2008)) and predictive control (Richards and How (2003); Rossi and Lovera (2002); Asawa et al. (2006)) are emerging. However, typically these methods allow the control signal (thrust) to take any value in an allowable range. This type of control signal is usually referred to as Pulse-Amplitude Modulated (PAM). A more realistic modeling of spacecraft thrusters would take into account that typically, thrusters are ON-OFF actuators, i.e., the thrusters are not able to produce any value of force, but can be only switched on (producing the maximum amount of force) or off (producing no force). Thus only the times where the thrusters are turned ON or OFF (the switching times) can be controlled. This type of control signal is usually referred to as Pulse-Width Modulated (PWM). ?The authors would like to acknowledge finantial support of the Spanish Ministry of Science and Innovation and of the European Commission for funding part of this work under grants DPI200805818 and project FPT ICT HD-MPC. The objective of this work is to present a trajectory planning algorithm for spacecraft rendezvous that is able to incorporate PWM control signals. The algorithm is based on the fact that, even though in the rendezvous model the PWM signals appear nonlinearly, this dependence on the switching times is only weakly nonlinear. Thus, linearization around reference values can faithfully capture the dependence of the system on the switching times, at least close to the reference values. To obtain a first solution around which linearize the systems, the trajectory planning problem is solved using PAM control signals; these are then converted to PWM signals, which are used as initialization for our algorithm. Iterating, the solution is refined until an optimal value is reached. The resulting algorithm is simple and fast and could be implemented online, or with a Model Predictive Controller (such as Gavilan et al. (2009)), allowing the MPC to include PWM control signals. While the idea of linearization to compute optimal PWM control signals is, as far as we know, original, we note that local linearization techniques have been used for optimal trajectory problems in other contexts (see e.g. Kim et al. (2002)). The structure of this paper is as follows. In Section 2 we introduce the mathematical model for rendezvous spacecraft used for trajectory planning, both in the PAM and PWM case. We follow in Section 3 where we formulate the planning problem, detailing the equality and inequality constraints, and the cost function. Section 4 describes the method we propose to solve the planning problem. In Section 5 we show simulations of the proposed method. We close the paper with some remarks in Section 6. 2. MODEL OF SPACECRAFT RENDEZVOUS In this work we use the the linear Hill-Clohessy-Wiltshire (HCW) equations, as introduced in Hill (1878) and CloProceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011 978-3-902661-93-7/11/$20.00 © 2011 IFAC 8473 10.3182/20110828-6-IT-1002.02445 hessy and Wiltshire (1960). Those describe the relative position of the spacecraft if the target is orbiting in a circular keplerian orbit and approaching vehicle is close enough to the target. The HCW model also assumes that the target vehicle is passive and moving along a circular orbit of radius R. Thus the angular speed of the target through its orbit is n=pµ R3, where µis the gravitation parameter of the Earth, µ= 398600.4 km3/s2. Denote by kthe time instant t=t0+kT , where Tis an adequately chosen sampling time (the sampling time chosen in this work is 60 s). We formulate two versions of the HCW equations in discrete time. In the first version, the control inputs are constant during a whole sampling time. This is referred to as the PAM discrete model. In the second version, it is assumed that thrusters are not able to produce any force value, but can be only switched on (producing the maximum force) or off (producing no force), and only once during each sampling time. This is referred to as the PWM discrete model. 2.1 PAM discrete model Assuming the control signal constant through the sampling time, it can be shown that x(k+ 1) = ATx(k) + BTu(k).(1) In (1), x(k), u(k) denote respectively the state at time k, and the input amplitude (from time kto k+ 1). Also, x= [x y z ˙x˙y˙z]T,u= [uxuyuz]T.(2) In these definitions, x,y, and zdenote the position of the chaser in a local–vertical/local–horizontal (LVLH) frame of reference fixed on the center of gravity of the target vehicle. In the LVLH frame, xrefers to the radial position, yto the in-track position, and zto the cross-track position. The velocity of the chaser in the LVLH frame is given by ˙x, ˙y, and ˙z. The variables ux,uy, and uzare the inputs (thrust actuation) acting on the chaser vehicle, referred to the LVLH axes as indicated by their respective subscripts. The matrices ATand BTappearing in (1) are given by AT=             4−3C0 0 S n 2(1 −C) n0 6(S−nT) 1 0 −2(1 −C) n 4S−3nT n0 0 0 C0 0 S n 3nS 0 0 C2S0 −6n(1 −C) 0 0 −2S4C−3 0 0 0 −nS 0 0 C             ,(3) BT=                 1−C n2 2nT −2S n20 2(S−nT) n2−3T2 2+ 41−C n20 0 0 1−C n2 S n21−C n0 2(C−1) n−3T+ 4S n0 0 0 S n                 ,(4) where S= sin nT and C= cos nT. Compact formulation Next we develop a compact formulation that simplifies the notation of the problem. The state at time k+ 1, given the initial state at time 0 (which is denoted as x0) and the input signals from time 0 to time k, is computed by applying recursively Equation (1): x(k) = Akx0+ k−1 X j=0 Ak−1−jBu(j).(5) Define now xSand uSas a stack of Npstates and input signals, respectively, spanning from time 1 to time Np, where Npis the planning horizon: xS=    x(1) x(2) . . . x(Np)     ,uS=    u(0) u(1) . . . u(Np−1)     . Then, xS=             ATx0+BTu(0) A2 Tx0+ 1 X j=0 A1−j TBTu(j) . . . ANp Tx0+ Np−1 X j=0 ANp−1−j TBTu(j)             ,(6) which can be written as xS=Fx0+GuS,(7) where Gis a block lower triangular matrix with its nonnull elements defined by (G)ij =Ai−j TBTand the matrix Fis defined as: F=     AT A2 T . . . ANp T      .(8) 2.2 PWM discrete formulation Considering now the case of ON-OFF thrusters, it is assumed that the thrusters are not able to produce any value of force, but can be only switched on or off. Also, it is assumed that there is an aligned pair of thrusters for each direction i= 1,2,3 with opposing orientation. To distinguish between the positive and negative they are denoted as u+ iand u− i, whereas the respective maximum thrust is referred to as u+ imax and u− imax , respectively. Finally, during each sample time each thruster is allowed to fire only once. Thus, the PWM output for each time interval kis completely described by two new control variables for each pair of thrusters: the pulse width κ+ i(k) and the pulse start time τ+ i(k) (for the positively oriented thruster in the direction i) and similarly κ− i(k) and τ− i(k) for the negatively oriented thruster in the direction i. Then, for t∈[kT, (k+ 1)T], we have: u+ i(t) =    0, t ∈kT, kT +τ+ i(k), u+ imax , t ∈kT +τ+ i(k), kT +τ+ i(k) + κ+ i(k), 0, t ∈kT +τ+ i(k) + κ+ i(k),(k+ 1)T, (9) 18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011 8474 and similarly for the negatively oriented thrusters. The new variables control variables verify κ+ i(k)>0, τ+ i(k)> 0 and τ+ i(k) + κ+ i(k)< T, and similarly for the negatively oriented thrusters. The last constraint prevent the PWM signal to spill over to the next time interval. Call the PWM control variables as uP(k): uP(k) =                    τ+ 1(k) κ+ 1(k) τ− 1(k) κ− 1(k) τ+ 2(k) κ+ 2(k) τ− 2(k) κ− 2(k) τ+ 3(k) κ+ 3(k) τ− 3(k) κ− 3(k)                    .(10) To find the HCW equations for PWM inputs, define b1 t=              1−C n2 2(S−nt) n2 0 S n 2(C−1) n 0              , b2 t=              2nt −2S n2 −3t2 2+ 41−C n2 0 21−C n −3t+ 4S n 0              , b3 t=           0 0 1−C n2 0 0 S n           ,(11) which are the three columns of the matrix Btin (4), but defined for different t. Then the system evolution equation (1) is replaced in the PWM case by x(k+ 1) = Ax(k) + BPWM(uP(k))umax,(12) where BPWM(uP(k)) =            AT−τ+ 1(k)−κ+ 1(k)b1 κ+ 1(k) AT−τ− 1(k)−κ− 1(k)b1 κ− 1(k) AT−τ+ 2(k)−κ+ 2(k)b2 κ+ 2(k) AT−τ− 2(k)−κ− 2(k)b2 κ− 2(k) AT−τ+ 3(k)−κ+ 3(k)b3 κ+ 3(k) AT−τ− 3(k)−κ− 3(k)b3 κ− 3(k)            T ,(13) umax =         u+ 1max −u− 1max u+ 2max −u− 2max u+ 3max −u− 3max         .(14) Notice that system (12) is nonlinear in the PWM variables uP(k). Compact formulation The compact formulation developed before can be readily adapted to PWM inputs. Equation (7) is now written as xS=Fx0+GPWM(uPS)umaxS,(15) where uPSis a stack vector with all the PWM signals, GPWM is a block lower triangular matrix with its non-null elements defined by (GPWM(uPS))ij = Ai−jBPWM(uP(j−1)), and the constant vector umaxS is defined as: umaxS=    umax umax . . . umax            N times.(16) 3. FORMULATION OF THE PLANNING PROBLEM Next we formulate our planning problem, introducing the constraints (both equalityand inequality-type constraints) and the objective function that has to be minimized. The formulation is done for both PAM and PWM control signals. 3.1 Constraints on the problem Inequality constraints on the state For sensing purposes (see Breger and How (2008)), during rendezvous it is required that the chaser vehicle remains inside a line of sight (LOS) area. To simplify the constraint, we consider a 2-D LOS area as shown in Figure 1. This LOS region is the intersection of a cone, given by the equations y≥cLOS(x−x0) and y≥ −cLOS(x+x0), and the region y≥0. Fig. 1. Line of Sight region. The LOS constraint is ALOS x(k)≤bLOS, where ALOS ="0−10000 cLOS −10000 −cLOS −10000#, bLOS ="0 cLOSx0 cLOSx0#. (17) Using the compact formulation that was developed in Section 2.1, the constraints equations for the state can be rewritten as: AcxS≤bc,(18) 18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011 8475 where Acand bcare given by: Ac=    ALOS ALOS ... ALOS     ,bc=    bLOS bLOS . . . bLOS     .(19) Then, for the case of PAM control and using equation (7), one can reformulate the LOS constraints as constraints for the control signals in the following way: AcGuS≤bc−AcFx0,(20) and similarly for the case of PWM control. Equality constraints on the state Equality constraints are formulated to ensure that the chaser spacecraft arrives at the origin with zero velocity at the end of the planning horizon. Thus, these constraints can be written as x(Np) = 0. Defining Aeq as Aeq =    0 ... 0 Id6×6     ,(21) the equality constraint can be written as: AeqxS=0.(22) Then, for the case of PAM control and using equation (7), one can reformulate the arrival conditions as equality constraints for the control signals in the following way: AeqGuS=−AeqFx0,(23) and similarly for the case of PWM control. Input constraints For the case of PAM control, one would have limitations on the magnitude of the control, given by umin ≤uS≤umax.(24) For the case of PWM control, the constraints are given by κ± i(k)>0, τ± i(k)>0 and τ± i(k) + κ± i(k)< T. 3.2 Objective function The objective function to be minimized in the planning problem is the 1-norm of the control signal, which is proportional to fuel consumption. PAM control inputs For the case of PAM control inputs, the selected control function is given by: JP AM = Np−1 X k=0 kuT(k)k1=kuSk1.(25) PWM control inputs For the case of PAM control inputs, using (9) it can be seen that the selected objective function is given by: JPWM = Np−1 X k=0 3 X i=1 u+ imax κ+ i(k) + u− imax κ− i(k).(26) 4. COMPUTATION OF THE OPTIMAL CONTROL INPUT As shown in Section 2.2, the discrete HWC equations in the PWM case are nonlinear in the switching times. Thus, it is neither easy nor fast to solve an optimization problem using the switching times as control variables. To find the optimal control input without needing to solve a nonlinear planning problem, the following scheme is proposed: Step 1. A PAM linear optimization problem is solved. Step 2. A PAM/PWM filter is used to convert the PAM input signals to PWM signals (i.e., switching times) that produce a very similar system output. Step 3. The plant equations, in the PWM formulation, are linearized around the previous step solution, thus obtaining a linear plant with respect to the switching time. A linear optimization problem is then posed and solved. The resulting solution is taken as a better approximation towards the real solution. Step 4. Repeat the linearization process of Step 3 around the new solution. Optimize again to find a better refinement. The process is iterated until the solution converges. Next, we describe all the step in our scheme. 4.1 Computation of PAM control input To compute the optimal control plan (with PAM control signals), one seeks the control signal that minimizes the cost function over the planning horizon, satisfying at the same time the PAM constraints: min uS JP AM (uS) (27) subject to: AcGuS≤bc−AcFx0 AeqGuS=−AeqFx0 umin ≤uS≤umax. Since the cost function and the constraints are linear, then (27) can be readily solved. 4.2 Initial PWM solution: A PAM/PWM filter The PAM solution found when solving (27) is transformed to an equivalent PWM solution using a PAM/PWM filter. This filter is formulated in the literature (for instance in Shieh et al. (1996); Ieko et al. (1999)) where several methods are proposed. These methods allow to, given the PAM inputs of a system, compute equivalent PWM inputs that produce a system output optimally approximating the output of the system when driven by the PAM signals. Following these references, a PAM control signal can be optimally approximated by a PWM control signal by using the following rules for each time instant kand direction i: (1) Use the positive or negative thruster according to the sign of the PAM signal ui(k). (2) The pulse width must be computed using the Principle of Equivalent Areas: κ± i(k) = T|ui(k)| u± imax , where u± imax is the maximum level of the (positive or negative) thruster i. 18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011 8476 (3) Since only one impulse per sample time is considered, it must be allocated at the center of the sampling interval, that is τ± i(k) = 1 2(T−κ± i(k)). The PWM signals uP(k) constructed by this method produce an almost identical output to the system driven by PAM signals. However, the PWM are not necessarily optimal since their constraints are quite different; in fact they might even not verify the constraints. Thus, this solution is only used as an initialization for the optimization algorithm proposed next. 4.3 Refined PWM solution: An optimization algorithm Noting that most nonlinearities appearing in the system equations (12) are of the form cos nt or sin nt, where n is the orbital angular velocity, and noting nt 1, we approximate (12) by linearizing BPWM(uP(k)) around uP(k). Then, the system equations are: x(k+1) = Ax(k)+BPWM(uP(k))umax+B∆(uP(k))∆(k), (28) where B∆=                              −A0 T−τ+ 1−κ+ 1 b1 κ+ 1 u+ 1max −A0 T−τ+ 1−κ+ 1 b1 κ+ 1 +AT−τ+ 1−κ+ 1b10 κ+ 1u+ 1max A0 T−τ− 1−κ− 1 b1 κ− 1 u− 1max A0 T−τ− 1−κ− 1 b1 κ− 1 −AT−τ− 1−κ− 1b10 κ− 1u− 1max −A0 T−τ+ 2−κ+ 2 b2 κ+ 2 u+ 2max −A0 T−τ+ 2−κ+ 2 b2 κ+ 2 +AT−τ+ 2−κ+ 2b20 κ+ 2u+ 2max A0 T−τ− 2−κ− 2 b2 κ− 2 u− 2max A0 T−τ− 2−κ− 2 b2 κ− 2 −AT−τ− 2−κ− 2b20 κ− 2u− 2max −A0 T−τ+ 3−κ+ 3 b3 κ+ 3 u+ 3max −A0 T−τ+ 3−κ+ 3 b3 κ+ 3 +AT−τ+ 3−κ+ 3b30 κ+ 3u+ 3max A0 T−τ− 3−κ− 3 b3 κ− 3 u− 3max A0 T−τ− 3−κ− 3 b3 κ− 3 −AT−τ− 3−κ− 3b30 κ− 3u− 3max                              T ,(29) and where the value of κand τvariables at time kis used. The variable ∆(k) represents the increments or decrements with respect to uP(k): ∆(k)=                    ∆τ+ 1(k) ∆κ+ 1(k) ∆τ− 1(k) ∆κ− 1(k) ∆τ+ 2(k) ∆κ+ 2(k) ∆τ− 2(k) ∆κ− 2(k) ∆τ+ 3(k) ∆κ+ 3(k) ∆τ− 3(k) ∆κ− 3(k)                    .(30) The matrices used in (29) are defined as: A0 t=       3nS 0 0 C2S0 6n(C−1) 0 0 −2S4C−3 0 0 0 −nS 0 0 C 3n2C0 0 −nS 2nC 0 −6n2S0 0 −2nC −4nS 0 0 0 −n2C0 0 −nS        ,(31) and b10 t=           S n 2(C−1) n 0 C −2S 0           , b20 t=           2−2C n −3t+ 4S n 0 2S −3+4C 0           , b30 t=         0 0 S n 0 0 C         .(32) Equation (33) is now compactly written as xS=Fx0+G∆(uPS)∆S+GPWM(uPS)umaxS,(33) where G∆(uPS) is a block lower triangular matrix with its non-null elements defined by (G∆(uPS))ij = Ai−jB∆(uP(Np−j−1)) and ∆Sis a stack vector of the increment in the PWM variables ∆(k). The LOS inequality constraints (20) can be reformulated as constraints for the PWM variables in the following way: AcG∆∆S≤bc−AcFx0−AcGPWMumaxS,(34) where the dependence of GPWM(uPS) and G∆(uPS) on uPShas been omited for simplicity. Similarly, the equality constraints become: AeqG∆∆S=−AeqFx0−AeqGPWMumaxS.(35) The constraints on the ∆(k) are as follows: −∆κ± i(k)≤κ± i(k),−∆τ± i(k)≤τ± i(k) (36) ∆τ± i(k)+∆κ± i(k)≤T−τ± i(k)−κ± i(k),(37) |∆(k)| ≤ ∆MAX ,(38) where (38) is used to avoid too large variations in each iteration step. These constraints can be summarized as A∆(k)∆S(k)≤b∆(k).(39) Finally, the objective function can be rewritten as a function of the PWM variables and their increments as J(uPS,∆S) = JPWM(uPS) + J∆(∆S), where J∆(∆S) = Np−1 X k=0 3 X i=1 u+ imax ∆κ+ i(k) + u− imax ∆κ− i(k).(40) Thus, a linear programming problem with PWM outputs can be posed as follows: min ∆S J∆(∆S) (41) s. t.: AcG∆∆S≤bc−AcFx0−AcGPWMumaxS, AeqG∆∆S=−AeqFx0−AeqGPWMumaxS, A∆∆S≤b∆. The solution of (41), ∆S, can be used to recompute new PWM variables: uNEW PS=uPS+∆S.(42) Then, the values of uNEW PSis used to recompute the various matrices appearing in (41), and the optimization problem is solved again. Iterating, it is expected that the solution keeps improving until a solution close to the global optimum is reached. Next, this is shown in simulations. 5. SIMULATION RESULTS We next show simulations of the trajectory planning problem for spacecraft rendezvous using our algorithm. 18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011 8477 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x (km) y (km) PWM PAM Fig. 2. Solution trajectory, showing the initial trajectory with PAM inputs (dashed) and the final trajectory with PWM inputs (solid). ! " # $ % & %'!( %'" %'"" %'"$ %'"& %'"( %'# %'#" %'#$)*!+!# ,-./01*23*45016542,7 !"#$%&'( Fig. 3. Value of objective function depending on iteration number. ! "!! #!!! #"!! $!!! $"!! %!!! !& !$ ! $'(#!!" )*+, -'*./+$, ! "!! #!!! #"!! $!!! $"!! %!!! !" ! "'(#!!" )*+, -0*./+$, Fig. 4. Initial PAM control signal. The planning horizon Npwas chosen as 50. The value of u± imax was chosen as 10−4km/s2. Initials conditions were r0= [0.2 0.6 0.2]Tkm,v0= [0.0015 0.002 −0.001]Tkm/s. The solution trajectory is shown in Fig. 2 (solid). The PAM solution that was used to initialize the algorithm is shown in the same figure (dashed). In this and other figures only the xand ycoordinates are shown for the sake of brevity. As it can be seen in Fig. 3, the algorithm converged after 6 iterations to a solution with a cost function approximately 2.5% smaller. Each iteration took less than a second on a conventional computer, using MATLAB’s linprog to calculate the solution. Comparing the initial PAM (Fig. 4) and PWM (Fig. 5) control signals, it can be seen that for both the control effort is concentrated at the beginning, with some mid-course and final maneuvers that are different for both cases. ! "!! #!!! #"!! $!!! $"!! %!!! !$ !# ! # $&'#!!( )*+, -&*./+$, ! "!! #!!! #"!! $!!! $"!! %!!! !$ ! $&'#!!( )*+, -0*./+$, Fig. 5. Final PWM control signals. 6. CONCLUDING REMARKS We have presented an algorithm to compute the optimal PWM control signals applied to rendezvous of spacecraft. The algorithm uses the HCW model with LOS and equality constraints, but can be generalized to more complicated models and constraints. Since the algorithm is very fast, it could be used online to implement an MPC scheme (e.g. Gavilan et al. (2009)) including PWM control signals. REFERENCES Asawa, S., Nagashio, T., and Kida, T. (2006). Formation flight of spacecraft in earth orbit via MPC. In SICEICASE International Join Conference. Breger, L. and How, J.P. (2008). Safe trajectories for autonomous rendezvous of spacecraft. J Guid Contr Dynam., 31(5), 1478–1489. Clohessy, W.H. and Wiltshire, R.S. (1960). Terminal guidance systems for satellite rendezvous. Journal of the Aerospace Sciences, 27(9), 653–658. Fehse, W. (2003). Automated Rendezvous and Docking of Spacecraft. Cambridge University Press. Gavilan, F., Vazquez, R., and Camacho, E.R. (2009). Robust model predictive control for spacecraft rendezvous with online prediction of disturbance bound. In Proceedings of AGNFCS’09, Samara, Russia,. Geller, D. (2006). Linear covariance techniques for orbital rendezvous analysis and autonomous onboard mission planning. J Guid Contr Dynam., 29(6), 1404–1414. Hill, G. (1878). Researches in lunar theory. American Journal of Mathematics, 1(3), 5–26, 129–147, 245–260. Ieko, T., Ochi, Y., and Kanai, K. (1999). New design method for pulse-width modulation control systems via digital redesign. J Guid Contr Dynam., 22(1), 123–128. Kim, H.J., Shim, D.H., and Sastry, S. (2002). Nonlinear model predictive tracking control for rotorcraft-based unmaned aerial vehicles. In Proceedings of ACC 2002. Richards, A.G. and How, J. (2003). Performance evaluation of rendezvous using model predictive control. AIAA Paper 2003-5507. Rossi, M. and Lovera, M. (2002). A multirate predictive approach to orbit control of small spacecraft. In Proceedings of ACC 2002. Shieh, L.S., Wang, W.M., and Sunkel, J. (1996). Design of PAM and PWM controllers for sampled-data interval systems. J Dyn Syst Meas Contr., 118(4), 673–681. Tong, C., Shijie, X., and Songxia, W. (2007). Relative motion control for autonomous rendezvous based on classical orbit element differences. J Guid Contr Dynam., 30(4), 1003–1014. Wie, B. (1998). Space vehicle dynamics and control. AIAA. Woffinden, D.C. and Geller, D.K. (2008). Navigating the road to autonomous orbital rendezvous. Journal of Spacecraft and Rockets, 44(4), 898–909. 18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011 8478