Modified Infinite-Time State-Dependent Riccati Equation Method for Nonlinear Affine Systems: Quadrotor Control

This paper presents modeling and infinite-time suboptimal control of a quadcopter device using the state-dependent Riccati equation (SDRE) method. It establishes a solution to the control problem using SDRE and proposes a new procedure for solving the problem. As a new contribution, the paper proposes a modified SDRE-based suboptimal control technique for affine nonlinear systems. The method uses a pseudolinearization of the closed-loop system employing Moore–Penrose pseudoinverse. Then, the algebraic Riccati equation (ARE), related to the feedback compensator gain, is reduced to state-independent form, and the solution can be computed only once in the whole control process. The ARE equation is applied to the problem reported in this study that provides general formulation and stability analysis. The effectiveness of the proposed control technique is demonstrated through the use of simulation results for a quadrotor device.


Introduction
The state-dependent Riccati equation (SDRE) method has become more popular due to the possibility of designing more effective nonlinear controllers and providing flexibility through state-dependent weighting matrices, which is well presented in [1]. This method scheme originates from the LQR method [2,3] and is also similar to the MPC (or its nonlinear version, NMPC) [4][5][6]. The main differences between SDRE and those two methods are, for example, unlike the LQR method, the SDRE sets up its control laws at every time-step. The MPC method requires a convex optimization technique, while SDRE and LQR methods utilize optimal control laws for linear systems [7,8].
The SDRE control method reformulates the nonlinear dynamics using parameterization to a linear structure with state-dependent coefficient (SDC) matrices, which are essential to solve an algebraic Riccati equation (ARE) and obtain the suboptimal control law. The coefficients of the ARE equation differ in the state space. Thus, at a given point in state space, it is essential to solve an algebraic state-dependent Riccati equation. The nonuniqueness of the state-dependent parameterization (SDP) allows us to obtain extra degrees of freedom, which may help to enhance controller performance. Minimization of the nonlinear performance index gives a quadratic-like structure [9][10][11][12].
The main disadvantage of this classic SDRE strategy is that Riccati equations are state dependent. Therefore, it is necessary to solve SDRE multiple times during the control process [25,26]. Nowadays, the solution can be successfully implemented in real control systems due to developments of digital controllers. However, the solution can be found faster, reducing computational effort and introducing linearization of closed-loop form. As a new contribution, the paper presents a possible solution for suboptimal control of nonlinear systems by solving the algebraic Riccati equation only once in the whole control process. Using linearization of the closed-loop system and introducing Moore-Penrose inversion [27], it is possible to modify SDRE into ARE. For the modified SDRE method, it is possible to obtain the same solution as for classic SDRE method. The advantage of the modified case is that there is no need to solve the state-dependent Riccati equation for each state of the control process. This allows us to save space for implementation in real-time control systems and simplifies computations. The specific differences between classical and modified methods are presented in the sections below.

