Fixed ground-target tracking control of satellites using a nonlinear model predictive control

Fixed ground-target tracking control of satellites using a nonlinear model predictive control

Abd El Hady M. ElbeltagyAhmed M. Youssef Amgad M. Bayoumy Yehia Z. Elhalwagy 

Aircraft Mechanics Dep., MTC, Cairo, Egypt

Aircraft Electric Equipment Dep., MTC, Cairo, Egypt

Mechanical Engineering Dep., MSA, Cairo, Egypt

Navigation Dep., MTC, Cairo, Egypt

Corresponding Author Email:
20 Feburary 2018
15 March 2018
31 March 2018
| Citation



This paper proposes a novel solution to the fixed-ground target tracking control problem of satellites utilizing a Nonlinear Model Predictive Control approach (NMPC). The Continuation / Generalized Minimal Residual (C/GMRES) algorithm is selected as a promising fast solver to an optimal control problem in real time. The algorithm could perfectly deal with the huge computational load of this approach, represented in solving Riccati differential equation, by simple and efficient approximations.  A new control-oriented model converting the main tracking problem into a simple regulation problem is developed. This simple and easy traceable reformulated model has an advantage in dealing with modeling errors and unplanned external environmental disturbances. The update of the control input is obtained by integrating a deduced time-dependent inputs and Lagrange multipliers vector; representing the solution of a set of linear equations and corresponding to the optimality conditions. The proposed algorithm is simulated using real satellite parameters to track a fixed-ground target for reconnaissance purposes. The simulation results show that the algorithm of C/GMRES method can track a desired fixed ground target robustly, with precise tracking error and guaranteed safe stability limits for shooting activities throughout the overpass flight.


C/GMRES method, ground-target tracking, image quality, optimization, predictive control

1. Introduction

The improvement of different strategic satellite attitude control systems by commissioning state-of-the-art intelligent control methods is an attractive arena of research recently. The extremely variant nature of satellite attitude control indicates that the elected effective intelligent control technique must be able to deal with a large Multi-Input Multi-Output (MIMO) system with nonlinear complex and constrained dynamics.  Furthermore, in many conditions, it may also be necessary to avoid the different actuators from long operating near their thresholds to keep power consumption within some specified limits or to keep the rate of wear as low as possible.

Satellite Ground-Target Tracking (GTT) controllers must also guarantee a precise tracking error and fast response as well. NMPC, as a unique optimal model-based approach to intelligent control systems design, is a talented candidate that covers all of these features. The capability to handle Two-Point Boundary Value Problems (TPBVP) with such approach has made this technique superior. However, the main challenge faces NMPC controller developers is guaranteeing method implementation in a real-time in spite of the needed huge computational loads. Hence, fast Real-Time Optimization (RTO) methods are carefully searched to overcome this weakness. Optimization methods based on C/GMRES method are candidates of these RTO algorithms that have revealed great potential for real-time applications in general. Using the feature of the solution of a nonlinear receding horizon control problem that it varies smoothly related to the time at certain conditions is the main core idea in Ref. [1]. Thus, tracing this solution is possible without the need for iterative optimization methods. This is known by continuation method. Ohtsuka use this continuation method combined with a fast solver algorithm for linear equations instead of solving Riccati differential equation to develop his real-time algorithm of nonlinear receding horizon control problem. As he succeeded to manipulate the equations of nonlinear optimality conditions such that all the unknown variables which can’t obtained recursively are put in one vector. Solving the nonlinear equations that are a function of this defined vector is replaced by searching this vector which can satisfy the features of any stable dynamic system. Hence, a set of linear equations are obtained and solved in real-time by employing the GMRES method [2].

Explicit MPC (eMPC) has been recently proposed as an alternative to online approaches in a trial to overcome the implementation challenges of a huge computational load for NMPC in each sampling time [3]. The nonlinear optimization problem in eMPC approach is completely solved offline to cancel the necessity of having an online optimization solver. However, the optimal solution is handled as a look-up table of linear gains that must be stored, the problem converted to storage challenges rather than computational ones to fulfill the limitations of the on-board processor.

The C/GMRES method is reapplied in a different application with using a developed modular code generation tool that can be tailored to be used in other applications [4]. However, some code changes should be employed to allow covering the problems which comprise time-variant parameters that cannot considered as a system states.

The GMRES-based NMPC method has been used in designing ecological adaptive cruise controller (E-ACC) for plug-in Hybrid electric vehicles (HEV) [5]. In the developed E-ACC the upcoming trip data and an on-board vehicle radar reading are used to optimize energy cost of the trip and to preserve safety and comfort as well.

This paper investigates the potential of novel employing C/GMRES based RTO algorithms to design intelligent attitude control systems, in particular, intelligent fixed GTT controllers.  The paper is organized to present a model predictive control and a control problem solution in section (‎2). In section (‎3) C/GMRES method is introduced in details. Section (‎4) is devoted to creating control-oriented model. Section (‎5) gives a brief about image quality requirements. In section (‎6) a case study with its simulation results is presented followed by the final conclusion statements in section (‎7).

2. Control Problem Solution Using MPC

In the approach of MPC, or alternatively receding horizon control, the optimal solution of the control problem, called the control input is found by solving a constraint optimization problem at each time instant with a moving horizon. This results in a control law with a feedback nature. Actually, optimization problem of the constraint objective function can be converted into a root finding problem of a set of nonlinear equations. These nonlinear equations are obtained by means of applying necessary conditions for optimality [1]. 

A general nonlinear dynamic system can be represented as:

$\begin{aligned} \dot{x}(t)=& f(x(t), u(t), p(t)) \\ & x\left(t_{0}\right)=x_{o} \end{aligned} \}$ (1)

where, $x(t) \in R^{n}$ is a state vector with states, $u(t) \in R^{m_{u}}$ is an input vector with mu inputs and $p(t) \in R^{m_{p}}$ is a time-varying vector. The optimal control problem for this system where a cost function is minimized while satisfying ng-dimensional equality constraints g and ng-dimensional inequality constraints h are represented as [1]:

