CN105425589A - 提高航天器惯性参数辨识精度的输入信号设计方法 - Google Patents

Abstract

本发明提供了一种提高航天器惯性参数辨识精度的输入信号设计方法，该方法针对航天器惯性参数可利用回归矩阵精确线性化的特点，提出一种基于直接配点法的最优输入设计方法。

Description

提高航天器惯性参数辨识精度的输入信号设计方法

技术领域

本发明涉及系统辨识领域，具体的涉及一种提高航天器惯性参数辨识精度的输入信号设计方法。

背景技术

系统参数辨识最主要的两个步骤分别是输入信号设计和输出信号分析处理，其中输入信号设计至关重要，因为如果采用不恰当的输入信号，即使采用最先进的分析处理算法，也不能从输出信号中获得准确的参数估计值。最优输入设计(OptimalInputDesign，OID)，是指输入信号在满足具体约束的条件下，设计使参数辨识误差最小的输入信号。最早的OID思想可追溯到1935年Fisher等人的研究中，但人们还是公认Levin是首次系统地提出OID概念的学者。

OID研究的一个重要问题是:如何提高连续时间动态系统的参数估计精度。这类系统以含有待估计参数的常微分方程(OrdinaryDifferentialEquation)或微分代数方程(DifferentialAlgebraicEquation)作为状态方程。早在20世纪70年代，解决该问题的基本思想就已成形，主要包括两点：一是采用Fisher信息矩阵(FisherInformationMatrix，FIM)作为优化依据；二是采用最优控制理论的方法寻找OID的解，从而获得高精度的参数估计值。时至今日，许多OID的实际应用依然遵循这一思想。

Morelli研究了飞行试验中飞机舵面运动设计，给出了一种状态有界约束下的设计方法，该方法以FIM逆矩阵的迹作为最优性能指标(即A-准则)，采用动态规划方法求解最优控制问题，获得分段连续的输入信号。Jauberthie在Morelli的基础上用连续可微函数对舵面运动指令进行近似，增加了梯度法优化步骤，获得了连续的最优输入信号，更适于实际情况下舵机的执行。HanYongsu研究了飞机多操纵面同时激励的情况，特点是各通道的输入信号由多个正弦基函数叠加而成，通过调整正弦函数的频率使各通道输入信号相互正交，从而实现观测数据与输入信号完全去相关。但是由于FIM含有参数灵敏度函数，复杂系统的灵敏度函数难以获得。限制了该方法在航天器领域的运用。

除了以FIM作为优化依据外，对于一类参数可线性化的系统，还可以用参数回归矩阵(RegressionMatrix)作为优化依据，现有技术中的这类方法常用于机器人参数辨识的研究，通常采用傅立叶级数构造输入信号，以傅立叶系数作为优化设计变量，以回归矩阵的条件数作为具体的代价函数，利用启发式优化算法(如遗传算法或粒子群算法等)进行求解。由于回归矩阵是状态轨迹和输入轨迹的函数，每次计算代价函数都要进行轨迹积分，因此该类方法的优化效率较低。

发明内容

本发明的目的在于提供一种提高航天器惯性参数辨识精度的最优输入设计方法，该发明解决了现有技术中已有的OID方法运用于航天器惯性参数辨识时，输入信号优化效率低，待估参数灵敏度函数难以获得的技术问题。

本发明提供一种提高航天器惯性参数辨识精度的输入信号设计方法，包括以下步骤：

步骤S100：选取在轨运行航天器中的未知惯性参数作为待估参数，将含待估参数的动力学微分方程以线性回归方程形式表示；

步骤S200：根据线性回归方程中的回归矩阵构造最小二乘法矩阵，通过步骤S300～步骤S400获得能够使最小二乘法矩阵的条件数接近1的输入信号；

步骤S300：以最小二乘法矩阵中各不相同的独立元素作为增广状态变量，扩充到航天器动力学微分方程中，得到增广状态微分方程；

步骤S400：将最优输入设计问题归纳为最优控制问题，并对最优控制问题进行求解，得到用于参数辨识实验所需的输入信号。

进一步地，步骤S400中最优控制问题的解决方法包括以下步骤：

步骤S401：以增广状态微分方程作为描述航天器惯性系统运动的状态方程；

步骤S402：根据航天器的姿控执行机构的实际物理局限，确定各通道输入信号的允许取值范围；

步骤S403：设定最小二乘法矩阵独立元素的初值为零，结合航天器角速度的真实初始状态，共同构成增广状态微分方程的初始状态，对增广状态微分方程的终端状态不设约束；

步骤S404：用静态的参数对动态的输入信号进行描述；

设所有通道的输入信号均为不同频率的正弦基函数之和，构成各通道输入信号的正弦基函数具有独立的相位、频率和幅值，设第i通道输入信号的正弦基函数的频率系数取值集合为：

S_i＝{(i+1),(i+1)+m,(i+1)+2m,(i+1)+3m,...}

