CN105404609A - 一种新型多目标动力系统参数辨识方法 - Google Patents

Abstract

本发明公开了一种新型多目标动力系统参数辨识方法，包括以下步骤：(1)针对待辨识动力系统建立运动微分方程；(2)选取多个待辨识动力系统输出响应的特征参数作为评价指标，并通过试验获得评价指标的值及待辨识动力系统的输入数据；(3)应用改进的分块增维精细积分法和待辨识动力系统的输入数据，通过待辨识动力系统的运动微分方程计算评价指标的值，将计算获得的评价指标的值与试验获得的评价指标的值相减获得差值，以差值的平方和作为待辨识动力系统的目标函数；(4)应用NSGA-II算法对目标函数进行多次寻优计算，直至输出结果稳定，获得最终的非劣解集；(5)选择结果。该方法效率高，适用面广，稳定性高。

Description

一种新型多目标动力系统参数辨识方法

技术领域

本发明涉及系统参数辨识领域，更具体的说，是涉及一种新型多目标动力系统参数辨识方法。

背景技术

参数辨识是根据实验数据和建立的模型来确定一组参数值，使得由模型计算得到的数值结果能最好的拟合测试数据，从而可以为生产过程进行预测，提供一定的理论指导，在工程实际中拥有非常广泛的应用。传统参数辨识一般采用单目标辨识，或者将多目标辨识进行线性加权转化成单目标辨识。单目标辨识方法简便，但准确性较差，由于目标函数少，辨识结果很难全面反映被辨识模型的全部特性，而线性加权的多目标参数辨识方法在使用过程中必须首先确定每个目标函数的权重值，而权重值的确定十分困难，因此在实际使用过程中通常采用经验值，而这又会极大地降低辨识结果的准确性。所以，开发一种不需要线性加权的新型多目标参数辨识方法是十分必要的。

发明内容

本发明的目的是为了克服现有技术中的不足，提供一种新型多目标动力系统参数辨识方法，该方法从功能上可以分为两部分，第一部分是计算目标函数的值，第二部分是对目标函数进行寻优计算以获得最优解。为了实现上述功能，本发明方法主要由两种算法构成，分别为改进的分块增维精细积分法和NSGA-II算法；其中，改进的分块增维精细积分法主要用来计算目标函数的值，NSGA-II算法主要用来对目标函数进行寻优计算。

本发明的目的是通过以下技术方案实现的：

一种新型多目标动力系统参数辨识方法，包括以下步骤：

(1)选定要进行参数识别的待辨识动力系统，根据待辨识动力系统的结构建立运动微分方程，并给出运动微分方程中已知参数的值和待辨识参数的范围及约束条件；

(2)选取多个待辨识动力系统输出响应的特征参数作为评价指标，并通过试验获得评价指标的值及待辨识动力系统的输入数据；

(3)应用改进的分块增维精细积分法和待辨识动力系统的输入数据，通过待辨识动力系统的运动微分方程计算评价指标的值，将计算获得的评价指标的值与试验获得的评价指标的值相减获得差值，以差值的平方和作为待辨识动力系统的目标函数；

(4)寻优计算，对NSGA-II算法进行相应的参数设定，应用NSGA-II算法对目标函数进行多次寻优计算，直至输出结果稳定，获得最终的非劣解集；

(5)选择结果，结合工程实际，在非劣解集中选取合适的点为最优点，以最优点对应的参数为最终辨识结果。

步骤(5)中所述最优点的选取原则是最优点使不同的目标函数值都达到最小。

与现有技术相比，本发明的技术方案所带来的有益效果是：

1.本发明方法中同时采用了改进的分块增维精细积分法和NSGA-II算法，其中改进的分块增维精细积分法是在精细积分法的基础之上发展而来的，它不但具有精细积分法求解精度高，稳定性和收敛性好的特点，而且克服了精细积分法不能求解非线性系统的缺点，提高了求解效率；NSGA-II算法具有计算量小、搜索效率高的特点，在保持种群多样性和防止优秀个体流失方面表现突出，而且不需要设定初值。

2.本发明方法的优点在于目标函数的个数不受限制，求解时不需要对多个目标函数进行线性加权，并且由于改进的分块增维精细积分法可以方便地求解线性/非线性偏微分方程任意精度的数值解，使得本发明的多目标动力系统参数辨识方法适用范围非常广泛，只要系统的模型可以转化为偏微分方程，均可以使用本方法进行系统辨识，尤其适合动力系统的参数辨识，无论是线性还是非线性动力系统，无论有无外部激励，无论激励是谐波激励还是随机激励，本方法都能很方便的应用。

附图说明

图1是履带急救车非线性减振系统结构示意图。

图2是履带急救车非线性减振系统四自由度振动模型。

图3是本发明方法的流程图。

图4是本发明方法的计算机操作流程图。

图5是应用本发明方法对履带急救车非线性减振系统进行参数辨识所获得的第一前端个体分布示意图。

图6是由辨识结果计算获得的担架台垂向振动加速度功率谱密度和由试验数据计算获得的担架台垂向振动加速度功率谱密度的对比示意图。

图7是由辨识结果计算获得的担架台垂向振动加速度概率分布和由试验数据计算获得的担架台垂向振动加速度概率分布的对比示意图。

附图标记：1-大板车厢2-担架台3-零刚度减振器4-橡胶阻尼减振器5-履带底盘M₁-担架抬质量M₂-大板车厢质量Ks-零刚度减振器的线性刚度Kz-零刚度减振器的三次非线性刚度C₁-零刚度减振器的阻尼K₂-橡胶阻尼减振器的刚度C₂-橡胶阻尼减振器的阻尼J₁、J₂-担架台和大板车厢的质量在几何中心处绕z轴的转动惯量l₁、l₂-担架台和大板车厢的几何中心距离各自减振器沿y轴的距离q₁、q₂-振动输入fun1-目标函数1fun2-目标函数2

具体实施方式

为了更加清楚的叙述本发明的实施过程，下面结合工程实例及说明书附图对本发明作进一步的详细描述：

图3为本发明方法的流程图，应用本发明所提出的新型多目标动力系统参数辨识方法对某履带急救车非线性减振系统进行参数辨识，履带急救车非线性减振系统结构示意图如图1所示。忽略其他两个方向，仅保留对成员舒适性影响最大的垂直方向，建立减振系统四自由度振动模型如图2所示，其中Ks、Kz、C₁、K₂和C₂作为待辨识参数，M₁、M₂、J₁、J₂、l₁和l₂为已知参数，取值如下：M₁＝180kg、M₂＝2000kg、J₁＝58kg·m2、J₂＝3665kg·m2、l₁＝0.8m、l₂＝1.539m。

本发明方法中所提出的新型多目标动力系统参数辨识方法由改进的分块增维精细积分法和NSGA-II算法构成，具体使用时需要通过计算机编程来实现，本实施例采用MATLAB编程，NSGA-II算法采用MATLAB内置的gamultiobj函数来实现，gamultiobj函数所采用的是NSGA-II算法的一种变形。而改进的分块增维精细积分法需要自行编程实现，下面结合履带急救车非线性减振系统的振动模型对改进的分块增维精细积分法的计算过程进行介绍。

对图2所示振动模型建立运动微分方程如下：

$M\stackrel{\·\·}{X}+C\stackrel{\·}{X}+KX=f(X,\stackrel{\·}{X},t)---\left(1\right)$

其中， $X=\left[\begin{array}{c}{x}_1\\ {\mathrm{\φ}}_1\\ x_2\\ {\mathrm{\φ}}_2\end{array}\right],M=\left[\begin{array}{cccc}{M}_1& 0& 0& 0\\ 0& J_1& 0& 0\\ 0& 0& M_2& 0\\ 0& 0& 0& J_2\end{array}\right],K=\left[\begin{array}{cccc}2K_s& 0& -2K_s& 0\\ 0& 2K_s{l}_1^2& 0& -2K_s{l}_1^2\\ -2K_s& 0& 2(K_s+K_2)& 0\\ 0& -2K_s{l}_1^2& 0& 2(K_s{l}_1^2+K_2{l}_2^2)\end{array}\right],$

$C=\left[\begin{array}{cccc}2C_1& 0& -2C_1& 0\\ 0& 2C_1{l}_1^2& 0& -2C_1{l}_1^2\\ -2C_1& 0& 2(C_1+C_2)& 0\\ 0& 2C_1{l}_1^2& 0& 2(C_1{l}_1^2+C_2{l}_2^2)\end{array}\right],$

