Main Content

This example shows how to convert a fixed-wing aircraft to a linear time invariant (LTI) state-space model for linear analysis.

This example describes:

Importing and filling data from a DATCOM file.

Constructing a fixed-wing aircraft from DATCOM data.

Calculating static stability of the fixed-wing aircraft.

Linearizing the fixed-wing aircraft around an initial state.

Validating the static stability analysis with a dynamic response.

Isolating the elevator-to-pitch transfer function and designing a feedback controller for the elevator.

This example uses a DATCOM file created for the Sky Hogg aircraft.

First, import the DATCOM output file using datcomimport.

```
allData = datcomimport('astSkyHoggDatcom.out', false, 0);
skyHoggData = allData{1}
```

`skyHoggData = `*struct with fields:*
case: 'SKYHOGG BODY-WING-HORIZONTAL TAIL-VERTICAL TAIL CONFIG'
mach: [0.1000 0.2000 0.3000 0.3500]
alt: [1000 3000 5000 7000 9000 11000 13000 15000]
alpha: [-16 -12 -8 -4 -2 0 2 4 8 12]
nmach: 4
nalt: 8
nalpha: 10
rnnub: []
hypers: 0
loop: 2
sref: 225.8000
cbar: 5.7500
blref: 41.1500
dim: 'ft'
deriv: 'deg'
stmach: 0.6000
tsmach: 1.4000
save: 0
stype: []
trim: 0
damp: 1
build: 1
part: 0
highsym: 1
highasy: 0
highcon: 0
tjet: 0
hypeff: 0
lb: 0
pwr: 0
grnd: 0
wsspn: 18.7000
hsspn: 5.7000
ndelta: 5
delta: [-20 -10 0 10 20]
deltal: []
deltar: []
ngh: 0
grndht: []
config: [1x1 struct]
version: 1976
cd: [10x4x8 double]
cl: [10x4x8 double]
cm: [10x4x8 double]
cn: [10x4x8 double]
ca: [10x4x8 double]
xcp: [10x4x8 double]
cma: [10x4x8 double]
cyb: [10x4x8 double]
cnb: [10x4x8 double]
clb: [10x4x8 double]
cla: [10x4x8 double]
qqinf: [10x4x8 double]
eps: [10x4x8 double]
depsdalp: [10x4x8 double]
clq: [10x4x8 double]
cmq: [10x4x8 double]
clad: [10x4x8 double]
cmad: [10x4x8 double]
clp: [10x4x8 double]
cyp: [10x4x8 double]
cnp: [10x4x8 double]
cnr: [10x4x8 double]
clr: [10x4x8 double]
dcl_sym: [5x4x8 double]
dcm_sym: [5x4x8 double]
dclmax_sym: [5x4x8 double]
dcdmin_sym: [5x4x8 double]
clad_sym: [5x4x8 double]
cha_sym: [5x4x8 double]
chd_sym: [5x4x8 double]
dcdi_sym: [10x5x4x8 double]

Next, prepare the DATCOM lookup tables.

DATCOM lookup tables might have missing values due to the tables only filling one value for the whole column.

This missing data is represented in the lookup tables as the value 99999 and can be filled using the "previous" method of fillmissing.

In this example, ${\mathrm{C}}_{\mathrm{y}\beta}$, ${\mathrm{C}}_{\mathit{n}\beta}$, ${\mathrm{C}}_{\mathrm{Lq}}$, and ${\mathrm{C}}_{\mathrm{mq}}$ have missing data.

skyHoggData.cyb = fillmissing(skyHoggData.cyb, "previous", "MissingLocations", skyHoggData.cyb == 99999); skyHoggData.cnb = fillmissing(skyHoggData.cnb, "previous", "MissingLocations", skyHoggData.cnb == 99999); skyHoggData.clq = fillmissing(skyHoggData.clq, "previous", "MissingLocations", skyHoggData.clq == 99999); skyHoggData.cmq = fillmissing(skyHoggData.cmq, "previous", "MissingLocations", skyHoggData.cmq == 99999);

With the missing data filled, the fixed-wing aircraft can be constructed.

First, the fixed-wing aircraft is prepared with the desired aircraft name.

Optionally, the aircraft name can be extracted from the "case" field on the DATCOM struct by passing an empty fixed-wing object.

skyHogg = Aero.FixedWing(); skyHogg.Properties.Name = "Sky_Hogg"; skyHogg.DegreesOfFreedom = "3DOF"; [skyHogg, cruiseState] = datcomToFixedWing(skyHogg, skyHoggData);

The datcomToFixedWing will convert all compatible data from the datcom struct into the fixed-wing object and its state. However, the returned state still needs processing to get the desired initial conditions of the aircraft.

In this example, the environment, mass, inertia, airspeed, and center of pressure need adjusting.

