UNIVERSITY OF CALIFORNIA, SANTA CRUZ

BOARD OF STUDIES IN COMPUTER ENGINEERING

CMPE 240: INTRODUCTION TO LINEAR DYNAMICAL SYSTEMS
Gabriel Hugh Elkaim Fall 2005
Aircraft Dynamics Example
longitudinal aircraft dynamics
wind gust & control inputs
linearized dynamics
steady-state analysis
eigenvalues & modes
impulse matrices
1 Longitudinal aircraft dynamics
PSfrag replacements
body axis
horizontal

variables are (small) deviations from operating point or trim conditions

state (components):
u: velocity of aircraft along body axis
v: velocity of aircraft perpendicular to body axis
(down is positive)
: angle between body axis and horizontal
(up is positive)
q =

: angular velocity of aircraft (pitch rate)
1
2 Inputs
disturbance inputs:
u
w
: velocity of wind along body axis
v
w
: velocity of wind perpendicular to body axis
control or actuator inputs:

e
: elevator angle (
e
> 0 is down)

t
: thrust
3 Linearized dynamics
for 747, level ight, 40000 ft, 774 ft/sec,

u
v
q

.003 .039 0 .322

.065 .319 7.74 0
.020 .101 .429 0
0 0 1 0

u u
w
v v
w
q

.01 1
.18 .04
1.16 .598
0 0

units: ft, sec, crad (= 0.01rad 0.57

)
matrix coecients are called stability derivatives
outputs of interest:
aircraft speed u (deviation from trim)
climb rate

h = v + 7.74
4 Steady-state analysis
DC gain from (u
w
, v
w
,
e
,
t
) to (u,

h):
H(0) = CA
1
B =

1 0 27.2 15.0
0 1 1.34 24.9

gives steady-state change in speed & climb rate due to wind, elevator & thrust changes
solve for control variables in terms of wind velocities, desired speed & climb rate

.0379 0.0229
.0020 .0413

u u
w

h + v
w

2
level ight, increase in speed is obtained mostly by increasing elevator (i.e., downwards)
constant speed, increase in climb rate is obtained by increasing thrust and increasing
elevator (i.e., downwards)
(thrust on 747 gives strong pitch up torque)
5 Eigenvalues and modes
eigenvalues are
0.3750 0.8818j, 0.0005 0.0674j
two complex modes, called short-period and phugoid, respectively
system is stable (but lightly damped)
hence step responses converge (eventually) to DC gain matrix
eigenvectors are
x
short
=

0.0005
0.5433
0.0899
0.0283

0.0135
0.8235
0.0677
0.1140

,
x
phug
=

0.7510
0.0962
0.0111
0.1225

0.6130
0.0941
0.0082
0.1637

6 Short-period mode
y(t) = Ce
tA
(x
short
) (pure short-period mode motion)
0 2 4 6 8 10 12 14 16 18 20
1
0.5
0
0.5
1
0 2 4 6 8 10 12 14 16 18 20
1
0.5
0
0.5
1
PSfrag replacements
u
(
t
)

h
(
t
)
3
only small eect on speed u
period 7 sec, decays in 10 sec
6.1 Phugoid mode
y(t) = Ce
tA
(x
phug
) (pure phugoid mode motion)
0 200 400 600 800 1000 1200 1400 1600 1800 2000
2
1
0
1
2
0 200 400 600 800 1000 1200 1400 1600 1800 2000
2
1
0
1
2
PSfrag replacements
u
(
t
)

