柔性翼飞机

在这个例子中，我们研究了从在明尼苏达大学无人飞行器实验室建造的小型柔性飞翼飞机收集的数据（1）。

飞机机翼在载荷作用下会发生较大变形。柔性模态频率低于机翼刚性较大的普通飞机。这种柔性设计降低了材料成本，增加了飞机的灵活性和飞行范围。但是，如果不加以控制，柔性模态可能导致灾难性的气动弹性失稳设计有效的控制律来抑制这些不稳定性需要精确地确定机翼的各种弯曲模式。

实验装置

实验的目的是收集飞机在不同位置对外部激励的振动响应。飞机在其重心处使用一个弹簧悬挂在木制框架上。弹簧具有足够的柔性，因此弹簧质量振动的固有频率不会干扰飞机的基本频率。输入力通过飞机中心附近的无螺栓Dickie 20型电动振动筛施加。

二十个PCB-353B16加速计沿机翼跨度放置，以收集振动响应，如下图所示。

振动筛输入命令指定为形式的恒定振幅啁啾输入 $\mathit{A.}\mathrm{罪}\left(\omega \left(\mathit{T}\right)\mathit{T}\right)$.啁啾频率随时间线性变化，即， $\omega \left(\mathit{T}\right)={\mathit{C}}_0+{\mathit{C}}_{1.}\mathit{T}$.输入信号所覆盖的频率范围为3–35 Hz。数据由两个加速计（前缘加速计和后缘加速计同时采集）采集x-位置）一次。因此，进行了10次实验，以收集所有20个加速度计响应。加速度计和力传感器测量值均以2000 Hz的频率采样。

数据准备

数据由10组两个输出/一个输入信号表示，每个信号包含600K样本。数据可在MathWorks支持文件网站上获得。请参阅金宝app免责声明. 下载数据文件并将数据加载到MATLAB®工作区。

网址='//www.tatmou.com/金宝appsupportfiles/ident/flexible-wing-GVT-data/FlexibleWingData.mat'；韦伯斯韦(“FlexibleWingData.mat”，url）；加载FlexibleWingData.mat测量数据

变量测量数据是一个具有字段的结构表1,Exp2, ...,Exp10.每个字段都是包含字段的结构Y和u包含两个响应和相应的输入力数据。绘制第一个实验的数据。

fs=2000；%数据采样频率Ts=1/fs；%采样时间y=测量数据Exp1.y；%输出数据（2列，每个加速计一列）u=测量数据Exp1.u；%输入力数据t=（0:length（u）-1）'*Ts；图形子地块（211）绘图（t，y）ylabel(“输出（m/s^2）”)传奇(“前沿”,“后缘”)子地块（212）地块（t，u）ylabel(“输入”)xlabel('时间（秒）')

为了准备用于模型识别的数据，将数据打包到iddata对象iddata对象是在系统标识工具箱中打包时域数据的标准方法™. 输入信号被视为带限信号。

ExpNames=字段名（MeasuredData）；Data=单元格（1,10）；对于k=1:10y=MeasuredData（ExpNames{k}）.y；u=测量数据（ExpNames{k}）.u；Data{k}=iddata（y，u，Ts，“样本间”,“bl”);终止

将数据集对象合并到一个多实验数据对象中。

数据=合并（数据{：}）

Data=包含10个实验的时域数据集。实验样本采样时间Exp1 600001 0.0005 Exp2 600001 0.0005 Exp3 600001 0.0005 Exp4 600001 0.0005 Exp5 600001 0.0005 Exp6 600001 0.0005 Exp7 600001 0.0005 Exp8 600001 0.0005 Exp9 600001 0.0005 Exp10 600001 0.0005输出单元（如果指定）y1 y2输入单元（如果指定）u1

模型识别

我们希望确定一个动态模型，其频率响应与实际飞机的频率响应尽可能接近。动态模型将系统输入和输出之间的数学关系封装为微分方程或差分方程。动态模型的例子有传递函数和状态空间模型。在系统识别工具箱中，动态模型由以下对象封装：idtf（对于传递函数），Ipoly（适用于AR、ARMA车型）和智能决策支持系统（对于状态空间模型）。动态模型可以通过运行估计命令来创建，例如tfest和ssest时域或频域中测量数据的命令。

对于这个例子，我们首先使用埃特菲命令然后使用估计的FRF识别飞机振动响应的状态空间模型。可以直接使用时域数据进行动态模型识别。然而，FRF形式的数据允许将大型数据集压缩为较少的样本，并且更容易将估计行为调整到相关的频率范围。FRF由以色列国防军物体。

估计每个数据实验的两个输出/一个输入频率响应函数（FRF）。不要使用窗口。使用24000个频率点进行响应计算。

G=单元（1，10）；N=24000；对于k=1:10%使用ETFE命令将时域数据转换为FRFData_k=getexp（Data，k）；G{k}=etfe（Data_k，[]，N）；%G{k}是一个@idfrd对象终止

将所有FRF串联成一个20输出/一输入FRF。

G=cat（1，G{:}）；%连接所有估计FRF的输出G.OutputName=“是的”;%命名输出'y1'、'y2'、…、'y20'G.样本间=“bl”;

为了获得估计频率响应的感觉，绘制输出1和15的振幅（任意拾取）。放大感兴趣的频率范围（4–35 Hz）。

opt=预兆；%绘图选项选择频率单位=“赫兹”;%以Hz为单位显示频率可选择的=“关”;%不显示相位OutputNum=[1 15]；%选择输出1和15进行打印clf bodeplot（G（OutputNum，：），opt）%绘制频率响应曲线xlim（[4 35]）网格在…上

FRF显示了至少9个共振频率。为了进行分析，我们希望将重点放在6-35 Hz的频率范围上，其中飞机最关键的柔性弯曲模式位于该频率范围内。因此，将FRF降低到该频率范围内。

f=G.频率/2/pi；%以Hz为单位提取频率向量（G以rad/s为单位存储频率）Gs=fselect（G，f>6&f<=32）%“fselect”选择请求范围内的FRF（6.5-35 Hz）

Gs=IDFRD模型。包含20个输出和1个输入的频率响应数据。624个频率点的响应数据可用，范围从37.96 rad/s到201.1 rad/s。采样时间：0.0005秒输出通道：“y（1）”，“y（2）”，“y（3）”，“y（4）”，“y（5）”，“y（6）”，“y（7）”，“y（8）”，“y（9）”，“y（10）”，“y（11）”，“y（12）”，“y（13）”，“y（14）”，“y（15）‘y（16）’、‘y（17）’、‘y（18）’、‘y（19）’、‘y（20）’输入通道：‘u1’状态：估计模型

Gs因此包含20个测量位置的频率响应测量值。接下来，确定一个状态空间模型以进行拟合Gs.子空间估计n4sid提供快速的非迭代估计。状态空间模型结构配置如下：

我们设置了以下各项的估算选项：n4sid使用N4SID选项.

