Нелинейная идентификация системы

В этом примере используются:

В этом примере показано, как выполнить динамическую идентификацию системы с использованием линейной модели ARX и нелинейной модели ANFIS.

Загрузить данные

Набор данных, используемый в этом примере для моделирования ANFIS и ARX, взят из лабораторного устройства «Feedback's Process Trainer PT 326» [1]. Устройство функционирует как фен: воздух раздувается через трубку и нагревается на входе. Термопара измеряет температуру воздуха. Вход $u (k$ ) представляет собой напряжение на сетке резисторных проводов для нагрева поступающего воздуха, а выход $y ($ k) представляет собой выходную температуру воздуха.

Загрузите тестовые данные и постройте график ввода и вывода.

load dryerdata
data_n = length(y);
output = y;
input = [[0; y(1:data_n-1)] ...
        [0; 0; y(1:data_n-2)] ...
        [0; 0; 0; y(1:data_n-3)] ...
        [0; 0; 0; 0; y(1:data_n-4)] ...
        [0; u(1:data_n-1)] ...
        [0; 0; u(1:data_n-2)] ...
        [0; 0; 0; u(1:data_n-3)] ...
        [0; 0; 0; 0; u(1:data_n-4)] ...
        [0; 0; 0; 0; 0; u(1:data_n-5)] ...
        [0; 0; 0; 0; 0; 0; u(1:data_n-6)]];
