A Unified Approach to Nonlinear Dynamic Inversion Control with Parameter Determination by Eigenvalue Assignment

This paper presents a unified approach to nonlinear dynamic inversion control algorithmwith the parameters for desired dynamics determined by using an eigenvalue assignment method, which may be applied in a very straightforward and convenient way. By using this method, it is not necessary to transform the nonlinear equations into linear equations by feedback linearization before beginning control designs.The applications of thismethod are not limited to affinenonlinear control systems or limited tominimum phase problems if the eigenvalues of error dynamics are carefully assigned so that the desired dynamics is stable.The control design by using this method is shown to be robust to modeling uncertainties. To validate the theory, the design of a UAV control system is presented as an example. Numerical simulations show the performance of the design to be quite remarkable.


Introduction
In the development of high-performance aircraft, control difficulties may be encountered over some parts of flight envelope.These difficulties arise from highly nonlinear aerodynamic properties [1] in some flight conditions.In order to solve these control difficulties, nonlinear controllers are required for high-performance aircraft.
Among many control methods, nonlinear dynamic inversion (NDI) is very popular and has been widely studied for flight control designs (e.g., [2,3]).NDI-based control systems are usually divided into fast-and slow-loop control subsystems according to the multiple time-scale method [4].In each subsystem, Lie derivatives [5] are used to transform the nonlinear equations into linear equations.Then, linear control design methods can be employed and the control inputs are obtained by converting the linear system control variables into the original coordinates.However, the control systems obtained by feedback linearization [6] may have nonminimum phase problems for affine or nonaffine nonlinear system [7] and robust issues in case of model mismatch.A typical nonminimum phase problem may be found in flight dynamics where the altitude-elevator transfer function usually has a right-half zero.The internal state control [8] is often used to overcome these nonminimum phase problems.In addition, the fuzzy logic control [9] was also applied for solving these kinds of problems.Furthermore, to overcome the robust problems, -analysis [10] and  ∞ method [11] were applied.Specially, incremental NDI (INDI) [12] was used to increase the robustness to aerodynamic uncertainties by calculating the control surface deflection changes instead of giving inputs directly.
To circumvent some aforementioned robust problems, an adaptive nonlinear model inversion control [13] was introduced, in which the design concept is similar to the conventional NDI yet without linearizing the nonlinear system.The model inversion method replaces the motion rates with a P-form or PI-form desired dynamics to negate the original dynamics.The choice of parameters in the desired dynamics is based on the bandwidth of response and time scales.The effects of different types of desired dynamics on the resulted control system were discussed [14].
Although the aforementioned NDI approaches are successful in many flight control system designs [14][15][16] over a large part of the flight envelope, the systems of dynamics in general have to be separated into several subsystems according to the rates of response.There are many cases in which the fast rate and the slow rate might not be distinguished so clearly however.Also, although the pole assignment method had been introduced to determine the parameters in the desired dynamics, very often the control system must first be transformed into a standard feedback control form from which the standard eigenvalue assignment method can be applied.
In consideration of the aforementioned problems existing in the current literature on control designs, a unified approach to nonlinear dynamic inversion control is proposed in this paper.The equations of motion will not be necessary to be separated into fast rate and slow rate groups, nor will they be limited to an affine system.Feedback linearizations will not be required to transform the nonlinear equations into linear equations.Nonminimum phase problems are solved by eigenvalue assignments for error dynamics.An iterative method for determining the parameters of the desired dynamics from the assigned eigenvalues of error dynamics is proposed.Analysis of robustness to model uncertainties or disturbances is conducted.This method will be ready for design without simplifying the system of equations based on physical insights once the governing dynamic equations are established and the state variables to be tracked are selected.The theory is to be developed in detail in the following sections.A UAV is introduced and its control system is designed with the developed method.Numerical simulations are conducted to validate this method.

