分数阶微积分是研究任意阶次微分、积分算子的特性及应用的数学分析方法^[17].其中,Caputo定义的分数阶微积分,其初值条件是整数阶的,具有清晰的物理意义,常用于解决工程应用和物理问题.

Caputo积分的定义式^[17]为

^C_aD^-v_t f(t)=1/(Γ(ν))∫_a^t(t-τ)^v-1f(τ)dτ(1)

Caputo微分的定义式^[17]为

^C_aD^α_t f(t)=1/(Γ(n-α))∫_a^t(f⁽ⁿ⁾(τ))/((t-τ)^α-n+1)dτ(2)

其中, ^C_aD^-v_t和 ^C_aD^α_t分别为Caputo定义的积分算子和微分算子; t和a为算子运算的上下限; τ为积分变量; v和α为阶次, 0<v<1, n-1<α<n, n为正整数; Γ(·)为Gamma函数; f⁽ⁿ⁾(·)为函数f(·)的n阶导数.由定义可知,分数阶微积分与以前所有时刻点的信息相关,具有空间全域相关特性和记忆特性,能更准确地描述中间过程以及临界现象的各类系统.

Caputo微积分定义的Laplace变换为^[17]

L[^C₀D^-v_t f(t)]=s^-vF(s)(3)

L[^C₀D^α_t f(t)]=s^αF(s)-∑^n-1_k=0s^α-k-1f^(k)(0)(4)

其中, L[·]为拉普拉斯变换; s为复频率; F[s]是变量为s的函数; k≥0且为整数, n=α」+1.

式(3)和式(4)表明,Caputo定义的微积分的Laplace变换简单,只需得到函数f(t)的整数阶导数的初值(f^(k)(0), k=0, 1, …,n-1)即可,极大地降低了分数阶微积分方程的求解难度,更具工程实用性.

分数阶微积分运算的求解,工程上常用波特图频域近似算法^[18-19],先对分数阶微积分进行时域-频域转换,再在频域中应用分段线近似法进行计算.由式(3)可知,阶次为α(0<α<1)的分数阶积分的传递函数为H(s)=1/s^α, 用单极点分数幂表示为H(s)=1/(1+s/p_T)^α, 其波特图为一条斜率为-20α dB/dec的直线.如图1,在工程计算中,通常用斜率为0和-20 dB的线段组成的锯齿线来近似,把分数阶函数转化为求解系统的零极点对的问题^[18].

图1 1/(1+s/pT)α的波特图及其锯齿线逼近[18]
Fig.1 Bode plot of 1/(1+s/pT)α and its approximation with zigzag straight lines[18]

假设最大角频率为ω_max的频率范围内,工程实际所要求动力学系统变量的计算误差不超过y(y为正值,单位:dB),工作频段限制为N个,取0<α<1, 则1/s^α的近似传递函数^[19]为

H(s)=1/s^α=1/((1+s/p_T)^α)≈

(∏^N-1_i=0(1+s/z_i))/(∏^N_i=0(1+s/p_i))(5)

其中, s为复频率; p_T为波特图中-3α dB处对应的角频率; z_i=(kd)ⁱkp₀, i=0, 1, …, N-1, 是系统所有的零点值, N=int[(lg(ω_max/p₀))/(lg(kd))]+1; p_i=(kd)ⁱp₀, 为系统所有的极点值, i=1, 2, …, N; p₀=p_T10^y/(20α), k=10^y/[10(1-α)], d=10^y/(10α). 该算法的实质是把分数阶积分用N+1阶的线性传递函数来近似.

分数阶模拟电容可采用树型、链型和混合型的阻容分抗电路进行模拟^[6].在复频域中,电容量为C₀、 阶次为α的分数阶电容的传递函数为 F(s)=1/(C₀s^α). 对链型阻容分抗电路,其两端之间的复频域表达式为^[6]

H(s)=(1/C₁)/(s+1/(R₁C₁))+(1/C₂)/(s+1/(R₂C₂))+

(1/C₃)/(s+1/(R₃C₃))+…+(1/C_n)/(s+1/(R_nC_n))(6)

根据工程需求,选定α、 ω_max、 p_T和y的值并代入式(5)计算,得到分数阶传递函数H(s)=1/s^α(α=0.1～0.9)的频域近似表达式.设分数阶模拟电容的值为C₀、 阶次为α, 对比式(5)和式(6)的参数,计算得R₁, R₂, …, R_n和C₁, C₂, …, C_n的值.