h
(
t
)
aects both speed and climb rate
period 100 sec; decays in 5000 sec
7 Dynamic response to wind gusts
impulse response matrix from (u
w
, v
w
) to (u,

h) (gives response to short wind bursts)
over time period [0, 20]:
0 5 10 15 20
0.1
0
0.1
0 5 10 15 20
0.1
0
0.1
0 5 10 15 20
0.5
0
0.5
0 5 10 15 20
0.5
0
0.5
PSfrag replacements
h
1
1
h
1
2
h
2
1
h
2
2
4
over time period [0, 600]:
0 200 400 600
0.1
0
0.1
0 200 400 600
0.1
0
0.1
0 200 400 600
0.5
0
0.5
0 200 400 600
0.5
0
0.5
PSfrag replacements
h
1
1
h
1
2
h
2
1
h
2
2
8 Dynamic response to actuators
impulse response matrix from (
e
,
t
) to (u,

h)
over time period [0, 20]:
0 5 10 15 20
2
1
0
1
2
0 5 10 15 20
2
1
0
1
2
0 5 10 15 20
5
0
5
0 5 10 15 20
0
0.5
1
1.5
2
2.5
3
PSfrag replacements
h
1
1
h
1
2
h
2
1
h
2
2
over time period [0, 600]:
5
0 200 400 600
2
1
0
1
2
0 200 400 600
2
1
0
1
2
0 200 400 600
3
2
1
0
1
2
3
0 200 400 600
3
2
1
0
1
2
3
PSfrag replacements
h
1
1
h
1
2
h
2
1
h
2
2
9 MATLAB code
% 747 longitudinal axis example for 263
% 40000 feet steady, level flight, 774 ft/sec
% from bryson p151
% x = u, w, q, theta
A = [ -0.003 0.039 0 -.322; -0.065 -0.319 7.74 0; 0.020 -.101 -.429
0 0 0 1 0];
%input = u_w, w_w, delta_e, delta_t
Bw= -A(:,[1,2]); Bc= [0.01 1; -.18 -.04; -1.16 .598; 0 0]; B = [Bw,
Bc];
% output: u, climb rate = -w + 7.74 theta
C = [ 1 0 0 0; 0 -1 0 7.74]; H0 = -C*inv(A)*B;
H01 = H0(:,[3,4]); % DC gain matrix from delta_e delta_t to
% speed, climb rate
% modal analysis
[V,Gam]=eig(A); xshort = real(V(:,1)); xphug = real(V(:,3)); xshort
= xshort/norm(xshort); xphug = xphug/norm(xphug);
Nsamp = 100; %number of time samples
yshort=zeros(2,Nsamp); t=linspace(0,20,Nsamp); for i=1:Nsamp,
yshort(:,i)=C*expm(t(i)*A)*xshort; end
figure(3) subplot(2,1,1) plot(t,yshort(1,:)); axis([0 20 -1 1])
ylabel(u) subplot(2,1,2) plot(t,yshort(2,:)); axis([0 20 -1 1])
ylabel(hdot)
6
print -deps aircraft_short
Nsamp = 400; %number of time samples
yshort=zeros(2,Nsamp); t=linspace(0,2000,Nsamp); for i=1:Nsamp,
yphug(:,i)=C*expm(t(i)*A)*xphug; end
figure(4) subplot(2,1,1) plot(t,yphug(1,:)); axis([0 2000 -2 2])
ylabel(u) subplot(2,1,2) plot(t,yphug(2,:)); axis([0 2000 -2 2])
ylabel(hdot)
print -deps aircraft_phug
% now do responses to various impulses
figure(1)
Nsamp = 200; %number of time samples
h1=zeros(2,Nsamp); h2=zeros(2,Nsamp); t=linspace(0,20,Nsamp); for
i=1:Nsamp,
h1(:,i)=C*expm(t(i)*A)*B(:,1); % impulse response from u_w
h2(:,i)=C*expm(t(i)*A)*B(:,2); % imp resp from v_w
end
subplot(2,2,1) plot(t,h1(1,:)); axis([0 20 -.1 0.1]) ylabel(h11)
subplot(2,2,2) plot(t,h2(1,:)); axis([0 20 -.1 0.1]) ylabel(h12)
subplot(2,2,3) plot(t,h1(2,:)); axis([0 20 -.5 0.5]) ylabel(h21)
subplot(2,2,4) plot(t,h2(2,:)); axis([0 20 -.5 0.5]) ylabel(h22)
print -deps aircraft_gust1
% now do same plots over longer time scale
figure(2); t=linspace(0,600,Nsamp); for i=1:Nsamp,
h1(:,i)=C*expm(t(i)*A)*B(:,1); % impulse response from u_w
h2(:,i)=C*expm(t(i)*A)*B(:,2); % imp resp from v_w
end
subplot(2,2,1) plot(t,h1(1,:)); axis([0 600 -.1 0.1]) ylabel(h11)
subplot(2,2,2) plot(t,h2(1,:)); axis([0 600 -.1 0.1]) ylabel(h12)
subplot(2,2,3) plot(t,h1(2,:)); axis([0 600 -.5 0.5]) ylabel(h21)
subplot(2,2,4) plot(t,h2(2,:)); axis([0 600 -.5 0.5]) ylabel(h22)
print -deps aircraft_gust2
% now do same things, but for actuator inputs
figure(1)
7
Nsamp = 200; %number of time samples
h1=zeros(2,Nsamp); h2=zeros(2,Nsamp); t=linspace(0,20,Nsamp); for
i=1:Nsamp,
h1(:,i)=C*expm(t(i)*A)*B(:,3); % impulse response from delta_e
h2(:,i)=C*expm(t(i)*A)*B(:,4); % imp resp from delta_t
end
subplot(2,2,1) plot(t,h1(1,:)); axis([0 20 -2 2]) ylabel(h11)
subplot(2,2,2) plot(t,h2(1,:)); axis([0 20 -2 2]) ylabel(h12)
subplot(2,2,3) plot(t,h1(2,:)); axis([0 20 -5 5]) ylabel(h21)
subplot(2,2,4) plot(t,h2(2,:)); axis([0 20 0 3]) ylabel(h22)
print -deps aircraft_act1
% now do same plots over longer time scale
figure(2); t=linspace(0,600,Nsamp); for i=1:Nsamp,
h1(:,i)=C*expm(t(i)*A)*B(:,3); % impulse response from delta_e
h2(:,i)=C*expm(t(i)*A)*B(:,4); % imp resp from delta_t
end
subplot(2,2,1) plot(t,h1(1,:)); axis([0 600 -2 2]) ylabel(h11)
subplot(2,2,2) plot(t,h2(1,:)); axis([0 600 -2 2]) ylabel(h12)
subplot(2,2,3) plot(t,h1(2,:)); axis([0 600 -3 3]) ylabel(h21)
subplot(2,2,4) plot(t,h2(2,:)); axis([0 600 -3 3]) ylabel(h22)
print -deps aircraft_act2
courtesy of Stephen Boyd @ Stanford University
8