FULL-MATRIX INVERSE-BASED MIMO QFT CONTROL DESIGN FOR SPACECRAFT FLYING IN FORMATION

FULL-MATRIX INVERSE-BASED MIMO QFT CONTROL DESIGN FOR SPACECRAFT FLYING IN FORMATION M. Garcia-Sanz 1, I. Eguinoa 1, S. Bennani 2 1

Automatic Control and Computer Science Department. Public University of Navarra. Campus Arrosadia, 31006 Pamplona, Spain 2

Guidance, Navigation and Control Section, ESA/ESTEC Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

Abstract: This paper introduces a reformulation to design full-matrix Quantitative Feedback Theory (QFT) controllers for multi-input-multi-output (MIMO) plants with model uncertainty. It considers a double-step procedure: an inverse-based decoupling and a consecutive loop-by-loop quantitative robust control design. The method generalizes several previous non-diagonal MIMO QFT techniques, avoiding some required prior hypotheses of such former methods, and simplifies the design procedure. It deals with two classical control problems: reference tracking and disturbance rejection at plant output. The paper ends applying the new technique to the design of a MIMO controller for a spacecraft flying in a formation that is moving with respect to a central body in a circular orbit. Copyright © 2007 IFAC Keywords: Robust Multivariable Control, Formation Flying, Satellite control, QFT.

1. INTRODUCTION Multiple-Input-Multiple-Output (MIMO) plants with large channel-interaction [large elements in the Relative Gain Array –RGA- (Bristol, 1966)] involve serious difficulties to be controlled (Skogestad and Morari, 1987). Model uncertainty makes the problem much more difficult. In the last few years a significant amount of work in robust control for MIMO uncertain systems has been done (e.g. Maciejowski, 1989; Skogestad and Postlethwaite, 2005; etc.). The Quantitative Feedback Theory (QFT) is among the most promising techniques for the mentioned problem (Horowitz, 1991; Houpis, Rasmussen and Garcia-Sanz, 2006). It is a frequency domain technique that deals with both magnitude and phase, and makes a quantitative synthesis, taking into account quantitative bounds on the plant uncertainty and quantitative tolerances on the acceptable closedloop system response.

Recently, some new methods for non-diagonal (fullmatrix) MIMO QFT robust control system design have been introduced (Yaniv, 1995; Franchek and Nwokah, 1995; Franchek, Herman and Nwokah, 1997; Boje and Nwokah, 1997; Boje and Nwokah, 1999; De Bedout and Franchek, 2002; Lan, Kerr and Jayasuriya, 2004; Kerr and Jayasuriya, 2006). Similarly, Garcia-Sanz et al. introduced a methodology to extend the classical QFT diagonal controller design for MIMO plants with uncertainty (Horowitz, 1982) to a fully populated matrix controller design, studying three classical cases: the reference tracking and the external disturbance rejection at plant input and output (Garcia-Sanz and Egaña, 2002; Garcia-Sanz, Egaña, Barreras, 2005). This method has shown good results in many scenarios (Houpis, Rasmussen and Garcia-Sanz, 2006), even improving some well-known robust control techniques. However, it requires the prior fulfilment of a hypothesis and, in some cases, an advanced expertise could be needed during the design procedure.

In this context, the present paper improves and generalizes that method by introducing a new fullmatrix MIMO QFT methodology. It handles the design for robust stability, reference tracking and disturbance rejection at plant output cases, avoiding the required prior hypothesis and simplifying the design procedure.

do’ Pdo(s)

do r’

F(s)

r

u

Gα(s)

Gβ(s)

P(s)

y

TY/R (s)

Section 2 formulates the theory of the new method. Section 3 illustrates the application of the technique by designing a MIMO QFT controller for a spacecraft flying in formation with respect to a central body in circular orbit. The results are compared to those obtained with the previous nondiagonal MIMO QFT technique. The last Section summarises the most relevant ideas of the paper. 2. FULL-MATRIX MIMO QFT CONTROLLER DESIGN METHODOLOGY. REFORMULATION

2.2 Procedure Consider a nxn linear multivariable system (Fig. 1), composed of a plant P, a fully populated matrix controller G = Gα Gβ, a prefilter F, and a plant output disturbance transfer function Pdo, where P ∈ ℑP , ℑP is the set of possible plants due to uncertainty, and, P = [ pij ] ; F = [ f ij ] ; Gα = [ g ijα ] ; Gβ = diag ( g iiβ ) (1)

The reference vector r’ and the external disturbance vector at plant output do’ are the inputs of the system. The output vector y is the variable to be controlled. The plant inverse P∗ is,

The proposed controller design includes the following steps:

)

RGA ( jω ) = P ( jω ) ⊗ P −1 ( jω )

The first objective of the non-diagonal elements of the full-matrix controller is to minimize the loop interactions (coupling). They act as feedforward functions, cancelling some of the dynamics of the non-diagonal elements of the plant matrix. As a consequence, they are sensitive to the model uncertainty. By contrast, the diagonal elements of the controller act as feedback functions, minimizing the effect of the uncertainty (sensitivity) and fulfilling the performance and stability specifications. They also have to overcome the uncertainty of the nondiagonal plant elements, minimizing the residual coupling that the non-diagonal elements of the controller have not cancelled. Taking this as a starting point, the new full-matrix controller design methodology adopts the following procedure.

P -1 = P * = [ pij* ]

Step A. Controller Structure: By using the RGA (Bristol, 1966) over the frequencies of interest (3) (Skogestad and Postlethwaite, 2005), and taking into account the requirement of minimum complexity for the controller, the method identifies the most appropriate structure for the compensator: the required off-diagonal and diagonal elements, which correspond with the significant values of the RGA matrix.

