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

Этот пример использует:

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

Загрузка данных

Набор данных, используемый в этом примере для моделирования ANFIS и ARX, получен из лабораторного устройства «Feedback's Process Trainer PT 326» [1]. Устройство функционирует как фен: воздух раздувается через трубку и нагревается у входного отверстия. Термопара измеряет температуру воздуха. Вход $u (k)$ - напряжение на mesh резисторных проводов для нагрева поступающего воздуха и выхода $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) + a_{1} \cdot y (k - 1) + . . . + a_{m} \cdot y (k - m) = b_{1} \cdot u (k - d) + . . . + b_{n} \cdot u (k - d - n + 1)$

Здесь:

$y (k)$ и $u (k)$ являются средними вычитаемыми версиями исходных данных.
$a_{i}$ и $b_{j}$ являются линейными параметрами.
$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 (System Identification Toolbox) и selstruc (System Identification Toolbox) функции.

% 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{10}{3}) = 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] Ljung, Lennart. Система идентификации: теория для пользователя. Серия информации и системных наук Prentice Hall. Englewood Cliffs, Нью-Джерси: Prentice Hall, 1987.

См. также

anfis | evalfis | genfis

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

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

Загрузка данных

Идентифицируйте модель ARX

Идентифицируйте модель ANFIS

Ссылка

См. также

Похожие темы

Документация по Fuzzy Logic Toolbox

Поддержка