一种基于多项式拟合的变阶分数阶微积分频域分析方法与流程

本发明涉及电子信息技术领域，尤其涉及一种基于多项式拟合的变阶分数阶微积分频域分析方法。

背景技术：

分数阶微积分学的发展起源于19世纪60年代，但其理论在近年来才被应用到工程技术中。目前主要应用领域包括了：涵盖粘弹性分析、反馈放大器、电容理论、控制领域、生物系统电导率分析、分数阶神经系统建模等领域。尤其是在一些自动化生产线中，会采用分数阶微积分优化自动控制流程，从而提高生产效率以及控制精度。

在现阶段的方案中，通常通过oustaloup滤波器来实现分数阶微积分，可以获得复频域s离散传递函数，但是在这种计算传递函数的方式，实时性较低，在对于实时性要求较高的应用领域中(比如自动控制领域)，限制了该方案的使用范围。

技术实现要素：

本发明的实施例提供一种基于多项式拟合的变阶分数阶微积分频域分析方法，能够提高传递函数的计算效率。

为达到上述目的，本发明的实施例采用如下技术方案：

所述方法用于一种分析系统，所述分析系统用于为外部系统提供分数阶微积分计算，所述分析系统包括：计算机设备、利用pci接口连接计算机的quanser板卡(qpide)、二级倒立摆；其中辅助控制二级倒立摆的硬件设备还有多功能数据采集卡、模拟传感器适配器电源和信号路由器。所述计算机设备上运行的分析系统是依靠quarc通用快速控制实时仿真系统，其通过matalb中simulink设计的控制机，直接生成实时代码，并将工业级的实时应用程序下载到电脑操作系统中。；所述计算机设备上运行仿真建模程序，所述仿真建模程序用于运行模拟了所述外部系统的架构和运行状态的仿真模型；

所述方法包括：

步骤(1)：对于分数阶积分算子所述计算机设备令阶次λ在(01)范围内取m个均匀点，并记为λi，其中i＝1,2,…,m；并且，对于分数阶微分算子所述计算机设备令阶次λ在(-10)范围内取m个均匀点，记为λi，其中i＝1,2,…,m；

骤(2)：所述计算机设备对于各阶次的λi，得到分数阶积分的频域整数阶近似传递函数；

步骤(3)：根据所述频域整数阶近似传递函数，构造分母系数矩阵a；

步骤(4)：根据所述频域整数阶近似传递函数，构造分子系数矩阵b；

步骤(5)：以离散点阶次x＝[λ0,λ1,…λm]为输入，以分母系数矩阵a的列向量转置为输出，采用最小二乘方法进行拟合计算，得到多项式fan(λ)，其中，n表示大于等于0的整数，a[n,1],a[n,2],…a[n,m]代表矩阵a第一列到最后一列的数值；

步骤(6)：以离散点阶次x＝[λ0,λ1,…λm]为输入，以分子系数矩阵b的列向量转置为输出，采用最小二乘方法进行拟合计算，得到多项式fbn(λ)；

步骤(7)：由多项式fan(λ)和fbn(λ)得到分数阶微积分的分子系数fb2n+1(λ),fb2n(λ),…,fb0(λ)和分母系数fa2n+1(λ),fa2n(λ),…,fa0(λ)，并通过所述分子系数和所述分母系数，构建分数阶积分的频域整数阶传递函数。

本发明实施例提供的基于多项式拟合的变阶分数阶微积分频域分析方法，通过频域拟合获得多个等间隔阶次时分数阶微积分算子的频域连续传递函数；然后构造频域传递函数同位置系数矩阵，并以系数矩阵的列向量进行多项式拟合，获得阶次为输入、系数为输出的拟合函数组；最后对于任意阶次的分数阶微积分算子，可快速解算并构建出其连续传递函数，从而实现变阶次分数阶微积分的频域数值计算。相对于现有的ousta_fod方法，本实施例提出的方法能够更快地运算出其传递函数。从而可以克服ousta_fod方法中无法实时获得传递函数的缺点，提高了传递函数的计算效率。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本发明实施例提供的系统架构示意图；

图2为本发明实施例提供的逻辑原理示意图；

图3至图14为本发明实施例提供的y相对于x的散点图。

具体实施方式

为使本领域技术人员更好地理解本发明的技术方案，下面结合附图和具体实施方式对本发明作进一步详细描述。下文中将详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的任一单元和全部组合。本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

本发明实施例提供一种变阶分数阶微积分频域分析方法，所述方法用于一种如图1所示的分析系统所述分析系统用于为外部系统提供分数阶微积分计算，所述分析系统包括：计算机设备、利用pci接口连接所述计算机设备的quanser板卡(qpide，一种多功能数据采集卡)、所述quanser板卡与模拟传感器适配器电源和信号路由器组成的通信模块连接二级倒立摆；其中辅助控制二级倒立摆的硬件设备还有多功能数据采集卡、模拟传感器适配器电源和信号路由器。所述计算机设备上运行的分析系统是依靠quarc通用快速控制实时仿真系统，其通过matalb中simulink设计的控制机，直接生成实时代码，并将工业级的实时应用程序下载到电脑操作系统中。其中，具体可以在计算机设备利用matlab仿真建模，将外部系统全部仿真到计算机设备中，建立控制系统的仿真模型。在matlab平台上编写步骤1到步骤7的程序(一种用于外部系统的控制算法)，在给定一个输入变量后，通过quanser板卡控制二级倒立摆。quanser板卡具体可以采用quanser公司的qpide卡。