$\begin{aligned} \underset{u}{\operatorname{minimize}} : J=\varphi &\left(x\left(t_{f}\right), p\left(t_{f}\right)\right) \\ &+\int_{t_{0}}^{t_{f}} L(x(\dot{t}), u(\dot{t}), p(\dot{t})) d \dot{t} \end{aligned}$ $\left\{\begin{aligned} \dot{x}(t)=& f(x(t), u(t), p(t)) \\ & x\left(t_{0}\right)=x_{o} \\ g(x(t), u(t), p(t)) &=0 \\ h(x(t), u(t), p(t)) & \leq 0 \end{aligned}\right.$ (2)

where, $T=t_{f}-t_{0}$ is the prediction horizon,  $\varphi( .)$ represents the terminal cost at the end of the prediction horizon and $L( .)$ denotes the trajectory cost over the horizon. The actual state x(t) is used as the initial state for the optimal control problem over the horizon $[t, t+T]$ .

Although the control input is determined over this horizon, the actual submitted input as a feedback to the system is the only optimum value at the initial time of this current horizon [1,4,5].

Solving an optimization problem with inequality constraints is not an easy task because of the complexity of the optimality conditions known as Karush–Kuhn–Tucker (KKT) conditions.  In this case, a particular interior point method known as barrier method is used to convert the inequality constraints to equality constraints. So, if $h_{i} \leq 0$ it can be written as $h_{i}+\alpha_{i}^{2}=0$ where $\alpha_{i}$ is a dummy input. If the problem is feasible, in the end, $h_{i}=-\alpha_{i}^{2}$ and for any given real value for $\alpha_{i}$, $h_{i} \leq 0$ The integrant part in cost function presented in Equation (2) has been also modified by adding a term $-r^{T}$ that is to give incentive to the solver to stay away from the boundaries of the feasible set [1,4].  Thus, the optimal control problem can be modified as:

$\begin{aligned} & \underset{u, \alpha}{\operatorname{minimize}} \cdot J=\varphi\left(x\left(t_{f}\right), p\left(t_{f}\right)\right) \\+& \int_{t_{0}}^{t_{f}}\left(L(x(\dot{t}), u(\dot{t}), p(\dot{t}))-r^{T} \alpha(\dot{t})\right) d \dot{t} \end{aligned}$ $\left\{\begin{aligned} \dot{x}(t)=& f(x(t), u(t), p(t)) \\ & x\left(t_{0}\right)=x_{o} \\ g(x(t), u(t), p(t)) &=0 \\ h(x(t), u(t), p(t))+\alpha(t)^{2} &=0 \end{aligned}\right.$ (3)

So, if the entries of are small positive numbers, the solution of this modified optimization problem is close to the solution of the original one with inequality constraints. Expanding the input vector u(t) to contain not only the real inputs but also the dummy inputs and considering  is the final equality constraints with dimension $m_{c}=m_{h} \cdot+m_{g}$ including the original equality constraints and the converted inequality constraints to equality constraints, the optimal problem can be written in more useful and generalized form as following [1,4]:

$\begin{aligned} \underset{u}{\operatorname{minimize} : J}=\varphi &(x(t+T), p(t+T)) \\ &+\int_{t}^{t+T} L(x(\dot{t}), u(\dot{t}), p(\dot{t})) d \dot{t} \end{aligned}$

Subjected to: $\left\{\begin{array}{c}{\dot{x}=f(x(t), u(t), p(t))} \\ {x(t)=x_{o}} \\ {C(x(t), u(t), p(t))=0}\end{array}\right.$ (4)

where, the optimal control $u^{*}(\dot{t} ; t, x(t))$ that minimize J is calculated over $\dot{t} \in[t, t+T]$, and its initial value at initial time t, $u^{*}(t ; t, x(t))$ is used as actual control input u(t) to the system.  The functions $\varphi$ and L

 are commonly selected for such case as:

$\varphi=\frac{1}{2}\left(x(t+T)-x_{f}\right)^{T} S_{f}\left(x(t+T)-x_{f}\right)$ (5)

$L=\frac{1}{2}\left(\left(x-x_{f}\right)^{T} Q\left(x-x_{f}\right)+u^{T} R \cdot u-r^{T} \alpha\right) )$ (6)

where Sf, Q, and R are free parameters' weighting matrices.

The NMPC problem is basically a family of finite horizon optimal control problems (FHOCPs) along fictitious time $\tau$ as follows [1,4]:

$\begin{aligned} \text { Minimize } & : J=\varphi\left(x^{*}(T, t), p(t+T)\right) \\ &+\int_{0}^{T} L\left(x^{*}(\tau, t), u^{*}(\tau, t), p(t+\tau)\right) d \tau \end{aligned}$ $\left\{\begin{array}{c}{x_{\tau}^{*}(\tau, t)=f\left(x^{*}(\tau, t), u^{*}(\tau, t), p(t+\tau)\right)} \\ {x^{*}(0, t)=x(t)} \\ {C\left(x^{*}(\tau, t), u^{*}(\tau, t), p(t+\tau)\right)=0 \cdots}\end{array}\right.$ (7)

where subscript $\tau$ represents partial differentiation with respect to $\tau$. The new state vector, $x^{*}(\tau, t)$, denotes a trajectory along the $\tau$ axis starting from $x(t) \cdot \operatorname{at} \tau=\cdot 0$. The optimal control input, $u^{*}(\tau, t)$, is obtained on the $\tau$ axis as the solution to FHOCP for each t, and actual control input is set by $u(t)=u^{*}(0, t)$. Horizon T is commonly a function of time, $T=T(t)$ .

To solve an optimization problem subject to equality constraints Lagrange multipliers method is used. The way it works is to write the Hamiltonian function as a sum of the integrant part of the cost function plus the dynamic equations multiplied by Lagrange multiplier $\lambda^{\prime} s$  known as costates plus the constraints equations multiplied by Lagrange multiplier $\boldsymbol{\mu}^{\prime} \boldsymbol{s}$ as seen in the following equation [1].

$\mathcal{H}(x, \lambda, u, \mu, p)=L(x, u, p)+\lambda^{T} f(x, u, p)+\mu^{T} C(x, u, p)$ (8)

Once the Hamiltonian function is defined, the first-order conditions necessary for optimal control is getting by calculus of variation as Euler Lagrange Equations (ELEs) [4].

$x_{\tau}^{*}=f\left(x^{*}, u^{*}, p\right)$ (9)

$x^{*}(0, t)=x(t)$ (10)

$\lambda_{\tau}^{*}=-\mathcal{H}_{x}^{T}\left(x^{*}, \lambda^{*}, u^{*}, \mu^{*}, p\right)$ (11)

$\lambda^{*}(T, t)=\varphi_{x}^{T}\left(x^{*}(T, t), p(t+T)\right)$ (12)

$\mathcal{H}_{u}\left(x^{*}, \lambda^{*}, u^{*}, \mu^{*}, p\right)=0$ (13)

$C\left(x^{*}, u^{*}, p\right)=0$ (14)

where the argument of p  is $t+\tau$ when absent. The ELEs describe a TPBVP, in which the value of initial state is known whereas the terminal costate is a function of the terminal state. Nonlinear TPBVP has to be solved within the sampling period for measured state x(t) at each sampling time, which represents a hard-computational load and is one of the major problems in NMPC. Now with discretizing these conditions for optimality on $\tau$-axis with forward difference by dividing the horizon into N time steps as follows [4].

$x_{i+1}^{*}(t)=x_{i}^{*}(t)+f\left(x_{i}^{*}(t), u_{i}^{*}(t), p(t+i \Delta \tau)\right) \Delta \tau$ (15)

$x_{0}^{*}(t)=x(t)$ (16)

$\lambda_{i}^{*}=\lambda_{i+1}^{*}+\mathcal{H}_{x}^{T}\left(x_{i}^{*}(t), \lambda_{i+1}^{*}(t), u_{i}^{*}(t), \mu_{i}^{*}(t), p(t+i \Delta \tau)\right) \Delta \tau$ (17)

$\lambda_{N}^{*}(t)=\varphi_{x}^{T}\left(x_{N}^{*}(t), p(t+T)\right)$ (18)

$\mathcal{H}_{u}^{T}\left(x_{i}^{*}(t), \lambda_{i+1}^{*}(t), u_{i}^{*}(t), \mu_{i}^{*}(t), p(t+i \Delta \tau)\right)=0$ (19)

$C\left(x_{i}^{*}(t), u_{i}^{*}(t), p(t+i \Delta \tau)\right)=0$ (20)

where time step $\Delta \tau=T / N$ and the real-time t is kept continuous while fictitious time $\mathrm{E} \tau$ is discretized. The computed sequences of the state, costate, input, and Lagrange multiplier on the discretized horizon are represented by $\left\{x_{i}^{*}(t)\right\}_{i=0}^{N}, \quad\left\{\lambda_{i}^{*}(t)\right\}_{i=0}^{N} \quad, \quad\left\{u_{i}^{*}(t)\right\}_{i=0}^{N} \quad$ and $\quad\left\{\mu_{i}^{*}(t)\right\}_{i=0}^{N}$, respectively.

Consequently, NMPC is formulated as TPBVP in Equation (15-20) for measured state x(t)  at real time t. The Equations (15,16) are the system dynamic equations with its initial conditions. The next two Equations (17,18) are dynamic equations corresponding to the Lagrange multiplier costates $\lambda$ with its final conditions. The rest of equations are equality constraints on all the variables that include the partial derivative of the Hamiltonian with respect to vector u, the vector of real and dummy inputs, and also the derivative of the Hamiltonian with respect to vector $\mu$ which lead to total equality constraints equation C=0 .

MPC is to solve an optimal control problem in each time step of a discrete time controller.  While the system is working the current states xn are measured or estimated then the optimal control problem is solved by computing the input ui for the next N steps and only the first step of the sequence is returned and that is the computed input which should be applied to the system then this loop is repeated again and again. Therefore, MPC is an implementation of an optimal controller in a real time.

3. C/GMRES Method

The C/GMRES method proposed by Ohtsuka in [1] allows for solving optimization problem in a real-time. The method consists of two steps. In the first step, the dynamic of the system Equations (15,16)  are solved forward in time and then the dynamic equations of the costates Equations (17,18)  are solved backward in time in order to get all the states xi and the costates $\lambda_{i}$ in terms of other variables $u_{i}^{*}$ and $u_{i}^{*}$. In the second step, all the remaining variables are put together in vector called U that includes the inputs u's , which contains the real inputs, dummy inputs and Lagrange multipliers $\mu^{\prime} s$ for the time step from 0  to N-1 as:

$U=\left[u_{0}^{T}, \ldots, u_{N-1}^{T}, \mu_{0}^{T}, \ldots, \mu_{N-1}^{T}\right]^{T}$ (21)

Then rewriting the remaining optimality conditions as:

$F(U(t), x(t), t)=0$ (22)


$F(U(t), x(t), t)=$ $\left[\begin{array}{c}{\mathcal{H}_{u}^{T}\left(x_{0}^{*}, \lambda_{1}^{*}, u_{0}^{*}, \mu_{0}^{*}, p(t)\right)} \\ {C\left(x_{0}^{*}, u_{0}^{*}, p(t)\right)} \\ {\cdot} \\ {\cdot} \\ {\mathcal{H}_{u}^{T}\left(x_{N-1}^{*}, \lambda_{N}^{*}, u_{N-1}^{*}, \mu_{N-1}^{*}, p(t+(N-1) \Delta \tau)\right)} \\ {C\left(x_{N-1}^{*}, u_{N-1}^{*}, p(t+(N-1) \Delta \tau)\right)}\end{array}\right]$ (23)

So, Equation (23) should be solved for U. The problem is that F is nonlinear function and this equation cannot be solved easily in general.  C/GMRES method is used to solve this problem. The continuation method is simply instead of solving F(U(0), x(0), t)=0, it selects  U(0) such that F(U(0), x(0), t)=0 and then find U so that [1]:

$\dot{F}(U(t), x(t), t)=-\xi F(U(t), x(t), t), \xi>0$ (24)

where $\dot{F}$ is the derivative of F with respect to time and the right-hand side is actually added to stabilize F=0. If Equation (24) is satisfied, regardless of where F is started which can be nonzero, it will converge to zero because this is the dynamics of the stable linear system. The derivative of F with respect to time can be computed using partial derivative with respect to x  and U  as:

$F_{x} \dot{x}+F_{U} \dot{U}+F_{t}=-\xi F(U, x, t)$ (25)

Equation (25) has a very interesting property which is its linearity in $\dot{U}$. Rewrite Equation (25) as:

$\dot{U}=F_{U}^{-1}\left(-\xi F-F_{x} \dot{x}-F_{t}\right)$ (26)

Given the current states x, the current value of vector U and the vector field of the system $f(x, u, p)$ or in other words $\dot{x}$,  everything can be computed except $\dot{U}$. Since Equation (26) is linear in $\dot{U}$, a linear solver to compute U(t) is needed. Once $\dot{U}$ is computed U(t) , update can be done to compute the new U and in turn, compute the control signal.

Unfortunately, Equation (26) still needs high computational operations represented in Jacobians FU, Fand Ft in addition to linear equation accompanied with $F_{U}^{-1}$. To reduce these computational loads, two main devices are employed. First the forward difference approximation for products of Jacobians and vectors, and second the GMRES method for the linear equation. The approximation for the products of the Jacobian and some $W \in R^{m N}, w \in R^{n}$ and $\omega \in R$ as

$F_{U}(U, x, t) W+F_{x}(U, x, t) w+F_{t}(U, x, t) \simeq$

$D_{h} F(U, x, t : W, w, \omega)=\frac{F\left(U+h_{d i r} W, x+h_{d i r} w, t+h \omega\right)-F(U, x, t)}{h_{d i r}}$ (27)

where hdir  is a positive real number. So, Equation (25) can be approximated using forward difference approximation method as [1,4]:

$D_{h} F(U, x, t : \dot{U}, \dot{x}, 1)=-\xi F(U, x, t)$ (28)

which can be written in an equivalent form as:

$D_{h} F(U, x+h \dot{x}, t+h : \dot{U}, 0,0)=b(U, x, \dot{x}, t)$ (29)


$b(U, x, \dot{x}, t)==-\xi F(U, x, t)-D_{h} F(U, x, t : 0, \dot{x}, 1)$ (30)

It is clear that the forward difference approximation of the products of the Jacobians and vectors can be estimated with only an additional calculation of the function, which needs particularly less computational load than approximation of the Jacobians themselves.

As Equation (29) approximates a linear equation in $\dot{U}$, a linear solver algorithm GMRES is applied with slight modification to include a certain initial guess of  $\hat{\dot{U}}$ [1]. The algorithm is employed with applying forward difference approximation so it is called FDGMRES [2] as follows:

Algorithm 1. $\dot{U} :=$ FDGMRES $\left(U, x, \dot{x}, t, \hat{U}, h, k_{\max }\right)$

(1) $\quad \hat{r} :=b(U, x, \dot{x}, t)-D_{h} F(U, x+h \dot{x}, t+h : \dot{U}, 0,0)$

$v_{1} :=\hat{r} /\|\hat{r}\|$

$\rho :=\|\hat{r}\|$

$\beta :=\rho$

$k :=0$

$\begin{array}{ll}{\text { (2) }} & {\text { While } k<k_{\max }, \mathrm{d} o} \\ {\text { a. }} & {k=k+1} \\ {\text { b. }} & {v_{k+1} :=D_{h} F\left(U, x+h \dot{x}, t+h : v_{k}, 0,0\right)} \\ {\text { c. }} & {\text { For } j=1 : k \quad h_{j k} :=v_{k+1}^{T} v_{j}} \\ {\text { d. }} & {v_{k+1} :=v_{k+1}-h_{j k} v_{j}}\end{array}$

$\begin{array}{ll}{e .} & {h_{k+1, k} :=\left\|v_{k+1}\right\|^{j}} \\ {f} & {v_{k+1} :=v_{k+1} /\left\|v_{k+1}\right\|}\end{array}$

$g . \quad$ For $e_{1}=[10 \ldots 0]^{T} \in \mathbb{R}^{k+1}$ and $H_{k}=\left(h_{i j}\right) \in$

$\mathbb{R}^{(k+1) \times k}$

$\left(h_{i j}=0 \text { for } i>j+1\right)$

Minimize $\left\|\beta e_{1}-H_{k} y^{k}\right\|$

Determine $y^{k} \in \mathbb{R}^{k} \rho :=\left\|\beta e_{1}-H_{k} y^{k}\right\|$

(3) $\quad \dot{U} :=\hat{U}+V_{k} y^{k}$ where $V_{k}=\left[v_{1} \ldots v_{k}\right] \in \mathbb{R}^{m N \times k}$

For a linear equation Ax=b with a nonsymmetric matrix A, the GMRES algorithm minimizes the residual $\rho :=\|b-A x\|$ at the kth  iteration with $x \in x_{0}+\mathcal{K}_{k}$, where x0 vector is the initial guess and $\mathcal{X}_{k}$ denotes the Krylov subspace expressed by $\mathcal{K}_{k} :=\operatorname{span}\left\{r_{0}, A r_{0}, \ldots, A^{k-1} r_{0}\right\}$ and $r_{0} :=b-A x_{0}$. GMRES creates an orthonormal basis $\left\{v_{j}\right\}_{j=1}^{k}$ for $\mathcal{X}_{k}$. Givens rotations is used for efficient minimization in step 2e shown in algorithm steps. Theoretically, GMRES has the ability to reduce the residual $\rho$ for a large linear equation monotonically and converges rapidly to an accurate solution with specified error tolerance within number of iterations equal to the dimension of the equation as maximum [1].

The product of the Jacobian matrix and a vector, $F_{U}(U, x, t) v_{j}$ in FDGMRES algorithm is computed using forward difference approximation, $D_{h} F\left(U, x+h \dot{x}, t+h : v_{k}, 0,0\right)$. Since the orthonormal basis$\left\{v_{j}\right\}_{j=1}^{k}, k$ vectors in $\mathbb{R}^{m N}$, must be saved during the implementation of FDGMRES, a massive amount of data storage may be required for a large-scale problem. Additionally, several iterations may need larger execution time which cannot be available in a real-time implementation. Therefore, the maximum number of iterations $k_{\max }$ should be carefully selected as small as possible to avoid excessive processing. Fortunately, both warm start at which the solution $\dot{U}$ at the previous sampling time is used as a consecutive initial guess for FDGMRES and selection of a small $k_{\max }$ often guarantee an accurate performance in a real-time application. On the other hand, in case of the need for a large number of iterations $k_{\max }$ while a limited amount of data storage is available, restart FDGMRES from the current iterate can efficiently handle this case. In restart technique, the orthonormal basis is reset with the latest obtained basis.

Once $\dot{U}$ is computed using FDGMRES algorithm, the update of the solution curve $U(t)$ of $F(U, x, t)=0$ is traced simply without iterative optimization methods by integrating this $\dot{U}$ with the Euler method in real time [6]. This approach is typical continuation method. Although the computational load of updating $U(t)$ for the explicit Euler technique related to only one iteration in Newton’s method, it realizes higher accuracy by taking the time dependence of the equation into account [4]. As the continuation method is combined with GMRES, this algorithm is named C /GMRES. The Ohtsuka’s suggestion of warm start in [1] which considers the solution for $\dot{U}$ in previous sampling time as the initial guess for current timestep's first inner iteration is used in this paper. The vector $\dot{x}$ in FDGMRES algorithm is computed through vector field equations $f(x(t), u(t), p(t))$ or approximated by finite difference in practical implementation as:

$\dot{x}(t) \approx(x(t+\Delta t)-x(t)) / \Delta t$ (31)

For designed sampling period $\Delta t$ and using the suffix(l) to express different variables at time $l \Delta t$, the C/GMRES algorithm for receding horizon control can be concluded as follows [1].

Algorithm $2 . \mathrm{C} / \mathrm{GMRES}$

(1) Initialization

$t :=0$

$l :=0$

Measure the initial state and set $x_{(0)} :=x(0)$

find $U_{(0)}$ analytically such that $\left\|F\left(U_{(0)}, x_{(0)}, 0\right)\right\| \leq \delta$ for

some positive $\delta$

(2) For $\dot{t} \in[t, t+\Delta t], \operatorname{set} U(\dot{t}) :=P_{0}\left(U_{(l)}\right)$

(3) At time $t+\Delta t$

set the state $x_{(l+1)} :=x(t+\Delta t)$

set $\Delta x_{(l)} :=x_{(l+1)}-x_{(l)}$

(4) compute $\dot{U}_{(l)} :=$ FDGMRES $({{U}_{(l)}},{{x}_{(l)}},\Delta {{x}_{(l)}}/\Delta t,t,{{\hat{\dot{U}}}_{(l)}},h,{{k}_{\max }})$ where the initial guess of ${{\hat{\dot{U}}}_{(l)}}$ is appropriately chosen; e.g., ${{\hat{\dot{U}}}_{(l)}}$:=0 or ${{\hat{\dot{U}}}_{(l)}}$:=$\dot{U}_{(l-1)}$ with ${{\hat{\dot{U}}}_{(-1)}}$:=0

(5) $\operatorname{Set} \Delta U_{(l)} :=\dot{U}_{(l)} \Delta t$

(6) Set $U_{(l+1)} :=U_{(l)}+\Delta U_{(l)}$

(7) Update $t :=t+\Delta t$$l :=l+1$

Go back to step 2.

It is clear from the C/GMRES algorithm that the usage of this iterative methodology is limited only in solving the linear Equation (26) with respect to $\dot{U}$. Then using this solution, the solution of the original nonlinear equation, $F(U, x, t)=0$, is traced easily with no line search or further Newton iteration. It means that the C/GMRES solves the linear Equation (26) only once per sampling time and, consequently, needs much less computational load than other standard iterative methods such as Newton’s method which needs to solve the linear equation many times to find search directions. Additionally, the C/GMRES algorithm doesn’t need line search, which is also a major difference from other regular optimization methods.

A simple way of initializing U(0) is to select a time-dependent horizon, T(t), as a smooth function satisfying $(T(0)=0)$ and $\left(T(t) \rightarrow T_{f}\right)$ at $(t=\infty)$ as [1]:

$T(t)=T_{f}\left(1-e^{-\alpha t}\right)$ (32)

where $T_{f}$ and $\alpha$ are positive numbers. Since the horizon has a zero length at (t=0), this leads to the following initial optimal parameters:

$\begin{array}{ll}{u_{i}^{*}(0)=u(0)} & {(i=0 : N-1)} \\ {\mu_{i}^{*}(0)=\mu(0)} & {(i=0 : N-1)} \\ {x_{i}^{*}(0)=x(0)} & {(i=0 : N)} \\ {\lambda_{i}^{*}(0)=\varphi_{x}^{T}(x(0))} & {(i=0 : N)}\end{array}$ (33)

The U(0) initialization reduces to finding only mu+mc unknown quantities, $u(0)$ and $\mu(0)$, such that [1]:

$\left\|\left[\begin{array}{c}{\mathcal{H}_{u}^{T}\left(x(0), \varphi_{x}^{T}(x(0)), u(0), \mu(0), p(0)\right)} \\ {C(x(0)), u(0), p(0) )}\end{array}\right]\right\| \leq \frac{\delta}{\sqrt{N}}$ (34)

where $\delta$ is a positive number. Equation (34) has a much-reduced size than the original Equation (23) and can be solved with Newton’s method off line or within the sampling period directly after initial state x(0) has been estimated.  The error in optimality, $\|F\|$, is guaranteed to be limited under several assumptions if the parameter $\xi$ is selected such that $(0<\xi \Delta t \leq \overline{\xi} \Delta t)$ for assumed sampling period $\Delta t$ and some $\bar{\xi }$ such that $(1<\overline{\xi} \Delta t \leq 2)$ [1].

Even if the problem of real-time optimization becomes a realistic option in implementing NMPC, another vital problem for practical application remains how to properly tune an NMPC performance index. Since it is problematic to have an explicit solution for a nonlinear optimal control problem, it is really difficult to relate the closed-loop response clearly to design parameters setting in a performance index such as weighting matrices and horizon length. Thus, tuning the performance index to realize a desired closed-loop performance is very effortful, not systematic and is considered as one of the major drawbacks of this approach.

An extensive search involving performance impact of numerous trials and modification of the free parameters is an alternative solution to overcome this drawback. This search is used here by randomly select the tuned parameters and hence their effect on the performance of the proposed NMPC controller is studied. This effect is used as a guide for optimum parameters selection.

4. Control Orinted Model

In fact, the efficient control performance of the NMPC scheme is always obstructed by the limited model precision. As for the complicated satellite dynamics with time-varying parameters and nonlinear characteristics, the used model may have some errors in addition to unexpected disturbance. Moreover, for the optimal tracking problem, the error dynamic model is crucial for the robust control design. Therefore, a new control-oriented model based on the attitude and rate errors is designed in this paper. Doing that allows handling the tracking task as a simple regulation task.

In this section, the main tracking problem is converted into a well-formulated global and simple regulation problem. The simple regulation problem is a more traceable problem and has a great potential dealing with disturbance and uncertainty. The attitude kinematic equation using quaternion error representation is described in [7] considering the inertial frame as a reference frame. The satellite is assumed to start GTT task from nadir pointing attitude and hence the orbital frame is used as a reference frame instead of the inertial frame for simplicity purposes. Thus, the attitude kinematic equation using quaternion error representation can be reformulated as:

$\delta \dot{\overline{q}}=0.5 \Omega(\delta \omega) \delta \overline{q}$ (35)


$\delta \omega=\omega_{S B C}^{O R C}-A_{d}^{S B C} \omega_{d}^{O R C}$ (36)

$\delta \overline{q}=\left[\begin{array}{c}{\boldsymbol{q}_{\boldsymbol{e}}} \\ {q_{4 e}}\end{array}\right]=\left[\begin{array}{cccc}{q_{4 d}} & {q_{3 d}} & {-q_{2 d}} & {-q_{1 d}} \\ {-q_{3 d}} & {q_{4 d}} & {q_{1 d}} & {-q_{2 d}} \\ {q_{2 d}} & {-q_{1 d}} & {q_{4 d}} & {-q_{3 d}} \\ {q_{1 d}} & {q_{2 d}} & {q_{3 d}} & {q_{4 d}}\end{array}\right]\left[\begin{array}{c}{q_{1}} \\ {q_{2}} \\ {q_{3}} \\ {q_{4}}\end{array}\right]$ (37)

The subscript d represents the desired frame, $\boldsymbol{\Omega}(\delta \omega)$ is 4*4 skew matrix function, the $\delta \overline{q}$  represents the quaternion error from the desired to the body frame, $\omega_{S B C}^{O R C}$ represents desired quaternion, q represents the current attitude, $\delta \omega$ is the error in angular rate represented in the body frame, $\omega_{S B C}^{O R C}$ is the orbital referenced body angular rate in body frame, $\omega_{d}^{O R C}$ is the orbital referenced desired angular rate in the desired frame and $A_{d}^{S B C}$ is the transformation matrix from desired to body frame.

The dynamic model of an Earth-pointing satellite using reaction wheels (RWs) only as active actuators is given by:

$J \dot{\omega}_{S B C}^{O R C}=N_{\omega o}+N_{G G}-\omega_{S B C}^{E I C} \times\left(J \omega_{S B C}^{E I C}+h_{w}\right)-\dot{h}_{w}$ (38)

Using Equation (36) to express the body angular rate $\omega_{S B C}^{O R C}$ as:

$\omega_{S B C}^{O R C}=\delta \omega+A_{d}^{S B C} \omega_{d}^{O R C}$ (39)

Differentiating Equation (39) to get

$\dot{\omega}_{B}^{O R C}=\delta \dot{\omega}-\Omega(\delta \omega) A_{d}^{S B C} \omega_{d}^{O R C}+A_{d}^{S B C} \dot{\omega}_{d}^{O R C}$ (40)

where $\Omega(\delta \omega)$ is $3 \times 3$ is skew matrix. Multiplying Equation (40) by J leads to

$J \dot{\omega}_{S B C}^{O R C}=J \delta \dot{\omega}-J\left[\Omega(\delta \omega) A_{d}^{S B C} \omega_{d}^{O R C}-A_{d}^{B} \dot{\omega}_{d}^{O R C}\right]$ (41)

Substituting $J \dot{\omega}_{B}^{O R C}$ from Equation (41) in Equation (38) gives:

$\begin{aligned} J \delta \dot{\omega}=& N_{\omega o}+N_{G G}+J &[ & \Omega(\delta \omega) A_{d}^{S B C} \omega_{d}^{O R C}-\\ & A_{d}^{S B C} \dot{\omega}_{d}^{O R C} ]-\omega_{S B C}^{E I C} \times\left(J \omega_{S B C}^{E I C}+h_{w}\right)-\dot{h}_{w} \end{aligned}$ (42)

Generally, the terms NGG  and $N_{\omega O}$ and $\omega_{S B C}^{E I C} \times\left(J \omega_{S B C}^{E I C}+h_{w}\right)$, are not significant for many practical rotation maneuvers [8]. The term $\dot{\omega}_{d}^{O R C}$ can also be neglected as it is lower several orders of magnitude than other terms for nearly circular orbit. All the omitted terms can be considered as some types of uncertainty which can be handled effectively using such presented model. Robustness check in the simulation is also used to confirm the validity of NMPC with the presented model in general. Rewriting Equation (42) in more simple and efficient form after neglecting the smaller terms yields:

$J \delta \dot{\omega}=J \Omega(\delta \omega) A_{d}^{S B C} \omega_{d}^{O R C}+u_{w}$

$\dot{h}_{w}=-u_{w}$ (43)

where uw  represents the RWs torque. Equation (35) and Equation (43) are the control-oriented model or in other words, the tracking error dynamic model used in this paper. This model has 10 states represented in attitude quaternion error, rate error, and RWs momentum. The model also includes time-varying parameter represented in $A_{d}^{S B C} \omega_{d}^{O R C}$ which add a little bit difficulty in implementation. The desired rate $\omega_{d}^{O R C}$ is [8]:

$\omega_{d}^{O R C}=u_{S / T}^{O R C} \times \dot{u}_{S / T}^{O R C}$ (44)

where $u_{S / T}^{O R C}$ is unit direction vector of the direction from satellite to the target represented in orbit frame, $X_{S / T}^{O R C}$ computed by [8]:

$X_{S / T}^{O R C}=A_{E I C}^{O R C}\left(A_{E F C}^{E I C} X_{T}^{E F C}-r_{s a t}\right)$ (45)

as $r_{s a t}$ is the satellite position in Earth-centered inertial frame, $X_{T}^{E F C}$ is the target location in Earth-centered fixed frame and the transformation matrix $A_{E F C}^{E I C}$ from the fixed frame to the inertial is time-dependent matrix can be estimated considering constant Earth angular rate around its rotation axis, $\omega_{E}$ in the tracking maneuver period and knowing the initial phase , $\alpha_{G}$ between the x-axes of both EFC and EIC frames at t=0 as:

$A_{E F C}^{E I C}=\left[\begin{array}{ccc}{\cos \left(\omega_{E} t+\alpha_{G}\right)} & {-\sin \left(\omega_{E} t+\alpha_{G}\right)} & {0} \\ {\sin \left(\omega_{E} t+\alpha_{G}\right)} & {\cos \left(\omega_{E} t+\alpha_{G}\right)} & {0} \\ {0} & {0} & {1}\end{array}\right]$ (46)

The transformation matrix $A_{E I C}^{O R C}$ from the inertial frame to orbital frame can be deduced by using the satellite position rsat  and velocity vsat  in the inertial frame as [8]:

$A_{E I C}^{O R C}=\left[\begin{array}{lll}{\hat{u}} & {\hat{v}} & {\widehat{w}}\end{array}\right]^{T}$ (47)

where the unit vectors $\widehat{u}$,  $\hat{v}$ and $\widehat{w}$ are:

$\begin{aligned} \widehat{w} &=-\frac{r_{s a t}}{\left\|r_{s a t}\right\|} \\ \hat{v}=&-\frac{r_{s a t} \times s a t}{\left\|r_{s a t} \times v_{s a t}\right\|} \\ \hat{u} &=\hat{v} \times \widehat{w} \end{aligned} \}$ (48)

The vector $\dot{u}_{S / T}^{O R C}$ can be calculated as following [8]:

$\dot{u}_{S / T}^{O R C}=\frac{1}{\left\|X_{S / T}^{O R C}\right\|}\left[I_{3}-u_{S / T}^{O R C}\left(u_{S / T}^{O R C}\right)^{T}\right] \dot{X}_{S / T}^{O R C}$ (49)

where the time derivative $\dot{X}_{S / T}^{O R C}$  is calculated by differentiating Equation (45) to get the following [8]:

$\dot{X}_{S / T}^{O R C}=\dot{A}_{E I C}^{O R C}\left(A_{E F C}^{E I C} X_{T}^{E F C}-r_{s a t}\right)+A_{E I C}^{O R C}\left(\dot{A}_{E F C}^{E I C} X_{T}^{E F C}-\right.$

$v_{s a t} )$ (50)

The time derivative of transformation matrices $A_{EIC}^{ORC}$ and $A_{\text{EFC}}^{\text{EIC}}$ in Equation (50) are calculated as following [8]:

$\dot{A}_{EIC}^{ORC}={{[\hat{\dot{u}}\text{ \hat{\dot{v}} \hat{\dot{w}}}]}^{T}}$ (51)


$\left. \begin{matrix}   \hat{\dot{w}}=-\frac{1}{\left\| {{r}_{sat}} \right\|}({{I}_{3}}-\bar{w}{{{\bar{w}}}^{T}}){{v}_{sat}} \\\hat{\dot{v}}=0\\   \hat{\dot{u}}=\frac{1}{\left\| {{r}_{sat}} \right\|}\hat{v}\times [({{I}_{3}}-\bar{w}{{{\bar{w}}}^{T}}){{v}_{sat}}]  \\\end{matrix} \right\}$  (52)

