Нелинейный System Identification

Этот пример показывает, как использовать команду anfis для нелинейной идентификации динамической системы.

Этот пример требует System Identification Toolbox™, когда сравнение сделано между нелинейным ANFIS и линейной моделью ARX.

Setup задач

Выйдите, если System Identification Toolbox не доступен.

if ~fuzzychecktoolboxinstalled('ident')
    errordlg('DRYDEMO needs the System Identification Toolbox.');
    return;
end

Набор данных для ANFIS и моделирования ARX был получен от преподавателя PT 326 Процесса названной Обратной связи лабораторного устройства, как описано в Главе 17 книги профессора Леннарта Лджанга "System Identification, Теория для Пользователя", Prentice Hall, 1987. Устройство функционирует как фен: воздух вентилируется через трубу и нагревается во входном отверстии. Температура воздуха измеряется термопарой при выходе. Вход u (k) является напряжением по сетке проводов резистора, чтобы нагреть входящий воздух; вывод y (k) является температурой воздуха выхода.

Вот результаты теста.

load drydemodata
data_n = length(y2);
output = y2;
input = [[0; y2(1:data_n-1)] ...
        [0; 0; y2(1:data_n-2)] ...
        [0; 0; 0; y2(1:data_n-3)] ...
        [0; 0; 0; 0; y2(1:data_n-4)] ...
        [0; u2(1:data_n-1)] ...
        [0; 0; u2(1:data_n-2)] ...
        [0; 0; 0; u2(1:data_n-3)] ...
        [0; 0; 0; 0; u2(1:data_n-4)] ...
        [0; 0; 0; 0; 0; u2(1:data_n-5)] ...
        [0; 0; 0; 0; 0; 0; u2(1:data_n-6)]];
data = [input output];
data(1:6, :) = [];
input_name = char('y(k-1)','y(k-2)','y(k-3)','y(k-4)','u(k-1)','u(k-2)','u(k-3)','u(k-4)','u(k-5)','u(k-6)');
index = 1:100;
subplot(2,1,1)
plot(index,y2(index),'-',index,y2(index),'o')
ylabel('y(k)','fontsize',10)
subplot(2,1,2)
plot(index,u2(index),'-',index,u2(index),'o')
ylabel('u(k)','fontsize',10)

Точки данных были собраны во время выборки 0,08 секунд. Одна тысяча точек данных ввода - вывода была собрана из процесса, когда вход u (k) был выбран, чтобы быть бинарной случайной переменой сигнала между 3.41 и 6,41 В. Вероятность сдвига входа на каждой выборке была 0.2. Набор данных доступен от System Identification Toolbox, и вышеупомянутые графики показывают выходную температуру y (k) и входное напряжение u (t) для первых 100 временных шагов.

Идентификация модели ARX

Условный метод состоит в том, чтобы удалить средние значения из данных и принять линейную модель формы:

y (k) +a1*y (k-1) +... +am*y (k-m) =b1*u (k-d) +... +bn*u (k-d-n+1)

где ай (i = 1 к m) и bj (j = 1 к n) линейные параметры, которые будут определены методами наименьших квадратов. Эта структура называется моделью ARX, и это точно задано тремя целыми числами [m, n, d]. Чтобы найти модель ARX для устройства сушилки, набор данных был разделен на обучение (k = 1 - 300) и проверка (k = 301 - 600) набор. Исчерпывающий поиск выполнялся, чтобы найти лучшую комбинацию [m, n, d], где каждое целое число позволено измененному от 1 до 10 независимо. Лучшая модель ARX, таким образом найденная, задана [m, n, d] = [5, 10, 2], с учебным RMSE 0,1122 и проверка RMSE 0,0749. Вышеупомянутые данные показывают подходящие результаты лучшей модели ARX.

trn_data_n = 300;
total_data_n = 600;
z = [y2 u2];
z = dtrend(z);
ave = mean(y2);
ze = z(1:trn_data_n,:);
zv = z(trn_data_n+1:total_data_n,:);
T = 0.08;

% Run through all different models
V = arxstruc(ze,zv,struc(1:10,1:10,1:10));
% Find the best model
nn = selstruc(V,0);
% Time domain plot
th = arx(ze,nn);
th.Ts = 0.08;
u = z(:,2);
y = z(:,1)+ave;
yp = sim(u,th)+ave;

xlbl = 'Time Steps';

