# A Novel Concept for Guidance and Control of Spacecraft Orbital Maneuvers

Authors: Matteo Dentis; Elisa Capello; Giorgio Guglieri
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

The purpose of this paper is the design of guidance and control algorithms for orbital space maneuvers. A 6-dof orbital simulator, based on Clohessy-Wiltshire-Hill equations, is developed in C language, considering cold gas reaction thrusters and reaction wheels as actuation system. The computational limitations of on-board computers are also included. A combination of guidance and control algorithms for an orbital maneuver is proposed: (i) a suitably designed Zero-Effort-Miss/Zero-Effort-Velocity (ZEM/ZEV) algorithm is adopted for the guidance and (ii) a linear quadratic regulator (LQR) is used for the attitude control. The proposed approach is verified for different cases, including external environment disturbances and errors on the actuation system.


---

## Body

## 1. Introduction

Space systems often need to be controlled by actuators with limited output level, enforcing strict requirements in terms of relative position and attitude. In the past years, autonomous spacecraft rendezvous and docking (RVD) maneuvers have been extensively studied, in order to obtain controlled trajectories during which the Chaser (active) tries to dock a passive Target spacecraft. The goal of this mission is to safely and efficiently approach the Target vehicle to within a few centimeters (surface-to-surface), following a predefined path, which is generated by the guidance algorithm. In these maneuvers, the control system has to maintain the strict requirements, despite the limitations imposed in the actuation system and the external environment disturbances.In order to overcome inherent difficulties to achieve the desired objectives, it can be seen as a cooperation of several complex subsystems, with their own dynamics. In this paper, a combination of guidance and control algorithms is proposed, to obtain a potentially flight compliant algorithm in which some on-board limitations and errors/disturbances are taken into account. This research is based on the previous works of [1, 2], in which a GNC system of a ground test-bed for spacecraft rendezvous and docking experiments is designed. The novelty of this approach is the design and validation of a suitably designed guidance system based on Zero-Effort-Miss/Zero-Effort-Velocity (ZEM/ZEV) theory for a complete orbital maneuver, reducing the fuel consumption.Guidance algorithms should be divided into two classes: (i) a predictive guidance scheme and (ii) a feedback-based guidance law, which also uses on-off pulses (see Algorithm1). In the first class, the following guidance laws are investigated: (i) the Lambert guidance [3] and (ii) time-varying state transition matrix (STM) guidance [4]. Other guidance laws, based on the same theory, have been derived to follow defined path to the Target, for example [5], or to intercept an asteroid with terminal velocity direction constraints [6]. Considering the results obtained in [3], in which a simple 3 DOF simulator is considered for the validation of the proposed guidance, it is demonstrated that the computational effort is too high for the implementation on-board. For this reason, in this paper, the second class of guidance laws is considered, in which the well-known Proportional Navigation (PN) algorithm and the ZEM/ZEV theory are included. The PN law issues acceleration commands, perpendicular to the instantaneous Chaser-Target Line-Of-Sight (LOS), which are proportional to the LOS rate and closing velocity [7].Algorithm 1:ZEM/ZEV Guidance Algorithm.
[Waypoint Initialization] (i) setS0 position with current position (ii) setS1,  S2 and S3 positions as Table 1 (iii) start tracking Target position and setS4 position (iv) computeΔXS1S2 [Initialization and Coefficients evaluation] (i) sett0,man=t s, Δt=1 Hz ifS1-S2 computeΔVx,S1S2 ifS2-S3 computeΔVz,S2S3 (ii) computeA1,B1,C1,D1,A2,B2,A3,B3,C3 according to (19) with initial conditions according to the current Si-Si+1 phase End. [Initial Acceleration] ifS0-S1 not accelerate ifS1-S2 while (x˙act<x˙0,S1S2) accelerate to x˙0,S1S2 ifS2-S3 while (z˙act<z˙0,S2S3) accelerate to z˙0,S2S3 ifS3-S4 while (x˙act<0.15) accelerate to 0.15 m/s after acceleration switch to[ZEM/ZEV Guidance] [ZEM/ZEV Guidance] ifS0-S1 (i) null acceleration command:a=0,0,0T ifS1-S2 or S2-S3 (i) computeA1p,B1p,C1p,D1p,A2p,B2p,A3p,B3p,C3p by (19), setting the initial conditions as the actual position and velocity (ii) compute ZEM and ZEV according to (16)-(17)-(18)-(20)-(21) (iii) compute command acceleration according to (15) (iv) ifS2 or S3 are approaching goto [Stationkeeping] ifS3-S4 (i) computeA1p,B1p,C1p,D1p,A2p,B2p,A3p,B3p,C3p by (19), setting the initial conditions as the actual position and velocity (ii) computer~f and v~f according to (20)-(21) (iii) setrf=-500+x0,S3S4t-t0,man,S3S4,0,0T (iv) setvf=0.15,0,0T (v) compute ZEM and ZEV by (16) (vi) compute command acceleration by (15) (vii) ifxact<-4 stop controlling thrust onx-axis: a=[0,ZEMy,ZEMz] ifxact<-1 stop controlling thrust:a=[0,0,0] End. [Stationkeeping] ifS1-S2 (i) break to maintain position within a 10 × 10 m box aroundS2 and velocity lower than 10 cm/s (three axes) ifS2-S3 (i) break to maintain position within a 5 × 5 m box aroundS3 and velocity lower than 5 cm/s (three axes) when Stationkeeping is reached goto[ZEM/ZEV Guidance] End.Table 1
Waypoint LVLH position.
Waypoint Description Position[m] x , y , z L V L H S 0 Initial simulation waypoint [ - 12000,0 , 3000 ] S 1 Initial waypoint for Hohmann [ - 12000 - Δ X S 1 S 2 , 0,3000 ] S 2 Terminal waypoint for Hohmann/initial waypoint for Radial Boost [ - 3000,0 , 0 ] S 3 Terminal waypoint for radial boost/initial waypoint for straight line [ - 500,0 , 0 ] S 4 Terminal position [ x t a r g e t , y t a r g e t , z t a r g e t ]Bryson and Ho [6] also discussed optimal feedback control laws for a simple rendezvous problem and the relationship between optimal feedback control and PN guidance algorithm. Ebrahimi et al. [8] proposed the new concept of the Zero-Effort-Velocity (ZEV) error, analogous to the Zero-Effort-Miss (ZEM) distance. Both these feedback guidance laws can be implemented with on-off thrusters using a simple Schmitt trigger or other pulse-modulation devices as described in [9]. As clearly explained in [10], these algorithms are usually used for asteroid intercept and rendezvous missions, even if different mission requirements and spacecraft capabilities require continued research on terminal-phase guidance laws. Motivated by this issue, the ZEM/ZEV algorithm is selected and deeply analyzed. This algorithm evaluates directly the control commands in terms of accelerations starting from the positions of the Chaser, measurable at each time step. In detail, as said before, the idea is to evaluate the performance of the chosen algorithm, combining it with a linear quadratic regulator (LQR) for the attitude dynamics. Extensive theory about LQR can be easily found in the literature and application of LQR for attitude control of spacecraft [11–14]. LQR controllers are very effective in controlling linear systems, but, applying a model linearization, LQR can be also used in controlling nonlinear systems, as in the case of spacecraft attitude dynamics. Since LQR belongs to the class of optimal controller, attitude control could be optimized to reach and maintain its goal orientation with minimum time or energy used more easily than PID controller, even though this aspect is not covered by the scope of this paper.The advantage of this proposed approach lies in its simplicity, and, in our case, the information in terms of velocity and position of the target are taken into account, despite the classical implementation. Using a short term prediction of terminal condition instead of the classical terminal conditions of the complete maneuver guarantees a more precise trajectory control and earlier trajectory error recovery. Moreover, a reduction of the fuel consumption is guaranteed, in particular on curvilinear maneuvers like Hohmann transfer.The paper is organized as follows. In Section2, a detailed spacecraft model is analyzed, including actuator models, position, and attitude dynamics. The proposed GNC strategy is developed in Section 3. The simulation results are in Section 4. Conclusions are drawn in Section 5.