Nonlinear Dynamic Inversion Control.
In general, dynamic equations of motion with control inputs can be expressed by where x ∈   is the state vector, u ∈   ( < ) the control vector, and f the nonlinear function representing the model of dynamics with controls.By extending the concept of dynamic inversion (DI) [4], the control vector u can be assumed to be computed from where x  ∈   is a desired state vector with its rate of change being designated.In this paper, the desired dynamics is designated as a set of stable first-order differential equations: where Ω ∈  × represents a constant matrix with  independent parameters which can be chosen.Substituting (3) into (2) yields which constitutes a set of  algebraic equations.Since  < , obviously, u cannot satisfy (4) if all elements of x  are to be designated.It means that only part of x  can be designated.
So let x  be divided into two groups, say x  ∈   and x  ∈  − , where x  contains the state variables which are to be controlled or designated and x  the residual ones.Both u and x  constitute  unknown variables which are to be determined from (4).To solve a set of nonlinear algebraic equations, the Newton-Raphson iteration method can be employed as follows: where f err ≜ f(x  , u) − Ω(x  − x) and  = 0, 1, 2, . . .represents the iteration number.

Parameter
Determination.Now, a question arises whether the state vector x will asymptotically follow the desired vector x  if u is determined from (4) and substituted in (1).
In order to answer it, let (1) and ( 4) be examined more carefully as follows.Subtracting (4) from (1) yields If x is very close to x  , then the above equation can be linearized as follows: With the concept that the desired variables x  are near constant, say ẋ  ≈ 0, (7) can be rewritten as Defining an error vector e = x − x  and replacing the approximate sign with the equal sign lead (8) to Equation ( 9) is a set of error dynamics in which the error vector e will vanish eventually if all the real parts of the eigenvalues of [−Ω + f  (x  , u)] are negative.This is possible if Ω is chosen appropriately.It means that the state vector x can approach the desired vector x  once e approaches 0. Recall that x  contains x  , a state vector to be tracked.Notice that the eigenvalue  in (9) can be determined by where I is an identity matrix.For simplicity, if Ω is diagonal and all elements are the same, say, Ω = I, then (10) can be rewritten as In order to make e bounded,  can be so chosen that real (  + ) < 0 (12) for all   ( = 1, 2, . . ., ).Although this method is simple, the resulting  may be unnecessarily large.
For more general cases, Ω contains a set of parameters,   , ( = 1, 2, 3, . . ., ), which may be chosen to determine the eigenvalues and eigenvectors of −Ω + f x (x  , u).Recall that (3) must be a stable model and, therefore, the simplest way of constructing Ω is to let all its elements vanish except those at the diagonal, which are assumed to be Ω  =   > 0. In fact, some off-diagonal elements can also be allowed to exist.For example, let Ω  = Ω  =   > 0, and Ω 2  = Ω 2  =  2  <  2  .It is trivial to prove that the latter case is also a stable model.Now, assume that the th ( = 1, 2, 3, . . ., ) eigenvalue and eigenvector are   and e  , respectively.Accordingly, In general, the eigenvector e  can be normalized so that However, if the eigenvalue   is desired rather than   , then it is necessary to adjust the parameters in Ω. Assume that an increment (Ω/  )Δ  to Ω is required.Then where Accordingly, the derivative of (13) with respect to   results in and the derivative of ( 14) becomes Equations ( 17) and ( 18) can be rearranged to from which   /  can be determined.Note that   /  represents the variation of   with respect to   .With all   /  , ( = 1, 2, 3, . . ., ), the th revised eigenvalue should be Recall that, in this equation,   represents the th eigenvalue of where Ω contains the parameters   , ( = 1, 2, 3, . . ., ).Now since   is the th desired eigenvalue, (20) in fact becomes a set of  simultaneous equations from which Δ  can be solved.With   being the initial guess, revisions for parameters can be made by Since   , ( = 1, 2, 3, . . ., ) are nonlinear function of   , ( = 1, 2, 3, . . ., ), iterative computations of ( 9) and ( 19)-( 21) are required in order to obtain a set of convergent   .To make it clear, the iteration procedures are summarized as follows.
(3) Determine the eigenvalues and eigenvectors of ( 9).Denote the th eigenvalue and eigenvector as   and e  , respectively.
(7) Replace   with   and go to step 3.
At this point, it must be mentioned that the iterations will converge only if the initial guesses are very close to the true answers.It is also emphasized that the desired eigenvalues   must be chosen to lie on the left-half -plane and the resulting parameters in Ω must satisfy the stable model requirements in (3).

