Time-varying Weights in Multi-Objective Optimal Control for Flexible Wing Aircraft

Time-Varying Weights in Multi-Objective Optimal Control for

Flexible Wing Aircraft

∗ Kelley E. Hashemi Stinger Ghaffarian Technologies, Inc., Moffett Field, CA, 94035 † Nhan T. Nguyen NASA Ames Research Center, Moffett Field, CA, 94035 ‡ Michael Drew Stinger Ghaffarian Technologies, Inc., Moffett Field, CA, 94035 A multi-objective optimal control technique is modified to accommodate changing cost function weights and is used to control a flexible wing aircraft model. Variation of the weights is used to adjust the relative importance of each objective according to either a prescribed function of time or of the state. Several techniques for obtaining a practical approximation to the optimal control solution are presented, and stability of a specific weight structure with the optimal controller is demonstrated. Functionality of the multi-objective control design with weight variation is demonstrated in simulation of a flexible wing transport aircraft and is shown to improve performance over the fixed weight version both at a constant flight condition and across changing flight conditions.

Nomenclature u = nominal control law n u = multi-objective control law m M = Wing root bending moment, ft-lb y ∆ C = change in drag coefficient D y = Accelerometer output, ft/s w = Gust disturbance vector γ = flight path angle, deg ψ = heading angle, deg Q , R = multi-objective weighting matrices q , q , q = multi-objective scalar weights e M D ˆ · = estimated quantity W = Riccati equation matrix I. Introduction ulti-objective optimal control relies on the careful selection of weights to balance the relative importance of M each objective and generate a useful effect. Poorly chosen weights can lead to issues such as dominance of a single objective or even competition among the objectives that results in little control effort. This paper explores the possibility of using time-varying weights to improve performance of the optimal control design as well as to respond appropriately to changes in environment that motivate renewed weight selection.

The topic of time variation in the context of optimal control has been widely addressed [ 1 – 5 ]. Here the focus is instead on handling variation of several weights in a cost function that combines multiple performance goals. Specifically, weight variation is used to improve upon the performance of a multi-objective optimal control design for a flexible ∗ Research Engineer, Intelligent Systems Division, AIAA Member † Technical Group Lead and Research Scientist, Intelligent Systems Division, AIAA Associate Fellow ‡ Research Engineer, Intelligent Systems Division, AIAA Member

wing aircraft that attempts to simultaneously provide flexible motion suppression, load alleviation, and drag reduction.

Fixed-weight versions of related multi-objective controllers have been previously considered [ 6 , 7 ]. Techniques for both time-dependent and state-dependent variation of the weights are presented. Implementation practicalities are discussed as each technique is demonstrated in a flexible aircraft control simulation.

The structure of the paper is as follows: Section II described the aircraft used for simulation of the weight variation techniques. The structure of the problem is presented in Section III, including a description of the plant dynamics, performance metrics, and control design. Section IV discuses time-based weight variation and contains simulation results demonstrating control of the aircraft across changing flight conditions. Section V considers state-based weight variation and discusses several techniques for approximating the optimal solution including polynomial weight structure that lends itself to closed-form gain expressions. Simulation results demonstrating improved performance at a single flight condition are presented for each. Stability of a controller based on a specific polynomial weight selection is demonstrated via Lyapunov. Finally, some concluding remarks are made in Section VI.

II. Aircraft Description All simulations are carried out on linear, 6-DoF, aeroservoelastic models of the aspect ratio 13.5 Common Research Model (CRM). The model mimics a modern, transonic, wide-body commercial transport aircraft. The aircraft wings are modeled as less structurally rigid than currently operational transport aircraft to investigate the impact of construction from new materials. The model is equipped with a suite of mini flaps along the trailing edge of the wing which can be used for both aeroelastic stabilization and to potentially exploit wing flexibly for performance gain through wing shaping. The model is also equipped with accelerometers along the leading edge of the wing. The wing sensor and actuator layout is detailed in Fig. 1. Further details on an earlier iteration of this aircraft model can be found in [8].

Fig. 1 Control surface and accelerometer layout for CRM 13.5 wing (not drawn to scale) The CRM simulations include some realistic effects such as linear actuator dynamics and gust disturbance’s vertical velocity and acceleration. The ability to calculate control gains on rigid body dynamics only or a reduced version of the model that excludes any user-defined subset of flexible modes is also available. An observer is used to provide estimates of all unmeasured states, and the gust disturbance is estimated from accelerometer measurements.

III. Problem Formulation The goal of the control design here is to reject a discrete gust disturbance and maintain steady level flight at cruise flight conditions well below the flutter boundary. Secondary performance goals include flexible mode suppression, load alleviation, and drag reduction. Two separate control laws are calculated to achieve these objectives. The nominal control law, u , provides heading and flight path angle tracking. The multi-objective control law, u , attempts to deliver n m balanced performance of secondary tasks and is calculated after the nominal loop has been closed.

