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.

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; and a 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: $v\left(t\right)=\frac{{\hat{}}_1}{{\hat{}}_2-u_d\left(t\right)}$ where {circumflex over ()}.sub.1 and {circumflex over ()}.sub.2 are constant estimates of uncertain parameters, and u.sub.d(t) denotes an auxiliary control term that incorporates the operating condition of the UAV.

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:
u.sub.d(t)={circumflex over ()}.sup.#(.sub.0(t).sub.1(t)) where {circumflex over ()}.sup.# denotes the pseudoinverse of {circumflex over ()}, and {circumflex over ()} is a constant estimate of an uncertain matrix, and where .sub.0(t) and .sub.1(t) are functions of roll, pitch, and yaw rate tracking errors.

5. The unmanned aerial vehicle according to claim 4, wherein:
.sub.0=(k.sub.s+I.sub.mx)e(t)(k.sub.s+I.sub.mn(0).sub.0.sup.(k.sub.s+I.sub.nn)e()d
.sub.1=.sub.0.sup.tsgn(e())d where: e(t) is a trajectory tracking error expressed in degrees per second, {circumflex over ()} is a constant estimate of the parametric uncertainty, and k.sub.s, , are positive, constant, user-defined amplifier settings.

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, wherein the plurality of sensors measure at least one of pitching rate and plunging rate as the operating condition, and wherein the control law comprises:
{dot over (u)}={circumflex over (B)}.sup.1((k.sub.s+I.sub.22)rsgn(e.sub.2(t)) where k.sub.s,.sup.22 denote constant, positive definite, diagonal control gain matrices, I.sub.22 denotes a 22 identity matrix, r and e.sub.2 denote auxiliary regulation error measurements, and {circumflex over (B)}.sup.22 denotes a constant estimate of the uncertain input gain matrix B.

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: $v\left(t\right)=\frac{{\hat{}}_1}{{\hat{}}_2-u_d\left(t\right)}$ where {circumflex over ()}.sub.1 and {circumflex over ()}.sub.2, are constant estimates and u.sub.d (t)={circumflex over ()}.sup.#(.sub.0(t).sub.1 (t)) where {circumflex over ()}.sup.# denotes the pseudoinverse of {circumflex over ()}, and {circumflex over ()} is a constant estimate of an uncertain matrix, and where .sub.0(t) and .sub.1(t) are functions of roll, pitch, and yaw rate tracking errors.

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: $v\left(t\right)=\frac{\hat{}}{{\hat{}}_2-u_d\left(t\right)}$ where {circumflex over ()}.sub.1 and {circumflex over ()}.sub.2 are constant estimates and u.sub.d(t) is based on the error; 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.

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.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

(2) 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.

(3) FIG. 3 is a schematic example of a SJA.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(29) FIG. 29 illustrates example details of a numerical grid.

(30) FIG. 30a illustrates flat plate pitching LCO amplitudes.

(31) FIG. 30b illustrates time histories of airfoil pitching responses.

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

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

(34) 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.

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

DISCUSSION OF INVENTION

(36) 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.

(37) 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.

(38) 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.sup.nn denotes the uncertain state matrix, B.sup.nn represents the uncertain input matrix, and (x,t).sup.n 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

(39) $u\left(t\right)\stackrel{}{=}{\left[u_1\left(t\right)u_2\left(t\right)\mathrm{.Math.}{u}_m\left(t\right)\right]}^{\mathrm{��}}{}^m$
represents the virtual surface deflections due to m arrays of synthetic jet, actuators (SJA).

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

(41) $\begin{array}{cc}{u}_i={}_{2i}^{*}-\frac{{}_{1i}^{*}}{v_i}i=1,2,\mathrm{.Math.},m& \left(2\right)\end{array}$
where .sub.i()=() denotes the peak-to-peak voltage acting, on the i.sup.th SJA array, and .sub.1i*, .sub.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.sub.i(t) is the peak-to-peak control input voltage to the SJA array expressed in volts for each SJA i=1n. .sub.1i*, .sub.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 .sub.1i*=33.33 volt-deg and .sub.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.

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

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

(44) In Equation (3), {circumflex over ()}.sub.1, and .sub.2* denote constant, best-guess estimates of the uncertain parameters .sub.1* and .sub.2*, 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 ()}.sub.1 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 .sub.1* is

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

(46) Also in Equation (3), u.sub.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.sub.d(t) is pre-multiplied by an uncertain term (a matrix in general) containing the uncertain parameter(s) .sub.1*, the constant estimate(s) {circumflex over ()}.sub.1, and the elements of the uncertain input matrix B. It can further be shown that the elements of the uncertain matrix contain the ratios .sub.1*/{circumflex over ()}.sub.1 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).

