**3. Example of PIL simulation application**

As an example of PIL simulation scheme, we present the application of this strategy on the attitude control problem of an artificial satellite. This case study consists of an embedded digital attitude control system (ACS) design for three axis stabilization [7]. The information obtained by the sensors is processed by a Digital Signal Processor (DSP). The controller design is based on Linear-quadratic Gaussian (LQG) optimal control approach, a well known control theory in terms of space technology applications.

### **3.1. The satellite attitude control problem**

executable code for the target computer from the model built in Simulink. The goal is to achieve simulation and real-time implementation of an algorithm integrating inertial navigation and GPS, using the validation and testing of the proposed algorithm in the FlightGear flight

In order to exemplify the possibility of use of different development environment to perform SIL and PIL simulations of a satellite control system; in the work of [6], a SIL simulation is made in MATRIXx environment where is obtained and tested the control algorithm, simulated the space environment and the satellite dynamics. This same software is used to generate the controller code, that after is downloaded in the satellite processor while its model is ported to a DSP that is controlled by a computer. This computer runs the real-time simulation together with the satellite processor, controlling the actuators and sensors signals too, in a PIL simula‐ tion. An interface allows the user to interact with the simulation, monitoring and changing

In the cited works, is observed as the simulations SIL, PIL and HIL can be useful in controllers project and test to aerospace applications and as this simulations can be implemented in different ways, different virtual development environments and together with several hardware platforms. The results are not only obtained much more quickly but come from tests very cheap that those made in a real platform, which is often not accessible to the researchers. Furthermore, the interaction with the user in real-time and the analysis instruments, available by co-simulation tools, make the tests much more dynamic and enables the analysis of the

Today there are many co-simulation tools and many are the possibilities of combination between these tools and the several processors hardware platforms. Such tools are in contin‐ uous development and making possible the coupling with several hardware platforms from different manufacturers, including additional functional facilities to become easier the task of

In these studies, the Simulink is used as a friendly graphical tool for systems development, design and generation of executable code for various devices and applications. The following section presents more details of a control system design procedure using the PIL approach,

As an example of PIL simulation scheme, we present the application of this strategy on the attitude control problem of an artificial satellite. This case study consists of an embedded digital attitude control system (ACS) design for three axis stabilization [7]. The information obtained by the sensors is processed by a Digital Signal Processor (DSP). The controller design is based on Linear-quadratic Gaussian (LQG) optimal control approach, a well known control

simulator (inside MATLAB), running on the host computer.

simulation parameters in real-time.

160 Computational and Numerical Simulations

system response in situations more realistic.

Simulink tools and a DSP processor.

**3. Example of PIL simulation application**

theory in terms of space technology applications.

development.

The mathematical model of the satellite attitude considers two reference systems: the orbital frame and the satellite body frame. The orbital frame moves jointly with the satellite. Its origin coincides with the satellite's center of the mass. The axis *z*<sup>0</sup> of the orbital frame is defined by the satellite radius vector, the axis *y*<sup>0</sup> by the orbital normal, and the axis *x*<sup>0</sup> completes a right handed coordinate frame.

The body frame (*x*, *y*, *z*) has its origin in the center of mass of the satellite. Their axes coincide with the satellite's principal axes of inertia. We consider the case where the nominal attitude has axes of the body frame aligned with the axes of the orbital frame [8]. The Euler angles are adopted as parameters of attitude, considering the sequence of rotations 3-2-1. In this case, the rotation matrix is given by [9].

$$R\_x^X = \begin{bmatrix} \cos\psi\cos\theta & \cos\psi\cos\theta & -\sin\theta \\ -\cos\phi\sin\psi + \sin\phi\sin\theta\cos\psi & \cos\phi\cos\psi + \sin\phi\sin\theta\sin\psi & \sin\phi\cos\theta \\ \sin\phi\sin\psi + \cos\phi\sin\theta\cos\psi & -\sin\phi\cos\psi + \cos\phi\sin\theta\sin\psi & \cos\phi\cos\theta \end{bmatrix} \tag{1}$$