\[\dot{A}_{EFC}^{EIC}={{\omega }_{E}}\left[ \begin{matrix}   -\sin ({{\omega }_{E}}t+{{\alpha }_{G}})-\cos ({{\omega }_{E}}t+{{\alpha }_{G}})\text{ 0}  \\   \cos ({{\omega }_{E}}t+{{\alpha }_{G}})-\sin ({{\omega }_{E}}t+{{\alpha }_{G}})\text{   0}  \\   0\text{  0   0}  \\\end{matrix} \right]\] (53)

The transformation matrix  AdSBC in Equation (43) can be extracted easily from the quaternion error $\delta \text{\bar{q}}$ given in Equation (37). The desired quaternion qd can be estimated as following [8]:

${{q}_{d}}=\left[ \begin{matrix}{{u}_{c}}\sin \frac{\delta }{2} \\\cos \frac{\delta }{2}  \\\end{matrix} \right]$ (54)

where δ is the angle between the time-dependent vector $u_{S/T}^{ORC}$ and commissioned axis or payload axis represented in body frame $u_{com}^{\text{SB}C}$ while uc define the axis of rotation needed for constructing qd given by [8]:

${{u}_{c}}=\frac{u_{com}^{SBC}\times u_{S/T}^{ORC}}{\left\| u_{com}^{SBC}\times u_{S/T}^{ORC} \right\|}$ (55)