data = [input output];
data(1:6,:) = [];
input_name = ["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;

figure
subplot(2,1,1)
plot(index,y(index),'-',index,y(index),'o')
title('Output Data')
ylabel('y(k)')
subplot(2,1,2)
plot(index,u(index),'-',index,u(index),'o')
title('Input Data')
ylabel('u(k)')

Figure contains 2 axes. Axes 1 with title Output Data contains 2 objects of type line. Axes 2 with title Input Data contains 2 objects of type line.

Точки данных отражают время выборки 0,08 секунды. Вход $u (k$ ) является двоичным случайным сигналом, сдвигающимся между 3,41 и 6,41. Вероятность сдвига входного сигнала в каждой выборке составляет 0,2. Набор данных содержит 1000 точек данных ввода/вывода. Эти графики показывают выходную температуру $y ($ k) и входное $напряжение u$ (k) для первых 100 временных шагов.

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

Модель ARX - это линейная модель следующего вида:

$y (k)_{} +a1⋅y (k-1) + ._{.} . +am\cdoty (_{} k-m) =b1\cdotu (_{k-d}) + . . .$ +bn⋅u (k-d-n + 1)

Здесь:

$y (k$ ) и $u ($ k) являются усредненными версиями исходных данных.
$_{ai}$ и $_{bj}$ - линейные параметры.
$m$ , $n$ и $d$ - три целых числа, точно определяющие модель ARX.

Чтобы найти модель ARX для сушильного устройства, сначала разделите набор данных на обучающий ( $k$ = 1-300) и проверочный ( $k$ = 301-600) набор.

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

Выполните исчерпывающий поиск, чтобы найти лучшую комбинацию $m$ , $n$ и $d$ , позволяя каждому целому числу изменяться от 1 до 10 независимо. Для выполнения поиска и выбора параметров ARX используйте arxstruc(Панель инструментов идентификации системы) и selstruc(Панель инструментов идентификации системы).

% Run through all different models.
V = arxstruc(ze,zv,struc(1:10,1:10,1:10));
% Find the best model.
nn = selstruc(V,0);
% Display model parameters
disp(['[m n d] = ' num2str(nn)])

[m n d] = 5  10   2

Лучшая модель ARX имеет $m$ = 5, $n$ = 10 и $d$ = 2. Создайте с обучающей корневой среднеквадратичной ошибкой (RMSE) 0,1122 и проверкой RMSE 0,0749. Постройте график исходной модели y $(k)$ вместе с этой моделью ARX.

Создайте и смоделируйте модель ARX с этими параметрами.

th = arx(ze,nn);
th.Ts = 0.08;
u = z(:,2);
y = z(:,1) + ave;
yp = sim(u,th) + ave;

Постройте график вывода модели ARX по данным обучения и проверки. Среднеквадратическая ошибка обучающего корня (RMSE) равна 0,1121, а RMSE проверки - 0,0748.

figure
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("Training Data (solid), ARX Prediction (dots)" + newline + "RMSE = " + num2str(rmse))
xlabel('Time Steps')

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("Validation Data (solid), ARX Prediction (dots)" + newline + "RMSE = " + num2str(rmse))
xlabel('Time Steps')

Figure contains 2 axes. Axes 1 with title Training Data (solid), ARX Prediction (dots) RMSE = 0.11208 contains 2 objects of type line. Axes 2 with title Validation Data (solid), ARX Prediction (dots) RMSE = 0.07483 contains 2 objects of type line.

Определение модели ANFIS

Модель ARX является линейной и может быстро выполнять определение структуры модели и параметров. Результаты работы на предыдущих участках представляются удовлетворительными. Однако, если вы хотите улучшить производительность, вы можете попробовать нелинейную модель, такую как адаптивная нейро-нечеткая система вывода (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).

Выполните последовательный прямой поиск входов с помощью функции sequentialSearch. Эта функция последовательно выбирает каждую входную переменную для оптимизации RMSE.

trn_data_n = 300;
trn_data = data(1:trn_data_n,:);
val_data = data(trn_data_n+1:trn_data_n+300,:);
[~,elapsed_time] = sequentialSearch(3,trn_data,val_data,input_name);

Selecting input 1 ...
ANFIS model 1: y(k-1) --> trn=0.2043, val=0.1888
ANFIS model 2: y(k-2) --> trn=0.3819, val=0.3541
ANFIS model 3: y(k-3) --> trn=0.5245, val=0.4903
ANFIS model 4: y(k-4) --> trn=0.6308, val=0.5977
ANFIS model 5: u(k-1) --> trn=0.8271, val=0.8434
ANFIS model 6: u(k-2) --> trn=0.7976, val=0.8087
ANFIS model 7: u(k-3) --> trn=0.7266, val=0.7349
ANFIS model 8: u(k-4) --> trn=0.6215, val=0.6346
ANFIS model 9: u(k-5) --> trn=0.5419, val=0.5650
ANFIS model 10: u(k-6) --> trn=0.5304, val=0.5601
Currently selected inputs: y(k-1)

Selecting input 2 ...
ANFIS model 11: y(k-1) y(k-2) --> trn=0.1085, val=0.1024
ANFIS model 12: y(k-1) y(k-3) --> trn=0.1339, val=0.1283
ANFIS model 13: y(k-1) y(k-4) --> trn=0.1542, val=0.1461
ANFIS model 14: y(k-1) u(k-1) --> trn=0.1892, val=0.1734
ANFIS model 15: y(k-1) u(k-2) --> trn=0.1663, val=0.1574
ANFIS model 16: y(k-1) u(k-3) --> trn=0.1082, val=0.1077
ANFIS model 17: y(k-1) u(k-4) --> trn=0.0925, val=0.0948
ANFIS model 18: y(k-1) u(k-5) --> trn=0.1533, val=0.1531
ANFIS model 19: y(k-1) u(k-6) --> trn=0.1952, val=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, val=0.0822
ANFIS model 21: y(k-1) u(k-4) y(k-3) --> trn=0.0806, val=0.0836
ANFIS model 22: y(k-1) u(k-4) y(k-4) --> trn=0.0817, val=0.0855
ANFIS model 23: y(k-1) u(k-4) u(k-1) --> trn=0.0886, val=0.0912
ANFIS model 24: y(k-1) u(k-4) u(k-2) --> trn=0.0835, val=0.0843
ANFIS model 25: y(k-1) u(k-4) u(k-3) --> trn=0.0609, val=0.0604
ANFIS model 26: y(k-1) u(k-4) u(k-5) --> trn=0.0848, val=0.0867
ANFIS model 27: y(k-1) u(k-4) u(k-6) --> trn=0.0890, val=0.0894
Currently selected inputs: y(k-1) u(k-3) u(k-4)

Figure contains an axes. The axes with title Error for Corresponding Inputs contains 3 objects of type line. These objects represent Training, Validation.

На этом графике показаны все комбинации входных данных, проверенных sequentialSearch. Поиск выбирает $y (k-1$ ), $u ($ k-3) $и u ($ k-4) в качестве входных данных, так как модель с этими входными данными имеет самый низкий RMSE обучения (0,609) и RMSE подтверждения (0,604).

В качестве альтернативы используйте исчерпывающий поиск по всем возможным комбинациям входных кандидатов. Как и прежде, ищите три входа из 10 кандидатов. Можно использовать функцию exhaustiveSearch для такого поиска; однако эта функция пробует все возможные комбинации кандидатов, $(\binom{}{103}) =$ 120 в данном случае.

Вместо exhaustiveSearchиспользуйте пользовательский код для поиска по подмножеству этих комбинаций. В этом примере не выбирайте какую-либо комбинацию входных данных исключительно из входных данных или исключительно из выходных данных.

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

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 u(k-6)

Укажите параметры и параметры обучения.

anfis_n = 6*length(group3);
index = zeros(anfis_n,3);
trn_error = zeros(anfis_n,1);
val_error = zeros(anfis_n,1);

% Create option set for generating initial FIS.
genOpt = genfisOptions('GridPartition','NumMembershipFunctions',2, ...
                       'InputMembershipFunctionType','gbellmf');
% Create option set for anfis function 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);

Модель ANFIS для каждой входной комбинации.

model = 1;
for i = 1:length(group1)
    for j = i+1:length(group2)
        for k = 1:length(group3)
            % Create input combinations.
            in1 = input_name(group1(i));
            in2 = input_name(group2(j));
            in3 = 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)]);
            val_data = data(trn_data_n+1:trn_data_n+300, ...
                [group1(i) group2(j) group3(k) size(data,2)]);
            
            % Create the initial FIS structure.
            in_fis = genfis(trn_data(:,1:end-1),trn_data(:,end),genOpt);
            
            % Set the initial FIS and validation data for ANFIS training.
            anfisOpt.InitialFIS = in_fis;
            anfisOpt.ValidationData = val_data;
            
            % Train the ANFIS system.
            [~,t_err,~,~,c_err] = anfis(trn_data,anfisOpt);
            trn_error(model) = min(t_err);
            val_error(model) = min(c_err);
            fprintf('ANFIS model = %d: %s %s %s',model,in1,in2,in3);
            fprintf(' --> trn=%.4f,',trn_error(model));
            fprintf(' val=%.4f',val_error(model));
            fprintf('\n');
            model = model+1;
        end
    end
end

ANFIS model = 1: y(k-1) y(k-2) u(k-1)

 --> trn=0.0990,

 val=0.0962

ANFIS model = 2: y(k-1) y(k-2) u(k-2)

 --> trn=0.0852,

 val=0.0862

ANFIS model = 3: y(k-1) y(k-2) u(k-3)

 --> trn=0.0474,

 val=0.0485

ANFIS model = 4: y(k-1) y(k-2) u(k-4)

 --> trn=0.0808,

 val=0.0822

ANFIS model = 5: y(k-1) y(k-2) u(k-5)

 --> trn=0.1023,

 val=0.0991

ANFIS model = 6: y(k-1) y(k-2) u(k-6)

 --> trn=0.1021,

 val=0.0974

ANFIS model = 7: y(k-1) y(k-3) u(k-1)

 --> trn=0.1231,

 val=0.1206

ANFIS model = 8: y(k-1) y(k-3) u(k-2)

 --> trn=0.1047,

 val=0.1085

ANFIS model = 9: y(k-1) y(k-3) u(k-3)

 --> trn=0.0587,

 val=0.0626

ANFIS model = 10: y(k-1) y(k-3) u(k-4)

 --> trn=0.0806,

 val=0.0836

ANFIS model = 11: y(k-1) y(k-3) u(k-5)

 --> trn=0.1261,

 val=0.1311

ANFIS model = 12: y(k-1) y(k-3) u(k-6)

 --> trn=0.1210,

 val=0.1151

ANFIS model = 13: y(k-1) y(k-4) u(k-1)

 --> trn=0.1420,

 val=0.1353

ANFIS model = 14: y(k-1) y(k-4) u(k-2)

 --> trn=0.1224,

 val=0.1229

ANFIS model = 15: y(k-1) y(k-4) u(k-3)

 --> trn=0.0700,

 val=0.0765

ANFIS model = 16: y(k-1) y(k-4) u(k-4)

 --> trn=0.0817,

 val=0.0855

ANFIS model = 17: y(k-1) y(k-4) u(k-5)

 --> trn=0.1337,

 val=0.1405

ANFIS model = 18: y(k-1) y(k-4) u(k-6)

 --> trn=0.1421,

 val=0.1333

ANFIS model = 19: y(k-2) y(k-3) u(k-1)

 --> trn=0.2393,

 val=0.2264

ANFIS model = 20: y(k-2) y(k-3) u(k-2)

 --> trn=0.2104,

 val=0.2077

ANFIS model = 21: y(k-2) y(k-3) u(k-3)

 --> trn=0.1452,

 val=0.1497

ANFIS model = 22: y(k-2) y(k-3) u(k-4)

 --> trn=0.0958,

 val=0.1047

ANFIS model = 23: y(k-2) y(k-3) u(k-5)

 --> trn=0.2048,

 val=0.2135

ANFIS model = 24: y(k-2) y(k-3) u(k-6)

 --> trn=0.2388,

 val=0.2326

ANFIS model = 25: y(k-2) y(k-4) u(k-1)

 --> trn=0.2756,

 val=0.2574

ANFIS model = 26: y(k-2) y(k-4) u(k-2)

 --> trn=0.2455,

 val=0.2400