h = 2000; [T, a, p, rho] = atmoscoesa(convlength(h, "ft", "m")); cruiseState.AltitudeMSL = h; cruiseState.Environment.Temperature = convtemp(T,"K", "R"); cruiseState.Environment.SpeedOfSound = convvel(a, "m/s", "ft/s"); cruiseState.Environment.Pressure = convpres(p, "pa", "psf"); cruiseState.Environment.Density = convdensity(rho, "kg/m^3", "slug/ft^3"); cruiseState.U = 169.42; cruiseState.Mass = 1299.214; cruiseState.Inertia.Variables = [5787.969 0 117.64;0 6928.93 0;-117.64 0 11578.329]; cruiseState.CenterOfPressure = [0.183, 0, 0];

Performing a static stability analysis helps determine the dynamic stability of the system without calculating a dynamic system response.

stability = staticStability(skyHogg, cruiseState)

`stability=`*6×8 table*
U V W Alpha Beta P Q R
________ _________ ________ ________ _________ _________ ________ _________
FX "Stable" "" "" "" "" "" "" ""
FY "" "Neutral" "" "" "" "" "" ""
FZ "" "" "Stable" "" "" "" "" ""
L "" "" "" "" "Neutral" "Neutral" "" ""
M "Stable" "" "" "Stable" "" "" "Stable" ""
N "" "" "" "" "Neutral" "" "" "Neutral"

With statically stable forces and moments, with perturbations to forward and vertical speeds, and perturbations to angle of attack, the dynamic stability of the system tends towards an oscillating steady-state when perturbing the forward and vertical speeds.

To verify this behavior, use the Control System Toolbox™.

To use the tools within the Control System Toolbox™, linearize the aircraft around a state.

This is done by using the linearize method with the same cruise state as before.

linSys = linearize(skyHogg, cruiseState)

linSys = A = XN XD U W Q XN 0 0 1 0 0 XD 0 0 0 1 0 U 0 1.319e-06 -0.001714 -0.0007608 0 W 0 0 -0.002705 -0.3319 2.557 Q 0 0 0.02379 -2.495 -26.41 Theta 0 0 0 0 1 Theta XN -2.586e-07 XD -2.957 U -0.5617 W -4.867e-08 Q 0 Theta 0 B = Delta XN 0 XD 0 U -0.0004314 W -0.02084 Q -2.321 Theta 0 C = XN XD U W Q Theta XN 1 0 0 0 0 0 XD 0 1 0 0 0 0 U 0 0 1 0 0 0 W 0 0 0 1 0 0 Q 0 0 0 0 1 0 Theta 0 0 0 0 0 1 D = Delta XN 0 XD 0 U 0 W 0 Q 0 Theta 0 Continuous-time state-space model.

With the linear state-space model constructed, you can plot the dynamic behavior of the system.

To verify the static stability results with the dynamic behavior of the system, plot the states space model against the initial conditions.

To induce a perturbation in the system, a 5 degree step for 1 second is added to the elevator signal.

x0 = getState(cruiseState, linSys.OutputName); t = linspace(0, 50, 500); u = zeros(size(t)); u(t > 1 & t < 2) = 5; lsim(linSys,u,t,x0)

As expected from the static stability analysis, the airspeed and pitch-rate is stable when responding to a small perturbation in the elevator.

In addition to the static stability verification, isolating the control surfaces to their intended dynamic response can help design controllers specific to the individual surfaces.

In this case, there is only a single control surface, the elevator.

The elevator controls the pitch response of the aircraft. To show a pitch response, isolate the elevator input to the pitch angle to elevator transfer function.

linSysElevatorTF = tf(linSys(6,1))

linSysElevatorTF = From input "Delta" to output "Theta": -2.321 s^3 - 0.7222 s^2 - 0.001232 s - 8.935e-09 ------------------------------------------------------------------ s^5 + 26.75 s^4 + 15.19 s^3 + 0.03932 s^2 + 0.008227 s + 5.712e-08 Continuous-time transfer function.

step(linSysElevatorTF)

As can be seen by the step plot, the pitch response to elevator input has an undesirable oscillatory nature and large steady-state error.

By adding a PID feedback controller to the elevator input, a much more desirable pitch response can be achieved.

`C = pidtune(linSysElevatorTF, "PID")`

C = 1 Ki * --- s with Ki = -0.000524 Continuous-time I-only controller.

elevatorFeedback = feedback(linSysElevatorTF * C, 1)

elevatorFeedback = From input to output "Theta": 0.001217 s^3 + 0.0003786 s^2 + 6.459e-07 s + 4.684e-12 ------------------------------------------------------------------------ s^6 + 26.75 s^5 + 15.19 s^4 + 0.04053 s^3 + 0.008605 s^2 + 7.03e-07 s + 4.684e-12 Continuous-time transfer function.

step(elevatorFeedback)