Boundary Stabilization of a Bending and Twisting Wing by Linear Quadratic Gaussian Theory

TL;DR

This paper develops a boundary stabilization method for a cantilever wing's bending and twisting using an infinite-dimensional LQR approach combined with a Kalman filter, incorporating aerodynamic forces based on Wagner's function.

Contribution

It introduces a novel boundary control and estimation framework for wing stabilization using infinite-dimensional LQR and Kalman filtering with aerodynamic modeling.

Findings

Effective boundary stabilization achieved for bending and twisting modes.

Integration of aerodynamic forces into the control model.

Demonstrated feasibility of the approach for high aspect ratio wings.

Abstract

We first consider the stabilization of the bending and twisting of a rectangular cantilever beam of moderate to high aspect ratio using full state feedback boundary control. Our approach is an infinite dimesnional extension of Linear Quadratic Regulation (LQR). The we develop an infinite dimensional Kalman filter that processes two point measurements and returns an estimate of the full state. The Linear Quadartuc Gaussian approach is to use this estimate in the place of the full state in the LQR feedback Then we add aerodynamic forces to obtain a model of a wing. The aerodyamic model is based on a two dimensional state space approximate realiztion of Wagner's indicial function by R.~T.~Jones.

Figures

Equations

$$\displaystyle\left[\begin{array}[]{ccccccccc}\mu&-S_{y}\\ -S_{y}&I_{y}\end{array}\right]\left[\begin{array}[]{ccccccccc}\frac{\partial^{2}w}{\partial t^{2}}(y,t)\\ \frac{\partial^{2}\theta}{\partial t^{2}}(y,t)\end{array}\right]$$

$$\displaystyle\begin{array}[]{llllll}\mu&\mbox{ mass per unit span}\\ EI&\mbox{ bending rigidity}\\ GJ&\mbox{ torsion rigidity}\\ S_{y}&\mbox{ static moment per unit span about $y$ axis}\\ I_{y}&\mbox{ moment of inertia per unit span about $y$ axis}\end{array}$$

$$\frac{{\partial }^{2}w}{\partial y^2}(L,t)=0,$$

$$w(0,t)=0,$$

$$\frac{\partial \theta }{\partial y}(L,t)$$

$$\frac{\partial \theta }{\partial y}(0,t)$$

$$\displaystyle B=\left[\begin{array}[]{ccccccccc}B_{1}&0\\ 0&B_{2}\end{array}\right]$$

$$z(y,t)$$

$$M\frac{\partial z}{\partial t}$$

$$D$$

$$M$$

$$z_1(0,t)=0,$$

$$\frac{{\partial }^2{z}_1}{\partial y^2}(L,t)=0,$$

$$\frac{\partial z_3}{\partial y}(0,t)=B_2{u}_2(t),$$

$$u(t)$$

$${\int }_0^{\infty }{\iint }_{\mathcal{S}}{z}^{\mathrm{T}}(y_1,t)Q(y_1,y_2)z(y_2,t)+u^{\mathrm{T}}(t)Ru(t)dAdt$$

$$\begin{array}[]{rllllllllllllllllllllllllllllllll}P_{2,j}(0,y_{2})=0,&&P_{i,2}(y_{1},0)=0\\ \frac{\partial^{2}P_{2,j}}{\partial y^{2}_{1}}(0,y_{2})=0,&&\frac{\partial^{2}P_{i,2}}{\partial y^{2}_{2}}(y_{1},0)=0\\ \frac{\partial^{2}P_{2,j}}{\partial y^{2}_{1}}(L,y_{2})=0,&&\frac{\partial^{2}P_{i,2}}{\partial y^{2}_{2}}(y_{1},L)=0\\ \frac{\partial^{3}P_{2,j}}{\partial y^{3}_{1}}(L,y_{2})=0,&&\frac{\partial^{3}P_{i,2}}{\partial y^{3}_{2}}(y_{1},L)=0\\ \frac{\partial P_{4,j}}{\partial y_{1}}(0,y_{2})=0,&&\frac{\partial P_{i,4}}{\partial y_{2}}(y_{1},0)=0\\ \frac{\partial P_{4,j}}{\partial y_{1}}(L,y_{2})=0,&&\frac{\partial P_{i,4}}{\partial y_{2}}(y_{1},L)=0\end{array}$$

$$0$$

$${\int }_0^{\infty }{\iint }_{\mathcal{S}}{z}^{\mathrm{T}}(y_1,t)Q(y_1,y_2)z(y_2,t)dAdt$$

$$0={\iint }_{\mathcal{S}}{z}^{\mathrm{T}}(y_1,0)MP(y_1,y_2)Mz(y_2,0)dA$$

$$\displaystyle+\int_{0}^{\infty}\iint_{\cal S}z^{\rm T}(y_{1},t)\left[\begin{array}[]{ccccccccc}-EI\frac{\partial^{4}P_{:,2}}{\partial y_{1}^{4}}(y_{1},y_{2})\\ P_{:,1}(y_{1},y_{2})\\ GJ\frac{\partial^{2}P_{:,4}}{\partial y_{1}^{2}}(y_{1},y_{2})\\ P_{:,3}(y_{1},y_{2})\end{array}\right]^{\rm T}Mz(y_{2},t)\ dA\ dt$$

$$\displaystyle+\int_{0}^{\infty}\iint_{\cal S}z^{\rm T}(y_{1},t)M\left[\begin{array}[]{ccccccccc}-EI\frac{\partial^{4}P_{:,2}}{\partial y_{2}^{4}}(y_{1},y_{2})\\ P_{:,1}(y_{1},y_{2})\\ GJ\frac{\partial^{2}P_{:,4}}{\partial y_{2}^{2}}(y_{1},y_{2})\\ P_{:,3}(y_{1},y_{2})\end{array}\right]z(y_{2},t)\ dA\ dt\$$

$$\displaystyle+\int_{0}^{\infty}\int_{0}^{L}u^{\rm T}(t)B\left[\begin{array}[]{ccccccccc}-EI\frac{\partial P_{:,2}}{\partial y_{1}}(0,y_{2}\\ -GJP_{:,4}(0,y_{2})\end{array}\right]^{\rm T}Mz(y_{2},t)\ dy_{2}\ dt$$

$$\displaystyle+\int_{0}^{\infty}\int_{0}^{L}z^{\rm T}(y_{1},t)M\left[\begin{array}[]{ccccccccc}-EI\frac{\partial P_{:,2}}{\partial y_{2}}(y_{1},0\\ -GJP_{:,4}(y_{1},0)\end{array}\right]Bu(t)\ dy_{1}\ dt$$

$${\iint }_{\mathcal{S}}{\left(u(t)-K(y_1)z(y_1,t)\right)}^{\mathrm{T}}R\left(u(t)-K(y_2)z(y_2,t)\right)dA$$

$$-RK(y_2)$$

$$K(y_2)$$

$$K(y_1)$$

$$\displaystyle\left[\begin{array}[]{ccccccccc}-EI\frac{\partial^{4}P_{2,:}}{\partial y_{1}^{4}}(y_{1},y_{2})\\ P_{1,:}(y_{1},y_{2})\\ GJ\frac{\partial^{2}P_{4,:}}{\partial y_{1}^{2}}(y_{1},y_{2})\\ P_{3,:}(y_{1},y_{2})\end{array}\right]^{\rm T}M$$

$$\displaystyle+M\left[\begin{array}[]{ccccccccc}-EI\frac{\partial^{4}P_{2,:}}{\partial y_{2}^{4}}(y_{1},y_{2})\\ P_{1,:}(y_{1},y_{2})\\ GJ\frac{\partial^{2}P_{4,:}}{\partial y_{2}^{2}}(y_{1},y_{2})\\ P_{3,:}(y_{1},y_{2})\end{array}\right]+Q(y_{1},y_{2})$$

$$\displaystyle=M\left[\begin{array}[]{ccccccccc}EI\frac{\partial P_{2,:}}{\partial y_{2}}(y_{1},0)\\ GJP_{4,:}(y_{1},0)\end{array}\right]^{\rm T}\Gamma\left[\begin{array}[]{ccccccccc}EI\frac{\partial P_{2,:}}{\partial y_{1}}(0,y_{2})\\ GJP_{4,:}(0,y_{2})\end{array}\right]M$$

$$\frac{\partial \theta }{\partial y}(0,t)=0,$$

$$\frac{{\partial }^2\theta }{\partial t^2}$$

$$\displaystyle\left[\begin{array}[]{ccccccccc}\cos{m\pi\over L}y&\pm\sqrt{GJ\ }\ {m\pi i\over L}\cos{m\pi\over L}y\end{array}\right]^{\rm T}$$

$$\varphi (0)=0,\frac{{\partial }^2\varphi }{\partial y^2}(0)=0,\frac{{\partial }^2\varphi }{\partial y^2}(L)=0,\frac{{\partial }^3\varphi }{\partial y^3}(L)=0$$

Full text

I Introduction

High Altitude Long Endurance (HALE) aircraft have highly flexible wings. Helicopter blades are also highly flexible. One would like to use active control to dampen their oscillations. One can model a wing by the beam equation or the plate equation. In [7] it is stated that ”even for relatively low aspect ratios $(4<AR<5)$ of thin rectangular cantilever plates, the model based on the beam theory is closer to (our) experimental model (as) compared to the plate theory.” So we choose a model based on the beam equation.

We assume that the beam can bend and twist and that there are two point actuators, one can deliver a bending moment to the root of the beam and the other can deliver a torque also at the root of the beam. We also assume that there are two point sensors, one measures the vertical velocity of the tip of the beam and other measures the angular velocity of the tip of the beam. We wish to design a compensator that processes the measurements and delivers a bending moment and torque that dampens the bending and torsion oscillations. We design this compensator using Linear Quadratic Gaussian (LQG) theory. After accomplishing this we convert the model of the beam into a model of a wing by adding the aerodynamic forces that the motions of the beam generate. We model the aerodynamic forces by a state space approximate realization of the classical indicial function of Wagner. Wagner’s indicial function is for an airfoil, a cross section of a wing, but we extend it to the whole wing.

Linear Quadratic Gaussian (LQG) theory is a well-known way of designing a dynamic compensator for a controlled and observed finite dimensional dynamical system. We extend LQG to designing a dynamic compensator for a controlled and observed infinite dimensional dynamical system under point actuation and point sensing. The LQG approach breaks the problem into two parts. The first part uses Linear Quadratic Regulator theory to find a stabilizing linear feedback. This feedback assumes that exact measuremnts of the full state are available. In most situations this assumption is unrealistic. The second part is to design a Kalman filter that processes the noisey measurements and returns an estimate of the full state. This estimate is used in place of actual state in the linear feedback. The result is a dynamic compenstor. The closed loop eigevalues of the combined system are the union of the eigenvalues of system under full state feedbck and the eigenvalues of the error dynamics of the Kalman filter. If all these eigevalues lie in the open left plane then the combined plant and Kalman filter is asymptotically stabe.

We start with the model presented in Section 3.6 of the classic treatise of Bisplinghoff, Ashley and Halfman, [2]. This is a linear model which ignores the nonlinear interactions between the bending and the torsion of the rectangular beam. We add point actuation and point sensing to the boundary of the model. Point actuation and point sensing are idealizations of actuation and sensing over small domains on the boundary and the idealization simplifies the mathematical analysis.

The linear model is neutrally stable so the bending and the torsion oscillations do not decay. We seek a compensator to asymptotically stabilize these oscillations. A model for designing a stabilizing compensator does not need to be as accurate as a model for simulation because any errors will decay during the stabilization process.

To the model of [2] we add two actuators located where the horizontal beam is joined to its support. One actuator affects the bending of the beam at its support and other affects the torsion of the beam at its support. We first seek a state feedback control law which drives the oscillations to zero. To find such a feedback we set up and solve a Linear Quadratic Regulator (LQR). We are certainly not the first to use LQR to find a state feedback that stabilizes a bending and twisting beam. Edwards, Breakwell and Bryson [5] used LQR to find a state feedback that stabilizes an airfoil, a cross section of a wing, using leading edge and trailing edge flaps. We stabilize the whole beam with two actuators at the root of the beam.

This work was presented at 2025 American Control Conference in Denver, CO and is repeated here because it is necessary background for the additinal material, the Kalman filter and aerodynamic model, that we present here. We add two sensors to the model, one which measures the vertical velocity at the tip of the wing and other that measures the rate of torsion also at the tip of the wing. We design a Kalman filter to process these measurements and return an estimate of the full state.

Finally we add a model of the aerodynamics forces that a bending and twisting beam generates to obtain a model of a wing. The model is based on the classic Wagner model for the aerodynamic forces exprienced by an airfoil in flight as described in [3]. We extend a standard finite dimensional model of an airfoil to a model of the complete wing and expand it using the same families of eigenfunctions that were used to describe the bending and twisting of the beam.

We present several simulations of the beam and wing to demonstrate the effectiveness of our approach.

II Dynamical System

Let the $y$ axis be the axis of rotation of the beam and suppose it extends from $y=0$ where it is attached to its support and its free end is at $y=L$. Let $w(y,t)$ be the vertical deflection of beam at location $y$ and time $t$ and let $\theta(y,t)$ be the angle of rotation of the beam around the $y$ axis at location $y$ and time $t$. According to [2], equations (3-155) and (3-156), the free vibrations of a uniform beam are governed by the two inertially coupled linear PDEs

[TABLE]

where

[TABLE]

The inertial coupling coefficient $S_{y}$ is zero if all the centers of gravity of the cross sections lie on the elastic axis.

The bending boundary conditions at the free end of the beam are

[TABLE]

and at the fixed end of the beam we assume that there is an actuator that can deliver a bending moment

[TABLE]

The torsion boundary condition at the free end of the beam is

[TABLE]

and at the fixed end of the beam we assume that there is an actuator that can deliver a torque

[TABLE]

We define

[TABLE]

We wish to express the dynamics as a first order system so we introduce a four vector valued variable

[TABLE]

then (7) becomes

[TABLE]

where

[TABLE]

The boundary conditions on $z(y,t)$ are

[TABLE]

We seek a feedback law of the form

[TABLE]

to stabilize the bending and torsion oscillations so we set up a Linear Quadratic Regulator (LQR). We choose a $4\times 4$ nonnegative definite matrix valued function $Q(y_{1},y_{2})$ that is symmetric in its arguments, $Q(y_{1},y_{2})=Q(y_{2},y_{1})$ and a positive definite $2\times 2$ matrix $R$. For a given initial condition $z(y,0)$ we seek to minimize by choice of $u(t)$ the quantity

[TABLE]

where ${\cal S}$ is the square $[0,L]^{2}$ and $dA=dy_{1}dy_{2}$.

Let $P(y_{1},y_{2})$ be an $4\times 4$ nonnegative definite matrix valued $C^{3}$ function that is symmetric in its arguments, $P(y_{1},y_{2})=P(y_{2},y_{1})$ and satisfies these homogeneous boundary conditions for $i,j=1,\ldots,4$

[TABLE]

These boundary conditions are analogous to the those on $z(y,t)$ when the the controls are zero and will be important when we integrate by parts below. We assume that the $P_{i_{1},i_{2}}(x_{1},x_{2})$ are $C^{3}$ so that these boundary conditions make sense.

LQR requires a stabilizability condition so we assume that for each $z(y,0)$ there is a $u(t)$ such that $w(y,t)\to 0$ and $\theta(y,t)\to 0$ as $t\to\infty$. Then by the Fundamental Theorem of Calculus

[TABLE]

LQR also requires a detectability condition so we assume that if $z(y,t)$ is such that

[TABLE]

then $z(y,t)\to 0$ as $t\to\infty$.

We bring the time differentiation inside the spatial integrals and integrate by parts several times to obtain

[TABLE]

where $P_{i,:}(y_{1},y_{2})$ and $P_{:,j}(y_{1},y_{2})$ denote the $i^{th}$ row and the $j^{th}$ column of $P(y_{1},y_{2})$ respectively.

We add the right side of (233) to the criterion (22) to get an equivalent criterion to be minimized. We wish to find a $4\times 2$ matrix valued function $K(y)$ such that the time integrand of equivalent criterion is equal to a perfect square of the form

[TABLE]

The terms quadratic in $u(t)$ match so we equate terms bilinear in $u^{\rm T}(t)$ and $z(y_{2},t)$. This yields

[TABLE]

so we assume that

[TABLE]

By symmetry

[TABLE]

Then by equating terms bilinear in $z^{\rm T}(y_{1},t)$ and $z(y_{2},t)$ we obtain the Riccati PDE for quadratic Fredholm kernel $P(y_{1},y_{2})$ of the optimal cost,

[TABLE]

where $\Gamma=BR^{-1}B$. This is an elliptic PDE with a quadratic nonlinearity.

III Fourier Analysis

We solve the Riccati PDE (54) by Fourier series. The partial differential operator $\frac{\partial^{2}}{\partial y^{2}}$ is self adjoint under the zero controlled boundary conditions

[TABLE]

so all the eigenvalues are real and the eigenfunctions are orthogonal. The eigenfunctions are $\Theta_{m}(y)=\cos{m\pi\over L}y$ and the eigenvalues are $\eta_{m}=-({m\pi\over L})^{2}$ for $n=0,1,2,\ldots$.

These are not the eigenpairs of the wave equation under the boundary conditions (60) but are related to them. The eigenpairs of the wave equation

[TABLE]

when written as a first order system are $\pm\sqrt{GJ\ }\ {n\pi i\over L}$ which are strictly imaginary if $m>0$. The corresponding eigenfunctions are

[TABLE]

The partial differential operator $-\frac{\partial^{4}}{\partial y^{4}}$ is self-adjoint when subject to the boundary conditions

[TABLE]

Note that these are not the boundary conditions of a cantilever beam. For a cantilever beam, the second boundary condition is $\frac{\partial\phi}{\partial y}(0)=0$, not $\frac{\partial^{2}\phi}{\partial y^{2}}(0)=0$

We look for eigenpairs $\nu,\Phi(x)$ that satisfy

[TABLE]

and the boundary conditions (62).

Also note that an eigenpair $\nu,\Phi(x)$ satisfying (63) and (62) is not an eigenpair of beam equation

[TABLE]

but it is related to them. When written as a first order system the beam eigenpairs are the eigenpairs of the matrix differential operator

[TABLE]

The beam eigenvalues are $\pm\sqrt{EI\nu\ \ }$ and since $EI\nu<0$ the beam eigenvalues are strictly imaginary. The corresponding beam eigenfunctions are vector valued,

[TABLE]

Because the differential operator $-\frac{\partial^{4}}{\partial y^{4}}(x)$ is self-adjoint under the zero controlled boundary conditions (62) it follows that all its eigenvalues $\nu_{n}$ are real and there is an orthogonal family of eigenfunctions, $\Phi_{n}(y)$. Regardless of the boundary conditions, any eigenfunction of the partial differential operator $-\frac{\partial^{4}}{\partial y^{4}}(x)$ must be of the form

[TABLE]

for some real $\nu,a,b,c,d$. Then the eigenvalue is $\lambda=-\nu^{4}$. We can express the boundary equations as a set of homogeneous linear equations depending on $\nu$,

[TABLE]

This system has a nontrivial solution iff $\nu$ is a root of the determinant of this matrix. The determinant is

[TABLE]

By adjusting the signs of $b$ and $d$ we can restrict our attention to nonnegative roots. There is no nonzero eigenfunction corresponding to $\nu=0$ so it is not an eigenvalue. The positive roots occur when $\tan\nu L=\tanh\nu L$. For each $n=1,2,\ldots$ there is exactly one root $\nu_{n}L\in[n\pi,(n+1/2)\pi)$. The first four roots are $\nu_{1}L=3.9266$, $\nu_{2}L=7.0686$, $\nu_{3}L=10.2102$, $\nu_{4}L=13.3518$. As $n\to\infty$ the $n^{th}$ root $\nu_{n}L$ is quickly converging to $(n\pi+{\pi\over 4})$.

Because of the boundary conditions (62) at $y=0$ we look for eigenfunctions are of the form

[TABLE]

Without loss of generality we can take $b=1$. Then the boundary conditions (62) at $y=L$ imply

[TABLE]

or equivalently

[TABLE]

For large $\nu_{n}L$ we have $\sinh\nu_{n}L\approx\cosh\nu_{n}L$ so we conclude that for large $n$, $\sin\nu_{n}L\approx\cos\nu_{n}L$. This happens when $\nu_{n}L\approx(n\pi+{\pi\over 4})$.

So the eigenfunctions $\Phi_{n}(y)$ are converging to .

[TABLE]

and

[TABLE]

We conjecture that the set of vectors $\{\Theta_{m}(y):m=0,1,2,3,\ldots\}\cup\{\Phi_{n}(y):n=1,2,3,\ldots\}$ is a Riesz basis for $L^{2}[0,L]$.

IV Series Solution of the Riccati PDE

Suppose $Q(y_{1},y_{2})$ has an expansion of the form

[TABLE]

We could consider more general $Q(y_{1},y_{2})$ but to keep the exposition relatively simple we do not. Notice that the ranges of the sums are different.

We also assume that the solution $P(y_{1},y_{2})$ of the Riccati PDE (54) has a similar but more complicated expansion.

[TABLE]

Note we are abusing notation, $P^{n_{1},n_{2}}$ is not necessarily the same as $P^{m_{1},m_{2}}$ even when $n_{1}=m_{1}$ and $n_{2}=m_{2}$. We use $n$ as a superscipt to indicate a coefficient of $\Phi_{n}(y)$ and we use $m$ as a superscipt to indicate a coefficient of $\Theta_{m}(y)$.

We plug these expansions into Riccati PDE and collect similar terms to obtain an infinite dimensional algebraic Riccati equation which has four uncoupled components. The $\Phi_{n_{1}}(y_{1})\Phi_{n_{2}}(y_{2})$ component is

[TABLE]

The $\Phi_{n_{1}}(y_{1})\Theta_{m_{2}}(y_{2})$ component is

[TABLE]

A simple solution to this equation is to take $P^{n,m}=0$ for all $n$ and $m$. Then by symmetry $P^{m,n}=0$. This would not be possible if we had chosen $Q^{m,n}\neq 0$.

Finally the $\Theta_{m_{1}}(y_{1})\Theta_{m_{2}}(y_{2})$ component is

[TABLE]

V Policy Iteration

We can approximately solve the algebraic Riccati equations (93) and (128) by policy iteration. The method can also be seen as value iteration. To find an initial iterate $\left(P^{n_{1},n_{2}}_{i,j}\right)^{(0)}$ and $\left(P^{m_{1},m_{2}}_{i,j}\right)^{(0)}$ we make the following assumptions.

$\left(P^{n_{1},n_{2}}_{i,j}\right)^{(0)}=0$ if $i=3,4$ or $j=3,4$. 2. 2.

$\left(P^{n_{1},n_{2}}_{i,j}\right)^{(0)}=0$ if $n_{1}\neq n_{2}$ when $i,j=1,2$. 3. 3.

$\left(P^{m_{1},m_{2}}_{i,j}\right)^{(0)}=0$ if $i=1,2$ or $j=1,2$. 4. 4.

$\left(P^{m_{1},m_{2}}_{i,j}\right)^{(0)}=0$ if $m_{1}\neq m_{2}$ when $i,j=3,4$.

These assumptions decouple bending from twisting so the four dimensional equations (93) and (128) reduce to two independent two dimensional equations. The Fourier expansion of the bending model is solely in terms $\Phi_{n}(y)$ and it satisfies

[TABLE]

Recall $\Phi^{\prime}_{n}(0)\to\nu_{n}$ as $n\to\infty$.

The Fourier expansion of the twisting model is solely in terms of $\Theta_{m}(y)$ and it satisfies

[TABLE]

The above assumptions on the initial iterate simplifies (137) to

[TABLE]

and simplifies (151) to

[TABLE]

Notice the $1,1$ component of (165) is a quadratic equation in $\left(P_{1,2}^{n,n}\right)^{(0)}=\left(P_{2,1}^{n,n}\right)^{(0)}$,

[TABLE]

Since it is harder to stabilize when $z_{1}(y,0)$ and $z_{2}(y,0)$ have the same sign we take the positive root of this quadratic to be the initial estimate

[TABLE]

The $2,2$ component of (165) is a quadratic equation in $\left(P^{n,n}_{2,2}\right)^{(0)}$

[TABLE]

and since $\left(P^{n,n}_{2,2}\right)^{(0)}$ must be nonnegative we define

[TABLE]

We can use the $1,2$ or the $2,1$ component of (165) to define $\left(P_{1,1}^{n,n}\right)^{(0)}$. Both lead to

[TABLE]

Then our initial iterate of the kernel of the bending part of the optimal cost is

[TABLE]

Theorem 1: The series (198) converges to a continuous $2\times 2$ matrix valued function. if there exist numbers $q>0$ and $r>8$ such that

[TABLE]

or $i,j=1,2$.

In addition if for some integer $\rho\geq 0$

[TABLE]

or $i,j=1,2$. then the series (198) converges to a $C^{\rho}$ $2\times 2$ matrix valued function.

Proof: The Mean Value Theorem applied to the function $g(s)=\sqrt{s}$ on the interval $[a,b]$ implies that there is an $s$ in $(a,b)$ such that

[TABLE]

Let $a=(\lambda_{n}EI)^{2}$ and $b=(\lambda_{n}EI)^{2}+Q_{1,1}(EI\Phi^{\prime}_{n}(0))^{2}\Gamma_{1,1}$. Then $\lambda_{n}EI=-\sqrt{(\lambda_{n}EI)^{2}}$ so there exists an $s$ between $(\lambda_{n}EI)^{2}$ and $(\lambda_{n}EI)^{2}+Q_{1,1}(EI\Phi^{\prime}_{n}(0))^{2}\Gamma_{1,1}$ such that

[TABLE]

The maximum of the right side of this equation between $(\lambda_{n}EI)^{2}$ and $(\lambda_{n}EI)^{2}+Q_{1,1}(EI\Phi_{n}(0)^{2}\Gamma_{1,1}$ occurs at $s=(\lambda_{n}EI)^{2}$ so

[TABLE]

Since $\left|\lambda_{n}\right|=\nu_{n}^{4}\approx\left({n\pi+\pi/4\over L}\right)^{4}$, $Q^{n,n}_{1,1}\leq{q\over n^{r}}$ with $r>8$ and $\left|\Phi_{n}(y)\right|\leq 2$ the series

[TABLE]

converges uniformly to a continuous function $P^{(0)}_{1,2}(y_{1},y_{2})$.

Now by (70) $\Phi_{n}^{\prime}(0)=\nu_{n}=O(n)$ so

[TABLE]

which is $O(n^{-1-r/2})$. Since $r>8$ the series

[TABLE]

converges uniformly to a continuous function $P^{(0)}_{2,2}(y_{1},y_{2})$.

Finally

[TABLE]

so $\left(P_{1,1}^{n,n}\right)^{(0)}=O(n^{3-r/2})$ and the series

[TABLE]

converges uniformly to a continuous function $P^{(0)}_{1,1}(y_{1},y_{2})$ because $3-r/2<-1$.

The assertion (200) is obtained by term by term differentiation of the series (198). QED

The way we interpret this result is that if we want the cost of moving the infinite number of imaginary eigenvalues into the open left half plane to be finite we must make $Q_{i,i}^{n,n}$ go to zero quite rapidly.

Notice the $3,3$ component of (179) is a quadratic equation in $P_{3,4}^{n,n}=P_{4,3}^{n,n}$

[TABLE]

Since it is harder to stabilize when $z_{3}(y,0)$ and $z_{4}(y,0)$ have the same sign we take the positive root of this quadratic to be the initial estimate

[TABLE]

The $4,4$ component of (179) is a quadratic equation in $P^{m,m}_{4,4}$

[TABLE]

and since $P^{m,m}_{4,4}$ must be nonnegative we define

[TABLE]

We can use the $3,4$ or the $4,3$ component of (179) to define $\left(P_{3,3}^{m,m}\right)^{(0)}$. Both lead to

[TABLE]

Our initial iterate of the kernel of the twisting part of the optimal cost is

[TABLE]

Theorem 2: The series (208) converges to a continuous $2\times 2$ matrix valued function if there exist numbers $q>0$ and $r>6$ such that for $m>0$

[TABLE]

$i,j=3,4$.

In addition if for some integer $\rho\geq 0$

[TABLE]

or $i,j=1,2$. then the series (208) converges to a $C^{\rho}$ $2\times 2$ matrix valued function.

Proof: Since $\eta_{m}=-\left({m\pi\over L}\right)^{2}$ the Mean Value Theorem applied to (202) implies that there is an $s$ between $\left(\eta_{m}GJ\right)^{2}$ and $\left(\eta_{m}GJ\right)^{2}+Q_{3,3}^{m,m}\Gamma_{2,2}\left(GJ\right)^{2}$ such that

[TABLE]

The maximum of the right side of this equation occurs at $s=\left(\eta_{m}GJ\right)^{2}$ so we conclude that

[TABLE]

and $\left(P_{3,4}^{m,m}\right)^{(0)}$ is of order ${1\over m^{8}}$. Since $\left|\Theta_{m}(y)\right|\leq 1$ the series

[TABLE]

converges uniformly to a continuous function.

Now from (204) we know that $\left(P_{4,4}^{m,m}\right)^{(0)}$ is of order ${1\over m^{r/2}}$ so the series

[TABLE]

converges uniformly to a continuous function.

Finally from (205) we see that $\left(P_{4,4}^{m,m}\right)^{(0)}$ is of order ${1\over m^{r/2-2}}$ and since $r>6$ we have ${r\over 2}-2>1$ so the series

[TABLE]

converges uniformly to a continuous function.

The statement (210) is obtained by term by term differentiation of the series (208). QED

Succesive iterates are found by plugging $\left(P^{n_{1},n_{2}}_{i,j}\right)^{(k)}$ and $\left(P^{m_{1},m_{2}}_{i,j}\right)^{(k)}$ into the right sides of the algebraic Riccati equations (93) and (128) and solving for $\left(P^{n_{1},n_{2}}_{i,j}\right)^{(k+1)}$ and $\left(P^{m_{1},m_{2}}_{i,j}\right)^{(k+1)}$ on the left side. A complication arises when $m_{1}=m_{2}=0$ because then $\nu_{0}=0$ and so the left sides is not an invertible expression in $\left(P^{0,0}_{2,:}\right)^{(k+1)}$ and $\left(P^{0,0}_{4,:}\right)^{(k+1)}$. A way around this is for all $m_{1}^{2}+m_{2}^{2}>0$ to plug $\left(P^{n_{1},n_{2}}_{i,j}\right)^{(k)}$ and $\left(P^{m_{1},m_{2}}_{i,j}\right)^{(k)}$ into the right sides of the algebraic Riccati equations (93) and (128) and solve for $\left(P^{n_{1},n_{2}}_{i,j}\right)^{(k+1)}$, $\left(P^{m_{1},m_{2}}_{2,:}\right)^{(k+1)}$ and $\left(P^{m_{1},m_{2}}_{4,:}\right)^{(k+1)}$ for all $n_{1},n_{2}$ and all $m_{1},m_{2}$ such that $m_{1}^{2}+m_{2}^{2}>0$. Then plug these $k+1$ iterates into (93) and (128) and solve the resulting quadratics for $\left(P^{m,m}_{2,:}\right)^{(k+1)}$ and $\left(P^{m,m}_{4,:}\right)^{(k+1)}$ when $m_{1}=m_{2}=0$. If the resulting quadratics do not determine $\left(P^{m,m}_{2,:}\right)^{(k+1)}$ and $\left(P^{m,m}_{4,:}\right)^{(k+1)}$ when $m_{1}=m_{2}=0$ just set them to zero. It should be noted that if $S_{y}\neq 0$ then the $k^{th}$ approximations to the bending and twisting kernels are not $2\times 2$ like (198) and (208) but are $4\times 4$ when $k>0$.

This is value (or policy) iteration so we know that the iterates are nonincreasing. If hypotheses of Theorems 1 and 2 are satisfied then the initial approximations (198) and (208) are continuous hence bounded on $[0,L]\times[0,L]$. Hence the initial estimate of the optimal cost is bounded as are all successive estimates. From this we conclude that all of the approximations to the optimal feedback move all the closed loop eigenvalues into the open left half plane. But we cannot conclude that the closed loop system is exponentially stable because the closed loop eigenvalues could (and do) approach the imaginary axis as the wave numbers $n$ and $m$ increase.

VI Finite Dimensional Approximating LQR

We construct a finite dimensional LQR whose algebraic Riccati equation is a truncation of (93) and (128). We choose an $N>0$ and construct a linear system with state $\zeta=[\zeta_{1},\zeta_{2},\zeta_{3},\zeta_{4}]^{\rm T}$ where $\zeta_{1},\zeta_{2},\zeta_{3},\zeta_{3}$ is each of dimension $N\times 1$. So the finite dimensional state $\zeta$ is of dimension $4N\times 1$. The dynamics is

[TABLE]

where

[TABLE]

This finite dimensional system approximates the infinite dimensional system in the following manner

[TABLE]

VII Example

We consider a $N=4$ example which leads to a $16$ dimensional approximation. For the time being we take all constants equal $1$ except $S_{y}=1/2$, $Q$ a $16\times 16$ identity matrix and $R$ a $2\times 2$ identity matrix. The $16$ open and closed loop poles are

[TABLE]

Notice how close the imaginary parts of the closed loop poles are to the open loop poles, All the closed loop poles have real parts less than $-1$.

Figure 1 is a simulation of the $16$ dimesional state under full state feedabck

Figure 2 shows the vertical displacement converging to zero under full state feedback. The control input at the root of the beam causes the ripples while stabilizing the vertical displacement of the beam.

Figure 3 shows the stabilization of the angle of attack under full state feedback.

Again the control input at the root of the beam causes the ripples while stabilizing the torsion of the beam.

VIII Kalman Filter

In the above sections we have constructed a full state feedback to asymptotically stabilize the linear system. But we don’t have the ability to measure the full state so need to construct a Kalman filter to deliver an estimate of the full state from a finite number of point measurements. We assume that we can measure $z_{2}(L,t)=\frac{\partial h}{\partial t}(L,t)$ and $z_{4}(L,t)=\frac{\partial\alpha}{\partial t}(L,t)$. We could use more sensors and/or different types of sensors but to keep thing simple we do not do so.

To develop the Kalman filter it is simpler to express the inertially coupled bending and twisting dynamics in a more standard form by multiplying the dynamics by the inverse of the coupling matrix $M$. We also add a distributed driving noise input $v(t)=[v_{1}(t),v_{2}(t)]^{\rm T}$ and obtain

[TABLE]

where the matrix partial differential operator is ${\cal A}=M^{-1}D$ and

[TABLE]

The boundary conditions on $z(y,t)$ are

[TABLE]

The measurements are the vertical velocity $\psi_{1}(t)=z_{2}(L,t)$ and angular velocity $\psi_{2}(t)=z_{4}(L,t)$ of the tip of the beam corrupted by white Gaussian noise,

[TABLE]

where

[TABLE]

and ${\cal D}$ is a $2\times 2$ invertible matrix with $w_{1}(t),w_{2}(t)$ standard independent white Gaussian noises. One way of obtaining such measurements is to integrate a pair of accelerometers as was done by Banks et al.[1].

Because of the Gaussian and linear assumptions we expect that when $u_{1}(t)=u_{2}(t)=0$ the optimal least squares estimate $\hat{z}(y,t)$ of $z(y,t)$ is a linear functional of the past observations,

[TABLE]

for some $4\times 2$ matrix valued function ${\cal L}(y,s)$. Notice that ${\cal L}(y,s)$ depends on the location $y$ because we are trying to estimate $z(y,t)$.

Given such a ${\cal L}(y,s)$ for each $y\in[0,L]$ we define a $4\times 4$ matrix valued function ${\cal H}(y,y_{1},s)$ by the differential equation

[TABLE]

The subscript on ${\cal A}_{1}$ indicates that it is involves partial differentiation with respect to $y_{1}$ not $y$.

This PDE is subject to the homogeneous boundary conditions

[TABLE]

and the initial condition

[TABLE]

Then

[TABLE]

Note the parenthesis around $\left({\cal H}(y,y_{1},s){\cal A}_{1}\right)$. The parenthesis implies that we apply the matrix partial differential operator ${\cal A}_{1}$ to ${\cal H}(y,y_{1},s)$ before doing anything else with this term.

Now we integrate by parts with respect to $s$

[TABLE]

We expect that the state $z(y_{1},t-s)$ in the far past $s>>0$ to be irrelevant to the estimate $\hat{z}(y_{1},t)$ of the state at time $t$ so we assume that ${\cal H}(y,y_{1},s)z(y_{1},t-s)\to 0$ as $s\to\infty$. Hence

[TABLE]

Then we integrate by parts with respect to $y_{1}$ using the boundary conditions

[TABLE]

Putting this all together yields

[TABLE]

So the estimation error $\tilde{z}(y,t)=z(y,t)-\hat{z}(y,t)$ is given by

[TABLE]

Because $v(t)$ and $w(t)$ are standard white Gaussian noises, the covariance of the error is

[TABLE]

For each $y\in[0,L]$ and for each pair of corresponding rows of ${\cal H}(y,y_{1},s)$ and ${\cal L}(y,s)$ we have an LQR in adjoint form with state ${\cal H}_{i,:}(y,y_{1},s)$, control ${\cal L}_{i,:}(y,s)$, linear dynamics (229) and a quadratic criterion (232). We can leave this adjoint LQR in matrix form as the optmal feedback gain ${\cal K}(y,y_{1})$ is the same for all rows,

[TABLE]

That is not to say that the rows of ${\cal L}(y,s)$ are all the same because the rows of ${\cal H}(y,y_{1},s)$ can be different, We are trying to estimate $z(y,t)$ which is why ${\cal H}(y,y_{1},s)$ and ${\cal K}(y,y_{1})$ might depand on $y$.

For each $y\in[0,L]$ let ${\cal P}(y,y_{1},y_{2})$ be a $6\times 6$ continuous matrix value function of $y_{1},y_{2}$ satisfying the homogeneous boundary conditions

[TABLE]

Using the Fundamental Theorem of Calculus assuming ${\cal H}(y,y_{1},s)\to 0$ as $s\to\infty$ we have

[TABLE]

We formally integrate by parts as if ${\cal P}(y,y_{1},y_{2})$ were $C^{2}$ in $y_{1},y_{2}$

[TABLE]

which reduces to

[TABLE]

We add the right side of (234) to the criterion (VIII) to be minimized to get an equivalent criterion

[TABLE]

We seek a ${\cal K}(y,y_{1})$ such that the time integrand in (VIII) is a perfect square of the form

[TABLE]

The terms quadratic in ${\cal L}(y,s)$ match so we equate terms bilinear in ${\cal L}(y,s)$ and ${\cal H}(y,y_{2},s)$ and obtain

[TABLE]

so we set

[TABLE]

Finally we match terms bilinear in ${\cal H}(y,y_{1},s)$ and ${\cal H}(y,y_{2},s)$ which yields the Riccati PDE of the infinite dimensional Kalman filter

[TABLE]

Actually this is a family of Riccati PDEs paramterized by $y$,

The optimal estimate is given by

[TABLE]

We make the substitution $s=t-\tau$ then

[TABLE]

We differentiate this with respect $t$ and obtain

[TABLE]

Then

[TABLE]

where ${\cal F}(y,y_{1})$ is the differential operator

[TABLE]

Put another way

[TABLE]

and so for each $y,y_{1}$

[TABLE]

We differentiate this with respect to $\tau$ and obtain

[TABLE]

so we conclude that

[TABLE]

and

[TABLE]

We differentiate the optimal estimate with respect $t$ and obtain

[TABLE]

so we conclude that

[TABLE]

or

[TABLE]

The boundary conditions on $\hat{z}(y,t)$ are the same as those (225) on ${z}(y,t)$,

[TABLE]

The so called innovations process is

[TABLE]

where

[TABLE]

This is the difference between the actual observation $\psi(t)$ and what we thing it should be $\hat{\psi}(t)$ given our optimal estimate of the state.

We can express the optimal filter as a copy of the original dynamics driven by the innovations

[TABLE]

subject to the boundary conditions (259).

If we use the estimate $\hat{z}(y,t)$ in place of the full state $z(y,t)$ in the feedback law we found before then the control input is

[TABLE]

We insert $u(t)$ into the plant boundary conditions (20). We also insert it into the boundary conditions (259) of the Kalman filter. This is dynamic compensation.

The error $\tilde{z}(y,t)=z(y,t)-\hat{z}(y,t)$ dynamics is given by

[TABLE]

subject to the homogeneous boundary conditions

[TABLE]

Notice that the error dynamics and its boundary conditions only depend on the error $\tilde{z}(y,t)$ so we can express the combined plant and Kalman filter in the coordinates $z(y,t)$ and $\tilde{z}(y,t)$. In these coordinates the combined system is upper triangular so if the full state feedback aymptotically stabilizes the system and if the error dynamics of the Kalman filter is asymptotically stable then the dynamic compensator stabilizes the system.

If the coefficient of the driving noise does not depend on $y$, ${\cal B}(y)={\cal B}$ then the family of LQRs does not depend on $y$. Therefore ${\cal L}(y,t)={\cal L}(t)$, ${\cal H}(y,y_{1},t)={\cal H}(y_{1},t)$, ${\cal K}(y,y_{1})={\cal K}(y_{1})$ and ${\cal P}(y,y_{1},y_{2})={\cal P}(y_{1},y_{2})$. The optimal filter is given by

[TABLE]

IX Filter Example

We apply the Kalman filter to the $16^{th}$ dimensional example that we treated above. The $16$ error variables are shown in Figure 4.

The vertical displacement error is shown in Figure 5 and the rotational displacement error is shown in Figure 6

Notice that the estimation errors are smaller near the sensors at the tip, $y=L$.

X Aerodynamic Model

The next step is to add a model of the aerodynamical forces that are generated by the bending and twisting beam to get a model of a wing.The classical aerodynamical models are those of Wagner and Theodoresen. Wagner’s model is in the time domain and Theodoresen’s is in the frequency domain. They are related by the Laplace transform.

These models are not exactly realizable by a linear time invariant system. But over the years considerable effort has gone into finding linear time invariant systems whose impulse response or transfer function approximates Wagner’s or Theodoresn’s model. We shall use the approximation to Wagner’s model developed by R. T. Jones in 1938 as reported in [3]. Wagner’s model is valid for an airfoil, a wing section. We shall extend it to model a wing by introducing spanwise dependence.

The state space realization of Wagner’s model adds two additional states to our four dimensional model. Of course each of these six state variables is a function of $y$ and $t$.

Aerodynamical models are parameterized by the free stream air velocity $U_{\infty}$ which we could treat as a static state. But then the model would become nonlinear so we do not do so..

Here is R.T. Jones’ approximate state space realization of Wagner’s impulse response. The aerodynamics state is $\xi=[\xi_{1}(y,t),\xi_{2}(y,t)]^{\rm T}$ and its dynamics is

[TABLE]

where

[TABLE]

The aerodynamic input is

[TABLE]

The aerodynamics output is

[TABLE]

where

[TABLE]

The lift per unit span is

[TABLE]

The moment per unit span is

[TABLE]

where $\psi(y,t)$ is the aerodynamic output.

Typical parameter values are

[TABLE]

The model is from [6] and the parameters are from [3].

We conjectured that the combination of $\Phi_{n}(y),\ n=1,2,3,\ldots$ and $\Theta_{m}(y),\ m=0.1,2,\ldots$ constitutes a Riesz basis for $L^{2}[0,L]$. We expand the aerodynamic state in these functions

[TABLE]

We differentiate with respect to time

[TABLE]

On the other hand

[TABLE]

[TABLE]

We multiply by the elements of the dual basis and get the modal dynamics

[TABLE]

The modal decomposition of the aerodynamic output is

[TABLE]

Our approximate model of beam has a sixteen dimensional state

[TABLE]

where $\zeta_{i}(t)$ is $4\times 1$. The approximate model of aerodynamics is also sixteen dimensional.

[TABLE]

where $\xi_{i}(t)$ is $4\times 1$. So the combined approximate model is $32$ dimensional.

Figure 7 shows the $16$ beam variables of $32$ variable combined model being stabilized by the $32$ dimensional full state feedback.

Figure 8 shows the $16$ aerodynamic variables of $32$ variable combined model being stabilized by the $32$ dimensional full state feedback. Notice how slow the stabilization is. The bending and twisting controls have very limited authority if any over the aerodynamical states. We can still use LQR because the eigenvalues of $A$ are $-0.0457$ and $-0.2998$ so aero model is lightly damped.

Figure 9 shows the $32$ dimensional full state feedback stabilizing the vertical displacement.

Figure 10 shows the $32$ dimensional full state feedback stabilizing the rotational displacement.

The combined beam and aero model is $32$ dimensional. The closed loop spectrum consists of $16$ complex eigenvalues and $16$ real eigenvalues. The complex eigenvalues are loosely associated to the closed loop dynamics of the beam while the real eigenvalues are loosely associated to the closed loop dynamics of the aero model. The real eigenvalues are very lightly damped, the largest one is $-0.04569$.

XI Conclusion

We started by considering the bending and twisting of a beam and showed how its oscillations can be damped by an LQR derived full state feedback using two point actuators located at the base of the beam. Then we derived a Kalman filter that processed two point measurements at the tip of the beam to obtain an estimate of the full state of the beam. This yields an Linear Quadratic Gaussian (LQG) synthesis. By adding a model of the aerodynamic forces generated by the bending and twisting of the beam we converted the model of the beam into a model of a wing and we showed that the LQR full state feedback stabilized the model of the wing.