Hzblf=Filt(Numblf, Denblf,Ts);

pidNOISE.m

clear all

Kd = 2.0;

Kp = 5.0;

Ki = 0.25;

Ts = 0.00025;

close all

NumBLF = [1 1];

DenBLF = [43.433 -41.433];

HzBLF=filt(NumBLF, DenBLF,Ts);

%freres(Hzf,'HzBLF')

NumZpid = [Kp+Kd+Ki -Kp-2*Kd Kd];

DenZpid=[1 -1];

Zpid = TF(NumZpid,DenZpid,Ts);

numZ=[1 0];

denomZ=[1];

Zp = TF(numZ,denomZ,Ts);

[ba,aa] = butter(7,0.1);

HzARF=filt(ba, aa,Ts);

%freres(HzARF,'HzARF')

Tz1=HzBLF*Zpid;

Tz=1/(Zp+Tz1)

freres(Tz,'Tz')

[mag,phase] = bode(Tz,2.0*pi*20.0);

mdf = 20*log10(mag)

phase

Signal transfer function with feedback

Tz1=HzBLF*Zpid;

Tz=Tz1/(Zp+Tz1)

freres(Tz,'Tz')

[z,p,k] = zpkdata(Tz,'v')

[b,a] = zp2tf(z,p,k)

[h,t] = stepz(b,a,200,1/Ts)

figure

plot(t,h)

pid1.m

clear all

Kd = 0.0;

Kp = 2.5;

Ki = 0.25;

Ts = 0.00025;

close all

NumPS = [1.0 -2.87463 0.716485 0.325489 -0.166391];

DenPS = [1.36339e-02 -7.6617e-04 -1.091e-02];

HzPS = filt(NumPS,DenPS,Ts);

%freres(HzPS,'HzPS')

NumCF = [1.36339e-02 -7.6617e-04 -1.091e-02];

DenCF = [1.0 -2.87463 0.716485 0.325489 -0.166391];

HzCF = filt(NumCF,DenCF,Ts);

%freres(HzCF,'HzCF')

NumBLF = [0.023023967 0.023023967];

DenBLF = [1.0 -0.953952064];

HzBLF=filt(NumBLF, DenBLF,Ts);

%freres(Hzf,'HzBLF')

[ba,aa] = butter(7,0.2);

HzARF=filt(ba, aa,Ts);

freres(HzARF,'HzARF');

NumZpid = [Kp+Kd+Ki -Kp-2*Kd Kd];

DenZpid=[1 -1];

Zpid = TF(NumZpid,DenZpid,Ts);

NumZ=[1 0];

DenZ=[1];

Zp = TF(NumZ,DenZ,Ts);

Tz1=HzPS*HzARF*HzCF*Zpid*HzBLF;

Tz= 1/(Zp+Tz1)

freres(Tz,'Tz')

[mag,phase] = bode(Tz,2.0*pi*20.0);

mdf = 20*log10(mag)

phase

Noise Transfer function attenuation at 20 HzMaximum Noise Transfer function attenuation

for Kp versus Ki when Kd=0.5for Kp versus Ki when Kd=0.5

pid2.m

clear all

Ts = 0.00025;

Fs=1.0/Ts;

close all

Kd=0.5;ir=0;

for Kp = 0.0:0.1:10

ir=ir+1;

jr=0;

for Ki = 0.0:0.1:2.0

jr=jr+1;

NumPS = [1.0 -2.87463 0.716485 0.325489 -0.166391];

DenPS = [1.36339e-02 -7.6617e-04 -1.091e-02];

HzPS = filt(NumPS,DenPS,Ts);

%freres(HzPS,'HzPS')

NumCF = [1.36339e-02 -7.6617e-04 -1.091e-02];

DenCF = [1.0 -2.87463 0.716485 0.325489 -0.166391];

HzCF = filt(NumCF,DenCF,Ts);

%freres(HzCF,'HzCF')

NumBLF = [0.023023967 0.023023967];

DenBLF = [1.0 -0.953952064];

HzBLF=filt(NumBLF, DenBLF,Ts);

%freres(Hzf,'HzBLF')

[ba,aa] = butter(7,0.1);

HzARF=filt(ba, aa,Ts);

%freres(HzARF,'HzARF');

NumZpid = [Kp+Kd -Kp-2*Kd+Ki Kd];

DenZpid=[1 -1 0];

Zpid = TF(NumZpid,DenZpid,Ts);

NumZ=[1 0];

DenZ=[1];

Zp = TF(NumZ,DenZ,Ts);

%Tz1=HzPS*HzARF*HzCF*Zpid*HzBLF;

Tz1= Zpid*HzBLF;

Tz= 1/(Zp+Tz1);

%freres(Tz,'Signal TF with FB')

[db20, dbmax] = maxnoise(Tz);

dbm20(ir,jr)=db20;

dbmmax(ir,jr)=dbmax;

end

end

Kp = 0.0:0.1:10.0;

Ki = 0.0:0.1:2.0;

figure(1);

mesh(Ki,Kp,dbm20);

xlabel('Ki');

ylabel('Kp');

zlabel('Gain(dB)');

title('Noise Attenuation Kd = 0.5');

figure(2);

mesh(Ki,Kp,dbmmax);

xlabel('Ki');

ylabel('Kp');

zlabel('Gain(dB)');

title('Maximum Noise Gain Kd = 0.5');

maxnoise.m

function [db20, dbmax] = maxnoise(Tz)

fs=0.1;

f=0;

i=0;

dbmax=-1000.0;

for j=1:100

f=f+fs;

i=i+1;

if (i==10)

i=1;

fs=fs*10;

end

[mag,phase] = bode(Tz,2.0*3.1415927*f);

md(j) = 20*log10(mag);

fr(j)=f;

ph(j)=phase;

if (md(j)>=dbmax)

dbmax=md(j); fmax=f;

end

if(f==20)

db20=md(j);

% fprintf(1,'%s %f %s %f %s \n','Noise transfer function attenuation Gain ',md(j),'dB and Phase ',ph(j),'Deg at 20 Hz')

end

if (f>= 1000)

break

end

end

% fprintf(1,'%s %f %s %f %s \n','Noise transfer function Maximum Gain ',dbmax,'dB at ',fmax,'Hz')

Noise Transfer function attenuation at 20 HzMaximum Noise Transfer function attenuation

for Kp versus Kd when Ki=0.1for Kp versus Kd when Ki=0.1

Noise Transfer function attenuation at 20 HzMaximum Noise Transfer function attenuation

for Ki versus Kd when Kp=3.0for Ki versus Kd when Kp=3.0

Noise Transfer function attenuation at 20 HzMaximum Noise Transfer function attenuation

for Kp versus Ki when Kd=0.0for Kp versus Ki when Kd=0.0

Noise Transfer function attenuation at 20 HzMaximum Noise Transfer function attenuation

for Kp versus Ki when Kd=0.0for Kp versus Ki when Kd=0.0

Noise Transfer function attenuation at 20 HzMaximum Noise Transfer function attenuation

for Kp versus Ki when Kd=0.0for Kp versus Ki when Kd=0.0

pid21.m

clear all

fid = fopen('pidpro.out', 'w');

Ts = 0.00025;

Fs=1.0/Ts;

close all

NumZ=[1 0];

DenZ=[1];

Zp = TF(NumZ,DenZ,Ts);

NumPS = [1.0 -2.87463 0.716485 0.325489 -0.166391];

DenPS = [1.36339e-02 -7.6617e-04 -1.091e-02];

HzPS = filt(NumPS,DenPS,Ts);

%freres(HzPS,'HzPS')

NumCF = [1.36339e-02 -7.6617e-04 -1.091e-02];

DenCF = [1.0 -2.87463 0.716485 0.325489 -0.166391];

HzCF = filt(NumCF,DenCF,Ts);

%freres(HzCF,'HzCF')

[ba,aa] = butter(7,0.1);

HzARF=filt(ba, aa,Ts);

%freres(HzARF,'HzARF');

NumBLF = [0.023023967 0.023023967];

DenBLF = [1.0 -0.953952064];

HzBLF=filt(NumBLF, DenBLF,Ts);

%freres(Hzf,'HzBLF')

Kd=0.0;ir=0;

for Kp = 0.0:0.1:12.0

ir=ir+1;

jr=0;

for Ki = 0.0:0.05:1.0

jr=jr+1;

NumZpid = [Kp+Kd -Kp-2*Kd+Ki Kd];

DenZpid=[1 -1 0];

Zpid = TF(NumZpid,DenZpid,Ts);

%Tz1=HzPS*HzARF*HzCF*Zpid*HzBLF;

Tz1= Zpid*HzBLF;

Tz= 1/(Zp+Tz1);

%freres(Tz,'Signal TF with FB')

[db20, dbmax] = maxnoise(Tz);

dbm20(ir,jr)=db20;

dbmmax(ir,jr)=dbmax;

if ((dbmax<2.0)&(db20<-12.0))

stb(ir,jr)=1;

else

stb(ir,jr)=0;

end;

count=fprintf(fid,'%05.2f %05.2f %+6.2f %+6.2f %1d \n',Kp,Ki,db20,dbmax,stb(ir,jr));

end

end

Kp = 0.0:0.1:12.0;

Ki = 0.0:0.05:1.0;

figure(1);

mesh(Ki,Kp,dbm20);

xlabel('Ki');

ylabel('Kp');

zlabel('Gain(dB)');

title('Noise Attenuation at 20Hz, Kd = 0.0');

figure(2);

mesh(Ki,Kp,dbmmax);

xlabel('Ki');

ylabel('Kp');

zlabel('Gain(dB)');

title('Maximum Noise Gain Kd = 0.0');

figure(4);

contour(Ki,Kp,stb);

xlabel('Ki');

ylabel('Kp');

title('Acceptable range of Kp and Ki with Kd = 0.0');

fclose(fid);

pidOPEN.m

clear all;

fid = fopen('step.out', 'w');

Ts = 0.00025;

Fs=1.0/Ts;

close all;

NumZ=[1 0];

DenZ=[0 1];

Zp = TF(NumZ,DenZ,Ts);

NumBLF = [0.023023967 0.023023967];

DenBLF = [1.0 -0.953952064];

HzBLF=filt(NumBLF, DenBLF,Ts);

%freres(Hzf,'HzBLF')

Kp=2.5;

Ki=0.02;

Kd=0.0;

NumZpid = [Kp+Kd -Kp-2*Kd+Ki Kd];

DenZpid=[1 -1 0];

Zpid = TF(NumZpid,DenZpid,Ts);

Tz1= Zpid*HzBLF;

Tz=Tz1/(Zp+Tz1);

[b,a] = tfdata(Tz,'v');

[h,t] = stepz(b,a,100,Fs);

for i=1:1:100

error(i)=abs(1.0-h(i));

t(i)=t(i)/Ts;

count=fprintf(fid,'%5.0f %8.5f %8.5f \n',t(i),h(i),error(i));

error(i)=error(i)*t(i);

end;

ITAE=sum(error)

figure(1)

plot(t,h)

grid on;

axis([0 100 0 1.2]);

xlabel('Time iteration');

ylabel('y(t)');

text(80,0.30,['Kp = ' num2str(Kp,'%6.3f')])

text(80,0.25,['Ki = ' num2str(Ki,'%6.3f')])

text(80,0.20,['Kd = ' num2str(Kd,'%6.3f')])

title(['ITAE = ' num2str(ITAE,'%8.2f') ' Ts = ' num2str(Ts,'%8.5f')]);

fclose(fid);

pidLoop.m

clear all;

close all;

%fid = fopen('step.out', 'w')

Ts = 0.00025;

Fs=1.0/Ts;

NumZ=[1 0];

DenZ=[0 1];

Zp = TF(NumZ,DenZ,Ts);

NumBLF = [0.023023967 0.023023967];

DenBLF = [1.0 -0.953952064];

HzBLF=filt(NumBLF, DenBLF,Ts);

%freres(Hzf,'HzBLF')

Ki=0.1;

ir=0;

for Kp=0:0.05:5

ir=ir+1;

jr=0;

for Kd=0:0.05:2.0

jr=jr+1;

NumZpid = [Kp+Kd -Kp-2*Kd+Ki Kd];

DenZpid=[1 -1 0];

Zpid = TF(NumZpid,DenZpid,Ts);

Tz1= Zpid*HzBLF;

Tz=Tz1/(Zp+Tz1);

[b,a] = tfdata(Tz,'v');

[h,t] = stepz(b,a,100,Fs);

for i=1:1:100

error(i)=abs(1.0-h(i));

t(i)=t(i)/Ts;

error(i)=error(i)*t(i);

end;

ITAE=sum(error);

if (ITAE > 2500) ITAE=2500; end;

y(ir,jr)=ITAE;

%count=fprintf(fid,'%+8.5f %8.5f %8.5f \n',Kp,Kd,y(ir,jr));

end;

end;

Kp = 0.0:0.05:5.0;

Kd = 0.0:0.05:2.0;

figure(1);

mesh(Kd,Kp,y);

xlabel('Kd');

ylabel('Kp');

zlabel('ITAE');

title('Integral Time Absolute Error, Ki = 0.1');

%fclose(fid);