Adaptive Fault-Tolerant Spacecraft Pose Tracking

IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY

1

Adaptive Fault-Tolerant Spacecraft Pose Tracking with Control Allocation Haichao Gui and Anton H. J. de Ruiter, Member, IEEE

Abstract—The concurrent position and attitude (pose) tracking of a rigid spacecraft is addressed in the presence of actuator faults, mass and inertia uncertainties, and unknown external disturbances. The control design relies on a novel hybrid dualquaternion integral sliding mode that incorporates a hysteretic switching to avoid the quaternion unwinding problem. On the sliding mode, the pose tracking error is globally finite-time convergent. The resultant control law has a simple structure and consists of a nominal control input, which realizes the sliding mode in the fault-free uncertainty-free case, and an adaptive part, which dynamically compensates the perturbations due to actuator faults and other uncertainties. Moreover, a real-time control allocation algorithm is devised to deliver the control command to position and attitude actuators in proportion to the effectiveness degree of each actuator. Rigorous proof shows that the proposed strategies can stabilize the spacecraft pose trajectory to a small neighborhood of the sliding mode in finite time. Formation flying of Earth-orbiting spacecraft and asteroid body-fixed hovering are simulated to demonstrate the efficacy and broad applicability of the proposed method. Index Terms—Adaptive control, control allocation, faulttolerant control, integral sliding mode, spacecraft attitude and position tracking.

I. I NTRODUCTION ONCURRENT position and attitude (pose) tracking of a spacecraft is gaining increasing interest because it provides an enabling technology for many space close proximity operations such as on-orbit monitoring and/or assembly, refueling, rendezvous and docking, space-based interferometry, small-body hovering, descent and landing, etc. The coupled 6degree-of-freedom (6DOF) spacecraft rigid-body motion must be fully taken into account for control design and stability analysis, unlike traditional translation- or rotation-only control issues. Many 6DOF spacecraft motion controllers were designed with asymptotic [1]–[5] or finite-time [6]–[8] convergence, using full-state feedback [1]–[7] or output feedback [8]–[11]. Additionally, techniques range from indirect or direct adaptive control [1], [3]–[5], [10], terminal sliding mode (TSM) [6], [7], immersion and invariance theory [12], and disturbance observer [13] were utilized to deal with uncertain mass and inertia as well as unknown disturbances. Among many algebraic methodologies such as quaternions or modified Rodrigues parameters (MRP) plus the position vector [1], [2], [9], [10],

C

H. Gui is with the School of Astronautics, Beihang University, Beijing, P. R. China (email: [email protected]) A. H. J. de Ruiter is with the Department of Aerospace Engineering, Ryerson University, Toronto, M5B 2K3, Canada (email: [email protected]). Manuscript received on March XX, 2017. This work is supported by the National Natural Science Foundation of China under projects 11702010 and 11432001, and by the Natural Sciences and Engineering Research Council of Canada through a Discovery Accelerator Supplement.

[12], [13], special Euclidean group SE(3) [14], [15] and twistors [16], dual quaternions, as an extension of quaternions, have emerged as an efficient tool for pose control design because they minimally globally describe the rigid-body pose while capturing the natural coupling between translation and rotation [3]–[8], [17]. Nonetheless, they inherit the ambiguities intrinsic in quaternions and form a double covering of SE(3). This property, if inappropriately treated, can cause unwinding problems like some quaternion-based controllers. The above studies all assumed healthy position and attitude control actuators that can produce the exact command force and torque as required. On-orbit spacecraft, however, can encounter unexpected actuator faults or failures. For example, the Hayabusa spacecraft carried four ion engines and three reaction wheels (RWs) but two RWs failed in 2005 and one ion engine became deteriorated in 2009 [18], [19]. More incidents of spacecraft actuators were reported in [20]. Hence, it is highly desirable to develop robust control strategies that can retain satisfactory performance when actuator faults or failures occur. Various fault-tolerant control (FTC) methods, either passive or active types, were obtained for attitude tracking [21]–[25]. In contrast, little attention has been paid to the equally important 6DOF FTC problem. Recently, Dong et al. [26] derived a dual-quaterion-based FTC algorithm using time-varying TSM. The design depends on Euler-Lagrange formulation of the pose dynamics with generalized coordinates being the vector part of the unit dual quaternion. The resultant controller thus is not globally effective since it does not allow the scalar part of the attitude quaternion to be zero. When the spacecraft is over-actuated with redundant actuators, control allocation (CA) is a useful approach to improve reliability and power efficiency. It aims to distribute the virtual control efforts to individual actuators in a proper manner, e.g., minimizing certain cost functions [27]. Model-dependent sliding mode control allocation algorithms were constructed in [28], [29] for FTC of linear systems and later extended by Shen et al. [30] for CA-based inertia-free fault-tolerant attitude tracking. Recently, torque allocation was formulated in [31] as a min-max optimization problem and solved numerically. This paper studies simultaneous position and attitude tracking of a rigid spacecraft in the framework of dual quaternion algebra. A novel combination of finite-time control, hybrid system techniques, integral sliding mode (ISM), adaptive control and CA is proposed for pose tracking with actuator faults, mass and inertia uncertainties, and unknown external disturbances. The design procedures start from a hybrid dualquaternion ISM (HDQISM) which encapsulates a globally finite-time stable nominal system. A virtual adaptive ISM controller (AISMC) is then derived and consists of a nominal


2

control input, which realizes the sliding mode in the fault-free uncertainty-free case, and an adaptive part, which dynamically compensates the total uncertainties according to the deviation from the sliding motion. The virtual controller is then coupled with a CA algorithm, which uses the measured or estimated effectiveness degree of the actuators to redistribute the control to individual actuators. The main contributions are threefold: (1) the HDQISM incorporates a hysteretically switched variable to dynamically select the shortest rotation toward the target attitude and thus avoids the unwinding problem; (2) The adaptive law ensures a priori bounded adaptive variables and thus prevents them from drifting arbitrarily away from the true values; (3) Rigorous proof shows that the proposed method can stabilize the pose tracking error to a small neighborhood of the HDQISM in finite time, despite actuator faults and multiple uncertainties. Thereafter, the real pose trajectory behaves closely to the nominal closed-loop performance without reconfiguring the controller. Spacecraft 6DOF proximity operations such as monitoring and circumnavigation of Earth-orbiting target and body-fixed hovering about asteroids are presented to illustrate the broad applicability of the proposed method. II. P RELIMINARIES Throughout this paper, denote by I n the n × n identity matrix, 0n = [0, ..., 0]T ∈ Rn , and In = {1, ..., n}. For any x ∈ R and α ≥ 0, let sgnα (x) = sgn(x)|x|α and satα (x) = sgn(x)min{|x|α , 1}, where sgn(·) is the standard sign function. Clearly, sgnα (x) is a continuous nonsmooth function if 0 < α < 1, while satα (x) becomes the standard saturation function sat(x) if α = 1. For x ∈ Rn and α ≥ 0, let sgnα (x) = [sgnα (x1 ), ..., sgnα (xn )] and satα (x) = [satα (x1 ), ..., satα (xn )]. Denote by k·k the Euclidean norm of a vector and the induced 2-norm of a matrix. For A ∈ Rm×n , let σ ¯ (A) and σ(A) be its maximum and minimum singular values respectively. Given x ∈ R3 , x× is the skew-symmetric matrix satisfying x× y = x × y, ∀y ∈ R3 , where × is the cross product on R3 . Given square matrices Ai , i ∈ In , of any dimension, denote by D = diag(A1 , · · · , An ) and A = adiag(A1 , · · · , An ) respectively the block diagonal and block anti-diagonal matrices generated from {Ai : i ∈ In }: 

0

A1 ..

 D=

0

  ,

.



0 ..

 A=

An

A1

An

1 0 (q ⊗ p − p∗ ⊗ q ∗ ) = ¯ + p0 q ¯+q ¯×p ¯ q0 p 2