## 2. Spacecraft Model

The spacecraft model, analyzed in this work, is a 6-degree-of-freedom (dof) model in which actuator errors, nonlinearities, and external disturbances are taken into account. As discussed in the Introduction, a complete rendezvous maneuver is considered (see Figure1), in which a Chaser vehicle has to reach the Target spacecraft, initially moving along a different altitude orbit [15]. A detailed description of maneuver phases and waypoint positions will be provided in Section 4.Figure 1
Orbital maneuver.


### 2.1. Position Mathematical Model

For circular orbits, the translational dynamics is written in local-vertical local-horizontal (LVLH) frame (see Figure2). The x-axis Vbar is in the direction of the orbital velocity vector and the y-axis Hbar is in the opposite direction of the angular momentum vector of the orbit, while the z-axis Rbar is radial from the spacecraft center of mass to the Earth center of mass.Figure 2
LVLH frame definition.Clohessy-Wiltshire-Hill equations [15] are implemented to describe the relative motion of the two bodies in neighboring orbits:(1)x¨=Fxmc+2ω0z˙,y¨=Fymc-ω02y,z¨=Fzmc-2ω0x˙+3ω02z,where x=x,y,zT∈R3 is the position vector, mc∈R is the Chaser mass (known and varying with time), ω0=μ/rHF∈R is the angular velocity of the LVLH frame at a distance rHF from center of Earth, μ is the gravitational parameter of Earth (known and constant), and F=Fx,Fy,FzT∈R3 is the total force vector, which is the sum of the forces due to the thrusters and the forces due to the action of the external environment disturbances affecting the system; that is,(2)F=Fthr+ΔFex,where Fthr∈R3 is the time-varying vector of the forces by the thrusters (it will be described later) and ΔFex∈R3 is the time-varying vector of the forces related to external disturbances. This last contribution includes the environment aerodynamics force and the J2 effect. The Chaser mass varies depending on the decrease of fuel mass, according to the model in [16] and its function of the total force acting on the system (2), of the gravity gradient vector at sea level, and of the thruster specific impulse.The forces obtained from the thrusters and the external disturbances are transformed from body frame to LVLH frame. Hence, we have(3)F=Robϕ,θ,ψF∗,where F∗∈R3 represents the vector of force in body frame (FB frame, see Figure 3), F∈R3 represents the force vector in LVLH frame (see (2)), and Rob(ϕ,θ,ψ)∈R(3,3) is the transformation matrix between these two reference frames. This transformation matrix couples the position and the attitude dynamics.Figure 3
Body frame definition.


### 2.2. Euler Equations

For the evaluation of the spacecraft attitude angular velocities and angles, the body reference frame is considered, the origin of which is the center of mass of the spacecraft (see Figure3).In this reference frame, the body tensorI∈R3,3 of inertia is considered diagonal and it is updated taking into account the center of mass (CoM) change of position. The angular velocities are defined from the Euler equation as(4)ω˙B=I-1MB-ωB×IωB+IRWωRW,where ωB=pB,qB,rBT∈R3 is the vector of Chaser angular velocities, IRW is the RW moment of inertia (known and constant), and ωRW=ωx,RW,ωy,RW,ωz,RWT∈R3 is the reaction wheel (RW) angular velocity. MB∈R3 is the total moment acting on the Chaser and is given by the sum(5)MB=Mthr+ΔMex+MRW,where ΔMex∈R3 is the external disturbance moment, due to gravity gradient effect and the solar radiation pressure.


### 2.3. Kinematic Equations

The quaternions are used for the attitude evaluation:(6)q˙=12ΣωBq,where q=q0,q1,q2,q3T∈R4 is the vector of quaternions and ΣωB∈R4,4 is defined as (7)ΣωB=0-ωBTωB-Ω,where Ω is the skew-symmetric matrix,(8)Ω=0-ωzωyωz0-ωx-ωyωx0.The attitude is obtained by the time integration of (6). The attitude is propagated with respect to the Earth Centered Inertial (ECI) frame. The ECI frame is centered in the Earth’s CoM and the x-axis lies in the equatorial plane, pointing toward the mean of the vernal equinox. The z-axis is normal to the equatorial plane pointing north and the y-axis is defined to form a right-handed system. A quaternion in the form q=1,0,0,0T represents the body frame aligned with ECI one.


### 2.4. Actuation System

The actuation system for position control exploits thrusters, which can exert monodirectional actions along fixed directions and with fixed given magnitude. As will be detailed in the following, in each required control direction a pair of actuators is used and these thrusters coupled by direction are always switched on together. This precise choice of design for the actuation system guarantees a nominal zero moment due to the thrusters in the ideal case, when no thrusters errors occur. Each thruster is characterized by a fixed output and the individual thruster can provide either the maximum amount of thrust when switched on or no force when switched off.Further characteristics of each thruster are given by the time duration of the pulse widthτoni and the thruster zero time before turning on again τoffi. In fact, if the controller switches on the ith thruster at time t0, this actuator provides the maximum thrust Tmaxi for a time τoni and then the same thruster cannot be turned on again for a time τoffi, as shown in Figure 4 and according to the following:(9)Ti+=Tmaxi,ift∈t0,t0+τoni,0,ift∈t0+τoni,t0+τoni+τoffi,where τoni and τoffi are known and constant for each thruster. The constant τi=τoni+τoffi can be computed and its inverse constitutes the maximum allowed frequency at which the ith thruster can be switched on.Figure 4
Thrust provided by theith thruster switched on at time t0.To model the different sources of uncertainty, the total thruster force of theith thruster is described as(10)Fthri=βiTmagidthri,i=1,…,Nthr,where Nthr is the total number of thrusters. The vector β∈RNthr is a Boolean vector related to the thruster switching on/off. The magnitude Tmagi∈R of the force applied by the ith thruster is modeled considering the maximum thrust and a magnitude error applied by each thruster. The vector dthri∈R3 is the vector representing the shoot direction of the ith thruster: that is,(11)dthri=dthri0+Δdthri,where dthri0∈R3 is the vector of the nominal direction of thrust for each individual thruster and Δdthri∈R3 is the shooting direction error due to nonperfect assembly. The design of both shooting and magnitude thruster errors is based on the theory of [17] and it includes bias and random components.The resulting total force applied by thrusters is given by(12)Fthr=∑i=1NthrFthri.The total moment due to the thrusters is evaluated as(13)Mthr=∑i=1NthrFthri×li,where li∈R3 is the ith thruster position vector for i=1,…,Nthr. The vector li changes during the maneuver, due to the variation of the spacecraft CoM position, and it is defined in terms of x,y,z position with respect to the CoM in body reference frame. To model this variation, we write(14)li=li0+Δli,where li0∈R3 represents the vector of initial position of the thrusters and Δli∈R3 is due to CoM position variation.See Figure5 for the thruster distribution on the spacecraft and for the scheme of thruster positions.Figure 5
Thrust provided by theith thruster switched on at time t0.The attitude control is based on three axial reaction wheels with inertiaIRW, driven by electric motors powered by the spacecraft electrical power supply. They are managed and controlled by the on-board attitude control computer. A reaction wheel actuator produces a moment MRW, causing its angular momentum to increase. A saturation on the maximum-minimum torque assigned by the RWs is also included. The angular velocities of the reaction wheels ωRW and the moment MRW are calculated in body reference frame (described in Section 2.2).


## 2.1. Position Mathematical Model

For circular orbits, the translational dynamics is written in local-vertical local-horizontal (LVLH) frame (see Figure2). The x-axis Vbar is in the direction of the orbital velocity vector and the y-axis Hbar is in the opposite direction of the angular momentum vector of the orbit, while the z-axis Rbar is radial from the spacecraft center of mass to the Earth center of mass.Figure 2
LVLH frame definition.Clohessy-Wiltshire-Hill equations [15] are implemented to describe the relative motion of the two bodies in neighboring orbits:(1)x¨=Fxmc+2ω0z˙,y¨=Fymc-ω02y,z¨=Fzmc-2ω0x˙+3ω02z,where x=x,y,zT∈R3 is the position vector, mc∈R is the Chaser mass (known and varying with time), ω0=μ/rHF∈R is the angular velocity of the LVLH frame at a distance rHF from center of Earth, μ is the gravitational parameter of Earth (known and constant), and F=Fx,Fy,FzT∈R3 is the total force vector, which is the sum of the forces due to the thrusters and the forces due to the action of the external environment disturbances affecting the system; that is,(2)F=Fthr+ΔFex,where Fthr∈R3 is the time-varying vector of the forces by the thrusters (it will be described later) and ΔFex∈R3 is the time-varying vector of the forces related to external disturbances. This last contribution includes the environment aerodynamics force and the J2 effect. The Chaser mass varies depending on the decrease of fuel mass, according to the model in [16] and its function of the total force acting on the system (2), of the gravity gradient vector at sea level, and of the thruster specific impulse.The forces obtained from the thrusters and the external disturbances are transformed from body frame to LVLH frame. Hence, we have(3)F=Robϕ,θ,ψF∗,where F∗∈R3 represents the vector of force in body frame (FB frame, see Figure 3), F∈R3 represents the force vector in LVLH frame (see (2)), and Rob(ϕ,θ,ψ)∈R(3,3) is the transformation matrix between these two reference frames. This transformation matrix couples the position and the attitude dynamics.Figure 3
Body frame definition.


## 2.2. Euler Equations

For the evaluation of the spacecraft attitude angular velocities and angles, the body reference frame is considered, the origin of which is the center of mass of the spacecraft (see Figure3).In this reference frame, the body tensorI∈R3,3 of inertia is considered diagonal and it is updated taking into account the center of mass (CoM) change of position. The angular velocities are defined from the Euler equation as(4)ω˙B=I-1MB-ωB×IωB+IRWωRW,where ωB=pB,qB,rBT∈R3 is the vector of Chaser angular velocities, IRW is the RW moment of inertia (known and constant), and ωRW=ωx,RW,ωy,RW,ωz,RWT∈R3 is the reaction wheel (RW) angular velocity. MB∈R3 is the total moment acting on the Chaser and is given by the sum(5)MB=Mthr+ΔMex+MRW,where ΔMex∈R3 is the external disturbance moment, due to gravity gradient effect and the solar radiation pressure.


## 2.3. Kinematic Equations

