System and method for robust nonlinear regulation control of unmanned aerial vehicles synthetic jet actuators

Abstract

An unmanned aerial vehicle (UAV) is provided with a plurality of synthetic jet actuators and a nonlinear robust controller. The controller compensates for uncertainty in a mathematic model that describes the function of the synthetic jet actuators. Compensation is provided by the use of constant feedforward best guess estimates that eliminate the need for more highly computationally burdensome approaches such as the use of time-varying adaptive parameter estimation algorithms.

Description

FIELD OF THE INVENTION

The present invention generally relates to a flight control system for selective operation of an array of synthetic jet micro-actuators employed on an unmanned aerial vehicle.

BACKGROUND

There has been a surge of interest in the design and application of unmanned aerial vehicles (UAVs), particularly fixed wing micro-air vehicles (MAV). UAVs can be used in numerous civilian applications, such as urban reconnaissance, package delivery, and area mapping. UAVs are utilized in various military applications as well. While some UAVs require human operators remotely controlling the operation and/or flight path of the UAV, other UAV's have been designed for autonomous operation. One of the biggest challenges involved in the autonomous operation of UAVs is in the design of flight tracking controllers for UAVs, which often operate in uncertain and possibly adverse conditions.

Flight tracking controllers should provide the UAV, operating autonomously, with the ability to remain within an acceptable range of their intended flight path while operating in an uncertain environment caused by, for example, changing wind velocity and direction resulting in gusts, vortical structures and turbulence.

Additionally, suppression of limit cycle oscillations (LCO) (or flutter) is another concern in UAV tracking control design. This is especially true for applications involving smaller, lightweight UAVs, where the aircraft wings are more susceptible to LCO. LCO refers to “flutter” behaviors in UAV wings that manifest themselves as constant-amplitude oscillations, which result from nonlinearities inherent in the aeroelastic dynamics of the UAV system. Note that LCOs may significantly affect the aerodynamic properties of an aircraft, and can be especially devastating for small UAVs interacting with impinging gusts with amplitudes comparable to the vehicle speed. Indeed, the gust-induced wing dynamics results in lift, drag, and pitching-moment oscillations deteriorating the aircraft aerodynamic performance and posing severe challenges to aircraft flight stability. Due to these behaviors, the LCO could surpass the safe flight boundaries of an aircraft and could potentially lead to structural damage causing the UAV to crash.

These engineering challenges necessitate the utilization of UAV flight controllers, which achieve accurate flight tracking in the presence of dynamic uncertainty while simultaneously suppressing LCO. Control applications for LCO suppression have been developed using mechanical deflection surfaces (e.g., flaps, ailerons, rudders, and elevators). However, when dealing with smaller UAVs, practical engineering considerations and physical constraints can preclude the addition of the large, heavy moving parts that generally are required for installation of deflection surfaces.

Due to their small size, ease of operation, and low cost, synthetic jet actuators (SJA) are promising tools for aircraft tracking control and flow control applications. SJA's transfer linear momentum to a flow system by using a vibrating diaphragm, which creates trains of vortices through the alternating ejection and suction of fluid through a small orifice. Since these vortices (i.e., jets) are formed entirely from the fluid (i.e., air) of the flow system, a key benefit of SJA is that they achieve this transfer of momentum with zero net mass injection across the flow boundary. Thus, SJAs do not require space for a fuel supply. SJAs can be utilized to modify the boundary layer flow field near the surface of an aircraft wing, which can improve aerodynamic performance. Moreover, synthetic jet actuators can expand the usable range of angle of attack, which can improve maneuverability. In addition to flow control applications, arrays consisting of several SJAs can be employed to achieve tracking control of aircraft, possibly eliminating the need for mechanical control surfaces. The benefits of utilizing SJAs on aircraft as opposed to mechanical control surfaces include reduced cost, weight, and mechanical complexity.

SJAs have been developed with the capability to achieve momentum transfer with zero-net mass-flux. This beneficial feature eliminates the need for an external fuel supply, since the working substance is simply the gas (i.e., air) that is already present in the environment of operation. This makes SJAs an attractive option in UAV applications, because of the significant reduction in the size of the required equipment. The SJAs synthesize the jet flow through the alternating suction and ejection of fluid through an aperture, which is produced via pressure oscillations in a cavity. The pressure oscillations can be generated using various methods, including pistons in the SJA's orifices or piezoelectric diaphragms. SJA's can achieve boundary-layer flow control near the surface of a UAV wing since they can provide instant actuation, unlike conventional mechanical control surfaces. In addition, SJA's can expand the usable range of angle of attack, resulting in improved UAV maneuverability.

FIGS. 1 and 2 show an example UAV 100 with a seamless fixed wing construction. The UAV 100 includes an airframe 102 and an array of synthetic jet actuators 104 disposed at selected locations around the periphery of the wing of the airframe 102. FIG. 3 schematically illustrates a synthetic jet actuator 104. The SJA includes a cavity 106, a membrane 108, also called a diaphragm, and an outlet orifice 110. The SJA 104 is provided with a voltage to cause the membrane 108 to oscillate, resulting in a time averaged pulsed jet of air being emitted from the orifice 110. The magnitude of the voltage applied to the SJA 104 may impact the flex or force of air expelled from the orifice 110. Under certain conditions, a series of vortex rings 112 formed during the expulsion phase are able to escape the reverse injection and convect away from the orifice as a time-averaged jet capable of manipulating the flight path of the airframe 102 and/or the flow of air affecting the airframe 102.

Use of SJA present challenges in control design due to uncertainties inherent in the dynamics of their operation. Specifically, the input-output (i.e., the control voltage to force delivered) characteristic of each SJA is nonlinear and contains parametric uncertainty. To compensate for this uncertainty, recently developed SJA-based control systems utilize online adaptive control algorithms, neural networks, and/or complex fluid dynamics computations in the feedback loop. FIG. 4 schematically illustrates one control system 200 that the present disclosure seeks to improve. Particularly, the control system 200 includes a processor 202 configured to generate a control command to a set of UAV system dynamics 204, such as an array of synthetic jet actuators. The processor 202 determines the controller command based on part on sensor measurements 206 which form a first feedback loop 208 to the processor 202. The instant control system 200 is characterized by the use of an adaptive parameter update law 210, effectively requiring a second feedback loop 212 such that the processor 202 generates the controller command based on a control law which is required to adapt or update itself with each iteration. However, implementation of these computationally burdensome control techniques, i.e., the adaptive parameter update law 210 and second feedback loop 212 requires additional computational resources (e.g., microprocessors, timer circuits, and/or time-consuming processing algorithms), which might not be available for smaller UAV applications with limited onboard space and computing power. Moreover, the computational complexity in these control system designs results in a slower control response rate, which can hinder control performance or even cause catastrophic failures.

SUMMARY

In one embodiment, the present disclosure provides a robust flight controller to achieve accurate trajectory tracking for an unmanned aerial vehicle equipped with an array of synthetic jet actuators over a wide envelope of actuator uncertainty and aircraft operating conditions. This disclosure suggests the use of a robust, continuous Lyapunov-based controller that includes a novel implicit learning characteristic.

The present disclosure includes an improved control system configured to achieve UAV tracking control performance that is comparable to adaptive or neural network-based methods with a reduction in the required computational resources (i.e. less expensive) and a faster response time. The control system may require only a single sensor feedback measurement loop. By eliminating an adaptive parameter feedback generally found in comparable systems, the control system of the present disclosure is more suitable for applications involving small UAVs with limited onboard space and computing power. Further, the simplicity in the disclosed control design endows the system with a faster control response rate, which results in more reliable UAV flight in unpredictable, time-varying flight conditions.

The present disclosure includes an unmanned aerial vehicle (UAV), comprising: an airframe, a plurality of synthetic jet actuators (SJA) affixed to the airframe, a plurality of sensors configured to determine an operating condition of the UAV, and a control system. Each SJA is configured to selectively produce a stream of air in response to a control command input thereto, each SJA being represented by a mathematical model having at least two uncertain parameters. The control system is configured to provide the control command to each of the plurality of synthetic jet actuators based on the operating condition of the UAV and a control law, the control law including a constant estimate for each of the uncertain parameters in the mathematical model of each corresponding SJA. The control system controls at least one or both of the trajectory of the UAV and the vibration of the UAV by activating at least one of the plurality of SJA with the control command.

The present disclosure also includes a control system for determining a control command usable by a synthetic jet actuator to at least one of, adjust a trajectory of an unmanned aerial vehicle, or suppress limit cycle oscillations of the unmanned aerial vehicle. The control system comprises a processor programmed to run a control law, the control law comprising a feedback term determined based on measurements from at least one sensor, and a constant estimation of at least two uncertain parameters necessitated by use of the synthetic jet actuator, which is represented by a mathematical model having at least two uncertain parameters. The control system is designed such that the control law is substantially free from time-varying adaptive parameter estimation algorithms.

The present disclosure also includes a method of controlling a micro air vehicle. The method comprises measuring at least one of roll rate, pitch rate, yaw rate, pitching rate, and plunging rate. The method may comprise comparing the measured value with a desired reference value to determine an error. The method may further comprise determining a control command voltage based on a control law having constant estimates and an error variable term. The method may include applying the control command voltage to a synthetic jet actuator of a plurality of synthetic jet actuators, and emitting a flow of air from the synthetic jet actuator to adjust the trajectory of the micro air vehicle.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a perspective schematic of an example unmanned aerial vehicle (UAV) equipped with an array of synthetic jet actuators (SJAs).

