В этом примере показано, как оценить параметры нелинейной модели серого прямоугольника с использованием нескольких данных эксперимента. В качестве основы для обсуждения будет использована система, демонстрирующая сухое трение между двумя твердыми телами. В этой системе одно тело фиксируется, в то время как другое тело перемещается вперед и назад над неподвижным телом под действием экзогенной силы согласно фиг.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, отражающего ситуацию моделирования. Также добавляется информация о наименованиях и единицах ввода, вывода, состояниях и параметрах модели структуры модели. Обратите внимание, что значения Parameters и InitityStates здесь указаны как векторы, что по умолчанию означает, что при вызове метода 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 и так далее.
Дополнительные сведения об идентификации динамических систем с помощью Toolbox™ System Identification см. на информационной странице System Identification Toolbox.