В этом примере показано, как использовать анализ чувствительности для сокращения числа параметров, которые необходимо оценить для соответствия модели. В этом примере используется модель вестибуло-глазного рефлекса, который генерирует компенсаторные движения глаз.
Вестибулокулярный рефлекс (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);
В выборке, приведенной выше, использовались опции по умолчанию, которые отражены на графиках выше. Значения параметров выбирались случайным образом из распределений, которые были однородными в диапазоне каждого параметра. Следовательно, графики гистограммы вдоль диагонали выглядят приблизительно однородными. При наличии Toolbox™ статистики и машинного обучения может использоваться множество других распределений, и выборка может быть выполнена с использованием последовательностей с низким несоответствием по Соболю или Халтону.
Внедиагональные графики выше показывают графики рассеяния между парами различных переменных. Поскольку мы не указали матрицу RureCorrelation в пс, графики рассеяния не указывают корреляции. Однако, если считалось, что параметры коррелируются, это может быть указано с помощью свойства RureCorrelation 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);
К другим типам анализа относятся корреляция и частичная корреляция (при наличии инструментов статистики и машинного обучения).
Мы можем просмотреть результаты анализа.
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)