(47) By incorporating the estimate {circumflex over ()}, the control term u.sub.d(t) is designed as
u.sub.d(t)={circumflex over ()}.sup.#(.sub.0(t).sub.1(t))(4)
where {circumflex over ()}.sup.# denotes the pseudoinverse of {circumflex over ()}, .sub.0(t) and .sub.1(t) are functions of the UAV roll, pitch, and yaw rate tracking errors defined as
.sub.0=(k.sub.s+I.sub.nn)e(t)(k.sub.s+I.sub.nn)e(0).sub.0.sup.t(k.sub.s+I.sub.nn)e()d
.sub.1=.sub.0.sup.tsgn(e())d.(5)

(48) 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).fwdarw.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.sub.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.

(49) 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.

(50) Example 1: Tracking Control Embodiments

(51) Using the expressions in (2) and (1), the dynamics can be expressed in terms of the i.sup.th SJA array as
{dot over (x)}=Ax+.sub.t=1.sup.m.sub.i.sub.i+(x,t).(6)

(52) In (6),

(53) $b_i\stackrel{}{=}{\left[b_{1i}{b}_{2i}\mathrm{.Math.}{b}_{\mathrm{ni}}\right]}^{\mathrm{��}}{}^{n}i=1,2,\mathrm{.Math.},m,$
where b.sub.ij by denotes the (i, j).sup.th element of the matrix B in (1).

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

(55) 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.sub.i(t) is utilized, which contains constant feedforward best-guess estimates of the uncertain parameters .sub.1i* and .sub.2i*. The robust inverse that compensates for the uncertain, jet array nonlinearities in (2) can be expressed as

(56) $\begin{array}{cc}{v}_i\left(t\right)=\frac{{\hat{}}_{1i}}{{\hat{}}_{2i}-u_{2i}\left(t\right)},i=1,\mathrm{.Math.},m& \left(7\right)\end{array}$
where {circumflex over ()}.sub.1i, {circumflex over ()}.sub.2i.sup.+ are constant feedforward estimates of .sub.1i* and .sub.2i*, respectively, and u.sub.di(t)i=1, . . . , m are subsequently defined auxiliary control signals.

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

(58) $\begin{array}{cc}{u}_{\mathrm{di}}\left(t\right)=\{\begin{array}{cc}{\hat{}}_{2i}-& \mathrm{if}g\left({}_0\left(t\right),{}_1\left(t\right)\right){\hat{}}_{2i}-\\ g\left({}_0\left(t\right),{}_1\left(t\right)\right)& \mathrm{otherwise}\end{array}& \left(8\right)\end{array}$
where .sup.+ is a small parameter, g(.) is a subsequently defined function, and .sub.0(t), .sub.1(t).sup.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.sub.i(t).

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

(60) Remark 2. The auxiliary control signal u.sub.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 ()}.sub.ji.sub.ji*, j=1, 2.

(61) 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)}.sub.m=A.sub.mx.sub.m+B.sub.m(9)
where x.sub.m(t).sup.n denotes the reference state, A.sub.m.sup.nn is a Hurwitz state matrix, B.sub.m.sup.n denotes the reference input matrix, and (t) is the reference input signal. The reference model in (9) is designed to exhibit desirable performance characteristics.

(62) Assumption 2: The model reference state x.sub.m(t) is bounded and sufficiently smooth in the sense that x.sub.m(t), {dot over (x)}.sub.m(t), {umlaut over (x)}.sub.m(t), (t).sub.t0.

(63) To quantify the control objective, the e(t).sup.n is defined as
e=xx.sub.m(10)

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

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

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

(67) 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.sub.i(t) is generated using (2) and (7); thus, the simulation model includes the complete actuator dynamics.

(68) The expression in (12) can be rewritten as
{dot over (r)}=+N.sub.d+{dot over ()}.sub.d(t)Se(13)
where .sup.nm denotes a constant uncertain matrix, S.sup.nn is a subsequently defined uncertain matrix and the control vector

(69) $u_d\left(t\right)\stackrel{}{=}{\left[u_{d1}\left(t\right),u_{d2}\left(t\right)\mathrm{.Math.}{u}_{dm}\left(t\right)\right]}^{}{}^m.$

(70) In (13), the unknown, immeasurable auxiliary functions (t) and N.sub.d(t) are defined as

(71) 0$\begin{array}{cc}\tilde{N}\stackrel{}{=}A\stackrel{.}{e}+\stackrel{.}{e}+\mathrm{��}e+\left(\stackrel{.}{f}\left(x,t\right)-\stackrel{.}{f}\left(x_m,t\right)\right)& \left(14\right)\\ N_d\stackrel{}{=}A{\stackrel{.}{x}}_m-{\stackrel{\mathrm{.Math.}}{x}}_m+f\left(x_m,t\right)& \left(15\right)\end{array}$

(72) 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:
.sub.0z,N.sub.dN.sub.d,{dot over (N)}.sub.dN.sub.d(16)
where .sub.0, .sub.N.sub.1, .sub.N.sub.d.sup.+ are known bounding constants and z(t).sup.2n is defined as

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

(74) Based on the open-loop error system in (13), the auxiliary control u.sub.d, (t) is designed as
u.sub.d(t)={circumflex over ()}.sup.#(.sub.0.sub.1)(18)
where {circumflex over ()}.sup.nm is a constant, best-guess estimate of the uncertain matrix , and [.].sup.# denotes the pseudoinverse of a matrix. In (18), .sub.0(t), .sub.1(t).sup.n 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.sub.d+{tilde over ()}({dot over ()}.sub.0{dot over ()}.sub.1)Se(19)
where the constant uncertain matrix {tilde over ()}.sup.nn is defined as
{tilde over ()}={circumflex over ()}.sup.#.(20)

(75) 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.sup.nn is a positive definite symmetric matrix; and T.sup.nn is a unity upper triangular matrix, which is diagonally dominant in the sense that
[T.sub.ii].sub.k=i+1.sup.n[T.sub.ik]Q,i=1, . . . ,n1.(22)

(76) In inequalities (22), (0,1) and Q.sup.+ are known bounding constants, and T.sub.ik denotes the (i, k).sup.th element of the matrix T.

(77) 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.sub.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 ()}.sub.j1 and {circumflex over ()}.sub.j2j=1, . . . , 6 differ from the mean values of the actual parameters .sub.j1*, and .sub.j2*, j=1, . . . , 6 as mean {circumflex over ()}.sub.j1=14.05, mean .sub.j1*, =25.28 mean {circumflex over ()}.sub.j2=7.25, mean .sub.j2*=12.08

(78) The values for {circumflex over ()}.sub.j1 and {circumflex over ()}.sub.j2j=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.

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

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

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

(82) Since S is positive definite, .sub.1(t) and N.sub.d1(t) satisfy the inequalities
.sub.1.sub.1z,N.sub.d1.sub.N.sub.d1{dot over (N)}.sub.d1.sub.N.sub.d1(25)
where .sub.1, .sub.N.sub.d1, .sub.N.sub.d1.sup.+ 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.sup.1{dot over (r)}=.sub.1+N.sub.d1+{dot over ()}.sub.0+T{dot over ()}.sub.oT{dot over ()}.sub.1e(26)
where

(83) $\stackrel{\_}{T}\stackrel{}{=}T-I_{\mathrm{nxn}}$
is a strictly upper triangular matrix, and I.sub.nn denotes an nn identity matrix. Based on the open-loop error system in (26), the auxiliary control terms .sub.0(t) and .sub.1(t) are designed as

(84) $\begin{array}{cc}{}_0=-\left(k_s+I_{\mathrm{nxn}}\right)e\left(t\right)-\left(k_s+I_{\mathrm{nxn}}\right)e\left(0\right)-{}_0^t\left(k_s+I_{\mathrm{nxn}}\right)e\left(r\right)\mathrm{dr}& \left(27\right)\\ {}_1={}_0^t\mathrm{sgn}\left(e\left(r\right)\right)\mathrm{dr}& \left(28\right)\end{array}$
where is a constant, positive control gain, k.sub.s.sup.nn 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 as
S.sup.1{dot over (r)}=.sub.1+T{dot over ()}.sub.o+N.sub.d1(k.sub.s+I.sub.nn)rT{dot over ()}.sub.1e.(29)

(85) After taking the time derivative of (27), the term T{dot over ()}.sub.o can be expressed as

(86) $\begin{array}{cc}\stackrel{\_}{T}{\stackrel{.}{}}_0=\left[\begin{array}{c}{.Math.}_{j=2}^n{\stackrel{\_}{T}}_{1j\stackrel{.}{}0j}\\ {.Math.}_{j=3}^n{\stackrel{\_}{T}}_{2j\stackrel{.}{}0j}\\ \mathrm{.Math.}\\ {\stackrel{\_}{T}}_{\left(n-1\right)n\stackrel{.}{}\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

(87) $A_p\stackrel{}{=}{\left[A_{p1}{A}_{p2}\mathrm{.Math.}{A}_{p\left(n-1\right)}\right]}^T{}^{n-1},$
with the individual elements defined as

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

(89) i1, . . . , n1 where the subscript j indicates the jth element of the vector. Based on the definitions in (27) and (30), A.sub.p satisfies, the inequality
A.sub.p.sub.1z(32)

(90) Remark 5: Note that based on (30) and (31), the bounding constant .sub.1 depends only on elements i+1 to n of the control gain matrix k.sub.s due to the strictly upper triangular nature of T. Thus, the element {dot over ()}.sub.01(t) of the control vector {dot over ()}.sub.01(t) does not appear in the term A.sub.p. This fact will be utilized in the subsequent stability proof.

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

(92) $\begin{array}{cc}{S}^{-1}\stackrel{.}{r}={\tilde{N}}_2+N_{d1}-\left(k_s+I_{\mathrm{nxn}}\right)-T{\stackrel{.}{}}_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}$

(93) Based on (25), (32), and (34), N.sub.2(t) satisfies the inequality
N.sub.2.sub.2z(35)

(94) where .sub.2 is a known bounding constant.

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

(96) $\begin{array}{cc}>\frac{1}{\mathrm{.Math.}}\left({}_{N_{d1}}+\frac{1}{}{}_{{\stackrel{.}{N}}_{d1}}\right)& \left(36\right)\end{array}$
where .sub.N.sub.d1 and .sub.N.sub.d1 are introduced in (25), and is introduced in (22).

(97) To facilitate the subsequent stability analysis, let D.sup.2n+1 be a domain containing (t)=0, where (t).sup.2n+1 is defined as

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

(99) 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.sup.T(0)N.sub.d1.sup.T(0)(38)
where the auxiliary function L(t) is defined as
L(t)=r.sup.T(N.sub.d1(t)T{dot over ()}.sub.1).(39)

(100) Lemma 1: Provided the sufficient gain condition in (36) is satisfied, the following inequality can be obtained:
.sub.o.sup.tL()drQ|e(0)|e.sup.T(o)N.sub.d1(0)(40)

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

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

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

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

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

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

(107) ${}_0\stackrel{}{=}\min \left\{a,1\right\},$
and .sub.min(.) denotes the minimum eigenvalue of the argument.

(108) Inequality, (44) can be used to show that V(w,t)L.sub.; h hence, e(t), r(t), P(t)L.sub.. Given that e(t), r(t)L.sub., (8) can be utilized to show that e(t)L.sub.. Since e(t)L.sub., (7) can be used along with the assumption that x.sub.m(t)L.sub. to prove that x(t), {dot over (x)}(t)L.sub.. Based on the fact that x(t)L.sub.. Assumption 1 can be utilized to show that (x,t)L.sub.. Given that x(t), {dot over (x)}(t), (x,t)L.sub., (1) can be used to show that u(t)L.sub.. Since e(t), r(t)L.sub., the time derivative of (27) and (28) can be used to show that {dot over ()}.sub.0(t), {dot over ()}.sub.1(t)L.sub.. Given that e(t), r(t){dot over ()}.sub.1(t)L.sub., (33) can be used along with (35) to show that r(t)L.sub.. Since e(t), r(t)L.sub., (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.sub.L.sub.2. Barbalat's Lemma can now be invoked to state that
z(t).fwdarw.0, as t.fwdarw.w(0).sup.2n+1(44)

(109) Based on the definition of z(t), (44) can be used to show that
e(t).fwdarw.0, as t.fwdarw.w(0).sup.2n+1.

(110) 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 .sub.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 as
x=[x.sub.1x.sub.2x.sub.3].sup.T

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

(112) 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, .sub.1i* and .sub.2i*, i=1, 2, . . . , 6) and estimated (i.e. {circumflex over ()}.sub.1i and {circumflex over ()}.sub.2i, i=1, 2, . . . , 6) values of the SJA parameters (see (2) and (7)) are shown in Table 1 above.

(113) The external disturbance used in the simulation is given by

(114) $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]$

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

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

(117) 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.sub.s=diag {0.10, 0.15, 2.3}

(118) 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.sub.1(t) and e.sub.3(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.

(119) Example 2: Management of Limit Cycle Oscillations (LCO)

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

(121) $\begin{array}{cc}{M}_s\stackrel{\mathrm{.Math.}}{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.sub.s, C.sub.s.sup.22 denote the structural mass and damping matrices, F(p).sup.22 is a nonlinear stiffness matrix, and p(t).sup.2 denotes the state vector. In Equation (1a), p(t) is explicitly defined as

(122) $\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.sub.s

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

(124) $\begin{array}{cc}{C}_s=2\left[\begin{array}{cc}{}_{h\sqrt{k_{h}m}}& 0\\ 0& {}_a\sqrt{k_{}{I}_{}}\end{array}\right]& \left(4a\right)\end{array}$
where the parameters .sub.h, .sub.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

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

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

(127) $\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_{}\stackrel{\mathrm{.Math.}}{p}+C_{}\stackrel{.}{p}+K_{}p+{}_{}+B_u& \left(6a\right)\end{array}$
where .sub.vj(t), M.sub.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).sup.2 denotes the aerodynamic state vector that relates the moment and lift to the structural modes. Also in Equation (6a), u(t).sup.2 denotes the SJA-based control input (e.g., the SJA air velocity or acceleration), and B.sup.22 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.sub.a, C.sub.a, K.sub.a, L.sub..sup.22 are described as

(128) $\begin{array}{cc}{M}_a={\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={\mathrm{pb}}^2\left[\begin{array}{cc}0& -U\\ 0& -\mathrm{Ub}\left(\frac{1}{2}-a\right)\end{array}\right]+2\mathrm{pUb}\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\mathrm{pUb}\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_{}=2\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.sub.1, b.sub.1, a.sub.2, b.sub.2 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.sub.n{dot over (p)}+K.sub.p+S.sub.(11a)

(129) The aerodynamic state matrices in Equation (11a), C.sub., K.sub., S.sub..sup.22, are explicitly defined as

(130) $\begin{array}{cc}{C}_{}=\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_{}=\frac{U}{b}\left[\begin{array}{cc}0& -U\\ 0& -U\end{array}\right]& \left(13a\right)\\ S_{}=\frac{U}{b}\left[\begin{array}{cc}-b_1& 0\\ 0& -b_2\end{array}\right]& \left(14a\right)\end{array}$

(131) By substituting Equation (6a) into Equation (1a), the LCO dynamics can be expressed as
{umlaut over (p)}=M.sup.1C{dot over (p)}M.sup.1Kp+M.sup.1L.sub.+M.sup.1Bu(15a)
where C=C.sub.sC.sub.a, K=F(p)K.sub.a, and M=M.sub.sM.sub.a. By making the definitions x.sub.1(t)=h(t), x.sub.2(t)=(t), x.sub.3(t)={dot over (h)}(t), x.sub.4(t)={dot over ()}(t), x.sub.5(t)=.sub.1(t), and x.sub.6(t)=.sub.2(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).sup.6 is the state vector, A(x).sup.66 is the state matrix (state-dependent). In Equation (16a), the input gain matrix B.sup.62 is defined as

(132) $\begin{array}{cc}\stackrel{\_}{B}=\left[\begin{array}{c}{0}_{22}\\ M^{-1}B\\ 0_{22}\end{array}\right]& \left(17a\right)\end{array}$
where 0.sub.22 denotes a 22 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).

(133) 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 as
M{umlaut over (p)}=g(h,,)+Bu(18a)
where g(h, , ) is an unknown, unmeasurable auxiliary function.

(134) 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.

(135) To quantify the control objective, a regulation error e.sub.1(t).sup.2 and auxiliary tracking error variables e.sub.2(t), r(t).sup.2 are defined as
e.sub.1=pp.sub.d(19a)
e.sub.2=.sub.1+.sub.1e.sub.1(20a)
r=.sub.2+.sub.2e.sub.2(21a)
where .sub.1, .sub.2>0 are user-defined control gains, and the desired plunging and pitching states p.sub.d=[h,a].sup.T=[0,0].sup.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 as
M{dot over (r)}=M.sub.2+.sub.2M.sub.2(22a)

(136) After using Equations (18a)-(21a), the open-loop error dynamics are obtained as
M{dot over (r)}=+N.sub.d+B{dot over ()}e.sub.2(23a)
where the unknown, unmeasurable auxiliary functions (e.sub.1, e.sub.2, r), N.sub.d (p.sub.d, .sub.d).sup.2 are defined as

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

(138) 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.sub.0z,N.sub.d.sub.N.sub.d,{dot over (N)}.sub.d.sub.N.sub.d(26a)
where p.sub.0, .sub.N.sub.d, .sub.N.sub.d.sup.+ are known bounding constants, and z(t).sup.6 is defined as

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

(140) Based on the open-loop error dynamics in Equation (23a), the control input is designed via
{dot over (u)}={circumflex over (B)}.sup.1((k.sub.s+I.sub.22)rsgn(e.sub.2(t))(28a)
where k.sub.s, .sup.22 denote constant, positive definite, diagonal control gain matrices, and I.sub.22 denotes a 22 identity matrix. In Equation (28a), {circumflex over (B)}.sup.22 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.sub.1(t) and e.sub.2(t) only which are the error of pitching and plunging respectively.

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

(142) $\begin{array}{cc}{}_{\min }\left(\right)>{}_{N_d}+\frac{1}{{}_2}{}_{{\stackrel{.}{N}}_d}& \left(29a\right)\end{array}$
where .sub.min(.) denotes the minimum eigenvalue of the argument. After substituting Equation (28a) into Equation (23a), the closed-loop error dynamics are obtained as
M{dot over (r)}=+N.sub.d(k.sub.s+I.sub.nn)r+sgn(e.sub.2(t))e.sub.2(30a)

(143) To reduce the complexity of the following stability analysis, it is assumed that the product B{circumflex over (B)}.sup.1 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.

(144) 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.

(145) Theorem 1. The controller given in Equation (28a) ensures asymptotic regulation of pitching and plunging displacements in the sense that
e.sub.1(t).fwdarw.0 as t.fwdarw.(31a)
provided the control gain k.sub.s is selected sufficiently large, and is selected according to the sufficient condition in Equation (29a).

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

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

(148) 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.sub.2(0)N.sub.d.sup.T(0)e.sub.2(0)(34a)
where the auxiliary function L(t) is defined as
L(t)=r.sup.T(N.sub.d(t)sgn(e.sub.2))(35a)

(149) Provided the sufficient condition in Equation (29a) is satisfied, the following inequality can be obtained:
.sub.0.sup.tL(r)dre.sub.2(0)N.sub.d.sup.T(0)e.sub.2(0)(36a)

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

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

(152) $\begin{array}{cc}V\left(w,t\right)\stackrel{}{=}\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.sub.1(t), e.sub.2(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 inequality
U.sub.1(w)V(w,t)U.sub.2(w)(38a)
provided the sufficient condition introduced in Equation (29a) is satisfied, where U.sub.1(w), U.sub.2(w) denote the positive definite functions

(153) 0$\begin{array}{cc}{U}_1\stackrel{}{=}{}_1{.Math.w.Math.}^2,U_2\stackrel{}{=}{}_2{.Math.w.Math.}^2& \left(39a\right)\end{array}$
where

(154) ${}_1=\min \left\{\frac{1}{2}{}_{\min }\left(M\right)\right\}$
and .sub.2=max{1,.sub.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

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

(156) $e_1^T{e}_2\frac{1}{2}{.Math.e_1.Math.}^2+\frac{1}{2}{.Math.e_2.Math.}^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

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

(158) Since k.sub.3>0, the upper bound in Equation (41a) can be expressed as

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

(160) The following expression can be obtained from Equation (42a):
{dot over (V)}(w,t)U(w)(43a)
where U(w)=cz.sup.2, for some positive constant c, is a continuous, positive semi-definite function.

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

(162) 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:

(163) TABLE-US-00002 TABLE 2 Constant simulation parameters. $=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}$ S.sub.a = 10.4 10.sup.3 kg .Math. m a.sub.1 = 0.1650 a.sub.2 = 0.0455 I.sub.a = 2.51 10.sup.3 kg .Math. m b.sub.1 = 0.3350 b.sub.2 = 0.3000 $k_h=450\frac{N}{m}$ 0$k_a=9.3\frac{N}{m}$ $k_{a^3}=55\frac{N}{m}$ .sub.h = 5.5 10.sup.3 .sub.a = 1.8 10.sup.2

(164) The following simulation results were achieved using control gains defined as

(165) $\begin{array}{cc}=\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}$

(166) 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 22 identity matrix, but the constant feedforward estimate {circumflex over (B)} used in the control law is given by

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

(168) 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.

(169) 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.

(170) 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.

(171) Example 3: Numerical Approach

(172) 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:

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

(174) 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.

(175) 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.

(176) 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.

(177) 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.

(178) 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.sub.g(x,t)=0) in Eqn. (1b below).

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

(180) In numerical simulations, such gust is generated with prescribed duration T.sub.g and the gust amplitude .sub.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).

