一种线性周期时变系统的模态振型分析方法及装置与流程

1.本发明属于线性周期时变系统动力学分析领域，更具体地，涉及一种线性周期时变系统的模态振型分析方法及装置。


背景技术：

2.在系统动力学分析领域，线性系统由于不存在非线性关系，系统状态变量之间满足叠加原理，使得线性系统解的性态得到了较为深入的研究。以状态空间的形式列写线性系统微分代数方程如下：
[0003][0004]
其中x(t)＝[x1(t),
…
,xn(t)]
t
为系统的n维状态向量；u(t)＝[u1(t),
…
,u
p
(t)]
t
为系统的p维输入向量；y(t)＝[y1(t),
…
,yq(t)]
t
为系统的q维输出向量；n为系统的阶次；a(t)、b(t)、c(t)和d(t)分别为n
×
n维、n
×
p维、q
×
n维、q
×
p维的系数矩阵。
[0005]
依据时变系数矩阵a(t)、b(t)、c(t)、d(t)所满足的性质，线性系统可进一步划分。其中，当a(t)、b(t)、c(t)、d(t)为定常矩阵时，系统为线性时不变系统。当a(t)、b(t)、c(t)、d(t)为一般的时变矩阵时，系统为线性时变系统。当时变系数矩阵a(t)满足式(2)关系时，系统为线性周期时变系统。
[0006]
a(t)＝a(t t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0007]
线性周期时变系统相比线性时不变系统更为一般，诸多物理、工程问题往往归结于线性周期时变系统的解的分析，如天体运动、离子加速器设计、姿态控制等。对于线性时变系统，其难以统一刻画的时变特性造成了一般分析方法的缺乏，往往需要进行一些假设后转化为线性周期时变系统进行研究。因此线性周期时变系统的研究对于线性系统动力学分析具有独特而重要的地位。
[0008]
现有线性周期时变系统分析方法主要以定性和解析两种思想作为理论基础，定性理论又以g.floquet与lyapunov的工作为代表。不考虑输入u(t)，floquet首先揭示了线性周期时变系统的解具有式(3)所示特征结构。
[0009]
x(t)＝p(t)e
qt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0010]
其中p(t)为周期矩阵，q为定常矩阵。
[0011]
该特征结构表明线性周期时变系统虽然相比线性时不变系统多引入了周期性因素，但在方程的解中仍然存在状态模式关联特性与稳定性相互解耦的特点。进一步地，lyapunov揭示了线性周期时变系统的可约性，即存在周期矩阵p(t)将系统矩阵a(t)转换为定常矩阵q：
[0012][0013]
lyapunov证明了时不变化变换矩阵p(t)的存在性，但对于一般线性周期时变系统，由于缺少解析构造p(t)的方法，因此基于定性理论发展的分析工具，如基于状态转移矩
阵(state transition matrix,stm)计算系统特征指数矩阵q，面临着无法揭示系统振荡模式频率的问题。进一步结合floquet所指出的线性周期时变系统状态模式关联特性与稳定性之间的解耦特点，系统状态模式之间的关联信息存在于时不变化变换矩阵p(t)之中，p(t)的缺失将造成无法定量刻画线性周期时变系统中的模态振型信息。
[0014]
线性周期时变系统解析理论基于级数逼近与谐波平衡在无穷阶收敛的情况下将a(t)转换为无穷阶定常矩阵进行分析。
[0015][0016]
将式(5)代入式(4)可得：
[0017][0018]
将式(6)写成式(7)所示矩阵形式，并通过无穷阶定常矩阵特征向量求解出pn(n为整数)，从而可以还原得到时不变化变换矩阵p(t)。
[0019][0020]
由于实际计算时需要对无穷阶矩阵进行截断处理，截断阶数直接影响数学关系的准确性与计算结果的收敛性，目前的理论尚无法对截断阶数的合理性提供有力支撑。因此在线性周期时变系统中，基于解析理论进行模态振型分析也尚不完备。
[0021]
因此，亟需提供一种有效分析计算线性周期时变系统模态振型的方法，从而保证能准确进行系统模态辨识、相对振荡状态量划分等系统动态稳定特性的分析。


技术实现要素：

