Этот пример показывает, как оценить параметры нелинейной серой модели поля, использующей несколько данных об эксперименте. Система, показывающая сухое трение между двумя твердыми телами, будет использоваться в качестве основания для обсуждения. В этой системе фиксируется одно тело, в то время как другое тело двигается вперед и назад по фиксированному телу из-за внешней силы согласно рисунку 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.