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

Этот пример использует передаточную функцию четвертого порядка со смешанной быстрой и медленной динамикой и умеренным затуханием. Умеренное затухание вызывает системную динамику к ослаблению за более длительный период времени и показывает способность сети 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 дБ, его значения DC.

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