Передаточная функция

Этот пример использует передаточную функцию четвертого порядка со смешанной быстрой и медленной динамикой и умеренным демпфированием. Умеренное демпфирование заставляет динамику системы демпфироваться в течение более длительного временного горизонта и показывает способность сети LSTM захватывать смешанную динамику без некоторой важной динамики отклика демпфирования. Создайте передаточную функцию путем определения полюсов и нулей системы.

fourthOrderMdl = zpk(-4,[-9+5i;-9-5i;-2+50i;-2-50i],5e5);
[stepResponse,stepTime] = step(fourthOrderMdl);

Постройте график переходной характеристики передаточной функции.

plot(stepTime,stepResponse)
grid on
axis tight
title('Fourth-order mixed step response')

Диаграмма Боде показывает полосу пропускания системы, которая измеряется как первая частота, где коэффициент усиления падает ниже 70,8% или около 3 дБ от своего значения постоянного тока.

bodeplot(fourthOrderMdl)

fb = bandwidth(fourthOrderMdl)

fb = 62.8858

Сгенерируйте обучающие данные

Создайте набор данных из входа и выходных сигналов, которые можно использовать для обучения сети LSTM. Для входов сгенерируйте случайный Гауссов сигнал. Симулируйте ответ передаточной функции fourthOrderMdl на этот вход для получения сигнала выхода.

Данные Гауссова шума

Задайте свойства для случайного сигнала Гауссова обучения шуму.

signalType = 'rgs'; % Gaussian 
signalLength = 5000; % Number of points in the signal
fs = 100; % Sampling frequency
signalAmplitude = 1; % Maximum signal amplitude

Сгенерируйте сигнал Гауссова шума, используя idinput и масштабируйте результат.

urgs = idinput(signalLength,signalType);
urgs = (signalAmplitude/max(urgs))*urgs';

Сгенерируйте сигнал времени на основе частоты дискретизации.

trgs = 0:1/fs:length(urgs)/fs-1/fs;

Используйте lsim функция для генерации отклика системы и хранения результата в yrgs. Транспонируйте моделируемый выход так, чтобы она соответствовала структуре данных LSTM, которая требует векторов-строк, а не векторов-столбцов.

yrgs = lsim(fourthOrderMdl,urgs,trgs);
yrgs = yrgs';

Точно так же создайте более короткий сигнал валидации для использования во время сетевого обучения.

xval = idinput(100,signalType);
yval = lsim(fourthOrderMdl,xval,trgs(1:100));

Создание и обучение сети

Следующая сетевая архитектура была определена при помощи байесовской стандартной программы оптимизации, где функция оптимизации Байеса использует независимые данные валидации (см. bayesianOptimizationForLSTM.mlx для получения дополнительной информации). Хотя несколько архитектур могут работать, эта оптимизация обеспечивает наиболее вычислительно эффективную. Процесс оптимизации также показал, что как только сложность передаточной функции увеличивается при применении LSTM к другим линейным передаточным функциям, архитектура сети существенно не меняется. Скорее число эпох, необходимых для обучения сети, увеличивается. Количество скрытых модулей, необходимых для моделирования системы, связано с тем, сколько времени требуется динамике для увлажнения. В этом случае существует две отдельные части отклика: Высокая частота реакция и низкая частотная характеристика. Для захвата низкочастотной характеристики требуется большее количество скрытых модулей. Если выбрано меньшее количество модулей, высокая частота реакция все еще моделируется. Однако оценка низкочастотной характеристики ухудшается.

Создайте сетевую архитектуру.

numResponses = 1;
featureDimension = 1;
numHiddenUnits = 100;
maxEpochs = 1000;
miniBatchSize = 200;

Networklayers = [sequenceInputLayer(featureDimension) ...
    lstmLayer(numHiddenUnits) ...
    lstmLayer(numHiddenUnits) ...
    fullyConnectedLayer(numResponses) ...
    regressionLayer];

Начальная скорость обучения влияет на успех сети. Использование начальной скорости обучения, которая является слишком высокой, приводит к высоким градиентам, которые приводят к более длительному времени обучения. Более длительное время обучения может привести к насыщению полностью подключенного уровня сети. Когда сеть насыщается, выходы расходятся, и сеть выводит NaN значение. Следовательно, используйте значение по умолчанию 0,01, которое является относительно низкой начальной скоростью обучения. Это приводит к монотонному уменьшению кривых невязки и потерь. Используйте кусочное расписание скорости, чтобы алгоритм оптимизации не попал в ловушку локальных минимумов в начале стандартной программы оптимизации.

