Объект

Объект может быть представлен как:

Здесь, G является передаточной функцией, и e является белым шумом. Объект имеет два рабочих режима. В первом рабочем режиме передаточная функция:

Слегка ослабленные полюса в G1 (s) имеют затухание 0.02 и собственная частота 30rad/s. Во втором рабочем режиме собственная частота этих полюсов является 60rad/s. В этом режиме передаточная функция:

Объект действует в первом режиме до t=10s, и затем переключается на второй режим.

Диаграммы Боде G1 и G2:

wn = 30; % natural frequency of the lightly damped poles
ksi = 0.02; % damping of the poles
G1 = tf(1,conv([1/5 1],[1/wn^2 2*ksi/wn 1])); % plant in mode 1
wn = 60; % natural frequency in the second operating mode
G2 = tf(1,conv([1/5 1],[1/wn^2 2*ksi/wn 1])); % plant in mode 2
bode(G1,G2,{1,125});
legend('G1','G2','Location','Best');

Онлайновая оценка параметра ARX

Цель состоит в том, чтобы оценить динамику объекта во время его операции. ARX является общей структурой модели, используемой с этой целью. Модели ARX имеют форму:

Здесь, $$q^{-1}$$ оператор сдвига времени. Отношение полиномов B (q)/A (q) получает модель ввода - вывода (u (t) к y (t)), и 1/A (q) получает шумовую модель (e (t) к y (t)). Вы оцениваете коэффициенты (q) и B (q) полиномы. e (t) является белым шумом.

Структура модели ARX является хорошим первым кандидатом на оценку линейных моделей. Связанные методы оценки имеют низкую вычислительную нагрузку, численно устойчивы, и имеют свойство выпуклости. Свойство выпуклости гарантирует, что нет никакого риска оценки параметра, застревающей в локальные оптимумы. Однако структура модели ARX не обеспечивает гибкость для шумовых моделей.

Отсутствие гибкости в шумовом моделировании может представлять трудности, если структура объекта не соответствует со структурой модели ARX, или если шум не является белым. Два подхода, чтобы исправить эту проблему:

Выберите Sample Time

Выбор шага расчета важен для хорошего приближения объекта непрерывного времени моделью дискретного времени. Эмпирическое правило должно выбрать частоту дискретизации как 20 раз доминирующая динамика системы. Объект имеет самую быструю доминирующую динамику в 60rad/s, или приблизительно 10 Гц. Частота дискретизации - поэтому 200 Гц.

Ts = 0.005; % [s], Sample time, Ts=1/200

Системное возбуждение

Для успешной оценки входные параметры объекта должны постоянно волновать его динамику. Простые входные параметры, такие как один вход шага обычно не достаточны. В этом примере объект управляется импульсом с амплитудой 10 и период 1 секунды. Ширина импульса составляет 50% своего периода.

Сгенерируйте сигналы ввода и вывода объекта:

t = 0:Ts:20; % Time vector
u = double(rem(t/1,1)-0.5 < 0); % pulse
y = zeros(size(u));
% Store random number generator's states for reproducible results.
sRNG = rng;
rng('default');
% Simulate the mode-switching plant with a zero-order hold.
G1d = c2d(G1,Ts,'zoh');
B1 = G1d.num{1}.'; 
A1 = G1d.den{1}.'; % B1 and A1 corresponds to G1.
G2d = c2d(G2,Ts,'zoh');
B2 = G2d.num{1}.';
A2 = G2d.den{1}.'; % B2 and A2 corresponds to G2.
idx = numel(B2):-1:1;
for ct=(1+numel(B2)):numel(t)
    idx = idx + 1;
    if t(ct)<10 % switch mode after t=10s
        y(ct) = u(idx)*B1-y(idx(2:end))*A1(2:end);
    else
        y(ct) = u(idx)*B2-y(idx(2:end))*A2(2:end);
    end
