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

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

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

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

Набор данных, используемый в этом примере для ANFIS и моделирования ARX, является от преподавателя 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 objects. Axes object 1 with title Output Data contains 2 objects of type line. Axes object 2 with title Input Data contains 2 objects of type line.

Точки данных отражают шаг расчета 0,08 секунд. Вход $u (k)$ бинарная случайная перемена сигнала между 3,41 и 6.41. Вероятность сдвига входа на каждой выборке 0.2. Набор данных содержит 1 000 точек данных ввода/вывода. Эти графики показывают выходную температуру $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 objects. Axes object 1 with title Training Data (solid), ARX Prediction (dots) RMSE = 0.11208 contains 2 objects of type line. Axes object 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 object. The axes object 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 object. The axes object 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 objects. Axes object 1 with title Training Data (solid), ANFIS Prediction (dots) RMSE = 0.047411 contains 2 objects of type line. Axes object 2 with title Validation Data (solid), ANFIS Prediction (dots) RMSE = 0.048532 contains 2 objects of type line.

Предсказания модели ANFIS соответствуют данным намного более тесно, чем предсказания модели ARX.

Ссылка

[1] Ljung, Lennart. System Identification: теория для пользователя. Информация о Prentice Hall и системный научный ряд. Englewood Cliffs, NJ: Prentice Hall, 1987.

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

anfis | genfis | evalfis

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

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

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

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

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

Ссылка

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

Похожие темы

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

Поддержка