Corso di Laurea Magistrale in Ingegneria Elettrica e dell

Corso di Laurea Magistrale in Ingegneria Elettrica e
dell’Automazione
note alle lezioni di Modellistica e Controllo di Sistemi
Meccanici (6CFU)
bozza del 10 dicembre 2010
Benedetto Allotta
1
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
2
Modellistica e Controllo di Sistemi Meccanici
Contenuti
Copertina
1
Indice
3
1 Sistemi anolonomi
1.1 Richiami di dinamica dei sistemi meccanici . . . . . . . . . . . . .
1.1.1 Equazioni di d’Alembert-Lagrange . . . . . . . . . . . . . .
1.1.2 Equazioni di Lagrange del II tipo . . . . . . . . . . . . . . .
1.2 Sistemi anolonomi . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Elementi di geometria differenziale . . . . . . . . . . . . . .
1.2.1.1 Varietà differenziabili . . . . . . . . . . . . . . . .
1.2.1.2 Spazio tangente . . . . . . . . . . . . . . . . . . .
1.2.1.3 Spazio duale . . . . . . . . . . . . . . . . . . . . .
1.2.1.4 Spazio cotangente . . . . . . . . . . . . . . . . . .
1.2.1.5 Campo vettoriale . . . . . . . . . . . . . . . . . . .
1.2.1.6 Campo covettoriale . . . . . . . . . . . . . . . . .
1.2.1.7 Distribuzione . . . . . . . . . . . . . . . . . . . . .
1.2.1.8 Codistribuzione . . . . . . . . . . . . . . . . . . .
1.2.1.9 Parentesi di Lie (Lie bracket) . . . . . . . . . . . .
1.2.1.10 Vincoli olonomi . . . . . . . . . . . . . . . . . . . .
1.2.1.11 Vincoli pfaffiani . . . . . . . . . . . . . . . . . . . .
1.2.1.12 Vincoli pfaffiani non integrabili (vincoli anolonomi)
1.2.1.13 Derivata di Lie . . . . . . . . . . . . . . . . . . . .
1.2.1.14 Distribuzioni completamente integrabili . . . . . .
1.2.1.15 Distribuzioni involutive . . . . . . . . . . . . . . . .
1.2.2 Importanti teoremi sui sistemi anolonomi . . . . . . . . . .
1.2.2.1 Teorema di Frobenius (uno dei tanti!) . . . . . . .
1.2.2.2 Teorema di Chow sui sistemi driftless . . . . . . .
1.2.2.3 Teorema di Brockett . . . . . . . . . . . . . . . . .
1.2.3 Vincoli anolonomi . . . . . . . . . . . . . . . . . . . . . . .
1.3 Esempi di sistemi anolonomi . . . . . . . . . . . . . . . . . . . . .
1.3.1 Modello cinematico dell’uniciclo su piano . . . . . . . . . .
1.3.2 Controllo dell’uniciclo . . . . . . . . . . . . . . . . . . . . .
1.3.3 Sfera che rotola su un piano . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
5
6
6
6
6
7
7
7
7
7
7
7
8
8
8
9
9
9
9
9
9
10
10
10
10
10
12
14
2 Visual servoing
2.1 Sistemi di visione . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 L’immagine digitale . . . . . . . . . . . . . . . . . . .
2.2 Telecamere . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Modelli di telecamera . . . . . . . . . . . . . . . . .
2.2.1.1 Il modello Basic Pinhole . . . . . . . . . . .
2.2.1.2 il modello CCD Full Perspective . . . . . .
2.2.1.3 il modello CCD con distorsione della lente
2.2.2 Calibrazione della telecamera . . . . . . . . . . . . .
2.3 Elementi di Geometria Proiettiva . . . . . . . . . . . . . . .
2.3.1 Primitive nello spazio immagine . . . . . . . . . . . .
2.3.1.1 Rappresentazione omogenea delle rette .
2.3.1.2 Rappresentazione omogenea dei punti . .
2.3.1.3 Gradi di libertà . . . . . . . . . . . . . . . .
2.3.1.4 Intersezione tra rette . . . . . . . . . . . .
2.3.1.5 Intersezione tre due rette parallele . . . . .
2.3.1.6 Punti ideali e linea all’infinito . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
17
17
18
19
20
21
23
24
25
27
27
27
28
28
28
28
28
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
3
Modellistica e Controllo di Sistemi Meccanici
2.3.1.7 Un modello per il piano proiettivo . . . . . .
2.3.2 Geometria di due immagini . . . . . . . . . . . . . .
2.3.2.1 Omografie . . . . . . . . . . . . . . . . . .
2.3.2.2 Trasformazioni di rette . . . . . . . . . . . .
2.3.3 Classificazione delle proiettività . . . . . . . . . . . .
2.3.3.1 Isometrie . . . . . . . . . . . . . . . . . . .
2.3.3.2 Similitudini . . . . . . . . . . . . . . . . . .
2.3.3.3 Affinità . . . . . . . . . . . . . . . . . . . .
2.3.3.4 Proiettività . . . . . . . . . . . . . . . . . .
2.3.3.5 Decomposizione di una proiettività . . . . .
2.3.3.6 Topologia del piano proiettivo . . . . . . . .
2.3.3.7 Punti fissi e linee fisse . . . . . . . . . . . .
2.3.3.8 Calcolo Omografia tra 2 immagini . . . . .
2.3.3.9 Structure from Motion . . . . . . . . . . . .
2.4 Asservimento visivo di robot . . . . . . . . . . . . . . . . . .
2.4.1 Classificazione di sistemi con asservimento visivo .
2.4.1.1 Position Based Visual Servoing (PBVS) . .
2.4.1.2 Image Based Visual Servoing (IBVS) . . .
2.4.1.3 Hybrid Visual Servoing . . . . . . . . . . .
2.4.2 Controllo dei sistemi ad Asservimento Visivo . . . .
2.4.3 Stabilità del Sistema . . . . . . . . . . . . . . . . . .
2.4.4 I Punti Immagine per l’identificazione dello Stato . .
2.4.4.1 Matrice di Interazione del punto immagine
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
29
29
29
30
31
31
31
32
32
33
33
34
35
36
37
38
38
39
39
40
41
42
42
3 Controllo sliding mode di sistemi in forma canonica di controllo
3.1 Controllo sliding mode di sistemi ad un solo ingresso . . . . . . . . .
3.1.1 Definizione della superficie di sliding . . . . . . . . . . . . . .
3.1.2 Sintesi della legge di controllo . . . . . . . . . . . . . . . . .
3.1.2.1 Parte model-based . . . . . . . . . . . . . . . . . .
3.1.2.2 Parte discontinua su S(t) . . . . . . . . . . . . . . .
3.1.2.3 Verifica della condizione di sliding . . . . . . . . . .
3.2 Controllo sliding mode di sistemi multivariabile . . . . . . . . . . . .
3.2.1 Variabili di stato ed ingressi . . . . . . . . . . . . . . . . . . .
3.2.2 Modello dinamico . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Modelli di incertezza . . . . . . . . . . . . . . . . . . . . . . .
3.2.4 Variabili di sliding . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.5 Superficie di sliding . . . . . . . . . . . . . . . . . . . . . . .
3.2.6 Progettazione del controllo . . . . . . . . . . . . . . . . . . .
3.2.7 Controllo basato sul modello . . . . . . . . . . . . . . . . . .
3.2.7.1 Definizione delle “variabili di riferimento” . . . . . .
3.2.7.2 Condizione di permanenza sulla superficie di sliding
3.2.8 Controllo robusto (parte discontinua) . . . . . . . . . . . . . .
3.2.9 Verifica delle condizioni di sliding . . . . . . . . . . . . . . . .
3.3 Un altro approccio al progetto del controllore SMC . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
45
45
45
46
46
47
47
49
49
49
50
50
50
50
51
51
51
51
52
53
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
4
Lista delle Figure
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Punto materiale in moto verticale sotto l’effetto della gravità. . . . . . . . . . . . .
Uniciclo che rotola su un piano. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Parcheggio dell’uniciclo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Possibili diverse manovre di parcheggio di un uniciclo. . . . . . . . . . . . . . . .
Esempio di sistema meccanico con vincoli anolonomi: sfera che rotola su un piano.
Sfera che rotola su un piano: proiezione ortogonale sul piano xz. . . . . . . . . .
Sfera che rotola su un piano: proiezione ortogonale sul piano yz. . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
Modello “pinhole”. L’immagine di un punto sullo schermo si forma lasciando passare un unico raggio attraverso un foro di piccole dimensioni praticato in uno schermo
opaco che funge da otturatore. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Modello “pinhole” di telecamera. Sono evidenziati: il centro della telecamera oc , il
punto principale dell’immagine, l’origine del sistema di riferimento immagine . . . .
Piano immagine con non perfetta ortogonalità tra gli assi xim e yim . . . . . . . . .
Piano immagine con non perfetta ortogonalità tra gli assi xim e yim . . . . . . . . .
Un modello del piano proiettivo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Proiezione centrale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Un altro esempio di proiettività . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
trasformazioni nell’affinità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Inversione dell’orientazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Punti fissi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Variabile di scivolamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
11
13
14
14
15
15
21
22
24
25
29
30
31
33
34
34
54
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
1
1.1
1.1.1
Modellistica e Controllo di Sistemi Meccanici
5
Sistemi anolonomi
Richiami di dinamica dei sistemi meccanici
Equazioni di d’Alembert-Lagrange
Per un sistema meccanico a n gradi di libertà (DOF) è possibile scrivere l’equazione di d’AlembertLagrange nel modo seguente:
n X
d ∂T
∂T
−
− Qk δqk
dt ∂ q̇k
∂qk
=
0 ,
(1)
k=1
T
dove q = [q1 q2 . . . qn ] è il vettore delle coordinate lagrangiane, T è, come di consueto, l’energia
cinetica del sistema, Qk è la k-esima componente lagrangiana delle forze attive (interne ed esterne) agenti sul sistema, e δqk è la k-esima componente di uno spostamento virtuale δq =
[δq1 δq2 . . . δqn ]T . Uno spostamento virtuale è uno spostamento del sistema nello spazio di
configurazione che sia compatibile con i vincoli esistenti.
Dato che nella (1) compare l’energia cinetica T invece della funzione lagrangiana L = T − V
(V è l’energia potenziale elastica e/o gravitazionale del sistema), sarà necessario tener conto
esplicitamente delle forze attive conservative nel termine Qk , ovvero:
Qk
=
c
Qnc
k + Qk
(2)
,
c
dove Qnc
k è la k-esima componente lagrangiana delle forze attive non conservative e Qk è la kesima componente lagrangiana delle forze attive conservative. La componente conservativa sarà
∂V
. Tenendo conto delle azioni conservative all’interno della funzione
calcolabile come: Qck = − ∂q
k
lagrangiana, è possibile, in alternativa, riscrivere la (1) come segue:
n X
d ∂L
∂L
nc
−
− Qk δqk
dt ∂ q̇k
∂qk
=
0 .
(3)
k=1
Esempio 1.1: Per un punto materiale di massa m soggetto al campo gravitazionale terrestre e vincolato a muoversi
verticalmente senza attrito, possiamo scegliere come sua unica coordinata lagrangiana q la quota del punto rispetto ad un
d ∂T
= mq̈. L’unica
riferimento fisso, con verso positivo verso l’alto. L’energia cinetica del punto vale T = 21 mq̇ 2 , da cui dt
∂ q̇
q
m
g
Figura 1: Punto materiale in moto verticale sotto l’effetto della gravità.
forza attiva agente sul sistema è la gravità, associata al campo (conservativo) gravitazionale, quindi avremo Q = Qc =
∂
∂
− ∂V
= − ∂q
V = − ∂q
(mgq) = −mg. L’equazione di d’Alembert-Lagrange diventa allora: (mq̈ + mg) δq = 0, da cui,
∂q
dovendo essere quest’ultima equazione valida ∀δq, deriva l’equazione del moto del sistema mq̈ = −mg.
In alternativa, essendo la gravità l’unica forza attiva agente sul sistema, se ne può tenere conto utilizzando l’Eq. (3).
d ∂L
Avremo quindi L = T − V = 12 mq̇ 2 − mgq, dt
= mq̈, ∂L
= mg e infine l’equazione del moto del sistema risulterà
∂ q̇
∂q
uguale a quella già ricavata, ovvero: mq̈ + mg = 0 cioè q̈ = −g.
Allotta, Fioravanti, Malvezzi
1.1.2
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
6
Equazioni di Lagrange del II tipo
Nei sistemi olonomi le componenti di un quasiasi spostamento virtuale δq = [δq1 δq2 . . . δqn ]T
sono tra loro indipendenti. Per ogni termine della sommatoria che compare nell’Eq. (1), deve
valere:
d ∂T
∂T
−
− Qk
dt ∂ q̇k
∂qk
=
0 ∀k = 1, . . . , n
,
(4)
∂L
d ∂L
−
− Qnc
k
dt ∂ q̇k
∂qk
=
0
∀k = 1, . . . , n
.
(5)
o, in alternativa:
T T
Esempio 1.2: Nel caso di un manipolatore a cui sia applicato1 un wrench di forze −γ i = [−f T
i − µi ] agente sul
nc
link i-esimo, il calcolo delle componenti lagrangiane delle forze attive Qk si effettua calcolando lo jacobiano J (li ) ed
(l )
(l )
i T
i
utilizzandolo come proiettore nel modo che segue: Qnc
è la k-esima colonna di J (li ) e
k = −(J k ) γ i + τk , dove J k
τk è, come di consueto, l’azione degli attuatori sul giunto k.
1.2
Sistemi anolonomi
In questo paragrafo, dopo aver introdotto alcuni concetti di geometria differenziale, opereremo
una classificazione dei vincoli, definendo in modo rigoroso i vincoli olonomi e quelli anolonomi.
1.2.1
Elementi di geometria differenziale
1.2.1.1 Varietà differenziabili Una varietà differenziabile (in Inglese differentiable manifold)
m-dimensionale è uno spazio topologico localmente diffeomorfo allo spazio euclideo Rm 2 In breve
diremo varietà. Per gli scopi di questo corso una Varietà M è un sottoinsieme di Rn (embedded
submanifold) definito dalla eguaglianza a zero di una funzione vettoriale h : Rn −→ Rp con p < n.
La funzione vettoriale h(x) è definita come segue
h(x)
=
T
[h1 (x) h1 (x) . . . hp (x)]
,
(6)
e la varietà differenziabile M , di dimensione m = n − p, è definita dalle equazioni:
h1 (x1 , x2 , . . . , xn )
=
0
h2 (x1 , x2 , . . . , xn )
=
0
...
hp (x1 , x2 , . . . , xn )
=
0 .
(7)
Ipotizzeremo che i p gradienti delle funzioni h1 (x), . . . , hp (x) siano tra loro linearmente indipenp×n
, definita come segue
denti, ovvero che la matrice jacobiana ∂∂ h
x ∈R
 ∂h

∂h1
∂h1
1
. . . ∂x
∂x1
∂x2
n




 ∂h2 ∂h2 . . . ∂h2 


∂x
∂x
∂x
2
2
n
∂h


= 
(8)



∂x
 ... ... ... ... 