(181) 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.sub.SJA(x,t)=A cos(.sub.SJA.sub.t)(43)

(182) The present example addresses SJA-based robust control of gust-induced LCO in NACA0012 airfoil. The near-airfoil region of the baseline 6493953 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.

(183) A fixed time step of t510.sup.5 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.sup.6 time steps.

(184) TABLE-US-00003 TABLE 3 = 1.1 kg/m.sup.3 b = 0.11 m a = 0.024 m m = 2.55 kg a.sub.1 = 0.165 a.sub.2 = 0.0455 S.sub.a = 1.04 10.sup.2 kg .Math. m b.sub.1 = 0.335 b.sub.2 = 0.300 I.sub.a = 2.51 10.sup.3 kg .Math. m k.sub.a = 9.3 N/m k.sub.a.sup.3 = 55 N/m k.sub.a = 450 N/m .sub.h = 5.5 10.sup.3 .sub.a = 1.8 10.sup.3

(185) 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 .sub.LCO 0.6 that abruptly switched to the lower-amplitude LCO at .sub.LCO0.2, whereas the nonlinear model reached such long-term response through continuous spectral shifts with lower-amplitude transient modes.

(186) 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.

(187) 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 .sub.l={dot over (p)}(0)=(0.002, 0.005), .sub.m={dot over (p)}(0)=(0.1, 0.2) and .sub.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.

(188) 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:

(189) $\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=\left[\begin{array}{cc}{10}^{-7}& 0\\ 0& 9{10}^{-4}\end{array}\right],=[\begin{array}{cc}{10}^{-3}& 0\\ 0& 25\end{array}]& \left(3b\right)\end{array}$

(190) 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.

(191) 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.