识别关键参数估计(代码)

打开脚本

这个例子展示了如何使用敏感性分析来缩小参数的数量,你需要估计适合一个模型。这个示例使用vestibulo-ocular反射模型,生成补偿眼球运动。

模型描述

vestibulo-ocular反射(伏尔)使眼睛以同样的速度移动,在相反的方向,所以这一愿景不是模糊当头部移动在正常活动。例如,如果头朝一个方向转,眼睛转向相反的方向,用同样的速度。这种情况甚至在黑暗中。事实上,VOR最容易的特点是测量在黑暗中,以确保眼球运动主要是由于刑事和解。

该文件sdoVOR_Data.mat包含均匀抽样数据的刺激和眼球运动。如果刑事和解完全补偿,那么一个阴谋的眼动数据,当垂直翻转,将完全覆盖一块头顶的运动数据。这样一个系统将被增加1和180度的阶段。然而,当我们把文件中的数据sdoVOR_Data.mat眼球运动是接近,但不完全补偿。

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

眼动数据并不完全覆盖头部的运动数据,这可以通过几个因素建模。头部的转动由内耳器官,感觉到的半规管。这些检测头运动和头部运动到大脑传递信号,电机命令发送给眼部肌肉,所以眼部运动补偿运动。我们想使用这个眼动数据估计模型中的参数这些不同的阶段。我们将使用如下所示的模型。有四个参数模型:延迟,获得,Tc,Tp。

model_name =“sdoVOR”;open_system (model_name)

延迟参数模型的事实有一些延迟通信信号从内耳的大脑和眼睛。这延迟是由于化学神经递质所需的时间穿越神经细胞之间的突触结晶。基于参与vestibulo-ocular反射的神经突触的数量,这延迟预计是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参数模型相关的动力学的半规管,以及一些额外的神经处理。运河是高通滤波器,因为一个主题后放入旋转运动,活跃的膜在运河慢慢放松休息回到位置,所以运河停止运动传感。因此在上面的图中,刺激经历过渡边缘后,眼球运动倾向于离开刺激。运河力学特性的基础上,加上额外的神经处理延长这个时间常数提高刑事和解的准确性,我们将估计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);

比较测量数据初步模拟输出

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

经验= sdo.Experiment (model_name);Exp.InputData = timeseries (HeadData、时间);

将眼动数据与模型输出。

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

添加EyeMotion的实验。

Exp.OutputData = EyeMotion;

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

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

创建一个模拟场景中使用实验,获得模拟输出。

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

寻找model_residual信号记录的仿真数据。

SimLog =找到(Simulator.LoggedData,…get_param (model_name“SignalLoggingName”));EyeSignal =找到(SimLog,“EyeMotion”);

模型输出不匹配的数据很好,残余如图所示,我们可以计算通过调用目标函数。

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

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

类型sdoVOR_Objective.m

函数vals = sdoVOR_Objective (v,模拟器、经验、方法)% % %比较模型输出和数据输入:% v -向量的参数和/或州%模拟器,用于模拟模型%实验——实验对象%方法——“SSE”标量输出残差的残差向量的% 2014 - 2015版权MathWorks, Inc . %要求设置点播= sdo.requirements.SignalTracking;要求的事情。类型=“= =”;要求的事情。方法=方法;%如果残差要求,继续与信号相同的规模,为策划开关方法“残差”要求。正常化=“关闭”;结束%模拟模型实验= setEstimatedValues (Exp, v);%使用向量的参数/州模拟器= createSimulator (Exp模拟器);模拟器= sim(模拟器);%比较模型输出数据SimLog =找到(模拟器。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);

上面的采样使用默认选项,这些反映在上面的情节。从分布参数值是随机选择,统一在每个参数的范围。因此,直方图地块沿对角线出现大约一致。如果统计和机器学习工具箱™是可用的,可以使用许多其他分布,采样可以通过使用Sobol或荷low-discrepancy序列。

上面的非对角的情节显示成对的不同变量之间的散点图。因为我们没有指定一个RankCorrelation矩阵在ps,散点图不显示的相关性。但是,如果参数被认为是相关的,这可以使用RankCorrelation属性指定的ps。

敏感性分析,简单的使用一个标量的目标,所以我们将指定的平方误差的总和,“SSE”:

estFcn = @ (v) sdoVOR_Objective (v,模拟器,经验,上交所的);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秒的最小值是大两秒钟的时间相比头运动刺激是在恒定速度。因此,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 9.23607 0.0843 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信息:[1 x1 struct] 1 x4 param.Continuous

可视化结果的优化

估计后获得模型的响应。寻找model_residual信号记录的仿真数据。

经验= setEstimatedValues (Exp、vOpt);模拟器= createSimulator (Exp、模拟器);模拟器= sim(模拟器);SimLog =找到(Simulator.LoggedData,…get_param (model_name“SignalLoggingName”));EyeSignal =找到(SimLog,“EyeMotion”);

眼睛比较测量数据与优化模型响应表明残差要小得多。

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