Studio in ambiente Matlab/Simulink del sistema di controllo di ...

Universita di Padova

FACOLTA DI INGEGNERIA

Corso di Laurea in Informazione

Studio in ambiente Matlab/Simulink

del sistema di controllo di sospensioni

attive automobilistiche

Tesi di Laurea in Automazione

Relatore:ALESSANDRO BEGHI

Presentata da:LUCA FABIETTI

Sessione 1Anno Accademico 2010/11

Indice

Elenco delle figure 1

1 Introduzione 3

2 Impostazione del Problema 42.1 Descrizione del sistema . . . . . . . . . . . . . . . . . . . . . . 42.2 Specifiche di controllo . . . . . . . . . . . . . . . . . . . . . . 52.3 Derivazione del modello matematico . . . . . . . . . . . . . . 5

2.3.1 Equazioni del moto . . . . . . . . . . . . . . . . . . . . 52.3.2 Funzioni di trasferimento . . . . . . . . . . . . . . . . 6

3 Risposta in catena aperta ai segnali canonici 73.1 Inizializzazione Matlab . . . . . . . . . . . . . . . . . . . . . . 73.2 Risposta al gradino . . . . . . . . . . . . . . . . . . . . . . . . 8

4 Controllo in catena chiusa 124.1 Controllore PID . . . . . . . . . . . . . . . . . . . . . . . . . . 12

4.1.1 Schema generale . . . . . . . . . . . . . . . . . . . . . 124.1.2 Implementazione Simulink . . . . . . . . . . . . . . . . 14

4.2 Luogo della Radici . . . . . . . . . . . . . . . . . . . . . . . . 194.3 Reti Correttrici . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.3.1 Errore a regime . . . . . . . . . . . . . . . . . . . . . . 264.3.2 Sovraelongazione . . . . . . . . . . . . . . . . . . . . . 274.3.3 Pulsazione di attraversamento . . . . . . . . . . . . . . 284.3.4 Determinazione della rete corretrice . . . . . . . . . . 28

5 Conclusione 36

1

Elenco delle figure

2.1 Modello fisico . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.1 Risposta al disturbo in catena aperta . . . . . . . . . . . . . . 83.2 Particolari della risposta al disturbo . . . . . . . . . . . . . . 93.3 Risposta al gradino da ingresso di controllo . . . . . . . . . . 93.4 Schema del sistema di spospensioni in catena chiusa . . . . . 103.5 Schema del sistema di spospensioni in catena chiusa . . . . . 11

4.1 Realizzazione Simulink del sistema di sospensioni automobi-listiche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

4.2 Schema generale PID . . . . . . . . . . . . . . . . . . . . . . 154.3 Implementazione Simulink PID controller . . . . . . . . . . . 154.4 Schema in catena chiusa con PID controller . . . . . . . . . . 164.5 Risposta ad un disturbo W = 0,1m . . . . . . . . . . . . . . . 174.6 Risposta ad un disturbo W = 0,1m . . . . . . . . . . . . . . . 184.7 Variazione dell’errore rispetto al riferimento r=0 . . . . . . . 184.8 Forza generata dal controllore PID in risposta al disturbo . . 194.9 Luogo radici in assenza di controllore . . . . . . . . . . . . . . 214.10 Luogo radici in presenza di controllore . . . . . . . . . . . . . 224.11 Risposta al disturbo 0.1m . . . . . . . . . . . . . . . . . . . . 234.12 Luogo con controllore definitivo . . . . . . . . . . . . . . . . . 244.13 Risposta al disturbo 0.1m . . . . . . . . . . . . . . . . . . . . 254.14 Schema a blocchi di un sistema con retroazione negativa . . . 264.15 Diagramma di Bode di G1(s) . . . . . . . . . . . . . . . . . . 274.16 Diagramma di Bode e margini di W (s) . . . . . . . . . . . . . 294.17 Diagramma di Bode di C(s) . . . . . . . . . . . . . . . . . . . 314.18 Diagramma di Bode e margini dopo l’aggiunta del controllore 324.19 Risposta al disturbo W=0.1m . . . . . . . . . . . . . . . . . . 324.20 Diagramma di Bode e margini di W(s) dopo l’aggiunta del

controllore definitivo . . . . . . . . . . . . . . . . . . . . . . . 344.21 Diagramma di Bode di C(s) . . . . . . . . . . . . . . . . . . . 344.22 Risposta al disturbo in presenza di controllore definitivo . . . 35

2

Capitolo 1

Introduzione

L’obiettivo di questa tesi e quello di usufruire delle conoscenze acquisite nelcorso di Controlli Automatici per effettuare lo studio ed il controllo di unsistema fisico semplificato. In prima istanza si procedera con la derivazioneanalitica delle relazioni che descrivono il sistema dato. Da queste si potrain seguito ricavare la risposta del sistema in catena aperta ad un ingressoa gradino, sia dal disturbo che dalla forza di controllo. Dopo aver ricavatoda questi studi alcune informazioni fondamentali sul sistema da controlla-re, verranno proposti diversi approcci per soddisfare determinate specificheprogettuali. In particolare le tecniche di controllo utilizzate risultano esseretre : controllore PID, sintesi attraverso il luogo delle radici, sintesi in fre-quenza. In ciascuno dei tre passi si utilizzeranno diversi tipi di conoscenzee metodi per giungere pero allo stesso risultato. Fondamentali nel corso ditutto il lavoro svolto sono stati gli strumenti Matlab e Simulink che hannoreso possibile, attraverso notevoli passi iterativi, il perfezionamento dei ri-sultati ottenuti. Per questo motivo si e ritenuto indispensabile dare largospazio, nell’elaborato, ai codici Matlab e agli schemi Simulink che hannopermesso di agevolare il processo realizzativo del controllore. Dopo aververificato per ogni metodo la piena aderenza con le specifiche di progetto,si cerchera di trarre alcune conclusioni riguardanti i risultati ottenuti, ledifficolta incontrate e le conoscenze acquisite.

