Motion-planning and pose-tracking based rendezvous and docking with a tumbling target

Journal Pre-proofs Motion-Planning and Pose-Tracking Based Rendezvous and Docking with a Tumbling Target Bang-Zhao Zhou, Xiao-Feng Liu, Guo-Ping Cai PII: DOI: Reference:

S0273-1177(19)30804-X https://doi.org/10.1016/j.asr.2019.11.013 JASR 14537

To appear in:

Advances in Space Research

Received Date: Revised Date: Accepted Date:

2 March 2019 25 October 2019 11 November 2019

Please cite this article as: Zhou, B-Z., Liu, X-F., Cai, G-P., Motion-Planning and Pose-Tracking Based Rendezvous and Docking with a Tumbling Target, Advances in Space Research (2019), doi: https://doi.org/ 10.1016/j.asr.2019.11.013

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.

Motion-Planning and Pose-Tracking Based Rendezvous and Docking with a Tumbling Target Bang-Zhao Zhou1,

Xiao-Feng Liu2,

Guo-Ping Cai3*

Department of Engineering Mechanics, State Key Laboratory of Ocean Engineering, Shanghai Jiaotong University, Shanghai 200240, China

Abstract: Rendezvous and docking (RVD) with a tumbling target is challenging. In this paper, a novel control scheme based on motion planning and pose (position and attitude) tracking is proposed to solve the pose control of a chaser docking with a tumbling target in the phase of close range rendezvous. Firstly, the current desired motion of the chaser is planned according to the motion of the target. In planning the desired motion, the “approach path constraint” is considered to avoid collisions between the chaser and the target, and the “field-of-view constraint” is considered to make sure the vision sensors on the chaser to obtain tight relative pose knowledge of the target with respect to the chaser. Then, the difference between the chaser’s motion and the desired motion is gradually reduced by a pose tracking controller. This controller is based on the non-singular terminal sliding mode (NTSM) method to make the tracking error converge to zero in finite time. Since the chaser nearly moves along the desired motion and the motion is reasonable, (1) it could safely arrive at the docking port of the target with a suitable relative attitude, (2) it will be always suitably oriented to observe the target well, and (3) the magnitude of the needed control inputs are less than that in existing literatures. The numerical results demonstrate the above three advantages of the proposed method. Keywords: Rendezvous and docking; Tumbling target; Motion planning; Pose tracking

1

Introduction The spacecraft rendezvous and docking (RVD) with a space object is the premise of on-orbit

servicing (OOS) missions such as berthing, refueling, repairing, upgrading, transporting, rescuing,

1

PhD student. Email address: [email protected].

2

Post-doctor. Email address: [email protected].

3*

Professor. Email address: [email protected] (Corresponding Author). 1 / 41

and orbital debris removal, here the space object refers to a malfunctioning satellite or, more generally, a space debris, both of which will be called “target” in this paper. All of the current and past RVD missions focus only on a cooperative target which is supposed to move smoothly in its orbit without rapid attitude changes (Flores-Abad et al., 2014; Ma et al., 2007; Dalla and Pathak, 2015, 2019a, 2019b), this steady attitude results in an easier docking. In practice, when a target satellite malfunctions, its attitude control system may be disabled. Furthermore, the angular momentum stored in the gyros and wheels will start migrating to its body as soon as failure occurs. As a result, the defunct satellite tumbles (rotates about time-varying axis). This phenomenon had been confirmed by the ground observations (Kawamoto et al., 2003). So far, RVD missions with a tumbling target have not been done yet (Flores-Abad et al., 2014; Ma et al., 2007). RVD missions with a tumbling target is still an open research area facing many technical challenges and has been a hot area of research for several years. In the final phase (where the distance between these two crafts is less than 1000 m) of a RVD mission, a chaser spacecraft should synchronously track both the time-varying desired relative positions and attitude trajectories accurately with respect to a target spacecraft (Dong et al., 2018). Therefore, pose (position and attitude) tracking is critical in the final phase. Kjellberg and Lightsey (2013) developed a general purpose guidance, navigation, and control algorithm for satellites with three degree-of-freedom rotation maneuverability. Pong and Miller (2015) developed a guidance and a control law to maneuver the boresight axis of a spacecraft relative to a reference vector, using knowledge of the reference vector and the angular rates of the spacecraft. Tekinalp and Tekinalp (2016) proposed a spacecraft attitude control which uses the quaternion parametrization. Cruz et al. (2012) applied retrospective cost adaptive control (RCAC) to spacecraft attitude control. Agarwal et al. (2012) studied the inertia-free spacecraft attitude control with control moment gyroscope actuation. Sanyal and Chaturvedi (2008) addressed the practical problem of tracking the attitude and angular velocity of a spacecraft in the presence of gravity and disturbance moments. Yoon and Tsiotras (2005) developed an adaptive control algorithm for the spacecraft attitude tracking problem for the case when the spin axis directions of the flywheel actuators are uncertain. Cruz et al. (2011) considered torque-saturated control laws for spacecraft attitude control. Bandyopadhyay et al. (2016) presented an attitude control strategy and a new nonlinear tracking controller for a spacecraft carrying a large object. Gao et al. (2016) addressed the high-precision coordinated control problem of spacecraft autonomous rendezvous and docking. Lee et al. (2015) 2 / 41

presented a tracking control scheme for spacecraft formation flying. The pose tracking control problem for a chaser approaching to a space target is researched in (Lee et al., 2014; Lee and Vukovich, 2016; Sun et al., 2017; Sun and Zheng, 2018). And these proposed control schemes are robust in the presence of external disturbances and uncertain system parameters. Sun and Zheng (2017) proposed an adaptive relative pose control strategy for a chaser in the final phase on a tumbling target. In their strategy, relative position vector between these two crafts is required to direct towards the docking port of the target and the attitude of them must be synchronized. However, in the final phase, the chaser also needs to follow some important motion constraints. On the one hand, to avoid collisions, the chaser should approach the docking port of the target from a safe direction. On the other hand, the attitude of the chaser also needs to be properly controlled to allow the vision sensors on the chaser to obtain tight relative pose knowledge, especially when the target is non-cooperative. These two constraints are called “approach path constraint” and “field-of-view constraint” respectively in (Dong et al., 2018), see Fig. 1. As shown in Fig. 1 (a), rC/T represents the position vector of the chaser with respect to the target, pT represents the docking axis of the target, and the approach path constraint requires that the angle β between rC/T and pT be less than βm. As shown in Fig. 1 (b), rT / C   rC / T , cb represents the optical axis of the camera on the chaser, and the field-of-view constraint requires that the angle α between rT/C and cb be less than αm. target

target rC/T βm

rT/C

β pT

α chaser

(a) Approach path constraint Fig. 1.

chaser αm cb

(b) Field-of-view constraint

Approach path constraint and field-of-view constraint (Dong et al., 2018).

In the previously mentioned researches (Gao et al., 2016; Lee et al., 2014; Lee et al., 2015; Lee and Vukovich, 2016; Sun et al., 2017; Sun and Zheng, 2017, 2018), the position of the chaser is controlled to track the docking port of the target, and the attitude of the chaser is controlled to 3 / 41

track that of the target. However, this control objective is a little insufficient. Firstly, the field-ofview constraint does not mean tracking the attitude of the target, on the contrary, tracking the attitude of the target may violate the field-of-view constraint in some situation. For example, see Fig. 2, in a situation where the optical axis of chaser’s camera is along its docking axis, if the chaser is asked to track the attitude of the target, the target can be observed when the chaser is at the position P1, while the observation will be fail when the chaser is at P2. Secondly, if the chaser is controlled in this way, since the tracking process is not constrained, the motion of the chaser may violate the above constraints during the process of reducing the tracking error. target

chaser

approach corridor

P1

Fig. 2.

P2

chaser

Two situations of chaser’s relative position with respect to target.

In order to avoid collision between the chaser and the target, some researchers studied the position control in the phase of close range rendezvous. In (Li et al., 2018), the relative distance between the two crafts is constrained to be larger than the radius of the danger zone during close proximity phase. However, if there is any solar panel attached to the target or chaser, this constraint will prevent the docking between them since the danger zone is too big. In (Li et al., 2017), the approach path constraint is simplified as a situation where βm = π/2, where βm can be find in Fig. 1. However, this constraint is not suitable in some situation where the shape of the target is complicated, such as the International Space Station (ISS), see Fig. 3. Li and Zhu (2018) proposed a position controller using model predictive control. The controlled chaser strictly complies with the approaching path constraint. However, they only studied the case where the target’s attitude is stable or spinning. Breger et al. (2008) presents a method for online generation of safe, fueloptimized rendezvous trajectories that guarantee collision avoidance for a large class of anomalous system behaviors. In (Dong et al., 2017; Richards et al., 2002; Weiss et al., 2015), the controlled chaser not only strictly complies with the approaching path constraint, but also has obstacle4 / 41

avoidance ability. And the found trajectories in (Richards et al., 2002) are also fuel-optimal.

Fig. 3.

ISS size information and axes definition (Dong et al., 2017).

