Pitch Control of a Rocket with a Novel LQG/LTR Control Algorithm

Corresponding Author: Bhar Kisabo Aliyu Centre for Space Transport and Propulsion (CSTP) Epe Lagos-State Nigeria E-mail: aliyukisabo@yahoo.com Abstract: This paper presents the design, simulation and analysis of a novel LQG and LQG/LTR control algorithm for the pitch angle of a sounding rocket. These improved LQG and LQG/LTR control algorithms stem from the fact that a Riccati Differential Equation (RDE) rather than the popular Algebraic Riccati Equation (ARE) is used to obtaining the Kalman gain in the observer of the traditional Linear quadratic Gaussian (LQG) control algorithm. Thus, eight (8) different controllers were design, simulated and analysed, three (3) of such controllers are novel and two out of these novel controllers were able to recover completely the robustness lost in the traditional LQG controller. All controllers synthesized were analysed using time response characteristics of closed-loop system and compared with the LQR and LQG control system. Using the LQR controller as the benchmark for best performance and the LQG as the worst. This study shows an application option that demonstrates optimal control system design in MATLAB/Simulink® and the approach put forward here proves to be very effective.


Introduction
Classical control system design is generally a trialand error process in which various methods of analysis are used iteratively to determine the design parameters of an "acceptable" system. Acceptable performance is genrally defined in terms of time and frequency domain criteria such as rise time, settling time, peak overshoot, gain and phase margin, and band width. To meet the demands of modern technology, different performance criteria must be satisfied, in a complex multiple-input multiple-out systems requirement. For example, the design of a spacecraft attitude control system that minimizes fuel expenditure is not amenable to solution by classical methods. A new and direct approach to the synthesis of these complex systems, called optimal control theory, has been made feasible by the development of digital computer.
The objective of the optimal control theory is to determine the control signals that will cause a process to satisfy the physical constraints and at the same time minimize (or maximize) some performance criterion. In certain cases, the problem statement may clearly indicate what to select for a performance measure (Kirk, 1998).
A lot of work has been done in the area of optimal controller design for aerospace vehicles (Jianqiao et al., 2011;Das and Halder, 2014;Moshen Ahmed et al., 2011;Liu, 2017;Nair and Harikumar, 2015) but non derives it Kalman gain from a diffferential reccati equation as we will demonstrate in this study. Also, in previous works, the synthesied LQG/LTR (Barzanooni, 2015;Barbosa et al., 2016;Ishihara and Zheng, 2017) did not restore completely the robustenses of the LQR or that of the Kalman filter.
The first stage of any control system theory is to obtain or formulate the system dynamics. There refers to modeling in terms of dynamical equations such as differential or difference equations. With such equations, the system is called the plant. This aspect has been addressed fully in Chapter One. In this Chapter, the realized plant will be used to design and analyse Optimal control algorithms (Brian and Moore, 1989) This chapter is divided into five sections. Section two presents the adotped mathematical model. In section three, stability analysis was done for the rocket mathematical model. Optimal Control theory was introduced in section four hence, the design of LQR, LQG and their two novel variants. Also, LQG/LTR was indroduced here with a design for a novel variant of it. Results were discussed in section five before the conclusion section.

Rocket Model
The modern control theory concerned with Multiple Inputs and Multiple Outputs (MIMO) is based on state variable representation in terms of a set of first order differential (or difference) equations. Here, the system (plant) is characterized by state variables, say, in linear, time invariant form as: where, dot denotes differentiation with respect to (w.r.t.) t, x(t), u(t) and y(t) are n, r and m dimensional state, control and output vectors respectively and A is n×n state, B is n×r input, Cis m×n output and Dis m×r transfer matrices. For the rocket system in this study, it can be shown to be represented mathematically as: where, θ is the pitch angle, q is the pitch rate and w is the velocity in z-axis.

Stability Analysis
After developing a model for the rocket as given in (2), it is imparetive to investigate the system properties as it relates to certain characteristics of control system design. First, we want to view the trajectorry of the system without a controller associated with it. This is called the open-loop response.