其中，m为总通道数，这样使得任意两通道的输入信号相互正交；

设构成第i通道输入信号的第k个正弦基函数的幅值和相位分别为A_k和ψ_k，这样动态的输入信号就可以用静态的参数{A_k,ψ_k}进行描述，比如第i个通道的输入信号被描述为：

$U_i\left(t\right)=\underset{k\∈S_i}{\Σ}{A}_{k}sin(\frac{2\mathrm{\π}kt}{T}+{\mathrm{\ψ}}_k),i=1,2,...,m$

步骤405：采用直接配点法，对动态的增广状态轨迹进行参数化，将动态的最优输入设计问题转化为静态的非线性规划问题；

步骤S406：采用串行优化策略，对静态非线性规划问题进行求解，获得优化后的描述输入信号的静态参数；

步骤S407：使用优化后的输入信号静态参数，重构得到动态的最优输入信号。

进一步地，步骤S405为：将输入持续时间区间进行分段，每段的两个端点称为节点，记节点处对应的状态变量为{X₁,X₂,...,X_N}，子区间[t_i,t_i+1]上的状态变量X(t)用三次Hermite插值多项式近似表示。

进一步地，步骤S404中还包括设置幅值和相位满足等式使得初始时刻t₀和终端时刻t_f输入信号值为0。

进一步地，约束条件还包括状态轨迹的边界条件约束和不等式约束，不等式约束由实际系统的物理约束确定。

本发明的技术效果：

1、本发明提供方法的优化依据有别于传统的FIM依据和回归矩阵依据，采用最小二乘法矩阵作为优化依据，性能指标函数只与终端时刻的最小二乘法矩阵有关，减轻了现有OID方法性能指标函数的计算复杂度。

2、本发明提供的方法采用正弦函数作为参数化输入信号的基函数，通过设置幅值和相位约束确保初/末时刻输入取值为零，从而使终端时刻状态轨迹回到初始位置附近，解决了物理实验中惯性航天器在实验结束时，难以及时恢复原始运动状态的问题，为系统辨识物理实验提供了便利。

3、本发明提供的方法采用直接配点法离散状态变量，设定了的串行优化策略中各步骤的操作内容，不仅保证了静态非线性规划问题的解为原动态输入信号优化问题的可行解，从而极大地提高了输入信号的优化效率。

4、本发明提供的方法已通过受控刚体航天器惯量矩辨识算例的证明，该方法能够快速获得优化的输入信号，使用所述的优化输入信号对待测系统进行激励，可以获得有用信息最多的观测数据，使用所述观测数据进行辨识计算，能够加快参数估计值的收敛速度，并极大的提高参数估计的精度。

具体请参考根据本发明的提高航天器惯性参数辨识精度的输入信号设计方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。

附图说明

图1是本发明优选实施例的最优输入求解策略流程示意图；

图2是本发明优选实施例提供的具有OID环节的参数估计流程示意图；

图3是本发明优选实施例算例中优化前后的输入信号对比示意图；

图4是本发明实施例算例中优化前后的输入信号在这两种输入信号激励下产生的航天器角速度轨迹曲线示意图；

图5是本发明优选实施例算例中作为对比的传统3211型输入信号和doublet型输入信号的变化曲线示意图；

图6是本发明优选实施例算例中在对比的传统3211型输入信号和doublet型输入信号激励下产生的航天器角速度轨迹曲线示意图；

图7是本发明优选实施例算例中分别在优化输入信号、未优化输入信号、3211型输入信号和doublet输入信号作用下进行参数辨识计算所得的待估参数估计值变化曲线示意图；

图8是本发明优选实施例算例中分别在优化输入信号、未优化输入信号、3211型输入信号和doublet输入信号作用下进行参数辨识计算所得的待估参数总体估计残差变化曲线示意图；

图9是本发明优选实施例提高航天器惯性参数辨识精度的最优输入设计方法流程示意图。

具体实施方式

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

本发明针对现有方法的缺点，研究一类参数可线性化系统辨识的最优输入问题，首先对OID问题的一般数学模型进行描述，并给出了限定问题性质的基本假设；然后提出了一种新型优化准则，该准则可以作为待估参数总体可辨识度的度量指标；之后提出输入信号的参数化新方法，这种输入信号的参数化方法可以使多维输入信号相互正交，并通过增加约束条件使输入信号在初始时刻和终端时刻为零；紧接着，基于直接配点法提出了一种求解OID问题的计算方法，该方法能以较少的时间获得静态的输入信号参数的优化解，使用该优化后的静态输入参数，可以重构出辨识实验所需的动态的最优输入信号；最后，将所提出的整套设计流程，应用于一个简单的空间受控刚体惯量矩参数辨识问题，给出了最优输入设计结果，通过与其他类型输入信号的辨识结果对比，证明采用本发明提供方法能够明显提高参数估计的效率和精度。

本文中所用到的数值积分器或非线性规划求解器均可以为市售产品。

参见图9，本发明提供的提高航天器惯性参数辨识精度的输入信号设计方法包括以下步骤：

步骤S100：选取在轨运行航天器中的未知惯性参数作为待估参数，将含待估参数的动力学微分方程以线性回归方程形式表示；

步骤S200：根据线性回归方程中的回归矩阵构造最小二乘法矩阵，通过步骤S300～步骤S400获得能够使最小二乘法矩阵的条件数接近1的输入信号；

步骤S300：以最小二乘法矩阵中各不相同的独立元素作为增广状态变量，扩充到航天器动力学微分方程中，得到增广状态微分方程；

步骤S400：将最优输入设计问题归纳为最优控制问题，并对最优控制问题进行求解，得到用于参数辨识实验所需的输入信号。

此处所用参数辨识算法可以为任何现有的参数辨识算法均可。对于最优控制问题的求解也可以用现有的控制问题进行求解。

优选的，为了避免使用变分法求解，最优控制问题的解决包括以下步骤：

步骤S401：以增广状态微分方程作为描述航天器惯性系统运动的状态方程；

步骤S402：根据航天器的姿控执行机构的实际物理局限，确定各通道输入信号的允许取值范围；

步骤S403：设定最小二乘法矩阵独立元素的初值为零，结合航天器角速度的真实初始状态，共同构成增广状态微分方程的初始状态，对增广状态微分方程的终端状态不设约束；

步骤S404：用静态的参数对动态的输入信号进行描述。

设所有通道的输入信号均为不同频率的正弦基函数之和，构成各通道输入信号的正弦基函数具有独立的相位、频率和幅值，设第i通道输入信号的正弦基函数的频率系数取值集合为：

S_i＝{(i+1),(i+1)+m,(i+1)+2m,(i+1)+3m,...}

其中，m为总通道数，这样使得任意两通道的输入信号相互正交；

设构成第i通道输入信号的第k个正弦基函数的幅值和相位分别为A_k和ψ_k，这样动态的输入信号就可以用静态的参数{A_k,ψ_k}进行描述，比如第i个通道的输入信号就被描述为

$U_i\left(t\right)=\underset{k\∈S_i}{\Σ}{A}_{k}sin(\frac{2\mathrm{\π}kt}{T}+{\mathrm{\ψ}}_k),i=1,2,...,m$

至此，输入变量可以用独立的静态设计变量{A_k,ψ_k}参数化表示；

优选的，设置幅值和相位满足等式使得初始时刻t₀和终端时刻t_f输入信号值为0。在该等式约束的作用下，系统在激励结束后将回到初始状态，这种性质在进行实际物理试验时是非常有利。输入信号初始时刻和终端时刻取值均为零，也有利于物理执行机构准确复现该曲线。

步骤405：采用直接配点法，对动态的增广状态轨迹进行参数化，将动态的最优输入设计问题转化为静态的非线性规划问题。具体为，将输入持续时间区间进行分段，每段的两个端点称为节点，记节点处对应的状态变量为{X₁,X₂,...,X_N}，子区间[t_i,t_i+1]上的状态变量X(t)用三次Hermite插值多项式近似表示。

至此通过S404和405将最优控制问题转化为静态最优化问题，优化变量为输入信号参数和节点处的状态变量。最优控制问题多用变分法解决，但我们通过2个参数化输入信号参数，得到设计变量和节点处的状态变量，从而将最优控制问题转化为静态非线性规划问题，就可以用处理该问题的算法来。

步骤S406：采用串行优化策略，对静态非线性规划问题进行求解，获得优化后的描述输入信号的静态参数。具体为该静态非线性规划问题通过SQP算法解决，采用串行优化策略得到SQP算法的解。该解中即包含了优化输入信号的设计变量，即{A_k,ψ_k}。

步骤S407：使用优化后的输入信号设计变量{A_k,ψ_k}，可重构出动态的最优输入信号。

该方法的使用过程：给定待估参数先验值θ₀。该先验值为航天器系统出厂时设备自带的标称参数，在轨使用一段时间后，该数值多会发生改变，为了能精确掌握运行一段时间后航天器的运动状态，需要设定能获得该参数的输入信号，对航天器进行激励，并对激励过程中获得的数据进行参数辨识，从而获得当前情况下，该航天器的实时参数估计值。

根据该先验值设计系统的最优输入信号u(t)并对真实系统进行激励，同样以该先验值作为辨识算法(该算法为任一现有的参数辨识算法)的初始估值，与采集得到的测量数据y(t)一起代入参数辨识算法中，计算得到参数的估计值即在轨运行一段时间后，航天器的具体估计参数值。

最优性能指标函数为关于最小二乘法矩阵的单调递增函数。

本发明提供的方法具体如下：

1本发明所适用的系统限制

考虑如下非线性系统

$\left\{\begin{array}{c}\stackrel{\·}{x}=f\left(t,x,\mathrm{\θ},u\right)\\ y=h\left(t,x,\mathrm{\θ},u\right)+\mathrm{\ϵ}\left(t\right)\end{array}\right.---\left(1.1\right)$

其中，表示时间；表示状态变量；表示输出变量；为未知参数；表示输入变量；为观测噪声，服从高斯分布且均值为零。

假设1:设ε(t)≡0，系统可改写为关于参数的回归方程的形式

$\mathrm{\φ}(t,\stackrel{\·}{x},x,y,u)\·\mathrm{\η}\left(\mathrm{\θ}\right)=\mathrm{\ν}(t,\stackrel{\·}{x},x,y,u)---\left(1.2\right)$

式中，即回归矩阵函数，对时间t一阶连续可微；为未知参数的线性或非线性函数，其反函数存在且唯一；且t∈[0,+∞)上不恒等于零.

系统参数辨识的OID问题可以描述为:设计有限输入变量u(t)，使得在u(t)激励下根据系统响应y(t)能以最小估计误差获得未知参数θ的估计值即并满足变量有界约束

$\begin{array}{c}{\stackrel{\·}{x}}_{lb}\≤\stackrel{\·}{x}\≤{\stackrel{\·}{x}}_{ub},x_{lb}\≤x\≤x_{ub}\\ y_{lb}\≤y\≤y_{ub},u_{lb}\≤u\≤u_{ub}\end{array},x_{lb}\≤x\≤x_{ub}$

y_lb≤y≤y_ub，u_lb≤u≤u_ub

其中下标lb和ub分别表示变量的下界与上界。

2最优输入设计方法

2.1最优输入性能指标

将输入持续时间区间t∈[t₀,t_f]等间隔化分为N-1个子区间

t₀＝t₁＜t₂＜...＜t_N＝t_f(1.3)

区间宽度Δt＝t_j-t_i,i＜j，设

${\mathrm{\φ}}_i=\mathrm{\φ}(t_i,\stackrel{\·}{x}(t_i),x\left(t_i\right),y\left(t_i\right),u\left(t_i\right))$

$v_i=\mathrm{\ν}(t_i,\stackrel{\·}{x}(t_i),x\left(t_i\right),y\left(t_i\right),u\left(t_i\right))$

$v_i=\mathrm{\ν}(t_i,\stackrel{\·}{x}(t_i),x\left(t_i\right),y\left(t_i\right),u\left(t_i\right))$

根据最小二乘原理，未知参数η(θ_i)的估计为

$\widehat{\mathrm{\η}}={\left(\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\φ}}_k^T{\mathrm{\φ}}_k\right)}^{-1}\left(\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\φ}}_k^T{\mathrm{\ν}}_k\right)---\left(1.4\right)$

当Δt→0时，式(1.4)改写为

$\widehat{\mathrm{\η}}={\left({\∫}_{t_0}^{t_f}{\mathrm{\φ}}^T\left(t\right)\mathrm{\φ}\left(t\right)dt\right)}^{-1}\left({\∫}_{t_0}^{t_f}{\mathrm{\φ}}^T\left(t\right)\mathrm{\ν}\left(t\right)dt\right)---\left(1.5\right)$

设Φ的条件数为cond(Φ)∈[1,+∞)。根据矩阵分析原理可知，cond(Φ)越接近于1，估计值受矩阵摄动ΔΦ的影响越小，即参数估计抗噪声的能力越强，估计值也就更接近于真实值η。因此，可以将cond(Φ)作为衡量参数估计精度的依据。注意到Φ为正定对称矩阵，引入记号

${\stackrel{\·}{\mathrm{\Φ}}}_{i,j}={\mathrm{\φ}}^T\mathrm{\φ}(i,j),i\≤j,j=1,2,...,N---\left(1.6\right)$

其中，Φ_i,j和φ^Tφ(i,j)分别表示矩阵和φ^Tφ的第i行第j列元素。

将Φ_i,j视为新的状态，设扩维状态矢量为X＝[x^T,Φ₁₁,...,Φ_pp]^T，最优输入设计问题可以被归纳成如下最优控制问题

min

κ＝L(X(t_f),t_f)(1.7)a

s.t.

$\stackrel{\·}{X}=F\left(X\right(t),U(t),t)---\left(1.7b\right)$

μ(X(t₀),t₀,X(t_f),t_f)＝0(1.7c)

g(X(t),U(t),t)≤0(1.7d)

γ(X(t),U(t),t)＝0(1.7e)

式中，表示最优性能指标，泛函L(·)为关于cond(Φ)的单调递增函数，这里定义为

L(X(t_f),t_f)＝lg(cond[Φ(t_f)])

式(1.7b)为最优控制问题的状态方程约束，扩维状态方程右函数定义为F(·)＝[f^T(·),φ^Tφ(i,j)]^T,i≤j≤p；式(1.7c)为状态轨迹的边界条件约束，注意在最优输入问题中Φ_i,j(t₀)＝0；式(1.7d)为不等式约束，由各变量的上下界约束变形得到；式(1.7e)为等式约束，根据观测方程得到。

对于式(1.7)所示的最优控制问题，很难使用变分法求得解析解，通常将其转化为非线性规划问题进行数值寻优求解，进行数值求解的首先要考虑的就是输入变量的参数化，即将输入变量用一组静态的设计变量表示。

2.2具有多输入正交性质的输入变量参数化

采用以下方法对输入变量U(t)进行参数化，首先设每一维输入为不同频率的正弦基函数之和。考虑到输入变量是多维的，如果各个维度的信号相互正交，将使得参数估计完全与输入无关，有利于提高模型参数的辨识精度。两个信号和在区间[t₀,t_f]内正交的定义为

$<U_1,U_2>={\&Integral;}_{t_0}^{t_f}{U}_1\left(t\right)U_2\left(t\right)dt=0$

其中，＜U₁,U₂＞表示两个信号的内积。

在系统辨识过程中，输入信号的频带越宽，对系统模型不确定性的鲁棒性就越强，因而希望输入信号频带不仅能够覆盖系统的模态频率，而且具有较宽的频带。设计构成输入变量的正弦基函数具有独立的相位，频率和幅值，U(t)的第i维信号可以写为

$U_i\left(t\right)=\underset{k\&Element;S_i}{\&Sigma;}{A}_{k}sin(\frac{2\mathrm{\&pi;}kt}{T}+{\mathrm{\&psi;}}_k),i=1,2,...,m---\left(1.8\right)$

式中，S_i是构成U_i的正弦基函数频率系数的集合；m是输入信号的维数，即控制执行机构的个数；T＝t_f-t₀表示输入的持续时间；A_k是第k个正弦函数的幅值；ψ_k是第k个正弦函数的相位。

考虑具有形如式(1.8)的两个输入信号

$U_1\left(t\right)=sin(\frac{2{\mathrm{\&pi;k}}_{1}t}{T}+{\mathrm{\&psi;}}_1),U_2\left(t\right)=sin(\frac{2{\mathrm{\&pi;k}}_{2}t}{T}+{\mathrm{\&psi;}}_2),k_1\&NotEqual;k_2$

其内积为

