System Identification and LQR Controller Design with Incomplete State Observation for Aircraft Trajectory Tracking

This paper presents a controller design process for an aircraft tracking problem when not all states are available. In the study, a nonlinear-transport aircraft simulation model was used and identified through Maximum Likelihood Principle and Extended Kalman Filter. The obtained mathematical model was used to design a Linear–Quadratic Regulator (LQR) with optimal weighting matrices when not all states are measured. The nonlinear aircraft simulation model with LQR controller tracking abilities were analyzed for multiple experiments with various noise levels. It was shown that the designed controller is robust and allows for accurate trajectory tracking. It was found that, in ideal atmospheric conditions, the tracking errors are small, even for unmeasured variables. In wind presence, the tracking errors were proportional to the wind velocity and acceptable for small and moderate disturbances. When turbulence was present in the experiment, state variable oscillations occurred that were proportional to the turbulence intensity and acceptable for small and moderate disturbances.


Introduction
The growth of commercial aviation results in increased requirements for the guidance systems as they need to provide better capabilities of the system while not lowering the safety level. Disturbances that act on the object (e.g., wind) can cause departure from the planned trajectory and thus lead to potentially dangerous situations. Moreover, if aircraft follows the prescribed path accurately, it is possible to optimize certain flight segments in order to save time and cost e.g., by using continuous descent operations [1]. Thus, the aircraft trajectory tracking problem is of great interest [2][3][4][5][6][7].
The full-state Linear-Quadratic Regulator (LQR) controller was successfully applied for trajectory tracking for various aircraft: fixed-wing [6,8,9], helicopters [10,11] and UAVs [12,13]. In case of incomplete state information, the LQG controller can be used instead. In [14] the LQG controller was designed to handle the tracking task for a linear business jet aircraft. In [15] a controller was designed for a nonlinear system that was sequentially linearized for missile path-following. In [16] the system stability is regained by using a LQG controller. A different approach is used in [17] to handle with the partial information-reduced-order observer is used to estimate unmeasured states.
The tracking problem can be also solved by using model predictive control (MPC). In this approach an object representation is used to predict the dynamic response of the object and basing on this information optimal control strategy is found. In [18] a nonlinear MPC was used to deal with tracking of a small-scale helicopter with unmodelled dynamics for a box manoeuvre in turbulent conditions. In [19] a nonlinear MPC was used to track fixed wing UAV in extreme manoeuvres (fast turn, sharp turn, flip), in case of actuators failures and for obstacle avoidance. The MPC was also used for off-road unmanned ground vehicle path tracking with taking into account sliding and steering constrains [20].
In the paper we propose another, new approach-to use LQR in which unmeasured states were directly omitted.
More precisely, this is a LQR problem with static output feedback, where the gain matrix is derived from the optimal gain matrix of the standard LQR problem by simply deleting columns corresponding to these state vector coordinates which are not present in the output vector. This simple construction proved to be very effective. We show conditions when it can be applied, that is, when the resulting closed-loop system is asymptotically stable. This approach requires an accurate model representing the aircraft dynamics, however, that can be obtained using system identification techniques. Several methods can be used for that purpose and those related to aerodynamic modelling are presented, e.g., in [21]. The most popular approach is the output-error method, which was used in numerous applications including icing modelling [22], missile microspoiler parameters estimation [23], aerodynamic coefficients identification [6,24] and performing near real time data compatibility check [25]. However, the output-error approach does not allow for process noise presence and thus aerodynamic turbulence strongly affects the identification quality. The problem is solved by including Kalman filter. In [26] the stability and control derivatives of the linear model were found by using this approach. In [27] the linear Kalman filter was used around the trim points to find the nonlinear model. A more common approach when dealing with nonlinear aircraft models is to use the Extended Kalman Filter. In [28] this was done for kite flight path reconstruction. Aerodynamic coefficients were estimated by using this approach in [29]. The nonlinear model can be also considered when Unscented Kalman Filter is introduced as in [30] for dealing with upset flight conditions. Unscented Kalman filter was also applied to identify a fixed-wing aircraft and a rotorcraft in [31]. Those results were compared with Extended Kalman Filter outcomes and time domain output error approach. It was found that the computational load of the Unscented Kalman Filter is at least twice of the Extended Kalman Filter and it does not result in significant gain in convergence or accuracy when aerodynamic coefficients are estimated [32]. Another possibility to deal with environmental uncertainties was presented in [33] where a nonlinear observer supplied by a Kalman Bucy filter was used to adjust tire-ground parameters in near-real time.
This paper is organized as follows. A nonlinear aircraft simulation model is presented in Section 2. In the remaining part of the paper this model is a treated like a real aircraft. This approach is often used in flight dynamics, as it allows to reduce costs when new designs or system modifications are tested [34,35]. In Section 3 the algorithm used for identifying the aircraft aerodynamic coefficients is shown. Estimation results and their validation is presented in this part as well. Problem formulation (trajectory tracking) is shown in Section 4. To solve this task a LQR controller was designed thought the approach shown in Section 5. The case scenario is presented in Section 6. The final outcomes are shown in Section 7. This is done for ideal atmospheric conditions and under atmospheric disturbances (wind, turbulence). The paper ends with short conclusions summary.
The novelty is twofold. First, the system identification approach with Kalman filtering is used to obtain an accurate nonlinear aircraft model for LQR design to track the aircraft trajectory. The second aspect of the novelty is related to the state measurements unavailability. In the case of incomplete state measurements a common practice is to use a LQG controller. In this paper we show that it is possible to solve our problem also by LQR with static output feedback-that is, without the restoration of the full state vector by an observer or a filter-applying gain matrices calculated during standard LQR controller synthesis. We explain when this method can be applied (in Section 5).

Nonlinear Aircraft Model
In the study, the Basic Aircraft Model (BACM) was used to represent the real aircraft. This is a nonlinear aircraft simulation model developed in the Matlab/Simulink environment. Besides aircraft dynamics, BACM also contains environment and sensors models. The environment model is based on International Standard Atmosphere model and WGS-84 system. Throttle position, aileron, elevator and rudder deflections are the inputs for the model. Model outputs consist of aircraft position and attitude, linear velocity components and angular rates, linear acceleration components and aerodynamic angles. Those elements are obtained from a IMU model (consisting of three-axis accelerometer and three-axis rate gyro) and pressure tube model.
The reason for BACM selection was twofold. First, using simulation models in flight dynamics is a common practice as it allows us to greatly lower the costs when new modifications or updated algorithms are to be assessed. Secondly, this would allow for direct comparison with [6] where BACM model was used for trajectory tracking when all states measurements were available.
The BACM is a nonlinear simulation model described in body axes. Equations of motion are obtained from momentum Π and angular momentum K O change theorems [36]: where F and M O are external force and moment, Ω is the angular rate and the dot symbol denotes time derivative. The object is treated as rigid body, which centre of gravity coincides with the origin of the body axes frame (Oxyz). The coordinate frame is shown in Figure 1. As the system identification experiment lasts for a short period of time, fuel consumption is neglected i.e., mass is assumed constant. Moreover, the airplane has a vertical symmetry plane in terms of geometry and mass. Due to those assumptions momentum Π and angular momentum K O are given as [36]: where V O is the linear velocity, m stands for mass and I is the inertia matrix in which I xy = 0 and I yz = 0. This leads to the following equation set [36]: where u, v and w are linear velocity components, p, q and r are angular velocity components, φ and θ are roll and bank angles,q is the dynamic pressure, S is the wing area andc and b stand for mean aerodynamic chord and wingspan, respectively.
The aerodynamic coefficients C () are given as [37]: C l = C l 0 + C l β β + (C l p + C l pα α) + (C l r + C l rα α)r + C l δA δ A + C l δR δ R C n = C n 0 + C n β β + (C n p + C n pα α) + (C n r + C n rα α)r + C n δA δ A + C n δR δ R where α and β are angle of attack and sideslip angle, δ A , δ R and δ R and ailerons, elevator and rudder deflections, k is the drag polar coefficient, x LH is the horizontal tail arm WB denotes wing-body and H index stands for horizontal tail. The aerodynamic angles were obtained directly from linear velocity components: α = arctan(w/u) and β = arcsin(v/V O ). In the aerodynamic coefficients expressions, C L 0 , C Y 0 , C D 0 stand for lift, side and drag force aerodynamic biases and C l 0 , C m 0 , C n 0 are rolling, pitching and yawing aerodynamic moment biases. The nondimensional coefficients C i j express the i-th force or moment coefficient change due to j-th flight parameter or control, e.g., C l β is the rolling moment coefficient change due to sideslip angle. The cross-coupling derivatives C i jk represent the i-th force or moment coefficient change due to k-th flight parameter effect on j, e.g., C n rα is the yawing moment coefficient change due to yaw rate resulting from angle of attack change.
The model was excited with D-optimal multi-step inputs presented in [37]. In that study, it was shown that designed inputs can provide accurate results for this aircraft if the flight is performed in calm atmosphere. The inputs are shown in Figure 2.
While the model response was evaluated the states (u, β, α, p, q, r, φ, θ) were corrupted by the process noise inclusion. After the evaluations, measurement noise was added to the outputs (states and acceleration components, i.e., a x , a y and a z ). From this point, the nondimensional aerodynamic parameters were considered as unknowns that must be estimated.