Robust Analysis.
The nonlinear dynamic inversion control design developed so far is based on a nominal dynamical model.In the real world, however, there always exist some modeling uncertainties which cannot be determined in advance.In order to check if the control design is robust, let the actual dynamical model be represented by With this assumption, now (6) can be modified to which can be rewritten as follows: where Δf  (x, x  , u) denotes the nonlinear part of f(x, u) − f(x  , u).The solution of the above equation can be represented as follows: Note that from ( 9), the eigenvalues of −Ω + f  (x  , u) are all in the left-half -plane since Ω has been determined with that assumption.Also, at this point, it is not unreasonable to assume that both Δf  (x(), x  , u()) and the modeling differences |f  (x, u) − f(x, u)| are bounded.Therefore, if the model in ( 3) is stable enough, the integral part in (25) will vanish along with  −Ω .Accordingly, as the time  gets large enough, the state variables x will approach the desired values x  as can be found from (25) even if there are some modeling uncertainties or disturbances.

Flight Dynamics Equations of Motion.
To illustrate the theory, a design of flight control system with the method developed is presented.Before the flight control design proceeds, a set of flight dynamics equations of motion must be formulated.Note that all aerodynamic forces and moments result from the relative motions between aircraft and the air.The aircraft is assumed to have a ground velocity: where (i, j, k) are unit vectors in the aircraft body moving frame and (I, J, K) the unit vectors in a fixed flat earth frame.The air is assumed to have a velocity: which is also known as the wind velocity.Then, the velocity of the aircraft relative to the air can be represented by from which, the aircraft total speed relative to the air, the angle of attack, and the sideslip angle can be determined, respectively, by the following equations: With the assumptions of fixed flat earth and winds being present, the motions of aircraft with six degrees of freedom can be represented by a set of nonlinear first-order differential equations as follows: where In (44),    ,    , and    are the three components of the force exerted by winds.At this point, it is worthy to mention that although there is no explicit term relating wind  38), but also have implicit effects on them through , , and  which obviously depend on   , , and .
To validate the method developed in this paper, a UAV as shown in Figure 1 is introduced.The parameters used for analysis are listed in Table 1.
The aerodynamic forces and moments of the UAV are computed by the following equations:  where  is the dynamic pressure and  the Mach number.
The aerodynamic coefficients and stability derivatives are determined by a computer code dubbed "VORSTAB" [17].
The computation results are in the form of discrete data which are then interpolated as continuous functions in the computer program.These functions are represented by   (, ,   ),    (, ), and so forth, in (45).The nonlinear database includes  ranging from −20 ∘ to 40 ∘ ,   from −24 ∘ to 24 ∘ ,   from −25 ∘ to 25 ∘ , and   from −30 ∘ to 30 ∘ .Also, the flight conditions include altitude ℎ ranging from sea level to 40, 000 ft and Mach number  from 0.1 to 0.9 and from 1.1 to 1.9.The thrust model of the UAV is assumed as where   is the thrust command.The thrust is assumed to have a limitation, say 0 ≤  ≤ 1.1.

