单自由度系统谐波平衡法装置的制造方法

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

【发明内容】

[0003] 针对谐波平衡法分析非线性系统时实现复杂、分析谐波项数少、耗时及精度低等 问题，本发明的目的在于提供一种单自由度非线性系统谐波平衡法装置，W解决上述问题。
[0004] 本发明的目的是该样实现的；一种单自由度系统谐波平衡法装置，采用一种谐波 平衡法，其在求解单自由度非线性系统的响应时，仅考虑响应中的有限次谐波，将所述有限 次谐波分为低频谐波部分与高频谐波部分， 从所述低频部分只包括均值与一次谐波开始，通过优化方法求解均值、所述低频谐波 部分的谐波幅值与相位，估计所述高频部分的谐波幅值与相位，并且通过迭代提高精度， 不断增加所述低频部分包含的最高次谐波的次数，每次增加一项，同时去掉高频谐波 部分的最低次谐波项，W所述低频谐波部分计算得到的结果作为初始值，最终求解出响应 中均值、所述有限次谐波的幅值与相位， 所述单自由度谐波平衡法装置包括W下部件： 非线性系统配置部件，用于定义非线性系统的表达式与激励， 符号变量定义部件，用于定义非线性系统响应中均值、所述有限次谐波的幅值与相位 的符号变量， 计算部件，用于通过所述谐波平衡法计算响应中均值、所述有限次谐波的幅值与相位， 输出部件，用于存储所述计算部件得到的结果、W图像化的方式显示结果， 所述计算部件包括；计算模块准备模块、计算模块主体模块、计算模块后处理模块， 其中，所述计算部件准备模块，用于设置迭代的初始值、一个基波周期内时间的离散数 目， 所述计算部件主体模块，用于通过迭代的方法计算响应中均值、所述有限次谐波的幅 值与相位，W及评价参数，所述计算主体部件由两层循环组成，外层循环控制所述低频谐波 部分最高次谐波的次数，内层循环通过优化方法求解响应中所述低频部分中均值、各次谐 波的幅值与相位； 所述计算部件主体模块中设置了一个目标函数构造单元，用来构造优化求解所述低频 谐波部分中均值、各次谐波的幅值与相位时用的目标函数，在迭代过程中，目标函数的表达 式形式、待优化的所述低频部分中均值、各次谐波的幅值与相位的个数在动态变化，所述目 标函数构造单元就是用来构造该个目标函数。
[0005]进一步地，所述计算部件中，求解响应中均值、所述有限次谐波的幅值与相位采用 的方法包括W下步骤： 51 ;代表所述低频谐波部分的最高次谐波次数的循环变量A、内层循环变量i置初值： 仁1，i=l，所述高频谐波部分置初值设置优化求解的初值X0 ; 52 ;优化方法计算所述低频谐波部分中均值、各次谐波的幅值与相位，包括： 521 ;计算所述低频谐波部分乃；p(f)的表达式
(4) 式中，4 .''2为均值、4|与巧t分别为第。次谐波的幅值与相化幼所述低频谐波部 分中最高次谐波的次数；巧为基波的角频率，即一次谐波的角频率， 522 ;计算非线性系统的响应_7(幻的表达式化）
(6) 523 ;求方程误差序列的符号表达式，求其FFT的符号表达式
(8) 式中，为离散时间点，^1、2、3、…、化A为所方程误差序列离散的点数； 524 ;求目标函数的表达式
(9) S25;W式巧）最小为优化目标，得到所述低频谐波部分中均值、各次谐波的幅值与相 位為，-?与，巧,A，巧，…，為，巧^的值； S26;将所述低频谐波部分中的各次谐波的幅值-韦為,转化为非负值，相位 (10) 巧，巧:,…，巧转化到区间[O,2?!)上，方法是； 53 ;更新迭代初值、高频谐波部分，包括； 531 ;计算所述方程误差序列句的数值表达式，计算方程误差均方根值e。,。； 532 ;判断内层循环是否满足退出条件：I所述方程误差均方根值e,,。、减小到设定的 精度，{數内层循环的迭代次数i达到设定的上限， 如果满足其中任一条件，则转S4增加所述低频谐波部分最高次谐波的次数A、继续迭 代， 如果都不满足，则转S33; 533 ;更新所述高频谐波部分，设置非线性优化的初值X0,循环变量增加1，i=i+l，转 S22继续迭代，其中，用到如下公式与步骤： 5331 ;计算式（11)的具体数值，
(11) 5332 ;用FFT计算式（11)包含的谐波的幅值与相位，将其转化为余弦函数与正弦函数 之和的形式，得到幅值0,,与，求所述高频谐波部分的估计值的表达式（12)
(12) 54 ;判断代表所述低频谐波部分的最高次谐波次数的循环变量A是否不大于设定的上 限次数，或未达到精度条件， 如果是，则设置所述非线性优化的初值X0,循环变量姑曽加1 ;仁知1，转步骤S2 ; 否则，计算结束。
[0006] 进一步地，步骤S4中，上限次数不大于步骤S2中所述方程误差序列离散的点数# 的一半￥2。
[0007] 进一步地，步骤S4中，精度条件为外层循环当前循环，即第每％循环，与上一轮循 环得到的所述方程误差均方根值相等或之差小于设定的误差限。
[0008] 在一些实施方式中，步骤S3中，所述非线性优化的初值X0的设置方法为；内层循 环从1到i的该J轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相 位， 步骤S4中，所述非线性优化的初值X0包含2知3个元素，为外层循环从1到铺勺该A轮循环中，第巧欠循环下，对应的内层循环中，第i次循环时，步骤S2与S3中得到的结果， 其中2知1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为步骤S3中FFT计算得到的对应谐波的幅值与相位。
[0009]在一些实施方式中，步骤S3中，所述非线性优化的初值X0的设置方法为；内层循 环从1到i的该J轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相 位， 步骤S4中，所述非线性优化的初值X0包含2知3个元素，为外层循环从1到铺勺该A轮循环中，第域C循环下，对应的内层循环从1到i的该J轮循环中，最小的误差均方根值 €。,;、对应的结果，其中2知1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另 外2个元素为步骤S3中FFT计算得到的对应谐波的幅值与相位。
[0010]在一些实施方式中，步骤S33的详细步骤为；判断所述方程误差均方根值是否 小于上一轮迭代得到的值，如果不小于，则非线性优化的初值X0在原来的基础上加上随机 扰动，循环变量增加1 ;i=i+l，转S2继续迭代；否则，执行S331与S332,迭代初值xO设置为 步骤S2中通过优化得到的均值、谐波幅值与相位，如果所述方程误差均方根值与上一 轮迭代得到的误差均方根值之差小于设定的误差限，则在xO再加上一个随机扰动， 步骤S4中，所述非线性优化的初值X0包含2知3个元素，为外层循环从1到铺勺该A轮循环中，第巧欠循环下，对应的内层循环中，第i次循环时，步骤S2与S3中得到的结果， 其中2知1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为步骤 S3中FFT计算得到的对应谐波的幅值与相位。
[0011] 进一步地，步骤S21~S24采用符号运算，步骤S3采用数值运算。
[0012]更进一步地，所述符号变量定义部件中，当采用数学软件MTLAB实现时，代表 均值、幅值与相位的符号变量由前缀与数字两部分组成，第一部分中均值、幅值符号变量 的前缀相同，相位符号变量的前缀相同，第二部分数字为正整数数字，数字位数的确定方 法是；设响应中所述有限次谐波的最高谐波次数为必则第二部分数字的位数的最小值为 _i〇ga./)j+1，如果位数不够，则前面用数字0补足。
[0013]再进一步地，所述目标函数构造单元采用的方法包括如下步骤：a;用char函数将表达式（4)转换为字符串myObjRmStrl; b;然后用字符串替换函数strrep将字符串myWj化nStrl的均值替换为字符串；x(l)，幅值替换为字符串；x(2)、x(3)、…、x(0、x(知1)，相位替换为字符串；x(知2)、 x(左+3)、???、x(ik)、x(ik+l)，得到代表目标函数的字符串型表达式myObjF'unStrS; C;在字符串mW)bj化nStr2前面加上字符串' @(x)'，得到新的字符串mW)bj化nStr3 ;d;用eval函数将字符串mW)bj化nStr3转换为目标函数。
[0014]本发明与现有技术相比具有显著的优点和有益效果： (1)在本发明中，采用符号运算与数值运算相结合，计算精度高、效率高。符号运算保证 了精度，数值运算提高运行效率、减少计算时间。
[0015](2)本发明采用的方法将非线性系统的响应中考虑的有限次谐波分为低频谐波 部分与高频谐波部分，从低频谐波部分只包含均值与一次谐波开始，只有=个待优化的未 知量，非线性优化简单易行；逐次增加低频谐波部分包含谐波的项数，每次增加一个高次谐 波，引入二个未知量，非线性优化从简单到复杂，保证了优化精度与效率。
[0016] (3)本发明提供的装置中设置了一个目标函数构造单元，使非线性目标函数始终 只包含待优化的未知量，即每次优化只有2知1个未知量，没有冗余的未知量，提高了优化 效率、大大减小了运行时间。
[0017] (4)步骤S3中，初值X0设置时加上随机扰动，避免了后续优化迭代时重复已经历 的计算过程，并且运用随机化增加了迭代收敛的速度，能够显著提高装置的运算精度、减少 整体迭代的次数。
[0018] (5)步骤S4中，初值X0设置时，对于新增的二个待优化的未知量采用步骤S3中 FFT计算得到的对应次谐波的幅值与相位，在没有先验信息的情况下，有充分的理由认为其 最接近真实值，因此可W提高下一轮非线性优化的效率，提高计算精度。
[0019] (6)代表均值、幅值与相位的符号变量第二部分数字由位数规定的正整数数字组 成，并且位数与要考虑的有限谐波次数关联，保证了目标函数构造单元实现的可读性、简洁 性与方便性，在使用本发明提供的装置时，只需要改变最高谐波次数为础勺值即可，方便使 用。
【附图说明】
[0020] 图1为本发明提供的装置的结构组成图； 图2为本发明装置采用的谐波平衡法的原理图； 图3为本发明装置采用的谐波平衡法的流程图； 图4为目标函数构造单元方法图； 图5为本发明实施例1的均值与谐波幅值图； 图6为本发明实施例2的均值与谐波幅值图； 图7为本发明实施例3的方程误差均方根值随迭代次数变化图； 图8为本发明实施例3的均值与谐波幅值图。
【具体实施方式】
[0021] W下结合附图及较佳实施例，对依据本发明提供的结构特征、【具体实施方式】及其 功效，详细说明如下。
[0022] 本发明提供了一种单自由度系统谐波平衡法装置，结构如图1所示。此装置适用 的单自由度非线性系统可表示为如下方程
(1) 式中，巧3时间；_fU)为系统的响应，即输出，下面为了表述方便，有些地地方省略了U);却1、^/)分别为响应的一阶、二阶导数；F(_v(/)J.'(〇)为关于抑:1和>如）的函 数；/U)为激励，即输入。
[0023] 用谐波平衡法求解非线性系统的响应时，将输出/响应表示为傅里叶级数的形 式，仅考虑有限次谐波，将响应中最高次谐波的次数设为一固定值#。
[0024] 则输出表示为 (2) 式中，4:,:' 2为均值，4- .1:-. 4.….?义f..,，.'k;与私梦;、斯…、斯fr斯f分别为肿 谐波的幅值与相位，巧为基波的角频率。
[00巧]如果不是精确解，则方程（1)等号左边与右边的就有一个误差eU)
(3)。
[0026] 谐波平衡法装置需要求解式（2)中的均值、必个谐波幅值、必个谐波相位，一共2iM 个未知量。本发明提供的所述单自由度系统谐波平衡法装置采用如图2所示的方法，将响 应式（2)中的所有谐波分为所述低频谐波部分与所述高频谐波部分。所述低频谐波部分包 含的谐波项数，也就是最高次谐波的次数，用循环变量匆5控制。
[0027] 所述低频谐波部分表示为
(4) 式中，为均值；4,与巧，分别为第"次谐波的幅值与相位；巧/所述低频谐波部分 中最高次谐波的次数；为基波的角频率，即一次谐波的角频率。
[002引式（4)所描述的所述低频谐波部分对应着1个均值、外谐波幅值、外相化所W所述低频谐波部分一共包含2知1个未知量。
[0029] 则所述高频谐波部分表示为
巧)。
[0030] 式（5)所描述的所述高次次谐波部分包含的最低次谐波的次数为知1，则所述高 频谐波部分包含的谐波项数为i^A。该i^仅m皆波对应于i^A个幅值与i^A个相位，所W所 述高频谐波部分一共包含2(i^0个未知量。
[003。开始时，仁1，所述低频谐波部分只包含均值与一次谐波，此时有均值、一次谐波的 幅值与相位3个未知量，通过非线性优化求该3个未知量，优化的初值包含=个元素，可W随机生成或根据非线性系统的先验知识设定。用式（3)计算方程误差的具体数值，通过FFT 计算得到所述高频谐波部分的幅值与相位的估计值，一共2 (JM)个量。根据上述结果，反 复迭代，提高均值、一次谐波幅值与相位3个未知量的精度。
[0032] 不断增加所述低频谐波部分包含的最高次谐波的次数，對I每次增加1。也就是说 低频谐波部分每次增加一项，同时去掉所述高频谐波部分的最低次谐波项，即知1次谐波 项。W外层循环为1，2,3,…，如寸的迭代结果作为第知1次迭代的初始值，反复迭代，最终 求解出响应中均值、皆波的幅值与相位。
[0033] 如图1所示，所述单自由度系统谐波平衡法装置包括W下部件： 非线性系统配置部件，用于定义非线性系统的表达式与激励， 符号变量定义部件，用于定义非线性系统响应中均值、所述有限次谐波的幅值与相位 的符号变量， 计算部件，用于通过所述谐波平衡法计算响应中均值、所述有限次谐波的幅值与相位， 输出部件，用于存储所述计算部件得到的结果、W图像化的方式显示结果， 所述计算部件包括；计算模块准备模块、计算模块主体模块、计算模块后处理模块， 其中，所述计算部件准备模块，用于设置迭代的初始值、一个基波周期内时间的离散数 目， 所述计算部件主体模块，用于通过迭代的方法计算响应中均值、所述有限次谐波的幅 值与相位，W及评价参数，所述计算主体部件由两层循环组成，外层循环控制所述低频谐波 部分最高次谐波的次数，内层循环通过优化方法求解响应中所述低频谐波部分中均值、各 次次谐波的幅值与相位； 所述计算部件主体模块中设置了一个目标函数构造单元，用来构造优化求解所述低频 谐波部分中均值、各次谐波的幅值与相位时用的目标函数，在迭代过程中，目标函数的表达 式形式、待优化的所述低频谐波部分中均值、各次谐波的幅值与相位的个数2知1在动态变 化，所述目标函数构造单元就是用来构造该个目标函数。
[0034]具体来说，本发明提供的所述单自由度系统谐波平衡法装置主要通过两层循环来 实现的。外层循环变量雄!制所述低频谐波部分包含的谐波的项数；内存循环变量i控制 在^定的情况下，通过迭代提高所述低频谐波部分包含的均值、谐波的相位与幅值的精 度，同时也提高所述高频谐波部分包含的谐波的幅值与相位估计值的精度。具体来所，主要 包括W下步骤。
[00巧]S1 ;代表所述低频谐波部分的最高次谐波次数的循环变量A、内层循环变量i置初 值：仁1，i=l，所述高频谐波部分九",(〇置初值；扣,0，设置优化求解的初值X0，初始X0 可随机生成，也可根据非线系统的先验知识确定。
[0036]S2;通过优化方法计算所述低频谐波部分中均值、各次谐波的幅值与相位，包括： 521 ;计算所述低频谐波部分ia"的表达式（4); 522 ;计算非线性系统的响应_7(幻的表达式化）； -1'(/) = _1'自',:(/) + _1'能（0 化） 523 ;并用式（3)求离散型方程误差序列 地,,,)=扣'化,.)牛-％'化,,） (7) 式中，为离散时间点，^1、2、3、…、化励所方程误差序列离散的点数， 并求方程误差序列e侣，)快速傅里叶变换的表达式，并求FFT的表达式
(8) 524 ;求非线性目标函数的表达式
(9) 525 式（9)最小为优化目标，得到所述低频谐波部分中均值、各次谐波的幅值与相 位4 ? 4。巧，f;,巧g…-韦.巧的值； S26 ;将所述低频谐波部分中的谐波幅值為，4，--，4转化为非负值，相位 巧?巧，…?巧转化到区间W上，方法是；
(10)。
[0037]S3;更新迭代初值、高频谐波部分，包括： 531 ;计算所述方程误差序列e化。》，计算其所述误差均方根值C'。;，； 532 ;判断内层循环是否满足退出条件； ,j;所述误差均方根值减小到设定的下限， 义迭代次数i达到设定的上限， 如果满足其中一个条件，则转S4增加第次谐波的次数A、继续迭代， 如果都不满足，则转S33; 533 ;更新所述高频谐波部分，S35 ;设置非线性优化的初值X0,循环变量增加1 ;i=i+l， 转S22继续迭代，其中，用到如下公式与步骤： 5331 ;计算式（11)的值，
(11) 5332 ;则用FFT计算式（11)包含的谐波的幅值与相位，将其转化为余弦函数与正弦函 数之和的形式，得到幅值与6,,，求高频谐波部分的估计值的表达式
(12)。
[0038]S4;判断表所述低频谐波部分的最高次谐波次数的循环变量A是否不大于设定的 上限次数，或未达到精度条件， 如果是 ，则设置所述非线性优化的初值X0， 循环变量姑曽加1 ;仁知1，转步骤S2 ;否则，计算结束。
[0039] 步骤S4中，上限次数不大于步骤S2中所述方程误差序列离散的点数碱勺一半 m。
[0040]随着對I的增加，所述自由度系统谐波平衡法装置得到的结果精度越来越高。当 姑曽加一定程度时，所述误差均方根值保持不变或减小的很慢，该时，可W认为迭代终 止。
[0041] 步骤S3和S4分别为内层与外层循环的下一轮循环设置迭代初值。如图2所示， 步骤S3中X0有2知1个元素，对应于均值、颂H氏次谐波的幅值与相化为内层的下一轮循 环准备优化初值；步骤S4中X0有2知3个元素，对应于均值、知1项谐波的幅值与相位，为 外层的下一轮循环准备优化初值。外层循环增加1，意味着所述低频谐波部分加进去一项 谐波。步骤S4中的XO比步骤S3中的XO多出的2个元素为所述低频谐波部分新增的该一 项谐波对应的幅值与相位。对于步骤S3与S4中迭代初值X0的设定，有W下3种都切实可 行、效果良好的技术方案。
[0042] 技术方案1;步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1 到i的该J轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位。
[0043] 步骤S4中，所述非线性优化的初值X0包含2知3个元素，为外层循环从1到左的 该每爸循环中，第巧欠循环下，对应的内层循环中，第i次循环时，步骤S2与S3中得到的结 果，其中2知1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为 步骤S3中FFT计算得到的对应谐波的幅值与相位。
[0044] 技术方案2;步骤S3中，所述非线性优化的初值X0的设置方法为：内层循环从1 到i的该J轮循环中，第i次循环时，步骤S2中通过优化得到的均值、谐波幅值与相位。
[0045] 步骤S4中，所述非线性优化的初值X0包含2知3个元素，为外层循环从1到左的 该每％循环中，第域C循环下，对应的内层循环从1到i的该J轮循环中，最小的误差均方根 值旬W对应的结果，其中2知1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位； 另外2个元素为步骤S3中FFT计算得到的对应谐波的幅值与相位。
[0046]技术方案3 ;步骤S33的详细步骤为；判断所述方程误差均方根值是否小于上 一轮迭代得到的值，如果不小于，则非线性优化的初值X0在原来的基础上加上随机扰动， 循环变量增加1 ;i=i+l，转S2继续迭代；否则，执行S331与S332,迭代初值xO设置为步骤 S2中通过优化得到的均值、谐波幅值与相位，如果所述方程误差均方根值<^"。，与上一轮迭 代得到的误差均方根值之差小于设定的误差限，则在xO再加上一个随机扰动，可用随机函 数随机生成。
[0047] 步骤S4中，所述非线性优化的初值X0包含2知3个元素，为外层循环从1到左的 该每爸循环中，第巧欠循环下，对应的内层循环中，第i次循环时，步骤S2与S3中得到的结 果，其中2知1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为 步骤S3中FFT计算得到的对应谐波的幅值与相位。
[0048] 在步骤S2中，式（8)与巧）中都包含了 2知1为待求未知量，而该些式子本身的形 式随着非线性系统方程（1)、娜日铺勺变化而变化，并且差异巨大；因此，在步骤S21~S24采 用符号运算，可充分利用工程软件或数学软件的强大的符号运算功能、求导等数学运算功 能，极大地减小了编程的复杂度。步骤S3是为了估计所述高频谐波部分的各次谐波的幅值 与相位，采用数值运算可W极大地提高运算速度。
[0049] 本发明提供的单自由度系统谐波平衡法装置可W在很多高级编程语言、工程软 件、数学软件中实施，如MTLAB。具体实施的过程中，式（2)中的均值、谐波幅值与相位一共 2ifrl个待求解的量，需要定义对应的变量来处理。为了构造目标函数的方便，代表均值、幅 值与相位的符号变量由前缀与数字两部分组成。均值、幅值符号变量的前缀相同，相位符号 变量的前缀相同；第二部分数字为正整数数字，数字位数的确定方法是；设需要分析的响 应中最高谐波的次数为必则第二部分数字的位数的最小值为U〇ga.f)j+i，如果位数不 够，则高位由数字0补足。
[0050] 例如，均值、幅值符号变量的前缀为An,相位符号变量的前缀为化in。当必：100时， _!Og(M)j牛！ = 2，则在MATLAB中定义符号变量的方法如下；
[0051] 如果M=12,则上面方法定义了 13个代表均值与谐波幅值的符号变量AnOl、An02、 An03、…、41111、41112、41113，还定义了13个代表与上述幅值对应的相位的符号变量化111〇1、 化in02、Phin03、…、Phinll、化inl2、Phinl3,并且还定义了两个组号数组An与化in,都包 含了 13个元素，分别为幅值及相位的符号变量。因为AnOl与均值对应，符号变量化inOl 只是为了配合程序的方便，使对应次谐波的幅值与相位变量下标一致而定义的，在此并不 使用。
[0052] 基于上面定义的符号变量，所述目标函数构造单元采用的方法包括如下步骤： a;用char函数将表达式（4)转换为字符串myObjRmStrl; b;然后用字符串替换函数strrep将字符串myWj化nStrl的均值替换为字符串；x(l)，幅值替换为字符串；x(2)、x(3)、…、X化)、x化+1)，相位替换为字符串；x化+2)、 x化+3)、…、x(化）、x(化+l)，得到代表目标函数的字符串型表达式mパ)bj化nStr2; C;在字符串mW)bj化nStr2前面加上字符串' @(x)'，得到新的字符串mW)bj化nStr3 ;d;用eval函数将字符串mW)bj化nStr3转换为目标函数。
[0053] 步骤b主要采用strr巧函数，具体方法是：
[0054] 代表均值、幅值与相位第二部分数字的位数的最小值为|>g(MH+ 】|，高位补0是 为了保证在strrep函数替换的过程中不出现错误替换。如果不固定位数，则会出现错误。 例如；若将变量定义为Anl、An2、An3、…、Anil、Anl2、Anl3,贝ij函数strrep在A1替换为 义(1)时，会将41110、41111、41112、41113错误替换成义(1)0、义(1)1、义(1)2、义(1)3。
[0055] 在数学软件MATLAB中实施本发明。下面W 3个实施例来说明本发明的有益效果。
[0056] 实施例1;单自由度非线性系统的方程为 r+OJ.i'+v+0.75/ =COS/。
[0057] 假设响应中一次谐波的角频率为1，设置迭代初值x0=[0, 0, 0]，仅考虑32次W下 的谐波，即i^32,内层循环的最大次数为5,采样点数^64,步骤S4中误差限设为！X10-K' ，随机扰动的范围为实施技术方案3。随着對I的增加，所述方程误差 均方根值迅速减小，前6次循环时的方程误差均方值依次为： 5.656854249492381e+000,2. 4968591790980876-003,8. 156122475628342e-008， 4.360844843951736e-012,4. 360844843951736e-012,4. 360844843951736e-012。
[005引各次谐波的幅值与相位如图5所示，其中0次谐波对应于均值，1~6次谐波由步骤S2的非线性优化得到，7~32次谐波由步骤S3中的数值FFT得到。
[0059] 从上述的数据可W明显看出，当必4W后，迭代结果趋于恒定，精度级别为 lOe-12,即 1()。。
[0060] 实施例2;将实施例1中的激励修改为COS2 假设响应中一次谐波的角频率为2， 设置迭代初值xO为随机生成，其他参数同实施例1。实施技术方案3。各次谐波的幅值与 相位如图6所示，其中0次谐波对应均值，1~6次谐波由步骤S2的非线性优化得到，7~32次 谐波由步骤S3中的数值FFT得到。当必2W后，迭代结果趋于恒定，精度级别为10-^气
[0061] 实施例3;将实施例1中的激励修改为cos3f，假设响应中一次谐波的角频率为3， 设置迭代初值xO为随机生成，仅考虑32次W下的谐波，即i^32,内层循环的最大次数为4， 采样点数^64，。实施技术方案3。随着對I的增加，所述方程误差均方根值迅速减小，如 图7所示，迭代结果精度级别为！〇U。各次谐波的幅值与相位如图8所示，其中0次谐波 对应均值，1~6次谐波由步骤S2的非线性优化得到，7~32次谐波由步骤S3中的数值FFT得 到。
[0062] 本发明提供的装置只需要改变础勺值即可控制谐波的次数，应用方便。上述的3 个实施例的结果，说明了本发明提供的装置具有很快的收敛速度与很高的精度。
[0063] 如果要分析亚谐波，只需要根据要亚谐波的角频率来设置一次谐波的频率即可。
[0064] 也可W将目标函数改为方程误差的谐波均方根值，同样可达到很好的效果。
[0065] 需要注意的是，上述具体实施例仅仅是示例性的，在本发明的上述教导下，本领域 技术人员可W在上述实施例的基础上进行各种改进和变形，而该些改进或者变形均落在本 发明的保护范围内。本领域技术人员应该明白，上面的具体描述只是 为了解释本发明的目 的，并非用于限制本发明。本发明的保护范围由权利要求及其等同物限定。
【主权项】
1. 一种单自由度系统谐波平衡法装置，其特征在于，采用一种谐波平衡法，其在求解单 自由度非线性系统的响应时，仅考虑响应中的有限次谐波，将所述有限次谐波分为低频谐 波部分与高频谐波部分， 从所述低频部分只包括均值与一次谐波开始，通过优化方法求解均值、所述低频谐波 部分的谐波幅值与相位，估计所述高频部分的谐波幅值与相位，并且通过迭代提高精度， 不断增加所述低频部分包含的最高次谐波的次数，每次增加一项，同时去掉高频谐波 部分的最低次谐波项，以所述低频谐波部分计算得到的结果作为初始值，最终求解出响应 中均值、所述有限次谐波的幅值与相位， 所述单自由度谐波平衡法装置包括以下部件： 非线性系统配置部件，用于定义非线性系统的表达式与激励， 符号变量定义部件，用于定义非线性系统响应中均值、所述有限次谐波的幅值与相位 的符号变量， 计算部件，用于通过所述谐波平衡法计算响应中均值、所述有限次谐波的幅值与相位， 输出部件，用于存储所述计算部件得到的结果、以图像化的方式显示结果， 所述计算部件包括：计算模块准备模块、计算模块主体模块、计算模块后处理模块， 其中，所述计算部件准备模块，用于设置迭代的初始值、一个基波周期内时间的离散数 目， 所述计算部件主体模块，用于通过迭代的方法计算响应中均值、所述有限次谐波的幅 值与相位，以及评价参数，所述计算主体部件由两层循环组成，外层循环控制所述低频谐波 部分最高次谐波的次数，内层循环通过优化方法求解响应中所述低频部分中均值、各次谐 波的幅值与相位； 所述计算部件主体模块中设置了一个目标函数构造单元，用来构造优化求解所述低频 谐波部分中均值、各次谐波的幅值与相位时用的目标函数，在迭代过程中，目标函数的表达 式形式、待优化的所述低频部分中均值、各次谐波的幅值与相位的个数在动态变化，所述目 标函数构造单元就是用来构造这个目标函数。2. 根据权利要求1所述的自由度系统谐波平衡法装置，其特征在于，所述计算部件 中，求解响应中均值、所述有限次谐波的幅值与相位采用的方法包括以下步骤： 51 :代表所述低频谐波部分的最高次谐波次数的循环变量h内层循环变量i置初值： 知1，i=l，所述高频谐波部分置初值:::+ 0,设置优化求解的初值XO ; 52 :优化方法计算所述低频谐波部分中均值、各次谐波的幅值与相位，包括： 521 :计算所述低频谐波部分《/)的表达式式中，J0 Z 2为均值、為,与奶I分别为第《次谐波的幅值与相位，々为所述低频谐波部 分中最高次谐波的次数；￡+?为基波的角频率，即一次谐波的角频率； 522 :计算非线性系统的响应_7(()的表达式（6) ,HO=VlAO+}w(〇 (6) 523 :求方程误差序列e(U的符号表达式，求其FFT的符号表达式式中，Q为离散时间点，&1、2、3、···、#，Λ为所方程误差序列离散的点数； 524 :求目标函数的表达式S25:以式（9)最小为优化目标，得到所述低频谐波部分中均值、各次谐波的幅值与相 位為,-s，'4:，％…，為,％ 的值； S26 :将所述低频谐波部分中的各次谐波的幅值為,為.….為转化为非负值，相位 的,朽，·' ·,%转化到区间2?丨上，方法是：S3 :更新迭代初值、所述高频谐波部分，包括： 531:计算所述方程误差序列《^,"；)的数值表达式，计算方程误差均方根值€_; 532 :判断内层循环是否满足退出条件：η所述方程误差均方根值€_减小到设定的精 度，内层循环的迭代次数i达到设定的上限， 如果满足其中任一条件，则转S4增加所述低频谐波部分最高次谐波的次数A继续迭 代， 如果都不满足，则转S33 ; 533 :更新所述高频谐波部分，设置非线性优化的初值X0,循环变量增加I :i=i+l，转 S22继续迭代，其中，用到如下公式与步骤： 5331 :计算式（11)的具体数值，5332 :用FFT计算式（11)包含的谐波的幅值与相位，将其转化为余弦函数与正弦函数 之和的形式，得到幅值U w与求所述高频谐波部分的估计值的表达式（12)S4 :判断代表所述低频谐波部分的最高次谐波次数的循环变量A是否不大于设定的上 限次数，或未达到精度条件， 如果是，则设置所述非线性优化的初值X0,循环变量^增加1 :知奸1，转步骤S2 ; 否则，计算结束。3. 根据权利要求2所述的自由度系统谐波平衡法装置，其特征在于，步骤S4中，上限次 数不大于步骤S2中所述方程误差序列离散的点数施]一半#/2。4. 根据权利要求3所述的自由度系统谐波平衡法装置，其特征在于，步骤S4中，精度条 件为外层循环当前循环，即第妇％循环，与上一轮循环得到的所述方程误差均方根值相等 或之差小于设定的误差限。5. 根据权利要求4所述的自由度系统谐波平衡法装置，其特征在于， 步骤S3中，所述非线性优化的初值XO的设置方法为：内层循环从1到i的这2轮循环 中，第i7欠循环时，步骤S2中通过优化得到的均值、谐波幅值与相位， 步骤S4中，所述非线性优化的初值XO包含2奸3个元素，为外层循环从1到础勺这左 轮循环中，第々次循环下，对应的内层循环中，第i7欠循环时，步骤S2与S3中得到的结果， 其中2奸1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为步骤 S3中FFT计算得到的对应谐波的幅值与相位。6. 根据权利要求4所述的自由度系统谐波平衡法装置，其特征在于， 步骤S3中，所述非线性优化的初值XO的设置方法为：内层循环从1到i的这2轮循环 中，第i7欠循环时，步骤S2中通过优化得到的均值、谐波幅值与相位， 步骤S4中，所述非线性优化的初值XO包含2奸3个元素，为外层循环从1到础勺这左 轮循环中，第砂：循环下，对应的内层循环从1到i的这i轮循环中，最小的误差均方根值 对应的结果，其中2奸1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另 外2个元素为步骤S3中FFT计算得到的对应谐波的幅值与相位。7. 根据权利要求4所述的自由度系统谐波平衡法装置，其特征在于， 步骤S33的详细步骤为：判断所述方程误差均方根值是否小于上一轮迭代得到的 值，如果不小于，则非线性优化的初值XO在原来的基础上加上随机扰动，循环变量增加1 : i=i+l，转S2继续迭代；否则，执行S331与S332,迭代初值x0设置为步骤S2中通过优化得 到的均值、谐波幅值与相位，如果所述方程误差均方根值与上一轮迭代得到的误差均 方根值之差小于设定的误差限，则在x〇再加上一个随机扰动， 步骤S4中，所述非线性优化的初值XO包含2奸3个元素，为外层循环从1到础勺这左 轮循环中，第々次循环下，对应的内层循环中，第i7欠循环时，步骤S2与S3中得到的结果， 其中2奸1个元素为步骤S2中通过优化得到的均值、谐波幅值与相位；另外2个元素为步骤 S3中FFT计算得到的对应谐波的幅值与相位。8. 根据权利要求5、6或7所述自由度系统谐波平衡法装置，其特征在于，步骤S21~S24 采用符号运算，步骤S3采用数值运算。9. 根据权利要求8所述的自由度系统谐波平衡法装置，其特征在于，所述符号变量定 义部件中，当采用数学软件MATLAB实现时，代表均值、幅值与相位的符号变量由前缀与数 字两部分组成，第一部分中均值、幅值符号变量的前缀相同，相位符号变量的前缀相同，第 二部分数字为正整数数字，数字位数的确定方法是：设响应中所述有限次谐波的最高谐波 次数为私则第二部分数字的位数的最小值为j +1，如果位数不够，则前面用数 字〇补足。10. 根据权利要求9所述的自由度系统谐波平衡法装置，其特征在于，所述目标函数构 造单元采用的方法包括如下步骤： a :用char函数将表达式（4)转换为字符串myObjFunStrl ; b :然后用字符串替换函数strrep将字符串myObjFunStrl的均值替换为字符串： x(l)，幅值替换为字符串：x(2)、x(3)、…、xU)、xU+Ι)，相位替换为字符串：xU+2)、 x〇+3)、…、x(i^)、χ(·^_+1)，得到代表目标函数的字符串型表达式my0bjFunStr2 ; c :在字符串my0bjFunStr2前面加上字符串' @(χ)'，得到新的字符串my0bjFunStr3 ; d :用eval函数将字符串myObjFunStr3转换为目标函数。
【专利摘要】本发明公开了一种单自由度系统谐波平衡法装置，属于计算机工程应用领域。该装置包括非线性系统配置部件、符号变量定义部件、计算部件、输出部件。计算部件中设置了一个目标函数构造单元，用来动态构造优化求解均值、各次谐波的幅值与相位时需要的目标函数。该装置采用的谐波平衡法，仅考虑响应中的有限次谐波，将其分为低次谐波部分与高频谐波部分。从低频部分只包括均值与一次谐波开始优化、迭代提高求解精度，不断增加低频部分包含的最高次谐波的次数，同时减小高频谐波部分的最低次谐波项，以低次谐波部分计算得到的结果作为初始值，最终求解出响应中均值、有限次谐波的幅值与相位。本发明提供的装置具有很快的收敛速度与很高的精度。
【IPC分类】G06F17/14
【公开号】CN104881394
【申请号】CN201510295236
【发明人】陈国强, 王耿, 康件丽
【申请人】河南理工大学
【公开日】2015年9月2日
【申请日】2015年6月3日