FIG. 2 is a side view schematic of an example arrangement of SJAs on the trailing edge of a wing of the UAV of FIG. 1.

FIG. 3 is a schematic example of a SJA.

FIG. 4 is a diagrammatic representation of a conventional control system for the UAV of FIG. 1.

FIG. 5 is a diagrammatic representation of a control system for the UAV of FIG. 1 according to embodiments of the present disclosure.

FIG. 6 shows plots of tracking error achieved during closed-loop controller operation according to one simulation.

FIG. 7 shows plots of control signals v_i (t), i=1, . . . 6 used during closed-loop controller operation according to the simulation of FIG. 6.

FIG. 8 plots uncontrolled pitching response for a free stream velocity U=18.25 m/s according to one example.

FIG. 9 plots uncontrolled plunging response for a free stream velocity U=18.25 m/s according to one example.

FIG. 10 plots uncontrolled pitching response for a free stream velocity U=19.5 m/s.

FIG. 11 plots uncontrolled plunging response for a free stream velocity U=19.5 m/s.

FIG. 12 plots the closed-loop pitching response for a free stream velocity U=18.25 m/s.

FIG. 13 plots the closed-loop plunging response for a free stream velocity U=18.25 m/s.

FIG. 14 plots actuator force required for free stream velocity U=18.25 m/s.

FIG. 15 plots actuator moment required for free stream velocity U=18.25 m/s.

FIG. 16 plots the closed-loop pitching response for a free stream velocity U=19.5 m/s.

FIG. 17 plots the closed-loop plunging response for a free stream velocity U=19.5 m/s.

FIG. 18 plots the actuator force required for free stream velocity U=19.5 m/s.

FIG. 19 plots the actuator moment required for free stream velocity U=19.5 m/s.

FIG. 20 plots the closed-loop pitching response for free stream velocity of U=20.5 m/s.

FIG. 21 plots the closed-loop plunging response for free stream velocity of U=20.5 m/s.

FIG. 22 plots actuator force required for free stream velocity of U=20.5 m/s.

FIG. 23 plots the actuator moment required for free stream velocity of U=20.5 m/s.

FIG. 24 plots a closed-loop pitching response in the presence of a disturbance.

FIG. 25 plots the closed-loop plunging response in the presence of a disturbance.

FIG. 26 plots actuator force used in the presence of the disturbance of FIGS. 24 and 24.

FIG. 27 plots actuator moment in response to the disturbance of FIGS. 24 and 25.

FIG. 28 illustrates a sharp-edge gust-airfoil interaction model.

FIG. 29 illustrates example details of a numerical grid.

FIG. 30a illustrates flat plate pitching LCO amplitudes.

FIG. 30b illustrates time histories of airfoil pitching responses.

FIGS. 31a and 31b illustrate example vorticity and streamline contours with suction-side actuation.

FIGS. 32a-c illustrate a transition to uncontrolled pitching LCO for increasing amplitudes of initial excitation.

FIGS. 33a-c illustrate an example of suppression of pitching LCO achieved by the feedback-loop robust control system with increasing amplitudes of initial excitation.

FIGS. 34a-f illustrate an actuator's control authority requirements for aerodynamic lift and moment with increasing amplitudes of initial excitation.

DISCUSSION OF INVENTION

The present invention is directed to a nonlinear robust controller method of control for operating synthetic jet actuators (SJA) used as part of unmanned aerial vehicles (UAV). In one embodiment, the SJA are configured for operation to assist with tracking of the UAV. In another embodiment, the use of the SJA may be configured to substantially suppress limit cycle oscillations (LCO) in unmanned aerial vehicle (UAV) systems with uncertain dynamics. In each embodiment, the controller method seeks to achieve accurate flight control in the presence of SJA non-linearities, parametric uncertainty and external disturbances (e.g., wind gusts). Particularly, the control method utilizes a control law that is continuous, making the method amenable to practical applications of small UAV, sometime referred to as micro air vehicles (MAV). Moreover, the control method presented herein is designed to be inexpensively implemented, requiring no online adaptive laws, function approximators, or complex fluid dynamics computations in the feedback loop. A matrix decomposition technique is utilized along with innovative manipulation in the error system development to compensate for the dynamic SJA uncertainty. The robust controller is designed with an implicit learning law, which is shown to compensate for bounded disturbances. A Lyapunov-based stability analysis is utilized to prove global asymptotic trajectory tracking in the presence of external disturbances, actuator nonlinearities, and parametric uncertainty in the system and actuator dynamics. Numerical simulation results are provided to complement the theoretical development. A salient feature of the robust control laws disclosed by this description is that their structure is continuous, and enables bounded disturbances to be asymptotically rejected without the need for infinite bandwidth. In one embodiment, a rigorous Lyapunov-based stability analysis is utilized to prove asymptotic pitching and plunging regulation, considering a detailed dynamic model of the pitching and plunging dynamics.

FIG. 5 schematically illustrates an aircraft control system 300 according to an embodiment of the present disclosure. The aircraft control system 300 may be configured for tracking along a predetermined flight path. In another embodiment, the aircraft control system 300 may be used to prevent limit cycle oscillations. The control system 300 is intended for use with a micro air vehicle (e.g. an unmanned aerial vehicle) having an airframe (see FIG. 1) and an actuation system 302, for example an array of synthetic jet actuators (see FIGS. 1-3). The airframe may include a plurality of sensors 304 configured to measure parameters of the airframe that impact the operation of the UAV. In one example, the sensors 304 may determine at least one, and potentially all three of the roll rate, pitch rate, and yaw rate of the airframe to be used to control tracking along a predetermined desired flight path. The sensors 304 may function based on displacement or acceleration in alternative to rates. Where LCO is being controlled, the motion of the airframe to be sensed may be the pitching and/or plunging action of the wings. The control system 300 also includes a controller 306 configured to provide a controller command, such as a peak-to-peak voltage, to selectively excite operation of the actuation system 302. The controller 306 may include a processor or a microprocessor configured to execute a control implementation and analysis law for determining the appropriate controller command. The control law may be configured to require input from the sensors 304 as the primary dynamic parameters. The sensed parameters provided by the sensors 304 may be compared to a predetermined desired reference condition. When attempting to control LCO, the desired reference condition could be a steady state, i.e. zero pitch or plunge rate. When tracking, the predetermined desired reference conditions may vary as the UAV is expected to follow a predetermined set of maneuvers. The remainder of the inputs used to execute the control law are intended to be constant, best-guess estimates of uncertain parameters. In other words, the control law does not function based upon the use of any time-varying adaptive parameter estimation algorithms.

Reiterating from above, use of SJA present challenges in control design due to uncertainties inherent in the dynamics of their operation. Specifically, the input-output (i.e., the control voltage to force delivered) characteristic of each SJA is nonlinear and contains parametric uncertainty. The uncertain aircraft dynamic model detailed herein contains parametric uncertainty due to linearization errors and un-modeled nonlinearities. Specifically, the aircraft system can be modeled via a linear time-invariant system as:{dot over (x)}=Ax+B+ƒ(x,t)  (1)where Aε^n×n denotes the uncertain state matrix, Bε^n×n represents the uncertain input matrix, and ƒ(x,t)εⁿ is a state- and time-dependent unknown, nonlinear disturbance. For example, ƒ(x,t) could include exogenous disturbances (e.g., due to wind gusts) or nonlinearities not captured in the linearized dynamic model. The state vector x(t) may contain roll, pitch, and yaw rate measurements. Also in (1), the control input

$u\left(t\right)\stackrel{\Delta }{=}{\left[u_1\left(t\right)u_2\left(t\right)\dots u_m\left(t\right)\right]}^{𝒯}\in {ℝ}^m$ represents the virtual surface deflections due to m arrays of synthetic jet, actuators (SJA).

Based on experimental data, the dynamics of the SJA can be modeled as

$\begin{array}{cc}{u}_i={\theta }_{2i}^{*}-\frac{{\theta }_{1i}^{*}}{v_i}i=1,2,\dots ,m& \left(2\right)\end{array}$ where _i()=()ε denotes the peak-to-peak voltage acting, on the i^th SJA array, and θ_1i*, θ_2i*ε are unknown positive physical parameters. u is the virtual control surface deflection angle expressed in degrees. In a standard UAV, this would represent the angles of the elevator, aileron, and rudder. In this SJA application, it is the equivalent deflection angle that is generated by an SJA array. v_i(t) is the peak-to-peak control input voltage to the SJA array expressed in volts for each SJA i=1−n. θ_1i*, θ_2i* are uncertain constant parameters, which are inherent in the SJA actuator dynamics. Nominal values for these parameters have been obtained in previous research by Deb et al as θ_1i*=33.33 volt-deg and θ_2i*=15 degrees. These parameter values can fluctuate significantly depending on the UAV operating conditions, and the resulting parametric uncertainty is the primary challenge that must be addressed in SJA-based UAV control design.

Based on the uncertain SJA actuator model given in (2), the control input voltage is designed via the feedback control law

$\begin{array}{cc}v\left(t\right)=\frac{{\hat{\theta }}_1}{{\hat{\theta }}_2-u_d\left(t\right)}& \left(3\right)\end{array}$

In Equation (3), {circumflex over (θ)}₁, and θ₂* denote constant, best-guess estimates of the uncertain parameters θ₁* and θ₂*, respectively. The use of constant estimates as opposed to time-varying adaptive estimates facilitates the improvements of the control system 300 of the present invention over the prior control system 200. Use of best guess estimates of the uncertain parameters enables the tracking control system 300 to eliminate the additional sensor feedback measurements and processing effort (i.e., the adaptive parameter update law) required by control system 200. {circumflex over (θ)}₁ depends on physical and aerodynamic parameters, including wing chord, the freestream velocity, and additional physical parameters that lead to one level of uncertainty. An empirically determined expression for θ₁* is