Flight Control Design.
In flight control design, only nominal flight dynamics models are considered since modeling uncertainties or wind disturbances, and so forth, cannot be determined in advance.Hence, all terms related to the wind disturbances in the flight dynamics equations of motion are neglected in this stage.For this flight dynamics model, the parameters involved in the control analysis are further elaborated as follows.
Intuitively,   ( = 1, 2, 3, . . ., 9) ∈  may just be arbitrarily assigned as long as they satisfy the stable model requirements.However, the resulted error dynamics in ( 9) may not just be stable.As will be shown in simulations, the desired eigenvalues must not only be able to make the system follow the desired dynamics, but also be able to generate a good set of  which makes  −Ω decrease so quickly that the system is also robust if modeling uncertainties exist or wind disturbances are encountered.The choice of these eigenvalues is not only just based on the control analysis alone but must also be simultaneously based on simulations.
Note that, in this case, longitudinal dynamics analysis also reveals a right-half zero  = 10.85 in the ℎ-  transfer function.It can therefore be identified as a nonminimum phase problem.For this kind of problem, the desired eigenvalues must be very carefully assigned lest the resulting  does not satisfy the stable model requirements.To determine , a twoway approach (  ) is proposed as follows.
(1) If the desired eigenvalues are equal to the open-loop eigenvalues as listed in Table 2, obviously,  = 0.
Choose small ( ̸ = 0) which satisfies the stable model requirements.
(3) If the error dynamics is not stable, modify  so that it makes the error dynamics stable.Use the modified values as the desired eigenvalues.
(4) Determine  by following the 7 steps of iteration procedure described in Section 2.2.
(5) If  does not satisfy the stable model requirements, modify  to make it satisfy.Go to step 2.
(6) In each step of simulations, determine x  and u from (4) by using the iteration method described in (5).In this step, only the nominal dynamics model is used.
(7) Use u in (1) for simulations.The equations of motion may include modeling uncertainties or wind disturbances.
(8) If the simulation results are not satisfactory, modify  and then go to step 2, or modify  and then go to step 4. By using these 8 steps of procedure, the designated eigenvalues and the resulting parameters for desired dynamics are obtained and listed in Tables 3 and 4, respectively.
These parameters indeed satisfy the stable model requirements.Note that initially  5 is small and  1,2 and  3,4 are not so deviated from their counterparts of open-loop eigenvalues.However, the simulation results show that the performance in tracking ℎ  is not good enough.Enlarging  5 does improve the performance but makes  1,2 and  3,4 deviate their counterparts of open-loop eigenvalues a lot.Also note that the large deviations of  6 ,  7,8 , and  9 from their counterparts of open-loop eigenvalues are due to the iteration procedure in determining a set of  6 ,  7 ,  8 , and  9 for satisfying the stable model requirements.
At each instant, ℎ  ,   ,   , and   are given, the controls   ,   ,   , and   along with the residual state variables   ,   ,   ,   , and   can be determined from (4) by using the iteration method described in (5).Note that, in the computation of the control, only nominal flight dynamics equations of motion are used since modeling uncertainties or wind disturbances are not known.Also note that, in using the Newton-Raphson method, the convergence can be guaranteed if the initial guess is sufficiently close to the solution.Since the solution changes only very little in each step of integration, numerical practices show that it takes only 4 iterations to converge to within 0.0001% of the correct value if the solution in the previous step of integration is taken as an initial guess.In the first step of integration, it may be necessary to take a few more iterations, say, 10 iterations, since the solution is different angle is closely associated with the heading angle as the bank angle is computed based on (48).Although the wind disturbances do not explicitly affect φ in (35), they do implicitly affect it through  and  which in turn affect the aerodynamic moments.Fluctuations in the bank angle and the heading angle are remarkable when the UAV passes through the wind zone but are still tolerable.
Sideslip angle seems to be affected by winds more heavily as shown in Figure 4(e).Fortunately, the control system is robust enough as being able to keep it small.As shown in Figure 5(a), all control surface deflection angles remain small.The elevator needs only to deflect slightly to rotate the UAV to climb.As mentioned early,   = −0.733∘ at the instant when ℎ  and   are given.When the UAV reaches the commanded state conditions, the elevator trim angle approaches −0.510 ∘ .Small deflection in aileron angle is enough to make the UAV bank turn and the rudder just keeps very small as does the sideslip angle.All control surface deflection angles are not remarkably affected as the UAV passes through the wind zone.
As shown in Figure 5(b), the thrust command increases very sharply as a demand to increase both the altitude and the Mach number simultaneously is given.Also, when the UAV passes the wind zone, the thrust command fluctuates very violently.Obviously, the thrust is very closely interacted with the Mach number.A low pass filter with time constant 10 sec alleviates the sharp and violent responses a lot for the actual thrust at the expense of delaying its reaction time.
In Figures 6(a)-6(e), responses of all residual state variables along with the computed commands for them are presented.In fact, these commands are generated in the inner system, not input from outer designations.It is interesting to observe how the actual state variables track the commands.
In Figure 6(a), it is observed that, in comparison to  trim = 6.04 ∘ , the computed command   = 4.85 ∘ at the instant when ℎ  and   are given.Then the angle of attack  follows   closely without apparent short-period mode oscillations when the designated short-period mode eigenvalues are as high as  1,2 = −4.22180± 1.57684.Winds do have remark effect on both   and .Finally, both approach a new value,  trim = 2.88 ∘ , for flight conditions at higher altitude and faster speed.
In Figures 6(b) and 6(c), it is observed that although both   and   increase in response to the requirement for increasing the cruise altitude,  and  decrease remarkably.As revealed in explaining ℎ response, the elevator angle computed is not enough to hold q = 0 and therefore q < 0 which in turn makes θ < 0. The minimum pitch rate is  = −0.717deg./s at  = 0.78 sec.In this case, the vertical wind does have more remarkable effects on  and  than the horizontal wind.Finally, both approach  trim =  trim = 2.88 ∘ and  trim = 0, respectively, for the new flight conditions.
It is very interesting to find that, in Figures 6(d) and 6(e), the shapes of  and  responses resemble, respectively, those of  and  responses as shown in Figures 4(d) and 4(e).The effects of horizontal wind on both responses are more remarkable than those of vertical wind in this case.The decay of  seems to be very slow.