subplot(2,1,1)
index = 1:trn_data_n;
plot(index, y(index), index, yp(index), '.')
rmse = norm(y(index)-yp(index))/sqrt(length(index));
title(sprintf(['(a) Training Data (Solid Line) and ARX Prediction (Dots)\nwith RMSE = ' num2str(rmse)]),'fontsize',10)
disp(['[na nb d] = ' num2str(nn)])
xlabel(xlbl,'fontsize',10)

subplot(2,1,2)
index = (trn_data_n+1):(total_data_n);
plot(index,y(index),index,yp(index),'.')
rmse = norm(y(index)-yp(index))/sqrt(length(index));
title(sprintf(['(b) Checking Data (Solid Line) and ARX Prediction (Dots)\nwith RMSE = ' num2str(rmse)]),'fontsize',10)
xlabel(xlbl,'fontsize',10)
[na nb d] = 5  10   2

Идентификация модели ANFIS

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

Чтобы использовать ANFIS для системы идентификации, первая вещь, которую мы должны сделать, выбрать вход. Таким образом, чтобы определить, какие переменные должны быть входными параметрами к модели ANFIS. Для простоты мы предполагаем, что существует 10 входных кандидатов (y (k-1), y (k-2), y (k-3), y (k-4), u (k-1), u (k-2), u (k-3), u (k-4), u (k-5), u (k-6)), и вывод, который будет предсказан, является y (k). Эвристический подход, чтобы ввести выбор называется последовательным прямым поиском, в котором каждый вход выбран последовательно, чтобы оптимизировать общую квадратичную невязку. Это может быть сделано функцией seqsrch; результат показывают в вышеупомянутом графике, где 3 входных параметров (y (k-1), u (k-3), и u (k-4)) выбраны с учебным RMSE 0,0609 и проверяющий RMSE 0,0604.

trn_data_n = 300;
trn_data = data(1:trn_data_n,:);
chk_data = data(trn_data_n+1:trn_data_n+300,:);
[~,elapsed_time] = seqsrch(3,trn_data,chk_data,input_name); % #ok<*ASGLU>
fprintf('\nElapsed time = %f\n',elapsed_time);
winH1 = gcf;
Selecting input 1 ...
ANFIS model 1: y(k-1) --> trn=0.2043, chk=0.1888
ANFIS model 2: y(k-2) --> trn=0.3819, chk=0.3541
ANFIS model 3: y(k-3) --> trn=0.5245, chk=0.4903
ANFIS model 4: y(k-4) --> trn=0.6308, chk=0.5977
ANFIS model 5: u(k-1) --> trn=0.8271, chk=0.8434
ANFIS model 6: u(k-2) --> trn=0.7976, chk=0.8087
ANFIS model 7: u(k-3) --> trn=0.7266, chk=0.7349
ANFIS model 8: u(k-4) --> trn=0.6215, chk=0.6346
ANFIS model 9: u(k-5) --> trn=0.5419, chk=0.5650
ANFIS model 10: u(k-6) --> trn=0.5304, chk=0.5601
Currently selected inputs: y(k-1)

Selecting input 2 ...
ANFIS model 11: y(k-1) y(k-2) --> trn=0.1085, chk=0.1024
ANFIS model 12: y(k-1) y(k-3) --> trn=0.1339, chk=0.1283
ANFIS model 13: y(k-1) y(k-4) --> trn=0.1542, chk=0.1461
ANFIS model 14: y(k-1) u(k-1) --> trn=0.1892, chk=0.1734
ANFIS model 15: y(k-1) u(k-2) --> trn=0.1663, chk=0.1574
ANFIS model 16: y(k-1) u(k-3) --> trn=0.1082, chk=0.1077
ANFIS model 17: y(k-1) u(k-4) --> trn=0.0925, chk=0.0948
ANFIS model 18: y(k-1) u(k-5) --> trn=0.1533, chk=0.1531
ANFIS model 19: y(k-1) u(k-6) --> trn=0.1952, chk=0.1853
Currently selected inputs: y(k-1) u(k-4)

Selecting input 3 ...
ANFIS model 20: y(k-1) u(k-4) y(k-2) --> trn=0.0808, chk=0.0822
ANFIS model 21: y(k-1) u(k-4) y(k-3) --> trn=0.0806, chk=0.0836
ANFIS model 22: y(k-1) u(k-4) y(k-4) --> trn=0.0817, chk=0.0855
ANFIS model 23: y(k-1) u(k-4) u(k-1) --> trn=0.0886, chk=0.0912
ANFIS model 24: y(k-1) u(k-4) u(k-2) --> trn=0.0835, chk=0.0843
ANFIS model 25: y(k-1) u(k-4) u(k-3) --> trn=0.0609, chk=0.0604
ANFIS model 26: y(k-1) u(k-4) u(k-5) --> trn=0.0848, chk=0.0867
ANFIS model 27: y(k-1) u(k-4) u(k-6) --> trn=0.0890, chk=0.0894
Currently selected inputs: y(k-1) u(k-3) u(k-4)

Elapsed time = 10.917078

Для входного выбора другой более в вычислительном отношении интенсивный подход должен сделать исчерпывающий поиск на всех возможных комбинациях входных кандидатов. Функция, которая выполняет исчерпывающий поиск, является exhsrch, который выбирает 3 входных параметров от 10 кандидатов. Однако exhsrch обычно включает существенное количество вычисления, если все комбинации пробуют. Например, если 3 выбран из 10, общее количество моделей ANFIS является C (10, 3) = 120.

К счастью, для идентификации динамической системы, мы знаем, что входные параметры не должны прибывать ни из одной из следующих двух групп входных кандидатов исключительно:

Y = {y (k-1), y (k-2), y (k-3), y (k-4)}

U = {u (k-1), u (k-2), u (k-3), u (k-4), u (k-5), u (k-6)}

Разумное предположение должно было бы взять два входных параметров из Y и один от U, чтобы сформировать входные параметры к ANFIS; общее количество моделей ANFIS затем C (4,2) *6=36, который является намного меньше. Вышеупомянутый график показывает, что выбранные входные параметры являются y (k-1), y (k-2) и u (k-3), с учебным RMSE 0,0474 и проверяющий RMSE 0,0485, которые лучше, чем модели ARX и ANFIS через последовательный прямой поиск.

group1 = [1 2 3 4];	% y(k-1), y(k-2), y(k-3), y(k-4)
group2 = [1 2 3 4];	% y(k-1), y(k-2), y(k-3), y(k-4)
group3 = [5 6 7 8 9 10];	% u(k-1) through y(k-6)

anfis_n = 6*length(group3);
index = zeros(anfis_n,3);
trn_error = zeros(anfis_n,1);
chk_error = zeros(anfis_n,1);
% ======= Training options
% Create option set for generating initial FIS.
genOpt = genfisOptions('GridPartition','NumMembershipFunctions',2, ...
                       'InputMembershipFunctionType','gbellmf');
% Create option set for |anfis| command and set options that remain constant
% for different training scenarios.
anfisOpt = anfisOptions('EpochNumber',1,...
                        'InitialStepSize',0.1,...
                        'StepSizeDecreaseRate',0.5,...
                        'StepSizeIncreaseRate',1.5,...
                        'DisplayANFISInformation',0,...
                        'DisplayErrorValues',0,...
                        'DisplayStepSize',0,...
                        'DisplayFinalResults',0);
% ====== Train ANFIS with different input variables
fprintf('\nTrain %d ANFIS models, each with 3 inputs selected from 10 candidates...\n\n',...
    anfis_n);
model = 1;
for i = 1:length(group1)
    for j = i+1:length(group2)
        for k = 1:length(group3)
            in1 = deblank(input_name(group1(i),:));
            in2 = deblank(input_name(group2(j),:));
            in3 = deblank(input_name(group3(k),:));
            index(model, :) = [group1(i) group2(j) group3(k)];
            trn_data = data(1:trn_data_n, [group1(i) group2(j) group3(k) size(data,2)]);
            chk_data = data(trn_data_n+1:trn_data_n+300, [group1(i) group2(j) group3(k) size(data,2)]);
            in_fismat = genfis(trn_data(:,1:end-1),trn_data(:,end),genOpt);
            % Set initial FIS and validation data in option set for ANFIS training.
            anfisOpt.InitialFIS = in_fismat;
            anfisOpt.ValidationData = chk_data;
            [~, t_err, ~, ~, c_err] = anfis(trn_data,anfisOpt);
            trn_error(model) = min(t_err);
            chk_error(model) = min(c_err);
            fprintf('ANFIS model = %d: %s %s %s',model,in1,in2,in3);
            fprintf(' --> trn=%.4f,',trn_error(model));
            fprintf(' chk=%.4f',chk_error(model));
            fprintf('\n');
            model = model+1;
        end
    end
end

% ====== Reordering according to training error
[~, b] = sort(trn_error);
b = flipud(b);		% List according to decreasing trn error
trn_error = trn_error(b);
chk_error = chk_error(b);
index = index(b,:);

% ====== Display training and checking errors
x = (1:anfis_n)';
subplot(2,1,1)
plot(x, trn_error,'-',x,chk_error,'-', ...
     x,trn_error,'o',x,chk_error,'*')
tmp = x(:, ones(1,3))';
X = tmp(:);
tmp = [zeros(anfis_n,1) max(trn_error,chk_error) nan*ones(anfis_n,1)]';
Y = tmp(:);
hold on
plot(X,Y,'g')
hold off
axis([1 anfis_n -inf inf])
h_gca = gca;
h_gca.XTickLabel = [];

% ====== Add text of input variables
for k = 1:anfis_n
    text(x(k), 0, ...
        [input_name(index(k,1),:) ' ' ...
         input_name(index(k,2),:) ' ' ...
         input_name(index(k,3),:)]);
end
h = findobj(gcf,'type','text');
set(h,'rot',90,'fontsize',11,'hori','right');

drawnow

% ====== Generate input_index for bjtrain.m
[a, b] = min(trn_error);
input_index = index(b,:);
title('Training (Circles) and Checking (Asterisks) Errors','fontsize',10)
ylabel('RMSE','fontsize',10)
Train 36 ANFIS models, each with 3 inputs selected from 10 candidates...

ANFIS model = 1: y(k-1) y(k-2) u(k-1) --> trn=0.0990, chk=0.0962
ANFIS model = 2: y(k-1) y(k-2) u(k-2) --> trn=0.0852, chk=0.0862
ANFIS model = 3: y(k-1) y(k-2) u(k-3) --> trn=0.0474, chk=0.0485
ANFIS model = 4: y(k-1) y(k-2) u(k-4) --> trn=0.0808, chk=0.0822
ANFIS model = 5: y(k-1) y(k-2) u(k-5) --> trn=0.1023, chk=0.0991
ANFIS model = 6: y(k-1) y(k-2) u(k-6) --> trn=0.1021, chk=0.0974
ANFIS model = 7: y(k-1) y(k-3) u(k-1) --> trn=0.1231, chk=0.1206
ANFIS model = 8: y(k-1) y(k-3) u(k-2) --> trn=0.1047, chk=0.1085
ANFIS model = 9: y(k-1) y(k-3) u(k-3) --> trn=0.0587, chk=0.0626
ANFIS model = 10: y(k-1) y(k-3) u(k-4) --> trn=0.0806, chk=0.0836
ANFIS model = 11: y(k-1) y(k-3) u(k-5) --> trn=0.1261, chk=0.1311
ANFIS model = 12: y(k-1) y(k-3) u(k-6) --> trn=0.1210, chk=0.1151
ANFIS model = 13: y(k-1) y(k-4) u(k-1) --> trn=0.1420, chk=0.1353
ANFIS model = 14: y(k-1) y(k-4) u(k-2) --> trn=0.1224, chk=0.1229
ANFIS model = 15: y(k-1) y(k-4) u(k-3) --> trn=0.0700, chk=0.0765
ANFIS model = 16: y(k-1) y(k-4) u(k-4) --> trn=0.0817, chk=0.0855
ANFIS model = 17: y(k-1) y(k-4) u(k-5) --> trn=0.1337, chk=0.1405
ANFIS model = 18: y(k-1) y(k-4) u(k-6) --> trn=0.1421, chk=0.1333
ANFIS model = 19: y(k-2) y(k-3) u(k-1) --> trn=0.2393, chk=0.2264
ANFIS model = 20: y(k-2) y(k-3) u(k-2) --> trn=0.2104, chk=0.2077
ANFIS model = 21: y(k-2) y(k-3) u(k-3) --> trn=0.1452, chk=0.1497
ANFIS model = 22: y(k-2) y(k-3) u(k-4) --> trn=0.0958, chk=0.1047
ANFIS model = 23: y(k-2) y(k-3) u(k-5) --> trn=0.2048, chk=0.2135
ANFIS model = 24: y(k-2) y(k-3) u(k-6) --> trn=0.2388, chk=0.2326
ANFIS model = 25: y(k-2) y(k-4) u(k-1) --> trn=0.2756, chk=0.2574
ANFIS model = 26: y(k-2) y(k-4) u(k-2) --> trn=0.2455, chk=0.2400
ANFIS model = 27: y(k-2) y(k-4) u(k-3) --> trn=0.1726, chk=0.1797
ANFIS model = 28: y(k-2) y(k-4) u(k-4) --> trn=0.1074, chk=0.1157
ANFIS model = 29: y(k-2) y(k-4) u(k-5) --> trn=0.2061, chk=0.2133
ANFIS model = 30: y(k-2) y(k-4) u(k-6) --> trn=0.2737, chk=0.2836
ANFIS model = 31: y(k-3) y(k-4) u(k-1) --> trn=0.3842, chk=0.3605
ANFIS model = 32: y(k-3) y(k-4) u(k-2) --> trn=0.3561, chk=0.3358
ANFIS model = 33: y(k-3) y(k-4) u(k-3) --> trn=0.2719, chk=0.2714
ANFIS model = 34: y(k-3) y(k-4) u(k-4) --> trn=0.1763, chk=0.1808
ANFIS model = 35: y(k-3) y(k-4) u(k-5) --> trn=0.2132, chk=0.2240
ANFIS model = 36: y(k-3) y(k-4) u(k-6) --> trn=0.3460, chk=0.3601

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