3

Capitolo 2

Impostazione del Problema

2.1 Descrizione del sistema

Il problema consiste nel controllo di un sistema di sospensioni di un autobus.Questo problema, che in generale coinvolgerebbe il sistema nella sua totalita,puo essere ridotto allo studio di un’unica sospensione, in modo da renderlopiu semplice. D’ora in avanti si considerera quindi una sola delle quattrosospensioni. Da un punto di vista matematico il sistema puo essere pensatocome composto da due masse, due molle e due smorzatori ideali come apparein Figura(2.1)

Figura 2.1: Modello fisico

4

2.2 Specifiche di controllo 5

I dati assegnati sono:

• massa principale(m1) = 2500 kg

• massa della sospensione(m2) = 320 kg

• costante elastica del sistema di sospensione(k1) = 80000 N/m

• costante elastica di ruota e pneumatici(k2) = 500000 N/m

• costante di smorzamento del sistema di sospensione(b1) = 350 Ns/m

• costante di smorzamento di ruota e pneumatico(b2) = 15020 Ns/m

• forza di controllo(u) = generata dal controllore che andremo a definire

2.2 Specifiche di controllo

Un buon sistema di sospensioni richiede che, di fronte ad una qualsiasi im-perfezione dell’asfalto, il sitema reagisca riducendo al di sotto di una certasoglia le possibili oscillazioni, e che queste vengano dissipate in un tem-po breve. Solo cosı si puo assicurare agli utenti a bordo una permanenzaconfortevole.

Ci si rende conto immediatamente della difficolta di valutazione delladistanza X1 - W; per questo, e per il fatto che X2-W e trascurabile, sipreferisce considerare come uscita del sistema la distanza X1-X2. Da unpunto di vista numerico le richieste sono quindi che, in presenza di un distur-bo della strada che verra simulato come un segnale a gradino, la distanzaX1-X2 abbia una sovraelongazione del 5%, e che il tempo di assestamentonon superi i 5s.

2.3 Derivazione del modello matematico

2.3.1 Equazioni del moto

Osservando la Figura 2.1 e considerando la legge di Newton, si ottengono leseguenti due equazioni, che descrivono il sistema:

−M1X1 − b1(X1 − X2)− k1(X1 − X2) + U = 0

−M2X2 + b1(X1 − X2) + k1(X1 − X2) + b2(W − X2) + k2(W − X2)− U = 0

2.3 Derivazione del modello matematico 6

2.3.2 Funzioni di trasferimento

Per i nostri scopi possiamo porre tutte le condizioni iniziali a zero, scelta lo-gica che semplifica la funzione di trasferimento, e che equivale a considerarel’autobus come appena uscito da una buca nel terreno. Grazie a questa con-siderazione e, applicando la trasformata di Laplace, si ottengono le seguenti:

−M1s2X1(s)− b1s(X1(s)− X2(s))− k1(X1(s)− X2(s)) + U(s) = 0

−M2s2X2(s) + b1s(X1(s)− X2(s)) + k1(X1(s)− X2(s)) + b2s(X1(s)− X2(s)) + k2(W(s)− X2(s))− U(s) = 0

Raccogliendo i termini simili si ottiene poi :

(M1s2 + b1s + k1)X1(s)− (b1s + K1)X2(s) = U(s)

−(b1s + k1)X1(s) + (M2s2 + (b1 + b2) + (k1 + k2))x2(s) = (b2s + k2)W(s)− U(s)

[M1s2 + b1s + k1 −(b1s + k1)

−(b1s + k1) ((M2s2 + (b1 + b2)s + (k1 + k2))

].

[X1(s)X2(s)

]=

[U(s)

(b2s + k2)W(s)− U(s)

]

A =

[M1s2 + b1s + k1 −(b1s + k1)

−(b1s + k1) ((M2s2 + (b1 + b2)s + (k1 + k2))

]Per trovare la soluzione nel dominio di Laplace basta invertire la matrice A,il cui determinante risulta :

∆ = (M1s2 + b1s + k1)(M2s2 + (b1 + b2)s + (k1 + k2))− (b1s + k1)(b1s + k1)

[X1(s)X2(s)

]=

1

∆

[(M2s2 + (b1 + b2)s + (k1 + k2)) (b1s + k1)

(b1s + k1) (M1s2 + b1s + k1)

] [U(s)

(b2s + k2)W(s)− U(s)

]

[X1(s)X2(s)

]=

1

∆

[(M2s2 + b2s + k2) b1b2s2 + (b1k2 + b2k1)s + k1k2)

(−M1s2) (M1b2)s3 + (M1k2 + b1b2)s

2 + (b1k2 + b2k1)s + k1k2)

] [U(s)W(s)

]

Da questa rappresentazione si possono infine ricavare le funzioni di trasfe-rimento G1(s) e G2(s) per l’uscita X1(s)−X2(s) relative rispettivamenteagli ingressi U(s) e W(s). Sotto la condizione W(s) = 0 si ottiene la fun-zione di trasferimento dall’ingresso di controllo U(s):