$f\left(X,\stackrel{\·}{X},t\right)=\left[\begin{array}{c}-K_z{\[\left(x_1-x_2\right)-l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3-K_z{\[\left(x_1-x_2\right)+l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3\\ K_z{l}_1{\[\left(x_1-x_2\right)-l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3-K_z{l}_1{\[\left(x_1-x_2\right)+l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3\\ K_z\left(q_1+q_2\right)+C_2\left({\stackrel{\·}{q}}_1+{\stackrel{\·}{q}}_2\right)+K_z{\[\left(x_1-x_2\right)-l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3+K_z{\[\left(x_1-x_2\right)+l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3\\ -K_z{l}_2\left(q_1-q_2\right)-C_2{l}_2\left({\stackrel{\·}{q}}_1-{\stackrel{\·}{q}}_2\right)-K_z{l}_1{\[\left(x_1-x_2\right)-l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3+K_z{l}_1{\[\left(x_1-x_2\right)+l_1\left({\mathrm{\φ}}_1-{\mathrm{\φ}}_2\right)\]}^3\end{array}\right].$

对运动微分方程中的K_s、K_z、C₁、K₂和C₂进行参数辨识，首先需要将式(1)改写如下形式

$\stackrel{\·}{\stackrel{\‾}{X}}=A\stackrel{\‾}{X}+F---\left(2\right)$

式中 $\stackrel{\‾}{X}={\left[\begin{array}{c}X\\ \stackrel{\·}{X}\end{array}\right]}_{2n\×1},A={\left[\begin{array}{cc}0& I_n\\ -M^{-1}K& -M^{-1}C\end{array}\right]}_{2n\×2n},F={\left[\begin{array}{c}0\\ M^{-1}f\end{array}\right]}_{2n\×1},$ I_n为n阶单位矩阵。

引入新的变量X_2n+1＝1，则式(2)增加一维后变成齐次方程形式，即

$\stackrel{\·}{Y}=\left[\begin{array}{c}\stackrel{\·}{\stackrel{\‾}{X}}\\ {\stackrel{\·}{X}}_{2n+1}\end{array}\right]=\left[\begin{array}{cc}A& F\\ 0& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\‾}{X}\\ X_{2n+1}\end{array}\right]=HY---\left(3\right)$

记时间步长为τ。令Δt＝τ/m，其中m＝2^L，L＝20。经过推导可得如下计算过程

$\left\{\begin{array}{c}{A}_{\mathrm{\α}}=A\mathrm{\Δ}t+{\left(A\mathrm{\Δ}t\right)}^2/2+{\left(A\mathrm{\Δ}t\right)}^3/6\\ B=I\mathrm{\Δ}t+{\mathrm{A\Δt}}^2/2+A^2{\mathrm{\Δt}}^3/6\end{array}\right.---\left(4\right)$

$\begin{array}{c}for(iter=0;iter<L;iter++)\\ \left\{B=\left(2I+A_{\mathrm{\&alpha;}}\right)B;A_{\mathrm{\&alpha;}}=2A_{\mathrm{\&alpha;}}+A_{\mathrm{\&alpha;}}^2\right\}\end{array}---\left(5\right)$

当循环(5)结束后，可得

$T_{k-1}=I+\left[\begin{array}{cc}{A}_{\mathrm{\&alpha;}}& BF\\ 0& 0\end{array}\right]=\left[\begin{array}{cc}{A}_{\mathrm{\&alpha;}}+I_{2n\&times;2n}& BF\\ 0& 0+I_{1\&times;1}\end{array}\right]=\left[\begin{array}{cc}{A}_{\mathrm{\&alpha;}}^I& G\\ 0& I_{1\&times;1}\end{array}\right]---\left(6\right)$

Y(kτ)＝T_k-1Y[(k-1)τ](k＝1,2,…)(7)

对矩阵B和矩阵F进行如下分块：

$B_{2n\&times;2n}=\begin{array}{cc}\left[\begin{array}{cc}{B}_{2n\&times;n}^1& B_{2n\&times;n}^2\end{array}\right]& F_{2n\&times;1}\end{array}={\left[\begin{array}{cc}{F}_{n\&times;1}^1& F_{n\&times;1}^2\end{array}\right]}^T---\left(8\right)$

由于矩阵F_2n×1的前n行为0，分块以后得到所以矩阵G的计算公式变为

$G=B_{2n\&times;2n}\&times;F_{2n\&times;1}=B_{2n\&times;n}^2\&times;F_{n\&times;1}^2---\left(9\right)$

考虑增维后H的特点，Y的最后一行实际上不会对所求结果产生影响，则式(7)可变成如下形式

$Y_{2n\&times;1}\left(k\mathrm{\&tau;}\right)=A_{\mathrm{\&alpha;}}^I{Y}_{2n\&times;1}\&lsqb;(k-1)\mathrm{\&tau;}\&rsqb;+G---\left(10\right)$

即 $\stackrel{\&OverBar;}{X}\left(k\mathrm{\&tau;}\right)=A_{\mathrm{\&alpha;}}^I\stackrel{\&OverBar;}{X}\&lsqb;(k-1)\mathrm{\&tau;}\&rsqb;+G---\left(11\right)$

公式(2)-(11)即为改进的分块增维精细积分法的计算过程，具体应用时需要将公式(2)-(11)转换为计算机程序。

以上即为构成新型多目标动力系统参数辨识方法的两种基本算法的实现方式。

在参数辨识之前需要对履带急救车进行野外行驶振动试验以获得履带底盘5、大板车厢1、担架台3等部位的振动加速度数据信号，为本方法提供数据基础。

如图4所示，是利用本发明的实施方法编制出的计算机流程图，具体实施步骤包括如下：

(1)应用MATLAB的load函数导入由履带急救车野外行驶振动试验获得的非线性减振系统的振动数据，并存储在相应的向量内；

(2)对履带急救车非线性减振系统的已知参数、改进的分块增维精细积分法的积分步长、计算过程中所需的向量空间等进行设定；

(3)应用MATLAB的pwelch函数和ksdensity函数计算由试验获得的振动数据中担架台垂向振动加速度的功率谱密度pxx和概率分布hx，并存储在相应的向量之中；

(4)根据式(2)的形式建立履带急救车非线性减振系统振动模型的运动微分方程，在计算程序中只需要给出矩阵A和矩阵F即可；

(5)将公式(2)-(11)编写成MATLAB程序段，并应用该程序段对履带急救车非线性减振系统振动模型的运动微分方程进行求解，获得履带急救车非线性减振系统的振动数据，并存储在相应的向量内；

(6)应用MATLAB的pwelch函数和ksdensity函数计算由仿真获得的振动数据中担架台垂向振动加速度的功率谱密度pxx1和概率分布hx1，并存储在相应的向量之中；

(7)以pxx与pxx1、hx与hx1差值的平方和建立目标函数1和目标函数2，分别表示为fun1和fun2；

(8)应用gaoptimset函数对NSGA-II算法的最优前端个体系数、种群大小、最大进化代数、停止代数、适应度函数值偏差等进行设置；

(9)给出待辨识参数的范围以及约束条件；

(10)应用gamultiobj函数对目标函数进行优化计算，获得第一前端个体分布以及第一前端个体分布所对应的非劣解集；

(11)从非劣解集中选取最优点，并以最优点对应数据为最终辨识结果。

以上即为本发明方法的计算机实现过程，具体编程时，第(1)-(7)步编写在一个M文件里，第(8)-(10)步编写在另一个M文件里，第(11)步为选择参数辨识结果，选择的原则是尽可能使两个目标函数值都达到最小。

应用本发明所提出的新型多目标动力系统参数辨识方法对工程实例中的履带急救车非线性减振系统进行参数辨识，第一前端个体分布(ParetoFront)如图5所示，第一前端个体分布对应的非劣解集如表1所示，其中加粗栏为参数辨识结果。

表1

将参数辨识结果重新带入图2所示模型的运动微分方程中计算担架台垂向振动加速度的功率谱密度和概率分布(以“仿真结果”来表示)，并与由试验数据计算获得的功率谱密度和概率分布(以“试验结果”来表示)进行对比，如图6和图7所示；观察图6和图7可以发现，试验结果与仿真结果基本吻合，尤其是在概率分布方面，吻合程度较高，这充分证明了本发明所提出的新型多目标动力系统参数辨识方法的有效性。

Claims (2)

1.一种新型多目标动力系统参数辨识方法，其特征在于，包括以下步骤：

(1)选定要进行参数识别的待辨识动力系统，根据待辨识动力系统的结构建立运动微分方程，并给出运动微分方程中已知参数的值和待辨识参数的范围及约束条件；

(2)选取多个待辨识动力系统输出响应的特征参数作为评价指标，并通过试验获得评价指标的值及待辨识动力系统的输入数据；

(3)应用改进的分块增维精细积分法和待辨识动力系统的输入数据，通过待辨识动力系统的运动微分方程计算评价指标的值，将计算获得的评价指标的值与试验获得的评价指标的值相减获得差值，以差值的平方和作为待辨识动力系统的目标函数；

(4)寻优计算，对NSGA-II算法进行相应的参数设定，应用NSGA-II算法对目标函数进行多次寻优计算，直至输出结果稳定，获得最终的非劣解集；

(5)选择结果，结合工程实际，在非劣解集中选取合适的点为最优点，以最优点对应的参数为最终辨识结果。

2.根据权利要求1所述的一种新型多目标动力系统参数辨识方法，其特征在于，步骤(5)中所述最优点的选取原则是最优点使不同的目标函数值都达到最小。