Sistema di controllo per il pendolo inverso -...

Sistema di controllo per il pendolo

inverso

Luca Rancati 644292

Servosistemi Aerospaziali

prof. Franco Bernelli Zazzera

Indice

1 Descrizione del sistema fisico 31.1 Equazioni di moto del sistema . . . . . . . . . . . . . . . . . . . . 31.2 Equilibri del sistema e linearizzazione delle equazioni differenziali 51.3 Modello agli stati . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Modello dell’attuatore . . . . . . . . . . . . . . . . . . . . . . . . 61.5 I sensori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Simulazione numerica 92.1 Il sistema non controllato . . . . . . . . . . . . . . . . . . . . . . 92.2 Il regolatore PID . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2.1 PID per il controllo dell’angolo θ . . . . . . . . . . . . . . 102.2.2 PID per il controllo della posizione x del carrello . . . . . 11

2.3 Modello numerico in Simulink e risultati . . . . . . . . . . . . . . 14

3 Il modello in laboratorio 18

4 Appendice 204.1 determinazione dell’attrito . . . . . . . . . . . . . . . . . . . . . . 204.2 Metodo di Newton . . . . . . . . . . . . . . . . . . . . . . . . . . 204.3 Funzioni matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

1

Elenco delle figure

1.1 sistema fisico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1 angolo (gradi) per 100 s di simulazione, θ(0) = 1.7 gradi . . . . . 152.2 transitorio iniziale per l’angolo teta . . . . . . . . . . . . . . . . . 152.3 posizione del carrello (cm) per 100 s di simulazione . . . . . . . . 162.4 posizione del carrello (cm) per 300 s di simulazione . . . . . . . . 162.5 angolo (gradi) per 300 s di simulazione, θ(0) = 1.7 gradi . . . . . 17

3.1 angolo (gradi), ottenuto in laboratorio . . . . . . . . . . . . . . . 183.2 posizione del carrello (cm) ottenuta in laboratorio . . . . . . . . . 19

2

Capitolo 1

Descrizione del sistemafisico

Si vuole descrivere il sistema di controllo per il pendolo inverso; per prima cosaviene sviluppato il modello numerico, poi quello di laboratorio una volta cheil modello numerico soddisfi le specifiche richieste. In questo caso le specifichesono:

• massa aggiunta al carrello di 37 g

• controllo PID

• massimo errore di 2 gradi per l’angolo dell’asta

Il sistema fisico consiste in un carrello che scorre su una guida orizzontale sulquale e incernierata un’asta; al carrello e aggiunta una massa di 0.037 kg; obiet-tivo del progetto e di controllare che l’asta resti verticale con un massimo erroreconsentito di due gradi.

1.1 Equazioni di moto del sistema

Per ricavare le equazioni di moto del sistema si individuano le coordinate libere,si ricava l’energia cinetica del sistema e l’energia potenziale, infine si utilizza laformula di Lagrange. Le coordinete libere sono la posizione x del carrello e l’an-golo θ dell’asta; l’angolo e considerato nullo sulla verticale. Ricaviamo l’energiacinetica del sistema; essa avra un contributo traslazionale ed uno rotazionale ele generiche espressioni per entrambi i contributi sono1:

Ttrasl =12m v2

Trot =12I ω2

L’energia potenziale e in generale2:

U = −m g y

1m e la massa del sistema e v la velocita, I e il momento di inerzia del sistema rispettoall’asse di rotazione e ω l’accelerazione angolare

2m e la massa e y la sua distanza normale al riferimento orizzontale

3

Servosistemi Aerospaziali Pendolo inverso

Figura 1.1: sistema fisico

Nel nostro caso l’energia cinetica e potenziale sono:

T =12(mc + mp) x2 +

12(mp/3) l2 θ2 +

12

mp l x θ cos(θ) (1.1)

U = −mp

2g l cos(θ) (1.2)

Ricaviamo ora la Lagrangiana del sistema:

L = T + U

e sostituendo otteniamo:

L =12(mc + mp) x2 +

12(mp/3) l2 θ2 +

12

mp l x θ cos(θ)− mp

