В этом примере показано, как использовать анализ чувствительности, чтобы сузить количество параметров, которые вы должны оценить, чтобы соответствовать модели. Этот пример использует модель вестибулоглазного рефлекса, который генерирует компенсаторные движения глаз.
Вестибулоокулярный рефлекс (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)
Параметр Delay моделирует тот факт, что существует некоторая задержка в передаче сигналов от внутреннего уха к мозгу и глазам. Эта задержка связана со временем, необходимым для прохождения химическими нейротрансмиттерами синаптических клефтов между нервными камерами. Исходя из количества синапсов, вовлеченных в вестибулоглазный рефлекс, эта задержка, как ожидается, составит около 5 мс. Для целей оценки мы примем, что это между 2 и 9 мс.
Delay = sdo.getParameterFromModel(model_name, 'Delay'); Delay.Value = 0.005; % seconds Delay.Minimum = 0.002; Delay.Maximum = 0.009;
Параметр Gain моделирует тот факт, что глаза движутся не так сильно, как голова. Мы будем использовать 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 моделирует динамику глазодвигательного объекта, т.е. глаз и прикрепленные к нему мышцы и ткани. Объект может быть смоделирована двумя полюсами, однако считается, что шест с большей временной константой отменяется предварительной компенсацией в мозге, чтобы позволить глазу совершать быстрые движения. Таким образом, на графике, когда стимуляция подвергается переходным ребрам, движения глаз следуют лишь с небольшой задержкой. Для параметра 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 или Halton низких несоответствий последовательностей.
Расположенные выше не по диагонали графики показывают графики поля точек между парами различных переменных. Поскольку мы не задали матрицу 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)