Сухое трение между двумя телами: оценка параметра Используя несколько данных об эксперименте
Этот пример показывает, как оценить параметры нелинейной серой модели поля, использующей несколько данных об эксперименте. Система, показывающая сухое трение между двумя твердыми телами, будет использоваться в качестве основания для обсуждения. В этой системе фиксируется одно тело, в то время как другое тело двигается вперед и назад по фиксированному телу из-за внешней силы согласно рисунку 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 является неизвестным скользящим параметром трения. Сухое трение, с другой стороны, является довольно комплексным явлением. В газете:
А. Клэвель, М. Сорайн и Ц. Чжан. "Моделирование и идентификация системы листовой рессоры". В третьем Семинаре IFAC по Усовершенствованиям в Автомобильном Управлении, 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'});
Данные ввода - вывода
На данном этапе мы загружаем доступные (моделируемые) данные ввода - вывода. Файл содержит данные из трех различных (моделируемых) тестовых прогонов каждое содержание, 1 000 поврежденных шумом выборок ввода - вывода сгенерировали использование уровня выборки (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'); ptrue = [400; 2e3; 0.0001; 1700]; fprintf(' %10.5f %10.5f\n', [ptrue'; getpvec(nlgr)']);
True Estimated parameter vector
400.00000 402.11783
2000.00000 1986.02824
0.00010 0.00011
1700.00000 1705.29327
Путем итогового использования команды 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.118 (estimated) in [-Inf, Inf]
p2 Sliding friction force coefficient [kg/s] 1986.03 (estimated) in [-Inf, Inf]
p3 Distance-related dry friction parameter [m] 0.000105164 (estimated) in [-Inf, Inf]
p4 Force-related dry friction parameter [N] 1705.29 (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.