微机保护算法仿真示例

1、l超高速方向继电器新算法的研究超高速方向继电器新算法的研究(II)实现与实现与仿真（参阅第八章：综合仿真）仿真（参阅第八章：综合仿真）l(以下引自：电力系统自动化，以下引自：电力系统自动化，99/10，p23)l数字滤波器的设计数字滤波器的设计l 基于工频电气量原理的继电保护算法，其严重基于工频电气量原理的继电保护算法，其严重的不足之处在于受谐波分量的影响。数字滤波的不足之处在于受谐波分量的影响。数字滤波器的性能好坏直接影响保护装置的性能。器的性能好坏直接影响保护装置的性能。l选择以下数字滤波器：选择以下数字滤波器：H1 (z) =（1 Z-m/4) 2H2 (z) =（1 Z7) / （1

2、Z1) 2 lm为每周期的采样点数，取为每周期的采样点数，取 m =96；lH1(z)是有限种激响应数字滤波器，为有梳是有限种激响应数字滤波器，为有梳状频谱特性的简单数字滤波器，主要用来初状频谱特性的简单数字滤波器，主要用来初步消去整次谐波及直流分量的影响；步消去整次谐波及直流分量的影响；lH2(z)为无限冲激响应数字滤波器，采用递为无限冲激响应数字滤波器，采用递归形式实现，是低通数字滤波器，能较好地归形式实现，是低通数字滤波器，能较好地消去高频分量的影响。消去高频分量的影响。(?！) 分析：分析：H1 (z) =（1 Z-m/4) 2 结构：相减滤波单元的级联结构：相减滤波单元的级联;阶次阶

3、次 : k = 96/4=24作用：使用相减滤波单元消除直流分量和作用：使用相减滤波单元消除直流分量和 4， 8，12，l4 次谐波分次谐波分量量类型：有限冲激响应数字滤波器类型：有限冲激响应数字滤波器 （FIR滤波器）滤波器） 分析：分析：H2 (z) =（1 Z7) / （1 Z1) 2 结构：积分滤波单元的级联结构：积分滤波单元的级联;阶次阶次 : k = 8作用：使用积分滤波单元增加对高次谐波作用：使用积分滤波单元增加对高次谐波 的衰耗的衰耗类型：有限冲激响应数字滤波器类型：有限冲激响应数字滤波器 （FIR滤波器；论文误为滤波器；论文误为IIR滤波器）滤波器）lN=96;f0=50;f

4、s=N*f0;la1=1;b1=1 zeros(1,23) -1; l%设置差分滤波单元的传递函数系数设置差分滤波单元的传递函数系数 %k=N/m=24; (谐波阶次谐波阶次 m=4)l% 文献给出的积分滤波单元的传递函数文献给出的积分滤波单元的传递函数la2=1;b2=ones(1,9); l%设置积分滤波单元的传递函数系数设置积分滤波单元的传递函数系数 k=8,l%另外注意：设置系数另外注意：设置系数b时，时，b0=1;lf=0:1:fs/2;lh1=abs(freqz(b1,a1,f,fs);lH1=h1/max(h1); l %由传递函数系数确定传递函数的幅频特性由传递函数系数确定传递

5、函数的幅频特性lH11=H1.2;l%级联级联差分差分滤波器滤波器lh2=abs(freqz(b2,a2,f,fs);lH2=h2/max(h2);lH22=H2.2;l%级联级联积分积分滤波器滤波器lh=h1.*h2;H=h./max(h); l %计算级联滤波器传递函数的幅频特性计算级联滤波器传递函数的幅频特性lHH=H.2;lFigurel%级联级联差分、积分差分、积分滤波器幅频特性滤波器幅频特性lsubplot(221);plot(f,H11);lxlabel(f/Hz);ylabel(H1);lsubplot(223);plot(f,H22);lxlabel(f/Hz);ylabel

6、(H2);lsubplot(222);plot(f,HH);lxlabel(f/Hz);ylabel(H);lsubplot(224);plot(f,20*log10(HH+10e-10);laxis(0,fs/2,-60,0);xlabel(f/Hz);ylabel(H/dB);分析：分析：lH2 (z) =（1 Z7) / （1 Z1) 2 并未针对特定谐波滤波进行设置，最终结并未针对特定谐波滤波进行设置，最终结果非常不理想：实际上二次谐波有最大输果非常不理想：实际上二次谐波有最大输出；且出；且3、5、6、7等各次滤波也有很大输等各次滤波也有很大输出。出。结论：结论：在在N=96时，论文给

7、出的时，论文给出的滤波器不能作为工频滤波器不能作为工频选频滤波器。选频滤波器。 改进：改进：考虑积分滤波单元考虑积分滤波单元1消除消除l3次谐波的影响次谐波的影响:m=3阶次阶次 : k=N/m1=31考虑积分滤波单元考虑积分滤波单元2消除消除l2 次谐波的影响次谐波的影响:m=2阶次阶次 : k=N/m 1=47H2 (z) =（1 Z31) （1 Z47) / （1 Z1) 2l仿真仿真m文件作以下修改：文件作以下修改：N=96;f0=50;fs=N*f0;a1=1;b1=1 zeros(1,23) -1; %设置差分滤波单元的传递函数系数设置差分滤波单元的传递函数系数 k=24 a2=1