The quaternions are used for the attitude evaluation:(6)q˙=12ΣωBq,where q=q0,q1,q2,q3T∈R4 is the vector of quaternions and ΣωB∈R4,4 is defined as (7)ΣωB=0-ωBTωB-Ω,where Ω is the skew-symmetric matrix,(8)Ω=0-ωzωyωz0-ωx-ωyωx0.The attitude is obtained by the time integration of (6). The attitude is propagated with respect to the Earth Centered Inertial (ECI) frame. The ECI frame is centered in the Earth’s CoM and the x-axis lies in the equatorial plane, pointing toward the mean of the vernal equinox. The z-axis is normal to the equatorial plane pointing north and the y-axis is defined to form a right-handed system. A quaternion in the form q=1,0,0,0T represents the body frame aligned with ECI one.


## 2.4. Actuation System

The actuation system for position control exploits thrusters, which can exert monodirectional actions along fixed directions and with fixed given magnitude. As will be detailed in the following, in each required control direction a pair of actuators is used and these thrusters coupled by direction are always switched on together. This precise choice of design for the actuation system guarantees a nominal zero moment due to the thrusters in the ideal case, when no thrusters errors occur. Each thruster is characterized by a fixed output and the individual thruster can provide either the maximum amount of thrust when switched on or no force when switched off.Further characteristics of each thruster are given by the time duration of the pulse widthτoni and the thruster zero time before turning on again τoffi. In fact, if the controller switches on the ith thruster at time t0, this actuator provides the maximum thrust Tmaxi for a time τoni and then the same thruster cannot be turned on again for a time τoffi, as shown in Figure 4 and according to the following:(9)Ti+=Tmaxi,ift∈t0,t0+τoni,0,ift∈t0+τoni,t0+τoni+τoffi,where τoni and τoffi are known and constant for each thruster. The constant τi=τoni+τoffi can be computed and its inverse constitutes the maximum allowed frequency at which the ith thruster can be switched on.Figure 4
Thrust provided by theith thruster switched on at time t0.To model the different sources of uncertainty, the total thruster force of theith thruster is described as(10)Fthri=βiTmagidthri,i=1,…,Nthr,where Nthr is the total number of thrusters. The vector β∈RNthr is a Boolean vector related to the thruster switching on/off. The magnitude Tmagi∈R of the force applied by the ith thruster is modeled considering the maximum thrust and a magnitude error applied by each thruster. The vector dthri∈R3 is the vector representing the shoot direction of the ith thruster: that is,(11)dthri=dthri0+Δdthri,where dthri0∈R3 is the vector of the nominal direction of thrust for each individual thruster and Δdthri∈R3 is the shooting direction error due to nonperfect assembly. The design of both shooting and magnitude thruster errors is based on the theory of [17] and it includes bias and random components.The resulting total force applied by thrusters is given by(12)Fthr=∑i=1NthrFthri.The total moment due to the thrusters is evaluated as(13)Mthr=∑i=1NthrFthri×li,where li∈R3 is the ith thruster position vector for i=1,…,Nthr. The vector li changes during the maneuver, due to the variation of the spacecraft CoM position, and it is defined in terms of x,y,z position with respect to the CoM in body reference frame. To model this variation, we write(14)li=li0+Δli,where li0∈R3 represents the vector of initial position of the thrusters and Δli∈R3 is due to CoM position variation.See Figure5 for the thruster distribution on the spacecraft and for the scheme of thruster positions.Figure 5
Thrust provided by theith thruster switched on at time t0.The attitude control is based on three axial reaction wheels with inertiaIRW, driven by electric motors powered by the spacecraft electrical power supply. They are managed and controlled by the on-board attitude control computer. A reaction wheel actuator produces a moment MRW, causing its angular momentum to increase. A saturation on the maximum-minimum torque assigned by the RWs is also included. The angular velocities of the reaction wheels ωRW and the moment MRW are calculated in body reference frame (described in Section 2.2).


## 3. Guidance, Navigation, and Control Algorithms

Proximity operations and docking require extremely delicate and precise translational and rotational maneuvering. During the orbital maneuver and, in particular, in the final approach of the proximity operations phase, the relative position, velocity, attitude, and angular rates between the two spacecrafts must be precisely controlled in order to obtain the required docking interface conditions. Thus, we propose an on-board Guidance, Navigation, and Control (GNC) concept to perform various rendezvous functions including translational and rotational control, targeting, and relative navigation, automatically and with extreme precision. As first approximation, we assume that an extended Kalman filter is implemented that filters all the sensor noise and all the measurement disturbances. The role of each element of the GNC system is(i)
to receive and elaborate signals from sensors (Navigation),(ii)
to compute the required acceleration or force needed to change, if required, the actual orbit and to follow a predefined path (Guidance),(iii)
to control the actuators enforcing the correct trajectory of the Chaser spacecraft pointing to the Target vehicle (Control).In the next section each function of the GNC software, developed in the present work, is deeply described.


### 3.1. Guidance Algorithm