The tracking error can be examined by calculating θerr  to define the angle between $u_{com}^{\text{SB}C}$ and $u_{S/T}^{ORC}$ vectors as [8]:

${{\theta }_{err}}={{\cos }^{-1}}(u_{S/T}^{SBC}\centerdot u_{S/T}^{ORC})$ (56)

where the unit vector $u_{S/T}^{\text{SB}C}$ defined for the current attitude matrix $A_{ORC}^{\text{SB}C}$ as:

$u_{S/T}^{SBC}=A_{ORC}^{SBC}\centerdot u_{S/T}^{ORC}$ (57)

5. Image Quality Requirements

The image quality enhancement is directly linked to accurate and low distortion images. To obtain low distortion images during tracking mode, the shooting conditions should guarantee stability requirements. Moreover, to ensure accurate shooting of definitely targeted area, minimum permissible geo-location error should be kept. For a given GSD (i.e. one pixel) capability of a certain payload its corresponding angular distance θpixel  should be obtained as [9,10].

${{\theta }_{pixel}}={{\tan }^{-1}}(\frac{GSD}{h})$ (58)


 is orbit altitude. The exposure time texp  of a single pixel (or a single row of pixels) is a function of the effective ground speed vgeff  as follows [9]:

${{t}_{\exp }}=\frac{GSD}{{{v}_{geff}}}$ (59)

${{v}_{geff}}=\frac{{{v}_{g}}}{{{n}_{fmc}}}$ (60)

As nfmc is the forward motion compensation (FMC) factor; reducing the effective ground speed of the image sensor boresight, vg  is the satellite ground speed defined, μE is the Earth gravitational constant and RE is the Earth radius.

The permissible safe variation in angular rate around the image sensor's boresight,ωsz, during shooting can be determined [9]:

${{\omega }_{sz}}\approx \frac{{{\tan }_{-1}}(\frac{\max \_smear\times GSD}{0.5\times {{L}_{swath}}})}{{{t}_{\exp }}}$ (61)

Where Lswath  is the image swath width while max_smear is the maximum pixel smear in percentage. The permissible safe variation in angular rate across the image sensor’s boresight, ωsxy, can similarly be computed via [9]:

${{\omega }_{sxy}}\approx \frac{\max \_smear\times GSD}{{{t}_{\exp }}}$ (62)

The stability requirement is selected to be the minimum of the calculated angular rate variation around the image sensor’s boresight,ωsz, and the calculated angular rate variation across the image sensor’s boresight,ωsxy.On the other hand, the payload integration time during which data is collected by the image sensor is another issue that directly affects the image quality. The higher the integration time, the more data can be acquired and consequently the more enhancement in the image quality.

6. Case Study

In order to evaluate the performance of using algorithm C/GMRES in solving a ground target tracking problem, the configuration UoSAT-12 satellite with higher RWs capabilities is used in the simulation. The inertia matrix J of the satellite in kgm2 is diag(40,40,32).

Table 1.  Controller and simulation parameters

Setting parameters





Dimension of state vector



Number of real input (actuators)



Number of constraints



Dimension of control input vector



Dimension of Time-Variant Parameters



final horizon length (sec.)



For defining variable prediction horizon



For Stabilization of Continuation Method



Step in forward difference approximation



Error tolerance in control input vector



Number of iteration in GMRES



Number of grids on the horizon


r1 to r6

For keeping boundary limit of feasible set



Simulation time (sec.)



Simulation step (sec.)



time step for saving data (sec.)


Weighting matrices and initial input vector


Description and setting


Initial guess for initial control input vector



Weighting values for terminal state



Weighting values for trajectory impact

diag(50,50,50,50,50,50,50,0.01,0.01, 0.01)


Weighting values for control input


Three reaction wheels installed in a three-axis configuration to allow full control of the satellite attitude orientation with maximum torque $\pm 0.2\text{Nm}$ and momentum capacity $\pm 6\text{Nmsec}$. The working orbit is almost circular with altitude 650km,  inclination angle 64.5deg, right ascension of ascending node 10deg, initial mean anomaly  0deg and initial mean phase 0deg. The tracked fixed ground-target in EFC frame is [4021.9  -35.1   4933.6]The satellite commissioned axis is [1  -1  9]and the satellite attitude is assumed initially aligned with orbit frame.

During simulation perfect attitude knowledge and measurements of the vectors rsat  and vsat are assumed. The developed code of C/GMRES method is outputted as an ANSI/C standard source code for simulation purposes in this paper. In Table  1 a list of all the controller variable and simulation parameters used during simulation for GTT task.

Since the clear desired states for reformulated regulation problem represented in the attitude quaternion and angular rate, it is planned to lightly punish the deviation between the desired zero angular momentum and realized momentum during optimum solution searching.

Actually, the restrict requirements concern the trajectory of angular momentum change is considered in the limits put on the input vector real elements (the first three elements in the input vector). All parameters which have no real physical limits are allowed to be changed freely that is why their corresponding weighting values are relatively too low.

UoSAT -12 conveys two redundant 32m ground sampling distance (GSD) multi-spectral cameras to provide 64Km wide swath width [11]. Furthermore, UoSAT-12 also carries a high-resolution narrow-angle panchromatic camera with 10m GSD and 10Km swath [11]. The stability requirement during imaging to obtain acceptable image quality with relaxed maximum pixel smear in percentage 17% as explained in section (‎5) is computed and found to be better than  0.1 deg/sec.  It can be realized from Figure 1 and Figure 2 that the satellite needs only (≈ 200 sec) to be kept in a stable tracking mode with a very small tracking error defined in Equation (56), ${{\theta }_{err}}\approx 0.0045\deg $.

