一种基于模态参数实时辨识方法的参数时变航天器控制方法*

2019-12-03贾贵鹏赵育善

飞控与探测 2019年6期

贾贵鹏，赵欣，赵育善，师鹏

(北京航空航天大学宇航学院·北京·100191)

0 引言

对于存在刚柔耦合航天器的姿态运动，为实现对姿态的控制，通常需要对系统状态量进行估计，获取系统尽可能精确的数学模型，但多数情况下，系统的数学精确模型不能通过理论计算获得，甚至系统的参数存在时变，系统参数辨识就是确定系统数学模型的理论方法[1-2]。通过实验分析进行参数辨识即通过选取结构上测试点，采集测试点的动态数据，通过特定的估计和辨识方法来获取系统模型及参数。在辨识获取系统姿态和模态状态量后，即可利用无模型控制方法对航天器姿态进行控制，而不需要受控对象的其他模型参数。

在航天领域，为解决大型挠性航天器具有低频密集的特点，时域辨识成为航天器参数辨识的主要方法，时域法无需激励数据，直接利用响应数据辨识系统参数，适合对系统进行在线参数辨识，获得系统动态性能[3-7]。 Liu[8]利用基于奇异值分解(Singular Value Decomposition，SVD)方法获取时变系统的参数估计，并提出伪模态方法来解决时变系统模态分析问题。Shao等[9]将仿射投影算法(Affine Projection Algorithm，APA)和符号算法(Sign Algorithms，SAS)的优点结合，通过使用多个投影同时根据L1范数优化标准更新其权重，提出了仿射投影符号算法(Affine Projection Sign Algorithm，APSA)，该方法计算复杂度较低，抗脉冲干扰性强，收敛速度快。Houtzager等[10]提出基于递归预测器的子空间辨识法(Recursive Predictor-Based Subspace Identification，RPBSID)，通过自适应滤波，利用外生向量自回归(Vector Autoregressive with Exogenous，VARX)预测器对系统参数进行辨识。针对RPBSID方法在对状态量进行估计时需要构建广义Hankel矩阵而导致计算量较大的问题，倪智宇等[11]将RPBSID方法和APSA方法进行结合，利用APSA方法对系统状态量进行估计，提高了计算效率。无模型控制是侯忠生提出的一种控制方法[12-14]，其主要特点是不需要系统的动力学模型，只需要系统的激励和响应量即可对系统实现控制，上述研究对象集中在机械、化工等工业领域，对于复杂航天器控制问题研究较少，本文针对一种因挠性结构转动引起模态参数变化的航天器的模态辨识与控制问题进行了研究。

本文首先建立参数时变刚柔耦合航天器的姿态动力学模型，其次利用文献[11]的参数递推辨识方法估算系统状态参数，最后利用无模型控制方法对该系统进行控制，在不需要获得系统其他模型参数的情况下可实现对系统的控制。

1 数学模型的建立

1.1 刚柔耦合航天器姿态时变动力学方程

对于带有大型柔性帆板或柔性天线的航天器的建模问题，常用的方式为将其简化为一个中心刚体附加若干柔性附件的模型结构，假设航天器带有两个可转动的太阳能帆板和两个具有桁架结构的天线，如图1所示。

图1 航天器本体和附件坐标系示意图Fig.1 Body coordinate system of the spacecraft

航天器本体为刚体，附件1、2为桁架结构的天线，附件3、4为太阳能帆板，帆板旋转角标记为δ，当帆板面向xb轴方向时δ=0。假设帆板以很小的角速度旋转以保持供能，可得到简化的动力学方程[12]：

(1)

(2)

(3)

式中，r=[rx,ry,rz]T为 J2000惯性坐标系下系统质心的位置向量，θ=[θroll,θpitch,θyaw]为姿态角矢量，ηi∈Rni为第i个挠性附件的模态坐标，Zi=diag[ζai,1,…,ζai,2]∈Rni×ni为对应ηi的阻尼比矩阵，Ωi=diag[ωai,1,…,ωai,2]∈Rni×ni为对应ηi的模态频率矩阵，F∈R3为作用在航天器系统上的力，T∈R3为作用在航天器系统上的力矩，M,J,Υ,Γi,Pi的定义如下：

由于本文只研究刚柔耦合航天器的姿态控制，因此认为力F很小，航天器的线性运动对振动的影响也可以忽略不计，忽略航天器的线性运动，由式(2)、(3)可得以下航天器动力学方程：

(4)

(5)

其中，Λi=2ZiΩi。假设帆板角度δ是时变的，即δ(t)，且J和P的变化依赖于帆板转角δ(t)。故可用J(t)和Pi(t)取代J和P。T为由反作用飞轮和扰动输入产生的三轴控制力矩输入。

简化上述动力学方程，得到：

(6)

(7)

其中，

其中，u∈R3为三轴控制力矩输入，w∈R3为干扰力矩输入。

将上述动力学方程式写成状态空间方程形式：

(8)

(9)

(10)

量测方程可以写为：

y(t)=Cx(t)

(11)

(12)

其中，C∈Rm×n为输出矩阵；Φ为相应的模态矩阵。

2 递推模态参数辨识方法

2.1 系统描述和假设

动态系统进行离散化并考虑噪声对系统的影响，可以进一步将其写为如下的新息形式：

xk+1=Akxk+Bkuk+Kkek

(13)

yk=Ckxk+ek

(14)

其中，xk∈Rn×1为状态向量；uk∈Rr×1为输入向量；yk∈Rm×1为输出向量；下标k为离散后的时间点，ek表示新息白噪声序列，而矩阵Kk表示卡尔曼增益矩阵；A∈Rn×n为状态矩阵；B∈Rn×r为输入矩阵；C∈Rm×n为输出矩阵。那么对式(13)和(14)做进一步的变换，可以写为：

(15)

yk=Ckxk+ek

(16)

定义一个如下形式的VARX预测器，则k时刻的输出yk可以表述为

(17)

yk=Ξkφk

(18)

式中：

那么利用自适应滤波技术，可以得到各个采样时刻的矩阵Ξk的最小二乘递推形式为：

(19)

(20)

式中：β1为遗忘因子，需要保证0<β1≤1；Ζk为定义的临时变量矩阵，通常初值给为Z0=I。通过式(19)的递推计算，可以得到各个时刻的Markov参数矩阵Ξk的值，接下来利用该矩阵递推估计系统的状态量xk。

2.2 基于仿射投影算法的系统状态量估计方法

对于式(15)、(16)，如果暂不考虑量测噪声ek，那么从k-p到k时刻的输出信号{yk-p,yk-p+1,…,yk}可以分别写为：

yk-p=Ck-pxk-p

(21)

(22)

(23)

则式(21)～(23)可以用矩阵形式表示为：

(24)

式中：

τk-p=Γk-pxk-p

(25)

那么对于k-1时刻，式(25)可以表示为如下的最小二乘形式：

τk-1=Γk-1xk-1

(26)

其中

(27)

对于式(26)，这里引入APA算法，从而得到各个时刻的状态量xk的递推形式为：

(28)

Φk-1=τk-1-Γk-1xk-1

(29)

式中：μ为递推过程中的仿射投影因子，一般应满足0<μ≤1。在式(28)的递推迭代中，矩阵Γk-1的初值Γ0应保证为满秩矩阵，如果没有其他先验知识，则Γ0一般可以选择为

(30)

在式(28)中，在递推得到xk后，矩阵Γk的更新同样可以通过如下的递推最小二乘算法得到：

(31)

(32)

Γk=Γk-1+(τk-Γk-1xk)Wk

(33)

式中：β2为遗忘因子，而矩阵Lk的初值选取方式和Γ0相似。根据式(28)～(33)，可以递推计算得到各个时刻的状态量xk。

3 基于辨识状态量的无模型控制方法

无模型控制方法是一种数据驱动控制方法，根据辨识得到的状态量xk可以利用无模型控制器对航天器进行控制，而不需要其他系统模型。考虑如下MIMO离散时间非线性系统

y(k+1)=

f(y(k),…,y(k-ny),u(k),…,u(k-nu))

(34)

uT(k),…,uT(k-Lu+1)]T

(35)

对形如式(34)的MIMO离散时间非线性系统，提出如下两个假设。

假设 3.1除有限时刻点外，fi(·),i=1,2，…,m，对各个变量的分量都存在连续偏导数。

(36)

且对任意时刻k，Φf,Ly,Lu(k)=[Φ1(k)…ΦLy+Lu(k)] 是有界的，其中Φi(k)∈Rm×m,i=1,…Ly+Lu。

考虑如下输入准则函数

J(u(k))=||y*(k+1)-y(k+1)||2+

λ||u(k)-u(k-1)||2

(37)

其中，λ>0是一个权重因子。

将式(36)代入准则函数(37)中，对u(k)求差分并令其等于零，得

Δu(k-i+1))

(38)

考虑到该控制算法中包含了矩阵求逆运算，在矩阵维数较大时运算量大，运算速度较慢，根据无模型自适应控制(Model Free Adaptive Control，MFAC)理论，可以得到针对离散时间MIMO非线性系统的MFAC控制方案如下

(39)

(40)

其中，i=1,2，…,m;j=1,2，…,m;ρ1,ρ2,…ρLy+Lu∈(0,1]；η∈(0,2]；λ>0，μ>0。

辨识与控制流程为：

(1)初始化

(2)辨识计算系统状态量xk

(5)由(40)计算控制输入u(k)

(6)回到第(2)步进行下一步循环

4 仿真验证

本文以文献[11]中航天器的模态参数为初始值进行仿真实验，本体初始姿态角为[0.2 0.3 0.3]rad，角速度为[0 0 0](rad/s)；太阳能帆板旋转角速度为7.278×10-5(rad/s)(即360°/天)，四个附件分别选取的模态阶数为8、8、7、7。对姿态角和四个附件的一阶模态坐标的辨识结果如图2所示。