${\theta }_1^{*}=\frac{p_2{p_3\left(U_{\infty }\right)}^3}{{\mathrm{fcp}}_4},$ where U_∞ denotes the freestream velocity, ƒ is the frequency of the input voltage v(t), c is the local wing chord, and p₂, p₃, and p₄ denote uncertain constant physical parameters. In the original experimental research of Deb et al, the values of the parameters p₂, p₃, and p₄ were selected arbitrarily as a baseline model resulting in θ₁* values consistent with those used in present simulations as found in Table 1, {circumflex over (θ)}₂ may be estimated as the maximum deflection angle that can be achieved by an SJA. Presently, this is suggested to be about 15 degrees.

Also in Equation (3), u_d(t) denotes an auxiliary control term, which incorporates the sensor feedback measurements (e.g., the roll, pitch, and yaw rates defining the UAV flight trajectory as shown in FIG. 5). After substituting Equation (3) into Equation (2) and considering Equation (1), it can be shown that the auxiliary control term u_d(t) is pre-multiplied by an uncertain term Ω (a matrix in general) containing the uncertain parameter(s) θ₁*, the constant estimate(s) {circumflex over (θ)}₁, and the elements of the uncertain input matrix B. It can further be shown that the elements of the uncertain matrix Ω contain the ratios θ₁*/{circumflex over (θ)}₁ of the uncertain parameters to their estimates. To compensate for the uncertainty in Ω, an estimate {circumflex over (Ω)} of Ω can be designed to be equivalent to the best guess estimate of the uncertain input matrix B in Equation (1).

By incorporating the estimate {circumflex over (Ω)}, the control term u_d(t) is designed asu_d(t)={circumflex over (Ω)}^#(μ₀(t)−μ₁(t))  (4)where {circumflex over (Ω)}^# denotes the pseudoinverse of {circumflex over (Ω)}, μ₀(t) and μ₁(t) are functions of the UAV roll, pitch, and yaw rate tracking errors defined asμ₀=−(k_s+I_n×n)e(t)−(k_s+I_n×n)e(0)−∫₀^tα(k_s+I_n×n)e(τ)dτμ₁=∫₀^tβsgn(e(τ))dτ.  (5)

The variables in Equation (5) are defined as follows:

• • e(t) is the UAV trajectory tracking error expressed in degrees per second. Physically, this represents the difference between the actual sensor measurements from sensors 304 and the reference (desired) UAV trajectory, for example based on a predetermined flight plan. The actual sensor measurements are the roll, pitch, and yaw rate measurements as shown in FIG. 5. The UAV tracking control objective can be stated as e(t)→0.
  • {circumflex over (Ω)} is a constant estimate of the parametric uncertainty resulting from the mismatch between the actual and estimated values of the uncertain SJA parameters.
  • k_s, α, β are positive, constant, user-defined control gains (i.e., amplifier settings). The values of these parameters are utilized to finely ‘tune’ the performance of the closed-loop system in terms of control response rate and trajectory tracking accuracy, for example.

Thus by configuring the controller 306 to operate in accordance with the control law consistent with equations (3-5) above, and particularly consistent with equations (3) and (4) above, the control command (e.g. v(t)) may be determined based on A) constant values combined with B) the measured rates of motion (such as pitch, roll and yaw for flight path control) available from sensors 304 fed back to the controller 306 during operation of the tracking control system 300 and compared to C) a predetermined motion state. The predetermined motion state may be a predetermined expect flight path in the example of tracking control, or the predetermined motion state may be steady state in the example of limiting oscillations.

Example 1: Tracking Control Embodiments

Using the expressions in (2) and (1), the dynamics can be expressed in terms of the i^th SJA array as{dot over (x)}=Ax+Σ_t=1^m_ii+ƒ(x,t).  (6)

In (6),

$b_i\stackrel{\Delta }{=}{\left[b_{1i}{b}_{2i}\dots b_{\mathrm{ni}}\right]}^{𝒯}\in {ℝ}^n\forall i=1,2,\dots ,m,$ where b_ij by denotes the (i, j)^th element of the matrix B in (1).

Assumption 1: If x(t)ε_∞, then ƒ(x,t) is bounded. Moreover, if x(t)ε_∞, then the first and second partial derivatives of the elements of ƒ(x,t) with respect to x(t) exist and are bounded.

A purely robust feedback control strategy can be utilized to compensate for the control input nonlinearity and input parametric uncertainty in (2). To this end, a robust inverse v_i(t) is utilized, which contains constant feedforward “best-guess” estimates of the uncertain parameters θ_1i* and θ_2i*. The robust inverse that compensates for the uncertain, jet array nonlinearities in (2) can be expressed as

$\begin{array}{cc}{v}_i\left(t\right)=\frac{{\hat{\theta }}_{1i}}{{\hat{\theta }}_{2i}-u_{2i}\left(t\right)},i=1,\dots ,m& \left(7\right)\end{array}$ where {circumflex over (θ)}_1i, {circumflex over (θ)}_2iε⁺ are constant feedforward estimates of θ_1i* and θ_2i*, respectively, and u_di(t)ε∀i=1, . . . , m are subsequently defined auxiliary control signals.

Remark 1. Singularity Issues Based on (7), the control signal v_i(t) will encounter singularities when u_di(t)={circumflex over (θ)}_2i. To ensure that the control law in (7) is singularity-free, the control signals u_di(t) for i=1, 2, . . . , m are designed using the following algorithm:

$\begin{array}{cc}{u}_{\mathrm{di}}\left(t\right)=\{\begin{array}{cc}{\hat{\theta }}_{2i}-\delta & \mathrm{if}g\left({\mu }_0\left(t\right),{\mu }_1\left(t\right)\right)\ge {\hat{\theta }}_{2i}-\delta \\ g\left({\mu }_0\left(t\right),{\mu }_1\left(t\right)\right)& \mathrm{otherwise}\end{array}& \left(8\right)\end{array}$ where δε⁺ is a small parameter, g(.) is a subsequently defined function, and μ₀(t), μ₁(t)ε^m are subsequently defined feedback control terms. Note that the parameter δ can be selected arbitrarily small such that the subsequent stability analysis remains valid for an arbitrarily large range of positive control voltage signals v_i(t).

In addition, the control terms u_i(t) in (2) will encounter singularities when v_i(t)=0, which occurs when {circumflex over (θ)}_1i=0 for any i. Since {circumflex over (θ)}_1i is a constant, user-defined feedforward estimate of the uncertain parameter θ_1i*, the singularity at v_i(t)=0 can be easily avoided by selecting {circumflex over (θ)}_1i>0 for i=1, 2, . . . , m,

Remark 2. The auxiliary control signal u_di(t) in (7) can be designed to achieve asymptotic tracking control and disturbance rejection for the uncertain dynamic model in (1) and (2) over a wide range of feedforward estimates {circumflex over (θ)}_ji≠θ_ji*, j=1, 2.

The control objective is to force the system state x(t) to track the state of a model reference system. Based on (1), reference model is selected as:{dot over (x)}_m=A_mx_m+B_mδ  (9)where x_m(t)εⁿ denotes the reference state, A_mε^n×n is a Hurwitz state matrix, B_mεⁿ denotes the reference input matrix, and δ(t)ε is the reference input signal. The reference model in (9) is designed to exhibit desirable performance characteristics.

Assumption 2: The model reference state x_m(t) is bounded and sufficiently smooth in the sense that x_m(t), {dot over (x)}_m(t), {umlaut over (x)}_m(t), (t)ε_∞∀t≧0.

To quantify the control objective, the e(t)εⁿ is defined ase=x−x_m  (10)

To facilitate the subsequent analysis, a filtered tracking error r(t) is also defined asr=ė+αe  (11)

After taking the time derivative of (11) and using (1) and (10), the open loop tracking error dynamics are obtained as

$\begin{array}{cc}\stackrel{.}{r}=A\stackrel{.}{e}+A{\stackrel{.}{x}}_m+{\sum }_{i=1}^m{b}_i\left(\frac{{\Theta }_{1i}^{*}}{{\Theta }_{2i}^{*}}{\stackrel{.}{u}}_{\mathrm{di}}\left(t\right)\right)+\stackrel{.}{f}\left(x,t\right)-{\ddot{x}}_m+\alpha \stackrel{.}{e}& \left(12\right)\end{array}$

Remark 3: Although the instant portion of the control input term vanishes upon taking the time derivative of the dynamics as in (12), the plant model: used in the subsequent numerical simulation retains the complete actuator dynamics. In the simulation, the control input u_i(t) is generated using (2) and (7); thus, the simulation model includes the complete actuator dynamics.

The expression in (12) can be rewritten as{dot over (r)}=Ñ+N_d+Ω{dot over (μ)}_d(t)−Se  (13)where Ωε^n×m denotes a constant uncertain matrix, Sε^n×n is a subsequently defined uncertain matrix and the control vector

$u_d\left(t\right)\stackrel{\Delta }{=}{\left[u_{d1}\left(t\right),u_{d2}\left(t\right)\dots u_{dm}\left(t\right)\right]}^{\tau }\in {ℝ}^m.$

In (13), the unknown, immeasurable auxiliary functions Ñ(t) and N_d(t) are defined as

$\begin{array}{cc}\tilde{N}\stackrel{\Delta }{=}A\stackrel{.}{e}+\alpha \stackrel{.}{e}+𝒮e+\left(\stackrel{.}{f}\left(x,t\right)-\stackrel{.}{f}\left(x_m,t\right)\right)& \left(14\right)\\ N_d\stackrel{\Delta }{=}A{\stackrel{.}{x}}_m-{\ddot{x}}_m+f\left(x_m,t\right)& \left(15\right)\end{array}$

The selective grouping of terms in (14) and (15) is motivated by the fact that Assumptions 1 and 2 can be utilized to develop the following inequalities:∥Ñ∥≦ρ₀∥z∥,∥N_d∥≦ζN_d,∥{dot over (N)}_d∥≦ζN_d  (16)where ρ₀, ζ_N1, ζ_Ndε⁺ are known bounding constants and z(t)ε²ⁿ is defined as

$\begin{array}{cc}z\stackrel{\Delta }{=}[\begin{array}{cc}{e}^T& {r^T]}^T\end{array}& \left(17\right)\end{array}$

Based on the open-loop error system in (13), the auxiliary control u_d, (t) is designed asu_d(t)={circumflex over (Ω)}^#(μ₀−μ₁)  (18)where {circumflex over (Ω)}ε^n×m is a constant, best-guess estimate of the uncertain matrix Ω, and [.]^# denotes the pseudoinverse of a matrix. In (18), μ₀(t), μ₁(t)εⁿ are subsequently defined feedback control terms. After substituting the time derivative of (18) into (13), the error dynamics can be expressed as{dot over (r)}=Ñ+N_d+{tilde over (Ω)}({dot over (μ)}₀−{dot over (μ)}₁)−Se  (19)where the constant uncertain matrix {tilde over (Ω)}ε^n×n is defined as{tilde over (Ω)}=Ω{circumflex over (Ω)}^#.  (20)

Assumption 3: Bounds, on the uncertain matrix Ω are known such that the feedforward estimate Ω can be selected, such that the product {tilde over (Ω)} can be decomposed as{tilde over (Ω)}=ST  (21)where Sε^n×n is a positive definite symmetric matrix; and Tε^n×n is a unity upper triangular matrix, which is diagonally dominant in the sense thatε≦[T_ii]−Σ_k=i+1ⁿ[T_ik]≦Q,i=1, . . . ,n−1.  (22)

In inequalities (22), εε(0,1) and Qε⁺ are known bounding constants, and T_ikε denotes the (i, k)^th element of the matrix T.

Remark 4: Assumption 3 is mild in the sense that (22) is satisfied over a wide range of {circumflex over (Ω)}≠Ω. Specifically, the auxiliary control signal u_di(t) in (7) and (18) can be designed to achieve asymptotic tracking control and disturbance rejection for the uncertain dynamic model in (1) and (2) when the mean values of the constant feedforward estimates {circumflex over (Θ)}_j1 and {circumflex over (Θ)}_j2∀j=1, . . . , 6 differ from the mean values of the actual parameters θ_j1*, and θ_j2*, ∀j=1, . . . , 6 as

• • mean {circumflex over (Θ)}_j1=14.05, mean θ_j1*, =25.28
  • mean {circumflex over (Θ)}_j2=7.25, mean θ_j2*=12.08

The values for {circumflex over (Θ)}_j1 and {circumflex over (Θ)}_j2∀j=1, . . . , 6 used in the simulation can be found in Table 1. This result demonstrates the capability of, the robust control design to compensate for significant dynamic uncertainty using, only a simple feedback controller structure.

TABLE 1
Array i	1	2	3	4	5	6
Θ1i* [volt-deg]	32.9	29.8	26.7	24.0	20.5	17.8
Θ2i* [deg]	14.7	13.8	12.8	11.7	10.0	9.5
{circumflex over (Θ)}1i* [volt-deg]	16.5	15.9	14.5	13.4	12.1	11.9
{circumflex over (Θ)}2i* [deg]	9.1	8.3	7.2	6.8	6.5	5.6

After utilizing the decomposition in (21), the error dynamics in (19) can be expressed as

$\begin{array}{cc}{𝒮}^{-1}\stackrel{.}{r}={\tilde{N}}_1+N_{d1}+\stackrel{\_}{T}\left({\stackrel{.}{\mu }}_0-{\stackrel{.}{\mu }}_1\right)-e\mathrm{where}& \left(23\right)\\ {\tilde{N}}_1\stackrel{\Delta }{=}{𝒮}^{-1}\tilde{N},N_{d1}\stackrel{\Delta }{=}{𝒮}^{-1}{N}_d.& \left(24\right)\end{array}$

Since S is positive definite, Ñ₁(t) and N_d1(t) satisfy the inequalities∥Ñ₁∥≦ρ₁∥z∥,∥N_d1∥≦ζ_Nd1∥{dot over (N)}_d1∥≦ζ_Nd1  (25)where ρ₁, ζ_Nd1, ζ_Nd1ε⁺ are known bounding constants. By utilizing the fact that the uncertain matrix T is unity upper triangular, the error dynamics in (23) can be rewritten:S⁻¹{dot over (r)}=Ñ₁+N_d1+{dot over (μ)}₀+T{dot over (μ)}_o−T{dot over (μ)}₁−e  (26)where

$\stackrel{\_}{T}\stackrel{\Delta }{=}T-I_{\mathrm{nxn}}$ is a strictly upper triangular matrix, and I_n×n denotes an n×n identity matrix. Based on the open-loop error system in (26), the auxiliary control terms μ₀(t) and μ₁(t) are designed as

$\begin{array}{cc}{\mu }_0=-\left(k_s+I_{\mathrm{nxn}}\right)e\left(t\right)-\left(k_s+I_{\mathrm{nxn}}\right)e\left(0\right)-{\int }_0^t\alpha \left(k_s+I_{\mathrm{nxn}}\right)e\left(r\right)\mathrm{dr}& \left(27\right)\\ {\mu }_1={\int }_0^t\mathrm{\beta sgn}\left(e\left(r\right)\right)\mathrm{dr}& \left(28\right)\end{array}$ where βε is a constant, positive control gain, k_sε^n×n is a constant, positive definite, diagonal control gain matrix, and α is introduced in (11). After substituting the time derivative of (27) into (26), the closed-loop error system is obtained asS⁻¹{dot over (r)}=Ñ₁+T{dot over (μ)}_o+N_d1−(k_s+I_n×n)r−T{dot over (μ)}₁−e.  (29)

After taking the time derivative of (27), the term T{dot over (μ)}_o can be expressed as

$\begin{array}{cc}\stackrel{\_}{T}{\stackrel{.}{\mu }}_0=\left[\begin{array}{c}{\sum }_{j=2}^n{\stackrel{\_}{T}}_{1j\stackrel{.}{\mu }0j}\\ {\sum }_{j=3}^n{\stackrel{\_}{T}}_{2j\stackrel{.}{\mu }0j}\\ ⋮\\ {\stackrel{\_}{T}}_{\left(n-1\right)n\stackrel{.}{\mu }\stackrel{.}{o}n}\\ 0\end{array}\right]=\left[\begin{array}{c}{A}_p\\ 0\end{array}\right]& \left(30\right)\end{array}$ where the auxiliary signal

$A_p\stackrel{\Delta }{=}{\left[A_{p1}{A}_{p2}\dots A_{p\left(n-1\right)}\right]}^T\in {ℝ}^{n-1},$ with the individual elements defined as

$\begin{array}{cc}{A}_{\mathrm{pi}}\stackrel{\Delta }{=}-{\sum }_{j=t+1}^n{\stackrel{\_}{T}}_{\mathrm{ij}}\left(k_{\mathrm{sj}}+1\right)r_j& \left(31\right)\end{array}$

∀i1, . . . , n−1 where the subscript j indicates the jth element of the vector. Based on the definitions in (27) and (30), A_p satisfies, the inequality∥A_p∥≦ρ_Λ1∥z∥  (32)

Remark 5: Note that based on (30) and (31), the bounding constant ρ_Λ1 depends only on elements i+1 to n of the control gain matrix k_s due to the strictly upper triangular nature of T. Thus, the element {dot over (μ)}₀₁(t) of the control vector {dot over (μ)}₀₁(t) does not appear in the term A_p. This fact will be utilized in the subsequent stability proof.

By utilizing (30), the error dynamics in (29) can be expressed as

$\begin{array}{cc}{S}^{-1}\stackrel{.}{r}={\tilde{N}}_2+N_{d1}-\left(k_s+I_{\mathrm{nxn}}\right)\tau -T{\stackrel{.}{\mu }}_1-e\mathrm{where}& \left(33\right)\\ {\stackrel{\_}{N}}_2={\tilde{N}}_1=\left[\begin{array}{c}{A}_p\\ 0\end{array}\right]& \left(34\right)\end{array}$

Based on (25), (32), and (34), N₂(t) satisfies the inequality∥N₂∥≦ρ₂∥z∥  (35)

where ρ₂ε is a known bounding constant.

To facilitate the subsequent stability analysis, the control gain β introduced in (28) is selected to satisfy

$\begin{array}{cc}\beta >\frac{1}{\varepsilon }\left({\zeta }_{N_{d1}}+\frac{1}{\alpha }{\zeta }_{{\stackrel{.}{N}}_{d1}}\right)& \left(36\right)\end{array}$ where ζ_Nd1 and ζ_Nd1 are introduced in (25), and ε is introduced in (22).

To facilitate the subsequent stability analysis, let D⊂²ⁿ⁺¹ be a domain containing (t)=0, where (t)ε²ⁿ⁺¹ is defined as

$\begin{array}{cc}\omega \left(t\right)\stackrel{\Delta }{=}{\left[z^T\left(t\right)\sqrt{P\left(t\right)}\right]}^T.& \left(37\right)\end{array}$

In (37), the auxiliary function P(t)ε is defined as the generalized solution to the differential equation{dot over (P)}(t)=−L(t), P(0)=βQ|e(0)−e^T(0)N_d1^T(0)  (38)where the auxiliary function L(t)ε is defined asL(t)=r^T(N_d1(t)−T{dot over (μ)}₁).  (39)

Lemma 1: Provided the sufficient gain condition in (36) is satisfied, the following inequality can be obtained:∫_o^tL(τ)dr≦βQ|e(0)|−e^T(o)N_d1(0)  (40)

Hence, (40) can be used to conclude that P(t)≧0.

Theorem 1: The robust control law given by (7), (18), (27), and (28) ensures asymptotic trajectory tracking in the sense that∥e(t)∥→0, as t→∞  (41)provided the control gain matrix k_s, introduced in (27) is selected sufficiently large (see the subsequent proof), and β is selected to satisfy the sufficient condition in (36).

Proof: Let V(w,t): Dx[0,∞)→ be a continuously differentiable, radially unbounded, positive definite function defined as

$\begin{array}{cc}V=\frac{1}{2}{e}^{T}e+\frac{1}{2}{r}^T{S}^{-1}r+P.& \left(42\right)\end{array}$

After taking the time derivative of (42) utilizing (11), (33) (33), and (39), {dot over (V)}(t) can be expressed as:

$\begin{array}{cc}\stackrel{.}{V}\le -\left({\lambda }_0-\frac{p_2^2}{4{\lambda }_{\min }\left(k_s\right)}\right){z}^2& \left(43\right)\end{array}$ where

${\lambda }_0\stackrel{\Delta }{=}\min \left\{a,1\right\},$ and λ_min(.) denotes the minimum eigenvalue of the argument.

Inequality, (44) can be used to show that V(w,t)εL_∞; h hence, e(t), r(t), P(t)εL_∞. Given that e(t), r(t)εL_∞, (8) can be utilized to show that e(t)εL_∞. Since e(t)εL_∞, (7) can be used along with the assumption that x_m(t)εL_∞ to prove that x(t), {dot over (x)}(t)εL_∞. Based on the fact that x(t)εL_∞. Assumption 1 can be utilized to show that ƒ(x,t)εL_∞. Given that x(t), {dot over (x)}(t), ƒ(x,t)εL_∞, (1) can be used to show that u(t)εL_∞. Since e(t), r(t)εL_∞, the time derivative of (27) and (28) can be used to show that {dot over (μ)}₀(t), {dot over (μ)}₁(t)εL_∞. Given that e(t), r(t){dot over (μ)}₁(t)εL_∞, (33) can be used along with (35) to show that r(t)εL_∞. Since e(t), r(t)εL_∞, (17) can be used to show that z(t) is uniformly continuous. Since z(t) is uniformly continuous, V(w,t) is radially unbounded, and (42) and (43) can be used to show that z(t)εL_∞∩L₂. Barbalat's Lemma can now be invoked to state that∥z(t)∥→0, as t→∞∀w(0)ε²ⁿ⁺¹  (44)

Based on the definition of z(t), (44) can be used to show that∥e(t)∥→0, as t→∞∀w(0)ε²ⁿ⁺¹.

A numerical simulation was created to verify the performance of the control law developed in (2), (7), (18), (27), and (28). The simulation is based on the dynamic model given in (1) and (2), where n=3 and m=6. Specifically, the control input μ_i(t), i=1, 2, . . . , 6 synthetic jet arrays, and the 3-DOF state vector is defined in terms, of the roll, pitch, and yaw rates asx=[x₁x₂x₃]^T

The state and input matrices A and B and reference system matrices A_m and B_m are defined based on the Barron Associates nonlinear tailless aircraft model (BANTAM).

The 3-DOF linearized model for the BANTAM was obtained analytically at the trim condition: Mach number M.=0.455; angle of attack α=2.7 deg, and side slip angle β=0. The linearized dynamic model does not produce the same result as the full nonlinear system with mechanical control surfaces, but the angular, accelerations caused by the virtual surface deflections are predicted accurately using the matrix β. The actual (i.e, Θ_1i* and Θ_2i*, i=1, 2, . . . , 6) and estimated (i.e. {circumflex over (θ)}_1i and {circumflex over (θ)}_2i, i=1, 2, . . . , 6) values of the SJA parameters (see (2) and (7)) are shown in Table 1 above.

The external disturbance used in the simulation is given by

$d=\left[\begin{array}{c}0.2\sin \left(0.5t\right)\\ 0.1\sin \left(0.5t\right)+0.2\sin \left(0.5t\right)\\ 0.15\sin \left(0.3t\right)\end{array}\right]$

The reference input δ(t) used in the simulation is given by

$\begin{array}{cc}\delta \left(t\right)=\left\{\begin{array}{cc}0.5\sin \left(t\right),& 0\le t\le 50\\ 0.5\sin \left(t\right)+\sin \left(2t\right)& 50<t\le 100\end{array}\right\}& \left(45\right)\end{array}$

FIGS. 6 and 7 illustrate the performance of the proposed troller with control gains selected as

• • α=0.3, β=diag {3.3, 0.3, 0.8}
    • k_s=diag {0.10, 0.15, 2.3}

FIG. 6 shows the tracking error response, which demonstrate the rapid convergence of the closed-loop system, and FIG. 7 shows the control commands during the closed-loop operation. The effects of the instantaneous change in the reference input δ(t) (i e., see (45)) at t=50s are apparent in the e₁(t) and e₃(t) plots in FIGS. 6 and 7, but the tracking errors converge to zero quickly in spite of this sudden switch in the reference command. The control command can be maintained within reasonable limits, and the level of chattering is minimal.

Example 2: Management of Limit Cycle Oscillations (LCO)

The equation describing LCO in an airfoil approximated as a 2-dimensional thin plate can be expressed as

$\begin{array}{cc}{M}_s\ddot{p}+C_s\stackrel{.}{p}+F\left(p\right)p=\left[\begin{array}{c}-\mathrm{Lift}\\ \mathrm{Moment}\end{array}\right]& \left(1a\right)\end{array}$ where the coefficients M_s, C_sε^2×2 denote the structural mass and damping matrices, F(p)ε^2×2 is a nonlinear stiffness matrix, and p(t)ε² denotes the state vector. In Equation (1a), p(t) is explicitly defined as

$\begin{array}{cc}p=\left[\begin{array}{c}h\\ a\end{array}\right]& \left(2a\right)\end{array}$ where h(t), α(t)ε denote the plunging [meters] and pitching [radians] displacements describing the LCO effects. Also in Equation (1a), the structural linear mass matrix M_s

$\begin{array}{cc}{M}_s=\left[\begin{array}{cc}m& S_{\alpha }\\ S_{\alpha }& I_{\alpha }\end{array}\right]& \left(3a\right)\end{array}$ where the parameters S_α, I_αε are the static moment and moment of inertia, respectively. The structural linear damping matrix is described as

$\begin{array}{cc}{C}_s=2\left[\begin{array}{cc}{\varsigma }_{h\sqrt{k_{h}m}}& 0\\ 0& {\varsigma }_a\sqrt{k_{\alpha }{I}_{\alpha }}\end{array}\right]& \left(4a\right)\end{array}$ where the parameters _h, _aε are the damping logarithmic decrements for plunging and pitching, and mε is the mass of the wing, or in this case, a flat plate. The nonlinear stiffness matrix utilized is

$\begin{array}{cc}F\left(p\right)=\left[\begin{array}{cc}{k}_h& 0\\ 0& k_{\alpha }+k_{{\alpha }^3}{\alpha }^2\end{array}\right]& \left(5a\right)\end{array}$ where k_α, k_α3ε denote structural resistances to pitching (linear and nonlinear) and k_hε structural resistance to plunging.

In Equation (1a), the total lift and moment are explicitly defined as

$\begin{array}{cc}\left[\begin{array}{c}-\mathrm{Lift}\\ \mathrm{Moment}\end{array}\right]=\left[\begin{array}{c}-\left(L+L_{v_j}\right)\\ \left(M+M_{v_j}\right)\end{array}\right]=M_{\alpha }\ddot{p}+C_{\alpha }\stackrel{.}{p}+K_{\alpha }p+{ℒ}_{\eta }\eta +B_u& \left(6a\right)\end{array}$ where _vj(t), M_vj(t)ε denote the equivalent control force and moment generated by the jth SJA, and (t), M(t)ε are the aerodynamic lift and moment due to the 2-degree-of-freedom motion. In Equation (6a), η(t)ε² denotes the aerodynamic state vector that relates the moment and lift to the structural modes. Also in Equation (6a), u(t)ε² denotes the SJA-based control input (e.g., the SJA air velocity or acceleration), and Bε^2×2 is an uncertain constant input gain matrix that relates the control input u(t) to the equivalent force and moment generated by the SJA. Also in Equation (6a), the aerodynamic and mode matrices M_a, C_a, K_a, L_ηε^2×2 are described as

$\begin{array}{cc}{M}_a=\pi {\mathrm{pb}}^2\left[\begin{array}{cc}-1& \mathrm{ba}\\ \mathrm{ba}& -b^2\left(\frac{1}{8}+a^2\right)\end{array}\right]& \left(7a\right)\\ C_a=\pi {\mathrm{pb}}^2\left[\begin{array}{cc}0& -U\\ 0& -\mathrm{Ub}\left(\frac{1}{2}-a\right)\end{array}\right]+2\pi \mathrm{pUb}\varphi \left(0\right)\left[\begin{array}{cc}-1& -b\left(\frac{1}{2}-a\right)\\ b\left(\frac{1}{2}+a\right)& b^2\left(\frac{1}{2}+a\right)\left(\frac{1}{2}-a\right)\end{array}\right]& \left(8a\right)\\ k_a=2\pi \mathrm{pUb}\varphi \left(0\right)\left[\begin{array}{cc}0& -U\\ 0& b\left(\frac{1}{2}+a\right)U\end{array}\right]& \left(9a\right)\\ L_{\eta }=2\pi \mathrm{pUb}\left[\begin{array}{cc}{a}_1{b}_1& a_2{b}_2\\ -b\left(\frac{1}{2}+a\right)a_1{b}_1& -b\left(\frac{1}{2}+a\right)a_2{b}_2\end{array}\right]& \left(10a\right)\end{array}$ where φ(0) is the Wagner solution function at 0, and the parameters a₁, b₁, a₂, b₂ε are the Wagner coefficients. In addition, a, bε denote the relative locations of the rotational axis from the mid-chord and the semi-chord, respectively. The aerodynamic state variables are governed by{dot over (η)}=C_n{dot over (p)}+K_ηp+S_ηη  (11a)

The aerodynamic state matrices in Equation (11a), C_η, K_η, S_ηε^2×2, are explicitly defined as

$\begin{array}{cc}{C}_{\eta }=\frac{U}{b}\left[\begin{array}{cc}-1& -b\left(\frac{1}{2}-a\right)\\ -1& -b\left(\frac{1}{2}-a\right)\end{array}\right]& \left(12a\right)\\ K_{\eta }=\frac{U}{b}\left[\begin{array}{cc}0& -U\\ 0& -U\end{array}\right]& \left(13a\right)\\ S_{\eta }=\frac{U}{b}\left[\begin{array}{cc}-b_1& 0\\ 0& -b_2\end{array}\right]& \left(14a\right)\end{array}$

By substituting Equation (6a) into Equation (1a), the LCO dynamics can be expressed as{umlaut over (p)}=−M⁻¹C{dot over (p)}−M⁻¹Kp+M⁻¹L_ηη+M⁻¹Bu  (15a)where C=C_s−C_a, K=F(p)−K_a, and M=M_s−M_a. By making the definitions x₁(t)=h(t), x₂(t)=α(t), x₃(t)={dot over (h)}(t), x₄(t)={dot over (α)}(t), x₅(t)=η₁(t), and x₆(t)=η₂(t), the dynamic equation in Equation (15a) can be expressed in state form as{dot over (x)}=A(x)x+Bu  (16a)where x(t)ε⁶ is the state vector, A(x)ε^6×6 is the state matrix (state-dependent). In Equation (16a), the input gain matrix Bε^6×2 is defined as

$\begin{array}{cc}\stackrel{\_}{B}=\left[\begin{array}{c}{0}_{2\times 2}\\ M^{-1}B\\ 0_{2\times 2}\end{array}\right]& \left(17a\right)\end{array}$ where 0_2×2 denotes a 2×2 matrix of zeros. The structure of the input gain matrix in Equation (17a) results from the fact that the control input u(t) only directly affects {umlaut over (h)}(t) and ä(t).

In some embodiments, an objective is to design a control signal u(t) to regulate the plunge and pitching dynamics (i.e., h(t), a(t)) resulting from LCO) to zero. To facilitate the control design, the expression in Equation (15a) is rewritten asM{umlaut over (p)}=g(h,α,η)+Bu  (18a)where g(h, α, η) is an unknown, unmeasurable auxiliary function.

Remark 1. Based on the open-loop error dynamics in Equation (18a), one of the control design challenges is that the control input u(t) is pre-multiplied by the uncertain matrix B. In the following control development and stability analysis, it will be assumed that the matrix B is uncertain, and the robust control law will be designed with a constant feedforward estimate of the uncertain matrix. The simulation results demonstrate the capability of the robust control law to compensate for the input matrix uncertainty without the need for online parameter estimation or function approximators.

To quantify the control objective, a regulation error e₁(t)ε² and auxiliary tracking error variables e₂(t), r(t)ε² are defined ase₁=p−p_d  (19a)e₂=ė₁+α₁e₁  (20a)r=ė₂+α₂e₂  (21a)where α₁, α₂>0ε are user-defined control gains, and the desired plunging and pitching states p_d=[h,a]^T=[0,0]^T for the plunging and pitching suppression objective. To facilitate the following analysis, Equation (21a) is pre-multiplied by M and the time derivative is calculated asM{dot over (r)}=Më₂+α₂Mė₂  (22a)

After using Equations (18a)-(21a), the open-loop error dynamics are obtained asM{dot over (r)}=Ñ+N_d+B{dot over (μ)}−e₂  (23a)where the unknown, unmeasurable auxiliary functions Ñ(e₁, e₂, r), N_d (p_d, _d)ε² are defined as