SDRE Method
Let Equation (1) be the continuous time nonlinear systeṁ where x ∈ R n and u ∈ R m are the state and input vectors, respectively. F(x) is a continuous (C k class) nonlinear function of x, and B ∈ R n×m is the constant control matrix. The problem consists of finding an admissible control u ∈ R m for t 0 ≤ t, which minimizes the performance index [28][29][30][31] where Q ∈ R n×n is a symmetric positive semidefinite, and R ∈ R m×m is a symmetric positive-definite matrix. Rewriting the nonlinear Equation (1), the state-dependent coefficient form (SDC) followsẋ where A(x) ∈ R n×n . There are many methods for SDC parameterization. Each parameterization should be true for all 0 ≤ α ≤ 1 [32] αA The control problem for the nonlinear continuous-time system Equation (1) can be formulated as follows. Given nonlinear functions F(x) ∈ R n , B ∈ R n×m , Q ∈ R n×n , and R ∈ R m×m , we attempt to find a control vector u ∈ R m for t ∈ [t 0 , ∞] that allows us to control the system state vector from x 0 to x ∞ while the performance index is minimized (Equation (19)).
If the integral functions x T Qx + u T Ru and A(x)x + Bu are continuously differentiable functions of each of their arguments, then we may imply that u ∈ C[t 0 , ∞] is a control that minimizes the functional J(u) : C[t 0 , ∞] → R. To solve the problem, consider the Hamiltonian where p ∈ R n is the co-state vector.
To be consistent with the Pontryagin minimum principle, if u ∈ C[t 0 , ∞] is a control for the functional Equation (19), subject to state Equation (3), and if x denotes the corresponding state, then there exists a p ∈ C[t 0 , ∞] such that Equations (6) and (7) are conditions allowing us to obtain suboptimal control while minimizing Equation (19). It follows that any optimal input u ∈ R m and the corresponding state x ∈ R n satisfies Equation (6); that is, From Equation (8), the optimal control vector is and subsequently, the adjoint differential equation iṡ Let p be a nonlinear combination of the state of the system where K(x) ∈ R n×n denotes feedback compensator gain, and let x be the solution of the nonlinear state equationẋ for t ∈ [t 0 , ∞], x(t 0 ) = x 0 . Then, feedback control is Let K(x) solve the state-dependent Riccati equation (SDRE), which is then the equation below satisfies the optimality condition [25,32] The suboptimal control for performance index (objective function) Equation (19) subject to dynamics Equation (1) can be found solving the state-dependent Riccati Equation (15). Often, it is difficult to find the solution of Equation (15) analytically. One approach allows us to solve the SDRE via symbolic software packages such as Matlab. However, for complex systems, the solution may become complicated, and then, it is necessary to approximate the solution.
To approximate, an interpolation method or Taylor series method can be used, for instance, ref. [25].

Modified SDRE Method
Let Equation (18) be the continuous time nonlinear systeṁ where x ∈ R n and u ∈ R m are the state and input vectors. F(x) is a continuous (C k class) nonlinear function of x, and B ∈ R n×m is the constant control matrix. The problem consists of finding an admissible control u ∈ R m for t 0 ≤ t that minimizes the performance index [28][29][30][31] where Q ∈ R n×n is a symmetric positive semidefinite, and R ∈ R m×m is a symmetric positive-definite matrix.
The first modification (compared to the classic approach) is to rewrite the system Equation (18) as a sum of state-independent and nonlinear state-dependent coefficient (SDC) forms Equation [25]: where is a sum of a constant matrix and a state-dependent matrix [12,33]. The SDC parameterization of Equation (20) can be performed the same way as for Equation (4), for all 0 ≤ α ≤ 1 [32]: To easily formulate controllability criteria, the state-dependent controllability matrix is introduced W(x) [25] If W(x) (state-dependent) has the full rank, then the system is controllable for all x ∈ R n [34]. As in the previous case of the SDRE method, an admissible control u(t) ∈ R m , that minimizes the performance index can be found for [28][29][30][31] where Q ∈ R n×n is a symmetric positive semi-definite matrix, and R ∈ R m×m is a symmetric positive definite matrix. The control problem for the nonlinear continuous-time system Equation (18) can be formulated as follows. Given nonlinear functions F(x) ∈ R n , B ∈ R n×m , Q ∈ R n×n and R ∈ R m×m , we attempt to find a control vector u ∈ R m for t ∈ [t 0 , ∞] that controls the system state vector from x 0 to x ∞ while minimizing the performance index (Equation (19)).
If the integrand of the cost functions x T Qx + u T Ru and A(x)x + Bu are continuously differentiable functions of each of their arguments, then we may imply that u ∈ C[t 0 , ∞] is a control that minimizes the functional J(u) : C[t 0 , ∞] → R. Note that, in this approach, two feedback compensators are defined A(x) = A 1 + A 2 (x). To solve the problem, consider the Hamiltonian where p ∈ R n is the co-state vector.
If u ∈ C[t 0 , ∞] is a control that minimizes Equation (19) subject to the state Equation (20), and if x denotes the corresponding state, then there exists a p ∈ C[t 0 , ∞] such that where p satisfies the condition for optimal control that minimizes Equation (19). It follows that any optimal input u ∈ R m and the corresponding state x ∈ R n satisfy Equation (26); that is, Thus, the optimal control is and subsequently, the adjoint differential equation iṡ Consider p as a linear and nonlinear combination of the state of the system Equation (18), which is the second modification of the classic approach where K(x), K 1 , K 2 (x) ∈ R n×n . Consider x as the solution of the nonlinear state equatioṅ Note that, in this approach, two feedback compensators are defined. Rearranging Equation (33), it is possible to obtaiṅ The first bracket of Equation (34) is state-independent, and the second one is statedependent; thus, it is possible to linearize it and solve the state-dependent gain matrix- as follows: Matrix BR −1 B T is singular; thus, the state-dependent matrix gain K 2 (x) may be computed only by the pseudoinverse operation [. . .] + .
In order to perform this procedure, a Moore-Penrose pseudoinverse is used, and the matrix is unique [27,35,36].
Equating the adjoint differential Equation (30) and the differential form of Equation (32), we have Consequently, employing Equations (22) and (29), it is possible to obtain and rearranging terms in Equation (39), we find Based on linearization Equations (34)-(36), let the matrix K 1 solve Equation (39), the algebraic Riccati equation (ARE): then, the optimality condition is given bẏ for K 2 (x) obtained from (36). So, the suboptimal control For index Equation (19), subject to Equation (18), it is possible to find the solution for the state-independent Riccati Equation (39). In general, it is difficult to find the solution of Equation (39) analytically, especially for complex systems [22,25,33,37]. Symbolic software packages (such Matlab or Mathematica) are commonly used to solve such problems. However, for complex systems, the solution may become complicated, and then, it is necessary to approximate the solution. To approximate, an interpolation method or Taylor series method seems to be most suitable [15].

