Этот пример показывает, как оценить параметры нелинейной модели серого ящика с помощью нескольких данных эксперимента. В качестве базиса для обсуждения будет использована система сухого трения между двумя твердыми телами. В этой системе одно тело фиксируется, в то время как другое тело движется вперед и назад по неподвижному телу от экзогенной силы согласно фиг.1.
Фигура 1: Схематический вид системы двух тел.
Используя третий закон движения Ньютона, движение движущееся тело описывается:
F_tot(t) = m*a(t) = m*d/dt v(t) = m*d^2/dt^2 s(t)
где F_tot (t) равна экзогенной силе F (t) минус сила трения, вызванная контактом между двумя телами. Сила трения принята как сумма силы трения скольжения F_s (t) и силы трения сухого F_d (t). Первый обычно моделируется как линейная функция скорости, т.е. F_s (t) = k * v (t), где k - неизвестный параметр скользящего трения. Сухое трение, с другой стороны, - довольно сложное явление. В статье:
А. Клавель, М. Сорине и К. Чжан. «Моделирование и идентификация листовой пружинной системы». На третьем семинаре ИФАК по Усовершенствования в области автомобильного контроля, 2001 год.
оно моделируется обыкновенным дифференциальным уравнением:
d/dt F_d(t) = -1/e*|v(t)|*F_d(t) + f/e*v(t)
где e и f - два неизвестных параметра с размерностями расстояние и сила, соответственно.
Обозначая входной сигнал u (t) = F (t) [N], вводя состояния как:
x1(t) = s(t) Position of the moving body [m]. x2(t) = v(t) Velocity of the moving body [m/s]. x3(t) = F_d(t) Dry friction force between the bodies [N].
и моделируйте параметры как:
m Mass of the moving body [m]. k Sliding friction force coefficient [kg/s]. e Distance-related dry friction parameter [m]. f Force-related dry friction parameter [N].
мы получаем следующую структуру модели пространства состояний:
d/dt x1(t) = x2(t) d/dt x2(t) = 1/m*(u(t) - k*x2(t) - x3(t)) d/dt x3(t) = 1/e*(-|x2(t)|*x3(t) + f*x2(t))
y(t) = x1(t)
Эти уравнения вводятся в файл модели C-MEX, twobodies_c.c. Его состояние и выходные уравнения обновления, compute_dx и compute_y, следующие:
/* State equations. */
void compute_dx(double *dx, double t, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* Retrieve model parameters. */
double *m, *k, *e, *f;
m = p[0]; /* Mass of the moving body. */
k = p[1]; /* Sliding friction force coefficient. */
e = p[2]; /* Distance-related dry friction parameter. */
f = p[3]; /* Force-related dry friction parameter. */
/* x[0]: Position. */
/* x[1]: Velocity. */
/* x[2]: Dry friction force. */
dx[0] = x[1];
dx[1] = (u[0]-k[0]*x[1]-x[2])/m[0];
dx[2] = (-fabs(x[1])*x[2]+f[0]*x[1])/e[0];
}
/* Output equation. */
void compute_y(double *y, double t, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* y[0]: Position. */
y[0] = x[0];
}
Написав файл, описывающий структуру модели, следующим шагом является создание объекта IDNLGREY, отражающего ситуацию моделирования. Также добавляется информация об именах и модулях входов, выходах, состояниях и параметрах модели структуры модели. Заметьте, что Параметры и InitialStates здесь заданы как векторы, что по умолчанию означает, что все параметры модели и ни один вектор начального состояния не будут оценены при вызове NLGREYEST.
FileName = 'twobodies_c'; % File describing the model structure. Order = [1 1 3]; % Model orders [ny nu nx]. Parameters = [380; 2200; 0.00012; 1900]; % Initial parameter vector. InitialStates = [0; 0; 0]; % Initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Two body system', ... 'InputName', 'Exogenous force', ... 'InputUnit', 'N', ... 'OutputName', 'Position of moving body', ... 'OutputUnit', 'm', ... 'TimeUnit', 's'); nlgr = setinit(nlgr, 'Name', {'Position of moving body' ... 'Velocity of moving body' ... 'Dry friction force between the bodies'}); nlgr = setinit(nlgr, 'Unit', {'m' 'm/s' 'N'}); nlgr = setpar(nlgr, 'Name', {'Mass of the moving body' ... 'Sliding friction force coefficient', ... 'Distance-related dry friction parameter' ... 'Force-related dry friction parameter'}); nlgr = setpar(nlgr, 'Unit', {'m' 'kg/s' 'm' 'N'});
На данной точке мы загружаем доступные (моделируемые) данные ввода-вывода. Файл содержит данные трех различных (моделируемых) тестовых запусков, каждый из которых содержит 1000 поврежденных шумом входных и выходных выборок, сгенерированных с использованием частоты дискретизации (Ts) 0,005 секунд. Вход u (t) является экзогенной силой [N], действующей на движущееся тело. В экспериментах вход был симметричным пилообразным формируемым сигналом, где частота повторения формы волны была одинаковой для всех экспериментов, но где максимальная амплитуда сигнала изменялась в период между тестовыми запусками. Выходом y (t) является положение [m] подвижного тела (относительно фиксированного). В наших целях моделирования мы создадим один объект IDDATA, содержащий три различных эксперимента:
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'twobodiesdata')); z = merge(iddata(y1, u1, 0.005), iddata(y2, u2, 0.005), iddata(y3, u3, 0.005)); z.Name = 'Two body system'; z.InputName = nlgr.InputName; z.InputUnit = nlgr.InputUnit; z.OutputName = nlgr.OutputName; z.OutputUnit = nlgr.OutputUnit; z.Tstart = 0; z.TimeUnit = nlgr.TimeUnit;
Данные ввода-вывода, используемые для экспериментов по идентификации вперед, показаны в окне графика.
figure('Name', [z.Name ': input-output data']); for i = 1:z.Ne zi = getexp(z, i); subplot(z.Ne, 2, 2*i-1); % Input. plot(zi.SamplingInstants, zi.InputData); title([z.ExperimentName{i} ': ' zi.InputName{1}]); if (i < z.Ne) xlabel(''); else xlabel([z.Domain ' (' zi.TimeUnit ')']); end axis('tight'); subplot(z.Ne, 2, 2*i); % Output. plot(zi.SamplingInstants, zi.OutputData); title([z.ExperimentName{i} ': ' zi.OutputName{1}]); if (i < z.Ne) xlabel(''); else xlabel([z.Domain ' (' zi.TimeUnit ')']); end axis('tight'); end
Фигура 2. Входные-выходные данные из системы с двумя телами.
Перед оценкой четырех неизвестных параметров мы моделируем систему с помощью начальных значений параметров. Мы используем решатель для дифференциальных уравнений по умолчанию (ode45) с опциями решателя по умолчанию. При вызове только с двумя входными параметрами COMPARE будет оценивать полные векторы начального состояния, в этом случае один на эксперимент, т.е. 3 эксперимента каждый с вектором состояния 3 на 1 подразумевает в общей сложности 9 предполагаемых начальных состояний. Моделируемые и истинные выходы показаны в окне графика, и, как видно, подгонка достойна, но не так хороша, как желательно.
clf compare(z, nlgr);
Фигура 3: Сравнение истинных выходов и моделируемых выходов исходной модели двух тел.
В порядок улучшения подгонки далее оцениваются четыре параметра. Мы используем данные первого и последнего экспериментов на фазе оценки и храним данные второго эксперимента в чистых целях валидации.
opt = nlgreyestOptions('Display', 'on'); nlgr = nlgreyest(getexp(z, [1 3]), nlgr, opt);
В порядок исследования эффективности предполагаемой модели, моделирование ее наконец выполняется. Путем настройки массива начальной структуры состояния nlgr можно полностью определить, какие состояния оценивать в каждом эксперименте, например, COMPARE. Давайте здесь определим и используем структуру, где начальные состояния x1 (0) и x2 (0) оценены для эксперимента 1, x2 (0) для эксперимента 2 и x3 (0) для эксперимента 3. С помощью этой модификации в окне графика показано сравнение между выходами измеренных и моделей.
nlgr.InitialStates = struct('Name', getinit(nlgr, 'Name'), ... 'Unit', getinit(nlgr, 'Unit'), ... 'Value' , zeros(1, 3), 'Minimum', -Inf(1, 3), ... 'Maximum', Inf(1, 3), 'Fixed', ... {[false false true]; [true false true]; [true true false]}); compare(z, nlgr, compareOptions('InitialCondition', 'model'));
Фигура 4: Сравнение истинных выходов и моделируемых выходов предполагаемой модели двух тел.
Особый интерес представляет результат с данными второго эксперимента, которые не использовались для оценки параметра. Динамика истинной системы четко смоделирована достаточно хорошо для всех экспериментов. Предполагаемые параметры также довольно близки к таковым, используемой для генерации экспериментальных данных:
disp(' True Estimated parameter vector');True Estimated parameter vector
ptrue = [400; 2e3; 0.0001; 1700];
fprintf(' %10.5f %10.5f\n', [ptrue'; getpvec(nlgr)']); 400.00000 402.16624
2000.00000 1987.96042
0.00010 0.00010
1700.00000 1704.51209
Наконец, используя команду PRESENT, мы можем получить дополнительную информацию о предполагаемой модели:
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'twobodies_c' (MEX-file):
dx/dt = F(t, u(t), x(t), p1, ..., p4)
y(t) = H(t, u(t), x(t), p1, ..., p4) + e(t)
with 1 input(s), 3 state(s), 1 output(s), and 4 free parameter(s) (out of 4).
Inputs:
u(1) Exogenous force(t) [N]
States: Initial value
x(1) Position of moving body(t) [m] xinit@exp1 0 (estimated) in [-Inf, Inf]
xinit@exp2 0 (estimated) in [-Inf, Inf]
xinit@exp3 0 (fixed) in [-Inf, Inf]
x(2) Velocity of moving body(t) [m/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
xinit@exp2 0 (estimated) in [-Inf, Inf]
xinit@exp3 0 (fixed) in [-Inf, Inf]
x(3) Dry friction force between the bodies(t) [N] xinit@exp1 0 (fixed) in [-Inf, Inf]
xinit@exp2 0 (fixed) in [-Inf, Inf]
xinit@exp3 0 (estimated) in [-Inf, Inf]
Outputs:
y(1) Position of moving body(t) [m]
Parameters: Value
p1 Mass of the moving body [m] 402.166 (estimated) in [-Inf, Inf]
p2 Sliding friction force coefficient [kg/s] 1987.96 (estimated) in [-Inf, Inf]
p3 Distance-related dry friction parameter [m] 0.000104534 (estimated) in [-Inf, Inf]
p4 Force-related dry friction parameter [N] 1704.51 (estimated) in [-Inf, Inf]
Name: Two body system
Status:
Model modified after estimation.
More information in model's "Report" property.
Этот пример описал, как использовать несколько данных эксперимента при выполнении моделирования IDNLGREY. Можно использовать любое количество экспериментов, и для каждого такого эксперимента можно полностью определить, какое начальное состояние или состояния оценивать в NLGREYEST, COMPARE, PREDICT и так далее.
Для получения дополнительной информации об идентификации динамических систем с помощью System Identification Toolbox™ посетите информационную страницу System Identification Toolbox.