FRF=挤压（一般响应数据）；加权=平均值（1./sqrt（绝对值（FRF）））。；n4Opt=N4SIDFOPTIONS(“估计协方差”错误的...“加权过滤器”，加权，...“输出重量”眼（20）；sys1=n4sid（Gs，24，“Ts”,0,“馈通”，对，n4Opt）；Fit=sys1.Report.Fit.FitPercent'

合身=1×2057.0200 57.9879 85.0160 86.3815 80.4879 80.4430 58.2216 45.2692 61.5057 76.7612 84.7305 86.2600 86.4266 85.0251 76.9208 82.1191 74.7982 79.6837 67.9078 76.7249

“拟合”数字显示数据之间的拟合百分比(Gs)模型的(系统1)使用归一化均方根误差（NRMSE）拟合优度度量的频率响应。下一步将绘制最差和最佳拟合。

[~，Imax]=min（Fit）；[~，Imax]=max（Fit）；clf bodeplot（Gs（[Imin，Imax]，：），sys1（[Imin，Imax]，：），opt）；xlim（[6 32]）标题(“数据和模型之间的最差（顶部）和最佳（底部）拟合”)网格在…上传奇(“数据”,“模型”)

模型实现的拟合系统1可以通过迭代非线性最小二乘法对模型参数进行细化来进一步改进ssest命令我们设置了以下各项的估算选项：ssest使用停止这一次还估计了参数协方差。

ssOpt=SSEOPTIONS(“估计协方差”符合事实的...“加权过滤器”，n4Opt.WeightingFilter，...“搜索方法”,“lm”,...%使用Levenberg-Marquardt搜索方法“显示”,“开”,...“输出重量”眼（20）；sys2=ssest（Gs、sys1、ssOpt）；%估计状态空间模型（这需要几分钟）Fit=sys2.Report.Fit.FitPercent'

合身=1×2089.7225 89.5185 89.7260 90.4986 88.5522 88.8727 81.3225 83.5975 75.9215 83.1763 91.1358 89.7310 90.6844 89.8444 89.6685 89.1467 87.8532 88.0385 84.2898 83.3578

和以前一样，我们设计了最差和最适合的方案。我们还可视化了模型频率响应的1-sd置信域。

[~，Imax]=min（Fit）；[~，Imax]=max（Fit）；clf h=bodeplot（Gs（[Imin，Imax]，：），sys2（[Imin，Imax]，：），opt）；showConfidence（h，1）xlim（[6.5 35]）标题(“数据和优化模型之间的最差（顶部）和最佳（底部）拟合”)网格在…上传奇(“数据”,'模型（sys2）')

精化状态空间模型系统2在7–20 Hz范围内，FRF非常适合。大多数共振位置的不确定度范围都很紧。我们估计了24阶模型，这意味着系统中最多有12个振荡模式系统2.使用莫达尔菲特命令获取这些模式的固有频率（以Hz为单位）。

f=总频率/2/pi；%数据频率（Hz）fn=modalfit（sys2，f，12）；%固有频率（Hz）disp（fn）

7.7721 7.7953 9.3147 9.4009 9.4910 15.3463 19.3291 23.0219 27.4164 28.7256 31.7014 63.3034

中的值fn建议两个频率非常接近7.8 Hz，三个频率接近9.4 Hz。对这些频率附近的频率响应的检查表明，峰值位置在输出端稍有偏移。可以通过更好地控制实验过程、执行直接时域识别（输入带宽限制在以这些频率为中心的窄范围内）或将单个振荡模式与这些频率周围的频率响应相匹配来消除这些差异。本例中未探讨这些备选方案。

模态参数计算

我们现在可以使用这个模型了系统2提取模态参数。对FRF的检查表明大约有10个模态频率，大致在频率[5 7 10 15 17 23 27 30]Hz左右。我们可以使用莫达尔斯德命令，当基础模型的顺序改变时，检查模态参数的稳定性。此操作需要很长时间。因此，生成的绘图直接作为图像插入。在下面的注释块中执行代码以再现图形。

FRF=排列（一般响应数据[3 1 2]）；f=总频率/2/pi；%{图形pf=modalsd（FRF，f，fs，'MaxModes'，20，'FitMethod'，'lsrf'）；%}

检查地块和pf这些值表示真实固有频率的精确列表：

频率=[7.8 9.4 15.3 19.3 23.0 27.3 29.2 31.7]；

我们使用中的值频率矢量作为从模型中选择最主要模式的指南系统2. 这是使用莫达尔菲特命令

[fn，dr，ms]=modalfit（sys2，f，length（Freq），“PhysFreq”，Freq）；

fn是固有频率（单位为Hz），博士相应的阻尼系数和太太归一化残差作为列向量（每个固有频率对应一列）。在提取这些模态参数的过程中，仅使用模型的稳定、欠阻尼极点太太列仅包含具有正虚部的极点的数据。

振型可视化

为了可视化各种弯曲模式，我们需要上面提取的模态参数。此外，我们需要有关测量点位置的信息。这些位置（x-y值）记录在矩阵中的每个加速计上阿克普斯:

AccelPos=[...%参见图216.63 18.48;%中心最近右16.63 24.48; 27.90 22.22; 27.90 28.22; 37.17 25.97; 37.17 31.97; 46.44 29.71; 46.44 35.71; 55.71 33.46; 55.71 39.46;%最右边-16.63 18.48;%离中心最近的左边-16.63 24.18; -27.90 22.22; -27.90 28.22; -37.17 25.97; -37.17 31.97; -46.44 29.71; -46.44 35.71; -55.71 33.46; -55.71 39.46];%最左边

振型包含在矩阵中太太其中每一列对应于一个频率的形状，通过在传感器坐标上叠加振型振幅并改变模式固有频率上的振幅来激活模式。动画显示没有阻尼的弯曲。作为例子，考虑在15.3赫兹的模式。

动画表情（3，fn，ms，AccelPos）；

结论

该示例显示了一种基于参数模型识别的模态参数估计方法。使用24阶状态空间模型，提取了频率范围为6–32 Hz的8个稳定振荡模式。将模态信息与加速计位置知识相结合，以可视化各种弯曲模式。一些这些模式如下图所示。

工具书类

[1] 古普塔、阿比内特、彼得·J·塞勒和布赖恩·P·达诺夫斯基。“柔性飞翼飞机的地面振动试验。”AIAA大气飞行力学会议，AIAA科学技术论坛。（AIAA 2016-1753）。

作用动画表情（ModeNum、fn、ModeShapes、AccelPos）%设置一个模式的动画。%ModeNum：模式的索引。%重新排列用于打印的传感器位置，以便，%非相交曲线围绕机身进行跟踪。PlotOrder=[19:-2:11,1:2:9,10:-2:2,12:2:20,19]；Fwd=PlotOrder（1:10）；Aft=PlotOrder（20:-1:11）；x=AccelPos（PlotOrder，1）；y=AccelPos（PlotOrder，2）；xAft=AccelPos（Aft，1）；yAft=AccelPos（Aft，2）；xFwd=AccelPos（Fwd，1）；yFwd=AccelPos（Fwd，2）；wn=fn（ModeNum）*）*2*pi；%模式频率，单位为rad/secT=1/fn（ModeNum）；%模态振荡周期Np=2.25；%要模拟的周期数tmax=Np*T；%模拟Np周期Nt=100；%动画的时间步数t=linspace（0，tmax，Nt）；ThisMode=ModeShapes（：，ModeNum）/max（abs（ModeShapes（：，ModeNum））；%标准化绘图z0=该模式（绘图顺序）；zFwd=该模式（Fwd）；zAft=该模式（后）；clf col1=[1]*0.5；xx=重塑（[[xAft，xFwd]'；NaN（2,10）]，[2,20]）；yy=重塑（[[yAft，yFwd]'；NaN（2,10）]，[2,20]）；图3（x，y，0*z0，'-'，x，y，0*z0，'.'，xx（：），yy（：），零（40,1），'--',...“颜色”，col1，“线宽”1.“MarkerSize”,10,“MarkerEdgeColor”，col1）；xlabel(“x”)伊拉贝尔(“是的”)兹拉贝尔(“振幅”)ht=最大值（绝对值（z0））；轴（[-100 100 10 40-ht]）视图（5,55）标题（sprintf('模式%d.FN=%s Hz'，ModeNum，num2str（fn（ModeNum），4））；%通过更新z坐标设置动画。持有在…上col2=[0.870.50]；h1=图3（x，y，0*z0，'-'，x，y，0*z0，'.'，xx（：），yy（：），零（40,1），'--',...“颜色”，col2，“线宽”1.“MarkerSize”,10,“MarkerEdgeColor”，col2）；h2=fill3（x，y，0*z0，col2，“FaceAlpha”,0.6); 持有关对于k=1:Nt Rot1=cos（wn*t（k））；Rot2=sin（wn*t（k））；z_t=real（z0）*Rot2；zAft_t=real（zAft）*Rot1-imag（zAft）*Rot2；zFwd_t=real（zFwd）*Rot1-imag（zFwd Rot2；zz=重塑（[[zAft_t，zFwd t）]；NaN（2,10，[220]）；set（h1（1），“ZData”，z_t）集（h1（2），“ZData”，z_t）集（h1（3），“ZData”，zz（：）h2.顶点（：，3）=z_t；暂停（0.1）终止终止