if ishghandle(winH1),delete(winH1);
end

trn_data = data(1:trn_data_n,[input_index, size(data,2)]);
chk_data = data(trn_data_n+1:600,[input_index, size(data,2)]);

% generate FIS matrix
in_fismat = genfis(trn_data(:,1:end-1),trn_data(:,end));
anfisOpt = anfisOptions('InitialFIS',in_fismat,...
                        'EpochNumber',1,...
                        'InitialStepSize',0.01,...
                        'StepSizeDecreaseRate',0.5,...
                        'StepSizeIncreaseRate',1.5,...
                        'ValidationData',chk_data);
[trn_out_fismat,trn_error,step_size,chk_out_fismat,chk_error] = ...
    anfis(trn_data,anfisOpt);

subplot(2,1,1)
index = 1:trn_data_n;
plot(index,y(index),index,yp(index),'.')
rmse = norm(y(index)-yp(index))/sqrt(length(index));
title(sprintf(['(a) Training Data (Solid Line) and ARX Prediction (Dots)\nwith RMSE = ' num2str(rmse)]),'fontsize',10)
disp(['[na nb d] = ' num2str(nn)])
xlabel('Time Steps','fontsize',10)
subplot(2,1,2)
index = (trn_data_n+1):(total_data_n);
plot(index, y(index),index,yp(index),'.')
rmse = norm(y(index)-yp(index))/sqrt(length(index));
title(sprintf(['(b) Checking Data (Solid Line) and ARX Prediction (Dots)\nwith RMSE = ' num2str(rmse)]),'fontsize',10)
xlabel('Time Steps','fontsize',10)
ANFIS info: 
	Number of nodes: 34
	Number of linear parameters: 32
	Number of nonlinear parameters: 18
	Total number of parameters: 50
	Number of training data pairs: 300
	Number of checking data pairs: 300
	Number of fuzzy rules: 8


Start training ANFIS ...

   1 	 0.0474113 	 0.0485325

Designated epoch number reached --> ANFIS training completed at epoch 1.

Minimal training RMSE = 0.047411
Minimal checking RMSE = 0.0485325
[na nb d] = 5  10   2

y_hat = evalfis(chk_out_fismat,data(1:600,input_index));

subplot(2,1,1)
index = 1:trn_data_n;
plot(index,data(index,size(data,2)),'-', ...
     index,y_hat(index),'.')
rmse = norm(y_hat(index)-data(index,size(data,2)))/sqrt(length(index));
title(sprintf(['Training Data (Solid Line) and ANFIS Prediction (Dots)\nwith RMSE = ' num2str(rmse)]),'fontsize',10)
xlabel('Time Index','fontsize',10)
ylabel('')

subplot(2,1,2)
index = trn_data_n+1:600;
plot(index,data(index,size(data,2)),'-',index,y_hat(index),'.')
rmse = norm(y_hat(index)-data(index,size(data,2)))/sqrt(length(index));
title(sprintf(['Checking Data (Solid Line) and ANFIS Prediction (Dots)\nwith RMSE = ' num2str(rmse)]),'fontsize',10)
xlabel('Time Index','fontsize',10)
ylabel('')

Заключение

Приведенная выше таблица является сравнением среди различных подходов моделирования. Моделирование ARX проводит наименьшее количество количества времени, чтобы достигнуть худшей точности, и ANFIS, моделирующий через исчерпывающий поиск, занимает большую часть количества времени, чтобы достигнуть лучшей точности. Другими словами, если быстрое моделирование является целью, то ARX является правильным выбором. Но если точность является предельным беспокойством, то мы должны пойти с ANFIS, который разработан для нелинейного моделирования и более высокой точности.

Смотрите также

| |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте