2. RELATIVE MOTION DYNAMICS

In the primary satellite's Local Vertical Local Horizontal (LVLH) frame, the position vector r _o of the secondary orbiting body is given by the following,

(1) $${\bf r}_o = {\bf r}+\vec {\rho }=\lpar {r+x} \rpar i+yi+zk, \quad {\bf r}_o \cdot {\bf r}_o =\lpar {r+x} \rpar ^2+y^2+z^2,$$

where r is the position vector of the primary satellite, r is the magnitude of r and $\vec {\rho }=xi+yi+zk$ is the position vector of the secondary orbiting body relative to the primary satellite. The angular velocity and acceleration of the LVLH frame, which is moving with the primary satellite, are given by $\omega =\dot{f}k=\lpar h /h {r^{2}} \rpar \!k$ and $\alpha =\ddot {f}k=-\lpar {2\dot {r}}/ r \rpar \omega =-\lpar 2 \dot{r} / r \rpar \lpar {h / {r^{2}}} \rpar \!k$ respectively, where f is the true anomaly and $h=\vert {\bf r} \times \dot{\bf r} \vert $ is the magnitude of the orbital angular momentum. Using Kepler's two body equations of motion for the primary satellite and for the secondary orbiting body both orbiting the Earth, $\ddot{r}=-{\mu_{e}}/{r^{2}}$ and $\ddot{\bf r}_{o} =-\mu_{e} {\bf r}_{o} / {r_{o}^{3}}$ respectively, the vector equations of relative motion dynamics defining the temporal evolution of $\vec{\rho}$ , under the influence of net perturbation acceleration vector ${\bf{a}}_{p} $ are,

(2) $$\displaystyle{{d^2\vec \rho } \over {dt^2}} = - \mu _e\left( {\displaystyle{{{{\rm (}{\bi r} + \vec \rho {\rm )}}^3} \over {\mid{\bi r} + \vec \rho \mid^3}} - \displaystyle{{\bi r} \over {r^3}}} \right) + {\bi a}_p,$$

where ${\bf a}_{p} = {\bf a}_{po} -{\bf{a}}_{ps} $ , a _ps and ${\bf{a}}_{po} $ are the perturbation accelerations due to the third body effects and Earth oblateness acting on the primary satellite and on the secondary orbiting body, which can be obtained from Vepa (Reference Vepa and Zhahir2017). From kinematics, the equation of motion for the secondary orbiting body in the primary satellite's LVLH frame is given by the following:

(3) $$\ddot{\bf r}_o = \ddot{\bf r} + \ddot{\rho}+2\omega \times \dot {\rho } + \omega \times \lpar {\omega \times \rho} \rpar +\alpha \times \rho.$$

Thus, in component form we have,

(4a) $$\ddot{x}-2\dot {f}\dot {y}-\dot {f}^2x-\ddot{f}y=a_{px} +\lpar {{\mu_e }/{r^2}} \rpar -{\mu _e \lpar {r+x} \rpar }/{\vert { {\bf r}+\vec {\rho }} \vert ^3},$$

(4b) $$\ddot {y}+2\dot {x}\dot {f}-\dot {f}^2y+\ddot {f}x=a_{py} -{\mu _e y}/{\vert { {\bf r}+\vec{\rho }} \vert ^3},$$

(4c) $$\ddot {z}=a_{pz} -{\mu _e z}/{\vert { {\bf r}+\vec {\rho }} \vert ^3}.$$

In matrix notion, after including the linear terms of the right-hand sides, in the left-hand side of the equation,

(5) $$\matrix{ & {\left[ {\matrix{ {\ddot x} \cr {\ddot y} \cr {\ddot z} \cr } } \right] + \dot f\left[ {\matrix{ 0 & { - 2} & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \cr } } \right]\left[ {\matrix{ {\dot x} \cr {\dot y} \cr {\dot z} \cr } } \right] + \left[ {\matrix{ { - {\dot f}^2 - 2{\rm (}\mu _e/r^3{\rm )}} & { - \ddot f} & 0 \cr {\ddot f} & { - {\dot f}^2 + {\rm (}\mu _e/r^3{\rm )}} & 0 \cr 0 & 0 & {{\rm (}\mu _e/r^3{\rm )}} \cr } } \right]\left[ {\matrix{ x \cr y \cr z \cr } } \right]{\rm }} \cr & {\quad = \left[ {\matrix{ {a_{px} - {\rm (}\mu _e){\rm }r + x{\rm (}/\mid{\bf r} + \vec \rho \mid^3 - {\rm (}\mu _e/r^2{\rm ) (}1 - 2x/r{\rm ))}} \cr {a_{py} - {\rm (}\mu _ey/\mid{\bf r} + \vec \rho \mid^3{\rm )} + {\rm (}\mu _e/r^2{\rm ) (}y/r{\rm )}} \cr {a_{pz} - {\rm (}\mu _ez/\mid{\bf r} + \vec \rho \mid^3{\rm )} + {\rm (}\mu _e/r^2{\rm ) (}z/r{\rm )}} \cr } } \right]} \cr } $$

Ignoring the right-hand sides' results in the equations which are called the Linearized Equations of Relative Motion (LERM). Equations (5) are the Nonlinear Equations of Relative Motion (NERM), which can be expressed as,

(6) $$\matrix{ & {\left[ {\matrix{ {\ddot x} \cr {\ddot y} \cr {\ddot z} \cr } } \right] + \dot f\left[ {\matrix{ 0 & { - 2} & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \cr } } \right]\left[ {\matrix{ {\dot x} \cr {\dot y} \cr {\dot z} \cr } } \right] + \left[ {\matrix{ { - {\dot f}^2 - 2{\rm (}\mu _e/r^3{\rm )}} & { - \ddot f} & 0 \cr {\ddot f} & { - {\dot f}^2 + {\rm (}\mu _e/r^3{\rm )}} & 0 \cr 0 & 0 & {{\rm (}\mu _e/r^3{\rm )}} \cr } } \right]\left[ {\matrix{ x \cr y \cr z \cr } } \right]{\rm }} \cr & {\quad = \left[ {\matrix{ {a_{px} + \displaystyle{{\mu _er} \over {r^3}}\left( {1 - 2\displaystyle{x \over r} - \left( {1 + \displaystyle{x \over r}{\rm }} \right){\rm }{\left( {1 + 2\displaystyle{x \over r} + \displaystyle{{x^2 + y^2 + z^2} \over {r^2}}} \right)}^{ - 3/2}} \right)} \cr {a_{py} + \displaystyle{{\mu _ey} \over {r^3}}\left( {1 - {\left( {1 + 2\displaystyle{x \over r} + \displaystyle{{x^2 + y^2 + z^2} \over {r^2}}} \right)}^{ - 3/2}} \right)} \cr {a_{pz} + \displaystyle{{\mu _ez} \over {r^3}}\left( {1 - {\left( {1 + 2\displaystyle{x \over r} + \displaystyle{{x^2 + y^2 + z^2} \over {r^2}}} \right)}^{ - 3/2}} \right)} \cr } } \right].} \cr } $$