end
% Add measurement noise
y = y + 0.02*randn(size(y));
% Restore the random number generator's states.
rng(sRNG);

Отобразите данные ввода - вывода на графике:

figure();
subplot(2,1,1);
plot(t,u);
ylabel('Input (u)');
subplot(2,1,2);
plot(t,y);
ylim([-0.2 1.2])
ylabel('Output (y)');
xlabel('Time [s]');

Отфильтруйте данные

Объект имеет форму:

где e (t) является белым шумом. В отличие от этого модели ARX имеют форму

Средство оценки будет использовать B (q) и (q), чтобы аппроксимировать G (q). Однако отметьте различие в шумовых моделях. Объект имеет белый шум e (t) непосредственно влияющий y (t), но модель ARX принимает, что термин белого шума, отфильтрованный 1/A (q), влияет на y (t). Это несоответствие будет негативно влиять на оценку.

Когда шумовая модель не представляет интереса, один метод, чтобы уменьшать удар этого несоответствия должен использовать фильтр данных. Используйте фильтр $$F(q)$$ и на u (t) и на y (t), чтобы получить $$u_f(t)=F(q)u(t)$$ и $$y_f(t)=F(q)y(t)$$. Затем используйте отфильтрованные сигналы $$u_f(t)$$ и $$y_f(t)$$ в средстве оценки вместо входа объекта u (t) и выход y (t). Выбор фильтра данных позволяет вам уменьшать влияние e (t) на оценке.

Данные фильтруют F (q), обычно lowpass или полосовой фильтр на основе важного частотного диапазона для приложения и характеристики e (t). Здесь, 4-й порядок Баттерворт фильтр lowpass с частотой среза 10 Гц используется. Это - приблизительно частота самой быстрой доминирующей динамики на объекте (60rad/s). Фильтр lowpass достаточен здесь, потому что шумовая часть не имеет низкочастотного содержимого.

% Filter coefficients
Fa = [1 -3.1806 3.8612 -2.1122 0.4383]; % denominator
Fb = [4.1660e-04 1.6664e-03 2.4996e-03 1.6664e-03 4.1660e-04]; % numerator
% Filter the plant input for estimation
uf = filter(Fb,Fa,u);
% Filter the plant output
yf = filter(Fb,Fa,y);

Настройте оценку

Используйте recursiveARX команду для онлайновой оценки параметра. Команда создает Систему object™ для онлайновой оценки параметра структуры модели ARX. Задайте следующие свойства объекта:

X = recursiveARX([3 1 0]); % [na nb nk]
X.ForgettingFactor = 0.995;

Создайте массивы, чтобы сохранить результаты оценки. Они полезны для проверки алгоритмов.

np = size(X.InitialParameterCovariance,1);
PHat = zeros(numel(u),np,np);
A = zeros(numel(u),numel(X.InitialA));
B = zeros(numel(u),numel(X.InitialB));
yHat = zeros(1,numel(u));

Используйте команду шага, чтобы обновить значения параметров с помощью одного набора данных ввода - вывода на каждом временном шаге. Это иллюстрирует онлайновую операцию средства оценки.

for ct=1:numel(t)
    % Use the filtered output and input signals in the estimator
    [A(ct,:),B(ct,:),yHat(ct)] = step(X,yf(ct),uf(ct));
    PHat(ct,:,:) = X.ParameterCovariance;
end

Просмотрите Диаграмму Боде предполагаемых передаточных функций:

G1Hat = idpoly(A(1000,:),B(1000,:),1,1,1,[],Ts); % Model snapshot at t=10s
G2Hat = idpoly(X); % Snapshot of the latest model, at t=20s
G2Hat.Ts = G1d.Ts; % Set the sample time of the snapshot
figure();
bode(G1,G1Hat);
xlim([0.5 120]);
legend('G1','Identified model at t=10s','Location','Best');

figure();
bode(G2,G2Hat);
xlim([0.5 120]);
legend('G2','Identified model at t=20s','Location','Best');

Проверка предполагаемой модели

Используйте следующие методы, чтобы подтвердить оценку параметра:

Абсолютная погрешность yf (t)-yhat (t) и относительная погрешность (yf (t)-yhat (t))/yf (t):

figure();
subplot(2,1,1);
plot(t,yf-yHat);
ylabel('Abs. Error');
subplot(2,1,2);
plot(t,(yf-yHat)./yf);
ylim([-0.05 0.05]);
ylabel('Rel. Error');
xlabel('Time [s]');

Абсолютные погрешности находятся порядка 1e-3, который мал по сравнению с самим измеренным выходным сигналом. График относительной погрешности в нижней части подтверждает это с ошибками, бывшими меньше 5% кроме в начале симуляции.

Диагонали ковариационной матрицы параметра, масштабируемой отклонением остаточных значений y (t)-yhat (t), получают отклонения оценок параметра. Квадратный корень отклонений является стандартными отклонениями оценок параметра. Первыми тремя элементами на диагоналях являются эти три параметра, оцененные в (q) полином. Последним элементом является один параметр в B (q) полином. Давайте посмотрим на первый предполагаемый параметр в (q)

noiseVariance = var(yf-yHat);
X.A(2) % The first estimated parameter. X.A(1) is fixed to 1

ans = -2.8635

sqrt(X.ParameterCovariance(1,1)*noiseVariance)

ans = 0.0175

Стандартное отклонение 0.0175 мало относительно абсолютного значения значения параметров 2.86. Это указывает на хорошее доверие к предполагаемому параметру.

figure();
plot(t,sqrt(PHat(:,1,1)*noiseVariance));
ylabel('Standard deviation estimate for the parameter A(2)')
xlabel('Time [s]');

Неопределенность мала и ограничена в течение оценки. Однако обратите внимание, что стандартные отклонения параметра являются также оценками. Они основаны на предположении, что остаточные значения y (t)-yhat (t) являются белыми. Это зависит от метода оценки, его связанных параметров, структуры предполагаемой модели и входного сигнала u. Различия между принятым и фактической структурой модели, отсутствием персистентного входного возбуждения или нереалистичных настроек метода оценки могут привести к чрезмерно оптимистическим или пессимистическим оценкам неопределенности.

Наконец, симулируйте предполагаемую модель ARX с помощью сохраненной истории предполагаемых параметров. Эта симуляция может также быть сделана одновременно с циклом оценки для валидации во время онлайновой операции.

ysim = zeros(size(y));
idx = numel(B2):-1:1;
for ct=(1+numel(B2)):numel(t)
    idx = idx + 1;
    ysim(ct) = u(idx(1))*B(idx(1),:)-ysim(idx(2:end))*A(ct,2:end)';
end
figure();
subplot(2,1,1);
plot(t,y,t,ysim);
ylabel('System Output');
legend('Measured','Estimated','Location','Best');
subplot(2,1,2);
plot(t,y-ysim);
ylim([-0.5 0.5]);
ylabel('Error, y(t)-ysim(t)');
xlabel('Time [s]');

Ошибка является большой первоначально, но она обосновывается к меньшему значению вокруг t=5s для первого рабочего режима. Большая начальная ошибка может уменьшаться путем обеспечения средству оценки исходного предположения для параметров модели и начальной ковариации параметра. Когда объект переключается на второй режим, ошибки растут первоначально, но успокаиваются со временем также. Это вселяет веру, что предполагаемые параметры модели способны получать поведение модели для данного входного сигнала.

Сводные данные

Вы выполнили онлайновую оценку параметра для модели ARX. Эта модель получила динамику объекта переключения режимов. Вы подтвердили предполагаемую модель путем рассмотрения ошибки между предполагаемой, симулированной, измеренной системой выходные параметры, а также оценки ковариации параметра.