识别评估的关键参数(代码)

打开脚本

这个示例展示了如何使用敏感性分析来缩小为适合模型而需要估计的参数数量。这个例子使用了前庭-眼反射的模型，它会产生代偿性的眼球运动。

模型描述

前庭眼反射(vestibulo-ocular reflex, VOR)可以使眼睛与头部以相同的速度和相反的方向移动，使头部在正常活动时移动时，视力不会模糊。例如，如果头部朝一个方向转动，眼睛就会以同样的速度朝相反的方向转动。即使在黑暗中也会发生这种情况。事实上，在黑暗中测量最容易识别VOR的特征，以确保眼球运动主要是由VOR驱动的。

该文件sdoVOR_Data.mat包含均匀采样的刺激和眼球运动数据。如果VOR是完美的补偿，那么当垂直翻转时，眼动数据图就会完全覆盖在头部运动数据图上。这样一个系统可以用增益为1，相位为180度来描述。然而，当我们在文件中绘制数据时sdoVOR_Data.mat在美国，眼球运动是闭合的，但并不是完全补偿性的。

负载sdoVOR_Data.mat；%列向量:Time HeadData EyeData图绘制(时间、HeadData“b”、时间、EyeData“g”)包含(的时间(秒)) ylabel (的角速度(度/秒)) ylim([-110 110])“头数据”，“眼数据”）

眼动数据并不完全覆盖头部运动数据，这可以由几个因素建模。头部的旋转是由内耳的器官，即半规管感知的。这些传感器检测头部运动，并将头部运动的信号发送给大脑，大脑再将运动指令发送给眼部肌肉，这样眼部运动就可以补偿头部运动。我们希望利用眼动数据来估计模型中各个阶段的参数。我们将使用的模型如下所示。模型中有四个参数:延迟，获得，Tc,Tp．

model_name =“sdoVOR”；open_system (model_name)

延迟参数模拟了这样一个事实:从内耳到大脑和眼睛的信号传递存在一定的延迟。这种延迟是由于化学神经递质穿过神经细胞之间的突触间隙所需要的时间。根据前庭眼反射所涉及的突触数量，这种延迟预计在5毫秒左右。出于估计的目的，我们假设它在2到9毫秒之间。