The kinematics equation is obtained by the time derivative of the rotation matrix equation [8], then:

$$
\dot{R}(t) = S(\omega(t))R(t) \tag{2}
$$

The kinematics equation can be simplified if one considers small angle maneuvers, in this case we can approximate: sin*ϕ* ≈*ϕ*, sin*θ* ≈*θ*, sin*ψ*≈*ψ*, cos*ϕ* ≈1, cos*θ* ≈1 and cos*ψ*≈1. Moreover, the non-linear terms (*ψ*˙ *θ*, *ψ*˙ *ϕ*, *θ*˙*ϕ*) are very small compared to the linear ones. The velocities are also small compared to the orbital velocity [10]. Considering those approximations the kinematic equation is given by:

$$
\omega\_{ib} = \begin{bmatrix} \phi \\ \theta \\ \psi \end{bmatrix} + \omega\_0 \begin{bmatrix} -\psi \\ -1 \\ 0 \end{bmatrix} \tag{3}
$$

In this work the satellite is modeled as a rigid body. Therefore the dynamic model can be obtained with the Euler equation that describes the rotation of a rigid body [9, 11]. The dynamic equation is given by:

$$J\,\dot{\omega}^{b}\_{\dot{\imath}b} + S\,(\omega^{b}\_{\dot{\imath}b})J\,\omega^{b}\_{\dot{\imath}b} = \tau^{b}\_{d} + \tau^{b}\_{p} \tag{4}$$

where *J* is the inertia matrix of the satellite, *τ<sup>d</sup> b* represents the torques from external perturba‐ tions acting on the satellite, and *τ<sup>p</sup> <sup>b</sup>* represent the control torques (*τx*, *<sup>τ</sup>y*, *<sup>τ</sup>z*), all vectors described in the body frame, *ωib <sup>b</sup>* is the angular velocity of the body with respect to the inertial frame written in the body frame. The gravity gradient is the only external perturbation considered in this work. Its effect cannot be neglected in a low-orbit satellite; an asymmetric body subject to the Earth gravitational field will experience a torque tending to align the axis of minor inertia with the field direction. The gravity gradient torque is modeled as:

$$\begin{aligned} \tau\_{\mathcal{S}}^b = 3\alpha\_0^2 \begin{bmatrix} (J\_z - J\_y)\phi \\ (J\_x - J\_z)\theta \\ 0 \end{bmatrix} \end{aligned} \tag{5}$$

Due to the presence of noises in sensors measurements and uncertainties in models, a filter is needed to obtain more reliable information about the states. For the inclusion of signals filtering, we adopted the controller design based on the theory of Linear Quadratic Gaussian control [12]. We consider in this work the presence of white noise in the observations, and that the system is observable. In the LQG problem, we want to minimize the cost function:

Processor-in-the-Loop Simulations Applied to the Design and Evaluation of a Satellite Attitude Control

( )

where *Qf* is the covariance matrix of the measurements noise, and *Rf* is the covariance matrix of the dynamics noise (or model uncertainties). The LQG control law determination will be given by solving two problems: setting the controller for the linear quadratic deterministic

^ <sup>+</sup> *Bu* <sup>+</sup> *<sup>L</sup>* (*<sup>y</sup>* <sup>−</sup>*Cx*

where *L* is the Kalman filter gain that is obtained by solving the algebraic Riccati matrix

−1

The linear quadratic controller action can be implemented using actuators such as reaction wheels and magnetic actuators for a continuous control command. However, during orbital operations, such as rapid detumbling maneuvers, the required torqueses are usually too high for reaction wheels. Therefore, on-off propulsion strategies are used for such operations [9]. The choice of cold gas jets actuation in the present study aims to test the control in a most

The firing of the gas jets is controlled by a pulse-width pulse-frequency modulator (PWPF) [13]. The PWPF is an interesting option for the thrusters control system due to its advantages over other types of pulse modulators. PWPF is designed to provide propulsion output proportional to the input command. The modulator optimizes the use of propellants; it

*B* is the transfer function of the attitude dynamics.

*lqg f f J x Q x u R u dt* (10)

<sup>1</sup> <sup>0</sup> - +- + = *T T A S SA SCR C S Q f f* (12)

<sup>1</sup> () ( ) () - *K G s K sI A BK LC LG s lqg* = -+ + (13)

^ of state *<sup>x</sup>*. The formulation

http://dx.doi.org/10.5772/57219

163

^) (11)