ANFIS model = 27: y(k-2) y(k-4) u(k-3)

 --> trn=0.1726,

 val=0.1797

ANFIS model = 28: y(k-2) y(k-4) u(k-4)

 --> trn=0.1074,

 val=0.1157

ANFIS model = 29: y(k-2) y(k-4) u(k-5)

 --> trn=0.2061,

 val=0.2133

ANFIS model = 30: y(k-2) y(k-4) u(k-6)

 --> trn=0.2737,

 val=0.2836

ANFIS model = 31: y(k-3) y(k-4) u(k-1)

 --> trn=0.3842,

 val=0.3605

ANFIS model = 32: y(k-3) y(k-4) u(k-2)

 --> trn=0.3561,

 val=0.3358

ANFIS model = 33: y(k-3) y(k-4) u(k-3)

 --> trn=0.2719,

 val=0.2714

ANFIS model = 34: y(k-3) y(k-4) u(k-4)

 --> trn=0.1763,

 val=0.1808

ANFIS model = 35: y(k-3) y(k-4) u(k-5)

 --> trn=0.2132,

 val=0.2240

ANFIS model = 36: y(k-3) y(k-4) u(k-6)

 --> trn=0.3460,

 val=0.3601

Постройте график ошибок обучения и проверки для каждой комбинации ввода в порядке убывания.

% Reorder according to training error.
[~, b] = sort(trn_error);
b = flipud(b);
trn_error = trn_error(b);
val_error = val_error(b);
index = index(b,:);

% Plot training and validation error.
x = (1:anfis_n)';
tmp = x(:, ones(1,3))';
X = tmp(:);
tmp = [zeros(anfis_n,1) max(trn_error,val_error) nan*ones(anfis_n,1)]';
Y = tmp(:);
figure
plot(x,trn_error,'-o',x,val_error,'-*',X,Y,'g')
title('Error for Corresponding Inputs')
ylabel('RMSE')
legend('Training','Validation','Location','northeast')

% Add ticks and labels.
labels = string(zeros(anfis_n,1));
for k = 1:anfis_n
    labels(k) = input_name(index(k,1))+ " & " + ...
	            input_name(index(k,2))+ " & " + ...
	            input_name(index(k,3));
end
xticks(x)
xticklabels(labels)
xtickangle(90)

Figure contains an axes. The axes with title Error for Corresponding Inputs contains 3 objects of type line. These objects represent Training, Validation.

Алгоритм выбирает входы $y (k-1$ ), $y ($ k-2) $и u ($ k-3) с обучающим RMSE 0,0474 и проверочным RMSE 0,0485. Эти значения RMSE улучшаются по сравнению с моделями ARX и моделью ANFIS, найденной путем последовательного прямого поиска.

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

Для этого сначала создайте набор данных.

[~,b] = min(trn_error);
input_index = index(b,:);
trn_data = data(1:trn_data_n,[input_index, size(data,2)]);
val_data = data(trn_data_n+1:600,[input_index, size(data,2)]);

Создание и обучение ANFIS.

in_fis = genfis(trn_data(:,1:end-1),trn_data(:,end));
anfisOpt = anfisOptions('InitialFIS',in_fis,...
                        'EpochNumber',1,...
                        'InitialStepSize',0.01,...
                        'StepSizeDecreaseRate',0.5,...
                        'StepSizeIncreaseRate',1.5,...
                        'ValidationData',val_data);
[trn_out_fis,trn_error,step_size,val_out_fis,val_error] = ...
    anfis(trn_data,anfisOpt);

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.0474113
Minimal checking RMSE = 0.0485325

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

figure
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("Training Data (solid), ANFIS Prediction (dots)" + newline + "RMSE = " + num2str(rmse))
xlabel('Time Steps')

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("Validation Data (solid), ANFIS Prediction (dots)" + newline + "RMSE = " + num2str(rmse))
xlabel('Time Steps')

Figure contains 2 axes. Axes 1 with title Training Data (solid), ANFIS Prediction (dots) RMSE = 0.047411 contains 2 objects of type line. Axes 2 with title Validation Data (solid), ANFIS Prediction (dots) RMSE = 0.048532 contains 2 objects of type line.

Прогнозы модели ANFIS подходят к данным гораздо ближе, чем прогнозы модели ARX.

Ссылка

[1] Люнг, Леннарт. Идентификация системы: теория для пользователя. Серия информационных и системных наук Прентис-Холла. Энглвуд Клиффс, Нью-Джерси: Прентис-Холл, 1987.

См. также

anfis | evalfis | genfis

Связанные темы

Нейро-адаптивное обучение и ANFIS

Документация