CN103914628A - 一种空间遥操作系统输出状态预测方法 - Google Patents

Abstract

本发明提出一种空间遥操作系统输出状态预测方法，首先建立空间遥操作系统的空间状态方程，其次根据建立的空间遥操作系统的空间状态方程，建立辅助模型；随后进行辅助模型状态估计，得到空间遥操作系统的空间状态方程在任意最小采样周期时刻的系统输出预测状态值，最后循环进行估计，直到完成全部状态输出预测过程。采用本发明的方法，能够保证在反馈信号频率较低的情况下，根据历史信息较为准确的预测系统输出状态。

Description

一种空间遥操作系统输出状态预测方法

技术领域

本发明属于系统状态预测控制领域，具体涉及一种空间遥操作系统输出状态预测方法。

背景技术

卡尔曼预测目前已经广泛的应用于航空/航天器的惯性导航系统、组合导航系统、目标预测跟踪和估计系统等，目前的研究的对象一般是把研究对象系统看为单率系统，这样不仅限制了方法的使用范围，也增加了很多的冗余设备，然而现实应用中存在一些被动的反馈信号，或者传感器和处理器等都具有不同的采样频率而造成的输入输出的采样频率不同步的情况，例如空间机器人的控制信号和关节反馈信号具有不同步的情况的，又如在遥操作控制中，地面的输入控制信号往往是连续的，相比于反馈的信号频率也要高很多。这种输入输出频率不同步的系统，我们称之为双率系统。目前关于双率系统采用的模型方法主要有以下几点：

1、基于提升技术的状态空间模型

2、基于多项式变换的差分方程模型

3、基于小波多分辨率分析的小波模型

以上的方法各有利弊，以提升技术为例，提升技术通过构造一个框架周期的提升系统，从而将双率系统转化为一个多变量的单率系统，然后再利用利用基本卡尔曼预测对状态值进行预测，但是该方法会造成计算量的增加，预测状态变量缺失的问题，若在遥操作过程中要求对系统的输出状态进行实时预测，该算法就不能也不能满足预测的要求，为此需要找到一个新的方法对双率系统输出状态值进行预测。

发明内容

技术方案

本发明的技术方案为：

所述一种空间遥操作系统输出状态预测方法，其特征在于：采用以下步骤：

步骤1：建立空间遥操作系统的空间状态方程：

$\left\{\begin{array}{c}x({\mathrm{KT}}_1+T_1)=\mathrm{Ax}\left({\mathrm{KT}}_1\right)+\mathrm{Bu}\left({\mathrm{KT}}_1\right)+w\left({\mathrm{KT}}_1\right)\\ y\left({\mathrm{KT}}_2\right)=\mathrm{Cx}\left(K{T}_2\right)+\mathrm{Du}\left({\mathrm{KT}}_2\right)+v\left({\mathrm{KT}}_2\right)\end{array}\right.$

其中K表示采样点，T₁为KT₁时刻空间遥操作系统的控制输入信号u(KT₁)的采样周期，x(KT₁)为KT₁时刻空间遥操作系统的系统状态值，A,B,C,D分别为状态方程参数；T₂为KT₂时刻空间遥操作系统反馈信号y(KT₂)的采样周期；w(KT₁)和v(KT₂)分别为过程噪声和测量噪声，都是平稳的零均值白噪声序列；

步骤2：根据步骤1建立的空间遥操作系统的空间状态方程，建立辅助模型：

$\left\{\begin{array}{c}x(\mathrm{kh}+h)=\mathrm{ax}\left(\mathrm{kh}\right)\mathrm{bu}\left(\mathrm{kh}\right)+r\left(\mathrm{kh}\right)\\ y\left(\mathrm{kh}\right)=\mathrm{cx}\left(\mathrm{kh}\right)+\mathrm{du}\left(\mathrm{kh}\right)+o\left(\mathrm{kh}\right)\end{array}\right.$

其中h为T₁和T₂的最小采样周期，满足条件T₁＝ph，T₂＝qh且p和q两个数互质，其中在同一采样周期内空间遥操作系统的控制输入信号保持不变；k表示辅助系统的采样点，x(kh)表示kh时刻的辅助系统中间状态值，u(kh)表示kh时刻的辅助系统输入值，y(kh)表示kh时刻的辅助系统目标状态值，a,b,c,d分别为辅助系统状态方程参数，o(kh)和v(kh)为辅助系统的过程噪声和测量噪声，都是平稳的零均值白噪声序列；辅助系统状态方程参数为：

$a=A^{\frac{1}{p}},b=B/\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{A}^{\frac{i}{p}},c=C,d=D,r\left(\mathrm{kh}\right)=w\left(\mathrm{kh}\right)/\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{A}^{\frac{i}{p}}$

步骤3：辅助模型状态估计：计算辅助模型的残差值：

根据公式

$o\left(\mathrm{kh}\right)=\frac{y\left(\mathrm{kh}\right)-{\mathrm{ca}}^{q}x\left(\mathrm{kh}\right)-(\mathrm{cb}\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}{a}^i+d)u\left(\mathrm{kh}\right)}{\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}{a}^i\left(\mathrm{cL}1(\mathrm{kh}+I)\right)}$

计算o(kh)，其中L1(kh)＝Kc，然后根据公式

$\mathrm{Ro}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}o(\mathrm{kh}+\mathrm{ih})o^T(\mathrm{kh}+\mathrm{ih})$

计算Ro，其中S表示残差序列长度；取误差值且为kh时刻x(kh)的估值，误差方差矩阵P(kh)和误差协方差矩阵M(kh)分别为P(kh)＝x(kh)x(kh)，M(kh)＝x(kh)x(kh-h)，计算kh+h时刻辅助系统中间状态的估计值x(kh+h)：

$\hat{x}(\mathrm{kh}+h)=a\hat{x}\left(\mathrm{kh}\right)+\mathrm{bu}\left(\mathrm{kh}\right)+L\left(\mathrm{kh}\right)[y\left(\mathrm{kh}\right)-y(\mathrm{kh}-h)]+L1\left(\mathrm{kh}\right)o\left(\mathrm{kh}\right)$

其中L(kh)、P(kh)、M(kh)和Rr的表达式为：

L(kh)＝[ac(M(kh)-P(kh))+L1(kh)Ro][2Ro+2c²P(kh)-2c²M(kh)]^-1

P(kh+h)＝a²P(kh)+Rr+L(kh)[ac(P(kh)-M(kh))+L1(kh)Ro]

+L1(kh)(L1(kh)+2I)^TRo

M(kh+h)＝aP(kh)+Rr+L(kh)(P(kh)-M(kh))c

$\mathrm{Rr}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}r(\mathrm{kh}+\mathrm{ih})r^T(\mathrm{kh}+\mathrm{ih})$

步骤4：估计步骤1中空间遥操作系统的空间状态方程在任意最小采样周期时刻T＝mh的系统输出预测状态值m为任意正整数：

$\hat{y}\left(\mathrm{kh}\right)=c\hat{x}\left(\mathrm{kh}\right)+\mathrm{du}\left(\mathrm{kh}\right)+o\left(\mathrm{kh}\right)$

步骤5：循环步骤3和步骤4直到完成全部状态输出预测过程。

有益效果

采用本发明的方法，能够保证在反馈信号频率较低的情况下，根据历史信息较为准确的预测系统输出状态。

附图说明

图1针对双率采样系统的目标状态预测方法的流程图

图2辅助模型的卡尔曼预测状态估计更新流程图

图3目标实际状态值与预测值的比较结果

具体实施方式

下面结合具体实施例描述本发明：

本实施例中的一种空间遥操作系统输出状态预测方法，采用以下步骤，其流程图如图1所示：

步骤1：对空间遥操作系统的回路过程建模，KT₁时刻系统的控制输入信号u(KT₁)频率较高，其采样周期为T₁，y(KT₂)为KT₂时刻遥操作系统的反馈信号，这种信号可以是力反馈或者其他反馈信号，且其采样周期为T₂，将遥操作过程简化为一个双率系统，其空间状态方程形式如下：

$\left\{\begin{array}{c}x({\mathrm{KT}}_1+T_1)=\mathrm{Ax}\left({\mathrm{KT}}_1\right)+\mathrm{Bu}\left({\mathrm{KT}}_1\right)+w\left({\mathrm{KT}}_1\right)\\ y\left({\mathrm{KT}}_2\right)=\mathrm{Cx}\left(K{T}_2\right)+\mathrm{Du}\left({\mathrm{KT}}_2\right)+v\left({\mathrm{KT}}_2\right)\end{array}\right.$

其中K表示采样点，x(KT₁)为KT₁时刻空间遥操作系统的系统状态值，A,B,C,D分别为状态方程参数；w(KT₁)和v(KT₂)分别为过程噪声和测量噪声，都是平稳的零均值白噪声序列；满足E(w(i))＝0，E(v(i))＝0cov(w(i),w(j)^T)＝Rwδ_ij，cov(v(i),v(j)^T)＝Rvδ_ij，其中δ_ij为kronecker-δ函数。

步骤2：根据步骤1建立的空间遥操作系统的空间状态方程，建立辅助模型：

$\left\{\begin{array}{c}x(\mathrm{kh}+h)=\mathrm{ax}\left(\mathrm{kh}\right)\mathrm{bu}\left(\mathrm{kh}\right)+r\left(\mathrm{kh}\right)\\ y\left(\mathrm{kh}\right)=\mathrm{cx}\left(\mathrm{kh}\right)+\mathrm{du}\left(\mathrm{kh}\right)+o\left(\mathrm{kh}\right)\end{array}\right.$

其中h为T₁和T₂的最小采样周期，满足条件T₁＝ph，T₂＝qh且p和q两个数互质，其中在同一采样周期内空间遥操作系统的控制输入信号保持不变；k表示辅助系统的采样点，x(kh)表示kh时刻的辅助系统中间状态值，u(kh)表示kh时刻的辅助系统输入值，y(kh)表示kh时刻的辅助系统目标状态值，a,b,c,d分别为辅助系统状态方程参数，o(kh)和v(kh)为辅助系统的过程噪声和测量噪声，都是平稳的零均值白噪声序列；辅助系统状态方程参数为：

$a=A^{\frac{1}{p}},b=B/\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{A}^{\frac{i}{p}},c=C,d=D,r\left(\mathrm{kh}\right)=w\left(\mathrm{kh}\right)/\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{A}^{\frac{i}{p}}$

步骤3：辅助模型状态估计：计算辅助模型的残差值：

根据公式

$o\left(\mathrm{kh}\right)=\frac{y\left(\mathrm{kh}\right)-{\mathrm{ca}}^{q}x\left(\mathrm{kh}\right)-(\mathrm{cb}\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}{a}^i+d)u\left(\mathrm{kh}\right)}{\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}{a}^i\left(\mathrm{cL}1(\mathrm{kh}+I)\right)}$

计算o(kh)，其中L1(kh)＝Kc，然后根据公式

$\mathrm{Ro}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}o(\mathrm{kh}+\mathrm{ih})o^T(\mathrm{kh}+\mathrm{ih})$

计算Ro，其中S表示残差序列长度；取误差值且为kh时刻x(kh)的估值，误差方差矩阵P(kh)和误差协方差矩阵M(kh)分别为P(kh)＝x(kh)x(kh)，M(kh)＝x(kh)x(kh-h)，计算kh+h时刻辅助系统中间状态的估计值x(kh+h)：

$\hat{x}(\mathrm{kh}+h)=a\hat{x}\left(\mathrm{kh}\right)+\mathrm{bu}\left(\mathrm{kh}\right)+L\left(\mathrm{kh}\right)[y\left(\mathrm{kh}\right)-y(\mathrm{kh}-h)]+L1\left(\mathrm{kh}\right)o\left(\mathrm{kh}\right)$

其中L(kh)、P(kh)、M(kh)和Rr的表达式为：

L(kh)＝[ac(M(kh)-P(kh))+L1(kh)Ro][2Ro+2c²P(kh)-2c²M(kh)]^-1

P(kh+h)＝a²P(kh)+Rr+L(kh)[ac(P(kh)-M(kh))+L1(kh)Ro]

+L1(kh)(L1(kh)+2I)^TRo

M(kh+h)＝aP(kh)+Rr+L(kh)(P(kh)-M(kh))c

$\mathrm{Rr}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}r(\mathrm{kh}+\mathrm{ih})r^T(\mathrm{kh}+\mathrm{ih})$

步骤4：估计步骤1中空间遥操作系统的空间状态方程在任意最小采样周期时刻T＝mh的系统输出预测状态值m为任意正整数：

$\hat{y}\left(\mathrm{kh}\right)=c\hat{x}\left(\mathrm{kh}\right)+\mathrm{du}\left(\mathrm{kh}\right)+o\left(\mathrm{kh}\right)$

步骤5：循环步骤3和步骤4直到完成全部状态输出预测过程。

本实施例中的考虑双率采样导航跟踪系统双率采样模型：