8、;b2=ones(1,32); %设置积分滤波单元的传递函数系数设置积分滤波单元的传递函数系数 k=31a3=1;b3=ones(1,48); %设置积分滤波单元的传递函数系数设置积分滤波单元的传递函数系数 k=47l重新运行重新运行m文件，结果如下：文件，结果如下：（参考（参考c4e2.m）l在在N值较大时，列写系数值较大时，列写系数a,b较繁琐较繁琐l考虑自动生成滤波器系数考虑自动生成滤波器系数l可使用可使用MATLAB提供的提供的WINDOWS标准对话标准对话框输入可选参数框输入可选参数l问题：问题：N 72后难以达到预期设计效果，后难以达到预期设计效果， 原因是什么？原因是什么？ （可

9、作课程考试题目进行讨论分析）（可作课程考试题目进行讨论分析）clearprompt = 采样点数采样点数(N),滤波次数滤波次数(m);title = 全零点滤波器设计输入参数全零点滤波器设计输入参数;lines= 1;def = 12,1;m = inputdlg(prompt,title,lines,def);N1=m(1,:);m1=m(2,:);N=str2num(char(N1);m=str2num(char(m1);r1=1:N/2;r2=-r1;r=r1 r2; q=r(:); %按列拼接矩阵元素按列拼接矩阵元素qq=0 q(1:N-1);pp=qq(1:2*m-1) qq(2*

10、m+2):N); %去除去除m次谐波对应元素次谐波对应元素z=exp(j*2*pi/N)*pp);p= ; b,a=zp2tf(z,p,1); f0=50; fs=N* f0; f=0:1: fs /2;h=freqz(b,a,f, fs);amp=abs(h)./max(abs(h); pha=angle(h);figuresubplot(121); zplane(b,a);subplot(222)plot(f,amp);xlabel(f/Hz);ylabel(H);subplot(224)plot(f,pha);xlabel(f/Hz);ylabel(pha);l含衰减分量的原始信号为含衰

11、减分量的原始信号为 i(t)=Aexp(-t/ ) + cos( t)l要减少衰减分量对算法的影响，应对原始信要减少衰减分量对算法的影响，应对原始信号修正号修正l修正方法：利用一个周波的数据估算衰减分修正方法：利用一个周波的数据估算衰减分量的幅值量的幅值A及时间常数及时间常数l对原始信号在对原始信号在0，T内积分，有内积分，有 I*(0)= i(t)dt/T =A 1- exp(-t/ ) / T = i(k)/N (1) l相隔一个周期的原始信号之差为相隔一个周期的原始信号之差为 = i(1) - i(N+1) = A1- exp(-T/ ) (2) （考虑其中周期分量的周期性，（考虑其中周

12、期分量的周期性， 由由衰减衰减分量确定）分量确定）l由（由（1）/（2），有），有 /T= I*(0)/ l求得求得 = T I*(0)/ = T i(k) / N (i(1) - i(N+1) 及及 A= / 1- exp(-T/ ) = i(1) - i(N+1)/ 1- exp(-T/ ) l对原始信号修正，有：对原始信号修正，有： i*(t)= i(t) Aexp(-t/ )l对修正后的信号运用付立叶算法，则可对修正后的信号运用付立叶算法，则可大大降低衰减分量带来的误差大大降低衰减分量带来的误差l消除非周期分量的交流采样数据修正法，消除非周期分量的交流采样数据修正法，不但可降低付立叶算

13、法的误差，也可降不但可降低付立叶算法的误差，也可降低其它算法的误差低其它算法的误差l缺点是运算量较大，适宜于非实时运算缺点是运算量较大，适宜于非实时运算l*自行完成相关算法仿真的自行完成相关算法仿真的m文件，与其文件，与其 它消除非周期分量的算法进行比较它消除非周期分量的算法进行比较clearT=0.02;t1=0.02; % t1：衰减分量时间常数：衰减分量时间常数N=24;Ts=T/N;t=0:Ts:4*T;y=(exp(-t/t1)-cos(2*pi*50*t); %原始信号原始信号a=1;i=1:N; bs=sin(2*pi*i/N); bc=cos(2*pi*i/N);ys=filt

14、er(bs,a,y); yc=filter(bc,a,y);ym=2*abs(complex(ys,yc)./N;%付立叶算法计算原始信号幅值付立叶算法计算原始信号幅值y1=y(1,1:N);%截取一周波原始信号原始信号I=sum(y1)/N;tt1 = T* I/(y(1)-y(N+1);A=(y(1)-y(N+1)/(1-exp(-T/tt1);y0=A*exp(-t./tt1); %计算衰减分量计算衰减分量y1=y-y0; %对原始信号修正对原始信号修正figureplot(t,y,t,y1,t,y2,t,exp(-t/t1);legend(y(原始信号原始信号),y0,y-y1(修正信

15、号修正信号),exp(-t/t1)ys2=filter(bs,a,y2); yc2=filter(bc,a,y2);ym2=2*abs(complex(ys2,yc2)./N;%付立叶算法计算修正信号幅值付立叶算法计算修正信号幅值figureplot(t,ym,t,ym2);f=50;N=24;t=0:0.02/N:0.03;va=sin(2*pi*f*t); x,y=size(va); % y:序列长度序列长度vb=sin(2*pi*f*t-2*pi/3);vc=sin(2*pi*f*t+2*pi/3); 输入正序电压输入正序电压va1=0,0,0,0, va(1:y-4); 延迟延迟NTs/6vb1=0,0,0,0, vb(1:y-4); vc1=0,0,0,0, vc(1:y-4); v1=(vb-va)-(vc1-va1)/3;v2=(vc-va)-(vb1-va1)/3;figuresubplot(221)plot(t,va,t,vb,t,vc)legend(va,vb,vc)title(