$\begin{array}{c}<U_1\left(t\right),U_2\left(t\right)>\\ \&ap;\underset{i=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\sin \left(\frac{2{\mathrm{\&pi;k}}_1{t}_i}{T}+{\mathrm{\&psi;}}_1\right)\sin \left(\frac{2{\mathrm{\&pi;k}}_2{t}_i}{T}+{\mathrm{\&psi;}}_2\right)\\ =\underset{i=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\sin \left(\frac{2{\mathrm{\&pi;k}}_{1}i}{N}+{\mathrm{\&psi;}}_1\right)\sin \left(\frac{2{\mathrm{\&pi;k}}_{2}i}{N}+{\mathrm{\&psi;}}_2\right)\end{array}$

根据三角恒等式，有

$\begin{array}{c}<U_1\left(t\right),U_2\left(t\right)>\&ap;\frac{1}{2}\underset{i=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\cos \left(\frac{2\mathrm{\&pi;}\left(k_1-k_2\right)i}{N}+{\mathrm{\&psi;}}_1-{\mathrm{\&psi;}}_2\right)\\ -\frac{1}{2}\underset{i=0}{\overset{N-1}{\mathrm{\&Sigma;}}}cos\left(\frac{2\mathrm{\&pi;}\left(k_1+k_2\right)i}{N}+{\mathrm{\&psi;}}_1+{\mathrm{\&psi;}}_2\right)\end{array}$

注意，对于任意整数k和常值相角ψ上式才等于0。因此，输入信号中的各正弦基函数取整数k可以保证各维信号是正交的。为此，设计频率参数的集合为

$\begin{array}{c}S=\{S_1,S_2,...,S_m\}\\ S_1=\{\mathrm{2,2}+m,2+2m,...\}\\ S_2=\{\mathrm{3,3}+m,3+2m,...\}\\ .\\ .\\ .\\ S_m=\{(m+1),(m+1)+m,(m+1)+2m,...\}\end{array}---\left(1.9\right)$

至此，输入变量可以用独立设计变量{A_k,ψ_k}参数化表示，再考虑到输入信号的工程可实现性，期望在t₀和t_f时刻输入变量为零，即控制参数要满足约束

$\underset{k\&Element;S_i}{\&Sigma;}{A}_k{\mathrm{sin\&psi;}}_k=0---\left(1.10\right)$

2.3基于直接配点法的状态变量参数化

控制变量参数化后，可以用显示积分处理微分等式约束(1.7b)，即所谓直接打靶法求解，但该方法需要在优化的每步迭代中数值积分状态轨迹，计算效率较低，因此本发明引入直接配点法思想，对状态轨迹也进行参数化。如式(1.3)所示，将输入持续时间历程分段，每段的两个端点称为节点(node)记节点处对应的状态变量和控制变量为{X₁,X₂,...,X_N}，对时间分段不一定是等分，但为简化问题，这里采用等分方式。子区间[t_i,t_i+1]上的状态变量X(t)用三次Hermite插值多项式近似表示状态，即

X(s)＝c₀+c₁s+c₂s²+c₃s³

T_i＝t_i+1-t_i,i＝1,2,...,N(1.11)

s＝(t-t_i)/T_is∈[0,1]

状态量在子区间上应满足边界条件

X(0)＝X_i

$\frac{\mathrm{dX}\left(s\right)}{\mathrm{ds}}{|}_{s=0}={\stackrel{\&CenterDot;}{X}}_i=f(X_i,U_i,t_i)---\left(1.12\right)$

X(1)＝X_i+1

$\frac{dX\left(s\right)}{ds}{|}_{s=1}={\stackrel{\&CenterDot;}{X}}_{i+1}=f(X_{i+1},U_{i+1},t_{i+1})$

根据边界条件，可求得Hermite多项式中的系数为

$c_0=X_i,c_2=-3X_i-2{\stackrel{\&CenterDot;}{X}}_i+3X_{i+1}-{\stackrel{\&CenterDot;}{X}}_{i+1}$

$c_1={\stackrel{\&CenterDot;}{X}}_i,c_3=2X_i+{\stackrel{\&CenterDot;}{X}}_i-2X_{i+1}+{\stackrel{\&CenterDot;}{X}}_{i+1}$

从而给定一组节点处的状态变量，便可求解相应的表示状态的多项式。

子区间的中点即s＝0.5处为配点(collocation)，可得配点处的状态为

$X_{c,i}=\frac{X_i+X_{i+1}}{2}+\frac{T_i({\stackrel{\&CenterDot;}{X}}_{i+1}-{\stackrel{\&CenterDot;}{X}}_i)}{8}---\left(1.13\right)$

${\stackrel{\&CenterDot;}{X}}_{c,i}=\frac{3(X_i+X_{i+1})}{2T_i}-\frac{{\stackrel{\&CenterDot;}{X}}_i+{\stackrel{\&CenterDot;}{X}}_{i+1}}{4}$

为确保多项式能更好地拟合状态变量，应使配点处由状态方程计算得到的状态变量导数

${\stackrel{\&CenterDot;}{X}}_{EOM,i}=f(X_{c,i},U(t_{c,i}),t_{c,i})$

与多项式求得的状态变量导数相等，即期望向量

$d_i={\stackrel{\&CenterDot;}{X}}_{EOM,i}-{\stackrel{\&CenterDot;}{X}}_{c,i}---\left(1.14\right)$

趋于零，d_i称为Defect向量，d_i＝0将问题(1.7)的微分等式约束转变为非线性等式约束。

2.4优化策略

经过状态变量和输入变量的参数化，以节点处的状态变量{X₁,X₂,...,X_N}以及输入信号参数(i＝1,2,...,m,k_i∈S_i)作为设计变量

$z=(X_1^T,X_2^T,X_N^T,A_{1,k_1},A_{2,k_2},...,A_{m,k_m},{\mathrm{\&psi;}}_{1,k_1},{\mathrm{\&psi;}}_{2,k_2},...,{\mathrm{\&psi;}}_{m,k_m})---\left(1.15\right)$

动态优化问题(1.7)转化为静态非线性规划问题

min

κ＝L(X_N)(1.16)a

s.t.

μ(X₁,X_N)＝0(1.16)b

g(z)≤0(1.16)c

$\left\{\begin{array}{c}\mathrm{\&gamma;}\left(z\right)=0\\ d\left(z\right)=0\\ A_{k_i}{\mathrm{sin\&psi;}}_{k_i}=0\end{array}\right.---\left(1.16\right)d$

对(1.17)这类高维非线性规划问题，可以采用序列二次规划(SQP)算法进行求解，SQP算法启动需要设计变量初值猜测值(InitialGuess，IG)，为了确保优化结果是原动态优化问题(1.7)的解，引入串行优化策略：首先，随机生成输入变量参数初值，显式积分状态方程(1.7b)获得状态轨迹；然后，取较少的节点数，并以静态优化问题的等式约束为目标函数进行优化求解，获得满足等式约束的可行解；之后，根据计算得到的可行解，插值获得更多节点，以此为优化初值计算静态优化问题(1.17)，获得优化的控制参数；最后，将优化后的控制参数再代入状态方程(1.7b)显示积分，验证优化结果是否为最优控制的解，以及最优输入性能指标是否得到改善。整个串行优化过程如图1所示。

引入OID环节后，动态系统的未知参数估计流程如图2所示。给定待估参数先验值θ₀，根据该先验值设计系统的最优输入信号u(t)并对真实系统进行激励，同样以该先验值作为辨识算法的初始估值，与采集得到的响应一起代入辨识算法计算参数的估计值

3仿真算例

3.1算例模型

为了验证本发明提出的最优设计方法对参数估值精度的提高能力及其计算效率，考虑一个简单的空间受控刚体惯量矩估计问题，空间受控刚体动力学方程为

$J\stackrel{\&CenterDot;}{\mathrm{\&omega;}}+{\mathrm{\&omega;}}^{\&times;}J\mathrm{\&omega;}=\mathrm{\&tau;}---\left(1.17\right)$

式中，J＝diag(J₁,J₂,J₃)为未知惯量矩参数；ω＝[ω₁,ω₂,ω₃]^T为旋转角速度矢量；τ＝[τ₁,τ₂,τ₃]^T为输入力矩。该系统的状态变量和输入变量分别为ω和τ，同时，假设角速度可以用陀螺测量，则输出变量也为ω。将(1.17)改写为回归方程的形式，有

$\mathrm{\&phi;}(t,\stackrel{\&CenterDot;}{x},x)=\left[\begin{array}{ccc}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_1& -{\mathrm{\&omega;}}_2{\mathrm{\&omega;}}_3& {\mathrm{\&omega;}}_2{\mathrm{\&omega;}}_3\\ {\mathrm{\&omega;}}_1{\mathrm{\&omega;}}_3& {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_2& -{\mathrm{\&omega;}}_1{\mathrm{\&omega;}}_3\\ -{\mathrm{\&omega;}}_1{\mathrm{\&omega;}}_2& {\mathrm{\&omega;}}_1{\mathrm{\&omega;}}_2& {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_3\end{array}\right],\mathrm{\&eta;}=\left[\begin{array}{c}{J}_1\\ J_2\\ J_3\end{array}\right],\mathrm{\&nu;}\left(t\right)=\left[\begin{array}{c}{\mathrm{\&tau;}}_1\\ {\mathrm{\&tau;}}_2\\ {\mathrm{\&tau;}}_3\end{array}\right]$

待估计参数的真实值为η＝[200,400,300]^Tkg×m²，最优输入设计和参数辨识算法用到的的待估参数先验值为η₀＝[100,200,150]^Tkg×m.设输入持续时间T＝60s，输入变量τ_i∈[0,4]N，设计正弦基函数频率参数为：

τ₁:k＝2,5,8,12

τ₂:k＝3,6,9,13

τ₃:k＝4,7,10,14

3.2最优输入设计结果

将本发明提出的最优输入设计方法在MATLAB2012b环境中开发相应的计算程序，SQP算法使用MATLAB优化工具箱中fmincon函数中的算法，初始可行解生成阶段取节点数为10，优化解生成阶段取节点数为60。共进行10次优化计算，平均每次计算耗时435.96sec。表1中列出了未优化的输入信号设计变量和优化后的输入信号设计变量取值，表中最后一行为法矩阵条件数，表征了不同输入信号激励下系统待估参数的可辨识度。从表中可见，优化前初始输入变量对应矩阵Φ的条件数为6.0696，优化后为1.7365，说明采用优化输入信号进行激励，可以提高待估参数的可辨识度，获得有用信息最多的测量数据。

表1优化前后的输入信号设计变量取值对比表

将初始输入和优化输入分别代入真实系统(即η＝[200,400,300]^Tkg×m²的系统)进行激励，使用表1中的设计变量重构出的输入信号随时间变化曲线及其作用下得到的角速度轨迹如图3～4所示。从图中可以看出，由于输入信号使用了正弦基函数，状态轨迹初始时刻与终端时刻取值接近，即系统在激励结束后将回到初始状态，这种性质在进行实际物理试验时是非常有利的。输入轨迹初末时刻均为零，有利于物理执行机构准确复现该曲线。

3.3参数辨识结果对比

为了进行对比，引入两种航空航天领域经常采用的输入信号：doublet型输入信号与3211型输入信号。两种信号自身随时间变化的曲线及其作用下产生的角速度轨迹如图5～6所示，与图3～4对比可见，doublet型和3211型输入信号的曲线变化范围更大，其作用产生的角速度轨迹变化范围也更大。从直观上看，系统状态的大范围变化有利于进行高精度的参数辨识，但后面的参数辨识结果将显示，采用本发明提出的方法设计的最优输入能够获得精度更高的辨识结果。

算例中采用增广无迹卡尔曼滤波器作为具体的参数辨识算法，增广一词指的是将待估参数视为常值状态量增加到原状态方程中，实现状态量与参数同时估计，仿真获得的角速度轨迹中增加了均方差为1e-3rad/s的高斯白噪声作为观测数据，不同输入信号作用下的参数估计值变化曲线如图7～8所示。从图中可以看出，除了doublet输入外，在其他三种输入信号激励下估计值均能很快收敛到真实值(200kgm²)，在优化后的输入信号激励下，的收敛速度明显快于其他三种输入。对于优化后的输入信号和未优化的输入信号均可以使估计值收敛向真值(300kgm²)。从总体的参数估计误差||ΔJ||曲线来看，在优化后的输入信号激励下，参数估计算法的收敛速度和估计精度均优于其他三种类型的输入信号。

为更充分的说明优化后的输入信号可以提高参数估计的精度，对参数辨识过程进行多次蒙特卡洛仿真，并计算参数估计值的均值和方差。共进行了100次蒙特卡洛仿真，得到不同输入信号作用下的参数估计均值和方差列于表2中。可以看出，根据本发明方法设计的优化输入从总体上提高了参数估计的准确度。

表2不同输入信号激励下的参数估计均值，方差和相对误差表

4结论

最优输入设计的目的是为了提高系统辨识的效率和准确度。本发明针对一类参数可线性化的非线性系统，提出了一种基于直接配点法的最优输入设计方法。该方法的优化依据有别于传统的FIM依据和回归矩阵依据，采用最小二乘信息矩阵作为优化依据，构造出的性能指标函数只与终端时刻的最小二乘信息矩阵有关，减轻了OID问题对应最优控制问题性能指标函数值的计算复杂度。此外还采用了具有幅值和相位参数约束的正弦基函数参数化输入信号，构造的多输入正交函数初/末时刻取值均为零，还可以使终端时刻状态轨迹回到初始位置附近，为系统辨识物理实验提供了便利。采用直接配点法离散状态变量，设计的串行优化策略不仅保证了非线性规划问题的解为原优化问题的解，还极大地提高了计算效率。通过一个简单的受控刚体惯量矩辨识算例，证明了通过本发明提出的设计方法，可以加速估计的收敛速度，并有效提高参数估计的精度。

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。

Claims (5)

1.一种提高航天器惯性参数辨识精度的输入信号设计方法，其特征在于，包括以下步骤：

步骤S100：选取在轨运行航天器中的未知惯性参数作为待估参数，将含所述待估参数的动力学微分方程以线性回归方程形式表示；

步骤S200：根据所述线性回归方程中的回归矩阵构造最小二乘法矩阵，通过步骤S300～步骤S400获得能够使所述最小二乘法矩阵的条件数接近1的输入信号；

步骤S300：以所述最小二乘法矩阵中各不相同的独立元素作为增广状态变量，扩充到所述航天器动力学微分方程中，得到增广状态微分方程；

步骤S400：将最优输入设计问题归纳为最优控制问题，并对所述最优控制问题进行求解，得到用于参数辨识实验所需的所述输入信号。

2.根据权利要求1所述的提高航天器惯性参数辨识精度的输入信号设计方法，其特征在于，所述步骤S400中所述最优控制问题的解决方法包括以下步骤：

步骤S401：以所述增广状态微分方程作为描述航天器惯性系统运动的状态方程；

步骤S402：根据所述航天器的姿控执行机构的实际物理局限，确定各通道输入信号的允许取值范围；

步骤S403：设定最小二乘法矩阵独立元素的初值为零，结合航天器角速度的真实初始状态，共同构成增广状态微分方程的初始状态，对增广状态微分方程的终端状态不设约束；

步骤S404：用静态的参数对动态的所述输入信号进行描述；

设所有通道的输入信号均为不同频率的正弦基函数之和，构成各通道输入信号的正弦基函数具有独立的相位、频率和幅值，设第i通道输入信号的所述正弦基函数的频率系数取值集合为：

S_i＝{(i+1),(i+1)+m,(i+1)+2m,(i+1)+3m,...}

其中，m为总通道数，这样使得任意两通道的输入信号相互正交；

设构成所述第i通道输入信号的第k个正弦基函数的幅值和相位分别为A_k和ψ_k，这样动态的输入信号就可以用静态的参数{A_k,ψ_k}进行描述，比如第i个通道的输入信号被描述为：

$U_i\left(t\right)=\underset{k\&Element;S_i}{\&Sigma;}{A}_{k}sin(\frac{2\mathrm{\&pi;}kt}{T}+{\mathrm{\&psi;}}_k),i=1,2,...,m$

步骤405：采用直接配点法，对动态的增广状态轨迹进行参数化，将动态的最优输入设计问题转化为静态的非线性规划问题；

步骤S406：采用串行优化策略，对所述静态非线性规划问题进行求解，获得优化后的描述输入信号的静态参数；

步骤S407：使用所述优化后的输入信号静态参数，重构得到动态的最优输入信号。

3.根据权利要求2所述的提高航天器惯性参数辨识精度的输入信号设计方法，其特征在于，

所述步骤S405为：将输入持续时间区间进行分段，每段的两个端点称为节点，记节点处对应的状态变量为{X₁,X₂,...,X_N}，子区间[t_i,t_i+1]上的状态变量X(t)用三次Hermite插值多项式近似表示。

4.根据权利要求3所述的提高航天器惯性参数辨识精度的输入信号设计方法，其特征在于，

所述步骤S404中还包括设置幅值和相位满足等式使得初始时刻t₀和终端时刻t_f输入信号值为0。

5.根据权利要求2所述的提高航天器惯性参数辨识精度的输入信号设计方法，其特征在于，所述约束条件还包括状态轨迹的边界条件约束和不等式约束，所述不等式约束由实际系统的物理约束确定。