$\left\{\begin{array}{c}x(k+1)=0.8x\left(k\right)+1.5u\left(k\right)+w\left(k\right)\\ y\left(k\right)=1.3x\left(k\right)-0.7u\left(k\right)+v\left(k\right)\end{array}\right.$

其中，输入采样周期T₁＝3h，输出采样周期为T₂＝4h，h＝1s，框架周期pqh＝12h，系统噪声和测量噪声分别为Rw＝0.5²，Rv＝0.25²，仿真结果如图3所示，其中虚线表示实际的系统输出结果，实线表示估计的输出结果，从图3可以看出，输入时刻对应的实际目标状态值(y)与利用差分卡尔曼预测估计值(ye)较为相近，可以计算平均误差值为4.74％，说明该方法能很好的跟踪和预测输出状态y(k)的变化。

Claims (1)

1.一种空间遥操作系统输出状态预测方法，其特征在于：采用以下步骤：

步骤1：建立空间遥操作系统的空间状态方程：

$\left\{\begin{array}{c}x({\mathrm{KT}}_1+T_1)=\mathrm{Ax}\left({\mathrm{KT}}_1\right)+\mathrm{Bu}\left({\mathrm{KT}}_1\right)+w\left({\mathrm{KT}}_1\right)\\ y\left({\mathrm{KT}}_2\right)=\mathrm{Cx}\left(K{T}_2\right)+\mathrm{Du}\left({\mathrm{KT}}_2\right)+v\left({\mathrm{KT}}_2\right)\end{array}\right.$

其中K表示采样点，T₁为KT₁时刻空间遥操作系统的控制输入信号u(KT₁)的采样周期，x(KT₁)为KT₁时刻空间遥操作系统的系统状态值，A,B,C,D分别为状态方程参数；T₂为KT₂时刻空间遥操作系统反馈信号y(KT₂)的采样周期；w(KT₁)和v(KT₂)分别为过程噪声和测量噪声，都是平稳的零均值白噪声序列；

步骤2：根据步骤1建立的空间遥操作系统的空间状态方程，建立辅助模型：

$\left\{\begin{array}{c}x(\mathrm{kh}+h)=\mathrm{ax}\left(\mathrm{kh}\right)\mathrm{bu}\left(\mathrm{kh}\right)+r\left(\mathrm{kh}\right)\\ y\left(\mathrm{kh}\right)=\mathrm{cx}\left(\mathrm{kh}\right)+\mathrm{du}\left(\mathrm{kh}\right)+o\left(\mathrm{kh}\right)\end{array}\right.$

其中h为T₁和T₂的最小采样周期，满足条件T₁＝ph，T₂＝qh且p和q两个数互质，其中在同一采样周期内空间遥操作系统的控制输入信号保持不变；k表示辅助系统的采样点，x(kh)表示kh时刻的辅助系统中间状态值，u(kh)表示kh时刻的辅助系统输入值，y(kh)表示kh时刻的辅助系统目标状态值，a,b,c,d分别为辅助系统状态方程参数，o(kh)和v(kh)为辅助系统的过程噪声和测量噪声，都是平稳的零均值白噪声序列；辅助系统状态方程参数为：

$a=A^{\frac{1}{p}},b=B/\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{A}^{\frac{i}{p}},c=C,d=D,r\left(\mathrm{kh}\right)=w\left(\mathrm{kh}\right)/\underset{i=0}{\overset{p-1}{\mathrm{\Σ}}}{A}^{\frac{i}{p}}$

步骤3：辅助模型状态估计：计算辅助模型的残差值：

根据公式

$o\left(\mathrm{kh}\right)=\frac{y\left(\mathrm{kh}\right)-{\mathrm{ca}}^{q}x\left(\mathrm{kh}\right)-(\mathrm{cb}\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}{a}^i+d)u\left(\mathrm{kh}\right)}{\underset{i=0}{\overset{q-1}{\mathrm{\Σ}}}{a}^i\left(\mathrm{cL}1(\mathrm{kh}+I)\right)}$

计算o(kh)，其中L1(kh)＝Kc，然后根据公式

$\mathrm{Ro}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}o(\mathrm{kh}+\mathrm{ih})o^T(\mathrm{kh}+\mathrm{ih})$

计算Ro，其中S表示残差序列长度；取误差值且为kh时刻x(kh)的估值，误差方差矩阵P(kh)和误差协方差矩阵M(kh)分别为P(kh)＝x(kh)x(kh)，M(kh)＝x(kh)x(kh-h)，计算kh+h时刻辅助系统中间状态的估计值x(kh+h)：

$\hat{x}(\mathrm{kh}+h)=a\hat{x}\left(\mathrm{kh}\right)+\mathrm{bu}\left(\mathrm{kh}\right)+L\left(\mathrm{kh}\right)[y\left(\mathrm{kh}\right)-y(\mathrm{kh}-h)]+L1\left(\mathrm{kh}\right)o\left(\mathrm{kh}\right)$

其中L(kh)、P(kh)、M(kh)和Rr的表达式为：

L(kh)＝[ac(M(kh)-P(kh))+L1(kh)Ro][2Ro+2c²P(kh)-2c²M(kh)]^-1

P(kh+h)＝a²P(kh)+Rr+L(kh)[ac(P(kh)-M(kh))+L1(kh)Ro]

+L1(kh)(L1(kh)+2I)^TRo

M(kh+h)＝aP(kh)+Rr+L(kh)(P(kh)-M(kh))c

$\mathrm{Rr}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}r(\mathrm{kh}+\mathrm{ih})r^T(\mathrm{kh}+\mathrm{ih})$

步骤4：估计步骤1中空间遥操作系统的空间状态方程在任意最小采样周期时刻T＝mh的系统输出预测状态值m为任意正整数：

$\hat{y}\left(\mathrm{kh}\right)=c\hat{x}\left(\mathrm{kh}\right)+\mathrm{du}\left(\mathrm{kh}\right)+o\left(\mathrm{kh}\right)$

步骤5：循环步骤3和步骤4直到完成全部状态输出预测过程。