The supreme advantage here is that the proposed controller succeeded to maintain this precise tracking error during the overpass flight simulation time. Normally this accuracy is deteriorated whenever the satellite position becomes closer to the target. The resultant bell-shaped figure with maximum satellite angular slew rate at a closer position to the target violate the stability conditions needed for shooting with accepted image quality. The proposed tracking controller here can react as fast as the kinematic of the target. Keeping the tracking error within this precise limit during the overpass flight of the target allows shooting this target at any time. Hence, an additional shooting mode covering stereo triplets imaging and video surveillance can be carried out.

Figure 1. Tracking error

Figure 2. Tracking error zoom in

The Euler angles representing the satellite attitude change related to the desired frame is shown in Figure 3. The initial attitude mismatch is diminished fast.

Figure 3. Current attitude in Euler angles

The angular rate error of the satellite body is shown in Figure 4 during the tracking control. It is seen that the that the orbit-referenced angular velocity of the satellite body $\omega _{B}^{ORC}$ follows the desired angular rate command $\omega _{\text{d}}^{ORC}$ without a noticeable error during the tracking period.

The stability requirement (better than 0.1deg/sec) during imaging to obtain acceptable image quality is initially guaranteed after only 57sec as cleared from Figure 5. Fortunately, this stability guarantee is kept through over the target overpass flight simulation time.  From  Figure 2 and Figure 5,  the shooting conditions with minimum geo-location error is guaranteed after 200sec of starting tracking process with tracking error 0.0045 deg.

Figure 4. Satellite angular rate error

Figure 5. Current angular rate

The geo-location error due to this error will be 0.05km. Therefore, the high-resolution narrow-angle panchromatic camera which offers 10Km swath (i.e. image size 10*10) can be used efficiently while tracking with the proposed intelligent tracking controller. The steady-state quaternion error qe=[0 0 0 1]T  represented in Figure 6 means that body attitude is totally aligned with the desired attitude.

Figure 6. Quaternion error

Since the UoSAT-12 satellite is equipped with body-fixed imager, aligning the body attitude with the desired attitude means that this imager (payload) pointing at the target.

The activity of the three-axis RWs is shown in Figure 7  and Figure 8. Figure 7 confirms that the torques of the three-axis RWs are generally far from saturation limits (0.2Nm) during the duration needed for tracking the ground target. A wheel torque transient, in general, happens at the beginning of the tracking maneuver due to initial attitude error between the satellite payload axis and the desired target direction. The max torque provided by the RW of Y-axis and it is almost (0.19Nm) with very limited duration time. Once the controller starts to respond to the desired command, this attitude error starts to be reduced and hence the torques needed from RWs are decreased.

The angular momentums of the three-axis RWs in Figure 8 are far from angular momentum capacity hwmax  (6.0Nm sec) during tracking task with the proposed controller. This safe working momentum region allows the RWs to further compensate for any disturbance easily and increase the robustness of the system as whole. It is also noticed that for this specific target and selected commissioned axis the effort needed from RW of Z-axis is minor however the other RWs provide the main effort for carrying out the required tracking maneuver. Since the disturbance is not considered during simulation, it is logic to note that the angular momentums of the RWs are kept unchangeable after achieving the tracking task.

Figure 7. Wheel torque

Figure 8. Wheel momenta

Figure 9. Robustness check against uncertainty in MOI

Because of the feedback nature of the proposed C/GMRES controller, the robustness of this controller against model uncertainties is guaranteed. Simulations were also done to examine the robust performance of the proposed tracking controller against a ±20% error in the moment of inertia of the satellite. These simulations demonstrate that the precise tracking error can still be preserved to designed values ± 0.003 degree as shown in Figure 9. However, in practice, the TE will ultimately be impacted by the attitude determination accuracy as well as the orbital measurements of the position and velocity vectors from the GPS receiver which have been not considered in this study.

7. Conclusion

In this paper, a feasible solution to enhance the UoSAT-12 satellite shooting capabilities in GTT task is proposed. Employing the C/GMRES, as fast RTO methods of NMPC approach, along with a proposed control-oriented model, that converts the main tracking problem into a well-formulated simple regulation problem, is presented in details.  The simulation results show that proposed algorithm succeeds to track a fixed-ground target robustly through a large angle maneuver with a precise tracking error of 0.0045 deg, while maintaining the stability requirements within 0.1deg/sec throughout the overpass flight. Therefore, the shooting mode capabilities are significantly enhanced.


[1] Ohtsukan T. (2004). A continuation/GMRES method for fast computation of nonlinear receding horizon control. Journal of Automatica 40 (4): 563–574. 10.1016/j.automatica.2003.11.005

[2] Kelley CT. (1995). GMRES Iteration, Iterative methods for linear and nonlinear equations, (16), Philadelphia, Pa.: Society for Industrial and Applied Mathematics, North Carolina, pp. 33 – 60.

[3] Kouramas KI, Panos C, Faisca NP, Pistikopoulos EN. (2013). An algorithm for robust explicit/multi-parametric model predictive control. Journal of Automatica 49 (2): 381-389. 10.1016/j.automatica.2012.11.035

[4] Ohtsuka T. (2015). A tutorial on C/GMRES and automatic code generation for nonlinear model. European Control Conference, Linz, Austria, pp. 73-86. 10.1109/ECC.2015.7330528

[5] Tajeddin S, Vajedi M, Azad NL. (2016). A Newton/GMRES approach to predictive ecological adaptive cruise control of a plug-in hybrid electric vehicle in car-following scenarios, International Federation of Automatic Control 49 (21): 059–065. 

[6] Richter SL,  Decarlo RA. (1983). Continuation methods: Theory and applications, IEE. Transactions on Circuits and Systems 28(6): 660–665.

[7] Weiss H. (1993). Quaternion-based rate/attitude tracking system with application to gimbal attitude control. Journal of Guidance, Control, and Dynamics 16 (4) : 609-616. 10.2514/3.21057

[8] Chen XJ, Willem HS, Hashida Y. (2000). Ground-target tracking control of earth-pointing satellites. AIAA Guidance, Navigation And Control Conference, Denver, USA. 10.2514/6.2000-4547

[9] Vuuren GH, Willem SH. (2015). The design and simulation analysis of an attitude determination and control system for a small earth satellite. M.S thesis, Department of Electrical and Electronic Engineering, Stellenbosch University, Stellenbosch.

[10] Wahballah W, Bazan TM, El-tohamy F, Fathy M. (2016). Analysis of smear in high-resolution remote sensing satellites. SPIE, Edinburgh. UK, 1-11.

[11] Sun W. (2000). In-Orbit Results From UOSAT-12 Earth Observation Minisatellite Mission, IAA-B3-0305P, Guildford, Surrey, 1-4.