Scaricando questo file matlab è possibile mandare in esecusione un esercizio completo sulla sintesi del controllore digitale.
Eseguendo digdem.m in ambiente Matlab si apre un menu a finestra. A questo punto è possibile scegliere tra 4 opzioni:
Il programma fornisce anche delle spiegazioni sintetiche sul tipo di operazioni che esegue.
Per quanto riguarda il controllo della velocità dell'hard disk. Il codice esegue i seguenti passi:
Matlab code - drag and drop:
echo off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DEMO su Controllo Digitale %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% risp=menu('Quale argomento ?','C 2 D','Disegno dal Continuo','Metodo Analitico', 'PID with Z-N','Progetto Hard--Disk') if risp==1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Dal continuo al discreto C --> D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t=0:.01:3; k=1:30; dt=0:(3/30):3; y=sin(t).*exp(-t); dy=sin(dt).*exp(-dt); plot(t,y,dt,dy,'*') disp(' %%%%%%%%%%%%%%%%%%%%') disp(' Conversioni Tempo Continuo--Tempo Discreto') %%%%%%%%%%%%%%%%%%%%') disp(' ') disp(' La conversione di un SEGNALE Tempo Continuo in un segnale Tempo') disp(' Discreto su intervalli di tempo T e'' senz''altro unica e facilmente ') disp(' realizzabile con la tecnica del CAMPIONAMENTO:') disp(' nel grafico e'' mostrata una funzione del tempo, ') disp(' y(t) = sin(t)*exp(-t), 0 < t < 3 sec ') disp(' ed il corrispondente segnale campionato ogni T = .1 sec ') disp('Any key to continue....') pause disp(' ') disp(' E'' possibile trovare la relazione tra la L-trasformata di y') disp(' e la Z-trasformata di dy. Per esempio, il gradino') disp(' 1/s ---campionamento---> Sum[exp(-skT)] ---def---> z/(z-1) ') disp(' [Sum e'' la Sommatoria per k=0,inf ]') disp(' [z e'' def. = exp(sT)]') disp(' Il passaggio sequenze numeriche <---> polinomi in z e'' chiarissimo') disp(' Quindi da una equazione alle differenze si passa automaticamente ad') disp(' una FdT in z') disp('') disp(' La conversione di un SISTEMA TC in TD e'' invece impossibile in modo') disp(' esatto. SISTEMA e'' una relazione dinamica, rappresentata da una equazione') disp(' differenziale o alle differenze, tra lo spazio dei segnali di ingresso e') disp(' quello dei segnali di uscita. Se si realizza un sistema TD la cui') disp(' uscita, in risposta ad un particolare ingresso campionato,') disp(' corrisponde esattamente alla uscita (campionata) del sistema continuo,') disp(' questo non avverra'' per nessun altra coppia ingresso-uscita.') disp('') disp(' Esempio: G(s) = 10/(s + 5), risposta impulsiva: y=10*exp(-5t);') disp(' Il campionamento di questa risposta ha z-trasformata 10z/(z-exp(-5T))') yc=impulse(10,[1 5],dt); yd=dimpulse([10 0],[1 -exp(-5*0.1)],length(dt)) ; stairs(dt,yd');hold on;plot(dt,yc'); hold off; disp('Any key to continue....') pause disp(' ') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step response %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% yc=step(10,[1 5],dt); yd=dstep([10 0],[1 -exp(-5*0.1)],length(dt)) ; stairs(dt,yd');hold on;plot(dt,yc'); hold off; disp('Any key to continue....') pause disp(' ') disp(' Esempio: G(s) = 10/(s + 5), risposta gradino: y=2-2*exp(-5t);') disp(' Il campionamento di questa risposta ha z-trasformata ') disp(' 2*z/(z-1) - 2*z/(z-exp(-5T))') disp(' Devo trovare G(z): G(z)*z/(z-1) eguagli quella risposta') disp(' Risulta G(z) = 2(1-exp(-5T))/(z-exp(-5T))') %%%%%%%%%%%%%%%%%%%% % step %%%%%%%%%%%%%%%%%%%% yc=step(10,[1 5],dt); yd=dstep([2*(1-exp(-5*0.1))],[1 -exp(-5*0.1)],length(dt)) ; stairs(dt,yd');hold on;plot(dt,yc'); hold off; disp('Any key to continue....') pause disp(' ') %%%%%%%%%%%%%%%%%%%% % impulse %%%%%%%%%%%%%%%%%%%% yc=impulse(10,[1 5],dt); yd=dimpulse([2*(1-exp(-5*0.1))],[1 -exp(-5*0.1)],length(dt)) ; stairs(dt,yd');hold on;plot(dt,yc'); hold off; disp('Any key to continue....') pause disp(' ') elseif risp==2 %%%%%%%%%%%%%%%%%%%% % Dal Continuo %%%%%%%%%%%%%%%%%%%% disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%% ') disp(' Trasposizione di un controllore TC in TD') disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%% ') disp(' E'' data la f.d.t di un motore elettrico') num=[1]; den=[1 1 0]; printsys(num,den); disp(' per il quale, in base a specifiche sul margine di fase e sulla ') disp(' velocita'' di risposta, e'' stato disegnato il compensatore') nc=[4 4]; dc=[1 2]; printsys(nc,dc); disp(' Trovare l''equivalente C(z) scegliendo Ts con la regola pratica') disp(' Ts = 0.3 / Wco') disp(' ') disp(' ') disp(' Wco e'' la pulsazione di crossover per il sistema compensato in anello aperto') disp(' Cominciamo dai diagrammi di Bode per l''impianto ed il compensatore') bode(num,den); disp('Any key to continue....') pause disp(' ') bode(nc,dc); disp('Any key to continue....') pause disp(' ') disp(' Il compensatore e'' un lead con massimo anticipo in') disp(' Wm = sqrt(1*2) , dove vale Phim = arcsin(((4+2)/2 - 2)/((4+2)/2)) = 20 deg') disp(' Diagramma compensato: ') [ncg,dcg]=series(nc,dc,num,den); bode(ncg,dcg) disp('Any key to continue....') pause disp(' ') disp(' Il taglio e'' circa in Wco = 2 rad/sec => Ts = 0.15') ts = 0.15; disp(' Applichiamo la regola di Tustin sostituendo') disp(' s con a (z-1)/(z+1) ') disp(' cercando di non distorcere alla pulsazione Wm, per la quale il') disp(' compensatore e'' stato progettato') Wm=sqrt(2); a=Wm/tan(Wm*ts/2), disp(' da cui si ottiene') ncz=3.74*[1 - 0.86] dcz=[1 -0.74] disp(' Confrontiamo le risposte in frequenza dei compensatori') w = logspace(0,30); [mag,phase] = bode(nc,dc,w); [mtus,ptus,wtus] = dbode(ncz,dcz,ts); % Now plot the results as a comparison. Press any key after the plot ... subplot(211) semilogx(w,20*log10(mag),wtus,20*log10(mtus)) xlabel('Frequency (rad/sec)'), ylabel('Gain db') title('c2d comparison plot') subplot(212) semilogx(w,phase,wtus,ptus) xlabel('Frequency (rad/sec)'), ylabel('Phase deg') pause % Press any key to continue ... disp('Ora confrontiamo i risultati finali, ad esempio in termini di') disp(' risposte al gradino degli anelli chiusi') [ncgl,dcgl]=cloop(ncg,dcg); yc = step(ncgl,dcgl,(1:50)*ts); [nz,dz]=c2dm(num,den,ts); %ZOH equiv. [ncgz,dcgz]=series(ncz,dcz,nz,dz); [ncglz,dcglz]=cloop(ncgz,dcgz); yd=dstep(ncglz,dcglz,50); stairs(1:50,yd');hold on;plot(1:50,yc'); hold off; disp('Any key to continue....') pause disp(' ') elseif risp==3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Metodo analitico %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Questa opzione non `e per il momento disponibile !!! ') elseif risp==4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PID con ZN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp(' G(s) = 10/(s(s+2)) '); ts=0.25; ng=10; dg=[1 2 0]; [ng_z,dg_z]=c2dm(ng,dg,ts) %% Usa ZOH axis('square'),zgrid('new'), rlocus(ng_z, dg_z) [km,pole]=rlocfind(ng_z,dg_z) %% trova il guadagno di oscillazione wm=angle(pole(1))/ts kp=0.6*km % parametri di Zieger--Nichols kd=kp*pi/(4*wm) ki=kp*wm/pi nc_z = [kp*ts+kd+ki*ts*ts, -kp*ts-2*kd, kd]; dc_z=[ts -ts 0] [ncg,dcg] = series(nc_z,dc_z,ng_z,dg_z) [ncl,dcl] = cloop(ncg,dcg,-1) ddamp(dcl,ts) hold off axis('square'),zgrid('new'), rlocus(ncg,dcg) [km,pole]=rlocfind(ncg,dcg) %% trova il guadagno di oscillazione elseif risp==5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % diskdemo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diskdemo end
Il file diskdemo.m che è richiamato da digdemo.m è qui di seguito riportato.
Matlab code - drag and drop:
%DISKDEMO.M Demonstration design of harddisk digital controller. echo on % This file demonstrates MATLAB's ability for classical digital control % system design by going through the design of a computer HARDDISK % read/write head position controller. echo on pause % Press any key to continue ... % Using Newton's law, the simplest model for the read/write head has the % following differential equation: % % I*theta_ddot + C*theta_dot + K*theta = Ki * i % % where I is the inertia of the head assembly % C is the viscous damping coefficient of the bearings % K is the return spring constant % Ki is the motor torque constant % theta_ddot, theta_dot, and theta are the angular acceleration, % angular velocity and position of the head % i is the input current % pause % Press any key to continue ... % Taking the laplace transform, the transfer function is: % Ki % H(s) = ---------------- % I s^2 + C s + K % % Using the values I=.01 Kg m^2, C=.004 Nm/(rad/sed), K=10 Nm/rad, and % Ki=.05 Nm/rad form the transfer function description of this system I = .01; C = 0.004; K = 10; Ki = .05; NUM = [Ki]; DEN = [I C K]; printsys(NUM,DEN,'s'); pause % Press any key to continue ... % Our task is to design a digital controller that can be used to provide % accurate positioning of the read/write head. We will do the design in the % digital domain. % First we must discretize our plant since it is continuous. Since our % plant will have a digital-to-analog-converter (with a zero-order hold) % connected to its input, use the 'zoh' discretization method % of the function C2DM. Use sample time Ts = 0.005 (5 ms) Ts = 0.005; w = logspace(0,3); [mag,phase] = bode(NUM,DEN,w); [num,den] = c2dm(NUM,DEN,Ts,'zoh'); [mzoh,pzoh,wzoh] = dbode(num,den,Ts); % Now plot the results as a comparison. Press any key after the plot ... echo off subplot(211) semilogx(w,20*log10(mag),wzoh,20*log10(mzoh)) xlabel('Frequency (rad/sec)'), ylabel('Gain db') title('c2d comparison plot') subplot(212) semilogx(w,phase,wzoh,pzoh) xlabel('Frequency (rad/sec)'), ylabel('Phase deg') pause % Press any key to continue ... echo on pause % Press any key to continue ... % Now analyze the discrete system. disp('Discrete system') printsys(num,den,'z') % Plot step response subplot(111) dstep(num,den); pause % Press any key after the plot ... % The system oscillates quite a bit. This is probably due to very light % damping. We can check this by computing the open loop eigenvalues. disp('Open loop discrete eigenvalues'), ddamp(den,Ts); zgrid('new'), pzmap(1,den); pause % Press any key after the plot ... % Note that the poles are very lightly damped and near the unit circle. % We need to design a compensator that increases the damping of this system. % Let's try to design a compensator. The simplest compensator is a simple gain. rlocus(num,den); hold off; pause % Press any key after the plot ... % As shown in the root locus, the poles quickly leave the unit circle and go % unstable. We need to introduce some lead or a compensator with some zeros. % Try the compensator: K(z + a) % D(z) = -------- where a < b % (z + b) pause % Press any key to continue ... % Form compensator and connect in series with our plant % Use a = -.85 and b = 0. [numc,denc] = zp2tf([.85 ]',[0]',1); [num2,den2] = series(numc,denc,num,den); % Lets see how the frequency response has changed. [mag,phase,w] = dbode(num,den,1); % Use normalized frequency [mag2,phase2] = dbode(num2,den2,1,w); % Now plot a comparison plot. Press any key after the plot ... echo off subplot(211), semilogx(w,20*log10(mag),w,20*log10(mag2)) xlabel('Frequency (rad/sec)'), ylabel('Gain dB') subplot(212), semilogx(w,phase,w,phase2) xlabel('Frequency (rad/sec)'), ylabel('Phase deg') % Plot -180 degree line hold on; semilogx([min(w(:)),max(w(:))],[-180,-180],'w--'); hold off; pause % Press any key to continue ... echo on % So our compensator has indeed added lead. % Now let's try the root locus again with our compensator subplot(111) zgrid('new'), rlocus(num2,den2); hold off; pause % Press any key after plot ... % This time the poles stay within the unit circle for some time. % Now its your turn, Using RLOCFIND chose the poles with greatest damping % ratio. (The lines drawn by ZGRID show the damping ratios from z=0 to 1 % in steps of .1) pause % Press any key and then choose a point on the plot [k,poles] = rlocfind(num2,den2); disp(['You chose gain: ',num2str(k)]), ddamp(poles,Ts); % Let's form the closed loop system so that we can analyze the design. [numc,denc] = feedback(num2,den2,k,1); % These eigenvalues should match the ones you chose. disp('Closed loop eigenvalues'), ddamp(denc,Ts); pause % Press any key to continue ... % Closed loop time response dstep(numc,denc); pause % Press any key after the plot ... % So the response looks pretty good and settles in about 14 samples % which is 14*Ts secs. disp(['Our disc drive will have a seek time > ',num2str(14*Ts),' seconds.']) pause % Press any key to continue ... % Let's now look at the robustness of our design. The most common classical % robustness criteria is the gain and phase margin. The criteria is determined % by forming a unity feedback system, calculating the Bode response and looking % for the phase and gain crossovers. MATLAB contains a function MARGIN that % determines the phase and gain margin given the Bode response. % Form unity feedback system by connecting our design system with the gain % we chose. Leave the loop open so we can compute the open loop Bode response. [num2,den2] = series(num2,den2,k,1); % Compute Bode response and margins [mag,phase,w] = dbode(num2,den2,Ts); [Gm,Pm,Wcg,Wcp] = margin(mag,phase,w); % Plot Bode plot with margins margin(mag,phase,w); pause % Press any key after the plot ... % Gain margin db, @ frequency, Phase margin, @ frequency Margins = [20*log10(Gm),Wcg,Pm,Wcp] % Our design is robust and can tolerate a 10 db gain increase and a 40 degree % phase lag without going unstable. By continuing this design process we may % be able to find a compensator that will stabilize the open loop system and % allow us to reduce the seek time (more damping would allow us to reduce the % settling time in the step response). echo off