The guidance function can be conceptualized as the trajectory generator. If a simple example is considered, the guidance law defines the trajectory that the spacecraft has to follow to reach the desired final position. The selected algorithm is the ZEM/ZEV, extensively discussed in [10, 18], and it corresponds to an optimal control only in the presence of uniform gravitational field, as in the case of Low Earth Orbit environment. The general formulation of the ZEM/ZEV control law is(15)a=6tgo2ZEM-2tgoZEV,where a∈R3 is the required control acceleration, tgo=tf-t is the time-to-go, and tf is the final maneuver time. The time tgo represents a parameter related to the expected time, the maneuver is concluded, and it is chosen a priori. The time t is the elapsing simulation time.Expanding theZEM and ZEV terms we obtain(16)ZEM=rf-r~f,ZEV=vf-v~f,where rf∈R3 and vf∈R3 are the desired final position and velocity and r~f∈R3 and v~f∈R3 are the predicted position and velocity at t=tgo.Choosing properly the final maneuver timetf is not a trivial task; indeed, a bad choice of tf could cause extremely high values of the command acceleration, since 1/tgo tends to infinity as the time t increases. To fix this behaviour it is possible to set a limiting value to tgo, once it becomes smaller than a specific limit value.The approach proposed in the present paper is as follows. Instead of computing the final position and velocity vectors required as terminal conditions of the entire maneuver, the ZEM/ZEV terms in (16) are implemented as positions and velocities that the spacecraft should have while it is following an ideal maneuver. As a consequence, rf and vf are the position and velocity vectors propagated foreword by a defined and constant Δt using the Clohessy-Wiltshire equations [15]. Considering rf=rfx,rfy,rfzT and vf=vfx,vfy,vfzT, the new required final position is(17)rfx=A1sin⁡ω0τpred-B1cos⁡ω0τpred+C1τpred+D1,rfy=A2cos⁡ω0τpred+B2sin⁡ω0τpred,rfz=A3cos⁡ω0τpred+B3sin⁡ω0τpred+C3,and, computing the time derivative, the new required final velocity is obtained(18)vfx=A1ω0cos⁡ω0τpred+B1ω0sin⁡ω0τpred+C1,vfy=-A2ω0sin⁡ω0τpred+B2ω0cos⁡ω0τpred,vfz=-A3ω0sin⁡ω0τpred+B3ω0cos⁡ω0τpred.The coefficients of (17)-(18) can be easily computed as(19)A1=4x˙0ω0-6z0,B1=2z˙0ω0,C1=6ω0z0-3x˙0,D1=x0+2z˙0ω0,A2=y0,B2=y˙0ω0,A3=2x˙0ω0-3z0,B3=z˙0ω0,C3=4z0-2x˙0ω0,where x0=x0,y0,z0T∈R3 is the initial position vector of the spacecraft and x˙0=x˙0,y˙0,z˙0T∈R3 is the initial linear velocity vector. Note that τpred=t-t0,man+Δt is the elapsed time from the beginning of the maneuver (t0,man) and Δt is the prediction horizon.The ideal trajectory and speed profiles are computed from (17) to (18). To compute the predicted position and velocity, respectively, r~f and v~f, Clohessy-Wiltshire equations are propagated over Δt, so we have(20)r~fx=A1psin⁡ω0Δt-B1pcos⁡ω0Δt+C1pΔt+D1p,r~fy=A2pcos⁡ω0Δt+B2psin⁡ω0Δt,r~fz=A3pcos⁡ω0Δt+B3psin⁡ω0Δt+C3p,(21)v~fx=A1pcos⁡ω0Δt+B1psin⁡ω0Δt+C1p,v~fy=-A2psin⁡ω0Δt+B2pcos⁡ω0Δt,v~fz=-A3psin⁡ω0Δt+B3pcos⁡ω0Δt,where all the coefficients are computed by (19) using the actual position and velocity as initial conditions.In order to generate the ideal trajectory that the Chaser shall follow for the phase of rendezvous and docking maneuver, Clohessy-Wiltshire equations have to be propagated forward with time from a starting timet0 to a finish time tf which varies according to the specific maneuver. In the ideal case, an impulsive maneuver of orbit altitude raising ends after half an orbital period, while a continuous thrust maneuver (with the same initial and terminal orbit altitude) ends after one orbital period. As depicted in Figure 1 the whole maneuver is composed of different phases (Hohmann, R-bar pulse (radial boost), etc.). This means that the initial conditions of the generated ideal maneuver have to be changed when a different phase is considered [15]. As a consequence, the GNC software has to compute the coefficients of (19) using different initial conditions for each phase before starting running the ZEM/ZEV algorithm and to generate the guidance law. So, t0,man is initialized at each waypoint Si after the station keeping conditions are reached. Finally, it is possible to generate the trajectory and the velocity profile that the Chaser should have Δt steps forward the current elapsed time t-t0,man from the beginning of the maneuver.To computeZEM and ZEV errors, it is required to estimate the predicted terminal conditions, which, in this proposed approach, are the position r~f and the velocity v~f computed forward in time by Δt. In this case, differently from the generation of the ideal trajectory, initial conditions used to compute coefficients of (19) are position and velocity vectors of the Chaser at the current simulation time t. Thus, it is possible to predict positions and velocities of the spacecraft after Δt in free dynamics (no control action and disturbances are applied). The coefficients of (20) are computed continuously during each phase of the simulation.


### 3.2. Attitude Control

Attitude control is performed by both activating reaction wheels and reaction control thrusters or Reaction Control System (RCS). Usually, fine and continuous attitude control is executed by RW, while the use of RCS, due to its impulsive behaviour, is required when the control torque is high (compared to the maximum torque generated RW) or to desaturate RW.Our idea, to satisfy pointing requirements and constraints, is to implement a LQR control system. For the definition of the LQR gains, a standard continuous time-invariant state-space formulation is evaluated:(22)x˙=Ax+Bu,where A∈R6,6 is the state matrix, B∈R6,3 is the control matrix, and the state is defined by x=q1,q2,q3,ω1,ω2,ω3T∈R6, with q1, q2, and q3 being the vectors of the first three components of the vector q (defined in Section 2.3). ω1,ω2,ω3 are, respectively, the x, y, and z components of the angular speed vector ωB of the spacecraft. The control vector u∈R3 is the control torque. Note that qrel and ωrel are relative quaternion and angular rates written in LVLH frame. Since the attitude dynamics model of the spacecraft is computed relative to the ECI frame, qrel and ωrel have to be derived in the inertial frame. The angular speed of the spacecraft relative to the ECI frame is(23)ωEBb=ωEOb+ωOBb,where ωEOb is the angular speed of the LVLH frame relative to ECI and ωOBb is the angular speed of the spacecraft with respect to LVLH frame. Rewriting the equation and expressing the rotation matrix from LVLH to body frame Rob we obtain(24)ωOBb=ωEBb-RobqωEOo,where, under the assumption of circular orbit, angular velocity of the LVLH frame is(25)ωEOo=0,ω0,0T,where ω0 is the angular velocity of the LVLH frame, as defined in Section 3.1. The rotation matrix Rob(q) can be defined using the quaternion qrel(26)Robq=q02+q12-q22-q332q1q2+q0q32q1q3-q0q22q1q2-q0q3q02-q12+q22-q332q2q3+q0q12q1q3+q0q22q2q3-q0q1q02-q12-q22+q33.In order to design the LQR controller, a linearization of the dynamic equations has to be applied. Bryson’s rule [19] is considered for the evaluation of the controller gains. The desired quaternion qrelid is then(27)qrelid=1,0,0,0Tand the resulting rotation matrix is(28)Robqrelid=12q3-2q2-2q312q12q2-2q11.The linearized angular velocities ωEOo and ωOBb are(29)ωEOo=-2q3ω0,-ω0,2q1ω0T,(30)ωOBb=ωEBb--2q3ω0-ω02q1ω0,with time derivative(31)ω˙OBb=ω˙EBb--2q˙3ω002q˙1ω0.The time derivative of the relative quaternion is(32)q˙relid=120-ω1-ω2ω3ω10ω3-ω2ω2-ω30ω1ω3ω2-ω10q0q1q2q3and, after linearization, the result is(33)q˙rel=12ω1,ω2,ω3,0T.The linearized dynamic model is obtained substituting (30) and (31) in (4). Then, it is possible to use the state-space representation (22) defining A∈R6,6 and B∈R6,3 matrix by(34)A14=A25=A36=12,A41=-2ω02Iyy-IzzIxx,A46=-Iyy-IzzIxx+1ω0,A63=-2ω02Ixx-IyyIzz,A64=-1+Ixx-IyyIzzω0,B41=1Ixx,B52=1Iyy,B63=1Izzin which gyroscopic torque due to RW is neglected due to limitation in software implementation. The other elements of the state and control matrices are zero. The design of the LQR controller consists in generating a control torque equal to(35)u=-Kx,where x is the state and K=R-1BTP is the static gain obtained by the solution P of the associated Algebraic Riccati Equation (ARE):(36)PA+ATP-PBR-1BTP+Q=0,where Q and R are weighting matrix of suitable dimensions.


