CN113189866B - 一种倾转旋翼无人机状态空间模型辨识方法 - Google Patents

一种倾转旋翼无人机状态空间模型辨识方法 Download PDF

Info

Publication number
CN113189866B
CN113189866B CN202110148131.4A CN202110148131A CN113189866B CN 113189866 B CN113189866 B CN 113189866B CN 202110148131 A CN202110148131 A CN 202110148131A CN 113189866 B CN113189866 B CN 113189866B
Authority
CN
China
Prior art keywords
space model
state space
equal
parameter
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110148131.4A
Other languages
English (en)
Other versions
CN113189866A (zh
Inventor
赵江
赵艳琪
蔡志浩
王英勋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202110148131.4A priority Critical patent/CN113189866B/zh
Publication of CN113189866A publication Critical patent/CN113189866A/zh
Application granted granted Critical
Publication of CN113189866B publication Critical patent/CN113189866B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开了一种倾转旋翼无人机状态空间模型辨识方法，基于倾转旋翼无人机试验数据，通过频域辨识，建立倾转旋翼无人机的状态空间模型，并通过对模型参数精确度进行分析以对模型结构进行优化，从而提高模型的精确度和可靠性，增强模型的鲁棒性。为了使模型的辨识结果能够更好地匹配频率响应数据，需要对辨识参数精确度或相对置信度指标进行分析，以便改进模型结构，通过模型优化使得辨识得到的模型在满足输入输出最大拟合的基础上没有被过度参数化，从而增加模型的精确度和可靠性。利用上述方法，可以得到倾转旋翼无人机不同飞行模式和飞行状态下的状态空间模型，从而应用于飞行控制系统设计，操纵品质评估和仿真验证。

Description

一种倾转旋翼无人机状态空间模型辨识方法
技术领域
本发明涉及无人机飞行控制与建模技术领域，尤其涉及一种倾转旋翼无人 机状态空间模型辨识方法。
背景技术
倾转旋翼无人机结合了直升机和固定翼飞机各自的优势，兼具起降要求低 和续航能力强的特点。倾转旋翼无人机与传统布局飞行器相比，面临的飞行力 学和飞行控制等问题更加复杂，因此，客观准确的数学模型对于研究其飞行动 态特性从而进一步研发以及设计飞行控制系统而言非常重要。
对于涉及到高阶甚至耦合的复杂MIMO系统，状态空间模型比传递函数 模型更能反映MIMO系统内部的动态特性，并且，很多应用中要求MIMO系 统辨识的最终结果是状态空间模型，其中，状态空间模型参数包括稳定导数和 控制导数，甚至一些物理系统参数，此外，对于很多现代MIMO飞行控制系 统的设计方法，如线性二次型调节器(LQR)和状态反馈H∞最优控制或动态 逆(NDI)，都是基于状态空间理论，要求得到精确的状态空间模型。
现有的状态空间模型辨识方法的研究重点主要在于辨识方法的选择，致力 于提高辨识精度和辨识算法收敛的速度。现有的状态空间模型辨识方法虽然能 够得到与输入输出数据拟合度很高的辨识结果，但是不能保证状态空间模型具 有可预测能力，且模型参数的物理意义也不明确。
发明内容
有鉴于此，本发明提供了一种倾转旋翼无人机状态空间模型辨识方法，用 以得到精度符合要求的状态空间模型，并保证模型的可预测能力，明确模型参 数的物理意义，以便于控制器设计。
本发明提供的一种倾转旋翼无人机状态空间模型辨识方法，包括如下步骤：
S1：建立倾转旋翼无人机的运动微分方程；
S2：将所述运动微分方程转换为状态空间模型，并根据先验知识确定所述 状态空间模型的参数初值；
S3：针对倾转旋翼无人机的飞行特性，设计控制器和辨识方案；
S4：根据所述辨识方案，在所述控制器控制下进行试飞试验，采集试验数 据，从中选取与所述状态空间模型的待辨识参数有关的时间历程数据，将选取 的时间历程数据通过线性调节Z变换计算频率响应，对所述频率响应进行处理， 生成频率响应数据；
S5：利用所述频率响应数据，通过优化算法辨识所述状态空间模型的参数； 判断所述状态空间模型的单输入单输出代价函数是否小于或等于第一阈值，且 所述状态空间模型的平均代价函数是否小于或等于第二阈值；若是，则执行步 骤S6；若否，则返回步骤S2，重复执行步骤S2～步骤S5；
S6：计算待辨识参数的克拉美-劳边界和不敏感度，判断待辨识参数的克 拉美-劳边界是否小于或等于第三阈值，且待辨识参数的不敏感度是否小于或 等于第四阈值；若是，则结束操作，辨识出参数，得到状态空间模型；若否， 则执行步骤S7；
S7：移除具有最大不敏感度的参数，返回步骤S5，重新进行辨识和判断； 若所述状态空间模型的单输入单输出代价函数大于第一阈值，且所述状态空间 模型的平均代价函数大于第二阈值，则将移除的参数重新插入所述状态空间模 型，执行步骤S8；若所述状态空间模型的单输入单输出代价函数小于或等于第 一阈值，且所述状态空间模型的平均代价函数小于或等于第二阈值，则进一步 判断所述状态空间模型的单输入单输出代价函数的变化量是否大于或等于第 五阈值，且所述状态空间模型的平均代价函数的变化量是否大于或等于第六阈 值，若是，则将移除的参数重新插入所述状态空间模型，执行步骤S8，若否， 则返回步骤S6，重新计算剩余参数的克拉美-劳边界和不敏感度，若剩余参数 的不敏感度小于或等于第四阈值，则执行步骤S8，若剩余参数的不敏感度大于 第四阈值，则返回步骤S7，继续移除具有最大不敏感度的参数；
S8：移除具有最大克拉美-劳边界的参数，返回步骤S5，重新进行辨识和 判断；若所述状态空间模型的单输入单输出代价函数大于第一阈值，且所述状 态空间模型的平均代价函数大于第二阈值，则将移除的参数重新插入所述状态 空间模型，结束操作；若所述状态空间模型的单输入单输出代价函数小于或等 于第一阈值，且所述状态空间模型的平均代价函数小于或等于第二阈值，则进 一步判断所述状态空间模型的单输入单输出代价函数的变化量是否大于或等 于第五阈值，且所述状态空间模型的平均代价函数的变化量是否大于或等于第 六阈值，若是，则将移除的参数重新插入所述状态空间模型，结束操作，若否， 则返回步骤S6，重新计算剩余参数的克拉美-劳边界和不敏感度，若剩余参数 的克拉美-劳边界小于或等于第三阈值，则结束操作，若剩余参数的克拉美-劳 边界大于第三阈值，则返回步骤S8，继续移除具有最大克拉美-劳边界的参数。
在一种可能的实现方式中，在本发明提供的上述倾转旋翼无人机状态空间 模型辨识方法中，步骤S3中，辨识方案，具体包括：
在控制空速和短舱倾角保持不变的情况下，分别进行俯仰通道、滚转通道 及偏航通道的扫频试验2～3次，待辨识通道称为主通道，其余称为副通道；采 用扫频信号作为主通道的输入信号，扫频幅值引起的飞机角度响应范围为-a～a， a的取值范围为5°～15°，频率范围为0.5ωBW≤ω≤2.5ω180，ωBW表示带宽频率，ω180表示相角为-180°时所对应的频率；扫频试验在倾转旋翼无人机处于配平状态 3s～5s后开始，扫频结束后在配平状态停留至少3s后结束一次扫频试验，每次 扫频试验的时间
Figure BDA0002931015760000031
ωmin表示最小频率，ωmin＝0.5ωBW；在副通道输 入幅值为主通道输入幅值20％的白噪声，记录控制器的输出量和所述状态空间 模型的输出量进行辨识。
在一种可能的实现方式中，在本发明提供的上述倾转旋翼无人机状态空间 模型辨识方法中，步骤S4中，从中选取与所述状态空间模型的待辨识参数有 关的时间历程数据，将选取的时间历程数据通过线性调节Z变换计算频率响应， 对所述频率响应进行处理，生成频率响应数据，具体包括：
截取扫频试验以及扫频试验前后各3s的时间历程数据，将截取的时间历 程数据通过线性调节Z变换计算频率响应，经过多输入规整和复合窗口组合将 所述频率响应转换为频率响应数据，剔除所述频率响应数据中频域响应相干值
Figure BDA0002931015760000044
或相干值曲线振荡的数据，剩余的频率响应数据用于辨识。
在一种可能的实现方式中，在本发明提供的上述倾转旋翼无人机状态空间 模型辨识方法中，步骤S5，利用所述频率响应数据，通过优化算法辨识所述状 态空间模型的参数；判断所述状态空间模型的单输入单输出代价函数是否小于 或等于第一阈值，且所述状态空间模型的平均代价函数是否小于或等于第二阈 值；若是，则执行步骤S6；若否，则返回步骤S2，重复执行步骤S2～步骤S5， 具体包括：
利用所述频率响应数据，采用正割法结合所述状态空间模型对频域数据进 行参数拟合；
单输入单输出代价函数定义如下：
Figure BDA0002931015760000041
其中，||表示每个频率ω处的幅值，∠表示每个频率ω处的相位，nω表示 频率采样点的数量，ω1表示拟合的起始频率，
Figure BDA0002931015760000045
表示拟合的终止频率；
Figure BDA0002931015760000042
Wγ表示加权函数；Wg表示幅值平方误差的相对权重，Wp表 示相位平方误差的相对权重，
Figure BDA0002931015760000043
表示复合频率响应估计矩阵，T表示辨识得到 的状态空间模型的频率响应矩阵；
平均代价函数为：
Figure BDA0002931015760000051
Figure BDA0002931015760000052
若所述状态空间模型的单输入单输出代价函数满足J1≤200且平均代价函 数满足Jave≤100，则执行步骤S6；否则，返回步骤S2，重复执行步骤S2～步骤 S5。
在一种可能的实现方式中，在本发明提供的上述倾转旋翼无人机状态空间 模型辨识方法中，步骤S6，计算待辨识参数的克拉美-劳边界和不敏感度，判 断待辨识参数的克拉美-劳边界是否小于或等于第三阈值，且待辨识参数的不 敏感度是否小于或等于第四阈值；若是，则结束操作，辨识出参数，得到状态 空间模型；若否，则执行步骤S7，具体包括：
第i个辨识参数的克拉美-劳边界CRi由Hessian矩阵H的逆阵中的对角线 元素确定：
Figure BDA0002931015760000053
其中，np×np阶矩阵H定义为：
Figure BDA0002931015760000054
第i个辨识参数的克拉美-劳边界CRi由收敛辨识参数值的百分比形式表示 为：
Figure BDA0002931015760000055
第i个辨识参数的不敏感度Ii由Hessian矩阵H中的对角线元素确定：
Figure BDA0002931015760000061
第i个辨识参数的不敏感度Ii由收敛辨识参数值的化百分比形式表示为：
Figure BDA0002931015760000062
若第i个辨识参数的克拉美-劳边界满足
Figure BDA0002931015760000063
且不敏感度满足
Figure BDA0002931015760000064
则结束操作，辨识出参数，得到状态空间模型；否则，执行步骤S7。
在一种可能的实现方式中，在本发明提供的上述倾转旋翼无人机状态空间 模型辨识方法中，步骤S7和步骤S8中，所述第五阈值的取值范围为1～2，所 述第六阈值的取值范围为10～20。
本发明提供的上述倾转旋翼无人机状态空间模型辨识方法，基于倾转旋翼 无人机试验数据，通过频域辨识，建立倾转旋翼无人机的状态空间模型，并通 过对模型参数精确度进行分析以对状态空间模型的结构进行优化，从而提高状 态空间模型的精确度和可靠性，增强状态空间模型的鲁棒性。为了使状态空间 模型的辨识结果能够更好地匹配频率响应数据，需要对辨识参数精确度或相对 置信度指标进行分析，以便改进状态空间模型的结构，通过模型优化使得辨识 得到的状态空间模型在满足输入输出最大拟合的基础上没有被过度参数化，从 而增加状态空间模型的精确度和可靠性。利用本发明提供的上述倾转旋翼无人 机状态空间模型辨识方法，可以得到倾转旋翼无人机不同飞行模式和飞行状态下的状态空间模型，从而应用于飞行控制系统设计，操纵品质评估和仿真验证。
附图说明
图1为本发明实施例1中倾转旋翼无人机状态空间模型辨识方法的流程示 意图；
图2为本发明实施例1中PID控制器的结构框图；
图3为本发明实施例1中滚转通道辨识结果曲线图；
图4为本发明实施例1中俯仰通道辨识结果曲线图。
具体实施方式
下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进 行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限 制本发明。
本发明提供的一种倾转旋翼无人机状态空间模型辨识方法，包括如下步骤：
S1：建立倾转旋翼无人机的运动微分方程；
S2：将运动微分方程转换为状态空间模型，并根据先验知识确定状态空间 模型的参数初值；
S3：针对倾转旋翼无人机的飞行特性，设计控制器和辨识方案；
S4：根据辨识方案，在控制器控制下进行试飞试验，采集试验数据，从中 选取与状态空间模型的待辨识参数有关的时间历程数据，将选取的时间历程数 据通过线性调节Z变换计算频率响应，对频率响应进行处理，生成频率响应数 据；
S5：利用频率响应数据，通过优化算法辨识状态空间模型的参数；判断状 态空间模型的单输入单输出代价函数是否小于或等于第一阈值，且状态空间模 型的平均代价函数是否小于或等于第二阈值；若是，则执行步骤S6；若否，则 返回步骤S2，重复执行步骤S2～步骤S5；
S6：计算待辨识参数的克拉美-劳边界和不敏感度，判断待辨识参数的克 拉美-劳边界是否小于或等于第三阈值，且待辨识参数的不敏感度是否小于或 等于第四阈值；若是，则结束操作，辨识出参数，得到状态空间模型；若否， 则执行步骤S7；
S7：移除具有最大不敏感度的参数，返回步骤S5，重新进行辨识和判断； 若状态空间模型的单输入单输出代价函数大于第一阈值，且状态空间模型的平 均代价函数大于第二阈值，则将移除的参数重新插入状态空间模型，执行步骤S8；若状态空间模型的单输入单输出代价函数小于或等于第一阈值，且状态空 间模型的平均代价函数小于或等于第二阈值，则进一步判断状态空间模型的单 输入单输出代价函数的变化量是否大于或等于第五阈值，且状态空间模型的平 均代价函数的变化量是否大于或等于第六阈值，若是，则将移除的参数重新插 入状态空间模型，执行步骤S8，若否，则返回步骤S6，重新计算剩余参数的克拉美-劳边界和不敏感度，若剩余参数的不敏感度小于或等于第四阈值，则 执行步骤S8，若剩余参数的不敏感度大于第四阈值，则返回步骤S7，继续移 除具有最大不敏感度的参数；
S8：移除具有最大克拉美-劳边界的参数，返回步骤S5，重新进行辨识和 判断；若状态空间模型的单输入单输出代价函数大于第一阈值，且状态空间模 型的平均代价函数大于第二阈值，则将移除的参数重新插入状态空间模型，结 束操作；若状态空间模型的单输入单输出代价函数小于或等于第一阈值，且状 态空间模型的平均代价函数小于或等于第二阈值，则进一步判断状态空间模型 的单输入单输出代价函数的变化量是否大于或等于第五阈值，且状态空间模型 的平均代价函数的变化量是否大于或等于第六阈值，若是，则将移除的参数重 新插入状态空间模型，结束操作，若否，则返回步骤S6，重新计算剩余参数的 克拉美-劳边界和不敏感度，若剩余参数的克拉美-劳边界小于或等于第三阈值， 则结束操作，若剩余参数的克拉美-劳边界大于第三阈值，则返回步骤S8，继 续移除具有最大克拉美-劳边界的参数。
下面通过一个具体的实施例结合图1对本发明提供的上述倾转旋翼无人机 状态空间模型辨识方法的具体实施进行详细说明。
实施例1：
第一步，建立倾转旋翼无人机的运动微分方程。
对于倾转旋翼无人机配平飞行状态，可以对运动方程进行小扰动线性化， 得到运动方程为：
Figure BDA0002931015760000091
其中，欧拉角速度为：
Figure BDA0002931015760000092
其中，u、v、w分别表示体轴系下x、y、z轴方向的扰动速度，p、q、r 分别表示扰动滚转角速度、扰动俯仰角速度和扰动偏航角速度，φ、θ、ψ分别 表示扰动滚转角、扰动俯仰角和扰动偏航角；
Figure BDA0002931015760000093
分别表示体轴系下x、 y、z轴方向的扰动速度对时间的导数，
Figure BDA0002931015760000094
分别表示扰动滚转角速度、 扰动俯仰角速度和扰动偏航角速度对时间的导数，
Figure BDA0002931015760000095
分别表示扰动滚转 角、扰动俯仰角和扰动偏航角对时间的导数；δa、δe为控制量，分别表示副翼 偏转角和升降舵偏转角；Xu、Yu、Zu、Lu、Mu、Nu分别表示纵向力X、横向 力Y、侧向力Z、滚转力矩L、俯仰力矩M和偏航力矩N对u的偏导数，Xv、Yv、Zv、Lv、Mv、Nv分别为X、Y、Z、L、M和N对v的偏导数，Xq、Yq、Zq、Mq分别为X、Y、Z、M对q的偏导数，
Figure BDA0002931015760000096
分别为X对δa、δe的偏导数，
Figure BDA0002931015760000097
分别为Y对δa、δe的偏导数，
Figure BDA0002931015760000098
分别为Z对δa、δe的偏导数,Lp、
Figure BDA0002931015760000099
分别为L对p、δa、δe的偏导数，Mw、
Figure BDA00029310157600000912
分别为M对w、δa、 δe的偏导数，Nr、
Figure BDA00029310157600000913
分别为N对r、δa、δe的偏导数；U0、V0、W0、Θ0分别表示配平状态下体轴系下x、y、z轴方向速度和俯仰角；g表示重力加速 度。
第二步，将运动微分方程转换为状态空间模型，并根据先验知识确定状态 空间模型的参数初值。
整理成状态空间的形式如下：
Figure BDA0002931015760000101
其中，
Figure BDA0002931015760000102
为控制量，分别为旋翼总距差动、旋翼纵 向周期变距、旋翼纵向周期变距差动、旋翼总距和方向舵偏转角；Xp、Yp、Zp、 Mp、Np分别为X、Y、Z、M和N对p的偏导数；Xw、Yw、Zw、Lw、Nw分别 为X、Y、Z、L和N对w的偏导数，Lq、Nq分别为L和N对q的偏导数，Xr、 Yr、Zr、Lr、Mr分别为X、Y、Z、L、M对r的偏导数，
Figure BDA0002931015760000103
Figure BDA0002931015760000104
分别为X对
Figure BDA0002931015760000105
δr的偏导数，
Figure BDA0002931015760000106
分 别为Y对
Figure BDA0002931015760000107
δr的偏导数，
Figure BDA0002931015760000108
分别 为Z对
Figure BDA0002931015760000109
δr的偏导数,
Figure BDA00029310157600001010
分别为L 对
Figure BDA00029310157600001011
δr的偏导数，
Figure BDA00029310157600001012
分别为M 对
Figure BDA00029310157600001013
δr的偏导数，
Figure BDA00029310157600001014
分别为N 对
Figure BDA0002931015760000111
δr的偏导数。
第三步，针对倾转旋翼无人机的飞行特性，设计控制器和辨识方案。其中， 辨识方案包括试飞次数、辨识方法以及试飞输入信号的类型、幅值、频段和长 度等。
辨识方案具体如下：在控制空速和短舱倾角保持不变的情况下，分别进行 俯仰通道、滚转通道及偏航通道的扫频试验2～3次，待辨识通道称为主通道， 其余称为副通道；采用扫频信号作为主通道的输入信号，扫频幅值引起的飞机 角度响应范围为-a～a，a的取值范围为5°～15°，频率范围为0.5ωBW≤ω≤2.5ω180， ωBW表示带宽频率，ω180表示相角为-180°时所对应的频率；扫频试验在倾转旋 翼无人机处于配平状态3s～5s后开始，扫频结束后在配平状态停留至少3s后结 束一次扫频试验，每次扫频试验的时间
Figure BDA0002931015760000112
ωmin表示最小频率， ωmin＝0.5ωBW；在副通道输入幅值为主通道输入幅值20％的白噪声，记录控制器 的输出量和状态空间模型的输出量进行辨识。
由于倾转旋翼无人机的某些飞行模态发散，并且计算机生成的扫频输入会 导致飞行状态相对于配平状态发生漂移，因此，扫频试验必须在有控制器参与 的闭环状态下进行，本发明采用PID控制器进行三轴姿态控制，如图2所示。
本发明实施例1中，针对倾转旋翼无人机的直升机模式滚转通道和俯仰通 道的扫频试验在定高配平状态下进行，配平后等待3s，姿态角指令输入幅值为 15°的扫频信号15s，频段1rad/s～30rad/s，对副通道输入小幅值(幅值为主 通道输入幅值的20％)的白噪声以避免矩阵奇异，扫频试验结束后在配平状态 再等待3s结束一次扫频试验。分别针对滚转通道和俯仰通道进行3次试验， 采用直接辨识方法进行辨识，即直接测量控制器输出量和状态空间模型输出量 进行辨识。
第四步，根据辨识方案，在控制器控制下进行试飞试验，采集试验数据， 从中选取与状态空间模型的待辨识参数有关的时间历程数据，将选取的时间历 程数据通过线性调节Z变换计算频率响应，对频率响应进行处理，生成频率响 应数据。
具体地，在采集的试验数据中，截取扫频试验以及扫频试验前后各3s的 时间历程数据，并记录在pixhawk中，通过MATLAB读取记录数据。本发明 实施例1针对直升机模式的辨识主要关注滚转通道和俯仰通道的稳定导数和操 纵导数，因此，在确定待辨识参数之后，选取状态量
Figure BDA0002931015760000121
和控制量
Figure BDA0002931015760000122
则待辨识的状态空间模型如下：
Figure BDA0002931015760000123
将截取的时间历程数据通过线性调节Z变换计算频率响应，经过多输入规 整和复合窗口组合将频率响应转换为频率响应数据，剔除频率响应数据中频域响应相干值
Figure BDA0002931015760000124
或相干值曲线振荡的数据，剩余的频率响应数据用于辨识。
第五步，利用频率响应数据，通过优化算法辨识状态空间模型的参数；判 断状态空间模型的单输入单输出代价函数是否小于或等于第一阈值，且状态空 间模型的平均代价函数是否小于或等于第二阈值；若是，则执行第六步；若否， 则返回第二步，重复执行第二步～第五步。
具体地，利用频率响应数据，采用正割法结合状态空间模型对频域数据进 行参数拟合；
单输入单输出代价函数定义如下：
Figure BDA0002931015760000125
其中，||表示每个频率ω处的幅值，∠表示每个频率ω处的相位，nω表示 频率采样点的数量，一般取nω＝20，ω1表示拟合的起始频率，
Figure BDA0002931015760000126
表示拟合的终 止频率；
Figure BDA0002931015760000131
Wγ表示加权函数，当
Figure BDA0002931015760000132
时，Wγ将平方误差权 重降低至50％以下；Wg表示幅值平方误差的相对权重，Wp表示相位平方误差 的相对权重，
Figure BDA0002931015760000133
表示复合频率响应估计矩阵，T表示辨识得到的状态空间模型 的频率响应矩阵；本发明实施例1取Wg＝1，
Figure BDA0002931015760000134
即每1dB幅值误差和7.75° 的相位误差相当；
平均代价函数为：
Figure BDA0002931015760000135
Figure BDA0002931015760000136
当状态空间模型的单输入单输出代价函数满足J1≤200且平均代价函数满 足Jave≤100时，执行第六步；当状态空间模型的单输入单输出代价函数不满足 J1≤200或平均代价函数不满足Jave≤100时，返回第二步，重复执行第二步～第五 步。
本发明实施例1中，初次辨识得到的频率响应数据的代价值，如表1所示。
表1
Figure BDA0002931015760000137
从表1可以看出，所有频率响应数据的单输入单输出代价函数J1均满足 J1≤200，平均代价Jave满足Jave≤100，可以继续第六步。
第六步，计算待辨识参数的克拉美-劳边界和不敏感度，判断待辨识参数 的克拉美-劳边界是否小于或等于第三阈值，且待辨识参数的不敏感度是否小 于或等于第四阈值；若是，则结束操作，辨识出参数，得到状态空间模型；若 否，则执行第七步。
具体地，第i个辨识参数的克拉美-劳边界CRi由Hessian矩阵H的逆阵中 的对角线元素确定：
Figure BDA0002931015760000141
其中，np×np阶矩阵H定义为：
Figure BDA0002931015760000142
第i个辨识参数的克拉美-劳边界CRi由收敛辨识参数值的百分比形式表示 为：
Figure BDA0002931015760000143
第i个辨识参数的不敏感度Ii由Hessian矩阵H中的对角线元素确定：
Figure BDA0002931015760000144
第i个辨识参数的不敏感度Ii由收敛辨识参数值的化百分比形式表示为：
Figure BDA0002931015760000145
当第i个辨识参数的克拉美-劳边界满足
Figure BDA0002931015760000146
且不敏感度满足
Figure BDA0002931015760000147
时，说明状态空间模型辨识具有很高的可靠性，且具有良好的预测能力，此时 可以结束操作，辨识出参数，得到状态空间模型；当第i个辨识参数的克拉美- 劳边界不满足
Figure BDA0002931015760000148
或不敏感度满足
Figure BDA00029310157600001410
时，执行第七步。
需要说明的是，在不损失可靠性的前提下，可以允许一些最大的
Figure BDA0002931015760000149
达到20％～40％，最大的
Figure BDA0002931015760000151
达到12％～20％。
本发明实施例1中，计算克拉美-劳边界和不敏感度结果如表2所示。
表2
Figure BDA0002931015760000152
从表2可以看出，
Figure BDA0002931015760000153
和
Figure BDA0002931015760000154
的
Figure BDA0002931015760000155
远远超出10％，Lq、Mp、
Figure BDA0002931015760000156
和
Figure BDA0002931015760000157
的
Figure BDA0002931015760000158
远远超出20％，说明目前的状态空间模型存在多余的参数，模型结构需要进一 步优化，即执行第七步。
第七步，移除具有最大不敏感度的参数，返回第五步，重新进行辨识和判 断：
1)若状态空间模型的单输入单输出代价函数大于第一阈值，且状态空间模 型的平均代价函数大于第二阈值，即若不满足J1≤200或Jave≤100，则将移除的 参数重新插入状态空间模型，执行第八步；
2)若状态空间模型的单输入单输出代价函数小于或等于第一阈值，且状态 空间模型的平均代价函数小于或等于第二阈值，即若满足J1≤200且Jave≤100， 则进一步判断状态空间模型的单输入单输出代价函数的变化量是否大于或等 于第五阈值(第五阈值的取值范围为1～2)，且状态空间模型的平均代价函数的 变化量是否大于或等于第六阈值(第六阈值的取值范围为10～20)，即判断两个 代价函数是否发生跳跃性变化，即判断是否满足ΔJave＜1～2和ΔJl＜10～20；
(1)若不满足ΔJave＜1～2和ΔJl＜10～20，则将移除的参数重新插入状态空间模型，执行第八步；
(2)若满足ΔJave＜1～2和ΔJl＜10～20，则返回第六步，重新计算剩余参数的 克拉美-劳边界和不敏感度；若剩余参数的不敏感度小于或等于第四阈值，即 若剩余参数的不敏感度满足
Figure BDA0002931015760000165
，则执行第八步；若剩余参数的不敏感度大 于第四阈值，即若剩余参数的不敏感不满足
Figure BDA0002931015760000166
，则返回第七步，继续移除 具有最大不敏感度的参数。
本发明实施例1中，依次移除具有最大不敏感度的参数
Figure BDA0002931015760000161
和
Figure BDA0002931015760000162
将参数
Figure BDA0002931015760000163
和
Figure BDA0002931015760000164
固定并设为0，重新辨识，结果满足J1≤200和Jave≤100，ΔJave＜1～2 和ΔJl＜10～20，此时需要返回第六步，重新计算剩余参数的克拉美-劳边界和不 敏感度，结果剩余参数的不敏感度小于第四阈值，即剩余参数的不敏感度满足 Ii≤10％，此时执行第八步。
第八步，移除具有最大克拉美-劳边界的参数，返回第五步，重新进行辨 识和判断：
1)若状态空间模型的单输入单输出代价函数大于第一阈值，且状态空间模 型的平均代价函数大于第二阈值，即若不满足J1≤200或Jave≤100，则将移除的 参数重新插入状态空间模型，结束操作；
2)若状态空间模型的单输入单输出代价函数小于或等于第一阈值，且状态 空间模型的平均代价函数小于或等于第二阈值，即若满足J1≤200且Jave≤100， 则进一步判断状态空间模型的单输入单输出代价函数的变化量是否大于或等 于第五阈值(第五阈值的取值范围为1～2)，且状态空间模型的平均代价函数的 变化量是否大于或等于第六阈值(第六阈值的取值范围为10～20)，即判断两个 代价函数是否发生跳跃性变化，即判断是否满足ΔJave＜1～2和ΔJl＜10～20；
(1)若不满足ΔJave＜1～2和ΔJl＜10～20，则将移除的参数重新插入状态空间模型，结束操作；
(2)若满足ΔJave＜1～2和ΔJl＜10～20，则返回第六步，重新计算剩余参数的 克拉美-劳边界和不敏感度；若剩余参数的克拉美-劳边界小于或等于第三阈值， 即若剩余参数的克拉美-劳边界满足
Figure BDA0002931015760000171
则结束操作；若剩余参数的克拉 美-劳边界大于第三阈值，即若剩余参数的克拉美-劳边界不满足
Figure BDA0002931015760000172
则 返回第八步，继续移除具有最大克拉美-劳边界的参数。
本发明实施例1中，依次移除具有最大克拉美-劳边界的参数Lq和Mp，将 参数Lq和Mp固定并设为0，重新辨识，结果满足J1≤200且Jave≤100，但不满足 ΔJave＜1～2和ΔJl＜10～20，但由于克拉美-劳边界值能接受的最大值为40％，因 此，模型结构还是将参数Lq和Mp去除，辨识结果如下：
Figure BDA0002931015760000173
辨识结果曲线如图3和图4所示，图3中的(a)为滚转角速度对差动总距的 频率响应辨识曲线，图3中的(b)为滚转角对差动总距的频率响应辨识曲线，图 4中的(a)为俯仰角速度对纵向周期变距的频率响应辨识曲线，图4中的(b)为俯 仰角对纵向周期变距的频率响应辨识曲线，其中实线为仿真数据计算得到的频 率响应曲线，虚线为辨识拟合曲线，从曲线拟合结果来看，滚转通道和俯仰通 道的幅相曲线整体拟合程度较高。
频率响应数据的最终辨识代价表如表3所示，最终辨识参数分析表如表4 所示，均满足J1≤200且Jave≤100，
Figure BDA0002931015760000174
且
Figure BDA0002931015760000175
说明状态空间模型辨识 具有很高的可靠性，且具有良好的预测能力。
表3
Figure BDA0002931015760000176
Figure BDA0002931015760000181
表4
Figure BDA0002931015760000182
本发明提供的上述倾转旋翼无人机状态空间模型辨识方法，基于倾转旋翼 无人机试验数据，通过频域辨识，建立倾转旋翼无人机的状态空间模型，并通 过对模型参数精确度进行分析以对状态空间模型的结构进行优化，从而提高状 态空间模型的精确度和可靠性，增强状态空间模型的鲁棒性。为了使状态空间 模型的辨识结果能够更好地匹配频率响应数据，需要对辨识参数精确度或相对 置信度指标进行分析，以便改进状态空间模型的结构，通过模型优化使得辨识 得到的状态空间模型在满足输入输出最大拟合的基础上没有被过度参数化，从 而增加状态空间模型的精确度和可靠性。利用本发明提供的上述倾转旋翼无人 机状态空间模型辨识方法，可以得到倾转旋翼无人机不同飞行模式和飞行状态下的状态空间模型，从而应用于飞行控制系统设计，操纵品质评估和仿真验证。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发 明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及 其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.一种倾转旋翼无人机状态空间模型辨识方法，其特征在于，包括如下步骤：
S1：建立倾转旋翼无人机的运动微分方程；
S2：将所述运动微分方程转换为状态空间模型，并根据先验知识确定所述状态空间模型的参数初值；
S3：针对倾转旋翼无人机的飞行特性，设计控制器和辨识方案；
S4：根据所述辨识方案，在所述控制器控制下进行试飞试验，采集试验数据，从中选取与所述状态空间模型的待辨识参数有关的时间历程数据，将选取的时间历程数据通过线性调节Z变换计算频率响应，对所述频率响应进行处理，生成频率响应数据；
S5：利用所述频率响应数据，通过优化算法辨识所述状态空间模型的参数；判断所述状态空间模型的单输入单输出代价函数是否小于或等于第一阈值，且所述状态空间模型的平均代价函数是否小于或等于第二阈值；若是，则执行步骤S6；若否，则返回步骤S2，重复执行步骤S2～步骤S5；
S6：计算待辨识参数的克拉美-劳边界和不敏感度，判断待辨识参数的克拉美-劳边界是否小于或等于第三阈值，且待辨识参数的不敏感度是否小于或等于第四阈值；若是，则结束操作，辨识出参数，得到状态空间模型；若否，则执行步骤S7；
S7：移除具有最大不敏感度的参数，返回步骤S5，重新进行辨识和判断；若所述状态空间模型的单输入单输出代价函数大于第一阈值，且所述状态空间模型的平均代价函数大于第二阈值，则将移除的参数重新插入所述状态空间模型，执行步骤S8；若所述状态空间模型的单输入单输出代价函数小于或等于第一阈值，且所述状态空间模型的平均代价函数小于或等于第二阈值，则进一步判断所述状态空间模型的单输入单输出代价函数的变化量是否大于或等于第五阈值，且所述状态空间模型的平均代价函数的变化量是否大于或等于第六阈值，若是，则将移除的参数重新插入所述状态空间模型，执行步骤S8，若否，则返回步骤S6，重新计算剩余参数的克拉美-劳边界和不敏感度，若剩余参数的不敏感度小于或等于第四阈值，则执行步骤S8，若剩余参数的不敏感度大于第四阈值，则返回步骤S7，继续移除具有最大不敏感度的参数；
S8：移除具有最大克拉美-劳边界的参数，返回步骤S5，重新进行辨识和判断；若所述状态空间模型的单输入单输出代价函数大于第一阈值，且所述状态空间模型的平均代价函数大于第二阈值，则将移除的参数重新插入所述状态空间模型，结束操作；若所述状态空间模型的单输入单输出代价函数小于或等于第一阈值，且所述状态空间模型的平均代价函数小于或等于第二阈值，则进一步判断所述状态空间模型的单输入单输出代价函数的变化量是否大于或等于第五阈值，且所述状态空间模型的平均代价函数的变化量是否大于或等于第六阈值，若是，则将移除的参数重新插入所述状态空间模型，结束操作，若否，则返回步骤S6，重新计算剩余参数的克拉美-劳边界和不敏感度，若剩余参数的克拉美-劳边界小于或等于第三阈值，则结束操作，若剩余参数的克拉美-劳边界大于第三阈值，则返回步骤S8，继续移除具有最大克拉美-劳边界的参数。
2.如权利要求1所述的倾转旋翼无人机状态空间模型辨识方法，其特征在于，步骤S3中，辨识方案，具体包括：
在控制空速和短舱倾角保持不变的情况下，分别进行俯仰通道、滚转通道及偏航通道的扫频试验2～3次，待辨识通道称为主通道，其余称为副通道；采用扫频信号作为主通道的输入信号，扫频幅值引起的飞机角度响应范围为-a～a，a的取值范围为5°～15°，频率范围为0.5ωBW≤ω≤2.5ω180，ωBW表示带宽频率，ω180表示相角为-180°时所对应的频率；扫频试验在倾转旋翼无人机处于配平状态3s～5s后开始，扫频结束后在配平状态停留至少3s后结束一次扫频试验，每次扫频试验的时间
Figure FDA0002931015750000021
ωmin表示最小频率，ωmin＝0.5ωBW；在副通道输入幅值为主通道输入幅值20％的白噪声，记录控制器的输出量和所述状态空间模型的输出量进行辨识。
3.如权利要求2所述的倾转旋翼无人机状态空间模型辨识方法，其特征在于，步骤S4中，从中选取与所述状态空间模型的待辨识参数有关的时间历程数据，将选取的时间历程数据通过线性调节Z变换计算频率响应，对所述频率响应进行处理，生成频率响应数据，具体包括：
截取扫频试验以及扫频试验前后各3s的时间历程数据，将截取的时间历程数据通过线性调节Z变换计算频率响应，经过多输入规整和复合窗口组合将所述频率响应转换为频率响应数据，剔除所述频率响应数据中频域响应相干值
Figure FDA0002931015750000035
或相干值曲线振荡的数据，剩余的频率响应数据用于辨识。
4.如权利要求3所述的倾转旋翼无人机状态空间模型辨识方法，其特征在于，步骤S5，利用所述频率响应数据，通过优化算法辨识所述状态空间模型的参数；判断所述状态空间模型的单输入单输出代价函数是否小于或等于第一阈值，且所述状态空间模型的平均代价函数是否小于或等于第二阈值；若是，则执行步骤S6；若否，则返回步骤S2，重复执行步骤S2～步骤S5，具体包括：
利用所述频率响应数据，采用正割法结合所述状态空间模型对频域数据进行参数拟合；
单输入单输出代价函数定义如下：
Figure FDA0002931015750000031
其中，||表示每个频率ω处的幅值，∠表示每个频率ω处的相位，nω表示频率采样点的数量，ω1表示拟合的起始频率，
Figure FDA0002931015750000032
表示拟合的终止频率；
Figure FDA0002931015750000033
Wγ表示加权函数；Wg表示幅值平方误差的相对权重，Wp表示相位平方误差的相对权重，
Figure FDA0002931015750000034
表示复合频率响应估计矩阵，T表示辨识得到的状态空间模型的频率响应矩阵；
平均代价函数为：
Figure FDA0002931015750000041
Figure FDA0002931015750000042
若所述状态空间模型的单输入单输出代价函数满足J1≤200且平均代价函数满足Jave≤100，则执行步骤S6；否则，返回步骤S2，重复执行步骤S2～步骤S5。
5.如权利要求4所述的倾转旋翼无人机状态空间模型辨识方法，其特征在于，步骤S6，计算待辨识参数的克拉美-劳边界和不敏感度，判断待辨识参数的克拉美-劳边界是否小于或等于第三阈值，且待辨识参数的不敏感度是否小于或等于第四阈值；若是，则结束操作，辨识出参数，得到状态空间模型；若否，则执行步骤S7，具体包括：
第i个辨识参数的克拉美-劳边界CRi由Hessian矩阵H的逆阵中的对角线元素确定：
Figure FDA0002931015750000043
其中，np×np阶矩阵H定义为：
Figure FDA0002931015750000044
第i个辨识参数的克拉美-劳边界CRi由收敛辨识参数值的百分比形式表示为：
Figure FDA0002931015750000045
第i个辨识参数的不敏感度Ii由Hessian矩阵H中的对角线元素确定：
Figure FDA0002931015750000051
第i个辨识参数的不敏感度Ii由收敛辨识参数值的化百分比形式表示为：
Figure FDA0002931015750000052
若第i个辨识参数的克拉美-劳边界满足
Figure FDA0002931015750000053
且不敏感度满足
Figure FDA0002931015750000054
则结束操作，辨识出参数，得到状态空间模型；否则，执行步骤S7。
6.如权利要求5所述的倾转旋翼无人机状态空间模型辨识方法，其特征在于，步骤S7和步骤S8中，所述第五阈值的取值范围为1～2，所述第六阈值的取值范围为10～20。
CN202110148131.4A 2021-02-03 2021-02-03 一种倾转旋翼无人机状态空间模型辨识方法 Active CN113189866B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110148131.4A CN113189866B (zh) 2021-02-03 2021-02-03 一种倾转旋翼无人机状态空间模型辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110148131.4A CN113189866B (zh) 2021-02-03 2021-02-03 一种倾转旋翼无人机状态空间模型辨识方法