A. Plant description The gust-disturbed dynamics of the aircraft subject to the two aforementioned control laws are taken to be of the form ˙ x = Ax + B r + B u + B u + E w (1) i n n m m The vector w ∈ R represents the the gust disturbance via its vertical velocity and acceleration. The r term is only included if a servomechanism is necessary for command tracking, and integrator states are appended to the state vector n T T T T as appropriate. All matrices are fully known. The the state x ∈ R is partitioned as x = [ x x x ] where the rigid r e i body states x are available directly, but the elastic states x must be estimated if they are to be used in the control r e implementation. Any other necessary integrator or filter states are contained in x . Accelerometer measurements i n o y ∈ R are available for estimation of the elastic state ˆ x . The additional performance metrics are captured by the wing e root bending moment M ∈ R y M = M x + M u + M u + M w (2) y x u n u m w n m as well as an approximate expression for the change in the drag coefficient ∆ C ∈ R from an optimized configuration D T T T T ∆ C = C x + C u + C u + x C x + u C u + u C u + 2 u C u . (3) D D D n D m D n D n m D m n D m x un um 2 2 2 2 x u nn u mm u nm Note that ∆ C is a positive quantity as defined in this study. An estimate of the gust vector ˆ w is computed from D accelerometer measurements as described in Section III of [ 7 ]. Accelerometers L14 , L19 , R14 , and R19 are used for gust estimation and as output signals for the observer.

B. Nominal Control Design The nominal controller u is designed to provide 6-DoF tracking of commands. It utilizes the elevator, rudder, and n ailerons RAILIN and LAILIN as illustrated in Fig. 1. Gains are calculated using only the rigid body dynamics of the aircraft. The structure of the controller is summarized here, while a detailed description can be found in [9].

The nominal controller consists of three components: servomechanism LQR for roll control, servomechanism LQR for flight path control, and a washout filter for yaw. A flight path angle command γ and a heading angle command ψ c c serve as the reference r ∈ R . All necessary error and filter states are appended to the plant state and are collectively denoted as x ∈ R . The final structure of the controller can be summarized as i         ~  δ   K 0 K   x   K 0  elevator ex ei r er            γ          c ~ ~   u =  δ  =  K 0 K   x  +  0 K  (4) n aileron ax ai e ar             ψ c         ~           δ K 0 K x 0 0 rudder r x ri i         or expressed in a more compact form as u = K x + K r . (5) n nx nr A block diagram of the nominal controller is provided in Fig. 2.

The design of the multi-objective controller requires the nominal loop to be closed. Equation (5) is used to obtain the nominal-loop-closed expressions for the plant dynamics and performance metrics presented in Section III.A. The plant dynamics become ¯ ¯ ˙ x = Ax + B r + B u + E w (6) i m m ¯ A = A + B K (7) n nx ¯ B = B + B K . (8) i i n nr Fig. 2 6-DoF nominal controller block diagram Wing root bending moment becomes ¯ ¯ M = M x + M r + M u + M w (9) y x r u m w m ¯ M = M + M K (10) x x u nx n ¯ M = M K . (11) r u nr n and the change in drag coefficient becomes T T T T T T ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ∆ C = C x + C r + C u + x C x + u C u + x C u + r C u + r C x + r C r (12) D D D D m D D m D m D m D D x r um 2 m 2 x um r um r x 2 x u mm r ¯ C = C + C K (13) D D D nx x x un ¯ C = C K (14) D D nr r un ¯ C = C (15) D D um um T ¯ C = C + K C K (16) D D D nx 2 2 nx 2 x x u nn ¯ C = C (17) D D 2 2 u mm u mm T ¯ C = 2 K C (18) D D x um nx 2 u nm T ¯ C = 2 K C (19) D D r um nr 2 u nm T ¯ C = 2 K C K (20) D D nx r x nr 2 u nn T ¯ C = K C K . (21) D D nr 2 nr 2 r u nn C. Multi-objective Control Design The multi-objective controller u is used to provide flexible mode suppression, load alleviation, and drag reduction m through the use of the optimal solution to a carefully weighted cost function. It uses the elevator and ailerons ROBFMPF1ROBFMPF8 , RIBFMPF1RIBFMPF3 , LOBFMPF1LOBFMPF8 , and LIBFMPF1:LIBFMPF3 . The infinite time horizon cost function balancing these performance objectives is given by ∫ t f ( ) T T 2 J = lim x Qx + u Ru + q M + q ∆ C dt . (22) m M D D m y t →∞ f where the first term in is intended to provide flexible motion stabilization, the third load alleviation, and the fourth drag reduction. The matrix Q ≥ 0 is chosen to be a diagonal matrix whose only nonzero entries are those that correspond to the modal displacement sates in the x portion of the plant’s state vector, or simply e   ~ ~ ~  0 0 0      ~ ~ Q =   . (23) 0 q I 0 e     ~ ~ ~  0 0 0    The matrix R > 0 as well as scalars q >0, q ≥ 0 , q ≥ 0 are weights available for the user to assign.

e M D Proceeding to solve the optimal control problem for u , the Hamiltonian is given by m ( ) 1 1 1 2 1 ( T T T ¯ ¯ ¯ ¯ ¯ ¯ H = x Qx + u Ru + q M x + M r + M u + M w + q C x + C r + C u + x C x m M x r u m w D D D D m D m m x r um 2 x 2 2 2 2 [ ] ) T T T T T T ¯ ¯ ¯ ¯ ¯ ¯ ¯ + u C u + x C u + r C u + r C x + r C r + λ Ax + B r + B u + E w . (24) D m D m D m D D i m m m 2 x um r um r x 2 u mm r The necessary conditions of optimality are T ( ) ( ) − ∂ H 1 T T T T ˙ ¯ ¯ ¯ ¯ ¯ ¯ ¯ λ = = − Qx − q M M x + M r + M u + M w − q C + 2 ¯ C x + C r + C u − A λ M x r u m w D D D m x m D 2 D x um x r x x ∂ x 2 (25) and T ( ) ( ) ∂ H 1 T T T T T ¯ ¯ ¯ ¯ ¯ = Ru + q M M x + M r + M u + M w + q C + C x + C r + 2 ¯ C u + B λ = 0 .

m M x r u m w D D m u m D D D 2 m m um x um r um u mm ∂ u 2 m (26) Taking the adjoint to be of the form λ = W x + V r + V w + V and combining with the optimality conditions results in r w 0 the optimal control solution u = K x + K r + K w + Λ (27) m x r w The gains are defined as ( ) − 1 T T T ¯ ¯ ¯ K = − R q C + q M M + B W (28) x D M x u m D m x um ( ) − 1 T T T ¯ ¯ ¯ K = − R q C + q M M + B V (29) r D M r r u m D m r um ( ) − 1 T T ¯ K = − R q M M + B V (30) w M w w u m m ( ) − 1 T T ¯ ¯ Λ = − R q C + B V (31) D 0 D m um Additional definitions are given by ( ) [ 1 1 − 1 T T − 1 T T ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ V = − V W B + q C + q M M − W B R q C + q M M r i D M r m D M r D x D u r x r um m 2 2 ( ) ( ) ] 1 1 T − 1 T T ¯ ¯ ¯ ¯ ¯ − q C + q M M R q C + q M M (32) D D M u D M r x um x m D u r um m 2 2 ( ) [ ( ) ( ) ] − 1 T − 1 T T − 1 T ¯ ¯ ¯ ¯ ¯ ¯ ¯ V = − V W E + q M M − W B R q M M − q C + q M M R q M M (33) w M w m M w D D M u M w x u x um x m u m m ( ) ( ) ( ) [ ] 1 1 1 1 ¯ − 1 T − 1 T T − 1 T ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ V = − V q C − W B R q C − q C + q M M R q C (34) 0 D m D D D M u D D D x um x m D x um um 2 2 2 2 T ¯ ¯ R = R + q C + q M M (35) D D M u u 2 m m u mm ¯ ¯ T − 1 ¯ ¯ ¯ V = A − W B R B (36) m m

where W is the solution to the Riccati equation

¯ ¯ T − 1 ¯

˙ ¯ ¯ ¯ ¯

W + W A + A W − W B R B W + Q = 0 (37)

m m

with

− 1 T − 1 T

¯

¯ ¯ ¯ ¯ ¯ ¯

A = A − q B R C − q B R M M (38)

D m M m x D u x um m

( ) ( )

1 1

¯ T T − 1 T T

¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯

Q = Q + q C + q M M − q C + q M M R q C + q M M . (39)

D D M x D D M u D M x x x um x m u 2 D m x x um

2 2

Typically a differential Riccati equation such as the one in Eq. (51) requires a backwards-in-time solution procedure to

˙

solve for W ( t ) . Here an infinite time horizon is assumed so that W = 0 and a constant W can be obtained from the algebraic Riccati equation.

IV. Weights Vary as Function of Time

Since weight selection is a significant factor in the performance of the multi-objective controller, ways to expand the permissible type of weights are explored next. The most straightforward extension is to replace the constant weights with weights that can vary as a function of time. Variable weights could then be used to shift the relative importance of the performance tasks in the cost function during operation, especially as the aircraft moves through the flight envelope.

The variable weight concept is demonstrated in a time-varying simulation of the CRM as it moves from one full fuel trim point at Mach 0.8 to another at Mach 0.9. The control system is tasked with simultaneously rejecting a discrete disturbance and pursuing the previously discussed performance tasks. Performance of the multi-objective controller with the same constant weights throughout the simulation is compared to the performance with the weights changing according to a predefined function of time. The time-varying simulation is executed by creating time-varying plant matrices from a convex combination of the existing constant plant matrices available at the beginning and ending trim points, e.g.

t

A ( t ) = ( 1 − c ( t )) A + c ( t ) A where c ( t ) = for 0 ≤ t ≤ t . (40)

1 2 f

t

f The nominal control and observer gains are computed using the plant models at the beginning trim point. They are not changed during the simulation. The plant matrices used in the multi-objective gain calculation are also held constant at their beginning trim point values. The Q and q weights are made time-varying using the same approach M

q ( t ) = ( 1 − c ( t )) q + c ( t ) q q ( t ) = ( 1 − c ( t )) q + c ( t ) q . (41)

e e 1 e 2 M M 1 M 2 It is important to note that the time dependence of the weights causes the Riccati equation that arises in the multi-objective

¯ ¯

¯ ¯ ˙

gain calculation, shown in Eq. (51) , to include time dependency in A and Q . The previous W = 0 justification no longer applies. Instead, the weights can be frozen and the algebraic Riccati equation resolved at each time step to obtain new multi-objective gains in the absence of an explicit closed-form solution.

◦ The simulation results of the aircraft attempting to maintain level flight ( γ = 0 ) when a five second gust disturbance is applied are shown in Figs. 4–8. The components of the gust disturbance vector w are shown in Fig. 3. Note that the same disturbance is used for all simulations in this paper. The load alleviation weight is decreased from q = 1 e − 10 M 1 to q = 1 e − 13 while the flexible motion suppression weight is increased from q = 1 e − 3 to q = 1 e − 2 since M 2 e 1 e 2 the changing flight conditions are moving towards the flutter boundary. Drag reduction is not considered and q = 0 .

D From the figures it is clear that maintaining constant weights during the time-varying simulation results in instability.

However, use of the time-varying weights recovers stable performance. The elevator and flap deflections shown in Fig. 7 and Fig. 8 remain within limits for the varying weights case.

V. Weights Vary as Function of State

Next consider the case where the weights vary as a function of the state. Similar to the function of time case, many previously constant parameters become state dependent. Additional terms are also introduced in the optimality conditions due to use of the partial derivative with respect to the state, leading to inconvenient or intractable closed-form solutions to the optimal control problem. Several suboptimal approximation methods omitting consideration of these terms are suggested here and their functionality demonstrated in simulation.

10 [ft/s] w 0 5 10 15 5 ] 0 [ft/s w -5 -10 -15 0 5 10 15 time [s]

Fig. 3 Discrete gust disturbance used in all simulations

0.5 command fixed weights time-varying weights 0.4 0.3 0.2 [deg] 0.1 -0.1 0 5 10 15 20 25 time [s]

Fig. 4 Gust-disturbed flight path angle tracking for fixed vs. time-varying weights

State-based variation of weights permits a range of functionality not possible with time-based weight variation.

Smooth transitions between objectives can be executed without having to preplan the switch and schedule it as a function of time. For example, if it is desired that the cost function include only load alleviation below a certain altitude and then transition to include only drag reduction above that altitude then the following pair of weight variation functions permits use of the single unifying cost function in Eq. (22):

[ ]

q ( x ) = q 1 − g ( x ) (42)

M M

q ( x ) = q g ( x ) (43)

D D

The subfunction g ( x ) is given by



 0 if h ≤ h

lb







h − h l b



g ( x ) = if h < h < h (44)

lb ub h − h

 u b l b







1 if h ≥ h

ub



where h is the aircraft’s altitude as taken from the state, h is the lower bound on the altitude region for the objective lb transition, and h is the upper bound. Weights can also be made to change based on sensor measurements or estimated ub state values. The result is a scheme that provides sensor-driven prioritization of objectives which will be used to demonstrate the techniques discussed in this section.

6 50 fixed weights fixed weights time-varying weights time-varying weights 3 10 || [lb-ft] y ||M mode 1 displacement -10 -20 0 -30 0 5 10 15 20 25 0 5 10 15 20 25 time [s] time [s] Fig. 5 Wing root bending moment for fixed Fig. 6 First mode’s displacement for fixed vs.

vs. time-varying weights time-varying weights

fixed weights 0.6 time-varying weights 1.5 0.4 0.5 0.2 flap deflection [deg] -0.5 elevator deflection [deg] -0.2 -1 -1.5 -2 5 0 5 10 15 20 25 10 time [s] flap number time [s] Fig. 7 Elevator deflection for fixed vs.

Fig. 8 Multi-objective controller flap

time-varying weights

deflection for time-varying weights

A. State Dependent Riccati Equation

The State Dependent Riccati Equation (SDRE) is one well-established method for handling the case of state dependency in the optimal control setting. Here the choice is made to replace Q and q with Q ( x ) and q ( x ) in the M M cost function. State dependency is introduced into the Riccati equation in Eq. (51), now known as an SDRE,

¯ ¯ T − 1 T ¯

¯ ¯ ¯ ¯

W ( x ) A ( x ) + A ( x ) W ( x ) − W ( x ) B R B W ( x ) + Q ( x ) = 0 (45)

m m and finding an explicit solution for W ( x ) may not be possible. Further, the optimality conditions will change and the choice of u in Eq. (27) will not be optimal unless an additional SDRE optimality condition is satisfied. However, an m adequate suboptimal approximation may be to repeatedly freeze the state-dependent parameters and solve the SDRE at a sufficiently high rate. Further information regarding the SDRE and its optimality criteria can be found in [10, 11].

To illustrate use of the SDRE, consider a time invariant simulation is conduced at a single trim point to demonstrate that use of variable weights can simply improve performance of the multi-objective controller. Again the five second gust disturbance is applied while the control system attempts to maintain level flight. Select state-dependent weights of

the form

[ ]

q ( x ) = η ‖ Sx ‖ + 1 q (46)

e e 0

[ ]

− 1

q ( x ) = η ‖ Sx ‖ + 1 q . (47)

M M 0 The matrix S is used to select the displacement of the first mode out of the full state vector such that the flexible motion suppression weight increases and the load alleviation weight decreases with larger displacements of the mode. The constant components of the weights are selected as q = 1 e − 10 and q = 1 e − 1 and the scaling constant η = 1 e − 2 .

M 0 e 0 Drag reduction is omitted by setting q = 0 . Results comparing use of the time-varying,state dependent weights with D fixed weights are presented in Figs. 9–13. Performance with nominal control only (multi-objective control u off) are m also included for reference. Figure 9 demonstrates how the flight path angle tracking is altered now that the control design is subject to pursuit of the other performance tasks. Figure 10 shows significant improvement in peak load alleviation with use of the time-varying weights. Similarly, Fig. 11 indicates significant peak reduction of the first mode’s displacement with varying weights. The elevator defelection for all controllers is shown in Fig. 12 while flap deflections for the varying weights controller only are shown in Fig. 13. The motion of all control surfaces remains within reasonable limits.

0.45 command u off 0.4 m tv weights fixed weights 0.35 0.3 0.25 0.2 [deg] 0.15 0.1 0.05 -0.05 0 5 10 15 time [s]

Fig. 9 Gust-disturbed flight path angle tracking for fixed vs. time-varying weights

Resolving the Riccati equation and forming new control gains at every step can become computationally burdensome as the dimension of the plant model increases. One simple alternative is to reduce the frequency of determining new gains to every few time steps. For comparison, the previous simulation is performed again with gains recomputed every 0.2 s. The wing root bending moment and flap deflections are shown in Figs. 14 and 15. Performance is clearly degraded with transients from the intermittent control gain recompilation disrupting the previously smooth control action, though the closed loop system remains stable. Other computationally tractable alternatives to the repeated Riccati solution are presented subsequently.

6 50 u off u off m m tv weights tv weights fixed weights fixed weights 3 0 || [lb-ft] y ||M -10 mode 1 displacement -20 -30 -40 0 -50 0 5 10 15 0 5 10 15 time [s] time [s] Fig. 10 Wing root bending moment for fixed Fig. 11 First mode’s displacement for fixed vs.

vs. time-varying weights time-varying weights

u off m tv weights fixed weights 0.5 -0.5 elevator deflection [deg] -1 -1.5 0 5 10 15 time [s] Fig. 12 Elevator deflection for fixed vs.

Fig. 13 Multi-objective controller flap

time-varying weights

deflection for time-varying weights

2.5 1.5 || [lb-ft] y ||M 0.5 0 5 10 15 time [s]

Fig. 15 Multi-objective controller flap

Fig. 14 Wing root bending moment for

deflection for reduced frequency gain

reduced frequency gain computation

computation

B. Polynomial Approximation

While repeatedly solving the SDRE as was done in the previous section can provide good results, it is a time consuming operation. A single expression for the the Riccati solution would be more viable for implementation even if it only approximates the true solution and results in a suboptimal controller. One such reasonable approximation could be a polynomial expansion for the Riccati solution when the mode suppression weight Q is also restricted to a

polynomial structure. For example, say Q is given by

Q ( x ) = Q f ( x ) + Q f ( x ) + Q (48)

2 1 0 with each of the constant Q matrices selected by the user. Then assume that W also has the form i

W ( x ) = W f ( x ) + W f ( x ) + W . (49)

2 1 0 For this structure it is possible to determine approximate expressions for the constant W matrices that can be solved i ahead of time. The resulting expression for W in Eq. (49) will not be exact and results from dropping higher order terms of f ( x ) . Thus it is necessary to select f ( x ) as a function of the state that decays. A study of the error involved in using the polynomial approximation and further investigation of its performance can be found in [12].

In the case where Q is the only state dependent weight the multi-objective cost function is given by

∫

t f

( )

T T 2

J = lim x Q ( x ) x + u Ru + q M + q ∆ C dt (50)

m m M y D D t →∞

f 2

where Q ( x ) is given by Eq. (48) . The W that appears in the multi-objective gains is assumed to have the form given in Eq. (49) . It is straightforward to obtain expressions for the W matrices by substituting the assumed polynomial i

expressions into the algebraic Riccati equation used to solve for W

¯ ¯ T − 1 T ¯

¯ ¯ ¯ ¯

W A + A W − W B R B W + Q = 0 . (51)

m m

¯

¯ ¯

The definitions for R in Eq. (35) and A in Eq. (38) demonstrate that they have no state dependence and thus remain

¯

¯

constant matrices. The definition for Q in Eq. (39) indicates state dependence through the presence of Q and can be

rewritten as

¯ 2

¯ ¯

Q ( x ) = Q f ( x ) + Q f ( x ) + Q (52)

2 1 0

¯

where Q = Q + T and

0 0

( ) ( )

1 1

T T − 1 T T

¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯

T = q C + q M M − q C + q M M R q C + q M M (53)

D D M x D D M u D M x 2 x x um x m D u x um m x

2 2

with the user-selected constant matrices Q , Q , Q > 0 . The Riccati equation becomes

2 1 0

[ ] [ ]

2 ¯ ¯ T 2

¯ ¯

W f ( x ) + W f ( x ) + W A + A W f ( x ) + W f ( x ) + W

2 1 0 2 1 0

[ ] [ ] [ ]

2 − 1 T 2 2

¯ ¯

− W f ( x ) + W f ( x ) + W B R B W f ( x ) + W f ( x ) + W + Q f ( x ) + Q f ( x ) + Q = 0 (54) 2 1 0 m 2 1 0 2 1 0 m where terms higher than f ( x ) have been dropped. Grouping constant terms gives the equation

¯ ¯ T − 1 T

¯ ¯ ¯ ¯

W A + A W − W B R B W + Q = 0 (55)

0 0 0 m 0 0 m which is a Riccati equation that can be solved offline for W . Grouping f ( x ) terms gives

( ) ( )

T

¯ − 1 ¯ − 1

¯ ¯ ¯ ¯

W A − B R B W + A − B R B W W + Q = 0 (56)

1 m m 0 m m 0 1 1 which is a Lyapunov equation that can be solved offline for W since W is now available. Similarly, the f ( x ) Lyapunov 1 0

equation

( ) ( ) ( )

T

¯ − 1 T ¯ − 1 T − 1 T

¯ ¯ ¯ ¯ ¯

W A − B R B W + A − B R B W W + Q − W B R B W = 0 (57)

2 m 0 m 0 2 2 1 m 1 m m m can be solved offline for W . Thus, W ( x ) can be constructed online without repeated online Riccati solutions.

Remark. The polynomial approximation cannot be easily implemented with every weight in the multi-objective cost function. Note that if q or q were to be polynomial expressions of the state the rewritten Riccati equation would M D − 1

¯

become more complex due to their involvement in R .

The polynomial approximation is used to create an expression for W during the same half second gust disturbance simulation used in the previous section. Here f ( x ) = ‖ Sx ‖ + 1 where again S selects the displacement of the first mode.

The polynomial terms are chosen using q = 1 , q = 2 , and q = 1 e − 3 to construct Q , Q , and Q according to e 0 e 1 e 2 0 1 2 Eq. (23) . The results are shown in Figs. 16–20 where performance using the time-varying controller based on the polynomial approximation is compared to performance when the Riccati equation is continually resolved. Performance without use of the multi-objective controller is also shown for reference. Note that the polynomial approximation controller does not fully match the Riccati resolve controller but performance trends are similar. Both provide notable peak reduction of wing root bending moment and displacement of the first mode as seen in Figs. 17 and 18 without significantly disrupting flight path angle tracking performance shown in Fig. 16. Control surface deflections are again reasonable as seen in Figs. 19 and 20.

0.5 command u off m polynomial approx 0.4 resolve 0.3 0.2 [deg] 0.1 -0.1 0 5 10 15 time [s]

Fig. 16 Gust-disturbed flight path angle tracking for fixed vs. time-varying weights

C. Quartic Optimal Control

For the polynomial approach discussed in the previous section it was noted that selecting a decaying f ( x ) would result in a better approximation to the true state dependent Riccati solution. However, when the multi-objective control 6 50 u off u off m m polynomial approx polynomial approx resolve resolve 3 0 || [lb-ft] y ||M -10 mode 1 displacement -20 -30 -40 0 -50 0 5 10 15 0 5 10 15 time [s] time [s] Fig. 17 Wing root bending moment for fixed Fig. 18 First mode’s displacement for fixed vs.

vs. time-varying weights time-varying weights

1.5 u off m polynomial approx resolve 0.5 -0.5 elevator deflection [deg] -1 -1.5 -2 0 5 10 15 time [s] Fig. 19 Elevator deflection for fixed vs.

Fig. 20 Multi-objective controller flap

time-varying weights

deflection for time-varying weights

gains can be bounded there exists a specific case of the polynomial structure that permits formal stability analysis while still providing simple, offline means of calculating the polynomial coefficient matrices.

Consider the multi-objective cost function where again only Q is permitted to vary as given in Eq. (50) with u as m

derived in Section III.C. Q is explicitly restricted to the form

Q ( x ) = Q ‖ x ‖ + Q (58)

2 0 where the constant weights Q , Q > 0 are user choice. The plant dynamics using this Q ( x ) with nominal and 2 0 multi-objective control loops closed can be shown to be bounded with respect to the reference dynamics

¯ ¯

˙ x = Ax + B r (59)

m m i that correspond to the undisturbed plant dynamics with nominal control loop closed and no multi-objective control term.

This optimal control problem is referred to as quartic optimal control due to the structure of the first term in the cost

function

∫

t f

[ ( ) ]

T 2 T 2

J = lim x Q ‖ x ‖ + Q x + u Ru + q M + q ∆ C dt . (60)

2 0 m M D D m y t →∞

f 2

The Riccati solution W is replaced with the similar polynomial

W ( x ) = W ‖ x ‖ + W (61)

2 0 and equations providing W and W can be determined through Lyapunov analysis as summarized in Theorem 1.

0 2 Theorem 1. Consider stable state dynamics of the form in Eq. (6) subject to a bounded disturbance w with bounded command signal r . Inclusion of the multi-objective controller in Eq. (27) when using state dependent weights of the form provided in Eq. (58) results in bounded tracking of the undisturbed, nominal loop closed dynamics when

W ( x ) = W ‖ x ‖ + W (62)

2 0 is used in place of W in the control gains. Additionally, W can be obtained from the Riccati equation

¯ ¯ T − 1 T

¯ ¯ ¯

W A + A W − W B R B W = − Q (63)

0 0 0 m 0 0 m

and W from the Lyapunov equation

( ) ( )

T − 1 T − 1 T

¯ ¯

¯ ¯ ¯ ¯

W A − 2 B R B W + A − 2 B R B W W = − Q (64)

2 m 0 m 0 2 3 m m where Q , Q > 0 .

0 3

Proof. Defining the error as e = x − x , the error dynamics are given by

m

¯

˙ e = Ae + B u + E w . (65)

m m

Express u in terms of e

m

u = K e + K x + K r + K w + Λ (66)

m x x m r w

and substitute it along with the K definition into the error dynamics

x

[ ( ) ]

− 1 T T − 1 T

¯ ¯ ¯ ¯ ¯

˙ e = A − B R q C + q M M − B R B W ( x ) e + B ( K r + K w + Λ + K x ) + E w . (67)

m D M x m m r w x m D u m x um m

This is more compactly expressed as

[ ]

¯ − 1 T

¯ ¯

˙ e = A − B R B W ( x ) e + η (68)

m m where η = B ( K r + K w + Λ + K x ) + E w . Similarly, including the K definition in u and using the resulting m r w x m x m

expression in ˙ x leads to

[ ]

¯ − 1 T

¯ ¯

˙ x = A − B R B W ( x ) x + ˜ η (69)

m m

( )

¯

where ˜ η = B + B K r + ( E + B K ) w + B Λ .

i m r m w m

Select the Lyapunov function

T

V ( e ) = e W ( x ) e (70)

and take its derivative

( )

T T T T T

˙

V = ˙ e W ( x ) e + e W ( x ) ˙ e + e W ˙ x x + W x ˙ x e . (71)

2 2 Substituting the error and state dynamics and manipulating scalar quantities where possible leads to

[ ] ( ) [ ( ) ]

T ¯ T ¯ − 1 T T T ¯ T ¯

˙ ¯ ¯ ¯ ¯ ¯

V = e A W ( x ) + W ( x ) A − 2 W ( x ) B R B W ( x ) e + e W e x A + A x

m 2 m

( ) [ ] ( )

T T − 1 T T T T

¯

− 2 e W e x B R B W ( x ) x + 2 e W ( x ) η + 2 e W e x ˜ η (72)

2 m 2 m Next, allow Q and Q to be the positive definite solutions of the algebraic Riccati equations 0 2

¯ ¯ T − 1 T

¯ ¯ ¯

W A + A W − W B R B W = − Q (73)

0 0 0 m 0 0 m

( ) ( )

T − 1 T − 1 T − 1 T

¯ ¯

¯ ¯ ¯ ¯ ¯

W A − 2 B R B W + A − 2 B R B W W − W B R B W = − Q . (74)

2 m 0 m 0 2 2 m 2 2 m m m

and include them in Eq. (72) to reach

[ ( ) ]

T − 1 T 2 − 1 2 2

˙ ¯ ¯

V = e − Q − W B R B W − Q ‖ x ‖ + W B R B W ‖ x ‖ 1 − 2 ‖ x ‖ e

0 0 m 0 2 2 m m 2 m

( ) [ ( ) ] ( ) [ ( ) ]

T T ¯ T ¯ T T − 1 T 2

¯ ¯ ¯

+ e W e x A + A x − 2 e W e x B R B W + W ‖ x ‖ x

2 2 m 0 2 m

( ) ( )

T 2 T T

+ 2 e W + W ‖ x ‖ η + 2 e W e x ˜ η. (75)

0 2 2

Selecting

− 1 T 2

¯

Q = Q + W B R B W ‖ x ‖ (76)

2 3 2 m 2 m

with Q > 0 leads to

[ ] ( ) [ ( ) ]

T − 1 T 2 − 1 4 T T ¯ T ¯

˙ ¯ ¯ ¯ ¯

V = e − Q − W B R B W − Q ‖ x ‖ − 2 W B R B W ‖ x ‖ e + e W e x A + A x

0 0 m 0 3 2 m m 2 2 m

( ) [ ( ) ] ( ) ( )

T T − 1 T 2 T 2 T T

¯

− 2 e W e x B R B W + W ‖ x ‖ x + 2 e W + W ‖ x ‖ η + 2 e W e x ˜ η. (77)

2 m 0 2 0 2 2 m Including this choice of Q in Eq. (76) also means that W can be found from the Lyapunov equation 2 2

( ) ( )

T − 1 T − 1 T

¯ ¯

¯ ¯ ¯ ¯

W A − 2 B R B W + A − 2 B R B W W = − Q (78)

2 m 0 m 0 2 3 m m instead of the Riccati equation. Define the constants c = λ ( Q ) , c = λ ( Q ) , c = λ ( W ) , c = λ ( W ) , 1 min 0 2 min 3 3 min 0 4 min 2

¯ ¯ ¯

− 1 T T

¯ ¯ ¯ ¯

and c = λ ( B R B ) . Also note that for appropriately chosen weights A is Hurwitz such that A + A ≤ 0 and 5 min m m

( ) ( )

¯ T ¯ − 1 T 2

¯ ¯ ¯

λ A + A = 0 serves as an upper bound. Additionally, B R B W + W ‖ x ‖ > 0 . The derivative can then be min m 1 2 m

bounded as

[ ( ) ] ( )

2 2 2 4 2 2 2

˙

V ≤ − c + c c + ( c + 2 c c c ) ‖ x ‖ + 2 c c + c ‖ x ‖ ‖ e ‖ + 2 c + c ‖ x ‖ ‖ e ‖‖ η ‖ + 2 c ‖ e ‖ ‖ x ‖‖ ˜ η ‖ (79) 1 5 2 4 5 3 5 4 3 4 4 3 4 such that the error dynamics are asymptotically stable when η = 0 and ˜ η = 0 as would be the case when r , w , and q D are zeroed.

˙

When η , 0 and ˜ η , 0 , the last two terms in V can be bounded by noting

2 2 2 2

‖ x ‖ = ‖ e + x ‖ ≤ ‖ e ‖ + 2 ‖ e ‖‖ x ‖ + ‖ x ‖ (80)

m m m 2 2

and recalling that ‖ a ‖‖ b ‖ ≤ ( ‖ a ‖ + ‖ b ‖ ) to reach

( ) ( )

2 2 2 2 2 2 2 2 c + c ‖ x ‖ ‖ e ‖‖ η ‖ + 2 c ‖ e ‖ ‖ x ‖‖ ˜ η ‖ ≤ 2 c ‖ e ‖ ‖ x ‖ + c ‖ ˜ η ‖ + ‖ η ‖ ‖ e ‖ 3 4 4 4 4

( )

2 2 2

+ 2 c ‖ x ‖‖ η ‖ + c ‖ η ‖ ‖ e ‖ + c ‖ x ‖ ‖ η ‖ (81)

4 m 3 4 m

2 2 ˙ Using the − c c ‖ e ‖ term to complete the square with the ‖ e ‖ term in the previous expression leads to the V bound [ ( ) ] c ‖ x ‖‖ η ‖ + c ‖ η ‖ 4 m 3 2 2 2 2 ˙ * + V ≤ − c − c ‖ ˜ η ‖ + ‖ η ‖ ‖ e ‖ − c c ‖ e ‖ − 1 4 5 c c , 3 - ( ) 2 2 2 2 2 4 2 − [ c − 2 c ] ‖ e ‖ ‖ x ‖ − 2 c c c ‖ e ‖ ‖ x ‖ − 2 c c + c ‖ x ‖ ‖ e ‖ + c (82) 2 4 4 5 3 5 4 6 where the constant ( ) c ‖ x ‖‖ η ‖ + c ‖ η ‖ 4 m 3 2 2 c ≥ c ‖ x ‖ ‖ η ‖ + (83) 6 4 m c c has been defined since η and ˜ η will be bounded quantities so long as the applied control gains are bounded, the reference ( ) 2 2 dynamics are stable, and r and w are bounded. Then, selecting weights such that c > 2 c and c > c ‖ η ‖ + ‖ ˜ η ‖ 2 4 1 4 ˙ means the only nonnegative term remaining in V is c . The tracking error will therefore converge to and remain in a set whose size is related to this quantity.

VI. Conclusion

This paper presents several approaches to addressing time or state-based variation of weights in an optimal control problem for multi-objective performance improvement of a flexible wing aircraft. The optimal gains balancing modal suppression, load alleviation, and drag reduction were determined as functions of the user-selected weights in the cost function. Some of the weights were then permitted to change online. Variation according to a prescribed function of time was shown to better accommodate changing flight conditions. Variation according to selected state values demonstrated improved gust rejection in level flight. A polynomial approximation to a necessary Riccati solution was proposed as a way to avoid intensive online gain calculations, and a special form of the polynomial approximation was shown to provide bounded tracking performance using Lyapunov stability analysis.

In the future existence of other closed form versions of the time-varying multi-objective controller will be investigated.

Other performance objectives have been proposed for combination with the existing performance terms in the cost function. A version of the controller with disturbance-based weight variation has also been proposed. The framework of the quartic optimal controller will be used to explore stability of more general weight expressions. Additionally, a thorough analysis of the state functions driving state-dependent weight variation must also be performed.

Acknowledgments

This work was supported by the Advanced Air Transport Technology Project of the NASA Aeronautics Research Mission Directorate. Additional funding was provided by an appointment to the NASA Postdoctoral Program, administered by Universities Space Research Association.

References

[1] Wernli, A. and Cook, G., “Suboptimal control for the nonlinear quadratic regulator problem," Automatica , Vol. 11, No. 1, pp.

75–84, 1975.

[2] Krener, A.J., Aguilar, C.O., and Hunt, T.W., “Series Solutions of HJB Equations," Festschrift in Honor of Uwe Helmke , ISBN 978-1470044008, 2013.

[3] Radhoush, S., Samavat, M., and Vali, M.A., “Optimal control of linear time-varying systems using the Chebyshev wavelets", Systems Science and Control Engineering , Vol. 2, No. 1, 2014, pp. 691–698.

[4] Sinha, S.C. and Joseph, P., “Control, of General Dynamic Systems with Periodically Varying Parameters via Lyapunov-Floquet Transformation", Journal of Dynamics, Measurement, and Control , Vol. 116, No. 4, 1994, pp. 650–658.

[5] Kanarachos, A.E. and Geramanis, K.T., “Optimal Control of Time-Varying Dynamic Systems", Progress in System and Robot Analysis and Control Design , Vol. 243, Springer, London, 1999, pp. 127–136.

[6] Hashemi, K.E., Nguyen, N.T., Drew, M., Chaparro, D., and Ting, E., “Performance Optimizing Gust Load Alleviation Control of Flexible Wing Aircraft,” AIAA SciTech, Kissimmee, FL, January 2018.

[7] Nguyen, N., Hashemi, K.E., and Drew, M., “Multi-Objective Adaptive Control for Load Alleviation and Drag Minimization of Flexible Aircraft," AIAA SciTech, Kissimmee, FL, January 2018.

[8] Vassberg, J.C. and Rivers, M.S.,“Development of a Common Research Model for Applied CFD Validation Studies," 26th AIAA Applied Aerodynamics Conference, Honolulu, HI, August 2008.

[9] Drew, M., Hashemi, K.E., Cramer, N., and Nguyen, N.T.,“Multi-objective Optimal Control of the 6-DoF Aeroservoelastic Common Research Model with Aspect Ratio 13.5 Wing," AIAA SciTech, San Diego, CA, January 2019 (submitted).

[10] Cloutier, J.R., D’Souza, C.N., and Mracek, C.P., “Nonlinear Regulation and Nonlinear H ∞ Control Via the State-Dependent Riccati Equation Technique: Part 1, Theory," International Conference on Nonlinear Problems in Aviation and Aerospace, Daytona Beach, FL, May 1996.

[11] Cloutier, J.R., “State-dependent Riccati Equation Techniques," American Control Conference, Albuquerque, NM, June 1997.

th [12] Lau, E. and Krener, A.J., “LPV Control of Two Dimensional Wing Flutter," 38 IEEE CDC, Phoenix, AZ, December 1999.