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

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)