Store 
System Identification Toolbox 7.3.1
A Guided Missile
In this case study, we investigate a rather large, quite involved, and nonlinear system. The purpose is to demonstrate IDNLGREY's ability to estimate quite a few parameters (16) in a system having many inputs (10) and outputs (5). The system considered is a guided missile, which was originally a generic missile studied and modeled by Saab Bofors Dynamics AB, Sweden. The following figure illustrates the modeling situation.
Figure 1: A guided missile.
Contents
Input-Output Data
As a first action we load the available missile data:
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'missiledat a'));
This file contains one data set with 501 samples, where data was equidistantly sampled. The sampling period was 0.02 seconds. This data set was generated from an even more elaborate model of the true missile system. With this as a background, we next create an IDDATA object to represent and hold the data: (y and u are available in missiledata.mat)
z = iddata(y, u, 0.02, 'Name', 'Missile'); set(z, 'InputName', {'Aileron angle' 'Elevator angle' ... 'Rudder angle' 'Dynamic pressure' ... 'Missile velocity' ... 'Measured angular velocity around x-axis' ... 'Measured angular velocity around y-axis' ... 'Measured angular velocity around z-axis' ... 'Measured angle of attack' ... 'Measured angle of sideslip'}, ... 'InputUnit', {'rad' 'rad' 'rad' 'kg/(m*s^2)' 'm/s' ... 'rad/s' 'rad/s' 'rad/s' 'rad' 'rad'}); set(z, 'OutputName', {'Angular velocity around x-axis' ... 'Angular velocity around y-axis' ... 'Angular velocity around z-axis' ... 'Acceleration in y-direction' ... 'Acceleration in z-direction'}, ... 'OutputUnit', {'rad/s' 'rad/s' 'rad/s' 'm/s^2' 'm/s^2'}); set(z, 'Tstart', 0, 'TimeUnit', 's');
The input data is shown in a figure.
figure('Name', [z.Name ': input data']); for i = 1:z.Nu subplot(z.Nu/2, 2, i); plot(z(:, [], i)); title(['Input #' num2str(i) ': ' z.InputName{i}]); xlabel(''); axis('tight'); if (i > z.Nu-2) xlabel([z.Domain ' (' z.TimeUnit ')']); end end
Figure 2: Measured inputs to the guided missile.
The output data is shown in another figure.
figure('Name', [z.Name ': output data']); for i = 1:z.Ny subplot(z.Ny, 1, i); plot(z(:, i, [])); title(['Output #' num2str(i) ': ' z.OutputName{i}]); xlabel(''); axis('tight'); end xlabel([z.Domain ' (' z.TimeUnit ')']);
Figure 3: Measured outputs of the guided missile.
At a first glance, it might here look strange to have measured variants of some of the outputs in the input vector. However, the used missile model involves several integrators, which often cause an unstable simulation behavior. To avoid this, the measurements are fed back via a nonlinear observer. Hence the notion of measured input signals.
Modeling the Guided Missile
The next step is to define an IDNLGREY model object, i.e., a state space structure with the ability to model the missile dynamics. A reasonable structure is obtained by the use of Newton's basic force and momentum laws (balance equations), which, to form a complete model structure, are complemented with a few basic aerodynamic relations (constitutive relations).
Omitting the modeling details we just display the resulting state and output equations, encoded in the C MEX-file missile_c.c, and make the observations that the structure is quite complex as well as nonlinear.
/* State equations. */
void compute_dx(double *dx, double *x, double *u, double **p)
{
/* Retrieve model parameters. */
double *F, *M, *C, *d, *A, *I, *m, *K;
F = p[0]; /* Aerodynamic force coefficient. */
M = p[1]; /* Aerodynamic momentum coefficient. */
C = p[2]; /* Aerodynamic compensation factor. */
d = p[3]; /* Missile diameter. */
A = p[4]; /* Missile reference area. */
I = p[5]; /* Moment of inertia, x-y-z. */
m = p[6]; /* Missile mass. */
K = p[7]; /* Feedback gain. */ /* x[0]: Angular velocity around x-axis. */
/* x[1]: Angular velocity around y-axis. */
/* x[2]: Angular velocity around z-axis. */
/* x[3]: Angle of attack. */
/* x[4]: Angle of sideslip. */
dx[0] = 1/I[0]*(d[0]*A[0]*(M[0]*x[4]+0.5*M[1]*d[0]*x[0]/u[4]+M[2]*u[0]
)*u[3]-(I[2]-I[1])*x[1]*x[2])+K[0]*(u[5]-x[0]);
dx[1] = 1/I[1]*(d[0]*A[0]*(M[3]*x[3]+0.5*M[4]*d[0]*x[1]/u[4]+M[5]*u[1]
)*u[3]-(I[0]-I[2])*x[0]*x[2])+K[0]*(u[6]-x[1]);
dx[2] = 1/I[2]*(d[0]*A[0]*(M[6]*x[4]+M[7]*x[3]*x[4]+0.5*M[8]*d[0]*x[2]
/u[4]+M[9]*u[0]+M[10]*u[2])*u[3]-(I[1]-I[0])*x[0]*x[1])+K[0]*(u[7]-x[2]);
dx[3] = (-A[0]*u[3]*(F[2]*x[3]+F[3]*u[1]))/(m[0]*u[4])-x[0]*x[4]+x[1]+
K[0]*(u[8]/u[4]-x[3])+C[0]*pow(x[4],2);
dx[4] = (-A[0]*u[3]*(F[0]*x[4]+F[1]*u[2]))/(m[0]*u[4])-x[2]+x[0]*x[3]+
K[0]*(u[9]/u[4]-x[4]);
} /* Output equations. */
void compute_y(double *y, double *x, double *u, double **p)
{
/* Retrieve model parameters. */
double *F, *A, *m;
F = p[0]; /* Aerodynamic force coefficient. */
A = p[4]; /* Missile reference area. */
m = p[6]; /* Missile mass. */ /* y[0]: Angular velocity around x-axis. */
/* y[1]: Angular velocity around y-axis. */
/* y[2]: Angular velocity around z-axis. */
/* y[3]: Acceleration in y-direction. */
/* y[4]: Acceleration in z-direction. */
y[0] = x[0];
y[1] = x[1];
y[2] = x[2];
y[3] = -A[0]*u[3]*(F[0]*x[4]+F[1]*u[2])/m[0];
y[4] = -A[0]*u[3]*(F[2]*x[3]+F[3]*u[1])/m[0];
}Before being able to create an IDNLGREY object we must also provide values of the believed (initial) parameters, in this case 23 different ones. We lump some of the parameters (aerodynamic force coefficients, aerodynamic momentum coefficients, and moment of inertia factors) into vectors, so that 8 different parameter objects are specified. The initial parameter values were obtained partly by physical reasoning and partly by quantitative guessing. The last 5 parameters (A, I, m and K) are more or less constants, so by fixing these we get a model structure with 16 free parameters, lumped into F, M, C and d).
Parameters = struct('Name', {'Aerodynamic force coefficient' ... % F, 4-by-1 vector. 'Aerodynamic momentum coefficient' ... % M, 11-by-1 vector. 'Aerodynamic compensation factor' ... % C, scalar. 'Missile diameter' ... % d, scalar. 'Missile reference area' ... % A, scalar. 'Moment of inertia, x-y-z' ... % I, 3-by-1 vector. 'Missile mass' ... % m, scalar. 'Feedback gain'}, ... % K, scalar. 'Unit', {'1/(rad*m^2), 1/(rad*m^2), 1/(rad*m^2), 1/ (rad*m^2)' ... ['1/rad, 1/(s*rad), 1/rad, 1/rad, ' ... '1/(s*rad), 1/rad, 1/rad, 1/rad^2, ' ... '1/(s*rad), 1/rad, 1/rad'] ... '1/(s*rad)' 'm' 'm^2' ... 'kg*m^2, kg*m^2,kg*m^2' 'kg' '-'}, ... 'Value', {[20.0; -6.0; 35.0; 13.0] ... [-1.0; 15; 3.0; -16.0; -1800; -50; 23.0; -20 0; -2000; -17.0; -50.0] ... -5.0 0.17 0.0227 ... [0.5; 110; 110] 107 6}, ... 'Minimum', {-Inf(4, 1) -Inf(11, 1) -Inf -Inf -Inf - Inf(3, 1) -Inf -Inf}, ... % Ignore constraints. 'Maximum', {Inf(4, 1) Inf(11, 1) Inf Inf Inf Inf(3, 1) Inf Inf}, ... % Ignore constraints. 'Fixed', {false(4, 1) false(11, 1) false true tru e true(3, 1) true true});
We also define the 5 states of the missile model structure in the same manner:
InitialStates = struct('Name', {'Angular velocity around x-axis' ... 'Angular velocity around y-axis' . .. 'Angular velocity around z-axis' . .. 'Angle of attack' 'Angle of sideslip'}, . .. 'Unit', {'rad/s' 'rad/s' 'rad/s' 'rad' 'rad'} , ... 'Value', {0 0 0 0 0}, ... 'Minimum', {-Inf -Inf -Inf -Inf -Inf}, ... % Ignore constraints. 'Maximum', {Inf Inf Inf Inf Inf}, ... % Ignore constraints. 'Fixed', {true true true true true});
The model file along with order, parameter and initial states data are now used to create an IDNLGREY object describing the missile system:
FileName = 'missile_c'; % File describing the model structur e. Order = [5 10 5]; % Model orders [ny nu nx]. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Missile', 'TimeUnit', 's');
The input and output signals of the missile system are next specified by simply copying the corresponding bookkeeping data from the IDDATA object:
set(nlgr, 'InputName', z.InputName, 'InputUnit', z.InputUnit); set(nlgr, 'OutputName', z.OutputName, 'OutputUnit', z.OutputUnit);
We thus have an IDNLGREY object involving 10 input signals, 5 states, and 5 output signals. As mentioned before the model also contains 23 parameters, 7 of which are fixed and 16 of which are free:
nlgr
Time-continuous nonlinear state-space model defined by 'missile_c' (MEX-file
):
dx/dt = F(t, u(t), x(t), p1, ..., p8)
y(t) = H(t, u(t), x(t), p1, ..., p8) + e(t)
with 10 inputs, 5 states, 5 outputs, and 16 free parameters (out of 23).
Performance of the Initial Guided Missile Model
Before estimating the 16 free parameters we next simulate the system using the initial parameter vector. By this we can gain useful information about the quality of the initial model.
figure; compare(z, nlgr); set(gcf,'Position',get(gcf,'Position')+[0 -287 0 287]) % enlarge figure to s how axes clearly
Figure 4: Comparison between measured outputs and the simulated outputs of the initial guided missile model.
From the plot we see that the agreement between the true and the simulated signals is not that bad, except maybe for the time span 4 to 6 seconds. This fact is clearly revealed in a plot of the prediction errors:
figure; pe(z, nlgr);
Figure 5: Prediction errors of the initial guided missile model.
Parameter Estimation
It is not that far-fetched to believe that this initial model is a plausible starting point for parameter estimation. The prediction error estimate of the 16 free parameters is next computed. This computation will take some time.
duration = clock; nlgr = pem(z, nlgr, 'Display', 'on'); duration = etime(clock, duration);
Criterion: Trace minimization
Scheme: Trust-Region Reflective Newton (LSQNONLIN, LargeScale = 'On')
--------------------------------------------------------------
Norm of First-order
Iteration Cost step optimality
--------------------------------------------------------------
0 1389.23 - -
1 576.647 10 4.41e+005
2 409.863 10 3.82e+004
3 369.655 10 9.3e+003
4 345.151 20 1.38e+004
5 323.765 40 1.08e+004
6 294.091 80 1.48e+004
7 241.874 160 1.98e+004
8 100.93 320 2.25e+005
9 34.6073 640 9.26e+005
10 7.3802 1.28e+003 1.14e+005
11 7.3802 1.52e+003 1.14e+005
12 4.04828 381 6.95e+004
13 1.26636 761 5.16e+004
14 0.375938 301 6.27e+003
15 0.356682 431 3.59e+003
16 0.352838 105 441
17 0.352269 71.5 256
18 0.352251 37.4 166
19 0.352214 9.34 45
20 0.352212 4.88 24.5
--------------------------------------------------------------
Performance of the Estimated Guided Missile Model
On the used computer, estimation of the parameters took the following amount of time to complete.
disp(['Estimation time : ' num2str(duration, 4) ' seconds']); disp(['Time per iteration: ' num2str(duration/nlgr.EstimationInfo.Iterations , 4) ' seconds.']);
Estimation time : 12.89 seconds Time per iteration: 0.6445 seconds.
To evaluate the quality of the estimated model and to illustrate the improvement compared to the initial model we next perform a simulation of the estimated model. A comparison between the measured and the simulated outputs are shown in a figure.
figure; compare(z, nlgr); set(gcf,'Position',get(gcf,'Position')+[0 -287 0 287])
Figure 6: Comparison between measured outputs and the simulated outputs of the estimated guided missile model.
The figure clearly indicates the improvement compared to the simulation result obtained with the initial model. The previously badly modeled region spanning the interval 4 to 6 seconds is without doubt captured with much higher accuracy than before. This is best displayed by looking at the prediction errors:
figure; pe(z, nlgr);
Figure 7: Prediction errors of the estimated guided missile model.
Let us finally conclude the case study by textually presenting some details about the estimated missile model:
present(nlgr);
Time-continuous nonlinear state-space model defined by 'missile_c' (MEX-file
):
dx/dt = F(t, u(t), x(t), p1, ..., p8)
y(t) = H(t, u(t), x(t), p1, ..., p8) + e(t)
with 10 inputs, 5 states, 5 outputs, and 16 free parameters (out of 23).
Inputs:
u(1) Aileron angle(t) [rad]
u(2) Elevator angle(t) [rad]
u(3) Rudder angle(t) [rad]
u(4) Dynamic pressure(t) [kg/(m*s^2)]
u(5) Missile velocity(t) [m/s]
u(6) Measured angular velocity around x-axis(t) [rad/s]
u(7) Measured angular velocity around y-axis(t) [rad/s]
u(8) Measured angular velocity around z-axis(t) [rad/s]
u(9) Measured angle of attack(t) [rad]
u(10) Measured angle of sideslip(t) [rad]
States: initial value
x(1) Angular velocity around x-axis(t) [rad/s] xinit@exp1 0 (fix)
in [-Inf, Inf]
x(2) Angular velocity around y-axis(t) [rad/s] xinit@exp1 0 (fix)
in [-Inf, Inf]
x(3) Angular velocity around z-axis(t) [rad/s] xinit@exp1 0 (fix)
in [-Inf, Inf]
x(4) Angle of attack(t) [rad] xinit@exp1 0 (fix)
in [-Inf, Inf]
x(5) Angle of sideslip(t) [rad] xinit@exp1 0 (fix)
in [-Inf, Inf]
Outputs:
y(1) Angular velocity around x-axis(t) [rad/s]
y(2) Angular velocity around y-axis(t) [rad/s]
y(3) Angular velocity around z-axis(t) [rad/s]
y(4) Acceleration in y-direction(t) [m/s^2]
y(5) Acceleration in z-direction(t) [m/s^2]
Parameters: value st
andard dev
p1(1) Aerodynamic force coefficient [1/(rad*m^2..] 21.3246 0
.348672 (est) in [-Inf, Inf]
p1(2) -7.78593 0
.186189 (est) in [-Inf, Inf]
p1(3) 36.2971 0
.752835 (est) in [-Inf, Inf]
p1(4) 13.0201
1.41627 (est) in [-Inf, Inf]
p2(1) Aerodynamic momentum coefficient [1/rad, 1/(..] -1.03446 0.
0750181 (est) in [-Inf, Inf]
p2(2) 15.337 0
.882165 (est) in [-Inf, Inf]
p2(3) 2.96768 0
.202111 (est) in [-Inf, Inf]
p2(4) -18.7409 0
.356223 (est) in [-Inf, Inf]
p2(5) -1160.03
234.007 (est) in [-Inf, Inf]
p2(6) -57.0413
1.36136 (est) in [-Inf, Inf]
p2(7) 34.2341
1.45059 (est) in [-Inf, Inf]
p2(8) -208.39
8.23738 (est) in [-Inf, Inf]
p2(9) -2480.44
278.64 (est) in [-Inf, Inf]
p2(10) -32.5455
3.165 (est) in [-Inf, Inf]
p2(11) -50.9551
1.6996 (est) in [-Inf, Inf]
p3 Aerodynamic compensation factor [1/(s*rad)] -0.651994 0
.715564 (est) in [-Inf, Inf]
p4 Missile diameter [m] 0.17
0 (fix) in [-Inf, Inf]
p5 Missile reference area [m^2] 0.0227
0 (fix) in [-Inf, Inf]
p6(1) Moment of inertia, x-y-z [kg*m^2, kg..] 0.5
0 (fix) in [-Inf, Inf]
p6(2) 110
0 (fix) in [-Inf, Inf]
p6(3) 110
0 (fix) in [-Inf, Inf]
p7 Missile mass [kg] 107
0 (fix) in [-Inf, Inf]
p8 Feedback gain [-] 6
0 (fix) in [-Inf, Inf]
The model was estimated from the data set 'Missile', which
contains 501 data samples.
Loss function 3.45246e-010 and Akaike's FPE 3.67297e-010
Created: 30-Jun-2009 16:10:34
Last modified: 30-Jun-2009 16:10:55
Concluding Remarks
The estimated guided missile model constitutes a decent basis for investigating the fundamental performance of different missile guidance control strategies. High-quality missile models, that preferably have a physical interpretation, are, e.g., vital components of so-called model predictive control systems.
Additional Information
For more information on identification of dynamic systems with System Identification Toolbox™ visit the System Identification Toolbox product information page.