Open-loop Response of a Plant
Typically, this is done by stimulating the system with a step, impulse or sine signal. The open-loop step response for the rocket system as described in (2) was simulated in MATLAB/Simulink. All the state of the system were view but only that of the pitch angle is presented. This was due to the fact that we are majorly interested in the control of the pitch angle of the rocket only.
The open-loop step response from Fig. 1 was obtained by the following code in MATLAB ® : sys=ss(A,B,C,D); step(sys,1) From Fig. 1, pitch angle keeps increasing with increase in time. This suggest instability, though we need to further investigates its eigenvalues.

Eigenvalues of a Plant
The eigenvalurs of a matrix are its roots. To further validate that the rocket system is unstable, we computed its eigenvalues. This was done easily in MATLAB ® with the command:

eig(A)
The above command gave S 1,2,3 = [3.0383, -3.8517, -0.1122] T , because only one real root is negative the rocket is unstable. Now that we are certain that the system as described in (2), it unstable we need to establish the fact that our mathematical model for the rocket is controllable-meaning, the formulated mathematical equation for the rocket is suitable for Optimal Control theory synthesis.

Controllability of a Plant
The concept of controllability was introduced by Kalman (1960) and plays an important role in the control of multivariable systems. Considering a system as described by (1), the system is said to be controllable if a control vector u(t) exist that will transfer the system from any initial state x(t 0 ) to some final state x(t) in a finite time interval. Sufficient condition for complete state controllability is defined by (3) such that M is of the rank n: The following MATLAB code was used to determine wether (2) is controllable: The result of the above code gave, the rank of (2) to be 3 and the number of uncontrollable states as 0. Hence,

Optimal Control Theory
In modern control theory, the optimal control problem is to find a control which causes the dynamical system to reach a target or follow a state variable (or trajectory) and at the same time extremize a performance index (Naidu, 2003). Thus, the design is usually with respect to a performance index-this is Optimization in its general form.
The control u(t) and state x(t) vectors are either unconstrained or constrained depending upon the physical situation. The unconstrainedproblem is less involved and gives rise to some elegant results. From thephysical considerations, often we have the controls and states, such as currents and voltages in an electrical circuit, speed of a motor, thrust of a rocket, constrained as: where, +, and -indicate the maximum and minimum values the variables can attain.
Hence we can state that the optimal control problemis about finding the optimal control u*(t) which causes the linear time-invariant plant to give the trajectory x*(t) that optimizes or extremizes (minimizes or maximizes) a performance indexwith some constraints on the control variables u(t) and/or the state variables x(t). This is summarized in Fig. 2. Thus, the Optimal control problembasically interested in finding the control u*(t) which when applied to the plant described, gives an optimal performance index J*. A lot of control algorithms fall in this class, starting with the Linear Quadratic Regulator.

Linear Quadratic Control (LQR)
Linear Quadratic Regulator seeks to minimize a cost function, in general terms, the optimal control problem is to find u which causes the system: To follow an optimal trajectory that minimizes the performance criterion (Burns, 2001), or cost function: The Hamilton-Jacobi equation is usually solved for the case of linear time-invariant plant with quadratic performance criterion (called the performance), which takes the form of the matrix Riccati equation. Thus, producing an optimal control law as a linear function of the state vector components which is stable provided the system is controllable.
Define a cost function of the form: where, over the time interval t 0 to t 1 : For a linear, time invariant plant, (7) becomes: And if (7) is a quadratic performance index, then: where, Q is an l×l symmetric positive-definite matrix and R an m×m symmetric positive-definite matrix. A first selection for the matrices Q and R is given by Bryson's rule as: Bryson's rule scales the variables that appear in (10) so that the maximum acceptable value for each term is one. Bryson's rule sometimes gives good results, often it is just the starting point to a trail-and-error iterative design procedure aimed at obtaining desirable properties for the closed-loop system (Franklin et al., 2006;Williams and Lawrence, 2007;Anderso and Moore, 1989). For this study, a combination of Bryson's rule with simulation was used to obtain the following values area used: The optimal control law can be shown to be given as: Where:

Fig. 2: LQR control schematic
Note that P is the stabilizing solotion of Algebriac Riccati Equation (ARE) in (16): The MATLAB in-built function [K LQR , P LQR , E LQR ] = lqr (sys, Q, R) was used to compute P LQR , of the associated algebraic Riccati equation as stated in (16), the state feedback gain matrix K LQR , as given in (15) and E = eig (A -BK LQR ), the closed-loop eigenvalues (User's Guide, 2014) for the associated rocket. These computed values are: The LQR controller as dipicted in Fig 3.3257 0.8509 0.0000 0.8509 0.2197 0.0003 0.0000 0.0003 0.0000 The state feedback controller with (18) as the feedback gain, was design and simulated in Simulink and the simulation result is given in Fig. 3.
It is clearly seen from Fig. 3 that the LQR has appreciable step-response characteristics and very good tracking of the pitch angle.LQR formulation has the drawback that all states must be available for feedback. This assumption is not realistic, hence the need for an observer to measure all or some of the state is inevitable.
It is impractical to expect in real-world applications that all state variables are measurable hence, the ability to implement a state feedback control law is in jeopardy.

Fig. 3: Set-point tracking LQR control of the rocket
This brings in the need for an estimate of the state vector derived from measurements of the input and output. Before an estimator is design the plant must be checked for observeberbility.

Observability of a Plant
An obsever in its simplest term is an algorithm that eatimate and predicts the state of a dynamic system. One of the most popular observer with wide range of application is Kalman Filter. A necessary condition for the Kalman Filter to work correctly is that the system for which the states are to be estimated, is observable. Therefore, you should check for observability before applying the Kalman Filter. (There may still be other problems that prevent the Kalman Filter from producing accurate state estimates, as a faulty or inaccurate mathematical model).
The system described by (2) is completely observable if the matrix N as given in (20) is of the rank n: Also, we use MATLAB ® to compute (20) and compred its rank with the system n×n matrix: If the determinant of Nis zero, the system is nonobservable. Non-observability has several concequences: • The transfer function from the input variable y to the output variabley has an order that is less than the number of state variables (n) • There are state variables or linear combinations of state variables that do not show any response • The steady-state value of the Kalman Filter gain cannot be computed. This gain is used to update the state estimates from measurements of the (real) system

Linear Kaman Filter
In real world system there would be noise (disturbance) on these measurements, and in order to analyse this problem realistically noise need to be added to the model. In general the output y is affected by measurement noise and the process dynamics are also affected by disturbance, specifically turbulence of atmosphere.
In light of this, a more reasonable model for the plant as given in (1) is: where, w(t) and v(t) are zero-mean Gaussian noise processes (uncorrelated from each other) with power spectrum: Q N is the process noise covarian and R N is the measurement noise covariance matrix. The following process (Q N ) and measurement (R N ) covariance matrices hold namely: A linear state observer is an n-dimensional linear state equation that accepts u and y as input and whose state represents the estimate of x. Suppose we construct the estimate x by replacing the process dynamics in (23) to be: To see if this would generate a good estimate for x, we can define the state estimation error in (24) and study its dynamics by: From (22) and (24), we can conclude that: This shows that when the matrix A is asymptotically stable the error e converges to zero for any input u, which is good news because it means that x eventually converges to x as t → ∞ . However, when A is not stable e is unbounded and x grow further and further apart from x as t → ∞ . To avoid this, one includes a correction term in (24) which yields: where, ŷ should be viewed as an estimate of y and K LQE a given n×k matrix. When x is equal (or very close) to x, then ŷ will be equal (or very close) to y and the correction term ( ) LQE K y y − plays no role. However, when x grows away from x, this term will (hopefully) correct the error. To see how this can be done, we rewrite the estimation error dynamics for (22) and (24) as: Now e converges to zero as long as A -K LQE C is asymptotically stable. It turns out that, even when A is unstable, in general we will be able to select K LQE so that A -K LQE C is asymptotically stable. The system in (27) can be re-written as: Now (29) is called a full-order observer for the process. Full-order observers have two inputs-the process' control input u and its measured output y (via sensors) -and a single output-the state estimatex . Schematically, the full-order observe is represented as shown in Fig. 4: Any choice of K LQE in (29) for which A-K LQE C is asymptotically stable will make the estimated state vector to converge to the actual state vector.
The Kalman gain, K LQE , for the optimal observer can be shown to be as given in (30), where P LQE is the stabilizing solution P LQE of Algebriac Riccati Equation (ARE) given as (31): With these values decleared in a MATLAB m-file, the comamd [K LQE , P LQE , E LQE ] = lqe(A,G N ,C,Q N ,R N ,P 0 ) was used to compute the kalman gain K LQE for our observer and was used to estimate the pitch angle of the rocket which gave the result in Fig. 5. The Kalman filter design has properties as given in (36), (37) and (38) The closer to zero, the elemets of the diagonal matrix of P LQE (covariance matrix) is, suggests that the error is minimal. Hence, the value of 0.0006 in (38) is associated with the kalman filter estimate of the pitch angle. Thus, it validates the result displayed in Fig. 6 (kalman filter estimate of the pitch angle is very close to the actual trajectory). The eigenvalues as displayed in (38) are all nagative, this suggests a stable observer.

LQG Control Design
Basically, the LQG approach addresses the problem where we consider a system dynamic model perturbed by a dynamical noise w and a state observation corrupted by measurement noise υ affecting the sensors data acquisition. It is essentially a combination of LQR and a Kalman based estimator (Aliyu, 2011). The estimator or observer addresses the short fall of the LQR's assumption that all states are available for feedback (via sensors). Here, all the required states will be estimated by the renowned kalman filter. From Fig. 6, z is the controlled system output and K LQR is the controller feedback gain in our case is the LQR gain. Consequently, the LQG control law is: where, x the estimated state by Kalman filter.
The designed LQG controller in MATLAB/simulink was simulated to track a pitch angle of 3 degrees (0.053 radians) and the result is dipicted in Fig. 7.
To check the stability of the LQG system, it is more convenient to consider the dynamics of the estimation error. In matrix form, the related equation can be shown as: For the LQG designed in this study, the eigenvalues of (40) are as given in (41). Notice that (34) is a combination of the roots of the LQR controller and those of the Kalman as earlier designed in (17) and (38):

Improved LQG Control Design (LQG i,1 )
Base on the dual principle, (16) is the algebriac riccatti equation for the LQR controller and (31) is its counterpart for the kalman filter. Hence, we can write (42) as the non-steady form of (31). This is a riccati differential equation (Aliyu et al., 2012): We propose an improved LQG controll design base on the solution of (42). This Riccati Differential Equation (RDE) is an initial value differential problem with initial conditon P 0 = [6.25 × 10 −5 6.25 × 10 −5 6.25 × 10 −5 ] T . Thus, we put forward an observer gain as given in (43). This kalman gain values change all through the regime of flight: State estimation was done using (43) as our Kalman gain and the result for pitch angle estimation is shown in Fig. 8a (zoom of Fig. 8b). Notice that the actual and estimated state is exactly the same, this is why both plots are super-imposed on each other as seen in Fig. 8b.
In MATLAB/Simulink ® , (43) as observer gain was combined with the LQR control designed earlier for the rocket to form the first improved LQG controller (LQG i,1 ). Simualation result of such controller is dipicted in Fig. 9.
Kalman gain for simulation in Fig. 9 gave is 916 values of Kalman gain which were updated throughout the simulation time of 5 sec.

Improved LQG Control Design (LQG i,2 )
To establish a fare comparison with the traditional LQG design, we decided to take the 916th value of the Kalman gain from the solution of (43) as given in (44) was used to designed an estimator which was able to estimate the pitch angle of the rocket as shown in Fig. 10. The estimator has properties as given in (45) and (46) (46), the estimator is stable (nagative eigenvalues). Comparing the P 11 element in (45) with that of (37), one could clearly say that the latter will give a better estimate of the pitch angle compared to the former. This fact is validated with Fig. 10, it can be seen that the algorithm poorly estimated the actual state for about 0.8 sec. When this variant of the Kalman gain was used as an observerto realise our second novel LQG algorithm (LQG i,2 ), simulation result (Fig. 11) gave us a controller with better time response characteristics than the traditional LQG. Our quest for a better controller will lead us to designing LQG/LTR hoping we can get one with better properties than that of LQG i,2 . Hence, we proceed with the traditional LQG/LTR design.

LQG/LTR Control Design
Based on the LQR solution of the above nominal optimal control, we can design a Linear Quadratic Gaussian (LQG) with Loop Transfer Recovery (LTR) controller. The LQG/LTR method is rooted in optimal control theory and in spite of systematic design procedure, shows some useful properties of robustness and good performance. In this method, the desired time response characteristics of the closed-loop system (plant and controller) must be designed in an LQG problem and then these time response characteristics are recovered at the input or output of the real plant by successive tuning of a gain in an LQR problem or reduction of the measurement noise as the case may be. Hence, there are two main methods for LQG/LTR design (Kulcsar, 2000), these are: • LQG/LTR at input of the plant by tuning the kalman gain matrix KLQE • LQG/LTR at output of the plant by tuning gain matrix KLQR

LQG/LTR Design at Plant Output
Now, breaking the LQG loop at the plant output as shown in Fig. 12, the return ratio at the output of the plant is: The associated sensitivity function of plant is defined as: Investigation the effect of decreasing the intensity Q of the LQR problem. Suppose: where, M is a fixed symmetric positive-definite weighting matrix and ߩ a positive number. Investigating the asymptotic behaviour of the closed-loop system as ߩ↓0. Also we assume that the plant in this case is minimum-phase. The feedback gain for the LQ problem is given by: Examining the situation while ߩ → ∞: And considering the fact that: From (59) it is inferred that with a q very high we can approach the return ratio of the Kalman filter. Also, during Loop Transfer Recovery, the poles of: Thus, the design procedure is: • First, we design a Kalman filter whose transfer function CΦ(s) K LQE is desirable. By choosing the power spectral density matrices w and υ so that the minimum singular value of CΦ(s) K LQE is large enough at low frequencies for good performance and its maximum singular value is small enough at high frequencies for robust stability; • When the singular values of CΦ(s) K LQE are thought to be satisfactory, loop transfer recovery is achieved by designing K LQR in an LQR problem with B=C in (59), Q = Q 0 + ρM, where ρ is a scalar.

LQG/LTR Design at Plant Input
Loop Tranfer Recovery at input of the plant (as shown in Fig. 13) is achieved by primarilily tuning the kalman gain matrix KLQE,its return ratio isgiven in (62). Notice that (62) is just a form of (47) being casted for LQG/LTR design at Plant input. Thus, extending the dual principle of the LQG problem to the LTR definition: Since the K LQR controller is ready, we only need to design the Kalman filter which can help recover the robustness of LQR controller. Notice, we are simply taking the converse steps of LQG/LTR at the plant output design as earler outlined. The Loop Transfer Recovery design is done by increasing the spectral density S w of the plant noise w(t). This method is popular in the research community (Ishihara et al., 2005): We choose different values of S w to design the Kalman filter and compare the results to the LQR results. We choose different values of S w by doing: where, I is an identity matrix and we want to find the K LQE feedback while ρ → ∞. The approach suggested in (64) is similar to that expressed in (49). Here, we were able to prove with simulation that the same approach can be implemented for LTR at plant input.
The Kalman filter gain for the LQG/LTR at the plant input for this study is obtained with MATLAB ® command [K LQE ,P LQE ,E LQE ] = lqe(A,GQnG,C,Q NN ,R N ,P 0 ). Next we show the system response of LQG/LTR control set-point tracking response for different values of ρ in Fig. 14.  Next we will design LQG/LTR using (44), as the Kalman gain for the observed. Remember this Kalman gain was harvested from a DRE. Simulation result for such design are presented in Fig. 15. We call this controller LQG/LTR i .

Discussion of Results
It can be clearly seen from Table 1, that the LQR controller has appreciable time ressponse characteristics with zero steady-state-error (see Fig. 4). This is viewed as an ideal situation and its position is well supported in literature (Aschepkov, 2016;ViorelBadescu, 2017;. The combination of LQR controller with Kalman filter as an observer caused some of the fine properties of both the LQR controller and kalman filter to be lost. The result is a closed-loop system with an overshoot (0.9%) as shown in Fig. 8. This result is also in line with popular litereture position (Wang Jian, 2016;Zhu et al., 2017). Our first novel controller, LQG i,1 , used updated values of kalman gain all through the regime of simulation to arrive at the result dipicted in Fig. 9 and Table 1. Notice that this controller was capable of restoring almost all the robustness of the of the LQR. In practical terms, LQG i,1 might be challenging to implement oweing to the fact that if the rocket will fly for a long time, Kalman gain values might increase beyond the 916steps that we have here for a simulation time of 5 second. This means large memory alocation on a hardware and at the same time prolonged computation throughout the regime of flight. Hence, we introduced the second novel controller LQG i,2 as a practical substitute to LQG i,1 . Here we used the last integrated value of Kalman gain as our preferred values of estimator gain to design the novel LQG i,2 controller. The simulation result of LQG i,2 are as dipicted in Table 1 and Fig. 11. This controller is also very promising in the sense that it was also capable of restoring the robustness of the LQR controller with even a better settling time than that of the LQR. The traditional LQG/LTR at plant input was designed, and the results for three scenerio are dipicted Fig. 14. The tunned Kalman gain with the largest value of ρ = 10 ( Fig. 14a) was able to restore all the robustnes and is comparable with the LQR, LQG i,1 and LQG i,2 controllers. For LQG/LTR with ρ = 0.1, in Fig. 15c, this controller is comparable with the LQG and LQG/LTR with ρ = 1 in Fig. 14b has better time response characteristics than the LQG but less, compared to LQR, LQG i,1 and LQG i,2 . For close perusal of the three LQG/LTR controllers synthesied, Fig. 14d is perspective. The novel LQG/LTR at plant input was designed using the Kalman gain from the solution of RDE and simulation results are as dipicted in Fig. 15. Notice, that LQG/LTR i with ρ = 0.1in Fig. 15a was able to restore all the robustness needed and has better time response characteristics compares to LQG i,1 and LQG i,2 as shown in Table 1. We tried to improve on this controller by making ρ = 1 and ρ = 10 but got the same result ( Fig. 15b and 15c). Notice, for this 3 LQG/LTR i controllers to be shown to have the same time response characteristics, Fig. 15d was plotted. Hence, we consider all the controllers in Fig. 15d to be just one-LQG/LTRi.

Conclusion
To design an optimal controller for a rocket, a mathematical model was developed and simulated in MATLAB/Simulink ® . The cascade combination of a matthematical model and a control algorithm (with feedback) gives a closed-loop control system (controller). Thus, eight controllers were designed in this study: LQR, LQG, LQG i,1 , LQG i,2 , LQG/LTR (ρ = 10), LQG/LTR (ρ = 1), LQG/LTR (ρ = 0.1) and LQG/LTR i . Simulation for all controllers were carried out in MATLAB/Simulink ® . Results of all controllers were compared based on time response characteristics and the LQR controller was used as the bench mark, LQG i,1 , LQG i,2 and LQG/LTR i are novel in nature based on the origin of the Kalman gain used the observer design. LQG/LTR i outperformed the other two coontrollers in terms of time response characteristics (better than the LQR).

Future Work
It is imparative to formulation an analytic mathematical proof for the simulation concept and result of the LQG/LTR control presented in this study.