2g l cos(θ) (1.3)

Sfruttiamo l’equazione di Lagrange per generare le due equazioni differenzialiche descrivono il moto del sistema.

d

d t(∂ L

∂ qk)− ∂ L

∂ qk= fk

Inserendo L trovata in precedenza e ricordando che qk sono le derivate temporalidelle coordinate libere e qk le coordinate libere stesse otteniamo:

(mc + mp) x + 12 mp l [cos(θ) θ − sin(θ) θ2] = fx

mp

3 l2 θ + mp

2 l cos(θ) x− mp

2 g l sin(θ) = fθ

(1.4)

Va evidenziato che la forzante fθ e nulla visto che non c’e la possibilita di con-trollare direttamente l’angolo θ. L’unica possibilita e di introdurre una forzadi controllo direttamente sul carrello, forza in direzione x. Inoltre e stato ag-giunto un termine di attrito proporzionale alla velocita del carrello per simularel’attrito tra carrello e guida;3 inoltre si considera:

mtot = mc + mp dove ad mc si e aggiunta la massa di 0.037 kg

3successivamente verra spiegato come e stato ricavato il coefficiente di attrito c

Luca Rancati644292

4


Il sistema di equazioni diventa:

mtot x + 12 mp l [cos(θ) θ − sin(θ) θ2] + c x = fx

mp

3 l2 θ + mp

2 l cos(θ) x− mp

2 g l sin(θ) = 0(1.5)

1.2 Equilibri del sistema e linearizzazione delleequazioni differenziali

Il sistema trovato e un generico sistema non lineare a tempo continuo del tipo:˙z(t) = f(z(t), u(t))

in cui z(t) ∈ Rn e il vettore di stato (nel nostro caso sara scelto z(t) ={θ(t), x(t), ˙θ(t), ˙x(t)}T ), u(t) ∈ Rm e il vettore degli ingressi e t ∈ R e il tempo;si definisce stato di equilibrio per u(t) = u ∀ t costante uno stato z(t) = z ∀ t

costante da cui discende che ˙z(t) = 0. Nel nostro caso troviamo:

˙x(t) = 0˙θ(t) = 0¨x(t) = 0¨θ(t) = 0

(1.6)

