柔性机翼飞机

在这个例子中，我们研究从明尼苏达大学无人飞行器实验室[1]建造的小型柔性飞翼飞机收集的数据。飞机的几何结构如下图所示。

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

实验装置

实验的目的是收集飞机在外部激励下不同位置的振动响应。这架飞机用一个位于其重心的弹簧悬挂在一个木架上。弹簧具有足够的灵活性，因此弹簧-质量振荡的固有频率不会干扰飞机的基频。输入力通过靠近飞机中心的Unholtz-Dickie 20型电动激振器施加。

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

激振器输入命令被指定为形式的恒幅啁啾输入 $\mathit{一个}\mathrm{罪}（\omega （\mathit{t}）\mathit{t}）$．啁啾频率随时间线性变化，即: $\omega （\mathit{t}）＝{\mathit{c}}_0+{\mathit{c}}_1\mathit{t}$．输入信号覆盖的频率范围为3 - 35hz。数据由两个加速度计收集(前置和后缘加速度计在一个位置)x-location)。因此，我们进行了10次实验来收集所有20个加速度计的响应。加速度计和力传感器的测量都在2000赫兹采样。

数据准备

数据由10组双输出/一输入信号表示，每组包含600K样本。数据可以在MathWorks支持文件站点上获得。万博1manbetx看到免责声明．下载数据文件并将数据加载到MATLAB®工作空间中。

url =“//www.tianjin-qmedu.com/万博1manbetxsupportfiles/ident/flexible-wing-GVT-data/FlexibleWingData.mat”；websave (“FlexibleWingData.mat”url);负载FlexibleWingData.matMeasuredData

的变量MeasuredData结构是否包含字段Exp1，Exp2、……Exp10．每个字段都是一个包含字段的结构yu包含两个响应和相应的输入力数据。绘制第一个实验的数据图。

Fs = 2000;%数据采样频率Ts = 1/fs;%采样时间y = MeasuredData.Exp1.y;%输出数据(2列，每个加速度计一列)u = MeasuredData.Exp1.u;%输入力数据t =(0:长度(u)-1)' * Ts;图subplot(211) plot(t,y) ylabel(“输出(m / s ^ 2)”)传说(“前沿”，“后缘”) subplot(212) plot(t,u)“输入”)包含(的时间(秒)）

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