Stability Analysis
Asymptotic stability of the closed-loop system ensures the possibility of controlling the states from the initial values to the final ones. The controlled system with the SDRE compensator-based feedback is locally asymptotically stable [25]. Let r > 0 be the largest radius, such that B r (0) ⊆ Ω. Assuming that the system is stabilizable at x = 0, it is possible to use optimal control theory to define matrix K 1 , such that all eigenvalues of (A 1 − BK 1 ) have negative real parts. There exist β > 0, such that Re(p) < −β for all eigenvalues p of (A 1 − BK 1 ), having the systeṁ where (A 1 + A 2 (x))x and ∂(A 1 x + A 2 (x)x) ∂x are continuous in x for ||x|| < r, where r > 0 is the largest radius around the origin x = 0.
and h(x) = g(x)x; then, the system (34) is described bẏ Evaluation of h(x) shows that it is an almost linear system: from the inequality Thus, ||g(x)|| → 0 when ||x|| → 0; then, h(x) satisfies the condition Equation (46). The theoretical considerations of almost linear systems leads to the conclusion that, if eigenvalues (A 1 − BK 1 ) have negative real parts, h(x) is continuous around the origin, and the condition of Equation (46) holds, then x = 0 is asymptotically stable [38].
For the proposed method, a global uniform stability can also be proved using the Lyapunov function. The proof can be obtained more easily than for classic SDRE, because feedback compensator gains are state-independent, and for a closed-loop system, there is no problem with generalization and definition of an attraction region.
Assuming that SDC parameterization is stabilizable and detectable for all x(t), then the closed-loop matrix A CL is symmetric for all x(t), and the solution of SDRE is asymptotically stable when t ∈ [t 0 , ∞]. Let the control law be in the form and K 1 satisfies ARE Then, the system in closed-loop forṁ is globally asymptotically stable if K 1 is the unique symmetric positive-defined matrix. Let us define the as Lyapunov function and for x ∈ U, where U in R n has the origin, and constants k 1 and k 2 are obtained from and for i = 1, 2, ..., n . Using Equations (54) and (56), it is possible to obtaiṅ Taking Equation (56), assuming that K 2 (x) assures linearization of Equations (45) and (46), and that the solution of ARE, Equation (54), is always symmetric, theṅ Then, taking into consideration (54), it is possible to obtaiṅ Let us assume thatV where for x ∈ U and i = 1, 2, ..., n.
Generally, in an SDRE control problem with infinite time horizon, the existence of constants k 1 , k 2 , k 3 do not assure global stability. The closed-loop matrix is state-dependent but can be generalized to global stability by defining the region of attraction [39].
In the proposed method, the existence of constants k 1 , k 2 , k 3 proposed in Equations (59), (60), and (65) assures global stability, because the matrix K 1 is state-independent, so there is no reason to define the region of attraction for the Riccati Equation (54).