∂hp
∂hp
∂hp
.
.
.
∂x1
∂x2
∂xn
abbia rango pieno (di riga), ovvero rango p.
1 La convenzione utilizzata anche nell’ambito del corso di Robotica Industriale prevede che γ sia il wrench applicato
dall’end-effector del manipolatore sull’esterno, quindi il wrench agente sull’end-effector del manipolatore sarà dato da
−γ. L’equazione (vettoriale) del moto, in presenza di un wrench di forze di interazione γ i tra il link i-esimo e l’ambiente,
risulterà: M (q)q̈ + h(q, q̇) + J (li) γ = τ ovvero M (q)q̈ + h(q, q̇) = −J (li) γ + τ . Il secondo membro membro di
quest’ultima equazione equivale al vettore delle componenti lagrangiane delle forze attive non conservative.
2 Un diffeomorfismo è una funzione differenziabile ed invertibile la cui inversa è anch’essa differenziabile.
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
7
1.2.1.2 Spazio tangente Data una varietà differenziabile M è possibile definire in ogni punto
x uno spazio tangente in x che indicheremo con il simbolo Tx M . Lo spazio tangente in x è
l’insieme dei vettori lungo cui x si può muovere rimanendo appartenente ad M , ovvero l’insieme
dei vettori tangenti alla varietà o, in altri termini, il complemento ortogonale dei gradienti definiti in
(8).
Esempio 1.3: La sfera unitaria S 2 definita dall’unica equazione h(x, y, z) = x2 + y 2 + z 2 − 1 = 0 è una (sotto)varietà
bidimensionale
dello spazio R3 e scriveremo che S 2 ⊂ R3 . Nella semisfera superiore (z > 0) possiamo scrivere:
p
2
z = 1 − x − y 2 . Lo spazio tangente alla sfera in un suo punto può essere descritto, ad esempio, dal range dei vettori
∂z T
∂z T
] = [1 0 √ −x2 2 ]T , v 2 = [0 1 ∂y
] = [0 1 √ −y2 2 ]T . Osserviamo che
v 1 , v 2 definiti come segue: v 1 = [1 0 ∂x
1−x −y
1−x −y
p
∂h ∂h
il gradiente di h, ovvero [ ∂h
] = [2x 2y 2z] = [2x 2y 2 1 − x2 − y 2 ] è ortogonale sia a v 1 v 2 .
∂x ∂x ∂z
1.2.1.3 Spazio duale Lo spazio duale U di uno spazio vettoriale dato V è, per i nostri scopi,
uno spazio vettoriale della stessa dimensione dello spazio vettoriale dato e per il quale possiamo
individuare una base ui = v Ti , ∀i = 1, . . . , m, dove v i è il generico elemento di una base
di V . Ad esempio, lo spazio duale di Rn (ricordiamo che Rn è lo spazio vettoriale dei vettori
colonna formati da n numeri reali) è lo spazio vettoriale dei vettori riga formati da n numeri reali.
Indicheremo il duale di Rn con Rn∗
1.2.1.4 Spazio cotangente Lo spazio cotangente ad una varietà M in suo punto x è lo spazio
∗
duale dello spazio tangente Tx M . Indicheremo lo spazio cotangente con il simbolo Tx
M . Per gli
scopi del corso assumeremo che Tx M sia lo spazio vettoriale formato dai vettori colonna tangenti
∗
M sia lo spazio vettoriale formato dai vettori riga tangenti a M in x.
a M in x mentre Tx
1.2.1.5 Campo vettoriale Un campo vettoriale su una varietà M è una applicazione f : M −→
Tx M (ovviamente con valori in Rn ):


f1 (x)
 f2 (x) 

(9)
f (x) = 
 ... 
fn (x)
1.2.1.6 Campo covettoriale Un campo covettoriale su una varietà M è una applicazione ω T :
∗
M −→ Tx
M (quindi con valori in R∗ n ):
ω T (x)
=
[ω1 (x) ω2 (x) . . . ωn (x)]
(10)
1.2.1.7 Distribuzione Siano X 1 (x), X 2 (x), . . . , X k (x) dei campi vettoriali definiti su M ∀x ∈
M . Definiamo distribuzione, e la indicheremo col simbolo ∆, l’insieme delle possibili combinazioni
lineari dei campi vettoriali X i , ∀i = 1, . . . , k:
∆
span {X 1 , X 2 , . . . , X k }
X1 X2 . . . Xk
= range
.
=
(11)
In sintesi una distribuzione assegna ad ogni punto x ∈ M della varietà uno spazio vettoriale
∆(x). ∆(x) è un sottospazio k-dimensionale dello spazio tangente m-dimensionale Tx M .
1.2.1.8 Codistribuzione Siano ω T1 (x), ω T2 (x), . . . , ω Tk (x) dei campi covettoriali definiti su M ∀x ∈
M . Definiamo codistribuzione, e la indicheremo col simbolo Ω, l’insieme delle possibili combinazioni lineari dei campi covettoriali X i , ∀i = 1, . . . , k:
Ω = span ω T1 , ω T2 , . . . , ω Tk
 T 
ω1
 ω T2 


= rowrange 
(12)
 . . .  .
T
ωk
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
8
Modellistica e Controllo di Sistemi Meccanici
In sintesi una codistribuzione assegna ad ogni punto x ∈ M della varietà uno spazio covettoriale
∗
Ω(x). Ω(x) è un sottospazio k-dimensionale dello spazio cotangente m-dimensionale Tx
M.
1.2.1.9 Parentesi di Lie (Lie bracket) Siano f (x) e g(x) due campi vettoriali definiti su Rn .
La Lie bracket di f (x) e g(x), che indicheremo con [f , g] è a sua volta un campo vettoriale definito
come segue:
[f , g]
∂f
∂g
f−
g
∂x
∂x
=
(13)
,
∂f
∂g
dove ∂ x e ∂ x sono gli jacobiani di f e g rispetto a x, ovvero:
∂f
∂fi (x)
∂g
∂gi (x)
=
=
∂x ij
∂xj
∂x ij
∂xj

Esempio 1.4: Calcoliamo la Lie bracket dei campi vettoriali, definiti su
[f , g]
=
=
=
R3 ,
f (x) = 
∂g
∂f
f−
g=
∂x
∂x


 
x2
0
0
0
0
1
 0 2x2 0   sin x1  −  cos x1 0
0
0
0
1
0
x1 + x23

  2  

0
x2
−x22
 2x2 sin x1  −  0  =  2x2 sin x1 
0
2x3
−2x3



0
x2
2


sin x1
x2  .
g(x) =
1
x1 + x23


0
0
0   x22 
2x3
1
.
Notiamo che [f , g] = −[g, f ], infatti:
[g, f ]
∂g
∂f
g−
f=
∂x
∂x
−[f , g].
=
=
(14)
1.2.1.10 Vincoli olonomi Supponiamo di avere un sistema meccanico la cui configurazione è
descritta da n coordinate lagrangiane soggette a p < n vincoli del tipo:
hi (q1 , q2 , . . . , qn )
=
0 i = 1, . . . , p
.
(15)
Un vincolo del tipo (15) si dice olonomo. Facciamo l’ipotesi che i gradienti delle funzioni hi siano
tra loro linearmente indipendenti. Derivando la (15) rispetto al tempo otterremo:
dhi
dt
∂hi
q̇ = 0 i = 1, . . . , p
∂q
=
,
(16)
ovvero i gradienti delle funzioni hi che definiscono i vincoli olonomi del sistema sono ortogonali al
vettore delle velocità q̇.
1.2.1.11
Vincoli pfaffiani Un vincolo sulle velocità del tipo:
ω Ti q̇
=
(17)
0
i
si dice pfaffiano. Se esiste una funzione hi tale che ω Ti = ∂h
∂ q , allora la (17) potrà essere scritta:
ω Ti q̇
=
∂hi
dhi
q̇ =
=0 .
∂q
dt
(18)
Integrando la (18) rispetto al tempo otterremo
hi
=
costante
,
(19)
ovvero un vincolo di tipo olonomo. Il vincolo pfaffiano sarà integrabile e, quindi, a tutti gli effetti
equivalente ad un vincolo olonomo.
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
9
Modellistica e Controllo di Sistemi Meccanici
1.2.1.12 Vincoli pfaffiani non integrabili (vincoli anolonomi) Se invece non esiste alcuna
T
i
funzione hi (q) tale che ω Ti = ∂h
∂ q allora il vincolo pfaffiano ω i q̇ = 0 si definisce non integrabile o
anolonomo.
Esempio 1.5: Come vedremo nel Par. 1.3.3, una sfera vincolata a ruotare senza strisciare sul piano xy è soggetta ai
vincoli:
ẋ − Rωy
=
0
ẏ + Rωx
=
0
,
dove R è il raggio della sfera, ẋ e ẏ sono le componenti non nulle della velocità di traslazione del centro della sfera e
ωx , ωy sono le componenti lungo x e lungo y della velocità angolare della sfera. Mentre è banale trovare delle funzioni le
cui derivate temporali siano rispettivamente uguali a ẋ e ẏ, è impossibile trovare delle funzioni le cui derivate siano pari
alle componenti di ω in quanto, a meno del caso particolare di moto piano di un corpo rigido, la velocità angolare non è la
derivata temporale di alcun vettore.
1.2.1.13 Derivata di Lie Sia f (x) un campo vettoriale f : Rn → Rn e sia h(x) una funzione
scalare. Si definisce derivata di Lie di h(x) lungo f (x) la seguente espressione:
∂h(x)
f (x) .
∂x
Lf h =
(20)
In altri termini, la derivata di Lie è la derivata direzionale della funzione h(x) lungo la direzione di
f (moltiplicata per il modulo di f ) e si ottiene calcolando il prodotto scalare tra il gradiente di h e
f.
1.2.1.14 Distribuzioni completamente integrabili Una distribuzione
∆ = span {X 1 , X 2 , . . . , X m } su Rn (m < n) si dice completamente integrabile se esistono n − m
funzioni scalari hj , j = 1, . . . , n − m, tra loro indipendenti, che soddisfano il seguente sistema
di equazioni differenziali:
LXi hj
=
∂hj (x)
X i = 0 ∀i = 1, . . . , m;
∂x
∀j = 1, . . . , n − m
.
(21)
In un sistema dinamico driftless in cui gli ingressi ui , i = 1, . . . , m sviluppano la loro azione
lungo le direzioni dei campi vettoriali d’ingresso X 1 , X 2 , . . . , X m come segue:
ẋ
=
m
X
(x ∈ Rn ,
X i ui
m < n) ,
(22)
i=1
il fatto che la distribuzione ∆ = span {X 1 , X 2 , . . . , X m } sia completamente integrabile comporta
che le direzioni di movimento dello stato rimangono confinate in una sottovarietà di Rn di dimensione m in quanto l’esistenza delle n − m funzioni hj , j = 1, . . . , n − m comporta dei vincoli di
tipo olonomo sulle variabili di stato. Infatti le funzioni hj rimangono costanti lungo le direzioni di
possibile evoluzione dello stato (cioè hj (x = cj , cj = cost) e possiamo agevolmente costruire
delle funzioni h0j = hj − cj tali che h0j = 0, cioè come nel Par. 1.2.1.10.
1.2.1.15 Distribuzioni involutive Si dice involutiva una distribuzione
∆ = span {X 1 , X 2 , . . . , X m } che è chiusa rispetto all’operazione parentesi di Lie, ovvero se
∀X i ∈ ∆, ∀X j ∈ ∆ risulta [X i , X j ] ∈ ∆.
Detto in altri termini, la distribuzione ∆ si dice involutiva se esistono delle funzioni scalari
αijk (x) tali che, ∀i, j vale:
[X i , X j ]
=
m
X
αijk X k
.
(23)
k=1
1.2.2
Importanti teoremi sui sistemi anolonomi
1.2.2.1 Teorema di Frobenius (uno dei tanti!) Condizione necessaria e sufficiente affinché
una distribuzione ∆ sia integrabile è che la distribuzione sia involutiva.
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
10
1.2.2.2 Teorema di Chow sui sistemi driftless Il teorema di Chow fornisce la condizione
necessaria e sufficiente per la controllabilità dei sistemi driftless. Dato il sistema (22) esso risulta
completamente controllabile sulla varietà Rn se la chiusura involutiva della distribuzione ∆ degli m
(m < n) campi vettoriali in ingresso ha dimensione n. In altri termini il sistema sarà controllabile
se esso è completamente anolonomo. In un’ottica controllistica, quindi, la presenza di vincoli
anolonomi è un fatto positivo in quanto consente di controllare completamente un sistema che
si muove di uno spazio di configurazione n-dimensionale tramite un numero m (inferiore a n) di
controlli.
Il teorema di Chow ci dice quando un sistema driftless è controllabile ma non ci dice come
controllarlo. Esistono in letteratura varie classi di metodi che sono stati sviluppati per il controllo
di sistemi anolonomi, tra questi:
• controllo ottimo;
• utilizzo di funzioni periodiche;
• utilizzo di funzioni costanti a tratti.
1.2.2.3
1.2.3
Teorema di Brockett
Vincoli anolonomi
Se nel sistema esistono n − m vincoli non integrabili del tipo:
n
X
aik q̇k
=
0 ∀i = 1, . . . , n − m
,
(24)
k=1
allora le componenti di un qualsiasi spostamento virtuale δq sono tra loro legate dal sistema
lineare:
n
X
aik δqk
=
0 ∀i = 1, . . . , n − m
.
(25)
k=1
Per ottenere le equazioni del moto del sistema a partire dalla (1) si può scegliere un sottoinsieme
di m coordinate lagrangiane indipendenti e procedere poi per eliminazione diretta delle altre n−m
coordinate lagrangiane dipendenti, sfruttando le Eq. (25). Il sistema ha m gradi di mobilità.
Esiste anche una tecnica alternativa per ricavare le equazioni del moto del sistema facendo
uso del metodo dei moltiplicatori di Lagrange.
1.3
1.3.1
Esempi di sistemi anolonomi
Modello cinematico dell’uniciclo su piano
L’uniciclo è un disco di raggio R che rotola senza strisciare su una superficie senza inclinarsi
lateralmente, come mostrato in Fig. 2. Per descriverne la configurazione possiamo usare, come
coordinate lagrangiane, la posizione x, y del punto di contatto con la superficie (che supponiamo
planare) su cui esso rotola, l’orientazione ϑ rispetto ad uno degli assi coordinati del piano e
l’angolo di rotolamento ϕ. Il sistema è soggetto al vincolo di tipo olonomo z = R, dove z è
la coordinata z del centro del disco. I vincoli anolonomi cui il sistema è soggetto sono il puro
rotolamento in avanti e l’impedita traslazione laterale. Il puro rotolamento impone che il modulo
T
della velocità di traslazione del centro del disco ẋ ẏ 0
, sia pari a Rϕ̇, ovvero:


cos ϑ
rϕ̇ − ẋ ẏ 0  sin ϑ  =
0
= rϕ̇ − ẋ cos ϑ − ẏ sin ϑ = 0 .
(26)
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
z
R
O
φ
x
θ
y
O
x
Figura 2: Uniciclo che rotola su un piano.
11
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
12
L’impedita traslazione laterale del disco può essere imposta tramite l’ortogonalità tra la velocità di
T
traslazione del centro del disco ed il versore sin ϑ − cos ϑ 0
, ovvero:
ẋ sin ϑ − ẏ cos ϑ = 0 .
(27)
I vincoli pfaffiani (26) e (27) vincolano la mobilità dell’uniciclo ma, come vedremo, sono di tipo
non integrabile. Quindi il sistema risulta essere anolonomo e, di conseguenza, controllabile,
utilizzando come ingressi u1 = ϕ̇ e u2 = ϑ̇.
Come ingressi del sistema utilizzeremo invece u1 = ϕ̇ e u2 = ϑ̇. Combinando la (26) e la (27)
otteniamo le equazioni di stato del sistema:
ẋ
= R cos ϑ u1
ẏ
= R sin ϑ u1
ϕ̇ = u1
ϑ̇
(28)
= u2
Per semplicità supponiamo adesso di non essere interessati a distinguere l’angolo di rotolamento
dell’uniciclo ϕ. In tal caso la dimensione dello spazio di configurazione si restringe a 3. Come
coordinate useremo x, y, e ϑ. Le equazioni di stato del sistema si riducono come segue:
Definendo il vettore di stato x =
nella forma:
x1
ẋ =
R cos ϑ u1
ẏ
=
R sin ϑ u1
ϑ̇
=
u2
T
x2
x3
(29)
=
x
y
ϑ
T
, possiamo scrivere il sistema
= g 1 u1 + g 1 u2 ,
(30)
T
dove g 1 = R cos ϑ R sin ϑ 0
e g2 = 0 0 1
. Verifichiamo che il sistema è controllabile calcolando la Lie Bracket [g 1 , g 2 ]:
ẋ
T
[g 1 , g 2 ]
∂g 1
∂g 2
g1 −
g =
∂x
∂x 2

 
0 0 0
R cos ϑ
0 0
=  0 0 0   R sin ϑ  −  0 0
0 0 0
0
0 0


−R sin ϑ
=  R cos ϑ  .
0
=


−R sin ϑ
0
R cos ϑ   0  =
0
1
(31)
¯
Come si vede dalla (31), [g 1 , g 2 ] è ortogonale sia a g 1 che a g 2 e quindi la chiusura involutiva ∆
della distribuzione generata dai campi vettoriali d’ingresso ∆ = span{g 1 , g 2 } ha dimensione 3 ed
il sistema risulta controllabile per il teorema di Chow. Il teorema di Chow non ci dice però come
controllare il sistema. Un possibile metodo è illustrato nel seguito.
1.3.2
Controllo dell’uniciclo
Per facilitare la sintesi di segnali di controllo idonei al controllo dell’uniciclo, operiamo un cambiamento delle variabili di stato e di quelle d’ingresso. Definiamo le variabili di stato x1 = x,
x2 = tan ϑ, x3 = y, e i segnali di controllo v1 = R cos ϑ u1 = R cos x2 u1 , v2 = cosu22 ϑ . Ricavando
v1
2
u1 e u2 in funzione delle nuove variabili otteniamo u1 = R cos
ϑ e u2 = v2 cos ϑ.
Con le nuove variabili le equazioni di stato del sistema si modificano come segue:
ẋ1
ẋ2
ẋ3
= v1
d
1
1
=
(tan ϑ) =
ϑ̇ =
u2 = v2
2
dt
cos ϑ
cos2 ϑ
R sin ϑ
= R sin ϑu1 =
v1 = x 2 v1 ,
R cos ϑ
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
13
Modellistica e Controllo di Sistemi Meccanici
ovvero:
ẋ1
= v1
ẋ2
= v2
ẋ3
= x 2 v1
(32)
.
Dalle prime due delle Eq. (32) risulta evidente che le x1 e x2 possono essere condotte con facilità
a qualsiasi valore desiderato. Una possibile strategia di controllo può allora prevedere una prima
fase in cui vengono raggiunti i valori desiderati di x1 e x2 , senza curarsi di x3 . Successivamente,
stabilire un tempo di completamento del compito T e sottoporre il sistema ad opportuni ingressi
v1 e v2 di tipo periodico, ad esempio sinusoidale, con periodo T (o un sottomultiplo di T ) in modo
che, dopo un tempo t = T , x1 e x2 tornino ai valori desiderati mentre anche x3 raggiunge il
setpoint.
y
x0 y0
O
xd yd
x
Figura 3: Parcheggio dell’uniciclo.
• configurazione iniziale: x(0) = 0, y(0) = y0 > 0, ϑ(0) = ϑ0 = 0;
• configurazione desiderata: xd = 0, yd = 0, ϑd = 0.
Applichiamo al sistema (32) dei segnali di controllo di tipo sinusoidale con pulsazione ω =
v1 (t)
= a sin ωt
v2 (t)
= b cos ωt
2π
T :
(33)
,
da cui, integrando ẋ2 = b cos ωt si ottiene x2 = ωb sin ωt + c. Dalle condizioni iniziali deriva che
c = 0. Sostituendo il valore trovato di x2 nell’ultima delle (32) otteniamo:
ẋ3
=
b
ab
sin ωt a sin ωt =
sin2 ωt
ω
ω
(34)
.
Integrando sul periodo T la (34) otteniamo:
x3 (T ) − x3 (0)
T
ab t
1
abT
ab2π
abπ
= yd − y0 =
−
sin 2ωt =
=
= 2
2
ω 2 4ω
2ω
2ω
ω
0
,
(35)
da cui:
ab =
4π
(yd − y0 ) .
T2
(36)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
14
Applichiamo la (36) al caso in cui si voglia eseguire un parcheggio come mostrato in Fig. 3,
ovvero yd = 0, y0 > 0. Dalla figura si evince che a e b devono avere segno opposto. Inoltre si
capisce che lo stesso obiettivo (annullamento di x3 = y) può essere perseguito con una sterzata
di piccola entità (|b|) ed un grosso avanzamento (|a|) o con una sterzata più decisa ed un piccolo
avanzamento, come mostrato in Fig. 4. D’altro canto, l’esperienza nella guida di un’automobile ci
porta alle stesse conclusioni.
12
10
8
x3
b=4
b=2
6
b=3
4
2
0
0
2
4
6
8
10
x1
Figura 4: Possibili diverse manovre di parcheggio di un uniciclo.
1.3.3
Sfera che rotola su un piano
Una sfera che rotola senza strisciare su un piano, come mostrato in Fig. 5, si muove in uno spazio
di configurazione di dimensione 5. Ciò perché l’unico vincolo olonomo esistente è che la quota
del centro della sfera rispetto al piano deve essere costante e pari al raggio R. Per il resto, il punto
di contatto può muoversi liberamente sul piano (2 DOF) e l’orientazione della sfera può essere
qualsiasi (3 DOF), quindi il sistema ha 5 DOF. Fissiamo una terna O0 x0 y0 z0 solidale al piano di
z0
o1z
z1
y0
o1y
O0
x0
y1
O1
x1
o1x
Figura 5: Esempio di sistema meccanico con vincoli anolonomi: sfera che rotola su un piano.
rotolamento con l’asse z ad esso ortogonale e una terna O1 x1 y1 z1 solidale alla sfera e avente
l’origine nel centro della sfera stessa. Come coordinate lagrangiane si possono scegliere:
• le due coordinate x = o1x e y = o1y del centro della sfera rispetto al piano, essendo la quota
z = o1z costante e pari al raggio R della sfera;
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
15
• tre angoli (ad esempio i tre angoli di Eulero ZYZ φ, θ, ψ) che descrivono
l’orientazione della sfera rispetto alla terna solidale al piano.
I gradi di mobilità della sfera, vincolata a rotolare senza strisciare sul piano, sono invece:
caso 1
3 gradi di mobilità, se permettiamo alla sfera di ruotare attorno
alla normale al piano passante per il centro della sfera stessa
(ωz 6= 0);
caso 2
2 gradi di mobilità, se impediamo alla sfera di ruotare attorno
alla normale al piano (ωz = 0).
Uno spostamento virtuale della sfera ha solo tre componenti indipendenti nel caso 1 e solo due
componenti indipendenti nel caso 2. Infatti, per il vincolo di puro rotolamento, le rotazioni infinitesime della sfera parallelamente agli assi x0 e y0 sono legate agli spostamenti δx e δy del centro
della sfera dalle seguenti relazioni desumibili osservando le proiezioni ortogonali sui piani xz e yz
mostrate nelle Fig. 6 e 7:
δx =
Rωy δt
δy
−Rωx δt
=
(37)
,
dove R è il raggio della sfera. Ricordiamo che il legame esistente tra la velocità angolare ω e le
ωy δt
δx
z0
O1
o1z
o1x
x0
y0
Figura 6: Sfera che rotola su un piano: proiezione ortogonale sul piano xz.
ωx δt
z0
o1z
−δy
O1
o1y
x0
y0
Figura 7: Sfera che rotola su un piano: proiezione ortogonale sul piano yz.
derivate temporali degli angoli di Eulero è:
ω
=
=
φ̇k + θ̇j 0 + ψ̇k1 =
 




0
−sφ
sθ cφ
 0  φ̇ +  cφ  θ̇ +  sθ sφ  ψ̇
1
0
cθ
,
dove k è il versore dell’asse z0 , j 0 è il versore della linea dei nodi orientata (asse y 0 ) e k1 è il
versore dell’asse z1 . da cui:
ωx
= −sφ θ̇ + cφ sθ ψ̇
(38)
ωy
= cφ θ̇ + sφ sθ ψ̇
(39)
ωz
= φ̇ + cθ ψ̇
(40)
.
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
16
Sostituendo le Eq. (38) e (39) nelle Eq. (37) si ottiene:
δx =
Rωy δt = R(cφ δθ + sφ sθ δψ)
(41)
δy
−Rωx δt = R(sφ δθ − cφ sθ δψ) .
(42)
=
Le Eq. (41) e (42), rappresentano due vincoli anolonomi, da aggiungere al vincolo olonomo z = R,
cui è soggetto il sistema sfera-piano.
• Nel caso 1, in cui è permesso alla sfera avere una componente di velocità angolare ωz ,
si possono scegliere come variabili indipendenti del generico spostamento virtuale le tre
rotazioni infinitesime δφ, δθ, δψ.
• Nel caso 2 invece, in cui non è permesso alla sfera avere una componente di velocità
angolare ωz , esiste un ulteriore vincolo anolonomo che otteniamo eguagliando a zero il II
membro dell’Eq. (40), ovvero:
φ̇ + cθ ψ̇
=
0 .
(43)
In questo caso, come variabili indipendenti del generico spostamento virtuale, si possono
scegliere, ad esempio, due delle tre rotazioni infinitesime δφ, δθ, δψ.
Allotta, Fioravanti, Malvezzi
2
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
17
Visual servoing
Nell’ambito del “Sensor-based Control” (metodologie generali che sfruttano le informazioni provenienti da diversi sensori per misurare/osservare lo stato di un sistema dinamico da controllare)
può essere inserita una serie di tecniche di controllo nota in letteratura come “Visual Servoing,”
cioè Asservimento Visivo.
Nel Visual Servoing si sfruttano in anello chiuso le informazioni visive provenienti da una o
più telecamere (i sensori di visione) per il controllo di un generico sistema robotico: le immagini
provenienti dalle telecamere vengono elaborate in “real time” (il tempo di campionamento del
sistema di controllo dipende dall’hardware utilizzato) e le informazioni da queste ricavate vengono
utilizzate per guidare il sistema verso la configurazione desiderata.
L’asservimento Visivo quindi rappresenta un’applicazione multidisciplinare che presenta differenti aspetti:
• i sensori di visione e la formazione dell’ immagine digitale;
• l’ Image Processing (tecniche di elaborazione dell’immagine digitale);
• la Geometria Proiettiva (studio delle proprietà geometriche dell’immagine);
• le tecniche di controllo non lineare per sistemi tempo-varianti;
• i sistemi operativi real-time;
• lo studio dei sistemi meccatronici (manipolatori industriali, robot mobili ecc...).
In questo capitolo tratteremo le principali tecniche di Visual Servoing soffermandoci inizialmente sulla modellazione dei sensori di visione, sull’elaborazione delle immagini digitali e sulle
principali nozioni di geometria proiettiva utili nei sistemi di asservimento visivo.
2.1
Sistemi di visione
L’architettura di un sistema di visione può essere molto differente in base al dispositivo meccatronico da controllare, all’ hardware utilizzato ed alle prestazioni richieste. In particolare i sistemi
di visione differiscono per:
• il numero di telecamere utilizzate
– si parla di Monocular System quando si sfrutta 1 telecamera;
– Stereo System quando si utilizzano 2 telecamere;
– Multi Camera System con più di 2 telecamere;
• la posizione delle telecamere nel sistema
– si parla di sistema Eye in Hand se il sensore visuale è posto sul sistema meccatronico
mobile (end-effector di un manipolatore, sopra un robot mobile etc...);
– nei sistemi Eye to Hand le telecamere sono solidali al sistema fisso (es. la base del
manipolatore);
• architettura del sistema di controllo
– Nel Direct Visual Servoing è presente un solo controllore, il visual controller, che genera direttamente i comandi agli attuatori e, per tale motivo, lavora ad una frequenza
non inferiore a 100 Hz;
– Nell’ Indirect Visual Servoing, invece, il visual controller lavora ad una frequenza più
bassa (circa 20÷30 Hz) ed assegna il riferimento per il controllore del robot, un controllore di basso livello che opera in un ulteriore anello di retroazione interno con una
frequenza di circa 100 Hz (tale architettura è conosciuta come Dynamic Look and
Move)
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
18
• conoscenza della scena osservata dalle telecamere
– Nel Model Based Approach si conoscono alcune caratteristiche geometriche della
scena osservata
– Nel Model Free Approach non si conoscono le dimensioni e le caratteristiche della
scena osservata
• conoscenza dei modelli matematici di telecamere utilizzate
– Nel Calibrated Approach si ha una conoscenza precisa dei parametri geometrici delle
telecamere utilizzate;
– Nell’Uncalibrated Approach non si ha una conoscenza precisa delle telecamere utilizzate.
2.1.1
L’immagine digitale
Nel Visual Servoing lo stato attuale del sistema controllato è ottenuto (anche) da informazioni
provenienti da immagini. Soffermiamoci quindi sul concetto di immagine:
Un’immagine è una quantità fisica misurata dipendente (anche) dalla posa relativa tra il
sensore e la scena osservata.
Ad esempio, l’intensità di radiazione luminosa (raggi X, raggi UV, spettro visibile e infrarosso),
un’onda acustica, l’altezza di una superficie misurate in funzione dello spazio costituiscono delle
immagini. In generale le immagini possono essere:
• Bidimensionali se rappresentano la proiezione di una scena 3D su un sensore 2D, in questo
caso avremo che un immagine può essere rappresentata da una funzione I = f (x, y);
• Tridimensionali si rappresenta una grandezza fisica misurata su uno spazio 3D (come nella
Tomografia Assiale Computerizzata)
Per trattare le immagini con un calcolatore occorre convertire l’immagine continua nella sua rappresentazione digitale: le telecamere digitali svolgono questa funzione e portano alla formazione
dell’ immagine digitale. Qualunque tipo di immagine digitale, indipendentemente dalla codifica, è
rappresentabile mediante una matrice di numeri.
• L’esatta corrispondenza fra un’immagine digitale ed il mondo fisico è determinata dal processo di acquisizione dipendente dal sensore utilizzato;
• ogni informazione reperibile nelle immagini (forme, distanze, identificazione di oggetti) deve
essere, quindi, estratta da una matrice di numeri che la codificano.
L’informazione numerica viene associata ad ogni elemento della matrice rappresentante l’immagine:
tale elemento viene chiamato pixel (ovvero PICture ELement). Maggiore è la grandezza della
matrice, maggiore risulterà la risoluzione dell’immagine (corrispondente al numero di pixel) ed
anche lo spazio di memoria richiesto per la memorizzazione nell’elaboratore. Le immagini digitali possono essere suddivise in due famiglie principali: le immagini mono-canale e le immagini
multi-canale. Tra le immagini monocanale ricordiamo:
• le immagini binarie, in cui ogni pixel può assumere soltanto 2 valori discreti, il valore 0
(corrispondente al bianco) ed il valore 1 corrispondente al nero; questo tipo di immagine
per essere memorizzata nel calcolatore necessita, quindi, di uno spazio di memoria assai
limitato;
• le immagini in scala di grigi, in cui ogni pixel può assumere un valore intero tra 1 e 256
compresi (ovvero tra 0 e 255 in convenzione zero-based), dove 1 corrisponde al nero, 256
corrisponde al bianco e gli altri valori corrispondono a dei grigi scalati dal nero al bianco; in
questo tipo di immagine ogni pixel occupa, quindi, 8 bit di memoria (normalmente si parla di
“immagini ad 8 bit”);
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
19
• le immagini indicizzate in cui ogni pixel assume un valore intero che viene poi codificato in
un ben determinato colore secondo una predefinita mappa colorimetrica (la colormap); il
valore del del pixel non è altro che il puntatore ad un determinato elemento della mappa
colorimetrica.
Nelle immagini multicanale, invece, ad ogni pixel è associato più di un valore: tali immagini possono quindi esser codificate da una matrice tridimensionale in cui gli indici di riga e di colonna
identificano la posizione del pixel mentre la terza dimensione della matrice (pari a 3 o più) identifica il numero di canali dell’immagine (la profondità dell’immagine). Fra le immagini multicanale
ricordiamo:
• le immagini truecolor a 3 canali (rappresentabili quindi da una matrice di dimensione I m×n×3 ),
nelle quali ad ogni pixel vengono associati 3 numeri interi compresi tra 1 e 256 solitamente
corrispondenti ai livelli dei tre colori primari rosso, verde e blu da cui deriva il nome di RGB
images oppure immagini a 32 bit (normalmente oltre alla colorimetria RGB possono essere
utilizzate anche altre codifiche come ad esempio l’HSV (Hue Saturation Value, letteralmente
“colore,” “saturazione,” “intensità” );
• le immagini multispettrali, in cui ad ogni pixel sono associati più di 3 canali (ed esempio, i
3 livelli di colore più 2 canali corrispondenti ai livelli di radiazione infrarossa e ultravioletta
emesse dalla scena); questo tipo di immagine viene utilizzato per lo più in fotogrammetria e
nelle foto satellitari (in cui si può arrivare anche ad una quindicina di canali per pixel!).
Solitamente un file rappresentante un’immagine digitale contiene oltre alla codifica dei pixel anche
alcune stringhe contenente informazioni aggiuntive di carattere generale chiamate “metadata”
(come ad esempio l’ora e la data della generazione dell’immagine ecc..).
2.2
Telecamere
Prima di addentrarci nello studio del Visual Servoing, occorre analizzare più da vicino un componente fondamentale del sistema di visione che rende possibile la realizzazione della retroazione
nello schema di controllo attraverso la generazione del flusso video (video-stream): il sistema di
acquisizione dell’immagine. Tale complesso sistema comprende diversi dispositivi ma, didatticamente, può esser pensato come costitutito da tre componenti principali:
• una lente (ottica della telecamera) utilizzata per far convergere la radiazione luminosa proveniente dalla scena esterna in un unico punto detto fuoco;
• un sensore, interposto fra il fuoco e la lente, che misura l’intensità luminosa e la codifica in
informazione digitale (attraverso anche un processore numerico);
• un dispositivo elettronico, chiamato frame-grabber, per prelevare l’immagine digitale con la
cadenza temporale desiderata ed inviarla al calcolatore (in passato tale dispositivo non era
integrato nella telecamera).
Per misurare l’intensità luminosa riflessa dalla scena il sensore impiega un elemento fotosensibile, detto pixel del sensore, in grado di trasformare l’energia luminosa in energia elettrica. Sono
disponibili sul mercato diversi tipi di sensore che differiscono per il principio fisico utilizzato per realizzare la trasformazione dell’energia. I dispositivi più diffusi sono i sensori CCD e CMOS basati
sull’effetto fotoelettrico dei semiconduttori.
Risulta chiaro che sensori CCD e CMOS assumono il ruolo di campionatori spaziali della radiazione luminosa (spazialmente continua) focalizzata dalla lente, mentre il frame-grabber ha invece il ruolo di campionatore temporale.
CCD
Un sensore CCD (Charge Coupled Device) è costituito da una matrice rettangolare di pixel. In
virtù dell’effetto fotoelettrico, quando un fotone colpisce la superficie del semiconduttore, si crea
un numero di elettroni liberi, in maniera tale che ogni elemento accumula una carica dipendente
dall’integrale temporale dell’illuminazione incidente sull’elemento fotosensibile. Tale carica viene
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
20
passata da un meccanismo di trasporto (simile a un registro a scorrimento analogico) ad un
amplificatore in uscita, mentre allo stesso tempo si scarica il pixel. Il segnale elettrico viene ulteriormente elaborato in maniera tale da produrre il segnale video reale. I sensori CCD sono quelli
che comunemente vengono utilizzati sia nelle fotocamere che nelle videocamere digitali e nelle
webcam a basso costo. Un effetto indesiderato comune a tutti i CCD è il fenomeno del “trabocco”:
questo si verifica quando un pixel saturato di elettroni trabocca influenzando i pixel vicini e quindi
alterando il loro livello di colore; tale effetto provoca una sfocatura dell’immagine nella zona di
trabocco.
CMOS
Un sensore CMOS Complementary Metal Oxide Semiconductor è costituito da una matrice rettangolare di fotodiodi. La giunzione di ciascun fotodiodo è precaricata e viene scaricata allorché
viene colpita da fotoni. Un amplificatore integrato in ciascun pixel è in grado di trasformare questa
carica in un livello di tensione o di corrente. La principale differenza con il sensore CCD è che i
pixel di un sensore CMOS sono dispositivi che non integrano; dopo esser stati attivati misurano
la quantità non il volume. In questo modo si impedisce l’effetto del trabocco dell’immagine che
talvolta può verificarsi nei sensori CCD. Tali dispositivi sono in generale più costosi e meno diffusi
dei CCD.
Alcuni parametri caratteristici dei sensori CCD sono:
• Dimensione del sensore (i più comuni dell’ordine di grandezza del cm);
• Grandezza della matrice di pixel, solitamente 320 × 240 o 640 × 480 per webcam più economiche; per le fotocamere digitali si arriva a risoluzioni anche molto maggiori (ad es.
2048 × 1536);
• dimensione del pixel, dell’ordine di una decina di µm;
• capacità del pixel, circa 700 elettroni per µm2 ;
• fattore di riempimento del sensore (il rapporto percentuale tra l’area sensibile del sensore e
l’area totale);
• la sensibilità assoluta intesa come il segnale più debole identificabile (più alto del rumore di
lettura);
• la sensibilità relativa intesa come il più piccolo incremento di segnale rilevabile
Come ogni sensore fisico anche il CCD è soggetto a diversi tipi di rumore che altera la completa
ripetitibilità e degrada la qualità della misura; i tipi di rumore più importanti sono:
• il Photon Noise, rumore dovuto alla natura particellare della luce; tale rumore è sempre
presente e può essere caratterizzato come un rumore bianco;
• il rumore di quantizzazione dovuto alla conversione discreta del segnale continuo;
• il rumore termico, una componente di rumore che incrementa all’aumentare della temperatura di lavoro del CCD.
2.2.1
Modelli di telecamera
Mentre in precedenza abbiamo illustrato le caratteristiche tecniche, analizziamo adesso la telecamera digitale dal punto di vista geometrico: una telecamera infatti non realizza altro che una
mappatura fra lo spazio Euclideo 3D e un’immagine appartenente a un piano 2D. Caratterizzare
con precisione la mappatura della telecamera è di vitale importanza nell’asservimento visivo: le
informazioni visive infatti sono, solitamente, le uniche che vengono sfruttate per il controllo del
sistema.
Ricordiamo d’altra parte che, come in ogni altro tipo di sensore (termometro, encoder ecc..), anche nella telecamera occorre da definire una funzione caratteristica, cioè una relazione che lega
gli input da misurare agli output del sensore. Per la definizione di tale funzione occorre innanzi
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
21
tutto definire il tipo di funzione caratteristica (ad esempio nel termometro si ipotizza una funzione
lineare tra innalzamento del livello di mercurio e la misura di temperatura) e poi stabilire fissare i
parametri che definiscono questa funzione.
In questo paragrafo verranno quindi illustrati:
• quali modelli geometrici possono essere utilizzati per mettere in relazione una scena 3D
(input da misurare) alla sua proiezione in una immagine 2D (output del sensore di misura)
che meglio si adattano ad una telecamera digitale (definizione della funzione caratteristica);
• come determinare i parametri specifici del modello utilizzato per una data telecamera, ovvero
come effettuare la calibrazione della telecamera.
2.2.1.1 Il modello Basic Pinhole Il modello di telecamera più semplice che possa essere
utilizzato è il cosiddetto modello “pinhole” (vedi Fig.8) in cui l’immagine di un punto P dello spazio
si forma su uno schermo grazie ad un unico raggio luminoso che attraversa un forellino di piccole
dimensioni. Seguendo le leggi dell’ottica geometrica, tanto più piccolo è il foro, tanto migliore
l’approssimazione che “un unico raggio” proveniente da un punto P formi un’immagine puntuale
p sullo schermo. In realtà ci sono due problemi: il primo problema è che tanto più piccolo è il foro,
tanta meno energia luminosa è disponibile per formare l’immagine (ad esempio su una pellicola
fotosensibile o su un CCD); il secondo problema è che con un foro di dimensioni molto piccole le
leggi dell’ottica geometrica non sono più idonee a rappresentare il fenomeno in quanto gli effetti
di bordo diventano apprezzabili.
schermo
f
P
punto
principale
asse ottico
p
foro
Figura 8: Modello “pinhole”. L’immagine di un punto sullo schermo si forma lasciando passare un
unico raggio attraverso un foro di piccole dimensioni praticato in uno schermo opaco che funge
da otturatore.
Sebbene le telecamere siano, ovviamente, costruite in modo molto diverso dallo schema di
principio di Fig.8, tale modello è una buona approssimazione del funzionamento reale della telecamera. Si faccia riferimento alla Fig.9 (modello di proiezione prospettica centrale), del tutto
equivalente a quello mostrato in Fig.8. Sono distinguibili:
oc
centro della telecamera
zc
asse ottico della telecamera
f
distanza focale
oc xc yc zc
oim
oim xim yim
sistema di riferimento solidale alla telecamera
origine del sistema di riferimento dell’immagine
sistema di riferimento dell’immagine
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
22
Modellistica e Controllo di Sistemi Meccanici
f
oim
origine immagine
x0
y0
centro camera
piano immagine
punto principale
zc
oc
xim
xc
yc
p
yim
P
Figura 9: Modello “pinhole” di telecamera. Sono evidenziati: il centro della telecamera oc , il punto
principale dell’immagine, l’origine del sistema di riferimento immagine
x0 y0
coordinate del punto principale nel piano immagine
P = [Xc Yc Zc ]T
coordinate di un punto dello spazio rispetto al sistema di riferimento oc xc yc zc
p = [x y]T
coordinate della proiezione di P sul piano immagine rispetto
al sistema immagine
Le coordinate del punto immagine p nel sistema di riferimento oc xc yc zc sono funzione delle coordinate del punto proiettato come segue:
Xc
Zc
Yc
= f .
Zc
x =
y
f
(44)
Le (44) esprimono il meccanismo della “proiezione prospettica”; si noti come queste relazioni
non siano lineari rispetto alle coordinate 3D del punto. Se per rappresentare i punti immagine
si utilizzano le coordinate omogenee (si veda a tal riguardo il prossimo paragrafo), la (44) può
essere trasformata in una mappatura lineare tra ogni punto dello spazio 3D (SE(3), ovvero Special
Euclidean di dimensione 3) e la sua proiezione sul piano immagine, infatti si ha:

 


f Xc
f 0 0
Xc
p̃ =  f Xc  =  0 f 0   Yc  = KP
(45)
Zc
0 0 1
Zc
La matrice K è chiamata matrice di calibrazione della telecamera e contiene tutti i parametri che
servono a caratterizzare la mappatura del sensore: i parametri interni a tale matrice sono detti
parametri intrinseci, caratteristici della singola telecamera.
Si può notare come l’unico parametro che caratterizza il modello di mappatura del sensore sia in
questo caso la distanza focale f (f variabile consente di modellare una telecamera con zoom).
Nel caso particolare in cui si abbia f = 1 K si semplifica nella matrice identica I e la (45) diventa:




mx
Xc /Zc
m̃ =  my  =  Yc /ZC  = (1/Zc ) P
(46)
1
1
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
23
dove m̃ è chiamato punto immagine normalizzato.
NOTE:
• Due punti sullo stesso raggio hanno stessa proiezione sull’immagine: si ha quindi una
perdita d’informazione nella mappatura di una scena 3D sull’immagine (ad esempio le profondità dei punti); la mappatura non è biunivoca;
• Se si conosce K è possibile risalire agli angoli reali che i raggi proiettivi formano tra di loro:
è possibile cioè misurare le direzioni reali degli oggetti visualizzati nell’immagine rispetto
alla direzione dell’asse focale. Si ha infatti
d = K −1 p̃
,
(47)
dove d è, infatti, un vettore nella direzione del punto P ;
• Il punto principale, il piano immagine, ed il sistema solidale alla telecamera sono entità ideali
che modellano il comportamento delle telecamere reali con una certa approssimazione (non
sono quindi direttamente identificabili sulla telecamera reale).
2.2.1.2 il modello CCD Full Perspective Sfortunatamente il modello basic pinhole non è sufficiente (data la sua eccessiva semplicità) a riprodurre con fedeltà la mappatura operata dalle
telecamere digitali reali. Per raffinare il modello di telecamera si introducono, quindi, dei parametri
che caratterizzano il sensore CCD considerando che l’immagine digitale risulta definita in pixel:
nasce quindi il modello CCD Full Perspective. Considerando questi cambiamenti la (45) può
essere rielaborata nella seguente:





f αx f αx cot(φ) x0
αx αx cot(φ) x0
f 0 0
f αy / sin(φ) y0  P =  0 αy / sin(φ) y0   0 f 0  = KP ;
(48)
p̃ =  0
0 0 1
0
0
1
0
0
1
La matrice di calibrazione K questa volta contiene i seguenti parametri:
• αx ed αy rappresentano i pixel per unità di lunghezza rispettivamente lungo la direzione di
base e lungo l’altezza del pixel stesso: solitamente il pixel non è perfettamente quadrato,
quindi questi due valori non coincidono;
• φ è l’angolo formato dai lati del pixel (di solito molto vicino a π/2);
• x0 ed y0 rappresentano le coordinate in pixel del punto principale: mentre nel modello basic pinhole l’origine del sistema immagine coincideva col punto principale, in questo caso
l’origine coincide con l’angolo in alto a sinistra (top-left coordinate system);
NOTE:
• il termine αx cot(φ) è detto skew parameter ed risulta 0 se la base e l’altezza del pixel sono
ortogonali
• il parametro αy /αx è noto come aspect ratio della telecamera;
• il modello full perspective CCD è ancora lineare se si rappresentano i punti immagine in
coordinate omogenee e la conoscenza di K implica ancora quella delle direzioni dei raggi
nello spazio 3D come in precedenza;
• questo modello ben si adatta alle telecamere professionali di alto costo;
Utilizzando questo modello un punto immagine in pixel di coordinate u e v è espresso dalla
seguente relazione:


u
p =  v  = K m̃ .
(49)
1
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
oim
Modellistica e Controllo di Sistemi Meccanici
24
xim
x0
φ
x'c
o'c
y0
θ
p
yim
y'c
Figura 10: Piano immagine con non perfetta ortogonalità tra gli assi xim e yim
2.2.1.3 il modello CCD con distorsione della lente Sebbene il modello CCD Full Perspective risulti nella maggior parte dei casi sufficientemente raffinato, occorre adottare un modello di
telecamera ancora più complesso quando si voglia utilizzare nell’asservimento visivo telecamere
a basso costo come le webcam commerciali. In questi casi infatti l’angolo di apertura della telecamera risulta molto ampio e l’ottica è di scarso pregio: la lente quindi produce una distorsione
nella mappatura sull’immagine che non può essere trascurata. Sia m un generico punto immagine normalizzato di coordinate mx = Xc /Zc ed my = Yc /Zc ; per inserire l’effetto distorcente della
lente occorre a questo punto considerare la seguente relazione:
md x
md =
= dr m + dt ,
(50)
md y
dove
• md rappresenta il punto distorto a causa della lente;
• dr è il fattore di distorsione radiale (scalare);
• dt è il vettore di distorsione tangenziale.
Un modello abbastanza accurato di distorsione calcola dr e dt nel seguente modo:
dr
dt
(1 + c1 r2 + c2 r4 + c3 r6 )
2p1 mx my + p2 (r2 + 2m2x )
=
2p2 mx my + p1 (r2 + 2m2y )
=
.
(51)
Nella (51), r2 = m2x + m2y rappresenta la distanza del punto normalizzato dal punto principale, i
coefficienti c1 ,c2 , c3 e p1 p2 vengono chiamati rispettivamente coefficienti di distorsione radiale e
tangenziale. A questo punto, per trovare il punto immagine espresso in pixel si applica al punto
distorto dalla lente md il modello CCD full projective (si considera l’effetto proiettivo del sensore):





u
f αx f αx cot(φ) x0
md x
f αy / sin(φ) y0   md y  .
(52)
p= v  = 0
1
0
0
1
1
NOTE :
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
y0
x0
αx xc
αy yc tan θ
xim
oim
φ
xc
o'c
αy y c
c os θ
25
Modellistica e Controllo di Sistemi Meccanici
/
θ
p
yim
yc
Figura 11: Piano immagine con non perfetta ortogonalità tra gli assi xim e yim
• il modello CCD con distorsione della lente non è esprimibile linearmente neppure se si
utilizzano per i punti immagine le coordinate omogenee;
• tale modello è funzione di 10 parametri (5 della matrice di calibrazione K più 5 parametri
di distorsione): il processo di calibrazione che porta alla determinazione di tali parametri
risulta quindi abbastanza complesso;
• l’effetto della distorsione, se molto accentuato come nelle webcam, provoca nell’immagine
il fenomeno del barrelling (le rette dello spazio 3D non si mappano più in rette ma in curve);
• per riuscire a calcolare correttamente le direzioni dei raggi proiettanti occorre eseguire
l’undistortion dell’immagine (si raddrizza l’immagine) conoscendo i coefficienti di distorsione
e K.
Il processo di undstortion è stato ideato da ricercatori dell’Università finlandese di Oulu e consiste
in un processo iterativo che consente di invertire la relazione (50) con il modello di distorsione
della (51); dall’immagine acquisita dalla telecamera sono noti i punti in pixel: occorre invece
calcolare i punti normalizzati corrispondenti per individuare le reali direzioni dei raggi proiettanti:
nota K i punti distorti possono essere individuati semplicemente invertendo la (52):




md x
u
 md y  = K −1 p =  v  .
1
1
Una volta noti i punti distorti, la procedura iterativa sfrutta l’espressione inversa della (??):
m=
md − dt (m)
dr (m)
,
(53)
dove si valuta il punto normalizzato di primo tentativo m1 calcolando d1t (md) e d1r (md ) e si raffina
poi la stima ai passi successivi calcolando al passo i-esimo mi con d1t (mi−1 ) e d1r (mi−1 ) (noti i
coefficienti di distorsione). Solitamente il ciclo converge dopo una decina di iterazioni.
2.2.2
Calibrazione della telecamera
Una volta scelto il modello di telecamera da utilizzare per misurare correttamente le direzioni e
le relazioni geometriche degli oggetti della scena proiettati sull’immagine, occorre individuare il
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
26
Modellistica e Controllo di Sistemi Meccanici
valore dei parametri che lo caratterizzano (un parametro per il basic pinhole, 5 per la CCD full
projective camera e 10 per la CCD con distorsione delle lenti): i metodi che vengono utilizzati a
tale scopo sono chiamati tecniche di calibrazione. La calibrazione deve essere effettuata su ogni
telecamera e può essere effettuata sulla stessa telecamera anche più di una volta a distanza di
tempo. I metodi classici di calibrazione sfruttano le corrispondenze tra i punti di una griglia (griglia
di calibrazione) di coordinate note e la loro proiezione nell’immagine. Consideriamo allora il punto
3D P o della griglia di calibrazione solidale ad al sistema di riferimento < o >: vediamo qual’è la
relazione fra questo punto e la sua immagine in pixel. Dapprima abbiamo:
c c Ro to
P
Po
=
(54)
1
1
0T 1
dove Roc e tco sono rispettivamente la matrice di rotazione ed il vettore traslazione della terna
griglia < o > rispetto alla terna telecamera < c >. Sfruttando il modello adottando il modello di
telecamera CCD full projective e la (54) abbiamo:
c c Po
Po
p = K Ro to
=Q
,
(55)
1
1
dove p è la proiezione
del punto P o sull’immagine espressa in coordinate omogenee. La matrice
Q = K Roc tco ∈ R3×4 è chiamata camera projection
matrix: mentre gli elementi di K sono
dei parametri intrinseci della telecamera gli elementi di Roc tco sono detti parametri estrinseci.
Nota Q è possibile ricavare la matrice di calibrazione K effettuando la decomposizione RQ del
suo primo blocco 3 × 3 KRoc . Il problema della calibrazione si è a questo punto ridotto al calcolo
della matrice Q. Da una generica corrispondenza punto immagine - punto della griglia, sfruttando
la (55) possiamo scrivere:


 

 Xo
q11 q12 q13 q14
u


 v  =  q21 q22 q23 q24   Yo  ,
(56)
 Zo 
1
q31 q32 q33 q34
1
dove con qij si sono indicati gli elementi di Q, le incognite del nostro problema. Sviluppando la
(56) abbiamo:
u = q11 Xo + q12 Yo + q13 Zo + q14
(57)
v = q21 Xo + q22 Yo + q23 Zo + q24
(58)
1 = q31 Xo + q32 Yo + q33 Zo + q34
(59)
.
Inserendo dapprima il secondo membro della (59) (uguale ad 1) nei primi membri delle (57) e
(58), riorganizzando infine dli elementi qi j in un vettore delle incognite a, si ricava il seguente
sistema lineare omogeneo:
L1 a = 0 ,
(60)
dove la matrice dei coefficienti (nota)L1 ∈ R2×12 ed il vettore a ∈ R12×1 risultano rispettivamente:
Xo Yo Zo 1 0
0
0 0 −uXo −uYo −uZo −u
L1 =
,
(61)
0
0
0 0 Xo Yo Zo 1 −vXo −vYo −vZo −v
e
a=
q11
q12
q13
q14
q21
q22
q23
q24
q31
q32
q33
q34
T
.
(62)
Da una corrispondenza punto griglia-punto immagine è dunque possibile ricavare 2 vincoli lineari
nelle incognite qij . Se si ha una griglia con N punti con N >> 6, sfruttando tutte le corrispondenze, si arriva ad un sistema lineare omogeneo sovradimensionato del tipo:
La = 0
,
(63)
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
con
Modellistica e Controllo di Sistemi Meccanici



L=


L1
...
Li
...
LN
27






,
(64)
matrice di coefficienti noti ∈ R2N ×12 . Il sistema lineare (63) non è risolvibile esattamente ma
si può trovare la soluzione ottima nel senso dei minimi quadrati: si cerca cioè la soluzione che
minimizza ||La|| imponendo ||a|| = 1. La soluzione a questo problema si calcola effettuando la
decomposizione in valori singolari (SVD)della matrice L che calcola le matrici U D e V tali che:
L = U DV T
(65)
con D matrice dei valori singolari diagonale ∈ R12×12 (ordinati sulla diagonale in valore decrescente), U matrice con colonne ortogonali ∈ R2N ×12 e V matrice ortogonale ∈ R12×12 . La
soluzione al problema è rappresentata dall’ultima colonna della matrice V (corrispondente al minimo valore singolare). Una volta trovato il vettore a ottimo (ovvero la matrice Q), è possibile
ricavare la matrice di calibrazione K come spiegato sopra .
NOTE:
• Una volta trovati gli elementi K non interessa la stima dei singoli parametri che li compongono (ad esempio αx , αy f etc..): la loro stima infatti risulterebbe molto difficoltosa e di
scarsa utilità dato che nota K è noto in questo caso tutto ciò che occorre per misurare con
accuratezza le direzioni dei raggi proiettivi;
• il metodo di calibrazione proposto è uno dei più semplici: solitamente questa tecnica lineare
viene utilizzata per inizializzare un ulteriore raffinamento degli elementi di K che sfrutta la
minimizzazione di un funzionale non-lineare (metodi di ottimizzazione Gauss-Newton o LM);
• ci sono anche tecniche di calibrazione che non prevedono l’utilizzo di alcuna griglia ma
sfruttano alcune proprietà geometriche note della scena visualizzata dalla telecamera: si
parla in tal caso di metodi di autocalibrazione. L’autocalibrazione solitamente non raggiunge
la precisione delle tecniche di calibrazione con griglia ma può essere effettuata on-line prima
di eseguire il controllo del sistema asservito.
2.3
Elementi di Geometria Proiettiva
Affrontiamo adesso lo studio dell’immagine da un punto di vista geometrico: ci interessano in
particolare dei modelli geometrici adeguati per interpretare le relazioni esistenti tra due o più
immagini differenti di uno stesso scena e caratterizzare gli elementi geometrici salienti all’interno
dell’immagine.
2.3.1
Primitive nello spazio immagine
Come sicuramente noto, un punto nel piano può essere rappresentato dalla coppia di coordinate
(x, y) ∈ R2 : di conseguenza è immediato identificare il piano con R2 . Considerando R2 come uno
spazio vettoriale, la coppia di coordinate (x, y) rappresenta un vettore; quindi identifichiamo un
punto come un vettore. Introduciamo adesso la notazione omogenea per i punti e le linee di un
piano.
2.3.1.1 Rappresentazione omogenea delle rette
da un’ equazione del tipo:
ax + by + c = 0
Una linea retta nel piano è rappresentata
;
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
28
differenti scelte dei parametri a, b e c generano rette differenti. Una retta quindi può essere naturalmente rappresentata dal vettore (a, b, c)T . La corrispondenza tra rette e vettori però non è biunivoca poiché i vettori differenti (a, b, c)T e (ka, kb, kc)T con k 6= 0 rappresentano la stessa linea.
Per tale ragione in geometria proiettiva due vettori che differiscano per una scalatura costante
sono considerati equivalenti. Da tale affermazione si deduce che ogni vettore (a, b, c)T è rappresentativo di una ben precisa classe di equivalenza.
Il set delle classi di equivalenza dei vettori ∈ R3 − (0, 0, 0)T forma lo spazio proiettivo P2 .
2.3.1.2 Rappresentazione omogenea dei punti Un punto x = (x, y)T giace sulla retta l =
(a, b, c)T se e solo se ax + by + c = 0. Questa relazione può essere scritta con un prodotto scalare
del tipo (x, y, 1) • (a, b, c)T = (x, y, 1) • l = 0, dove il punto (x, y)T in R2 è rappresentato come un
vettore ∈ R3×1 con l’aggiunta della terza coordinata pari ad 1. Osserviamo che per ogni costante
k 6= 0 abbiamo (kx, ky, k) • l = 0 se e solo se (x, y, 1) • l = 0. Risulta allora naturale considerare
il set di vettori (kx, ky, k) come una rappresentazione del punto (x, y)T in R2 . Cosı̀ come le linee,
i punti sono rappresentati da vettori omogenei.
Un generico vettore omogeneo della forma x = (x1 , x2 , x3 )T rappresenta il punto x = (x1 /x3 , x2 /x3 )T
in R2 .
Anche i punti quindi, come vettori omogenei sono elementi di P2 .
Risultato 1 Il punto x appartiene alla retta l se e solo se xT l = 0.
2.3.1.3 Gradi di libertà Come sappiamo, per specificare un punto occorrono due valori, cioè
le due coordinate x e y. Analogamente una retta è definita da due parametri indipendenti (i due
rapporti indipendenti a : b : c) e quindi ha due gradi di libertà.
0
0
0
0
2.3.1.4 Intersezione tra rette Date due rette l = (a, b, c)T ed l = (a , b , c )T vogliamo trovare
0
la loro intersezione. Definiamo il vettore x = l ∧ l dove si indica con ∧ il prodotto vettoriale. Dalle
0
0
0
proprietà dei vettori (doppio prodotto misto) abbiamo: lT (l∧l ) = (l ∧l) = 0; cioè lT x = l T x = 0.
0
Se quindi x è la rappresentazione di un punto, x giace sia su l che su l e perciò è l’intersezione
delle due rette. Allora:
0
0
Risultato 2 L’intersezione di due rette l e l è il punto x = l ∧ l .
Con una dimostrazione analoga si ricava che:
0
0
Risultato 3 La linea passante per due punti x e x è l = x ∧ x .
2.3.1.5 Intersezione tre due rette parallele Consideriamo due rette parallele ax + by + c = 0
0
0
0
e ax + by + c = 0. Queste sono rappresentate dai vettori l = (a, b, c)T e l = (a, b, c )T che hanno
0
0
le prime due coordinate uguali. L’intersezione tra le rette risulta il punto l ∧ l = (c − c)(b, −a, 0)T
0
che, ignorando il fattore di scala (c − c), diventa (b, −a, 0)T .
Adesso, se cerchiamo di ricavare la rappresentazione non omogenea del punto, otteniamo (b/0, −a/0, 0)T .
Tale scrittura non ha senso ma ci indica che il punto ha coordinate infinitamente grandi in R2 . In
generale i punti con coordinate omogenee (x, y, 0)T non corrispondono ad alcun punto finito di R2 .
Questa osservazione concorda con l’idea comune che due rette parallele si incontrano all’infinito.
Più rigorosamente si dice che due rette parallele si incontrano in un punto e si incontrano in un
punto improprio o punto ideale.
2.3.1.6 Punti ideali e linea all’infinito I vettori omogenei (x1 , x2 , x3 )T tali che x3 6= 0 corrispondono a punti finiti di R2 . Possiamo ampliare R2 aggiungendo i punti con coordinata x3 = 0.
Lo spazio risultante comprende il set di tutti i vettori omogenei cioè lo spazio proiettivo P2 .
I punti con x3 = 0 sono detti punti ideali o punti all’infinito o punti impropri. Tale set di punti può
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
29
essere scritto (x1 , x2 , 0)T con ogni punto particolare specificato dal rapporto x1 : x2 . Possiamo
notare che questo set di punti appartiene tutto ad una singola retta, la retta all’infinito rappresentata dal vettore l∞ = (0, 0, 1)T . Utilizzando il risultato (2), si trova che la linea l = (a, b, c)T
0
0
interseca l∞ nel punto ideale (b, −a, 0)T , poiché (b, −a, 0)l = 0. Una retta l = (a, b, c )T , parallela
ad l, interseca l∞ nello stesso punto ideale (b, −a, 0)T . In notazione non omogenea (b, −a)T è un
vettore tangente alla retta ed ortogonale alla normale alla retta (a, b)T , e quindi ne rappresenta la
direzione.
Col variare della direzione della retta, varia il punto ideale (b, −a, 0)T su l∞ ; per questi motivi la
linea all’infinito può essere pensata come il set delle direzioni delle rette nel piano. Nel piano
proiettivo P2 possiamo trattare punti al finito ed all’infinito nello stesso modo; questo non è vero
nella geometria euclidea standard di R2 .
Lo studio della geometria di P2 è conosciuto come geometria proiettiva.
2.3.1.7 Un modello per il piano proiettivo Un modo utile di immaginare il piano proiettivo
P2 , è di pensarlo come un set di raggi in R3 (vedi Figura 12). Il set di tutti i vettori k(x1 , x2 , x3 )T
al variare di k, forma un raggio attraverso l’origine. Questo raggio può essere pensato come la
rappresentazione di un singolo punto di P2 . In questo modello ogni retta in P2 corrisponde ad un
piano passante per l’origine in R3 .
Si può infatti verificare che due raggi non coincidenti appartengono ad un solo piano e due piani
distinti si intersecano in un solo raggio (questo è il caso analogo di due punti distinti che definiscono una sola retta e di due rette che si intersecano in un solo punto).
Punti e rette possono essere ottenuti intersecando tale set di raggi e piani col piano x3 = 1. I
raggi che rappresentano i punti ideali sono paralleli a tale piano.
punto
ideale
l
O
x
Figura 12: Un modello del piano proiettivo.
2.3.2
Geometria di due immagini
Vediamo come è possibile modellare le relazioni che intercorrere fra due immagini. Fra tutte
le possibili trasformazioni dell’immagine le più utili in Geometria Proiettiva sono le omografie (o
trasformazioni proiettive): noi ci limiteremo a trattare questo tipo di mappature.
2.3.2.1 Omografie Una omografia o trasformazione proiettiva è una mappatura invertibile da
punti in P2 a punti in P2 che mappa linee in linee. Più precisamente:
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
30
Definizione 1 Una proiettività è una mappatura invertibile h da P2 in se stesso, tale che tre punti
x1 , x2 e x3 appartengono alla stessa retta se e solo se h(x1 ), h(x2 ) e h(x3 ) giacciono su una
stessa retta.
Le proiettività formano un gruppo, poiché sia l’inversa che la composizione di due proiettività
risultano essere sempre una proiettività. Sinonimo di proiettività è anche omografia.
Diamo adesso una equivalente definizione algebrica di proiettività:
Teorema 1 Una mappatura h : P2 → P2 è una proiettività se e solo se esiste una matrice 3 × 3
non singolare H tale che per ogni punto in P2 rappresentato dal vettore x vale: h(x) = Hx.
H dunque costituisce una mappatura lineare tra coordinate omogenee rappresentanti punti in P2 .
Non daremo qui dimostrazione del teorema, per maggiori chiarimenti si consulti [?].
Una trasformazione planare proiettiva è dunque una trasformazione lineare di vettori omogenei
3 × 1 rappresentati da una matrice non singolare 3 × 1:
 0  


x1
h11 h12 h13
x1
0
 x2  =  h21 h22 h23   x2  ,
(66)
0
h31 h32 h33
x3
x3
0
o più brevemente: x = Hx.
Notiamo che moltiplicare la matrice H per un fattore di scala arbitrario λ 6= 0 non altera la trasformazione proiettiva. Diciamo allora che H è una matrice omogenea poiché, in una rappresentazione omogenea di punti, è importante solo il rapporto tra gli elementi della matrice.
Ci sono 8 rapporti indipendenti tra gli elementi di H, quindi un’omografia ha in generale 8 gradi di
libertà.
Un’omografia allora proietta ogni figura in un’altra proiettivamente equivalente; nel modello dei
raggi, precedentemente proposto per P2 , una trasformazione proiettiva è semplicemente una
trasformazione lineare di R3 . Degli esempi di omografia sono mostrati nelle Figure 13 e 14.
x’
x
y’
y
x’
O
x
Figura 13: La proiezione centrale da un punto dello spazio O mappa punti e linee del piano π in punti e
0
linee del piano π ; quindi è una proiettività .
2.3.2.2 Trasformazioni di rette Si può dimostrare che, se un punto x appartiene ad una retta
0
l, allora il trasformato del punto x = Hx appartiene alla linea:
0
l = (H −1 )T l
.
(67)
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
immagine 1
Modellistica e Controllo di Sistemi Meccanici
R,T
31
immagine 2
piano oggetto
Figura 14: La relazione tra le immagini generate da due proiezioni centrali differenti di una stessa figura
planare è ancora una proiettività .
2.3.3
Classificazione delle proiettività
Vediamo adesso come possono essere classificate le proiettività. Le matrici di omografia nel
piano proiettivo P2 costituiscono un gruppo chiamato PL(3).
Questo gruppo contiene alcuni sottogruppi molto importanti:
• il gruppo affine, costituito dalle omografie che hanno come ultima riga (0, 0, 1);
• il gruppo Euclideo, che è a sua volta un sottogruppo del gruppo affine, avente in più la
sottomatrice superiore sinistra 2 × 2 ortogonale.
Introduciamo queste trasformazioni partendo dalle più specializzate.
2.3.3.1 Isometrie Le isometrie sono trasformazioni del piano R2 che non alterano la distanza
euclidea. Un’ isometria è rappresentata da una matrice del tipo:
 0  


x1
cos(θ) − sin(θ) tx
x1
 x0  =  sin(θ) cos(θ) ty   x2  ,
(68)
2
0
0
1
1
1
con = ±1. Se = 1 l’isometria è detta “orientation-preserving” ed è una trasformazione Euclidea (corrispondente ad una composizione di una traslazione ed una rotazione dell’immagine);
se invece = −1 viene chiamata “orientation-reversing”: un classico esempio di questo tipo di
isometrie è la riflessione con H = diag(−1, 1, 1).
Una trasformazione Euclidea, come ben sappiamo, modella il moto rigido di un oggetto e può
essere scritta in maniera più concisa:
0
R t
x = HE x =
x ,
(69)
0T 1
dove R è una matrice ortogonale 2 × 2 e t rappresenta il vettore traslazione 2 × 1. Una trasformazione Euclidea del piano ha 3 gradi di libertà, uno per la rotazione ed due per la traslazione e
può essere valutata conoscendo la corrispondenza tra due punti prima e dopo la trasformazione.
Tale trasformazione non altera le lunghezze, le aree e gli angoli.
2.3.3.2 Similitudini Una similitudine è un’isometria composta con una scalatura isotropa. Se
l’isometria è una trasformazione Euclidea abbiamo:
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici



0 
x1
s cos(θ) −s sin(θ) tx
x1
0
 x  =  s sin(θ) s cos(θ) ty   x2 
2
0
0
1
1
1
32

,
(70)
oppure, in modo più conciso:
0
x = Hs x =
sR
0T
t
1
x
,
(71)
dove s rappresenta il parametro di scalatura isotropo. Tale trasformazione non altera la forma
dell’immagine (quindi mantiene i rapporti tra le lunghezze e tra le aree) ed ha ovviamente 4 gradi
di libertà.
2.3.3.3 Affinità Una trasformazione affine, o più semplicemente affinità, è una trasformazione
non singolare lineare seguita da una traslazione ed ha tale rappresentazione matriciale:
 0  


x1
a11 a12 tx
x1
 x0  =  a21 a22 ty   x2  ,
(72)
2
0
0
1
1
1
oppure, in modo più conciso:
0
x = Ha x =
A
0T
t
1
x
,
(73)
con A matrice 2 × 2 non singolare. Una affinità planare ha 6 gradi di libertà corrispondenti ai sei
elementi della matrice e può quindi essere ricavata dalle corrispondenze di tre punti.
Per capire gli effetti geometrici che la matrice A provoca sull’immagine, scomponiamola in tal
modo:
A = R(θ)R(−ϕ)DR(ϕ) ,
(74)
con R(θ) e R(ϕ) rotazioni rispettivamente di θ e di ϕ (Figura 15)e D matrice diagonale del tipo:
λ1 0
D=
.
0 λ2
Tale scomposizione deriva direttamente dalla SVD di Ha .
La matrice Ha è quindi una trasformazione composta da:
• una rotazione di ϕ;
• una scalatura di λ1 e λ2 nelle direzioni x e y ruotate;
• una controrotazione di ϕ;
• una rotazione di θ.
La nuova trasformazione rispetto alla similitudine è la scalatura non isotropa su direzioni comunque ortogonali. Tale trasformazione non altera il parallelismo, il rapporto tra segmenti paralleli ed il rapporto tra aree. L’affinità può preservare o no l’orientazione, a seconda del segno del
determinante di A (det(A) ≤ 0 inverte l’orientazione).
2.3.3.4 Proiettività La proiettività è già stata definita nella (66). Si tratta di una trasformazione
lineare non singolare generale di coordinate omogenee. Riscriviamo una proiettività nella forma:
0
A t
x = Hp x =
x ,
(75)
vT v
con il vettore v = (v1 , v2 )T . Tale trasformazione ha 8 gradi di libertà e può essere ricavata dalle
corrispondenze di 4 punti, tre a tre non allineati. Essa lascia invariato il rapporto tra i rapporti di
due linee parallele detto anche “cross-ratio”.
La differenza fondamentale tra una proiettività ed un’ affinità è data dal vettore v. Questo vettore
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
33
rotazione
deformazione
Figura 15: Rotazione di R(θ) e deformazione di D dopo una rotazione di R(ϕ).
è responsabile degli effetti non lineari della proiettività.
Nell’ affinità i punti ideali rimangono all’infinito:

 " #
x
x1
A t  1 
A
x2 =
x2
0T 1
0
0
,
mentre nella proiettività i punti all’infinito sono mappati in punti al finito:

 "
#
x
x1
A t  1 
A
x2 =
x2
.
vT v
0
v1 x 1 + v2 x 2
(76)
(77)
2.3.3.5 Decomposizione di una proiettività Una proiettivà può essere scomposta in una sequenza di trasformazioni più particolari:
sR t
K 0
I 0
A t
H = Hs Ha Hp =
=
,
(78)
vT v
vT v
0T 1
0T 1
con A matrice non singolare data da A = sRK + tv T , K matrice triangolare superiore normalizzata poiché det(K) = 1. Questa scomposizione è valida se v 6= 0, ed è unica se s è scelto
positivo. La matrice Hp è chiamata “elation” e sposta la linea all’infinito.
2.3.3.6 Topologia del piano proiettivo Abbiamo visto che il piano proiettivo P2 può essere
pensato come il set di tutti i vettori omogenei 3 × 1. Un vettore di questo tipo x = (x1 , x2 , x3 )T può
sempre esser normalizzato in modo che x21 + x22 + x23 = 1. Tale punto appartiene ad una sfera di
raggio unitario in R3 . Comunque qualsiasi vettore x ed il suo opposto −x rappresentano lo stesso
punto in P2 , dal momento che differiscono per il fattore di scala −1. Esiste quindi una corrispondenza due-a-uno tra la sfera unitaria S 2 in R3 ed il piano proiettivo P2 . Il piano proiettivo potrebbe
perciò esser immaginato come una sfera unitaria con i punti diametralmente opposti sulla sfera
coincidenti. Topologicamente, per rappresentare i punti di P2 , è sufficiente una sola semisfera
(poiché l’altra parte di sfera originerebbe gli stessi punti) con la circonferenza diametrale che ha
i punti opposti coincidenti. Ma poiché una semisfera è topologicamente equivalente ad un disco
piano, possiamo immaginare P2 semplicemente come un disco avente i punti diametralmente opposti sul bordo coincidenti. Questa superficie non è fisicamente realizzabile. Si dimostra che P2
non è semplicemente connesso.
Il piano proiettivo inoltre non è orientabile:
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
34
Figura 16: Rappresentazione del piano proiettivo; l’orientazione locale (rappresentata dalla L nei vari riquadri), si inverte quando attraversa i confini del piano proiettivo .
questo vuol dire che non è possibile definire un’orientazione locale in P2 che sia consistente
sull’intera superficie (vedi Figura 16). Questo concetto ci sarà di aiuto in seguito per meglio comprendere il comportamento apparentemente “anomalo” del controllo in alcune situazioni.
2.3.3.7 Punti fissi e linee fisse Ci sono certi punti del piano proiettivo P2 che rimangono fissi
dopo la trasformazione proiettiva. Questi punti sono chiamati punti fissi della trasformazione.
Analizziamo adesso questo argomento in maniera più approfondita.
0
Affinché il punto di partenza x ed il punto di arrivo x siano fissi, devono coincidere, cioè deve
valere:
He = λe
(79)
;
dunque il punto fisso deve coincidere con l’autovettore della trasformazione relativo all’autovalore
λ (infatti e e λe rappresentano nel piano proiettivo lo stesso punto). Spesso gli autovettori o gli
autovalori di una trasformazione hanno un significato fisico o geometrico particolare. Una matrice
immagine 1
immagine 2
H
e
e
e
e
e
e
Figura 17: Un’omografia planare ha in generale tre punti fissi distinti che individuano tre rette fisse come
set ma non puntualmente.
3 × 3 ha tre autovalori e conseguentemente una trasformazione proiettiva ha fino a tre punti fissi,
come in Figura 17. Poiché l’equazione caratteristica della trasformazione è in questo caso una
cubica, uno o tre autovalori (ed i corrispondenti autovettori) sono reali. Dato che le linee si trasformano secondo la (67), una simile trattazione può essere fatta per le rette fisse che corrispondono
agli autovettori di (H −1 )T .
Si può notare che tali rette, individuate dai punti fissi, sono in generale fisse come set ma non
puntualmente: ed esempio un punto generico su una linea fissa viene rimappato sulla stessa
linea ma non sullo stesso punto. Si possono poi verificare dei casi in cui si hanno due autovalori
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
35
uguali. Supponiamo perciò che i primi due autovalori λ1 e λ2 della trasformazione H siano uguali
e che ci siano due autovettori distinti (e1 , e2 ) corrispondenti; allora la retta passante per e1 ed e2
sarà fissa puntualmente.
Se infatti supponiamo x = αe1 + βe2 , risulta:
Hx = λ1 αe1 + λ1 βe2 = λ1 x
(80)
.
Un’altra possibilità è che valga λ1 = λ2 ma ci sia un solo autovettore corrispondente. In questa
situazione, l’autovettore ha molteplicità geometrica inferiore alla molteplicità algebrica del corrispettivo autovalore nell’equazione caratteristica: c’è un punto fisso in meno (due invece di tre). Si
fa notare che il gruppo affine ed i gruppi più specializzati hanno due punti fissi all’infinito (dunque
la retta all’infinito si mantiene in tal caso fissa ) che corrispondono agli autovettori della sottomatrice superiore sinistra 2 × 2 estratta da H; l’altro punto fisso è in genere finito.
Queste considerazioni ci saranno utili nella progettazione della pianificazione del controllo basato
su visione.
2.3.3.8 Calcolo Omografia tra 2 immagini In questo Paragrafo si affronta il problema del calcolo dell’omografia esistente nel piano immagine tra due viste differenti di uno stesso oggetto
piano.
Il metodo esposto è lo stesso che utilizziamo nell’algoritmo di controllo. Una omografia, come
già specificato, non è altro che una mappatura invertibile fra punti del piano proiettivo P2 , in punti
dello stesso piano, fatta in modo da trasformare linee in linee . I punti del piano proiettivo sono
rappresentati da vettori omogenei ∈ R3×1 .
Ad ogni omografia, o proiettività, corrisponde una matrice omogenea non singolare ∈ R3×3 . Tale
matrice si dice omogenea poiché non conta il valore che assumono i suoi elementi, bensı̀ il loro
rapporto, quindi ha in totale 8 gradi di libertà.
Se abbiamo a disposizione le coordinate di 4 punti (3 a 3 non allineati) in due immagini differenti
e conosciamo le corrispondenze tra i punti stessi, allora possiamo procedere cosı̀:

h11
sia
Hd =  h21
h31


p̃i1
siano p̃i =  p̃i2 
p̃i3
h12
h22
h32
e

h13
h23 
h33

l’omografia cercata,

0
p̃i1
0
p̃i =  p̃i2 
0
p̃i3
0
per i = 1, ..., 4
le coordinate omogenee dei quattro punti noti nelle due immagini. Per ogni coppia di punti corrispondenti vale la relazione:
p̃0i = Hd p̃i

 
p˜0i1
h11
 ˜0  
h21
=
 pi2 
0
˜
h31
pi3
h12
h22
h32
ovvero


h13
p˜i1
h23   p˜i2 
h33
p˜i3
i = 1, ..., 4.
Sviluppiamo la relazione appena descritta per un generico punto i :
p˜0i1 = h11 p˜i1 + h12 p˜i2 + h13 p˜i3
p˜0i2 = h21 p˜i1 + h22 p˜i2 + h23 p˜i3
p˜0i3 = h31 p˜i1 + h32 p˜i2 + h33 p˜i3
(81)
(82)
.
(83)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
36
Modellistica e Controllo di Sistemi Meccanici
Ricaviamo quindi 3 nuove equazioni ottenute moltiplicando in maniera incrociata i membri delle
precedenti. Otteniamo rispettivamente da (81) e (82), (81) e (83) e da (82) e (83) :
h21 p˜i1 p˜0i1 + h22 p˜i2 p˜0i1 + h23 p˜i3 p˜0i1 = h11 p˜i1 p˜0i2 + h12 p˜i2 p˜0i2 + h13 p˜i3 p˜0i2
h31 p˜i1 p˜0i1 + h32 p˜i2 p˜0i1 + h33 p˜i3 p˜0i1 = h11 p˜i1 p˜0i3 + h12 p˜i2 p˜0i3 + h13 p˜i3 p˜0i3
h31 p˜i p˜0 + h32 p˜i p˜0 + h33 p˜i p˜0 = h21 p˜i p˜0 + h22 p˜i p˜0 + h23 p˜i p˜0
1
i2
2
i2
3
i2
1
i3
2
i3
3
i3
.
Portando tutto al primo membro otteniamo 3 equazioni omogenee nelle 9 incognite hii .
Ripetendo questo procedimento per tutte le 4 corrispondenze note, alla fine si giunge ad un sistema omogeneo lineare di 12 equazioni nei 9 elementi di H che, espresso in forma matriciale,
risulta:
P h = 0 P ∈ R12×9
h = [h11
h12
h13
h21
h22
h23
h31
(84)
h32
h33 ]T .
(85)
Se tutte le equazioni del sistema fossero linearmente indipendenti (cioè se il rango di P fosse 9),
il sistema ammetterebbe per h la sola soluzione nulla; tuttavia il rango di P è 8 ed essa possiede
un kernel di dimensione 1. Le reali incognite del problema sono gli 8 rapporti h11 : h12 : h13 :
h21 : h22 : h23 : h31 : h32 : h33 (Hd è definita a meno di una costante moltiplicativa). Esistono
perciò solo 8 equazioni linearmente indipendenti. La soluzione per h è una qualsiasi base del
N (P ) ∈ R9×1 .
2.3.3.9 Structure from Motion Nel precedente paragrafo abbiamo affrontato il problema del
calcolo dell’omografia fra due immagini dalle corrispondenze di 4 punti; vediamo adesso come
sono in relazione fra loro due viste di un oggetto piano acquisite da una stessa telecamera in
posizioni differenti (nell’asservimento visivo infatti risulta interessante valutare le relazioni che
legano la geometria delle immagini prodotte agli spostamenti relativi telecamera-scena).
Il problema della determinazione dei parametri cinematici di spostamento della telecamera corrispondente a due differenti immagini è chiamato Structure From Motion problem (SFM).
n
di
Pi
Pf
t
C
f
i
i
C
R
f
i
f
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
37
Si considerino due viste Initiale Ii and Finale If di un bersaglio planare definito da punti; abbiamo
P f = Rf i P i + t f i
⇒
Zf
1
mf = Rf i mi + tf i
Zi
Zi
,
(86)
dove abbiamo usato la (46). Poiché il punto ∈ al piano π, nTi P i = di ⇒ 1/Zi = (nTi mi )/di da cui:
Zf
nTi mi
Zf
nTi
mf = Rf i mi + tf i
,⇒
mf = Rf i + tf i
mi = Hf i mi ,
(87)
Zi
di
Zi
di
Hf i è chiamata Omografia Euclidea; questa definisce la mappature tra due immagini di un bersaglio
piano con i punti espressi in coordinate normalizzate. Risulta dunque:
nT
=
Rf i + tsf i nTi
,
(88)
Hf i = Rf i + tf i i
di
T
n
dove tsf i = dii è la traslazione scalata tra la terna telecamera nella posa finale Cf e la terna
telecamera nella posa iniziale Ci . Hf i può anche essere espressa da:
Hf i = Rif T (I − tsif ni T ) ,
(89)
dove tutti i termini sono espressi in terna Ci .
L’eq.(87) utilizzando le coordinate in pixel diventa:
pf ∝ Gf i pi
(90)
,
dove (assumendo la stessa matrice di calibrazione K nelle 2 viste) :
Gf i = KHf i K −1
.
(91)
G è l’omografia tra due viste in pixel di un oggetto piano.
Da due immagini calibrate (cioè nota la matrice K) di un oggetto piano è quindi possibile ricavare:
• L’omografia Gf i fra le due viste utilzzando il metodo precedentemente illstrato, se si conoscono
almeno 4 corrispondenze di punti;
• Dalla Gf i nota K si può facilmente calcolare l’omografia euclidea Hf i
• Esistono delle procedure geometriche per decomporre Hf i e stimare i principali parametri
cinematici di spostamento della telecamera fra le due viste: in particolare è possibile ricavare Rif , tsif ed ni (si noti che conoscendo tsif il modulo della traslazione reale non è
noto ma solo la sua direzione). I parametri cinematici possono quindi essere stimati dalle
sole informazioni visive a meno di una fattore di scala.
NOTE:
• se Ci e π sono assegnate H possiede solo 6 DOF, i gradi di libertà della telecamera (mentre
come abbiamo visto un’omografia generica possiede 8 DOF);
• anche se il bersaglio non è planare esistono tecniche che calcolano la SFM riferendosi ad
un piano virtuale che passa per 3 punti non allineati caratteristici del bersaglio.
2.4
Asservimento visivo di robot
Dopo aver chiarito alcuni i concetti fondamentali che sono propedeutici alla comprensione del
Visual Servoing, andiamo a presentare più in dettaglio le tecniche di controllo ad asservimento
visivo, considerandone pregi e difetti.
Il concetto di asservimento visivo è quello di guidare il sistema verso la configurazione desiderata
attraverso un loop in retrazione che sfrutti il lo stream video fornito da una o più telecamere.
Tutti i sistemi ad asservimento visivo, durante il controllo del sistema, sfruttano alcune primitive
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
38
estratte dall’immagine in real-time : le Feature Immagine. Le Feature immagine possono essere
costituite da punti, rette, coniche o altre primitive geometriche rilevabili sul bersaglio inquadrato
dalla telecamera. Uno dei compiti richiesti più comunemente ad un sistema asservito con telecamere è il positioning task che prevede il controllo del posizionamento del sistema rispetto ad
un bersaglio fisso.
Per l’esecuzione di
• Registrare un’immagine desiderata del bersaglio (immagine di riferimento) corrispondente
alla posizione desiderata del sistema;
• Il sistema in generale si troverà in una configurazione iniziale diversa da quella desiderata
a cui corrisponderà un’immagine attuale del bersaglio;
• Dalle immagini iniziale e desiderata occorre a questo punto estrarre delle feature immagine
(ad esempio dei punti ben rilevabili sul bersaglio): per fare ciò si utilizzano tecniche di elaborazione digitale delle immagini (Image Processing) che portano generalmente all’identificazioni
delle coordinate di corner o di egde sull’immagine (fase di Feature Extraction);
• Una volta note le coordinate delle feature nelle due immagini si passa alla fase del Feature
Matching: si cerca di risolvere il problema delle corrispondenze fra le feature delle due immagini (ad esempio si cerca di capire qual’è il punto corrispondenti dell’immagine desiderata
ad un determinato punto nell’immagine iniziale); questo problema risulta molto complesso
da automatizzare e solitamente è l’utente che manualmente individua le corrispondenze tra
le feature;
• A questo punto si definisce un errore e di sistema funzione delle feature immagine attuali
e di riferimento che andrà minimizzato attraverso un’opportuna legge di controllo; l’input di
controllo dovrà eccitare gli attuatori del sistema dinamico per guidarlo nella configurazione
corrispondente all’immagine desiderata in cui e = 0;
• Nella fase di controllo ad ogni passo di integrazione intervallo di campionamento del visual
controller andrà aggiornato il valore corrente dell’errore, andranno cioè stimate in real-time
le successive coordinate delle feature immagine correnti: tale processo di stima prende il
nome di Tracking delle feature. Il Tracking delle feature è una procedura automatica (la più
dispendiosa di tutte le fasi dal punto di vista computazionale) che può essere paragonata
ad un Matching locale delle feature: nota la posizione delle feature nell’immagine all’istante
precedente si cercano in un intorno abbastanza piccolo (solitamente una finestra quadrata
di 10 pixel se si tratta di punti) le feature nell’immagine corrente.
2.4.1
Classificazione di sistemi con asservimento visivo
A seconda di come è definito l’errore di sistema e i sistemi ad asservimento diviso possono essere
suddivisi in 3 tipologie:
• Se e è definito da primitive cartesiane 3D (come la rotazione e la traslazione tra due terne)
si parla diPosition Based Visual Servoing (PBVS) o di 3D Visual Servoing;
• Se l’errore e è definito da primitive immagine 2D (attraverso esempio punti e rette nell’immagine)
si esempio punti e rette nell’immagine) si parla di Image Based Visual Servoing (IBVS) o di
2D Visual Servoing;
• Se e è definito sia da primitive 3D che da primitive immagine 2D si parla di Hybrid Visual
Servoing o di 2 − 1/2D Visual Servoing.
Analizziamo adesso i principali vantaggi e svantaggi di questi tre tipi di asservimento visivo.
2.4.1.1 Position Based Visual Servoing (PBVS) Nel PBVS la legge di controllo utilizza un
errore direttamente legato alla posa 3D ella telecamera.
VANTAGGI:
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
39
• Il controllo assegna direttamente la traiettoria 3D del sistema: si ha quindi un pieno controllo
della gestione della traiettoria cartesiana e la legge di controllo risulta di più comune utilizzo
(ad esempio il controllo nello spazio operativo del robot).
SVANTAGGI:
• Il bersaglio potrebbe uscire dal campo visivo della telecamera durante il controllo (failure del
sistema);
• Se la telecamera non è perfettamente calibrata si ha un grasso errore nella stima dei
parametri Cartesiani della scena con un rapido degrado delle prestazioni del sistema (il
sistema raggiunge una configurazione diversa da quella desiderata);
• è necessario il modello 3D del bersaglio.
2.4.1.2 Image Based Visual Servoing (IBVS) Nell’IBVS la funzione errore è espressa tutta
nello spazio immagine e quindi il loop di retroazione si chiude tutto nell’immagine.
VANTAGGI:
• Il sistema presenta un’elevata robustezza rispetto ad errori di calibrazione della telecamera
(funziona anche utilizzando una stima approssimativa della matrice K);
• Non è necessario il modello 3D del bersaglio (approccio model-free);
SVANTAGGI:
• la stabilità del sistema è assicurata analiticamente dalle attuali leggi di controllo soltanto in
un intorno della posizione di equilibrio (stabilità locale): il sistema risulta infatti accoppiato e
non lineare (nascono inoltre problemi di singolarità e di minimi locali della legge di controllo);
• la traiettoria 3D del sistema non è facilmente individuabile: potrebbe anche risultare irrealizzabile.
• Il bersaglio potrebbe uscire dal campo visivo della telecamera durante il controllo (failure del
sistema);
2.4.1.3 Hybrid Visual Servoing Nell’ Hybrid Visual Servoing l’errore è calcolata in parte nell’immagine
ed in parte nello spazio 3D.
VANTAGGI:
• Controlla in maniera disaccoppiata rotazione e traslazione della telecamera;
• Sono state analiticamente provate condizioni di stabilità asintotica per la stabilità locale ed
è stato individuato il dominio di convergenza del sistema rispetto ad errori di calibrazione.
• è un approccio model-free
SVANTAGGI:
• presenta alcuni problemi di runore quando il sistema si porta vicino alla convergenza;
• Non controlla il moto 3D del sistema.
• Il bersaglio potrebbe uscire dal campo visivo della telecamea durante il controllo (failure del
sistema);
Nel seguito analizzeremo a fondo i sistemi IBVS che sono i più robusti rispetto ad errori di
calibrazione della telecamera e del sistema meccatronico: l’input del visual controller infatti sarà
diverso da finché rimane presente un errore nell’immagine.
Allotta, Fioravanti, Malvezzi
2.4.2
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
40
Controllo dei sistemi ad Asservimento Visivo
In questo paragrafo faremo riferimento all’ architettura più comune di sistema asservito: si utilizza
una singola telecamera con configurazione eye-in-hand e bersaglio fisso. Determiniamo nel sistema k primitive indipendenti xi con i = 1...k, 2D o 3D, funzione delle feature immagine con k ≥
dei DOF del sistema da controllare (per il controllo di un manipolatore nello spazio libero avremo
k ≥ 6).
Se le feature immagine coincidessero con tutte le k primitive il sistema asservito sarebbe di tipo
Image Based: l’esempio più semplice è assumere come primitive indipendenti le coordinate mx
ed my di una serie di punti normalizzati sull’immagine (vedremo che comunque saranno necessari almeno 4 punti immagine → k = 8).
Il vettore formato dalle k primitive identifica lo stato del sistema da controllare definito dal vettore
x ∈ Rk×1 . Sia e l’errore del sistema da minimizzare espresso da:
e = (xd − x) ∈ Rk×1
,
(92)
dove xd rappresenta lo stato di riferimento del sistema ed x lo stato corrente.
Sia r(t) il vettore che definisce posizione ed attitudine della telecamera ∈ R6×1 (ad esempio
l’origine della terna telecamera ed il suo vettore di rotazione rispetto ad al sistema fisso), risulta
che:
x = x(r(t)) ,
(93)
cioè, come è naturale intuire, lo stato del sistemaè funzione della posa della telecamera.
Sia w il Twist-Screw della telecamera composto dalla velocità di traslazione e dalla velocità angolare della terna telecamera c (rispetto ad una terna fissa) espresse in terna c :
v
w=
= [vx vy vz ωx ωy ωz ]T .
(94)
ω
La (93) tradotta in termini differenziali diventa:
ẋ = J(x, r(t))w
.
(95)
La matrice J(x, r(t)) ∈ Rk×6 è chiamata Matrice di Interazione o Jacobiano Immagine. J rappresenta una mappatura lineare locale fra la variazione delle primitive che identificano lo stato del
sistema ed il twist-screw della telecamera. Normalmente J dipende da:
1. lo stato attuale del sistema (è funzione della configurazione);
2. la posa della telecamera.
Per il controllo del sistema è naturale utilizzare una legge di controllo del tipo:
w = Jˆ† (ke + ẋd ) ,
(96)
ˆ −1 JˆT è la pseudoinversa sinistra dello Jacobiano immagine (stimato) e k è un
dove Jˆ† = (JˆT J)
guadagno scalare positivo.
La (96) presenta un termine proporzionale all’errore di sistema ed un eventuale termine di feedˆ = 6 ovvero Jˆ ha rango
forward (se xd 6= cost ): questa legge può essere utilizzata solo se rank(J)
massimo).
xd rappresenta la traiettoria desiderata per le feature immagine ed è il riferimento per il sistema
di controllo:
1. senza nessuna pianificazione ⇒ xd = xf = cost (il riferimento coincide con lo stato finale
desiderato);
2. con pianificazione del riferimento xd 6= cost varia dallo stato iniziale xi a quello finale
desiderato xf .
Allotta, Fioravanti, Malvezzi
2.4.3
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
41
Stabilità del Sistema
Per indagare sulla stabilità del sistema si usa il metodo di Lyapunov; si definisce quindi una
funzione dell’errore di sistema V = V (e) definita positiva detta Candidata di Lyapunov:
V =
1 T
e e
2
(97)
.
• Se la derivata della V rispetto al tempo risulta una funzione definita negativa ovvero V̇ < 0
∀ configurazione del sistema con e 6= 0, il sistema controllato risulta globalmente stabile;
• Se V̇ risulta semidefinita negativa ovvero V̇ ≤ 0 ∀ configurazione del sistema con e 6= 0 la
convergenza è assicurata solo in un intorno della posizione di equilibrio (sistema localmente
stabile).
Sfruttando (95), (97) e la (96) nella espressione di V̇ otteniamo:
V̇
= eT ė
V̇
= eT (ẋd − ẋ)
V̇
V̇
= eT (ẋd − Jw)
= eT (ẋd − J Jˆ† (ke + ẋd ))...
V̇
= −keT J Jˆ† e + eT (I − J Jˆ† )ẋd
(98)
.
• il termine −keT J Jˆ† e della (99) risulta semidefinito negativo, infatti risulta 0 sia per e = 0
che per e ∈ N (J T ).
Se J è quadrato tale termine diventa invece definito negativo;
• sul segno di eT (I − J Jˆ† )ẋd niente si può dire; questo termine comunque si annulla in tre
casi:
– se ẋd = 0 cioè a riferimento costante (se non si esegue alcuna pianificazione);
– se ẋd oppure e appartengono allo spazio immagine di J, poiché il termine (I − J Jˆ† )
rappresenta un proiettore nello spazio N (J T );
– se J è quadrato (in tal caso avremmo N (J) = 0.
Per quanto detto finora è possibile affermare che:
1. se non si esegue la pianificazione del riferimento (ẋd = 0 ⇒ xd = xf ) il sistema è localmente
stabile: la convergenza è assicurata analiticamente solo in un intorno della posizione di
equilibrio); si avrebbe stabilità asintotica in caso di J quadrato (infatti in tal caso avremmo
N (J) = 0);
2. se si esegue una pianificazione del riferimento per il controllo (ẋd 6= 0) e ẋd 6= N (J)T il
sistema risulta localmente stabilà Tale condizione in realtà risulta difficile da realizzare in
quanto dipende dal tipo di pianificazione effettuata: vedremo in seguito alcuni esempi. Se
J fosse quadrato avremmo anche in questo caso stabilità asintotica.
NOTE:
• l’analisi di stabilità svolta precedentemente è valida solo se J ha rango massimo (nel caso
di matrice quadrata, quando lo Jacobiano non è singolare): le configurazini in cui J perde
rango sono dette singolarità del sistema asservito. Se il sistema passa vicino ad una configurazione singolare la (96) perde significato, il sistema diventa instabile e si ha uno stallo
del controllo: tali configurazioni sono quindi da evitare.
Le caratteristiche di singolarità dello Jacobiano dipendono direttamente dalle primitive scelte
per conporre il vettore x: bisogna quindi porre particolare attenzione nella scelta del numero
e del tipo di primitive che compongono lo stato x.
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
42
Modellistica e Controllo di Sistemi Meccanici
• se non si esegue la pianificazione del riferimento (ẋd = 0) e lo Jacobiano non è quadrato
ci sono configurazioni in cui il sistema può arrestarsi in una posizione differente da quella
desiderata (ovvero si ha che w = 0 con e 6= 0) tali situazioni prendono sono chiamati
minimi locali del sistema e si verificano quando e ∈ N (J)T . Tali situazioni si verificano
quando per arrivare nella configurazione desiderata la legge di controllo spinge il sistema
in configurazioni intermedie irrealizzabili. Tali situazione dovrebbero essere evitate se si
pianifica il riferimento opportunamente.
• Nei ragionamenti fatti in precedenza si ipotizza che Jˆ = J ovvero che lo Jacobiano stimato
per implementare la legge di controllo (96) corrisponda allo Jacobiano reale: questa condizione risulata molto difficile da realizzare in quanto Jˆ dipende, oltre che dal vettore di stato,
anche dalla posa della telecamera r(t) che di solito non è totalmente nota; d’altra parte una
calibrazione troppo imprecisa può portare ad una stima di x e quindi di Jˆ grossolana.
Per tenere conto dell’incertezza sulla stima dei parametri occorrerebbe eseguire un’analisi
di robustezza ed individuare il dominio di convergenza del sistema in funzione dell’errore di
stima: questa analisi talvolta risulta molto complessa.
Nella pratica comunque i sistemi ad asservimento visivo di tipo IBVS risultano essere molto robusti rispetto ad errori di calibrazione sia della telecamera che del sistema meccatronico controllato.
2.4.4
I Punti Immagine per l’identificazione dello Stato
Il caso più semplice suggerisce di scegliere per il nostro sista asservito dei semplici punti immagine ben rilevabili sul bersaglio come image feature e di utilizzarli anche come primitive per
costruire il vettore di stato x.
Supponiamo allora di scegliere 4 punti complanari sul nostro bersaglio rispetto al quale si desidera
eseguire dei positioning task (la scelta di 4 punti sarà giustificata in seguito). Siano
mxi
Xci /Zci
mi =
=
,
(99)
myi
Yci /Zci
con i = 1...4 le coordinate normalizzate dei 4 punti in analisi definiti come nella (46); questi
possono eser ricavati dai rispettivi punti in pixel invertendo la (49), supponendo di assumere un
modello di telecamera CCD full projective con matrice K nota e costante durante il compito.
Costruiamo allora il vettore di stato x utilizzando come primitive direttamente le coordinate dei 4
punti normalizzati:
x(t) =
mx1
my1
mx2
my2
mx3
my3
mx4
my4
T
.
(100)
L’errore di sistema (92) risulta:


md1 − m1
 md2 − m2 

e(t) = 
 md3 − m3 
md4 − m4
,
(101)
cioè la differenza tra le coordinate di riferimento e le coordinate attuali dei punti normalizzati. Con
tale scelta il vettore di stato e, quindi, l’errore risultano tutti definiti nell’immagine: si parla quindi
di Image based Visual Servoing (IBVS).
2.4.4.1 Matrice di Interazione del punto immagine Per l’implementazione della legge di controllo (96) è necessario ricavare la matrice J: dobbiamo cioè trovare qual’è la relazione che lega
il twist-screw di telecamera w alla velocità dei punti nell’immagine normalizzata. Descriviamo
quindi il moto relativo della telecamera rispetto all’bersaglio supposto in quiete, ricordando che il
camera-twist-screw
v
w=
,
ω
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
43
è composto dalla velocità di traslazione v dell’origine della terna telecamera < C > e dalla velocità angolare di < C > ω , entrambi espressi nel sistema di riferimento < C >.
Innanzitutto, vediamo come esprimere la velocità di un punto del bersaglio in funzione dello screw
di velocità w definito sopra.
T
Tale punto ha al solito coordinate in terna < C > date da: P c = Xci Yci Zci
.
Dalla cinematica relativa ricaviamo:
Ṗ o(c) = Ċ o(c) + ω co(c) ∧ (P − C)(c) + Ṗ c(c)
(102)
,
avendo indicato con C l’origine della terna telecamera ed inoltre con:
• Ẋ i(j) la velocità del punto X dello spazio 3D rispetto alla terna I, espresso nelle coordinate
della terna J;
• ω ij(k) la velocità angolare della terna I rispetto alla terna J, espressa nelle coordinate della
terna K;
• (X − Y )i il vettore orientato Y al punto X espresso in terna I.
poiché l’oggetto è solidale appunto ad < O >, avremo:
Ṗ o(c) = Ṗ o(o) = 0
(103)
;
sostituendo la (103) nella (102) si ricava:
Ṗ c(c) = −Ċ o(c) − ω co(c) ∧ (P − C)(c)
(104)
.
Possiamo osservare che, per come abbiamo definito w, risulta:
e
Ċ o(c) = v
ω co(c) = ω
quindi la (104) diventa:
Ṗ c = −v − ω ∧ (P − C)
Sviluppando l’equazione vettoriale (105), risulta:


 

Ẋci
vx
0
 Ẏci  = −  vy  −  ωz
−ωy
vz
Żci
−ωz
0
ωx
(105)
.


Xci
ωy
−ωx   Yci 
Zci
0
;
(106)
da cui si ricava:
Ẋci
= −vx + ωz Yci − ωy Zci
(107)
Ẏci
= −vy − ωz Xci + ωx Zci
(108)
Żci
= −vz + ωy Xci − ωx Yci
.
(109)
Ricordando la (??) :
mi =
mxi
myi
=
Xc /Zc
Yc /Zc
e derivando tale relazione rispetto al tempo risulta:
#
"
1
[(Ẋc Zc ) − Żc Xc ]
Zc 2
ṁi =
1
[(Ẏc Zc ) − Żc Yc ]
Zc 2
(110)
.
.
(111)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
44
Modellistica e Controllo di Sistemi Meccanici
A questo punto, sostituendo le equazioni (107), (108) e (109) nella (111), ricaviamo:

ṁi


= 

−
1
Zi
0
0
1
−
Zi
xi
Zi
yi
Zi
− 1 + x2i
xi yi
1+
yi2
−xi yi
yi
−xi



 w = Jp w

,
(112)
La matrice ottenuta nell’equazione precedente è lo jacobiano immagine al punto immagine normalizzato mi .
Allotta, Fioravanti, Malvezzi
3
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
45
Controllo sliding mode di sistemi in forma canonica di controllo
3.1
Controllo sliding mode di sistemi ad un solo ingresso
Il controllo sliding mode (o a struttura variabile) è una tecnica di controllo nonlineare con notevoli
proprietà di robustezza. Lo esamineremo inizialmente con particolare riferimento al controllo di
un sistema nonlineare del II ordine ad un ingresso con incertezza di modellazione.
3.1.1
Definizione della superficie di sliding
Sia dato un sistema nonlineare (lineare rispetto all’unico ingresso u ovvero input-affine) del tipo:
x(n) = f (x) + b(x)u
,
(113)
dove x = [x, ẋ, . . . , x(n−1) ]T è il vettore di stato, f (x) sia approssimata con incertezza limitata
da una funzione fˆ(x), ovvero |f (x) − fˆ(x)| ≤ F (x), con F (x) nota, b(x) sia una funzione di
segno noto, e sia assegnata una traiettoria desiderata del sistema nello spazio di stato. Quindi
l’obbiettivo del sistema di controllo (problema di asservimento) può essere definito come segue:
x
= xd (t)
ẋ
= ẋd (t)
ẍ
= ẍd (t)
...
= ...
x(n−1)
= xd
x(n)
(n−1)
(t)
(n)
= xd (t)
Sia x̃ = x − xd (pari all’errore e = xd − x cambiato di segno); la sliding surface è una superficie
tempovariante S(t) ∈ Rn definita dall’equazione:
(114)
s=0 ,
dove s è uno scalare, detto sliding variable, definito come segue:
s = (λ +
d n−1
)
x̃
dt
,
(115)
Dove λ > 0 è una costante fissata dal progettista. Per un sistema del II ordine n = 2, quindi, la
varibile di sliding vale:
s = λx̃ + x̃˙
(116)
.
Se lo stato del sistema è tale da rimanere sulla superficie di sliding), la dinamica d’errore è
governata dall’equazione s = 0 ossia, nel caso del sistema del II ordine, dalla seguente equazione
differenziale:
λx̃ + x̃˙ = 0 .
(117)
La (117) è lineare tempo invariante e quindi x̃, da una condizione iniziale x̃(t0 ) = x̃0 , convergerà
esponenzialmente a zero con legge temporale:
x̃ = x̃0 e−λ(t−t0 )
.
(118)
L’appartenenza dello stato alla superficie di sliding è quindi un punto cruciale del controllo sliding
mode. Possiamo riassumere gli obbiettivi del controllo come segue:
• portare lo stato del sistema il più rapidamente possibile sulla superficie di sliding;
• mantenere lo stato del sistema sulla superficie.
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
46
La tecnica sliding mode consiste nell’applicare un controllo in retroazione composto da due parti:
• una prima parte del segnale di controllo è basata sul modello disponibile del sistema e è
tale da permettere la permanenza dello stato sulla sliding surface in caso di conoscenza
esatta del modello;
• La seconda parte del segnale di controllo serve a fornire robustezza rispetto alle incertezze
del modello ed al rumore di misura e, in pratica, a riportare velocemente il sistema sulla
sliding surface, qualora lo stato tendesse a lasciarla.
Il sistema di controllo deve assicurare che la superficie di sliding venga raggiunta in un tempo
finito e non venga abbandonata. Queste condizioni sono entrambe soddisfatte se è verificata la
relazione:
1 d 2
s ≤ −η|s| ,
2 dt
(119)
dove η è una costante strettamente positiva. L’equazione precedente può essere manipolata
ottenendo:
ṡs ≤ −η|s| ,
(120)
o anche:
ṡ sign(s) ≤ −η
(121)
.
Il significato della (121) è che la derivata temporale di s ha segno opposto a s ed è in modulo
maggiore di un numero η. Quindi ṡ è tale da attrarre lo stato verso S(t) su entrambe le due
porzioni dello spazio di stato caratterizzate rispettivamente da s > 0 ed s < 0.
Riguardo al tempo di raggiungimento di S(t) supponiamo, per fissare le idee, che s(t0 ) > 0
e che sia tf l’istante in cui lo stato raggiunge la superficie di sliding. Quindi, integrando la (121),
otteniamo:
Z tf
Z tf
ṡdt ≤ −η
dt
t0
t0
t
s|tf0
≤
−η(tf − t0 )
s(tf ) − s(t0 ) ≤
−η(tf − t0 )
s(t0 )
.
η
(tf − t0 ) ≤
(122)
Ripetendo il ragionamento a partire da una condizione iniziale s(t0 ) < 0, si ottiene la condizione
(tf − t0 ) ≤ − s(tη0 ) , quindi in generale:
(tf − t0 ) ≤
3.1.2
|s(t0 )|
η
.
(123)
Sintesi della legge di controllo
Al fine di assicurare il rispetto della condizione di sliding (120), il controllo viene costruito come
segue.
3.1.2.1 Parte model-based La permanenza dello stato sulla superficie di sliding può essere
espressa dall’equazione ṡ = 0:
¨ + λx̃˙ = 0 .
ṡ = x̃
(124)
Sostituendo nella (124) l’espressione della dinamica del sistema (113 che, per il sistema del II
ordine, risulta ẍ = f (x) + b(x)u), si ottiene:
¨ + λx̃˙ = ẍ − ẍd + λx̃˙ = f (x) + b(x)u − ẍd + λx̃˙ = 0 .
x̃
(125)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
47
La parte del segnale di controllo model-based û si costruisce in modo che l’equazione precedente
sia verificata in caso di conoscenza esatta del modello del sistema, ovvero:
1 ˆ
û =
−f (x) + ẍd − λx̃˙
.
(126)
b̂(x)
3.1.2.2 Parte discontinua su S(t) La parte discontinua su S(t) serve a riportare il sistema
sulla superficie di sliding nel caso che lo stato se ne allontanasse per qualche ragione. Essa ha
la forma:
−k
1
b̂(x)
sign(s) ,
(127)
dove k è una costante (k > 0 se b(x) > 0 e k < 0 se b(x) < 0).
3.1.2.3 Verifica della condizione di sliding Il segnale di controllo è quindi la somma di û
della (126) e del termine discontinuo della (127):
1 ˆ
−f (x) + ẍd − λx̃˙ − k sign(s)
.
(128)
u=
b̂(x)
Affinché la condizione di sliding sia verificata è necessario dimensionare in modo opportuno il
guadagno k sulla base del modello d’incertezza. Esaminiamo i vari casi.
Modello noto con esattezza In tal caso sarà fˆ ≡ f e b̂ ≡ b. Sostituendo la (128) nella (120)
si ottiene:
ṡs =
=
¨ + λx̃)s
˙
(x̃
˙
(ẍ − ẍd + λx̃)s
˙
(f (x) + b(x)u − ẍd + λx̃)s
"
#
b(x) ˆ
˙
˙
=
f (x) +
−f (x) + ẍd − λx̃ − k sign(s) − ẍd + λx̃ s
b̂(x)
h
i
= f (x) − fˆ(x) + ẍd − λx̃˙ − k sign(s) − ẍd + λx̃˙ s
=
= −k sign(s)s = −k|s| .
(129)
La condizione di sliding è quindi rispettata se −k|s| ≤ −η|s|, ovvero k ≥ η.
Incertezza su f
sua stima:
Supponiamo adesso che b sia ancora nota con esattezza e valga per f e la
|f (x) − fˆ(x)| < F (x) ,
(130)
con F (x) nota. Ripetendo i conti della (129) nel caso di stima incerta di f (x) si ottiene:
"
#
b(x) ˆ
−f (x) + ẍd − λx̃˙ − k sign(s) − ẍd + λx̃˙ s
ṡs =
f (x) +
b̂(x)
h
i
= f (x) − fˆ(x) + ẍd − λx̃˙ − k sign(s) − ẍd + λx̃˙ s
=
(f (x) − fˆ(x))s − k sign(s)s
.
(131)
Vogliamo che sia ṡs ≤ η|s| e possiamo soddisfare tale condizione imponendo che:
F (x)|s| − k|s| =
(F (x) − k)|s| ≤ −η|s| .
La condizione di sliding è allora verificata se (F (x) − k) ≤ −η, ovvero k ≥ F (x) + η
(132)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
48
Incertezza su b Supponiamo adesso che f sia nota con esattezza e il guadagno b sia invece
incerto. Assumiamo per b un’incertezza di tipo moltiplicativo:
ξ1 ≤ 1 −
In tal modo la quantità 1 −
b(x)
b̂(x)
b(x)
b̂(x)
≤ ξ2
(133)
.
risulterà limitata in modulo da max{ξ1 , ξ2 }. Se, ad esempio:
r
bm ≤ b ≤ bM
bM
=β>1 ,
bm
(134)
possiamo utilizzare, come stima del guadagno b:
p
bm bM
b̂ =
(135)
.
In tal modo, risulta:
bm
b̂
r
=
b
√ m
=
bm bM
r
=
b
√ M
=
bm bM
bm
1
=
bM
β
.
(136)
.
(137)
analogamente:
bM
b̂
bM
=β
bm
Come è possibile verificare, risulta:
1−
bM
b̂
≤1−
b
b̂
≤1−
bm
b̂
1
1−β ≤1− ≤1−
β
b̂
β−1
b
1−β ≤1− ≤
β
b̂
b β − 1 ≤ 1 − ≤ β − 1
β
b̂
b
(138)
La condizione di sliding risulta:
"
ṡs =
=
#
b(x) ˆ
f (x) +
−f (x) + ẍd − λx̃˙ − k sign(s) − ẍd + λx̃˙ s
b̂(x)
!
b
b(x)
f (x) − ẍd + λx̃˙ s − k sign(s)s ≤ −η|s|
1−
b̂(x)
b̂
La (139) può essere modificata come segue:
b
b
− k |s| ≤ −η|s| − 1 −
f − ẍd + λx̃˙ s
b̂
b̂
b̂
b
˙
k|s| ≥
η|s| + 1 −
f − ẍd + λx̃ s
b
b̂
La (139) può essere soddisfatta (condizione sufficiente) scegliendo k come segue:
h
i
k ≥ β η|s| + (β − 1) f − ẍd + λx̃˙ (139)
(140)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
49
Modellistica e Controllo di Sistemi Meccanici
Incertezza su f e b È possibile trattare anche il caso di incertezze sia su f che su b. Assumiamo che per f e b sussistano rispettivamente le condizioni (130) e (138). Calcoliamo ṡ:
ṡ =
=
=
=
=
ẍ − ẍd + λx̃˙ =
i
bh ˆ
f+
−f + ẍd − λx̃˙ − ksign(s) − ẍd + λx̃˙ =
b̂
i
bh ˆ
f − fˆ + fˆ +
−f + ẍd − λx̃˙ − ksign(s) − ẍd + λx̃˙ =
b̂
b
b
b ˆ
ˆ
f+
− 1 ẍd − λx̃˙ − ksign(s) =
f −f + 1−
b̂
b̂
b̂
h
i b
b
f − fˆ + 1 −
fˆ − ẍd + λx̃˙ − ksign(s) .
b̂
b̂
(141)
dalla (141) risulta che per il soddisfacimento della condizione di scivolamento deve essere:
i
b
b hˆ
f − ẍd + λx̃˙ s − ksign(s)s ≤ −η|s|
(142)
f − fˆ s + 1 −
b̂
b̂
da cui:
b
b̂
i
b hˆ
ˆ
≥ η|s| + f − f s + 1 −
f − ẍd + λx̃˙ s
b̂
k|s|
che sarà sicuramente soddisfatta se:
"
k
≥
β η + F + (β − 1) fˆ − ẍd + λx̃˙ (143)
,
#
(144)
.
Chattering La parte discontinua del segnale di controllo (128) fa sı̀ che, una volta raggiunta
la superficie di sliding, venga richiesta agli attuatori una continua ed istantanea commutazione tra
± kb̂ . Questo fenomeno, denominato chattering, oltre ad essere praticamente irrealizzabile, sottopone gli attuatori a un regime di solito non tollerabile. L’effetto del ritardo di commutazione degli
attuatori sull’evoluzione dello stato è caratterizzato dalla permanenza dello stato nelle vicinanze
della sliding surface, in una zona definita boundary layer.
3.2
3.2.1
Controllo sliding mode di sistemi multivariabile
Variabili di stato ed ingressi
Supponiamo di avere un sistema la cui configurazione sia definita da un vettore di m coordinate
q = [q1 , q2 , . . . , qm ]T ed esista un vettore di m ingressi u = [u1 , u2 , . . . , um ]T . Un sistema di questo
tipo si dice “quadrato.” Ipotizziamo che lo stato del sistema sia descrivibile con un vettore x in cui
compaia ognuna delle qi , insieme alle sue derivate fino all’ordine ni − 1:
(n1 −1)
x = [q1 , q̇1 , . . . , q1
(n2 −1)
, q2 , q̇2 . . . . , q2
(nm −1) T
, . . . qm , q̇m , . . . , qm
]
.
(145)
Definiamo inoltre un vettore q (n) in cui compaiano le derivate ni -esime delle qi :
(n1 )
q (n) = [q1
3.2.2
(n2 )
, q2
(nm ) T
, . . . , qm
]
(146)
Modello dinamico
Il sistema è governato dalle equazioni:
q (n) = f (x) + B(x)u
(147)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
3.2.3
Modellistica e Controllo di Sistemi Meccanici
50
Modelli di incertezza
Facciamo l’ipotesi che B(x) sia invertibile e che anche la sua stima B̂(x) sia invertibile. Supponiamo inoltre che:
B(x) = I + ∆(x) B̂(x) ,
(148)
dove ∆(x) è una matrice di incertezza tale che:
|∆ij | ≤ Dij
, ∀i = 1, . . . , m,
∀j = 1, . . . , m
(149)
con Dij ≥ 0. Notiamo che D è una matrice di numeri non negativi.
Esaminiamo la struttura di D. Il teorema di Frobenius-Perron ci assicura che una matrice con
elementi non negativi ha un autovalore reale e non negativo σF P (detto autovalore di FrobeniusPerron) che rappresenta il raggio spettrale della matrice D:
|σ| < σF P
∀ σ = σ(D),
σ 6= σF P
.
(150)
Al fine di contenere il livello di incertezza della stima di B(x), facciamo l’ipotesi che l’autovalore
di Frobenius-Perron σF P della matrice D sia strettamente minore di 1:
σF P < 1
3.2.4
(151)
.
Variabili di sliding
Definiamo ∀i = 1, . . . , m una variabile di sliding nel modo consueto, già visto per il sistemi SISO:
si = (
d
+ λ)(ni −1) q̃i
dt
,
(152)
dove q̃i = qi − qdi è l’errore tra qi (t) ed il suo valore desiderato qdi (t), e λ > 0 è il parametro
che definisce la velocità di convergenza a zero dell’errore una volta che il sistema si trovi sulla
superficie di sliding.
Esempio 3.1: Se il sistema è un robot a m DOF, allora ni = n = 2 ∀i = 1, . . . , m, e le variabili di sliding sono definite
da:
d
(153)
si = ( + λ)(2−1) q̃i = q̃˙i + λq̃i .
dt
3.2.5
Superficie di sliding
La superficie di sliding, come si può intuire, è definita dall’uguaglianza a zero di tutte le variabili di
sliding si . In sintesi, definendo il vettore delle variabili di sliding:
s = [s1 , s2 , . . . , sm ]T
,
(154)
la superficie di sliding è definita dall’equazione:
s=0 .
3.2.6
(155)
Progettazione del controllo
Come visto nel caso SISO, il controllo sarà costituito da una parte basata sul modello û ed una
parte ur discontinua e destinata ad irrobustire il controllo, annullando l’effetto delle incertezze del
modello:
u = û + ur .
(156)
Il contributo di controllo basato sul modello deve essere tale da mantenere il sistema sulla superficie di sliding una volta che esso vi si trovi (ṡ = 0) in assenza di incertezze nel modello
(B̂ = B, fˆ = f ).
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
3.2.7
51
Modellistica e Controllo di Sistemi Meccanici
Controllo basato sul modello
3.2.7.1 Definizione delle “variabili di riferimento” Per comodità nell’esposizione sucessiva,
definiamo ∀i = 1, . . . , m, una “variabile di riferimento” qri nel modo che vedremo. Sviluppando la
potenza del binomio che compare nell’Eq. (152) otteniamo:
si
=
=
(ni −1)
(ni −2)
+ · · · + (ni − 1)λ(ni −2) q̃˙i + λ(ni −1) q̃i =
(n −1)
(n −1)
(n −2)
qi i
−q i
+ (ni − 1)λq̃i i
+ · · · + (ni − 1)λ(ni −2) q̃˙i + λ(ni −1) q̃i =
{z
}
| di
q̃i
+ (ni − 1)λq̃i
(n −1)
−qri i
=
(ni −1)
qi
(n −1)
− qri i
dove è stata definita la variabile
(n −1)
qri i
(n −1)
= qdi i
(n )
qri i
(ni −2)
− (ni − 1)λq̃i
(n −1)
La derivata temporale di qri i
(157)
,
(n −1)
qri i :
− · · · − (ni − 1)λ(ni −2) q̃˙i − λ(ni −1) q̃i
.
(158)
, che utilizzeremo in appresso per compattezza di notazione, vale:
d (ni −1)
q
=
dt ri
(n )
(n −1)
= qdi i − (ni − 1)λq̃i i
− · · · − (ni − 1)λ(ni −2) q̃¨i − λ(ni −1) q̃˙i
=
.
(159)
3.2.7.2 Condizione di permanenza sulla superficie di sliding Affinché il sistema permanga
sulla superficie di sliding, una volta che vi sia giunto, per la generica variabile di sliding si deve
essere verificata la condizione:
(n )
(n )
ṡi = qi i − qri i = 0 .
(160)
Se consideriamo l’intero vettore delle variabili di sliding, potremo scrivere:
ṡ = q (n) − q (n)
=0 ,
r
(n)
dove, in analogia alla (146), è stato definito il vettore q r
q (n)
r
=
(161)
come segue:
(n ) (n )
(nm ) T
[qr1 1 , qr2 2 , . . . , qrm
]
.
(162)
Sostituendo nella (161) l’espressione di q̇ (n) ricavata dalla (147), otteniamo:
ṡ = q (n) − q (n)
=
r
= f (x) + B(x)u − q (n)
r
.
(163)
Eguagliando a zero il secondo membro della (163) e risolvendo in u, ricaviamo la seguente
espressione (in cui è stata omessa per brevità la dipendenza funzionale di B e f da x):
u = B −1 q (n)
−
f
.
(164)
r
La (164) rappresenta il controllo che, conoscendo esattamente il modello, consentirebbe al sistema di rimanere sulla superficie di sliding. In presenza di incertezza, la (164) ci permette di
calcolare la parte del segnale di controllo basato sul modello:
ˆ
û = B̂ −1 q (n)
−
f
.
(165)
r
3.2.8
Controllo robusto (parte discontinua)
La parte discontinua ur del controllo, in analogia con quanto fatto per i sistemi SISO, si costruisce
nel seguente modo:
ur
= −B̂ −1 Ksign(s) ,
(166)
dove K = diag{ki } è una matrice diagonale di guadagni non negativi, s è il vettore delle variabili
di sliding definito nella (154) e sign(•) è un vettore di funzioni segno.
Il controllo risulta complessivamente:
ˆ
u = B̂ −1 q (n)
−
f
−
Ksign(s)
.
(167)
r
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
3.2.9
52
Modellistica e Controllo di Sistemi Meccanici
Verifica delle condizioni di sliding
Cerchiamo adesso di determinare i guadagni ki che permettono di verificare la condizione di
sliding per ogni singola variabile di sliding si .
Utilizzando la (167) Ricaviamo l’espressione della derivata del vettore delle variabili di sliding:
ṡ = q (n) − qr(n) =
= f + Bu − qr(n) =
−1
(n)
(n)
= f + I + ∆ B̂
| B̂
{z } qr − f̂ − Ksign(s) − qr =
Im×m
f+
− f̂ − Ksign(s) + ∆ qr(n) − f̂ − ∆Ksign(s) − qr(n) =
=
f − f̂ − I + ∆ Ksign(s) + ∆ qr(n) − f̂
.
qr(n)
=
(168)
Dalla (168) ricaviamo adesso l’espressione di ṡi :
ṡi
=
=
(n )
(n )
qi i − qri i =
m
m
X
X
(n)
fi − fˆi − ki sign(si ) −
∆ij kj sign(sj ) +
∆ij qrj − fˆj
j=1
.
(169)
j=1
Il rispetto delle condizioni di sliding per ogni si richiede che:
ṡi si ≤ −ηi |si | ,
(170)
ovvero, utilizzando la (169),
m
m
X
X
(n)
fi − fˆi si − ki sign(si )si −si
∆ij kj sign(sj ) + si
∆ij qrj − fˆj
| {z }
j=1
|si |
≤
−ηi |si |
j=1
(171)
m
X
(n )
ηi + sign(si ) fi − fˆi + sign(si )
∆ij qrj i − fˆj
≤
ki + sign(si )
j=1
m
X
∆ij kj sign(sj )
j=1
(172)
.
Dalla (172), si deduce che dei valori adeguati di ki possono essere ricavati se è soddisfatta la
seguente relazione:
ki
≥
m
X
m
X
(n )
Dij kj + ηi + fi − fˆi +
Dij qrj i − fˆj j=1
,
(173)
j=1
in particolare col segno uguale. Possiamo sintetizzare le (173) con il segno uguale come segue:
I −D k = z ,
(174)
dove sono stati implicitamente definiti il vettore dei guadagni k = [k1 , k2 , . . . , km ]T ed il vettore di
numeri positivi z = [z1 , z2 , . . . , zm ]T .
Il teorema di Frobenius-Perron ci assicura che una matrice non negativa D ha un autovalore
reale positivo σF P di molteplicità 1 che costituisce anche il suo raggio spettrale (quindi tutti gli
autovalori di D hanno modulo inferiore a σF P ). Se, nel modello d’incertezza, ammettiamo che
σF P < 1 per la matrice D, allora è possibile dimostrare che la matrice discriminante del sistema
lineare (174), ovvero I −D ha un’inversa che è positiva nel senso “element-wise,” ovvero elemento
per elemento:
−1
(I − D)
> 0 ∀i = 1, . . . , m j = 1, . . . , m .
(175)
ij
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
53
Teorema 1 Sia D una matrice positiva con autovalore di Frobenius-Perron σF P e sia λ > σF P un
numero positivo. La serie geometrica:
1
1
1
1
I + 2 D + 3 D2 + 4 D3 + . . .
λ
λ
λ
λ
(176)
converge a (λI − D)−1 se |λ| è maggiore del massimo autovalore di D, cioè σF P . In particolare
ciò è vero se λ = 1 e σF P < 1 e, in questo caso, risulta:
(I − D)−1
= I + D + D2 + D3 + · · · =
∞
X
Di
(177)
i=0
(PROVARE PER CREDERE: )
(I − D)(I + D + D2 + D3 + . . . )
(I − D) + (D − D2 ) + (D2 − D3 ) + · · · = I
=
(178)
Tutti gli elementi della serie sono ovviamente matrici positive (elemento per elemento) e quindi
si può concludere che anche la soluzione del sistema (174) è fatta di elementi tutti positivi.
3.3
Un altro approccio al progetto del controllore SMC
È possibile progettare un controllore che assicura il raggiungimento della superficie di scivolamente anche con ipotesi meno gravose di quelle fatte nel paragrafo 3.2. Invece di trovare dei
guadagni che impongano il rispetto delle singole condizioni di scivolamento
si ṡi
≤
−η|si |
∀i = 1, . . . , m
(179)
,
si impone il rispetto di una condizione di scivolamento “cumulativa”:
sT ṡ
≤
−ηksk
(180)
.
Dimostriamo adesso che la (180) assicura il raggiungimento della superficie di scivolamente in un
tempo finito a partire da una condizione iniziale:
s(0)
T
= s0 = [s10 . . . si0 . . . sm0 ]
.
(181)
La (180) si può scrivere anche come segue:
1 dksk2
2 dt
≤
−ηksk
(182)
,
ovvero:
ηdt ≤
−
kskdksk
= dksk
ksk
(183)
.
Supponiamo che esista un tempo t = tf di raggiungimento della superficie di scivolamento. La
(183), integrata tra le condizioni iniziali e quelle di raggiungimento della superficie di scivolamento
fornisce:
Z tf
ηtf ≤ −
dksk = − 0 − ks0 k = ks0 k ,
(184)
0
da cui segue:
tf
≤
ks0 k
η
(185)
,
Invece del controllo (167) utilizziamo il seguente controllo:
ˆ
u = B̂ −1 q (n)
−
f
−
ksign(s)
r
,
(186)
BOZZA 10 dicembre 2010
Allotta, Fioravanti, Malvezzi
Modellistica e Controllo di Sistemi Meccanici
54
Figura 18: Variabile di scivolamento
formalmente identica alla (167) dove k però è uno scalare al contrario del K che compare nella
(167).
Imponiamo adesso il rispetto della condizione di scivolamento (182). Calcoliamo la derivata
temporale di s:
ṡ =
=
=
q (n) − qr(n) =
f + Bu − qr(n) =
−1
(n)
f + I + ∆ B̂
B̂
q
−
f̂
−
ksign(s)
− qr(n) =
| {z } r
Im×m
=
=
f + qr(n) − f̂ − ksign(s) + ∆ qr(n) − f̂ − k∆sign(s) − qr(n) =
f − f̂ − k sign(s) − k ∆ sign(s) + ∆ qr(n) − f̂
.
(187)
Il rispetto della condizione di scivolamento (180) si ottiene se:
sT ṡ = sT f − f̂ − k sT sign(s) − k sT ∆ sign(s) + sT ∆ qr(n) − f̂
h
i
= −k ksk1 − k sT ∆ sign(s) + sT f − f̂ + ∆ qr(n) − f̂
≤
−ηksk
(188)
.
Dalla (188) si può derivare un valore di k che consente di rispettare la (188):
h
i
k ksk1 ≥ −k sT ∆ sign(s) + sT f − f̂ + ∆ qr(n) − f̂ + ηksk
k ksk1 ≥ k ksk k∆k + ksk kF k + k∆kkqr(n) − f̂ k + ηksk
k ksk1 ≥ ksk k k∆k + kF k + k∆kkqr(n) − f̂ k + η
.
(189)
Dato che ksk ≤ ksk1 , la (189) è sicuramente verificata se:
k ksk1 ≥ ksk1 k k∆k + kF k + k∆kkqr(n) − f̂ k + η
k ≥ k k∆k + kF k + k∆kkqr(n) − f̂ k + η
k 1 − k∆k
≥ kF k + k∆kkqr(n) − f̂ k + η
k
≥
kF k + k∆kkqr(n) − f̂ k + η
1 − k∆k
.
(190)
Allotta, Fioravanti, Malvezzi
BOZZA 10 dicembre 2010
Modellistica e Controllo di Sistemi Meccanici
Affinché sia possibile determinare un guadagno scalare positivo k è sufficiente che k∆k < 1.
55