在本实施例中，二级倒立摆具体可以采用目前的倒立摆控制系统。倒立摆控制系统是一个一个复杂的、不稳定的、非线性系统的实验平台，主要用于控制理论教学及开展各种控制实验。对倒立摆系统的研究能有效的反映控制中的许多典型问题:如非线性问题、鲁棒性问题、镇定问题、随动问题以及跟踪问题等。通过对倒立摆的控制，用来检验新的控制方法是否有较强的处理非线性和不稳定性问题的能力。同时，其控制方法在军工、航天、机器人和一般工业过程领域中都有着广泛的用途，如机器人行走过程中的平衡控制、火箭发射中的垂直度控制和卫星飞行中的姿态控制等。

所述方法的执行流程包括：

步骤(1)：对于分数阶积分算子所述计算机设备令阶次λ在(01)范围内取m个均匀点，并记为λi，其中i＝1,2,…,m；并且，对于分数阶微分算子所述计算机设备令阶次λ在(-10)范围内取m个均匀点，记为λi，其中i＝1,2,…,m，分数阶积分算子和分数阶微分算子的符号都为其中积分的运用范围是(0,1)，微分的运用范围是(-1,0)。

具体的，s表示拉普拉斯变换后的算子，例如设一个系统的输入函数为x(t)，输出函数为y(t)，则y(t)的拉氏变换y(s)与x(t)的拉氏变换x(s)的商:w(s)＝y(s)/x(s)称为这个系统的传递函数。，后面的斜体和正体是一样的。；

例如：通过quanser板卡的多功能采集器，从外部系统采集输入变量给到了外部二级倒立摆。

步骤(2)：所述计算机设备对于各阶次的λi，得到分数阶积分的频域整数阶近似传递函数；

其中，具体可以运用oustaloup滤波器方法得到分数阶积分的频域整数阶近似传递函数。

步骤(3)：根据所述频域整数阶近似传递函数，构造分母系数矩阵a；

步骤(4)：根据所述频域整数阶近似传递函数，构造分子系数矩阵b；

步骤(5)：以离散点阶次x＝[λ0,λ1,…λm]为输入，以分母系数矩阵a的列向量转置为输出，采用最小二乘方法进行拟合计算，得到多项式fan(λ)，其中，n表示大于等于0的整数，a[n,1],a[n,2],…a[n,m]代表矩阵a第一列到最后一列的数值；

步骤(6)：以离散点阶次x＝[λ0,λ1,…λm]为输入，以分子系数矩阵b的列向量转置为输出，采用最小二乘方法进行拟合计算，得到多项式fbn(λ)；

步骤(7)：由多项式fan(λ)和fbn(λ)得到分数阶微积分的分子系数fb2n+1(λ),fb2n(λ),…,fb0(λ)和分母系数fa2n+1(λ),fa2n(λ),…,fa0(λ)，并通过所述分子系数和所述分母系数，构建分数阶积分的频域整数阶传递函数。

由计算机设备中的程序，计算得到传递函数的完整表达式后，可以输入quanser板卡，计算机设备的显示模块(比如屏幕)显示quanser板卡在输入传递函数的运行状态，其中，步骤(1)到步骤(5)可以采用matlab平台编写，步骤(6)可以采用目前主流的程序语言编写。

在本实施例中，所述步骤(2)，具体包括：

所述计算机设备对于各阶次的λi，得到分数阶积分的频域整数阶近似传递函数

其中，对于任一阶次λi，可以运用目前oustaloup滤波器方法得到分数阶积分的频域整数阶近似传递函数：其中2n+1为滤波器阶次，b[0,i],…,b[2n+1,i]和a[0,i],…,a[2n+1,i]为系数,i＝1,…,m。

在本实施例中，所述步骤(3)，具体包括：根据m个离散的频域整数阶近似传递函数，获取分母系数矩阵a：其中a2n+1,…,a0为2n+2个m×1维列向量。

在本实施例中，所述步骤(4)，具体包括：根据m个离散的频域整数阶近似传递函数，获取分子系数矩阵b：其中b2n+1,…,b0为2n+2个m×1维列向量。

在本实施例中，所述步骤(5)，具体包括：以离散点阶次x＝[λ0,λ1,…λm]为输入，以分母系数矩阵a的列向量转置为输出，采用最小二乘方法进行拟合计算，得到多项式fan(λ)＝n[h,n]λ^h+n[h-1,n]λ^h-1+...+n[1,n]λ+n[0,n]，n＝0,1,…,2n+1，和多项式fbn(λ)＝m[h,n]λ⁵+m[h-1,n]λ⁴+...+m[1,n]λ+m[0,n]，n＝0,1,…,2n+1；

其中，λ为任意阶次输入，fbn(λ)表示离散传递函数分子多项式中系数bn，fan表示离散传递函数分母多项式中系数an。这里以系数矩阵a和b为基础构建多项式fan(λ)和fbn(λ)，具体的包括：以离散点阶次x＝[λ0,λ1,…λm]为输入，以系数矩阵a的列向量转置为输出，采用最小二乘方法进行拟合，可得到如下2n+2多项式：fan(λ)＝n[h,n]α^h+n[h-1,n]α^h-1+...+n[1,n]α+n[0,n]，n＝0,1,…,2n+1，其中h为多项式最高项阶次，可以根据精度需要选取适当数值。

以离散点阶次x＝[λ0,λ1,…λm]为输入，以系数矩阵b的列向量转置为输出，采用最小二乘方法进行拟合，可得到如下2n+2多项式：fbn(λ)＝m[h,n]λ^h+m[h-1,n]λ^h-1+...+m[1,n]λ+m[0,n]，n＝0,1,…,2n+1，其中h为多项式最高项阶次，可以根据精度需要选取适当数值。

在现有方案中，通常通过oustaloup滤波器直接来实现分数阶微积分，可以获得复频域s离散传递函数:其中，α为分数阶微积分的阶次，2n+1表示滤波器的阶次,a2n+1…a0为分母系数,b2n+1…b0为分子系数。可以看出现有技术所设计的积分仅仅能在知道α和n的情况下，得出传递函数，例如取α＝0.5，n＝2得传递函数：因此此方法无法做到实时获得传递函数。

具体的计算逻辑如图2所示，包括：本发明公开了一种快速得到任意阶次分数阶微积分的频域传递函数实现方法。该方法首先在(01)或(-10)区间内取m个均匀阶次点，然后运用oustaloup滤波器方法得到m个离散阶次分数阶微积分的频域传递函数组；然后利用上述传递函数组的分子系数、分母系数利用最小二乘方法构造拟合多项式；最后对于任意可变阶次，可迅速通过多项式函数获得该阶次的对应分数阶微积分的频域传递函数。在matlab平台上编写步骤(1)到步骤(5)的程序，在给定一个输入变量后，向quanser板卡传输实时获取的传递函数，并通过quanser板卡再输出给外设，从而控制二级倒立板。

本发明所要解决的技术问题是提供一种快速得到任意阶次的离散s传递函数的方法。在oustaloup递推滤波器方法的基础上，基于此提出并设计了一种变阶次微积分α(t)的数值实现方法，通过对比验证了基于oustaloup滤波器分数阶离散传递函数系数是可拟合的。并依据此对变阶次分数阶微积分数值计算进行频域二次拟合，方便了变阶次分数阶的数值方法实现。

在本实施例中，提出一种变阶次分数阶微积分数值计算的频域二次拟合法。该方法首先用oustaloup递推滤波器方法，通过频域拟合获得多个等间隔阶次时分数阶微积分算子的频域连续传递函数；然后构造频域传递函数同位置系数矩阵，并以系数矩阵的列向量进行多项式拟合，获得阶次为输入、系数为输出的拟合函数组；最后对于任意阶次的分数阶微积分算子，可快速解算并构建出其连续传递函数，从而实现变阶次分数阶微积分的频域数值计算。通过多项式计算的方法快速获得变阶分数阶微积分的频域传递函数,且具有很高的精度。

对于任意可变阶次λ，根据多项式fan(λ)和fbn(λ)可以快速得到分数阶微积分的分子系数fb2n+1(λ),fb2n(λ),…,fb0(λ)和分母系数fa2n+1(λ),fa2n(λ),…,fa0(λ),因此可以快速构建的频域整数阶传递函数：

参照图3至图14所示的y相对于x的散点图。其中，离散的阶次x＝[α1,α2,…,αm]为横坐标,以为纵坐标。从图3-图14中可以看出，阶次在(0,1)上时得到的s传递函数系数可以拟合成连续的曲线。

相对于现有方案中，通常通过oustaloup滤波器直接来实现分数阶微积分的ousta_fod方法，本实施例提出的方法能够更快地运算出其传递函数，从表1-3中任意取的三个阶次的运算时间如表4(计算环境：cpu主频2ghz)：

表1α＝0.23时系数分析表

表2α＝0.57时系数分析表

表3α＝0.72时系数分析表

表4运算时间分析表

由此可见，相对于ousta_fod方法，本实施例提出的方法能够更快地运算出其传递函数；且能够更快的实现变阶次分数阶微积分的数值计算；还有就是相对误差也较小。从而可以克服ousta_fod方法中无法实时获得传递函数的缺点，可快速解算并构建出其连续传递函数，从而实现变阶次分数阶微积分的频域数值计算。因此，本发明提供的方案实现了快速求解分数阶传递函数，提高了传递函数的计算效率。也间接提高了在工业控制时的精度以及效率。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于设备实施例而言，由于其基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。