Завод

Завод может быть представлен следующим образом:

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

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

Установка работает в первом режиме до 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 имеют вид:

Здесь $${}^{\mathrm{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, или если шум не белый. К решению этой проблемы относятся два подхода:

Выбор времени выборки

Выбор времени выборки важен для хорошей аппроксимации установки непрерывного времени с помощью дискретной модели времени. Правило большого пальца состоит в выборе частоты выборки, в 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) для $${}_{\mathrm{получения}}uf(t)=F$$(q) u $${}_{(}t)и\mathrm{yf}(t$$) = F (q) y (t). Затем используйте $${}_{}$$отфильтрованные $${}_{}\mathrm{сигналы}$$uf (t) и yf (t) в оценщике вместо входного сигнала u (t) установки и выходного сигнала y (t). Выбор фильтра данных позволяет уменьшить влияние e (t) на оценку.

Фильтр данных F (q) обычно является фильтром нижних частот или полосовым фильтром на основе частотного диапазона, важного для применения, и характеристик e (t). Здесь используется фильтр нижних частот Баттерворта 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);

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

Для оценки параметров в режиме реального времени используется команда recurexedARX. Команда создает системный 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 = 5 с. Большая начальная ошибка может быть уменьшена путем предоставления оценщику начального предположения для параметров модели и ковариации начального параметра. Когда завод переключается во второй режим, ошибки сначала растут, но с течением времени также уменьшаются. Это дает уверенность в том, что оцененные параметры модели хорошо улавливают поведение модели для данного входного сигнала.

Резюме

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