In the phase of close range rendezvous, a chaser should accurately track not only a suitable position trajectory to comply with the approaching path constraint, but also a desired attitude trajectory to satisfy the field-of-view constraint. The control problem of the chaser under both these two constraints is a challenging 6-DOF constrained control problem. Most of the existing literature can only solve a part of this complex problem (see, for example, in (Breger et al., 2008; Dong et al., 2017; Li and Zhu, 2018; Li et al., 2018; Li et al., 2017; Richards et al., 2002; Weiss et al., 2015), only the approaching path constraint is considered, while in (Ventura et al., 2016), only the field-of-view constraint is considered). To the best of the authors’ knowledge, there are only two studies (Dong et al., 2018; Lee and Mesbahi, 2014) in which this 6-DOF constrained control problem is solved. Both these two controllers in (Dong et al., 2018; Lee and Mesbahi, 2014) ensure the arrival of the chaser at the docking port of the target with a suitable relative attitude, and no constraint violation happens during the entire process. However, there is some imperfection in their studies: (1) the approach in (Lee and Mesbahi, 2014) can only deal with a specialized class of stabilization problems, but cannot be readily applied in more general tracking problems (Dong et al., 2018); (2) in both these two studies, the approaching path constraint is treated in the way where it demands that the chaser should be always in the infinite cone, see Fig. 1 (a), while we think this treatment is a little too rigorous. In fact, if the relative distance between the two crafts is large enough, the chaser will not collide with the target, so in this case it is not necessary to constrain the chaser in the infinite cone. Besides, this treatment brings unnecessary trouble. On the one hand, it requires the initial position of the chaser is in the infinite cone. On the other hand, it may need unnecessarily bigger control force to drive the chaser if this treatment is applied, 5 / 41

