Politecnico di Bari
Facoltà di Ingegneria
Dispensa per il Corso di
Controlli Automatici I
Uso del software di calcolo “Matlab 4.2” per il
tracciamento del luogo delle radici
Ing. Mariagrazia Dotoli
Uso del software di calcolo “Matlab 4.2” per il
tracciamento del luogo delle radici
1. Tracciamento del luogo delle radici
Consideriamo il sistema in figura e tracciamone il luogo
delle radici e il luogo delle radici complementare per
mezzo della procedura Matlab ‘rlocus’.
+
-
k
s ( s + 5)( s + 7)
» num=1;
» den=[1 12 35 0];
» roots(den)
ans =
0
-7
-5
»
»
»
»
rlocus(num,den), grid
title('Luogo delle radici')
figure(2),rlocus(-num,den), grid
title('Luogo delle radici complementare')
Determiniamo ora con la funzione ‘rlocfind’ i punti del luogo sull’asse immaginario e il punto
doppio del luogo sull’asse reale per il luogo delle radici propriamente detto.
» [k1,p1]=rlocfind(num,den)
Select a point in the graphics window
selected_point =
Luogo delle radici
Luogo delle radici complementare
10
8
8
6
6
4
4
2
Imag Axis
Imag Axis
2
0
-2
0
-2
-4
-4
-6
-6
-8
-10
-8
-6
-4
-2
0
-8
2
Real Axis
-8
-6
-4
-2
Real Axis
2
0
2
0.0138 + 5.9064i
k1 =
419.6004
p1 =
-11.9978
-0.0011 + 5.9138i
-0.0011 - 5.9138i
» [k2,p2]=rlocfind(num,den)
Select a point in the graphics window
selected_point =
-1.8940 - 0.0585i
k2 =
30.0590
p2 =
-8.1638
-1.9181 + 0.0535i
-1.9181 - 0.0535i
In effetti, si ha una equazione caratteristica
s 3 + 12 s 2 + 35s + k = 0
e quindi la tabella di Routh è la seguente.
3
1
35
2
12
k
1
420-k
0
k
Dunque il luogo interseca l’asse immaginario per k=420, ossia per
12s 2 + k = 12s 2 + 420 = 0
s = ± j 35 = ± j 5.92
3
come calcolato. Il punto doppio sull’asse reale corrisponde invece a
3s 2 + 24 s + 35
d
G( s) H (s) = −
=0
ds
s 2 ( s + 5) 2 ( s + 7) 2
ossia ad un valore
s=
− 12 + 39
= −1.92
3
con
k = −( s 3 + 12 s 2 + 35s ) = 30.04
come calcolato.
Consideriamo ora quattro diversi valori del guadagno, cui corrispondono diverse situazioni in anello
chiuso. In ciascun caso determiniamo la funzione di trasferimento in anello chiuso per mezzo della
funzione Matlab ‘feedback’ o della equivalente ‘cloop’. Fissiamo prima k=20. I tre poli in anello
chiuso sono reali. Si ha:
» k1=20;
» [numc1,denc1]=feedback(k1*num,den,1,1)
numc1 =
0
0
0
20
12
35
20
denc1 =
1
» roots(denc1)
ans =
-7.8809
-3.3649
-0.7542
Fissiamo quindi k=100. I tre poli in anello chiuso sono ora due complessi dominanti ed uno reale. Si
ha:
» k2=100;
» [numc2,denc2]=feedback(k2*num,den,1,1)
numc2 =
0
0
0
100
4
denc2 =
1
12
» roots(denc2)
35
100
ans =
-9.4099
-1.2951 + 2.9916i
-1.2951 - 2.9916i
Consideriamo ora il valore del guadagno critico k=420. A parte un polo lontano reale, gli altri due
sono immaginari puri.
» k3=420;
» [numc3,denc3]=feedback(k3*num,den,1,1)
numc3 =
0
0
0
420
12
35
420
denc3 =
1
» roots(denc3)
ans =
-12.0000
0 + 5.9161i
0 - 5.9161i
Infine, con k=500, si hanno due poli complessi instabili e un terzo polo reale non dominante.
» k4=500;
» [numc4,denc4]=feedback(k4*num,den,1,1)
numc4 =
0
0
0
500
12
35
500
denc4 =
1
» roots(denc4)
ans =
-12.4226
0.2113 + 6.3407i
5
0.2113 - 6.3407i
Risposte al gradino
2
1.5
» t1=[0:0.1:10]’;
» t2=[0:0.1:10]’;
» t3=[0:0.1:10]’;
» t4=[0:0.1:10]’;
» [y1,x1]=step(numc1,denc1,t1);
» [y2,x2]=step(numc2,denc2,t2);
» [y3,x3]=step(numc3,denc3,t3);
» [y4,x4]=step(numc4,denc4,t4);
» plot(t1,y1,t2,y2,t3,y3,t4,y4)
» grid
» axis([0 10 -2 2])
» title('Risposte al gradino')
» legend('Poli reali','Poli complessi stabili','Poli immaginari',
'Poli complessi instabili')
1
0.5
0
-0.5
-1
Poli reali
Poli complessi stabili
Poli immaginari
Poli complessi instabili
-1.5
-2
0
1
2
3
4
5
6
7
8
9
10
Come si nota, nel primo caso la risposta al gradino è analoga a quella di un sistema sovrasmorzato;
nel secondo a quella di un sistema sottosmorzato (il terzo polo è completamente dominato); nel
terzo caso si ha un polo reale dominato e due poli immaginari, per cui la risposta è in pratica
puramente oscillatoria; infine, nell’ultimo caso la risposta è divergente, poiché due dei poli sono
instabili.
Consideriamo ora in particolare il caso k=100 e analizziamo il tempo di assestamento e la
sovraelongazione percentuale del sistema.
Risposta al gradino per k=100
1.4
» yf2=y2(length(y2));
» yb1=(yf2*1.05)*ones(size(y2));
» yb2=(yf2*0.95)*ones(size(y2));
» plot(t2,y2,'k-',t2,yb1,'m',t2,yb2,'m-')
» title('Risposta al gradino per
k=100')
» grid
» zoom
1.2
1
0.8
0.6
0.4
0.2
Si determina graficamente un tempo di
assestamento al 5% di circa 2.4 secondi. Ciò è
coerente con la formula ricavata approssimando la
risposta al gradino con il suo inviluppo esponenziale:
t s5% ≅
0
0
1
2
3
4
5
6
7
8
9
10
3
δω n
dove δ e ω sono rispettivamente il coefficiente di smorzamento e la pulsazione naturale della coppia
di poli complessi e coniugati, ricavabili come segue:
» damp(denc2)
Eigenvalue
Damping
6
Freq. (rad/s)
-1.30e+000 + 2.99e+000i
-1.30e+000 - 2.99e+000i
-9.41e+000
3.97e-001
3.97e-001
1.00e+000
3.26e+000
3.26e+000
9.41e+000
Avendo memorizzato i valori delle ordinate yi, si può calcolare anche la sovraelongazione
percentuale.
» Mp2=(max(y2)-y2(length(y2)))/(y2(length(y2))-y2(1))*100
Mp2 =
23.8404
Sempre per questo sistema, possiamo calcolare l’errore di posizione, ossia l’errore a regime con
ingresso a gradino. Questo è nullo, poiché il tipo della funzione di trasferimento di anello è unitario.
Infatti:
» y2(length(t2))
ans =
1.0000
Del resto, la funzione di trasferimento dell’errore vale:
E ( s) =
s ( s + 5)( s + 7)
1
=
1 + G ( s ) s ( s + 5)( s + 7) + k
quindi l’errore di posizione è
e p = lim s→0 s ⋅ E ( s ) ⋅
1
1
=
=0.
s 1 + G (0)
ossia esso è nullo indipendentemente dal valore di k (se il sistema è stabile). Infine, verifichiamo
che per k=100 il terzo polo reale è dominato dagli altri due, calcolando l’espressione in fratti
semplici della funzione di trasferimento e i rapporti di dominanza, dati dal rapporto dei moduli dei
residui e dei poli.
» [r2,p2]=residue(numc2,denc2)
r2 =
1.3369
-0.6685 - 1.8132i
-0.6685 + 1.8132i
p2 =
-9.4099
-1.2951 + 2.9916i
-1.2951 - 2.9916i
7
» d21=abs(r2(1,1))/abs(p2(1,1))
d21 =
0.1421
» d22=abs(r2(2,1))/abs(p2(2,1))
d22 =
0.5928
» d23=abs(r2(3,1))/abs(p2(3,1))
d23 =
0.5928
Come si vede, il rapporto di dominanza relativo al polo reale è circa un quinto di quello
corrispondente ai poli complessi, che risultano quindi essere dominanti sul primo, come anticipato.
2 Bibliografia
R. C. Dorf, R. H. Bishop, Modern Control Systems, settima edizione, 1995, Addison-Wesley.
K. Ogata, Modern Control Engineering, terza edizione, 1997, Prentice Hall.
H. Saadat, Computational Aids in Control Systems Using Matlab, 1993, Mc Graw-Hill.
M. Tibaldi, Note introduttive a Matlab e Control System Toolbox, 1993, Progetto Leonardo.
8