options = trainingOptions('adam', ...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'GradientThreshold',10, ...
    'Shuffle','once', ...
    'Plots','training-progress',...
    'ExecutionEnvironment','gpu',...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropPeriod',100,...
    'Verbose',0,...
    'ValidationData',[{xval'} {yval'}]);

loadNetwork = true; % Set to true to train the network using a parallel pool.
if loadNetwork
    load('fourthOrderMdlnet','fourthOrderNet')
else
    poolobj = parpool;
    fourthOrderNet = trainNetwork(urgs,yrgs,Networklayers,options);
    delete(poolobj)
    save('fourthOrderMdlnet','fourthOrderNet','urgs','yrgs');
end

Оценка производительности модели

Сеть работает хорошо, если она успешна в захвате динамического поведения системы. Оцените эффективность сети путем измерения способности сети точно предсказать отклик системы к шаговому входу.

Создайте вход шага.

stepTime = 2; % In seconds
stepAmplitude = 0.1;
stepDuration = 4; % In seconds

% Construct step signal and system output.
time = (0:1/fs:stepDuration)';
stepSignal = [zeros(sum(time<=stepTime),1);stepAmplitude*ones(sum(time>stepTime),1)];
systemResponse = lsim(fourthOrderMdl,stepSignal,time);

% Transpose input and output signal for network inputs.
stepSignal = stepSignal';
systemResponse = systemResponse';

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

fourthOrderMixedNetStep = predict(fourthOrderNet,stepSignal);

figure
title('Step response estimation')
plot(time,systemResponse,'k', time, fourthOrderMixedNetStep)
grid on
legend('System','Estimated')
title('Fourth-Order Step')

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

Инициализация сети и настройка подгонки

Чтобы инициализировать состояние сети до правильного начального условия, необходимо обновить состояние сети, чтобы оно соответствовало состоянию системы в начале тестового сигнала.

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

Установите начальное состояние сети

Когда сеть выполняет оценку с использованием шага входа от 0 до 1, состояния сети LSTM (камер и скрытых состояний слоев LSTM) дрейфуют к правильному начальному условию. Чтобы визуализировать это, извлекайте камеру и скрытое состояние сети в каждый временной шаг с помощью predictAndUpdateState функция.

Используйте только значения камер и скрытых состояний перед шагом, который происходит через 2 секунды. Задайте маркер времени в течение 2 секунд и извлеките значения до этого маркера.

