[发明专利]一种基于COMI‑PSO算法的不相关多源频域载荷识别方法

权利要求书

1.一种基于COMI-PSO算法的多源频域载荷识别方法，适用于响应测点的个数大于等于激励的个数、线性不变的系统、平稳随机的载荷激励、各个载荷源近似不相关或独立，其特征在于包括如下步骤：

步骤1、设系统有m个载荷激励输入f_i(i=1,…,m)，n个测点输出y_j(j=1,…,n)，对应于每一个输出y_j，有m个脉冲响应函数h_ji(t),i=1,…,m，j=1,…n,测得多输入多输出系统在m个不相关载荷同时作用下n个响应测点的时域响应y_j(t)j=1,2,…,n，计算响应测点间的互功率谱密度矩阵S_yy(ω)，其中ω为频率；

步骤2、对系统依次分立对各个载荷点施加单个激励输入f_i(i=1,…,m)，计算其自功率谱为测得系统在该单独激励输入下各响应测点的响应输出并计算其自功率谱按照(15)式的一元线性回归模型识别出的传递函数模的平方

β^0=S‾yyj,i(ω)-|H^j,i(ω)|2S‾fi(ω)|H^j,i(ω)|2=Σl=1k[Syyj,il(ω)-S‾yyj,i(ω)][Sfil(ω)-S‾fi(ω)]Σl=1k[Sfil(ω)-S‾fi(ω)]2---(15)]]>

(15)式中，为载荷点i所施加的k次载荷的平均值，为响应点j在载荷点i所施加的k次响应的平均值；

步骤3、利用载荷源的不相关特性对载荷识别方程(1)式进行解耦得到（4）式：

Syy(ω)=12π∫-∞∞∫-∞∞∫-∞∞h(u)Cff(τ+u-v)·hT(v)e-jωτdudvuτ=H‾(ω)Sff(ω)HT(ω)---(1)]]>

(1)式中h(u)是系统的单位脉冲响应矩阵，h^T(u)是系统的单位脉冲响应矩阵的转置，是输入各激励之间的协方差函数矩阵，是系统频率特性矩阵，H^T(ω)是系统频率特性矩阵的转置，是系统频率特性矩阵的共轭，S_ff(ω)是输入各激励之间的互功率谱密度，

在m个输入载荷激励都是零均值的平稳随机过程，且互不相关的情况下，输入协方差函数矩阵为对角阵，即：输入功率谱矩阵S_ff(ω)也为对角阵此时，输出功率谱中主对角线上的任意一元素Syyjj(ω)(j=1,2,···,n)]]>满足：

Syyjj(ω)=H‾j1(ω)···H‾ji(ω)···H‾jn(ω)·diag[Sffii(ω)]·Hj1(ω)···Hji(ω)···Hjn(ω)T=Σi=1mH‾ji(ω)Sffii(ω)HjiT(ω)=Σi=1m|Hji(ω)|2Sffii(ω)---(3)]]>

(3)式写成矩阵后的形式为：

其中，|H_j,i(ω)|²是输入f_i对响应y_j的传递函数模的平方，是待识别的载荷源f_i的自功率谱，是响应y_j的自功率谱；

记S→Y(ω)=ΔSyy11(ω)···Syynn(ω)T]]>，

则(4)式可简写为：S→Y(ω)=B(ω)S→F(ω);]]>

步骤4、对(4)式识别出的S_yy(ω)和采用(16)式传统最小二乘广义逆计算(4)式载荷值：

S→F′(ω)=[B(ω)TB(ω)]-1B(ω)TS→Y(ω)---(16)]]>

步骤5、利用(27)式条件数的定义计算（4）式的条件数，若条件数小于阈值，直接使用（16）式的载荷值作为最终的识别载荷值,载荷识别结束；如果条件数大于阈值，则使用（16）式的计算结果控制粒子群的初始化位置，并随机初始化粒子的初试速度；

在求解方程Ax=b的过程中，引入的舍入误差总会导致其数值解x或多或少地不等于其理论解x'=A^-1B，

定义误差为：δx=Δx-x′---(25)]]>

定义剩余为：δb=ΔAx-b---(26)]]>

定义条件数为：cond(A)=Δ||A-1||||A||---(27);]]>

步骤6、以(20)式作为适应值函数，计算每个粒子的适应值：

minf(Sff′(ω))=min(maxj{|Syyjj(ω)-Σi=1m|Hj,i(ω)|2Sffii′(ω)||Syyjj(ω)|})---(20)]]>

步骤7、找到每个粒子的历史最优位置和最优适应值：

对每个粒子而言，如果它的当前适应度比其历史最优适应度还要小，则使用当前的适应度替代该粒子的最优适应度，并保存当前位置为该粒子历史最优位置；

步骤8、找到整个群体的最优位置和最优适应值；

对每个粒子而言，如果它的当前适应度比整个群体最优适应度还要小，则使用当前的适应度替代粒子群最优适应度，并保存当前位置为整个群体最优位置；

步骤9、根据COMI-PSO算法，使用(21)式、(22)式计算各个粒子最新速度和最新位置

vidt+1=w*vidt+c1r1(pid-xidt)+c2r2(pgd-xidt)---(21)]]>

xidt+1=xidt+vidt+1---(22)]]>

其中惯性参数w=max_W-(max_W-min_W)*ln(1+(e-1)*CT/TS)，max_W为惯性参数的上限，min_W为惯性参数的下限，e为自然对数的底数，上式能保证w的取值范围是[min_W,max_W]；r₁,r₂为(0,1)之间的随机数，p_id表示当前群体中第i个个体历史最优位置在第d维上数值，p_gd表示当前群体最优位置在第d维上的数值，CT为当前的迭代次数，TS为预设的总迭代次数，加速因子c₁,c₂分别为：

c₁=4*(w-min_W)²/(max_W-min_W)²

c₂=4*(max_W-w)²/(max_W-min_W)²

通过上述操作完成迭代过程中第i个个体从第t代的速度和位置更新为第(t+1)代新的速度和位置

步骤10、通过(20)式计算新位置下，每个粒子的适应度，根据遗传算法中的选择算子，保留种群中适应度较好的部分粒子直接进入下一次迭代，对另一部分适应度较小的粒子采用遗传算法中的交叉算子和突变算子进行预处理，计算预处理之后的子代的适应值，将适应度比父代好的部分相同数目的粒子取代原粒子群中的父辈进入下一次迭代；

步骤11、根据是否达到最大迭代次数或达到最好的适应值判断迭代终止条件，如果达到转入步骤12，否则转入步骤6；

步骤12、输出全局最优解，即群体最优位置，作为最终识别的多个不相关载荷源值,载荷识别算法结束。