especially when the target is rotating. If the target rotates, this body-fixed infinite cone also rotates about the target. If the chaser is constrained in this rotating cone, the speed and acceleration of the chaser will be rather big when the distance between these two crafts is big. As a result, the control force will be rather large (this phenomenon can be found in the simulation results of (Dong et al., 2018). Therefore, this treatment is not perfect. Since the above 6-DOF constrained control problem is really challenging, and it is difficult to suitably treat the approaching path constraint, we adopt another method to solve the pose control of a chaser docking with a tumbling target in the phase of close range rendezvous. This method is based on motion planning and pose tracking. Firstly, a method for planning the current desired pose motion of the chaser according to the motion of the target is proposed. If the chaser moves along the desired motion, its motion will be suitable in the phase of close range rendezvous. Namely it will arrive at the docking port of the target with a suitable relative attitude, and no constraint violation happens during the entire process. Then, a pose tracking controller is designed to gradually reduce the difference between the chaser’s motion and the desired motion. It should be noted that the desired motion also has the following characteristics: (1) its initial position and velocity is close to that of the chaser, (2) the linear and angular velocity of the desired motion is continuous, and (3) the linear and angular acceleration is within a small range. Therefore, it is easy to track the desired motion, namely the tracking error will be small. If the tracking error is small enough, the real motion of the chaser will be also suitable in the final phase. Besides, unlike in the existing literature (Dong et al., 2018; Lee and Mesbahi, 2014), in this paper, even if the initial direction of the chaser with respect to the target is arbitrary, our method always works. It should be noted that, Filipe and Tsiotras (2014) have applied this idea in their simulations. However, their method is not perfect enough: (1) they only considered a case where the target’s attitude is stable; (2) their method for motion planning did not achieve a high generality; (3) it is not easy to track their planned motion since the planned velocity is not continuous. In this paper, we avoid these shortcomings and improve their method. The contents of this paper are organized as follows. In Section 2, the relative translational and rotational dynamics of the chaser are presented. Then in Section 3, a pose tracking controller is designed. In Section 4, a method is proposed to plan the desired motion of the chaser. After that numerical simulation results are shown in Section 5. Finally, the conclusions are stated in Section 6. 6 / 41

2

Relative Motion Kinematics and Dynamics In this paper, the target and chaser are assumed to be rigid bodies, and some disturbances

including the Earth oblateness (J2) effect, the gravity gradient moment, and unknown bounded disturbances are considered in our dynamic model. The mass center of the target moves freely in an elliptical orbit, and the rotational motion is tumbling. The chaser need to track a desired motion to gradually dock with the target. The desired motion is obtained by specifying the relative pose motion with respect to the target, see Section 4. The frame ΣD is called “desired frame” in this paper, and its pose motion is the desired motion. The relative motion kinematics and dynamics of the chaser with respect to ΣD is introduced in this section. This is necessary in designing the pose tracking controller in Section 3.

2.1 Translational Motion Differential Equation In order to describe the chaser’s motion, the following frames are introduced. As shown in Fig. 4, the frame ΣN is an Earth-centered inertial frame. Its x-axis points to the vernal equinox, its z-axis is along the Earth’s axis of rotation. The frame ΣH built at the target’s mass center T is the local-level local-horizon Euler-Hill frame. Its x-axis is along the position vector of the target with respect to ΣN, and the z-axis is aligned with the orbit angular momentum of the target. The bodyfixed frame ΣC of the chaser is built at its mass center. The position vector of ΣD with respect to ΣH is ρD, ρC is that of ΣC, and ρr represents the position vector of ΣC with respect to ΣD.

Fig. 4.

Some frames in this paper.

7 / 41

If the target’s mass center moves freely in an elliptical orbit, the relative translational motion4 (Lee et al., 2014) of the chaser with respect to ΣH is

aCH  f  ( ρCH , vCH , t ) 

FCH (t )  a JH2 ( ρCH , t )  aH (t ) mC

(1)

H T where ρC  [ xC , yC , zC ] is the coordinate array of ρC in ΣH (in this paper, the superscript “T”

represents transposition, a vector with a superscript “H” denotes the coordinate array of this vector in ΣH, note that vT does not represent the transposition “vT” but the coordinate matrix of v in ΣT), H 2 d H d H T  v  ρC   [ xC , y C , zC ] , aC  2 ρCH   [ xC , yC , zC ]T , FCH is the applied control force, dt dt

H C

H

mC is the mass of the chaser, aJH2 is the relative effect of the Earth oblateness due to J2 (see Appendix A), aH is the sum of all other unknown disturbance accelerations, and f ( ρCH , vCH , t ) is given by

  ( rT  xC )    2 2 y C f  yC f  xC f  [( r  x ) 2  y 2  z 2 ]3 / 2  r 2  T C C C T     yC f  ( ρCH , vCH , t )    2 x C f  xC f  yC f 2   2 2 2 3/ 2 [( rT  xC )  yC  zC ]     zC    [( r  x ) 2  y 2  z 2 ]3 / 2 T C C C  

(2)

where rT is the distance between the target and the Earth, f is the true anomaly of the target, and μ = 3.986032×1014 m3/s2 is the gravitational constant. Assuming that the motion of the target is known, namely the variables (rT, f, f , and f ) are given. As is shown in Fig. 4, ρr = ρC – ρD, so H 2 d2 H d H a  a  a , here a  2 ρr  and aD  2 ρDH  . Therefore, the translational motion dt dt H r

H C

H D

H r

H

differential equation of the chaser with respect to ΣD can be written as

arH  G ( ρrH , vrH , t ) 

4

FCH (t )  aH (t ) mC

(3)

Strictly speaking, when higher harmonics in gravitational potential is considered, the orbit is no longer elliptical and no flat.

Consequently, Eqs. (1) and (2) is not valid for long time. Anyway for the considered case in this paper these two equations are appropriate. 8 / 41

where

G ( ρrH , vrH , t )  f ( ρCH , vCH , t )  aJH2 (t )  aDH (t )  H  ρC  ρrH  ρDH (t ), vCH  vrH  vDH (t )

(4)

Here ( ρ DH (t ) , v DH (t ) , a DH (t ) ) represents the known desired translational motion.

2.2 Rotational Motion Differential Equation The rotation vector and angular velocity of ΣC with respect to ΣD are denoted by Θr = θrpr and ωr respectively. Here the unit vector pr is the axis of rotation and  r  Θr  [0,  ] the angle of rotation. The kinematic equation (Lee and Vukovich, 2016) is

Θ r  A(Θr )ωrC

(5)

1~ 2

~

~

2 where A(Θr )  I 3  Θr  f ( r )(Θr ) , here Θr is the skew-symmetric matrix of Θr and

 1 1  cos  r   f  ( r )   r2 2 r sin  r  1 / 12

( r  0,  )

(6)

( r  0)

Therefore, lim f ( r )  1 / 12 , and the sole singular point of Eq. (5) is θr = π. On the one hand,  r 0

the maximum value of θr is π. On the other hand, from Eq. (B4), we know that θr is monotonically decreasing if our controller is applied. Therefore, only when the initial value of θr is π, the singular point could be encountered; namely Eq. (5) is almost completely applicable in this paper. Denote ωC as the angular velocity of ΣC with respect to ΣN. The Euler’s dynamic equation of the chaser is

~ C J ω C  M C (t )  M C (t )  d C (t ) J C ω CC  ω C C C C g 

(7)

C C where JC is the chaser moment of inertia, M C is the applied control torque, d is the unknown

C disturbance torque, and Mg 

3 C r  JC rCC is the gravity gradient moment. Meanwhile, ωr is 5 C rC

C C CD D the angular velocity of ΣC with respect to ΣD, so ωr  ωC  A ωD . Here ωDD is the angular

velocity of ΣD with respect to ΣN, and ACD is the direction cosine matrix from ΣD to ΣC. Therefore, 9 / 41

the rotational motion differential equation of the chaser with respect to ΣD can be written as:

Θ r  A(Θr )ωrC  C ω r  G (Θr , ωrC , t )  J C1 [ M CC (t )  dC (t )]

(8)

d C ωr  , and dt ~ C J ωC  M C (t )]  ACD ω D  ω ~ C ACD ω D G (Θr , ωrC , t )  J C1[ ω C C C g D r D

 rC  where ω

C

(9)

D C CD D C  DD (t ) ) represent the known desired Here ωC  A ωD  ωr , and the parameters ( ωD (t ) , ω

rotational motion.

3

Controller Design Sliding mode controller is well known for its robustness to system parameter variations and

external disturbances (Man and Xing Huo, 1997). In applying a sliding mode controller, firstly the controller drives and constrains the system states to lie within a neighborhood of the sliding manifold; then since when in the sliding manifold, the closed-loop response becomes totally insensitive to both internal parameter uncertainties and external disturbances, all the system states will gradually converge to zero. The non-singular terminal sliding mode (NTSM) controller (Feng et al., 2002) overcomes the singularity problem of the terminal sliding mode controller. Moreover, both the time taken to reach the manifold from any initial state and the time taken to reach the equilibrium point in the sliding mode can be guaranteed to be finite time. In this section, a NTSM controller is designed, so that the pose of the chaser (ΣC) will converge to that of ΣD in finite time. For the design of NTSM controller, the unknown disturbances are assumed to be bounded by some known constants ( d C,i ,max , a H, j ,max ), namely  d C,i  d C,i ,max , i  1, 2, 3   H H  a  , j  a  , j ,max , j  1, 2, 3

(10)

The NTSM model is defined as

s(t )  xd  Bxvp / q

(11)

T H ,T T where xd  [Θr , ρr ] , B = diag(b1, b2, b3, b4, b5, b6) is a positive definite matrix, and p / q is

10 / 41

a power of corresponding values, namely xvp / q  [rp, x/ q rp, y/ q rp, z/ q v rp, x/ q v rp, y/ q v rp, z/ q ]T . Here [ωr,x, ωr,y, ωr,z]T = ωrC , [vr,x, vr,y, vr,z]T = vrH , and in order to avoid singular point when applying the NTSM controller (Feng et al., 2002), the positive odd integers q, p are chosen such that 1 < p/q < 2. A standard NTSM model requires that x d  xv , so in the NTSM manifold (s(t) = 0), all the system states converge to zero in finite time (Feng et al., 2002). It should be noted that in the defined NTSM model Eq. (11), x d  x v . However, in the NTSM manifold of this paper, the system states xd and xv also converge to zero in finite time, and this is proved in Appendix B. The time derivative of s(t) is

Θ r  ω rC  C H s(t )   H   BG(ωr , vr ) H  vr  ar  where G(ωr , vr )  C

H

(12)

p diag(r(,pxq) / q , ... , vr(,pzq) / q ) . Substituting Eqs. (3) and (8) into Eq. (12), then q

G  J C1[ MCC (t )  dC (t )]  A(Θr )ωrC  s(t )      BG H H H G  FC (t )/mC  a (t ) vr 

(13)

Let [ BG1 , BG 2 ]  BG , where BG 1 , BG 2  R 63 . Then Eq. (13) can be written as C  d C ( t )   A(Θ ) ω C   B   M C (t ) 1 s( t )   H r r   BG 1G   BG 2 G    BG 1 J C1 G 2   H   BG 1 J C BG 2  H  m C   FC ( t )    a  ( t )   v r 





(14)

Let us consider the following control law for the systems described by Eqs. (3) and (8)  M CC (t )  H   Γ c , eq  Γ c , n  FC (t ) 

(15)

where 1 C    BG 2    A(Θr )ωr  1   B G  B G  Γ c ,eq    BG1 J C     G 1 G 2   mC   vrH      1   BG 2  1 diag( D(t )  η) sgn( s)  Γ c ,n    BG1 J C mC     d C   D(t )  [ D (t ), ..., D (t )]T  abs B J 1 B    ,max  1 6 G1 C G2 H   a ,max  





11 / 41

(16)

Here η = [η1, … , η6]T (ηi > 0, i = 1, 2, …, 6), sgn(s) = [sgn(s1), . . . , sgn(s6)]T is a column of signum functions, both dC,max  [d C,1,max , d C, 2,max , d C,3,max ]T and aH,max  [aH,1,max , aH,2,max , aH,3,max ]T

 A11  have been defined in Eq. (10), and the function abs(A) is defined as abs( A)     Am1 Theorem 1

  

  . Amn  A1n 

For systems described by Eqs. (3) and (8), with the NTSM model defined in Eq. (11),

if the control law is designed as in Eq. (15), then the NTSM manifold will be reached in finite time. Furthermore, the states xd and xv will converge to zero in finite time. Proof

Consider the following Lyapunov function candidate:

1 V  sT s 2

(17)

Taking time derivative of the Lyapunov function V results in V  s T s

(18)

Substituting Eq. (14) into the above equation, one can obtain   A(Θr )ωrC   M CC (t )  T 1 BG 2  1  V s  H  H   BG1 J C BG 2   BG1G  BG 2G    BG1 J C   v r mC   FC (t )    





dC (t )    H   a (t ) 

(19)

Substituting the control law (Eqs. (15) and (16)) results in

 V  s T  BG1 J C1 BG 2  



Let BG1 J

1 C

BG 2





 dC (t )     ) sgn( ) ) ( diag( s η D t  H    a (t ) 

(20)

dC (t )   H  = [d1, … , d6]T, then  a (t )

6

6

i 1

i 1

V    ( Di  i ) si  d i si     ( Di  d i )  i  si

(21)

From Eqs. (10) and (16), we know Di(t) > di(t), so 6

V    i si  0

(22)

i 1

Therefore, the state trajectory of the system will finally reach the NTSM manifold s(t) = 06×1. Moreover, the time for this reaching is no more than 2 6V0 /  j , here V0 is the initial value of V, 12 / 41

and j may be 1, 2, …, or 6. It is proved in Appendix C. As mentioned before, after the NTSM manifold is reached (in finite time), in this NTSM manifold, the system states xd and xv converge to zero also in finite time. Therefore, the tracking error will converge to zero in finite time if our pose tracking controller is applied.

4

Motion Planning In this section, the current desired pose motion of the chaser is planned in real time according

to the motion of the target. If the chaser moves along this desired motion, its motion will be suitable in the phase of close range rendezvous. Namely it will safely arrive at the docking port of the target with a suitable relative attitude, and it will be always suitably oriented to let its camera observe the target with an excellent line-of-sight angle α. Here α refers to the angle between the direction of the target with respect to the chaser and the optical axis of camera, see Fig. 1 (b). Once the desired motion is planned, then the chaser can track it using the NTSM controller proposed in the last section. If the tracking error is small enough, the chaser’s real motion will be also suitable in the final phase. In order to guarantee a small tracking error, the desired motion needs to be tracked well. Therefore, in planning the desired motion, the chaser’s initial conditions should be considered, the planned velocity should be continuous and the planned acceleration should be within a small range. The desired motion is planned below.

4.1 Planning of Translational Motion The desired translational motion of the chaser is planned in this subsection, and the rotational motion in the next subsection. Without loss of generality, assuming the docking axis of the target is on the y-axis of its body-fixed frame ΣT (shown in Fig. 5). The chaser’s attitude needs to be suitable with respect to ΣT when docking. In this paper, the chaser’s body-fixed frame ΣC is built at its mass center, and its orientation is chosen in the following way: ΣC is just parallel with ΣT when docking. Therefore, the pose of ΣC with respect to ΣT when docking can be represented in Fig. 5. In this figure, the sphere is the danger zone of the target, F is the docking port of the target, SF is the approach corridor, and S is the entrance of the approach corridor. The length of TF is rf, and rs for TS. The chaser will not collide with the target if its mass center is not in the danger zone.

13 / 41

danger zone (sphere) zT zC T xT

Fig. 5.

F

yC

S

xC

yT

The pose of ΣC when docking.

The planned translational motion in this paper can be divided into the following five phases: (1) ΣD moves from the initial position to a position D1, here the length of TD1 is krsrs, (krs > 1, namely ΣD is always out of the danger zone during this phase), (2) ΣD moves from D1 to S, (3) ΣD keeps staying at S for a timespan of τ3, (4) ΣD moves to F along the approach corridor SF, (5) ΣD keeps staying at F for a timespan of τ5. The details of each Phase are shown as follows. Phase 1: When planning the translational motion of Phases 1 and 2, a best value vopt of the current desired velocity vD is firstly determined, then the desired acceleration aD is calculated according to the difference between vopt and vD, so that vD can converge to vopt quickly. And aD is given by

 v am v aD   k a v v v 

k k

a

a

v  am v  am

 (23)



where Δv = vopt–vD, am denotes the maximum acceleration in motion planning, and ka is a parameter in motion planning. After aD is calculated, the planned velocity vD and position ρD can be obtained by integration. The way for determining vopt is shown below. During Phase 1, the distance between these two crafts is more than rs, so the chaser will not collide with the target. The vopt points to the mass center T, and its magnitude is vm, namely H opt

v

  vm

ρDH

(24)

D

Then the desired acceleration aDH can be calculated through Eq. (23). After that, the translational motion parameters ( ρDH , vDH , and aDH ) of D with respect to ΣH can be obtained through integration.

14 / 41

Phase 2: During this phase, the chaser should arrive at the entrance S with a zero-velocity, and it should not enter the danger zone before the arrival. Let ρDT  [ xDT , yDT , zDT ]T denote the position vector of D in ΣT. As shown in Fig. 6, if y DT  rs , the direction of vopt is along the tangent DP, or the direction is along DS. The process of calculating the direction of vopt can be found in Appendix D. zT

T xT

danger zone

yT (yTp)

F

S βT P

xTp D

Fig. 6.

The magnitude of vopt is vopt





 vm kv DS  vm   kv DS  vm kv DS



Direction of vopt in ΣT.



(25)

T  vopt evT , where kv is a parameter in motion planning. So far, the vopt in ΣT has been obtained as vopt

here e vT

denotes the direction of vopt. Then the desired acceleration aDT in ΣT, namely

d2 T T ( ρD ) , can be calculated by substituting vopt into Eq. (23). After that all the parameters ρ DT dt2 T

and v DT can be obtained through integration. Phase 3: During this phase, the chaser stays still at S. Therefore, ρDT = [0, rs, 0]T, v DT = aDT = 0. Phase 4: During this phase, the chaser moves to the docking port F along the approach corridor SF, and both the initial and terminal values of v DT are 0. Therefore, the magnitude of the desired velocity can be planned as

15 / 41

am (0     sf )  ( sf     f   sf ) v D  vm ,sf a (   ) (       ) f sf f  m f

(26)

where τ = t – ts4, ts4 is the initial moment of phase 4, τsf = vm,sf/am, and τf = (rs - rf)/vm,sf + τsf. Then all the parameters ( ρDT , v DT , and aDT ) can be obtained through Eq. (26). Phase 5: this phase is analogous to Phase 3. It should be noted that, during Phases 2–5, the obtained motion parameters ( ρDT , v DT , and aDT ) are all with respect to ΣT. However, the needed motion parameters ( ρDH , vDH , and aDH ) in the translational motion differential equations (3) and (4) are all with respect to ΣH. These needed parameters can be calculated in the following way. As mentioned before, if the reference frame is chosen as ΣH, the frame ΣT rotates about its origin T (which is also the origin of ΣH), and the angular velocity and acceleration are known as

ωTH/ H  AHT ωTT / N  ωHH / N  H ω T / H  AHT ω TT / N  ω HH / N

(27)

respectively. And the motion parameters of D with respect to ΣT are known. Therefore, the motion parameters of D with respect to ΣH can be calculated as  ρDH  A HT ρDT  H H H v D  vrel  ve a H  A HT a T  ω H  ρ H  2ω H  v H  ω H  v H D T /H D T /H rel T /H e  D

(28)

H where v rel  A HT v DT and v eH  ω TH / H  ρ DH .

4.2 Planning of Rotational Motion In this section, the desired rotational motion of the chaser is planned. If the chaser moves along the planned motion, the optical axis of chaser’s camera will always point to the target exactly, and the attitude of the chaser will be suitable when docking. It should be noted that, the optical axis of the camera on the chaser should be along the minus y-axis of ΣC. Otherwise the chaser probably could not observe the target when docking, because ΣC should be parallel with ΣT at that time, see Fig. 5. The planned rotational motion can be divided 16 / 41

into the following three stages. Stage 1: the stage 1 of the rotational motion corresponding to the Phases 1 and 2 of the planned translational motion. During this stage, the planned motion should guarantee that the optical axis will always point to the target exactly. Stage 2: at the end of Stage 1, ΣD just arrive at S, then it will stay here for some time. Now the y-axis of ΣD is along TS (for suitable observation), so it is parallel with the y-axis of ΣT, see Fig. 7. Therefore, during this stage, ΣD is planned to rotate about the y-axis of ΣT till ΣD is parallel with ΣT. This rotation is single axis rotation, and its angular velocity can be planned using Eq. (26). Here the planned maximum angular velocity is denoted by ωp,m, and the planned angular acceleration is denoted by αp. They correspond to vm,sf and am in Eq. (26) respectively. single axis rotation

danger zone zT

zD

xD T

F S

yD

xT

Fig. 7.

yT

The pose of ΣD when it just arrives at S.

Stage 3: ΣD is planned to be always parallel with ΣT. During Stage 3, D stays on the line FS, see Fig. 7. Therefore, this planned motion can guarantee that the optical axis will always point to the target exactly, and the attitude of the chaser will be suitable when docking. The planned rotational motion of Stages 2 and 3 has been explained clearly above. And the details of Stage 1 are shown below. The target’s mass center T is at the origin of ΣH, see Fig. 8. When planning the rotational motion, the y-axis of ΣD should be always along with ρDH exactly, so that the optical axis will always point to the target exactly. As is shown in Fig. 8, the rotation angle dθ of ΣD with respect to ΣH during a timespan of dt is γ dθ T

Fig. 8.

vDdt

D ρD

ΣH

Relationship between rotational and translational motions. 17 / 41

vDdt sin  D  vDdt cos

d 

(29)

where vD is the norm of v DH , γ is the angle between ρDH and vDH . Therefore,  equals to   v D sin  /  D

(30)

It is the norm of the angular velocity ωD/H of ΣD with respect to ΣH. As is shown in Fig. 8, the direction of ωD/H is along ρDH  vDH . Therefore, ωD/H can be written as ω

H D/H

 y z  y z  1   2 zx  zx   D   x y  x y 

(31)

where [x, y, z]T = ρDH , [ x, y , z]T  vDH , and  D2 = x2 + y2 + z2. Then the angular acceleration of ΣD with respect to ΣH can be calculated as ω

H D/H

 yz  yz  d H 1   2( xx  yy  zz ) H ωD / H   ωD / H  2 zx  zx     D2 dt D   xy  xy  H

(32)

So far, the rotational motion parameters of ΣD with respect to ΣH can be obtained as ΘD/H, ω DH / H , and ω DH / H by integration. Here ΘD/H is the rotation vector of ΣD with respect to ΣH. In the rotational motion differential equations (8) and (9), the rotational motion parameters ( ωDD , ω DD ) of ΣD with respect to ΣN are needed. These parameters can be calculated by

ωDD  ωDD/ H  ωHD / N  D DH H D D DH H ω D  A ω D / H  ωH / N  ωD / H  A ω H / N

(33)

where ADH is the direction cosine matrix from ΣH to ΣD,

ω DD / H  A DH ω DH / H , and

ω HD / N  A DH ω HH / N .

5

Numerical Simulations In this section, numerical simulation examples are presented to illustrate the effectiveness of

the proposed method. Assuming that the target satellite is running on a Molniya orbit, with the orbital elements given in Table 1.

18 / 41

Table 1 Orbital elements of the target satellite (Dong et al., 2018) Values

Orbital elements Semimajor axis, km Eccentricity Inclination, deg

26553.937 0.729677 63.4

Argument of perigee, deg

-90

Right ascension of the ascending node, deg

0

The mass of the chaser is mC = 15 kg, and its inertia matrix in body-fixed frame ΣC is JC = diag(3.0514, 2.6628, 2.1879) kg  m 2 . The target’s inertia matrix in body-fixed frame ΣT is diag(3.85, 4.36, 4.90) kg  m 2 . At the beginning, the target is at the perigee. The initial conditions for translational motion of the chaser are given by

 ρCH (t0 )  [200, 150,  250]T (m)  H vC (t0 )  [0.1, 0.2, 0.3]T (m/s)

(34)

As for the initial conditions for attitude motion, at the beginning, the rotation angle between the chaser’s attitude and the best attitude is a uniformly distributed pseudorandom numbers on [0 10] deg, and the elements of ω CC ( t 0 ) are all uniformly distributed pseudorandom numbers on [-0.25, 0.25] deg/s. Here the best attitude is a suitable attitude of the chaser so that the optical axis of the camera on the chaser points to the target exactly. The parameters of the NTSM controller are as follows: p = 17, q = 15, B = I6, and η = [1 1 1 1 1 1]T. The sampling period is chosen as 0.1 s. The common phenomenon of chattering when applying a sliding mode controller can be mitigated by replacing the signum function in the control scheme Eq. (16) with the continuous saturation function defined as:  si /  , si   sat ( si ,  )   sgn( si ), si  

(35)

where ε is set as 10 in our simulations. The control inputs are required to be bounded by FCH, i  mC ×0.5 m/s2 and M CC,i  1 N m . The unknown disturbances are given as H 4 4 T 2  a (t )  10  10 [sin(t1 )  sin(t2 ),  sin(t1 )  sin(t2 ), cos(t1 )  cos(t2 )] ( m/s )  C 5 5 T  d (t )  10  10 [sin(t1 )  sin(t2 ),  sin(t1 )  sin(t2 ), cos(t1 )  cos(t2 )] ( Nm)

19 / 41

(36)

where t1 = πt/125, and t2 = πt/200. Therefore, the upper bound of the disturbances can be denoted by aH,max = 3×10-5[1 1 1]T and dC, max = 3×10-5[1 1 1]T. The parameters for motion planning are as follows: krs = 2, rs = 10 m, rf = 3 m, τ3 = 300 s, τ5 = 200 s, am = 0.05 m/s2, vm = 0.5 m/s, vm,sf =0.2 m/s, ka = 0.1 m0.5s-1.5, kv = 0.1 m0.5s-1, αp = 0.01 rad/s2, and ωp,m = 1 deg/s. In theory, the initial conditions of desired motion should be the same as that of the chaser, so that the desired motion can be tracked exactly just at the beginning. However, the real motion of the chaser cannot be known exactly, so when planning the chaser’s desired translational motion, a slight difference between the initial conditions of desired motion and real motion is considered, and the differences are given by

 ρDH (t0 )  ρCH (t0 )  0.2[ r1 , r2 , r3 ]T  H vD (t0 )  vCH (t0 )  0.02[ r4 , r5 , r6 ]T

(37)

where r1, r2, …, r5, and r6 are standard normally distributed pseudorandom numbers.

5.1 An Uncontrolled Tumbling Target An uncontrolled tumbling target is common in OOS missions. This kind of target is considered in this section. The target’s mass center moves freely in its orbit. As for its attitude motion, the gravity gradient moment and small disturbance dT,T  dC are considered. And the initial conditions for its attitude motion are given by

ΘT/N (t0 )  [0,0,0]T  T ωT/N (t0 )  [0.025,0.025,0.05]T (rad / s)

(38)

T where ΘT/N (t0 ) and ωT / N (t0 ) are the initial value of the rotation vector and angular velocity of

ΣT with respect to ΣN respectively. In order to demonstrate the necessity of the motion planning in our method, two contrast simulations are shown below.

5.1.1

Simulation I

In this simulation, the position of the chaser is controlled to track that of the docking port F on the target, and the attitude of the chaser is controlled to track that of the target. Namely motion 20 / 41

planning is not used in this simulation. The position and attitude tracking errors are shown in Fig. 9 and Fig. 10 respectively. As shown in these two figures, the position and attitude tracking errors both converge to about zero finally. This means the chaser could finally arrive at the docking port with a suitable relative attitude. 150

400

(deg)

200

r

r

(m)

300

100 0

0

100

200

300

100 50 0

400

0

100

200

Time (s)

Fig. 9.

300

400

Time (s)

Position tracking error.

Fig. 10.

Attitude tracking error.

However, this method cannot make sure that no constraint violation happens during the entire process. Firstly, note that the initial distance between these two crafts is about 354 m, and the chaser arrived the docking port at about t = 320 s; namely in this example the chaser’s velocity in ΣH is too big (more than 1 m/s). Therefore, it’s necessary to constrain the chaser’s velocity in this example. Secondly, simulation results of the line-of-sight angle α(t) is given in Fig. 11. As shown in this figure, α could be up to 150 deg during the entire process. However, a reasonable maximum value of α is about 30 deg (Dong et al., 2018). Therefore, in this example, the field-of-view constraint is violated, namely during the entire process the chaser sometimes cannot observe the target.

(deg)

200 150 100 50 0

0

100

200

300

400

Time (s)

Fig. 11.

Line-of-sight angle α(t).

Finally, in order to clearly show the translational motion of the chaser with respect to the target, after the chaser is close enough to the target, the trajectory in the target’s body-fixed frame ΣT is shown in Fig. 12. In this figure, the sphere represents the danger zone, the point on the spherical surface is the entrance S of the approach corridor, and the approach corridor is along the straight line through S and the center T of the danger zone. As shown in Fig. 12, the chaser did not approach 21 / 41

the docking port along the approach corridor. Therefore, the chaser could collide with the target.

Fig. 12.

Translational motion of the chaser in ΣT.

In summary, if the chaser is controlled in this way in the phase of close range rendezvous, its velocity could be too big, and it may violate both the approach path constraint and the field-ofview constraint. It is not surprising that such a result appears. In this way, the chaser is controlled to track the pose of docking port F, while there is no constraint on how to reduce the tracking error. And the initial value of the position tracking error (namely the initial distance between the chaser and F) is usually up to hundreds of meters. Therefore, the chaser may violate some constraints during the process of reducing the tracking error. It should be noted that this simulation also demonstrates that only pose tracking could not well solve the pose control of a chaser docking with a tumbling target in the phase of close range rendezvous.

5.1.2

Simulation II

In the phase of close range rendezvous, the proposed method for the control of the chaser is making the chaser track the motion planned in Section 4. The translational motion of the chaser in ΣH is shown in Fig. 13. As illustrated in Section 4.1, the planned translational motion is divided into five phases. Here the five phases are separated by four vertical dotted lines in this figure. As shown in this figure: (1) during Phase 1, the chaser almost moves uniformly and in a straight line (with respect to ΣH, which is almost an inertial frame); (2) during Phase 2, the chaser’s velocity keeps changing to gradually track the motion of S (rotation about T); (3) during Phase 3, the chaser tracks the motion of S, so its position keeps changing; (4) during Phase 4, the chaser moves to the 22 / 41

docking port F along the approach corridor, so the amplitude of its position decreases gradually; (5) during Phase 5, the chaser tracks the motion of F (rotation about T), so its position keeps changing. 1200 1000

x y z

20

(m)

600

H C

800

400

0 -20

200

600

800

1000

1200

0 -200 0

500

1000

1500

Time (s)

Fig. 13.

Translational motion of the chaser in ΣH.

It should be noted that in this simulation, the desired motion on the timespan of [0, 1322] s of the chaser is calculated, and its required time is 45.9 s. Namely calculating 1s of the desired motion requires only 0.035 s. Note that it is calculated on Matlab 2014a, the CPU model is Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz, and the size of the RAM is 8.0 GB. If our algorithm is programmed by C, the required time will be much less. Therefore, our motion planning algorithm can calculate the desired motion of the chaser nearly in real time. After the chaser is close enough to the target, its translational trajectory in the target’s bodyfixed frame is shown in Fig. 14. As shown in this figure, the chaser is always outside the danger zone before it arrives at S (the color of the trajectory in the danger zone is lighter); then it approaches the docking port F along the approach corridor. It should be noted that the trajectory before the chaser arrives at S seems to be complicated. In fact, as shown in Fig. 13, the trajectory in ΣH before the chaser arrives at S is nearly a straight line. Here the trajectory is shown in the target’s body-fixed frame ΣT. Since the target tumbles, the trajectory looks complicated.

23 / 41

Fig. 14.

Translational motion of the chaser in ΣT.

The control inputs are shown in Fig. 15 and Fig. 16. The control force is directly relative to the translational motion in Fig. 13. As shown in Fig. 15, the control force is also divided into five phases as the translational motion in Fig. 13. At the beginning of phase 1, 2, and 4, and the end of phase 4, the chaser’s acceleration changes suddenly, so the control force at these moments changes suddenly too. As mentioned in Section 4.2, in this paper, the planned rotational motion is divided into three stages. These three stages are separated by two vertical dotted lines in Fig. 16. At the beginning of stage 1 and 2, and the end of stage 2, the chaser’s angular acceleration changes suddenly, so the control torque at these moments changes suddenly too. 0.15

0.08

MC (Nm)

F C/m C (m/s 2)

0.1

0.06 0.04

0.1

0.05

0.02 0

0

500

1000

0

1500

0

500

Time (s)

Fig. 15.

1000

1500

Time (s)

Control force.

Fig. 16.

Control torque.

The position and attitude tracking errors are shown in Fig. 17 and Fig. 18 respectively. These two figures are divided into several parts as in Fig. 15 and Fig. 16 respectively. As shown in these figures, the tracking errors decrease rapidly within dozens of seconds; after that, the position

24 / 41

tracking error is less than 2 mm, and 0.04 deg for the attitude tracking error. Also note that in these figures, there are tiny pinnacles at the moment when control inputs change suddenly. Since the tracking error is very small, the desired motion is tracked well, namely the arrival of the chaser at the docking port with a suitable relative attitude can be guaranteed. 10 8

0.4

6

0.3

4

0.2

2

0.1

0

10

0

500

0

0

0.08

8

(deg)

0.5

-3

r

r

(m)

0.6

1000

500

0.04 4

0.02

2

1500

1000

0.06

6

1500

0

0

0

500

Time (s)

Fig. 17.

500

0

1000

1500

1000

1500

Time (s)

Position tracking error.

Fig. 18.

Attitude tracking error.

It should be noted that, both the initial position and velocity are considered when planning the desired translational motion of the chaser, so the position tracking error is not large (no more than 0.6 m) at the beginning. Moreover, the NTSM controller is used to track the desired motion, so the tracking error rapidly converges to about zero, specifically the position tracking error rapidly converges within 2 mm during the early stage of Phase 1 (of the planned translational motion). During Phase 1, the distance between the chaser and the target is more than 2rs (20 m). Since the chaser will not collide with the target if the distance is more than rs, it is acceptable that even the position tracking error reaches 1 m during Phase 1. Therefore, combining the results in Fig. 14, it can be concluded that the collision between the chaser and the target is avoided. Simulation results of the line-of-sight angle α(t) is shown in Fig. 19. As shown in this figure, the α is controlled be less than 1 deg rapidly. As explained in the first paragraph of Section 4, α refers to the angle between the direction of the target with respect to the chaser and the optical axis of camera on the chaser, see Fig. 1 (b). Therefore, once α is within a range (e.g., α is less than 30 deg), the chaser can observe the target and measure the target’s pose. In this paper, the α is controlled be less than 1 deg rapidly, which means the target will be nearly at the center of visual field of the camera on the chaser. Therefore, the chaser can observe the target well and our method is excellent. Note that the α(t) does not converge to zero. The reason is as follows. As shown in

25 / 41

Fig. 20, the target is at T, the chaser’s desired pose is at C1 when t  t1 , the chaser’s desired pose is at C2 when t  t2 , and the arrows denote the desired direction of the chaser’s optical axis of camera. If the attitude error is zero and the position error is δ, (1) when t  t1 , the chaser may be at C1’, then its line-of-sight angle is α1; (2) when t  t2 , the chaser may be at C2’, then its line-ofsight angle is α2. Therefore, as shown in Fig. 19 and Fig. 20, the α(t) may become larger when the chaser is approaching the target; nevertheless, the α is controlled be less than 1 deg rapidly, which ensures good observation of the chaser to the target satellite.

(deg)

10

5

0

C1’

α1

C2’

δ 0

500

1000

C1

1500

Time (s)

Fig. 19.

Line-of-sight angle α(t).

Fig. 20.

C2

α2 T

The desired and real pose of the chaser.

In summary, on the one hand, the desired motion has a series of advantages: collision avoidance, guarantee of excellent line-of-sight angle, and reasonable velocity; on the other hand, the desired motion could be tracked well. Therefore, if our method is applied in the phase of close range rendezvous, it is not necessary to constrain the chaser’s motion when tracking the desired motion. In contrast, if the chaser directly tracks the position of the docking port, its motion also needs to be constrained; otherwise both the approach path constraint and the field-of-view constraint could be violated, and its velocity could be too large. Therefore, in this paper, (1) the challenging 6-DOF constrained control problem as in (Dong et al., 2018) is avoided, namely the design of the pose tracking controller is greatly simplified; (2) the safe arrival of the chaser at the docking port with a suitable relative attitude can be guaranteed; (3) the chaser will be suitably oriented to observe the target well; and (4) the chaser’s velocity will be reasonable (very close to the desired one).

5.2 A Controlled Tumbling Target In this section, an attitude-controlled tumbling target is considered. This case is more general. Both the translational and rotational motion of the target is the same as the target in (Dong et al.,

26 / 41

2018). The translational motion has been given at the beginning of Section 5. The rotational motion is given by T T T  ΘT/N (t0 )  [0,0,0] , ωT/N (t0 )  [0,0,0]  T T 2  ω T/N (t )  0.01  [ 0.45 sin(0.1t ), cos(0.2t ), 0.5 cos(0.1t )] (rad/s )

(39)

The control inputs are shown in Fig. 21. As shown in this figure, the maximum value of FC/mC and MC are about 0.18 m/s2 and 0.08 N m respectively. While in (Dong et al., 2018), these values are up to about 3.5 m/s2 and 1.5 N m respectively. It should be noted that, both the inertial parameter of the chaser and the pose motion of the target are the same as those in (Dong et al., 2018). Therefore, the amplitude of the control input is greatly slashed in this paper. The reason for this great improvement is shown below. In the method of (Dong et al., 2018), the chaser is always constrained within a body-fixed cone of the target, see Fig. 1 (a), even when the chaser is far (e.g., hundreds of meter) away from the target. If the target tumbles, all the points in the target’s bodyfixed frame ΣT rotates about the mass center of the target, so the farther the point is, the bigger the acceleration of this point is. Therefore, if the chaser is far from the tumbling target, it will need rather large control input to keep the chaser within the body-fixed cone. While in this paper, because the chaser will not collide with the target if the distance between them is more than rs (10 m), we think there is no need to keep the chaser in the cone during Phase 1 (where the distance is more than 2rs). Therefore, in this paper, on the one hand, the chaser almost moves uniformly and in a straight line during Phase 1, so the control input is rather small during this phase; on the other hand, after that, the chaser gradually keep still in ΣT, since the distance between the two crafts is 0.2

0.08

0.15

0.06

MC (Nm)

F C/m C (m/s 2)

less than 2rs, the control input will not be too big now.

0.1

0.02

0.05 0

0.04

0

500

1000

0

1500

0

500

1000

1500

Time (s)

Time (s)

(a) Control force

(b) Control torque Fig. 21.

Control inputs.

Although the needed control inputs are greatly slashed in this paper, the control objective can 27 / 41

still be fulfilled. After the chaser is close enough to the target, its translational trajectory in the target’s body-fixed frame is shown in Fig. 22. As shown in this figure, the chaser is always outside the danger zone before it arrives at S; then it approaches the docking port F along the approach corridor. Similarly, here the trajectory looks complicated because the reference frame ΣT is tumbling.

Fig. 22.

Translational motion of the chaser in ΣT.

The position and attitude tracking errors are shown in Fig. 23 and Fig. 24 respectively. As shown in these two figures, the tracking errors decrease rapidly within dozens of seconds, after that, the position tracking error is less than 4 mm, and 0.04 deg for the attitude tracking error. Namely the tracking error is very small, so the safe arrival of the chaser at the docking port with a suitable relative attitude can be guaranteed. 0.3

12 0.01

0.25

(deg)

0.005

0.15

0.08

8

0.06

6

0.04

4

0.02

2

0

r

r

(m)

0.2

10

0.1 0

0.05 0

0

0

500

1000

500

1500

1000

1500

0

0

0

500

Time (s)

Fig. 23.

500

1000

1500

1000

Time (s)

Position tracking error.

Fig. 24.

28 / 41

Attitude tracking error.

1500

Simulation results of the line-of-sight angle α(t) is shown in Fig. 25. As shown in this figure, the α is controlled be less than 1 deg rapidly, which means the target will be nearly at the center of visual field of the camera on the chaser. Therefore, even if the target is controlled and tumbling, the chaser can observe the target well.

(deg)

10

5

0

0

500

1000

1500

Time (s)

Fig. 25.

Line-of-sight angle α(t).

In summary, even in a case where the attitude of the target is controlled and tumbling, if our method is applied in the phase of close range rendezvous, (1) the safe arrival of the chaser at the docking port with a suitable relative attitude can be also guaranteed, (2) the chaser will be suitably oriented to observe the target well, and (3) the magnitude of the needed control input is relatively smaller.

5.3 Monte Carlo Simulations In order to verify the effectiveness of the method, the following simulation was designed. Let H the elements of ρC (t0 ) (the chaser’s initial position vector) be uniformly distributed random

H numbers on [15, 200] m, and the elements of vC (t0 ) (the chaser’s initial velocity vector) be

uniformly distributed random numbers on [-0.3, 0.3] m/s. In the case of uncontrolled tumbling target (in Section 5.1.1) and attitude-controlled tumbling target (in Section 5.1.2), 1000 simulations were performed respectively, and the simulation results were statistically analyzed. In this way, almost all initial direction of the chaser with respect to the target are considered. Therefore, the validity of our method can be checked by this simulation. Before the chaser arrives at the point S (the entrance of the approach corridor, see Fig. 5), the maximum depth of the chaser’s entering the spherical danger zone is denoted by dm. The frequency histogram of dm is shown in Fig. 26. As shown in this figure, dm is no more than 7 mm. It should be noted that the radius of the danger zone is 10 m, which is much larger than 7 mm. Therefore, if

29 / 41

our method is applied in the phase of close range rendezvous, it can be considered that the chaser almost will not enter the danger zone before arriving at the point S, and so the chaser will not 500

400

400

300

Frequency

Frequency

collide with the target during this phase.

300 200 100 0

0

2

4

200 100 0

6

2

0

4

6

d m (mm)

d m (mm)

(a) Uncontrolled tumbling target (b) Attitude-controlled tumbling target Fig. 26. Statistical results of dm.

After the chaser arrives at the point S, it approaches the target along the approach corridor SF. If the pose tracking errors during this phase are too large, the chaser may collide with the target, and the subsequent docking operation cannot be executed. The maximum pose tracking errors during this phase are denoted by r,m and r,m respectively, and the frequency histograms of them are shown in Fig. 27 and Fig. 28. As shown in these figures, r,m is mainly distributed on [1, 2] mm if the target is uncontrolled, while r,m is mainly distributed on [2, 3] mm if the target is attitude-controlled, and r,m is mainly distributed on [0, 0.4] deg no matter the target is controlled or not. Such tracking error is acceptable in close range rendezvous. Because generally speaking, if the position tracking error is less than 20 mm and the attitude tracking error is less than 1 deg, the subsequent docking operation can be executed smoothly. Therefore, the precision 200

250

150

200

Frequency

Frequency

of our method is acceptable.

100 50 0

150 100 50

1

1.2

1.4

1.6

0

1.8

(mm) r,m

0

0.2

0.4

0.6

(deg) r,m

(a) Position tracking error (b) Attitude tracking error Fig. 27. Statistical results of pose tracking errors (uncontrolled tumbling target).

30 / 41

250

300

200

Frequency

Frequency

400

200 100 0

150 100 50

2

2.5 r,m

0

3

0

0.2

(mm)

0.4 r,m

0.6

(deg)

(a) Position tracking error (b) Attitude tracking error Fig. 28. Statistical results of pose tracking errors (attitude-controlled tumbling target).

In this paper, the chaser’s initial attitude is set to a random number within a certain range, so that α(t0) (the initial angle between the direction of the target with respect to the chaser and the optical axis of camera) will be a uniformly distributed random numbers on [0 10] deg. And the elements of ω CC ( t 0 ) are all uniformly distributed pseudorandom numbers on [-0.25, 0.25] deg/s, see the paragraph after Eq. (34). The angle α may increase at the beginning because of inertia even if the chaser is controlled. If α increases to a certain value (e.g., 31 deg), the chaser could not observe the target. The maximum value of the increasing amount of α is denote by δα, and its frequency histogram is shown in Fig. 29. As shown in this figure, δα is mainly distributed on [0, 1] deg. Therefore, once α(t0) is less than 29 deg, α will be no more than 30 deg; i.e., once α(t0) is less than a certain value, the chaser will be always suitably oriented to observe the target. 1000

Frequency

Frequency

1000

500

0 0

0.2

0.4

0.6

0.8

500

0 0

1

0.2

0.4

0.6

0.8

1

(deg)

(deg)

(a) Uncontrolled tumbling target (b) Attitude-controlled tumbling target Fig. 29. Statistical results of δα.

After the chaser arrives at the point S, the maximum value of α is denoted by αsm, and its frequency histogram is shown in Fig. 30. As shown in this figure, αsm is mainly distributed on [0.2, 0.7] deg. It means that the optical axis of the chaser’s camera is almost point to the target, and this viewing angle is excellent. In summary, if our method is applied in the phase of close range rendezvous, the chaser will be always suitably oriented to observe the target, and the viewing angle will be excellent after the chaser arrives at the point S.

31 / 41

400

150

300

Frequency

Frequency

200

100 50 0 0.2

0.3

0.4 sm

0.5

200 100 0

0.6

(deg)

0.3

0.4 sm

0.5

0.6

(deg)

(a) Uncontrolled tumbling target (b) Attitude-controlled tumbling target Fig. 30. Statistical results of αsm.

6

Conclusion In this paper, a novel control scheme based on motion planning and pose tracking is proposed

to solve the pose control of a chaser docking with a tumbling target in the phase of close range rendezvous. Unlike the existing literature, in this paper, even if the initial direction of the chaser with respect to the target is arbitrary, our method always works. The proposed scheme is: firstly, the current desired motion of the chaser is planned in real time according to the motion of the target, then the difference between the chaser’s motion and the desired motion is gradually reduced by a pose tracking controller. It should be noted that, on the one hand, in planning the desired motion, the chaser’s initial conditions are considered, the planned velocity is continuous, and the planned acceleration is within a small range, namely the desired motion is easy to track; on the other hand, the NTSM controller is used to track the desired motion, so the tracking errors converge to about zero in finite time. Therefore, the desired motion can be tracked well by the chaser. Moreover, the desired motion has a series of advantages: collision avoidance, guarantee of excellent line-of-sight angle, and reasonable velocity and acceleration. Therefore, it is not necessary to constrain the chaser’s motion when tracking the desired motion. In summary, (1) the challenging 6-DOF constrained control problem as in (Dong et al., 2018) is avoided, namely the design of the pose tracking controller is greatly simplified; (2) the safe arrival of the chaser at the docking port of the target with a suitable relative attitude can be guaranteed; (3) the chaser will be always suitably oriented to observe the target well; (4) the chaser’s velocity will be reasonable (very close to the desired one); and (5) the magnitude of the needed control inputs is relatively small.

32 / 41

Appendix A: Earth Oblateness Effect In this section, the relative effect of the Earth oblateness due to J2 is calculated. For an Earth satellite, the perturbing acceleration due to the J2 effect in the inertial frame ΣN is described as  x 1  5( z / r ) 2   3J 2 R  aJN2 ( r N )   y 1  5( z / r ) 2   5 2r  z 3  5( z / r ) 2    2 e

(A1)

where rN = [x, y, z]T is the position vector of the satellite with respect to ΣN, J2 = 0.0010826267, Re = 6378.137 km, and μ = 3.986032×1014 m3/s2. Therefore, the aJH2 in Eq. (1) can be computed as

aJH  AHN [aJN (rCN )  aJN (rTN )] 2

2

(A2)

2

N where rC is the position vector of the chaser with respect to ΣN, rTN is that of the target, and AHN

represents the direction cosine matrix from ΣN to ΣH.

Appendix B: Convergence of System States in the NTSM Manifold In this section, it is proved that in the NTSM manifold of this paper, the system states converge to zero in finite time. In the NTSM manifold of this paper, the following equations should be held q/ p   x   x      x  b1 xp / q  0   b1      s( t )  0     q/ p p/q  x  x  b4 x  0    x    b    4   

(B1)

where [ x , y , z ]T  Θr , [ x ,  y , z ]T  ωrC , [ x, y , z ]T  ρrH , and [ x, y , z]T  vrH . Substituting Eq.(B1) into Eq.(5), one can obtain   [x , y , z ]T  A(Θr )   x  b1

  

q/ p

 y    b2

  

q/ p

 z    b3

  

q/ p

  

T

(B2)

Meanwhile,  r   x2   y2   z2 is the norm of Θ r , its time derivative is

r  (xx   yy  zz ) / r

(B3) 33 / 41

Substituting Eq.(B2) into the above equation, one can obtain q/ p q/ p q/ p  z    y  1   x    r   x      y      z      r   b1   b2   b3   

(B4)

Note that p, q are positive odd integers, and b1 > 0, so  x   x / b1 

q/ p

namely r  0 . Let 1  q / p b1   x  r   r  

1

r

=  b1 

q / p

q / p

( p  q) / p

 0,

i  kr , k [ 3 / 3, 1], here i  max{ x ,  y , z } , which gives

( pq) / p

 b2 

q / p

y

( pq) / p

 b3 

q / p

z

( pq) / p

 

bi  q / p  i ( p  q ) / p

 k( p  q ) / p bi 

x

(B5)

 rq / p

 ( 3 / 3) ( p  q ) / p bi 

q / p

 rq / p

It is known that, in a standard NTSM manifold, s(t )  x  bx p / q  0 , (b > 0), so

dx  kx q / p dt

(B6)

where k = (b)-q/p > 0. If

x0

, then

x  q / p d x  k d t





p t t t t d x ( p  q ) / p  k t t  t pq p x (t0  t f ) ( p  q ) / p  x (t0 ) ( p  q ) / p  kt f pq 0

0

f

0

f

(B7)

0





Since x(t0 + tf) = 0, tf could be expressed as

tf 

p x ( t0 ) ( p  q ) / p k ( p  q)

(B8)

Therefore, in a standard NTSM manifold, the state converges to zero in finite time tf. In this paper, the last three equations in Eq.(B1) belongs to the equations of standard NTSM manifold. Namely the last three states converge to zero in finite time. As for the first three, from Eq. (B5), q / p we know  r  [( 3 / 3) ( p  q ) / p bi  ]   rq / p , by contrasting it with Eqs. (B6) and (B8), it can be

concluded that the time for θr to converge to zero will be

t f , 

p  r ( t0 ) ( p  q ) / p k ( p  q)

(B9)

34 / 41

where k  ( 3 / 3) ( p  q ) / p bi 

q / p

. Namely the first three states converge to zero in finite time, too.

Appendix C: The Time for Reaching the NTSM Manifold In this section, it is proved that the value of the Lyapunov function (Eq. (17)) converges to zero (namely reaching the NTSM manifold s(t) = 06×1) in finite time. From Eq. (22), we know 6

V   i si   j s j , here s j  max(s1 , ... , s6 ) . Meanwhile, from Eq. (17), we know that i 1

V  ( s12  ...  s62 ) / 2 , so

s j  ks V , ks [ 3 / 3, 2 ]

(C1)

（ j V）. Reexamining Eqs. (B6) and (B8), let q/p = 0.5, it can concluded Therefore, V   3/3 that if x   k x , the time for x to converge to zero is 2/k ( x0 ) . Therefore, the time for V to converge to zero is no more than 2 3V0 / j , here V0 is the initial value of V, and j may be 1, 2, …, or 6.

Appendix D: Direction of the Best Value of the Desired Velocity In this section, the direction of the best value vopt of the desired velocity is calculated. As shown in Fig. D1 and Fig. D2, the line DP is the tangent of the danger zone, and the triangle TDP is on the plane xTp-T-yTp. zT

danger zone

yT (yTp)

T T xT

yT (yTp)

F

S

xTp

P θ0

D

Fig. D1.

rs

ρD

βT

θ1

D

Direction of vopt in ΣT.

Fig. D2.

35 / 41

P θ2

xTp

Direction of vopt in ΣTp.

As shown in Fig. D1, The frame ΣTp is defined below. If ΣT rotates an angle of βT about the axis yT, its attitude will be the same as that of ΣTp. So βT can be determined by the position of D as 0 ( x DT  0, z DT  0)  ( x DT  0, z DT  0)  / 2   T    / 2 ( x DT  0, z DT  0)  T T ( x DT  0) arctan(  z D / x D ) arctan(  z DT / x DT )   ( x DT  0) 

(D1)

where ρDT  [ xDT , yDT , zDT ]T is the position vector of D in ΣT. Then the direction cosine matrix from ΣTp to ΣT can be calculated as A

TTp

 cos  T  0    sin  T

0 sin  T  1 0   0 cos  T 

(D2)

T If y D  rs , the direction of vopt is along DS, or the direction is along the tangent DP. The

direction of DP can be solved below. As shown in Fig. D2, θ0 = arccos( y DT /ρD), θ1 = arcsin(rs/ρD), Tp so θ2 =   0  1 can be obtained now. Then unit vector of vopt can be obtained as: e v = [T TTp Tp sinθ2, cosθ2, 0]T, and e v can be computed as A ev .

Acknowledgments This work was supported by the Natural Science Foundation of China [grant numbers 11772187, 11802174], the China Postdoctoral Science Foundation [grant number 2018M632104], and the research project of the Key Laboratory of Infrared System Detection and Imaging Technology of Chinese Academy of Sciences [grant number CASIR201702]. The authors greatly appreciate their financial support.

References Agarwal, K., Weiss, A., Kolmanovsky, I., Bernstein, D., 2012. Inertia-Free Spacecraft Attitude Control with Control Moment Gyroscope Actuation. AIAA Guidance, Navigation, and Control Conference. American Institute of Aeronautics and Astronautics. Bandyopadhyay, S., Chung, S.-J., Hadaegh, F.Y., 2016. Nonlinear Attitude Control of Spacecraft with a Large Captured Object. Journal of Guidance, Control, and Dynamics. 39(4), 754-769.

36 / 41

Breger, L.S., How, J.P., 2008. Safe Trajectories for Autonomous Rendezvous of Spacecraft. Journal of Guidance, Control, and Dynamics. 31(5), 1478-89. Cruz, G., D'Amato, A., Bernstein, D., 2012. Retrospective Cost Adaptive Control of Spacecraft Attitude. AIAA Guidance, Navigation, and Control Conference. American Institute of Aeronautics and Astronautics. Cruz, G., Yang, X., Weiss, A., Kolmanovsky, I., Bernstein, D., 2011. Torque-saturated, InertiaFree Spacecraft Attitude Control. AIAA Guidance, Navigation, and Control Conference. American Institute of Aeronautics and Astronautics. Dalla, V. K., Pathak, P. M., 2015. Trajectory tracking control of a group of cooperative planar space robot systems. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering. 229(10), 885-901. Dalla, V. K., Pathak, P. M., 2019. Impedance control in multiple cooperative space robots pulling a flexible wire. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science. 233(6), 2190-2205. Dalla, V. K., Pathak, P. M., 2019. Power-optimized motion planning of reconfigured redundant space robot. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering. 233(8), 1030-1044. Dong, H., Hu, Q., Akella, M.R., 2017. Safety Control for Spacecraft Autonomous Rendezvous and Docking Under Motion Constraints. Journal of Guidance, Control, and Dynamics. 40(7), 1680-1692. Dong, H., Hu, Q., Akella, M.R., 2018. Dual-Quaternion-Based Spacecraft Autonomous Rendezvous and Docking Under Six-Degree-of-Freedom Motion Constraints. Journal of Guidance, Control, and Dynamics. 41(5), 1150-1162. Feng, Y., Bao, S., Yu, X., 2002. Design method of non-singular terminal sliding mode control systems. Control and Decision. 17(2), 194-198. Feng, Y., Yu, X., Man, Z., 2002. Non-singular terminal sliding mode control of rigid manipulators. Automatica. 38(12), 2159-2167. Filipe, N., Tsiotras, P., 2014. Adaptive Position and Attitude-Tracking Controller for Satellite Proximity Operations Using Dual Quaternions. Journal of Guidance, Control, and Dynamics. 38(4), 566-577. Flores-Abad, A., Ma, O., Pham, K., Ulrich, S., 2014. A review of space robotics technologies for on-orbit servicing. Progress in Aerospace Sciences. 68, 1-26. Gao, W., Zhu, X., Zhou, M., Jiang, Q., 2016. ADRC Law of Spacecraft Rendezvous and Docking in Final Approach Phase. Cybernetics and Systems. 47(3), 236-248. Kawamoto, S., Nishida, S., Kibe, S., 2003. Research on a space debris removal system. NAL Res. 37 / 41

Prog. (Nat. Aerosp. Lab. Jpn.). 2002/2003, 84-87. Kjellberg, H.C., Lightsey, E.G., 2013. Discretized Constrained Attitude Pathfinding and Control for Satellites. Journal of Guidance, Control, and Dynamics. 36(5), 1301-1309. Lee, D., Bang, H., Butcher, E.A., Sanyal, A.K., 2014. Nonlinear Output Tracking and Disturbance Rejection for Autonomous Close-Range Rendezvous and Docking of Spacecraft. Transactions of the Japan Society for Aeronautical and Space Sciences. 57(4), 225-237. Lee, D., Sanyal, A.K., Butcher, E.A., 2015. Asymptotic Tracking Control for Spacecraft Formation Flying with Decentralized Collision Avoidance. Journal of Guidance, Control, and Dynamics. 38(4), 587-600. Lee, D., Vukovich, G., 2016. Robust adaptive terminal sliding mode control on SE(3) for autonomous spacecraft rendezvous and docking. Nonlinear Dynamics. 83(4), 2263-2279. Lee, U., Mesbahi, M., 2014. Dual Quaternion based Spacecraft Rendezvous with Rotational and Translational Field of View Constraints. AIAA/AAS Astrodynamics Specialist Conference. American Institute of Aeronautics and Astronautics. Li, P., Zhu, Z.H., 2018. Model predictive control for spacecraft rendezvous in elliptical orbit. Acta Astronautica. 146, 339-348. Li, Q., Yuan, J., Wang, H., 2018. Sliding mode control for autonomous spacecraft rendezvous with collision avoidance. Acta Astronautica. 151, 743-751. Li, Q., Yuan, J., Zhang, B., Gao, C., 2017. Model predictive control for autonomous rendezvous and docking with a tumbling target. Aerospace Science and Technology. 69, 700-711. Ma, Z., Ma, O., Shashikanth, B.N., 2007. Optimal approach to and alignment with a rotating rigid body for capture. The Journal of the Astronautical Sciences. 55(4), 407-419. Man, Z., Xing Huo, Y., 1997. Terminal sliding mode control of MIMO linear systems. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications. 44(11), 10651070. Pong, C.M., Miller, D.W., 2015. Reduced-Attitude Boresight Guidance and Control on Spacecraft for Pointing, Tracking, and Searching. Journal of Guidance, Control, and Dynamics. 38(6), 1027-1035. Richards, A., Schouwenaars, T., How, J.P., Feron, E., 2002. Spacecraft Trajectory Planning with Avoidance Constraints Using Mixed-Integer Linear Programming. Journal of Guidance, Control, and Dynamics, 25(4), 755-64. Sanyal, A., Chaturvedi, N., 2008. Almost Global Robust Attitude Tracking Control of Spacecraft in Gravity. AIAA Guidance, Navigation and Control Conference and Exhibit. American Institute of Aeronautics and Astronautics. Sun, L., Huo, W., Jiao, Z., 2017. Adaptive nonlinear robust relative pose control of spacecraft 38 / 41

autonomous rendezvous and proximity operations. ISA Transactions. 67, 47-55. Sun, L., Zheng, Z., 2017. Adaptive relative pose control for autonomous spacecraft rendezvous and proximity operations with thrust misalignment and model uncertainties. Advances in Space Research. 59(7), 1861-1871. Sun, L., Zheng, Z., 2018. Adaptive relative pose control of spacecraft with model couplings and uncertainties. Acta Astronautica. 143, 29-36. Tekinalp, O., Tekinalp, A., 2016. Tracking Control of Spacecraft Attitude on Time Dependent Trajectories. AIAA/AAS Astrodynamics Specialist Conference. American Institute of Aeronautics and Astronautics. Ventura, J., Ciarcià, M., Romano, M., Walter, U. Fast and Near-Optimal Guidance for Docking to Uncontrolled Spacecraft. Journal of Guidance, Control, and Dynamics. 40(12), 3138-54. Weiss, A., Baldwin, M., Erwin, R.S., Kolmanovsky, I., 2015. Model Predictive Control for Spacecraft Rendezvous and Docking: Strategies for Handling Constraints and Case Studies. IEEE Transactions on Control Systems Technology. 23(4), 1638-47. Yoon, H., Tsiotras, P., 2005. Adaptive Spacecraft Attitude Tracking Control with Actuator Uncertainties. AIAA Guidance, Navigation, and Control Conference and Exhibit. American Institute of Aeronautics and Astronautics. Figure captions in a separate list Fig. 1.

Approach path constraint and field-of-view constraint (Dong et al., 2018).

Fig. 2.

Two situations of chaser’s relative position with respect to target.

Fig. 3.

ISS size information and axes definition (Dong et al., 2017).

Fig. 4.

Some frames in this paper.

Fig. 5.

The pose of ΣC when docking.

Fig. 6.

Direction of vopt in ΣT.

Fig. 7.

The pose of ΣD when it just arrives at S.

Fig. 8.

Relationship between rotational and translational motions.

Fig. 9.

Position tracking error.

Fig. 10.

Attitude tracking error.

Fig. 11.

Line-of-sight angle α(t).

Fig. 12.

Translational motion of the chaser in ΣT.

Fig. 13.

Translational motion of the chaser in ΣH.

Fig. 14.

Translational motion of the chaser in ΣT.

Fig. 15.

Control force. 39 / 41

Fig. 16.

Control torque.

Fig. 17.

Position tracking error.

Fig. 18.

Attitude tracking error.

Fig. 19.

Line-of-sight angle α(t).

Fig. 20.

The desired and real pose of the chaser.

Fig. 21.

Control inputs.

Fig. 22.

Translational motion of the chaser in ΣT.

Fig. 23.

Position tracking error.

Fig. 24.

Attitude tracking error.

Fig. 25.

Line-of-sight angle α(t).

Fig. 26.

Statistical results of dm.

Fig. 27.

Statistical results of pose tracking errors (uncontrolled tumbling target).

Fig. 28.

Statistical results of pose tracking errors (attitude-controlled tumbling target).

Fig. 29.

Statistical results of δα.

Fig. 30.

Statistical results of αsm.

Fig. D1.

Direction of vopt in ΣT.

Fig. D2.

Direction of vopt in ΣTp.

40 / 41

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

41 / 41