$\begin{array}{cc}\tilde{N}\stackrel{\Delta }{=}\stackrel{.}{g}\left(p,\eta \right)-\stackrel{.}{g}\left(p_d,\eta \right)+{\alpha }_1\left(r-{\alpha }_2{e}_2-{\alpha }_1{e}_2+{\alpha }_1^2{e}_1\right)+{\alpha }_{2}M\left(r-{\alpha }_2{e}_2\right)+e_2& \left(24a\right)\\ N_d\stackrel{\Delta }{=}-{\stackrel{⃛}{p}}_d+\stackrel{.}{g}\left(p_d,\eta \right)& \left(25a\right)\end{array}$

The motivation for defining the auxiliary functions in Equations (24a) and (25)a is based on the fact that the following inequalities can be developed:∥Ñ∥≦p₀∥z∥,∥N_d∥≦_Nd,∥{dot over (N)}_d∥≦_Nd  (26a)where p₀, _Nd, _Ndε⁺ are known bounding constants, and z(t)ε⁶ is defined as

$\begin{array}{cc}z\stackrel{\Delta }{=}{\left[\begin{array}{ccc}{e}_1^T& e_2^T& r^T\end{array}\right]}^T& \left(27a\right)\end{array}$

Based on the open-loop error dynamics in Equation (23a), the control input is designed via{dot over (u)}={circumflex over (B)}⁻¹(−(k_s+I_2×2)r−βsgn(e₂(t))  (28a)where k_s, βε^2×2 denote constant, positive definite, diagonal control gain matrices, and I_2×2 denotes a 2×2 identity matrix. In Equation (28a), {circumflex over (B)}ε^2×2 denotes a constant, feedforward “best guess” estimate of the uncertain input gain matrix B. The control input u(t) does not depend on the unmeasurable acceleration term r(t), since Equation (28a) can be directly integrated to show that u(t) requires measurements of e₁(t) and e₂(t) only which are the error of pitching and plunging respectively.

To facilitate the following stability proof, the control gain matrix β in Equation (28a) is selected to satisfy the sufficient condition

$\begin{array}{cc}{\lambda }_{\min }\left(\beta \right)>{\varsigma }_{N_d}+\frac{1}{{\alpha }_2}{\varsigma }_{{\stackrel{.}{N}}_d}& \left(29a\right)\end{array}$ where λ_min(.) denotes the minimum eigenvalue of the argument. After substituting Equation (28a) into Equation (23a), the closed-loop error dynamics are obtained asM{dot over (r)}=Ñ+N_d−(k_s+I_n×n)r+βsgn(e₂(t))−e₂  (30a)

To reduce the complexity of the following stability analysis, it is assumed that the product B{circumflex over (B)}⁻¹ is equal to identity. It can be proven that asymptotic regulation can be achieved for the case where the feedforward estimate {circumflex over (B)} is within some prescribed finite range of the actual matrix B.

The following simulation results demonstrate the performance of a controller according to embodiments of the present invention that seek to control LOC in the presence of uncertainty in the input gain matrix B.

Theorem 1. The controller given in Equation (28a) ensures asymptotic regulation of pitching and plunging displacements in the sense that∥e₁(t)∥→0 as t→∞  (31a)provided the control gain k_s is selected sufficiently large, and β is selected according to the sufficient condition in Equation (29a).

Lemma 1. To facilitate the following proof let ⊂⁷ be a domain containing w(t)=0, where w(t)ε⁷ is defined as

$\begin{array}{cc}w\left(t\right)\stackrel{\Delta }{=}{\left[z^T\sqrt{P\left(t\right)}\right]}^T& \left(32a\right)\end{array}$

In Equation (32a), the auxiliary function P(t)ε is the generalized solution to the differential equation{dot over (P)}(t)=−L(t)  (33a)P(0)=β∥e₂(0)∥−N_d^T(0)e₂(0)  (34a)where the auxiliary function L(t)ε is defined asL(t)=r^T(N_d(t)−βsgn(e₂))  (35a)

Provided the sufficient condition in Equation (29a) is satisfied, the following inequality can be obtained:∫₀^tL(r)dr≦β∥e₂(0)∥−N_d^T(0)e₂(0)  (36a)

Hence, Equation (36a) can be used to conclude that P(t)≧0.

Proof 1. (See Theorem 1) Let V(w,t): x[0, ∞)→+ defined as the nonnegative function

$\begin{array}{cc}V\left(w,t\right)\stackrel{\Delta }{=}\frac{1}{2}{e}_1^T{e}_1+\frac{1}{2}{e}_2^T{e}_2+\frac{1}{2}{r}^T\mathrm{Mr}+P& \left(37a\right)\end{array}$ where e₁(t), e₂(t), and r(t) are defined in Equations (19a)-(21a), respectively; and the positive definite function P(t) is defined in Equation (33a). The function V(w,t) satisfies the inequalityU₁(w)≦V(w,t)≦U₂(w)  (38a)provided the sufficient condition introduced in Equation (29a) is satisfied, where U₁(w), U₂(w)ε denote the positive definite functions

$\begin{array}{cc}{U}_1\stackrel{\Delta }{=}{\lambda }_1{w}^2,U_2\stackrel{\Delta }{=}{\lambda }_2{w}^2& \left(39a\right)\end{array}$ where

${\lambda }_1=\min \left\{\frac{1}{2}{\lambda }_{\min }\left(M\right)\right\}$ and λ₂=max{1,λ_max(M)} After taking the time derivative of Equation (37a) and utilizing Equation (20a), Equation (21a), Equation (30a), and Equation (33a), V(w,t) can be upper bounded as

$\begin{array}{cc}\stackrel{.}{V}\left(w,t\right)\le -\left({\alpha }_1-\frac{1}{2}\right){e_1}^2-\left({\alpha }_2-\frac{1}{2}\right){e_2}^2-{r}^2-k_s{r}^2+p_0zr& \left(40a\right)\end{array}$ where the bounds in Equation (26a) were used, and the fact that

$e_1^T{e}_2\le \frac{1}{2}{e_1}^2+\frac{1}{2}{e_2}^2$ (i.e., Young's inequality) was utilized. After completing the squares in Equation (40a), the upper bound on V(w,t) can be expressed as

$\begin{array}{cc}\stackrel{.}{V}\left(w,t\right)\le -\left({\alpha }_1-\frac{1}{2}\right){e_1}^2-\left({\alpha }_2-\frac{1}{2}\right){e_2}^2-{r}^2-{k_s\left(r-\frac{p_0}{2k_s}z\right)}^2+\frac{p_0^2}{4_{\mathrm{ks}}}{z}^2& \left(41a\right)\end{array}$

Since k₃>0, the upper bound in Equation (41a) can be expressed as

$\begin{array}{cc}\stackrel{.}{V}\left(w,t\right)\le -\left({\lambda }_0-\frac{p_0^2}{4k_s}\right){z}^2\mathrm{where}{\lambda }_0\stackrel{\Delta }{=}\min \left\{{\alpha }_1-\frac{1}{2},{\alpha }_2-\frac{1}{2},1\right\}.& \left(42a\right)\end{array}$

The following expression can be obtained from Equation (42a):{dot over (V)}(w,t)≦−U(w)  (43a)where U(w)=c∥z∥², for some positive constant cε, is a continuous, positive semi-definite function.

It follows directly from the Lyapunov analysis that e₁(t), e₂(t), r(t)ε_∞. This implies that ė₁(t), ė₂(t)ε_∞ from the definitions given in Equations (20a) and (21a). Given that ė₁(t), e₂(t), r(t)ε_∞, it follows that ë₁(t)ε_∞ from Equation (21a). Thus, Equation (19a) can be used to prove that p(t), {dot over (p)}(t), {umlaut over (p)}(t)ε_∞. Since p(t), {dot over (p)}(t), {umlaut over (p)}(t)ε_∞, Equation (18a) can be used to prove that u(t)ε_∞. Since r(t), u(t)εL_∞, Equation (28a) can be used to show that {dot over (u)}(t)ε_∞. Given that e₁(t), e₂(t), r(t), {dot over (u)}(t)ε_∞, Equation (30a) can be used along with Equation (26a) to prove that {dot over (r)}(t)ε_∞. Since ė₁(t), ė₂(t), {dot over (r)}(t)ε_∞, e₁(t), e₂(t), r(t), are uniformly continuous. Equation (27a) can then be used to show that z(t) is uniformly continuous. Given that e₁(t), e₂(t), r(t)ε_∞, Equation (37a) and Equation (42a) can be used to prove that z(t)ε_∞∩₂. Barbalat's lemma can now be invoked to prove that ∥z(t)∥→0 as t→∞. Hence, ∥e₁(t)∥→0 as t→∞ from Equation (27a). Further, given that V(w,t) in Equation (37a) is radially unbounded, convergence of e₁(t) is guaranteed, regardless of initial conditions—a global result.

A numerical simulation was created to demonstrate the performance of the control law developed in Equation (28a). In order to develop a realistic stepping stone to high-fidelity numerical simulation results using detailed computational fluid dynamics models, the following simulation results are based on detailed dynamic parameters and specifications. The simulation is based on the dynamic model given in Equation (1a) and Equation (11a). The dynamic parameters utilized in the simulation are summarized in Table 2:

TABLE 2
Constant simulation parameters.
$\rho =1.225\frac{\mathrm{kg}}{m^3}$	a = −0.24	$U=18.25\frac{m}{s},19.5\frac{m}{s},20.5\frac{m}{s}$
m = 2.55 kg	b = 0.11 m	$v=18\frac{m}{s}$
Sa = 10.4 × 10−3 kg · m	a1 = 0.1650	a2 = 0.0455
Ia = 2.51 × 10−3 kg · m	b1 = 0.3350	b2 = 0.3000
$k_h=450\frac{N}{m}$	$k_a=9.3\frac{N}{m}$	$k_{a^3}=55\frac{N}{m}$
ζh = 5.5 × 10−3	ζa = 1.8 × 10−2

The following simulation results were achieved using control gains defined as

$\begin{array}{cc}\beta =\left[\begin{array}{cc}0.0001& 0\\ 0& 25\end{array}\right],k_s=\left[\begin{array}{cc}0.00001& 0\\ 0& 0.11\end{array}\right]& \left(44a\right)\\ a_1=\left[\begin{array}{cc}1& 0\\ 0& 35\end{array}\right],a_2=\left[\begin{array}{cc}1& 0\\ 0& 35\end{array}\right]& \left(45a\right)\end{array}$

The control gains given in Equations (44a) and (45a) were selected based on achieving a desirable response in terms of settling time and required control effort. To test the case where the input gain matrix B is uncertain, it is assumed in the simulation that the actual value of B is the 2×2 identity matrix, but the constant feedforward estimate {circumflex over (B)} used in the control law is given by

$\hat{B}=\left[\begin{array}{cc}0.9& 0.1\\ -0.1& 1.1\end{array}\right]$

FIGS. 8-11 show the open-loop pitching and plunging response of the simulated system for the case where the free stream velocity U is 18.25 m/s and 19.5 m/s, respectively. These figures demonstrate that sustained LCO occur in the absence of feedback control. FIGS. 12 and 13 show the closed-loop pitching and plunging response for the case where U=18.25 m/s, and FIGS. 14 and 15 show the required control force and moment, respectively, in this case. These simulation results demonstrate the capability of the robust control law to asymptotically suppress LCO using a control force and moment that are within reasonable limits. This demonstrates the practical efficacy of the control design to suppress LCO using the limited control authority that is available from the SJA. It should be noted, however, that in actual application, arrays consisting of several SJAs can be utilized to realize a greater overall control force and moment.

FIGS. 16-19 show the closed-loop pitching and plunging response and the required control force and moment, respectively, for the case where the free stream velocity U=19.5 m/s. FIGS. 20-23 show the closed-loop pitching and plunging response and the required control moment and force for the case where U=20.5 m/s. As shown in the figures, an increased control effort is required to suppress the LCO in the presence of the increased air velocity.

FIGS. 24 and 25 show the robustness of the proposed controller in the event of sudden un-modeled disturbances. In these plots, a sudden increase (i.e., step function) in the pitching and plunging rates was programmed to occur at 10.9 s. Simultaneously, the free stream velocity U was increased from 19 m/s to 25 m/s for approximately 0.2 seconds to simulate a disturbance due to a wind gust. The pitching and plunging responses in FIGS. 24 and 25 demonstrate the capability of the control law to compensate for the sudden and unexpected disturbances. FIGS. 26 and 27 show the commanded control force and moment, respectively, of the closed-loop system in response to the disturbance.

Example 3: Numerical Approach

To demonstrate the developed SJA-based robust flight control technology, a high-fidelity numerical approach is examined which employs a modified version of the Implicit Large Eddy Simulation (ILES) Navier-Stokes solver. The following features of the original version of the code are particularly beneficial for the analysis of fluid-structure interaction and its control:

Implicit time marching algorithms (up to 4th-order accurate) are particularly suitable for the low-Reynolds number wall-bounded flows characteristic of MAV airfoils.

High-order spatial accuracy (up to 6th-order accurate) is achieved by use of implicit compact finite-difference schemes, thus making LES resolution attainable with minimum computational expense.

Robustness is achieved through a low-pass Pade-type non-dispersive spatial filter that regularizes the solution in flow regions where the computational mesh is not sufficient to fully resolve the smallest scales. Note that the governing equations are represented in the original unfiltered form used unchanged in laminar, transitional or fully turbulent regions of the flow. The highly efficient Implicit LES (ILES) procedure employs the high-order filter operator in lieu of the standard SGS and heat flux terms, with the filter selectively damping the evolving poorly-resolved high-frequency content of the solution.

Overset grid technique is adopted for geometrically complex configurations, with high-order interpolation maintaining spatial accuracy at overlapping mesh interfaces. The code employs an efficient MPI parallelization that has been successfully utilized on various Beowulf cluster platforms.

The present example employs the developed in the code and successfully validated capability to simulate the coupled aerodynamic and aeroelastic responses of 1-DOF and 2-DOF elastically-mounted airfoils (FIG. 1) and their transition to LCO induced by an impinging gust. The aeroelastic response module has been implemented within the framework of the viscous flow solver. In the numerical formulation, the equations governing the fluid dynamics are essentially coupled with equations governing 1-DOF or 2-DOF airfoil motion so that the fluid and structure are treated as a single dynamical system. In the time-marching procedure, the aerodynamic loads are supplied through Navier-Stokes simulations while the structural response module determines the displacement vector which, in turn, defines the grid motion. At each physical time step, the internal iterative loop achieves the balance of the new airfoil position and the corresponding unsteady flowfield. The loop is efficiently merged with the sub-iterative procedure implemented as part of the flow solver's implicit time marching scheme.

In this example, the airfoil's LCO is induced by an impinging sharp-edge gust, with details of the numerical implementation of the gust-airfoil interaction model (FIG. 28). The model is analytically described in terms of the upwash velocity profile (with the streamwise component u_g(x,t)=0) in Eqn. (1b below).

$\begin{array}{cc}{v}_8\left(x,t\right)=\{\begin{array}{cc}{\varepsilon }_{g}f(t-x\text{/}{u}_{\infty )},& u_{\infty }\left(t-T_g\right)\le x\le u_{\infty }t\\ 0,& \mathrm{otherwise}\end{array}& \left(1b\right)\end{array}$

In numerical simulations, such gust is generated with prescribed duration T_g and the gust amplitude ε_g in the momentum source region located upstream of the airfoil, and undergoes ramp-up and ramp-down phases similar to natural flows as represented by function ƒ in Eqn. (1b).

The present example addresses effectiveness of SJA (e.g., FIG. 3) for active control of unsteady flow over SD7003 low-Re airfoil in presence of a sharp-edge gust, where only the actuator's orifice with a properly defined fluctuating-velocity boundary condition specified by Eqn. (2b) is embedded into the airfoil surface to allow for the synthetic jet to freely interact with the grazing boundary-layer flow.v_SJA(x,t)=A cos(ω_SJAt)  (43)

The present example addresses SJA-based robust control of gust-induced LCO in NACA0012 airfoil. The near-airfoil region of the baseline 649×395×3 O-grid is illustrated in FIG. 29. The freestream conditions are imposed at the farfield boundary located more than 100 chords away from the airfoil, with the grid rapidly stretching towards the boundary to ensure effective elimination of spurious reflections achieved in conjunction with the low-pass spatial filtering. The synthetic-jet actuators embed the actuator's orifice mesh in the airfoil surface and provide an adequate overlap with the original airfoil mesh (FIG. 29). The proper implementation of the employed overset grid methodology involves 7 meshes generated using Pointwise© software in the near-orifice overlap region for each of the two actuators symmetrically positioned on the top and the bottom airfoil surfaces at 0.95 c chordwise locations on each side. The overset grid connectivity is established using NASA'S PEGASUS and AFRL's BELLERO software, with the connectivity data produced by the former serving as input for the latter handling grid decomposition and establishing the intra-grid communication required for the grid system subdivided into blocks for parallel processing.

A fixed time step of Δt≈5×10⁻⁵ is used in the code parallel simulations, with the baseline mesh efficiently partitioned into a set of 572 overlapped blocks assigned to different processors. Such computations require about 10 CPU hours on a DOD HPC system to establish a clearly-defined LCO in approximately 10⁶ time steps.

	TABLE 3
	ρ = 1.1 kg/m3	b = 0.11 m	a = −0.024 m
	m = 2.55 kg	a1 = 0.165	a2 = 0.0455
	Sa = 1.04 × 10−2 kg · m	b1 = 0.335	b2 = 0.300
	Ia = 2.51 × 10−3 kg · m	ka = 9.3 N/m	ka3 = 55 N/m
	ka = 450 N/m	ζh = 5.5 × 10−3	ζa = 1.8 × 10−3

A representative set of the aeroelastic model's parameters shown in Table 3 was selected to provide a realistic model of elastically-mounted NACA0012 wing section. The structural parameters were employed to match a critical (flutter) speed of about 16 m/s. FIG. 30a shows comparison of the pitching experimental results. The corresponding transition to pitching LCO induced by an impinging sinusoidal gust predicted based on the high-accuracy simulations for 18 and 19 m/s is shown in FIG. 30b. Note that the latter reveals the LCO frequency and amplitude shifts in the structural response. However, the comparison of the quasilinear and nonlinear structural response models indicates different transient processes towards establishing similar final LCO patterns with a phase shift. In particular, the response frequency analysis of the quasilinear model reveals a transient region characterized by a high-amplitude LCO at ω_LCO ˜0.6 that abruptly switched to the lower-amplitude LCO at ω_LCO˜0.2, whereas the nonlinear model reached such long-term response through continuous spectral shifts with lower-amplitude transient modes.

FIG. 31 illustrates vorticity and streamline contours for the suction-side actuation based on the numerical experiments performed with actuators embedded on both sides of the airfoil close to the trailing edge.

FIG. 32 first illustrates results of the reduced-order simulations for uncontrolled 2-DOF flat-plate pitching LCO obtained for the flow speed of 19 m/s for three initial airfoil excitation amplitudes corresponding to ε_l={dot over (p)}(0)=(0.002, 0.005), ε_m={dot over (p)}(0)=(0.1, 0.2) and ε_h={dot over (p)}(0)=(1, 2) (the state vector p(t) is defined in Eqn. (2)). Clearly, the final LCO amplitudes are the same in all cases but the transition process is different. The plunging LCO characteristics have very similar features and thus are not shown.

The required control authority of the actuators changes correspondingly depending on the initial excitation and the LCO amplitudes (i.e., the flow speed, as shown in FIG. 30a). Test computations are performed for the following selection of the control gains in the robust controller model:

$\begin{array}{cc}{a}_1=\left[\begin{array}{cc}1& 0\\ 0& 35\end{array}\right],a_2=\left[\begin{array}{cc}1& 0\\ 0& 35\end{array}\right],k_3=\hspace{1em}\left[\begin{array}{cc}{10}^{-7}& 0\\ 0& 9\times {10}^{-4}\end{array}\right],\beta =[\begin{array}{cc}{10}^{-3}& 0\\ 0& 25\end{array}]& \left(3b\right)\end{array}$

Successful suppression of the pitching LCO is demonstrated for the three initial excitation amplitudes in FIG. 33, whereas the corresponding time histories of the aerodynamic lift and moment produced by the actuator governed by the robust controller are shown in FIG. 34. As expected, both the time required to suppress LCO oscillations and the amplitudes of the forces and moments to be delivered by the actuator operating in the feedback-loop robust control system increase with higher initial excitation amplitudes. Similarly, the control authority requirements become more demanding at higher supercritical flight speeds, and generally new optimized sets of control gains should be determined.

A nonlinear robust control law for SJA-based LCO suppression in UAV wings is presented. The control law is rigorously proven to achieve global asymptotic regulation of the pitching and plunging displacements to zero in the presence of dynamic model uncertainty and parametric actuator uncertainty. Furthermore, the control law is shown via numerical simulation to compensate for un-modeled external disturbances (i.e., due to wind gusts and un-modeled effects). It is further shown that the robust control law can achieve suppression of LCO using minimal control effort.

Claims

1. An unmanned aerial vehicle (UAV), comprising: an airframe;a plurality of synthetic jet actuators (SJA) affixed to the airframe, each SJA configured to selectively produce a stream of air in response to a control command input thereto, each SJA being represented by a mathematical model having at least two uncertain parameters;a plurality of sensors configured to determine an operating condition of the UAV; anda control system configured to provide the control command to each of the plurality of synthetic jet actuators based on the operating condition of the UAV and a control law, the control law including a constant estimate for each of the uncertain parameters in the mathematical model of each corresponding SJA;wherein the control system controls at least one or both of the trajectory of the UAV and the vibration of the UAV by activating at least one of the plurality of SJA with the control command.

2. The unmanned aerial vehicle according to claim 1, wherein the plurality of SJA are selectively controlled by the control system to track a predetermined flight path, and wherein the control law comprises:

3. The unmanned aerial vehicle according to claim 2, wherein the plurality of sensors measure the roll rate, pitch rate and yaw rate of the airframe as the operating condition.

4. The unmanned aerial vehicle according to claim 2, wherein: ud(t)={circumflex over (Ω)}#(μ0(t)−μ1(t))

5. The unmanned aerial vehicle according to claim 4, wherein: μ0=−(ks+Im×x)e(t)−(ks+Im×n(0)−∫0ƒα(ks+In×n)e(τ)dτμ1=∫0tβsgn(e(τ))dτ

6. The unmanned aerial vehicle according to claim 1, wherein the plurality of SJA are selectively controlled by the control system to dampen limit cycle oscillations,

7. The control system according to claim 1, wherein the control law is configured to generate the control command for use in flight path tracking, the control law comprises:

8. A method of controlling a micro air vehicle, comprising: measuring at least one of roll rate, pitch rate, yaw rate, pitching rate, and plunging rate;comparing the measured value with a desired reference value to determine an error;determining a control command voltage based on a control law, the control law comprising:

9. The method of controlling a micro air vehicle according to claim 8, wherein the method is free from time-varying adaptive parameter estimation algorithms.