%Program 1 --- Expression of transfer function --- num=5*[1 3]; den=conv([1 6],[1 8]); printsys(num,den,'s'); %Program 2 --- Tranformation from frequency domain to space state domain ---- [aa,bb,cc,dd]=tf2ss(num,den) %Program 3 --- Expression of state space model --- a=[-3 -4 -12; 1 0 0; 0 1 0]; b=[1 0 0]'; c=[1 3 2]; d=0; printsys(a, b, c, d); %Program 4 --- Tranformation from state space domain to frequency domain --- [numm, denn]=ss2tf(a, b, c, d, 1); printsys(numm, denn, 's') %Program 5 --- Series and parallel connection of two model --- [nums,dens]=series(numm,denn,numm,denn); printsys(nums,dens,'s') [nump,denp]=parallel(numm,denn,numm,denn); printsys(nump,denp,'s') %Program 6 --- Closed loop transfer function for unit feedback system --- [numc,denc]=cloop(numm,denn); printsys(numc,denc,'s') %Program 7 --- Closed loop transfer function for nonunit feedback system --- [numf,denf]=feedback(num,den,numm,denn); printsys(numf,denf,'s') %Program 8 --- Reduction from subsystems to single transfer function --- n1=1;d1=1;n2=1;d2=[1 1];n3=1;d3=[1 2];n4=1;d4=[1 3];n5=4;d5=1;n6=8;d6=1; n7=12;d7=1;nblocks=7;blkbuild q=[2 -5 1 0 0; 3 2 -6 0 0; 4 3 -6 2 -7 5 3 0 0 0; 6 3 0 0 0; 7 4 0 0 0]; input=1; output=4;[aa,bb,cc,dd]=connect(a,b,c,d,q,input,output); [num,den]=ss2tf(aa,bb,cc,dd); printsys(num,den,'s') %Program 9 --- Computing poles and zeros --- [z, p, k]=tf2zp(num, den) [z, p, k]=ss2zp(aa, bb, cc, dd) %Program 10 --- Creating transfer function from poles and zeros --- [numcr, dencr]=zp2tf(z, p, k); printsys(numcr, dencr, 's') %Program 11 --- Step response and plot --- [num, den]=zp2tf([], [-2+5*i -2-5*i], 2); step(num,den) %Program 12 --- Analysis of step response --- [num,den]=zp2tf([],[-1+3*i -1-3*i],3); finalv=polyval(num,0)/polyval(den,0) [y,x,t]=step(num,den); [Y,k]=max(y); peak=t(k) percent=100*(Y-finalv)/finalv %Compute rise time n=1; while y(n)<0.1*finalv n=n+1;end m=1; while y(m)<0.9*finalv m=m+1;end rtime=t(m)-t(n) %Compute settling time l=length(t); while (y(l)>0.98*finalv)&(y(l)<1.02*finalv) l=l-1;end stime=t(l) pause %Program 13 --- Meshplot of second-order system --- clg t=[0:0.05:5]; number=10; y=zeros(length(t),number); n=1; while n<=number, [num,den]=zp2tf([],[-n/4+2.5*i -n/4-2.5*i],(n/4)^2+8); [y(1:length(t),n)x,tdumb]=step(num,den,t); n=n+1; end clg mesh(y') title('Mesh Plot Showing Step Response for Twelve Pole Locations') pause %Program 14 --- Plot of time history response and phase plane --- clg [num,den]=ord2(10,0.1);subplot(221) step(num,den) subplot(222) impulse(num,den) temp=fliplr(den);a=[0 1;-temp([1 1 0])]; b=[0;1]; c=[1 0]; d=0; [y,x,t]=step(a,b,c,d); subplot(223) plot(x(:,1),x(:,2)) xlabel('x1-axis'),ylabel('x2-axis'), title('x2 vs. x1 for step input'), [y1,x1,t1]=impulse(a,b,c,d); subplot(224), plot(x1(:,1),x1(:,2)) xlabel('x1-axis'),ylabel('x2-axis') title('x2 vs. x1 for impulse input') %Program 15 --- Root locus plot --- clg num=1; den=conv(conv([1 0.5], [1 1]), [0.5 1]); [r, kvector]=rlocus(num, den); v1=0.5; v2=2.5; h1=2.5; h2=0.2; axis([-h1 h2 -v1 v2]); plot(r, '-') xlabel('Real Axis') ylabel('Imag Axis') hold on plot(real(r(1, :)), imag(r(1, :)), 'x') plot([-h1 h2], [0 0], 'w:') plot([0 0],[-v1, v2], 'w:') hold off dmp=0.707; wn=1:1:4; sgrid(dmp, wn) %Program 16 --- Nyquist plot and bode plot --- clg [num,den]=zp2tf([],[-0.5 -1 -2],5); printsys(num,den,'s') [numgc,dengc]=cloop(num,den); printsys(numgc,dengc,'s') nyquist(num,den) pause clg bode(num,den) w=logspace(-1,1,100); [mag,phase,w]=bode(num,den); [gainm,phasem,wc,wp]=margin(mag,phase,w) pause; bode(numgc,dengc) [maggc,phasegc,w]=bode(numgc,dengc); [mp,k]=max(maggc) wp=w(k) n=1; while 20*log10(maggc(n))>=-3,n=n+1; end bandwidth=w(n) %Program 17 --- Nichols chart --- num=1; den=conv([1 10],conv([1 0.0001],[4 3])); axis([-360 0 40 80]) %[mag,phase,w]=nichols(num,den); [mag1,phase1,w1]=nichols(100*num,den); plot(phase,20*log10(mag),'-',phase1,20*log10(mag1),'--') title('Nichols plot') ngrid %Program 18 --- Control system design using phase compensator --- %********************************* %Input of control object clg w=logspace(-1,2,100); num=16; den=conv([1 0],conv([1 1],[0.2 1])); bode(num,d) [magg,phasee]=bode(num,den); [gain_m,phase_m,wgg,wpp]=margin(magg,phasee,w) title('Bode plot of third-order system') pause %******************************** %Specification for control object zetaspec=0.5 settling_time_spec=2; omegabw=4*sqrt(1-2*zetaspec^2+sqrt(4*zetaspec^4-4*zetaspec^2+2))/(zetaspec*settling_time_spec) pause %******************************** %Design of phase compensator using phase lead and lag clg %lag compensator design Tla=2*pi/(0.1*omegabw); beta=0.1; gamma=1/beta; [numla,denla]=zp2tf(-1/Tla,-1/(gamma*Tla),1); %lead compensator design; omegamax=omegabw; Tle=1/(omegamax*sqrt(beta)) [numle,denle]=zp2tf(-1/Tle,-gamma/Tle,1) [numc,denc]=series(numle,denle,numla,denla); bode(numc,denc) title('Bode plot of lead/lag compensator') pause %******************************** %Open loop transfer clg [numt,dent]=series(numc,denc,num,den) bode(numt,dent) title('Bode plot of compensated system') [mag,phase,w]=bode(numt,dent); pause %******************************** %Closed loop transfer function and evaluation clg [numcl,dencl]=cloop(numt,dent); bode(numcl,dencl) pause subplot(211) axis([0 10 0 2]) step(numcl,dencl) subplot(212) axis([0 10 0 100]) step(num,den) title('Step response of closed-loop system') %Program 19 --- Control system design using PID compensator %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% num=16; den=conv([1 0],conv([1 1],[0.2 1])); bode(num,den) pause step(num,den) pause %______________________________________________ r=abs(roots(den));i=find(r>0); rp=r(i);rmx=max(rp); rmn=min(rp); wst=0.1*round(rmn); wf=20*round(rmx); dw=wf/800; w=wst:dw:wf; clc KI=1 pm=30 [mag, phase] = bode(num, den, w); phase=180/pi*unwrap(phase*pi/180); if phase(1) > (-180+pm) i = find(phase < (-180 + pm)); else i = find(phase > (-180 + pm));end if length(i)==0 disp('Phase does not cross (-180 + P.M.).') return, else i2=i(1); i1 = i2 -1; if i1 ==0 wpm=w(i2); else wa = w(i1); wb = w(i2); p1 = phase(i1); p2 = phase(i2); wpm = wa + (-180+pm - p1)/(p2-p1)*(wb-wa); end end w1=wpm;[M,ph]=bode(num, den, w1); % Returns the mag. and phase of G(w)H(w1) thta=-180 + pm - ph; thtar=thta*pi/180; wmx=round(10*w1); dw1=wmx/100; wmn=w1/10; dw2=w1/10; stab=0; while sin(thtar)/(M*w1) +KI/(w1^2) < 0 | cos(thtar) < 0 & w1 < wmx w1=w1+dw1; [M,ph]=bode(num, den, w1); thta=-180 + pm - ph; thtar=thta*pi/180; end wmp=w1; while sin(thtar)/(M*w1)+KI/(w1^2) > 0 & cos(thtar) > 0 & w1>wmn stab=stab+1; w1=w1-dw2; [M,ph]=bode(num, den, w1); thta=-180 + pm - ph; thtar=thta*pi/180; end w1mn=w1+dw2; w1=wpm; [M,ph]=bode(num, den, w1); thta=-180 + pm -ph; thtar=thta*pi/180; while sin(thtar)/(M*w1)+KI/(w1^2) >0 & cos(thtar) > 0 & w1 m o=zeros(1,n-m); mk=[o,1]; num1=conv(num,mk); else, num1=num, end numgc=[KD,Kp,KI]; numopen=conv(numgc,num1); dengc=[0, 1, 0]; denopen=conv(dengc, den); denclsd=denopen+numopen; fprintf('Row vectors of polynomial coefficients of the compensated system:\n') fprintf('Open-loop num. '),disp(numopen) fprintf('Open-loop den. '),disp(denopen) fprintf('Closed-loop den'),disp(denclsd) [magp,phasep]=bode(numopen,denopen,w); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bode(numopen,denopen) pause [Gm, Pm, wpc,wgc]=margin(magp,phasep,w); if Pm>360 Pm=Pm-360; else, end fprintf('Gain Margin = %7.3g',Gm),fprintf(' GAIn crossover w = %7.3g\n',wgc) fprintf('Phase Margin = %7.3g',Pm),fprintf(' Phase crossover w = %7.3g\n',wpc) fprintf('\n') [M,ph]=bode(numopen,denclsd, w); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subplot(211) bode(num,den) subplot(212) bode(numopen,denclsd) pause subplot(211) step(num,den) subplot(212) step(numopen,denclsd) pause frqspecexample(w,M) discr2=[ 'Roots of the compensated characteristic equation: ']; disp(discr2) r=roots(denclsd); disp(r) rreal=real(r); for l=1:n-1 if rreal(l) >=0 fprintf(' Root of the RHP, system is unstable. Change KI or phase margin & repeat.\n\n') else,end end