图2是由运算放大器组成的一种广义阻抗变换器(general impedance converter, GIC)电路^[20],它既能模拟电容,也能模拟电感,电路所呈现的阻抗性质由Z₁～Z₄及Z_L所选择的电容或电阻来决定.

图2 一种GIC电路[20]
Fig.2 A general impedance converter circuit[20]

若运算放大器是理想的,则该电路输入阻抗为

Z_i(s)=(Z₁(s)Z₃(s))/(Z₂(s)Z₄(s))Z_L(s)(7)

在图2中,若Z₁、 Z₂、 Z₃和Z_L为电阻R₁、 R₂、 R₃和R_L, 且Z₄为阶次为α的电容C时,其输入阻抗Z_i(s)=s^α(R₁R₃C)/(R₂)R_L=s^αL_eq, L_eq为电感量,可模拟阶次为α电感量为L_eq的电感; 若Z₁和Z₃分别是阶次为α的电容C₁和阶次为β的电容C₃时, 则当Z₂、 Z₄和Z_L为电阻R₂、 R₄和R_L时,其输入阻抗Z_i(s)=1/(s^α+βC₁C₃R₂R₄)R_L=1/(s^α+βC_eq), C_eq为电容量,可模拟阶次为α+β、 电容量为C_eq的电容.因此,应用GIC电路,把α阶分数阶电容转换为α阶分数阶电感,把分数阶电容的阶次扩展为0～2阶,最终可实现0～2阶的分数阶电容和分数阶电感.

根据式(7)和图2的GIC电路,可把β阶的分数阶电容回转为β阶的模拟接地电感.图3是基于运放LF351D设计的GIC电路,取R₁=R₂=R₃=R_L=1 kΩ, C₄为β阶电容量为1 μF的分数阶电容,可得到β阶电感量为1 H的模拟接地电感,取不同的电阻值可得到不同的模拟电感值.分数阶电容C₄采用文献[8]的参数ω_max=100 rad/s、 p_T=0.01 rad/s、 y=2 dB,由式(5)和式(6)计算出β阶次、电容值为1 μF时链型分抗电路的阻容元件参数.

图3 β阶电容回转为β阶电感的GIC电路
Fig.3 A GIC circuit to transform β-order capacitor into β-order inductor

用正弦电压u_i=sin(wt)激励时,流过β阶电感的电流相位比电压超前βπ/2. 用Multisim软件对图3的电路进行仿真,得到β阶模拟接地电感的正弦电压激励的电流波形,在虚拟示波器上显示的波形如图4.波形①为电感两端的电压波形,波形②、③和④分别对应阶次β为0.4、0.8和1.0,电感值为1 H时的电感电流,电压和电流的转换关系为每500 mV表示1 mA.可见,电感电流与电压的相位关系与理论分析基本一致.

图4 分数阶电感的正弦电压激励响应
Fig.4 The sine wave response of the fractional order inductor

为验证所设计分数阶模拟电容和模拟电感的正确性,对分数阶L_βC_α串联电路进行仿真.使用Multisim分析自带的波特仪对串联电路的幅频特性,得到分数阶串联电路的谐振频率,再与理论计算得到的谐振频率ω_res对比. ω_res的计算公式为^[12]

ω_res={((α-β)cos[(α+β)π/2]+((α-β)²cos²[(α+β)π/2]+4αβ)^1/2)/(2βL_βC_α)}^1/(α+β)(8)

选取分数阶模拟电容和模拟电感的谐振参数时,应确保谐振频率位于其正常工作频域内.由式(8)可得,当α=0.4, β=0.8, C_α=1 μF,L_β=5 kH时, f_res=ω_res/(2π)=10.80 Hz,与图5(a)仿真得到的谐振频率基本一致; 当α=0.8, β=0.4, C_α=1 μF, L_β=6 kH时, f_res=ω_res/(2π)=13.78 Hz,与图5(b)中仿真得到的谐振频率基本一致.可见,电路仿真得到的谐振频率与理论计算得到的基本一致,证明所设计的分数阶模拟电容和模拟电感有效可行.

图5 分数阶LβCα串联电路仿真的波特图
Fig.5 Simulation Bode plot of fractional order LβCα serial circuit