1) G1(s) =Y(s)U(s) = (M1+M2)s2+b2s+K2

∆

Ponendo invece U(s) = 0 si ottiene la funzione di trasferimento rispettoall’ingresso di disturbo W(s):

2) G2(s) =Y(s)W(s) = −M1b2s3−M1K2s2

∆

Capitolo 3

Risposta in catena aperta aisegnali canonici

3.1 Inizializzazione Matlab

Per effettuare uno studio della risposta in catena aperta ai segnali canoni-ci, si inizializza una m-file Matlab con le specifiche del sistema da studiare.Prima si definiscono le variabili in gioco con i corrispettivi valori:

m1 = 2500;

m2 = 320;

k1 = 80000;

k2 = 500000;

b1 = 350;

b2 = 15020;

Attraverso il comando tf(num,den) si creano le funzioni di trasferimentoprecedentemente calcolate:

num1 = [(m1+m2) b2 k2];

den1 = [(m1*m2) (m1*(b1+b2)+ b1*m2) (m1*(k1+k2)+b1*b2+m2*k1) (b1*k2+k1*b2)

k1*k2];

G1 = tf(num1,den1);

num2 = [-(m1*b2) -(m1*k2) 0 0];

den2 = [(m1*m2) (m1*(b1+b2)+ b1*m2) (m1*(k1+k2)+b1*b2+m2*k1) (b1*k2+k1*b2)

k1*k2];

G2 = tf(0.1*(num2),den2);

7

3.2 Risposta al gradino 8

3.2 Risposta al gradino

Una volta creato il sistema in ambiente Matlab, si puo procedere analizzandoil comportamento del sistema in reazione ad ingressi a gradino; attraverso ilcomando step(0.1*G2) si ottiene la risposta al disturbo W. In particolaresi andra a simulare un disturbo di 10 cm e, come precedentemente esposto,l’obiettivo sara quello di ottenere delle oscillazioni che variano in un rangedi ±5mm, per poi dissiparsi completamente nell arco di 5 secondi.

0 5 10 15 20 25 30 35 40 45 50−0.15

−0.1

−0.05

0

0.05

0.1

Step Response

Time (sec)

Am

plitu

de

Figura 3.1: Risposta al disturbo in catena aperta

Da questo primo grafico si osserva come le oscillazioni perdurino per untempo inaccettabile (circa 60 secondi). Per quanto riguarda l’ ampiezzaconviene riscalare gli assi attraverso il comando axis([0 10 -0.1 0.1])ottenendo :


0 1 2 3 4 5 6 7 8 9 10−0.1

−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

0.06

0.08

0.1

Step Response

Time (sec)

Am

plitu

de

Figura 3.2: Particolari della risposta al disturbo

Analizzando la Figura 3.2 si nota come in seguito ad un disturbo di 0.1msi ottengano delle oscillazioni di ampiezza di circa 0.14m, ben oltre la sogliadelle specifiche. Puo essere utile anche studiare la reazione del sistema inrisposta ad un gradino da ingresso di controllo. A tal proposito, digitandostep(G1) si ottiene:

0 5 10 15 20 25 30 35 40 45 500

0.5

1

1.5

2

2.5x 10

−5Step Response

Time (sec)

Am

plitu

de

Figura 3.3: Risposta al gradino da ingresso di controllo


Anche in questo caso si osserva come il vincolo sul tempo di assestamentonon venga rispettato: le sospensioni continuano ad oscillare per circa 60 sec.Per quanto concerne l’ampiezza delle oscillazioni, invece il vincolo e rispet-tato.La soluzione risulta essere quindi l’inseririmento di un controllore in retroa-zione. Lo schema del sistema in catena chiusa e riportato in Figura 3.4.

Figura 3.4: Schema del sistema di spospensioni in catena chiusa

Considerando lo schema del sistema e necessario, innanzitutto, associare adogni blocchetto la funzione di trasferimento appropriata. Si osserva imme-diatamente che Plant = G1(s); per ricavare il valore di F(s) si deve imporreche la FdT in catena aperta tra il disturbo e l’uscita nello schema equivalentecorrisponda a G2(s). Attraverso semplici passaggi si ricava quindi:

F (s) ·G1(s) = G2(s)⇒numf

denf· num1

den1=num2

den2

Da cui:numf

denf=num2

den2· den1

num1

Osservando infine che den1 = den2 si ottiene:

numf

denf=num2

6 den2· 6 den1

num1⇒ F (s) =

num2

num1

Lo schema definitivo del sistema sara quindi quello riportato in Figura 3.4:


Figura 3.5: Schema del sistema di spospensioni in catena chiusa

Il passo successivo sara quindi quello di delineare un controllore C(s)capace di portare il sistema a soddisfare le specifiche di progetto.Per realizzare questo controllore si seguiranno diversi approcci; ognuno diessi portera a sintetizzare un compensatore con caratteristiche differenti da-gli altri, ma in grado di risolvere il problema dato con la stessa efficacia.

I metodi proposti risultano essere tre:

• Controllore PID

• Sintesi attraverso il Luogo delle radici

• Studio in frequenza e sintesi di Bode

Capitolo 4

Controllo in catena chiusa

4.1 Controllore PID

4.1.1 Schema generale

In questa sezione si cerchera di delineare la funzione di trasferimento delcontrollore attraverso un metodo di sintesi euristico e semplice da definirenella pratica. Il risultato sara un controllore PID (proporzionale, inte-grativo, derivativo).In generale un controllore PID e della forma:

KP + KDs +KI

s=

KDs2 + KPs + KI

s(4.1)

Dove:

• KPe e il parametro proporzionale all errore,

• KDe e il parametro proporzionale alla derivata dell errore,

• KIe/ e il parametro proporzionale all integrale dell errore.

I tre parametri sopra definiti realizzano diverse azioni, il cui risultato com-plessivo e quello di portare il sistema a soddisfare le specifiche progettuali.

L’azione proporzionale realizza la funzione piu semplice ed immediata, ov-vero tende a correggere l’uscita non appena l’errore, rispetto al riferimento,diventa diverso da zero. Il suo contributo risulta essere quindi:

uP = KP ∗ e (4.2)

L’azione integrativa invece va ad incidere sul comportamento a regime delsistema. Grazie ad essa il controllore ha memoria dei valori passati del se-gnale d’errore. Questa proprieta da al controllore PID la capacita di portareil processo ad inseguire perfettamente il riferimento richiesto.

uI =KI

s(4.3)

12

4.1 Controllore PID 13

L’azione derivativa, infine, aumenta le prestazioni del controllore. Il suo sco-po e quindi quello di migliorare il comportamento del sistema retroazionatodurante il transitorio. Questo contributo diventa preponderante, infatti, nelmomento in cui l’errore inizia ad aumentare. Attraverso l’azione derivativaviene generato un contributo della forma:

ud = KDde

dt(4.4)

Spesso si tende a tralasciare l’azione derivativa ponendo KD = 0. Dei trecontributi essa rappresenta quello che piu di tutti puo destabilizzare il si-stema. In presenza, ad esempio, di un rapido mutamento del riferimento,a causa dell’azione derivativa si potrebbe generare un termine correttivotroppo elevato indesiderato. Per quanto riguarda il problema specifico dellesospensioni, questo incoveniente non si presenta, poiche e ragionevole pen-sare che il riferimento si mantenga costante a zero.Grazie alle seguenti considerazioni, nel seguito si procedera con l’implemen-tazione Simulink del nostro processo. Essendo un metodo di sintesi euristica,questo tipo di approcio permettera di agire comodamente sui i tre parametria disposizione per delineare infine la FdT voluta.


4.1.2 Implementazione Simulink

Per prima cosa si riportano di seguito le equazioni che descrivono il sistema:

−M1X1 − b1(X1 − X2)− k1(X1 −X2) + U = 0 (4.5)

−M2X2+b1(X1−X2)+k1(X1−X2)+b2(W−X2)+k2(W−X2)−U = 0 (4.6)

L’implementazione relativa in ambiente Simulink e riportata in figura (4.1)

Figura 4.1: Realizzazione Simulink del sistema di sospensioni automobilisti-che

Una volta definito il processo, il passo successivo e rappresentato dalla crea-zione del controllore PID. Lo schema generale di questo tipo di controlloree riportato in figura (4.2)


Figura 4.2: Schema generale PID

Figura 4.3: Implementazione Simulink PID controller


Dopo aver creato in Simulink uno schema equivalente, si puo infine rea-lizzare la messa a punto del sistema in catena chiusa. Si osservi l’utilitadella funzione create subsystem attraverso cui si e potuto semplificarevisivamente il sistema. Figura(4.4).

Figura 4.4: Schema in catena chiusa con PID controller

Andando ad effettuare diverse simulazioni, ed osservando la risposta delsistema in catena chiusa, si possono modificare ed ottimizzare i tre parame-tri del controllore PID. In prima approssimazione si nota che per i valori:KP = 650000,KI = 400000,KD = 170000 il sistema reagisce ad un di-sturbo a gradino W = 0,1m, smorzando le oscillazioni entro 5s (Figura 4.5);la specifica riguardante il tempo di assestamento e quindi soddisfatta anchese l’ampiezza delle oscillazioni rimane troppo elevata (circa 18mm). Per ov-viare a questo problema bisogna quindi migliorare il controllo nel transitorioaumentando KD.


0 1 2 3 4 5−0.015

−0.01

−0.005

0

0.005

0.01

Time (sec)

Am

plitu

deStep response

Figura 4.5: Risposta ad un disturbo W = 0,1m

Sempre attraverso continui passi iterativi si giunge infine a delineare i treparametri che assicurano una soddisfacente risposta del sistema al disturboW. Impostando KP = 1120000,KI = 900000,KD = 440000, tutte le spe-cifiche progettuali vengono rispettate. Si riportano di seguito i grafici dellarisposta del sistema al disturbo W (Figura 4.6), della variazione dell’erro-re rispetto al riferimento (Figura 4.7) e della forza generata dal controllore(Figura 4.8).


0 1 2 3 4 5−5

−4

−3

−2

−1

0

1

2

3

4x 10

−3 Risposta al disturbo 0.1m

Time (sec)

Am

plitu

de

Figura 4.6: Risposta ad un disturbo W = 0,1m

0 1 2 3 4 5−4

−3

−2

−1

0

1

2

3

4

5x 10

−3

Time (sec)

Am

plitu

de

Errore dal riferimento

Figura 4.7: Variazione dell’errore rispetto al riferimento r=0

4.2 Luogo della Radici 19

0 1 2 3 4 5−3

−2

−1

0

1

2

3

4

5x 10

4

Time (sec)

Am

plitu

de

Forza di controllo

Figura 4.8: Forza generata dal controllore PID in risposta al disturbo