[0022]
针对现有技术的缺陷，本发明提供了一种线性周期时变系统的模态振型分析方法及装置，其目的在于，判断线性周期时变系统振荡频率，确定振荡存在于哪一部分系统状态量之间，为准确分析线性周期时变系统动态稳定特性提供理论支撑。
[0023]
按照本发明的一个方面，提供了一种线性周期时变系统的模态振型分析方法及装置，包括以下步骤：
[0024]
步骤1：基于线性周期时变系统状态空间模型，通过数值计算求解系统状态转移矩阵在初始时刻为零和观测时刻为最小周期t的值。根据状态转移矩阵的物理意义，在其作用下，系统状态由初始时刻的状态值x(t0)转移到观测时刻t的状态x(t)，因此，当初始时刻t0为零时刻，观测时刻t为最小周期t时，状态转移矩阵φ(t,0)满足：
[0025][0026]
由于状态转移矩阵φ(t,0)不受初始状态的影响，因此初始状态可以选择任意n维的极大线性无关组。本发明直接选择n维单位矩阵i的n个列向量作为系统的初始状态，可以
得到线性周期时变系统n个状态响应x(t)，即：
[0027][0028][0029][0030][0031]
通过上述过程可以得到状态转移矩阵φ(t,0)的n个列向量，于是状态转移矩阵可以表示为：
[0032]
φ(t,0)＝[φ1(t,0) φ2(t,0)
…
φn(t,0)]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0033]
步骤2：根据floquet-lyapunov理论计算线性周期时变系统时不变化得到的定常矩阵q，其中q阵与步骤1所得的状态转移矩阵φ(t,0)满足式(13)所示关系：
[0034][0035]
步骤3：基于线性周期时变系统状态空间模型，选择n维单位矩阵i的n个列向量分别作为系统的初始状态，通过数值计算获得系统状态转移矩阵在一个周期t内的数值解：
[0036][0037]
步骤4：根据线性周期时变系统解的结构，基于系统状态初值x(0)，定常矩阵q与状态转移矩阵φ(t,0)求解时不变化变换矩阵p(t)
[0038]
p(t)＝φ(t,0)x(0)e-qt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0039]
步骤5：根据时不变化变换矩阵p(t)的周期特性，对(15)计算得到的p(t)进行傅里叶分解：
[0040][0041]
将(16)回代入(3)中判断振型信息：
[0042]
[0043]
步骤6：通过比较矩阵pn中各元素p
n(rs)
模值大小进行系统模态辨识；
[0044]
步骤7：通过比较矩阵pn各元素的相位进行相对振荡状态量划分，对于不同子系统中具有相同物理意义的状态量，相位相差趋近0
°
被划分为共模振荡，相位相差趋近180
°
被划分为差模振荡。
[0045]
作为进一步改进，步骤1与步骤3可以同时进行，求解一个周期内的所有状态响应后，提取初始时刻与一个周期结束后时刻的状态响应作为计算定常矩阵q的基础。
[0046]
作为进一步改进，步骤1与步骤3所得到的状态转移矩阵会受到解法器、数值计算步长与误差控制的影响，需要根据实际需求进行合理选择。
[0047]
按照本发明的另一个方面，提供了一种线性周期时变系统的模态振型分析装置，包括：系统状态转移矩阵求解模块、定常矩阵求解模块、时不变化变换矩阵模块和模态振型分析模块；
[0048]
所述状态转移矩阵求解模块，用于基于线性周期时变系统状态空间模型，选择n维单位矩阵i的n个列向量分别作为系统的初始状态，通过数值计算求解系统状态转移矩阵在初始时刻为零和观测时刻为最小周期t的值φ(t,0)以及系统状态转移矩阵在一个周期t内的数值解φ(t,0)；
[0049]
定常矩阵求解模块，用于根据floquet-lyapunov理论计算线性周期时变系统时不变化得到的定常矩阵q；
[0050]
所述时不变化变换矩阵模块，用于根据线性周期时变系统解的结构，基于系统状态初值x(0)，定常矩阵q与状态转移矩阵φ(t,0)求解时不变化变换矩阵p(t)；
[0051]
所述模态振型分析模块，用于根据时不变化变换矩阵p(t)的周期特性，对p(t)进行傅里叶分解，利用傅里叶分解后的系数矩进行断模态振型分析。
[0052]
按照本发明的又一个方面，提供了一种计算机可读存储介质，包括存储的计算机程序。计算机程序被处理器执行时，控制计算机可读存储介质所在设备执行本发明提供的线性周期时变系统的模态振型分析方法。
[0053]
通过本发明所构思的以上技术方案，与现有技术相比，能够取得以下有益效果：
[0054]
1.本发明通过融合线性周期时变系统定性理论与数值计算，避免了构造时不变化变换矩阵p(t)没有一般解析方法的缺陷，将问题转化为数值计算系统状态转移矩阵φ(t,0)、φ(t,0)与定常矩阵q，进而准确求解时不变化变换矩阵p(t)。
[0055]
2.本发明基于线性周期时变系统定性理论作为求解时不变化变换矩阵p(t)的基础，避免了解析理论中无穷维定常矩阵所面临的截断问题，准确度更高。
[0056]
3.本发明提出的线性周期时变系统基于时不变化变换矩阵p(t)的模态振型分析方法，能够有效进行系统模态辨识并划分相对振荡的状态量，为进一步分析线性周期时变系统动态稳定特性提供支撑。
附图说明
[0057]
图1为本发明实施例提供的线性周期时变系统时不变化变换矩阵计算方法流程图；
[0058]
图2为本发明实施例提供的线性周期时变系统状态转移矩阵求解流程图。
具体实施方式
[0059]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间不构成冲突就可以相互组合。
[0060]
一种线性周期时变系统时不变化变换矩阵计算方法，如图1所示，包括：状态转移矩阵获取步骤、定常矩阵获取步骤以及时不变化变换矩阵获取步骤。
[0061]
状态转移矩阵φ(t,0)获取步骤(步骤1)：基于线性周期时变系统状态空间模型，通过数值计算求解系统状态转移矩阵在初始时刻为零和观测时刻为最小周期t的值。根据状态转移矩阵的物理意义，在其作用下，系统状态由初始时刻的状态值x(t0)转移到观测时刻t的状态x(t)，因此，当初始时刻t0为零时刻，观测时刻t为最小周期t时，矩阵φ(t,0)满足：
[0062][0063]
由于状态转移矩阵φ(t,0)不受初始状态的影响，因此初始状态可以选择任意n维的极大线性无关组。
[0064]
为简化相关计算，作为一种优选的实施方式，本实施例直接选择n维单位矩阵i的n个列向量作为系统的初始状态，可以得到线性周期时变系统n个状态响应x(t)，即：
[0065][0066]
状态变量初值赋予规则为xk(0)＝ek，计数器k遍历1到n的所有整数值，xk(0)表示第k次数值求解时的状态变量初值所形成的列向量。通过上述过程可以得到状态转移矩阵φ(t,0)的n个列向量，相应的流程图如图2所示，状态转移矩阵可以表示为：
[0067]
φ(t,0)＝[φ1(t,0) φ2(t,0)
…
φn(t,0)]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0068]
可选地，本实施例中，状态转移矩阵φ(t,0)可采用ode45计算，并将步长选择为t/1000，相对误差控制在1e-4；ode45为matlab提供的一种数值计算方法；在本发明其他的一些实施例中，基于实际的计算精度要求，也可以采用其他的数值计算方法、步长与误差控制设置完成状态转移矩阵的计算。
[0069]
定常矩阵获取步骤(步骤2)：根据floquet-lyapunov理论计算线性周期时变系统时不变化得到的定常矩阵q，其中q阵与步骤1所得的状态转移矩阵φ(t,0)满足式(21)所示关系：
[0070][0071]
状态转移矩阵φ(t,0)获取步骤(步骤3)：基于线性周期时变系统状态空间模型，选择n维单位矩阵i的n个列向量分别作为系统的初始状态，通过数值计算系统状态转移矩阵在一个周期t内的数值解，相应的流程图如图2所示。
[0072][0073]
时不变化变换矩阵获取步骤(步骤4)：根据线性周期时变系统解的结构，基于n组系统状态初值xk(0)(k＝1,2,
…
,n)，定常矩阵q与状态转移矩阵φ(t,0)求解时不变化变换矩阵p(t)：
[0074]
p(t)＝φ(t,0)[x1(0),x2(0),...,xn(0)]e-qt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0075]
振型信息获取步骤(步骤5)：根据时不变化变换矩阵p(t)的周期特性，对(25)计算得到的p(t)进行快速傅里叶变换(fast fourier transform,fft)，求取p(t)的傅里叶系数矩阵pn，本实施例中，快速傅里叶变换可采用matlab的fft函数实现。
[0076][0077]
将(26)回代入(3)中判断振型信息：
[0078][0079]
模态辨识步骤(步骤6)：通过比较矩阵pn中各元素p
n(rs)
模值大小进行系统模态辨识。
[0080]
相对振荡状态量划分步骤(步骤7)：通过比较pn中各元素的相位进行相对振荡状态量划分，对于不同子系统中具有相同物理意义的状态量，相位相差趋近0
°
被划分为共模振荡，相位相差趋近180
°
被划分为差模振荡。
[0081]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。