Publications (2)

Publication Number Publication Date
CN113189866A CN113189866A (zh) 2021-07-30
CN113189866B true CN113189866B (zh) 2022-04-26

Family

ID=76972846

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110148131.4A Active CN113189866B (zh) 2021-02-03 2021-02-03 一种倾转旋翼无人机状态空间模型辨识方法

Country Status (1)

Country Link
CN (1) CN113189866B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113901572B (zh) * 2021-12-09 2022-03-04 中航金城无人系统有限公司 基于频域辨识的多旋翼飞行器传递函数模型建模方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106507927B (zh) * 2005-05-20 2010-03-10 成都飞机工业（集团）有限责任公司 用时域气动等效确定飞机运动模态的方法
CN105843040A (zh) * 2016-04-06 2016-08-10 沈阳上博智拓科技有限公司 无人直升机动力学参数的辨识方法和装置
CN111538237A (zh) * 2020-03-20 2020-08-14 北京航空航天大学 一种倾转旋翼无人机非线性浅灰模型辨识与校正方法
CN112068582A (zh) * 2020-09-24 2020-12-11 北京航空航天大学 一种倾转旋翼无人机过渡模式模型辨识方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7623993B2 (en) * 2003-12-09 2009-11-24 Iowa State University Research Foundation, Inc. Method and system to perform energy-extraction based active noise control

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106507927B (zh) * 2005-05-20 2010-03-10 成都飞机工业（集团）有限责任公司 用时域气动等效确定飞机运动模态的方法
CN105843040A (zh) * 2016-04-06 2016-08-10 沈阳上博智拓科技有限公司 无人直升机动力学参数的辨识方法和装置
CN111538237A (zh) * 2020-03-20 2020-08-14 北京航空航天大学 一种倾转旋翼无人机非线性浅灰模型辨识与校正方法
CN112068582A (zh) * 2020-09-24 2020-12-11 北京航空航天大学 一种倾转旋翼无人机过渡模式模型辨识方法

