微机保护算法仿真示例

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

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

结构：相减滤波单元的级联;阶次:k=96/4=24作用：使用相减滤波单元消除直流分量和4，8，12，…l•4次谐波分量类型：有限冲激响应数字滤波器（FIR滤波器）分析：H2(z)=[（1–Z–7)/（1–Z–1)]2

结构：积分滤波单元的级联;阶次:k=8作用：使用积分滤波单元增加对高次谐波的衰耗类型：有限冲激响应数字滤波器（FIR滤波器；论文误为IIR滤波器）N=96;f0=50;fs=N*f0;a1=1;b1=[1zeros(1,23)-1];

%设置差分滤波单元的传递函数系数%k=N/m=24;(谐波阶次m=4)%文献给出的积分滤波单元的传递函数a2=1;b2=[ones(1,9)];

%设置积分滤波单元的传递函数系数k=8,%另外注意：设置系数b时，b0=1;f=0:1:fs/2;h1=abs(freqz(b1,a1,f,fs));H1=h1/max(h1);%由传递函数系数确定传递函数的幅频特性H11=H1.^2;%级联差分滤波器h2=abs(freqz(b2,a2,f,fs));H2=h2/max(h2);H22=H2.^2;%级联积分滤波器h=h1.*h2;H=h./max(h);%计算级联滤波器传递函数的幅频特性HH=H.^2;Figure%级联差分、积分滤波器幅频特性subplot(221);plot(f,H11);xlabel('f/Hz');ylabel('H1');subplot(223);plot(f,H22);xlabel('f/Hz');ylabel('H2');subplot(222);plot(f,HH);xlabel('f/Hz');ylabel('H');subplot(224);plot(f,20*log10(HH+10e-10));axis([0,fs/2,-60,0]);xlabel('f/Hz)');ylabel('H/dB');分析：H2(z)=[（1–Z–7)/（1–Z–1)]2

并未针对特定谐波滤波进行设置，最终结果非常不理想：实际上二次谐波有最大输出；且3、5、6、7等各次滤波也有很大输出。结论：在N=96时，论文给出的滤波器不能作为工频选频滤波器。改进：考虑积分滤波单元1消除l•3次谐波的影响:m=3阶次:k=N/m–1=31考虑积分滤波单元2消除l•2次谐波的影响:m=2阶次:k=N/m–1=47H2(z)=（1–Z–31)（1–Z–47)/（1–Z–1)2仿真m文件作以下修改：N=96;f0=50;fs=N*f0;a1=1;b1=[1zeros(1,23)-1];

%设置差分滤波单元的传递函数系数k=24a2=1;b2=[ones(1,32)];

%设置积分滤波单元的传递函数系数k=31a3=1;b3=[ones(1,48)];

%设置积分滤波单元的传递函数系数k=47重新运行m文件，结果如下：（参考c4e2.m）滤波效果比较

原始输入信号y0=sin(2*pi*f0*t)+0.5*sin(2*2*pi*f0*t)+0.5*sin(3*2*pi*f0*t)+0.3*sin(4*2*pi*f0*t)+0.3*sin(5*2*pi*f0*t);滤波输出信号1:应用论文给出的滤波器滤波输出信号2:应用改进的滤波器C4example2：

全零点滤波器设计在N值较大时，列写系数a,b较繁琐考虑自动生成滤波器系数可使用MATLAB提供的WINDOWS标准对话框输入可选参数问题：N72后难以达到预期设计效果，原因是什么？（可作课程考试题目进行讨论分析）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=[r1r2];q=r(:);%按列拼接矩阵元素qq=[0q(1:N-1)];pp=[qq(1:2*m-1)qq((2*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');C4example3：

消除非周期分量的交流采样数据修正法

含衰减分量的原始信号为

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

=

i(1)-i(N+1)=A[1-exp(-T/)](2)（考虑其中周期分量的周期性，由衰减分量确定）由（1）/（2），有

/T=I*(0)/求得=T·I*(0)/=T·[i(k)]/[N·(i(1)-i(N+1)]及

A=/[1-exp(-T/)]=[i(1)-i(N+1)]/[1-exp(-T/)]

对原始信号修正，有：i*(t)=i(t)–A·exp(-t/)对修正后的信号运用付立叶算法，则可大大降低衰减分量带来的误差消除非周期分量的交流采样数据修正法，不但可降低付立叶算法的误差，也可降低其它算法的误差缺点是运算量较大，适宜于非实时运算*自行完成相关算法仿真的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=filter(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(修正信号)','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);C4example4：

滤序算法

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('三相正序电压')subplot(223)pl