В этом примере показано серая модель большой и сложной нелинейной системы. Цель состоит в том, чтобы показать возможность использования модели 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 для представления структуры пространства состояний системы. Разумная структура получается при использовании основных законов силы и импульса Ньютона (уравнений баланса). Для полного описания структуры модели используются также основные аэродинамические соотношения (конститутивные соотношения).
C MEX-файл, 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, 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.
Предполагаемая модель является хорошей отправной точкой для изучения фундаментальных показателей различных стратегий контроля. Модели высокой точности, которые предпочтительно имеют физическую интерпретацию, являются, например, жизненно важными компонентами так называемых «модельных прогностических систем управления».
Дополнительные сведения об идентификации динамических систем с помощью Toolbox™ идентификации систем см. на странице инструментария идентификации систем.