which are both associative and distributive but are not commutative. Additionally, (q⊗p)∗ = p∗ ⊗q ∗ and (q×p)∗ = p∗ ×q ∗ . Letting 1 = [1, 0, 0, 0]T be the identity in Q, the set of unit quaternions is then defined by QU = {q ∈ Q : q ⊗ q ∗ = 1}. Dual quaternions are extensions of quaternions by means of a dual unit which satisfies 2 = 0 but 6= 0. A dual quaterˆ = q r + q d , where q r = [qr0 , q ¯ Tr ]T ∈ Q nion is denoted as q T T ¯ d ] ∈ Q are called the real and dual and q d = [qd0 , q ˆ are parts respectively [32]. The scalar and vector parts of q ¯ r + ¯ ˆ is given by qr0 + qd0 and q q d . The conjugate of q ˆ ∗ = q ∗r + q ∗d . Denote by DQ the set of dual quaternions q and define the associative, distributive, noncommutative dual ˆ and any quaternion product and cross product between q ˆ = pr + pd ∈ DQ as p ˆ⊗p ˆ = q r ⊗ pr + (q r ⊗ pd + q d ⊗ pr ) q ˆ ×p ˆ= q

1 ˆ −p ˆ∗ ⊗ q ˆ ∗ ) = q r × pr + (q r × pd + q d × pr ). (ˆ q ⊗p 2

ˆ ∗ ⊗ˆ The above definitions imply (ˆ q ⊗ˆ p)∗ = p q ∗ and (ˆ q ×ˆ p)∗ = ∗ ∗ ˆ ˆ ×q ˆ . Letting 1 = 1 + 04 be the identity in DQ, the set of p unit dual quaternions is then defined by DQU = {ˆ q ∈ DQ : ˆ ˆ⊗q ˆ ∗ = 1}. q Similarly, dual vectors, as an extension of 3-D vectors, can ˆ = q r + q d , where q r , q d ∈ R3 . Denote by be defined as q DQV the set of dual vectors. As a 3-D vector can be viewed as a quaternion with zero scalar part, dual vectors are dual quaternions with zero scalar parts. This is an implied fact when performing conjugation, dual quaternion product, and cross product for elements in DQV . By stacking its real part and dual part in order, each dual quaternion (vector) corresponds a unique column vector. Hence, DQ (DQV ) is isomorphic to R8 (R6 ). In this paper, ˆ = q r + q d (∈ DQ or ∈ DQV ) is also treated as q ˆ = [q Tr , q Td ]T for convenience in algebraic computations. q Following this, all operations on the Euclidean space (e.g., the Euclidean norm, matrix product, the functions sgn(·) and ˆ sat(·) previously defined, etc.) can also be applied to any q (∈ DQ or ∈ DQV ). In addition, the functions re(ˆ q ) and du(ˆ q) ˆ. respectively extract the real part and dual part of q

  .

.

q×p=

0

A. Quaternions and Dual Quaternions Denote by Q the set of quaternions, which is isomorphic to ¯ T ]T ∈ Q, R4 . A quaternion can be represented as q = [q0 , q ¯ are called the scalar and vector parts of q, where q0 and q respectively. The conjugate of q is q ∗ = [q0 , −¯ q T ]T . The quaternion product and cross product between q and any p = ¯ T ]T ∈ Q are defined as [p0 , p ¯T p ¯ q0 p0 − q q⊗p= ¯ + p0 q ¯+q ¯×p ¯ q0 p

B. Spacecraft Dynamics in Dual Quaternions Dual quaternions are a favorable tool to compactly encapsulate relative position and attitude, and thus provide a unified description of the spacecraft 6DOF motion. First, let r yz , v yz , and ω yz be the linear displacement, linear velocity, and angular velocity of an orthogonal coordinate frame, say a Y frame relative to a Z frame and their components in an arbitrary orthogonal coordinate frame, say the X frame, are given by ?xyz ∈ R3 , where ? ∈ {r, v, ω}. The attitude of the Y frame relative to the Z frame is described by q yz ∈ QU . Clearly, the Z frame can be transformed to the Y frame by a rotation q yz followed by a translation r yz or vice versa. The pose of the Y frame relative to the Z frame can then be represented by a unit dual quaternion as [7]


3

1 1 ˆ yz = q yz + q yz ⊗ r yyz = q yz + r zyz ⊗ q yz . (1) q 2 2 The dual velocity of the Y frame with respect to the Z frame expressed in the X frame is defined as ˆ xyz = ω xyz + (v xyz + ω xyz × r xxy ) ω

(2)

where v xyz = Rxy v yyz = Rxy (r˙ yyz + ω yyz × r yyz ) and Rxy is the rotation matrix from the Y frame to the X frame. By means of q xy ∈ QU , Rxy is computed as 2 ¯ Txy q ¯ xy )I 3 + 2¯ ¯ Txy − 2qxy0 q ¯× Rxy = (qxy0 −q q xy q xy . w ˆ xyz and ω ˆw Given yet another frame, say W , ω yz = ω yz + w w w (v yz + ω yz × r wy ) satisfy the following transformation formula [4]

ˆ ∗xw ⊗ ω ˆw ˆ xw . ˆ xyz = q ω yz ⊗ q

(3)

The above general definitions make it fairly easy to obtain dual quaternions and dual velocities for a spacecraft in different reference frames. Denote by FI an inertial frame, FB the spacecraft body frame with its origin attached to the center-ofmass (CM) of the spacecraft, and FD the frame representing the desired pose of the spacecraft. We can derive the unit ˆ bi and dual quaternions of FB and FD relative to FI as q ˆ di , and FB relative to FD as q ˆ bd by (1), as well as the q ˆ bi , ω ˆ di and ω ˆ bd by (2). Noting corresponding dual velocities ω q bd = q ∗di ⊗ q bi , ω bbd = ω bbi − ω bdi , and r bbd = r bbi − r bdi , it can ˆ bd = q ˆ ∗di ⊗ q ˆ bi and ω ˆ bbd = ω ˆ bbi − ω ˆ bdi , where be verified that q b ∗ d ˆ di = q ˆ bd ⊗ ω ˆ di ⊗ q ˆ bd according to (3). ω Let m > 0 and J = J T > 0 be the spacecraft mass and inertia matrix. They are utilized to define a 6 × 6 block antidiagonal matrix as 03×3 mI 3 M = adiag(mI 3 , J ) = (4) J 03×3 which can be viewed as a matrix counterpart of the dual inertia defined in [4], [7] and has lower dimension than the 8×8 dual inertia matrix introduced in [3], [17]. Assume that the total dual force acting on the spacecraft includes the dual control b b force fˆ c ∈ DQV plus the dual environment force fˆ e ∈ DQV . By means of M , the pose kinematics and dynamics of the spacecraft relative to FD developed in [4] can be rewritten as 1 ˆ ⊗ω ˆ bbd ˆq˙ bd = q 2 bd

(5)

ˆ bbd − ω ˆ bdi × M ω ˆ bdi Mω ˆ˙ bbd = S(ω bbd , ω bdi )ω b b ˆ ) + fˆ + fˆ ˆ˙ d ⊗ q −M (ˆ q∗ ⊗ ω bd

di

bd

c

(6)

e

where S(ω bbd , ω bdi ) takes the following form: S(ω bbd , ω bdi ) = adiag(S p (ω bbd , ω bdi ), S a (ω bbd , ω bdi )) S p (ω bbd , ω bdi ) = −m(ω bbd + ω bdi )× S a (ω bbd , ω bdi ) = (J (ω bbd + ω bdi ))× − (ω bdi )× J − J (ω bdi )× . (7)

The spacecraft dynamics relative to FI can be readily derived through (5) and (6) by setting FI as the desired frame. As such, it follows that b b 1 ˆ bi ⊗ ω ˆ bbi , M ω ˆ bbi × M ω ˆ bbi + fˆ c + fˆ e . ˆ˙ bbi = −ω ˆq˙ bi = q 2 b The dual control force fˆ c = f bc + τ bc consists of the b 3 control force f c ∈ R and torque τ bc ∈ R3 provided by b position and attitude control actuators respectively. fˆ e depends on the dynamical environment where the spacecraft operates, and usually includes the gravity effect of the central body, solar radiation pressure, etc. In the following, we denote b ˆ bdi × M ω ˆ bdi + M (ˆ ˆ bd ) fˆ D = ω q ∗bd ⊗ ω ˆ˙ ddi ⊗ q

(8)

which is the dual force induced by the desired velocity and acceleration. b In realistic scenarios, m, J , and fˆ e can barely be known precisely and only nominal values, denoted by m0 , J 0 , and b fˆ e0 respectively, are available from modeling and/or experimental test. The following assumptions are reasonable in practice and utilized in the subsequent developments. Assumption 2.1: The spacecraft mass and inertia, the environmental forces and torques, the velocity and acceleration of the desired pose trajectory are all uniformly bounded in terms of the Euclidean norm. This also implies the boundedness of b the nominal quantities m0 , J 0 , and fˆ e0 . C. Actuator Faults and Control Allocation Assume that there are np (≥ 3) position actuators and na (≥ 3) attitude actuators fixed to the spacecraft. Their alignment matrices are D p ∈ R3×np and D a ∈ R3×na , respectively. Let 0 ≤ e?i (t) ≤ 1, ? ∈ {p, a}, represent the effectiveness degree of the ith position (? = p, i ∈ Inp ) or attitude (? = a, i ∈ Ina ) actuator. The ith actuator is fully functional when e?i (t) = 1, and fails completely if e?i (t) = 0; otherwise it loses a portion of its efficacy. Denote by E(t) = diag(ep1 (t), · · · , epnp (t), ea1 (t), · · · , eana (t)) the actuator effectiveness matrix, and D = diag(D p , D a ). The dual control force from actuators can then be written as b fˆ c = DE(t)fˆ u , fˆ u = [f Tu , τ Tu ]T np

na

(9)

where f u ∈ R and τ u ∈ R are command force and torque fed to the corresponding actuator units. Denote by ˜ E(t) = diag(˜ ep1 (t), · · · , e˜pnp (t), eã1 (t), · · · , eãna (t)), where 0 ≤ e˜?i (t) ≤ 1, ? ∈ {p, a}, the measured or estimated (e.g., from a real-time fault diagnosis algorithm) actuator effectiveness matrix. The following assumption is fundamental and restricts our study to the scenario that the spacecraft is always fully actuated in 6DOF despite actuator faults. Assumption 2.2: The effectiveness matrix and its estimate ˜ 3 (t)D T } = 6. satisfy rank{DE(t)D T } = 6 and rank{D E b Given a virtual control input fˆ vc ∈ DQV from a high-level control law, a CA algorithm is required to obtain fˆ u that drives the actuators. When either the position or attitude actuators are b redundant (np > 3 or na > 3), it is desirable to allocate fˆ vc


4

in proportion to the effectiveness degree of each actuator. In the following, we consider b ˜ 2 (t)D T (D E ˜ 3 (t)D T )−1 fˆ u = D + fˆ vc , D + = E

(10)

which solves the optimization problem: T −1 ˜ (t)fˆ u subject to fˆ b = D E(t) ˜ fˆ u . min fˆ u E vc

Clearly, the actuator with weak effectiveness is weighted heavily in the above cost function. If actuator faults occur and there is enough redundancy, the CA scheme given in (10) tends to allocate more control effort to the actuators with large effectiveness degree while minimizing the use of less healthy actuators. In fact, (10) is an extension of the torqueonly allocation method derived in [30]. Substituting (10) into (9) yields the actual control produced by the actuators:

and sgn(x) = sgn(x) for x 6= 0. The following virtual control input is then proposed: ) b b b b fˆ vc0 = fˆ k + fˆ D0 − fˆ e0 x∈C b ˆ (ˆ ˆ bbd ) fˆ k = −K 1 κ q bd , 1 − α1 , h) − K 2 satα2 (ω x∈D ˆ bbd , sgn(qbd0 )) x+ = (ˆ q bd , ω C = {x ∈ X : hqbd0 ≥ −δ} with D = {x ∈ X : hqbd0 ≤ −δ} (12) where K i = adiag(kpi I 3 , kai I 3 ) for kpi , kai > 0, i = 1, 2, b ˆ bdi × M 0 ω ˆ bdi + 0 < α1 < 1, α2 = 2α1 /(1 − α1 ), fˆ D0 = ω ∗ b d ˆ bd ), and x = (ˆ ˆ bd , h) ∈ X , DQU × M 0 (ˆ q bd ⊗ ω ˆ˙ di ⊗ q q bd , ω DQV × H; 0 < δ < 1 is a constant ; x+ denotes the state value immediately after a discontinuous jump and

b ˆ (ˆ κ q bd , 1 − α1 , h) = κ a (hq bd , 1 − α1 ) + κp (r bd , 1 − α1 ) ¯ h q bd  √ = (11) 1−α1 , if hqbd0 6= 1 2(1−hqbd0 ) κa (hq bd , 1 − α1 ) =  ˜ Denote by E δ (t) = E(t) − E(t) the mismatch between 03 , if hqbd0 = 1 b ˜ r E(t) and E(t), and M δ = M − M 0 the mismatch between 1 κp (r bbd , 1 − α1 ) = √ 2 b 2 kr b kbd1−α1 1+η kr bd k bd M and M 0 , where M 0 = adiag(m0 I 3 , J 0 ). Note that (13) γ1 = sup kDE δ (t)D + k and γ2 = sup kM δ M −1 DE(t)D + k where η > 0 is a constant. Note that κ ˆ (hˆ q bd , 1 − α1 ) t≥0 t≥0 is a continuous function of hˆ q bd for 0 < α1 < 1 since are both bounded and represent the maximal multiplicative limr b →0 kr bbd kα1 −1 r bbd = 03 and κ(hq bd , 1 − α1 ) is also 3 bd b perturbing effect on fˆ vc , induced by E δ (t) and M δ respec- continuous at hqbd0 = 1 [33]. Controller (12) is a hybrid of two modes, namely, a contintively. Clearly, γ1 and γ2 are proportional to the uncertainties flow E δ (t) and M δ and are zeros when actuator faults and the dual uous mode and a discrete mode governed by the so-called S inertia are precisely known. This perturbing effect should not set C and jump set D respectively (note that C D = X ). exceed the desired control effect, as assumed below, in order When x ∈ C, it produces a continuous control command b to be able to establish the closed-loop stability. fˆ vc0 which consists of a nonsmooth proportional-derivative Assumption 2.3: There exists a constant 0 < γ0 < 1 such feedback plus feedforward compensation terms. When x ∈ D, that γ1 + γ2 ≤ 1 − γ0 . it induces a state jump which reverses the sign of h and thus the spacecraft rotation direction toward the desired attitude ˆ bd and ω ˆ bbd continuous. From (12), it can while maintaining q D. Problem Formulation ˆ and q ˆ represent the alignment be seen that h changes its sign only when the sign mismatch ˆ bd = 1 ˆ bd = −1 Since both q between h and qbd0 exceeds the specified amount δ. Since of the spacecraft in the desired pose, our objective is to design 0 < δ < 1, the sign switching of h is hysteretic and b a virtual dual control force fˆ vc such that the tracking error makes a tradeoff between robustness to measurement noise b ˆ 06 ) ˆ bd ) is stabilized to (a small neighborhood of) (±1, (ˆ q bd , ω and selection of the shortest rotation path [34]. b b under the CA scheme given in (10). To this end, the following ˜ Setting fˆ vc = fˆ vc0 and recalling E(t) = E(t) , I np +na , assumption is made to clarify the state variables required for b b ˆ ˆ M = M 0 and f e = f e0 , the closed-loop flow dynamics control design. ˆ bi and (x ∈ C) become Assumption 2.4: The spacecraft can obtain its pose q ˆ bbi from the navigation module. In addition, dual velocity ω b ˆ bbd + fˆ k M 0ω ˆ˙ bbd = S 0 (ω bbd , ω bdi )ω (14) ˆ di , velocity ω ˆ ddi , and acceleration the desired pose trajectory q ω ˆ˙ ddi , which are usually generated by the spacecraft guidance where S 0 (ω b , ω b ) is obtained from (7) by replacing m and bd di module, are also known to the spacecraft. J with m0 and J 0 . The stability of the closed-loop system b fˆ c

b DE(t)D + fˆ vc .

III. N OMINAL C ONTROLLER D ESIGN W ITHOUT U NCERTAINTY First, consider the case that all actuators are fully functional ˜ and the system model is precisely known, i.e., E(t) = E(t) , b b ˆ ˆ I np +na , M = M 0 , and f e = f e0 . To account for the double-covering feature of unit dual quaternions, we introduce a binary variable h ∈ H , {−1, 1} and define sgn : R → H as an outer semicontinuous set-valued map, where sgn(0) ∈ H

under controller (12) is stated as follows. Theorem 3.1: Consider the 6DOF spacecraft tracking system given in (5), (6) and (9), and the CA scheme given in (10). b By employing the hybrid controller (12) and setting fˆ vc = b ˆ bbd ) is globally fˆ vc0 for x ∈ C, the pose tracking error (ˆ q bd , ω ˆ 06 ) in finite time. stabilized to (h1, Proof: See Appendix A. Remark 3.1: The hybrid pose tracking law given by (12) reduces to the hybrid finite-time attitude controller in [33] after


5

removing its translational part. In addition, it ensures uniform ˆ 06 ) if α1 = α2 = 1. global asymptotic stability of (h1, IV. ROBUST A DAPTIVE C ONTROLLER D ESIGN Controller (12) depends on zero-uncertainty assumption for global stability. In this section, it is further integrated with ISM and adaptive techniques to derive a new pose controller with increased robustness to actuator faults, mass and inertia uncertainties, and unknown disturbances. A. Integral Sliding Mode Inspired by the closed-loop flow dynamics given by (14), a novel HDQISM is designed in the form ˆ(t) = M 0 [ω ˆ bbd (t) − ω ˆ bbd (0)] − s

Z

t

b ˆ bbd (ν) + fˆ k (ν)]dν [S 0 (ν)ω

0

(15) where S 0 (ν) , S 0 (ω bbd (ν), ω bdi (ν)). The above integration is b well defined since fˆ k (ν) is Lebesgue integrable despite the ˆ(0) = 06 and s ˆ(t) discontinuous switches of h. In addition, s remains continuous when h switches its sign. ˆ(t) ≡ 06 , (15) is equivalent to (14). In other words, When s the reduced dynamics on the HDQISM is identical to those resulting from the hybrid pose controller (12). The following result follows from Theorem 3.1. ˆ(t) ≡ 06 , (ˆ ˆ bbd ) Corollary 4.1: On the sliding mode s q bd , ω ˆ 06 ) in finite time. globally converges to (h1, ˆ(t) ≡ 06 is Remark 4.1: Since the system trajectory on s finite-time convergent, the HDQISM is also a nonsingular TSM. A terminal ISM for attitude tracking was previously devised in [25] but the sliding motion is merely almost globally finite-time stable due to the presence of multiple, undesirable, unstable equilibria. As a result, sluggish transients can occur for the initial conditions close to or exactly residing on any of the undesirable equilibria. In contrast, the HDQISM overcomes these drawbacks and also the unwinding problem due to its advantageous global finite-time convergence feature. ˆ(t) ≡ 06 is When α1 = α2 = 1, the system trajectory on s globally asymptotically convergent according to Remark 3.1. B. Adaptive Integral Sliding Mode Controller ˆ(t) ≡ 06 from onset with zero Controller (12) ensures s uncertainty but can no longer achieve the same objective in the presence of system uncertainties. An intuitive way to enhance robustness is to introduce additional control inputs that adaptively reject the unknown perturbations. Hence, we propose a virtual dual control force in the form b b b fˆ vc = fˆ vc0 + fˆ vc1 b

(16)

where fˆ vc0 is given by (12), which is the equivalent control b on the sliding mode, and fˆ vc1 ∈ DQV aims at disturbance rejection and is to be designed later. ˆ, applying (6), (11), (12) and For x ∈ C, differentiating s (16) and invoking M 0 M −1 DED + = I 6 + DE δ D + − M δ M −1 DED + yields

b ˆ bbd + fˆ k ˆs˙ = M 0 ω ˆ˙ bbd − S 0 ω b b b ˆ bbd − fˆ D + DE(t)D + (fˆ k + fˆ D0 = M 0 M −1 [S ω (17) b b b b ˆ b + fˆ −fˆ + fˆ ) + fˆ ] − S 0 ω e0

vc1

e

bd

b = δˆ S + δˆ D + δˆ e + δˆ k + δˆ ξ + fˆ vc1

k

where ˆ bbd − S 0 ω ˆ bbd δˆ S = M 0 M −1 S ω b −1 ˆ bd − M δ M S ω ˆ bbd with S δ = S − S 0 = Sδ ω b b δˆ D = M 0 M −1 [DED + fˆ D0 − fˆ D ] b b b = M 0 M −1 [fˆ − fˆ + DE δ D + fˆ ] D0

D

D0

b b δˆ e = M 0 M −1 [fˆ e − DED + fˆ e0 ] b b b = M 0 M −1 [fˆ e − fˆ e0 − DE δ D + fˆ e0 ] b b δˆ k = M 0 M −1 DED + fˆ − fˆ k

b

k

(18)

b

= DE δ D + fˆ k − M δ M −1 DED + fˆ k b δˆ ξ = [DE δ D + − M δ M −1 DED + ]fˆ vc1 . Assumption 2.1 implies that the uncertain terms δˆ S , δˆ D , δˆ e , and δˆ k satisfy the following growth condition. Lemma 4.1: There exist constants ξp ≥ 0 and ξa ≥ 0 such that kre(δˆ S + δˆ D + δˆ e + δˆ k )k ≤ ξp $p kdu(δˆ S + δˆ D + δˆ e + δˆ k )k ≤ ξa $a where $p = 1 + kv bbd k + kv bbd kkω bbd k and $a = 1 + kω bbd k + kω bbd k2 . Proof: See Appendix B. Evidently, ξp and ξa are unknown due to uncertainties, but an estimation of their upper bounds can usually be derived with the available modeling information. It is assumed that ξ¯p > 0 and ξ¯a > 0 are known such that ξp < ξ¯p and ξa < ξ¯a . b Letting sp = re(ˆ s) and sa = du(ˆ s), we then design fˆ vc1 as follows b ˆ0 = [sTp0 , sTa0 ]T ∈ DQV fˆ vc1 = −γ0−1 Ξˆ s0 , s ˜ (kp + 2ξp )$p I 3 0 Ξ= 0 (ka + 2ξã )$a I 3 (19)  s?  ks? k if ks? k > µ? , µ? = µ$?0? , ? ∈ {p, a} s?0 =  s? if ks? k ≤ µ? , µ?

where k? > 0 and µ?0 > 0, ? ∈ {p, a}, are constants to be properly selected; γ0 is the constant specified in Assumption 2.3. The variables ξ˜p and ξã are respectively the dynamic estimates of ξp and ξa and updated by ˙ ξ˜? = ρ?1

ξ˜? 1− ¯ ξ?

! $? ks? k − ρ?2 ξ˜? , ? ∈ {p, a}

(20)

where ρ?1 , ρ?2 > 0 are constants. The initial estimate is chosen as 0 ≤ ξ˜? (0) ≤ ξ¯? . It can be seen that the above adaptive law compensates for the total system uncertainties according to the deviation from the sliding motion, which is characterized by ˆ. the magnitude of s


6

AISMC adaptive law d di

sˆ

d di

( qˆdi , ωˆ , ωˆ )

HDQISM

(x p , x a ) spacecraft

actuator faults

fˆvc

sˆ

control allocation

fû = D + fˆvcb

virtual dual control force

E( E (t )

fû fu

position actuators

τu

attitude actuators

fˆcb

fault diagnosis

( qˆbi , ωˆ bib )

( qˆbi , ωˆ bib ) Fig. 1: Architecture of the AISMC and control allocation algorithm

˙ ˙ As inferred from (20), ξ˜? < 0 if ξ˜? ≥ ξ¯? , ξ˜? ≥ 0 if ξ˜? = 0, ˙˜ and ξ? > 0 if ξ˜? < 0, ? ∈ {p, a}. By setting 0 ≤ ξ˜? (0) ≤ ξ¯? , the proposed adaptive law guarantees 0 ≤ ξ˜? (t) ≤ ξ¯? for all time and hence prevents the adaptive variable ξ˜? from drifting arbitrarily away from its true value. In contrast, these notable properties are not possessed by the adaptive pose controllers developed in [5], [10], [12], [13], [26]. Figure 1 shows the architecture of the proposed AISMC and CA algorithm. The behavior of the entire closed-loop system is described as follows. Theorem 4.1: Consider the 6DOF spacecraft tracking system given in (5), (6) and (9), and the CA scheme given in (10). By employing the AISMC given in (16), (19), and (20), the closed-loop trajectory is stabilized in finite time to a ˆ = 06 given by neighborhood of s η0 ˆ ∈ DQV : lim sup kˆ s s(t)k ≤ (21) η1 t→∞ where η1 = min {kp , ka } η0 =

η ξ¯ (kp +2ξ¯p )µp0 ρ ¯2 √1 p + + ρp2 ξ ρp1 4 p1 p ¯ ¯ (ka +2ξa )µa0 ρa2 ¯2 η1 ξ a + ρa1 ξa + √ρa1 + 4

(22)

(23)

Proof: See Appendix C. Due to actuator faults and various uncertainties, the closedˆ = 06 and, if drifting too loop trajectory can drift from s far, the adaptive law will introduce increased compensation ˆ = 06 . Through such fˆ vc1 to pull the trajectory back toward s mechanism, the proposed AISMC stabilizes the closed-loop ˆ = 06 and then approximately trajectory to a neighborhood of s recovers the behavior of the nominal system. Therefore, the control parameters can be tuned in a two-step procedure: first tune the parameters of fˆ vc0 to set up the desired performance of the nominal system and then choose the parameters of fˆ vc1 to determine the speed to recover nominal system behavior. Large values of ρp1 and ρa1 can lead to fast adjustment and possibly increased overshoots, of the adaptive variables. As shown in (21)-(23), increasing kp and ka can reduce ˆ. In addition, decreasing µp0 and the ultimate bound on s

µa0 shrinks the boundary layer on the closed-loop trajectory and hence increases the pose tracking accuracy but also the likelihood of chattering. Therefore, the selection of µp0 and µa0 must trade off between the steady-state accuracy and chattering effect of the control. Remark 4.2: A stuck actuator usually gives rise to zero or constant output, which can be viewed as a constant disturbance. It is noteworthy that the proposed method can also handle actuator stuck as long as the spacecraft is still fully actuated, i.e., three control forces and three control torques remain available. A numerical example illustrating this type of fault will be shown in Section V.C. Unlike the faulttolerant scheme in [35], the methods in this paper does not take the actuator dynamics into account. In other words, we assume that the actuators are responsive enough and can provide the allocated control command almost immediately. How to incorporate actuator dynamics into the control design for spacecraft pose tracking with fault tolerance remains an interesting topic for future research. Remark 4.3: The above combination of CA and AISMC algorithms has some notable features: 1) the nominal control fˆ vc cannot be precisely restored unless the fault detection is ˜ exact (i.e., E(t) = E(t)), and 2) the nominal performance of controller (12) can only be approximated instead of being exactly recovered. These features are similar to the FTC methods in [28]–[30] but differ from the virtual actuator approach [36] for fault-tolerant control allocation, which can restore the nominal controller and the nominal performance exactly, in spite of actuator faults. The virtual actuator method in [36], however, is limited to linear systems and its extension to nonlinear systems, such as the spacecraft pose tracking system, is not trivial and needs further study. V. N UMERICAL E XAMPLES This section presents three sets of illustrative numerical simulations to demonstrate the application of the proposed methods. The first set analyzes the behavior of the nominal controller (12) as a function of the design parameters, aiming to provide a gain-tuning guide for the desirable sliding-mode dynamics behavior for the AISMC given in (16), (19), and (20). In the latter two sets, a circumnavigation operation in the


TABLE I: Nominal Orbital Elements of Target Spacecraft RAAN Inclination angle Argument of perigee at t = 0 s Eccentricity Perigee altitude

xT

yT

follower

rdt yD

xB

yB

target

zT

120 deg 60 deg 30 deg 0 1000 km

xD

rbd

zB

zI

zD

yI xI Fig. 2: Reference frames in the formation flight around Earth

Earth gravity field and a hovering operation about an asteroid are performed using the AISMC respectively.

A. Parametric Analysis of the Nominal Controller: Monitoring of an Earth-Orbiting Target This subsection investigates the performance of the nominal controller (12) as a function of control parameters via formation flight of two Earth-orbiting spacecraft, namely, a target and a follower. Figure 2 defines the inertial frame FI , the target frame FT , the desired frame FD , and the follower’s bodyfixed frame FB , where {x? , y ? , z ? }, ? ∈ {I, T , D, B}, represent the three orthogonal coordinate axes of the corresponding frame. Note that FI is set as the Earth-centered-inertial frame while FT is set as the local-vertical-local-horizontal (LVLH) frame of the target and is centered at its CM. Table I gives the the nominal orbit elements of the target under the Earth’s mass-point gravity. The motion states of the target spacecraft relative to FI , i.e., (r iti , v iti , v˙ iti , q ti , ω iti , ω˙ iti ), can be propagated by applying the classical orbital mechanics and integrating the mass-point gravitational acceleration plus the Earth’s J2 perturbation (oblateness effect). By means of (1) and (2), these quantities can then be used to derive the dual ˆ ti , dual velocity ω ˆ iti , and dual acceleration ω quaternion q ˆ˙ iti . The follower is required to achieve a stationary monitoring position relative to the target at r tdt = [−200, 0, 0] m while simultaneously aligning its attitude along FT such that xB ˆ tdt = ω always points to the target, i.e., q dt = 1 and ω ˆ˙ tdt = 06 . ˆ di = q ˆ ti ⊗ q ˆ dt and (2) and (3), the desired velocity Invoking q and acceleration are computed as [3]: ˆ ddi = q ˆ ∗dt ⊗ ω ˆ tdt ⊗ q ˆ dt + ω ˆ dti , ω

ˆ dti = q ˆ ∗di ⊗ ω ˆ iti ⊗ q ˆ di (24) ω

ˆ ∗dt ⊗ ω ˆ dt + q ˆ ∗di ⊗ ω ˆ di − ω ˆ ddi × ω ˆ dti (25) ω ˆ˙ ddi = q ˆ˙ tdt ⊗ q ˆ˙ iti ⊗ q

7

The true mass and inertia matrix of the follower are   60 −7 4 53 −1.5  kg · m2 . m = 150 kg J =  −7 4 −1.5 70 while the nominal values available for control implementation are assumed to be m0 = 160 kg and J 0 = diag{67, 55, 72} kg · m2 . The environmental force and torque acting on the follower take the following form b b b b b b fg ˆb = fd fˆ e = fˆ g + fˆ d , fˆ g = , f d τ bd τ bg where f bg and τ bg are the gravitational force and torque while f bd and τ bd are the external disturbance force and torque, all expressed in FB . Letting z I = [0, 0, 1]T and i i i r¯ ibi = r ibi /kr ibi k = [¯ rbi1 , r¯bi2 , r¯bi3 ]T , f bg and τ bg are then computed as r b ×J r b f bg = mRbi aig , τ bg = 3µE bikr b k5 bi bi i i ¯ ibi +2¯ µ rb 3µ J R2 [(1−5(¯ rbi3 )2 )2 r rbi3 zI ] aig = − krEb kbi3 − E 2 E i 4 2k r k bi bi

(26)

where µE = 398, 600.4418 km3 /s2 and J2 = 0.0010826267, and RE = 6378.137 km is the Earth’s mean equatorial radius. It is assumed that f bd and τ bd are given by f bd = 2[3 + s(ω0 t), 2 − 1.5c(ω0 t), 1 + 2s(ω0 t)]T × 10−3 N τ bd = [1 − c(ω0 t), 2 + 2s(ω0 t), 2 + 1.2s(ω0 t)]T × 10−3 N · m where ω0 = 0.004 rad/s, s(·) , sin(·) and c(·) , cos(·). Assume that the follower spacecraft is equipped with four pairs of thrusters for position control and four reaction wheels for attitude control. Their alignment matrices are set to √ √ √     1 0 0 √3/3 1 0 −√ 3/3 √3/3 D p =  0 0 √3/3 √3/3  , D a =  0 1 0 √3/3  . 3/3 3/3 0 1 0 0 1 3/3 (27) The maximum thruster force and wheel torque are set to 0.5 N and 0.15 N · m, respectively. The initial pose and velocity of the follower relative to the target are set to r tbt (0) = [−400, 0, 0]T m, q bt (0) = [0.2646, 0.8, −0.5, −0.2]T , v tbt (0) = [0.1, −0.2, −0.1]T m/s, and ω tbt (0) = [0.05, 0.05, 0.05]T rad/s. The control gains are kp1 = 0.03, kp2 = 1.2, η = 0.001, ka1 = 0.1, ka2 = 1, δ = 0.1, and h(0) = 1. Note that the simulations in this Sec˜ tion assume fully healthy actuators (i.e., E(t) = E(t) = I 8) b b and the nominal value of fˆ e is fˆ e0 = [(f g0 )T , (τ bg0 )T ]T , where f bg0 and τ bg0 are obtained by replacing m and J in b (26) with m0 and J 0 . In other words, fˆ d is assumed to be completely unknown. The influence of random measurement noise is also taken into account in the simulations. For most close proximity operations between spacecraft, the follower spacecraft can usually obtain its pose and velocity relative to the target in high precision. Hence, we assume that the position and velocity errors from the navigation module are given by r bbd,m = r bbd +r n and v bbd,m = v bbd +v n , where r n ∈ R3 and v n ∈ R3 are zero mean white Gaussian noise processes with variances (0.01 m)2 I 3


TABLE II: Steady-State Accuracy for Different α1 α1

k¯ q bd k

kω bbd k (rad/s)

kr dbd k (m)

kv dbd k (m/s)

0.6 0.7 0.8 0.9 1

3.8 × 10−3 8.3 × 10−3 1.5 × 10−2 2.4 × 10−2 3.5 × 10−2

2.2 × 10−5 4.6 × 10−5 8.3 × 10−5 1.3 × 10−4 2.2 × 10−4

0.29 0.34 0.37 0.4 0.41

4.2 × 10−4 5.1 × 10−4 5.5 × 10−4 5.7 × 10−4 2.6 × 10−3

and (0.001 m/s)2 I 3 , respectively. Similarly, the attitude and angular velocity errors from the navigation module are given by q bd,m = q bd ⊗ [cos(kφn k/2), φTn sin(kφn k/2)/kφn k]T and ω bbd,m = ω bbd + ω n , where φn ∈ R3 and ω n ∈ R3 are white Gaussian noise processes with variances (20 arcsec)2 I 3 and (1 × 10−6 rad/s)2 I 3 , respectively. Simulations are conducted first for α1 = 0.6, 0.7, 0.8, 0.9, 1 (the other parameters are fixed). Figure 3 compares the red sponses of the follower’s pose error {qbd1 , rbd1 } and velocity b d error {ωbd1 , vbd1 } while Fig. 4 compares the command force for Thruster 1 and command torque for RW 1. As α1 increases, the spacecraft position and attitude trajectories become more and more oscillatory. Besides, the transient overshoots and control input magnitudes also increase. Table II compares the ¯ bd , ω bbd , r dbd , and v dbd in terms of their norms. residual error of q Clearly, smaller α1 results in smaller steady-state pose tracking error. Hence, with other parameters fixed, increasing α1 can decrease the damping on the closed-loop pose motion, and also reduce the robustness against parametric uncertainties and unknown disturbances. When the damping effect is too weak, the transient response in the attitude and position trajectory exhibits slow convergence speed and significant oscillations. The latter then causes highly oscillatory and substantial command torque and force (Fig. 4). Next, simulations are conducted for η = 0.001, 0.005, 0.01, 0.02 and α1 = 0.6 (the other parameters are fixed). In this case, only the behavior of the position d d , and fu1 , vbd1 motion is affected. The responses of rbd1 are plotted in Fig. 5. As η increases, the magnitude of the command force decreases at the beginning phase. As a result, the maximum translation speed relative to the desired frame, occurred during transient phase, decreases and it takes longer and longer time for the follower to reach the desired position. As shown in Fig. 5b, for η = 0.2 the follower almost coasts d at a constant speed (vbd1 ≈ 0.1 m/s) along the xT direction toward the desired position during the time interval [0, 27] min. Hence, increasing η can reduce the damping on the close-loop position motion. Note that varying η only affects the transient behavior of the position trajectory and does not affect the position tracking error in steady phase. From (12), it is expected that increasing kp1 and ka1 enlarges the spring force imposed on the spacecraft while increasing kp2 and ka2 tends to lift the damping effect. These influences on closed-loop behavior are quite evident and are thus not shown numerically. B. Circumnavigation Around an Earth-Orbiting Target This section considers the same two-spacecraft formation flight system as the previous section and examines the per-

8

formance of the proposed AISMC in the presence of actuator faults and the other uncertainties. The follower is required to achieve and maintain a circular circumnavigation around the target with a radius of 100 m in the y T − z T plane and with −y B pointed continuously toward the target. More precisely, the desired position and attitude relative to FT are r tdt = [0, 100c(0.001t), 100s(0.001t)]T m and q dt = [c(0.001t/2), s(0.001t/2), 0, 0]T . It is direct to derive that ω tdt = [0.001, 0, 0]T rad/s, ω˙ tdt = 03 , v tdt = r˙ tdt and ˆ dt , ω ˆ tdt and ω v˙ tdt = r¨ tdt . After the construction of q ˆ˙ tdt , one d d ˙ ˆ di = q ˆ ti ⊗ q ˆ dt , ω ˆ di and ω can then compute q ˆ di following (24) and (25). The measurement noise remains the same as the previous section. The actuators are assumed to experience the following fault scenarios: ep1 (t) = 0.98[t − 300]− + 0.8[t − 300]+ + 0.008rand(0, 1) ep2 (t) = 0.95[t − 600]− + 0.6[t − 600]+ + 0.005rand(0, 1) ep3 (t) = 0.98 + 0.005rand(0, 1) ep4 (t) = 0.95[t − 800]− + 0.5[t − 800]+ + 0.002rand(0, 1) ea1 (t) = 0.98[t − 25]− + 0.4[t − 25]+ + 0.02s(0.01t) ea2 (t) = 0.92[t − 65]− + 0.5[t − 65]+ + 0.05s(0.02t) ea3 (t) = 0.99[t − 100]− + 0.7[t − 100]+ ea4 (t) = 0.8 − 0.03s(0.01t) where rand(x, y) generates real-time noise subject to the normal distribution with mean x ∈ R and standard deviation y ≥ 0; the functions [·]− and [·]+ are defined as 1 x µ? , ? ∈ {p, a}, it follows from (19) that ks?0 k = 1, while if ks? k ≤ µ? , it is derived that $? ks? k[1 − ks?0 k] = $? ks? k[1 − µ−1 ? ks? k] −1 = µ?0 $? ks? k[µ?0 − $? ks? k] ≤ µ4?0 Hence, (31) leads to V˙ 2 ≤ −kp $p ksp k − ka $a ksa k

A PPENDIX C P ROOF OF T HEOREM 4.1

T

p

= −kp $p ksp k − ka $a ksa k +(kp + 2ξ˜p )$p ksp k[1 − ksp0 k] +(ka + 2ξã )$a ksa k[1 − ksa0 k] −1 ¯2 ¯2 +ρ−1 p1 ρp2 ξp + ρa1 ρa2 ξa √

ˆT δˆ ξ s

The above analysis together with Lemma 4.1 amounts to V˙ 2 ≤ −kp $p ksp kksp0 k − ka $a ksa kksa0 k +2ξ˜p $p ksp k[1 − ksp0 k] + 2ξã $a ksa k[1 − ksa0 k] +ρ−1 ρp2 ξ¯2 + ρ−1 ρa2 ξ¯2

kre(δˆ S + δˆ D + δˆ e + δˆ k )k ≤ γp (t)$p kdu(δˆ S + δˆ D + δˆ e + δˆ k )k ≤ γa (t)$a

V˙ 2

?

3, we can

Summarizing the above analysis amounts to

where

˜˙ ¯−1 ˜ ˜ −ρ−1 ?1 [∆ξ? ]ξ? = −[∆ξ? ]$? ks? k + ξ? [ξ? − ξ? ]ξ? $? ks? k −1 ˜ +ρ?1 ρ?2 ξ? ∆ξ? ≤ [−ξ? + 2ξ˜? ]$? ks? k − ξ¯?−1 ξ˜?2 $? ks? k ¯2 +ρ−1 ?1 ρ?2 ξ? ≤ [−ξ? + 2ξ˜? ]$? ks? k + ρ−1 ρ?2 ξ¯2

a

It follows from (19) that sT? s?0 = ks? kks?0 k, ? ∈ {p, a}. In addition, Assumption 2.3 gives γ0 ≤ 1 − γ1 − γ2 . We can further obtain b ˆT δˆ ξ + s ˆT fˆ vc1 ≤ −[kp + 2ξ˜p ]$p ksp kksp0 k s −[ka + 2ξã ]$a ksa kksa0 k

Since the adaptive law provides 0 ≤ ξ˜? (t) ≤ ξ¯? , ? ∈ {p, a}, it follows that |∆ξ? | ≤ ξ¯? . This property and 0 < ξ? /ξ¯? < 1 can be used to show that

˜ (k +2ξ˜ )µ ρ ¯2 + p 4 p p0 + (ka +24ξa )µa0 + ρp2 ξ + ρρa2 ξ¯2 p1 p a1 a ≤ −kp ksp k − ka ksa k ¯ (k +2ξ¯ )µ ρ ¯2 + p 4 p p0 + (ka +24ξa )µa0 + ρp2 ξ + ρρa2 ξ¯2 p1 pi a1 a h |∆ξp | | √ a ≤ − η1 ksp k + ksa k + √ρp1 + |∆ξ ρa1 + η0

(32)

where $? ≥ 1 and ξ˜? ≤ ξ¯ are used when deriving the second inequality; η1 and η0 are given in (22) and (23) respectively. Applying the inequality [38] !α n n X X 0