Also Published As

Publication number Publication date
CN113189866A (zh) 2021-07-30

Similar Documents

Publication Publication Date Title
Wang et al. Quadrotor fault tolerant incremental sliding mode control driven by sliding mode disturbance observers
Wang et al. Flexible aircraft gust load alleviation with incremental nonlinear dynamic inversion
CN108717265B (zh) 一种基于控制变量参数化的无人飞行器巡航跟踪控制系统及控制方法
CN111538237B (zh) 一种倾转旋翼无人机非线性浅灰模型辨识与校正方法
CN111367182A (zh) 考虑输入受限的高超声速飞行器抗干扰反步控制方法
CN110794864A (zh) 基于姿态角速率与攻角测量的飞行器稳定控制方法
CN109164708B (zh) 一种高超声速飞行器神经网络自适应容错控制方法
CN113110511B (zh) 一种基于广义模糊双曲模型的智能船舶航向控制方法
CN111399530A (zh) 一种基于反解传函的小型飞行器攻角滑模跟踪方法
CN113189866B (zh) 一种倾转旋翼无人机状态空间模型辨识方法
CN107831653B (zh) 一种抑制参数摄动的高超声速飞行器指令跟踪控制方法
CN111695193A (zh) 一种全局相关三维气动力数学模型的建模方法及系统
CN110929339B (zh) 高精度亚音速固定翼飞行器非线性动力学建模方法
CN113128035A (zh) 民用飞机飞控传感器信号重构容错控制方法
Suresh et al. An on-line learning neural controller for helicopters performing highly nonlinear maneuvers
Chen Full state constrained stochastic adaptive integrated guidance and control for STT missiles with non-affine aerodynamic characteristics
CN111427267A (zh) 一种采用力与力矩自适应估计的高速飞行器攻角跟踪方法
CN112068582B (zh) 一种倾转旋翼无人机过渡模式模型辨识方法
CN111176117B (zh) 一种无人直升机的模糊自适应弹性控制方法
Grauer et al. Aeroelastic modeling of X-56A stiff-wing configuration flight test data
CN115793449B (zh) 直升机飞行姿态控制器设计方法及装置、设备、存储介质
Seher-Weiss et al. EC135 system identification for model following control and turbulence modeling
CN116360255A (zh) 一种非线性参数化高超声速飞行器的自适应调节控制方法
CN106651090A (zh) 一种归一化的人机系统飞行品质预测方法
CN117519257B (zh) 一种基于反步法的超声速巡航高度控制方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant