В этом примере показано, как оценить параметры нелинейной серой модели поля, использующей несколько данных об эксперименте. Система, показывающая сухое трение между двумя твердыми телами, будет использоваться в качестве основания для обсуждения. В этой системе фиксируется одно тело, в то время как другое тело двигается вперед и назад по фиксированному телу из-за внешней силы согласно рисунку 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');True Estimated parameter vector
ptrue = [400; 2e3; 0.0001; 1700];
fprintf(' %10.5f %10.5f\n', [ptrue'; getpvec(nlgr)']); 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.