延迟= sdo.getParameterFromModel (model_name,“延迟”）;延迟。值= 0.005;%秒延迟。最小值= 0.002;延迟。最大= 0.009;

增益参数模拟了眼睛不像头部那样移动的事实。我们将使用0.8作为我们最初的猜测，并假设它在0.6和1之间。

获得= sdo.getParameterFromModel (model_name,“获得”）;收益。值= 0.8;收益。最小值= 0.6;收益。最大= 1;

Tc参数模拟了与半规管相关的动力学，以及一些额外的神经处理。这些管道是高通过滤器，因为当一个实验对象进行旋转运动后，管道中的神经活动膜会慢慢放松回到休息的位置，所以管道停止感知运动。因此，在上图中，刺激经过过渡边缘后，随着时间的推移，眼球运动趋向于离开刺激。根据运河的力学特性，结合额外的神经处理，延长这个时间常数，以提高VOR的准确性，我们将Tc参数估计为15秒，并假设它在10 - 30秒之间。

Tc = sdo.getParameterFromModel (model_name,“Tc”）;Tc。值= 15;Tc。最小= 10;Tc。最大= 30;

最后，Tp参数模拟动眼植物的动力学，即眼睛和附着在它上面的肌肉和组织。这种植物可以用两个极点来模拟，然而，人们认为，大脑中的预先补偿抵消了时间常数较大的极点，使眼睛能够快速移动。因此，在情节中，当刺激经历过渡边缘时，眼球运动紧随其后，只有一点延迟。对于Tp参数，我们将使用0.01秒作为我们的初始猜测，并假设它在0.005到0.05秒之间。

Tp = sdo.getParameterFromModel (model_name,“Tp”）;Tp。值= 0.01;Tp。最小值= 0.005;Tp。最大= 0.05;

将这些参数收集到一个向量中。

v =[延迟增益Tc Tp];

比较测量数据和初始模拟输出

创建一个Experiment对象。指定HeadData作为输入。

经验= sdo.Experiment (model_name);Exp.InputData = timeseries(HeadData, Time);

将眼动数据与模型输出相关联。

EyeMotion = 金宝appSimulink.SimulationData.Signal;EyeMotion。Name =“EyeMotion”；EyeMotion。BlockPath = [model_name' /动眼神经的植物的];EyeMotion。PortType =“输出港”；EyeMotion。PortIndex = 1;EyeMotion。Values = timeseries(EyeData, Time);

添加EyeMotion的实验。

Exp.OutputData = EyeMotion;

在模型中使用数据的时间特性。

stop_time =时间(结束);set_param (gcs,“StopTime”num2str (stop_time));dt = Time(2) - Time(1);set_param (gcs,“FixedStep”num2str (dt))

利用实验创建仿真场景，并获得仿真输出。

Exp = setEstimatedValues(Exp, v);%使用参数/状态向量模拟器= createSimulator (Exp);模拟器= sim(模拟器);

在记录的模拟数据中搜索model_residual信号。

SimLog =找到(模拟器。LoggedData,．..get_param (model_name“SignalLoggingName”));EyeSignal =找到(SimLog,“EyeMotion”）;

模型输出与数据的匹配不是很好，如残差所示，我们可以调用目标函数来计算残差。

stfcn = @(v) sdoVOR_Objective(v，模拟器，Exp，)“残差”）;Model_Error = estFcn (v);情节(时间、EyeData“g”，．..EyeSignal.Values。时间、EyeSignal.Values.Data“——c”，．..时间,Model_Error。F,“- r”）;包含(的时间(秒)）;ylabel (的角速度(度/秒)）;传奇(“眼数据”，“模型”，“残留”）;

上面使用的目标函数在文件“sdoVOR_Objective.m”中定义。

类型sdoVOR_Objective.m

函数val = sdoVOR_Objective(v, Simulator, Exp, Method) %比较模型输出和数据% %输入:% v -参数和/或状态的向量% Simulator -用于模拟模型% Exp -实验对象% Method - 'SSE'用于标量输出，'Residuals'用于残差向量% Copyright 2014-2015 the MathWorks, Inc. % requirements setup req = sdo.requirements.SignalTracking;要求的事情。类型=“= =”;要求的事情。方法=方法;%如果要求残差，保持与信号相同的比例，用于标绘开关方法情况“残差”要求。正常化=“关闭”;end %模拟模型Exp = setEstimatedValues(Exp, v);% use vector of parameters/states Simulator = createsemulator (Exp,Simulator);模拟器= sim(模拟器); % 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([延迟;所得;Tc;Tp]);

从参数空间生成100个样本。

rng默认的；%的再现性x = sdo。样本(ps, 100);sdo.scatterPlot (x);

上面的抽样使用了默认选项，这些都反映在上面的图中。参数值是从每个参数范围内均匀的分布中随机选取的。因此，沿对角线的直方图图看起来近似一致。如果Statistics和Machine Learning Toolbox™可用，可以使用许多其他分布，并且可以使用Sobol或Halton低差异序列进行采样。

上面的非对角图显示了不同变量对之间的散点图。由于我们没有在ps中指定RankCorrelation矩阵，散点图并不表明相关性。然而，如果参数被认为是相关的，这可以使用RankCorrelation属性指定ps．

对于敏感性分析，使用标量目标比较简单，因此我们将指定误差的平方和“SSE”:

stfcn = @(v) sdoVOR_Objective(v，模拟器，Exp，)上交所的）;y = sdo。评估(estFcn, ps, x);

模型在100个样本下评估。

使用并行计算还可以加快计算速度。

得到标准化回归系数。

选择= sdo.AnalyzeOptions;选择。方法=“StandardizedRegression”；敏感性= sdo。分析(x, y,选择);

其他类型的分析包括相关性，如果统计学和机器学习工具箱可用的话，还有偏相关性。

我们可以查看分析结果。

disp(敏感)

F _________延迟0.01303增益0.90873 Tc -0.044395 Tp 0.19919

对于标准化回归，对模型输出影响较大的参数的灵敏度值接近于1。另一方面，影响较小的参数具有较小的灵敏度。我们看到，这个目标函数对增益和Tp参数的变化很敏感，但对时延和Tc参数的变化不敏感。

您可以通过重新取样和重新评估样本的目标函数来验证灵敏度分析结果。您还可以使用工程直觉进行快速分析。例如，在这个模型中，时间常数Tc范围从10秒到30秒。即使是10秒的最小值，与头部运动刺激以恒定速度保持的2秒持续时间相比也是很大的。因此,Tc预计不会对产量有很大影响。然而，即使这种直觉在其他模型中不容易得到，敏感性分析可以帮助突出哪些参数是有影响的。

根据灵敏度分析的结果，指定延迟和Tc参数作为固定时优化。自由参数数量的减少加速了优化。

延迟。自由= false;Tc。自由= false;

优化

我们可以用灵敏度分析得到的最小值作为优化的初始猜想。

[fval, idx_min] = min(y.F);延迟。值= x.Delay (idx_min);收益。值= x.Gain (idx_min);Tc。值= x.Tc (idx_min);Tp。值= x.Tp (idx_min);％v =[延迟增益Tc Tp];选择= sdo.OptimizeOptions;选择。方法=“fmincon”；

与灵敏度分析中模型评估的情况一样，并行计算可以加快优化速度。

vOpt = sdo。优化(estFcn, v,选择);disp (vOpt)

优化01 - 9月- 2021年开始16:55:59 max一阶Iter F-count f (x)约束步长最优18 0 5 13.4798 0 1 12.2052 11.1441 0.129 305 2 30 0 0.0648 10.0493 790 3 41 0 0.0843 9.23607 290 4 46 0 0 0.0183 10.1 0.0758 286 5 51 8.76122 6 56 8.75862 0 8.75862 0.00184 0.476 7 57 0 8.41 0.476 e-05局部最小值。约束满足。Fmincon停止，因为当前步长小于步长公差的值，约束满足到约束公差的值之内。(1,1) =名字:“延迟”值:0.0038最低:0.0020最高:0.0090免费:0比例:0.0078信息:[1 x1 struct](1、2)=名字:“获得”值:0.9012最低:0.6000最大:1自由:1规模:1信息:[1 x1 struct](1、3)=名字:Tc的值:16.6833最低:10最大:30自由:0规模:16信息:[1 x1 struct](1,4) =名字:Tp的值:0.0157最低:0.0050最大值:0.0500自由:1比例:0.0156信息:[1x1结构]1x4参数。连续

优化结果可视化

估计后得到模型响应。在记录的模拟数据中搜索model_residual信号。

Exp = setEstimatedValues(Exp, vOpt);模拟器= createSimulator (Exp、模拟器);模拟器= sim(模拟器);SimLog =找到(模拟器。LoggedData,．..get_param (model_name“SignalLoggingName”));EyeSignal =找到(SimLog,“EyeMotion”）;

将测量的眼数据与优化后的模型响应进行比较，发现残差要小得多。

stfcn = @(v) sdoVOR_Objective(v，模拟器，Exp，)“残差”）;Model_Error = estFcn (vOpt);情节(时间、EyeData“g”，．..EyeSignal.Values。时间、EyeSignal.Values.Data“——c”，．..时间,Model_Error。F,“- r”）;包含(的时间(秒)）;ylabel (的角速度(度/秒)）;传奇(“眼数据”，“模型”，“残留”）;