4.2 Luogo della Radici

Un altro possibile approccio e dato dal metodo del luogo delle radici. Attra-vero questo metodo grafico si puo studiare il comportamento di un sistemain catena chiusa al variare di un parametro reale K. Sebbene la sintesi at-traverso il luogo delle radici sia studiata per sistemi del primo e del secondoordine, questa puo essere applicata anche a sistemi di ordine superiore grazieall’approssimazione dei poli dominanti.

Per prima cosa si deve ricavare la posizione dei poli del sistema non

compensato. A tal proposito si utilizza il comando roots(C) con cui si

ottengono le radici del polinomio i cui coefficienti sono elementi del vettore

C. La risposta di Matlab all’istruzione R= roots(den1) risulta essere :

R =

-23.9758 +35.1869i

-23.9758 -35.1869i

-0.1098 + 5.2504i

-0.1098 - 5.2504i

Si vede subito che i poli dominanti sono le radici -0.1098+5.2504i e -

0.1098-5.2504i, che sono radici complesse coniugate.


Si consideri quindi un modello di secondo grado della forma:

s2 + 2ξωn + ω2n (4.7)

Assumendo che i poli complessi coniugati ad esso associati siano:

p1,2 = −σ ± jω σ, ω ≥ 0

Rielaborando l’equazioni date si possono ottenere le seguenti relazioni:

ω2n = σ2 + ω2 ξωn = σ

La costante ωn e detta pulsazione naturale mentre ξ e detto coefficientedi smorzamento.

Il legame tra sovraelongazione e parametri del modello (4.7) risulta essere:

S = e− ξπ√

1−ξ2 = e−π

tanθ

Mentre per quanto riguarda il tempo di assestamento il vincolo e dato da:

Ta =3

ξωn

Questi due vincoli possono essere quindi tradotti nel posizionamento dei poliin catena chiusa, in opportune zone del piano complesso. Per prima cosasi deve procedere studiando il luogo dei poli del sistema da controllare, inassenza di controllore:

rlocus (G1)

z =- log (0,05) / sqrt (pi ^ 2 + (log (0,05) ^ 2))

sgrid (z, 0)

Il comando rlocus(SYS) fornisce e traccia il luogo delle radici di SYS alvariare del parametro K reale tra 0 e inf. Per quanto riguarda il comandosgrid(Z,Wn) esso forma nel piano complesso una regione che dipende daisuoi due parametri. In particolare il primo rappresenta il legame con lasovraelongazione, mentre il secondo con la pulsazione naturale del sistema.Grazie a questa funzione Matlab, si potra scegliere il guadagno del control-lore, in modo da soddisfare entrambe le specifiche. In figura e riportato ilrisultato dei comandi sopra indicati:


−25 −20 −15 −10 −5 0 5−150

−100

−50

0

50

100

150

0.69

0.69

Luogo delle radici in assenza di compensatore e con regione di smorzamento al 5%

Real Axis

Imag

inar

y A

xis

Figura 4.9: Luogo radici in assenza di controllore

Da questo grafico si capisce immediatamente che, comunque scelto, nonesiste valore del parametro K che porti il sistema a soddisfare le specifiche;e quindi necessario aggiungere un filtro che modifichi il luogo delle radici,spostandolo nella regione delineata dal cono di piano che individua la posi-zione richiesta per i poli. Per raggiungere questo obiettivo si deve cercare diattrarre il luogo verso il piano complesso stabile, ed in particolare, all’internodella regione d’interesse. Per farlo ci si puo servire di un paio di zeri moltovicini all’asse immaginario, che andranno a bilanciare i poli dominanti inassenza di compensatore. D’altro canto, per avere un controllore realizzabi-le, si devono aggiungere un paio di poli, che verranno scelti sufficientementelontani dall’asse immaginario, al fine di rendere piu pronto il sistema. Inseguito a queste considerazioni si propone il controllore:

C(s) =s2 + 4s+ 4.25

s2 + 90s+ 1800(4.8)

Questo puo essere ottenuto attraverso i seguenti comandi:

z1=2+0.5i;

z2=2-0.5i;


p1=30;

p2=60;

numc=conv([1 z1],[1 z2]);

denc=conv([1 p1],[1 p2]);

contr=tf(numc,denc);

Decisa la posizione di poli e zeri, si e sfruttato il comando conv(A, B)che fornisce, nel caso in cui A e B siano coefficienti polinomiali, la loromoltiplicazione. Una volta definito il controllore, si verifica come e statomodificato il luogo delle radici:

rlocus(G1*contr)

axis([-30 10 -30 30])

z=-log(0.05)/sqrt(pi^2+(log(0.05)^2));

sgrid(z,0)

Il cui risultato e:

−30 −25 −20 −15 −10 −5 0 5 10−30

−20

−10

0

10

20

300.69

0.69

Luogo radici con primo controllore

Real Axis

Imag

inar

y A

xis

Figura 4.10: Luogo radici in presenza di controllore

Attraverso l’aggiunta del filtro si e spostato il luogo delle radici nella re-gione voluta; ora, digitando il comando [k,poles]= rlocfind(G1*contr),apparira, direttamente nella grafica, un cursore attraverso cui scegliere il


guadagno K del controllore. Matlab fornira dunque il valore del guadagnoscelto e la posizione dei poli del sistema:

Selected point = -1.8009+1.9565i

k = 3.6523e+007

poles = 1.0e+002*

-0.6381+3.6207i

-0.0339-0.1320i

