Моделирование аэродинамического тела
Этот пример показывает моделирование серого поля большого, и объедините нелинейную систему. Цель состоит в том, чтобы показать способность использовать модель IDNLGREY, чтобы оценить большое количество параметров (16) в системе, имеющей много входных параметров (10) и выходные параметры (5). Система является аэродинамическим телом. Мы создаем модель, которая предсказывает ускорение и скорость тела с помощью измерений его скоростей (переводный и угловой) и различные углы, связанные с его поверхностями управления.
Данные ввода - вывода
Мы загружаем измеренные скорости, углы и динамическое давление файла данных, названного aerodata.mat:
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'aerodata'));
Этот файл содержит однородно набор выборочных данных с 501 выборкой в переменных u и y. Шаг расчета составляет 0,02 секунды. Этот набор данных был сгенерирован из другой более тщательно продуманной модели аэродинамического тела.
Затем, мы создаем объект IDDATA представлять и хранить данные:
z = iddata(y, u, 0.02, 'Name', 'Data'); z.InputName = {'Aileron angle' 'Elevator angle' ... 'Rudder angle' 'Dynamic pressure' ... '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'}; z.InputUnit = {'rad' 'rad' 'rad' 'kg/(m*s^2)' 'm/s' ... 'rad/s' 'rad/s' 'rad/s' 'rad' 'rad'}; z.OutputName = {'V(x)' ... % Angular velocity around x-axis 'V(y)' ... % Angular velocity around y-axis 'V(z)' ... % Angular velocity around z-axis 'Accel.(y)' ... % Acceleration in y-direction 'Accel.(z)' ... % Acceleration in z-direction }; z.OutputUnit = {'rad/s' 'rad/s' 'rad/s' 'm/s^2' 'm/s^2'}; z.Tstart = 0; z.TimeUnit = 's';
Просмотрите входные данные:
figure('Name', [z.Name ': input data'],... 'DefaultAxesTitleFontSizeMultiplier',1,... 'DefaultAxesTitleFontWeight','normal',... 'Position',[50, 50, 850, 620]); for i = 1:z.Nu subplot(z.Nu/2, 2, i); plot(z.SamplingInstants, z.InputData(:,i)); title(['Input #' num2str(i) ': ' z.InputName{i}],'FontWeight','normal'); xlabel(''); axis tight; if (i > z.Nu-2) xlabel([z.Domain ' (' z.TimeUnit ')']); end end
Рисунок 1: Входные сигналы.
Просмотрите выходные данные:
figure('Name', [z.Name ': output data']); h_gcf = gcf; Pos = h_gcf.Position; h_gcf.Position = [Pos(1), Pos(2)-Pos(4)/2, Pos(3), Pos(4)*1.5]; for i = 1:z.Ny subplot(z.Ny, 1, i); plot(z.SamplingInstants, z.OutputData(:,i)); title(['Output #' num2str(i) ': ' z.OutputName{i}]); xlabel(''); axis tight; end xlabel([z.Domain ' (' z.TimeUnit ')']);
Рисунок 2: Выходные сигналы.
На первый взгляд может выглядеть странным измерить варианты некоторых выходных параметров во входном векторе. Однако модель, используемая для генерации данных, содержит несколько интеграторов, которые часто приводят к нестабильному поведению симуляции. Чтобы избежать этого, измерения некоторых выходных сигналов возвращены через нелинейного наблюдателя. Это вход 6 - 8 в наборе данных z. Таким образом это - система, действующая в замкнутом цикле, и цель осуществления моделирования состоит в том, чтобы предсказать "будущие" значения тех выходных параметров с помощью измерений текущего и прошлого поведения.
Моделирование системы
Чтобы смоделировать динамику интереса, мы используем объект модели IDNLGREY, чтобы представлять структуру пространства состояний системы. Разумная структура получена при помощи основной силы Ньютона и законов об импульсе (уравнения баланса). Чтобы полностью описать образцовую структуру, основные аэродинамические отношения (конститутивные отношения) также используются.
Файл MEX на C, aero_c.c, описывает систему с помощью состояния и выходных уравнений и начальных условий, как описано затем. Не использование деривации детализирует для уравнений движения, мы только отображаем конечное состояние и выводим уравнения и замечаем, что структура является довольно комплексной, а также нелинейной.
/* 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]; /* Body diameter. */
A = p[4]; /* Body reference area. */
I = p[5]; /* Moment of inertia, x-y-z. */
m = p[6]; /* 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]; /* Body reference area. */
m = p[6]; /* 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];
}
Мы должны также обеспечить начальные значения этих 23 параметров. Мы задаем некоторые параметры (аэродинамические коэффициенты силы, аэродинамические коэффициенты импульса, и момент факторов инерции) как векторы в 8 различных объектах параметра. Начальные значения параметров получены частично физическим обоснованием и частично количественным предположением. Последние 4 параметра (A, я, m и K) являются более или менее константами, таким образом, путем фиксации этих параметров мы получаем образцовую структуру с 16 свободными параметрами, распределенными среди объектов параметра F, M, C и 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. 'Body diameter' ... % d, scalar. 'Body reference area' ... % A, scalar. 'Moment of inertia, x-y-z' ... % I, 3-by-1 vector. '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; -200; -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 true true(3, 1) true true});
Мы также задаем 5 состояний образцовой структуры таким же образом:
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});
Образцовый файл наряду с порядком, параметром и данными о начальных состояниях теперь используется, чтобы создать объект IDNLGREY, описывающий систему:
FileName = 'aero_c'; % File describing the model structure. Order = [5 10 5]; % Model orders [ny nu nx]. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Model', 'TimeUnit', 's');
Затем, мы используем данные из объекта IDDATA задать сигналы ввода и вывода системы:
nlgr.InputName = z.InputName; nlgr.InputUnit = z.InputUnit; nlgr.OutputName = z.OutputName; nlgr.OutputUnit = z.OutputUnit;
Таким образом у нас есть объект IDNLGREY с 10 входными сигналами, 5 состояниями и 5 выходными сигналами. Как упомянуто ранее, модель также содержит 23 параметра, 7 из которых фиксируются и 16 из которых свободны:
nlgr
nlgr =
Continuous-time nonlinear grey-box model defined by 'aero_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 input(s), 5 state(s), 5 output(s), and 16 free parameter(s) (out of 23).
Name: Model
Status:
Created by direct construction or transformation. Not estimated.
Производительность первоначальной модели
Прежде, чем оценить 16 свободных параметров, мы моделируем систему с помощью начального вектора параметра. Симуляция предоставляет полезную информацию о качестве первоначальной модели:
clf
compare(z, nlgr); % simulate the model and compare the response to measured values
Рисунок 3: Сравнение между измеренными выходными параметрами и моделируемыми выходными параметрами первоначальной модели.
Из графика мы видим, что измеренные и моделируемые сигналы совпадают друг с другом тесно за исключением отрезка времени 4 - 6 секунд. Этот факт ясно показан в графике ошибок прогноза:
figure; h_gcf = gcf; Pos = h_gcf.Position; h_gcf.Position = [Pos(1), Pos(2)-Pos(4)/2, Pos(3), Pos(4)*1.5]; pe(z, nlgr);
Рисунок 4: ошибки Прогноза первоначальной модели.
Оценка параметра
Первоначальная модель, как задано выше является вероятной отправной точкой для оценки параметра. Затем, мы вычисляем ошибочную оценку прогноза 16 свободных параметров. Это вычисление займет время.
duration = clock; nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on')); duration = etime(clock, duration);
Производительность предполагаемой аэродинамической модели тела
На используемом компьютере оценка параметров заняла следующее количество времени, чтобы завершиться:
disp(['Estimation time : ' num2str(duration, 4) ' seconds']);
Estimation time : 18.29 seconds
disp(['Time per iteration: ' num2str(duration/nlgr.Report.Termination.Iterations, 4) ' seconds.']);
Time per iteration: 0.8709 seconds.
Чтобы оценить качество предполагаемой модели и проиллюстрировать улучшение по сравнению с первоначальной моделью, мы моделируем предполагаемую модель и сравниваем измеренные и моделируемые выходные параметры:
clf compare(z, nlgr);
Рисунок 5: Сравнение между измеренными выходными параметрами и моделируемыми выходными параметрами предполагаемой модели.
Фигура ясно указывает на улучшение по сравнению с результатом симуляции, полученным с первоначальной моделью. Системная динамика в отрезке времени 4 - 6 секунд теперь получена с намного более высокой точностью, чем прежде. Это лучше всего отображено путем рассмотрения ошибок прогноза:
figure; h_gcf = gcf; Pos = h_gcf.Position; h_gcf.Position = [Pos(1), Pos(2)-Pos(4)/2, Pos(3), Pos(4)*1.5]; pe(z, nlgr);
Рисунок 6: ошибки Прогноза предполагаемой модели.
Давайте завершим тематическое исследование путем отображения модели и оцененной неопределенности:
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'aero_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 input(s), 5 state(s), 5 output(s), and 16 free parameter(s) (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) 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 (fixed) in [-Inf, Inf]
x(2) Angular velocity around y-axis(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(3) Angular velocity around z-axis(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(4) Angle of attack(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(5) Angle of sideslip(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf]
Outputs:
y(1) V(x)(t) [rad/s]
y(2) V(y)(t) [rad/s]
y(3) V(z)(t) [rad/s]
y(4) Accel.(y)(t) [m/s^2]
y(5) Accel.(z)(t) [m/s^2]
Parameters: Value Standard Deviation
p1(1) Aerodynamic force coefficient [1/(rad*m^2..] 21.2863 0.339394 (estimated) in [-Inf, Inf]
p1(2) -7.62502 0.180264 (estimated) in [-Inf, Inf]
p1(3) 35.0799 0.657227 (estimated) in [-Inf, Inf]
p1(4) 8.58246 1.08611 (estimated) in [-Inf, Inf]
p2(1) Aerodynamic momentum coefficient [1/rad, 1/(..] -1.0476 0.0733533 (estimated) in [-Inf, Inf]
p2(2) 15.6854 0.883102 (estimated) in [-Inf, Inf]
p2(3) 3.00613 0.199227 (estimated) in [-Inf, Inf]
p2(4) -17.7963 0.324639 (estimated) in [-Inf, Inf]
p2(5) -1060.91 224.269 (estimated) in [-Inf, Inf]
p2(6) -53.5594 1.25436 (estimated) in [-Inf, Inf]
p2(7) 34.6095 1.37299 (estimated) in [-Inf, Inf]
p2(8) -210.237 7.95211 (estimated) in [-Inf, Inf]
p2(9) -2641.55 273.034 (estimated) in [-Inf, Inf]
p2(10) -33.6327 3.05742 (estimated) in [-Inf, Inf]
p2(11) -50.9269 1.64086 (estimated) in [-Inf, Inf]
p3 Aerodynamic compensation factor [1/(s*rad)] -0.640669 0.706338 (estimated) in [-Inf, Inf]
p4 Body diameter [m] 0.17 0 (fixed) in [-Inf, Inf]
p5 Body reference area [m^2] 0.0227 0 (fixed) in [-Inf, Inf]
p6(1) Moment of inertia, x-y-z [kg*m^2, kg..] 0.5 0 (fixed) in [-Inf, Inf]
p6(2) 110 0 (fixed) in [-Inf, Inf]
p6(3) 110 0 (fixed) in [-Inf, Inf]
p7 Mass [kg] 107 0 (fixed) in [-Inf, Inf]
p8 Feedback gain [-] 6 0 (fixed) in [-Inf, Inf]
Name: Model
Status:
Termination condition: Maximum number of iterations or number of function evaluations reached..
Number of iterations: 21, Number of function evaluations: 22
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Data".
Fit to estimation data: [52.93;94.91;91.4;96.07;98.84]%
FPE: 4.627e-10, MSE: 1.672
More information in model's "Report" property.
Заключительные замечания
Предполагаемая модель является хорошей отправной точкой для исследования основной производительности различных стратегий управления. Высокочастотные модели, которые предпочтительно имеют физическую интерпретацию, являются, например, жизненными компонентами так называемых "образцовых прогнозирующих систем управления".
Дополнительная информация
Для получения дополнительной информации об идентификации динамических систем с System Identification Toolbox™ см. страницу продукта System Identification Toolbox.