## 3.1. Guidance Algorithm

The guidance function can be conceptualized as the trajectory generator. If a simple example is considered, the guidance law defines the trajectory that the spacecraft has to follow to reach the desired final position. The selected algorithm is the ZEM/ZEV, extensively discussed in [10, 18], and it corresponds to an optimal control only in the presence of uniform gravitational field, as in the case of Low Earth Orbit environment. The general formulation of the ZEM/ZEV control law is(15)a=6tgo2ZEM-2tgoZEV,where a∈R3 is the required control acceleration, tgo=tf-t is the time-to-go, and tf is the final maneuver time. The time tgo represents a parameter related to the expected time, the maneuver is concluded, and it is chosen a priori. The time t is the elapsing simulation time.Expanding theZEM and ZEV terms we obtain(16)ZEM=rf-r~f,ZEV=vf-v~f,where rf∈R3 and vf∈R3 are the desired final position and velocity and r~f∈R3 and v~f∈R3 are the predicted position and velocity at t=tgo.Choosing properly the final maneuver timetf is not a trivial task; indeed, a bad choice of tf could cause extremely high values of the command acceleration, since 1/tgo tends to infinity as the time t increases. To fix this behaviour it is possible to set a limiting value to tgo, once it becomes smaller than a specific limit value.The approach proposed in the present paper is as follows. Instead of computing the final position and velocity vectors required as terminal conditions of the entire maneuver, the ZEM/ZEV terms in (16) are implemented as positions and velocities that the spacecraft should have while it is following an ideal maneuver. As a consequence, rf and vf are the position and velocity vectors propagated foreword by a defined and constant Δt using the Clohessy-Wiltshire equations [15]. Considering rf=rfx,rfy,rfzT and vf=vfx,vfy,vfzT, the new required final position is(17)rfx=A1sin⁡ω0τpred-B1cos⁡ω0τpred+C1τpred+D1,rfy=A2cos⁡ω0τpred+B2sin⁡ω0τpred,rfz=A3cos⁡ω0τpred+B3sin⁡ω0τpred+C3,and, computing the time derivative, the new required final velocity is obtained(18)vfx=A1ω0cos⁡ω0τpred+B1ω0sin⁡ω0τpred+C1,vfy=-A2ω0sin⁡ω0τpred+B2ω0cos⁡ω0τpred,vfz=-A3ω0sin⁡ω0τpred+B3ω0cos⁡ω0τpred.The coefficients of (17)-(18) can be easily computed as(19)A1=4x˙0ω0-6z0,B1=2z˙0ω0,C1=6ω0z0-3x˙0,D1=x0+2z˙0ω0,A2=y0,B2=y˙0ω0,A3=2x˙0ω0-3z0,B3=z˙0ω0,C3=4z0-2x˙0ω0,where x0=x0,y0,z0T∈R3 is the initial position vector of the spacecraft and x˙0=x˙0,y˙0,z˙0T∈R3 is the initial linear velocity vector. Note that τpred=t-t0,man+Δt is the elapsed time from the beginning of the maneuver (t0,man) and Δt is the prediction horizon.The ideal trajectory and speed profiles are computed from (17) to (18). To compute the predicted position and velocity, respectively, r~f and v~f, Clohessy-Wiltshire equations are propagated over Δt, so we have(20)r~fx=A1psin⁡ω0Δt-B1pcos⁡ω0Δt+C1pΔt+D1p,r~fy=A2pcos⁡ω0Δt+B2psin⁡ω0Δt,r~fz=A3pcos⁡ω0Δt+B3psin⁡ω0Δt+C3p,(21)v~fx=A1pcos⁡ω0Δt+B1psin⁡ω0Δt+C1p,v~fy=-A2psin⁡ω0Δt+B2pcos⁡ω0Δt,v~fz=-A3psin⁡ω0Δt+B3pcos⁡ω0Δt,where all the coefficients are computed by (19) using the actual position and velocity as initial conditions.In order to generate the ideal trajectory that the Chaser shall follow for the phase of rendezvous and docking maneuver, Clohessy-Wiltshire equations have to be propagated forward with time from a starting timet0 to a finish time tf which varies according to the specific maneuver. In the ideal case, an impulsive maneuver of orbit altitude raising ends after half an orbital period, while a continuous thrust maneuver (with the same initial and terminal orbit altitude) ends after one orbital period. As depicted in Figure 1 the whole maneuver is composed of different phases (Hohmann, R-bar pulse (radial boost), etc.). This means that the initial conditions of the generated ideal maneuver have to be changed when a different phase is considered [15]. As a consequence, the GNC software has to compute the coefficients of (19) using different initial conditions for each phase before starting running the ZEM/ZEV algorithm and to generate the guidance law. So, t0,man is initialized at each waypoint Si after the station keeping conditions are reached. Finally, it is possible to generate the trajectory and the velocity profile that the Chaser should have Δt steps forward the current elapsed time t-t0,man from the beginning of the maneuver.To computeZEM and ZEV errors, it is required to estimate the predicted terminal conditions, which, in this proposed approach, are the position r~f and the velocity v~f computed forward in time by Δt. In this case, differently from the generation of the ideal trajectory, initial conditions used to compute coefficients of (19) are position and velocity vectors of the Chaser at the current simulation time t. Thus, it is possible to predict positions and velocities of the spacecraft after Δt in free dynamics (no control action and disturbances are applied). The coefficients of (20) are computed continuously during each phase of the simulation.


## 3.2. Attitude Control

Attitude control is performed by both activating reaction wheels and reaction control thrusters or Reaction Control System (RCS). Usually, fine and continuous attitude control is executed by RW, while the use of RCS, due to its impulsive behaviour, is required when the control torque is high (compared to the maximum torque generated RW) or to desaturate RW.Our idea, to satisfy pointing requirements and constraints, is to implement a LQR control system. For the definition of the LQR gains, a standard continuous time-invariant state-space formulation is evaluated:(22)x˙=Ax+Bu,where A∈R6,6 is the state matrix, B∈R6,3 is the control matrix, and the state is defined by x=q1,q2,q3,ω1,ω2,ω3T∈R6, with q1, q2, and q3 being the vectors of the first three components of the vector q (defined in Section 2.3). ω1,ω2,ω3 are, respectively, the x, y, and z components of the angular speed vector ωB of the spacecraft. The control vector u∈R3 is the control torque. Note that qrel and ωrel are relative quaternion and angular rates written in LVLH frame. Since the attitude dynamics model of the spacecraft is computed relative to the ECI frame, qrel and ωrel have to be derived in the inertial frame. The angular speed of the spacecraft relative to the ECI frame is(23)ωEBb=ωEOb+ωOBb,where ωEOb is the angular speed of the LVLH frame relative to ECI and ωOBb is the angular speed of the spacecraft with respect to LVLH frame. Rewriting the equation and expressing the rotation matrix from LVLH to body frame Rob we obtain(24)ωOBb=ωEBb-RobqωEOo,where, under the assumption of circular orbit, angular velocity of the LVLH frame is(25)ωEOo=0,ω0,0T,where ω0 is the angular velocity of the LVLH frame, as defined in Section 3.1. The rotation matrix Rob(q) can be defined using the quaternion qrel(26)Robq=q02+q12-q22-q332q1q2+q0q32q1q3-q0q22q1q2-q0q3q02-q12+q22-q332q2q3+q0q12q1q3+q0q22q2q3-q0q1q02-q12-q22+q33.In order to design the LQR controller, a linearization of the dynamic equations has to be applied. Bryson’s rule [19] is considered for the evaluation of the controller gains. The desired quaternion qrelid is then(27)qrelid=1,0,0,0Tand the resulting rotation matrix is(28)Robqrelid=12q3-2q2-2q312q12q2-2q11.The linearized angular velocities ωEOo and ωOBb are(29)ωEOo=-2q3ω0,-ω0,2q1ω0T,(30)ωOBb=ωEBb--2q3ω0-ω02q1ω0,with time derivative(31)ω˙OBb=ω˙EBb--2q˙3ω002q˙1ω0.The time derivative of the relative quaternion is(32)q˙relid=120-ω1-ω2ω3ω10ω3-ω2ω2-ω30ω1ω3ω2-ω10q0q1q2q3and, after linearization, the result is(33)q˙rel=12ω1,ω2,ω3,0T.The linearized dynamic model is obtained substituting (30) and (31) in (4). Then, it is possible to use the state-space representation (22) defining A∈R6,6 and B∈R6,3 matrix by(34)A14=A25=A36=12,A41=-2ω02Iyy-IzzIxx,A46=-Iyy-IzzIxx+1ω0,A63=-2ω02Ixx-IyyIzz,A64=-1+Ixx-IyyIzzω0,B41=1Ixx,B52=1Iyy,B63=1Izzin which gyroscopic torque due to RW is neglected due to limitation in software implementation. The other elements of the state and control matrices are zero. The design of the LQR controller consists in generating a control torque equal to(35)u=-Kx,where x is the state and K=R-1BTP is the static gain obtained by the solution P of the associated Algebraic Riccati Equation (ARE):(36)PA+ATP-PBR-1BTP+Q=0,where Q and R are weighting matrix of suitable dimensions.


## 4. Results