-0.0339-0.1320i

-0.0189+0.0197i

-0.0189-0.0197i

Infine, si definisce la FdT dal disturbo all’uscita e se ne studia la rispo-sta ad un gradino 0.1m:

sys=F*feedback(G1,k*contr);

t=0:0.01:5;

step(sys,t)

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−12

−10

−8

−6

−4

−2

0

2

4x 10

−3Step Response

Time (sec)

Am

plitu

de

Figura 4.11: Risposta al disturbo 0.1m


Il tempo di assestamento e rispettato, le oscillazioni terminano entro 5s;per quanto riguarda invece la sovraelongazione, la risposta al disturbo nonpuo essere soddisfacente. Si introduce allora un filtro con FdT:

C(s) =s2 + 8s+ 25

s2 + 75s+ 1250(4.9)

Ovvero gli zeri e i poli sono rispettivamente:

z1=4+3i;

z2=4-3i;

p1=25;

p2=50;

Procedendo come prima si ricavano il luogo delle radici, il guadagno, e infinela risposta al disturbo.

−40 −35 −30 −25 −20 −15 −10 −5 0 5 10−30

−20

−10

0

10

20

300.69

0.69

Luogo con controllore definitivo

Real Axis

Imag

inar

y A

xis

Figura 4.12: Luogo con controllore definitivo

Il guadagno selezionato risulta essere K = 6.6593e+007. Grazie a que-sta scelta i poli dominanti del sistema sono compresi nella regione a 5% dismorzamento. Ricordando infine le specifiche progettuali ovvero: tempo diassestamento minore di 5s e ampiezza delle oscillazioni compresa tra +5mm


e -5mm, si osserva che le scelte effettuate soddisfano entrambe le specifiche.Il grafico e riportato in Figura 4.13:

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−4

−3

−2

−1

0

1

2

3x 10

−3Risposta al disturbo con controllore definitivo

Time (sec)

Am

plitu

de

Figura 4.13: Risposta al disturbo 0.1m

4.3 Reti Correttrici 26

4.3 Reti Correttrici

L’ultimo metodo di sintesi proposto sfrutta l’analisi in frequenza del sistema.La strategia su cui si basa questo metodo e quella di tradurre le specificheprogettuali assegnate in vincoli in frequenza, che il sistema in catena aperta,avente funzione di trasferimento W (s) = C(s)G1(s), deve soddisfare, e infi-ne, ottenere C(s) in modo che tali vincoli siano rispettati. Per prima cosa siriportano le specifiche di partenza:

• tempo di assestamento Ta = 5sec.,

• sovraelongazione massima S = 5% (S=0.05),

• errore a regime 1% (ε = 0.01).

Figura 4.14: Schema a blocchi di un sistema con retroazione negativa

4.3.1 Errore a regime

Per quanto riguarda l’errore a regime sappiamo che esso e legato al guadagnostatico di W(s) e al tipo del sistema. Nel caso in esame attraverso la specificasull’errore si trova che W(s) dovra essere di tipo 0, ovvero hW = 0, e cheKW = 1

ε = 10.01 = 100. Dalla relazioneW (s) = C(s)G1(s) e scrivendo tutte

le funzioni di trasferimento in forma di Bode, si puo passare a vincoli sulcontrollore C(s):

W (S) =KW

shWW (s), W (0) = 1

C(S) =KC

shCC(s), C(0) = 1


G1(S) =KG1

shG1G1(s), G1(0) = 1

Da queste ralazioni si vede facilmente che:

KW = KC ·KG1 , hW = hC + hG

Avendo gia ricavato vincoli su KW e hW , bastera studiare il diagramma diBode di G1(s) per ricavare hG1 e KG1 .Attraverso il comando margin(G1) si ottiene la rappresenzione di Bode delsistema (Fig. 4.15):

−160

−140

−120

−100

−80

Mag

nitu

de (

dB)

100

101

102

103

−180

−135

−90

−45

0

Pha

se (

deg)

Bode DiagramGm = Inf dB (at Inf rad/sec) , Pm = Inf

Frequency (rad/sec)

Figura 4.15: Diagramma di Bode di G1(s)

Da questo grafico si ricava immediatamente il valore del guadagno stati-co di G1(s) espresso in decibel, (KG1)dB = −97.78dB ⇒ KG1 = 1.29 ∗ 10−5

,da cui:

KC =KW

KG1

= 7744618, hC = hW − hG1 = 0

4.3.2 Sovraelongazione

Dal vincolo sulla sovraelongazione si puo ricavare il margine di fase di W(s).E’ infatti abbastanza intuitivo che un sistema con sovraelongazione moltogrande sia vicino all’instabilita e quindi ad un margine di fase piccolo. Da


queste considerazioni si puo ricavare una stima grossolana tra margine difase e sovraelongazione, ovvero:

S decresce ⇔ mφ cresce mφ = 1− 0.8S

Da questa relazione si ricava mφ ≈ 0.96rad(55◦)

4.3.3 Pulsazione di attraversamento

In generale, per sistemi abbastanza regolari, si sa che aumentando la bandapassante del sistema si ha un guadagno in termini di prontezza del sistema.Quindi, aumentando la banda, il tempo di salita Ts diminuira in modo pres-soche proporzionale. Ma poiche vale la relazione approssimata wA ≈ 5B, sinota che aumentando la banda del sistema anche la pulsazione di attraver-samento, cioe la pulsazione per cui vale |W (jwA)| = 1, aumentera.E’ possibile ricavare anche in questo caso una stima della relazione cheintercorre tra pulsazione wA e tempo di salita Ts:

wA ≈2

Ts

Nel caso in esame non e presente un valore preciso per il tempo salita; perdeterminarne uno di riferimento e conveniente andare a studiare i risultatiforniti dagli altri metodi di sintesi. Un buon riferimento puo essere sceltocome Ts = 0.35sec. da cui si ricava wA = 5.66rad/s.

4.3.4 Determinazione della rete corretrice

Da quanto studiato in precedenza, a partire dalle specifiche progettuali, si egiunti alle seguenti relazioni:

• margine di fase mϕ = 0.96rad.,

• guadagno statico e tipo del controllore KC = 7744618 e hC = 0,

• pulsazione di attraversamento wA = 5.66rad/sec

Tutti questi vincoli si riferiscono alla FdT W(s); l’obiettivo sara ora quello dideterminare la FdT del compensatore in modo che il sistema soddisfi tuttequeste richieste.Dalla decomposizione:

W (s) = C(s)G(s) =KC

shCC(s)G1(s)

Definendo:

W (s) =KC

shCG1(s)

si osseva che questa funzione di trasferimento risulta fissata univocamente, equindi possibile tracciarne il diagramma di Bode e definire alcuni parametrifondamentali:


• woA= pulsazione di attraversamento prima dell’aggiunta del controllore

cioe tale che |W (jwoA)| = 1,

• wA= pulsazione di attraversamento richiesta,

• mϕ= margine di fase richiesto,

• moϕ = π+∠W (jwA) ovvero il margine di fase presente, prima dell’ag-

giunta del controllore, nella pulsazione di attraversamento richiestawA.

Sempre attraverso il comando margin(K*G1) (K=7744618) si ottiene ildiagramma di Bode voluto con i margini di fase e d’ampiezza gia calcolati(Fig. 4.16):

Bode DiagramGm = Inf dB (at Inf rad/sec) , Pm = 15.3 deg (at 167 rad/sec)

Frequency (rad/sec)10

010

110

210

3−180

−135

−90

−45

0

System: untitled1Frequency (rad/sec): 5.66Phase (deg): −160

Pha

se (

deg)

−40

−20

0

20

40

60

80

System: untitled1Frequency (rad/sec): 5.66Magnitude (dB): 54

Mag

nitu

de (

dB)

Figura 4.16: Diagramma di Bode e margini di W (s)

Ottenendo:

• woA= 167 rad/sec,

• wA= 5.66 rad/sec,


• mϕ= 0.96 rad,

• moϕ = 60o.

Quello che si nota immediatamente, e che la pulsazione di attraversamentorichiesta e inferiore a quella presente, mentre il margine di fase e insufficiente.Valgono cioe le relazioni :

woA > wA, ∆ϕ = 40o > 0

La rete corretrice da utilizzare in questo caso sara quindi una rete a selladetta anche rete anticipatrice e ritardatrice. Si tratta di un sistema con FdTcon tre gradi di liberta:

C(s) =1 + sTa1 + saTa

· 1 + sTr1 + srTr

A questa si impone r=1/a, Ta > 0, 0 < a < 1, Tr > 0, r > 1 e inoltreTa � Tr. Per determinare i parametri di questa rete si sfrutteranno alcunerelazioni note di cui non verra fornita la dimostrazione formale. E’ impor-tante notare che ∆ϕ = φr + φa.Si utilizzera il grado di liberta per fissare φr che di solito e scelto piccoloin un range di −6o < φr < −3o. Imponendo quindi φr = −5o anche φarisultera fissata.Definendo:

q = tan(∆ϕ− φr) = 1, C =1

|W (jwA)|si puo procedere con la soluzione della seguente eq. di secondo grado in a:

(q2 − C2 + 1)a2 + 2C2q2a+ C2(q2C2 + C2 − 1) = 0

Da cui si ricava:

a = 0.0014, r = 1/a = 708.1096

Per determinare Ta si utilizzera la seguente:

T 2a =

1

wA2

C2 − a2

a2(1− C2)

Dovendo essere Ta > 0 si trova Ta = 0.1772. Infine per Tr vale:

(rwA2tanφr)T

2r − wa(1− r) + tanφr = 0

Il cui risultato e Tr = 2.0163.

La rete complessiva puo essere riscritta come :

C(s) =1 + s(0.1772)

1 + s(0.00025)· 1 + s(2.0163)

1 + s(1428.97)


−60

−50

−40

−30

−20

−10

0M

agni

tude

(dB

)

10−5

100

105

−90

−45

0

45

90

Pha

se (

deg)

Bode Diagram

Frequency (rad/sec)

Figura 4.17: Diagramma di Bode di C(s)

Si fornisce in Figura 4.17 il diagramma di Bode di C(s) da cui si nota chegrazie al vincolo r=1/a si ha C(∞) = 1.

Digitando il seguente codice sul prompt dei comando Matlab si ottiene larisposta al disturbo W=0.1m Fig.(4.18):

numc=conv([Ta 1],[Tr 1]);

denc=conv([a*Ta 1],[r*Tr 1]);

contr=tf(numc,denc);

margin(k*G1*contr)

sys=F*feedback(G1,k*contr);

step(0.1*sys)


Bode DiagramGm = Inf dB (at Inf rad/sec) , Pm = 59 deg (at 5.65 rad/sec)

Frequency (rad/sec)10

−510

010

5−180

−135

−90

−45

0

45

