В этом примере показано серое моделирование статической системы с одним входом и одним выходом с использованием функции MATLAB в качестве файла ODE.
Идентификация системы обычно связана с идентификацией параметров динамических моделей. Однако статические модели также представляют интерес, иногда сами по себе, а иногда в качестве субмоделей более крупных более вовлеченных моделей. Пример последнего рассматривается в тематическом исследовании «An Industrial Robot Arm» (idnlgreydemo13.m), где модель статического трения используется в качестве фиксированного (предварительно оцененного) компонента модели руки робота.
Мы обсуждаем, как представить явление статического трения как модель IDNLGREY и оценить соответствующие коэффициенты.
Модели прерывистого и кусочно-непрерывного трения часто проблематичны для высокопроизводительных контроллеров непрерывного действия. Этот факт побудил Маккара, Диксона, Сойера и Ху предложить новую непрерывно и дифференцированную модель трения, которая фиксирует наиболее распространенные явления трения, встречающиеся на практике. Новая структура модели трения была представлена в
C. Makkar, W. E. Dixon, W. G. Sawyer, and G.Hu "A New Continuously
Differentiable Friction Model for Control Systems Design", IEEE(R)/ASME
International Conference on Advanced Intelligent Mechatronics,
Monterey, CA, 2005, pages 600-605.
и послужит основой для наших экспериментов по статической идентификации.
Модель трения, предложенная Маккаром и др., связывает скорость скольжения v (t) тела, находящегося в контакте с другим телом, с силой трения f (t) через статическое соотношение.
f(t) = g(1)*(tanh(g(2)*v(t) - tanh(g(3)*v(t))
+ g(4)*tanh(g(5)*v(t)) + g(6)*v(t)
где g (1),..., g (6) - 6 неизвестных положительных параметров. Эта структура модели отображает ряд хороших свойств, возникающих в реальных приложениях:
Модель трения симметрична относительно начала координат
Коэффициент статического трения аппроксимируется g (1) + g (4)
Первый член уравнения, tanh (g (2) * v (t) - tanh (g (3) * v (t), фиксирует так называемый эффект Стрибека, где член трения показывает довольно внезапное падение с увеличением скорости скольжения вблизи начала координат.
Куломбный эффект трения моделируется термином g (4) * tanh (g (5) * v (t))
Рассеяние вязкого трения отражается последним членом, g (6) * v (t)
Дополнительные сведения о трении в целом и предлагаемой структуре модели в частности см. в вышеупомянутом документе.
Теперь создадим объект модели IDNLGREY, описывающий статическое трение. Как обычно, отправной точкой является запись файла моделирования IDNLGREY, и здесь мы создаем файл MATLAB, friction_m.m, с содержанием следующим.
function [dx, f] = friction_m(t, x, v, g, varargin)
%FRICTION_M Nonlinear friction model with Stribeck, Coulomb and viscous
% dissipation effects.
% Output equation.
f = g(1)*(tanh(g(2)*v(1))-tanh(g(3)*v(1))) ... % Stribeck effect.
+ g(4)*tanh(g(5)*v(1)) ... % Coulomb effect.
+ g(6)*v(1); % Viscous dissipation term.
% Static system; no states.
dx = [];
Обратите внимание, что файл модели всегда должен возвращать обновление состояния dx и что оно должно быть пустым ([]) в случаях статического моделирования.
Следующим шагом является передача файла модели, информации о порядке модели, предполагаемом векторе параметров и т.д. в качестве входных аргументов конструктору IDNLGREY. Мы также указываем имена и единицы ввода и вывода и заявляем, что все параметры модели должны быть положительными.
FileName = 'friction_m'; % File describing the model structure. Order = [1 1 0]; % Model orders [ny nu nx]. Parameters = {[0.20; 90; 11; ... 0.12; 110; 0.015]}; % Initial parameters. InitialStates = []; % Initial initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Static friction model', ... 'InputName', 'Slip speed', 'InputUnit', 'm/s', ... 'OutputName', 'Friction force', 'OutputUnit', 'N', ... 'TimeUnit', 's'); nlgr = setpar(nlgr, 'Minimum', {zeros(5, 1)}); % All parameters must be >= 0.
После этих действий мы имеем исходную модель трения со следующими свойствами.
present(nlgr);
nlgr =
Continuous-time nonlinear static model defined by 'friction_m' (MATLAB file):
y(t) = H(t, u(t), p1) + e(t)
with 1 input(s), 0 state(s), 1 output(s), and 6 free parameter(s) (out of 6).
Inputs:
u(1) Slip speed(t) [m/s]
Outputs:
y(1) Friction force(t) [N]
Parameters: Value
p1(1) p1 0.2 (estimated) in [0, Inf]
p1(2) 90 (estimated) in [0, Inf]
p1(3) 11 (estimated) in [0, Inf]
p1(4) 0.12 (estimated) in [0, Inf]
p1(5) 110 (estimated) in [0, Inf]
p1(6) 0.015 (estimated) in [0, Inf]
Name: Static friction model
Status:
Created by direct construction or transformation. Not estimated.
More information in model's "Report" property.
В наших экспериментах по идентификации мы заинтересованы не только в модели полного трения, но и в изучении того, как модель пониженного трения будет работать. Под уменьшенной мы здесь подразумеваем модель трения, которая содержит два из трех членов полной модели. Чтобы исследовать это, создаются три копии полной структуры модели, и в каждой копии мы фиксируем вектор параметров так, что только два члена будут способствовать:
nlgr1 = nlgr; nlgr1.Name = 'No Striebeck term'; nlgr1 = setpar(nlgr1, 'Value', {[zeros(3, 1); Parameters{1}(4:6)]}); nlgr1 = setpar(nlgr1, 'Fixed', {[true(3, 1); false(3, 1)]}); nlgr2 = nlgr; nlgr2.Name = 'No Coulombic term'; nlgr2 = setpar(nlgr2, 'Value', {[Parameters{1}(1:3); 0; 0; Parameters{1}(6)]}); nlgr2 = setpar(nlgr2, 'Fixed', {[false(3, 1); true(2, 1); false]}); nlgr3 = nlgr; nlgr3.Name = 'No dissipation term'; nlgr3 = setpar(nlgr3, 'Value', {[Parameters{1}(1:5); 0]}); nlgr3 = setpar(nlgr3, 'Fixed', {[false(5, 1); true]});
В нашем распоряжении 2 различных (смоделированных) набора данных, в которых входная скорость проскальзывания была увеличена от -10 м/с до 10 м/с по типу клина. Мы загружаем данные и создаем два объекта IDDATA для наших экспериментов идентификации, ze для оценки и zv для целей проверки.
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'frictiondata')); ze = iddata(f1, v, 1, 'Name', 'Static friction system'); set(ze, 'InputName', 'Slip speed', 'InputUnit', 'm/s', ... 'OutputName', 'Friction force', 'OutputUnit', 'N', ... 'Tstart', 0, 'TimeUnit', 's'); zv = iddata(f2, v, 1, 'Name', 'Static friction system'); set(zv, 'InputName', 'Slip speed', 'InputUnit', 'm/s', ... 'OutputName', 'Friction force', 'OutputUnit', 'N', ... 'Tstart', 0, 'TimeUnit', 's');
Данные ввода-вывода, которые будут использоваться для оценки, отображаются в окне графика.
figure('Name', ze.Name);
plot(ze);
Рис. 1: Данные ввода-вывода из системы, характеризующейся трением.
Что касается доступа к данным ввода-вывода и четырем различным исходным моделям трения, то теперь очевидный вопрос заключается в том, насколько хороши эти модели? Рассмотрим это для набора оценочных данных с помощью моделирования, выполняемого COMPARE:
set(gcf,'DefaultLegendLocation','southeast'); compare(ze, nlgr, nlgr1, nlgr2, nlgr3);
Рис. 2: Сравнение истинного выхода с моделируемыми выходами четырех исходных моделей трения.
Ни одна из начальных моделей не способна правильно описать истинные выходные данные. Для преодоления этого мы оцениваем параметры модели всех четырех структур модели. Мы настраиваем все оценки для выполнения не более 30 итераций и для остановки итераций только в том случае, если допуск почти равен нулю (что на практике никогда не будет для реальных данных). Эти вычисления займут некоторое время.
opt = nlgreyestOptions('Display', 'on'); opt.SearchOptions.MaxIterations = 30; opt.SearchOptions.FunctionTolerance = eps; opt.EstimateCovariance = false; nlgr = nlgreyest(nlgr, ze, opt); nlgr1 = nlgreyest(nlgr1, ze, opt); nlgr2 = nlgreyest(nlgr2, ze, opt); nlgr3 = nlgreyest(nlgr3, ze, opt);
Производительность моделей еще раз исследуют путем сравнения истинных выходных данных с моделируемыми выходными данными четырех моделей, полученными с помощью команды СРАВНИТЬ, но на этот раз сравнение основано на наборе проверочных данных zv.
compare(zv, nlgr, nlgr1, nlgr2, nlgr3);
Рис. 3: Сравнение истинного выхода с моделируемыми выходами четырех расчетных моделей трения.
Для этой системы мы ясно видим, что полная модель превосходит уменьшенные. Тем не менее, уменьшенные модели, по-видимому, способны довольно хорошо фиксировать эффекты, которые они моделируют, и в каждом случае оценка приводит к гораздо лучшей подгонке. Наихудшая посадка получается для модели, где термин вязкого рассеяния был исключен. Впечатляющая подгонка полной модели не является большим сюрпризом, поскольку ее структура совпадает с той, которая используется для генерации истинных выходных данных. Параметры полной модели также близки к параметрам, которые использовались для генерации истинных выходных данных модели:
disp(' True Estimated parameter vector');True Estimated parameter vector
ptrue = [0.25; 100; 10; 0.1; 100; 0.01];
fprintf(' %7.3f %7.3f\n', [ptrue'; getpvec(nlgr)']); 0.250 0.249
100.000 106.637
10.000 9.978
0.100 0.100
100.000 87.699
0.010 0.010
Критерий Final Prediction Error (FPE) (низкие значения хорошие), применяемый ко всем четырем моделям трения, подтверждает превосходство модели полного трения:
fpe(nlgr, nlgr1, nlgr2, nlgr3);
1.0e-03 *
0.0002 0.1665 0.1584 0.7931
Что касается динамических систем, мы также можем изучить ошибки прогнозирования статической модели с помощью PE. Мы делаем это для модели полного трения и делаем вывод, что остатки, по-видимому, имеют случайную природу:
pe(zv, nlgr);
Рис. 4: Ошибки прогнозирования, полученные при расчетной модели полного трения.
Далее мы проверяем случайность, просматривая остатки («остатки») модели полного трения:
figure('Name', [nlgr.Name ': residuals of estimated IDNLGREY model']); resid(zv, nlgr);
Рис. 5: Остатки, полученные с расчетной моделью полного трения.
Реакция статических моделей на шаг также может быть вычислена и нанесена на график. Применим единичный шаг и сделаем это для всех четырех расчетных моделей трения:
figure('Name', [nlgr.Name ': step responses of estimated models']); step(nlgr, nlgr1, nlgr2, nlgr3); legend(nlgr.Name, nlgr1.Name, nlgr2.Name, nlgr3.Name, 'location', 'SouthEast');
Рис. 6: Единичные ступенчатые характеристики четырех расчетных моделей трения.
Наконец, мы отображаем ряд свойств, таких как расчетные стандартные отклонения параметров, функция потерь и т.д. модели полного трения.
present(nlgr);
nlgr =
Continuous-time nonlinear static model defined by 'friction_m' (MATLAB file):
y(t) = H(t, u(t), p1) + e(t)
with 1 input(s), 0 state(s), 1 output(s), and 6 free parameter(s) (out of 6).
Inputs:
u(1) Slip speed(t) [m/s]
Outputs:
y(1) Friction force(t) [N]
Parameters: Value
p1(1) p1 0.249402 (estimated) in [0, Inf]
p1(2) 106.637 (estimated) in [0, Inf]
p1(3) 9.97835 (estimated) in [0, Inf]
p1(4) 0.0999916 (estimated) in [0, Inf]
p1(5) 87.6992 (estimated) in [0, Inf]
p1(6) 0.0100019 (estimated) in [0, Inf]
Name: Static friction model
Status:
Termination condition: Change in parameters was less than the specified tolerance..
Number of iterations: 9, Number of function evaluations: 10
Estimated using Solver: FixedStepDiscrete; Search: lsqnonlin on time domain data "Static friction system".
Fit to estimation data: 99.68%
FPE: 2.428e-07, MSE: 2.414e-07
More information in model's "Report" property.
В этом примере показано, как выполнять IDNLGREY-моделирование статической системы. Процедура выполнения этой операции в основном аналогична процедуре моделирования динамических систем.