stepMarker = time <= 2;
yhat = zeros(sum(stepMarker),1);
hiddenState = zeros(sum(stepMarker),200); % 200 LSTM units
cellState = zeros(sum(stepMarker),200);
for ntime = 1:sum(stepMarker)
    [fourthOrderNet,yhat(ntime)] = predictAndUpdateState(fourthOrderNet,stepSignal(ntime)');
    hiddenState(ntime,:) = fourthOrderNet.Layers(2,1).HiddenState;
    cellState(ntime,:) = fourthOrderNet.Layers(2,1).CellState;
end

Затем постройте график скрытого состояния и состояния камер за период до шага и подтвердите, что они сходятся к фиксированным значениям.

figure
subplot(2,1,1)
plot(time(1:200),hiddenState(1:200,:))
grid on
axis tight
title('Hidden State')
subplot(2,1,2)
plot(time(1:200),cellState(1:200,:))
grid on
axis tight
title('Cell State')

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

initializationSignalDuration = 10; % In seconds
initializationValue = 0;
initializationSignal = initializationValue*ones(1,initializationSignalDuration*fs);

fourthOrderNet = predictAndUpdateState(fourthOrderNet,initializationSignal);

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

figure
zeroMapping = predict(fourthOrderNet,initializationSignal);
plot(zeroMapping)
axis tight

Теперь, когда сеть правильно инициализирована, используйте сеть, чтобы предсказать переходную характеристику снова и построить график результатов. Начальное нарушение порядка исчезло.

fourthOrderMixedNetStep = predict(fourthOrderNet,stepSignal);

figure
title('Step response estimation')
plot(time,systemResponse,'k', ...
    time,fourthOrderMixedNetStep,'b')
grid on
legend('System','Estimated')
title('Fourth-Order Step - Adjusted State')

Настройка смещения сети

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

bias = mean(predict(fourthOrderNet,initializationSignal));
fourthOrderMixedNetStep = fourthOrderMixedNetStep-bias;

figure
title('Step response estimation')
plot(time,systemResponse,'k',time,fourthOrderMixedNetStep,'b-')
legend('System','Estimated')
title('Fourth-Order Step - Adjusted Offset')

Выезд за пределы области значений обучения

Все сигналы, используемые для обучения сети, имели максимальную амплитуду 1, а шаговая функция - 0,1. Теперь исследуйте поведение сети за пределами этих областей значений.

Сдвиг во времени

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

stepTime = 3; % In seconds
stepAmplitude = 0.1;
stepDuration = 5; % In seconds
[stepSignal,systemResponse,time] = generateStepResponse(fourthOrderMdl,stepTime,stepAmplitude,stepDuration);

fourthOrderMixedNetStep = predict(fourthOrderNet,stepSignal);
bias = fourthOrderMixedNetStep(1) - initializationValue;
fourthOrderMixedNetStep = fourthOrderMixedNetStep-bias;

figure
plot(time,systemResponse,'k', time,fourthOrderMixedNetStep,'b')
grid on
axis tight

Амплитудный сдвиг

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

figure
histogram(urgs,'Normalization','pdf')
grid on

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

pValues = [60:2:98, 90:1:98, 99:0.1:99.9 99.99];
stepAmps = prctile(urgs,pValues); % Amplitudes
stepTime = 3; % In seconds
stepDuration = 5; % In seconds

stepMSE = zeros(length(stepAmps),1);
fourthOrderMixedNetStep = cell(length(stepAmps),1);
steps = cell(length(stepAmps),1);

for nAmps = 1:length(stepAmps)
    % Fourth-order mixed
    [stepSignal,systemResponse,time] = generateStepResponse(fourthOrderMdl,stepTime,stepAmps(nAmps),stepDuration);
    
    fourthOrderMixedNetStep{nAmps} = predict(fourthOrderNet,stepSignal);
    bias = fourthOrderMixedNetStep{nAmps}(1) - initializationValue;
    fourthOrderMixedNetStep{nAmps} = fourthOrderMixedNetStep{nAmps}-bias;
    
    stepMSE(nAmps) = sqrt(sum((systemResponse-fourthOrderMixedNetStep{nAmps}).^2));
    steps{nAmps,1} = systemResponse;
end

figure
plot(pValues,stepMSE,'bo')
title('Prediction Error as a Function of Deviation from Training Rrange')
grid on
axis tight

subplot(2,1,1)
plot(time,steps{1},'k', time,fourthOrderMixedNetStep{1},'b')
grid on
axis tight
title('Best Performance')
xlabel('time')
ylabel('System Response')
subplot(2,1,2)
plot(time,steps{end},'k', time,fourthOrderMixedNetStep{end},'b')
grid on
axis tight
title('Worst Performance')
xlabel('time')
ylabel('System Response')

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

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

Изменение полосы пропускания системы

Исследуйте эффект полосы пропускания системы на количество скрытых модулей, выбранных для сети LSTM, путем моделирования передаточной функции смешанной динамики четвертого порядка с четырьмя различными сетями:

Загрузите обученные сети.

load('variousHiddenUnitNets.mat')

Сгенерируйте сигнал шага.

stepTime = 2; % In seconds
stepAmplitude = 0.1;
stepDuration = 4; % In seconds

% Construct step signal.
time = (0:1/fs:stepDuration)';

stepSignal = [zeros(sum(time<=stepTime),1);stepAmplitude*ones(sum(time>stepTime),1)];
systemResponse = lsim(fourthOrderMdl,stepSignal,time);

% Transpose input and output signal for network inputs.
stepSignal = stepSignal';
systemResponse = systemResponse';

Оцените отклик системы с помощью различных обученных сетей.

smallNetStep = predict(smallNet,stepSignal)-smallNetZeroMapping(end);
medNetStep = predict(medNet,stepSignal)-medNetZeroMapping(end);
fullnetStep = predict(fullNet,stepSignal) - fullNetZeroMapping(end);
doubleNetStep = predict(doubleNet,stepSignal) - doubleNetZeroMapping(end);
 

Постройте график расчетного отклика.

figure
title('Step response estimation')
plot(time,systemResponse,'k', ...
    time,doubleNetStep,'g', ...
    time,fullnetStep,'r', ...
    time,medNetStep,'c', ...
    time,smallNetStep,'b')
grid on
legend({'System','Double Net','Full Net','Med Net','Small Net'},'Location','northwest')
title('Fourth-Order Step')

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

figure
title('Slow dynamics component')
plot(time,movmean(systemResponse,50),'k')
hold on
plot(time,movmean(doubleNetStep,50),'g')
plot(time,movmean(fullnetStep,50),'r')
plot(time,movmean(medNetStep,50),'c')
plot(time,movmean(smallNetStep,50),'b')
grid on
legend('System','Double Net','Full net','Med Net','Small Net','Location','northwest')
title('Fourth Order Step')

Добавьте шум к измеренному Отклику системы

Добавьте случайный шум к выходу системы, чтобы исследовать эффект шума на эффективность LSTM. Для этого добавьте белый шум с уровнями 1%, 5% и 10% к измеряемым откликам системы. Используйте зашумленные данные для обучения сети LSTM. С теми же наборами зашумленных данных оцените линейные модели с помощью tfest. Симулируйте эти модели и используйте симулированные отклики в качестве базового уровня для сравнения эффективности.

Используйте ту же функцию шага, что и прежде:

stepTime = 2; % In seconds
stepAmplitude = 0.1;
stepDuration = 4; % In seconds

[stepSignal,systemResponse,time] = generateStepResponse(fourthOrderMdl,stepTime,stepAmplitude,stepDuration);

Загрузите обученные сети и оцените отклик системы.

load('noisyDataNetworks.mat')
netNoise1Step = predictAndAdjust(netNoise1,stepSignal,initializationSignal,initializationValue);
netNoise5Step = predictAndAdjust(netNoise5,stepSignal,initializationSignal,initializationValue);
netNoise10Step = predictAndAdjust(netNoise10,stepSignal,initializationSignal,initializationValue);

Оценщик передаточной функции (tfest) используется для оценки функции на вышеуказанных уровнях шума для сравнения устойчивости сетей к шуму (см. сопутствующие noiseLevelModels.m для получения дополнительной информации.

load('noisyDataTFs.mat')
tfStepNoise1 = lsim(tfNoise1,stepSignal,time);
tfStepNoise5 = lsim(tfNoise5,stepSignal,time);
tfStepNoise10 = lsim(tfNoise10,stepSignal,time);

Постройте график сгенерированных откликов.

figure
plot(time,systemResponse,'k', ...
    time,netNoise1Step, ...
    time,netNoise5Step, ...
    time,netNoise10Step)
grid on
legend('System Response','1% Noise','5% Noise','10% Noise')
title('Deep LSTM with noisy data')

Теперь постройте график предполагаемых передаточных функций.

figure
plot(time,systemResponse,'k', ...
    time,tfStepNoise1, ...
    time,tfStepNoise5, ...
    time,tfStepNoise10)
grid on
legend('System Response','1% Noise','5% Noise','10% Noise')
title('Transfer functions fitted to noisy data')

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

msefun = @(y,yhat) mean(sqrt((y-yhat).^2)/length(y));

% LSTM errors
lstmMSE(1,:) = msefun(systemResponse,netNoise1Step);
lstmMSE(2,:) = msefun(systemResponse,netNoise5Step);
lstmMSE(3,:) = msefun(systemResponse,netNoise10Step);

% Transfer function errors
tfMSE(1,:) = msefun(systemResponse,tfStepNoise1');
tfMSE(2,:) = msefun(systemResponse,tfStepNoise5');
tfMSE(3,:) = msefun(systemResponse,tfStepNoise10');

mseTbl = array2table([lstmMSE tfMSE],'VariableNames',{'LSTMMSE','TFMSE'})

mseTbl=3×2 table
     LSTMMSE        TFMSE   
    __________    __________

    1.0115e-05    8.8621e-07
    2.5577e-05    9.9064e-06
    5.1791e-05    3.6831e-05

Шум оказывает аналогичный эффект как на результаты LSTM, так и на оценку передаточной функции.

Вспомогательные функции

function [stepSignal,systemResponse,time] = generateStepResponse(model,stepTime,stepAmp,signalDuration)
%Generates a step response for the given model.
%

%Check model type
modelType = class(model);
if nargin < 2
    stepTime = 1;
end

if nargin < 3
    stepAmp = 1;
end

if nargin < 4
    signalDuration = 10;
end

% Constuct step signal
if model.Ts == 0
    Ts = 1e-2;
    time = (0:Ts:signalDuration)';
else
    time = (0:model.Ts:signalDuration)';
end
stepSignal = [zeros(sum(time<=stepTime),1);stepAmp*ones(sum(time>stepTime),1)];
switch modelType
    case {'tf', 'zpk'}
        systemResponse = lsim(model,stepSignal,time);
    case 'idpoly'
        systemResponse = sim(model,stepSignal,time);
    otherwise
        error('Model passed is not supported')
end

stepSignal = stepSignal';
systemResponse = systemResponse';
end