В этом примере показано, как использовать анализ чувствительности, чтобы сузить количество параметров, которые необходимо оценить, чтобы подобрать модель. Этот пример использует модель vestibulo-окулярного отражения, которое генерирует компенсационные движения глаз.
Vestibulo-окулярное отражение (VOR) позволяет глазам переместиться на той же скорости и в противоположном направлении как голова, так, чтобы видение не было размыто, когда голова перемещается во время нормального действия. Например, если голова поворачивается в одном направлении, повороте глаз в противоположном направлении, с той же скоростью. Это происходит даже в темноте. На самом деле VOR наиболее легко характеризуется измерениями в темноте, чтобы гарантировать, что движения глаз преимущественно управляются VOR.
Файл sdoVOR_Data.mat содержит однородно выборочные данные стимуляции и движений глаз. Если бы VOR были совершенно компенсационными, то график данных о движении глаз, когда инвертировано вертикально, наложил бы точно сверху графика главных данных о движении. Такая система была бы описана усилением 1 и фазой 180 градусов. Однако, когда мы отображаем данные на графике в файле sdoVOR_Data.mat, движения глаз являются близкими, но не совершенно компенсационными.
load sdoVOR_Data.mat; % Column vectors: Time HeadData EyeData figure plot(Time, HeadData, ':b', Time, EyeData, '-g') xlabel('Time (sec)') ylabel('Angular Velocity (deg/sec)') ylim([-110 110]) legend('Head Data', 'Eye Data')
Данные о движении глаз отлично не накладывают главные данные о движении, и это может быть смоделировано несколькими факторами. Главное вращение обнаруживается органами во внутренних ушах, известных как полукруглые каналы. Они обнаруживают главное движение и передают сигналы о главном движении к мозгу, который отправляет моторные команды в глазные мышцы, так, чтобы движения глаз компенсировали главное движение. Мы хотели бы использовать эти данные о движении глаз, чтобы оценить параметры в моделях для этих различных этапов. Модель, которую мы будем использовать, показывают ниже. В модели существует четыре параметра: Delay, Gain, Tc, и Tp.
model_name = 'sdoVOR';
open_system(model_name)
Параметр Задержки моделирует то, что существует некоторая задержка передачи сигналов от внутреннего уха до мозга и глаз. Эта задержка происходит из-за времени, необходимого для химических нейромедиаторов, чтобы пересечь синаптические расселины между нервными клетками. На основе количества синапсов, вовлеченных в vestibulo-окулярное отражение, эта задержка, как ожидают, составит приблизительно 5 мс. В целях оценки мы примем, что это между 2 и 9 мс.
Delay = sdo.getParameterFromModel(model_name, 'Delay'); Delay.Value = 0.005; % seconds Delay.Minimum = 0.002; Delay.Maximum = 0.009;
Параметр Усиления моделирует то, что глаза не перемещаются вполне так, как голова делает. Мы будем использовать 0.8 в качестве нашего исходного предположения и принимать, что это между 0,6 и 1.
Gain = sdo.getParameterFromModel(model_name, 'Gain');
Gain.Value = 0.8;
Gain.Minimum = 0.6;
Gain.Maximum = 1;
Параметр Tc моделирует динамику, сопоставленную с полукруглыми каналами, а также некоторой дополнительной нейронной обработкой. Каналы являются фильтрами высоких частот, потому что после того, как предмет был помещен во вращательное движение, нейронным образом активные мембраны в каналах медленно возвращаются к покоящемуся положению, таким образом, каналы прекращают обнаруживать движение. Таким образом в графике выше, после того, как стимуляция подвергается ребрам перехода, движения глаз имеют тенденцию вылетать от стимуляции в зависимости от времени. На основе механических характеристик каналов, объединенных с дополнительной нейронной обработкой, которая продлевает эту постоянную времени, чтобы улучшить точность VOR, мы оценим, что параметр Tc составляет 15 секунд и принимает, что это между 10 и 30 секундами.
Tc = sdo.getParameterFromModel(model_name, 'Tc');
Tc.Value = 15;
Tc.Minimum = 10;
Tc.Maximum = 30;
Наконец, параметр Tp моделирует динамику oculomotor объекта, i.e. глаз и мышцы и ткани, присоединенные к нему. Объект может быть смоделирован двумя полюсами, однако считается, что полюс с большей постоянной времени отменяется предварительной компенсацией в мозгу, чтобы позволить глазу сделать быстрые перемещения. Таким образом в графике, когда стимуляция подвергается ребрам перехода, движения глаз следуют только с небольшой задержкой. Для параметра Tp мы будем использовать 0,01 секунды в качестве нашего исходного предположения и принимать, что это между 0,005 и 0,05 секундами.
Tp = sdo.getParameterFromModel(model_name, 'Tp');
Tp.Value = 0.01;
Tp.Minimum = 0.005;
Tp.Maximum = 0.05;
Соберите эти параметры в вектор.
v = [Delay Gain Tc Tp];
Создайте объект Experiment. Задайте HeadData как введено.
Exp = sdo.Experiment(model_name); Exp.InputData = timeseries(HeadData, Time);
Объединенные данные о движении глаз с выходом модели.
EyeMotion = Simulink.SimulationData.Signal; EyeMotion.Name = 'EyeMotion'; EyeMotion.BlockPath = [model_name '/Oculomotor Plant']; EyeMotion.PortType = 'outport'; EyeMotion.PortIndex = 1; EyeMotion.Values = timeseries(EyeData, Time);
Добавьте EyeMotion к эксперименту.
Exp.OutputData = EyeMotion;
Используйте характеристики синхронизации данных в модели.
stop_time = Time(end); set_param(gcs, 'StopTime', num2str(stop_time)); dt = Time(2) - Time(1); set_param(gcs, 'FixedStep', num2str(dt))
Создайте сценарий симуляции с помощью эксперимента и получите симулированный выход.
Exp = setEstimatedValues(Exp, v); % use vector of parameters/states
Simulator = createSimulator(Exp);
Simulator = sim(Simulator);
Ищите сигнал model_residual в регистрируемых данных моделирования.
SimLog = find(Simulator.LoggedData, ... get_param(model_name, 'SignalLoggingName') ); EyeSignal = find(SimLog, 'EyeMotion');
Выход модели не совпадает с данными очень хорошо, как показано невязкой, которую мы можем вычислить путем вызывания целевой функции.
estFcn = @(v) sdoVOR_Objective(v, Simulator, Exp, 'Residuals'); Model_Error = estFcn(v); plot(Time, EyeData, '-g', ... EyeSignal.Values.Time, EyeSignal.Values.Data, '--c', ... Time, Model_Error.F, '-r'); xlabel('Time (sec)'); ylabel('Angular Velocity (deg/sec)'); legend('Eye Data', 'Model', 'Residual');
Целевая функция, используемая выше, задана в файле "sdoVOR_Objective.m".
type sdoVOR_Objective.m
function vals = sdoVOR_Objective(v, Simulator, Exp, Method) % Compare model output with data % % Inputs: % v - vector of parameters and/or states % Simulator - used to simulate the model % Exp - Experiment object % Method - 'SSE' for scalar output, 'Residuals' for vector of residuals % Copyright 2014-2015 The MathWorks, Inc. % Requirement setup req = sdo.requirements.SignalTracking; req.Type = '=='; req.Method = Method; % If Residuals requested, keep on same scale as signals, for plotting switch Method case 'Residuals' req.Normalize = 'off'; end % Simulate the model Exp = setEstimatedValues(Exp, v); % use vector of parameters/states Simulator = createSimulator(Exp,Simulator); Simulator = sim(Simulator); % Compare model output with data SimLog = find(Simulator.LoggedData, ... get_param(Exp.ModelName, 'SignalLoggingName') ); OutputModel = find(SimLog, 'EyeMotion'); Model_Error = evalRequirement(req, OutputModel.Values, Exp.OutputData.Values); vals.F = Model_Error;
Создайте объект произвести пространство параметров.
ps = sdo.ParameterSpace([Delay ; Gain ; Tc ; Tp]);
Сгенерируйте 100 выборок от пространства параметров.
rng default; % for reproducibility x = sdo.sample(ps, 100); sdo.scatterPlot(x);
Выборка выше используемых опций по умолчанию, и они отражаются в графиках выше. Значения параметров были выбраны наугад из распределений, которые были универсальны в области значений каждого параметра. Следовательно, графики гистограммы по диагонали кажутся приблизительно универсальными. Если Statistics and Machine Learning Toolbox™ доступен, много других распределений могут использоваться, и выборка может быть сделана с помощью последовательностей низкого несоответствия Sobol или Холтона.
Недиагональные графики выше графиков рассеивания показа между парами различных переменных. Поскольку мы не задавали матрицу RankCorrelation в PS, графики рассеивания не указывают на корреляции. Однако, если параметры, как полагали, коррелировались, это может быть задано с помощью свойства RankCorrelation ps.
Для анализа чувствительности более просто использовать скалярную цель, таким образом, мы зададим сумму квадратичных невязок, "SSE":
estFcn = @(v) sdoVOR_Objective(v, Simulator, Exp, 'SSE');
y = sdo.evaluate(estFcn, ps, x);
Model evaluated at 100 samples.
Оценка могла также быть ускорена с помощью параллельных вычислений.
Получите стандартизированные коэффициенты регрессии.
opts = sdo.AnalyzeOptions;
opts.Method = 'StandardizedRegression';
sensitivities = sdo.analyze(x, y, opts);
Другие типы анализа включают корреляцию и, если Statistics and Machine Learning Toolbox является доступной, частичной корреляцией.
Мы можем просмотреть результаты анализа.
disp(sensitivities)
F
_________
Delay 0.01303
Gain -0.90873
Tc -0.044395
Tp 0.19919
Для стандартизированной регрессии параметры, которые высоко влияют на выход модели, имеют величины чувствительности близко к 1. С другой стороны, менее влиятельные параметры имеют меньшие величины чувствительности. Мы видим, что эта целевая функция чувствительна к изменениям в параметрах Gain и Tp, но намного менее чувствительна к изменениям в параметрах Delay и Tc.
Можно подтвердить результаты анализа чувствительности путем передискретизации и переоценки целевой функции для выборок. Можно также использовать техническую интуицию для быстрого анализа. Например, в этой модели, постоянная времени Tc диапазоны с 10 до 30 секунд. Даже минимальное значение 10 секунд является большим по сравнению с 2 второй длительностью, на которую главная стимуляция движения сохранена при постоянной скорости. Поэтому Tc как ожидают, не будет влиять на выход значительно. Однако, даже когда этот вид интуиции не легко доступен в других моделях, анализ чувствительности может помочь подсветить, какие параметры влияют.
На основе результатов анализа чувствительности определяйте Delay и Tc параметры, как зафиксировано при оптимизации. Это сокращение количества свободных параметров ускоряет оптимизацию.
Delay.Free = false; Tc.Free = false;
Мы можем использовать минимум от анализа чувствительности как исходное предположение для оптимизации.
[fval, idx_min] = min(y.F); Delay.Value = x.Delay(idx_min); Gain.Value = x.Gain(idx_min); Tc.Value = x.Tc(idx_min); Tp.Value = x.Tp(idx_min); % v = [Delay Gain Tc Tp]; opts = sdo.OptimizeOptions; opts.Method = 'fmincon';
Как имел место с оценками модели в анализе чувствительности, параллельные вычисления могли использоваться, чтобы ускорить оптимизацию.
vOpt = sdo.optimize(estFcn, v, opts); disp(vOpt)
Optimization started 27-Jan-2021 14:05:07
max First-order
Iter F-count f(x) constraint Step-size optimality
0 5 13.4798 0
1 18 12.2052 0 0.129 305
2 30 11.1441 0 0.0648 790
3 41 10.0493 0 0.0843 290
4 46 9.23607 0 0.0758 286
5 51 8.76122 0 0.0183 10.1
6 56 8.75862 0 0.00184 0.476
7 57 8.75862 0 8.41e-05 0.476
Local minimum possible. Constraints satisfied.
fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are
satisfied to within the value of the constraint tolerance.
(1,1) =
Name: 'Delay'
Value: 0.0038
Minimum: 0.0020
Maximum: 0.0090
Free: 0
Scale: 0.0078
Info: [1x1 struct]
(1,2) =
Name: 'Gain'
Value: 0.9012
Minimum: 0.6000
Maximum: 1
Free: 1
Scale: 1
Info: [1x1 struct]
(1,3) =
Name: 'Tc'
Value: 16.6833
Minimum: 10
Maximum: 30
Free: 0
Scale: 16
Info: [1x1 struct]
(1,4) =
Name: 'Tp'
Value: 0.0157
Minimum: 0.0050
Maximum: 0.0500
Free: 1
Scale: 0.0156
Info: [1x1 struct]
1x4 param.Continuous
Получите ответ модели после оценки. Ищите сигнал model_residual в регистрируемых данных моделирования.
Exp = setEstimatedValues(Exp, vOpt); Simulator = createSimulator(Exp,Simulator); Simulator = sim(Simulator); SimLog = find(Simulator.LoggedData, ... get_param(model_name, 'SignalLoggingName') ); EyeSignal = find(SimLog, 'EyeMotion');
Сравнение измеренных данных о глазе с оптимизированным ответом модели показывает, что остаточные значения намного меньше.
estFcn = @(v) sdoVOR_Objective(v, Simulator, Exp, 'Residuals'); Model_Error = estFcn(vOpt); plot(Time, EyeData, '-g', ... EyeSignal.Values.Time, EyeSignal.Values.Data, '--c', ... Time, Model_Error.F, '-r'); xlabel('Time (sec)'); ylabel('Angular Velocity (deg/sec)'); legend('Eye Data', 'Model', 'Residual');
Закройте модель
bdclose(model_name)