Abstract
This paper presents a unified approach to nonlinear dynamic inversion control algorithm with 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 this method are not limited to affine nonlinear control systems or limited to minimum 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.
1. 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–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.
2. A Unified Approach to Nonlinear Dynamic Inversion Control
2.1. Nonlinear Dynamic Inversion Control
In general, dynamic equations of motion with control inputs can be expressed bywhere is the state vector, () the control vector, and the nonlinear function representing the model of dynamics with controls. By extending the concept of dynamic inversion (DI) [4], the control vector can be assumed to be computed fromwhere 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) yieldswhich constitutes a set of algebraic equations. Since , obviously, cannot satisfy (4) if all elements of are to be designated. It means that only part of can be designated. So let be divided into two groups, say and , where contains the state variables which are to be controlled or designated and the residual ones. Both and 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 and represents the iteration number.
2.2. Parameter Determination
Now, a question arises whether the state vector will asymptotically follow the desired vector if 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) yieldsIf is very close to , then the above equation can be linearized as follows:With the concept that the desired variables are near constant, say , (7) can be rewritten asDefining an error vector and replacing the approximate sign with the equal sign lead (8) toEquation (9) is a set of error dynamics in which the error vector will vanish eventually if all the real parts of the eigenvalues of are negative. This is possible if is chosen appropriately. It means that the state vector can approach the desired vector once approaches . Recall that contains , a state vector to be tracked.
Notice that the eigenvalue in (9) can be determined bywhere is an identity matrix. For simplicity, if is diagonal and all elements are the same, say, , then (10) can be rewritten asIn order to make bounded, can be so chosen thatfor all . Although this method is simple, the resulting may be unnecessarily large.
For more general cases, contains a set of parameters, , which may be chosen to determine the eigenvalues and eigenvectors of . 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 . In fact, some off-diagonal elements can also be allowed to exist. For example, let , and . It is trivial to prove that the latter case is also a stable model.
Now, assume that the th eigenvalue and eigenvector are and , respectively. Accordingly,In general, the eigenvector can be normalized so thatHowever, if the eigenvalue is desired rather than , then it is necessary to adjust the parameters in . Assume that an increment to is required. Thenwhere and are assumed to be approximated to and , respectively, and is also normalized so thatAccordingly, the derivative of (13) with respect to results inand the derivative of (14) becomesEquations (17) and (18) can be rearranged tofrom which can be determined. Note that represents the variation of with respect to . With all , , the th revised eigenvalue should beRecall that, in this equation, represents the th eigenvalue of , where contains the parameters , (). 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 bySince are nonlinear function of , 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.(1)Let the desired eigenvalues be which are all distinct.(2)Guess a set of parameters .(3)Determine the eigenvalues and eigenvectors of (9). Denote the th eigenvalue and eigenvector as and , respectively.(4)With each , determine from (19).(5)Determine from (20). If is less than a preset small value, then stop; else continue with the next step.(6)Determine from (21).(7)Replace with and go to step .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).
2.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 byWith this assumption, now (6) can be modified towhich can be rewritten as follows:where denotes the nonlinear part of . The solution of the above equation can be represented as follows:Note that from (9), the eigenvalues of 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 and the modeling differences 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 will approach the desired values as can be found from (25) even if there are some modeling uncertainties or disturbances.
3. Nonlinear Dynamic Inversion Flight Control System Design for a UAV: An Example
3.1. 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 are unit vectors in the aircraft body moving frame and 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 byfrom 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 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:where 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 , , and are, respectively, represented by , , and which are expressed as follows:whereIn (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 disturbances to , , and in (39)–(41), winds do have effects on , , and through which , , and are affected. Also, winds not only have explicit effects on , , and in (36)–(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 to , from to , from to , and from to . Also, the flight conditions include altitude ranging from sea level to ft and Mach number from to and from to . The thrust model of the UAV is assumed aswhere is the thrust command. The thrust is assumed to have a limitation, say .
3.2. 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.(1)In the flight dynamics equations of motion, only , , , , , , , , and are coupled. Therefore , , and .(2)The control input vector and . According to the theory derived in the above, only 4 state variables can be controlled or designated. In this paper, the state vector to be controlled is chosen and the residual vector is .(3)Determine the state variables and control variables for some trim conditions, on which the control design is based.(4)The model matrix in (3) is constructed aswhere the arrangement for the pair () is enlightened from the phugoid mode in which and are coupled. Similar arrangements are for the pairs () and (). In order to make the desired dynamics stable, the elements on the diagonal, , , , , , and , must be positive. Also, those off-diagonal elements must satisfy the conditions, , , and . For convenience of later usage, these necessary conditions are dubbed “the stable model requirements.”
In the design, the initial cruise flight conditions are ft and . The trim conditions are , lb, and . Based on these data, the open-loop eigenvalues of are determined and listed in Table 2.
Longitudinal dynamics analysis reveals that and are closely associated with the short period mode and the phugoid mode, respectively, and lateral dynamics analysis shows that , , and are closely associated with the pure roll mode, the Dutch roll mode, and the spiral mode, respectively. Finally, the rest one, , can be inferred to be associated with the altitude dynamics.
Intuitively, 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 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 two-way approach () is proposed as follows.(1)If the desired eigenvalues are equal to the open-loop eigenvalues as listed in Table 2, obviously, . Choose small which satisfies the stable model requirements.(2)Determine the eigenvalues of the error dynamics in (9).(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 .(6)In each step of simulations, determine and from (4) by using the iteration method described in (5). In this step, only the nominal dynamics model is used.(7)Use 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 , or modify and then go to step .
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 is small and and 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 does improve the performance but makes and deviate their counterparts of open-loop eigenvalues a lot. Also note that the large deviations of , , and from their counterparts of open-loop eigenvalues are due to the iteration procedure in determining a set of , , , and for satisfying the stable model requirements.
3.3. Flight Simulations
In simulations, the initial conditions are , ft, , , , , , , , , and . The states to be tracked are ft, , , and . Although the heading angle is to be tracked, for pilots, it would make more sense to regulate the bank angle rather than the heading angle by assumingwhere is a constant parameter. In this study, .
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 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 from that in the trim conditions, which is usually taken as an initial guess.
For illustration, simulation block diagrams including , , , , and are shown in Figure 2.
After the controls , , , and are determined, (30)–(41) are used for simulations. In order to ascertain if the control system is robust or not, a fictitious wind model is put enroute as follows:where ft, ft, and ft represent the center position of wind zone, and ft, ft, and ft represent the maximum range of wind zone from its center. From (49)–(51), it is also known that the three maximum wind velocity components in the earth fixed frame are knots, knots, and knots, respectively.
In this study, two cases of simulations are conducted; one is without wind disturbances and another with wind disturbances. The results are presented in Figures 3–6.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
As shown in Figure 3(a), there is no much difference in ground trajectories whether wind disturbances are present or not. While flying, the UAV first suffers vertical wind disturbances in between sec and sec and then horizontal wind disturbances in between sec and sec as shown in Figures 3(b) and 3(c).
The state variables to be tracked are all plotted in Figures 4(a)–4(e). From Figure 4(a), it is found that decreases initially when ft is commanded. A careful study reveals that the elevator angle computed is , which is not enough to hold the UAV level as is required. The lowest altitude reached is ft with in Table 4 being used. Increasing will increase the lowest altitude reached but at the expense of increasing overshoot. Fortunately, after descending to ft, the UAV begins to climb. During its climb, the UAV encounters an ascending wind and then a descending wind. But the influences on the altitude are almost negligible. In contrary, the horizontal wind seems to have more influences as it makes the altitude fluctuate. Whether wind disturbances are present or not, the UAV can always approach the commanded altitude asymptotically.
As shown in Figure 4(b), the Mach number does not seem to be so affected by the vertical wind as does by the horizontal wind. In this case, the horizontal wind causes the Mach number to fluctuate up to . The sharp angles shown in the figure are due to the discontinuity of wind acceleration at the edges of wind zone as shown in Figure 3(c).
When the UAV is commanded to make a turn, its heading angle gradually increases and asymptotically approaches the command as shown in Figure 4(c) whether winds are present or not. Figure 4(d) shows that the bank 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, at the instant when and are given. When the UAV reaches the commanded state conditions, the elevator trim angle approaches . 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 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 , the computed command 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 . Winds do have remark effect on both and . Finally, both approach a new value, , 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 and therefore which in turn makes . The minimum pitch rate is deg./s at sec. In this case, the vertical wind does have more remarkable effects on and than the horizontal wind. Finally, both approach and , 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.
4. Conclusions
In this paper, a unified nonlinear dynamic inversion control system is successfully developed. With this method, the 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors would like to express their gratitude to Chuan-Tau Edward Lan for his helps in computing the aerodynamic data of the UAV used as an example for the control design demonstrated in this paper.