To obtain the Hill Clohessy Wiltshire (HCW) equations, we let $\dot{f}-n=\ddot {f}=0$ , $\dot {f}^{2}={\mu_{e} / r^{3}}=n^{2}$ , where n represents the mean rate of motion of the primary satellite or the mean motion. The system can be described by the state vector $ {\bf x} = \lsqb x\ y\ z\ \dot{x}\ \dot{y}\ \dot{z} \rsqb ^{T}$ . In state-space form, the LERM are given by the equation, $\dot{\bf x} = {\bf A}\lpar t\rpar {\bf x}$ , where

(7) $${\bf A}(t) = \left[ {\matrix{ 0 & 0 & 0 & 1 & 0 & 0 \cr 0 & 0 & 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 0 & 0 & 1 \cr {{\dot f}^2 + 2\mu _e/r^3} & {\ddot f} & 0 & 0 & {2\dot f} & 0 \cr { - \ddot f} & {{\dot f}^2 - \mu _e/r^3} & 0 & { - 2\dot f} & 0 & 0 \cr 0 & 0 & { - \mu _e/r^3} & 0 & 0 & 0 \cr } } \right].$$

The matrix ${\bf{A}}\lpar t\rpar $ is time varying as r, $\dot {f}$ , $\ddot f$ vary with time. Thus ${\bf{A}}\lpar t\rpar $ is periodic with the period equal to T = 2π/n. The NERM may be expressed compactly as,

(8a) $$\displaystyle{{d^2} \over {dt^2}}\left[ {\matrix{ x \cr y \cr z \cr } } \right] = - \dot f^2{\bf K}_{lerm}\left[ {\matrix{ x \cr y \cr z \cr } } \right] - \dot f{\bf G}_{lerm}\displaystyle{d \over {dt}}\left[ {\matrix{ x \cr y \cr z \cr } } \right] + \left[ {\matrix{ {{\tilde a}_{nx}} \cr {{\tilde a}_{ny}} \cr {{\tilde a}_{nz}} \cr } } \right],$$

where,

(8b) $${\bf K}_{lerm} = \displaystyle{1 \over {{\dot f}^2}}\left[ {\matrix{ { - {\dot f}^2 - 2\displaystyle{{\mu _e} \over {r^3}}} & { - \ddot f} & 0 \cr {\ddot f} & { - {\dot f}^2 + \displaystyle{{\mu _e} \over {r^3}}} & 0 \cr 0 & 0 & {\displaystyle{{\mu _e} \over {r^3}}} \cr } } \right],\quad {\bf G}_{lerm} = \left[ {\matrix{ 0 & { - 2} & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \cr } } \right]$$

and

(8c) $$\left[ {\matrix{ {{\tilde a}_{nx}} \cr {{\tilde a}_{ny}} \cr {{\tilde a}_{nz}} \cr } } \right] = \left[ {\matrix{ {a_{px} - \left( {\mu _e{\rm (}r + x{\rm )}/\mid{\bf r} + \vec \rho \mid^3 - {\rm (}\mu _e/r^2{\rm ) (}1 - 2x/r{\rm )}} \right){\rm }} \cr {a_{py} - {\rm (}\mu _ey/\mid{\bf r} + \vec \rho \mid^3{\rm )} + {\rm (}\mu _e/r^2{\rm ) (}y/r{\rm )}} \cr {a_{pz} - {\rm (}\mu _ez/\mid{\bf r} + \vec \rho \mid^3{\rm )} + {\rm (}\mu _e/r^2{\rm ) (}z/r{\rm )}} \cr } } \right].$$

3. DERIVATION OF THE NONLINEAR CONTROLLED TH EQUATIONS

The TH equations are obtained by using the true anomaly f as the independent variable. Moreover, Tschauner and Hempel (Reference Tschauner and Hempel1965) also apply a coordinate scaling prior to transforming the equations. The scaled state variable $\tilde x$ is related to the original state by the relation, $\tilde {x}=\lpar {1+e\cos\, f} \rpar x$ . Further,

(9) $$d/df = d/\dot fdt = {\rm (}r^2/h{\rm )}d/dt,\quad r = p/{\rm (}1 + e\cos {\mkern 1mu} f{\rm )},$$

where p is the semi-latus rectum and $h=\sqrt {\mu p} $ . It follows that,

(10) $$ \displaystyle{{d\tilde {x}}\over{df}}=-ex\sin\, f+\lpar {1+e\cos\, f} \rpar \displaystyle{{dx}\over{df}}, \quad \displaystyle{{dx}\over{df}}=\displaystyle{{r^2}\over{h}}\displaystyle{{dx}\over{dt}}=\displaystyle{{p^2}\over{h\lpar {1+e\cos\, f} \rpar ^2}}\displaystyle{{dx}\over{dt}}.$$

Thus,

(11) $$\displaystyle{{d\tilde {x}}\over{df}}=-ex\sin\, f+\displaystyle{{p^2}\over{h\lpar {1+e\cos\, f} \rpar }}\displaystyle{{dx}\over{dt}},$$

It follows that,