(

2.1 Overview

G = [ g ij ] = Gα Gβ

Fig. 1. Two Degree of Freedom MIMO System

(2) methodology

T

(3)

where ⊗ denotes element-by-element multiplication (Schur product). Step B. Design of Gα: The fully populated matrix controller G is composed of two matrices: G = Gα Gβ. The main objective of the pre-compensator Gα is to diagonalize the plant P as much as possible. That is to say,  pˆ jj ∆ˆ ji   Gα = g ijα = Pˆ −1 Pˆdiag =   ∆ˆ 

[ ]

(4)

where Pˆ is a plant matrix selected within the uncertainty, ∆ˆ is the determinant of the Pˆ matrix and ∆ˆ ji the jith cofactor of the Pˆ matrix. The plant matrix Pˆ is selected so that the expression of the extended matrix Px = P Gα = [pijx] presents the closest form to a diagonal matrix, nulling the offdiagonal terms as much as posible by minimizing the cost function J, where, J = ∑ µij pijx , with (0 ≤ µ ij ≤ 1)

(5)

i≠ j

and where, provided the G controller does not introduce non-minimum phase transmission zeros (see section below), the pijx expression is, pijx =

pˆ jj ∆ˆ

(p

ˆ + p ∆ˆ + L + p ∆ˆ i 2 j2 in jn

i1 ∆ j1

)

(6)

Step C. Design of Gβ: After determining Gα, the method goes on with the design of a diagonal matrix Gβ that will fulfil the desired robust stability and

robust performance specifications for the extended plant Px = P Gα. The Px* matrix is reorganized so that (p11x*)−1 has the smallest bandwidth, (p22x*)−1 the next smallest bandwidth, and so on. The compensators gkkβ(k = 1 to n) are calculated by using a sequential (loop-by-loop) standard QFT loopshaping methodology (Houpis, Rasmussen and Garcia-Sanz, 2006) for the inverse of the equivalent plant [piix*e]i−1. It satisfies the recursive relationship of (7) – (Franchek, Herman and Nwokah, 1997) with gij = 0, i ≠ j-, and takes into account the compensator elements previously designed.  p x * e =  p x * e −  ij  k  ij  k −1

[

i, j ≥ k; P x*

e

]

k =1

[p ] [p ] [p ] + [g ]

= P x*

x* e i (k -1) k −1

x*e (k -1) (k -1) k −1

x* e (k -1) j k −1

β

(k −1) (k −1) k −1

(7)

The presence of model uncertainty (P ∈ ℑP) reduces the diagonalization effect of the pre-compensator Gα over P. That diagonalization is a real cancellation only when the plant P is exactly at Pˆ (working point): P Gα = Pˆ Pˆ −1 Pˆdiag = Pˆdiag . If the plant P is

working on a different point, then the extended plant that the compensator Gβ sees will include offdiagonal elements as well. Consequently, the performance specifications used to design the elements gkkβ will be modified in order to avoid that residual coupling, according to the classical methodology that takes into account the coupling effect loop by loop (Houpis, Rasmussen and GarciaSanz, 2006, pp. 177-181). Stability conditions: The closed-loop stability of the MIMO system with the full matrix controller G = [gij] is guaranteed by the following sufficient conditions (De Bedout and Franchek, 2002): i. each Li(s) = giiβ(s) [piix*e(s)]i−1, i = 1, ..., n, satisfies the Nyquist encirclement condition, ii. no RHP pole-zero cancellations occur between giiβ(s) and [piix*e(s)]i−1, i = 1,...,n, iii. no Smith-McMillan pole-zero cancellations occur between P(s) and G(s), and iv. no Smith-McMillan pole-zero cancellations occur in P*(s) + G(s).

Conditions i and ii are checked loop by loop when the compensators giiβ(s) are calculated in Step C. Conditions iii and iv are checked after Step C is finished. Avoiding non-minimum phase transmission zeros due to the controller: Although remote, there exists a theoretical possibility of introducing right-half plane (RHP) transmission zeros due to the controller design. This undesirable situation cannot be detected until the multivariable system design is completed. To avoid it, the proposed methodology –Steps A, B and C– is inserted in a more general procedure

(Garcia-Sanz and Eguinoa, 2005). Once the matrix compensator G(s) is designed, the transmission zeros of P(s) G(s) are determined by using the SmithMcMillan form and over the set of possible plants ℑP due to uncertainty. If there exist new RHP zeros apart from those initially present in P(s), they can be removed by modifying conveniently some nondiagonal elements of the matrix G(s). That involves an additional role for these controller elements. Case 2x2

By using [Eqs. (1), (2), (4)], the expression of the final controller G(s) for a 2x2 MIMO plant is: g12  g G = Gα Gβ =  11 =  g 21 g 22  − pˆ12 pˆ11 pˆ 22   β β g 22  pˆ pˆ − pˆ pˆ g11  ˆ ˆ p p pˆ11 − 12 21  11 22 12 21  pˆ 22  =  − pˆ 21 pˆ11 pˆ 22 β  β g11 g  pˆ pˆ pˆ11 pˆ 22 − pˆ12 pˆ 21 22   pˆ 22 − 12 21 pˆ11  

(8)

For tracking or disturbance rejection at plant output, the expression of the previous non-diagonal MIMO QFT methodology described by (Garcia-Sanz, Egaña and Barreras, 2005) for 2x2 MIMO plants was: g G =  11  g 21

 g11 g12    = g 22   − pˆ 21 g  pˆ 22 11

− pˆ 12  g 22  pˆ 11  g 22  

(9)

It is clear that the new approach (8) is a generalization of the previous one (9). Equation (8) tends to equation (9) if ( pˆ11 pˆ 22 >> pˆ12 pˆ 21 ) , which is just the diagonal dominance hypothesis of the previous method. Remarks

There are well-known difficulties (Skogestad and Morari, 1987) that the inversed-based controllers could introduce under model uncertainty, which yields for the new formulation P Gα ≠ Pdiag. However, that problem is solved by means of the Gβ controller, which deals with the model uncertainty within the QFT framework. 3. SPACECRAFT FLYING IN FORMATION IN A CIRCULAR ORBIT To illustrate its benefits, the MIMO QFT reformulation presented above is applied to control a spacecraft (denoted as s/c from now on) flying in a formation that is moving with respect to a central body in a circular orbit.

3.1 Model

The non-linear dynamic equation of motion of a spacecraft flying in formation, relative to the centre of mass of the formation (c.m.f), following a closed Keplerian reference orbit, is shown in Eq. (10) –Fig. 2- (Wiltshire and Clohessy, 1960; Ploen, Scharf, Hadaegh et al., 2004; Garcia-Sanz and Hadaegh, 2007). m &ρ& + 2 m [ω] ρ& + m ([ω]2 + [ω& ]) ρ + + m (b& + [ω] b) = Q + D − U

(10)

q

where the reference orbit defines an orbit reference frame (x,y,z) –Fig. 2-, which serves as the primary frame to analyze the dynamic of the s/c formation. The unit vector x(o1) points anti-nadir, z(o3) points in the direction of the orbit normal, and y(o2) completes the right-handed triad. The inertial frame of reference (n1,n2,n3) –Fig. 2-, is attached to the centre of the Earth. The unit vector n1 points toward the vernal equinox, n3 points toward the geographic North Pole, and n2 completes the right-handed triad. m is the mass of the s/c; Q = [Qx , Qy , Qz]T are the actuator forces (thrusters); D = [Dx , Dy , Dz]T are the disturbance forces acting on the satellite due to lunisolar attraction, solar radiation pressure, Earth’s oblatness, tidal effects, Earth’s geomagnetic field, ∂ − µ m etc.; U q =   is the gradient of the Earth ∂q  R  potential field; [ω]=[0 -ω0 0 ; ω0 0 0 ; 0 0 0]; b = [ R& 0 , R0 ω 0 , 0]T , R = [( x + R0 ) 2 + y 2 + z 2 ]1 / 2 . The system outputs to be controlled are the distances r r r r ρ = [ x, y, z ]T ; ρ = x O1 + y O2 + z O3 of the s/c to the c.m.f (orbit reference frame). R0 is the mean orbit radius, ω0 = µ / R03 is the mean orbit rate and µ = 3.9860x105 Km3/sec2 denotes the gravitational parameter of the Earth. Both parameters (R0 and ω0) present some uncertainty due to the disturbances. Alternatively, the reference orbit can also be described by the orbital elements a (semi-major axis), e (eccentricity), i1 (inclination), Ω1 (longitude of the ascending node), ω1 (argument of perigee), and T (time of perigee passage) –Fig. 2-. For circular orbits, the orbit radius R0 can present a slow time-variation due to uncertainty, which && ≈ 0 . Applying these conditions means R& ≈ 0 and R 0

0

in Eq. (10), the relative equations of motion of a s/c relative to the c.m.f in circular orbit are, &x& − 2 ω0 y& − 3 ω02 x =

&y& + 2 ω0 x& = &z& + ω02 z =

1 (Qx + Dx ) m

1 (Q y + D y ) m

1 (Qz + D z ) m

(11) (12) (13)

which are the so-called Clohessy-Wiltshire-Hill (CWH) equations (Wiltshire and Clohessy, 1960). Applying the Laplace Transform on these equations, the transfer function matrix P relating the actuator forces and the position of the s/c becomes, 1   2 + s ω02   x  y  = 1  − 2 ω0   m  s (s 2 + ω 2 ) 0   z   0  

2 ω0

    0   1  s 2 + ω02  0

s ( s 2 + ω02 ) s 2 − 3 ω02 s 2 ( s 2 + ω02 ) 0

Qx  Q   y Qz 

(14) which is a 2x2 MIMO system (axes x and y) plus a SISO (Single-Input-Single-Output) system (axis z). Orbit Normal

North Pole

n3 ω

z (o3)

y (o2)

Periapsis direction AntiNadir

i1

c.m.f

x (o1)

ρ ω1

n2

Earth

n1

R0

s/c

Ω1

Vernal equinox

Fig. 2. Spacecraft general orbital geometry. 3.2 Parameter values Equation (14) is particularised for a s/c flying in formation in a circular geostationary Earth orbit. Since the sidereal day is 23 h 56 min 4 sec, then ω0 = 2 π / 86164 sec = 7.29212 x 10-5 rad/sec. It corresponds to an orbit radius R0 = 42164 km. This radius presents some uncertainty because of external disturbances: R0 ∈ [41953 , 42375] km . So, the mean orbit rate is ω0 ∈ [7.2378, 7.3472]x10-5 rad/sec. Due to fuel consumption, the mass of the satellite also presents some uncertainty: m ∈ [1350, 1650] kg. 3.3 Objectives and Specifications The notation adopted in the next sections for transfer function expressions denotes the steady state gain as a constant without parenthesis; simple poles and zeros as (ω), which corresponds to (s/ω + 1); poles and zeros at the origin as (0); conjugate poles and zeros as [ζ ; ωn], with ((s/ωn)2 + (2 ζ/ωn) s + 1); nmultiplicity of poles and zeros as an exponent ( )n. The main objective of the controller is to keep the satellite stable, rejecting disturbances while

controlling the distance ρ between the s/c and the c.m.f. The main specifications are:

g11 =

Robust Stability

g21 =

[

g iiβ ( s ) piix*e ( s )

[

]

−1 i

1 + g iiβ ( s ) piix*e ( s )

]

−1 i

Reference Tracking

[

g iiβ ( s) piix*e ( s)

[

]

−1 i

1 + g iiβ ( s) piix*e ( s)

]

−1 i

≤ TRU (ω )

TRL (jω ) = (a) [ζ ; ω n ]

TRU (jω ) = 1 {(σ 1 ) (σ 2 ) (σ 3 )}

−4 2

(1) (4.5)(1.2633×10 )

; g22 =

5×10−4 (0.012)(1.1×10−4 ) (0) (1) (4.5)

(21) (15)

TRL (ω ) ≤

7.3122(0.014)(10−4 )

5×10−4 (0.012)(1.1×10−4 ) g33 = ; g13 = g23 = g31 = g32 = 0 (0)(1) (4.5)

∀ω , i = 1,2,3

≤ 1.05

8×10−4 (0.014)(10−4 ) − 7.2937×10−8 (0.012)(1.1×10−4 ) ; g12 = (0) (1) (4.5) (0)2 (1) (4.5)

(16)

(17)

where a = 0.7 ; ζ = 0.65 ; ωn = 0.3 ; σ1 = 2 ωn ; σ2 = 0.3 ωn ; σ3 = ωn , for ω = [5x10-6 , 1.0] rad/sec, i=1,2,3.

 (0.42) (0.048)   (0.17) (0.038) F = 0    0 

    0  (0.42) (0.048)   (0.22) (0.038) 

0

0

(0.42) (0.048) (0.17) (0.038) 0

(22)

Fig. 3 shows how the diagonal dominance hypothesis of the previous method is fulfilled for the main part of the frequencies of interest. Then Eq. (8) can be simplified to Eq. (9) and both methods (the previous one and the new one) are equivalent in this particular example, as can be seen in the following section. 400 350

3.4 Controller Design The controller Gα is designed according to Step B in Section 2 and Eqs. (4) and (14). Due to the specific characteristics of the plant –see Eq. (14)-, the Gα literally derived from Eq. (4) presents a right half plane (RHP) transmission zero at s = ωˆ 0 3 rad/sec and a pair of imaginary complex poles at s = ± j ωˆ 0 rad/sec, according to the Smith-McMillan form. The Gα finally designed avoids these elements by changing the RHP transmission zero for its corresponding all-pass filter and damping the imaginary poles. It yields,

( )

3 3α 2   [ζ ;α ]  Gα =  2 α −1 (0)   [ζ ;α ]   0 

− 6α

( 3α )

(0) [ζ ;α ]

( )

3 3α

[ζ ;α ] 0

2

2

 0    0   0 

(18)

α = ωˆ 0 = 7.2937 × 10 −5 rad/sec ; ζ = 0.8 The controller Gβ and the prefilter F are: Gβ =

5 ×10−4 (0.012) (1.1×10−4 ) I 3×3 (0) (1) (4.5)

(19)

F=

(0.42) (0.048) I 3×3 (0.22) (0.038)

(20)

The controller and the pre-filter designed with the previous full-matrix MIMO QFT method (GarciaSanz and Egaña, 2002; Garcia-Sanz, Egaña and Barreras, 2005) are,

Magnitude (dB)

300 250

| p11 p22 |

200 150 100 50

| p12 p21 |

0 10

-6

10

-5

-4

10

10

-3

10

-2

10

-1

Frequency (rad/sec)

Fig. 3. Bode Analysis of ( p11 p22 >> p12 p21 ) 3.5 Results Fig. 4 shows the results for the new full-matrix MIMO QFT controller design [Eqs. (18)-(20)], being fully equivalent to that obtained with the previous method [Eqs. (21)-(22)]. The figure presents the three outputs [x, y, z] for a reference tracking problem and for the plants defined by the contour of the parameter uncertainty space (ω0max, ω0min, mmax, mmin). These upper and lower outputs delimit the rest of the plant responses within the uncertainty. The coupling among loops (see t = 600 sec and t = 1100 sec in Fig. 4) is minimized and the reference tracking specifications are fulfilled. Besides, other classical multivariable specifications (Skogestad and Postlethwaite, 2005) can be checked to assess the stability and performance of the designed control system. The maxω[σmax(T)] is lower than 0.42 dB within the uncertainty. With respect to other matrices of interest, the maximum values within the uncertainty and over the frequency range are: maxω[σmax(SO)] = 1.90 dB and maxω[σmax(SI)] = 1.90 dB, where SO = (I + P G)-1 and SI = (I + G P)-1.

Christian Philippe for his valuable comments and contributions for the success of this work.

1.015 m =1650 kg , ω

,ω

m =1350 kg , ω

,ω

max

1.01

max

min min

REFERENCES

Position x [m]

1.005 1 0.995 0.99 0.985 0.98

0

200

400

600

800 1000 Time [sec]

1200

1400

1600

1.015 m =1650 kg , ωmax , ωmin 1.01

m =1350 kg , ωmax , ωmin

Position y [m]

1.005 1 0.995 0.99 0.985 0.98

0

200

400

600

800 1000 Time [sec]

1200

1400

1600

1.015 m =1650 kg , ωmax , ωmin

1.01

m =1350 kg , ω

max

,ω

min

Position z [m]

1.005 1 0.995 0.99 0.985 0.98 1500

1600

1700

1800 Time [sec]

1900

2000

2100

Fig. 4. Time domain results. 4. CONCLUSIONS This paper has presented a new methodology to design full-matrix QFT controllers for MIMO plants with model uncertainty. It generalizes a previous non-diagonal MIMO QFT method, solving simultaneously two classical problems: reference tracking and disturbance rejection at plant output. The new formulation does not need a previous hypothesis of diagonal dominance and simplifies the calculations for the off-diagonal elements, finding the solution directly. The paper illustrated the new technique with the design of a MIMO QFT controller for a spacecraft flying in a formation that is moving with respect to a central body in a circular orbit. ACKNOWLEDGMENT This work has been performed under ESA/ESTEC contract 19326/05/NL/LvH. Authors gratefully appreciate the support given by the Spanish Government (MEC) under grant CICYT DPI’200615522-C02-01. Acknowledgment is also given to Dr.

Boje, E. and O. Nwokah (1997). Quantitative Feedback Design Using Forward Path Decoupling. Symposium on Quantitative Feedback Theory and Other Frequency Domain Methods, University of Strathclyde (United Kingdom). Boje, E. and O. Nwokah (1999). Quantitative Multivariable Feedback Design for a Turbofan Engine with Forward Path Decoupling. Int. J. Robust Nonlinear Control 9(12): 857-882. Bristol, E. H. (1966). On a new measure of interaction for multivariable process control. IEEE Transactions on Automatic Control AC-11(1): 133-134. De Bedout, J. M. and M. A. Franchek (2002). Stability conditions for the sequential design of non-diagonal multivariable feedback controllers. Int. Journal of Control 75(12): 910-922. Franchek, M. A., P. Herman and O. D. I. Nwokah (1997). Robust nondiagonal controller design for uncertain multivariable regulating systems. Transactions of the ASME. Journal of Dynamic Systems, Measurement and Control 119(1): 80-85. Franchek, M. A. and O. Nwokah (1995). Robust Multivariable Control of Distillation Columns Using Non-Diagonal Controller Matrix. IMECE, ASME Dynamics Systems and Control Division 57(1): 257-264. Garcia-Sanz, M., I. Egana and M. Barreras (2005). Design of Quantitative Feedback Theory non-diagonal controllers for use in uncertain multiple-input multiple-output systems. IEE Proceedings-Control Theory and Applications 152(2): 177-187. Garcia-Sanz, M. and I. Egaña (2002). Quantitative non-diagonal controller design for multivariable systems with uncertainty. Int. Journal of Robust and Nonlinear Control 12(4): 321-333. Garcia-Sanz, M. and I. Eguinoa (2005). Improved non-diagonal MIMO QFT design technique considering non-minimum phase aspects. 7th Int. Symp. on QFT and Robust Frequency Domain Methods, Lawrence, Kansas, USA. Garcia-Sanz, M. and F. Y. Hadaegh (2007). Load-Sharing Robust Control of Spacecraft Formations: Deep Space and Low Earth Elliptic Orbits. IET Control Theory & Applications. 1(2): 475484. Horowitz, I. (1982). Improved design technique for uncertain multiple-input-multiple-output feedback systems. International Journal of Control 36(6): 977-988. Horowitz, I. (1991). Survey of quantitative feedback theory (QFT). International Journal of Control 53(2): 255-291. Houpis, C. H., S. J. Rasmussen and M. Garcia-Sanz (2006). Quantitative Feedback Theory Fundamentals and Applications. Florida, USA, a CRC Press Book. Taylor and Francis. Kerr, M. L. and S. Jayasuriya (2006). An improved non-sequential multi-input multi-output quantitative feedback theory design methodology. Int. J. Robust Nonlinear Control 16(8): 379-95. Lan, C., M. L. Kerr and S. Jayasuriya (2004). Synthesis of Controllers for Non-minimum Phase and Unstable Systems Using Non-sequential MIMO Quantitative Feedback Theory. American Control Conference, Boston, MA. Maciejowski, J. M. (1989). Multivariable Feedback Design, Addison-Wesley Publishing Company. Ploen, S. R., D. P. Scharf, F. Y. Hadaegh and A. B. Acikmese (2004). Dynamics of Earth orbiting formations. Collection of Technical Papers - AIAA Guidance, Navigation, and Control Conference, Aug 16-19 2004, Providence, RI, United States. Skogestad, S. and M. Morari (1987). Implications of large RGA elements on control performance. Industrial & Engineering Chemistry Research 26(11): 2323-2330. Skogestad, S. and I. Postlethwaite (2005). Multivariable feedback control. Analysis and design. Chichester, West Sussex, England, John Wiley & Sons Ltd. Wiltshire, R. S. and W. H. Clohessy (1960). Terminal guidance for rendezvous in space. J. Aerospace Science 27(9): 653-658. Yaniv, O. (1995). MIMO QFT using non-diagonal controllers. International Journal of Control 61(1): 245-253.