-
摘要: 本文研究了一类具有边界执行器动态特性的双曲线型偏微分方程(Partial differential equation, PDE)系统的输出调节问题. 特别地, 执行器由一组非线性常微分方程(Ordinary differential equation, ODE)描述, 控制输入出现在执行器的一端而非直接作用在PDE系统上, 这使得控制任务变得相当困难. 基于几何设计方法和有限维与无限维反步法, 本文提出了显式表达的输出调节器, 实现了该类系统的扰动补偿及跟踪控制. 并且我们采用Lyapunov稳定性理论严格证明了闭环系统及跟踪误差在范数意义上的指数稳定性. 仿真实例对比验证了所提出控制方法的有效性.Abstract: This paper investigates the output regulation problem for a class of hyperbolic partial differential equation (PDE) systems with boundary actuator dynamics. Particularly, the control input appears at one end of the actuator described by a set of ordinary differential equation (ODE) rather than directly in the PDE system, which makes the control task rather difficult. Based on the geometric design method as well as finite and infinite dimensional backstepping methods, an output regulator is explicitly provided in the paper so that the disturbance compensation and tracking control of this system are implemented. Moreover, we rigorously prove the exponential stability of both the closed-loop system and the tracking error in the norm by employing the Lyapunov stability theory. The simulation example comparatively demonstrates the effectiveness of the proposed control method.
-
近年来, 仿生机器人在机器人界得到越来越多的关注[1-4], 而扑翼飞行正是仿生机器人研究中的焦点之一[5-7]. 对扑翼飞行原理的探索不仅具有深刻而广泛的理论价值, 更具有较高的工程应用价值. 与此同时, 由于扑翼飞行的高效性、灵活性和隐蔽性, 其在地形探索、目标追踪、军事侦察、生态监测等领域均有广阔应用前景.
由于扑翼飞行的形态学特性, 弹性结构以及其特殊的往复运动模式, 在扑动过程中, 气体的粘性力和惯性力的作用均不可忽视, 进而造成了多种复杂的非定常气动力学现象[1, 8-10]. 具体而言, Ellington等在文献[11-12]中详细分析了由于前缘涡的存在, 产生失速延时(可以不失速地大攻角平动)的现象, 阐明了其在提升扑翼升力中所发挥的重要作用, 并在实际实验中对前缘涡的具体产生和演变过程进行了观测与分析. Dickinson等在文献[13-14]中强调了翅翼旋转产生的升力对昆虫飞行的重要作用, 并揭示了其在昆虫机动飞行中所扮演的重要角色. Weis-Fogh[15]在研究丽蚜小蜂的扑翼飞行时发现了著名的拍合−剥离特性(Clap-and-fling), 他认为在拍合过程中并没有明显的气动力学影响, 但是翅翼拍合保证了翅翼在剥离之前位置的确定性, 并且指出在剥离过程中在左右翼之间出现绕翼环流, 在完成剥离之后环流也随即迅速演化成翅翼边缘的稳定涡流, 整个过程均提升了翅翼所产生的升力. Miller等[16]和Percin等[17]分别研究了翅翼弹性在拍合−剥离过程中的具体作用与影响, 并指出拍合过程中向下的射流也有助于提升升力. Shyy等和Chin等分别在文献[18]和文献[19]中, 对扑翼的非定常气动力特性进行了系统总结, 具体包括前缘涡效应、翼面快速旋转、拍合−剥离、附加质量效应、尾迹捕捉等.
建立扑翼飞行的动力学模型需要对复杂的气动力特性进行必要的简化. 而求解不可压粘性流体运动对应的纳维−斯托克斯方程非常困难, 建立扑翼空气动力学的准稳态模型, 几乎是解析分析扑翼系统动力学特性的唯一方式. 根据不同的具体系统特性, 以及模型具体的应用场景, 可以建立不同的准稳态模型[14, 19-23]. 所建立的准稳态扑翼动力学模型, 既可以用于扑翼动力学系统特性的分析, 也可以用于实际扑翼飞行器的设计与控制. Sun等在文献[24-26]中详细分析了昆虫躯干−翅翼的多刚体模型后, 用“刚体假设”简化了翅翼的运动从而获得机体动力学模型, 并进一步用扑翼周期平均的方法分析了系统的稳定性和能控性. 他们还将系统在周期解附近线性化, 利用Floquet原理分析了系统周期解附近小扰动干扰下的稳定性. Taha等在文献[27]中强调了高阶平均化理论在分析扑翼稳定性中的重要性, 但是在系统简化过程中他们直接用三角函数模拟了翅翼的运动学行为, 忽略了翅翼周期运动的动力学特性. Cheng等在文献[28]中考虑了机体运动对翅翼运动学的影响, 并利用气动准稳态模型进一步建立了三维空间6自由度刚体的模型, 且分析了4种昆虫的被动稳定性. 在实际的扑翼飞行器控制中, 作用于机体的力与力矩通常和翅翼的运动学关联, 或者更直接地与决定翅翼运动的参数相关联. 通过把翅翼步态映射到周期平均力与力矩, 并以6自由度刚体运动模型作为系统模型, 从而设计相应的控制器来完成针对不同环境 或任务的机体姿态、位置或速度的控制方法[29-32]. 其中需要特别指出, Ramezani等在文献[33-34]中采用拉格朗日方法对仿蝙蝠扑翼飞行器进行多刚体建模, 并指出扑翼系统的零动态难以确定, 而采用反馈线性化的方法来设计控制器, 并且分析了闭环系统位置环和姿态环的时标分离.
在这些扑翼系统建模和控制的方法中, 极少考虑翅翼的动力学特性. 这也就默认了翅翼系统周期运动的稳定性, 且其具体运动特性不会被机体本身运动所干扰. 但这种稳定性并不是系统固有的, 所以这种单纯从运动学角度考虑翅翼行为的方式, 显然是对系统的过度简化, 难以全面描述系统的具体动力学特性. 因此, 在扑翼系统动力学系统中, 对翅翼动力学系统的周期稳定性分析是非常必要的. 此外, 扑翼动力学系统中往往包含着子系统之间的时标分离. 并且, Orlowski等在文献[35]中指出, 翅翼和机体之间的运动可能不处于同一时间尺度, 但由于翅翼系统没有确定的收敛点而可能只存在某个收敛轨迹, 所以无法使用常规的针对临界点的奇异摄动理论进行分析.
对于上述存在的问题, 本文针对扑翼飞行周期性系统进行奇异摄动分析. 具体而言, 我们在给定周期进行频闪采样, 利用这些采样构建相应的离散系统, 进而反映扑翼飞行的状态变化. 通过周期输入和简化模型共同确定扑翼周期的步态, 并在一定程度上忽略实际周期内的具体行为. 进而通过观测构建的离散系统, 综合周期输入和简化模型确定的步态信息, 共同估计扑翼飞行的周期状态, 并在这一基础上使用奇异摄动理论对其稳定性进行了分析. 本文主要贡献可以概括为以下两个方面:
1)建立尽可能简洁的扑翼多刚体动力学模型, 为扑翼飞行周期动力学的奇异摄动分析提供了基础;
2)基于奇异摄动理论, 提出了适用于实际扑翼飞行问题的系统周期稳定性分析方法.
本文其余部分组织如下: 在第1节中建立了扑翼多刚体模型, 在第2节中利用奇异摄动理论分析了扑翼飞行周期性系统的稳定性, 在第3节中利用实际实验验证了所提方法的有效性和可行性, 最后在第4节做出了总结.
1. 扑翼飞行的多刚体模型
尽管扑翼飞行的种类多种多样, 其系统在某种程度上依旧可以视为在不可压低粘性流体中运动的多刚体动力学模型. 在准定常分析气动力的前提下, 进而可以方便地利用系统运动状态计算出系统所受广义气动力, 从而可以使用经典的 拉格朗日动力学建模方法进行处理[34, 36-37]. 为了加以区别, 如图1所示, 我们将除躯干或者机体之外的刚体称为连杆, 进而将连接占主要质量和转动惯量的躯干的连杆称为初级连杆, 将连接初级连杆的连杆称为次级连杆, 以此类推.
图 1 不同种类扑翼飞行多刚体模型中的坐标系建立Fig. 1 The coordinate setup in the multi-rigid body dynamics of different flapping wing flights通常而言, 我们可以定义如下的广义坐标
$$ {\boldsymbol{q}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{q}}_p^{\rm{T}}} & {{\boldsymbol{q}}_a^{\rm{T}}} & {{\boldsymbol{q}}_A^{\rm{T}}} & {{\boldsymbol{q}}_B^{\rm{T}}} \\ \end{array}} \right]^{\rm{T}}} $$ (1) 其中,
$ {\boldsymbol{q}}_p\in {\bf{R}}^3 $ 为占主要质量和转动惯量的躯干或者机体的三维位置, 而$ {\boldsymbol{q}}_a\in {\cal{S}}^3 $ 为上述躯干或者机体的三维姿态,$ {\boldsymbol{q}}_A $ 为连接上述躯干或者机体与初级连杆对应关节的广义关节坐标, 而$ {\boldsymbol{q}}_B $ 为连接初级与次级连杆对应关节被动的广义关节坐标, 其具体维数由具体飞行器结构决定. 整体动力学可以由拉格朗日方程所描述$$ \frac{\rm{d}}{{{\rm{d}}t}}\left( {\frac{{\partial K}}{{\partial {{\dot { q}}}}}} \right) - \frac{{\partial K}}{{\partial {q}}} + \frac{{\partial P}}{{\partial {q}}} = {G} $$ (2) 其中,
$ K = K\left( {{\boldsymbol{q}}, \dot {{\boldsymbol{q}}}} \right) $ 为系统总功能,$ P = P\left( {\boldsymbol{q}} \right) $ 为系统总势能. 具体来看,$ K$ 包括多刚体系统各组分的转动和平动动能,$ P $ 包括多刚体系统各组分的重力势能和弹性势能(很多扑翼生物或飞行器会周期性地存储和释放弹性势能), 而$ G $ 则主要指主动关节的驱动力或者力矩和多刚体系统各组分所受的空气动力或者力矩(包括惯性力和粘性力). 由于运动的拓扑结构限制, 对于大多数扑翼机体, 某个特定翼面通常以一个特定关节或以一组集中的关节组与 占主要质量和转动惯量的躯干相连接. 为了统一, 这些关节或者关节组均被认为是具有转动幅度或者方向限制的三个自由度旋转轴. 在图2中, 我们具体描述了惯性坐标系$ {\cal{F}}_I $ 、体坐标系$ {\cal{F}}_B $ 、各级连杆刚体中心坐标系以及转轴坐标之间的关系. 请注意, 同一刚体的转轴坐标和刚体质心坐标之间只存在平移变换$ {\boldsymbol{D}}_i $ 和$ {\boldsymbol{D}}_j $ , 而无旋转变换.
图 2 翅翼两级连杆系统中各刚体坐标系之间的旋转平移示意图Fig. 2 The schematic for the rotations and the translations between different frames of a two-link wing以下我们以两层连杆构成的动力学系统为例进行分析, 此系统的动能可以公式化地规范为
$$ \begin{split} K\left( {{\boldsymbol{q}}, \dot {{\boldsymbol{q}}}} \right) =\;& \frac{1}{2}{m_1}\dot {{\boldsymbol{q}}}_p^{\rm{T}}{{\dot {{\boldsymbol{q}}}}_p} + \frac{1}{2}{{\boldsymbol{\omega}}_1 ^{\rm{T}}}({{\boldsymbol{q}}}_a, {{{\dot {{\boldsymbol{q}}}}_a}}){I_1}{\boldsymbol{\omega}}_1 ({{\boldsymbol{q}}}_a, {{{\dot {{\boldsymbol{q}}}}_a}} ) +\\ & \frac{1}{2}\sum\limits_{i \in A} {{m_i}{\boldsymbol{v}}_i^{\rm{T}}( {Q_j}\backslash \left\{ {{{\boldsymbol{q}}_j}, {{\dot {{\boldsymbol{q}}}}_j}} \right\} ){{\boldsymbol{v}}_i}( {Q_j}\backslash \left\{ {{{\boldsymbol{q}}_j}, {{\dot {{\boldsymbol{q}}}}_j}} \right\})}\;+\\ & \frac{1}{2}\sum\limits_{i \in A} {{\boldsymbol{\omega}} _i^{\rm{T}}({{\boldsymbol{q}}}_a, {{{\dot {{\boldsymbol{q}}}}_a}, {{\boldsymbol{q}}}_i, {{\dot {{\boldsymbol{q}}}}_i}} ){I_i}{{\boldsymbol{\omega}} _i}\left( {{\boldsymbol{q}}}_a, {{{\dot {{\boldsymbol{q}}}}_a}, {{\boldsymbol{q}}}_i, {{\dot {{\boldsymbol{q}}}}_i}} \right)}\;+ \\ & \frac{1}{2}\sum\limits_{j \in B} {{m_j}{\boldsymbol{v}}_j^{\rm{T}}\left( {{Q_j}} \right){{\boldsymbol{v}}_j}\left( {{Q_j}} \right)}\;+ \\ & \frac{1}{2}\sum\limits_{j \in B} {{\boldsymbol{\omega}} _j^{\rm{T}}\left( {Q_j}\backslash {{{\dot {{\boldsymbol{q}}}}_p}} \right){I_j}{{\boldsymbol{\omega}} _j}\left( {Q_j}\backslash {{{\dot {{\boldsymbol{q}}}}_p}} \right)} \\[-20pt]\end{split} $$ (3) 其中,
$ {Q_j} = \left\{ {{\dot {{\boldsymbol{q}}}_p, {{\boldsymbol{q}}}_a, {\dot {{\boldsymbol{q}}}}_a}, {{\boldsymbol{q}}_{i\left( j \right)}}, {{\dot {{\boldsymbol{q}}}}_{i\left( j \right)}}, {{\boldsymbol{q}}_j}, {{\dot {{\boldsymbol{q}}}}_j}} \right\}, $ ${\boldsymbol{q}}_p\in {\cal{S}}^3 ,$ 而${\boldsymbol{q}}_a\in {\cal{S}}^3 ,$ ${\boldsymbol{q}}_i\in {\cal{S}}^3 ,$ 以及${\boldsymbol{q}}_i\in {\cal{S}}^3 ,$ 则表示基于四元数表示的姿态或姿态的旋转,$ m \in {\bf{R}} $ 和$ I \in {\bf{R}}^{3\times3} $ 分别表示刚体的质量和基于其固连坐标系的转动惯量,$ i(j) $ 表示连杆的连接关系,$ A $ 表示初级连杆的集合, 而$ B $ 表示次级连杆的集合. 对应地, 本文中$ i $ 在不特别指定时表示初级连杆的标号, 而$ j $ 表示次级连杆的标号. 上述系统中的线速度和角速度可显式地表示为$$ \begin{split} &{{\boldsymbol{v}}_j} = \frac{\rm{d}}{{{\rm{d}}t}}\left( {{\boldsymbol{q}}_p} + T\left( {{{\boldsymbol{q}}_a}, {{\boldsymbol{L}}_a}} \right)T\left( {{{\boldsymbol{q}}_{i\left( j \right)}, {{\boldsymbol{L}}_{i(j)}}}} \right)T\left( {{{\boldsymbol{q}}_j}, {{\boldsymbol{L}}_j}} \right){{\boldsymbol{D}}_j} \right) \\ &{{\boldsymbol{\omega}} _j} = {{\boldsymbol{\omega}} _{{{\boldsymbol{q}}} ,\dot {{\boldsymbol{q}}}}}\left( {{{\boldsymbol{q}}_j} \otimes {{\boldsymbol{q}}_{i(j)}} \otimes {{\boldsymbol{q}}_a}, \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {{{\boldsymbol{q}}_j} \otimes {{\boldsymbol{q}}_{i(j)}} \otimes {{\boldsymbol{q}}_a}} \right)} \right) \\ &{{\boldsymbol{v}}_i} = \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {{{\boldsymbol{q}}_p} + T\left( {{{\boldsymbol{q}}_a} ,{{\boldsymbol{L}}_a}} \right)T\left( {{{\boldsymbol{q}}_i}, {{\boldsymbol{L}}_i}} \right)} {{\boldsymbol{D}}_i}\right) \\ &{{\boldsymbol{\omega}} _i} = {{\boldsymbol{\omega}} _{{\boldsymbol{q}}, \dot {{\boldsymbol{q}}}}}\left( {{{\boldsymbol{q}}_i} \otimes {{\boldsymbol{q}}_a}, \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {{{\boldsymbol{q}}_i} \otimes {{\boldsymbol{q}}_a}} \right)} \right) \\ &{{\boldsymbol{\omega}} _1} = {{\boldsymbol{\omega}} _{{\boldsymbol{q}} ,\dot {{\boldsymbol{q}}}}}\left( { {{\boldsymbol{q}}_a} ,{{\dot{{\boldsymbol{q}}}_a}} } \right)\\[-10pt] \end{split} $$ (4) 其中,
$ T({\boldsymbol{q}}, {\boldsymbol{L}})\in {\bf{R}}^{4\times4} $ 表示由旋转$ {\boldsymbol{q}} \in {\cal{S}}^3 $ 和关节参数$ {\boldsymbol{L}} \in {\bf{R}}^3 $ 决定的齐次矩阵,$ {\boldsymbol{D}}_i $ 与$ {\boldsymbol{D}}_j $ 均表示旋转中心到对应连杆质心的偏移量. 对于四元数${\boldsymbol{q}} = $ $ {[ \eta {{{\boldsymbol{\epsilon}} ^{\rm{T}}}} ]^{\rm{T}}}\in {\cal{S}}^3$ , 其中$ \eta $ 表示四元数标量部分, 而$ {\boldsymbol{\epsilon}} $ 表示四元数矢量部分. 操作符$ \otimes $ 具体含义表示为$$ \begin{split}{{\boldsymbol{q}}_1}\otimes{{\boldsymbol{q}}_2} = \;&\left[ {\begin{array}{*{20}{c}} {{\eta _1}{\eta _2} - {\boldsymbol{\epsilon}} _1^{\rm{T}}{{\boldsymbol{\epsilon}} _2}} \\ {{\eta _1}{{\boldsymbol{\epsilon}} _2} + {\eta _2}{{\boldsymbol{\epsilon}} _1}+S\left( {{{\boldsymbol{\epsilon}} _1}} \right){{\boldsymbol{\epsilon}} _2}} \\ \end{array}} \right]= \\ &M_L({\boldsymbol{q}}_1){\boldsymbol{q}}_2 = M_R({\boldsymbol{q}}_2){\boldsymbol{q}}_1 \end{split} $$ (5) 其中,
$ M_L $ 与$ M_D $ 表示利用矩阵和向量乘法代替四元数乘法的相应矩阵. 且考虑到四元数表述姿态的特性, 可知$$ {{\boldsymbol{\omega}} _{{\boldsymbol{q}}, \dot {{\boldsymbol{q}}}}\left( {\boldsymbol{q}}, \dot {{\boldsymbol{q}}}\right)}{\text{ = }}2{{{\boldsymbol{q}}}^*} \otimes \dot {{\boldsymbol{q}}}{\rm{ = }}2M_L^{ - 1}\left( {\boldsymbol{q}} \right)\dot {{\boldsymbol{q}}} $$ (6) 这里角速度和惯量矩阵均需扩展为四维, 并且 其中反对称矩阵
$ S\left( {\boldsymbol{x}} \right) $ 表示为$$ S\left( {\boldsymbol{x}} \right) = \left[ {\begin{array}{*{20}{c}} 0 & { - {x_3}} & {{x_2}} \\ {{x_3}} & 0 & {{x_1}} \\ { - {x_2}} & { - {x_1}} & 0 \end{array}} \right] $$ (7) 与动能相似, 同样我们也可以规范化地表示出系统的势能:
$$ \begin{split} P\left({\boldsymbol{q}} \right) = \;& {P_g}\left( {{{\boldsymbol{q}}_p}} \right) +\\ & \sum\limits_{i \in A} {{P_{gi}}\left( {{{\boldsymbol{q}}_p}, {{\boldsymbol{q}}_a}, {{\boldsymbol{q}}_i}} \right)} \; + \sum\limits_{i \in A} {{P_{ei}}\left( {{{\boldsymbol{q}}_i}} \right)}\; + \\ & \sum\limits_{j \in B} {{P_{gj}}\left( {{{\boldsymbol{q}}_p}, {{\boldsymbol{q}}_a}, {{\boldsymbol{q}}_{i\left( j \right)}}, {{\boldsymbol{q}}_j}} \right)} + \sum\limits_{j \in B} {{P_{ej}}\left( {{{\boldsymbol{q}}_j}} \right)} \end{split} $$ (8) 其中, Pg表示重力势能, 而Pe表示弹性势能. 另外, 为了简化问题, 我们忽略扑翼运动对周围空气流场的影响, 并将各刚体所受气动广义力分别独立地在惯性坐标系下表出:
$$\begin{split} &{{\boldsymbol{F}}_p} = {{\boldsymbol{F}}_{{\rm{aero}}, p}}\left( {{{\dot {{\boldsymbol{q}}}}_p}, {{\boldsymbol{q}}_a}, {{\dot {{\boldsymbol{q}}}}_a}} \right) \\ &{{\boldsymbol{F}}_a} = {{\boldsymbol{F}}_{{\rm{aero}}, a}}\left( {{{\dot {{\boldsymbol{q}}}}_p}, {{\boldsymbol{q}}_a}, {{\dot {{\boldsymbol{q}}}}_a}} \right) \\ &{{\boldsymbol{F}}_{pi}} = {{\boldsymbol{F}}_{{\rm{aero}}, pi}}\left( {{{\dot {{\boldsymbol{q}}}}_p}, {{\boldsymbol{q}}_a}, {{\dot {{\boldsymbol{q}}}}_a}, {{\boldsymbol{q}}_i}, {{\dot {{\boldsymbol{q}}}}_i}} \right) \\ &{{\boldsymbol{F}}_{ai}} = {{\boldsymbol{F}}_{{\rm{aero}}, ai}}\left( {{{\dot {{\boldsymbol{q}}}}_p}, {{\boldsymbol{q}}_a}, {{\dot {{\boldsymbol{q}}}}_a}, {{\boldsymbol{q}}_i}, {{\dot {{\boldsymbol{q}}}}_i}} \right) \\ &{{\boldsymbol{F}}_{pj}} = {{\boldsymbol{F}}_{{\rm{aero}}, pj}}\left( {Q_j} \right) \\ &{{\boldsymbol{F}}_{aj}} = {{\boldsymbol{F}}_{{\rm{aero}}, aj}}\left( {Q_j} \right) \end{split} $$ (9) 其中,
${\boldsymbol{F}}_{{\rm{aero}}, p}, \,{\boldsymbol{F}}_{{\rm{aero}}, pi},\, {\boldsymbol{F}}_{{\rm{aero}}, pj}\in {\bf{R}} ^3,$ ${\boldsymbol{F}}_{{\rm{aero}}, a}, {\boldsymbol{F}}_{{\rm{aero}}, ai}, $ $ {\boldsymbol{F}}_{{\rm{aero}}, aj}\in {\bf{R}} ^4$ 分别表示准稳态意义下气动广义力中气动力和气动力矩对应的项. 同时, 为简洁起见, 这里并不列出气动力具体形式. 进一步, 根据虚功原理可以得到广义力为$$ \begin{split} &{\boldsymbol{G}} = J\left({\left[ {{\boldsymbol{F}}_p^{\rm{T}}, {\boldsymbol{F}}_a^{\rm{T}}, {\boldsymbol{F}}_{pA}^{\rm{T}}, {\boldsymbol{F}}_{aA}^{\rm{T}}, {\boldsymbol{F}}_{pB}^{\rm{T}}, {\boldsymbol{F}}_{aB}^{\rm{T}}} \right]^{\rm{T}}}+ {{\boldsymbol{F}}_{u}}\left( t \right)\right) \\ &{{\boldsymbol{F}}_{u}}\left( t \right) = {\left[ {{{\bf{0}}^{\rm{T}}}, {{\bf{0}}^{\rm{T}}}, {{\bf{0}}^{\rm{T}}}, {\boldsymbol{F}}_{uA}^{\rm{T}}\left( t \right), {{\bf{0}}^{\rm{T}}}, {\boldsymbol{F}}_{uB}^{\rm{T}}\left( t \right)} \right]^{\rm{T}}}\\[-10pt] \end{split} $$ (10) 其中, 雅克比矩阵
$J = \dfrac{\partial {\boldsymbol{r}}} {\partial {\boldsymbol{q}}}$ ,$ {\boldsymbol{r}}({\boldsymbol{q}}) $ 表示各刚体位置姿态在惯性坐标系下表出与广义关节坐标的相对关系函数.$ {\boldsymbol{F}}_{ui}\in {\bf{R}} ^ 3, i \in \left\{A, B\right\} $ 表示关节主动输入力距. 从控制角度而言$ {{\boldsymbol{F}}_{u}} $ 可以解释为系统输入. 根据式(2) ~ (9), 推导可得次级关节动力学$$ \begin{split} &\frac{\rm{d}}{{{\rm{d}}t}}\left( {{m_j}{\left({\frac{{\partial {{\boldsymbol{v}}_j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_j}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_j} + {\left({\frac{{\partial {{\boldsymbol{\omega}} _j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_j}}}}\right)^{\rm{T}}}{I_j}{{\boldsymbol{\omega}} _j}} \right) - {m_j}{\left(\frac{{\partial {{\boldsymbol{v}}_j}}}{{\partial {{\boldsymbol{q}}_j}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_j} - \\ &\qquad {\left(\frac{{\partial {{\boldsymbol{\omega}} _j}}}{{\partial {{\boldsymbol{q}}_j}}}\right)^{\rm{T}}}{I_j}{{\boldsymbol{\omega}} _j} + {\left(\frac{{\partial {P_{gj}}}}{{\partial {{\boldsymbol{q}}_j}}}\right)^{\rm{T}}} + {\left(\frac{{\partial {P_{ej}}}}{{\partial {{\boldsymbol{q}}_j}}}\right)^{\rm{T}}} = {{\boldsymbol{G}}_{j}}\\[-15pt] \end{split} $$ (11) 类似地, 初级关节的动力学表示为
$$ \begin{split} &\frac{\rm{d}}{{{\rm{d}}t}}\left( {{m_i}{\left({\frac{{\partial {{\boldsymbol{v}}_i}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_i}{ + }{\left({\frac{{\partial {{\boldsymbol{\omega}} _i}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}}\right)^{\rm{T}}}{I_i}{{\boldsymbol{\omega}} _i} + {\boldsymbol{Sd}}{_i}} \right)- \\ & \quad\left( {{m_i}{\left({\frac{{\partial {{\boldsymbol{v}}_i}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_i}{ + }{\left({\frac{{\partial {{\boldsymbol{\omega}} _i}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}}}{I_i}{{\boldsymbol{\omega}} _i} + {{\boldsymbol{S}}_i}} \right)+ \left(\frac{{\partial {P_{gi}}}}{{\partial {{\boldsymbol{q}}_i}}}\right)^{\rm{T}}+ \\ & \quad\sum\limits_{j \in B, i\left( j \right) = i} \left({\frac{{\partial {P_{gj}}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}} + \left(\frac{{\partial {P_{ei}}}}{{\partial {{\boldsymbol{q}}_i}}}\right)^{\rm{T}} = {{\boldsymbol{G}}_{i}} \\ &{\boldsymbol{Sd}}{_i} = \sum\limits_{j \in B, i\left( j \right) = i} {\left( {{m_j}{\left({\frac{{\partial {{\boldsymbol{v}}_j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_j} + {\left({\frac{{\partial {{\boldsymbol{\omega}} _j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}}\right)^{\rm{T}}}{I_j}{{\boldsymbol{\omega}} _j}} \right)} \\ &{{\boldsymbol{S}}_i} = \sum\limits_{j \in B, i\left( j \right) = i} {\left( {{m_j}{\left({\frac{{\partial {{\boldsymbol{v}}_j}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_j} + {\left({\frac{{\partial {{\boldsymbol{\omega}} _j}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}}}{I_j}{{\boldsymbol{\omega}} _j}} \right)} \end{split} $$ (12) 对于占主要质量和转动惯量的躯干的关节有
$$ \begin{split} & \frac{\rm{d}}{{{\rm{d}}t}}\left( {{m_1}{{\dot {{\boldsymbol{q}}}}_p} + \sum\limits_{i \in A} {{m_i}{\left({\frac{{\partial {{\boldsymbol{v}}_i}}}{{\partial {{\dot {{\boldsymbol{q}}}}_p}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_i}} + \sum\limits_{j \in B} {{m_j}{\left({\frac{{\partial {{\boldsymbol{v}}_j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_p}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_j}} } \right)+ \\ & \qquad\left(\frac{{\partial {P_g}}}{{\partial {{\boldsymbol{q}}_p}}}\right)^{\rm{T}} + \sum\limits_{i \in A} \left({\frac{{\partial {P_{gi}}}}{{\partial {{\boldsymbol{q}}_p}}}}\right) ^{\rm{T}} + \sum\limits_{j \in B} \left({\frac{{\partial {P_{gj}}}}{{\partial {{\boldsymbol{q}}_p}}}}\right)^{\rm{T}} = {{\boldsymbol{G}}_{p}} \end{split} $$ (13) $$ \begin{split} &\frac{{\rm{d}}}{{{\rm{d}}t}}\left( {4{{\left( {M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right)} \right)}^{\rm{T}}}{I_1}M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right){{\dot {{\boldsymbol{q}}}}_a} + {\boldsymbol{S}}{{\boldsymbol{d}}_1}} \right)- \\ & \qquad\left( {4{{\left( {\frac{{\partial M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right){{\dot {{\boldsymbol{q}}}}_a}}}{{\partial {{\boldsymbol{q}}_a}}}} \right)}^{\rm{T}}}{I_1}M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right){{\dot {{\boldsymbol{q}}}}_a} + {{\boldsymbol{S}}_1}} \right)+ \\ &\qquad \sum\limits_{i \in A} {{\left({\frac{{\partial {P_{gi}}}}{{\partial {{\boldsymbol{q}}_a}}}}\right)^{\rm{T}}}} + \sum\limits_{j \in B} {{\left({\frac{{\partial {P_{gj}}}}{{\partial {{\boldsymbol{q}}_a}}}}\right)^{\rm{T}}}} = {{\boldsymbol{G}}_{a}} \\ &{\boldsymbol{S}}{{\boldsymbol{d}}_1} = \sum\limits_{i \in A, B} {\left( {{m_i}{\left({\frac{{\partial {{\boldsymbol{v}}_i}}}{{\partial {{\dot {{\boldsymbol{q}}}}_a}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_i} + {\left({\frac{{\partial {{\boldsymbol{\omega}} _i}}}{{\partial {{\dot {{\boldsymbol{q}}}}_a}}}}\right)^{\rm{T}}}{I_i}{{\boldsymbol{\omega}} _i}} \right)} \\ &{{\boldsymbol{S}}_1} = \sum\limits_{i \in A, B} {\left( {{m_i}{\left({\frac{{\partial {{\boldsymbol{v}}_i}}}{{\partial {{\boldsymbol{q}}_a}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_i} + {\left({\frac{{\partial {{\boldsymbol{\omega}} _i}}}{{\partial {{\boldsymbol{q}}_a}}}}\right)^{\rm{T}}}{I_i}{{\boldsymbol{\omega}} _i}} \right)}\\[-15pt] \end{split} $$ (14) 至此, 我们便基本建立了扑翼飞行的多刚体动力学模型.
2. 扑翼多刚体动力学的奇异摄动分析
本节主要研究由翅翼与躯干质量差异导致的奇异摄动现象. 首先考虑展开后的次级关节动力学
$$ \begin{split} & \frac{{\rm{d}}} {{{\rm{d}}t}}\Bigg\{ {{m_j}\left( {T_a}{T_{i\left( j \right)}}\frac{{\partial {{\dot T}_j}{D_j}}} {{\partial {{\dot {{\boldsymbol{q}}}}_j}}}\right)^{\rm{T}} {{{\boldsymbol{v}}}_j}}+\\ & \quad {2{{\left[ {{{\left( {{{{\boldsymbol{q}}}_j} \otimes {{{\boldsymbol{q}}}_{i\left( j \right)}} \otimes {{{\boldsymbol{q}}}_a}} \right)}^*} \otimes {M_R}\left( {{{{\boldsymbol{q}}}_{i\left( j \right)}} \otimes {{{\boldsymbol{q}}}_a}} \right)} \right]}^{\rm{T}}}{I_j}{{\boldsymbol\omega} _j}} \Bigg\}- \\ & \quad {m_j}{\left(\frac{{\partial {{{\boldsymbol{v}}}_j}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}}{{{\boldsymbol{v}}}_j} -{\left(\frac{{\partial {{\boldsymbol\omega} _j}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}}{I_j}{{\boldsymbol\omega} _j} + \\ &\quad{\left(\frac{{\partial {P_{gj}}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}} + {\left(\frac{{\partial {P_{ej}}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}} = {{\boldsymbol{G}}_j}\\[-10pt] \end{split} $$ (15) 其中,
$ T_a $ 是$ T\left( {{{\boldsymbol{q}}_a}, {{\boldsymbol{L}}_a}} \right) $ 的简写,$ T_{i(j)} $ 是$ T\left( {{{\boldsymbol{q}}_{i(j)}}, {{\boldsymbol{L}}_{i(j)}}} \right) $ 的简写, 而$ T_j $ 是$ T\left( {{{\boldsymbol{q}}_j}, {{\boldsymbol{L}}_j}} \right) $ 的简写. 考虑到次级关节动力学速度与角速度的具体形式:$$ {{\boldsymbol{v}}_j} = \frac{\rm{d}}{{{\rm{d}}t}}\left( {{\boldsymbol{q}}_p} + T_a T_{i(j)}T_j{\boldsymbol{D}}_j \right)\qquad \;\qquad\qquad$$ (16) $${{\boldsymbol{\omega}} _j} = 2{\left( {{{\boldsymbol{q}}_j} \otimes {{\boldsymbol{q}}_{i\left( j\right)}} \otimes {{\boldsymbol{q}}_a}} \right)^*} \otimes \frac{{\rm{d}}} {{{\rm{d}}t}}\left( {{{\boldsymbol{q}}_j} \otimes {{\boldsymbol{q}}_{i\left( j \right)}} \otimes {{\boldsymbol{q}}_a}} \right) $$ (17) 可知次级连杆动力学(11)杂糅了次级连杆的广义关节坐标加速度项
$ \ddot{{\boldsymbol{q}}}_j $ , 对应的初级连杆的广义关节坐标加速度项$ \ddot{{\boldsymbol{q}}}_{i\left( j\right)} $ , 和机体坐标角速度项$ \ddot{{\boldsymbol{q}}}_a $ 和$ \ddot{{\boldsymbol{q}}}_p $ . 与此同时, 该动力学方程中也受初级关节和机体运动所产生的惯性项影响. 由于翅翼质量的占比较少, 可以设置${m_j} = \varepsilon {m_{\varepsilon j}} ,$ ${I_j} = \varepsilon {I_{\varepsilon j}},\; j \in B,$ 使得$ m_{\varepsilon j} $ 和$m_{1},$ $ I_{\varepsilon j} $ 和$ I_{1} $ 近似同量级, 其中$ \varepsilon $ 可以视为表示时间尺度的正值常数. 这样实际上在慢子系统中隔绝了惯性力和重力的影响, 当设置$ \varepsilon = 0 $ 时, 有$$ \left(\frac{{\partial P_{ej}}} {{\partial {{\boldsymbol{q}}_j}}}\right)^{\rm{T}} = {\bar{{\boldsymbol{G}}}_j}$$ (18) 其中, 标称广义力
$ {\bar{{\boldsymbol{G}}}_j} $ 由广义摩擦力、广义气动力和广义输入构成, 具体而言, 在准定常气动假设下的广义气动力由系统运动学和外部环境共同决定, 而广义输入由关节输入力矩和系统运动学共同决定. 此外,${\bar{{\boldsymbol{G}}}_i} ,$ ${\bar{{\boldsymbol{G}}}_p} ,$ $ {\bar{{\boldsymbol{G}}}_a} $ 均以此种方式构成. 接下来考虑初级连杆的动力学方程:$$ \begin{split} &\frac{{\rm{d}}}{{{\rm{d}}t}}\left\{ {{m_i}{{\left( {{T_a}\frac{{\partial {{\dot T}_i}{{\boldsymbol{D}}_i}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}} \right)}^{\rm{T}}}{{\boldsymbol{v}}_i}}\right.+\\ &\quad\left. { 2{{\left[ {{{\left( {{{{\boldsymbol{q}}}_j} \otimes {{{\boldsymbol{q}}}_a}} \right)}^*} \otimes {M_R}\left( {{{{\boldsymbol{q}}}_a}} \right)} \right]}^{\rm{T}}}{I_j}{{\boldsymbol{\omega}} _j}}\right.+\\ &\quad\left. { \sum\limits_{j \in B,i\left( j \right) = i} {\left( {{m_j}{\left({\frac{{\partial {{\boldsymbol{v}}_j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_j} + {\left({\frac{{\partial {{\boldsymbol{\omega}} _j}}}{{\partial {{\dot {{\boldsymbol{q}}}}_i}}}}\right)^{\rm{T}}}{I_j}{{\boldsymbol{\omega}} _j}} \right)} } \right\} -\\ &\quad \left( {{m_i}{\left({\frac{{\partial {{\boldsymbol{v}}_i}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}}}{{\boldsymbol{v}}_i}{ + }{\left({\frac{{\partial {{\boldsymbol{\omega}} _i}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}}}{I_i}{{\boldsymbol{\omega}} _i} + {{\boldsymbol{S}}_i}} \right)+ \left(\frac{{\partial {P_{gi}}}}{{\partial {{\boldsymbol{q}}_i}}}\right)^{\rm{T}}+\\ & \quad \sum\limits_{j \in B, i\left( j \right) = i} \left({\frac{{\partial {P_{gj}}}}{{\partial {{\boldsymbol{q}}_i}}}}\right)^{\rm{T}} + \left(\frac{{\partial {P_{ei}}}}{{\partial {{\boldsymbol{q}}_i}}}\right)^{\rm{T}} = {{\boldsymbol{G}}_i} \\[-15pt]\end{split} $$ (19) 同样基于翅翼质量占比较少, 且同样设置
$ {m_i} = \varepsilon {m_{\varepsilon i}} $ ,$ {I_i} = \varepsilon {I_{\varepsilon i}}, i \in A $ , 当$ \varepsilon = 0 $ 时, 有$$\left(\frac{{\partial P_{ei}}} {{\partial {{\boldsymbol{q}}_i}}}\right)^{\rm{T}} = {\bar{{\boldsymbol{G}}}_i} $$ (20) 基于同样假设, 考虑躯干动力学方程简化模型
$$m_1 \ddot{{\boldsymbol{q}}}_p+ \left(\frac{{\partial {P_g}}}{{\partial {{\boldsymbol{q}}_p}}}\right)^{\rm{T}} = {{\boldsymbol{G}}_{p}} $$ (21) $$ \begin{split} &\frac{{\rm{d}}}{{{\rm{d}}t}}\left( {4{{\left( {M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right)} \right)}^{\rm{T}}}{I_1}M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right){{\dot {{\boldsymbol{q}}}}_a} } \right) -\\ & \qquad {4{{\left( {\frac{{\partial M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right){{\dot {{\boldsymbol{q}}}}_a}}}{{\partial {{\boldsymbol{q}}_a}}}} \right)}^{\rm{T}}}{I_1}M_L^{ - 1}\left( {{{\boldsymbol{q}}_a}} \right){{\dot {{\boldsymbol{q}}}}_a} } = {{\boldsymbol{G}}_{a}} \end{split} $$ (22) 事实上, 在躯干动力学层面的摄动问题并没有奇异性, 作为扰动的惯性项仅作为常规扰动出现. 值得指出的是,
$ {\boldsymbol{G}}_p $ 和$ {\boldsymbol{G}}_a $ 中的可以作为输入给定的部分为初级连杆根部的关节力矩, 而非作用于翅翼的气动力. 显然, 由简化模型(18), (20) ~ (22)可以共同构成扑翼飞行的动力学简化模型, 也就是慢子系统. 在给定稳定外部环境时, 假设可以在去除$ {\boldsymbol{q}}_p $ 的广义关节的状态空间(包括其广义关节速度项, 且包括$\dot{{\boldsymbol{q}}}_p )$ 依此动力学模型找到对应某闭合解对应的广义输入, 此闭合解可能是闭合轨线或者临界点. 直观而言, 闭合轨线对应着周期稳定的扑翼飞行, 而临界点对应着稳定的滑翔飞行. 这里着重考虑闭合轨线的情况, 并在相空间将此轨线描述为$ \gamma $ . 沿此轨线我们首先描述次级连杆的边界层动力学:$$ \begin{split} &\frac{{\rm{d}}} {{{\rm{d}}t}}\Bigg\{ { \varepsilon {m_{\varepsilon j}}\left( {T_a}{T_{i\left( j \right)}}\frac{{\partial {{\dot T}_j}{D_j}}} {{\partial {{\dot {{\boldsymbol{q}}}}_j}}}\right)^{\rm{T}} {{{\boldsymbol{v}}}_j} }+\\ &\quad { 2\varepsilon {{\left[ {{{\left( {{{{\boldsymbol{q}}}_j} \otimes {{{\boldsymbol{q}}}_{i\left( j \right)}} \otimes {{{\boldsymbol{q}}}_a}} \right)}^*} \otimes {M_R}\left( {{{{\boldsymbol{q}}}_{i\left( j \right)}} \otimes {{{\boldsymbol{q}}}_a}} \right)} \right]}^{\rm{T}}} {I_{\varepsilon j}}{{\boldsymbol\omega} _j}} \Bigg\}_{\bot \gamma } - \\ & \quad\varepsilon {m_{\varepsilon j}}{\left(\frac{{\partial {{{\boldsymbol{v}}}_j}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}}{{{\boldsymbol{v}}}_j} - \varepsilon {\left(\frac{{\partial {{\boldsymbol\omega} _j}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}} {I_{\varepsilon j}}{{\boldsymbol\omega} _j} + \\ &\quad \varepsilon{\left(\frac{{\partial {P_{\varepsilon gj}}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}} + {\left(\frac{{\partial {P_{ej}}}} {{\partial {{{\boldsymbol{q}}}_j}}}\right)^{\rm{T}}_{\bot \gamma }} = {{\boldsymbol{G}}_{\varepsilon j}} \\[-15pt]\end{split} $$ (23) 其中, 下标
$ \star_{\bot \gamma} $ 表示$ \star $ 此量相对于确定闭合轨线$ \gamma $ 动力学之外的部分, 广义力$ {{\boldsymbol{G}}_{\varepsilon j}} = {\boldsymbol{G}}_{j{\bot \gamma }} $ 则表示广义控制输入的高频分量, 和其他边界层系统的复杂耗散力之和. 显然, 对于初级连杆和躯干也可以获得其相对应的边界层动力学, 其通过线性运算可以进一步获得各广义坐标的动力学.因为飞行器多刚体结构的复杂性、弹性势能场的复杂性、以及各种复杂耗散力的综合作用, 探讨上述边界层系统的具体动力学特性极其复杂, 需要对问题进行抽象简化. 具体而言, 简化系统和边界层系统为
$$ \dot {\bar{{\boldsymbol{q}}}} = {{\boldsymbol{f}}_G}\left(t, {\bar {{\boldsymbol{q}}},\tilde{{\boldsymbol{q}}},\varepsilon } \right) $$ (24) $$ \varepsilon {{\dot {\tilde {{\boldsymbol{q}}}}}} = {{\boldsymbol{g}}_G}\left(t, {\bar {{\boldsymbol{q}}},\tilde {{\boldsymbol{q}}},\varepsilon } \right) $$ (25) 其中,
$ \bar{{\boldsymbol{q}}} $ 是对应于$ \gamma $ 的标称广义坐标状态光滑周期解, 而$ \tilde{{\boldsymbol{q}}} $ 是真实广义坐标状态$ {{\boldsymbol{q}}_s} $ 和轨线$ \gamma $ 之间在广义坐标状态空间的距离向量中模数最小者, 且可以由$ \tilde{{\boldsymbol{q}}} $ 和$ \gamma $ 重构出${{\boldsymbol{q}}} ,$ $ {{\boldsymbol{f}}_G} $ 和$ {{\boldsymbol{g}}_G} $ 则是对应于广义力$ G $ 的动力学. 因为研究扑翼飞行器的具体特性预设固定, 所以这里$ G $ 代表着不同控制策略与环境作用的综合. 值得指出的是, 这里所谓的广义坐标状态$ {{\boldsymbol{q}}_s} $ 既包括广义坐标也包括广义坐标的导数.综上所述, 这是一个典型的周期系统奇异摄动问题, 可以使用几何奇异摄动理论给出扰动系统的周期解存在性和稳定性的条件[38-39], 同时也可以参考文献[40]中定理6.6.1和定理6.6.2. 如要获得对应于整体奇异摄动系统(24)和(25)的且微分同胚于闭合轨线
$ \gamma $ 的稳定不变流形$ M_\varepsilon $ , 则需要轨线$ \gamma $ 上每一点关于摄动系统的线性化都需要其对应的雅克比矩阵特征值实部不为零, 且其关于简化系统至少是轨线$ \gamma $ 邻域上局部吸引的. 这显然需要闭合轨线附近的动力学满足: 对于所设计控制器和环境综合作用$ G $ , 可以在任一点检测上述条件. 然而考虑到实际情况, 这是极其困难的. 常见扑翼飞行器的控制问题中既不能完成周期内的精细调控, 也不可能获得准确的飞行器动力学.本文采用一种类似频闪的方式来观测扑翼飞行, 实际实现中, 即基于固定周期
$ T $ 的观测对系统进行奇异摄动分析. 值得注意的是, 这里$ T $ 不必要是闭合轨线的周期. 这种方法需要假设已知非摄动系统的先验信息, 也即需要从给定的扑翼飞行器控制策略到闭合轨线$ \gamma $ 的映射关系, 此映射关系可以由简化模型(18), (20) ~ (22)基于准定常气动力模型构建, 也可从实际实验的观测中归纳. 以下分析假设扑翼飞行器控制策略、控制目标和环境影响固定不变, 也即$ \gamma $ 固定. 那么对于广义坐标状态空间$ Q_S $ 上的每一点$ {{\boldsymbol{q}}_s} $ 可以定义其到$ \gamma $ 距离为$$ d\left( {{{\boldsymbol{q}}_s},\gamma } \right) = \mathop {\inf }\limits_{{\boldsymbol{x}} \in \gamma } d\left( {{{\boldsymbol{q}}_s},{\boldsymbol{x}}} \right) $$ (26) 其中,
$ d $ 为定义在$ Q_s $ 上的某确定度量. 根据在离散点$t = T,2T, \cdots, kT$ 的观测可以进一步获得离散点到$ \gamma $ 的距离$ d\left( {{{\boldsymbol{q}}_s}\left( {kT} \right),\gamma } \right) $ , 简记为$ d_k $ . 基于上述观测, 可以获得离散的动力学系统:$$ {d_{k + 1}} = h\left( {{d_k}} \right) $$ (27) 其中,
$ h $ 表示频闪观测的系统$ \gamma $ - 轨道距离状态变换函数, 是由摄动系统动力学、观测时间间隔和闭合轨线$ \gamma $ 共同决定的.定理 1. 对于上述奇异摄动问题, 如果简化系统的周期轨迹
$ \gamma $ 是渐近稳定且孤立的, 同时频闪观测间隔$ T $ 是$ {\cal{O}} (\varepsilon) $ 的, 且对于离散系统(27)有$ {d_k} \le $ $ W\left( {kT} \right) $ , 其中$ {W} $ 是正定连续可微函数, 其关于时间的导数$\dot W\left( {kT} \right) \le - {\lambda} {{d_k}} + {\beta _\gamma },$ 正值常数$ \lambda $ 是$ {\cal{O}} (1) $ 的, 而正值常数$ {\beta _\gamma } $ 是$ {\cal{O}} (\varepsilon) $ 的, 那么存在$ \varepsilon_0 $ 对于所有$ \varepsilon < \varepsilon_0 $ , 设若奇异摄动系统(24)和(25)有在$ \gamma $ 附近的的一个或一族周期解①, 取其中任意解的对应轨迹为$ \gamma_\varepsilon $ , 奇异摄动系统会收敛到$ \gamma_\varepsilon $ 的$ {\cal{O}} (\varepsilon) $ 邻域内.证明. 此证明重点考虑频闪观测和Poincare映射在反映连续系统行为上的相似性. 由所述条件可知在
$ k\to \infty $ ,$ {d_{k }} $ 最终收敛到$ {\cal{O}} (\varepsilon) $ 之内. 考虑到观测间隔$ T $ 是$ {\cal{O}} (\varepsilon) $ 的, 由系统的连续性可知在两次观测之间的系统状态, 在上述收敛的前提下也有类似性质, 具体而言, 对于$ t\to \infty $ , 有$ d\left( {{q_s}\left( t \right),\gamma } \right) $ 最终收敛到$ {\cal{O}} (\varepsilon) $ , 也即系统收敛到$ \gamma $ 的$ {\cal{O}} (\varepsilon) $ 邻域内. 因为此邻域包含系统的唯一前向不变集, 所以周期轨线$ \gamma_\varepsilon $ 必然包含于此集中. 进而闭合轨线$ \gamma $ 和$ \gamma_\varepsilon $ 之间的Hausdorff距离是$ {\cal{O}} (\varepsilon) $ 的, 由度量的定义, 进而可知奇异摄动系统会收敛到$ \gamma_\varepsilon $ 的$ {\cal{O}} (\varepsilon) $ 邻域内.对于所述奇异摄动问题, 当
$ \varepsilon\to 0 $ 时, 有$ \gamma_\varepsilon\to\gamma $ , 其对应的周期$ \tau(\gamma_\varepsilon)\to\tau(\gamma) $ . 上述定理阐述了频闪观测结果和系统动力学行为的关系, 说明在所述条件下基于简化模型设计周期输入控制效果的准确性, 同时也未对复杂的边界层系统提出除光滑性和有界性等物理系统常见性质之外的要求. 周期广义力输入和周期非摄动系统轨迹的关系可以用数值方法, 根据实验结果或者扑翼飞行准定常动力学模型给出近似, 而不用根据动力学模型进行在线解算. 我们将进一步结合图3, 纲要性地描述上述频闪观测等方法对应的控制策略: 在周期性广义力的作用下, 控制器依慢子系统动力学所确定的周期轨迹控制输入周期变化, 其系统在轨迹周围存在一个$ {\cal{O}} (\varepsilon) $ 误差容许范围. 当频闪观测到系统偏离此范围时, 控制器采取不同的周期性的广义力促使其恢复原状态. 值得指出的是, 控制器同样基于慢子系统动力学来预期这个周期性的广义力的作用效果.
图 3 扑动翅翼周期轨道控制示意图Fig. 3 The schematic for the control of flapping-wing periodic orbit本文所述的奇异摄动方法基于慢子系统动力学和频闪观测进行分析, 并在一定范围内确定整体摄动系统的动力学性质, 其具体优势如下:
1)采用频闪观测的方法, 固定观测时间间隔, 便于机载硬件实现;
2)即便广义力输入不具有精确的解析形式, 所采用的分析方法也可以通过数值方法来分析系统的稳定性;
3)相对于平均化等方法, 所设计方法考虑系统周期力或力矩输入对翅翼动力学和其周期运动步态的影响, 为完成更灵活的扑翼飞行和更复杂的运动控制任务铺平道路.
3. 扑翼飞行实验与分析
这一节我们分析自制的四自由度扑翼飞行器在悬停时的具体行为, 以验证所述基于奇异摄动理论分析方法的有效性. 在悬停时翅翼周期往复拍动, 躯干受翅翼运动所致在小范围内振动. 如图4所示, 系统共存在四个控制输入, 其中空心杯电机驱动X-型扑翼机构, 三个直线舵机分别驱动两个翅翼的转角和一个尾翼舵机转角. 至于飞行器的具体实验设置和控制方法, 可以参考文献[41]. 此外, 还对该机型在测力/力矩传感器上进行了对其偏航力矩的实际测量. 上述扑翼飞行器在实际实验中, 存在翅翼运动和机体运动时标分离现象. 如图5所示, 我们在Qualisys 动作捕捉系统下做了所述机型的悬停实验, 使用文献[41]中所述姿态控制器, 利用左右翅翼前后偏转的共同作用产生俯仰力矩, 利用尾翼偏转产生滚转力矩. 在控制器保持
$ z $ 轴竖直的情况下, 同时将左右翅翼转角的差值设置为最大幅值的30%, 使得飞行器获得绕$ z $ 轴的转矩. 在此控制器设置下, 对左翼靠前和右翼靠前两种情况分别设置了两组不同的实验. 根据系统输出的记录和动作捕捉系统的姿态反馈, 我们获得如图6所示的飞行数据. 其中$\phi,$ $ \theta $ 和$ \psi $ 分别表示飞行器的俯仰、滚转和偏航, FA表示扑动翅翼驱动电机的电压输入, LW, RW和VT分别表示左翼、右翼和尾翼的驱动舵机偏转.
图 4 四自由度扑翼飞行器示意图 ((a)惯性坐标系与机体坐标系设置; (b)翅翼驱动与控制机构; (c1)无偏转的翅翼姿态; (c2)偏转后的翅翼姿态; (d1) 无偏转的尾翼姿态; (d2) 偏转后的尾翼姿态)Fig. 4 The schematic for the 4-DoF flapping wing vehicle ((a) Inertial coordinate system and body coordinate system setup; (b) Flapping wings actuation and control mechanism; (c1) Flapping wings in default state; (c2) Flapping wings in deflection state; (d1) Tail in default state; (d2) Tail in deflection state)
图 5 运动捕捉系统辅助的四自由度扑翼飞行器实际飞行Fig. 5 The real flight of the developed 4-DoF flapping wing system with the help of the motion capture system
图 6 不同两翼转角差下悬停飞行时的控制输入与机体姿态Fig. 6 The control inputs and fuselage attitudes in the hovering flight with different two-wings deflection-angle deviations除此之外, 如图7所示, 利用测力平台测得飞行器在不同翅翼转角设置下所产生的对应气动力矩. 具体而言, 确定驱动扑翼机构的空心杯电机电压为 3 V, 将左右翅翼转角的差值设置为最大幅值的30%, 定义翅翼转角偏向机体系
$ x $ 轴方向为负, 反之为正. ATI Nano 17六维力/力矩传感器采用 SI-12-0.12 标定参数, 采样频率为 7 000 Hz, 其测力精度约为 1/320 N, 测力矩精度约为 1/64 N·mm. 这里, 由于测力平台坐标系$ {\cal{F}} _S $ 和飞行器机体坐标系$ {\cal{F}} _B $ 不重合, 需要利用平行轴定理进行修正, 其具体计算方法可以参考文献[42]. 实验分别独立测量了左右翅翼转角差值为最大转角幅值 0.3 和 −0.3倍时的$ z $ 轴力矩, 由于系统装配的误差, 在转角差值为 0 时依然有平均值为 1.6 左右的$ z $ 轴力矩, 所以两个翅翼转角差值所对应的周期整体平均转矩, 在减去偏置后实际上是近似大小相同而方向相反的.根据图6和 图8所示的实验结果进一步分析, 可以得到如下结论:
图 8 由翅翼偏转引起的绕z轴周期变化的力矩 ((a)在不同偏转角下力矩的周期变化; (b)转角差值为−0.3倍最大转角的力矩多周期统计特性; (c)转角差值为0.3倍最大转角时的力矩多周期统计特性)Fig. 8 The periodically changing torque around z-axis induced by the wings deflection ((a) Periodic variation of torque under different deflection angles; (b) Statistical characteristics of the torque with the deflection angle of −0.3 limit angle; (c) Statistical characteristics of the torque with the deflection angle of 0.3 limit angle)1)因为所设计控制器基于第2节中所述的慢子系统, 而飞行器整体动力学对应着包含快慢子系统的摄动系统, 实验完成的稳定飞行从一定程度上说明: 在快慢子系统时间尺度分离较为完整的情况下, 尽管如图8所示的快子系统振动中的不能被平缓周期力反映的高频波动部分会在一定程度上影响整体系统, 但以扑翼慢子系统稳定性为设计目标的控制器依旧适用于包含快慢子系统的摄动系统;
2)实验中的扑翼飞行通过动作捕捉系统的100 Hz频闪观测来确定系统状态, 尽管其不能观测翅翼状态, 同时可观测部分也存在一定误差, 但实验结果也足可说明当频闪观测足够快速时, 其可以有效反映包含快慢子系统的摄动系统整体运动状态, 并能够为基于慢子系统设计的控制器提供有效反馈信息;
3)通常而言, 扑翼飞行器常基于平均化理论进行其动力学分析和控制器设计, 简单而言, 上述两者多基于平滑的周期平均系统而非原有的复杂周期变化系统进行稳定性分析. 具体而言, 使用平均化方法认为系统使用图8所示的平均力或力矩而非周期变化的力或力矩作为控制输入, 并基于此做针对单刚体机体的控制器设计. 采用此种方法便意味着忽略系统的周期性特征, 从而无法完成周期内具体特性的控制, 比如本文所述的系统收敛到某一确定周期轨迹附近. 与此不同的是, 本文所述的系统将多刚体模型下的周期运动视为慢动态, 而将其他高频摄动视为快动态, 从而在一定程度反映了扑翼系统的周期特性, 较原有方法更为精细. 尽管视觉捕捉系统不能完全反应周期系统的轨迹, 但是其频闪采样的数据依旧可以在一定程度上刻画周期运动特性.
4. 结论
本文首先基于常见的扑翼飞行系统, 建立了扑翼飞行器多刚体模型, 描述了各刚体之间的动力学关系. 然后, 本文还通过应用奇异摄动理论来分析周期系统所对应的离散系统, 在翅翼和机体运动存在时标分离的前提下, 分析了扑翼系统的周期稳定性. 其中, 针对已有方法在动力学分析和设计扑翼控制器时忽略翅翼运动的问题, 以及其难以应用于翅翼周期运动控制的不足, 本文详细讨论了翅翼动力学的周期稳定性问题, 并说明了通过奇异摄动方法分析翅翼动力学的必要性. 最后, 结合上述理论对扑翼飞行实验结果进行了分析, 在一定程度上验证了文中所提方法的有效性. 在后续研究中, 我们将继续探究奇异摄动理论在扑翼周期问题中的应用, 并进一步发掘扑翼飞行器的机动性潜力.
-
图 1 一类典型的流程工业多反应器串联系统信号流图
Fig. 1 A typical signal flow graph of multi-reactor series system in process industry
图 2 输出反馈闭环控制系统结构框图
Fig. 2 The block diagram of output-feedback closed-loop control system
图 3 状态反馈控制下的被控输出$y\left(t\right)$及参考信号$y_r\left(t\right) = 2\cos \left(t\right)$
Fig. 3 The controlled output $y\left(t\right)$ under the state-feedback control and the reference signal $y_r\left(t\right) = 2\cos \left(t\right)$
图 4 状态反馈控制器$U\left(t\right)$轨迹
Fig. 4 The trajectory of state-feedback controller $U\left(t\right)$
图 5 输出反馈控制下的被控输出$y\left(t\right)$及参考信号$y_r\left(t\right) = 2\cos \left(t\right)$
Fig. 5 The controlled output $y\left(t\right)$ under the output-feedback control and the reference signal $y_r\left(t\right) = 2\cos \left(t\right)$
图 6 输出反馈控制器$\hat U\left(t\right)$轨迹
Fig. 6 The trajectory of output-feedback controller $\hat U\left(t\right)$
图 7 PDE子系统观测误差的范数$\left\| {{{\tilde u}_1}} \right\|$和$\left\| {{{\tilde u}_2}} \right\|$
Fig. 7 The norms $\left\| {{{\tilde u}_1}} \right\|$ and $\left\| {{{\tilde u}_2}} \right\|$ of observer errors of PDE subsystem
图 8 执行器状态观测误差${X_1} - {{\hat X}_1}$和${X_2} - {{\hat X}_2}$
Fig. 8 The observer errors ${X_1} - {{\hat X}_1}$ and ${X_2} - {{\hat X}_2}$ of actuator states
图 9 PDE子系统受到的常值扰动${D_i} = Q_i^\mathrm{T}{v_{d1}},i = 1, 2 ,3$及相应的扰动观测值${\hat D_i} = Q_i^\mathrm{T}{\hat v_{d1}},i = 1, 2 ,3$
Fig. 9 The constant perturbations ${D_i} = Q_i^\mathrm{T} {v_{d1}}$, $i = 1, 2 ,3$ to PDE subsystem and the corresponding disturbance observations ${\hat D_i} = Q_i^\mathrm{T}{\hat v_{d1}},i = 1, 2 ,3$
图 10 执行器受到的周期性扰动${d_i} = q_i^\mathrm{T}{v_{d2}},i = 0,1 ,2$以及相应的扰动观测值${{\hat d}_i} = q_i^\mathrm{T}{{\hat v}_{d2}},i = 0,1 ,2$
Fig. 10 The periodic perturbations ${d_i} = q_i^\mathrm{T} {v_{d2}}$, $i = 0, 1 ,2$ to actuator and the corresponding disturbance observations ${{\hat d}_i} = q_i^\mathrm{T}{{\hat v}_{d2}}$, $i = 0, 1 ,2$
图 11 采用控制器(79)的系统被控输出 $y\left(t\right)$及参考信号$y_r\left(t\right) = 2\cos \left(t\right)$
Fig. 11 The system controlled output $y\left(t\right)$ and reference signal $y_r\left(t\right) = 2\cos \left(t\right)$ using the controller (79)
-
[1] Xu C, Sallet G. Exponential stability and transfer functions of processes governed by symmetric hyperbolic systems. ESAIM: Control, Optimisation and Calculus of Variations, 2002, 7: 421-442 doi: 10.1051/cocv:2002062 [2] Landet I S, Pavlov A, Aamo O M. Modeling and control of heave-induced pressure fluctuations in managed pressure drilling. IEEE Transaction Control System Technology, 2012, 21(4): 1340-1351 [3] Boskovic D M, Krstic M. Backstepping control of chemical tubular reactors. Computers & Chemical Engineering, 2002, 26(7-8): 1077-1085 [4] Deutscher J. Output regulation for general linear heterodirectional hyperbolic systems with spatially-varying coefficients. Automatica, 2017, 85: 34-42 doi: 10.1016/j.automatica.2017.07.027 [5] Xu X, Dubljevic S. Output regulation boundary control of first-order coupled linear mimo hyperbolic pide systems. International Journal of Control, 2020, 93(3): 410-23 doi: 10.1080/00207179.2018.1475749 [6] Deutscher J, Gehring N, Kern R. Output feedback control of general linear heterodirectional hyperbolic ode-pde-ode systems. Automatica, 2018, 95: 472-480 doi: 10.1016/j.automatica.2018.06.021 [7] Xu X, Dubljevic S. Output regulation for a class of linear boundary controlled first-order hyperbolic pide systems. Automatica, 2017, 85: 43-52 doi: 10.1016/j.automatica.2017.07.036 [8] Deutscher J, Gabriel J. Minimum time output regulation for general linear heterodirectional hyperbolic systems. International Journal of Control, 2020, 93(8): 1826-1838 doi: 10.1080/00207179.2018.1533648 [9] Zhang J, Qi J, Dubljevic S, Bo S. Output regulation for a first-order hyperbolic PIDE with state and sensor delays. European Journal of Control, 2022, 65: Article No. 100643 [10] Irscheid A, Deutscher J, Gehring N, Joachim R. Output regulation for general heterodirectional linear hyperbolic PDEs coupled with nonlinear ODEs. Automatica, 2023, 148: Article No. 110748 [11] Deutscher J. Finite-time output regulation for linear 2×2 hyperbolic systems using backstepping. Automatica, 2017, 75: 54-62 doi: 10.1016/j.automatica.2016.09.020 [12] Owens B A, Mann B P. Linear and nonlinear electromagnetic coupling models in vibration-based energy harvesting. Journal of Sound and Vibration, 2012, 331(4): 922-937 doi: 10.1016/j.jsv.2011.10.026 [13] Susto G A, Krstic M. Control of pde-ode cascades with neumann interconnections. Journal of Franklin Institute, 2010, 347(1): 284-314 doi: 10.1016/j.jfranklin.2009.09.005 [14] Xu X, Tian Y, Yuan Y, Luan X, Liu F, Dubljevic S. Output regulation of linearized column froth flotation process. IEEE Transaction Control System Technology, 2020 29(1): 249-262 [15] Li J, Wu Z, Liu Y. Adaptive stabilization for an uncertain reaction-diffusion equation with dynamic boundary condition at control end. System & Control Letters, 2022, 162, 105-180 [16] Liu W J, Krstic M. Backstepping boundary control of burgers’ equation with actuator dynamics. System & Control Letters, 2000, 41(4): 291-303 [17] Wang J, Krstic M. Output-feedback control of an extended class of sandwiched hyperbolic pde-ode systems. IEEE Transaction on Automatic Control, 2020, 66(6): 2588-2603 [18] Wang J, Krstic M. Delay-compensated control of sandwiched ode-pde-ode hyperbolic systems for oil drilling and disaster relief. Automatica, 2020, 120: 109-131 [19] Wang J, Krstic M. Event-triggered output-feedback backstepping control of sandwich hyperbolic pde systems. IEEE Transaction Automatic Control, 2022, 67(1): 220-235 doi: 10.1109/TAC.2021.3050447 [20] Wang J, Krstic M. Output feedback boundary control of a heat pde sandwiched between two odes. IEEE Transaction Automatic Control, 2019, 64(11): 4653-4660 doi: 10.1109/TAC.2019.2901704 [21] Li J, Wu Z, Wen C. Adaptive stabilization for a reaction-diffusion equation with uncertain nonlinear actuator dynamics. Automatica, 2021, 128: 109-594 [22] Xiao Y, Yuan Y, Yang C, Luo B, Xu X, Dubljevic S. Adaptive neural tracking control of a class of hyperbolic PDE with uncertain actuator dynamics. IEEE Transactions on Cybernetics, DOI: 10.1109/TCYB.2022.3223168 [23] Zhang B, Yang C, Zhu H, Shi P, Gui W. Controllable-domain-based fuzzy rule extraction for copper removal process control. IEEE transactions on fuzzy systems, 2017 26(3): 1744-1756 [24] Lin Z, Stoorvogel A A, Saberi A. Output regulation for linear systems subject to input saturation. Automatica, 1996, 32(1): 29-47 doi: 10.1016/0005-1098(95)00110-7 [25] Raghavan S, Hedrick J K. Observer design for a class of nonlinear systems. International Journal of Control, 1994, 59(2): 515-528 doi: 10.1080/00207179408923090 [26] Vazquez R, Krstic M, Coron J M. Backstepping boundary stabilization and state estimation of a 2 × 2 linear hyperbolic system. In: Proceedings of Conference on Decision and Control and European Control Conference. Orlando, USA: 2011. 4937−4942 [27] Kailath T. Linear Systems. Englewood Cliffs, NJ: Prentice Hall, 1980. [28] Khalil H. Nonlinear systems (3rd edition). NJ: Prentice Hall, 2002. -
下载:
下载:
计量
- 文章访问数: 506
- HTML全文浏览量: 270
- PDF下载量: 179
- 被引次数: 0
