Моделирование аэродинамического тела
Этот пример показывает серое прямоугольное моделирование большой и сложной нелинейной системы. Цель состоит в том, чтобы показать способность использовать модель 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, чтобы представлять структуру пространства состояний системы. Разумная структура получается при помощи основных законов силы и импульса Ньютона (уравнения баланса). Для полного описания структуры модели используются также основные аэродинамические зависимости (конститутивные зависимости).
aero_c.c С файл MEX на 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, I, 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.4 seconds
disp(['Time per iteration: ' num2str(duration/nlgr.Report.Termination.Iterations, 4) ' seconds.']);
Time per iteration: 0.876 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.