Объект

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

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

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

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

Диаграммы Боде 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)). Вы оцениваете коэффициенты полиномов A (q) и B (q). e (t) - белый шум.

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

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

Выбор шага расчета

Выбор шага расчета важен для хорошего приближения объекта непрерывного времени к модели дискретного времени. Правило thumb состоит в том, чтобы выбрать частоту дискретизации, как в 20 раз больше доминирующей динамики системы. Самая быстрая доминирующая динамика у объекта - 60 рад/с, или около 10Hz. Поэтому частота дискретизации 200Hz.

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) и A (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). Здесь используется фильтр lowpass Баттерворта 4-го порядка с 10Hz частоты отсечения. Это примерно частота самой быстрой доминирующей динамики на объекте (60рад/с). Низкочастотного фильтра здесь достаточно, потому что шумовой термин не имеет низкочастотного содержимого.

% 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));

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

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), захватывают отклонения оценок параметров. Квадратный корень дисперсий является стандартными отклонениями от оценок параметра. Первые три элемента диагоналей являются тремя параметрами, оцененными в A (q) полиноме. Последний элемент является одним параметром в полиноме B (q). Рассмотрим первый оцененный параметр в A (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. Эта модель захватила динамику объекта переключения режимов. Вы подтвердили предполагаемую модель, посмотрев на ошибку между оцененными, моделируемыми, измеренными выходами системы, а также ковариационными оценками параметра.