В этом примере показан серый ящик статической системы с одним входом и одним выходом с использованием функции MATLAB в качестве файла ODE.
Система идентификации обычно касается идентификации параметров динамических моделей. Однако статические модели также представляют интерес, иногда сами по себе, а иногда как подмодели больших более вовлеченных моделей. Пример последнего обсуждается в тематическом исследовании «Промышленная рука робота» (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.
и послужит базисом для наших экспериментов по статической идентификации.
Модель трения, предложенная Makkar, et al, связывает скорость скольжения 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);
Эффективность моделей снова исследуется путем сравнения истинного результата с моделируемыми выходами четырех моделей, полученными с помощью COMPARE, но на этот раз сравнение основано на наборе данных валидации 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 статической системы. Процедура для этого в основном та же, что и для моделирования динамических систем.