Описание модели

Вестибулоокулярный рефлекс (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)