单自由度系统谐波平衡法装置

IPC分类号 : G06F17/14

申请号

CN201510295236.7

可选规格: 数量

库存1件

确认取消

￥30000; 库存1件

首页

立即咨询

看了又看

专利摘要

本发明公开了一种单自由度系统谐波平衡法装置，属于计算机工程应用领域。该装置包括非线性系统配置部件、符号变量定义部件、计算部件、输出部件。计算部件中设置了一个目标函数构造单元，用来动态构造优化求解均值、各次谐波的幅值与相位时需要的目标函数。该装置采用的谐波平衡法，仅考虑响应中的有限次谐波，将其分为低次谐波部分与高频谐波部分。从低频部分只包括均值与一次谐波开始优化、迭代提高求解精度，不断增加低频部分包含的最高次谐波的次数，同时减小高频谐波部分的最低次谐波项，以低次谐波部分计算得到的结果作为初始值，最终求解出响应中均值、有限次谐波的幅值与相位。本发明提供的装置具有很快的收敛速度与很高的精度。

权利要求

1.一种单自由度系统谐波平衡法装置，其特征在于，采用一种谐波平衡法，其在求解单自由度非线性系统的响应时，仅考虑响应中的有限次谐波，将所述有限次谐波分为低频谐波部分与高频谐波部分，

从所述低频部分只包括均值与一次谐波开始，通过优化方法求解均值、所述低频谐波部分的谐波幅值与相位，估计所述高频部分的谐波幅值与相位，并且通过迭代提高精度，

不断增加所述低频部分包含的最高次谐波的次数，每次增加一项，同时去掉高频谐波部分的最低次谐波项，以所述低频谐波部分计算得到的结果作为初始值，最终求解出响应中均值、所述有限次谐波的幅值与相位，

所述单自由度谐波平衡法装置包括以下部件：

非线性系统配置部件，用于定义非线性系统的表达式与激励，

符号变量定义部件，用于定义非线性系统响应中均值、所述有限次谐波的幅值与相位的符号变量，

计算部件，用于通过所述谐波平衡法计算响应中均值、所述有限次谐波的幅值与相位，

输出部件，用于存储所述计算部件得到的结果、以图像化的方式显示结果，

所述计算部件包括：计算模块准备模块、计算模块主体模块、计算模块后处理模块，

其中，所述计算部件准备模块，用于设置迭代的初始值、一个基波周期内时间的离散数目，

所述计算部件主体模块，用于通过迭代的方法计算响应中均值、所述有限次谐波的幅值与相位，以及评价参数，所述计算主体部件由两层循环组成，外层循环控制所述低频谐波部分最高次谐波的次数，内层循环通过优化方法求解响应中所述低频部分中均值、各次谐波的幅值与相位；

所述计算部件主体模块中设置了一个目标函数构造单元，用来构造优化求解所述低频谐波部分中均值、各次谐波的幅值与相位时用的目标函数，在迭代过程中，目标函数的表达式形式、待优化的所述低频部分中均值、各次谐波的幅值与相位的个数在动态变化，所述目标函数构造单元就是用来构造这个目标函数。

2. 根据权利要求1所述的自由度系统谐波平衡法装置，其特征在于，所述计算部件中，求解响应中均值、所述有限次谐波的幅值与相位采用的方法包括以下步骤：

S1：代表所述低频谐波部分的最高次谐波次数的循环变量k、内层循环变量i置初值：k=1，i=1，所述高频谐波部分置初值：，设置优化求解的初值X0；

S2：优化方法计算所述低频谐波部分中均值、各次谐波的幅值与相位，包括：

S21：计算所述低频谐波部分的表达式

(4)

式中，为均值、与分别为第n次谐波的幅值与相位，k为所述低频谐波部分中最高次谐波的次数；为基波的角频率，即一次谐波的角频率；

S22：计算非线性系统的响应y(t)的表达式(6)

(6)

S23：求方程误差序列的符号表达式，求其FFT的符号表达式

(8)

式中，为离散时间点，m=1、2、3、…、N，N为所方程误差序列离散的点数；

S24：求目标函数的表达式

(9)

S25：以式(9) 最小为优化目标，得到所述低频谐波部分中均值、各次谐波的幅值与相位的值；

S26：将所述低频谐波部分中的各次谐波的幅值转化为非负值，相位转化到区间上，方法是：

(10)

S3：更新迭代初值、所述高频谐波部分，包括：

S31：计算所述方程误差序列的数值表达式，计算方程误差均方根值；

S32：判断内层循环是否满足退出条件：所述方程误差均方根值减小到设定的精度，内层循环的迭代次数i达到设定的上限，

如果满足其中任一条件，则转S4增加所述低频谐波部分最高次谐波的次数k、继续迭代，

如果都不满足，则转S33；

S33：更新所述高频谐波部分，设置非线性优化的初值X0，循环变量增加1：i=i+1，转S22继续迭代，其中，用到如下公式与步骤：

S331：计算式(11)的具体数值，

(11)

S332：用FFT计算式(11)包含的谐波的幅值与相位，将其转化为余弦函数与正弦函数之和的形式，得到幅值与，求所述高频谐波部分的估计值的表达式(12)

(12)

S4：判断代表所述低频谐波部分的最高次谐波次数的循环变量k是否不大于设定的上限次数，或未达到精度条件，

如果是，则设置所述非线性优化的初值X0，循环变量k增加1：k=k+1，转步骤S2；

否则，计算结束。

3.根据权利要求2所述的自由度系统谐波平衡法装置，其特征在于，步骤S4中，上限次数不大于步骤S2中所述方程误差序列离散的点数N的一半N/2。

4.根据权利要求3所述的自由度系统谐波平衡法装置，其特征在于，步骤S4中，精度条件为外层循环当前循环，即第k轮循环，与上一轮循环得到的所述方程误差均方根值相等或之差小于设定的误差限。

5.根据权利要求4所述的自由度系统谐波平衡法装置，其特征在于，

步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1到i的这i轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位，

步骤S4中，所述非线性优化的初值X0包含2k+3个元素，为外层循环从1到k的这k轮循环中，第k次循环下，对应的内层循环中，第i次循环时，步骤S2与S3中得到的结果，其中2k+1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为步骤S3中FFT计算得到的对应谐波的幅值与相位。

6.根据权利要求4所述的自由度系统谐波平衡法装置，其特征在于，

步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1到i的这i轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位，

步骤S4中，所述非线性优化的初值X0包含2k+3个元素，为外层循环从1到k的这k轮循环中，第k次循环下，对应的内层循环从1到i的这i轮循环中，最小的误差均方根值对应的结果，其中2k+1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为步骤S3中FFT计算得到的对应谐波的幅值与相位。

7.根据权利要求4所述的自由度系统谐波平衡法装置，其特征在于，

步骤S33的详细步骤为：判断所述方程误差均方根值是否小于上一轮迭代得到的值，如果不小于，则非线性优化的初值X0在原来的基础上加上随机扰动，循环变量增加1：i=i+1，转S2继续迭代；否则，执行S331与S332，迭代初值x0设置为步骤S2中通过优化得到的均值、谐波幅值与相位，如果所述方程误差均方根值与上一轮迭代得到的误差均方根值之差小于设定的误差限，则在x0再加上一个随机扰动，

8.根据权利要求5、6或7所述自由度系统谐波平衡法装置，其特征在于，步骤S21~S24采用符号运算，步骤S3采用数值运算。

9.根据权利要求8所述的自由度系统谐波平衡法装置，其特征在于，所述符号变量定义部件中，当采用数学软件MATLAB实现时，代表均值、幅值与相位的符号变量由前缀与数字两部分组成，第一部分中均值、幅值符号变量的前缀相同，相位符号变量的前缀相同，第二部分数字为正整数数字，数字位数的确定方法是：设响应中所述有限次谐波的最高谐波次数为M，则第二部分数字的位数的最小值为，如果位数不够，则前面用数字0补足。

10.根据权利要求9所述的自由度系统谐波平衡法装置，其特征在于，所述目标函数构造单元采用的方法包括如下步骤：

a：用char函数将表达式(4)转换为字符串myObjFunStr1；

b：然后用字符串替换函数strrep将字符串myObjFunStr1的均值替换为字符串：x(1)，幅值替换为字符串：x(2)、x(3)、…、x(k)、x(k+1)，相位替换为字符串：x(k+2)、x(k+3)、…、x(2k)、x(2k+1)，得到代表目标函数的字符串型表达式myObjFunStr2；

c：在字符串myObjFunStr2前面加上字符串'(x)'，得到新的字符串myObjFunStr3；

d：用eval函数将字符串myObjFunStr3转换为目标函数。

说明书

技术领域

本发明属于计算机工程应用领域，具体而言，涉及一种单自由度系统谐波平衡法装置。

背景技术

谐波平衡法是处理非线性系统问题的一种非常有效的方法。很多实际的控制系统都满足应用这种方法的条件，因此在工程技术领域中获得了广泛的应用，例如齿轮传动、振动输送机、振动冷却机、振动磨机、发动机、汽车转向、电子电路等机器与系统的分析与设计中。其中单自由度非线性系统最为常见，也是多自由度非线性系统研究与应用的基础，得到了国内外专家学者、工程技术人员的广泛关注，并投入了大量的研究。它应用的基础及采用的处理方法是：将输出/响应表示为傅里叶级数的形式，仅考虑有限次谐波，通过数学方法求解谐波的幅值与相位，从而获得系统的频域特性。在非线性系统方程中代入输入/激励与输出/响应的傅里叶级数，比较两边系数来建立方程组，进而求解谐波的幅值与相位；但是当非线性系统变得复杂时，无法得到显式的方程组，此方法将面临着很多应用上得困难。直接应用快速傅里叶变换(FFT)，采用数值方法分析非常复杂的非线性系统，则具有很强的通用性。但当系统输入/激励包含的谐波项数或响应中要分析的谐波项数很多时，并考虑到系统输出的亚谐波与超谐波，则涉及到谐波幅值与相位等大量待解未知数，影响了分析的效率与精度，实现起来非常复杂、精度差，而且计算时间急剧增加，甚至无法有效地分析。

发明内容

针对谐波平衡法分析非线性系统时实现复杂、分析谐波项数少、耗时及精度低等问题，本发明的目的在于提供一种单自由度非线性系统谐波平衡法装置，以解决上述问题。

本发明的目的是这样实现的：一种单自由度系统谐波平衡法装置，采用一种谐波平衡法，其在求解单自由度非线性系统的响应时，仅考虑响应中的有限次谐波，将所述有限次谐波分为低频谐波部分与高频谐波部分，

所述单自由度谐波平衡法装置包括以下部件：

非线性系统配置部件，用于定义非线性系统的表达式与激励，

符号变量定义部件，用于定义非线性系统响应中均值、所述有限次谐波的幅值与相位的符号变量，

计算部件，用于通过所述谐波平衡法计算响应中均值、所述有限次谐波的幅值与相位，

输出部件，用于存储所述计算部件得到的结果、以图像化的方式显示结果，

所述计算部件包括：计算模块准备模块、计算模块主体模块、计算模块后处理模块，

其中，所述计算部件准备模块，用于设置迭代的初始值、一个基波周期内时间的离散数目，

进一步地，所述计算部件中，求解响应中均值、所述有限次谐波的幅值与相位采用的方法包括以下步骤：

S1：代表所述低频谐波部分的最高次谐波次数的循环变量k、内层循环变量i置初值：k=1，i=1，所述高频谐波部分置初值：，设置优化求解的初值X0；

S2：优化方法计算所述低频谐波部分中均值、各次谐波的幅值与相位，包括：

S21：计算所述低频谐波部分的表达式

(4)

式中，为均值、与分别为第n次谐波的幅值与相位，k为所述低频谐波部分中最高次谐波的次数；为基波的角频率，即一次谐波的角频率，

S22：计算非线性系统的响应y(t)的表达式(6)

(6)

S23：求方程误差序列的符号表达式，求其FFT的符号表达式

(8)

式中，为离散时间点，m=1、2、3、…、N，N为所方程误差序列离散的点数；

S24：求目标函数的表达式

(9)

S25：以式(9) 最小为优化目标，得到所述低频谐波部分中均值、各次谐波的幅值与相位的值；

S26：将所述低频谐波部分中的各次谐波的幅值转化为非负值，相位转化到区间上，方法是：

(10)

S3：更新迭代初值、高频谐波部分，包括：

S31：计算所述方程误差序列的数值表达式，计算方程误差均方根值；

S32：判断内层循环是否满足退出条件：所述方程误差均方根值减小到设定的精度，内层循环的迭代次数i达到设定的上限，

如果满足其中任一条件，则转S4增加所述低频谐波部分最高次谐波的次数k、继续迭代，

如果都不满足，则转S33；

S33：更新所述高频谐波部分，设置非线性优化的初值X0，循环变量增加1，i=i+1，转S22继续迭代，其中，用到如下公式与步骤：

S331：计算式(11)的具体数值，

(11)

S332：用FFT计算式(11)包含的谐波的幅值与相位，将其转化为余弦函数与正弦函数之和的形式，得到幅值与，求所述高频谐波部分的估计值的表达式(12)

(12)

S4：判断代表所述低频谐波部分的最高次谐波次数的循环变量k是否不大于设定的上限次数，或未达到精度条件，

如果是，则设置所述非线性优化的初值X0，循环变量k增加1：k=k+1，转步骤S2；

否则，计算结束。

进一步地，步骤S4中，上限次数不大于步骤S2中所述方程误差序列离散的点数N的一半N/2。

进一步地，步骤S4中，精度条件为外层循环当前循环，即第k轮循环，与上一轮循环得到的所述方程误差均方根值相等或之差小于设定的误差限。

在一些实施方式中，步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1到i的这i轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位，

在一些实施方式中，步骤S33的详细步骤为：判断所述方程误差均方根值是否小于上一轮迭代得到的值，如果不小于，则非线性优化的初值X0在原来的基础上加上随机扰动，循环变量增加1：i=i+1，转S2继续迭代；否则，执行S331与S332，迭代初值x0设置为步骤S2中通过优化得到的均值、谐波幅值与相位，如果所述方程误差均方根值与上一轮迭代得到的误差均方根值之差小于设定的误差限，则在x0再加上一个随机扰动，

进一步地，步骤S21~S24采用符号运算，步骤S3采用数值运算。

更进一步地，所述符号变量定义部件中，当采用数学软件MATLAB实现时，代表均值、幅值与相位的符号变量由前缀与数字两部分组成，第一部分中均值、幅值符号变量的前缀相同，相位符号变量的前缀相同，第二部分数字为正整数数字，数字位数的确定方法是：设响应中所述有限次谐波的最高谐波次数为M，则第二部分数字的位数的最小值为，如果位数不够，则前面用数字0补足。

再进一步地，所述目标函数构造单元采用的方法包括如下步骤：

a：用char函数将表达式(4)转换为字符串myObjFunStr1；

c：在字符串myObjFunStr2前面加上字符串'(x)'，得到新的字符串myObjFunStr3；

d：用eval函数将字符串myObjFunStr3转换为目标函数。

本发明与现有技术相比具有显著的优点和有益效果：

(1)在本发明中，采用符号运算与数值运算相结合，计算精度高、效率高。符号运算保证了精度，数值运算提高运行效率、减少计算时间。

(2)本发明采用的方法将非线性系统的响应中考虑的有限次谐波分为低频谐波部分与高频谐波部分，从低频谐波部分只包含均值与一次谐波开始，只有三个待优化的未知量，非线性优化简单易行；逐次增加低频谐波部分包含谐波的项数，每次增加一个高次谐波，引入二个未知量，非线性优化从简单到复杂，保证了优化精度与效率。

(3)本发明提供的装置中设置了一个目标函数构造单元，使非线性目标函数始终只包含待优化的未知量，即每次优化只有2k+1个未知量，没有冗余的未知量，提高了优化效率、大大减小了运行时间。

(4)步骤S3中，初值X0设置时加上随机扰动，避免了后续优化迭代时重复已经历的计算过程，并且运用随机化增加了迭代收敛的速度，能够显著提高装置的运算精度、减少整体迭代的次数。

(5)步骤S4中，初值X0设置时，对于新增的二个待优化的未知量采用步骤S3中FFT计算得到的对应次谐波的幅值与相位，在没有先验信息的情况下，有充分的理由认为其最接近真实值，因此可以提高下一轮非线性优化的效率，提高计算精度。

(6)代表均值、幅值与相位的符号变量第二部分数字由位数规定的正整数数字组成，并且位数与要考虑的有限谐波次数关联，保证了目标函数构造单元实现的可读性、简洁性与方便性，在使用本发明提供的装置时，只需要改变最高谐波次数为M的值即可，方便使用。

附图说明

图1为本发明提供的装置的结构组成图；

图2为本发明装置采用的谐波平衡法的原理图；

图3为本发明装置采用的谐波平衡法的流程图；

图4为目标函数构造单元方法图；

图5为本发明实施例1的均值与谐波幅值图；

图6为本发明实施例2的均值与谐波幅值图；

图7为本发明实施例3的方程误差均方根值随迭代次数变化图；

图8为本发明实施例3的均值与谐波幅值图。

具体实施方式

以下结合附图及较佳实施例，对依据本发明提供的结构特征、具体实施方式及其功效，详细说明如下。

本发明提供了一种单自由度系统谐波平衡法装置，结构如图1所示。此装置适用的单自由度非线性系统可表示为如下方程

(1)

式中，t为时间； y(t)为系统的响应，即输出，下面为了表述方便，有些地地方省略了(t)；、分别为响应的一阶、二阶导数；为关于和的函数； f(t)为激励，即输入。

用谐波平衡法求解非线性系统的响应时，将输出/响应表示为傅里叶级数的形式，仅考虑有限次谐波，将响应中最高次谐波的次数设为一固定值M。

则输出表示为

(2)

式中，为均值，与分别为M个谐波的幅值与相位，为基波的角频率。

如果y不是精确解，则方程(1)等号左边与右边的就有一个误差e(t)

(3)。

谐波平衡法装置需要求解式(2)中的均值、M个谐波幅值、M个谐波相位，一共2M+1个未知量。本发明提供的所述单自由度系统谐波平衡法装置采用如图2所示的方法，将响应式(2)中的所有谐波分为所述低频谐波部分与所述高频谐波部分。所述低频谐波部分包含的谐波项数，也就是最高次谐波的次数，用循环变量k来控制。

所述低频谐波部分表示为

(4)

式中，为均值；与分别为第n次谐波的幅值与相位；k为所述低频谐波部分中最高次谐波的次数；为基波的角频率，即一次谐波的角频率。

式(4)所描述的所述低频谐波部分对应着1个均值、k个谐波幅值、k个相位，所以所述低频谐波部分一共包含2k+1个未知量。

则所述高频谐波部分表示为

(5)。

式(5)所描述的所述高次次谐波部分包含的最低次谐波的次数为k+1，则所述高频谐波部分包含的谐波项数为M-k。这M-k项谐波对应于M-k个幅值与M-k个相位，所以所述高频谐波部分一共包含2(M-k) 个未知量。

开始时，k=1，所述低频谐波部分只包含均值与一次谐波，此时有均值、一次谐波的幅值与相位3个未知量，通过非线性优化求这3个未知量，优化的初值包含三个元素，可以随机生成或根据非线性系统的先验知识设定。用式(3)计算方程误差的具体数值，通过FFT计算得到所述高频谐波部分的幅值与相位的估计值，一共2(M-1)个量。根据上述结果，反复迭代，提高均值、一次谐波幅值与相位3个未知量的精度。

不断增加所述低频谐波部分包含的最高次谐波的次数，k值每次增加1。也就是说低频谐波部分每次增加一项，同时去掉所述高频谐波部分的最低次谐波项，即k+1次谐波项。以外层循环为1，2，3，…，k时的迭代结果作为第k+1次迭代的初始值，反复迭代，最终求解出响应中均值、M项谐波的幅值与相位。

如图1所示，所述单自由度系统谐波平衡法装置包括以下部件：

非线性系统配置部件，用于定义非线性系统的表达式与激励，

符号变量定义部件，用于定义非线性系统响应中均值、所述有限次谐波的幅值与相位的符号变量，

计算部件，用于通过所述谐波平衡法计算响应中均值、所述有限次谐波的幅值与相位，

输出部件，用于存储所述计算部件得到的结果、以图像化的方式显示结果，

所述计算部件包括：计算模块准备模块、计算模块主体模块、计算模块后处理模块，

其中，所述计算部件准备模块，用于设置迭代的初始值、一个基波周期内时间的离散数目，

所述计算部件主体模块，用于通过迭代的方法计算响应中均值、所述有限次谐波的幅值与相位，以及评价参数，所述计算主体部件由两层循环组成，外层循环控制所述低频谐波部分最高次谐波的次数，内层循环通过优化方法求解响应中所述低频谐波部分中均值、各次次谐波的幅值与相位；

所述计算部件主体模块中设置了一个目标函数构造单元，用来构造优化求解所述低频谐波部分中均值、各次谐波的幅值与相位时用的目标函数，在迭代过程中，目标函数的表达式形式、待优化的所述低频谐波部分中均值、各次谐波的幅值与相位的个数2k+1在动态变化，所述目标函数构造单元就是用来构造这个目标函数。

具体来说，本发明提供的所述单自由度系统谐波平衡法装置主要通过两层循环来实现的。外层循环变量k控制所述低频谐波部分包含的谐波的项数；内存循环变量i控制在k一定的情况下，通过迭代提高所述低频谐波部分包含的均值、谐波的相位与幅值的精度，同时也提高所述高频谐波部分包含的谐波的幅值与相位估计值的精度。具体来所，主要包括以下步骤。

S1：代表所述低频谐波部分的最高次谐波次数的循环变量k、内层循环变量i置初值：k=1，i=1，所述高频谐波部分置初值：，设置优化求解的初值X0，初始X0可随机生成，也可根据非线系统的先验知识确定。

S2：通过优化方法计算所述低频谐波部分中均值、各次谐波的幅值与相位，包括：

S21：计算所述低频谐波部分的表达式(4)；

S22：计算非线性系统的响应y(t)的表达式(6)；

(6)

S23：并用式(3)求离散型方程误差序列

(7)

式中，为离散时间点，m=1、2、3、…、N，N为所方程误差序列离散的点数，

并求方程误差序列快速傅里叶变换的表达式，并求FFT的表达式

(8)

S24：求非线性目标函数的表达式

(9)

S25：以式(9)最小为优化目标，得到所述低频谐波部分中均值、各次谐波的幅值与相位的值；

S26：将所述低频谐波部分中的谐波幅值转化为非负值，相位转化到区间上，方法是：

(10)。

S3：更新迭代初值、高频谐波部分，包括：

S31：计算所述方程误差序列，计算其所述误差均方根值；

S32：判断内层循环是否满足退出条件：所述误差均方根值减小到设定的下限，迭代次数i达到设定的上限，

如果满足其中一个条件，则转S4增加第次谐波的次数k、继续迭代，

如果都不满足，则转S33；

S33：更新所述高频谐波部分，S35：设置非线性优化的初值X0，循环变量增加1：i=i+1，转S22继续迭代，其中，用到如下公式与步骤：

S331：计算式(11)的值，

(11)

S332：则用FFT计算式(11)包含的谐波的幅值与相位，将其转化为余弦函数与正弦函数之和的形式，得到幅值与，求高频谐波部分的估计值的表达式

(12)。

S4：判断表所述低频谐波部分的最高次谐波次数的循环变量k是否不大于设定的上限次数，或未达到精度条件，

如果是，则设置所述非线性优化的初值X0，

循环变量k增加1：k=k+1，转步骤S2；否则，计算结束。

步骤S4中，上限次数不大于步骤S2中所述方程误差序列离散的点数N的一半N/2。

随着k值的增加，所述自由度系统谐波平衡法装置得到的结果精度越来越高。当k增加一定程度时，所述误差均方根值保持不变或减小的很慢，这时，可以认为迭代终止。

步骤S3和S4分别为内层与外层循环的下一轮循环设置迭代初值。如图2所示，步骤S3中X0有2k+1个元素，对应于均值、k项低次谐波的幅值与相位，为内层的下一轮循环准备优化初值；步骤S4中X0有2k+3个元素，对应于均值、k+1项谐波的幅值与相位，为外层的下一轮循环准备优化初值。外层循环增加1，意味着所述低频谐波部分加进去一项谐波。步骤S4中的X0比步骤S3中的X0多出的2个元素为所述低频谐波部分新增的这一项谐波对应的幅值与相位。对于步骤S3与S4中迭代初值X0的设定，有以下3种都切实可行、效果良好的技术方案。

技术方案1：步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1到i的这i轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位。

技术方案2：步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1到i的这i轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位。

技术方案3：步骤S33的详细步骤为：判断所述方程误差均方根值是否小于上一轮迭代得到的值，如果不小于，则非线性优化的初值X0在原来的基础上加上随机扰动，循环变量增加1：i=i+1，转S2继续迭代；否则，执行S331与S332，迭代初值x0设置为步骤S2中通过优化得到的均值、谐波幅值与相位，如果所述方程误差均方根值与上一轮迭代得到的误差均方根值之差小于设定的误差限，则在x0再加上一个随机扰动，可用随机函数随机生成。

在步骤S2中，式(8)与(9)中都包含了2k+1为待求未知量，而这些式子本身的形式随着非线性系统方程(1)、M和k的变化而变化，并且差异巨大；因此，在步骤S21~S24采用符号运算，可充分利用工程软件或数学软件的强大的符号运算功能、求导等数学运算功能，极大地减小了编程的复杂度。步骤S3是为了估计所述高频谐波部分的各次谐波的幅值与相位，采用数值运算可以极大地提高运算速度。

本发明提供的单自由度系统谐波平衡法装置可以在很多高级编程语言、工程软件、数学软件中实施，如MATLAB。具体实施的过程中，式(2)中的均值、谐波幅值与相位一共2M+1个待求解的量，需要定义对应的变量来处理。为了构造目标函数的方便，代表均值、幅值与相位的符号变量由前缀与数字两部分组成。均值、幅值符号变量的前缀相同，相位符号变量的前缀相同；第二部分数字为正整数数字，数字位数的确定方法是：设需要分析的响应中最高谐波的次数为M，则第二部分数字的位数的最小值为，如果位数不够，则高位由数字0补足。

例如，均值、幅值符号变量的前缀为An，相位符号变量的前缀为Phin。当M<100时，，则在MATLAB中定义符号变量的方法如下：

for myi=1:(2*M+1)

eval(['syms ','An',num2str(myi,'%02d')])

eval(['syms ','Phin',num2str(myi,'%02d')])

An(myi)=eval(['An',num2str(myi,'%02d')]);

Phin(myi)=sym(['Phin',num2str(myi,'%02d')]);

end。

如果M=12，则上面方法定义了13个代表均值与谐波幅值的符号变量An01、An02、An03、…、An11、An12、An13，还定义了13个代表与上述幅值对应的相位的符号变量Phin01、Phin02、Phin03、…、Phin11、Phin12、Phin13，并且还定义了两个组号数组An与Phin，都包含了13个元素，分别为幅值及相位的符号变量。因为An01与均值对应，符号变量Phin01只是为了配合程序的方便，使对应次谐波的幅值与相位变量下标一致而定义的，在此并不使用。

基于上面定义的符号变量，所述目标函数构造单元采用的方法包括如下步骤：

a：用char函数将表达式(4)转换为字符串myObjFunStr1；

b：然后用字符串替换函数strrep将字符串myObjFunStr1的均值替换为字符串：x(1)，幅值替换为字符串：x(2) 、x(3) 、…、x(k)、x(k+1)，相位替换为字符串：x(k+2) 、x(k+3) 、…、x(2k)、x(2k+1)，得到代表目标函数的字符串型表达式myObjFunStr2；

c：在字符串myObjFunStr2前面加上字符串'(x)'，得到新的字符串myObjFunStr3；

d：用eval函数将字符串myObjFunStr3转换为目标函数。

步骤b主要采用strrep函数，具体方法是：

myObjFunStr2=strrep(myObjFunStr1,'An01','x(1)')

for tempk=1:k

myObjFunStr2=strrep(myObjFunStr2,['An',num2str(tempk+1,'%02d')],

['x(',num2str(tempk+1,'%d'),')'])

myObjFunStr2=strrep(myObjFunStr2,['Phin',num2str(tempk+1,'%02d')],

['x(',num2str(tempk+k+1,'%d'),')'])

end。

代表均值、幅值与相位第二部分数字的位数的最小值为，高位补0是为了保证在strrep函数替换的过程中不出现错误替换。如果不固定位数，则会出现错误。例如：若将变量定义为An1、An2、An3、…、An11、An12、An13，则函数strrep在A1替换为x(1)时，会将An10、An11、An12、An13错误替换成x(1)0、x(1)1、x(1)2、x(1)3。

在数学软件MATLAB中实施本发明。下面以3个实施例来说明本发明的有益效果。

实施例1：单自由度非线性系统的方程为

。

假设响应中一次谐波的角频率为1，设置迭代初值x0=[0,0,0]，仅考虑32次以下的谐波，即M=32，内层循环的最大次数为5，采样点数N=64，步骤S4中误差限设为，随机扰动的范围为。实施技术方案3。随着k值的增加，所述方程误差均方根值迅速减小，前6次循环时的方程误差均方值依次为：

5.656854249492381e+000，2.496859179098087e-003，8.156122475628342e-008，

4.360844843951736e-012，4.360844843951736e-012，4.360844843951736e-012。

各次谐波的幅值与相位如图5所示，其中0次谐波对应于均值，1~6次谐波由步骤S2的非线性优化得到，7~32次谐波由步骤S3中的数值FFT得到。

从上述的数据可以明显看出，当k>4以后，迭代结果趋于恒定，精度级别为10e-12，即。

实施例2：将实施例1中的激励修改为cos2t，假设响应中一次谐波的角频率为2，设置迭代初值x0为随机生成，其他参数同实施例1。实施技术方案3。各次谐波的幅值与相位如图6所示，其中0次谐波对应均值，1~6次谐波由步骤S2的非线性优化得到，7~32次谐波由步骤S3中的数值FFT得到。当k>2以后，迭代结果趋于恒定，精度级别为。

实施例3：将实施例1中的激励修改为cos3t，假设响应中一次谐波的角频率为3，设置迭代初值x0为随机生成，仅考虑32次以下的谐波，即M=32，内层循环的最大次数为4，采样点数N=64，。实施技术方案3。随着k值的增加，所述方程误差均方根值迅速减小，如图7所示，迭代结果精度级别为。各次谐波的幅值与相位如图8所示，其中0次谐波对应均值，1~6次谐波由步骤S2的非线性优化得到，7~32次谐波由步骤S3中的数值FFT得到。

本发明提供的装置只需要改变M的值即可控制谐波的次数，应用方便。上述的3个实施例的结果，说明了本发明提供的装置具有很快的收敛速度与很高的精度。

如果要分析亚谐波，只需要根据要亚谐波的角频率来设置一次谐波的频率即可。

也可以将目标函数改为方程误差的谐波均方根值，同样可达到很好的效果。

需要注意的是，上述具体实施例仅仅是示例性的，在本发明的上述教导下，本领域技术人员可以在上述实施例的基础上进行各种改进和变形，而这些改进或者变形均落在本发明的保护范围内。本领域技术人员应该明白，上面的具体描述只是为了解释本发明的目的，并非用于限制本发明。本发明的保护范围由权利要求及其等同物限定。

单自由度系统谐波平衡法装置专利购买费用说明