System: untitled1Frequency (rad/sec): 7.73Phase (deg): −119

Pha

se (

deg)

−150

−100

−50

0

50System: untitled1Frequency (rad/sec): 7.73Magnitude (dB): −16.8

Mag

nitu

de (

dB)

Figura 4.18: Diagramma di Bode e margini dopo l’aggiunta del controllore

0 2 4 6 8 10 12 14−0.12

−0.1

−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

0.06

0.08

Step Response

Time (sec)

Am

plitu

de

Figura 4.19: Risposta al disturbo W=0.1m


Nonostante le prestazioni del sistema siano migliorate, si osserva che enecessaria un’ulteriore iterazione nella sintesi. Per rendere il sistema piurobusto, e quindi immune al disturbo, puo essere utile aumentare il marginedi fase a bassa frequenza, lasciando pressoche immutata la pulsazione diattraversamento wA = 7.66rad/sec. In questo caso la soluzione migliore eaggiungere una rete anticipatrice della forma :

A(s) =1 + sT

1 + saT

Come fatto in precedenza, anche qui si sfrutteranno alcune relazioni appros-simate che permetteranno di modificare il comportamento del sistema. Perprima cosa si definiscono:

∆ϕ = 81o, C =1

|W (j(7.66))|= 1.5, q = tan(81o)

Dalle relazioni:

C2(C2q2 + C2 − 1)a2 + 2C2q2 + (q2 − C2 + 1) = 0

T 2 =1

wA2

C2 − 1

1− a2C2

si ricavano i valori dei due parametri che sostituiti nell’equazione generaleforniscono:

A(s) =1 + s(0.9254)

1 + s(0.00166527)

La funzione di trasferimento complessiva del controllore sara quindi data da:

C(s) = KCC(s)A(s) = 7744618· 1 + s(0.9254)

1 + s(0.00166527)· 1 + s(0.1772)

1 + s(0.00025)· 1 + s(2.0163)

1 + s(1428.97)

Si riportano di seguito i diagramma di Bode in catena aperta di W(s) (Fig4.20) e quello di C(s) (Fig.4.21):


−60

−40

−20

0

20

40M

agni

tude

(dB

)

10−5

100

105

−180

−90

0

90

180

Pha

se (

deg)

Bode DiagramGm = Inf dB (at Inf rad/sec) , Pm = 65.5 deg (at 2.99e+003 rad/sec)

Frequency (rad/sec)

Figura 4.20: Diagramma di Bode e margini di W(s) dopo l’aggiunta delcontrollore definitivo

−60

−40

−20

0

20

40

60

Mag

nitu

de (

dB)

10−5

100

105

−90

−45

0

45

90

135

180

Pha

se (

deg)

Bode Diagram

Frequency (rad/sec)

Figura 4.21: Diagramma di Bode di C(s)


Dalla Figura 4.20 si puo notare come in realta vi siano molteplici pul-sazioni di attraversamento ma, essendo interessati allo studio in bassa fre-quenza, questo non influisce sul comportamento del sistema. Dalla Figura4.21, invece, si vede come il controllore definitivo porti un grande guadagnonella fase a basse frequenze e , allo stesso tempo, si nota che, rispetto alsistema originale, la pulsazione di attraversamento sia stata sensibilmenteridotta.

La risposta definitiva al disturbo W=0.1m e riportata in Figura 4.22:

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−5

−4

−3

−2

−1

0

1

2

3

4

5x 10

−3Step Response

Time (sec)

Am

plitu

de

Figura 4.22: Risposta al disturbo in presenza di controllore definitivo

Tutte le specifiche sono ora rispettate e non vi e la necessita di ulterioriiterazioni nel processo di sintesi.

Capitolo 5

Conclusione

Quest’ultimo capitolo e dedicato ad una breve conclusione riguardo il lavorosvolto nella stesura di questa tesi. Come ampiamente discusso nei capitoliprecedenti lo scopo dell’elaborato e stato quello di prendere in esame unproblema reale, semplificarlo, tradurlo in temini matematici ed elaborarediverse strategie per risolverlo. Le difficolta incontrate sono state innumere-voli: cercare di capire come impostare il problema avendo specifiche calibratesul disturbo anziche, come di norma nel corso di studi effettuato, sul riferi-mento; riuscire ad individuare per ogni strategia di sintesi la soluzione piuappropriata, sia in termini di prestazioni che di costi del controllore; cercaredi sfruttare al massimo le potenzialita degli strumenti a disposizione, pri-ma nella gestione pratica del problema, (Matlab/Simulink) poi nella stesuradella tesi stessa (LATEX); infine rapportarsi con il sito di riferimento e le so-luzioni da questo proposte, che molte volte si discostavano dalla proceduradidattica standard fornita dal corso di Controlli Automatici. Nell’affrontarele problematiche emerse nel corso del lavoro ho avuto modo di mettere inpratica le conoscenze acquisite durante il mio percorso universitario e, allostesso tempo, di approfondirle ed ampliarle.

36

Bibliografia

[1] S.Zampieri ,Dispensa di controlli Automatici, Edizioni LibreriaProgetto Padova, 2011;

[2] R.C. Dorf, R.H.Bishop, Controlli Automatici, Pearson, 2010;

[3] http://www.engin.umich.edu/class/ctms/, Control Tutorials for Ma-tlab and Simulink.

37

Studio in ambiente Matlab/Simulink del sistema di controllo di ...

Documents

Transcript of Studio in ambiente Matlab/Simulink del sistema di controllo di ...