(12) $$\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = {\rm (}1 + e\cos {\mkern 1mu} f{\rm )}\left[ {\matrix{ x \cr y \cr z \cr } } \right],\quad \displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = - e\sin {\mkern 1mu} f\left[ {\matrix{ x \cr y \cr z \cr } } \right] + \displaystyle{{p^2} \over {h{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}\displaystyle{d \over {dt}}\left[ {\matrix{ x \cr y \cr z \cr } } \right].$$

The inverse transformations are,

(13) $$\left[ {\matrix{ x \cr y \cr z \cr } } \right] = \displaystyle{1 \over {{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right],\quad \displaystyle{d \over {df}}\left[ {\matrix{ x \cr y \cr z \cr } } \right] = \displaystyle{{p^2} \over {h{{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}^2}}\displaystyle{d \over {dt}}\left[ {\matrix{ x \cr y \cr z \cr } } \right].$$

Thus

(14) $$\displaystyle{{{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}} \over {\dot f}}\displaystyle{d \over {dt}}\left[ {\matrix{ x \cr y \cr z \cr } } \right] = \displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{{e\sin {\mkern 1mu} f} \over {{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right].$$

Consequently, we obtain,

(15) $$\displaystyle{d \over {df}}\left[ {\matrix{ x \cr y \cr z \cr } } \right] = \displaystyle{{e\sin {\mkern 1mu} f} \over {{{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{1 \over {{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}\displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right].$$

Thus, the transformation relating the state vector takes the form, $\tilde{\bf{x}}=\bf{Tx}$ ,

(16) $$\matrix{ {\widetilde{{\bf x}}} & { = {\left[ {\matrix{ {\tilde x} & {\tilde y} & {\tilde z} & {{\tilde x}^{\prime}} & {{\tilde y}^{\prime}} & {{\tilde z}^{\prime}} \cr } } \right]}^T,{\bf x} = {\left[ {\matrix{ x & y & z & {\dot x} & {\dot y} & {\dot z} \cr } } \right]}^T,} \cr } $$

(17) $$\matrix{ {\bf T} & { = \left[ {\matrix{ {{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}{\bf I}} & {\bf 0} \cr { - {\rm (}e\sin {\mkern 1mu} f{\rm )}{\bf I}} & {\displaystyle{{p^2} \over {h{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}{\bf I}} \cr } } \right],\quad {\bf T}^{ - 1} = \left[ {\matrix{ {\displaystyle{1 \over {{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}{\bf I}} & {\bf 0} \cr {\displaystyle{{h{\rm (}e\sin {\mkern 1mu} f{\rm )}} \over {p^2}}{\bf I}} & {\displaystyle{{h{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}} \over {p^2}}{\bf I}} \cr } } \right].} \cr } $$

Furthermore,

(18) $$\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = - e\cos {\mkern 1mu} f\left[ {\matrix{ x \cr y \cr z \cr } } \right] + \displaystyle{{p^2} \over {h{\rm (}1 + e\cos {\mkern 1mu} f{\rm )}}}\displaystyle{d \over {dt}}\displaystyle{d \over {df}}\left[ {\matrix{ x \cr y \cr z \cr } } \right].$$

Hence,

(19) $$\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = - e\cos {\mkern 1mu} f\left[ {\matrix{ x \cr y \cr z \cr } } \right] + \displaystyle{{1 + e\cos {\mkern 1mu} f} \over {{\dot f}^2}}\displaystyle{{d^2} \over {dt^2}}\left[ {\matrix{ x \cr y \cr z \cr } } \right].$$

It follows that,

(20) $$\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = - \left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{1 \over {1 + e\cos {\mkern 1mu} f}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{{1 + e\cos {\mkern 1mu} f} \over {{\dot f}^2}}\displaystyle{{d^2} \over {dt^2}}\left[ {\matrix{ x \cr y \cr z \cr } } \right].$$

Noting that $p/r=\lpar {1+e\cos\, f} \rpar $ , from Equation (8a) we have,

(21) $$\matrix{ {\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right]} & { = - \left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{r \over p}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - {\bf K}_{lerm}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - {\bf G}_{lerm}\left\{ {\displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{{re\sin {\mkern 1mu} f} \over p}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right]} \right\}{\rm }} \cr {} & {\quad + \displaystyle{r \over {p{\dot f}^2}}\left[ {\matrix{ {{\tilde a}_{nx}} \cr {{\tilde a}_{ny}} \cr {{\tilde a}_{nz}} \cr } } \right].} \cr } $$

Thus,

(22) $$\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = - {\rm (}{\bi K}_{lerm} + {\bi I}{\rm )}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - \displaystyle{{re\sin f} \over p}{\bi G}_{lerm}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \displaystyle{{{\bf I}r} \over {p}} \left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - {\bi G}_{lerm}\displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right]{\rm }\quad + \displaystyle{r \over {p{\dot f}^2}}\left[ {\matrix{ {{\tilde a}_{nx}} \cr {{\tilde a}_{ny}} \cr {{\tilde a}_{nz}} \cr } } \right].{\rm }$$

But,

(23) $$\eqalign{&{{\bf K}_{lerm} + {\bf I} + \displaystyle{{r{\rm (}e\sin {\mkern 1mu} f{\rm )}} \over p}{\bf G}_{lerm} - \displaystyle{r \over p}{\bf I}} \cr & {\quad = \left( {\displaystyle{1 \over {\dot f}} - 1} \right)\displaystyle{{2re\sin {\mkern 1mu} f} \over p}\left[ {\matrix{ 0 & 2 & 0 \cr { - 2} & 0 & 0 \cr 0 & 0 & 0 \cr } } \right] + \displaystyle{r \over p}\left[ {\matrix{ { - 2 - 1} & 0 & 0 \cr 0 & {1 - 1} & 0 \cr 0 & 0 & {1 + 1 + e\cos {\mkern 1mu} f} \cr } } \right].}} $$

Hence,

(24) $$\matrix{ & {{\bf K}_{lerm} + {\bf I} + \displaystyle{{{\rm (}e\sin f{\rm )}} \over {{\rm (}1 + e\cos f{\rm )}}}{\bf G}_{lerm} - \displaystyle{1 \over {1 + e\cos f}}{\bf I}{\rm }} \cr & { = \left( {\displaystyle{{p^2} \over {h{{\rm (}1 + e\cos f{\rm )}}^2}} - 1} \right)\displaystyle{{2e\sin f} \over {{\rm (}1 + e\cos f{\rm )}}}\left[ {\matrix{ 0 & 2 & 0 \cr { - 2} & 0 & 0 \cr 0 & 0 & 0 \cr } } \right]{\rm }} \cr & {\quad + \displaystyle{1 \over {{\rm (}1 + e\cos f{\rm )}}}\left[ {\matrix{ { - 3} & 0 & 0 \cr 0 & 0 & 0 \cr 0 & 0 & {2 + e\cos f} \cr } } \right].} } $$

Hence if we let,

(25) $$\eqalign{ {{\widetilde{{\bf K}}}_{lerm}} & { = \displaystyle{1 \over {{\rm (}1 + e\cos f{\rm )}}}\left\{ {\left[ {\matrix{ { - 3} & 0 & 0 \cr 0 & 0 & 0 \cr 0 & 0 & {2 + e\cos f} \cr } } \right]} \right.{\rm }} \cr & {\quad \left. { + 2e\sin f\left( {\displaystyle{{p^2} \over {h{{\rm (}1 + e\cos f{\rm )}}^2}} - 1} \right)\left[ {\matrix{ 0 & 2 & 0 \cr { - 2} & 0 & 0 \cr 0 & 0 & 0 \cr } } \right]} \right\},{\rm }} \cr {\quad {\widetilde{{\bf G}}}_{lerm}} & { = \left[ {\matrix{ 0 & { - 2} & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \cr } } \right],} } $$

the nonlinear controlled TH equations with control inputs $ \tilde{\bf{u}} \lpar \tilde{\bf{x}} \rpar $ are given by,

(26) $$\matrix{ {\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right]} & { = - {\widetilde{{\bf K}}}_{lerm}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - {\widetilde{{\bf G}}}_{lerm}\displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right]{\rm }} \cr & {\quad + \displaystyle{{p^4} \over {h^2{{\rm (}1 + e\cos f{\rm )}}^3}}\left[ {\matrix{ {{\tilde a}_{nx}{\rm (}\tilde x,\tilde y,\tilde z{\rm )}} \cr {{\tilde a}_{ny}{\rm (}\tilde x,\tilde y,\tilde z{\rm )}} \cr {{\tilde a}_{nz}{\rm (}\tilde x,\tilde y,\tilde z{\rm )}} \cr } } \right] + \widetilde{{\bf u}}{\rm (}\widetilde{{\bf x}}{\rm )},} \cr } $$

where $\tilde {a}_{n^{\ast}} \lpar \tilde {x}\comma \; \tilde {y}\comma \; \tilde {z}\rpar $ are the perturbation accelerations expressed in terms of the scaled coordinates. Given a full state control input vector in the scaled TH frame as,

(27) $$\widetilde{{\bf u}}{\rm (}\widetilde{{\bf x}}{\rm )} = - \widetilde{{\bf K}}_f\widetilde{{\bf x}} = - \widetilde{{\bf K}}_f\left[ {\matrix{ {\tilde x} & {\tilde y} & {\tilde z} & {\displaystyle{{d\tilde x} \over {df}}} & {\displaystyle{{d\tilde y} \over {df}}} & {\displaystyle{{d\tilde z} \over {df}}}} } \right]^T,$$

since,

(28) $$\widetilde{{\bf x}} = {\bf Tx},\quad {\bf T} = \left[ {\matrix{ {{\rm (}1 + e\cos f{\rm )}{\bf I}} & {\bf 0} \cr { - {\rm (}e\sin f{\rm )}{\bf I}} & {\displaystyle{{p^2} \over {h{\rm (}1 + e\cos f{\rm )}}}{\bf I}}} } \right],$$

the physical control input vector is,

(29) $$\matrix{ & {{\bf u}{\rm (}{\bf x}{\rm )} = \displaystyle{{h^2{{\rm (}1 + e\cos f{\rm )}}^3} \over {p^4}}\widetilde{{\bf u}}{\rm (}{\bf Tx}{\rm )} = - \displaystyle{{h^2{{\rm (}1 + e\cos f{\rm )}}^3} \over {p^4}}{\widetilde{{\bf K}}}_f{\bf Tx}{\rm }} \cr & { = - \displaystyle{{h^2{{\rm (}1 + e\cos f{\rm )}}^3} \over {p^4}}{\widetilde{{\bf K}}}_f{\bf T}{\left[ {\matrix{ x & y & z & {\displaystyle{{dx} \over {df}}} & {\displaystyle{{dy} \over {df}}} & {\displaystyle{{dz} \over {df}}} \cr } } \right]}^T.} } $$

Thus, in terms of the original state variables the physical control input vector is,

(30) $${\bf u}{\rm (}{\bf x}{\rm )} = - {\bf K}_f{\bf x} = - {\bf K}_f\left[ {\matrix{ x & y & z & {\displaystyle{{dx} \over {df}}} & {\displaystyle{{dy} \over {df}}} & {\displaystyle{{dz} \over {df}}} \cr } } \right]^T,$$

with

(31) $$\matrix{ {{\bf K}_f} & { = - \displaystyle{{h^2{{\rm (}1 + e\cos f{\rm )}}^3} \over {p^4}}{\widetilde{{\bf K}}}_f{\bf T}{\rm }} \cr & { = - \displaystyle{{h^2{{\rm (}1 + e\cos f{\rm )}}^3} \over {p^4}}{\widetilde{{\bf K}}}_f\left[ {\matrix{ {{\rm (}1 + e\cos f{\rm )}{\bf I}} & {\bf 0} \cr { - {\rm (}e\sin f{\rm )}{\bf I}} & {\displaystyle{{p^2} \over {h{\rm (}1 + e\cos f{\rm )}}}{\bf I}} \cr } } \right].} } $$

4. LINEAR QUADRATIC NONLINEAR GAUSSIAN (LQNG) DESIGN

The strategy adopted to design the controller is to adopt the methodology of the linear quadratic regulator for the synthesis of the control laws. To apply the methodology of the linear quadratic regulator, Equations (26) must first be linearized. To linearize Equations (26) the control input vector is expressed as, ${\tilde{\bf{u}}}\!\lpar {{\tilde{\bf{x}}}} \rpar ={\tilde{\bf{u}}}_{0} \!\lpar {{\tilde{\bf{x}}}} \rpar +{\tilde{\bf{u}}}_{1}\! \lpar {{\tilde{\bf{x}}}} \rpar $ . The control input ${\tilde{\bf{u}}}_{0}\! \lpar {{\tilde{\bf{x}}}} \rpar $ is chosen in terms of $\hat {\tilde {x}}$ , $\hat {\tilde {y}}$ and $\hat {\tilde {z}}$ which are estimates of $\tilde x$ , $\tilde y$ and $\tilde z$ respectively as,

(32) $$\widetilde{{\bf u}}_0(\widetilde{{\bf x}}) = - \displaystyle{{p^4} \over {h^2{(1 + e\cos f)}^3}}\left[ \matrix{ {{\tilde a}_{nx}({\hat {\tilde x}},{\hat {\tilde y}},{\hat {\tilde z}})} \cr {{\tilde a}_{ny}({\hat {\tilde x}},{\hat {\tilde y}},{\hat {\tilde z}})} \cr {{\tilde a}_{nz}({\hat {\tilde x}},{\hat {\tilde y}},{\hat {\tilde z}})} \cr } \right],$$

Thus Equations (26) reduce to,

(33) $$\eqalign{ \displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] & = - {{\tilde {\bf K}}}_{lerm}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - {{\bf \tilde G}}_{lerm}\displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] \cr & {\quad + \displaystyle{{p^4} \over {h^2{(1 + e\cos f)}^3}}\left[ {\matrix{ {{\tilde a}_{nx}(\tilde x,\tilde y,\tilde z) - {\tilde a}_{nx}({\hat {\tilde x}},{\hat {\tilde y}},{\hat {\tilde z}})} \cr {{\tilde a}_{ny}(\tilde x,\tilde y,\tilde z) - {\tilde a}_{ny}({\hat {\tilde x}},{\hat {\tilde y}},{\hat {\tilde z}})} \cr {{\tilde a}_{nz}(\tilde x,\tilde y,\tilde z) - {\tilde a}_{nz}({\hat {\tilde x}},{\hat {\tilde y}},{\hat {\tilde z}})} \cr } } \right] + {{\tilde {\bf u}}}_1({\tilde {\bf x}}).} \cr } $$

The estimates of the perturbation accelerations components $\tilde {a}_{n\ast } \lpar {\hat {\tilde {x}}\comma \; \hat {\tilde {y}}\comma \; \hat {\tilde {z}}} \rpar $ are evaluated at the mid-point of the current integration step. Thus, the right-hand side of Equation (33) represents the error in the actual and estimated perturbation accelerations, which can be assumed to be a disturbance vector and ignored for purposes of designing the control input ${\tilde{\bf{u}}}_{1}\! \lpar {{\tilde{\bf{x}}}} \rpar $ .

However, it is often not possible to measure all the states of the system, which is in reality a nonlinear system. Moreover, the measurements are almost always noisy. For this reason, we adopt a nonlinear estimation methodology to reconstruct all the states from the measurements. The nonlinear estimation methodology is based on the unscented Kalman filter. Because the equations of motion are nonlinear, the separation theorem, which is often adopted to justify the independent design of the control laws and the solution of the state estimation problem, is no longer valid. Nevertheless, we will design the control input ${\tilde{\bf{u}}}_{1} \!\lpar {{\tilde{\bf{x}}}} \rpar $ based on the linearized true anomaly-dependent discretized equations for the relative motion and implement these inputs by feeding back the estimates of the states rather than the original states of the system.

5. THE LINEAR TH EQUATIONS AND THEIR DISCRETIZATION

The linearized TH equations in the presence of control accelerations are given by,

(34) $$\displaystyle{{d^2} \over {df^2}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] = - \widetilde{{\bf K}}_{lerm}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] - \widetilde{{\bf G}}_{lerm}\displaystyle{d \over {df}}\left[ {\matrix{ {\tilde x} \cr {\tilde y} \cr {\tilde z} \cr } } \right] + \left[ {\matrix{ 1 & 0 & 0 \cr 0 & 1 & 0 \cr 0 & 0 & 1 \cr } } \right]\left[ {\matrix{ {{\tilde u}_{1x}} \cr {{\tilde u}_{1y}} \cr {{\tilde u}_{1z}} \cr } } \right].$$

Equations (33) may be expressed as a set of first order equations in matrix form with the true anomaly as the independent variable,

(35a) $$\displaystyle{d \over {df}}\widetilde{{\bf x}} = {\bf A}{\rm (}f{\rm )}\widetilde{{\bf x}} + \left[ {\matrix{ {\bf 0} & {\bf I} \cr } } \right]^T\widetilde{{\bf u}}_1,$$

with,

(35b) $${\bf A}{\rm (}f{\rm )} \equiv \left[ {\matrix{ {\bf 0} & {\bf I} \cr { - {\widetilde{{\bf K}}}_{lerm}} & { - {\widetilde{{\bf G}}}_{lerm}} \cr } } \right].$$

The coefficients of the equation are periodic with a period T = 2π. Thus, there exists a non-singular matrix ${\bf P} \lpar f \rpar $ , periodic with period T = 2π with ${\bf P} \lpar 0 \rpar = {\bf I}$ , such that the change of the variables ${\tilde{\bf{x}}}={\bf P}\lpar f \rpar {\bf z}$ , transforms the system into a linear system with constant coefficients. The Floquet transition matrix which is also known as the monodromy matrix $\Phi \lpar {f\comma \; 0} \rpar $ , which determines the stability of the system with f = T, satisfies,

(36) $$\displaystyle{d \over {df}}\Phi {\rm (}f,0{\rm )} = {\bf A}{\rm (}f{\rm )}\Phi {\rm (}f,0{\rm )},\quad \Phi {\rm (}0,0{\rm )} = {\bf I}.$$

The solution ${\tilde{\bf{x}}}\!\lpar f \rpar $ then satisfies,

(37) $$\widetilde{{\bf x}}{\rm (}f{\rm )} = \Phi {\rm (}f,0{\rm )}\widetilde{{\bf x}}{\rm (}0{\rm )} + \int\limits_0^f {\Phi {\rm (}f,\varphi {\rm )}{\left[ {\matrix{ {\bf 0} & {\bf I} \cr } } \right]}^T{\widetilde{{\bf u}}}_1d\varphi } ,$$

which can also be expressed as,

(38) $$\widetilde{{\bf x}}{\rm (}f + \Delta f{\rm )} = \Phi {\rm (}f + \Delta f,f{\rm )}\widetilde{{\bf x}}{\rm (}f{\rm )} + \int\limits_0^{\Delta f} {\Phi {\rm (}f + \Delta f,f + \varphi {\rm )}{\left[ {\matrix{ {\bf 0} & {\bf I} \cr } } \right]}^T{\widetilde{{\bf u}}}_1{\rm (}f + \varphi {\rm )}d\varphi } .$$

Assuming that the control input could be approximated as a constant over each integration interval and given by $ \bar{\bf{u}} ={\tilde{\bf{u}}}_{1} $ , we have the discrete system,

(39) $${\tilde{\bf{x}}}\!\lpar {f+\Delta f} \rpar =\Phi \!\lpar {f+\Delta f,f} \rpar {\tilde{\bf{x}}}\!\lpar f \rpar +\Psi \!\lpar {f+\Delta f,f} \rpar \!\bar{bf{u}}\!\lpar f \rpar ,$$

with,

(40) $$\Psi {\rm (}f + \Delta f,f{\rm )} = \int\limits_0^{\Delta f} {\Phi {\rm (}f + \Delta f,f + \varphi {\rm )}{\left[ {\matrix{ {\bf 0} & {\bf I} \cr } } \right]}^Td\varphi } .$$

The discrete system given by Equation (39) may be expressed as,

(41) $${\tilde{\bf{x}}}\!\lpar {k+1} \rpar ={\tilde{bf{A}}}\!\lpar k \rpar \!{\tilde{\bf{x}}}\!\lpar k \rpar +{\tilde{\bf{B}}}\!\lpar k \rpar \bar{bf{u}}\!\lpar k \rpar .$$

Equation (41) may be used to synthesise an optimal control law.

6. LINEAR QUADRATIC REGULATOR DESIGN

The design process for the discrete-time Linear Quadratic Regulator (LQR) with a finite final time when magnitude constraints are imposed on the state and control variables is discussed in this section. A finite final time quadratic performance index could be used, so the ${\bf{Q}}$ and ${\bf{R}}$ weighting matrices will have to be chosen and an algebraic Riccati equation must be solved. In the absence of terminal constraints, the performance index is of the form:

(42) $$J = \displaystyle{1 \over 2}\sum\limits_{k = 0}^{N - 1} {{\rm (}\widetilde{{\bf x}}_k^T {\bf Q}{\widetilde{{\bf x}}}_k + \overline {\bf u} _k^T {\bf R}{\overline {\bf u} }_k{\rm )}} .$$

When a spacecraft is actuated upon by thrusters, the thrusts generated are small. Typically, a thruster generates thrusts from 1-2 N. Thus, it is essential that in the design of the control laws the elements of the matrix ${\bf{R}}$ are of the order of 10¹ or more. During a controlled rendezvous the position coordinates are reduced from 10-100 m to less than a millimetre. Thus, the elements of the matrix ${\bf{Q}}$ are of the order of 10⁵ or less. An iterative process is used to update the weighting matrices used in the performance index to meet the terminal constraints. An alternate but equivalent approach is to use the Model Predictive Control (MPC) approach. To briefly describe the synthesis of linear optimal control based on the MPC approach, consider a linear discrete time system in the form,

(43) $${\bf{x}}\!\lpar {k+1} \rpar ={\bf{A}}\!\lpar k \rpar {\bf{x}}\!\lpar k \rpar +{\bf{B}}\!\lpar k \rpar {\bf{u}}\!\lpar k \rpar , \quad {\bf{y}}\!\lpar k \rpar ={\bf{C}}\!\lpar k \rpar {bf{x}}\!\lpar k \rpar .$$

Our aim is to find an optimal control input sequence defined over a control prediction window, ${ \bf{u}}\!\lpar j \rpar $ , $j=0\comma \; 1\comma \; 2\ldots N-1$ or the vector ${\bf U} = \left[ {\matrix{ {{\bf u}^T{\rm (}0{\rm )}} & {{\bf u}^T{\rm (}1{\rm )}} & \cdots & {{\bf u}^T{\rm (}N - 1{\rm )}} \cr } } \right]^T$ so as to minimise the performance index, which now includes a terminal weighting matrix ${\bf{Q}}_{N} $ ,

(44) $$J\!\lpar {{\bf{x}}\!\lpar 0 \rpar ,{\bf{U}}} \rpar =\sum\limits_{k=0}^{N-1}\! {\lcub {{\bf{y}}^T\!\lpar k \rpar q{\bf{Qy}}\!\lpar k \rpar +{\bf{u}}^T\!\lpar k \rpar r{\bf{Ru}}\!\lpar k \rpar } \rcub +{\bf{y}}^T\!\lpar N \rpar q{\bf{Q}}_N {\bf{y}}\!\lpar N \rpar },$$

where q and r are scalar scaling parameters which are used to re-scale the relative contributions of the states and the control inputs to the cost function. When the magnitude of the control input is relatively large when compared with the magnitude of the state, r is assumed to be a small number relative to q and vice-versa.

Defining the vector, ${\bf X} = \left[ {\matrix{ {{\bf x}^T{\rm (}1{\rm )}} & {{\bf x}^T{\rm (}2{\rm )}} & \cdots & {{\bf x}^T{\rm (}N - 1{\rm )}} & {{\bf x}^T{\rm (}N{\rm )}} \cr } } \right]$ , we may write,

(45) $$ J\!\lpar {{\bf{x}}\!\lpar 0 \rpar ,{\bf{U}}} \rpar =q{\bf{x}}^T\lpar 0 \rpar {\bf{C}}^T\lpar 0 \rpar {\bf{QC}}\lpar 0 \rpar {\bf{x}}\lpar 0 \rpar +q{\bf{X}}^T \bar{\bf{Q}}\bf{X}+r{\bf{ U}}^T \bar{\bf{R}}\bf{U},$$

where $\bar{\bf{Q}}$ is a block diagonal matrix with matrix ${\bf{C}}^{T}\lpar k \rpar {\bf{QC}}\lpar k \rpar $ , $k=1\comma \; 2\ldots N-1$ along the diagonal except the last element which is ${\bf{C}}^{T}\lpar N \rpar {\bf{Q}}_{N} {\bf{C}}\lpar N \rpar $ ; $\bar{\bf{R}}$ is a block diagonal matrix with matrix $\bf{R}$ along the diagonal. Using the state space model Equation (43) recursively, we may construct a prediction model in the form,

(46) $${\bf{X}}={\bf{SU}}+{\bf{Tx}}\!\lpar 0 \rpar $$

where

(47) $${\bf S} = \left[ {\matrix{ {\bf B} & {\bf 0} & {\bf 0} & {\bf 0} \cr {{\bf A}{\rm (}{\rm 1}{\rm )}{\bf B}} & {\bf B} & {\bf 0} & {\bf 0} \cr \cdots & \cdots & \cdots & \cdots \cr {\left( {\prod\limits_{k = N - {\rm 1}}^{\rm 1} {\bf A} {\rm (}k{\rm )}} \right){\bf B}} & {\left( {\prod\limits_{k = N - {\rm 2}}^{\rm 1} {\bf A} (k{\rm )}} \right){\bf B}} & \cdots & {\bf B} \cr } } \right],{\bf T} = \left[ {\matrix{ {{\bf A}{\rm (}{\rm 0}{\rm )}} \cr {{\bf A}{\rm (}{\rm 1}{\rm )}{\bf A}{\rm (}{\rm 0}{\rm )}} \cr \cdots \cr {\prod\limits_{k = N - {\rm 1}}^{\rm 0} {\bf A} {\rm (}k{\rm )}} \cr } } \right].$$

We have assumed for simplicity that ${\bf{B}}$ is constant but not ${\bf{A}}$ . The above formula could be generalised for the case when ${\bf{B}}$ is not constant. Thus, the cost function may be expressed as,

(48) $$J\lpar {{\bf{x}}\!\lpar 0 \rpar ,{\bf{U}}} \rpar =\displaystyle{{1}\over{2}}{\bf{x}}\!\lpar 0 \rpar ^T{\bf{Gx}}\!\lpar 0 \rpar +\displaystyle{{1}\over{2}}{\bf{U}}^T{\bf{HU}}+{\bf{x}}^T\lpar 0 \rpar {\bf{FU}}, $$

with, ${\bf{G}}=2q\lpar {{\bf{T}}^{T} {\bar{\bf{Q}}\bf{T}}+{\bf{C}}^{T}\lpar 0 \rpar {\bf{QC}}\lpar 0 \rpar } \rpar $ , ${\bf{H}}=2r \bar{\bf{R}}+2q{\bf{S}}^{T} \bar{\bf{Q}}\bf{S}$ and $\bf{F}=q\bf{T}^{T} \bar{\bf{Q}}\bf{S}$ .

The optimum control sequence is obtained by setting the gradient of $J\!\lpar {{\bf{x}}\!\lpar 0 \rpar \comma \; {\bf{U}}} \rpar $ to zero. Minimising the cost function results in,

(49) $$\eqalign{ {dJ{\rm (}{\bf x}{\rm (}0{\rm )},{\bf U}{\rm )}/d{\bf U}} & { = {\bf U}^T{\bf H} + {\bf x}^T{\rm (}0{\rm )}{\bf F} = 0\quad \Rightarrow \quad {\bf HU} + {\bf F}^T{\bf x}{\rm (}0{\rm )} = 0\quad \Rightarrow {\bf U}} \cr {} & { = - {\bf H}^{ - 1}{\bf F}^T{\bf x}{\rm (}0{\rm )}.}} $$

The state ${\bf{x}}\!\lpar 0 \rpar $ , at the start of the prediction window, is assumed to represent the state at the next time instant, in real time. The control law based on the receding horizon is,

(50) $${\bf u}{\rm (}k{\rm )} = - \left[ {\matrix{ 1 & 0 & \cdots & 0 \cr } } \right]{\bf H}^{ - 1}{\bf F}^T{\bf x}{\rm (}k{\rm )}.$$

The control sequence is recursively calculated over successive control prediction windows. The product ${\bf{H}}^{-1}{\bf{F}}^{T}$ may be expressed as

(51) $${\bf H}^{ - 1}{\bf F}^T = \displaystyle{q \over 2}\left( {r\overline {\bf R} + q{\bf S}^T\overline {\bf Q} {\bf S}} \right)^{ - {\bf 1}}{\bf S}^T\overline {\bf Q} {\bf T} = \displaystyle{{q^{\bf 2}} \over 2}\left( {\displaystyle{r \over q}\overline {\bf R} + q{\bf S}^T\overline {\bf Q} {\bf S}} \right)^{ - {\bf 1}}{\bf S}^T\overline {\bf Q} {\bf T}.$$

If we let q = 1, Equation (51) reduces to,

(52) $${\bf H}^{ - {\rm 1}}{\bf F}^T = \displaystyle{{\rm 1} \over {\rm 2}}{\rm (}r\overline {\bf R} + {\bf S}^T\overline {\bf Q} {\bf S}{\rm )}^{ - {\rm 1}}{\bf S}^T\overline {\bf Q} {\bf T},$$

where r is treated as a free parameter to be chosen.

7. THE UNSCENTED KALMAN FILTER

The basic unscented Kalman filter is identical to the filter implemented in Vepa and Amzhari (Reference Vepa and Zhahir2011) and Vepa (Reference Vepa2017) which is briefly summarised in the following.

Consider a random variable $\bf{w}$ with dimension L which is going through the nonlinear transformation, $\bf{y} = \bf{f}\lpar \bf{w}\rpar $ . The initial conditions are that $\bf{w}$ has a mean $\bar{\bf{w}}$ and a covariance ${\bf{P}}_{ww} $ . To calculate the statistics of $\bf{y}$ , a matrix χ of 2L + 1 sigma vectors is formed. We have chosen to use the scaled unscented transformation proposed by Julier (Reference Julier2002), as this transformation gives the added flexibility of scaling the sigma points to ensure that the covariance matrices are always positive definite.

Given a general discrete-time nonlinear dynamic system in the form,

(53) $${\bf{x}}_{k+1} ={\bf{f}}_k \lpar {{\bf{x}}_k ,{\bf{u}}_k } \rpar +{\bf{w}}_k , \quad{\bf{y}}_k ={\bf{h}}_k \lpar {{\bf{x}}_k } \rpar +{\bf{v}}_k$$

where ${\bf{x}}_{k} \in R^{n}$ is the state vector, ${\bf{u}}_{k} \in R^{r}$ is the known input vector, ${\bf{y}}_{k} \in R^{m}$ is the output vector at time k. ${\bf{w}}_{k} $ and ${\bf{v}}_{k}$ are, respectively, the disturbance or process noise and sensor noise vectors, which are assumed to be Gaussian white noise with zero mean. Furthermore ${\bf{Q}}_{k}$ and ${\bf{R}}_{k}$ are assumed to be the covariance matrices of the process noise sequence, ${\bf{w}}_{\bf{k}}$ and the measurement noise sequence, ${\bf{v}}_{\bf{k}}$ respectively. The unscented transformations of the states corresponding to the functions ${\bf{f}}_{k} \lpar {{\bf{x}}_{k} \comma \; {\bf{u}}_{k} } \rpar $ and ${\bf{h}}_{k} \lpar {{\bf{x}}_{k} } \rpar $ are denoted respectively as,

(54) $$ {\bf{f}}_k^{UT} ={\bf{f}}_k^{UT} \lpar {{\bf{x}}_k ,{\bf{u}}_k } \rpar , \quad {\bf{h}}_k^{UT} ={\bf{h}}_k^{UT} \lpar {{\bf{x}}_k } \rpar $$

while the transformed covariance matrices and cross-covariance are respectively denoted as,

(55) $${\bf P}_k^{ff} = {\bf P}_k^{ff} {\rm (}\widehat{{\bf x}}_k,{\bf u}_k{\rm )},\quad {\bf P}_k^{hh - } = {\bf P}_k^{hh} {\rm (}\widehat{{\bf x}}_k^ - {\rm )and}{\bf P}_k^{xh - } = {\bf P}_k^{xh - } {\rm (}\widehat{{\bf x}}_k^ - ,{\bf u}_k{\rm )}.$$

The UKF estimator can then be expressed in a compact form. The state time-update equation, the propagated covariance, the Kalman gain ${\bf{K}}_{k}$ , the state estimate $\hat{\bf{x}}_{k}$ and the updated covariance $\hat{\bf{P}}_{\bf{k}}$ are respectively given by,

(56a) $$\hat{\bf{x}}_k^{-} ={\bf{f}}_{k-1}^{UT} \lpar {\hat{\bf{x}}}_{k-{\bf{1}}} \rpar $$

(56b) $$\hat{\bf{P}}_k^- ={\bf{P}}_{k-1}^{ff} +{\bf{Q}}_{k-1}$$

(56c) $${\bf K}_k = \widehat{{\bf P}}_k^{xh - } {\rm (}\widehat{{\bf P}}_k^{hh - } + {\bf R}_k{\rm )}^{ - 1}$$

(56d) $$\widehat{{\bf x}}_k = \widehat{{\bf x}}_k^ - + {\bf K}_k[{\bf z}_k - {\bf h}_k^{UT} {\rm (}\widehat{{\bf x}}_k^ - {\rm )}]$$

(56e) $$\widehat{{\bf P}}_{\rm k} = \widehat{{\bf P}}_k^ - - {\bf K}_k{\rm (}\widehat{{\bf P}}_k^{hh - } + {\bf R}_k{\rm )}^{ - {\rm 1}}{\bf K}_k^T .$$

Equations (56a)-(56e) are in the same form as the traditional Kalman filter and the extended Kalman filter. Thus, higher order non-linear models capturing significant aspects of the dynamics may be employed to ensure that the Kalman filter algorithm can be implemented to effectively estimate the states in practice. The control input $\tilde{\bf{u}}_{1} \lpar \tilde{\bf{x}} \rpar $ is assumed to be given by,

(57) $$\widetilde{{\bf u}}_1{\rm (}\widetilde{{\bf x}}{\rm )} = \overline {\bf u} {\rm (}k{\rm )} = - \left[ {\matrix{ 1 & 0 & \cdots & 0 \cr } } \right]{\bf H}^{ - {\rm 1}}{\bf F}^T\widehat{{\widetilde{{\bf x}}}}{\rm (}k{\rm )},$$

where the vector $\hat{\tilde{\bf{x}}}$ is an estimate of the vector $\hat{\tilde{\bf{x}}}$ obtained by the unscented Kalman filter.

8. TYPICAL SIMULATIONS AND RESULTS

The initial position and velocity of the primary satellite are obtained from the initial classical orbital elements defined in Table 1. From these orbital elements, the initial position and velocity vector of the satellite are obtained. The actual initial location of the secondary orbiting body is defined by the orbital elements set shown in Table 1 and converted to position and velocity vectors relative to the satellite. The orbital elements defining the reference Keplerian orbit are identical to the primary satellite's initial orbital elements.

Table 1.

Orbital elements defining the initial position of the Primary Satellite and orbiting body.

In all three cases the size and shape of the orbital ellipse and the inclination of the orbit plane are assumed to be the same. Only the argument of the ascending node (Ω) is assumed to be different by 0·00001^°. From the orbital elements, the initial position and velocity of the secondary orbiting body are defined relative to the satellite. The numerical integration was performed using the Runge-Kutta Dormand-Prince 4(5) (Dormand and Prince, Reference Dormand and Prince1980). The overall prediction and estimation window is in excess of 2π radians which is covered in 7,500 integration steps. The non-dimensionalising distance was chosen to be $r_{sc} = 35,786\,{\rm km}$ , which is the radius of the geostationary orbit, and the unit of time was t _sc  = 86,400 sec or a day. The position coordinates used for integration are the scaled TH coordinates. All the results are plotted in SI units.

A total of 13 noisy measurements are made: a pair of angles that the satellite makes with the Sun, and another pair with the Moon, the satellite's distance vector from the Earth which is assumed to be relatively very noisy or uncertain and the satellite's relative position and velocity vectors relative to the secondary orbiting body. The positions of the Sun and Moon are assumed to be known from the relevant ephemeris and obtained in the same manner as in Vepa (Reference Vepa2017).

Figure 1 presents the simulated and estimated position coordinates of the primary satellite, relative to the secondary orbiting body in the closed loop case, which are compared. The errors in the relative position coordinate estimates, relative to the secondary orbiting body are shown in Figure 2. Figure 3 presents the corresponding simulated and estimated velocity components of the satellite relative to the secondary orbiting body in the closed loop case, which are compared against each other. The errors in these relative velocity components are shown in Figure 4. The control laws were evaluated both by applying the discrete linear quadratic regulator algorithm (with an infinite horizon) at every true anomaly step as well as by the MPC algorithm over a finite control prediction horizon. The results in both cases were quite similar and for this reason only the latter results are presented here.

Figure 1.

The simulated and estimated [x  y  z] position coordinates of the primary satellite, relative to the secondary orbiting body in the closed loop case with the optimal control law.

Figure 2.

The errors in the simulated and estimated [x  y  z] position coordinates of the primary satellite, relative to the secondary orbiting body.

Figure 3.

The simulated and estimated [x  y  z] velocity components of the primary satellite, relative to the secondary orbiting body in the closed loop case.

Figure 4.

The errors in the simulated and estimated [x  y  z] velocity components of the primary satellite, relative to the secondary orbiting body.

Figures 5 and 6 show the estimates of the second body position and velocity in the TH coordinates, relative to a particle in a notional Keplerian orbit. Figures 7 and 8 show the corresponding errors in the estimates of the second body position and velocity. Although these errors are an order of magnitude greater than the errors in the satellites' relative position and velocity estimates they are still less than 0·16% of the secondary orbiting body's absolute position and velocity component magnitudes relative to the Earth's centre.

Figure 5.

The simulated and estimated [x  y  z] position coordinates of the secondary orbiting body relative to the position of a notional particle in a Keplerian orbit.

Figure 6.

The simulated and estimated [x  y  z] velocity components of the secondary orbiting body relative to the position of a notional particle in a Keplerian orbit.

Figure 7.

Errors in the simulated and estimated [x  y  z] position coordinates of the secondary orbiting body relative to the position of a notional particle in a Keplerian orbit.

Figure 8.

Errors in the simulated and estimated [x  y  z] velocity components of the secondary orbiting body relative to the position of a notional particle in a Keplerian orbit.

Examining the results in Figure 1 closely it is clear that there is an accumulation of position error in the radial direction which approaches 4 m and that the satellite is falling behind in its pursuit of the secondary orbiting object. Moreover, the response resembles the typical response of a first order system to a step input as the error between the pursuer, the primary satellite, and the secondary orbiting body approaches a steady state value, albeit slowly, over one orbit. This is a classic case that demands from a practical standpoint the use of integral control to drive the position error, due to the residual disturbance acceleration, to zero. The control is now assumed to be,

(58) $$ \tilde{\bf{u}}_1\! \lpar \tilde{\bf{x}} \rpar = \bar{\bf{u}}\!\lpar k \rpar +K_I \sum\limits_{j=1}^{k-1} \bar{\bf{u}} \lpar j \rpar \Delta f,$$

where the integral matrix gain K _I is nominally chosen to be equal to,

(59) $$K_I = \left[ {\matrix{ 2 & 0 & 0 \cr 0 & {0.5} & 0 \cr 0 & 0 & 0 \cr } } \right],$$

Δ f is the integration step in terms of the true anomaly and $\bar{\bf{u}}\!\lpar k \rpar $ is the previously designed optimal control law.

The results corresponding to Figure 1, with the additional term representing integral control in the control law defined in Equation (58), are shown in Figure 9. There is now a significant reduction in the position error of the pursuer, the primary satellite, relative to the secondary orbiting body which approaches an almost steady state value after one orbit. The position error components are well below 4 m and decreasing over one orbit. The relative positions and relative velocities are within acceptable limits and rendezvous could be established as soon as the error components are near zero and attitude synchronisation is achieved.

Figure 9.

The simulated and estimated [x  y  z] position coordinates of the primary satellite, relative to the secondary orbiting body with the additional integral term in the control law.