确定评估的关键参数(代码)
这个示例展示了如何使用灵敏度分析来缩小需要估计以拟合模型的参数数量。这个例子使用前庭-眼反射模型，产生代偿性眼球运动。
模型描述
前庭眼反射(VOR)使眼睛以与头部相同的速度和相反的方向移动，因此在正常活动中头部移动时，视觉不会模糊。例如，如果头部向一个方向转动，眼睛则以相同的速度向相反的方向转动。即使在黑暗中也会发生这种情况。事实上，VOR最容易在黑暗中测量，以确保眼球运动主要由VOR驱动。
该文件sdoVOR_Data.mat包含统一采样的刺激和眼球运动数据。如果VOR是完全补偿的，那么当垂直翻转时，眼动数据图将恰好覆盖在头部运动数据图上。这样的系统可以用增益1和相位180度来描述。但是，当我们绘制文件中的数据时sdoVOR_Data.mat时，眼球运动是接近的，但不是完全代偿的。
负载sdoVOR_Data.mat；%列向量:Time HeadData EyeData图表(Time, HeadData，“b”、时间、眼数据、“g”)包含(的时间(秒)) ylabel (角速度(度/秒)ylim([-110 110])“头数据”，“眼数据”）
眼球运动数据不能完美地覆盖头部运动数据，这可以由几个因素建模。头部转动是由内耳的半规管器官感知的。它们检测头部运动并将头部运动的信号传递给大脑，大脑将运动指令发送给眼部肌肉，因此眼部运动补偿了头部运动。我们想用这些眼动数据来估计这些不同阶段模型中的参数。我们将使用的模型如下所示。模型中有四个参数:延迟，获得，Tc,Tp．
model_name =“sdoVOR”；open_system (model_name)
延迟参数模拟了从内耳向大脑和眼睛传递信号时存在延迟的事实。这种延迟是由于化学神经递质穿过神经细胞之间的突触间隙所需的时间。根据参与前庭-眼反射的突触数量，这种延迟预计在5毫秒左右。出于估计的目的，我们将假设它在2到9毫秒之间。
Delay = 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。Value = 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作为输入。
Exp = sdo.Experiment(model_name);Exp.InputData =时间序列(HeadData，时间);
将眼动数据与模型输出联系起来。
EyeMotion = Simulink.SimulationData.Signal;EyeMotion。Name =“EyeMotion”；EyeMotion。BlockPath = [model_name .' /动眼神经的植物的];EyeMotion。PortType =“输出港”；EyeMotion。PortIndex = 1;EyeMotion。Values =时间序列(EyeData, Time);
添加EyeMotion为了实验。
Exp.OutputData = EyeMotion;
在模型中使用数据的时序特征。
stop_time =时间(结束);set_param (gcs,“StopTime”num2str (stop_time));dt =时间(2)-时间(1);set_param (gcs,“FixedStep”num2str (dt))
利用实验创建仿真场景，并获得仿真输出。
Exp = setEstimatedValues(Exp, v);%使用参数/状态向量模拟器= createSimulator(Exp);模拟器= sim(模拟器);
在记录的模拟数据中搜索model_residual信号。
SimLog = find(模拟器。LoggedData,.．.get_param (model_name“SignalLoggingName”));眼信号= find(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
function vals = sdoVOR_Objective(v, Simulator, Exp, Method) %比较模型输出与数据% %输入:% v -参数和/或状态的向量%模拟器-用于模拟模型% Exp -实验对象%方法- 'SSE'用于标量输出，'残差'用于残差向量%版权所有2014-2015 the MathWorks, Inc. %需求设置req = sdo.requirements.SignalTracking;要求的事情。Type =' ==';要求的事情。方法=方法;%如果要求残差，保持与信号相同的比例，以绘制开关方法案例的残差要求。Normalize = 'off';模拟模型Exp = setEstimatedValues(Exp, v);%使用参数/状态的矢量模拟器=创建模拟器(Exp，模拟器);模拟器= sim(模拟器);比较模型输出与数据SimLog = find(模拟器。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或Halton低差异序列进行抽样。
上面的非对角线图显示了不同变量对之间的散点图。由于我们没有在ps中指定秩相关矩阵，所以散点图并不表示相关性。但是，如果参数被认为是相关的，则可以使用的RankCorrelation属性指定ps．
对于灵敏度分析，使用标量物镜更简单，因此我们将指定误差平方和“SSE”:
stfcn = @(v) sdoVOR_Objective(v，模拟器，Exp，上交所的）;Y = sdo。evaluate(estFcn, ps, x);
模型在100个样本中进行评估。
使用并行计算也可以加快计算速度。
得到标准化回归系数。
opts = sdo.AnalyzeOptions;选择。方法=“StandardizedRegression”；敏感性= sdo。分析(x, y, opts);
其他类型的分析包括相关性，如果统计和机器学习工具箱可用，则包括部分相关性。
我们可以查看分析结果。
disp(敏感)
F _________延迟0.01303增益-0.90873 Tc -0.044395 Tp 0.19919
对于标准化回归，高度影响模型输出的参数的敏感性值接近1。另一方面，影响较小的参数具有较小的敏感量。我们看到这个目标函数对增益和Tp参数的变化很敏感，但对延迟和Tc参数的变化不太敏感。
您可以通过重新采样和重新评估样本的目标函数来验证灵敏度分析结果。您还可以使用工程直觉进行快速分析。例如，在这个模型中，时间常数Tc取值范围是10 ~ 30秒。与头部运动刺激以恒定速度保持的2秒持续时间相比，即使是10秒的最小值也很大。因此,Tc预计不会对产量产生很大影响。然而，即使这种直觉在其他模型中不容易获得，敏感性分析也可以帮助突出哪些参数具有影响。
根据灵敏度分析结果，指定延迟而且Tc参数为优化时固定的。自由参数数量的减少加快了优化速度。
延迟。Free = false;Tc。Free = 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];opts = sdo.OptimizeOptions;选择。方法=“fmincon”；
就像灵敏度分析中的模型评估一样，并行计算可以用来加速优化。
vOpt = sdo。优化(estFcn, v, opts);disp (vOpt)
优化开始31- 8月2022 08:13:18 max一阶Iter f -count f(x)约束步长优化05 13.4798 01 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局部最小可能值。约束满足。Fmincon停止，因为当前步长的大小小于步长公差的值，并且约束被满足到约束公差的值之内。(1,1) =名称:'Delay'值:0.0038最小值:0.0020最大值:0.0090空闲:0刻度:0.0078信息:[1x1 struct](1,2) =名称:'Gain'值:0.9012最小值:0.6000最大值:1空闲:1刻度:1信息:[1x1 struct](1,3) =名称:'Tc'值:16.6833最小值:10最大值:30空闲:0刻度:16信息:[1x1 struct](1,4) =名称:'Tp'值:0.0157最小值:0.0050最大值:0.0500空闲:1刻度:0.0156信息:[1x1 struct] 1x4参数。连续
优化结果可视化
得到估计后的模型响应。在记录的模拟数据中搜索model_residual信号。
Exp = setEstimatedValues(Exp, vOpt);模拟器= createsemulator (Exp，模拟器);模拟器= sim(模拟器);SimLog = find(模拟器。LoggedData,.．.get_param (model_name“SignalLoggingName”));眼信号= find(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 (角速度(度/秒)）;传奇(“眼数据”，“模型”，“残留”）;
关闭模型。
bdclose (model_name)
相关的话题
您也可以从以下列表中选择一个网站: