第5章
分数阶微积分算子与系统的近似
动态系统是描述很多物理现象的数学模型基础。从系统分析与描述角度看，通常可以将系统分为线性系统和非线性系统。从本章开始将引入系统的概念。5.1节首先介绍基于 MATLAB控制系统工具箱的整数阶线性系统建模与分析工具，为以后将介绍的系统行为的比较奠定基础。
在前面章节介绍分数阶微分 Dαy(t)数值计算的时候，一直假设函数 y(t)或其采样值是已知的，这样先构造出 y的采样值向量，然后才能计算出其分数阶的导数与积分。如果分数阶算子 Dα被一个事先未知的信号驱动，就像在控制系统中分数阶受控对象被控制器驱动那样，则前面讨论的方法是不能计算这样的分数阶动作的。这时，通常需要建立起一个能模拟分数阶行为的装置完成这样的任务。有了这样的装置，就可以将信号馈入该装置，而该装置的输出就是所需的 α阶微分信号。当然，这样的思想同样适用于分数阶积分信号的获得。
若想要设计这样一个装置，比较有效且常用的思路是设计一个整数阶的线性连续滤波器。如线性连续的传递函数模型 G(s)，使得其频域响应尽可能好地逼近原始的分数阶算子模型。
由理论分析可知，分数阶微分环节 sα的 Bode幅频特性就是斜率为 20α dB/dec的斜线，而相频特性就是值为 απ/2的水平直线，而整数阶连续传递函数模型的幅频特性的渐近线为斜率为 20k dB/dec的斜线，且 k只能为整数。所以不可能存在一个整数阶传递函数模型，使得在整个频率范围内都能对分数阶行为进行完全一致的逼近，所以只能退而求其次，得出在某个频率段内能逼近分数阶行为的滤波器。
在早期整数阶传递函数逼近的研究中，人们主要采用连分式类近似方法，包括低频与高频部分的连分式近似、Carlson近似与 Matsuda–Fujii近似等，5.2节将给出这类近似与滤波器设计的一般介绍。以现在的观点来看，连分式类近似方法效果很差，不太适合实用，因为没有办法指定感兴趣的频率段，而频率段的选择在分数阶系统的数值仿真中是至关重要的。
·114·分数阶微积分学——数值算法与实现
法国学者 Oustaloup教授及其同事提出的分数阶算子的滤波器近似开启了复杂分数阶系统仿真的新时代 [1]。这类滤波器允许用户自行选择感兴趣的频率段与阶次，利用整数阶传递函数模型逼近分数阶微积分算子。这类滤波器的设计将在
5.3节给出详细介绍，该节还介绍这类滤波器的改进形式。如果分数阶系统或传递函数中每个分数阶算子都被这类滤波器替代，则可以将整个分数阶系统用高阶整数阶传递函数逼近，也可以用低阶最优模型降阶技术对其逼近。 5.4节将给出这些近似的具体实现方法。
对于那些不能由分数阶传递函数标准形式描述的无理系统，例如，如果系统中含有 pγ(s)项，其中 p(s)为整数阶传递函数模型。 5.5节将介绍三类近似方法：第一类是 1/(τs + 1)ν隐式模型的时域近似；第二类是用频率响应拟合的方法设计滤波器；第三类是 Charef滤波器设计技术。该节还将介绍最优 Charef滤波器设计算法及其实现。
除了前面介绍的连续滤波器近似之外， 5.6节还介绍用离散滤波器或离散化的滤波器对分数阶算子进行近似[2，3]。不过，从频域响应拟合角度看，拟合模型在高频时常伴有极高的增益，从而使这类滤波器对噪声信号放大效果严重，不一定适合实际应用。
本书将各种通过频域响应拟合得出的高阶整数阶传递函数统称为滤波器设计，因为其作用是让某个函数通过拟合模型后得出所期望的分数阶运算。
5.1线性整数阶模型的表示与分析
MATLAB的控制系统工具箱对整数阶线性时不变系统的模型输入与分析提供了强大的支持。作为本章的底层支持，这里简要介绍在 MATLAB环境下线性整数阶系统的建模与分析工具。
5.1.1数学模型输入与处理
整数阶线性时不变（linear time invariant，LTI）连续系统一般由传递函数和状态方程表示。本节先介绍这两种模型的数学形式，再介绍利用 MATLAB输入这些模型的方法。
系统的传递函数模型有两种输入方法：
（1）先用 num=[b1,b2,··· ,bm+1]、den=[a1,a2,··· ,an+1]命令输入传递函数的分子、分母多项式，然后用 G=tf(num,den)命令构造传递函数模型。
（2）由 s=tf('s')命令定义 Laplace算子 s，然后在 MATLAB下用表达式输入传递函数模型。若已知状态方程模型，即已知 A、B、C和 D矩阵，则由 G=ss(A,B,C,D)命令可以输入系统的状态方程模型。
如果已知 G1(s)和 G2(s)两个模块串联连接，且这两个模块分别用 G1、G2两个变量表示，则串联系统总的模型可以由 G=G2 * G1直接求出；如果这两个模块并联连接，则可以由 G=G1 + G2直接求出总模型；如果这两个模块负反馈连接，其中前向通路为 G1，反向通路为 G2，则总系统可以由 G=feedback(G1,G2)命令直接求出，而求正反馈总模型应该用 G=feedback(G1,G2,1)命令。
5.
1.2时域与频域响应
若已知系统的模型变量 G，则不论是传递函数还是状态方程模型，均可以对该系统进行分析。例如，由 bode(G)、nyquist(G)、nichols(G)函数可以分别绘制系统的 Bode图、Nyquist图和 Nichols图；由 step(G)、impulse(G)语句可以绘制系统的阶跃或冲激响应；若已知时间向量 t和激励信号的采样点训练 u，则可以调用 lsim(G,u,t)直接绘制时域响应曲线，或由 y=lsim(G,u,t)获得时域响应数据向量 y。
5.
1.3分数阶线性系统的建模与分析
作者的 FOTF工具箱是在仿照上面的建模与分析语句调用格式的基础上开发的，系统分析函数与调用格式尽量与上面介绍的函数保持一致，可以直接调用。后面还将专门介绍 FOTF工具箱的使用方法。
·116·分数阶微积分学——数值算法与实现
5.2基于连分式的几种近似方法
在早期的研究中，人们经常采用连分式展开的形式设计滤波器，而连分式技术也经常被认为是逼近某些非线性函数的有效工具。
5.2.1连分式近似
从理论上看是不可能直接对分数阶算子 s.α进行连分式近似的，所以，研究者分别在高频和低频段对 G1(s) = 1/(1+Ts)α和 G2(s) = (1+1/s)α进行连分式近似。其中，当 ωT .1时用 G1(s)逼近 1/sα；而 ω .1时，用 G2(s)逼近该积分算子。
MATLAB的符号运算引擎 MuPAD提供了底层函数 contfrac()，作者编写了该函数的一个接口函数，对给定函数进行连分式近似。
function [F,r]=contfrac(f,s,n,a)
arguments, f(1,1), s(1,1), n(1,1), a(1,1), end
F=feval(symengine,'contfrac',f,[inputname(2) '=' num2str(a)],n);
if nargout==2
r=feval(symengine,'contfrac::rational',F); end, end该函数的调用格式为 [f1,r]=contfrac(f,s,n,a)，其中， f是被近似函数的符号表达式，s为自变量，n为连分式的级次，a为参考点。返回的 f1为连分式的符号表达式，r为有理式。遗憾的是，由于近期版本 MuPAD底层的 contfrac()函数似乎不再支持求取函数的连分式表达式，这里给出的接口函数不能正常工作。如果想得出连分式可以尝试早期的版本，如 MATLAB 2018b，或改用 Padé近似技术。
例 5-1试分别在高频与低频段用连分式近似分数阶积分算子 1/ √ s，并评价其逼近效果。 解选择参考点 a =2，连分式级次为 n =9，可以得出 4阶高阶近似模型。
>> syms s; T1=1/(1+s)^0.5; [c1,G2]=contfrac(T1,s,9,2)
得出的连分式模型为
√
3 s . 2
c1(s) ≈ +
√
s . 2
3 .6 3+ √
23 s . 2
. +
√
s . 2
9
.54 3+ √
23 s . 2
. + √
s . 2
45 .150 3+ √
23
. + ···
105
对应的有理函数近似模型为
s4 + 112s3 + 1464s2 + 4864s + 4240
G2(s)= √√ √ √√
3s4 + 288 3s3 + 1944 3s2 + 4032 3s + 2448 3
可以看出，这里得出9的近似模型与下面显示的文献 [4]中的模型不一致：
0.3513s4 +1.405s3 +0.8433s2 +0.1574s +0.008995
Gh(s)=
s4 +1.333s3 +0.478s2 +0.064s +0.002844
s4 +4s3 +2.4s2 +0.448s +0.0256
Gl(s)=
9s4 + 12s3 +4.32s2 +0.576s +0.0256
一个可能的原因是选择的参考点不同。近似模型和原始模型 1/ √ s的 Bode图比较可以由下面语句得出，如图 5-1所示，其中，图中标记的英文字符是自动生成的， Bode Diagram为 Bode图，Magnitude（dB）为幅值，单位为分贝（dB），Phase（deg）为相位，单位为度，Frequency（rad/s）为频率，单位为 rad/s。
图 5-1连分式拟合的 Bode图比较
>> w=logspace(-3,3); s=fotf('s'); G0=s^-0.5; s=tf('s');
G2=(s^4+112*s^3+1464*s^2+4864*s+4240)/(9*3^(1/2)*s^4+...
288*3^(1/2)*s^3+1944*3^(1/2)*s^2+4032*3^(1/2)*s+2448*3^(1/2));
·118·分数阶微积分学——数值算法与实现
n=[0.3513 1.405 0.8433 0.1574 0.008995];
d=[1 1.333 0.478 0.064 0.002844]; Gh=tf(n,d);
n=[1 4 2.4 0.448 0.0256]; d=[9 12 4.32 0.576 0.0256];
Gl=tf(n,d); H=bode(G0,w); bode(H,G2,Gh,Gl)
bode(H,'-',G2,'--',Gh,':',Gl,'-.')
可以看出，模型 G2(s)的拟合效果比 Gh(s)的稍有改善，但是在低频段较差，从整个频率段的拟合看并不是很好。 例 5-2试比较几个滤波器对 f(t)= e.t sin(3t + 1)信号的 0.5阶积分的近似效果，并与 Grünwald–Letnikov分数阶积分的计算结果进行比较。 解可以用下面的语句得出这几个滤波器的滤波结果，这些结果是由 MATLAB控制系统工具箱 lsim()函数计算出来的，滤波结果与理论值曲线如图 5-2所示。可以看出，这样的滤波结果与理论值的差距太大，所以这种滤波器不能真正用于分数阶系统的仿真。
图 5-2 0.5阶积分的比较
>> t=0:0.01:10; y=exp(-t).*sin(3*t+1);
y0=glfdiff9(y,t,-0.5,5); y1=lsim(G2,y,t);
y2=lsim(Gh,y,t); plot(t,y0,'-',t,y1,'--',t,y2,':')
5.2.2 Carlson近似
另一个基于连分式的近似方法是 Carlson方法 [4]。Carlson方法的目标是对给定整数阶基础模型 G(s)分数次幂的形式找到整数阶近似模型 H(s)，即
H(s) ≈ Gα(s) （5-2-2）
基于上述的算法，可以编写出下面的 MATLAB函数设计 Carlson滤波器:
function H=carlson_fod(alpha,G,iter,epsx)
arguments
alpha(1,1) double, G(1,1)
iter(1,1) {mustBePositiveInteger}=2, epsx(1,1)=1e-4;
end
q=alpha; m=q/2; H=1;
for i=1:iter
H=minreal(H*((q-m)*H^2+(q+m)*G)/((q+m)*H^2+(q-m)*G),epsx); end, end该函数的调用格式为 H=carlson_fod(α,G,n,.)，其中 G为整数阶基础模型，α是分数阶幂次， n是递推的次数， H为设计出的 Carlson滤波器。请注意，滤波器的实际阶次随着递推次数增加而迅速增加，一般情况下选择 n =2即可。用户还可以选择较大的误差容限 .，例如，取值 10.4，对得出的近似模型做最小实现处理。
例 5-3试用 Carlson算法给 1/ √ s设计一个滤波器，并观察频域响应的拟合效果，另外，求给定函数 f(t)= e.t sin(3t + 1)的 0.5阶积分，并评价其精度。 解对这个具体的问题而言，基础模型为 G(s) = 1/s，幂次为 α =0.5，选择递推次数为 2，这样可以设计出 Carlson滤波器：
>> s=tf('s'); G=1/s; alpha=0.5;
H2=carlson_fod(alpha,G,2,1e-4) %选择一个较大的误差容限
这里得出的模型如下，与文献 [4]给出的是等效的：
0.1111s4 +4s3 + 14s2 +9.333s +1
H2(s)=
s4 +9.333s3 + 14s2 +4s +0.1111
可以由下面的语句绘制出频域响应拟合曲线，如图 5-3所示。可以看出，如果选择感兴趣的频率区间 10.1 , 101 rad/s，则滤波器的逼近效果比较理想。
>> G3=tf([1 36 126 84 9],[9 84 126 36 1]); w=logspace(-3,3);
s=fotf('s'); H=bode(s^-0.5,w); bode(H,'-',H2,'--')
如果用这些滤波器计算给定信号的 0.5阶积分，得出的结果如图 5-4所示。可见，得出的滤波效果与 Grünwald–Letnikov定义得出的理论值几乎完全一致。
·120·分数阶微积分学——数值算法与实现
图 5-3两个滤波器的 Bode图比较
>> t=0:0.01:10;
y=exp(-t).*sin(3*t+1); y0=glfdiff9(y,t,-0.5,5);
y1=lsim(H2,y,t); plot(t,y0,'-',t,y1,'--')
图 5-4 0.5阶积分的计算效果比较
如果时间区间增加到 t ∈ (0, 100)，则可以得出滤波器与理论值的时域响应拟合比较，如图 5-5所示。
>> t=0:0.01:100;
y=exp(-t).*sin(3*t+1); y0=glfdiff9(y,t,-0.5,3);
y1=lsim(H2,y,t); plot(t,y0,'-',t,y1,'--')
不过，这时的拟合效果在 t较大时产生很大误差，而 t很大对应于低频段，由此可知，误差产生的原因是低频处拟合不佳。结合图 5-3中的频域响应拟合，可以看出，拟合区域下限选择为 10.1 rad/s是不够的，应该扩展区间下限到 ω< 10.2 rad/s甚至更低，不过这样的要求在 Carlson滤波器的设计中是做不到的，所以需要允许用户自行选择拟合区间的更好算法。
图 5-5更大时间范围的 0.5阶积分拟合效果
5.2.3 Matsuda–Fujii近似
Matsuda和 Fujii提出了一种利用连分式技术近似无理传递函数 G(s)的整数阶传递函数的逼近方法[4，5]。遗憾的是，其原始算法的数学描述有误，这里只给出改正后的算法。
下面给出 Matsuda–Fujii滤波器算法的 MATLAB实现：
function G=matsuda_fod(G0,n,wb,wh)
arguments, G0,n,wb(1,1),wh(1,1){mustBeGreaterThan(wh,wb)}, end
if nargin==2, f=G0; w=n; n=length(w);
else
if isa(G0,'double'), s=fotf('s'); G0=s^G0; end
w=logspace(log10(wb),log10(wh),n); f=mfrd(G0,w);
·122·分数阶微积分学——数值算法与实现
end
v=abs(f(:).'); s=tf('s'); n=length(w);
for k=1:n, a(k)=v(k); v=(w-w(k))./(v-v(k)); end
G=a(n); for k=n-1:-1:1, G=a(k)+(s-w(k))/G; end
end
说明 5-1 Matsuda–Fujii滤波器设计函数
（1）该函数有下面 3种调用格式： G=matsuda_fod(f,ω)，其中需要给定原系统的频域响应数据。 G=matsuda_fod(γ,n,ωb,ωh)，在频率 ωb,ωh内对 sγ算子进行 n点近似。 G=matsuda_fod(G0,n,ωb,ωh)，在频率 ωb,ωh内对 G0进行 n点近似。
（2）由于 ωb与 ωh是选择频率边界内的内点，所以实际拟合的频率段要比指定的 ωb,ωh稍大些。
（3）如果想得到较好的拟合，则应该选择奇数 n，这样，实际滤波器的阶次等于 (n . 1)/2。
（4）虽然拟合的目标是频域响应的幅值拟合，但相位的拟合效果也比较好。
例 5-4试为 1/ √ s设计一个 Matsuda–Fujii滤波器，并评价其拟合效果。 解如果选择频率段 10.1 , 101 rad/s和 10.2 , 102 rad/s，则可以直接设计出 Matsuda–Fujii滤波器，得出的 Bode图比较如图 5-6所示。
>> s=fotf('s'); w=logspace(-3,3); G1=matsuda_fod(-0.5,9,1e-1,1e1) G2=matsuda_fod(-0.5,9,1e-2,1e2) %更大频率段的滤波器设计 H=bode(s^-0.5); bode(H,G1,'--',G2,':',{1e-3,1e3})
图 5-6 Bode图比较
用两种调用格式设计出的滤波器是完全一致的，与文献 [4]中给出的结果也是完全一致的。 0.0855s4 +4.876s3 + 20.84s2 + 13s +1 G1(s)=
s4 + 13s3 + 20.84s2 +4.876s +0.0855
0.04401s4 +8.142s3 + 58.85s2 + 30.76s +1
G2(s)=
s4 + 30.76s3 + 58.85s2 +8.142s +0.04401
从频域响应比较可见，低频拟合效果与 Carlson滤波器接近，可以说，低频拟合的问题尚未解决。好在 Matsuda–Fujii滤波器的频率段与阶次都是可以人为选定的，所以实际的频域响应拟合效果有很大的改善余地。例如，可以增加滤波器的阶次，并增大拟合的频率范围，则由下面语句直接设计 8阶滤波器（取 n = 17），这时得出的 Bode图如图 5-7所示。可见，这样设计的滤波器的效果有明显的改善，其 Bode图拟合的频率范围更宽，拟合效果也更好。
>> w=logspace(-3,3); G1=matsuda_fod(-0.5,9,1e-2,1e2);
G2=matsuda_fod(-0.5,17,1e-2,1e2), H=bode(s^-0.5);
bode(H,G1,'--',G2,':',{1e-3,1e3})
图 5-7 Bode图比较
有了滤波器，就可以将已知信号馈入该滤波器，得出该信号的 0.5阶积分，并和数值计算的精确解相比较，得出如图 5-8所示的比较结果。可以看出， G2滤波器与数值计算所得的曲线几乎重合，而 G1滤波器的拟合效果（虚线）比较差。
>> t=0:0.01:100; y=exp(-t).*sin(3*t+1);
y0=glfdiff9(y,t,-0.5,3); y1=lsim(G1,y,t); y2=lsim(G2,y,t);
plot(t,y0,'-',t,y1,'--',t,y2,':')
5.2.4拟合效果与滤波器参数选择的关系
前面介绍了各种滤波器的设计方法，也在时域响应角度比较了滤波器阶次与感兴趣频段对逼近效果的影响，得出了一些定性结论。这里将结合 Laplace变换的性质，更准确地解释这些滤波器参数对时域响应拟合的影响。
·124·分数阶微积分学——数值算法与实现
图 5-8时域响应拟合效果
因为频率与 s的关系为 s = jω，所以，由定理 5-1可知，若想在初始时刻（即 t比较小时）得到很好的拟合效果，则应该改善高频段的拟合，即提高感兴趣频率区域的上限 ωh；如果想得到更好的时域响应精度，则应该减小感兴趣频率区间的下限 ωb；若想得到更好的整体的时域响应精度，则应该使得中频段频域响应曲线更好地逼近理论值的斜线，而增加滤波器的阶次是一种有效的改善方法。增加滤波器的阶次当然会增大计算量，不过这里增加的计算量对当前的计算机软硬件系统而言并不会增加很大的负担。
5.3 Oustaloup滤波器近似
前面已经通过例子演示过，由于基于连分式的滤波器设计方法在频域响应拟合上不甚理想，尤其是连分式类方法不允许使用者自由选择合适的拟合频率段，使得这样的滤波器应用大打折扣。本节将介绍更实用的 Oustaloup滤波器。
5.3.1常规的 Oustaloup近似
假设感兴趣的频率段为 ωb,ωh，可以考虑用图 5-9中给出的一组折线逼近分数阶微积分算子幅频响应的直线特性。法国学者 Oustaloup教授等基于这样的想法提出了滤波器设计的方法 [6]，本书中称之为 Oustaloup滤波器。所有这些折线都是由整数阶的零点与极点生成的，使得幅频特性渐近线的斜率在 0dB/dec与 .20 dB/dec之间交替变化，这样频域响应本身会很好地逼近一条斜线。
根据上面的公式，可以总结出设计 Oustaloup滤波器的算法。
基于上述算法可以直接编写出 Oustaloup滤波器的 MATLAB函数如下：
function G=ousta_fod(gam,N,wb,wh)
arguments
gam(1,1) double, N(1,1) {mustBePositiveInteger}=9
wb(1,1) double=1e-4;
wh(1,1) double {mustBeGreaterThan(wh,wb)}=1e4;
end
if round(gam)==gam, G=tf('s')^gam;
else
k=1:N; wu=sqrt(wh/wb);
wkp=wb*wu.^((2*k-1-gam)/N); wk=wu^(2*gam/N)*wkp;
G=zpk(-wkp,-wk,wh^gam);
幅值（dB）
f
ω ′ ω ′ ω ′ ω ′
O
ω1 ω212 ω（rad/sω） N.1 N.1 ωNN
图 5-9 Oustaloup滤波器的分段折线逼近
·126·分数阶微积分学——数值算法与实现
end, end
该函数的调用格式为 G=ousta_fod(γ,N,ωb,ωh)，其中， γ为分数阶算子的阶次，N为滤波器的阶次， ωb与 ωh为用户选择的感兴趣频率段的下限与上限。一般情况下，所设计的滤波器 G的 Bode图在指定频段效果很好，但在该频率段之外效果不好。
说明 5-2 Oustaloup滤波器
（1）这里给出的滤波器实现避免了 ωbωh =1限制，两端均可以独立选择。
（2）如果 γ为整数，可以直接构造整数阶微积分器 G(s)= sγ。
（3）这里的 γ既可为正也可为负，分别表示微分与积分。
（4）γ的绝对值允许大于 1，例如 γ =3.7仍可以直接拟合，不过在实际应用中，建议将阶次范围限定在 .1 <γ< 1，其余的部分用整数阶微积分算子表示，例如， s3.7 = s3s0.7，或 s3.7 = s4s.0.3。
（5）尽管从图 5-9上看近似的幅频特性像折线，这些折线实际上是 Bode图的
渐近线，如果参数选择合适，实际的 Bode图本身应该很接近斜线。将 y(t)信号馈入 Oustaloup滤波器，则滤波器的输出信号为 RLDγ y(t)，即输入
0 t
信号的 Riemann–Liouville分数阶导数或积分。
例 5-5假设感兴趣的频率段为 ωb =0.01，ωh = 1000 rad/s，可以选择不同的阶次设计 Oustaloup滤波器。如果输入函数为 f(t)= e.t sin(3t + 1)，试设计 0.5阶积分的滤波器，并评价滤波器的近似效果。 解可以用下面的语句直接设计 5阶 Oustaloup滤波器 G1(s)和 Matsuda–Fujii滤波器 G2(s)：
>> G1=ousta_fod(-0.5,5,0.01,1000), G1a=zpk(G1)
G2=matsuda_fod(-0.5,11,0.01,1000), G2a=zpk(G2)
s=fotf('s'); G=s^-0.5; w=logspace(-4,5);
H=bode(G,w); bode(H,'-',G1,'--',G2,':')
设计出来的两个滤波器的零极点形式分别为
0.031623(s + 562.3)(s + 56.23)(s +5.623)(s +0.5623)(s +0.05623)
G1a(s)=
(s + 177.8)(s + 17.78)(s +1.778)(s +0.1778)(s +0.01778)
0.013865(s + 1792)(s + 71.46)(s +5.81)(s +0.503)(s +0.03427)
G2a(s)=
(s + 291.8)(s + 19.88)(s +1.721)(s +0.1399)(s +0.00558)
可以绘制出该滤波器的 Bode图，如图 5-10所示，同时绘制出分数阶算子 s.0.5的 Bode图，分数阶算子由 fotf()函数定义，该函数将在第 6章详细讨论。可以看出，对本例而言，Matsuda–Fujii滤波器的拟合带宽要宽一些。
还可以得出两个滤波器的滤波结果与 Grünwald–Letnikov 0.5阶积分的数值解，如
图 5-10两个滤波器的 Bode图拟合
图 5-11所示。可以看出，这样得出的分数阶积分不是很精确，两个滤波器的近似效果差不多。若想改善近似效果，建议增大感兴趣频域区间，并增加阶次。
>> t=0:0.001:pi; y=exp(-t).*sin(3*t+1);
y1=lsim(G1,y,t); y2=lsim(G2,y,t);
y0=glfdiff9(y,t,-0.5,3); plot(t,y0,t,y1,'--',t,y2,':')
图 5-11分数阶积分的近似效果
例 5-6现在考虑为 Heaviside函数设计 0.5阶导数的 Oustaloup滤波器，并假设时间范围很大，试观察 Oustaloup滤波器参数对导数计算精度的影响。 解 Heaviside函数的频率为零。现在选择感兴趣的频率段 (0.01, 1000) rad/s，并假设仿真的时间区间为 (0,20)s，就可以使用下面的语句设计 Oustaloup滤波器。这样，用两种方法得出分数阶导数曲线，如图 5-12所示。可以看出，当 t =2或 t很大时用滤波器的计算结果很不精确。
>> t=0:0.01:20;
y=ones(size(t)); G=ousta_fod(0.5,5,0.01,1000),
y1=lsim(G,y,t); y2=glfdiff9(y,t,0.5,3);
·128·分数阶微积分学——数值算法与实现
plot(t,y1,t,y2,'--'), ylim([0,1])
图 5-12 Heaviside函数的分数阶导数
选择考虑增加 Oustaloup滤波器的阶次，则由下面的语句可以重新计算分数阶导数，虽然在局部拟合上有所改进，在 t较大时误差仍然很大。这意味着只改变滤波器的阶次并不能改进大时间区间的近似效果。
>> G=ousta_fod(0.5,9,0.01,1000); y1=lsim(G,y,t);
plot(t,y1,t,y2,'--'), ylim([0,1])
因为时间 t很大时存在大误差，这意味着低频段 Oustaloup滤波器的拟合不好，所以应该降低 ωb的值，例如选择 ωb =0.001 rad/s，这时在新滤波器下分数阶导数的计算结果如图 5-13所示。可以看出已经得到了很好的计算精度。
>> G=ousta_fod(0.5,9,0.001,1000); y1=lsim(G,y,t);
plot(t,y1,t,y2,'--'), ylim([0,1])
图 5-13改进的计算结果
如果想进一步改进计算精度，则可以考虑选择更小的频率下限 ωb，并选择更高的滤波器阶次。
例 5-7从前面的例子可以看出，若想提高计算精度，则可以增大感兴趣的频率段，那么如何选择滤波器的阶次呢？ 解例 5-5演示过，如果感兴趣的频率段为 10.2 , 103 rad/s，则 5阶 Oustaloup滤波器效果已经很好了。一般情况下，如果频率段每增加 1个十倍频程，则应该至少增加阶次 1到 2阶。例如，如果感兴趣频率段选择为 10.4 , 104 rad/s，即频率段增加了 3个十倍频程，则需要设计一个 11阶的滤波器。一般情况下，如果计算机资源允许，不妨选择一个更高的阶次。不同阶次滤波器的 Bode图如图 5-14所示。
>> s=fotf('s'); G1=s^0.5;
w=logspace(-5,5); H=bode(G1,w);
G=ousta_fod(0.5,5,1e-4,1e4); G1=ousta_fod(0.5,7,1e-4,1e4);
G2=ousta_fod(0.5,9,1e-4,1e4); G3=ousta_fod(0.5,11,1e-4,1e4);
bode(H,G,'--',G1,':',G2,'-.')
图 5-14不同阶次下的滤波器
为了更好地演示 11阶 Oustaloup滤波器，给出如图 5-15所示的拟合效果，对幅频特性曲线局部放大，可以看出，这时的频域响应曲线几乎处处都是直线，而不是折线。
说明 5-3滤波器设计参数选择的建议
（1）如果时间 t较大时滤波近似效果不佳，则应该相应地减小低频下限 ωb，而初始时刻近似不佳时则应增大频率上限 ωh；
（2）不要害怕使用很高的滤波器阶次和很宽的频段，例如 (10.6 , 106)，N = 30并不会显著地增加运算时间。
5.3.2一种改进的 Oustaloup滤波器
从前面的例子可以看出， Oustaloup滤波器的近似效果在选择的频率段边界 ωh，ωb处不太理想。另外，滤波器分子的阶次与分母的阶次相同，传递函数不是严格正则系统，所以可以考虑使用改进的滤波器模型。
·130·分数阶微积分学——数值算法与实现
图 5-15 11阶 Oustaloup滤波器的频率响应近似效果
这样，这一滤波器的阶次 γ必须满足 γ ∈ (0, 1)。正常情况下，加权参数可以选择为 b = 10，d =9。对算法 5-3稍加变化就可以设计出改进的 Oustaloup滤波器。下面的 MATLAB函数可以直接用于改进 Oustaloup滤波器的设计：
function G=new_fod(r,N,wb,wh,b,d)
arguments
r(1,1) double, N(1,1){mustBePositiveInteger}=9
wb(1,1)double=1e-4; wh(1,1){mustBeGreaterThan(wh,wb)}=1e4;
b(1,1) double=10; d(1,1) double=9;
end
G=(d/b)^r*tf([d,b*wh,0],[d*(1-r),b*wh,d*r]);
G=G*ousta_fod(r,N,wb,wh);
end
该函数的调用格式为 G=new_fod(γ,N,ωb,ωh,b,d)，其中，参数 b，d可以省略，直接使用默认值。
例 5-8考虑例 5-5中的问题，选择 ωb =0.01 rad/s，ωh = 1000 rad/s，试观测新滤波器的近似结果，并比较分数阶导数的近似效果。 解由下面的语句可以设计这两个滤波器，并得出 Bode图的近似效果，如图 5-16所示。可见改进算法的拟合带宽更宽。
>> G1=ousta_fod(0.5,5,0.01,1000);
G2=new_fod(0.5,5,0.01,1000); zpk(G2)
s=fotf('s'); G0=s^0.5; w=logspace(-5,5,100);
H=bode(G0,w); bode(G1,'-',G2,'--',H,':')
图 5-16 Bode图比较
得出滤波器 G2(s)的零极点表达式为 60s(s+0.01778)(s+0.1778)(s+1.778)(s+17.78)(s+177.8)(s+1111)
G2(s)= (s+0.056)(s+0.56)(s+5.623)(s+56.23)(s+562.3)(s+2222)(s+0.00045)
计算得出的分数阶导数曲线如图 5-17所示。可以看出，改进的滤波器得出的拟合结果稍有改善。
>> t=0:0.001:pi; y=exp(-t).*sin(3*t+1);
y1=lsim(G1,y,t); y2=lsim(G2,y,t);
y0=glfdiff9(y,t,0.5,3);
plot(t,y1,t,y2,'--',t,y0,':'), axis([0,pi,-2 8])
例 5-9考虑下面的分数阶传递函数模型：
s +1
G(s)=
10s3.2 + 185s2.5 + 288s0.7 +1
试比较两个滤波器对分数阶传递函数模型的拟合效果。 解如果设计 0.2阶导数的滤波器，则得出的两个滤波器的 Bode图如图 5-18所示。可以由 FOTF工具箱中提供的 bode()重载函数直接绘制分数阶传递函数的 Bode图。还可以由下面的代码绘制出两个滤波器的 Bode图，如图 5-19所示。可以看出，改进滤波器得出的拟合结果明显好于 Oustaloup滤波器。
>> b=[1 1]; a=[10,185,288,1]; nb=[1 0]; na=[3.2,2.5,0.7,0];
G0=fotf(a,na,b,nb); w=logspace(-4,4,200); H=bode(G0,w);
·132·分数阶微积分学——数值算法与实现
图 5-17不同滤波器下的分数阶导数计算
s=zpk('s'); N=5; w1=1e-3; w2=1e3; b=10; d=9;
g1=ousta_fod(0.2,N,w1,w2); a1=g1;
g2=ousta_fod(0.5,N,w1,w2); g3=ousta_fod(0.7,N,w1,w2);
G1=(s+1)/(10*s^3*g1+185*s^2*g2+288*g3+1);
g1=new_fod(0.2,N,w1,w2); g2=new_fod(0.5,N,w1,w2);
g3=new_fod(0.7,N,w1,w2); bode(g1,a1,'--');
figure, G2=(s+1)/(10*s^3*g1+185*s^2*g2+288*g3+1);
bode(H,G1,'--',G2,':',w)
图 5-18 s
0.2的滤波器拟合
5.4分数阶传递函数的整数阶近似
5.4.1分数阶传递函数的高阶近似
如果遇到例 5-9中给出的分数阶传递函数模型，可以用该例中给出的方法得出高阶整数阶近似模型。但是，如果每次都这样用底层命令建立高阶模型是很烦琐的，所以可以编写一个 MATLAB函数进行这样的自动转换，得出高阶整数阶传递函数
图 5-19分数阶传递函数的 Bode图比较
模型。
基于上述的算法编写出 MATLAB函数进行所期望的转换，其中用到了很多 FOTF工具箱的函数，该函数应该置于 @fotf文件夹，其清单如下：
function Ga=high_order(G0,filter,wb,wh,N)
arguments
G0, filter='ousta_fod', wb(1,1) double=1e-3
wh(1,1) double {mustBeGreaterThan(wh,wb)}=1e3
N(1,1) {mustBePositiveInteger}=5
end
[n,m]=size(G0); F=filter;
for i=1:n, for j=1:m
if G0(i,j)==fotf(0), Ga(i,j)=tf(0);
else, G=simplify(G0(i,j)); [a na b nb]=fotfdata(G);
G1=pseudo_poly(b,nb,F,wb,wh,N)/pseudo_poly(a,na,F,wb,wh,N);
Ga(i,j)=minreal(G1);
end, end, end, end
%伪多项式的近似
function G1=pseudo_poly(a,na,filter,wb,wh,N)
G1=0; s=tf('s');
for i=1:length(a), na0=na(i); n1=floor(na0);
·134·分数阶微积分学——数值算法与实现
if na0>n1, g1=eval([filter '(na0-n1,N,wb,wh)']);
else, g1=1; end
G1=G1+a(i)*s^n1*g1;
end, end
函数的调用格式为 Ga =high_order(G0,filter,ωb,ωh,N)，其中， G0为分数阶传递函数矩阵对象（ FOTF对象），可以直接处理多变量的问题；输入变量 filter可以为 'ousta_fod'，'new_fod'或'matsuda_fod'；ωb与 ωh用于描述感兴趣频率段； N为每个分数阶算子的近似阶次。可以看出，该函数的结构是开放的，如果有其他滤波器设计函数，也可以仿照这样的结构直接嵌入。
例 5-10考虑例 5-9中给出的 FOTF对象，试得出整数阶的传递函数近似模型。 解可以用下面的语句直接输入 FOTF对象模型，并求出用 Oustaloup滤波器或改进 Oustaloup滤波器替代分数阶算子后的高阶传递函数近似模型，它们的频域响应拟合结果与图 5-19中给出的完全一致。
>> b=[1 1]; a=[10,185,288,1]; nb=[1 0]; na=[3.2,2.5,0.7,0];
G0=fotf(a,na,b,nb); N=4; w1=1e-3; w2=1e3; b=10; d=9;
G10=high_order(G0,'ousta_fod',w1,w2,N)
G20=high_order(G0,'new_fod',w1,w2,N)
w=logspace(-3,3); bode(G0), hold on; bode(G10,'--',G20,':')
这样得出的近似模型分别为
0.025119(s + 595.7)(s + 421.7)(s + 251.2)(s + 18.84)× (s + 13.34)(s +7.943)(s + 1)(s +0.5957)(s +0.4217)× (s +0.2512)(s +0.01884)(s +0.01334)
G10(s)=
(s + 595.7)(s + 507)(s + 154.7)(s + 40.89)(s + 18.78)(s +7.377)×
(s +2.041)(s +0.4418)(s +0.2511)(s +0.05445)(s +0.01334)×
(s +0.002349)(s2 +0.3824s +1.613)
0.020523(s + 1)(s +0.5957)(s +0.4217)(s +0.2512)(s +7.943)×
(s + 13.34)(s + 18.84)(s + 251.2)(s + 421.7)(s + 595.7)(s + 1389)×
(s + 2222)(s + 3704)(s +0.01884)(s +0.01334)(s +0.007943)×
(s +0.00063)(s +0.00045)(s +0.00018)
G20(s)=
(s + 3704)(s + 2325)(s + 1111)(s + 595.7)(s + 488.2)(s + 152.6)×
(s + 40.07)(s + 18.78)(s +7.358)(s +2.041)(s +0.4422)(s +0.2511)×
(s +0.05453)(s +0.01334)(s +0.007943)(s +0.002197)(s +0.00045)×
(s +0.0002201)(s +0.00018)(s2 +0.3762s +1.574)
这两个滤波器模型得出的阶跃响应曲线可以由下面的语句得出，如图 5-20所示。可以看出，两个滤波器的近似效果都不理想，原因还是低频段拟合不好。
>> t=0:0.01:20; y=step(G0,t); y1=step(G10,t);
plot(t,y,t,y1,'--')
对这样的问题仍然需要降低 ωb的值，比如令其为 10.4 rad/s，并适当提高滤波器的阶次，这样可以得出如图 5-21所示的阶跃响应曲线。可以看出，即使这里采用了更大
图 5-20阶跃响应比较
的时间响应区域，计算精度仍然是令人满意的。
图 5-21更长时间段的阶跃响应曲线比较
>> N=7; w1=1e-4; w2=1e3; G10=high_order(G0,'ousta_fod',w1,w2,N)
t=0:0.01:100; y=step(G0,t); y1=step(G10,t);
plot(t,y,t,y1,'--')
5.4.2基于模型降阶技术的低阶近似方法
顾名思义，模型降阶就是用一个低阶模型逼近高阶模型的行为。由于使用滤波器逼近分数阶算子，经常会得到阶次特别高的整数阶模型，所以，这里主要研究用一个低阶整数阶模型逼近分数阶系统的行为。
若想得到一个较好的降阶模型，则应该定义一个最优准则。比如，如果原始模
·136·分数阶微积分学——数值算法与实现
型和降阶模型用同样的信号激励，则由两个模型的输出之差可以定义误差信号，由这个误差信号再定义目标函数。根据这样的目标函数将要求解的问题转换成数值最优化问题。
由于式（ 5-4-2）存在延迟项，可以考虑用 Padé近似逼近延迟项。这样目标函数可以变成下面的范数求解形式，而最优化问题变换成
G(s) . δ
J = min δGr/m(s)2（5-4-4）
θ
上述问题是没有解析解的，只能用数值最优化技术求解原始问题。由参考文献 [8]给出的思路，可以引出模型最优降阶算法。
可以看出，从原来分数阶传递函数模型得到的整数阶近似模型阶次非常高，所以可以用最优降阶的方法得出效果满意的低阶模型 [8，9]。最优模型降阶的 MAT-LAB函数可以写成
function Gr=opt_app(G,r,k,key,G0)
arguments
G, r(1,1) {mustBePositiveInteger}
k(1,1) {mustBePositiveInteger}
key(1,1){mustBeMember(key,[0,1])}=0, G0=0
end
GS=tf(G); Td=totaldelay(GS);
GS.ioDelay=0; s=tf('s'); GS.InputDelay=0; GS.OutputDelay=0;
if nargin<5, G0=(s+1)^r/(s+1)^k; end
beta=G0.num{1}(k+1-r:k+1); alph=G0.den{1};
Tau=1.5*Td; x=[beta(1:r),alph(2:k+1)];
if abs(Tau)<1e-5, Tau=0.5; end
dc=dcgain(GS); if key==1, x=[x,Tau]; end
y=opt_fun(x,GS,key,r,k,dc);
x=fminsearch(@opt_fun,x,[],GS,key,r,k,dc);
alph=[1,x(r+1:r+k)]; beta=x(1:r+1); beta(r+1)=alph(end)*dc;
if key==0, Td=0; end
if key==1, Tau=x(end)+Td; else, Tau=0; end
Gr=tf(beta,alph,'ioDelay',Tau);
end
%目标函数计算
function y=opt_fun(x,G,key,r,k,dc)
ff0=1e10; a=[1,x(r+1:r+k)]; b=x(1:r+1); b(end)=a(end)*dc;
if key==1, tau=x(end);
if tau<=0, tau=eps; end, [n,d]=pade(tau,3); gP=tf(n,d);
else, gP=1; end
Ge=G-tf(b,a)*gP;
Ge.num{1}=[0,Ge.num{1}(1:end-1)]; [y,ierr]=geth2(Ge);
if ierr==1, y=10*ff0; else, ff0=y; end
end
%计算 H2范数
function [v,iE]=geth2(G)
G=tf(G); num=G.num{1}; den=G.den{1}; iE=0; v=0; n=length(den);
if abs(num(1))>eps
disp('System not strictly proper'); iE=1; return
else, a1=den; b1=num(2:length(num)); end
for k=1:n-1
if (a1(k+1)<=eps), ierr=1; return
else
aa=a1(k)/a1(k+1); bb=b1(k)/a1(k+1); v=v+bb*bb/aa; k1=k+2;
for i=k1:2:n-1,
a1(i)=a1(i)-aa*a1(i+1); b1(i)=b1(i)-bb*a1(i+1);
end, end, end, v=sqrt(0.5*v);
end
该函数的调用格式为 Gr =opt_app(G,r,m,key,G0)，其中 r和 m为用户指定的降阶模型分子与分母多项式的阶次； key表示在降阶模型中是否需要延迟项；如果有必要，用户也可以提供自己的初始模型 G0。
·138·分数阶微积分学——数值算法与实现
例 5-11考虑如下给出的分数阶传递函数模型
0.63 . 4
.2s
G(s)=
2s3.501 +3.8s2.42 +2.6s1.798 +2.5s1.31 +1.5
试找出其最优降阶模型并评价降阶模型与原模型的匹配程度。 解可以用下面的语句直接得出高阶近似模型：
>> b=[-2 -4]; nb=[0.63 0]; a=[2 3.8 2.6 2.5 1.5];
na=[3.501 2.42 1.798 1.31 0]; G=fotf(a,na,b,nb);
G1=high_order(G,'ousta_fod',1e-3,1e3,5); order(G1)
可以看出，这样得出的模型为 45阶的整数阶模型。现在考虑找到一个降阶模型，例如选择 r =4，m =5，则可以直接得出降阶模型，并得出高阶与降阶模型的 Bode图，如图 5-22所示。这里的相位图看似有很大差异，不过相位图是平行的相差 360.，所以，可以认为它们是一致的。
>> G2=opt_app(G1,4,5,0), zpk(G2)
w=logspace(-4,4); H=bode(G,w); bode(H,G1,'--',G2,':')
图 5-22 Bode图比较
这样得出的低阶近似模型为
.4.7911(s +0.2452)(s +0.01657)(s2 . 8.983s + 31.65)
G2(s)= (s + 140.9)(s +0.3799)(s +0.0173)(s2 +0.2303s +0.2479)
可以看出，这样得出的降阶模型的频域响应还是很接近原始模型的，而高频部分的拟合看起来和原始模型有很大的误差。不过不必太担心，因为这里的 Bode幅频特性使用的单位是 dB（分贝），当增益小于 .20 dB时，其实际的倍数小于 .20 dB，即 10.1，而 .50 dB大约为 10.2.5，故二者都是很小的增益，并没有太显著的差异，至少没有图 5-22中看起来那么大的差异。相频特性拟合中，几条曲线基本上相差的是 360.和 720.，这里给出的是手工修改后的相频曲线。可见，实际分数阶系统的频域响应曲线与高阶整数阶拟合模型几乎完全一致。
还可以求出各个模型的阶跃响应，如图 5-23所示。可见，降阶模型和原始模型还是
很接近的。 >> t=0:0.01:10; y=step(G,t); y1=step(G1,t); y2=step(G2,t); plot(t,y,t,y1,'--',t,y2,':')
图 5-23阶跃响应比较
例 5-12考虑下面的分数阶传递函数模型，试得出降阶模型并比较其效果。
5
G(s)=
s2.3 +1.3s0.9 +1.25
解首先先得出原始模型的高阶整数阶近似模型，然后采用不同的降阶方法，可以得出如图 5-24所示的阶跃响应比较。这里的阶跃响应曲线是 step()函数自动绘制的。
>> w1=1e-3; w2=1e3; N=9;
s=fotf('s'); G=5/(s^2.3+1.3*s^0.9+1.25);
G0=high_order(G,'ousta_fod',w1,w2,N); zpk(G0)
G1=opt_app(G0,1,2,0), G2=opt_app(G0,2,3,0)
G3=opt_app(G0,3,4,0), step(G0,G1,'--',G2,':',G3,'-.')
图 5-24分数阶传递函数的阶跃响应比较
上面的语句首先得出一个 19阶的高阶近似模型（其数学形式从略），还可以得出最
·140·分数阶微积分学——数值算法与实现
优降阶模型为
.2.045s +7.654 .0.5414s2 +4.061s +2.945
G1(s)= ,G2(s)=
s2 +1.159s +1.917 s3 +0.9677s2 +1.989s +0.7378 .0.2592s3 +3.365s2 +4.9s +0.3911
G3(s)=
s4 +1.264s3 +2.25s2 +1.379s +0.09797可以看出，二阶模型 G1(s)与 3阶模型 G2(s)的拟合效果不甚理想， 4阶与 5阶近似模型的响应很接近分数阶模型。在一些特定的例子中，比如在例 5-9的问题中，如果降阶模型的阶次选择得过低，则逼近效果会很差，所以得到降阶模型后应该检验一下才能真正使用。
5.5无理分数阶模型的近似
在分数阶控制系统中，有的时候系统的某个组成部分不能很好地用有理分数阶传递函数直接描述。例如，若传递函数存在 (as + b)/(cs + d) α这样的项，而 α不是整数，则该项只能展开成无穷项的级数。如果存在 asβ + b α形式的项则更难以处理。本节将介绍这类无理传递函数的近似方法。
5.5.1隐式无理模型的近似
考虑一种简单的无理传递函数
11
G(s)= (τs + 1)ν = τ.ν (s + 1/τ)ν（5-5-1）
证明由传递函数模型可得 Y (s)= G(s)U(s)= G(s + 1/τ)U(s)，则 Y (s . 1/τ)= G(s . 1/τ)U(s . 1/τ)= τ.νU(s . 1/τ)/sν。取 Laplace反变换，则利用 Laplace变换的平移性质 L [e.atf(t)] = F (s + a)，取 a = .1/τ，得出
L [Y (s . 1/τ)] = τ.νD.ν[et/τ u(t)]（5-5-3）
再利用平移性质，则得出下式，定理得证。
et/τ y(t)= τ.ν D.ν[et/τ u(t)]（5-5-4）
t
× ×
（a）框图表示 （b）Simulink实现（c5mimp.slx）
图 5-25 隐式无理模型建模
5.5.2 频域响应近似方法
说明 5-4隐式无理微分方程
（1）这里给出了一种近似方法，但该方法不能得出传递函数近似模型，无法比较频域响应拟合效果。
（2）利用这里介绍的方法，可以用如图 5-25（a）所示的框图描述近似方法。基于这里给出的框图，可以直接建立 Simulink仿真模型，如图 5-25（b）所示。其中， 1/sν模块可以由 Oustaloup滤波器实现。
（3）由于 Simulink模型中存在 t.ν项，所以若 ν> 0，仿真过程中不宜将 t的初值设置为 0，可以将其设置成微小的值，如 0.0001。
当然，并不是所有的无理分数阶模型都能用式（ 5-5-1）的简单形式表示，所以需要探讨一般无理传递函数模型的近似方法。
如果测得或计算出某个模型的一些频率响应点的数据，则可以由 MATLAB提供的 invfreqs()函数获得整数阶近似模型。一般情况下，这样的模型可以很好地逼近给出的频域响应数据。该函数的调用格式为 G=invfreqs(H,w,m,n)，其中， H为频域响应数据，w为频率点向量，m和 n分别为预期的分子与分母的阶次。
·142·分数阶微积分学——数值算法与实现
例 5-13考虑文献 [11]给出的分数阶定量反馈理论（quantitative feedback theory， QFT）控制器模型 0.96 1.76
s +0.011 8.8×10.5s +1 1
Gc(s)=1.8393
2
s 8.096×10.5s +1 (1+ s/0.29)
试找出较好的有理近似模型。 解这样的控制器结构很复杂，只能借助于比较好的有理传递函数模型逼近。由于 MATLAB控制系统工具箱中的 frd()函数只能处理整数阶模型，该函数返回一个结构体变量，可以对其成员变量 ResponseData进行处理，获得其非整数次幂，最终得出无理传递函数的精确频域响应数据。有了频域响应数据，就可以调用 invfreqs()函数获得拟合模型。对本例而言，可以先选择感兴趣的频域段 ω ∈ 10.4 , 102 rad/s，再输入下面的 MATLAB命令：
>> w=logspace(-4,2); G1=tf([1 0.011],[1 0]); F1=frd(G1,w);
G2=tf([8.8e-5 1],[8.096e-5 1]); F2=frd(G2,w);
s=tf('s'); G3=1/(1+s/0.29)^2; F3=frd(G3,w); F=F1;
h1=F1.ResponseData; h2=F2.ResponseData; h3=F3.ResponseData;
h=1.8393*h1.^0.96.*h2.^1.76.*h3; F.ResponseData=h;
[n,d]=invfreqs(h(:),w,4,4); H1=zpk(tf(n,d))
这样得出的近似控制器为
.1.5×10.10(s.3.913×104)(s+2.635×104)(s+0.01071)(s+0.0006019)
H1(s)= (s +0.2922)(s +0.2878)(s +0.0007664)(s +1.284×10.7)
从结果可以看出，有的稳定零极点对相距很近，所以可以通过最小实现技术将它们对消，得到更低阶次的模型：
>> H1a=minreal(H1,1e-3), H1b=minreal(H1,1e-2)
这样得出的更低阶近似为
.1.5002×10.10(s . 3.913×104)(s +2.635×104)(s +0.01071)
H1a(s)=
(s +0.2922)(s +0.2878)(s +1.284×10.7) .1.5002×10.10(s . 3.913×104)(s +2.635×104)
H1b(s)=
(s +0.2922)(s +0.2878)
对于无理系统拟合的问题还可以尝试使用 Matsuda–Fujii滤波器。用下面语句直接设计这样的滤波器，不过该滤波器 H2(s)是不稳定的模型，说明 Matsuda–Fujii设计方法并不能保证滤波器的稳定性，使用时需要慎重。
>> w=logspace(-4,2,11); F1=frd(G1,w); F2=frd(G2,w); F3=frd(G3,w);
h1=F1.ResponseData; h2=F2.ResponseData; h3=F3.ResponseData;
h=1.8393*h1.^0.96.*h2.^1.76.*h3; H2=zpk(matsuda_fod(h,w))
选择更大一些的频率段 10.6 , 104 rad/s验证得出的模型，就可以得出原控制器与各个低阶近似模型的 Bode图，如图 5-26所示。理论值用点线表示，可见各个幅频特性的逼近很理想，相频特性相差 360.，
（这里做了手工平移）也是比较理想的。
>> w=logspace(-6,4,200); F1=frd(G1,w); h1=F1.ResponseData;
F2=frd(G2,w); h2=F2.ResponseData; F3=frd(G3,w);
h3=F3.ResponseData; h=1.8393*h1.^0.96.*h2.^1.76.*h3;
F=F1; F.ResponseData=h; bode(H1a,'--',H1b,F,':',H1,'-.',w)
图 5-26频域响应比较
从比较结果看，最小实现模型 H1a(s)与 4阶控制器 H1(s)甚至与原无理传递函数模型很接近，而二阶近似模型 H1b(s)的效果很差，这意味着原来的无理传递函数模型至少需要 3阶传递函数近似。上面得出的 H1a(s)是通过最小实现得出的，其实更应该由 invfreqs()函数直接得出：
>> w=logspace(-4,2); F1=frd(G1,w); F2=frd(G2,w); F3=frd(G3,w);
h1=F1.ResponseData; h2=F2.ResponseData; h3=F3.ResponseData;
h=1.8393*h1.^0.96.*h2.^1.76.*h3; [n,d]=invfreqs(h(:),w,3,3);
H3=zpk(tf(n,d))
直接得出的 3阶模型为
.1.4886×10.10(s . 3.93×104)(s +2.644×104)(s +0.01009)
H3(s)=
(s +0.3028)(s +0.2768)(s +1.103×10.5)
例 5-14试用隐式模块逼近例 5-13中给出的 QFT控制器。 解对 QFT控制器传递函数进行重组，可以写成如下形式：
0.96
1.8393/0.0110.96 s
Gc(s)= +1 (8.8×10.5 s+1)1.76(8.096×10.5 s+1).1.76
2
s0.95 (1+s/0.29)0.011
由这个数学形式可以直接构造出如图 5-27所示的 Simulink仿真模型，其整数阶部分由传递函数模型给出，其中，分子编辑框 num填写 1.8393/0.011^0.96，分母编辑框
·144·分数阶微积分学——数值算法与实现
den填写 conv([1/0.29 1],[1/0.29 1])。
图 5-27 QFT控制器的 Simulink模型（c5mqft.slx）
5.5.3 Charef近似
Charef近似技术是对下面无理传递函数模型的一种有效的近似方法[4，12]： H(s)= (1+ s1/pT)α（5-5-6）
可以建立起下面的算法设计 Charef滤波器。
根据上面的算法可以编写出 Charef滤波器设计的 MATLAB函数：
function G=charef_fod(alpha,pT,delta,wM)
arguments
alpha(1,1) double, pT(1,1) double
delta(1,1) double, wM(1,1) double
end
p0=pT*10^(delta/20/alpha);
a=10^(delta/10/(1-alpha)); b=10^(delta/10/alpha);
n=ceil(log(wM/p0)/log(a*b)); ii=1:n; p=p0*(a*b).^ii;
p=[p0 p]; z=a*p(1:n); K=prod(p)/prod(z); G=zpk(-z,-p,K);
end
该函数的调用格式为 G=charef_fod(α,pT,δ,ωM)，该函数的输入变量与对
应的数学公式是一致的，返回变量则是设计出来的 Charef滤波器。
例 5-15试设计无理传递函数 G(s) = 1/(1+0.2s)0.5的 Charef滤波器。
解可以选择误差容限 δ =1 dB与频率段上限 ωM = 1000 rad/s。这样可以由下面
的语句直接设计出 Charef滤波器：
>> G1=charef_fod(0.5,1/0.2,1,1000)
设计出来的 Charef滤波器为
99.763(s +9.976)(s + 25.06)(s + 62.95)(s + 158.1)(s + 397.2)(s + 997.6)
G1(s)= (s+6.295)(s+15.81)(s+39.72)(s+99.76)(s+250.6)(s+629.5)(s+1581)
此外，由前面介绍的频域响应拟合方法也可以求取原模型的整数阶近似模型：
>> w=logspace(-2,3); G0=tf(1,[0.2 1]);
H=frd(G0,w); h=H.ResponseData; h=h.^0.5;
[n,d]=invfreqs(h(:),w,5,5); G2=zpk(tf(n,d))
得到的近似模型为
0.01541(s + 8899)(s + 1121)(s + 340.2)(s + 91.78)(s + 17.16)
G2(s)= (s + 2412)(s + 609.2)(s + 183.6)(s + 41.34)(s +7.41)
事实上，前面介绍的 Carlson滤波器也适合于近似无理传递函数模型，可以通过下
面的语句直接设计该滤波器：
>> s=tf('s'); G=1/(1+0.2*s); H1=carlson_fod(0.5,G,2); zpk(H1)
设计出的滤波器模型为
0.11111(s+165.8)(s+20)(s+8.52)(s+6.667)(s+6.666)(s+5.662)(s+5)6
H1(s)=
(s + 42.74)(s + 12.1)(s +6.678)(s +4.887)(s2 + 13.32s + 44.37)× (s2 + 10.33s + 26.67)(s2 +9.859s + 24.31)(s2 + 10.08s + 25.42)
无理传递函数与 Charef滤波器的 Bode图比较如图 5-28所示。可以看出， Charef滤
波器的频域响应也很令人满意。 >> w=logspace(-3,5); M=tf(1,[0.2 1]); H=frd(M,w); h=H.ResponseData; H.ResponseData=h.^0.5; bode(H,G1,'--',G2,':',H1,'-.'), xlim([1e-3 1e5])
·146·分数阶微积分学——数值算法与实现
图 5-28频域响应比较
无理传递函数 G(s)的阶跃响应可以通过数值 Laplace反变换的方法求出，该阶跃响应与 Charef滤波器和频率响应拟合模型阶跃响应的比较如图 5-29所示，其中，函数 INVLAP_new()是基于数值 Laplace变换的求解函数，第 6章将详细介绍。
>> [t,y]=INVLAP_new('1/(1+0.2*s)^0.5',0,3,1001,0,'1/s');
y1=step(G1,t); y2=step(G2,t); y3=step(H1,t);
plot(t,y,t,y1,'--',t,y2,':',t,y3,'-.')
图 5-29阶跃响应比较
例 5-16考虑下面带有两个无理因子的模型，试得出较好的 Charef滤波器：
1
G(s)=
0.60.3
(1 + s/1.6)(1+ s/6.2)
解原始模型可以认为是两个无理传递函数的乘积，可以得出 Charef滤波器。
>> G11=charef_fod(0.6,1.6,1,1000);
G12=charef_fod(0.3,6.2,1,1000); G1=G11*G12; zpk(G1)
这样得出的近似模型为
12018(s +3.447)(s +8.997)(s + 12.64)(s + 23.48)× (s + 37.85)(s + 61.3)(s + 113.3)(s + 160)(s + 339.2)× (s + 417.6)(s + 1015)(s + 1090)
G1(s)=
(s +1.938)(s +5.06)(s +9.1)(s + 13.21)(s + 27.24)×
(s + 34.47)(s + 81.55)(s + 89.97)(s + 234.8)(s + 244.1)×
(s + 613)(s + 730.8)(s + 1600)(s + 2188)
选择感兴趣频率段为 ω ∈ 10.2 , 103 rad/s，还可以得出频域响应近似模型：
>> G01=tf(1,[1/1.6 1]); G02=tf(1,[1/6.2 1]);
w=logspace(-2,3); H1=frd(G01,w); h1=H1.ResponseData;
H2=frd(G02,w); h2=H2.ResponseData; H=H1;
h=h1.^0.6.*h2.^0.3; [n,d]=invfreqs(h(:),w,5,5);
H.ResponseData=h; G2=zpk(tf(n,d))
这样得出的模型为
0.00011629(s +4.456×104)(s + 1348)(s + 358.2)(s + 81.78)(s +7.508)
G2(s)= (s + 1579)(s + 407.5)(s + 97.05)(s + 10.91)(s +2.266)
如果用两个频域响应拟合模型串联的形式逼近原始模型，则可以得出第 3个整数阶传递函数的近似模型：
>> h1=h1.^0.6; h2=h2.^0.3; [n,d]=invfreqs(h1(:),w,5,5);
g=tf(n,d); [n,d]=invfreqs(h2(:),w,5,5); G3=zpk(tf(n,d)*g)
这样得出的模型为
0.00029552(s +1.066×104)(s + 6244)(s + 1122)(s + 1016)× (s + 321.9)(s + 319.3)(s + 87.85)(s + 76.86)(s + 17.36)(s +9.993)
G3(s)=
(s + 2975)(s + 2079)(s + 709.9)(s + 526.2)(s + 223.3)(s + 143.8)× (s + 55.69)(s + 24.43)(s + 10.66)(s +2.536)
这 3个模型与原始无理传递函数模型的 Bode图比较如图 5-30所示。可以看出，分段拟合模型 G1(s)的效果最差。
图 5-30频域响应拟合的比较
>> w=logspace(-4,5); H1=frd(G01,w); h1=H1.ResponseData;
·148·分数阶微积分学——数值算法与实现
H2=frd(G02,w); h2=H2.ResponseData; H=H1;
h=h1.^0.6.*h2.^0.3; H.ResponseData=h;
bode(H,G1,'--',G2,':',G3,'-.',w)
各个模型的阶跃响应曲线如图 5-31所示。可以看出， 3个近似模型的效果都比较好。当然，Oustaloup滤波器方法效果最好，且有进一步改进的余地。
>> sys='1/(1+s/1.6)^0.6/(1+s/6.2)^0.3';
[t,y]=INVLAP_new(sys,0,3,1001,0,'1/s');
y1=step(G1,t); y2=step(G2,t); y3=step(G3,t);
plot(t,y,t,y1,'--',t,y2,':',t,y3,'-.')
图 5-31阶跃响应的比较
例 5-17重新考虑例 5-16中的无理模型，试利用 Simulink近似该模型。 解显然，原模型可以由两个隐式模型串联搭建，由此可以构造如图 5-32所示的 Simulink模型。在该模型下可以得出如图 5-33所示的阶跃响应曲线，从曲线看不出与数值 Laplace变换结果的区别。
图 5-32 Simulink仿真模型（c5msim1.slx）
>> sys='1/(1+s/1.6)^0.6/(1+s/6.2)^0.3';
[t,y]=INVLAP_new(sys,0,3,1001,0,'1/s'); N=15;
ww=[1e-4 1e4]; [t1,~,y1]=sim('c4msim1',3);
plot(t,y,t1,y1,'--')
5.5.4复杂无理模型的最优 Charef滤波器设计
现在回顾一下定义 5-8中的 Charef滤波器，该滤波器中的比值为 zi/pi = a与 pi+1/zi = b，i =0, 1, ··· ,n . 1，其中 a与 b为固定的数值。如果比值不再固定，则它
图 5-33阶跃响应的比较
们可以在某种性能指标下寻优，引入最优的 Charef滤波器[13]的概念。
选择一组频率点 ω1,ω2, ··· ,ωN，并计算出无理传递函数的频域响应数据 ha，最优化的目标是得到一组决策变量 x = c, wu0 , k, y0,y1, ··· ,yn的值，使得选定的目标函数最小化。这里最优化的目标函数为
N
.
w2
f(x)= w1||.h.(jωi)||∞ +.h.(jωi)+
N
i=1
N
.
w3||.∠h.(jωi)||∞ + wN 4 .∠h.(jωi)（5-5-14）
i=1
式中，wi为加权系数；h(jω)为滤波器模型的频域响应，且 .h.(jωi)= |h(jωi)|.|ha(jωi)|（5-5-15）
·150·分数阶微积分学——数值算法与实现
.∠h.(jωi)= ∠h(jωi) . ∠ha(jωi)（5-5-16）
如果选择无穷范数，则表示求向量所有元素绝对值的最大值。
基于上述的算法，可以编写出最优 Charef滤波器设计的 MATLAB函数：
function Ga=charef_opt(wr,n,G,wt,a,wc)
arguments, wr(:,1), n(1,1){mustBePositiveInteger}
G(:,1), w(1,4)=[1 1 1 1], a(1,1)=0.5, wc(1,1)=1
end
e=1; x0(1)=10^(e/(10*a*(1-a)));
x0(2)=10^(e/(20*(1-a)))*wc*10^(e/(20*a));
x0(3:n+3)=ones(1,n+1)/10^(e/(20*(1-a)));
ff=optimset; ff.Display='iter'; A=[];B=[];Aeq=[];Beq=[];CF=[];
xm=[1 0 zeros(1,n+1)]; xM=[3 10 ones(1,n+1)];
x1=fmincon(@charef_obj,x0,A,B,Aeq,Beq,xm,xM,CF,ff,wr,G,wt);
c=x1(1); wu0=x1(2); wu=wu0*c.^[0:n]; p=wu.*x1(3:end);
z=wu.^2./p; k=prod(p)/prod(z); Ga=zpk(-z',-p',k);
end
%计算目标函数值的子函数
function f=charef_obj(x,wr,G,wt)
n=length(x)-2;
c=x(1); wu0=x(2); wu=wu0*c.^(0:(n-1)); p=wu.*x(3:end);
z=wu.^2./p; k=prod(p)/prod(z); Ga=zpk(-z',-p',k);
Ga1=frd(Ga,wr); h=Ga1.ResponseData; Ga_fr=h(:);
f=wt(1)*norm(abs(angle(Ga_fr)-angle(G)),inf)+...
wt(2)*sum(abs(angle(Ga_fr)-angle(G)))/length(wr)+...
wt(3)*norm(abs(abs(Ga_fr)-abs(G)),inf)+...
wt(4)*sum(abs(abs(Ga_fr)-abs(G)))/length(wr);
end
其中，这样的最优化目标函数可以用 MATLAB子函数描述出来， x为决策变量向量。除了 x之外，其他变量都是附加变量，都应该以附加变量的形式传递给求解函数 fmincon()。这里的附加变量包括： a为阶次的初值； ωc为一阶系统时间常数 ωT的初值； ωr为频率点构成的向量； G为频率点 ωr处原无理传递函数的频域响应数据；wt为加权系数 w1，w2，w3与 w4构成的向量，这些加权值的默认值都是 1。可以看出，在实际寻优过程中 ωc和 a的值并不是很重要，它们就是数值最优化的初值，使用默认值就可以了。该函数调用了 MATLAB函数 fmincon()求最优化问题的数值解。
该函数调用格式为 Ga =charef_opt(ωr,n,G,wt,a,ωc)，其中， ωr与 G提供无理传递函数的频域响应信息，n为期待的最优 Charef滤波器的阶次。 说明 5-5最优 Charef滤波器
（1）感兴趣频率段由频率向量 ωr给出，本算法兼顾幅值与相位数据的拟合。
（2）频率 ωm的选择对最终结果并没有太大影响，实际的频率下界由受控对象模型的转折点确定。
（3）原始模型的静态增益必须为 1，否则应该事先提取出来。
（4）在实际应用中，最好将 c的上限设置为 3，这样零极点对 (zi,pi)可能不会交叉重叠。
（5）选择阶次 n的建议：首先选择一个初始整数 n0，得出最优 Charef滤波器并显示拟合误差，再选择阶次 n0 +1，如果拟合误差显著降低则进一步选择更高的阶次，否则接受 n0。
例 5-18试为例 5-16的模型设计最优的 Charef滤波器。
解选择一个感兴趣的频率段 ω = 10.2 , 103 rad/s，就可以计算出无理传递函数的频域响应数据；再选择最优 Charef滤波器的预期阶次为 6，就可以由下面语句设计出最优 Charef滤波器：
>> G01=tf(1,[1/1.6 1]); G02=tf(1,[1/6.2 1]);
w=logspace(-2,3); H1=frd(G01,w); h1=H1.ResponseData;
H2=frd(G02,w); h2=H2.ResponseData; H=H1;
h=h1.^0.6.*h2.^0.3; h=h(:); G1=charef_opt(w,6,h)
得出的最优 Charef滤波器模型为
·152·分数阶微积分学——数值算法与实现
0.00096744(s +6.692)(s + 14.16)(s + 46.48)(s + 126.5)× (s + 299.6)(s + 929)(s + 3537)
G1(s)= (s+2.274)(s+7.892)(s+17.67)(s+47.71)(s+148)(s+350.5)(s+676.4)
两个模型的 Bode图比较如图 5-34所示。可以看出，这时得出的 6阶滤波器模型的逼近效果与例 5-16得出的 17阶滤波器的效果相仿。
>> H.ResponseData=h; bode(H,G1,'--',w)
图 5-34频域响应比较
还可以指定一个更大的频率段 10.2 , 106 rad/s，从而设计出一个 16阶的最优 Charef滤波器，可以得到指定范围的完美拟合，如图 5-35所示，在拟合结果上几乎看不出任何区别。
>> w=logspace(-2,7,100); H1=frd(G01,w); h1=H1.ResponseData;
H2=frd(G02,w); h2=H2.ResponseData; H=H1;
h=h1.^0.6.*h2.^0.3; h=h(:); H1.ResponseData=h;
G2=charef_opt(w,15,h), bode(H1,G2,'--',w)
图 5-35更大区间的频域响应拟合比较
原始无理传递函数与滤波器模型的阶跃响应曲线如图 5-36所示，两个滤波器的响应很接近，但与原始模型之间稍有差异。
>> f='1/(1+s/1.6)^0.6/(1+s/6.2)^0.3/s';
[t,y]=INVLAP_new(f,0,3,1001);
step(G1,G2,'--',t); hold on, plot(t,y,':'), hold off
图 5-36阶跃响应比较
例 5-19现在考虑一个更复杂的无理传递函数模型：
0.3
(1 + s/6.2)
G(s)=
0.6
(1 + s/1.6)(s0.7/3 + 1)0.2
试设计最优 Charef滤波器并分析其效果。 解在频率段 10.4 , 105 rad/s内产生一个向量，计算原始无理传递函数的频域响应，然后设计出最优 Charef滤波器：
>> G01=tf(1,[1/1.6 1]); G02=tf([1/6.2 1],1);
w=logspace(-4,5,200); H1=frd(G01,w); h1=H1.ResponseData;
H2=frd(G02,w); h2=H2.ResponseData; H=H1;
G3=fotf([1/3 1],[0.7 0],1,0); H3=bode(G3,w);
h3=H3.ResponseData; h=h1.^0.6.*h2.^0.3.*h3.^0.2;
H1.ResponseData=h; G1=charef_opt(w,11,h(:));
bode(H1,G1,{1e-4,1e5})
原无理传递函数与整数阶拟合模型的 Bode图如图 5-37所示。可以看出，拟合效果还是很理想的。利用数值 Laplace反变换技术可得出无理系统的阶跃响应曲线，同时也可绘制近似模型 G1(s)的阶跃响应曲线，如图 5-38所示。可以看出，两个系统的响应是很接近的。
>> f='(1+s/6.2)^0.3/(1+s/1.6)^0.6/(s^0.7/3+1)^0.2';
[t,y]=INVLAP_new(f,0,3,1000,0,'1/s'); step(G1,t); line(t,y)
·154·分数阶微积分学——数值算法与实现
图 5-37频域响应比较
图 5-38阶跃响应的比较
设计出来的 12阶最优 Charef滤波器为
0.0027751(s +3.591)(s +8.44)(s + 23.97)(s + 64.13)×
(s + 179.1)(s + 488.1)(s + 1351)(s + 3681)(s +1.019×104)×
(s +2.772×104)(s +7.561×104)(s +2.515×105)
G1(s)=
(s +1.757)(s +5.638)(s + 14.97)(s + 42.22)(s + 114)×
(s + 315.6)(s + 860.2)(s + 2380)(s + 6487)×
(s +1.799×104)(s +4.973×104)(s +1.128×105)
5.6离散滤波器近似
前面介绍的是分数阶与无理系统的连续滤波器逼近。本节将介绍这些环节的
离散滤波器近似方法。
假设输入信号为 u(n)，则经过滤波器后的输出信号可以由差分方程表示为 y(k)= .a1y(k . 1) . a2y(k . 1) .· ··. amy(k . m)+
b1u(k)+ b2u(k . 1) + ··· + bn+1u(k . n)（5-6-2）
在实际应用中，经常使用的离散滤波器有下面两种特殊形式：
5.6.1 FIR滤波器逼近
对可以近似分数阶系统的离散滤波器来说，也可以考虑引入 FIR滤波器 [15]或 IIR滤波器 [16]的形式对其近似。利用 MATLAB滤波器设计工具箱中的 filt()函数，由 Ivo Petrá.教授开发的分数阶微积分器的 FIR滤波器设计函数 [17]。其核心部分清单（经过与本书风格一致的改写）如下：
function H=dfod2(n,T,r) %这里根据需要对程序结构进行了修改
arguments
n(1,1) {mustBePositiveInteger}
T(1,1) {mustBePositive}, r(1,1) double
end
if r>0
bc=cumprod([1,1-((r+1)./[1:n])]); H=filt(bc,T^r,T);
elseif r<0
bc=cumprod([1,1-((-r+1)./[1:n])]); H=filt(T^(-r),bc,T);
end, end
其中，n为期望的滤波器阶次，T为滤波器的采样周期，r为所需的导数阶次，用该函
·156·分数阶微积分学——数值算法与实现
数可以直接设计出滤波器 H。若某信号经过滤波器 H，则可以得出该信号的 r阶导数，经过滤波器得出的输出信号可以由 lsim()函数计算出来。
例 5-20试用离散 FIR滤波器对例 5-5中给出的信号求取分数阶微分。
解用下面的语句可以直接设计出滤波器，并绘制出 Bode图，如图 5-39所示。由得出的频域响应看， n = 100明显优于 n = 50时的滤波器，但显然远不如前面介绍的用 Oustaloup算法设计的连续滤波器。
>> t=0:0.001:pi;
y=exp(-t).*sin(3*t+1); y3=glfdiff(y,t,0.5);
G1=dfod2(50,0.001,0.5); G2=dfod2(100,0.001,0.5);
bode(G1,G2,'--')
图 5-39不同 (n, a)组合下的数值微分器频域响应比较
在这些滤波器下还可以得出给定信号的分数阶微分，如图 5-40所示，然而计算精度都很不理想，不能真正使用。
图 5-40数值微分器的微分计算
>> y1=lsim(G1,y,t); y2=lsim(G2,y,t);
plot(t,y1,t,y2,'--',t,y3,':'), ylim([-1 8])
5.6.2 IIR滤波器逼近
FIR滤波器的阶次较高，所以可以考虑采用 IIR滤波器生成分数阶微积分信号。从连续系统角度看，分数阶微积分可以用 Laplace变换写成 s±γ。对其进行离散化，就可以引入变换函数 s = wz.1近似分数阶微分或积分运算。
文献[16]给出了一种选择变换函数的新方法，该方法兼顾了 Simpson积分算法和梯形积分算法，可以比较好地实现离散分数阶算子滤波器。
对高阶微积分来说，可以先进行整数阶微积分，对其结果再利用各种滤波器进行分数阶微积分。文献 [16]还给出了基于连分式的滤波器设计方法和求解函数。不过，该函数是基于 Maple的连分式函数计算的，在当前版本 MATLAB下已经不能使用；最新 MATLAB版本下的 MuPAD对连分式的支持也不正常，所以这里将其核心部分用 Padé近似取代，修改后的函数如下：
function H=iir_pade(r,a,n,T)
arguments, r(1,1) double, a(1,1){mustBeInRange(a,0,1)}=0
n(1,1){mustBePositiveInteger}=3
T(1,1){mustBePositive}=0.01
end
syms x, b=(3+a-2*sqrt(3*a))/(3-a);
k0=6*b/T/(3-a); f=((1-x^2)/(1+b*x)^2)^r;
c=taylor(f,x,'Order',2*n); c=sym2poly(c);
c=c(end:-1:1); [N,D]=padefcn(c,n-1,n);
H=k0^r*tf(N(end:-1:1),D(end:-1:1),'Variable','z^-1','Ts',T);
end
其中， r0为阶次， a ∈ [0, 1]为加权系数， n为滤波器阶次， T为采样周期，这样可以由
·158·分数阶微积分学——数值算法与实现
该函数直接设计出离散的滤波器模型 H。下面将通过例子演示 IIR滤波器的设计及分数阶微积分近似解法。
例 5-21试用不同的 n, a组合构造出 s0.5阶滤波器，并比较其频域响应效果。
解用前面给出的程序可以立即设计出如下的滤波器，并绘制出滤波器的 Bode图，如图 5-41所示。从得出的频率响应曲线看，若选择 a =0.5则相位曲线在给定的区域内不为恒值，所以拟合效果将不理想，但高频处赋值将有下降趋势，对噪声将有抑制作用。
>> H1=iir_pade(0.5,0,3,0.01); H2=iir_pade(0.5,0.5,3,0.01);
H3=iir_pade(0.5,0,4,0.01); H4=iir_pade(0.5,0,5,0.01);
bode(H1,'-',H2,'--',H3,':',H4,'-.')
图 5-41不同 (n, a)组合下的频域频域响应拟合
用 H4,a z.1表示滤波器，可以得出设计出的各类滤波器。例如，
.1 113.1 . 113.1z.1 . 56.57z.2 + 56.57z.3
H4,0 z =
8 . 8z.2 + z.4
如果提高滤波器的阶次，如，选择 n = 10，且 a分别选择为 0和 0.5，就可以设计离散滤波器，得出的时域拟合结果如图 5-42所示（由于 H1在高频增益过大，对噪声的放大很大，故在图中没有绘制）。可见，即使选取的 a =0.5已经抑制了高频增益，总体时域响应曲线仍然不是很理想。
>> H1=iir_pade(0.5,0,10,0.01); H2=iir_pade(0.5,0.5,10,0.01);
t=0:0.01:pi; y=exp(-t).*sin(3*t+1); y0=glfdiff(y,t,0.5);
y1=lsim(H1,y,t); y2=lsim(H2,y,t); plot(t,y0,t,y2,'--')
和前面介绍的连续滤波器（尤其是 Oustaloup滤波器）相比，离散滤波器拟合的频段远比连续滤波器窄，且往往伴随高频处的超高增益，不适合用于信号的分数阶导数的求取，不建议使用。
H2(z)
理论值
图 5-42时域响应比较
5.6.3基于阶跃或冲激响应不变性的离散滤波器
表 5-1中列出了一组 MATLAB函数，可以对分数阶微积分算子甚至分数阶无理传递函数设计离散滤波器。这些滤波器是在阶跃或冲激响应不变性基础上设计的[18]，FOTF工具箱也提供了这 3个函数，可以直接调用。
表 5-1分数阶环节的离散滤波器
irid_fod() G=irid_fod(α,T ,N)基于冲激响应不变性的 sα算子的滤波器 srid_fod() G=srid_fod(α,T ,N)基于阶跃响应不变性的 sα算子的滤波器 irid_folpf() G=irid_folpf(τ,α,T ,N)基于冲激响应不变性的 (τs + 1).α环节的滤波器
例 5-22选择采样周期 T =0.01 s，并令滤波器阶次为 5，使利用前两个函数设计
0.5阶积分器的离散滤波器，并比较其效果。 解由下面的命令可以直接设计这两个滤波器，滤波器的 Bode图如图 5-43所示。
>> G1=irid_fod(-0.5,0.01,5); G2=srid_fod(-0.5,0.01,5);
s=fotf('s'); H=bode(s^(-0.5),logspace(-3,3));
bode(H), hold on, bode(G1,'--',G2,':'), hold off
由上面语句设计的离散滤波器如下：
0.09354z5 . 0.2395z4 +0.2094z3 . 0.06764z2 +0.003523z +0.0008224
G1(z)=
z5 . 3.163z4 +3.72z3 . 1.966z2 +0.4369z . 0.02738
2.38×10.6z5 +0.1128z4 . 0.367z3 +0.4387z2 . 0.2269z +0.04241
G2(z)=
z5 . 3.671z4 +5.107z3 . 3.259z2 +0.882z . 0.05885
可以将已知的 y(t)= e.t sin(3t + 1)信号馈入这两个滤波器，可以得出函数的积分曲线，如图 5-44所示。从得出的效果看，G1(z)滤波器的效果不佳，G2(z)滤波器（即基于阶跃响应不变性设计的滤波器）对积分行为有比较好的逼近，如果提高滤波器的阶次，
·160·分数阶微积分学——数值算法与实现
图 5-43滤波器的 Bode图比较
可能得到更好的效果。
>> t=0:0.01:pi; y=exp(-t).*sin(3*t+1);
y0=glfdiff9(y,t,-0.5,5);
y1=lsim(G1,y,t); y2=lsim(G2,y,t);
plot(t,y0,t,y1,'--',t,y2,':')
理论值与 G2(z)
G1(z)
图 5-44滤波器的积分效果比较
例 5-23现在考虑例 4-9的基准测试问题。试设计一个积分算子的离散滤波器，并
和 Oustaloup滤波器进行比较，评价其近似精度。 解选择采样周期 T =0.005 s，可以设计两种离散滤波器，在这些滤波器下获得的分数阶积分效果与理论值的比较在图 5-45所示。图中还给出了由 caputo9()函数得出的计算结果，其最大误差为 0.0358，得出的曲线由点线表示，误差比较大；由阶跃响应不变性离散滤波器得出的最大误差为 0.0092，图中的虚线比较接近理论值；冲激响应不变性离散滤波器的最大误差为 0.0496，图中的点画线也存在较大的误差；而 Oustaloup滤波器得出的最大误差为 0.0010，且不能在曲线上分辨出来，该曲线与理论值完全重合。可见，离散滤波器的效果比 Oustaloup滤波器差很多。
>> f=@(t)1/gamma(1.3)*t.^0.3-2/gamma(2.3)*(t-1).^1.3.*(t>1);
y=@(t)t+1-(t-1).^2.*(t>1); T=0.005; t=[0:T:2]'; y0=f(t);
N=5; G=srid_fod(-0.7,T,N); ya=y(t); y1=lsim(G,y0,t)+1;
G=irid_fod(-0.7,T,N); ya=y(t); y4=lsim(G,y0,t)+1;
G1=ousta_fod(-0.7,15,1e-4,1e3); y3=lsim(G1,y0,t)+1;
t1=0:0.02:2; y01=f(t1); ya1=y(t1); y2=caputo9(y01,t1,-0.7,6)+1;
plot(t,ya,t,y1,'--',t1,y2,':',t,y4,'-.',t,y3)
max(abs(y1-ya)), max(abs(y2-ya1(:)))
max(abs(y3-ya)), max(abs(y4-ya))
图 5-45基准测试问题的求解比较
必须指出的是，在设计离散滤波器和控制器时，如果阶次取得过高，可能得出不稳定的结果。
本章习题
（1）试用函数 contfrac()近似分数阶积分算子 1/s0.8，应用 Bode图评价近似效果。
（2）试用函数 carlson_fod()近似分数阶积分算子 1/s0.51，用得到的传递函数和前面介绍的函数 glfdiff9()计算 RLDt .0.51 sin t2，检验得到的结果是否一致。
（3）试用函数 matsuda_fod()近似阶积分算子 1/s0.6，设置不同的参数可以得
分数0
RLD.0.6e.t
到不同的传递函数。用这些传递函数计算 0 t ，应用第 2章提供的函数 ml_func()检验哪种参数下的计算结果更精确。0.7RLD.0.7√
（4）试用函数 ousta_fod()近似分数阶积分算子 1/s，并计算 0 t sin t +3，设计数值算法检验之前的结果是否正确。
（5）试用函数 ousta_fod()近似分数阶微分算子 s0.3，设置不同的参数可以得到不同RLD0.3 √
的传递函数。用设计的这些传递函数计算 0 t cosh t，应用第 4章提供的函数 glfdiff9()检验哪种参数下计算结果更精确。
0.2RLD0.2
（6）试用函数 new_fod()近似分数阶积分算子 s，并计算 0 t (1/Γ(t))。另外设
·1
62·分数阶微积分学——数值算法与实现
计一种计算 RLDt 0.2 (1/Γ(t))的方法，检验之前的计算是否正确。
...
（7）本章介绍了种连续滤波器的设计方法，包括连分式方法、 Matsuda–Fujii方法、
多0
Carlson方法、 Oustaloup方法及 Charef方法等。同样选择 α = 1/2和 α = .1/2，试通过调节设计参数的方法，找出能够在 ω ∈ [10.4 , 104]频率段内得到的最好的拟合效果。同时，利用设计的滤波器评价分数阶导数与积分计算的精度。
（8）考虑例 5-6的问题，试找出合适的频率段与阶次，使得 Oustaloup滤波器可以在 t ∈ [0, 5000] s区间很好地拟合 s.0.5分数阶算子的行为。
（9）试将分数阶传递函数输入 MATLAB工作空间[8]。
s0.3 +3 2
G(s)= (s0.2 + 3)(s0.4 + 4)(s0.4 + 3)
（10）假设典型单位负反馈控制系统的模型为
0.8s1.2 +2 1.2s0.72 +1.5s0.33
G(s)= ,Gc(s)=
3s0.8
1.1s1.8 +0.8s1.3 +1.9s0.5 +0.4
..
试将此模型输入 MATLAB工作空间[8]。
（11）将下面的分数阶传递函数模型输入 MATLAB工作空间，并利用函数 ousta_fod()和 new_fod()对此模型进行拟合。
0.5s0.9 +1
G(s)= s3 + 23s2.5 + 19s0.7 +1
（12）已知多变量系统的分数阶传递函数矩阵模型为
12 1.35s1.2 +2.3s0.9 +1 4.13s0.7 +1 G(s)= 11
.
0.52s1.5 +2.03s0.7 +1 3.8s0.8 +1
试将其输入 MATLAB工作空间。
（13）考虑下面的无理分数阶模型，应用频域响应近似方法找出它的有理近似模型。
1
G(s)= √
s2 + s +1
（14）考虑无理分数阶模型 G(s) = 1/ √ s +1，试用函数 charef_fod()找出其有理近似模型。
√
（15）考虑无理模型 G(s)= e. s/( √ s + 1)，试绘制阶跃响应曲线。已知系统阶跃响应的解析解如下，试评价不同计算步长下阶跃响应的计算精度。
√
11
y(t)= erfc √. et+1erfc t + √
2 t 2 t
（16）考虑隐式无理系统模型 G(s) = 1/(5s + 1)0.8。可以推导出其阶跃响应的解析解为 y(t)=(t/5)0.8 E01.,81.8 (.t/5)，试比较数值 Laplace反变换、连续和离散滤波器，评价参数对求解精度的影响。
（17）下面是离子交换聚合金属材料的辨识模型[19]：
340
G(s)= s0.756(s2 +3.85s + 5880)1.15
试用本章提供的各种方法找出它的有理近似模型，并比较这些有理近似模型的频域特性。
（18）考虑下面的无理分数阶模型，试用 Implicit model模块搭建 Simulink框图。设输入信号为 u(t)= cos t/Γ(t)，利用搭建的 Simulink框图计算模型的输出，并设法验证输出是否正确。 1
G(s)= (1 + 2s)0.4(1 + 4s)0.2
（19）选择采样周期 T =0.01 s，滤波器阶次为 5，试用函数 irid_fod()和 srid_fod()设计 0.7阶积分器的离散滤波器。应用得到的滤波器计算 RLD.0.7 cosh t/Γ(t)，
0 t
并应用 glfdiff9()函数检验计算结果。
（20）试选择合适的整数阶传递函数，近似下面的分数阶模型，并比较频域响应拟合的效果[8]。
25 562920(s +1.0118)0.6774
G(s)= ,G(s)=
(s2 +8.5s + 25)0.2 (s2 + 54.7160s + 590570)0.8387
（21）已知分数阶模型[8]
5 5s0.6 +2
G1(s)= ,G2(s)=
s2.3 +1.3s0.9 +1.25 s3.3 +3.1s2.6 +2.89s1.9 +2.5s1.4 +1.2
试求出能够较好拟合原始模型的整数阶模型，讨论采用何种阶次组合能得出较好的效果。试从频域响应和阶跃响应角度比较系统降阶模型。
参考文献
[1] Oustaloup A. La dérivation non entière：Théorie，synthèse et applications[M]. Paris： Hermès，1995.
[2] Petrá. I. Fractional-order nonlinear systems：Modelling，analysis and simulation[M]. Beijing：Higher Education Press，2011.
[3] ChenYQ，Vinagre B M. A new IIR-type digital fractional order differentiator[J]. Signal Processing，2003，83（11）：2359–2365.
[4] Petrá. I，Podlubny I，O’Leary P. Analogue realization of fractional order controllers [R]. TU Ko.ice：Fakulta BERG，2002.
[5] Matsuda K， Fujii H. H∞-optimized wave-absorbing control： analytical and
experimental results[J]. Journal of Guidance，Control and Dynammics，1993，16（6）：1146–1153.
[6] Oustaloup A，Levron F，Nanot F，et al. Frequency band complex non integer differentiator：characterization and synthesis[J]. IEEE Transactions on Circuits and Systems I：Fundamental Theory and Applications，2000，47（1）：25–40.
·164·分数阶微积分学——数值算法与实现
[7] XueDY，Zhao C N，Chen Y Q. A modified approximation method of fractional order system[C]. Proceedings of IEEE Conference on Mechatronics and Automation. Luoyang, China，2006，1043–1048.
[8] 薛定宇 .控制系统计算机辅助设计—— MATLAB语言与应用 [M]. 4版.北京：清华大学出版社，2022.
[9] XueDY，Chen Y Q. Suboptimum H2 pseudo-rational approximations to fractional-order linear time invariant systems. // Sabatier J，Agrawal O P，Machado J A
T. Advances in fractional calculus — Theoretical developments and applications in physics and engineering[M]. Dordrecht：Springer，2007.
[10] Sabatier J，Farges C. Analysis of fractional models physical consistency[J]. Journal of Vibration and Control，2015，23（6）：895–908.
[11] Monje C A，Chen Y Q，Vinagre B M，et al. Fractional-order systems and controls： Fundamentals and applications[M]. London：Springer，2010.
[12] Charef A，Sun H H，Tsao Y Y，et al. Fractal system as represented by singularity function[J]. IEEE Transaction on Automatic Control，1992，37（9）：1465–1470.
[13] Meng L，Xue D Y. A new approximation algorithm of fractional order system models based optimization[J]. Journal of Dynamic Systems，Measurement and Control，2012， 134（4）：504–511.
[14] MathWorks. Signal processing toolbox user’s guide [Z]，2007.
[15] Tseng C-C，Pei S-C，Hsia S-C. Computation of fractional derivatives using Fourier transform and digital FIR differentiator[J]. Signal Processing，2000，80（1）：151–159.
[16] Chen Y Q，Vinagre B M. A new IIR-type digital fractional order differentiator[J]. Signal Processing，2003，83：2359–2365.
[17] Petrá. I. Digital fractional order differentiator/inte grator — FIR type[OL]. [2022-10-2]. https://www.mathworks.com/matlabcentral/fileexchange/3672-digital-fractional-order-differentiator-integrator-iir-type.
[18] Chen Y Q. Impulse response invariant discretization of fractional order integra-tors/differe ntiators[OL]. [2022-10-2]. http://www.mathworks.cn/matlabcentral/ fileexchange/21342-impulse-response-invariant-discretization-of-fracti onal-order-integrators-differentiators.
[19] Caponetto R，Dongola G，Fortuna L，et al. Fractional order systems：Modeling and control applications[M]. Singapore：World Scientific，2010.