Aircraft System Identification
Presented aircraft model can be described by the following equations [32]: where f and g are nonlinear functions, x, y and z are states, outputs and measurements, respectively, u is control, Θ m are model parameters, w and ν are process and measurement noise, whilst F and G are their distributions and t k is time at discrete k-th point.
The unknown parameters Θ can be found by using the maximum likelihood principle. In this approach, a set of parameters that maximize conditional probability of observing measurements z is found [32]: For the multidimensional normal distribution, the probability density function is given as: The same result is obtained by performing a simpler task-minimization of the negative log-likelihood [32,37]: where tilde symbol denotes predicted value. To estimate the residuals covariance matrix following formula can be used [32,37]: For uncorrelated measurements errors and after neglecting constant terms the cost function simplifies to [32,37]: In the study, the aircraft response was recorded in turbulent conditions. Due to process noise presence (turbulence), it was not possible to use the output error method, which is the most frequently used system identification method. To solve this issue, a Kalman filter was incorporated in the estimation algorithm accordingly to [32]. Resulting filter error method, which can be seen as the output error method extension, allows also for less local minima and faster convergence. To achieve this aim, the unknown parameters vector consists of model parameters and process noise distribution On the other hand, one has to bear in mind that the filter error method will allow us to obtain good match between the identified outputs and measurements even if the estimates are inaccurate. Thus, when performing the identification, special attention must be paid to statistical properties of the estimates. Moreover, using data from multiple experiments (for which process noise is different) requires estimating multiple process noise distributions.
Nonlinear model is incorporated into the Kalman filter prediction step to estimate the state variables x and outputs y [32]: The correction step is given as [32]: The steady-state Kalman gain matrix is used, as the nonlinearities are not large and the system is time invariant [32]: The output matrix is found through numerical linearization: State prediction error covariance matrix P is obtained by solving the Riccati equation [32]: The equation can be solved by eigenvector decomposition (Potter's method). First, the Hamiltonian matrix is build [32]: After evaluating its eigenvectors and eigenvalues, those eigenvectors are sorted in descending order with respect to real parts of their corresponding eigenvalues. The obtained eigenvectors rearranged matrix is further decomposed [32]: Finally, the state prediction error covariance matrix P is given as [32]: To perform cost function minimization Gauss-Newton algorithm was used to estimate the aerodynamic coefficients [32]:Θ where the Fisher Information matrix is given as [32]: and the gradient matrix is as follows [32]: Finite difference approximation is used to evaluate gradients of each model response j for each model parameter l: where δΘ m,l is the l-th parameter small perturbation and e is a column vector with 1 in l-th row and zeros elsewhere. Evaluation of the predicted outputs when the model parameters are perturbed y(t k , Θ m + δΘ m,l e l ) requires evaluation of the predicted statesx(t k , Θ m + δΘ m,l e l ) and corrected stateŝ x(t k , Θ m + δΘ m,l e l ) for the perturbed parameters as well. This is done by evaluating (11) and (13) for the perturbed parameters Θ m + δΘ m,l e l .
To provide that the measurement noise covariance matrix GG T is physically meaningful it must be positive-semi definitive. This is true when all eigenvalues of KC are real and less than one. The KC matrix tends to be diagonal-dominant, thus, the condition can be approximated by verifying if the diagonal elements satisfy the condition. The nonlinear inequality constraints can be approximated by [32]: If this is satisfied for all conditions, the iteration is complete. Otherwise, the solution is given by [32]: The gradient in constrained equation M is [32]: where a = 1, ..., n c and b = 1, ..., n p are element indices, n c is the violated constrains number and n p is the estimated parameters number. The auxiliary vector is [32]: In the method, the R is estimated when all other parameters are fixed, and their values are estimated for fixed R. This approach works well if there is no strong correlation between those two.
When the strong correlation exists the F matrix is updated to compensate for revised R. This is done by using a heuristic procedure [32]: where ξ, η, ζ are indices, n w is the number of process noise variables being estimated, while old and new denote former and updated terms. The described filter error method was used to identify nondimensional aerodynamic derivatives. We used the BACM model without assuming any knowledge of the aerodynamic derivatives, as in a case of dealing with true aircraft. To determine the aerodynamic coefficients, we identified the system by using the described Filter Error Method. When performing the system identification, we assumed initial values of the aerodynamic derivatives based on general flight dynamics knowledge (e.g., phugoid motion due to potential and kinetic energy conversion), empirical formulas [38] and aerodynamic derivatives values for similar aircraft. As estimating to many parameters at the same time may lead to local minima (and thus less accurate results), extends the estimation time and lowers convergence, a two-stage procedure was adopted. This approach was possible because cross-coupling had a secondary effect on aircraft dynamics.
Step 1a: first, while performing the identification it was assumed that u, α, q and θ have primary influence on longitudinal motion. Those components formed the state vector used to evaluate the outputs. Besides mentioned state variables, the output vector consisted of longitudinal and vertical acceleration (a x and a z ) for better convergence. The remaining motion variables (β, p, r, φ, a y ) were weighted by 0, thus lateral-dimensional dynamics did not influence the longitudinal motion outcomes. In this step only the coefficients related to longitudinal motion were identified (C L , C L WB , C L H αH , C L H δH , C D 0 , C m 0 , C m q ). The remaining aerodynamic derivatives were kept fixed at their initial values.
Step 1b: then, it was assumed that β, p, r and φ influence mainly the lateral-directional motion. Besides mentioned state variables, the output vector consisted of lateral acceleration a y for better convergence. The remaining motion variables (u, α, q, θ, a x and a z ) were weighted by 0, thus longitudinal dynamics did not influence the lateral-directional motion outcomes. In this step only the coefficients related to lateral directional motion were estimated ( C l β , C l p , C l r , C l δA , C l δR , C n 0 , C n β , C n p , C n r , C n δA , C n δR ). The remaining aerodynamic derivatives were kept fixed at their previous values (i.e., estimated from longitudinal motion data in Step 1a).
Step 2: in the second stage, all state variables and acceleration components formed the output vector. The cross-coupling derivatives (C Y pα , C Y rα , C l pα , C l rα , C n pα , C n rα ) were also included in the estimation, i.e., all aerodynamic derivatives were identified. The values of the aerodynamic derivatives in the first iteration were taken from the estimates obtained when only longitudinal and only lateral-directional data was used (Step 1a and Step 1b, respectively).
In each step (Steps 1a, 1b, 2) the estimation was stopped when the relative change of the cost function in (10) was smaller than 1 × 10 −4 .
Relative deviations of the estimated aerodynamic coefficients are shown in Table 1. A rule of thumb is that relative standard deviations below 10% denote accurate estimates [32]. From the obtained outcomes, it can be seen that all the parameters were estimated with high precision. The results are also comparable with the outcomes obtained in previous studies [6]-when process noise was not present in the data and thus output error method was used.
The estimated model was also validated by feeding it with the data not used during system identification and comparing obtained time histories with the measured ones. This is shown in Figure 3, where blue lines denote measured responses and red lines stand for the estimated model outputs. From Figure 3 it can be seen that a very good visual match was obtained for the estimated model. To express this match quantitatively J RMS was used [39]: where W is a diagonal weighting factor matrix. When the observations are in given in ft/s, ft/s 2 , rad or rad/s the J RMS ≤ 2 denotes high accuracy of the estimates if a six-DoF model with cross-coupling effects is used [39]. In the study, longitudinal velocity u and linear acceleration components a x , a y , a z were given in m/s or m/s 2 , respectively. Therefore, a factor of 0.3048 was selected for those outputs in the weighting matrix W. The remaining non-zero elements in the weighting matrix were set to 1. In result the J RMS = 0.9597, what indicates high accuracy of the estimates and supports previous outcomes. To assess the estimated model predicting capabilities Theil Inequality Coefficient can be used [39]: A rule of thumb is that TIC ≤ 0.3 indicates good predicting capabilities [39]. The TIC obtained for the validation data was TIC = 0.1387. Thus, it also confirms that the model was estimated with high precision.

Output Tracking Problem in Aircraft Modelling
Aircraft model development and verification leads to the output tracking problem, inherently. Availability of the aircraft state vector is crucial for it. However, the aircraft full state vector usually is not measured: side velocity, angular rates or WGS Cartesian coordinates are not used in standard flight control. Environmental state is known only partially: GPS can be used to obtain wind velocity and direction, but vertical air velocity is usually unknown. Usually, Pitot tube is used to measure angle of attack and sideslip angle, but it is also possible to obtain aerodynamic angles from measured velocity components or by using five-hole differential pressure probe. Another approach is to use Flush Air Data System, but this system needs to be designed, especially for each particular aircraft type. In practice, turbulence parameters are not available and can be identified for, e.g., Dryden or von Karman spectra. Air density ρ, crucial in flight mechanics applications is hard to measure and can be obtained only indirectly from ideal gas law.
Thus, the tracking problem in aviation is then essentially an output tracking problems in which number of measurements m is usually less than state dimension n, m < n.
The output tracking problem can be formulated in general as follows: Given the measurements and inputs of a real aircraft and registered atmospheric conditions, find the inputs that provide the computed aircraft model outputs to be close to the measurements. Optionally, controls providing the outputs and the measurements matching should be close to those recorded in flight.
The success in tracking problem handling essentially depends on data availability and accuracy. If all states of an aircraft are available, the computed control u(t) should assure the state x(t) local error minimization with regard to the recorded state x r (t). In output tracking problem one should require control u that uniformly minimizes local error of the output y(t) with regard to the measurements y r (t).
In case of incomplete state measurements, when p < n, much less accuracy of tracking error is to be expected, in general. The control that realizes output tracking is expected to be less regular than that obtained in tracking based on full state measurement. Issues that concern controlability and stabilizability will hardly influence control evaluation for output tracking problem.
The tracking problem is of an evolution type that makes use of recorded data given in long sampled time series. Thus, the time-domain approach are in favour over the frequency domain methods for its analysis.
There are several possibilities of handling the output tracking problems in time domain. The simplest one, based on classical PID controllers, is rather precluded because of large dimension of problem (greater than twelve) and strong cross dependencies between the system states.
The robust H-type (H 1 , H 2 , H ∞ ) method defined in the frequency domain (e.g., [40,41]) can be formulated in time domain, but this require selection of pre-and post-compensator modules for plant (aircraft) dynamics, which is very difficult for real systems.
The most common approach to the output tracking problem is to use the Linear Quadratic Gauss regulator (LQG). It makes use of Kalman filter to restore the state vector from measurements in some stochastic sense. It is widely known that LQG was very successful in a variety of aerospace applications [40]. However, the LQG control requires knowledge (or assumption) on probabilistic characteristics of model and measurements disturbances, which are assumed to be Gaussian. Such information is rarely known in practice, that makes the use of LQG problematic in the considered output tracking problem.
In the present work, the LQR with static output feedback is proposed for handling the output tracking problem. This version of LQR is not widely used, but it will be shown here that the LQR can be applied successfully if some conditions (concerning detectability and stabilizability of the system, but not only) are fulfilled [42,43]. It will be shown that the robustness of LQR seems to be enough for assuring good properties of the controller in output tracking problem considered here.

LQR with Static Output Feedback
The output tracking problem can be solved by a LQR optimal control. The LQR working principle with incomplete states observation is presented in Figure 4. The LQR control problem consists in the determination of the output feedback control law µ o (index o indicates the output feedback) [40,42] where u(t) ∈ R m and y(t) ∈ R p are output and control vectors at time t, and m, p are dimensions of control and output, respectively. The controller design process applied in the study is shown in Figure 5. The optimal control law of type (32) is found by solving the constrained optimization problem where J(µ o , x 0 ) ∈ R is a quadratic performance index which will be minimized with the initial state x 0 ∈ R n , where n is the dimension of the state space and weighting matrices Q ∈ R n×n and R ∈ R m×m for state and control, respectively. The constrains are given by the control law (32), together with the linear plant and output modelṡ with initial conditions where A ∈ R n×n , B ∈ R n×m and C ∈ R p×n are the state, control and output matrices, respectively. It is assumed that the initial state x 0 is a random variable of a known probability distribution.
If the full state is measured, p = n, C = I, y ≡ x, the initial state x 0 is given as a (deterministic) vector of numbers and we have the classical LQR problem: Index s of the control law indicates the state feedback. When [44]: in which the optimal gain matrix K * s ∈ R m×n is calculated as where P * s ∈ R n×n is the unique, positive semidefinite, symmetric solution of the stationary Riccati equation Moreover, the closed-loop dynamic system (38) and (40), now with the dynamicṡ is asymptotically stable and the criterion (37) takes the value: LQR control has favorable robustness properties, e.g., in SISO case: infinite gain and 60 • phase margin and thus is appealing in various applications.
Unfortunately, it is not the case when output is given by (35) and p < n. However, as we will show, the gain matrix K * s , that is the solution of the standard LQR problem, can still be used.
In [45] it is assumed, that control u(t) is generated also via output linear feedback with time-invariant feedback gains, i.e., in (32): and the optimization problem (33)-(36) can be rewritten as: In the case, when the state transformation matrix is stable, that is its eigenvalues lie in the negative halfplane, and, as in [45], the initial state x 0 is a random variable uniformly distributed on the surface of the n-dimensional unit sphere, in order for K * o to be optimal it is necessary that K * o , positive semidefinite P * o and positive definite L * are solutions of the system of nonlinear equations [42,45]: is given by (51). The optimal value of the criterion (48) is: Unfortunately, unlike the standard LQR problem with state feedback (37)-(40), Equations (52)-(54) are only necessary conditions, and there might be several solutions of these equations which are not optimal [45]. Moreover, conditions for the existence and global uniqueness of solutions to (52)-(54) such, that P * o is positive semidefinite, L * is positive definite and the matrix is stable are not known [46]. The latter problem itself is often cited as one of the difficult open problems in systems and control for which a satisfactory answer has yet to be found [43,47,48]. The algorithms solving LQR problem with output feedback were proposed by many authors [42]. The most popular solve many times the Riccati Equations (52) and (53) for properly changed, fixed matrices K o [42,45,49]. A basic issue for all of them is the selection of an initial stabilizing static output feedback gain K o [42]. Methods proposed in [42] are quite complicated. Since it is a very difficult problem, perhaps NP-hard [47], heuristic methods can be proposed, which proved to be effective while solving specific problems.
One of the specific class of problems is when the output vector y is formed by some of the coordinates of the state vector x. Let us denote the set of indices of x present in y as Z (it may be ordered, but it is not necessary). Hence, we may write: The matrix C will be then formed from the n-dimensional row versors e T i corresponding to subsequent elements of the set Z, that is Let us take as a stabilizing gain matrixK o the submatrix of the matrix K * s , being the optimal solution of the standard (that is with state feedback) LQR problem (37)- (40), which contains only columns corresponding to the coordinates of the vector y, that is The closed-loop system (49) will take the form: Let us notice that the n × n matrix is built of versors e i and n-dimensional vectors of zeros. A versor e i is in the i-th column of W when i ∈ Z, that is the coordinate x i of the state vector is one of the coordinates of the output vector y. The remaining columns contain zeros. We may write: where χ Z (i) is the characteristic function of the set Z: Let us denote by b T i , i = 1, . . . , n the i-th row of the matrix B, by k * sj the j-th column of the matrix K * s and by H the product of matrices B and K * s . We will have: Let us now apply the Gershgorin circle theorem [50] to the solution of the standard LQR problem (37)- (40). The LQR theory guarantees, that the system with matrix A * s = A − H is asymptotically stable. Assume that the sufficient condition resulting from the Gershgorin theorem is satisfied, that is for all rows i = 1, . . . , n of A * s there is: and Our system matrix in the case of the output feedback is The matrix HW is simply the matrix H with zeroed columns corresponding to those coordinates of the state vector x, which were not taken to the output vector y, that is having indices j / ∈ Z. Hence, the Gershgorin conditions for all the eigenvalues of the matrixÃ o to have negative real parts will take the form: Comparing (69) with (66) we obtain that when the matrices A and H satisfy the conditions: that is, when the matrices A, B and K * s , making use of (64), satisfy the conditions: and from (66) and (68) ∀i / ∈ Z, a ii < 0 (73) the closed-loop system with linear output feedback (60) for K o =K o = K * s C T is asymptotically stable. It is not a solution of the problem (48)-(50) itself, but of a different LQR problem with output feedback, where some matrices were modified [43].
In LQR problems there are no constraints either on the current control u(t), nor on the current state x(t). Such problems, including those with mixed functional constraints on current control and state are very important in many areas of engineering. The methods of their solutions are presented e.g., in [51,52].

Case Scenario
The task is to compute the control for a cargo aircraft resulting in its flight path close to a given trajectory with time synchronization. In other words, aircraft should track the trajectory at the same time time points as during the flight recording. Such a problem is thus the output tracking problem in the sense defined in Section 4. In practice, this is a problem of flight management that provides high safety level in dense airspace.
The key point in practical implementation of the output tracing problem is the number and quality of the available measurements.
The full state vector of an aircraft dimension of n = 12 is where u, v, w are the linear velocities defined in aircraft fixed coordinate system, p, q, r are the angular velocities defined in aircraft fixed coordinate system, X, Y, Z are the aircraft's positions in WGS coordinate system and φ, θ, ψ are orientation angles of an aircraft with respect to the WGS coordinate system. It is assumed that the basic instrumentation is available, including GPS. The five hole differential pressure probe is also available, thus the sideslip angle can be measured. The measurements are thus: total velocity V, sideslip angle β, latitude λ and longitude κ, attitude angles: pitch θ, roll φ and yaw (heading) ψ. The positions X and Y are computed from latitude and longitude taken from GPS. The side velocity v is computed from total velocity V and sideslip angle β. The output vector is then: (75) Its dimension p = 7 is considerably less than that of state dimension n = 12. The real problem described above are surely too hard to handle with PID controllers, because of its nonlinearity and couplings between states. Tuning the separated PID controllers seems to be rather problematic in such case. The extended LQR controller, in turn, gives a chance to find control K for the whole system at once states due to the use of weighting matrices Q and R.
The synthesis of the LQR controller requires selection of weighting matrices Q and R and deriving linearized state and control matrices A and B.
Matrices A and B are computed as the Jacobians of nonlinear model equations f (x, u): They depend on actual state x(t) and control u(t) and were computed at each time step using an approximate finite difference formula of fourth order. The four point formula was chosen since it turned out to be more accurate than two point formula and was stable enough. Jacobians were computed for non-equilibrium states in the same manner as it was done in many papers concerning generalization of the LQR control for nonlinear problems. Both Jacobians did not change significantly in time, so we tried also the use of constant Jacobians computed once for all at the beginning of simulated flight. No significant differences were observed in that case.
Weighting matrices Q and R selection (both diagonal) matrices Q and R is a key aspect of the LQR design. Poorly selected Q and R may cause inaccurate tracking or aircraft motion instability. To assure stability, the weighting matrices components maximum values were limited. This is because too-high values caused oscillations of control that lead to flight instability. On the other hand, too small contributions of the weights result in inaccurate tracking and departure from the prescribed trajectory. Designing of matrices Q and R is not trivial, because of implicit dependencies between their elements, which are not known. Therefore, experience based trial and error procedure was used to select them. First, low values for elements of Q and R are assumed and simulation runs are performed. In those simulations, tracking accuracy and stability of flight are examined. Next, the elements of Q and R are increased until instability appears. If it is observed, some elements of Q or R are decreased by a few percent until stability is recovered. Usually a few tens of runs are required for tuning matrices Q and R.
Note, that matrix Q has two zeros on the positions that correspond to the vertical velocity w and roll angle φ. Such a choice was found to fulfil detectability conditions.
The system is controllable since matrix A − BK is stable in each time step. This was checked by computing its eigenvalues and verifying that their real parts are negative. What is interesting, there was one unstable mode of A (the spiral one), but it was both observable and controllable. The system is observable for chosen outputs since they have been selected in such a manner to assure that observation matrix O = C CA CA 2 . . . CA (n−1) T has full rank. It has been verified in a standard manner by checking nonsingularity of Gramian matrix O T O. Due to the inherent nonlinearity of the considered output tracking problem the the LQR controller properties can be assessed only qualitatively through numerical simulations. Robustness of the LQR control was investigated experimentally by introducing disturbances induced by atmospheric phenomena: steady wind and turbulence.

Results
A number of tests were performed to assess quantitatively the stability and detectability of the nonlinear output tracking problem based on the LQR control. The simulation approach has been used for this purpose. The robustness of output LQR control with respect to atmospheric disturbances was examined by applying two types of such of disturbances: systematic (constant wind) and irregular disturbances with zero mean value (turbulence).
Assessing the stability and tracing accuracy of the extended LQR control was done for long lasting flights. Three cases were examined: flight with no atmospheric disturbances, flight with steady side wind and no turbulence and flight with omnidirectional turbulence.
In Figure 6 a reference trajectory (only the geometry) is shown. Figure 6. Exemplary trajectory.

Ideal Atmospheric Conditions
To investigate the extended LQR control stability and tracking accuracy during long lasting flights, a scenario with no environmental perturbations was evaluated first. For this purpose, the local errors of states e(t) were examined. In Figure 7 output tracking errors are shown. The flight is stable and it can be seen that the output tracking is very accurate in the considered case. The local errors e(t) are low (roughly a few percent of the nominal values) and roughly uniform for over all the flight time, even these of the unobserved states w, q, r, Z. It can be stated that linear velocity errors are quite small, being of order 0.1-0.5 m/s. The position errors are less than ten meters at the end of flight, what is possible only if velocity errors are small enough. Considerable oscillations, which are of order 5 deg/s, are observed in roll angular rate p which is not measured. However, other states that are not measured (w, q, r, Z) do not exhibit such behaviour. Oscillations of measured roll angle φ, attaining almost two degrees, are a consequence of roll rate oscillations. Pitch and yaw angle errors are very small, generally up to one degree over all the flight time. No significant oscillations of them are observed.
The LQR generated controls (control surfaces deflections) match good those that were recorded in flight, as can be seen in Figure 8. The elevator and rudder deflections match well the recorded values whereas computed deflections of aileron, although also close to recorded values, exhibits some oscillations. This coincides with oscillations of roll angular velocity p and roll angle φ. The reasons of roll rate oscillations are not clear. They cannot be attributed fully to the lack of its measurement, because other not measured states do not exhibit oscillations. It is thus a suspicion that its origin lies in improper tuning of Q matrix that results in middle instability of the nonlinear closed loop system (aircraft-LQR controller).
Comparing to the similar results obtained for full state tracking with LQR [6], one can conclude that the overall stability and sensitivity of the output tracking problem with assumed measurements, although degraded to some extent (which was expected) does not significantly loose its ability to accurately track the prescribed flight path.
In a conclusion, results described above indicate that the generalized LQR control used for output tracking problem is stable. Moreover, it provides accurate tracking for the given path in the absence of external perturbations e.g., due to atmospheric conditions.

Influence of Atmospheric Disturbances
To investigate the robustness of the extended LQR control in an output-tracking problem, the influence of disturbances was examined. There are two types of disturbances: modelling and physical ones. Modelling errors are considered as exogenous disturbances [40] and are assumed neglectable, because the aircraft model was identified with high fidelity. Physical disturbances are caused by atmospheric phenomena, such as wind or turbulence. The perturbations concern both state and environment and they are usually diverse. Typical errors of disturbances are presented in Table 2. The first test concerns checking the robustness of LQR control with respect to constant disturbances. Steady horizontal right-side wind with velocity 5 m/s was introduced. No turbulence has been added. Results are presented in Figure 9. There are no significant changes in LQR control efficiency when compared to outcomes presented in testing it with no disturbances. Increasing the wind velocity results in an approximately proportional decrease in tracking accuracy. The aircraft stability was affected only moderately. The selective oscillations of the roll rate and roll angle can be observed. Higher vertical velocity error and the Z coordinate error are caused by the changes in flight conditions introduced by wind that was not present in recorded flight. However, too high wind velocity or its adverse direction (tail wind), destroy stability and tracking completely and causes quick departure of an aircraft from the prescribed flight path. Maximal level of wind velocity that does not destroy tracking depends strongly on its direction and is roughly 20 m/s for head wind, 10 m/s for side wind and 5 m/s for tail wind.
The control errors are presented in Figure 10. It can be seen that the controls computed by LQR are less accurate when compared to the recorded ones. This is caused by the wind influence that was not present in recorded of flight. Higher level of aileron oscillations comparing to that of undisturbed case may also be observed. In conclusion, it can be stated that steady wind does not cause a significant change in tracking accuracy and aircraft flight stability. Therefore the LQR control used in the output-tracking problem is robust enough with respect to the constant disturbances of moderate level.
The second test concerns checking the LQR control robustness regarding the stochastic type perturbations. Random turbulence with zero mean and variance proportional to velocity components in Table 2 was added to the measurement vector y. Wind has not been added in this test. Results are presented in Figure 11. As in the previous test with steady wind, there are no significant changes in LQR control efficiency when compared to the outcomes shown in testing it with no disturbances. The main difference is that oscillations are present in all states. Oscillations of the roll rate and roll angle are higher than those of other states. Increasing the intensity of turbulence results in roughly proportional increasing oscillations of states, but tracking accuracy and aircraft flight stability was affected only moderately. However, too-high level of turbulence intensity destroys stability and tracking completely and causes quick departure of an aircraft from the flight path. The maximal level of turbulence intensity that did not destroy tracking was of order 3 m/s. The controls computed by LQR are presented in Figure 12. It can be seen that they are largely distorted comparing to the recorded ones, although their mean values are close to these recorded in flight with no turbulence. A higher level of aileron oscillations is present compared to that of the undisturbed case. It is worth of notice, however, that oscillations of control computed by LQR are usually quite significant, even for problems with smooth states. In conclusion, it can be stated that moderate turbulence does not cause a significant distortion of tracking accuracy as well as stability of aircraft flight. Therefore, it can be concluded that the LQR control used in output tracking problem is sufficiently robust with respect to stochastic type disturbances with zero mean values.

Conclusions
In this paper, an LQR controller was designed for an aircraft-tracking problem when some states were not observed. The controller design was based on a nonlinear aircraft model obtained from system identification when measurement and process noise was present in the data. The time domain filter error method based on the maximum likelihood principle allowed us to obtain accurate estimates of aerodynamic derivatives and the model had good predicting capabilities.
It has been shown that the LQR control can be used for handling of aircraft output tracking problem. The design was robust enough with respect to disturbances of moderate level. Even in the case of incomplete and inaccurate measurements, the output LQR controller has the ability of stable tracking the flight trajectory with quit good precision. Such abilities of the LQR depend essentially on the number and types of measured states. Too-poor measurements result in destabilizing the aircraft motion. The tracking accuracy decreases also because environment parameters were not estimated and thus not included in the LQR controller design.
In case of ideal atmospheric conditions, the aircraft tracking was sufficiently accurate and the local errors were small. This was true also for the unobserved states, with roll rate exception, for which considerable oscillations were observed. The LQR generated controls were in good agreement with their reference values. If the wind was added, the tracking accuracy dropped proportionally to the wind velocity. The LQR generated controls did not match reference values because the wind was not present in the reference flight. The controller was found to be robust if the wind was not too strong and did not act in diverse directions. Similar conclusions were drawn when turbulence was present in the experiment. Increasing turbulence had a proportional effect on all state variables oscillations, but again it was possible to provide accurate tracking for small and moderate disturbances.
It was possible to obtain accurate tracking, i.a., because of the high quality of the system identification model. For the model estimated with lower fidelity, the tracking quality would drop. This would be even more visible in case of fewer states present in the LQR static output feedback as less information would be available and thus the inaccuracies of the estimates would have greater overall contribution.
The theory concerning linear systems with static output feedback does not deliver conditions for the existence of stabilizing gain matrix that are easy to be checked-it is still a big challenge. However, it is possible to propose such a matrix for specific systems, e.g., as in our case, where the gain matrix being the solution of a standard LQR problem (that is with state feedback) together with the dynamic system matrices A, B, C satisfies 2n simple scalar inequalities.
A success in using LQR control for nonlinear output tracking problem depends essentially on using a high fidelity flight mechanics model that compensates unmodelled aircraft dynamics due to, e.g., uncertainties or excitations shape degradation. The output LQR control is a useful instrument in the estimation and tuning process of highly nonlinear aircraft models from flight data, e.g., for high fidelity full flight simulator certification.
Enhancing of accuracy or stability of the LQR output feedback control depends essentially on internal stability of the object and, even more, on the number and types of available measurements. The brand new idea in this article consists of using LQR weighting matrices Q and R in the first step of a truncated iteration used for computing the exact gain matrix for output feedback [43]. The possibility of using the more advanced algorithm of such type for further accuracy enhancing is currently investigated.
In future, we plan to investigate tracking accuracy degradation due to state unavailability, weighting matrices auto-tuning and compare the outcomes with the results obtained for nonlinear quadratic regulator. We expect that nonlinear quadratic regulator will improve the tracking accuracy and that decrease in information content stored in states measurements (e.g., related to their unavailability) will degrade the trajectory tracking. However, this needs to be verified through at least numerical simulations. The final step of this research is to verify the results through hardware in the loop simulations, followed by high-fidelity full flight simulator tests and, finally, experiments performed on true aircraft.
Funding: This research received no external funding.

Acknowledgments:
The authors would like to thank German Aerospace Center (DLR) for providing the BACM model.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following symbols were used in the manuscript: probability density function p, q, r roll, pitch and yaw ratē q dynamic pressure r reference value t, t k time and time at k-th point u, v, w longitudinal, side and vertical velocity u control vector x, y, z state, output and measurement vector x H horizontal tail arm x 0 initial state A, B, C state, control and output matrices C X , C Y , C Z longitudinal, side and vertical force coefficients C l , C m , C n rolling, pitching and yawing moment coefficients C L , C D lift and drag coefficients C L 0 , C Y 0 , C D 0 lift, side-force and drag coefficients bias C l 0 , C m 0 , C n 0 rolling, pitching and yawing moment coefficients bias C i j i-th force or moment coefficient change due to j-th flight parameter or control C i jk i-th force or moment coefficient change due to k-th flight parameter effect on j D triangular matrix F, G process and measurement noise distribution matrices H auxliary matrix I, I ij inertia matrix and its component about i axis when the object is rotated about j axis J cost function J RMS , TIC root mean square match and Theil inequality coefficient K s gain matrix in LQR problem with state feedback K * s optimal gain matrix in LQR problem with state feedback