ExpNames =字段名(MeasuredData);数据= cell(1,10);为k = 1:10 y = MeasuredData.(ExpNames{k}).y;u = MeasuredData.(ExpNames{k}).u;数据{k} = iddata(y, u, Ts，“InterSample”，“提单”）;结束

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

数据=合并(数据{:})

数据=包含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(对于传递函数)，idpoly(AR, ARMA模型)和中的难点(用于状态空间模型)。动态模型可以通过运行估计命令来创建，例如特遣部队而且党卫军在时域或频域对测量数据的命令。

对于这个例子，我们首先将测量的时域数据转换为频响数据的经验传递函数估计使用etfe命令。然后使用估计的频响函数来识别飞机振动响应的状态空间模型。可以直接利用时域数据进行动态模型识别。然而，频响函数形式的数据允许将大型数据集压缩到更少的样本中，并且更容易将估计行为调整到相关频率范围。对频响进行封装idfrd对象。

估计每个数据实验的双输出/单输入频响函数(FRF)。不要使用窗口。使用24000个频率点进行响应计算。

G = cell(1,10);N = 24000;为K = 1:10使用ETFE命令将时域数据转换为频响Data_k = getexp(数据，k);G{k} = etfe(Data_k， []， N);% G{k}是一个@idfrd对象结束

将所有频响连接成单个20输出/ 1输入频响。

G = cat(1, G{:});%连接所有估计频响的输出G.OutputName =“y”；% name输出'y1'， 'y2'，…,“日元”G.InterSample =“提单”；

为了感受估计的频率响应，绘制输出1和15的振幅(任意选择)。放大到感兴趣的频率范围(4-35 Hz)。

Opt = bodeoptions;%绘图选项opt.FreqUnits =“赫兹”；%表示频率，单位为Hzopt.PhaseVisible =“关闭”；%不显示相位OutputNum = [1 15];% pick输出1和15用于绘图clf bodeploy (G(OutputNum，:)， opt)%图频率响应Xlim([4 35])网格在

频响显示至少9个共振频率。为了进行分析，我们将重点放在6- 35hz的频率跨度上，这是飞机最关键的柔性弯曲模式所在。因此，将频响函数减小到该频率区域。

f = g .频率/2/pi;%提取频率矢量，单位为Hz (G存储频率，单位为rad/s)Gs = fselect(G, f>6 & f<=32)% "fselect"选择请求范围内的频响(6.5 - 35hz)

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使用n4sidOptions．

FRF =挤压(Gs.ResponseData);加权=平均值(1./平方根(abs(FRF))).';n4Opt = n4sidOptions(“EstimateCovariance”假的,.．.“WeightingFilter”权重,.．.“OutputWeight”、眼睛(20));sys1 = n4sid(Gs,24，“t”0,“引线”,真的,n4Opt);Fit = sys1.Report.Fit.FitPercent'

适合=1×2054.2689 55.3738 84.0188 85.7017 81.9259 80.2293 55.7448 40.6773 58.6429 76.2997 83.5370 84.6555 85.9240 84.8780 76.8937 81.0399 74.8321 79.4069 65.0803 76.3328

“拟合”数字显示数据(Gs)和模型的(sys1)频率响应采用归一化均方根误差(NRMSE)拟合优度测度。下面是最穷和最适合的情况。

[~，Imin] = min(Fit);[~，Imax] = max(Fit);clf bodeplot (Gs (Imin, Imax,:), sys1 ((Imin, Imax):),选择);Xlim ([6 32])“数据和模型之间的最差(上)和最佳(下)匹配”网格)在传奇(“数据”，“模型”）

拟合结果与模型吻合sys1通过对模型参数进行迭代非线性最小二乘细化，可以进一步提高模型的精度。可以使用党卫军命令。我们设置了评估选项党卫军使用ssestOptions命令。这一次也估计了参数协方差。

ssOpt = ssestOptions(“EstimateCovariance”,真的,.．.“WeightingFilter”, n4Opt。WeightingFilter,.．.“SearchMethod”，“lm”，.．.%使用Levenberg-Marquardt搜索法“显示”，“上”，.．.“OutputWeight”、眼睛(20));sys2 = ssest(Gs, sys1, ssOpt);%估计状态空间模型(这需要几分钟)Fit = sys2.Report.Fit.FitPercent'

适合=1×2088.8715 88.6587 88.9388 89.8630 87.8597 88.1235 79.7016 82.5930 74.2458 82.6593 90.2663 88.7739 89.9033 89.0571 88.8040 88.2251 86.9307 87.5250 83.2249 82.8676

和之前一样，我们画出了最坏和最合适的情况。我们还可视化了模型频率响应的1-sd置信区域。

[~， Imin] = min(Fit);[~， Imax] = max(Fit);clf h = bodeplot (Gs (Imin, Imax,:), sys2 ((Imin, Imax):),选择);showConfidence(h, 1) xlim([6.5 35]) title(“数据和改进模型之间的最差(上)和最佳(下)拟合”网格)在传奇(“数据”，“模型(sys2)”）

细化的状态空间模型sys2在7 - 20hz范围内很好地拟合了频响。大多数共振位置周围的不确定性界限很紧。我们估计了一个24阶模型，这意味着系统中最多有12个振荡模式sys2．使用modalfit命令获取这些模式的固有频率(以Hz为单位)。

f = Gs.Frequency/2/pi;%数据频率(Hz)Fn = modalfit(sys2, f, 12);%固有频率(Hz)disp (fn)

7.7720 7.7954 9.3176 9.4052 9.4919 15.3462 19.3281 23.0219 27.4035 28.6873 31.7036 64.2624

中的值fn建议有两个频率非常接近7.8 Hz，三个频率接近9.4 Hz。检查这些频率附近的频率响应，发现峰值位置在输出中有一点移动。这些差异可以通过更好地控制实验过程来消除，将输入带宽限制在以这些频率为中心的窄范围内进行直接时域识别，或将单个振荡模式拟合到这些频率周围的频率响应。本例中没有探讨这些替代方案。

模态参数计算

我们现在可以使用这个模型了sys2提取模态参数。对频响的检查表明大约有10个模态频率，大致在频率[5 7 10 15 17 23 27 30]Hz附近。的方法可以使这个评估更加具体modalsd命令，当底层模型的顺序发生变化时检查模态参数的稳定性。该操作耗时较长。因此，生成的图形直接作为图像插入。执行下面注释块中的代码以重现图。

FRF =置换(Gs。ResponseData，[3 1 2]);f = Gs.Frequency/2/pi;% {数字pf = modalsd(FRF,f,fs，'MaxModes'，20，'FitMethod'，'lsrf');%}

检查地块和pfValues给出了一个精确的真实固有频率列表:

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

我们使用的值频率向量作为从模型中选择最主要模式的指南sys2．这是使用modalfit命令。

[fn,dr,ms] = modalfit(sys2,f,length(Freq)，“PhysFreq”、频率);

fn为固有频率(Hz)，博士相应的阻尼系数和女士归一化残差作为列向量(每个固有频率一列)。在这些模态参数的提取过程中，只使用了模型的稳定欠阻尼极点。的女士列只包含虚部为正的极点的数据。

模态振型可视化

为了可视化各种弯曲模态，我们需要上面提取的模态参数。此外，我们还需要有关测量点位置的信息。这些位置(x-y值)被记录在矩阵中的每个加速度计AccePos：

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 Hz的模式。

AnimateOneMode(3, fn, ms, AccelPos);

结论

本例展示了一种基于参数模型辨识的模态参数估计方法。利用24阶状态空间模型，在频率6 ~ 32 Hz范围内提取了8个稳定的振荡模态。模态信息与加速度计位置知识相结合，以可视化各种弯曲模态。其中一些模式如下图所示。

参考文献

[1] Gupta, Abhineet, Peter J. Seiler, Brian P. Danowsky。柔性飞翼飞机的地面振动试验美国航空协会大气飞行力学会议，美国航空协会科技论坛。张仁(2016 - 1753)。

函数AnimateOneMode(ModeNum, fn, ModeShapes, AccelPos)动画一个模式。% ModeNum:模式索引。%重新排序传感器位置，以便连续，%不相交的曲线是围绕机体的。PlotOrder = [19: 2:11 1:2:9 10: 2:2, 12:2:20, 19);前进= PlotOrder(1:10);Aft = PlotOrder(20:-1:11);x = AccelPos(PlotOrder,1);y = AccelPos(PlotOrder,2);xAft = AccelPos(Aft,1);yAft = AccelPos(船尾，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 = ThisMode(PlotOrder);zFwd = ThisMode(Fwd);zAft = ThisMode(Aft);CLF col1 = [1 1 1]*0.5;xx =重塑([[xAft, xFwd]';南(2,10)],[2]20);yy =重塑([[yAft, yFwd]';南(2,10)],[2]20);plot3 (x, y, 0 * z0,“- - -”0 * z0, x, y,,“。”， xx(:)， yy(:)， 0 (40,1)，“——”，.．.“颜色”col1,“线宽”, 1“MarkerSize”10“MarkerEdgeColor”, col1);包含(“x”) ylabel (“y”) zlabel (“振幅”) ht = max(abs(z0));轴([-100 100 10 40 -ht ht])视图(5,55)% d’模式。FN = %s Hz'， ModeNum, num2str(fn(ModeNum)，4)));通过更新z坐标进行动画。持有在Col2 = [0.87 0.5 0];H1 = plot3(x,y,0*z0，“- - -”0 * z0, x, y,,“。”， xx(:)， yy(:)， 0 (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)*Rot1 - imag(z0)*Rot2;zAft_t = real(zAft)*Rot1 - imag(zAft)*Rot2;zFwd_t = real(zFwd)*Rot1 - imag(zFwd)*Rot2;zz =重塑([[zAft_t, zFwd_t]';南(2,10)],[2]20);集(h1 (1),“ZData”z_t)组(h1 (2),“ZData”z_t)组(h1 (3),“ZData”，zz(:)) h2.Vertices(:，3) = z_t;暂停(0.1)结束结束