Numerical Simulation
To present a numerical example, the mentioned quadcopter model was chosen. Typical SDRE control technique was compared with its modified version.
Considering a quadcopter (presented in Figure 1) as a vehicle with four independent drives and with an electric power system located at its center of gravity, it is important to define coordinate systems in order to determine the position. A model device has six degrees of freedom, where vertical movement in the global coordinate system and three Euler angles are controlled parameters. There is a need to assume a body fixed frame and a symmetrical structure of the model with the origin in the center of mass, where every drive is independently controlled and aerodynamic effects are negligible. The mathematical model is taken from [26]. As mentioned, this model needs to be described using Euler's angles: roll angle φ ∈ −π, π , pitch angle θ ∈ − π 2 , π 2 , and yaw angle ψ ∈ −π, π . These angles are a sequence of three rotations This leads to the rotation matrix R, which describes rotations of the local system in terms of the global system.
The representation of the actual model in the simulations requires us to determine differential equations that describe all relations in the model. The orientation of the vehicle is described by two vectors, r T = (x, y, z) and Ω T = (ψ, φ, θ). The quadrotor's acceleration is described by the differential equation where: g-the gravitational acceleration; R-the rotation matrix; b-the thrust factor; ω i -the speed of rotor i = 1, 2, 3, 4. The second differential equation is presented in Vector τ is defined as where: l-the length the of quadrotor's arm; d-the drag factor. The quadrotor's model has four control inputs (vertical movement in the global coordinate system and three Euler angles). It is assumed that the input values, which represent forces, are directly proportional to the squared angular velocity of the motor.
Let us create a new variable that describes relative motor speed Using Equations (68)-(72), a dynamic model is presented as a system of six equations It is possible to describe the model in the state-space form, where is the control vector, and x ∈ R 12 is the state vector Combining Equations (73) and (75), we obtaiṅ (cos x 7 sin x 9 cos x 11 + sin x 7 sin x 11 ) u 1 m x 4 (cos x 7 sin x 9 sin x 11 − sin x 7 cos x 11 ) u 1 m x 6 −g + (cos x 7 cos x 9 ) u 1 m x 8 x 10 x 12 x 8 I 2 + J R I y x 8 ω d + l I y u 3 x 12 where Considering the attitude of the control, and applying SDRE method, the state vector and the control vector After factorization, it is possible to obtaiṅ The above state-space model has SDC structure, where A(x I ) is state-dependent matrix, and B is constant, as assumed in Equation (18). Applying the proposed methodology and replacing A by a sum of two matrices as presented in Equation (20) (for proposed method), it is possible to obtaiṅ To present the influence on the system, weight matrices were changed during this experiment. There were three different configurations: The values of parameters used in modeling of the object are presented below in Table 1 [40]. (82) Simulation time was 5s and time step was 1 × 10 −3 s. The aim of the control process is stabilization at "zero" point. The comparison of the execution time for both methods ( Table 2) shows that the proposed method reduces computational effort, which can be seen in the shortened time of computation below. The computations were executed on Windows 7 Professional, 64-bit Intel® Core™ i5-5200M CPU @ 2.40GHz.

s
Simulation results for configuration R 1 , Q 100 are presented. Differences between classical and modified SDRE method for next two configurations are unnoticeable.
In every case of Q and R matrices configuration after a short transition phase, all angles are stabilized at a required point (Figures 2-4). The waveforms for both methods are convergent. Changing the value of Q and R matrices affects only the speed of regulation. Small differences are due to numerical errors and pseudoinversion inaccuracy.

Conclusions
The modeling and control of a quadcopter using the SDRE method and modified SDRE-based proposition were presented in this paper. A new method for solving a suboptimal control problem is established, and a procedure for solving the problem was presented and illustrated with a numerical example. As shown, the modified SDRE method gives the same solution as the classic SDRE method, but in the modified case, there is no need to solve a state-dependent Riccati equation for each state of the control process. This allows us to save space for implementation in real-time control systems and simplifies computations.