*C <sup>T</sup> S*. After obtaining *L*, it is possible to obtain

*T T*

0 = + ò *T*

problem, and setting a Kalman-Bucy filter for optimum estimation *x*

*x* ^˙ <sup>=</sup> *Ax*

the transfer function of open loop LQG controller according to [10]:

critical situation in terms of the small attitude adjustments difficulties.

of the Kalman-Bucy filter is given by:

The gain of the optimal filter is given by *L* =*Rf*

−1

equation:

where *G*(*s*)=*C*(*s*)(*sI* − *A*)

Substituting the applied torques and the kinematic equation (Eqs. 5 and 3) into Eq. 4, and representing this dynamics in the state space form, we have [10]:

$$\begin{aligned} \dot{x} &= A\chi + Bu\\ y &= C\chi + Du \end{aligned} \tag{6}$$

### *3.1.1. The LQG control approach*

The Linear Quadratic Regulator is designed upon the linearization of the dynamic model. The theory is developed for linear systems. However, the system simulation takes into considera‐ tion the complete non-linear model of the satellite. This optimal control approach consists of minimizing a quadratic cost function and computing a feedback gain matrix [12]. The optimi‐ zation problem aims to obtain a control law expressed by a linear relationship between the state variable *x* (expressing here the Euler's angles and the angular velocities) and the control variable *u* (the applied torques), i.e. *u* = − *Kx*. The matrix of gain *K* is obtained by the minimi‐ zation of a quadratic cost function formulated as follows:

$$J\_{lqr} = \bigcap\_{0}^{T} \left( \mathbf{x}^{T} \underline{Q}\_{c} \mathbf{x} + \mathbf{u}^{T} R\_{c} \mathbf{u} \right) dt \tag{7}$$

where *Qc* and *Rc* are the weight matrices of state and control vectors, respectively. The value of *K* results from solving the algebraic matrix Riccati equation for a time-invariant system and considering an infinite horizon context:

$$A^T P + PA - PBR\_c^{-1}B^T P + Q\_c = 0\tag{8}$$

where *A* and *B* are the matrices of the linearized attitude dynamic model. The optimal control gain is then given in terms of the solution *P* of the Riccati equation:

$$K = R^{-1}B^T P \tag{9}$$

Due to the presence of noises in sensors measurements and uncertainties in models, a filter is needed to obtain more reliable information about the states. For the inclusion of signals filtering, we adopted the controller design based on the theory of Linear Quadratic Gaussian control [12]. We consider in this work the presence of white noise in the observations, and that the system is observable. In the LQG problem, we want to minimize the cost function:

described in the body frame, *ωib*

162 Computational and Numerical Simulations

*3.1.1. The LQG control approach*

*<sup>b</sup>* is the angular velocity of the body with respect to the inertial

*J J* (5)

(6)

frame written in the body frame. The gravity gradient is the only external perturbation considered in this work. Its effect cannot be neglected in a low-orbit satellite; an asymmetric body subject to the Earth gravitational field will experience a torque tending to align the axis

> ( ) 3( ) 0

Substituting the applied torques and the kinematic equation (Eqs. 5 and 3) into Eq. 4, and

The Linear Quadratic Regulator is designed upon the linearization of the dynamic model. The theory is developed for linear systems. However, the system simulation takes into considera‐ tion the complete non-linear model of the satellite. This optimal control approach consists of minimizing a quadratic cost function and computing a feedback gain matrix [12]. The optimi‐ zation problem aims to obtain a control law expressed by a linear relationship between the state variable *x* (expressing here the Euler's angles and the angular velocities) and the control variable *u* (the applied torques), i.e. *u* = − *Kx*. The matrix of gain *K* is obtained by the minimi‐

( )

where *Qc* and *Rc* are the weight matrices of state and control vectors, respectively. The value of *K* results from solving the algebraic matrix Riccati equation for a time-invariant system and

where *A* and *B* are the matrices of the linearized attitude dynamic model. The optimal control

*lqr c c J x Q x u R u dt* (7)

<sup>1</sup> <sup>0</sup> - + - += *T T A P PA PBR B P Q c c* (8)


*T T*

0 = + ò *T*

gain is then given in terms of the solution *P* of the Riccati equation:

*J J*

é ù ê ú = - ê ú ê ú ë û

*z y*

f

 q

of minor inertia with the field direction. The gravity gradient torque is modeled as:

2 0

*g xz*

*x*˙ = *Ax* + *Bu y* =*Cx* + *Du*

tw

*b*

representing this dynamics in the state space form, we have [10]:

zation of a quadratic cost function formulated as follows:

considering an infinite horizon context:

$$J\_{\rm lqg} = \prod\_{0}^{T} \left( \mathbf{x}^{T} \mathbf{Q}\_{f} \mathbf{x} + \mathbf{u}^{T} \mathbf{R}\_{f} \mathbf{u} \right) dt \tag{10}$$

where *Qf* is the covariance matrix of the measurements noise, and *Rf* is the covariance matrix of the dynamics noise (or model uncertainties). The LQG control law determination will be given by solving two problems: setting the controller for the linear quadratic deterministic problem, and setting a Kalman-Bucy filter for optimum estimation *x* ^ of state *<sup>x</sup>*. The formulation of the Kalman-Bucy filter is given by:

$$\stackrel{\wedge}{\alpha} = A\stackrel{\wedge}{\alpha} + Bu + L \stackrel{\wedge}{\alpha} (y - C\stackrel{\wedge}{x}) \tag{11}$$

where *L* is the Kalman filter gain that is obtained by solving the algebraic Riccati matrix equation:

$$\boldsymbol{A}^{T}\boldsymbol{S} + \boldsymbol{S}\boldsymbol{A} - \boldsymbol{S}\boldsymbol{C}\boldsymbol{R}\_{f}^{-1}\boldsymbol{C}^{T}\boldsymbol{S} + \boldsymbol{Q}\_{f} = \boldsymbol{0} \tag{12}$$

The gain of the optimal filter is given by *L* =*Rf* −1 *C <sup>T</sup> S*. After obtaining *L*, it is possible to obtain the transfer function of open loop LQG controller according to [10]:

$$K\_{lqq}G(\mathbf{s}) = K\left(\mathbf{s}I - A + BK + LC\right)^{-1}LG(\mathbf{s})\tag{13}$$

where *G*(*s*)=*C*(*s*)(*sI* − *A*) −1 *B* is the transfer function of the attitude dynamics.

The linear quadratic controller action can be implemented using actuators such as reaction wheels and magnetic actuators for a continuous control command. However, during orbital operations, such as rapid detumbling maneuvers, the required torqueses are usually too high for reaction wheels. Therefore, on-off propulsion strategies are used for such operations [9]. The choice of cold gas jets actuation in the present study aims to test the control in a most critical situation in terms of the small attitude adjustments difficulties.

The firing of the gas jets is controlled by a pulse-width pulse-frequency modulator (PWPF) [13]. The PWPF is an interesting option for the thrusters control system due to its advantages over other types of pulse modulators. PWPF is designed to provide propulsion output proportional to the input command. The modulator optimizes the use of propellants; it provides a smoother control and increases the equipment life. The PWPF structure is shown in Fig. 1 [10].

sampling rate is 10 to 30 times the bandwidth of the system *ωB* in closed loop [15]. A suggestion of [16] is to adopt a sampling frequency greater than 20 *ωB*, in order to have a fairly smooth control response. In [17], the indication is to adopt the sampling period *T* =1 / (10 *f <sup>B</sup>*), where *f <sup>B</sup>* =*ω<sup>B</sup>* / 2*π* and *f <sup>B</sup>* is the bandwidth of the closed loop continuous system. System dynamics frequencies' analysis and recursive tuning adjustments led to the value of 10ms as a quite satisfactory value, considering the desired pointing accuracy. Furthermore, since the maximum speed of the adopted DSP core is 600MHz, a sampling period of 10ms allows the processing of the received signal and computing the control

Processor-in-the-Loop Simulations Applied to the Design and Evaluation of a Satellite Attitude Control

http://dx.doi.org/10.5772/57219

165

The design of control systems in the continuous domain is mathematically simpler and allows the use of a large set of tools. In the case of control system design in discrete domain, the mathematical problem is quite more complicated. In addition, in the continuous time domain, the visualization of the relationship between physical reality and mathematical representation of a control system is more evident. Therefore, the usual starting point for a discrete time control design is the continuous time control system study, followed by discretization

There are several methods of discretization of a given continuous time system, and they are basically divided into open loop methods and closed loop methods. In the open-loop methods the discretization essentially consists of turning the transfer function *G(s)* in *G(z)*. This transformation is performed by substituting the terms in *s* by *z* terms in order to satisfy some criterion. In the closed loop methods, the discretization of a controller function *G(s)* is obtained taking into account information from the operation of the closed-loop system, and also the knowledge of all the transfer functions involved in the system, including the plant that is is

Several methods of discretization have been proposed and evaluated [18]. We adopted the open loop method of transformation of Tustin, also called bilinear transformation, giving satisfactory results even when compared to closed loop methods (which often have better results by taking into account the whole system), when applied to systems of low order. This method is based on the approximation of the integral represented by the factor *1/s*. The integral

> 2 1 1 - » <sup>+</sup>

The dynamics of the satellite plant in the case of discrete control remains the same as in the continuous time case. Figure 2 illustrates the LQG controller scheme after its discretization. The main change observed in this scheme is the addition of a zero-order holder in the output signal which is sent to the satellite's block, as indicated by the output named "torque". The sensors signals are received by the "measured states" gateway. Finally, the controller subsys‐

*<sup>z</sup> <sup>s</sup> Tz* (14)

can be approximated by trapezoidal integration method in order to obtain:

tem was changed by the discretization in terms of its integration function.

signal in the co-simulation in the interval between samples.

procedures.

intrinsically continuous.

**Figure 1.** Scheme of the PWPF modulator.

When the positive input in the Schmitt trigger is greater than *Uon*, the trigger output is *Um*. If the input falls below *Uoff* , the trigger output becomes null. This response is also reflected for negative inputs. The error signal *e(t)* is the difference between the output of Schmitt trigger *Uon* and system input *r(t)*. This error is sent to a pre-filter whose input *f(t)* feeds the Schmitt trigger [10].

### *3.1.2. The digital LQG controller*

The LQG controller, originally designed for continuous time systems, must be adapted to discrete time application. The appropriate selection of the sampling period *T* is a crucial factor in digital controller design, since if this period is too large there are problems in the signal reconstruction, and if it is too small, system instability and processing capacity problems can occur. In principle, one can believed that smaller sampling period is the best digital approxi‐ mation. However, if the sampling period is too small, the controller poles approach the unit, causing instability. According to the equation mapping from *s* and *z* spaces, where, we can see that if *T* is very small, tending to zero, the poles of the controller in *z* tends to 1, which are the poles of a marginally unstable system, making the closed-loop system unstable. However, it is not necessary that *T* approaches to zero to start the problems, if it is small enough the poles at *z* cannot be anymore distinguished by the computer unit. Furthermore, small sampling periods can introduce significant distortions in the system dynamics behavior [14].

Very large sampling periods may result in violation of the rule established by the sam‐ pling theorem, which says that the sampling frequency *ωs* must be twice greater than the highest signal component frequency *ω<sup>M</sup>* [15]. If the condition of the theorem is not satisfied, there are information losses in the signal reconstruction. A reasonable choice for the sampling rate is 10 to 30 times the bandwidth of the system *ωB* in closed loop [15]. A suggestion of [16] is to adopt a sampling frequency greater than 20 *ωB*, in order to have a fairly smooth control response. In [17], the indication is to adopt the sampling period *T* =1 / (10 *f <sup>B</sup>*), where *f <sup>B</sup>* =*ω<sup>B</sup>* / 2*π* and *f <sup>B</sup>* is the bandwidth of the closed loop continuous system. System dynamics frequencies' analysis and recursive tuning adjustments led to the value of 10ms as a quite satisfactory value, considering the desired pointing accuracy. Furthermore, since the maximum speed of the adopted DSP core is 600MHz, a sampling period of 10ms allows the processing of the received signal and computing the control signal in the co-simulation in the interval between samples.

provides a smoother control and increases the equipment life. The PWPF structure is shown

When the positive input in the Schmitt trigger is greater than *Uon*, the trigger output is *Um*. If the input falls below *Uoff* , the trigger output becomes null. This response is also reflected for negative inputs. The error signal *e(t)* is the difference between the output of Schmitt trigger *Uon* and system input *r(t)*. This error is sent to a pre-filter whose input *f(t)* feeds the Schmitt

The LQG controller, originally designed for continuous time systems, must be adapted to discrete time application. The appropriate selection of the sampling period *T* is a crucial factor in digital controller design, since if this period is too large there are problems in the signal reconstruction, and if it is too small, system instability and processing capacity problems can occur. In principle, one can believed that smaller sampling period is the best digital approxi‐ mation. However, if the sampling period is too small, the controller poles approach the unit, causing instability. According to the equation mapping from *s* and *z* spaces, where, we can see that if *T* is very small, tending to zero, the poles of the controller in *z* tends to 1, which are the poles of a marginally unstable system, making the closed-loop system unstable. However, it is not necessary that *T* approaches to zero to start the problems, if it is small enough the poles at *z* cannot be anymore distinguished by the computer unit. Furthermore, small sampling

periods can introduce significant distortions in the system dynamics behavior [14].

Very large sampling periods may result in violation of the rule established by the sam‐ pling theorem, which says that the sampling frequency *ωs* must be twice greater than the highest signal component frequency *ω<sup>M</sup>* [15]. If the condition of the theorem is not satisfied, there are information losses in the signal reconstruction. A reasonable choice for the

in Fig. 1 [10].

164 Computational and Numerical Simulations

**Figure 1.** Scheme of the PWPF modulator.

*3.1.2. The digital LQG controller*

trigger [10].

The design of control systems in the continuous domain is mathematically simpler and allows the use of a large set of tools. In the case of control system design in discrete domain, the mathematical problem is quite more complicated. In addition, in the continuous time domain, the visualization of the relationship between physical reality and mathematical representation of a control system is more evident. Therefore, the usual starting point for a discrete time control design is the continuous time control system study, followed by discretization procedures.

There are several methods of discretization of a given continuous time system, and they are basically divided into open loop methods and closed loop methods. In the open-loop methods the discretization essentially consists of turning the transfer function *G(s)* in *G(z)*. This transformation is performed by substituting the terms in *s* by *z* terms in order to satisfy some criterion. In the closed loop methods, the discretization of a controller function *G(s)* is obtained taking into account information from the operation of the closed-loop system, and also the knowledge of all the transfer functions involved in the system, including the plant that is is intrinsically continuous.

Several methods of discretization have been proposed and evaluated [18]. We adopted the open loop method of transformation of Tustin, also called bilinear transformation, giving satisfactory results even when compared to closed loop methods (which often have better results by taking into account the whole system), when applied to systems of low order. This method is based on the approximation of the integral represented by the factor *1/s*. The integral can be approximated by trapezoidal integration method in order to obtain:

$$s \approx \frac{2z - 1}{Tz + 1} \tag{14}$$

The dynamics of the satellite plant in the case of discrete control remains the same as in the continuous time case. Figure 2 illustrates the LQG controller scheme after its discretization. The main change observed in this scheme is the addition of a zero-order holder in the output signal which is sent to the satellite's block, as indicated by the output named "torque". The sensors signals are received by the "measured states" gateway. Finally, the controller subsys‐ tem was changed by the discretization in terms of its integration function.

**Figure 2.** Scheme of the discret LQG controller.