da cui si ricava che: {fx = 0 ∀xmp

2 g l sin(θ) = 0Dalla seconda equazione si evince che vi sono due posizioni di equilibrio perl’angolo θ: {

θ = 0θ = π

(1.7)

L’unica posizione di equilibrio fisicamente plausibile e quella nulla ed e una po-sizione di equilibrio instabile. E possibile a questo punto linearizzare le equazioniattorno alla posizione di equilibrio instabile ottenendo:

mtot x + 12 mp l θ + c x = fx

mp

3 l2 θ + mp

2 l x− mp

2 g l θ = 0(1.8)

1.3 Modello agli stati

Il sistema lineare a tempo continuo trovato in precedenza e un sistema di dueequazioni differenziali ordinarie del secondo ordine; e possibile trasformarlo in unsistema di quattro equazioni differenziali ordinarie del primo ordine aggiungendodue equazioni in questo modo:

θ − θ = 0

x− x = 0

mtot x + 12 mp l θ + c x = fx

mp

3 l2 θ + mp

2 l x− mp

2 g l θ = 0

(1.9)

Luca Rancati644292

5


Il sistema in forma matriciale diventa:

1 0 0 00 1 0 00 0 mp

3 l2mp

2 l0 0 mp

2 l mt

θx

θx

+

0 0 −1 00 0 0 −1

−mp

2 g l 0 0 00 0 0 c

θx

θx

+

000−1

{fx} =

0000

Inoltre un generico sistema dinamico ha y ∈ Rp vettore delle uscite definitedalla relazione di trasformazione delle uscite g, cosı che il sistema completodiventa: { ˙z(t) = f(z(t), u(t))

y = g(z(t), u(t))

Nel caso di sistemi lineari otteniamo:{ ˙z(t) = Az(t) + B u(t)

y = C z(t) + D u(t)

Nel nostro caso le uscite sono due, date dalle due misure effettuabili per laposizione x e l’angolo θ. Il sistema matriciale scritto in precedenza e dellaforma:

[V ] { ˙z(t)}+ [A] {z(t)}+ [B] {u(t)} = {0}Si possono quindi ricavare le matrici del sistema ponendo:

A = −V −1 A

B = −V −1 B

Le matrici presenti nella trasformazione dell’uscita sono:

C =[1 0 0 00 1 0 0

]e D =

[00

]

A questo sistema viene ora aggiunto il modello dell’attuatore, costituito da unmotore elettrico.

1.4 Modello dell’attuatore

L’attuatore e un motore elettrico che ricevendo una tensione in ingresso fornisceuna coppia in uscita. Indicando con T e V rispettivamente la coppia e la tensioneapplicate, e con ω e i la velocita angolare del motore e la corrente circolante,sfruttando la legge di Kirchhoff alle maglie e un bilancio di coppie, possiamoscrivere:

i = −Rm

Lmi− hm

Lmω + 1

LmV

ω = km

Jmi− Bm

Jmω − 1

JmT

(1.10)

dove Rm e Lm sono la resistenza e l’induttanza interne al motore, hm e lacostante della forza elettromotrice inversa, Km e la costante di coppia, Jm ilmomento di inerzia del rotore, Bm il coefficiente di attrito viscoso. Il motoree anche dotato di un riduttore con rapporto τm il quale e collegato medianteuna ruota dentata di raggio rm ad una cremagliera, quindi la coppia in uscita

Luca Rancati644292

6


e convertita in forza, e la velocita angolare in velocit‘a lineare x del carrello.Traduciamo quindi il tutto con le seguenti relazioni:

T = fx rm

τm

ω = τm xrm

=⇒ ω = τm xrm

Vi e inoltre il rendimento del motore ηm. Il sistema diventa quindi:

i = −Rm

Lmi− hm τm

Lm rmx + ηm

LmV

x = km rm

Jm τmi− Bm

Jmx− r2

m

Jm τ2m

fx

(1.11)

Questo sistema va inserito nel modello agli stati per completare la descrizionedel sistema fisico; tuttavia possiamo trascurare il transitorio elettrico rispettoa quello meccanico, cioe considerare la corrente a regime durante il transitoriomeccanico. 4 Il tutto si traduce nel considerare:

i = 0

0 = −Rm

Lmi− hm τm

Lm rmx + ηm

LmV

x = km rm

Jm τmi− Bm

Jmx− r2

m

Jm τ2m

fx

Ricavando i dalla seconda equazione e sostituendo nella terza si ricava la forzasul carrello in funzione di x e V :

fx = α x + β V

dove α e β sono ⇒ α = −τ2m

r2m

(km hm

Rm+ Bm)

β =ηm km τm

Rm rm

Il modello lineare completo comprendente anche l’attuatore diventa:

1 0 0 00 1 0 00 0 mp

3 l2mp

2 l0 0 mp

2 l mt

θx

θx

+

0 0 −1 00 0 0 −1

−mp

2 g l 0 0 00 0 0 (c + α)

θx

θx

+

000−β

{V } =

0000

in cui α e β sono stati definiti prima. Analogamente a prima si possono ricavarele matrici A,B,C,D del sistema. Nella simulazione numerica verra utilizzatoquesto modello agli stati. Considerando che c’e una tensione minima al di sottodella quale non si ha forza prodotta ed una tensione massima al di sopra dellaquale si ha la saturazione dell’attuatore, nel modello si inseriranno un bloccodead zone5 ed un blocco saturation port.6

4il tempo caratteristico del transitorio elettrico e molto molto piccolo rispetto a quellomeccanico

5con zona morta tra ± 0, 4 [V ]6con tensione limitata tra ± 6 [V ]

Luca Rancati644292

7


1.5 I sensori

Come sensori sono presenti due potenziometri; uno misura la posizione dell’asta,l’altro quella del carrello. Per simulare correttamente il loro funzionamento si eseguito il seguente procedimento:

1. analisi del segnale misurato in laboratorio

2. esecuzione della densita spettrale di potenza

3. verifica dell’ergodicita del rumore presente

4. calcolo della sensibilita del potenziometro

Eseguendo la densita spettrale di potenza si nota che il rumore ha intensitapraticamente costante a tutte le frequenze e puo quindi essere modellato comeun rumore bianco;7 A questo punto e necessario filtrare il segnale; La funzionedi trasferimento8 del filtro e la seguente:

ξfiltrata = G(s) ξrumore dove G(s) =ω

s + ω

Si e scelta una frequenza di taglio pari a 10 [Hz]. Per il calcolo della sensibilitasono stati misurati quanti volt si registravano in uscita rispetto ai centimetripercorsi dal carrello e dall’angolo dell’asta9. Si sono ricavati i seguenti valori:

• asta 4 [V/rad]

• carrello 11.11 [V/m]

7si utilizzera il blocco simulink Band limited white noise8ξ e una generica variabile e ω e la pulsazione corrispondente alla frequenza di taglio:

ω = 2 π f9misure fatte con riga e goniometro per il carrello e l’asta; i volt sono stati registrati con

matlab in laboratorio

Luca Rancati644292

8

Capitolo 2

Simulazione numerica

Vedremo ora l’analisi in anello aperto del sistema non controllato e le sue propri-eta di stabilita. Vedremo poi come formulare una legge di controllo utilizzandoun regolatore PID lavorando sulla matrice delle funzioni di trasferimento delsistema.

2.1 Il sistema non controllato

Dopo aver ricavato il modello agli stati, descritto nel capitolo precedente, sihanno le matrici A,B,C,D del sistema. Utilizzando matlab la matrice A ed isuoi autovalori sono:

A =

0 0 1 00 0 0 1

30.6043 0 0 49.4236−2.6357 0 0 −20.0989

⇒ λ =

05.0421−4.7055−20.4355

Essendoci un autovalore a parte reale maggiore di zero si conclude che ilpunto di equilibrio θ = 0 e instabile. L’autovalore nullo rappresenta la possibilitadi avere una posizione di equilibrio per ogni x. E possibile calcolare la matricedelle funzioni di trasferimento date la matrici del sistema nel seguente modo:

N(s) = C (s I −A)−1 B + D

e si ottengono due funzioni di trasferimento, una che lega l’angolo dell’asta allaforza sul carrello ed una che lega la posizione del carrello alla forza su di esso;esse sono costituite da un numeratore e dallo stesso denominatore il quale e ilpolinomio caratteristico della matrice A:

θ(s)f(s)

=−6.489 s2

s4 + 20.1 s3 − 30.6 s2 − 484.8 s(2.1)

x(s)f(s)

=2.639 s2 − 63.66

s4 + 20.1 s3 − 30.6 s2 − 484.8 s(2.2)

ottenibili con i comandi matlab:

sist = ss(A,B,C, D) crea il sistema dinamico date le matrici

9


tf(sist) funzioni di trasferimento dato il sistema dinamico

Simbolicamente le due funzioni di trasferimento sono:

θ(s)f(s)

=n1 s2

s4 + d1 s3 + d2 s2 + d3 s(2.3)

x(s)f(s)

=z1 s2 + z2

s4 + d1 s3 + d2 s2 + d3 s(2.4)

E ora possibile lavorare separatamente su ciascuna funzione di trasferimento percontrollare il sistema, elaborando un regolatore PID per ogni funzione.

2.2 Il regolatore PID

Questo tipo di regolatore prevede un controllo proporzionale all’errore rispettoad un riferimento, derivativo ed integrale. L’effeto del controllo proporzionele equello di contrastare l’errore, quello derivativo e quello di anticipare la rispostadel sistema, quello integrale di annullare l’errore a regime. Questo si traduce nelcalcolare i coefficienti del regolatore in modo tale che il sistema controllato abbiatutti i poli a parte reale minore di zero; lavorando sulle funzioni di trasferimentooccorre quindi aggiungere al sistema non controllato la funzione di trasferimentodel regolatore PID e studiare la stabilita del sistema in anello chiuso. La funzionedi trasferimento del regolatore PID e:

RPID(s) =kp s + kd s2 + ki

s(2.5)

Questa verra aggiunta alle due funzioni di trasferimento per l’angolo teta e laposizione del carrello.

2.2.1 PID per il controllo dell’angolo θ

La funzione di trasferimento G(s) e quella ricavata prima per θ. La funzione ditrasferimento dell’anello chiuso e:

GCL(s) =RPID(s) G(s)

1 + RPID(s)G(s)

e cioe:GCL(s) =

numPID numG

denPID denG + numPID numG

quindi simbolicamente:

GCL(s) =(n1 s2)(kp s + kd s2 + ki)

s5 + d1 s4 + d2 s3 + d3 s2 + (n1 s2)(kp s + kd s2 + ki)

Raccogliendo i termini dello stesso grado otteniamo:

GCL(s) =(n1 s2)(kp s + kd s2 + ki)

s2 [s3 + (d1 + n1 kd) s2 + (d2 + n1 kp) s + (d3 + n1 ki)]

A questo punto al denominatore compaiono i coefficienti del regolatore in mo-do indipendente; e pero possibile solo controllare tre poli su quattro potendo

Luca Rancati644292

10


agire solo sul polinomio di grado tre ed avendo tre gradi di liberta ovvero itre coefficienti. per calcolarli analiticamente si confrota la parte di terzo gradodel denominatore della funzione dell’anello chiuso con un generico polinomio diterzo grado le cui radici sono −p1,−p2,−p3:

(s + p1) (s + p2) (s + p3)

che sviluppata da:

s3 + (p1 + p2 + p3) s2 + (p1 p2 + p1 p3 + p2 p3) s + (p1 p2 p3)

Uguagliando i coefficienti che moltiplicano i termini dello stesso grado otteni-amo:

kd = [(p1+p2+p3)−d1]n1

kp = [(p1 p2+p1 p3+p2 p3)−d2]n1

kd = [(p1 p2 p3)−d3]n1

(2.6)

Quindi scegliendo i tre poli in anello chiuso p1, p2, p3 si calcolano i coefficientidel regolatore. Il polo a parte reale nulla non si riesce a spostare, quindi in anellochiuso si riesce a rendere questa parte di sistema semplicemente stabile. 1Lafunzione matlab coefpid.m calcola i coefficienti del regolatore dati in ingresso itre poli in anello chiuso.

2.2.2 PID per il controllo della posizione x del carrello

PID sulla funzione di trasferimento x(s)/f(s)

Eseguendo lo stesso procedimento fatto per l’angolo teta si ottiene:

x(s)f(s)

=z1 s2 + z2

s4 + d1 s3 + d2 s2 + d3 s

aggiungendo il regolatore:

GCL(s) =(z1 s2 + z2)(kp s + kd s2 + ki)

s5 + d1 s4 + d2 s3 + d3 s2 + (z1 s2 + z2)(kp s + kd s2 + ki)

e risolvendo si ottiene:

GCL(s) =(z1 s2 + z2)(kp s + kd s2 + ki)

[s5 + (d1 + z1 kd) s4 + (d2 + z1 kp) s3 + (d3 + z1 ki + z2 kd) s2 + z2 kp s + z2 ki]

Come si nota, i coefficienti non sono piu indipendenti; inoltre se confrontassi ildenominatore con un generico polinomio di quinto grado del tipo (s + p1) (s +p2) (s + p3) (s + p4) (s + p5) otterei cinque condizioni da verificare con solo treparametri. Quindi non verra presa questa soluzione.

1parte di sistema corrispondente alla funzione di trasferimento per teta

Luca Rancati644292

11


Compensatore per la funzione di trasferimento x(s)/f(s)

Si e quindi pensato di aggiungere solo un compensatore che avesse uno zeroche cancellasse il polo instabile della funzione di trasferimento, un polo in altafrequenza ed un guadagno:

R(s) = µ(s + z)(s + p)

(2.7)

in cui z ha la funzione di cancellare il polo instabile e p e un polo in altafrequenza. Cosı facendo pero in anello chiuso si otterrebbe un sistema che haancora un polo a parte reale maggiore di zero anche se molto vicino allo zero. Ilsistema non diverge, le specifiche sull’angolo teta sono rispettate ma il carrellooscilla continuamente; il sistema e a regime su un ciclo per valori del guadagnoprossimi a 350.

PID sulla funzione di trasferimento (x + k θ)(s)/f(s)

Questa soluzione porta a definire una matrice C che abbia la seconda riga paria [k, 1, 0, 0] con k parametro libero. La funzione di trasferimento diventa:

(x + k θ)(s)f(s)

=(a k + z1) s2 + z2

s4 + d1 s3 + d2 s2 + d3 s

con a coefficiente costante ricavato dai calcoli con Matlab. Aggiungendo ilregolatore PID si ottiene in anello chiuso:

GCL(s) =[(a k + z1) s2 + z2](kp s + kd s2 + ki)

[s5 + (d1 + z1 kd + a k kd) s4 + (d2 + z1 kp + a k kp) s3 + (d3 + z1 ki + z2 kd + a k ki) s2 + z2 kp s + z2 ki]

Come si vede cosı facendo non si puo confrontare il denominatore con ungenerico polinimio di quinto grado e scegliere i poli in anello chiuso, in quantosi avrebbero cinque condizioni da rispettare con quattro parametri. Tuttavia epossibile rendere tutti i coefficienti del denominatore concordi, condizione perosolo necessaria perche le radici siano tutte negative. le condizioni sono:

d1 + z1 kd + a k kd > 0d2 + z1 kp + a k kp > 0

d3 + z1 ki + z2 kd + a k ki > 0kp < 0ki < 0

(2.8)

dove kp e ki devono essere negativi perche z2 che li moltiplica negli ultimi duecoefficienti e negativo. Ponendo per semplicita k = 1 si ricava:

kp > −7.9483kd > −7.3555

ki < −1(2.9)

I valori scelti per la simulazione sono:

• kp = −35

• kd = −8

Luca Rancati644292

12


• ki = −0.8

Tuttavia, essendo verificata solo la condizione nacessaria, se con questi coeffi-cienti si calcolano le radici del denominatore della funzione di trasferimento inanello chiuso si scopre che e presente ancora un polo positivo e vicino allo ze-ro, cadendo cosı nel problema trovato in precedenza. Comunque il sistema conquesti valori ha oscillazioni limitate e rispetta la condizione per l’errore su teta.

PID sulla funzione di trasferimento (k1x + k2 θ)(s)/f(s)

Questa soluzione porta a definire una matrice C che abbia la seconda riga paria [k1, k2, 0, 0]. Si introducono due parametri liberi per ottenere una funzionein anello chiuso con cinque parametri, questi due e i tre del regolatore PID.Le cinque condizioni sul confronto con un polinomio di grado quinto del tipo(s+p1) (s+p2) (s+p3) (s+p4) (s+p5) da luogo, scelti i poli ad un sistema nonlineare nelle incognite k1, k2, kp, kd, ki. A questo punto si e tentato di risolverlocon il metodo di Newton per diverse condizioni iniziali ma esso non converge2.La funzione di trasferimento e:

GCL(s) =[(k1 B3 + k2 B4) + k2 z2] (kp s + kd s2 + ki)

s5 + d1 s4 + d2 s3 + d3 s + (k1 B3 + k2 B4) + k2 z2] (kp s + kd s2 + ki)

dove i termini B3 e B4 sono i coefficienti della matrice B. Confrontando ildenominatore con il polinomio di quinto grado:

(s + p1) (s + p2) (s + p3) (s + p4) (s + p5) che diventas5 + a1 s4 + b1 s3 + c1 s2 + d s + e1

dove i coefficienti sono:

a1 = p1 + p2 + p3 + p4 + p5

b1 = p1p2 + p1p3 + p1p4 + p1p5 + p2p3 + p2p4 + p2p5 + p3p4 + p3p5 + p4p5

c1 = p1p2p3 + p1p2p4 + p1p2p5 + p1p3p4 + p1p3p5 + p1p4p5 + p2p3p4 + p2p3p5 + p2p4p5 + p3p4p5

d = p1p2p3p4 + p1p2p3p5 + p1p2p4p5 + p1p3p4p5 + p2p3p4p5

e1 = p1p2p3p4p5

Si arriva uguagliando i coefficienti ad un sistema non lineare come detto inprecedenza; il sistema e:

d1 + kd k1 B3 + kd k2 B4 − a1 = 0

d2 + kp k1 B3 + kp k2 B4 − b1 = 0

d3 + ki k1 B3 + ki k2 B4 + kd k2 z2 − c1 = 0

kp k2 z2 − d = 0

ki k2 z2 − e1 = 0

(2.10)

2bisognerebbe avere una stima a priori della soluzione per dare la condizione iniziale sensata

Luca Rancati644292

13


Con il metodo di Newton3 per i sistemi non lineari non si e arrivati alla soluzioneper il vettore delle incognite [k1, k2, kp, kd, ki]. Per cui si adotta nella simulazionela soluzione PID sulla funzione di trasferimento (x + k θ)(s)/f(s) che e quellache minimizza l’oscillazione del carrello.

2.3 Modello numerico in Simulink e risultati

Il sistema fisico e rappresentato dal sistema dinamico linearizzato attorno allaposizione θ = 0,il blocco state space, che da in uscita le variabili x e θ. Questepassano nel blocco rumore e filtraggio in cui viene aggiunto il rumore e viene poifiltrato.4 Le variabili filtrate entrano nel blocco controllo in cui viene generatol’errore rispetto al riferimento nullo per le variabili θ e x + k θ data la sceltadelle funzioni di trasferimento su cui implementare il PID. I coefficienti sonocalcolati per il PID su θ tramite la funzione coefpid.m mentre per l’altra sonocalcolati rispettando i vincoli (2.6) ponendo k = 1. Per entrambi i regolatori PIDi coefficienti sono moltiplicati per l’errore, la derivata dell’errore e l’integrale.Per ottenere la derivata e l’integrale dell’errore si sono utilizzate due funzioni ditrasferimento rappresentanti uno pseudoderivatore ed uno pseudointegratore inquesto modo:

s

1/(2π f) s + 1pseudoderivatore con f = 10 [Hz]

s

s2 + 0.5 spseudointegratore

Il comando finale e la somma dei singoli comandi elaborati dai PID. Questo co-mando passa dai due blocchi dead zone e saturation come spiegato in precedenzae poi viene nuovamente filtrato sempre con frequenza di taglio f = 10 [Hz], edinfine entra nel sistema dinamico costituito dal blocco state space. Vediamol’andamento delle due variabili teta ed x:

3vedere appendice4rumore e filtri sono come spiegato in precedenza

Luca Rancati644292

14


0 10 20 30 40 50 60 70 80 90 100−1.5

−1

−0.5

0

0.5

1

1.5

2

Figura 2.1: angolo (gradi) per 100 s di simulazione, θ(0) = 1.7 gradi

Le specifiche sull’angolo sono rispettate, partendo da una condizione inizialeperturbata; anche durante il transitorio iniziale le specifiche sono rispettate:

0 1 2 3 4 5 6 7 8 9

−1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

minimo non supera i −2 gradi

Figura 2.2: transitorio iniziale per l’angolo teta

Luca Rancati644292

15


0 10 20 30 40 50 60 70 80 90 100−0.015

−0.01

−0.005

0

0.005

0.01

0.015

Figura 2.3: posizione del carrello (cm) per 100 s di simulazione

Come si vede il carrello si mantiene limitato nella sua oscillazione.Nel caso in cui si utilizzasse la soluzione di compensare semplicemente la fun-zione di trasferimento per il carrello si otterrebbero i seguenti risultati con:

R(s) = µs + z

s + p

con:

• µ = 350 guadagno

• z = −5.041 zero che cancella il polo instabile

• p = 100 polo in alta frequenza

0 50 100 150 200 250 300−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

0.06

0.08

Figura 2.4: posizione del carrello (cm) per 300 s di simulazione

Luca Rancati644292

16


0 50 100 150 200 250 300−1

−0.5

0

0.5

1

1.5

2

Figura 2.5: angolo (gradi) per 300 s di simulazione, θ(0) = 1.7 gradi

Come si vede le specifiche sull’angolo sono rispettate ed il carrello ha oscillazionipiu marcate; questo perche in anello chiuso la presenza di uno zero a parte realemaggiore di zero nella funzione di trasferimento del carrello, fa si che nonostanteil regolatore cancelli il polo positivo, il denominatore in anello chiuso presentaancora un polo positivo ma molto vicino allo zero. Provando ad aggiungereun PID lasciando il guadagno come parametro libero si incappa nel problemadi avere cinque condizioni da rispettare con quattro parametri. Si riesce asoddisfare solo la condizione necessaria5.

5sempre confrontando il denominatore con un polinomio di grado quinto e uguagliando icoefficienti

Luca Rancati644292

17

Capitolo 3

Il modello in laboratorio

Il modello linearizzato rappresentato dal blocco state space e sostituito dall’es-perimento vero e proprio. Le variabili x e teta vengono misurate dai sensori;attraverso una scheda di acquisizione dati i segnali entrano nel modello simulinkche si interfaccia con gli ingressi della scheda con il blocco analog input ; questisono volt e vengono convertiti in radianti e metri moltiplicando per le sensibilitaricavate nei capitoli precedenti. Il modello e poi analogo a quello numerico; ilcomando calcolato dal controllo entra in un blocco analog output che attraver-so la scheda manda il comando al motore elettrico del carrello. E necessariointrodurre un blocco time based collegato ad un terminator. Il rumore e giapresente nel segnale quindi sono presenti solo i filtri. Il resto e uguale a prima.Si riportano gli andamenti delle variabili:

0 5 10 15 20 25 30 35 40−4

−2

0

2

4

6

8

Figura 3.1: angolo (gradi), ottenuto in laboratorio

Le specifiche sono rispettate anche perturbando periodicamente l’asta. 1

1dando un colpo all’asta

18


0 5 10 15 20 25 30 35 40−10

−5

0

5

10

15

20

Figura 3.2: posizione del carrello (cm) ottenuta in laboratorio

Luca Rancati644292

19

Capitolo 4

Appendice

4.1 determinazione dell’attrito

Si valuta il coefficiente c dovuto all’attrito. E stato valutato per confronto ecioe: in laboratorio si e dato in ingresso al carrello un segnale a scalino diampiezza 2 e si e registrata la posizione x. E poi stato fatto un modello numericoche simulasse l’esperimento in precedenza descritto. tale modello prevede chel’attrito sia proporzionale alla velocita del carrello ed e:

mtot x + c x = fx (4.1)

Sostituendo l’espressione trovata per la forza si ottiene:

mtot x + (c− α) x = β V

dove α e β sono i parametri del motore noti dal costruttore. si e variato ilcoefficiente c fino a che le risposte del sistema simulato ben rappresentasseroil sistema reale di laboratorio. Scegliendo c pari a 4.1 [N/(m/s)] si ottiene ilrisultato cercato.

4.2 Metodo di Newton

Si e utilizzata la funzione qssnewtonxs1 per i sistemi non lineari di ordine N.

4.3 Funzioni matlab

Nel floppy sono contenute le due funzioni matlab utilizzate per il progetto:

• datip.m carica le costanti utilizzate2, crea le matrici del sistema lineariz-zato, crea il sistema dinamico e calcola gli autovalori; effettua il test di os-servabilita e controllabilita e stampa a video gli autovalori ed un messaggioche indica se il sistema e osservabile e controllabile.

1funzione utilizzata nel corso di Calcolo Numerico2dati fisici del costruttore

20


• coefpid.m dati i tre poli in anello chiuso per la funzione di trasferimentoθ(s)/f(s) calcola i coefficienti del regolatore PID.

Sono anche presenti i files dei modelli simulink.

Luca Rancati644292

21

Sistema di controllo per il pendolo inverso -...

Documents

Transcript of Sistema di controllo per il pendolo inverso -...