Pulse-width predictive control for LTV systems with application to spacecraft rendezvous

Pulse-Width Predictive Control for LTV Systems with Application to Spacecraft Rendezvous R. Vazqueza,∗, F. Gavilana, E. F. Camachob aDepartamento de Ingenier´ıa Aeroespacial. Escuela Superior de Ingenieros. Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092,Sevilla, Spain bDepartamento de Ingenier´ıa de Sistemas y Autom´atica. Escuela Superior de Ingenieros. Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092,Sevilla, Spain Abstract This work presents a model predictive controller (MPC) that is able to handle linear time-varying (LTV) plants with PWM control. The MPC is based on a planner that employs a PAM or impulsive approximation as a hot-start and then uses explicit linearization around successive PWM solutions for rapidly improving the solution by means of linear programming. As an example, the problem of rendezvous of spacecraft for eccentric target orbits is considered. The problem is modeled by the LTV Tschauner-Hempel equations, whose transition matrix is explicit; this is exploited by the algorithm for rapid convergence. The efficacy of the method is shown in a simulation study. Keywords: Spacecraft autonomy, Space robotics, Pulse-width modulation, Trajectory planning, Optimal trajectory, Linear Time-Varying Systems. 1. Introduction Aerospace systems often need to be controlled by using pulse-width modulated (PWM) actuators, i.e., actuators whose output level is fixed and can only be turned on and off, such as spacecraft thrusters. It would be therefore desirable to use control design methods that directly take into account pulsed actuators. However, most feedback design and motion planning methods ignore ∗Corresponding author, Fax number +34954486041, Email address: rv[email protected] Email addresses: [email protected] (F. Gavilan), [email protected] (E. F. Camacho) Preprint submitted to Control Engineering Practice October 2, 2018 arXiv:1511.00869v1 [math.OC] 3 Nov 2015 variable width pulses and approximate the control variables either by impulses (which produce instantaneous changes in some combination of the states) or pulse-amplitude modulated (PAM) control. However, neither impulsive actuation nor PAM actuation capture with precision the behavior of pulsed actuators such as spacecraft thrusters. A more realistic model has to take into account that, typically, thrusters are ON-OFF actuators, i.e., the thrusters are not able to produce arbitrary forces, but instead can only be switched on (producing the maximum amount of force) or off (producing no force). These switching times are the only signals that can be controlled. This type of control signal is usually referred to as Pulse-Width Modulated (PWM). Control design with PWM actuation poses a challenge because the system becomes nonlinear in the switching times, even if the system is linear. One can find in the literature several procedures to find an equivalent PWM solution starting from a PAM solution (for instance in Shieh et al. (1996); Ieko et al. (1999); Bernelli-Zazzera et al. (1998)). These methods allow to, given the PAM inputs of a system, compute PWM inputs that produce a system output optimally approximating the output of the system when driven by the PAM signals. The results are based on the so-called Principle of Equivalent Areas, which computes the PWM signal so that it covers the same area as the PAM signal. However, while these procedures are quite effective in the sense that the output produced by the approximate PWM signals is very similar to the one produced by PAM signals, they assume that the plant is linear time-invariant. In this paper, Model Predictive Control (MPC) is used to directly find PWM signals to control the system. MPC (see, e.g., Camacho and Bordons (2004)) is a family of methods that originated in the late seventies and has developed considerably since then. In MPC, the process model is used to predict the future plant outputs, based on past and current values and on the proposed optimal future control actions. These actions are calculated by the optimizer taking into account the cost function as well as the constraints. Since the plant is nonlinear in the control signals (ON-OFF times), the underlying optimization problem is nonlinear and possibly non-convex. To solve the problem, the algorithm starts 2 from an initial guess computed by solving an optimal linear program with PAM or impulsive actuation, approximate the solution with ON-OFF thrusters, and then iteratively linearize around the obtained solutions to improve the PWM solution. While the idea of linearization to specifically compute optimal PWM control signals in the context of MPC is, to the best knowledge of the authors, original, it must be noted that local linearization techniques have been used for optimal trajectory problems in other contexts (see e.g. Kim et al. (2002)). As an application the problem of rendezvous of spacecraft is considered, i.e., the controlled close encounter of two space vehicles. Autonomous spacecraft rendezvous capabilities are becoming a necessity as access to space continues increasing. The field has become very active in recent years, with a rapidly growing literature. Among others, approaches based on trajectory planning and optimization (Breger and How (2008); Arzelier et al. (2013, 2011); Louembet et al. (2015); Deaconu et al. (2015, 2014); D’Amico et al. (2013); Gaias et al. (2014)) and predictive control (Richards and How (2003); Rossi and Lovera (2002); Asawa et al. (2006); Gavilan et al. (2009, 2012); Larsson et al. (2006); Hartley et al. (2012); Leomanni et al. (2014); Jewison et al. (2015); Weiss et al. (2012)) are emerging. Classically, in these approaches the problem of rendezvous is modeled by using impulsive maneuvers; one computes a sequence of (possibly optimal) impulses (usually referred to as ∆V’s) to achieve rendezvous. Recently, Vazquez et al. (2011, 2014) introduced a trajectory planning algorithm algorithm for spacecraft rendezvous that was able to incorporate PWM control signals. The former considered the linear time-invariant Clohessy-Wiltshire model (target orbiting in a circular Keplerian orbit, see Clohessy and Wiltshire (1960)). The latter extended the approach to elliptical target orbits by using the linear time-varying Tschauner-Hempel model (see Tschauner and Hempel (1965)). Both methods start from an initial guess computed by solving an optimal linear program with PAM or impulsive actuation, approximate the solution with ON-OFF thrusters, and then iteratively linearize around the obtained solutions to improve the PWM solution. For both circular and elliptical target 3 orbits the algorithms are simple and reasonably fast, and simulations favorably compare with an impulsive-only approach. These results were extended in Vazquez et al. (2015) to a decreasing-horizon model predictive controller able to take into account orbital perturbations, disturbances or model errors. In this paper, a receding horizon model predictive controller with PWM inputs is formulated for general LTV plants and both alternatives (PAM or impulsive starting guess) are discussed in detail, with an application to rendezvous given at the end of the paper. The structure of the paper is as follows. In Section 2 the plant model is introduced. Three types of inputs are considered: PWM, PAM and impulsive. Section 3 follows with a formulation of the underlying optimization problem. Section 4 describes a method that solves the planning problem using PWM signals. Section 5 develops the model predictive controller. Next, Section 6 describes the application to spacecraft rendezvous. Section 7 presents a simulation study of the method applied to spacecraft rendezvous. The paper finishes with some remarks in Section 8. 2. System Model Consider a linear time-varying system given as ˙x=A(t)x+B(t)u, (1) where x∈Rnis the state, u∈Rmis the input (control) vector, and A(t) and B(t) are, respectively, n×nand n×mmatrices depending on time t≥0. Considering that, for some time tk≥0, initial conditions x(tk)∈Rnare given and the input is known, the solution to (1) for t>tkis given by x(t) = Φ(t, tk)x(tk) + Zt tk Φ(t, s)B(s)u(s)ds, (2) where Φ(t, tk) is the system transition matrix, see for instance Rugh (1996). This matrix can be computed numerically (or analytically if possible) as the 4 unique solution to the linear matrix differential equation ˙ Φ(t, tk) = A(t)Φ(t, tk), t > tk(3) Φ(tk, tk)=I.(4) To obtain an unified notation in term of the inputs, denote by Bi(t) the i-th column of B(t), corresponding to the i-th input ui(t), for i= 1, . . . , m. In the paper, time intervals starting at some initial time tkand ending at tk+1 =tk+T are considered, where Twill be an adequate sample time. Then equation (2) can be written as x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 Ztk+1 tk Φ(tk+1, s)Bi(s)ui(s)ds, (5) The objective is solving the problem with PWM inputs. In addition, two other types of inputs are considered; they will be used as an intermediate step towards computing PWM inputs by the algorithm. All types of input are analyzed in the following sections. 2.1. Pulse width-modulated (PWM) control In the PWM case, each input uiis a pulse starting at time τk,i (relative to tk) with pulse width κk,i, with constant magnitude umax =uW k,i, as shown in Fig. 1, i.e., ui(t) =          0, t ∈[tk, tk+τk,i], uW k,i, t ∈[tk+τk,i, tk+τk,i +κk,i], 0, t ∈[tk+τk,i +κk,i, tk+1], (6) with κk,i >0, τk,i >0 and τk,i +κk,i < T , where the last constraint prevent the PWM signal to spill over to the next time interval. Then, substituting ui(t) in (5) one obtains x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 Ztk+τk,i+κk,i tk+τk,i Φ(tk+1, s)Bi(s)ds!uW k,i,(7) and denoting BW k,i(τk,i, κk,i) = Ztk+τk,i+κk,i tk+τk,i Φ(tk+1, s)Bi(s)ds, (8) 5 max u t T · ¿ Figure 1: PWM Variables. one can write the solution as x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 BW k,i(τk,i, κk,i)uW k,i.(9) There is an important difference between a PWM input and the PAM or impulsive inputs that will be subsequently introduced. While the latter can in principle take positive o negative values at different times, the former is fixed either as positive or negative for all time. Thus, typically a PWM model has twice number of inputs than a PAM/impulsive model. To make this explicit in the model (10), denote with a plus or minus super-index the positive or negative inputs, as follows x(tk+1) = Φ(tk+1, tk)x(tk)+ m X i=1 hBW k,i(τ+ k,i, κ+ k,i)uW+ k,i −BW k,i(τ− k,i, κ− k,i)uW− k,i i, (10) with uW+ k,i , τ+ k,i, κ+ k,i and uW− k,i , τ− k,i, κ− k,i denoting, respectively, the magnitude, start, and width of the positive and negative i-th input pulses. 2.2. Pulse amplitud-modulated (PAM) control In this case, each control ui(t) in (5) is constant inside the interval [tk, tk+1], and equal to uA k,i. Then, substituting ui(t) in (5) one obtains x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 Ztk+1 tk Φ(tk+1, s)Bi(s)dsuA k,i,(11) and denoting by BA k,i =Ztk+1 tk Φ(tk+1, s)Bi(s)ds, (12) 6 the solution can be written as x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 BA k,iuA k,i.(13) 2.3. Impulsive control In this case, ui(t) = uI k,iδ(t−(tk+τk,i)), where δ(t) is Dirac’s delta function, tk+τk,i is the instant at which the impulse is given, and uk,i is the magnitude of the impulse. Then, assuming 0 < τk,i < T for all i(all the impulses are given inside the considered time interval) and substituting ui(t) in (5) one obtains x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 Φ(tk+1, tk+τi)Bi(tk+τi)uI k,i,(14) and denoting by BI k,i(τk,i) = Φ(tk+1, τi)Bi(tk+τi), x(tk+1) = Φ(tk+1, tk)x(tk) + m X i=1 BI k,i(τk,i)uI k,i.(15) 2.4. Discretization and compact notation Consider now a sequence of time instants tk=t0+kT,k= 0, . . ., and denote xk=x(tk). Then, it is possible to write, for both PAM and impulsive control, xk+1 =Akxk+BkUk,(16) where Ak= Φ(tk+1, tk), and Bkand Ukdepend on the input type. In the PWM case, write xk+1 =Akxk+B+ kU+ k−B− kU− k,(17) where B+ kis a matrix whose i-th column is BW k,i(τ+ k,i, κ+ k,i) and U+ k,τ+ k, and κ+ kare column vectors whose i-th entries are, respectively, uW k,i,τ+ k,i and κ+ k,i. The same definitions (with minus super-index) are used for the negative pulses. Then, to reach model (16), define Bk=hB+ k−B− ki, Uk=  U+ k U− k , τk=  τ+ k τ− k , κk=  κ+ k κ− k .(18) The definitions of Bkare simpler for the other types of actuation. In the PAM case, Bkis a matrix whose i-th column is BA k,i and Uka column vector 7 whose i-th entry is uA k,i. In the impulsive case, Bkis a matrix whose i-th column is BI k,i(τk,i), and Uk,τkare column vector whose i-th entries are, respectively, uI k,i and τk,i. Next a compact formulation is developed to simplify the notation of the problem. The state at time tk+j+1, given the state xkat time tk, and the input signals from tkto time tk+j, is computed by applying recursively, in the PAM and impulsive cases, by applying Equation (16): xk+j+1 =Ak+j+1,kxk+ k+j X i=k Ak+j,iBk,iUk,i,(19) where the definition Ak,k = I, Ak+1,k =Ak, and if j > 0 then Ak+j+1,k = Ak+jAk+j−1. . . Akhas been used. Define now Xkand Ukas a stack vector of Npstate and input vectors, respectively, spanning from time tk+1 to time tk+Np for the state and from time tkto time tk+Np−1for the controls, where Npis the planning horizon: Xk=     xk+1 . . . xk+Np      ,Uk=     Uk . . . Uk+Np−1      . Similarly, for the impulsive and PWM cases, define Γk=     τk . . . τk+Np−1      ,Λk=     κk . . . κk+Np−1      ,Υk=  Γk Λk  Then one can write Xk=Fkxk+GkUk,(20) where Gkis a square, block lower triangular matrix of size mNp, defined as Gk=         Bk0··· 0 Ak+2,k+1BkBk+1 ··· 0 . . .. . ..... . . Ak+Np,k+1BkAk+Np,k+2Bk+1 . . . Bk+Np−1         ,(21) 8 this is, its non-null blocks are defined by (Gk)jl =Ak+j,k+lBk+l−1, and the matrix Fkis defined as: Fk=         Ak+1,k Ak+2,k . . . Ak+Np,k         .(22) It is important to note that, in the impulsive case, Gkis a (nonlinear) function of Γk, whereas in the PWM case it is a (nonlinear) function of Υk. To avoid lengthy expressions this dependence has been omitted. Another important remark is that, in the PWM case, Ukis fixed whereas in the other cases is the input variable 3. Formulation of the planning problem Next the planning problem is formulated, introducing the constraints and the objective function. The formulation is done for the three types of control signals. 3.1. Constraints on the problem First constraints on the state and input are introduced. While only inequality are considered, equality constraints would be treated similarly. 3.1.1. Inequality constraints on the state In this work it is assumed that the state is subject to inequality constraints along the planning horizon, which can vary as time advances. These can be formulated in general as AkXk≤bk, and using (20), one reaches a expression in term of inputs, namely AkGkUk≤bk−AkFkxk.(23) 9 compute numerically if an explicit solution for the integrals (10) is not known or possible. 4.4. Steps 4 and 5. Refined PWM solution: An optimization algorithm To linearize (20) around inputs Υi k, notice from (8) that ∂ ∂τk,i BW k,i(τk,i, κk,i) = ∂ ∂τk,i Ztk+τk,i+κk,i tk+τk,i Φ(tk+1, s)Bi(s)ds, = Φ(tk+1, tk+τk,i +κi)Bi(tk+τk,i +κk,i) −Φ(tk+1, tk+τk,i)Bi(tk+τk,i),(46) and ∂ ∂κk,i BW k,i(τk,i, κk,i) = ∂ ∂κk,i Ztk+τk,i+κk,i tk+τk,i Φ(tk+1, s)Bi(s)ds, = Φ(tk+1, tk+τk,i +κi)Bi(tk+τk,i +κk,i).(47) Thus, (16) can be explicitly linearized around some given ¯τ+ k, ¯τ− kand ¯κ+ k, ¯κ− k, reaching xk+1 =Akxk+BkUk+B∆τ kδτk+B∆κ kδκk,(48) where Bkis computed with ¯τ+ k, ¯τ− k, ¯κ+ k, and ¯κ− k, and define B∆τ k=hBδτ+ kBδτ− ki, B∆κ k=hBδκ+ kBδκ+ ki,(49) δτk=  ¯τ+ k−τ+ k ¯τ− k−τ− k , δκk=  ¯κ+ k−κ+ k ¯κ− k−κ− k ,(50) where the i-th entries of the Bδ kmatrices in (50) are given, respectively, by (Bδ¯τ+ k)i= Φ(tk+T, tk+ ¯τ+ k,i + ¯κ+ k,i)Bi(tk+ ¯τ+ k,i + ¯κ+ k,i)uW+ k,i −Φ(tk+T, tk+ ¯τ+ k,i)Bi(tk+ ¯τ+ k,i)uW+ k,i (51) (Bδ¯τ− k)i=−Φ(tk+T, tk+ ¯τ− k,i + ¯κ− k,i)Bi(tk+ ¯τ− k,i + ¯κ− k,i)uW− k,i +Φ(tk+T, tk+ ¯τ− k,i)Bi(tk+ ¯τ− k,i)uW− k,i (52) (Bδ¯κ+ k)i= Φ(tk+T, tk+ ¯τ+ k,i + ¯κ+ k,i)Bi(tk+ ¯τ+ k,i + ¯κ+ k,i)uW+ k,i ,(53) (Bδ¯κ− k)i=−Φ(tk+T, tk+ ¯τ− k,i + ¯κ− k,i)Bi(tk+ ¯τ− k,i + ¯κ− k,i)uW− k,i (54) 16 Thus, defining stack vectors with the increments in the PWM variables at step ias ∆Γi k=     δτk . . . δτk+Np−1      ,∆Λi k=     δκk . . . δκk+Np−1      , and grouping all increments as ∆i k=  ∆Γi k ∆Λi k , B∆ k=  B∆τ k B∆κ k , and defining G∆ kas in (55), i.e., G∆ k=         B∆ k0··· 0 Ak+2,k+1B∆ kB∆ k+1 ··· 0 . . .. . ..... . . Ak+Np,k+1B∆ kAk+Np,k+2B∆ k+1 . . . B∆ k+Np−1         ,(55) one can write Xi k≈Fkxk+Gk(Υi k)Uk+G∆ k(Υi k)∆i k,(56) The state constraints (85) become AkG∆ k(Υi k)∆i k≤bk−AkFkxk−Gk(Υi k)Uk.(57) The constraints on ∆Γi kand ∆Λi kare as follows: −AW∆i k≤bW−AWΥi k,(58) ∆max ≤∆i k≤∆max,(59) where the last constraint (59) is used to avoid large variations that might make the linearization approximation to fail. All these constraints are might be summarized as A∆∆i k≤b∆(60) 17 Finally, the objective function can be rewritten in terms of ∆i kas Ji k=Ji k(Υi k)+ J∆ k(Υi k,∆i k), where Ji k= 2xT kFT kQkGk(Υi k)Uk+UT kGT k(Υi k)QkGk(Υi k)Uk+αAJ kΥi k,(61) J∆ k= 2(xT kFT k+UT kGT k(Υi k))QkG∆ k(Υi k)∆i k+ (∆i k)TG∆T k(Υi k)QkG∆ k(Υi k)∆i k +αAJ k∆i k.(62) Noting that J∆ kis quadratic in ∆i k, a quadratic optimization program with linear restriction, formulated on the output increments, can be posed as follows: min ∆i k J∆ k(Υi k,∆i k) (63) s. t.: AkG∆ k(Υi k)∆i k≤bk−AkFkxk−Gk(Υi k)Uk., A∆∆i k≤b∆. The solution ∆kis used to recompute new PWM inputs, Υi+1 k= Υi k+ ∆i k. Then the linearization process can be repeated around the new Υi+1 k, refining the solution in each iteration. 5. Model Predictive Control with PWM inputs In this section, building upon the trajectory planning algorithm of Section 4, which is open-loop and has a finite time-horizon, a closed-loop algorithm is developed based on the ideas of model predictive control (also known as receding horizon control). Model predictive control closes the loop by simply re-planning the maneuver at each time step, after applying just the set of control inputs corresponding to the first time step, and keeps looking ahead Nptime steps. Thus, the algorithm starts at k= 0 and is repeated for each k. The re-planning is done from the actual position at each time step, which seldom coincides with the planned position due to disturbances, thus effectively closing the loop. However, except at the start, it is not necessary to repeat all the steps of Section 4. Since the new position should be close to the planned one, one can apply the linearization scheme of the planning algorithm starting from the last available linearization. The MPC algorithm is summarized next: 18 Step 1. At time step k= 0 and starting from x0apply the Planning algorithm of Section 4, obtaining a set of PWM inputs Υ0that would optimize the planning problem (38) for the next Nptime steps, if there were no disturbances. Step 2. Apply impulses corresponding to the first time instant; save the rest of impulses. Set k= 1 Step 3. One arrives at xk, which probably is not the intended value of the state at time kbut close. Thus re-planning is necessary. Step 4. For re-planning, apply the planning algorithm of Section 4. However, to avoid the initial step of having to use a PAM or impulsive model and compute an initial guess, use instead as an initial guess the impulses of Υk−1that were not used (all of them except those corresponding to time k−1) and guess the remaining impulses (at the end) as zeros. In this way, form an initial guess Υ0 k. Step 4. Apply the linearization algorithm of Section 4.4 using Υ0 kas initial guess to obtain, after iterating, a new set of impulses Υk. Apply the set of impulses corresponding to time k. Save the rest of impulses. Step 5. Repeat step 3. 6. Example application: Spacecraft Rendezvous Rendezvous of spacecraft is the controlled close encounter of two (or more) space vehicles. This work assumes just two vehicles, one of which is the target vehicle (which is in a known orbit, and considered passive) and the other is the chaser spacecraft, which begins from a known position and maneuvers until very close to target. Only close range rendezvous Fehse (2003) is considered, which starts at hundreds of meters and ends when the chaser is very close to target (a few meters with speeds of centimeters per second). 19 There are numerous mathematical models for spacecraft rendezvous; which one should be used depends on the parameters of the scenario. In Carter (1998) a survey of numerous mathematical models for spacecraft rendezvous can be found. For instance, if the target is orbiting in a circular Keplerian orbit, the general equations of the relative movement between an active chaser spacecraft close to a passive target vehicle are linear time-invariant Hill-Clohessy-Wiltshire (HCW) equations (introduced in Hill (1878) and Clohessy and Wiltshire (1960)). While these equations are frequently used in the literature, it must be noted that, in many situations, the HCW equations are not accurate. For instance, if the target vehicle is moving in a Keplerian eccentric orbit (see Inalhan et al. (2002)) or if some orbital perturbations are taken into account (see for example Humi and Carter (2008)). A more complex model, the Tschauner-Hempel model (see Tschauner and Hempel (1965) or Carter (1998)) assumes that the target vehicle is passive and moving along an elliptical orbit with semi-major axis a and eccentricity e. The system equations are linear time-varying and cannot be exactly integrated in time to obtain a discrete transition model; however, if one substitutes the time tby the eccentric anomaly of the target orbit, E, it is possible to obtain explicit expressions for the system evolution in the PWM, impulsive, and PAM actuation cases. This will be the model considered in this work. The model can be expressed in cartesian coordinates, but also in the so-called relative orbital elements (see, e.g., Gaias et al. (2014) or Sinclair et al. (2014)). The former has been chosen for simplicity. Let us first establish some notation. Define the orbital mean motion n= pµ a3, where µis the gravitational parameter of the central body around which the target spacecraft is orbiting. Now, note that tand Eare related in a one-to-one fashion by using Kepler’s equation: n(t−tp) = E−esin E, (64) where tpis the time at periapsis used as a starting point to measure the ec- 20 Figure 2: LVLH frame. centric anomaly E. The time tpis chosen such that it is equal or less than the starting time which is denoted as t0(subtracting, if necessary, any number of orbital periods). Kepler’s equation is not analytically invertible, but its inverse can be found numerically with any desired degree of precision (see any Orbital Mechanics reference, such as Wie (1998)). Denote its inverse by the function K, i.e. E=K(t). Denote by E0the true anomaly corresponding to t0, this is, E0=K(t0). Then, Ek=K(tk) = K(t0+kT), where Tis the sampling time (not to be confused with the orbital period). Call as rx,k,ry,k, and rz,k 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 at time tk. In the (elliptical) LVLH frame, xrefers to the radial position, zto the out-of-plane position (in the direction of the orbital angular momentum), and yis perpendicular to these coordinates (not necessarily aligned with the target velocity given that its orbit is not circular). The velocity and inputs of the chaser in the LVLH frame at time tkare denoted, respectively, by vx,k,vy,k, and vz,k, and by ux,k,uy,k, and uz,k. If there is no actuation (i.e. ux,k =uy,k =uz,k = 0), the resulting transition 21 equation was obtained in a simple form in Yamanaka and Ankersen (2002) as follows: xk+1 = Φ(tk+1, tk)xk(65) where xk= [rx,k ry,k rz,k vx,k vy,k vz,k]T,(66) and where Φ(tk+1, tk) = YK(tk+1)Y−1 K(tk),(67) with Ytkbeing the fundamental matrix solution of the Tschauner-Hempel model, which are expressed in Yamanaka and Ankersen (2002) as a function of true anomaly θ. However there is a one-to-one relation between Eand θgiven by tan θ 2=r1 + e 1−etan E 2,(68) which is exploited in the sequel. The explicit expression of the matrices1is found in (69) and (70), YE=         s0 0 2/ρ −3esJ −c0 c(1 + 1/ρ) 1/ρ 0−3ρJ s(1 + 1/ρ) 0 0 0 c/ρ 0 0 s/ρ αρ2c0 0 α(−es −3eρ2Jc)αρ2s0 αs(−1−ρ2)αes 0αρ(3esρJ −3) α(c+e+cρ2) 0 0 0 −sα 0 0 (c+e)α         ,(69) 1These expressions slightly differ from Yamanaka and Ankersen (2002) because the two transformation matrices that appear in that paper have been pre-multiplied; also, the reference axes are not the same. 22 Y−1 E=3J (1 −e2)         eρ2(1 + ρ)−e2ρ2s0e2s/α eρ/α 0 ρ2(1 + ρ)−eρ2s0es/α ρ/α 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         +1 (1 −e2) ×         −s(ρ2+ 2ρ+e2)es2(1 + ρ) 0 c−2e/ρ α−s(ρ+1) ρα 0 −es(1 + ρ)2ρ2(1 −ce) + e2s20ec−2/ρ α−es(ρ+1) ρα 0 0 0 (c+e)(1 −e2) 0 0 −s(1−e2) αρ ρ2(1 + ρ)−esρ20es α ρ α0 3ρ(c+e)−eρs2−esc(1 + ρ)−e2s0s α c(ρ+1)+e αρ 0 0 0 s(1 −e2) 0 0 c(1−e2) αρ         (70) where the following symbols are used (expressed in terms of E): ρ=1−e2 1−ecos E, s =√1−e2sin E 1−ecos E, c =cos E−e 1−ecos E,(71) J=αE−ˆ E−e(sin E−sin ˆ E) (1 −e2)3/2,α =n (1 −e2)3/2,(72) where ˆ Ein (72) can be substituted by zero or any other desired reference value of E. For instance, if when evaluating (67) one chooses ˆ E=Ek=K(tk), then for Y−1 K(tk)one gets J= 0 and the first matrix in (70) becomes zero. Using (67), one gets Akin (16) explicitly, as well as Bkin the impulsive case (explicitly defined in in terms of Φ(tk+1, tk)). To obtain the Bkmatrix in the PWM and PAM cases, one needs to solve (8) or (12), respectively, which involves an integral. Defining bi(r1, r2, r3) = Zr2 r1 Φ(r3, s)Bi(s)ds, (73) one has that, from (8), BW k,i(τk,i, κk,i) = bi(tk+τk,i, tk+τk,i +κk,i, tk+1),(74) and, from (12), BA k,i =bi(tk, tk+1, tk+1).(75) To compute the bi’s, the following integral is needed ZY−1 K(t)Ci+3dt, (76) 23 where Ciis a 6-element column vector of zeros with a value of one at row i. For the computation, define the functions fi(t), for i= 1,2,3, as the indefinite integrals of (76), in terms of eccentric anomaly fi(E) = ZY−1 ECi+3 1−ecos E ndE. (77) Once the fi’s are computed, one finds the bi’s as bi(r1, r2, r3) = YK(r3)(fi(K(r2)) −fi(K(r1))) (78) Inserting the expression of (70) in (77) and integrating, one obtains f1=1−e2−7/2 2α2         2(1 + 6e2)S−3e(2 + e2)E+e2Ch +e3ˆ S/2 2e(8 −e2)S−(4 + 7e2−2e4)E+eCh +e2ˆ S/2 0 −2eC 1−e23/2 −2C1−e23/2 0         ,(79) f2=1−e2−3 2α2           C(4(1 + e2)−eC)−eEh −3eE2 eC(10 −2e2−eC)−Eh −3E2 0 2E1−e23/2 √1−e2(4S−e(3E+ˆ S/2)) 0           ,(80) f3=1−e2−5/2 4α2         0 0 2√1−e2C(2 −eC) 0 0 4e2+ 1S−e(6E+ˆ S)         ,(81) where S= sin(E), ˆ S= sin(2E), C= cos E,h= 6α(ˆ E−E−esin( ˆ E)). Similar expressions for the Bmatrices can be found in Ankersen (2010), however using a slightly different definition of reference axes. Note that, using these formulas, it is possible to express (16) explicitly for all actuation types. This greatly speeds up the algorithms. 6.1. Constraints for the rendezvous problem Besides the input constraints (which were given in Section 3.1.2), the inequality state constraints which were generically specified in Section 3.1.1 are, 24 in general, related to safety and sensing considerations (see e.g. Breger and How (2008)). In this work, it is considered that during rendezvous the chaser vehicle has to remain inside a line of sight (LOS) area. To simplify the constraint2, in this work a 2-D LOS area is used as shown in Figure 3. This LOS region is the intersection of a cone, given by the equations ry≥cLOS(rx−rx0) and ry≥ −cLOS(rx+rx0), and the region ry≥0. Figure 3: Line of Sight region. The LOS constraint is ALOSxk≤bLOS, where ALOS =     0−10000 cLOS −10000 −cLOS −10000      , bLOS =     0 cLOSrx0 cLOSrx0      .(82) Using the compact formulation that was developed in Section 2, the con- 2More complicated constraints could be considered, see Gavilan et al. (2012) for examples including a rotating LOS constraint. 25 Ieko, T., Ochi, Y., Kanai, K., 1999. New design method for pulse-width modulation control systems via digital redesign. J. Guid. Contr. Dynam. 22 (1), 123–128. Inalhan, G., Tillerson, M., How, J. P., 2002. Relative dynamics and control of spacecraft formations in eccentric orbits. J. Guid. Contr. Dynam. 25 (1), 48–59. Jewison, C., Erwin, R. S., Saenz-Otero, A., 2015. Model predictive control with ellipsoid obstacle constraints for spacecraft rendezvous. In: ACNAAV 2015 IFAC workshop. Kim, H. J., Shim, D. H., Sastry, S., 2002. Nonlinear model predictive tracking control for rotorcraft-based unmaned aerial vehicles. In: Proceedings of ACC 2002. Larsson, R., Berge, S., Bodin, P., J¨onsson, U., 2006. Fuel efficient relative orbit control strategies for formation flying and rendezvous within prisma. In: Proceedings of the 29th AAS guidance and control conference. Leomanni, M., Rogers, E., Gabriel, S. B., 2014. Explicit model predictive control approach for low-thrust spacecraft proximity operations. Journal of Guidance, Control, and Dynamics 37 (6), 1780–1790. Louembet, C., Arzelier, A., Deaconu, G., 2015. Robust rendezvous planning under maneuvering errors. Journal of Guidance, Control, and Dynamics 38 (1), 76–93. Richards, A. G., How, J., 2003. Performance evaluation of rendezvous using model predictive control. AIAA Paper 2003-5507. Rossi, M., Lovera, M., 2002. A multirate predictive approach to orbit control of small spacecraft. In: Proceedings of ACC 2002. Rugh, W. J., 1996. Linear System Theory (2Nd Ed.). Prentice-Hall, Inc., Upper Saddle River, NJ, USA. 32 Shieh, L.-S., Wang, W.-M., Sunkel, J., 1996. Design of PAM and PWM controllers for sampled-data interval systems. J Dyn Syst Meas Contr. 118 (4), 673–681. Sinclair, A. J., Sherrill, R. E., Lovell, T. A., 2014. Calibration of linearized solutions for satellite relative motion. Journal of Guidance, Control, and Dynamics 37 (4), 1362–1367. Tschauner, J., Hempel, P., 1965. Rendevous zu einem in elliptischer bahn umlaufenden. Ziel. Acta Astronaut. II (2), 104–109. Vazquez, R., Gavilan, F., Camacho, E. F., 2011. Trajectory planning for spacecraft rendezvous with on/off thrusters. In: Proc. of IFAC World Congress 2011. Vazquez, R., Gavilan, F., Camacho, E. F., 2014. Trajectory planning for spacecraft rendezvous in elliptical orbits with On/Off thrusters. In: IFAC World Congress, Cape Town. Vazquez, R., Gavilan, F., Camacho, E. F., 2015. Model predictive control for spacecraft rendezvous in elliptical orbits with On/Off thrusters. In: ACNAAV 2015 IFAC workshop. Weiss, A., Kolmanovsky, I., Baldwin, M., Erwin, R. S., 2012. Model predictive control of three dimensional spacecraft relative motion. In: American Control Conference (ACC), 2012. IEEE, pp. 173–178. Wie, B., 1998. Space vehicle dynamics and control. AIAA. Yamanaka, K., Ankersen, F., 2002. New state transition matrix for relative motion on an arbitrary elliptical orbit. J. Guid. Contr. Dynam. 25 (1), 60–66. 33 List of Figure Captions 1 PWMVariables. ........................... 6 2 LVLHframe. ............................. 21 3 LineofSightregion. ......................... 25 4 System trajectories in the target orbital plane: open-loop PWM inputs computed from impulsive solution (dashed), closed-loop Model Predictive Control with PWM inputs using impulsive model (dot-dashed), and closed-loop Model Predictive Control with PWM inputs using the PWM planning algorithm (solid). . . . . . . . . 27 5 System trajectories in the target orbital plane, with inexact orbit model: open-loop PWM inputs computed with the planning algorithm (dashed), closed-loop Model Predictive Control with PWM inputs using impulsive model (dot-dashed), and closedloop Model Predictive Control with PWM inputs using the PWM planning algorithm (solid). . . . . . . . . . . . . . . . . . . . . . . 28 34