Conclusions
In this paper, a unified nonlinear dynamic inversion control system is successfully developed.With this method, the Pitch rate (deg./s)q (with wind) q r (with wind) q (no wind) q r (no wind) Mathematical Problems in Engineering 13 parameters for desired dynamics can be determined with a set of assigned eigenvalues for error dynamics.The control system has been proved to be robust to modeling uncertainties or wind disturbances.A nonaffine nonlinear flight dynamics system with right-half zero has been used as an example for the control design.In the design process, it is not necessary to use the feedback linearization to transform the nonlinear equations into linear equations.Numerical simulations of the control system show that the desired state variables can be successfully tracked whether winds are present or not.

Figure 1 :
Figure 1: The configuration of the UAV.
(no wind) r (no wind)  (with wind)
and e  are assumed to be approximated to   + (  /  )Δ  and e  +(e  /  )Δ  , respectively, and e 11 cos  cos  +  21 sin  +  31 sin  cos )   ( 12 cos  cos  +  22 sin  +  32 sin  cos ) = −  ( 13 cos  cos  +  23 sin  +  33 sin  cos )   +   −   )     −      +     −     −         −      +       −       −       −           −     +         −      +   (  −   +   )     −      −       −       +       −       , (41) where (  ,   ) represents the position in a fixed flat earth frame, ℎ the altitude,   the elements of direction cosine matrix for transferring a fixed flat earth frame to the aircraft body moving frame,  the angle of attack,  the sideslip angle,  the heading angle,  the pitch angle,  the bank angle,  the roll rate,  the pitch rate, and  the yaw rate.Also,   ,   , and   represent, respectively, three components of the total force in an aircraft body moving frame.The three force components are composed of the thrust , the lift , the drag , the side force , and the weight  ( is the aircraft mass and  the gravity acceleration) by the following equations:  is the angle between the thrust and the longitudinal axis.Moreover,   ,   , and   represent the roll moment, the pitch moment, and the yaw moment, respectively, about the center of gravity, and   ,   ,   , and   are the components of the moment-of-inertia tensor.Furthermore, the wind disturbances on V  , α , and β are, respectively, represented by    ,    , and    which are expressed as follows: =  =  sec  sin  +  sec  cos  (33) θ ≜   =  cos  −  sin  (34) φ ≜   =  +  tan  sin  +  tan  cos (35)= − cos  tan  +  −  sin  tan  =   =  cos   +  sin  −  cos  −  sin    =  +  sin  cos    = − sin   −  cos  −  sin  +  cos  cos ,(42)where

Table 2 :
The open-loop eigenvalues of f x (x trim , u trim ).

Table 3 :
The designated eigenvalues.

Table 4 :
The resulting parameters for desired dynamics.