Chance-constrained Model Predictive Control for Near Rectilinear Halo Orbit spacecraft rendezvous

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Elsevier in Aerospace Science and Technology, Vol. 100 on May 2020, available at: https://doi.org/10.1016/j.ast.2020.105827 © 2020 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND Chance-constrained Model Predictive Control for Near Rectilinear Halo Orbit Spacecraft Rendezvous Julio C. Sancheza, Francisco Gavilana, Rafael Vazqueza,∗ aDepartamento de Ingenier´ıa Aeroespacial, Escuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, 41092, Sevilla, Spain Abstract This work presents a robust Model Predictive Controller (MPC) to solve the problem of spacecraft rendezvous in the context of the restricted three-body problem (R3BP) as will be required to dock with space stations in cislunar space. The employed methodology is both valid for chemical and electric thrusters. By exploiting the state transition matrix and using a chance-constrained approach, the robust MPC assures constraints satisfaction under the presence of disturbances in a probabilistic sense. The perturbations parameters are computed on-line using a disturbance estimator. The robust controller is tested for a rendezvous scenario with a target placed in an Earth-Moon L2 Near-Rectilinear Halo Orbit. Numerical results are shown and discussed. Keywords: Spacecraft rendezvous, Three body problem, Model predictive control, Robust control. 1. Introduction Demonstrating rendezvous capabilities in the context of multi-body environments is becoming a growing and active field of research as International Space Station (ISS) partners have interest in building a space station in the cislunar space, named as the Lunar Orbital Platform Gateway (LOP-G), see [1]. Moreover, this lunar space station will greatly enhance scientific opportunities by allowing to return samples from the Moon, see [2]. ∗Corresponding author Email address: [email protected] (Rafael Vazquez) Accepted version in Aerospace Science and Technology May 1, 2020 Several options have been studied to place the LOP-G, see [3], being the Near Rectilinear Halo Orbits (NRHOs), around the L2Earth-Moon point, the most attractive candidates. NRHOs are members of the broader set of L1and L2families of Halo orbits existing in the circular restricted three-body problem (CR3BP), see [4] for more details about CR3BP orbits. The NRHOs also persist in higher-fidelity models since they present favourable stability properties, see [5]. Typically, far-rendezvous operations, where fuel consumption is the key driver instead of safety considerations, have been extensively studied in the literature. Reference [6] exploits the method of invariant manifolds connections whereas surrogate models, to ease the computational burden of global optimization, have been proposed by [7]. Finally, [8] compared the fuel efficiency of classical phasing strategies with invariant manifolds connections. On the other hand, close rendezvous operations (where safety is a main concern) are starting to gain more momentum. Reference [9] proposed a targeting law combined with a navigation filter for restricted three body problem (R3BP) rendezvous operations. Practical rendezvous scenarios for Earth-Moon Halo orbits were proposed in [10], whereas shooting methods to achieve rendezvous have been studied in [11]. The previous works have expressed the system dynamics in the Earth-Moon co-rotating reference frame. However, this frame is not very useful to describe state constraints attached to the target. This is the reason why local frames are widely preferred in close rendezvous operations, see [12]. In [13], a local frame of reference is proposed taking into account that the LOP-G will be orbiting the Moon in a practical sense. The purpose of this work is to develop a robust rendezvous controller for R3BP scenarios. The key idea behind robust control is to explicitly take into account disturbances and uncertainties in the optimization problem. In the case of Keplerian rendezvous operations, several robust techniques have been explored. Reference [14] employed the chance-constrained approach to guarantee constraints satisfaction probabilistically. A 2 worst-case scenario methodology, to minimize the size of the terminal arrival set, was proposed by [15]. Finally, a tube-based method, guaranteeing constraint satisfaction for bounded disturbances, has been experimentally validated in [16]. The main contribution of this work is the extension of the chance-constrained approach, developed in [14], to R3BP rendezvous. The proposed method explicitly considers the disturbances, affecting the state constraints, in a probabilistic sense. Then, these probabilistic constraints are bounded at a certain probability, which allows to compute control signals in a deterministic way. Since a priori knowledge of the disturbances statistical properties is required, an on-line estimator of stochastic parameters is also employed. The robust program is embedded into a Model Predictive Control (MPC) scheme, see [17], so the robust program is updated after each sampling time. Moreover, for this type of mission, the propulsive plant of the chaser can be either chemical or electrical. To extend the potential application of this work, both the impulsive and continuous thrust models are considered. For the continuous thrust case, it is assumed that the control signal can be linearly parameterized by some decision variables. As an additional contribution, basis splines (B-splines), typically employed for attitude control as in [18] and [19], are chosen to parameterize the control signal. The structure of this work is as follows. Section 2 describes motion in the restricted three body problem and the linearized relative model. Section 3 follows describing the rendezvous problem. Section 4 formulates the chance-constrained based MPC and the on-line disturbance estimator. Section 5 shows numerical results through a Monte Carlo comparison of the robust and non-robust controllers. Section 6 closes the paper with some final remarks. 2. Relative motion in the restricted three body problem This section studies the relative motion between two vehicles in the R3BP. Firstly, the motion of a particle, under R3BP assumptions, is described. Additionally, some 3 facts about NRHOs are given. Then, the local-vertical local-horizontal (LVLH) frame is introduced and the R3BP relative dynamics deduced. Finally, the relative motion is linearized assuming that the vehicles are close enough. 2.1. Restricted three body problem and NRHOs Under R3BP assumptions, where µ1≥µ2≫µ, being µ1and µ2the gravitational parameters of the two primaries and µthat of the vehicle, the spacecraft dynamics are conveniently expressed in the synodic frame, see [20]. Denote the inertial frame by I:{O,iI,jI,kI}where Ois the position of the system barycenter. Denote the synodic frame by S:{O,iS,jS,kS}, with iScoincident with the line uniting the two primaries and positive in the direction of the second primary, kSparallel to the system kinetic momentum and jScompleting a right-handed system, see Fig.1. The R3BP equations in Figure 1: Inertial, synodic and LVLH frames of reference for the Earth-Moon system. the Sframe are ¨r|S=−µ1(r−r1) ∥r−r1∥3 2 −µ2(r−r2) ∥r−r2∥3 2 −2ω ω ωS/I ×˙r|S−˙ω ˙ω ˙ωS/IS×r−ω ω ωS/I ×(ω ω ωS/I ×r) + u, (1) where ris the spacecraft position, r1and r2the primaries position, ω ω ωS/I the angular velocity of the synodic frame with respect to the inertial and uthe control acceleration. Eq.(1) allows primaries in elliptic orbits. To obtain the CR3BP equations (circular orbits), set ω ω ωS/I =nkSand ˙ω ˙ω ˙ωS/I =0in Eq.(1), obtaining ¨r|S=−µ1(r−r1) ∥r−r1∥3 2 −µ2(r−r2) ∥r−r2∥3 2 −2nkS×˙r|S−nkS×(nkS×r) + u,(2) 4 where n=p(µ1+µ2)/D3and Dis the distance between the two primaries. The CR3BP system (2) has five libration points, named as Lagrange points (Li, i = 1 . . . 5), with associated families of periodic orbits around them, see [4]. Amongst these periodic orbits, the ones receiving more attention, for practical purposes, are the Halo orbits around collinear equilibria. Since these are unstable, the Halo orbits are in turn inherently unstable, requiring station-keeping to be maintained. Amongst each set of L1 and L2Halo orbits, there exists a subset (NRHOs) with favourable stability properties. These properties have shown to persist in higher-fidelity models, and hence these orbits may support long-term missions near the Moon. Regarding scientific opportunities, the preferred Earth-Moon NRHOs are the ones associated to the Southern L2family. This family allows great coverage for both the lunar South pole and far side of the Moon, see [21]. Covering these areas is of great scientific interest due to the existence of water ice in the South pole, see [22], and the impossibility to observe the far side of the Moon from Earth. The Southern L2Halo family and their subset of NRHOs for the Earth-Moon system are shown in Fig.2 in a non-dimensional synodic frame. Note that they can be practically seen as lunar orbits, with the perilune at the North pole. Figure 2: Green: Southern L2Halo family; blue: Southern L2NRHOs; black: Sec.V NRHO. Parameter ais the Earth-Moon semimajor axis. To evaluate the stability properties of CR3BP periodic orbits, [23] proposed the 5 stability index parameter ν ν=1 2λmax +1 λmax ,(3) which is a function of λmax, the absolute value of the monodromy matrix (state transition matrix after one orbital period) maximum real eigenvalue (in absolute value). The monodromy matrix of an autonomous Hamiltonian system is symplectic, hence each eigenvalue λhas an opposite one λ−1, see [24] for the details. Since the orbit is periodic, two monodromy matrix eigenvalues are always equal to the unity. As a consequence ν≥1 and the periodic orbit is marginally stable if ν= 1 and unstable if ν > 1. Both the stability indexes and orbital periods for the Southern L2NRHOs are shown in Fig.3. It can be seen that at very close distances from the Moon surface, the NRHOs are almost marginally stable. As distance from the Moon increases, the stability indexes rise and decrease until they become almost marginally stable for altitudes ranging from 11500 km to 16750 km. Afterwards the stability indexes begin to increase quickly becoming highly unstable. Additionally, Fig.3 shows the period increases monotonically with respect to the perilune radius. As remarked by [23], some practical orbits exist within the EarthMoon NRHOs. A 9:2 resonance with the Moon synodic period (∼29.5 days) can be found at an altitude of ∼1500 km, whereas another 4:1 resonance arises at ∼4150 km, which are useful to avoid Earth eclipses at all times. 0 5000 10000 15000 0 1 2 3 4 5 6 7 8 9 10 11 Figure 3: Stability indexes and periods for Southern L2NRHOs 6 2.2. Relative motion in the R3BP For relative dynamics, following [13], a local frame (LVLH) is employed. The frame is denoted by L:{rt,iL,jL,kL}, where rtis the target position, kLis pointing towards the second primary, jLis in the opposite direction to the target kinetic momentum as view from the Sframe with respect to the second primary and iLcompletes the right-handed frame. Fig.1 shows the Lframe as well as the target position rt, the chaser position r and the relative position ρ ρ ρ=r−rt. The relative dynamics in the Lframe is given by ¨ρ ¨ρ ¨ρ|L=¨r|L−¨rt|L,which can be further developed by using Eq.(1), reaching ¨ρ ¨ρ ¨ρ|L=−2ω ω ωL/I ×˙ρ ˙ρ ˙ρ|L−ω ω ωL/I ×(ω ω ωL/I ×ρ ρ ρ) −˙ω ˙ω ˙ωL/I L×ρ ρ ρ−µ1ρ ρ ρ+r1t ∥ρ ρ ρ+r1t∥3 2 −r1t ∥r1t∥3 2 −µ2ρ ρ ρ+r2t ∥ρ ρ ρ+r2t∥3 2 −r2t ∥r2t∥3 2+u, (4) where r1t=rt−r1and r2t=rt−r2denote the relative position of the target with respect to the primaries. Note ω ω ωL/I =ω ω ωL/S +ω ω ωS/I,(5) ˙ω ˙ω ˙ωL/IL= ˙ω ˙ω ˙ωL/SL+ ˙ω ˙ω ˙ωS/IS−ω ω ωL/S ×ω ω ωS/I,(6) hence, rt,ω ω ωL/S and ˙ω ˙ω ˙ωL/SLdepend on the target motion with respect to the synodic frame whereas r1,r2,ω ω ωS/I and ˙ω ˙ω ˙ωS/ISdepend on the primaries motion (Eq.(4) is still valid if the primaries evolve in elliptic orbits). 2.3. Linearized relative motion in the R3BP Considering close-range rendezvous operations, that is, ∥r1t∥2,∥r2t∥2≫ ∥ρ ρ ρ∥2, one has r ∥r∥3 2 ≈r0 ∥r0∥3 2 −1 ∥r0∥3 2I−3r0rT 0 ∥r0∥2 2(r−r0),(7) 7 being r0the linearization point. Introducing the linearization of Eq.(7) into Eq.(4), one obtains ¨ρ ¨ρ ¨ρ=−˙ Ω ˙ Ω ˙ ΩL/I + Ω Ω Ω2 L/I −µ1 r3 1tI−3r1trT 1t r2 1t−µ2 r3 2tI−3r2trT 2t r2 2tρ ρ ρ−2Ω Ω ΩL/I ˙ρ ˙ρ ˙ρ+u, (8) where Ω Ω ΩL/I and ˙ Ω ˙ Ω ˙ ΩL/I are the cross-product matrices associated to ω ω ωL/I and ˙ ω ω ωL/I respectively, see [25]. This can be written as a linear time-varying (LTV) system: d dt     ρ ρ ρ ˙ρ ˙ρ ˙ρ     =    0 I A˙ρ ˙ρ ˙ρρ ρ ρ−2Ω Ω ΩL/I         ρ ρ ρ ˙ρ ˙ρ ˙ρ     +    0 I     u,(9) where A˙ρ ˙ρ ˙ρρ ρ ρ=−˙ Ω ˙ Ω ˙ ΩL/I + Ω Ω Ω2 L/I −µ1 r3 1tI−3r1trT 1t r2 1t−µ2 r3 2tI−3r2trT 2t r2 2t.(10) Defining x= [ρ ρ ρT,˙ρ ˙ρ ˙ρT]T, Eq.(9) is of the form ˙x(t) = A(t)x(t) + Bu(t), which has as general solution, see [26], x(t) = ϕ ϕ ϕ(t, t0)x0+Zt t0 ϕ ϕ ϕ(t, τ)Bu(τ)dτ, (11) with ϕ ϕ ϕ(t, t0) the state transition matrix, verifying ˙ ϕ ˙ ϕ ˙ ϕ(t, t0) = A(t)ϕ ϕ ϕ(t, t0), ϕ ϕ ϕ(t0, t0) = I.(12) 3. Rendezvous planning problem Next, the control inputs are described and parameterized; then, the objective function and the constraints are described. Finally, the rendezvous problem is stated. 3.1. Control input In this work, both chemical and electric thrusters are considered; thus, u=uC+uE, where uCand uEdenote the chemical and electric accelerations respectively. For the chemical thrusters, the control signal can be described by impulses (i.e. instantaneous changes of velocity) lim ∆t→0Ztk+∆t tk uC(t)dt = ∆V(tk)δ(t−tk),(13) 8 The probability of constraint satisfaction, by adding the bounding term bδ, should be near one. This guarantees that the chaser remains within the LOS region for almost all perturbations. Considering that the disturbances are normally distributed, δ δ δ∼N6(¯ δ ¯ δ ¯ δ,Σ Σ Σδ), with known mean, ¯ δ ¯ δ ¯ δ, and covariance matrix, Σ Σ Σδ= Σ Σ ΣT δ≻0, the following relation holds (see [32] for more details) δ δ δ∼N6(¯ δ ¯ δ ¯ δ,Σ Σ Σδ)−→ (δ δ δ−¯ δ ¯ δ ¯ δ)TΣ Σ Σ−1 δ(δ δ δ−¯ δ ¯ δ ¯ δ)∼χ2(6),(34) where χ2(6) is a chi-square probability distribution with six degrees of freedom. Making the hypothesis that the statistical properties of the disturbances are time-invariant (quasisteady approach), Eq.(34) is valid at all times (δ δ δk+j−¯ δ ¯ δ ¯ δ)TΣ Σ Σ−1 δ(δ δ δk+j−¯ δ ¯ δ ¯ δ)∼χ2(6), j = 0 . . . N, (35) hence the following probabilistic relations hold P(χ2(6) ≤α) = p−→ (δ δ δk+j−¯ δ ¯ δ ¯ δ)TΣ Σ Σ−1 δ(δ δ δk+j−¯ δ ¯ δ ¯ δ)≤α, where finding αfrom a given p, the right side inequality is guaranteed with probability p. Then, the parameter pis the probability of constraint satisfaction and should be as close to unity as possible. The bounding term bδ(k) can be found by solving the following minimization problem for each row iof −ALSGk,δ denoted as ai(k) min δ δ δS (bδ(k))i=ai(k)δ δ δS, s.t. (δ δ δk+j−¯ δ ¯ δ ¯ δ)T(αΣ Σ Σδ)−1(δ δ δk+j−¯ δ ¯ δ ¯ δ)≤1, (36) It can be proved, see [14], that the rows of bδ(k) are (bδ(k))i= N X j=0 −qaijH−1aij +aij¯ δ ¯ δ ¯ δ.(37) 15 Once the vector bδ(k) is computed through Eq.(37), the control input at time tkis obtained by solving the following robust program min ∆VS(k),ξ ξ ξS(k)J(xk,∆VS(k),ξ ξ ξS(k),¯ δ ¯ δ ¯ δS(k)), s.t. ALS(Gk,∆V∆VS(k) + Gk,ξξ ξ ξS(k)) ≤bLS −ALSFkxk+bδ(k), −∆VS,max ≤∆VS(k)≤∆VS,max, −uS,max ≤Bk,Sξξ ξ ξS(k)≤uS,max, (38) which is a quadratic programming (QP) problem. 4.4. Disturbance estimator The robust satisfaction of constraints, presented in the section 4.3, requires a priori knowledge of the perturbations statistical properties, ¯ δ ¯ δ ¯ δand Σ Σ Σδ. However, such properties are typically unknown and they have to be estimated on-line. Since the disturbances have been assumed as normally distributed such that δ δ δ∼N6(¯ δ ¯ δ ¯ δ,Σ Σ Σδ), the normal distribution parameters ¯ δ ¯ δ ¯ δand Σ Σ Σδare estimated a posteriori at each time kby taking into account all past disturbances δ δ δi=xi+1 −ϕ ϕ ϕ(ti+1, ti)xi−Zti+1 ti ϕ ϕ ϕ(ti+1, τ)Bu(τ)dτ, (39) with i= 1 . . . k −1. The estimates of ¯ δ ¯ δ ¯ δand Σ Σ Σδat time k, based on disturbances up to k−1, are named as ˆ δ ˆ δ ˆ δkand ˆ Σ ˆ Σ ˆ Σk,δ, and following [14] one can use recursive formulas for their estimation as follows ˆ δ ˆ δ ˆ δk=e−λ γk (γk−1ˆ δ ˆ δ ˆ δk−1+δ δ δk−1), ˆ Σ ˆ Σ ˆ Σk,δ =e−λ γkγk−1ˆ Σ ˆ Σ ˆ Σk−1+ (δ δ δk−1−ˆ δ ˆ δ ˆ δk)(δ δ δk−1−ˆ δ ˆ δ ˆ δk)T, with ˆ δ ˆ δ ˆ δ0=0and ˆ Σ ˆ Σ ˆ Σ0,δ =0. 5. Results In this section, an application case of rendezvous with a target located in an EarthMoon NRHO is considered. A comparison between the proposed chance-constrained 16 MPC algorithm against a non-robust MPC is carried out. 5.1. Simulation model The non-linear R3BP relative dynamics given by Eq.(4) are used to obtain the numerical results of this section. As reported in [13], the position error between the linear and non-linear models increases faster at the NRHO perilune (∼40 m in 1 h) compared to its apolune (∼2 m in 1 h). The minimum and maximum distance between Earth and Moon are taken as ∥r12∥= 363104 km and ∥r12∥= 405696 km, whereas the primaries gravitational parameters are µ1= 398600.4 km3/s2and µ2= 4904.869 km3/s2. The manoeuvre is considered to take place when the distance between Moon and Earth is minimal. Apart from model mismatch, numerical integration is required to obtain the LTV transition matrix with Eq.(12), hence cumulative integration errors are expected to arise. Another source of disturbances is the computation of the target NRHO which is done with the continuation software AUTO (see [33]), using the CR3BP model, see Eq.(2). The target L2Southern NRHO is taken as the one with ν= 1.0120, T= 10.35 days and closest distance to the Moon surface of 15674 km, see Fig.3. Regarding the thrusters performance, in the same sense as [14], the real control inputs ∆Vreal = [∆Vx,∆Vy,∆Vz]Tand ureal = [ux, uy, uz]Tdo not match the computed control signals ∆Vand uE ∆Vreal =R(δθ δθ δθ)(∆V+δ δ δV),(40) ureal =R(δθ δθ δθ)(uE+δ δ δuE),(41) where Ris a rotation matrix and δθ δθ δθ ∼N3(δ δ δ¯ θ ¯ θ ¯ θ,Σ Σ Σδθ) is a vector of small random angles modelling imperfect alignment of thrusters, whereas δ δ δV∼N3(δ δ δ¯ V,Σ Σ ΣδV ) and δ δ δuE∼ N3(δ δ δ¯ uE,Σ Σ ΣδuE) are additive random noises to the impulse or electric thrust amplitude respectively. Note that Nndenotes a n-dimensional gaussian distribution. 17 5.2. Simulation results In this section, the previously designed robust controller performance is evaluated for each one of the thrusters configurations. The initial manoeuvre time is chosen at the instant when the target is closest to the Moon (perilune), thus potentially representing a lunar sample return scenario, see [2]. The simulations are done in MATLAB with Gurobi as the QP solver (see [34]). The state transition matrices are computed numerically by solving the ODE system (12) with the ode45 routine of MATLAB which implements a 4th order Runge-Kutta method with a variable time step. For the continuous thrust case, the second term of the right hand-side of Eq.(11) is computed with a trapezoidal method integration. Although numerical integrations augments the computational burden, especially the one concerning the transition matrix, in practice these matrices can be computed on ground and uplinked to the probe before starting the manoeuvre. The common conditions for both scenarios are shown in Table 1. Note that a docking sensor has a cone half-angle of t01d 7h 9m 10s tf1d 19h 9m 10s cy1/tan(π/6) cz1/tan(π/6) y05 m z05 m δ δ δ¯ θ ¯ θ ¯ θ[2.5◦, 2.5◦, 2.5◦]TΣ Σ Σδθ (2.5◦)2I Table 1: Global simulation conditions 30◦. The controller tuning parameters are taken, for both scenarios, as N= 40, γ= 106, α= 0.95 and λ= 0.25. On the other hand, the specific continuous thrust parameters are chosen as nu= 400, q= 4 and nc= 44, hence assuring C4continuity. Since the disturbances evolve stochastically, 100 random realizations of them, see Eq.(40)-(41), are simulated. By doing this, the proposed robust controller can be effectively compared with a non-robust one (δ δ δS=0). 18 Figure 5: Chaser 3D trajectory for the first random realization using the robust controller. 5.2.1. Impulsive scenario Consider the impulsive scenario defined by Table 2. The thrust level bias could r0[400, 200, -200]Tm v0[0.1, -0.1, 0.1]Tm/s ∆Vmax [0.1, 0.1, 0.1]Tm/s umax [0, 0, 0]Tm/s2 max(δ δ δ¯ V) 5·10−4·[1, 1, 1]Tm/s Σ Σ ΣδV (5·10−4)2Im2/s2 Table 2: Impulsive scenario conditions potentially influence the results. As a consequence, it is considered to be in the interval δ δ δ¯ V∈[−max(δ δ δ¯ V),max(δ δ δ¯ V)] with constant probability. Note that the continuous thrusters are not operative since their bounds are null. The robust controller simulation results are shown in Fig.5-7. For the sake of clarity, only the trajectory on the XZ plane is shown, since that projection depicts the most critical LOS constraints. Although, a close range rendezvous scenario is considered, the trajectory start diverging from the target up to 1.5 km approximately, (see Fig.5-6), which highlights the capability of the linear model to provide fair accuracy at distances above the typical rendezvous ones (<1 km). The final in-track impulses ∆Vx, see Fig.7, 19 Figure 6: XZ plane trajectories for all the random realizations using the robust controller. 0 2 4 6 8 10 12 -0.01 0 0.01 0 2 4 6 8 10 12 -5 0 5 10 10-3 0 2 4 6 8 10 12 -0.02 -0.01 0 Figure 7: Impulses for the first random realization using the robust controller. Blue: computed impulses; red: applied impulses. Figure 8: XZ plane trajectories for all the random realizations using the non-robust controller. 20 012345 10-4 0.5 0.55 0.6 Figure 9: Mission cost against the thrust level bias for all the random realizations. Blue: LOS satisfaction; red: LOS violation. are positive to brake the chaser and avoid collision with the target. Two critical moments happen along the manoeuvre, the first one taking place just after the departure and the other one at the end of the rendezvous operation, see the zoomed areas of Fig.6 and Fig.8. Note that the non-robust controller is not capable of guaranteeing LOS constraint satisfaction in any case, see Fig.8, whereas the robust controller avoid the two arising conflicts for the 93% of the cases, see Fig.6. However, in exchange for the safeness increment, the mission cost also increases when comparing the robust approach with the non-robust one, see Fig.9. The computation times, for a i7-860 CPU at 2.80 GHz, are of 2.3451 s to compute the stack matrices whereas each MPC step requires 0.4606 s in average requiring the worst case 0.6911 s. 5.2.2. Continuous thrust scenario Consider the continuous thrust scenario characterized by Table 3. Again, the thrust level bias has been considered to vary with constant probability within the interval ¯ uE∈ [−max(¯ uE),max(¯ uE)]. The robust controller simulation results are shown in Fig.10-12. Again, the final control in the in-track direction, ux, is positive to avoid collision with the target, see 21 Figure 10: Chaser 3D trajectory for the first random realization using the robust controller. Figure 11: XZ plane trajectories for all the random realizations using the robust controller. 0 2 4 6 8 10 12 -1.5 -1 -0.5 0 0.5 110-5 Figure 12: Thrust acceleration for the first random realization using the robust controller. Solid: computed thrust; dotted: applied thrust. 22 r0[600, 300, -200]Tm v0[0.1, -0.1, 0]Tm/s ∆Vmax [0, 0, 0]Tm/s umax [10−4, 10−4, 10−4]Tm/s2 max(¯ uE) 5·10−7[1, 1, 1]Tm/s2 Σ Σ ΣδuE(5·10−7)2Im2/s4 Table 3: Continuous thrust scenario conditions 012345 10-7 2 4 6 8 10 12 10-6 Figure 13: Mission cost against the thrust level bias for all the random realizations. Blue: LOS satisfaction; red: LOS violation. Fig.12. Comparing the robust against the non-robust controller yields again the conclusion that the non-robust one shows worst performance in terms of constraints satisfaction when compared to the chance-constrained method, see Fig.11. As a matter of fact LOS constraint satisfaction is of 76% (appearing most of the violations at high bias levels, see Fig.13) for the robust controller whereas the non-robust controller achieves a LOS constraint satisfaction of 1%. Moreover, in this case, the chance-constrained method consumes less, in general, than the non-robust method, see Fig.13. Regarding computation times, 6.3504 s are required to compute the stack matrices at the beginning whereas in 23 average each robust MPC step takes 2.0369 s with the most severe computation requiring 3.5514 s. 6. Conclusions In this work, a chance-constrained MPC with disturbance estimation, for restricted three body problem rendezvous, is presented. Moreover, this robust controller is formulated to consider both chemical and electric thrusters, thus increasing the flexibility of the method. The chemical thrusters are modelled as impulsives and the electric ones are parameterized in terms of B-splines. The controller is limited to close rendezvous operations where the system dynamics can be linearized. The simulations have shown a great increase of mission success, sometimes at the expense of the cost, for the robust controller when compared to the non-robust one. The main drawback of the algorithm is the numerical integration of the state transition matrices since the dynamics is LTV. However these matrices can be computed by the ground control segment and loaded via uplink to the spacecraft. It is left as future work to evaluate the performance of this algorithm against other robust techniques such as worst-case methodologies, see [15], and tube-based MPC, see [16]. In conclusion, the presented chance-constrained model predictive controller describes an implementable, flexible and relatively fuel efficient algorithm for spacecraft close rendezvous operations in a complex dynamical system under the presence of disturbances. Acknowledgements The authors thank Jos´e Manuel Montilla and Jorge Gal´an-Vioque for discussions and help with NRHOs. The authors also gratefully acknowledge financial support from Universidad de Sevilla, through its V-PPI US, and from the Spanish Ministerio de Ciencia, Innovaci´on y Universidades under grant PGC2018-100680-B-C21. 24