The simulation model is tested for a generic Chaser-Target combination involved in sequential flight phases, as in Figure1. The Chaser is considered in stable initial conditions along an orbit of height h = 500 km and the Target center of mass is located at the origin of the LVLH reference frame, and it is 12 km far from the Chaser. A cubic-shape Chaser (1.2 m) is considered with an initial mass of 600 kg. The layout adopted for thrusters is as in Figure 5, as discussed in Section 2.4. The parameters and the initial conditions are selected according to [1] and to the approaching maneuver examples in [15], including safety requirements.Position of waypoints relative to LVLH frame is summarized in Table1, while the complete rendezvous and docking maneuver is divided into four leading phases as summarized in Table 2.Table 2
Rendezvous and docking maneuver phase description.
Phase Name Initial condition Units S 0-S1 Free drift x 0 = - 12,  y0=0,  z0=3 [km] x ˙ 0 = 3 2 ω 0 z 0,  y˙0=z˙0=0 [m/s] S 1-S2 Hohmann x 0 = - 3 - Δ X S 1 S 2,  y0=0,  z0=3 [km] x ˙ 0 = 3 2 ω 0 z 0 + Δ V x , S 1 S 2,  y˙0=z˙0=0 [m/s] S 2-S3 Radial boost x 0 = - 3,  y0=0,  z0=0 [km] x ˙ 0 = y ˙ 0 = 0,  z˙0=ΔVz,S2S3 [m/s] S 3-S4 Straight line x 0 = - 500,  y0=0,  z0=0 [m] x ˙ 0 = 0.15,  y˙0=0,  z˙0=0 [m/s]Note that, with the exception of the starting point of the Hohmann maneuver (S1), the initial points of free drift, radial boost, and straight line maneuvers correspond to the position of the waypoints, respectively, S0, S2, and S3. The position of S1 results from computing ΔVx,S1S2 required to raise the orbit of the Chaser to the Target orbit (3 km in this case). When ΔVx,S1S2 is known, we evaluate the distance along Vbar that the Chaser will cover during the maneuver, and Vbar position of S1 is computed as ΔXS1S2 behind S2. ΔVx,S1S2 and ΔXS1S2 are calculated in the following form:(37)ΔVx,S1S2=ω04z0,ΔXS1S2=3πω0ΔVx,S1S2,using the value of z0 corresponding to S1-S2. Once the distance between two consecutive waypoints S2 and S3 is fixed, it is possible to compute the required ΔVx,S2S3 to complete the radial boost maneuver, so(38)ΔVz,S3S4=ω04S2x-S3x.The following cases are analyzed:(1)
Ideal case: in the ideal case, no external disturbances due to the LEO environment and errors induced by RCS are considered (see Table3). The maneuver described in Table 2 is executed completely, including reaching station keeping conditions before switching to the next phase, as described in the guidance algorithm. The frequency of the guidance control loop is 1 Hz and both the attitude control and the simulation frequency are set to 100 Hz. The results obtained are presented in Figure 6. As shown in Figure 6(d) ideal trajectories (green and magenta dashed lines) are well followed. The performance of the control law can be evaluated by analyzing Figure 6(g): except for initial acceleration boosts, the trajectory is controlled with few pulses along Vbar (Fx) and Rbar (Fz), with the presence of disturbance thrusts along Hbar (Fy) due to attitude misalignment. The two-sided chattering is due to the magnitude and sign of the control acceleration command of (15). The deadband is managed by the PWM function: the thrusters are switched on only if the acceleration required by the guidance algorithm is higher than the minimum feasible acceleration pulse.(2)
Real case: in the second case both atmospheric drag disturbances and RCS errors are activated. The same maneuver as in the ideal case is analyzed. In the next considered cases, the frequency of the guidance control loop is 1 Hz. We analyze two subcases:
(2.A)
We set the frequency of the control loop at 100 Hz (as in case (1)), to evaluate the performance of GNC algorithms as in the ideal case but introducing leading disturbances. See Table 4 for the results. In Figures 7(a) and 7(c) the effort required by the attitude control to counteract disturbance torque due to thrust errors is shown. Such misalignment causes also a greater effort for the guidance loop with respect to the ideal case, mainly highlighted by the presence of pulsed thrust along Hbar in Figure 7(g). The request of a larger control activation is translated into a larger consumption of fuel (see Table 4). Better results are obtained in terms of Δy˙f as the control along Hbar is active, instead of the ideal case, in which ZEM/ZEV errors along Hbar are not large enough (free dynamics motion) to require the activation of the control law.(2.B)
The control loop frequency is 50 Hz. The simulation frequency is 100 Hz and the guidance loop frequency is 1 Hz. As it is shown in Figures8(a), 8(c), and 8(g) the required control effort by both guidance and control is greater than the two previous cases, even though the Chaser is able to complete the entire maneuver even with a degradation of performance. This is due to the “low” update frequency of the controller in comparison with the high dynamics variation occurring during the maneuver. This is also translated in a high consumption of fuel (see Table 5) that summarizes the numerical results.Table 3
Ideal case terminal conditions.
Position error Δ x f = 3.0 ⋅ 1 0 - 4 [m] Δ y f = 2.5 ⋅ 1 0 - 3 Δ z f = - 7.4 ⋅ 1 0 - 3  Velocity error Δ x ˙ f = 0.15 [m/s] Δ y ˙ f = 2.5 ⋅ 1 0 - 7 Δ z ˙ f = 2.2 ⋅ 1 0 - 3  Propellant m p = 15.52 [kg]Table 4
Real case terminal conditions: attitude control at 100 Hz.
Position error Δ x f = 4.8 ⋅ 1 0 - 4 [m] Δ y f = 7.2 ⋅ 1 0 - 4 Δ z f = 7.4 ⋅ 1 0 - 3  Velocity error Δ x ˙ f = 0.15 [m/s] Δ y ˙ f = 4.3 ⋅ 1 0 - 5 Δ z ˙ f = 2.1 ⋅ 1 0 - 3  Propellant m p = 15.73 [kg]Table 5
Real case terminal conditions: attitude control at 50 Hz.
Position error Δ x f = 2.9 ⋅ 1 0 - 4 [m] Δ y f = 1.8 ⋅ 1 0 - 3 Δ z f = 7.2 ⋅ 1 0 - 3  Velocity error Δ x ˙ f = 0.15 [m/s] Δ y ˙ f = 1.6 ⋅ 1 0 - 4 Δ z ˙ f = 2.2 ⋅ 1 0 - 3  Propellant m p = 16.17 [kg]Figure 6
Ideal case results: controller frequency 100 Hz.
(a)
Relative quaternionqrel (b)
Relative angular velocityωrel (c)
RW torque command (d)
V b a r - R b a r trajectory (e)
x , y , z L V L H versus time (f)
x ˙ , y ˙ , z ˙ L V L H versus time  (g)
Thrust force in LVLH frame (h)
Final approachVbar-Rbar highlightFigure 7
Real case results: controller frequency 100 Hz.
(a)
Relative quaternionqrel (b)
Relative angular velocityωrel (c)
RW torque command (d)
V b a r - R b a r trajectory (e)
x , y , z L V L H versus time (f)
x ˙ , y ˙ , z ˙ L V L H versus time  (g)
Thrust force in LVLH frame (h)
Final approachVbar-Rbar highlightFigure 8
Real case results: controller frequency 50 Hz.
(a)
Relative quaternionqrel (b)
Relative angular velocityωrel (c)
RW torque command (d)
V b a r - R b a r trajectory (e)
x , y , z L V L H versus time (f)
x ˙ , y ˙ , z ˙ L V L H versus time  (g)
Thrust force in LVLH frame (h)
Final approachVbar-Rbar highlight


## 5. Conclusions

In the present work, a new approach for the development of a GNC algorithm has been proposed. The implementation of a modified version of ZEM/ZEV algorithm is novel as this innovative guidance method still finds limited application among the available references. The guidance design is combined with a classical LQR attitude control. This solution is novel and also very effective in driving a Chaser spacecraft along a complete rendezvous and docking maneuver, even in presence of disturbances due to external perturbations and to the imperfections of the actuation system. The impact of these disturbances and modeling errors was investigated simulating a reference maneuver with different controller frequency. The results confirm that the proposed solution is promising in terms of disturbance rejection, robustness, and performance, providing strict compliance of terminal conditions with nominal reference states, as a major achievement of ZEM/ZEV algorithm implementation.



---

*Source: 7695257-2016-11-29.xml*