Università degli Studi di Siena
Facoltà di Ingegneria
Corso di Laurea specialistica in Ingegneria Informatica
Inseguimento ciclico di sistemi multi robot
con vincoli di connettività.
Relatore
Prof. Domenico Prattichizzo
Correlatore
Ing. Fabio Morbidi
Tesi di laurea di
Luca Lucherini
A. A. 2007/2008
2
Indice
1 Introduzione
1.0.1 Struttura della tesi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Controllo di una squadra di robot . . . . . . . . . . . . . . . . . . . . . . . . . .
2 Invarianza controllata decentralizzata
2.1 Introduzione all’approccio geometrico e definizioni
2.1.1 Note sull’algebra dei sottospazi . . . . . . .
2.2 Il problema dell’invarianza decentralizzata . . . . .
2.3 L’approccio geometrico . . . . . . . . . . . . . . . .
2.4 L’approccio algebrico . . . . . . . . . . . . . . . . .
2.4.1 Caso con soluzione non unica . . . . . . . .
2.4.2 Studio di stabilità della Friend Matrix . . .
3 Controllo decentralizzato in retroazione
3.1 Consensus . . . . . . . . . . . . . . . . . . . . . .
3.2 Sistema multi-agente: Inseguimento ciclico . . . .
3.2.1 Controllo retroazionato per mantenimento
3.2.2 Modifiche al problema di ottimizzazione .
3.3 Risultati delle simulazioni e scelte progettuali . .
3.3.1 Strumenti per le simulazioni . . . . . . . .
3.3.2 Risultati del problema di ottimizzazione .
3.3.3 Definizione dei parametri liberi c e d . . .
3.4 Conclusioni . . . . . . . . . . . . . . . . . . . . .
3
5
6
6
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
9
9
9
11
11
13
14
14
. . . . . . . . .
. . . . . . . . .
di connettività
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
17
17
19
25
29
31
31
35
37
41
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4
INDICE
Ringraziamenti
Durante il periodo in cui ho lavorato su questo progetto molte sono le persone che in svariati
modi mi hanno aiutato, alcune coscientemente, altre senza saperlo, un ringraziamento va
quindi a tutti coloro che mi hanno sostenuto in un modo o nell’altro.
In particolare vorrei ringraziare il Prof. Domenico Prattichizzo per avermi dato la possibilità di lavorare su questo progetto e l’Ing. Fabio Morbidi che mi ha aiutato durante una
consistente parte del mio lavoro. Senza queste due persone la realizzazione del qui presente
lavoro non sarebbe stata possibile.
Ultimi, ma più importanti di tutti, ringrazio i miei genitori che mi hanno permesso di
raggiungere questo obiettivo, credendo sempre nelle mie potenzialità e sostenendomi durante
tutti gli anni del mio corso di studi. Senza di loro non sarei potuto arrivare dove sono, grazie.
Capitolo 1
Introduzione
Negli ultimi anni abbiamo testimoniato un crescente interesse nella coordinazione e nel controllo dei sistemi multi agente. La ricerca in quest’area è stata stimolata dalle recenti migliorie
tecnologiche introdotte nell’ambito delle comunicazioni wireless e dalla constatazione che una
squadra di agenti può svolgere delle operazioni in maniera più efficace ed affidabile rispetto
ad un singolo robot. In particolare gli algoritmi di consenso hanno ricevuto una notevole
attenzione in questo campo. I Consensus problems hanno una lunga storia nell’informatica e
vanno di fatto a formare quello che potremmo chiamare il “fondamento” della computazione
distribuita. Lo studio formale di questi problemi trova le sue origini nel 1960 ad opera di
DeGroot. La sua idea di statistical consensus theory riappare due decenni dopo, in un lavoro
di aggregazione di informazioni con incertezza ottenute da sensori multipli. Lo sviluppo da
allora è stato notevole e distribuito su molti campi tra cui uno dei più noti è senza dubbio la
computazione distribuita su reti.
Negli stessi anni in cui si sviluppavano le prime idee sui problemi di consenso un altro
problema stava venendo intensivamente studiato ovvero la sintesi di regolatori che risultassero
invarianti rispetto ai disturbi. L’impostazione classica dello studio dei controlli invarianti si
basa sulla rappresentazione matematica dei sistemi dinamici mediante matrici di funzioni di
trasferimento e i criteri usualmente impiegati per giudicare della realizzabilità dell’invarianza
fanno riferimento a particolari proprietà di tali matrici. In accordo alla tendenza del tempo di
abbandonare nello studio dei sistemi dinamici complessi i metodi trasformazionali, che mal si
prestano all’analisi mediante dispositivi di calcolo automatici, si passò all’analisi nello spazio
degli stati ottenendo delle soluzioni al problema basate sull’approccio geometrico ovvero sullo
studio di sottospazi invarianti. Tra i pionieri di questo approccio troviamo Basile e Marro.
In questa tesi verrà preso in considerazione il framework dell’inseguimento ciclico, questo
problema ha già ricevuto una notevole attenzione nella letteratura in particolare nel caso con
dinamiche a singolo integratore. Più precisamente il caso qui presentato è quello di un sistema
multi agente con un comportamento globale coordinato che si basa su una semplice regola
locale. Il vantaggio dell’uso di una struttura di agenti decentralizzata (ovvero senza leader)
sono numerosi.
Nella letteratura preesistente la comunicazione tra i veicoli è sempre considerata possibile.
Questo però non è sempre vero in quanto i robot reali hanno un raggio di connettività limitato dovuto alle tecniche di propagazione del segnale. Si cercherà quindi una soluzione per
mantenere i veicoli sempre connessi, questa sarà attuata sottoponendo il sistema a limitazioni
sul controllo e sullo stato. La determinazione di un controllore che rispetti queste limitazioni
5
6
CAPITOLO 1. INTRODUZIONE
è in genere risolta risolvendo un problema di controllo non vincolato associato al sistema e
poi modificandone la soluzione affinchè i vincoli su stato e controllo vengano soddisfatti. Genericamente queste modifiche consistono nella saturazione del vettore di controllo, risultando
però in un sistema ad anello chiuso la cui stabilità non è sempre garantita. Degli approcci
più efficienti per la soluzione del problema in questione sono basati sulle proprietà degli insiemi positively invariant. Si definirà quindi un controllo in anello chiuso che mantenga la
connessione tra i veicoli sfruttando la teoria sugli invarianti.
Verrà inoltre analizzata la retroazione di un sistema e si cercherà di vincolare i controlli in
retroazione a dipendere da elementi specifici del vettore di stato. Per risolvere questo problema
si sfruttano elementi dalla teoria sugli invarianti, verrà introdotto il concetto di invarianza
controllata decentralizzata per sostituire un controllo centralizzato. Verranno analizzate due
soluzioni a questo problema, una geometrica ed una algebrica.
1.0.1
Struttura della tesi
La tesi è divisa in due parti principali:
• Nella prima parte viene trattato il problema dell’invarianza controllata decentralizzata
facendo un richiamo sulla teoria degli invarianti e passando poi a descrivere i due approcci
studiati per la soluzione del problema. Viene poi analizzato il sistema in retroazione
andando a vedere se questo è stabile.
• Nella seconda parte viene introdotto il problema del consenso ed il framework di lavoro.
Una volta mostrati i comportamenti ad anello aperto del sistema viene introdotto il
problema della comunicazione a raggio limitato. Si illustrano i vincoli che il sistema deve
rispettare e si và a definire il controllo in retroazione sfruttando la teoria sugli insiemi
positively invariant. Vengono poi mostrati i risultati delle simulazioni sul sistema ad
anello chiuso e si mostra come varia il comportamento del sistema al variare di alcuni
parametri.
1.1
Controllo di una squadra di robot
Parlando di robot la prima cosa che ci viene in mente sono senza dubbio quelle macchine con
sembianze umanoidi e volontà proprie a cui l’immaginario comune ci ha abituato a pensare.
Nonostante la realtà sia ancora distante da questa idea i progressi fatti negli ultimi anni sono
notevoli. Tra i tipi di robot attualmente impiegati nella ricerca troviamo:
• swarm robots;
• haptic interface robots;
• soft robot ed altri.
Il punto chiave nell’utilizzo di questi dispositivi è il loro controllo. Questo può essere attuato
su un solo robot oppure su una squadra di robot. Tra gli ambiti di ricerca che coinvolgono
squadre di robot troviamo:
• esplorazione di ambienti,
• pianificazione di percorsi,
1.1. CONTROLLO DI UNA SQUADRA DI ROBOT
7
• giochi di squadra.
Prendiamo ad esempio l’esplorazione di ambienti, il nostro controllo in questo caso sarà un
controllo di formazione. Con questo termine intendiamo il problema di controllare la posizione relativa e l’orientamento dei robot in un gruppo, permettendo al gruppo di muoversi in
formazione. Un approccio possibile per questo problema può essere quello del leader-follower,
in questo caso un robot della formazione, il leader, si muove seguendo una traiettoria predefinita, mentre gli altri robot, gli inseguitori (o follower), mantengono una postura desiderata
(distanza ed orientamento) dal leader. Il maggiore problema di questo approccio risiede nel
fatto che il ruolo del leader è essenziale per la formazione e questo và quindi ad influenzare
l’affidabilità del sistema. Va però detto che l’approccio presenta una particolare flessibilità e
modularità, particolarmente gradite nel caso in cui nuovi follower si uniscono alla formazione.
Questo tipo di approccio presenta una struttura degli agenti definita struttura centralizzata.
Nella prossima sezione useremo l’approccio geometrico per passare ad una struttura decentralizzata. I vantaggi che si ottengono sfruttando questo tipo di struttura sono molteplici:
• il sistema fornisce robustezza al fallimento di singoli agenti;
• il sistema fornisce maggiore modularità e flessibilità di una struttura tipica;
• la reattività e l’affidabilità del sistema sono maggiori.
8
CAPITOLO 1. INTRODUZIONE
Capitolo 2
Invarianza controllata decentralizzata
2.1
Introduzione all’approccio geometrico e definizioni
L’essenza dell’approccio geometrico consiste nello sviluppo di un supporto matematico in una
forma libera da coordinate, per poter così sfruttare risultati semplici ed eleganti, che facilitano
una comprensione del significato di istruzioni e procedure. L’elemento principale di questo
approccio è il concetto di invarianza di un sottospazio rispetto ad una trasformazione lineare.
2.1.1
Note sull’algebra dei sottospazi
Ci si riferirà di seguito a sottospazi X, Y, Z degli spazi vettoriali a dimensioni finite con prodotto interno Rn , Rm (C n , C m ) e si indicherà con A sia una matrice m × n reale (complessa)
sia la corrispondente trasformazione lineare da Rn a Rm (da C n a C m ). Le operazioni base
sui sottospazi sono:
1. Somma:
Z = X + Y = {z : z = x + y, x ∈ X, y ∈ Y }
(2.1)
2. Trasformazione lineare:
Y = AX = {y : y = Ax, x ∈ X}
(2.2)
3. Complementazione ortogonale:
Y = X ⊥ = {y : hy, xi = 0, x ∈ X}
(2.3)
Z = X ∩ Y = {z : z ∈ X, z ∈ Y }
(2.4)
4. Intersezione
5. Trasformazione lineare inversa
X = A−1 Y = {x : y = Ax, y ∈ Y }
(2.5)
Consideriamo adesso il seguente sistema lineare stazionario:
x(t)
˙
= Ax(t) + Bu(t)
(2.6)
y(t) = Cx(t)
(2.7)
9
10
CAPITOLO 2. INVARIANZA CONTROLLATA DECENTRALIZZATA
in cui u ∈ Rp , y ∈ Rq , x ∈ X = Rn e A, B, C indicano matrici reali di dimensioni opportune.
Nell’applicazione dell’approccio geometrico allo studio della struttura dei sistemi lineari stazionari risulta fondamentale poter stabilire se esistono traiettorie appartenenti a dati sottoinsiemi
dello spazio degli stati, specificatamente a sottospazi.
Lemma 1 (lemma fondamentale dell’approccio geometrico) Condizione necessaria e
sufficiente affinchè una generica traiettoria x|[t0 ,t1 ] del sistema 2.6 appartenga ad un dato
sottospazio L ⊆ X è che sia x(t)
˙
∈ L quasi ovunque in [t0 , t1 ]
Introduco adesso la definizione di sottospazio invariante.
Definizione 1 Sia V uno spazio vettoriale e A : V → V una funzione lineare. Un sottospazio
X ⊆ V si dice invariante in A o A-invariante se:
AX ⊆ X
Si consideri il sistema libero
x(t)
˙
= Ax(t)
(2.8)
vale la seguente proprietà.
Proprietà 1 Un sottospazio L ⊆ X è luogo di traiettorie del sistema 2.8 se e solo se è
invariante in A.
Il concetto di invariante risulta essenziale in rapporto alla controllabilità e osservabilità del
sistema a 3 mappe (2.6, 2.7). Con riferimento a questo sistema si indicano con Rt+1 (x) e Rt−1 (x)
rispettivamente gli insiemi degli stati raggiungibili dallo stato x nell’intervallo [0, t1 ] e degli
stati controllabili allo stato x nello stesso intervallo. Si indicano invece con ε−
t1 (u(·), y(·)) e
ε+
(u(·),
y(·))
rispettivamente
gli
insiemi
degli
stati
iniziali
e
finali
compatibili
con
la funzione
t1
di ingresso u(·) e quella di uscita y(·) nell’intervallo [0, t1 ] detti anche insieme degli stati
inosservabili e degli stati non ricostruibili in [0, t1 ] in relazione alle date funzioni di ingresso e
di uscita. Si citano adesso le seguenti proprietà:
Proprietà 2 Nel caso di sistemi lineari stazionari continui il sottospazio Rt+1 (0), degli stati
raggiungibili dall’origine nell’intervallo [0, t1 ] con t1 > 0, è min I(A, B), il minimo invariante
in A contenente il sottospazio delle azioni forzanti B.
Proprietà 3 Nel caso dei sistemi lineari stazionari continui il sottospazio ε−
t1 (0, 0), degli stati
inosservabili nell’intervallo [0, t1 ] con t1 > 0, è max I(A, C), il massimo invariante in A
contenuto nel sottospazio degli stati inaccessibili C.
Quanto finora visto è valido in assenza di controllo, nel caso in cui quest’ultimo sia presente
possiamo introdurre il concetto di invariante controllato.
Definizione 2 Dati una trasformazione lineare A : X → X e un sottospazio B ⊆ X, un
sottospazio V ⊆ X si dice invariante controllato in (A, B) se vale la relazione
AV ⊆ V + B
Si vuole adesso vedere qual’è il collegamento tra gli invarianti controllati ed i sistemi dinamici,
ci viene in aiuto la seguente proprietà:
2.2. IL PROBLEMA DELL’INVARIANZA DECENTRALIZZATA
11
Proprietà 4 Un sottospazio V ⊆ X è luogo di traiettorie del sistema 2.6 se e solo se è
invariante controllato in (A, B)
Citiamo adesso il seguente teorema
Teorema 1 Un sottospazio V ⊆ X è invariante controllato in (A, B) se e solo se esiste una
matrice F tale che (A + BF )V ⊆ V.
La F sopra definita è chiamata Friend Matrix. Sfruttando la F si ha che u = F x e quindi la
2.6 diventa:
x(t)
˙
= Ax(t) + BF x(t)
2.2
Il problema dell’invarianza decentralizzata
In genere la matrice F dipende da uno o più degli elementi del vettore di stato, noi vorremmo
vincolare F ad avere dipendenza da elementi precisi del vettore di stato, in modo da poter
decidere su quali elementi del vettore di stato agiscono gli ingressi.
Esempio
Consideriamo il caso in cui il vettore di stato x è composto da quattro elementi x = (x1 x2 x3 x4 )T ,
e la matrice B è tale per cui ho due ingressi. Vorremmo che gli ingressi siano vincolati alle
variabili di x nel modo seguente:
u1 = f unction of (x1 , x2 )
u2 = f unction of (x3 , x4 )
Per avere queste dipendenze la Friend Matrix deve avere la seguente struttura:


x1
 x2 
F1,1 F1,2
0
0
u1


=
u2
0
0
F2,1 F2,2  x3 
x4
Sono stati analizzati due metodi per ottenere il risultato desiderato:
• Approccio geometrico;
• Approccio algebrico.
2.3
L’approccio geometrico
In questo approccio si cerca di raggiungere il risultato della decentralizzazione sfruttando le
nozioni teoriche sui sottospazi invarianti. Introduco la notazione che useró in seguito applicata
al caso con vettore di stato a quattro elementi.
x1
x3
x2 =
x1 =
x2
x4
x1 (k + 1)
A11 A12
x1 (k)
B1
u1 (k)
=
+
x2 (k + 1)
x2 (k)
B2
A21 A22
u2 (k)
12
CAPITOLO 2. INVARIANZA CONTROLLATA DECENTRALIZZATA
In questo caso tutte le Aij sono matrici 2*2, cosí come B1 e B2 . Divido inoltre in due parti
di uguale dimensione il sottospazio invariante:
V1
V=
V2
Sfruttando questa notazione posso riscrivere la definizione di sottospazio invariante controllato
come segue:
A11 V1 + A12 V2 ⊆ V1 + B1
A21 V1 + A22 V2 ⊆ V2 + B2
(2.9)
(2.10)
Scrivo adesso delle condizioni, solo sufficienti, che soddisfano 2.9 e 2.10:
A11 V1 ⊆ V1 + B1
A12 V2 ⊆ V1 + B1
⇓
A11 V1 + A12 V2 ⊆ V1 + B1
A21 V1 ⊆ V2 + B2
A22 V2 ⊆ V2 + B2
⇓
A21 V1 + A22 V2 ⊆ V2 + B2
Riscrivendo in forma matriciale le precedenti equazioni ottengo:
B1
V1
A11
0
V1
+
⊆
B2
V2
V2
0 A22
B1
V1
0 A12
V1
+
⊆
B2
V2
A21
0
V2
(2.11)
(2.12)
Per ottenere la F con la struttura desiderata modifico l’equazione 2.12 nel seguente modo,
inserendo così un ulteriore vincolo di conservatività:
0 A12
V1
V1
⊆
(2.13)
A21
0
V2
V2
Con questa modifica vado a chiedere che il sottospazio sia invariante semplice, chiedo quindi
che non ci sia bisogno del controllo per rimanervi. Se questo è vero allora la F può essere nulla
dove desiderato. Un metodo per trovare il sottospazio per cui valgono le equazioni 2.11, 2.13
è il seguente:
%Definisco le matrici...
A=[1 2 1 1; 4 0 -1 0; 0 0 1 0; 0 0 1 1]
B=[1 1 0 0;0 0 1 1]’
E=[1 0 0 0; 0 0 1 0]
Adiag=[A(1:2,1:2) zeros(2,2);zeros(2,2) A(3:4,3:4)]
Aandi=[zeros(2,2) A(1:2,3:4);A(3:4,1:2) zeros(2,2)]
Q1=robcoin([Adiag Aandi],[B zeros(size(B))], [E E])
13
2.4. L’APPROCCIO ALGEBRICO
Una volta trovato il più grande sottospazio che verifica le equazioni si calcola su questo la
Friend Matrix sfruttando la seguente istruzione Matlab:
F=effe(A,B,Q1)
La matrice così trovata dovrebbe essere decentralizzata e verificare quindi le seguenti:
A11 A12
B1
F1 0
+
= A + BF
A21 A22
B2
0 F2
A11 + B1 F1
A12
V1
V1
⊆
A21
A22 + B2 F2
V2
V2
2.4
L’approccio algebrico
Si affronterà adesso lo stesso problema con un nuovo approccio che in qualche modo è “meno”
geometrico e più algebrico. Riprendendo la definizione di invariante controllato è possibile
vedere che:
AV ⊆ V + B
(2.14)
X
AV = V B
(2.15)
U
X
Se il ker[V B] = 0 allora esiste un’unica
che risolve 2.15. Esplicitando rispetto a
U
X
ottengo:
U
#
X
= V B
AV
(2.16)
U
dove # stá per pseudo-inversa.
Sapendo che:
U = −F x
x∈V
posso scrivere che U = −F V, da cui ricavo che −V T F T = U T . Esiste peró un kernel di V T ,
l’equazione da cui calcolerò F T è quindi la seguente:
F T = −(V T )# U T + ker(V T )γ
(2.17)
A questo punto si vuole vedere se tra le ∞ F ne esiste almeno una con la struttura desiderata.
In 2.17 abbiamo a nostra disposizione un termine libero (γ), il nostro problema è risolvibile
se ∃γ tale per cui le parti indesiderate di −(V T )# U T possono essere azzerate. È immediato
vedere che mentre prima lavoravamo a livello di sottospazio invariante V stavolta lavoriamo
sulla F sfruttando il sottospazio invariante trovato senza vincoli alcuni sfruttando l’istruzione
matlab:
mainco(A,B,ker(E)).
Mostreró adesso con un esempio il procedimento utilizzato nell’algoritmo di soluzione di questo
approccio supponendo di voler trovare una struttura decentralizzata come vista precedentemente.
14
CAPITOLO 2. INVARIANZA CONTROLLATA DECENTRALIZZATA
Esempio
Iniziamo osservando l’equazione 2.17 e vedendo da cosa sono composti i suoi termini, come
nel caso precedente abbiamo un sistema a 2 ingressi u1 ed u2 .
p q
(V T )# U T =
(2.18)
r s
c
T
γ1 γ2
ker(V )γ =
(2.19)
d
Entrambi i termini sopra riportati sono dimensionalmenti uguali ad F T . Osservando 2.18
è facile vedere che per avere la struttura decentralizzata i termini q ed r dovrebbero essere
azzerati. Questo è possibile solo se si verificano le seguenti condizioni:
• d è nello spazio colonna di r
• c è nello spazio colonna di q.
Se queste condizioni sono rispettate allora é possibile trovare degli opportuni moltiplicatori
tali per cui la Friend Matrix assume struttura decentralizzata.
2.4.1
Caso con soluzione non unica
Abbiamo fino a questo momento supposto che ker[V B] fosse nullo nell’equazione
2.15,
se
X
però questa condizione viene a mancare allora non esiste una soluzione unica per
. Per
U
implementare questo nuovo scenario l’equazione 2.16 viene modificata come segue:
#
X
= V B
AV + ker V B ξ
(2.20)
U
La nuova U estratta dall’equazione 2.20 viene quindi sostituita nella formula della Friend
Matrix (equazione 2.17) ottenendo cosí la seguente:
F T = −(V T )# (Upart + Uomog )T + ker(V T )γ
(2.21)
La nostra nuova equazione presenta due gradi di libertà ξ e γ che sfrutteremo ancora una volta
per azzerare i contributi indesiderati di F . L’algoritmo che risolve questo approccio (dove
possibile) sfrutta il Symbolic Math Toolbox. Grazie a questo strumento si calcola una F i cui
termini sono equazioni nelle quali compaiono delle variabili “libere”. Tramite un procedimento
iterativo si vanno a scegliere degli opportuni valori per ξ e γ in modo da azzerare le equazioni
corrispondenti a zero nella F desiderata. L’approccio algebrico presenta genericamente dei
calcoli più complessi se paragonati a quelli a dell’approccio geometrico.
2.4.2
Studio di stabilità della Friend Matrix
Abbiamo fino ad ora imposto dei vincoli solo sulla forma della matrice F, si vuole ora analizzare
le proprietá di stabilitá dell’invariante V. A partire dalla proprietá 6.5.5 riportata su [2] ci
è noto che un invariante controllato V è stabilizzabile sia internamente sia esternamente se e
solo se esiste almeno una matrice reale F tale che sia (A + BF )V ⊆ V con A + BF stabile.
15
2.4. L’APPROCCIO ALGEBRICO
Come ben noto dalla teoria sui sistemi a tempo continuo, un sistema è stabile quando tutti
i suoi autovalori hanno parte reale negativa. Usando l’approccio algebrico otteniamo una F
della forma desiderata con al suo interno uno o più termini liberi. Purtroppo la risoluzione del
calcolo degli autovalori in forma simbolica è genericamente un calcolo computazionalmente
molto complesso. Mostrerò adesso un esempio di quello su cui viene fatto il calcolo degli
autovalori.
Esempio
Supponiamo di lavorare con le seguenti matrici A, B, C.




1 0
1 2 1 1
 1 0 
 2 2 −1 1 
1 1 1 1




C=
B=
A=
0 −1 
1 1 1 1
1 0 0 0 
0 −1
0 0 1 1
Applicando l’approccio algebrico ottengo la seguente F :
−1.16666 + y20 −1.16666 + y20
0
F =
0
0
−0.3333 −
y0
2
0
−0.6666 −
y0
2
Come si vede questa è proprio della forma desiderata, però la presenza di y0 1 rende la computazione degli autovalori di A + BF non solubile dal calcolatore. Oltretutto questo non è
neppure uno dei casi peggiori che si possono presentare.
Alternative al calcolo diretto degli autovalori
Per poter ovviare alla presenza di termini liberi all’interno delle matrici, che risulta in un
problema computazione, si presentano 2 vie principali:
• Cercare di diminuire la complessitá della matrice in questione sfruttando opportuni
artifici matematici.
• Usufruire di opportuni teoremi per verificare le proprietá degli autovalori direttamente
dalla matrice A + B × F
Strada 1
Al fine di semplificare la computazione si cerca di ridurre la complessitá della matrice tramite
artefici matematici. Si può sfruttare ad esempio un cambiamento di base tramite la matrice
T = [V V ⊥ ]. Svolgendo il calcolo:
T −1 × (A + BF ) × T
nel caso in cui ker[V B] = 0 ottengo una matrice la cui prima colonna é composta da zeri,
questo perché c’é un autovalore nell’origine e posso quindi analizzare il blocco di dimensioni
n − 1 × n − 1 rimanente. Andiamo a osservare le modifiche apportate dal cambiamento di
base, possiamo riscrivere il cambiamento di base come segue:
T −1 × (A + BF ) × T ⇒ T −1 × A × T + T −1 × BF × T
1
y0 è il parametro libero a nostra disposizione in questo caso.
16
CAPITOLO 2. INVARIANZA CONTROLLATA DECENTRALIZZATA
Andando ad analizzare la moltiplicazione di F e T, sostituisco ad F l’equazione 2.17 opportunamente trasposta ottenendo cosí:
F = −U ((V T )# )T + γ T (ker(V T ))T γ
(2.22)
Tramite oppurtune considerazioni algebriche si vede che é possibile semplificare come segue:
((V T )# )T × T = [eye(size(V, 2)) zeros(rimanenticolonne)]
(ker(V T ))T × T = [zeros(size(V, 2)) eye(rimanenticolonne)]
Purtroppo questo metodo non porta a semplificazioni degne di nota, si decide quindi di vedere
se esistono teoremi matematici che possano aiutarci nell’analisi.
Strada 2
Tra i teoremi che consentono di analizzare le proprietá degli autovalori di una matrice a partire
dalla sua struttura troviamo il Teorema di Gershgorin. Questo afferma il seguente:
Teorema
P 2 Sia A una matrice di dimensioni n*n, con componenti aij . Per i ∈ {1, ...n} scrivi
Ri = i6=j |aij | dove |aij | é il valore assoluto dell’elemento in questione. Sia D(aii , Ri ) il disco
chiuso centrato in aii e di raggio Ri . Questi dischi sono chiamati Dischi di Gerschgorin.
Teorema 3 (di Gerschgorin) Ogni autovalore di A giace all’interno di almeno uno dei
dischi di Gerschgorin D(aii , Ri ).
Si puó a questo punto fare in modo che le nostre variabili simboliche assumano valori tali per
cui questo teorema é soddisfatto, se questo é possibile la stabilitá dell’invariante V è verificata.
Per soddisfare la nostra condizione sulla stabilità si devono definire i cerchi in modo che siano
interamente nella zona del piano a parte reale negativa. Purtroppo questo teorema presenta
una condizione solamente sufficiente per la stabilitá, puó quindi succedere che i cerchi non
siano interamente nella zona a parte reale negativa, ma che gli autovalori siano comunque
validi per la stabilitá. Questa è purtroppo una forte limitazione all’applicazione di questo
metodo in forma di algoritmo.
Considerazioni
Dovendo valutare i due approcci sopra illustrati possiamo dire che l’approccio algebrico trova
soluzione in un maggior numero di applicazioni, pur essendo meno elegante e compatto rispetto
all’approccio geometrico. Punto saldo nell’analisi dei due approcci è comunque il fatto che i
metodi sono risultati troppo conservativi.
Capitolo 3
Controllo decentralizzato in
retroazione
La decentralizzazione come introdotta precedentemente è principalmente una tematica teorica,
nella realtà il problema sopra introdotto deve essere modificato opportunamente. Il problema
dovrà quindi subire una generalizzazione o più che altro una caratterizzazione. Quest’ultima potrebbe consistere nell’inserimento dei legami tra i robot, o nella definizione dei loro
comportamenti, si vanno quindi ad aggiungere dei vincoli al nostro problema.
Di seguito verrà introdotto il lavoro effettuato per controllare una squadra di robot, questi
erano stati focalizzati su un obiettivo comune, il cui raggiungimento era ottenuto con l’ausilio
di un algoritmo di consenso, che provvedeva alla coordinazione dei loro movimenti.
3.1
Consensus
Nelle reti di agenti “consensus” significa raggiungere un accordo riguardante una certa quantità
di interesse che dipende dallo stato di tutti gli agenti. Un “algoritmo di consenso” (o protocollo)
è una regola di interazione che definisce lo scambio di informazioni tra un agente e tutti i suoi
“vicini” sulla rete.
Per poter lavorare con queste tematiche introduco di seguito alcuni dei concetti necessari.
La topologia di interazione di una rete di agenti è rappresentata usando un grafo direzionato
G = (V, E) in cui V={1, 2,· · · ,n} è il set di nodi ed E ⊆ V × V l’insieme degli archi. I vicini
di un agente sono denotati da Ni = {j ∈ V : (i, j) ∈ E}, sono quindi vicini due agenti i cui
rispettivi nodi sono connessi da un arco orientato da i verso j. Un semplice algoritmo per
raggiungere un accordo riguardante lo stato di n agenti integratori con dinamica x˙i = ui può
essere espresso come un sistema lineare su un grafo
x˙i (t) =
X
j∈Ni
(xj (t) − xi (t)), xi (0) = zi
questo protocollo prende il nome di “average consensus”. La dinamica collettiva del gruppo
di agenti che segue l’algoritmo può essere riscritta come:
x˙ = −Lx
17
18
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
dove L è il grafo laplaciano della rete ed i suoi elementi sono definiti come:
−1, j ∈ Ni
lij =
|Ni |, j = i
In questa formula |Ni | denota il numero di vicini del nodo i. Vista la definizione
di grafo
P
laplaciano la somma di tutte le righe di L dà zero come risultato visto che j lij = 0. Quindi
il grafo laplaciano ha sempre un autovalore nullo λ1 = 0. Questo autovalore nullo corrisponde
all’autovettore 1=(1, · · · , 1)T perchè 1 appartiene al null-space di L (L1 = 0). Si può quindi
trovare un equilibrio del sistema in uno stato della forma P
x∗ = (α, · · · , α)T = α1 in cui tutti
1
gli agenti concordano. Questo valore di consenso è α = n i zi o in altre parole la media dei
valori iniziali degli agenti.
Se si volessero analizzare delle proprietà di convergenza degli algoritmi di consenso dovremmo sfruttare le proprietà spettrali della matrice laplaciana. In accordo con il teorema di
Gerschgorin, tutti gli autovalori di L nel piano complesso sono situati in dischi chiusi centrati
in ∆ + 0j il cui raggio è ∆ = maxi di . Per grafi non orientati, L è una matrice simmetrica
con autovalori reali, il set di autovalori può essere quindi ordinato in ordine ascendente come:
0 = λ1 ≤ λ2 ≤ · · · ≤ λn ≤ 2∆
Come precedentemente citato l’autovalore nullo è noto come autovalore banale di L. Il secondo
più piccolo autovalore del laplaciano λ2 rappresenta una misura di prestazioni per l’algoritmo
di consenso ed è chiamato “algebraic connectivity” di un grafo.
Altri tipi di algoritmi di consenso
Oltre all’average consensus algorithm in letteratura troviamo molti altri algoritmi come ad
esempio il Weighted-Average Consensus. In questo protocollo viene sfruttato un vettore dei
pesi γ = (γ1 , · · · , γn ) per modificare la dinamica collettiva come segue:
K x˙ = −Lx
dove K = diag(γ1 , · · · , γn ). Questo equivale ad un nodo con rapporto variabile di integrazione
basato sul protocollo
X
γi x˙i =
aij (xj − xi )
j∈Ni
Negli esempi visti sopra gli algoritmi di consenso trovano un accordo sulla base delle condizioni
iniziali degli agenti, esistono però anche algoritmi che forniscono agli agenti degli “obiettivi”
time-varying. Un esempio di questo tipo di protocolli può essere quello studiato da Ren e
Sorensen in [6]. In questo articolo una squadra di agenti ha l’obiettivo di mantenere una certa
formazione mentre insegue un centroide tempo variante, l’elemento su cui si vuole raggiungere
un accordo è in questo caso il centro della formazione. La squadra di agenti è divisa tra leader
e follower, i primi conoscono le coordinate del centroide mentre i secondi devono ottenerle
tramite un algoritmo di consenso. In questo caso le dinamiche degli agenti sono più complesse
e dipendono dal centroide stimato da ogni veicolo ξ˙ivc . La stima per i leader viene fatta usando
la seguente:
P
vc [ξ˙vc − γ(ξ vc − ξ vc )]
ξ˙dvc − γ(ξivc − ξdvc ) + nj=1 gij
j
i
j
vc
Pn
,
ξ˙i =
vc
1 + j=1 gij
3.2. SISTEMA MULTI-AGENTE: INSEGUIMENTO CICLICO
19
mentre per i follower viene utilizzata:
Pn
gvc [ξ˙jvc − γ(ξivc − ξjvc )]
˙ξ vc = j=1 ij P
n
i
vc
j=1 gij
dove ξdvc sono le reali coordinate del centroide passate ai leader, γ è un parametro positivo ed
i termini gij definiscono se c’è comunicazione tra un veicolo e l’altro.
Questo protocollo è stato citato per mostrare le potenzialità di questo tipo di algoritmo,
queste sono tutt’altro che limitate e grazie alla loro versatilità gli algoritmi di consenso si
prestano per essere utilizzati in una vasta gamma di applicazioni.
3.2
Sistema multi-agente: Inseguimento ciclico
Andiamo adesso ad introdurre il framework di lavoro che è stato utilizzato. Il nostro sistema
multiagente si basa sul “cyclic pursuit”, ogni agente conosce quindi solo lo stato dell’agente a
lui successivo, la topologia di interazione tra gli agenti è rappresentata in Figura 3.1. In questa
struttura non sono presenti leader, il sistema di agenti può essere definito decentralizzato. I
vantaggi di questo tipo di struttura sono già stati elencati nella sezione 1.1 e non verranno
quindi nuovamente citati.
<
=
0
>
1
2
?
3
@
;4
5
6
/7
8
9
:.A
-B
+
,C
*#
D
$
%
)"E
(F
!&
G
'K
JH
L
M
IN
~
O
}P
|m
Q
u
w
v{zyxtR
sln
ko
p
q
rS
T
jU
g
h
V
ia
b
d
ce
W
f_
`^]X
\Y
Z
[
2
3
1
i
n
Figura 3.1: Topologia di interazione tra gli agenti nel cyclic pursuit
Una matrice circolante di ordine n è una matrice quadrata della forma:


c1 c2 · · ·
cn
 cn c1 · · · cn−1 


C= .
.. . .
..  .
 ..
.
.
. 
c2
c3 · · ·
c1
Gli elementi di ogni riga di C sono identici a quelli della precedente riga ma spostati di
un elemento verso destra (l’elemento che uscirebbe dalla matrice viene reinserito a sinistra).
L’intera matrice è quindi definita unicamente dalla prima riga e può essere riscritta come
C = circ[c1 , c2 , · · · , cn ].
L’utilizzo di questo tipo di matrici è comodo per poter determinare prontamente autovalori ed
autovettori della stessa (è possibile diagonalizzare C tramite una matrice di Fourier). Risultati
20
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
equivalenti sono validi anche nel caso di matrici circolanti a blocchi in cui l’intera matrice è
definita dal primo blocco di righe, una matrice di questo tipo si può riscrivere come:
A = circ[A1 , A2 , · · · , An ].
Il nostro framework è composto da n agenti mobili la cui posizione è definita ad ogni istante
t ≥ 0 da ξi (t) = [xi (t), yi (t)]T . La cinematica di ogni robot è descritta da un semplice
integratore:
ξ˙i (t) = ui (t)
Nella strategia di inseguimento ciclico classica l’ingresso è:
ui = k(ξi+1 − ξi ), k ∈ R+ .
(3.1)
Si ha quindi che l’agente i insegue il suo vicino, l’agente i + 1, lungo la linea di vista, con
velocità proporzionale al vettore che và dall’agente i all’agente i + 1. Il parametro k verrà
considerato da qui in avanti, senza perdità di generalità pari ad 1.
Inseriamo adesso una nuova ipotesi nel nostro problema, si suppone che ogni robot insegua
il suo vicino lungo una linea di vista ruotata di un angolo α ∈ [−π, π] comune a tutti gli agenti.
L’ingresso di controllo per ogni veicolo (equazione 3.1) viene quindi modificato come segue:
ui (t) = R(α)(ξi+1 (t) − ξi (t)),
dove R(α) è una matrice di rotazione così definita:
cos(α) sin(α)
R(α) =
− sin(α) cos(α)
Volendo scrivere il nostro nuovo sistema in una notazione compatta possiamo riformulare la
dinamica globale del sistema come segue:
˙ = Aξ(t),
ξ(t)
dove ξ = [ξ1T , ξ2T , · · · , ξnT ] ed A è una matrice circolante a blocchi descritta dalla seguente
equazione:
A = circ[−R(α), R(α), 02 , · · · , 02 ],
dove 02 è una matrice di zeri a dimensione 2 × 2.
Per poter conoscere i comportamenti del sistema si devono analizzare i suoi autovalori, si
sfrutta quindi la seguente formula per ottenere una diagonalizzazione a blocchi di A:
Di = A1 + ϕi A2 + ϕ2i A3 + · · · + ϕin−1 An , i = 1, 2, · · · , n.
(3.2)
Partendo dalla forma diagonale a blocchi di A è facile ricavare gli autovalori che risultano
essere una collezione di:
±jα
, i = 1, 2, · · · , n.
λ±
i = (ϕi − 1)e
Per poter esplicitare la struttura degli autovalori definisco:
i−1
±
θi =
π±α .
n
21
3.2. SISTEMA MULTI-AGENTE: INSEGUIMENTO CICLICO
Dopo aver effettuato alcune manipolazioni algebriche posso scrivere gli autovalori nel seguente
modo:
i−1
±
λi = 2 sin
π (− sin θi± + j cos θi± ).
n
Si può facilmente vedere che la matrice A ha un autovalore nullo a molteplicità 2 e che vale
la seguente proprietà:
+
λ+
i = 2, 3, · · · , n,
(3.3)
i = λn−i+2
dove λ denota il complesso coniugato di λ. Un’altro elemento necessario per lo studio dei possibili comportamenti è la conoscenza degli autovettori, si devono quindi studiare gli autovettori
della nostra matrice A i quali hanno la seguente struttura1 :
i
h
(n−1)
(n−1) T
,
wi+ = 1, j, ϕi , jϕi , · · · , ϕi
, jϕi
i
h
(n−1)
(n−1) T
,
wi− = 1, −j, ϕi , −jϕi , · · · , ϕi
, −jϕi
−
−
dove wi+ corrisponde all’autovalore λ+
i e wi corrisponde a λi .
Visto che A è diagonalizzabile la soluzione generale del sistema avrà la seguente forma:
ξ(t) =
q mλ
X
Xk
αkj etλk wkj ,
(3.4)
k=1 j=1
dove q è il numero di autovalori distinti, i parametri αkj sono costanti ed i termini wkj sono
gli autovettori corrispondenti al k-esimo autovalore. È possibile riscrivere l’equazione 3.4
estraendo i due autovettori corrispondenti all’autovalore nullo e sostituendoli con i seguenti:
1
1
= (w1+ + w1− ) = [1, 0, 1, 0, · · · , 1, 0, 1],
wG
2
1 +
2
wG = (w1 − w1− ) = [0, 1, 0, 1, · · · , 0, 1, 0],
2j
ottenendo così:
ξ(t) =
q mλ
X
Xk
1
2
+ yG wG
,
αkj etλk wkj + xG wG
(3.5)
k=2 j=1
dove xG e yG sono le coordinate del baricentro. Si possono adesso andare ad analizzare i
comportamenti del sistema al variare dell’angolo α.
Se |α| < πn tutti gli autovalori non nulli risiedono nel semipiano sinistro aperto del piano
complesso. Abbiamo quindi che:
1
2
lim ξ(t) = xG wG
+ yG wG
.
t→+∞
Da questo deriva che, per ogni condizione iniziale, tutti gli agenti convergono esponenzialmente
ad un singolo punto limite che risulta essere il baricentro, il comportamento è mostrato in
Figura 3.2 (a) e (b), come si può vedere più l’angolo è piccolo e più veloce è la convergenza.
Se α = nπ , allora oltre ai due autovalori nulli ci sono altri due autovalori non nulli che
risiedono sull’asse immaginario, mentre tutti gli altri autovalori risiedono nel semipiano sini−
stro aperto del piano complesso. Gli autovalori immaginari sono λ+
n e λ2 ed il loro valore
1
Per informazioni più dettagliate riguardo gli autovalori e gli autovettori del sistema si rimanda a [4].
22
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
è ∓j2 sin( πn ) = ∓jγ. Possiamo riscrivere l’equazione 3.5 estraendo anche i due autovettori
corrispondenti agli autovalori puramente immaginari ottenendo così:
ξ(t) =
q mλ
Xk
X
1
2
αkj etλk wkj + d1 e−jtγ wn+ + d2 ejtγ w2− + xG wG
+ yG wG
,
(3.6)
k=4 j=1
dove d1 e d2 sono due costanti. Si sostituiscono poi le due eigenfunction complesse con due
eigenfunction reali ed independenti ottenute come segue:
1 jtγ −
e w2 + e−jtγ wn+ ,
2
1 jtγ −
e w2 − e−jtγ wn+ .
=
2j
1
wdom
=
2
wdom
Otteniamo quindi, definendo δi =
2π(i−1)
n
i = 1, 2, · · · , n :
1
wdom
= [cos(δ1 + γt), sin(δ1 + γt), · · · , cos(δn + γt), sin(δn + γt)]T ,
2
wdom
= [sin(δ1 + γt), − cos(δ1 + γt), · · · , sin(δn + γt), − cos(δn + γt)]T .
Per t → +∞ si ha quindi:
1
2
2
1
+ xG wG
+ yG wG
.
+ c2 wdom
ξ(t) = c1 wdom
Dopo un transitorio, la traiettoria di ogni agente risulta essere (assumendo c1 6= 0 e c2 6= 0):
xi (t) ≃ c1 cos(δi + γt) + c2 sin(δi + γt) + xG ,
yi (t) ≃ c1 sin(δi + γt) − c2 cos(δi + γt) + yG ,
p
ovvero un moto circolare centrato nel baricentro avente raggio r = c21 + c22 e periodo 2π
γ .
Gli agenti risultano inoltre essere equispaziati visto che la distanza angolare tra l’agente i e
l’agente i + 1 è 2π
n . Il comportamento in questo caso è mostrato
in Figura 3.2(c).
L’ultimo caso rimasto da analizzare è quello con α ∈ nπ , 2π
n . Qui oltre ad avere la coppia di
−
autovalori nulli ho una coppia di autovalori complessi positivi (λ+
n e λ2 con valori β ∓jγ). Con
procedimenti analoghi al caso precedente, otteniamo che le traiettorie degli agenti risultano
essere (per c1 6= 0 e c2 6= 0):
xi (t) ≃ eβt (c1 cos(δi + γt) + c2 sin(δi + γt)) + xG ,
yi (t) ≃ eβt (c1 sin(δi + γt) − c2 cos(δi + γt)) + yG ,
ovvero una spirale logaritmica con tasso di crescita β, mostrata in Figura 3.2(d).
Comunicazione tra i veicoli
Come in tutti i problemi di consensus è necessario che i veicoli possano comunicare per poter
raggiungere un accordo, fino ad ora la comunicazione tra i veicoli era data per scontata,
andiamo adesso ad introdurre delle condizioni in presenza delle quali si ha che i veicoli sono
connessi.
23
3.2. SISTEMA MULTI-AGENTE: INSEGUIMENTO CICLICO
6
6
5.5
5.5
5
5
4.5
4.5
R1
4
y [m]
y [m]
R
2
R3
3.5
R
4
R1
R
2
4
R3
R
4
3.5
3
3
2.5
2.5
2
2
1.5
−3
−2
−1
0
x [m]
1
2
3
−2.5
−2
−1.5
−1
−0.5
(a)
0
x [m]
0.5
1
1.5
2
2.5
(b)
43
x 10
6
3
5
2
1
R1
R1
R2
y [m]
y [m]
4
R3
3
R2
0
R3
R4
R4
−1
2
−2
1
−3
−3
−2
−1
0
1
x [m]
2
3
4
5
−4
−3
−2
(c)
−1
0
x [m]
1
2
3
4
43
x 10
(d)
Figura 3.2: (a),(b) Traiettorie del sistema con |α| < πn rispettivamente con α =
(c)Traiettoria del sistema con α = nπ (d)Traiettoria del sistema per |α| > πn
<
=
0
>
1
2
?
3
@
;4
5
6
/7
8
9
:.A
-B
+
,C
*#
D
$
%
)"E
(F
!&
G
'K
JH
L
M
IN
~
O
}P
|m
Q
u
w
v{zyxtR
sln
ko
p
q
rS
T
jU
g
h
V
ia
b
d
ce
W
f_
`^]X
\Y
Z
[
veicolo i
veicolo i + 1
ri+1
Figura 3.3: Situazione di connessione tra 2 veicoli.
π
5
eα=
π
9
24
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
Come evidenziato anche in [3], il veicolo i + 1 si dice connesso al veicolo i al’istante di
tempo t se il veicolo i è all’interno del disco di comunicazione con raggio ri+1 2 in quel preciso
momento (vedi Figura 3.3).
Due veicoli (i ed i + 1) si dicono connessi per ogni istante di tempo t ≥ 0 se il veicolo i
è sempre all’interno del disco di comunicazione del veicolo i + 1 per ogni t ≥ 0. Un metodo
per mantenere la connettività potrebbe essere quello di cercare un raggio che me la preservi.
Supponiamo che il disco di comunicazione del veicolo i + 1 abbia raggio
ri+1 ≥ max ||(T [2i − 1, 2i] − T [2i + 1, 2i + 2])eAt T −1 ξ(0)||,
t≥0
dove i ∈ {1, 2, · · · , n}, la matrice di trasformazione T è
T = ℜ(w1+ ) ℑ(w1+ ) ℜ(w2+ ) · · · ℜ(wn+ ) ℑ(wn+ ) ,
ed eAt è una matrice diagonale a blocchi fatta come segue:
eAt = blkdiag[I2 , B2 (t), B3 (t), · · · , Bn (t)],
cos(ℑ(λ+
)t) sin(ℑ(λ+
)t)
ℜ(λ+
)t
k
k
e k ∈ {2, 3, · · · , n}, allora gli agenti
dove Bk (t) = e k
+
− sin(ℑ(λ+
k )t) cos(ℑ(λk )t)
raggiungono il rendezvous o convergono ad una formazione circolare equispaziata senza perdere
la connettività per ogni t ≥ 0.
Alternativamente si può provare a creare un controllo retroazionato per il sistema che
mi assicuri la connettività dei veicoli e quindi che ||ξi (t) − ξi+1 (t)|| < ri+1 . Osservando la
letteratura esistente in fatto di controllo vincolato si vede che è possibile creare dei controllori
che rispettano dei vincoli lineari sullo stato. Per poter usufruire di tali risultati si sostituisce
il disco di connettività dei robot con dei quadrati (vedi figura 3.4).
Al fine di scegliere un quadrato correlato al disco di comunicazione originale si presentano
due opzioni:
1. scegliere il più piccolo quadrato contenente il disco di comunicazione originale,
2. scegliere il più grande quadrato contenuto nel disco di comunicazione originale.
Il primo caso presenterà un’area di copertura maggiore e più precisamente di circa il 21%,
infatti:
area quadrato = 4r 2 , area triangolo = πr 2
surplusarea = (4 − π)r 2 ≈ 0.858r 2
0.858r 2
→
≈ 21% area in eccesso.
4r 2
Equivalentemente per l’altro caso abbiamo:
area triangolo = πr 2 , area quadrato interno = 2r 2
surplusarea2 = (π − 2)r 2 ≈ 1.141r 2
1.141r 2
≈ 36% area in meno.
→
πr 2
Sono state quindi mostrate le differenze di area nei due casi, qualora si volesse lavorare pensando ai dischi come riferimento, queste sono le differenze di cui si deve tenere conto. Da ora
in poi ci si riferirà al raggio pensando a metà del lato del quadrato cioè alla distanza (minima)
che separa il centro del quadrato da un suo lato.
2
Il disco di raggio ri+1 definisce l’area di trasmissione del veicolo i + 1.
25
3.2. SISTEMA MULTI-AGENTE: INSEGUIMENTO CICLICO
quadrato esterno
1
0.5
ri
0
−0.5
quadrato interno
−1
−1.5
−1
−0.5
0
0.5
1
1.5
Figura 3.4: Trasformazione del disco di comunicazione in un quadrato
3.2.1
Controllo retroazionato per mantenimento di connettività
Il nostro scopo è quello di creare un controllo retroazionato che impedisca al raggio di connettività tra i veicoli di superare il limite consentito, si vuole quindi imporre dei vincoli sul
vettore di stato. Il nostro problema rientra nella classe dei Linear Constrained Regulation
Problem. Al fine di risolvere questa tematica siamo partiti dai risultati presenti in [5]. In
questo articolo viene studiato come determinare una legge di controllo lineare u(t) = F x(t)
che deve rispettare tre tipi di vincoli:
1. vincoli sul vettore di controllo u(t);
2. vincoli sul vettore di stato;
3. vincoli sull’insieme degli stati iniziali.
Andiamo adesso ad analizzare i vincoli e come vogliamo che siano nel nostro caso di interesse.
Il vincolo sull’ampiezza dell’ingresso di controllo è del tipo −ρL ≤ u ≤ ρU , questo serve
soprattutto a dare un feeling reale al problema, mettendo delle limitazioni che rispecchiano
quelle meccaniche.
Il vincolo sul vettore di stato è quello per noi di maggiore interesse ed è del tipo Dx ≤ c
con D ∈ Rp×n , rank(D) = n, p ≥ n e c > 0. Il nostro scopo è come già detto contenere la
distanza tra i veicoli, si sfrutteranno a questo scopo dei vincoli del tipo |x1 − x2 | < r. Da
notare che nell’articolo [5] è richiesto che il vettore di stato sia vincolato ad appartenere ad
un poliedro chiuso. Per fare questo si inseriscono oltre ai vincoli necessari per il nostro scopo
anche dei vincoli ausiliari del tipo |x1 + x2 | < M , in cui M è un parametro opportunamente
grande tale per cui il problema non risulta realmente vincolato. Abbiamo così inserito dei
bound sulla distanza tra i veicoli ottenendo lo spazio di lavoro rappresentato in Figura 3.5 (è
rappresentato lo spazio che contiene una coordinata di un veicolo).
26
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
<
=
0
>
1
2
?
3
@
;4
5
6
/7
8
9
:.A
-B
+
,C
*#
D
$
%
)"E
(F
!&
G
'K
JH
L
M
IN
~
O
}P
|m
Q
u
w
v{zyxtR
sln
ko
p
q
rS
T
jU
g
h
V
ia
b
d
ce
W
f_
`^]X
\Y
Z
[
|x1 − x2 | < r
|x1 + x2 | < M
"
#
&$
$
%
!$
!#
!"
!!
!"
!#
!$
%
&'
$
#
"
!
Figura 3.5: Schematizzazione dei vincoli sul vettore di stato.
Per quanto concerne l’ultimo vincolo questo ha una forma analoga al precedente ed è
Gx ≤ w, con G ∈ Rq×n con rank(G) = n, q ≥ n e w > 0. Il nostro unico vincolo sugli stati
iniziali è che i veicoli siano in condizioni di connettività quindi i vincoli sono gli stessi che per
il vettore di stato.
Andiamo adesso ad analizzare le condizioni per l’esistenza del controllore lineare vincolato,
queste si basano sul concetto di positive invariance. Un insieme non vuoto Ω si dice insieme
positively invariant del sistema (S) se e solo se x0 ∈ Ω implica che:
x(t; x0 ; t0 ) ∈ Ω per ogni t0 ∈ T e t ≥ t0 ,
dove x(t; t0 ; x0 ) è la traiettoria di S con condizioni iniziali (t0 , x0 ). Definiamo adesso alcuni
insiemi che utilizzeremo successivamente:
Q(F, ρL , ρU ) = {x ∈ Rn : −ρL ≤ F x ≤ ρU }
è l’insieme poliedrale associato alla legge di controllo in retroazione, ovvero la regione degli
stati dove il vettore di controllo non viola il primo vincolo,
P (D, c) = {x ∈ Rn : Dx ≤ c}
è invece l’insieme degli stati che soddisfano il secondo vincolo, per la stessa notazione il set
degli stati iniziali che soddisfano il terzo vincolo è indicato da P (G, w). La legge di controllo
u = F x con F ∈ Rm×n è una soluzione dell’LCRP se e solo se:
1. Gli autovalori λi , i = 1, 2, · · · , n, della matrice A+BF hanno parte reale negativa.
2. Esiste un insieme positevely invariant Ω ⊆ Rn del sistema ad anello chiuso x˙ = (A +
BF )x tale che valgono le seguenti:
P (G, w) ⊆ Ω ⊆ P (D, c)
3.2. SISTEMA MULTI-AGENTE: INSEGUIMENTO CICLICO
27
Ω ⊆ Q(F, ρL , ρU )
I risultati sopra citati possono essere sfruttati in molti differenti algoritmi, l’approccio che noi
seguiremo si basa sul seguente corollario:
Corollario 1 Supponiamo che F sia una matrice m × n tale che:
• gli autovalori λi , i = 1, 2, · · · , n, della matrice A+BF abbiano parti reali negative3 ,
• P (G, w) sia un insieme positevely invariant del sistema ad anello chiuso,
• P (G, w) ⊆ Q(F, ρL , ρU ).
Se queste richieste sono verificate allora la legge di controllo in retroazione è una soluzione
dell’LCRP.
Al fine di soddisfare tutte le richieste sopra citate si inizia con la ricerca di condizioni
algebriche tali per cui P (G, w) è un insieme positevely invariant rispetto al sistema in anello
chiuso. È noto che un sistema asintoticamente stabile possiede insiemi hyperellipsoidal positevely invariant. Questi insiemi sono generati da funzioni di Lyapunov quadratiche xT P x,
dove P ∈ Rn×n è una funzione simmetrica definita positiva che soddisfa l’equazione matriciale
P A + AT P − Q = 0, e Q ∈ Rn×n è una matrice simmetrica definita positiva. Usando una
funzione di Lyapunov non quadratica possiamo anche generare degli insiemi positively invariant di tipo poliedrale come nel caso di P (G, w). Al fine di definire queste funzioni di tipo
Lyapunov si introducono le matrici reali n × n di tipo Mn .
Definizione 3 La matrice reale n × n B = (bij ) si dice appartenente alla classe Mn se e solo
se bij ≥ 0, per ogni i 6= j.
Visto che P (G, w) è un insieme chiuso può essere definito dall’espressione
P (G, w) = {x ∈ Rn : v(x) ≤ 1}
dove
v(x) = max
1≤i≤q
(Gx)i
wi
,
(3.7)
(Gx)i denota l’i-esima componente del vettore Gx. La seguente proposizione fornisce una
condizione necessaria e sufficiente per cui v(x) definita in 3.7 risulta essere una funzione di
Lyapunov e con essa P (G, w) un insieme positively invariant del sistema.
Proposizione 1 Il set poliedrale P (G, w) è un insieme positively invariant del sistema y˙ = Ay
se e solo se esiste una matrice H ∈ Rq×q tale che:
H ∈ Mq ,
GA − HG = 0
e Hw ≤ 0.
Una conseguenza diretta di quanto appena detto è la seguente:
Proposizione 2 Supponiamo F ∈ Rm×n e che esista una matrice H ∈ Rq×q tale che:
1. H ∈ Mq ,
3
La matrice B è stata supposta una matrice identità di dimensioni opportune, quindi la matrice F è
dimensionalmente uguale ad A.
28
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
2. Hw ≤ 0,
3. GA + GBF = HG,
4. gli autovalori λi , i = 1, 2, · · · , n, della matrice A+BF abbiano parte reale negativa,
5. P (G, w) ⊆ Q(F, ρL , ρU ).
Allora la legge di controllo u = F x è una soluzione dell’LCRP.
Alcune delle condizioni riportate sopra possono essere sostituite o modificate opportunamente4
per raggiungere il seguente risultato:
Proposizione 3 Se esiste una matrice stabile H ∈ Mq ed una matrice F ∈ Rm×n che soddisfa
le seguenti:
1. Hw ≤ 0,
2. GA + GBF = HG,
3. −ρL ≤ F xi ≤ ρU , i = 1, 2, · · · , r, 5
allora la legge di controllo u = F x è una soluzione dell’LCRP.
Come si può facilmente vedere l’applicazione del precedente metodo è facilmente attuabile
tramite il seguente problema di ottimizzazione vincolata:
funz. obiettivo: J(F, H, ε) = max ε
L
vincoli : GA + GBF = HG,
(3.8)
Hw ≤ −εw,
(3.9)
i
U
−ρ ≤ F x ≤ ρ , i = 1, 2, · · · , r,
hij ≥ 0, i 6= j,
ε ≥ 0,
(3.10)
(3.11)
(3.12)
(3.13)
dove ε è un numero reale.
Osservando il problema di ottimizzazione si vede che compare un nuovo parametro, ovvero
ε, questo viene inserito nel vincolo 3.9 per velocizzare la convergenza dell’algoritmo, infatti
maggiore sarà ε e minore sarà il valore del prodotto Hw. Al fine di massimizzare la velocità
di convergenza la funzione obiettivo viene dedicata a ε.
Risultati dopo l’applicazione del controllo trovato come soluzione dell’LCRP
Una volta risolto il problema di ottimizzazione definito nella precedente sezione abbiamo a
nostra disposizione una matrice F da utilizzare per il controllo in retroazione. Andiamo
a vedere il comportamento che assume il sistema ad anello chiuso utilizzando la suddetta
matrice.
4
5
per approfondimenti sulle modifiche si rimanda a [5]
Con xi si indicano i vertici del poliedro chiuso P (G, w)
29
3.2. SISTEMA MULTI-AGENTE: INSEGUIMENTO CICLICO
4
4
3.5
3.5
3
R1
3
R1
R2
R2
2.5
R3
y [m]
y [m]
R3
2.5
2
1.5
2
1
1.5
0.5
1
0.5
1
1.5
2
x [m]
2.5
3
3.5
0
−1
−0.5
0
(a)
0.5
1
1.5
x [m]
2
2.5
3
3.5
4
(b)
Figura 3.6: (a) Traiettoria del sistema ad anello aperto nel caso di rendezvous al baricentro;
(b) Traiettoria del sistema ad anello chiuso
Come si può vedere dalla Figura 3.6(b) si ha una convergenza diretta verso l’origine che
indubbiamente mantiene i veicoli connessi, ma al prezzo di annullare le dinamiche pre-esistenti
dei veicoli. Andando a vedere la struttura della matrice A+F si vede che questa è: −y ∗I, dove
y è un parametro positivo ed I è una matrice identità di dimensioni opportune. I problemi
causati da questa struttura sono i seguenti:
1. perdita delle dinamiche originali del sistema,
2. perdita del punto di equilibrio originale del sistema.
Si devono quindi cercare delle modifiche da apportare al problema di ottimizzazione che
abbiano un impatto più “conservativo” sulle dinamiche.
Nota: per quanto riguarda il parametro y questo risulta coincidere con il parametro ε
presente nella funzione obiettivo del mio problema di programmazione lineare.
3.2.2
Modifiche al problema di ottimizzazione
Il sistema ad anello aperto presentava delle particolari caratteristiche quali ad esempio:
• la circolarità,
• la dipendenza di ogni veicolo da un solo vicino6 .
Al fine di mantenere queste caratteristiche sono stati inseriti degli opportuni vincoli nel
problema di ottimizzazione. Questi sono riassumibili con:
F = circ[F1 , F2 , 02 , · · · , 02 ]
dove F1 , F2 ∈ R2×2 . Questi vincoli utilizzati da soli sono inutili e non risolvono nessuno dei
nostri problemi, infatti anche la struttura precedentemente trovata dal problema di ottimizzazione rispetta questi vincoli, è però importante inserirli al fine di evitare soluzioni con una
forma non desiderata.
6
Le dinamiche del veicolo i dipendono anche dalla posizione del veicolo i+1.
30
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
Analizziamo adesso il problema della convergenza al baricentro. Come visto precedentemente dalla teoria sul consenso il nostro sistema dovrebbe convergere ad un punto che è la
media delle condizioni iniziali di ogni veicolo. Le proprietà precedentemente riportate valevano per via della struttura del grafo laplaciano, tra le caratteristiche di quest’ultimo che ci
interessa replicare c’è il fatto che la somma degli elementi di ogni riga è nulla. Questa dà luogo
alla presenza dell’autovalore nullo, il quale ci permette di avere un equilibrio del sistema in un
punto del tipo β ∗(1, 1, · · · , 1). Questo andrebbe sicuramente bene nel caso monodimensionale
ma non nel nostro caso. Il punto di equilibrio che noi vogliamo raggiungere è del tipo:
x∗ = [β, γ, β, γ, · · · , β, γ]T ∈ R2 n
dove β e γ sono rispettivamente la media delle coordinate iniziali x e di quelle y. Si può
facilmente intuire che questo punto risulta essere un equilibrio del sistema se la somma di
tutti gli elementi che moltiplicano una coordinata x su una riga è 0, e se anche la somma di
tutti gli elementi che moltiplicano una coordinata y è a sua volta nulla. Il vincolo può essere
schematizzato come segue:
Pn
Fi,2j−1 = 0 Pj=1
(3.14)
n
i = 1, · · · , n,
j=1 Fi,2j = 0
dove n è il numero di veicoli che compongono il sistema. Un modo alternativo per vedere
questo vincolo è il seguente:
Se x∗ ∈ ker(A + BF ) allora il sistema in anello chiuso ha un equilibrio in x∗ .
Se il vincolo sul baricentro è soddisfatto allora il mio sistema in retroazione avrà per
costruzione una coppia di autovalori in 0, che sono sempre presenti anche nel sistema ad
anello aperto. Purtroppo questo vincolo non esclude la più banale delle soluzioni ovvero
F = −A, al fine di evitare questa evenienza indesiderata si impone un vincolo che non annulli
degli elementi di F. Dopo vari test è stato visto che gli elementi di cui è più facile evitare
l’annullamento sono quelli sulla diagonale dei blocchi F1 ed F2 . Imponendo dei vincoli del
tipo Fi,j ≥ Ai,j sugli elementi designati si riusciva ad evitare l’annullamento delle dinamiche
ad anello chiuso. Formalizziamo quest’ultimo vincolo. Le matrici F1 ed F2 appartengono ad
R2×2 per costruzione e si ripetono n volte nella matrice F. Si vuole quindi che:
F1,i,i 6= R(α)i,i i = 1, 2,
F2,i,i 6= −R(α)i,i in modo tale da poter evitare la soluzione banale del problema. Visto che la matrice R(α) è
nota si vanno a sostituire i rispettivi elementi nei vincoli:
F1,i,i 6= cos(α) i = 1, 2.
(3.15)
F2,i,i 6= − cos(α) Osservando i vincoli in 3.15 si notano facilmente due cose:
1. I vincoli sono speculari,
2. i vincoli non sono lineari,
3. i vincoli dipendono da α.
31
3.3. RISULTATI DELLE SIMULAZIONI E SCELTE PROGETTUALI
La specularità dei vincoli è dovuta a 3.14 ed al fatto che ogni agente dipende complessivamente
da soli 2 veicoli, si può quindi utilizzare solo uno dei due vincoli e soddisfare comunque anche
l’altro7 . Per quanto riguarda la non linearità dei vincoli, vista la natura del problema da
risolvere (ottimizzazione lineare), si deve provvedere ad una sostituzione. Il valore che si vuole
evitare divide la retta dei punti possibili in due parti, si deve quindi scegliere se nella soluzione
noi vogliamo:
F2,i,i > −κ · cos(α);
F2,i,i < −κ · cos(α),
(3.16)
(3.17)
dove i = 1, 2 e κ è un parametro che uso per evitare che la soluzione vada a cadere in un
punto che differisce di un infinitesimo dal punto che non desidero8 . Da vari test il vincolo 3.16
risulta essere il migliore, questo viene però applicato in una sua variante shiftata:
F2,i,i > κ · cos(α), i = 1, 2,
(3.18)
Il motivo dietro a questa scelta verrà spiegato nella sezione 3.3.1, dove si illustrano anche
alcune tematiche legate alla scelta del solver per l’ottimizzazione. Per quanto riguarda la
dipendenza dei vincoli da α questo potrebbe essere un problema visto che noi usiamo angoli
α ∈ [−π, π], andando però ad osservare la funzione coseno si vede che nell’intervallo [−π, 0]
ed in quello [0, π] il comportamento è uguale quindi si può lasciare il vincolo invariato senza
agire su eventuali cambi di segno. Ora che tutti i vincoli sono stati introdotti non ci resta che
vedere la formulazione finale del nostro problema di ottimizzazione:
3.3
funz.obiettivo : {















vincoli :














GA + GBF = HG,
Hw ≤ 0,
−ρL ≤ F xi ≤ ρU , i = 1, 2, · · · , r,
hij ≥ 0, i 6= j,
FP
= circ[F1 , F2 , 02 , · · · , 02 ],
n
Fi,2j−1 = 0 Pj=1
n
i = 1, · · · , n,
j=1 Fi,2j = 0
F2,i,i > κ · cos(α).
Risultati delle simulazioni e scelte progettuali
In questa sezione mostrerò alcuni risultati delle simulazioni ed illustrerò gli strumenti necessari
per svolgerle. Il nostro “workbench” è stato Matlab.
3.3.1
Strumenti per le simulazioni
Le nostre simulazioni si articolavano in due passi principali:
• ricerca di una soluzione per il problema di ottimizzazione vincolata,
7
8
Si sceglie di utilizzare il vincolo su F2
Ovviamente κ sarà minore di 1 se uso il vincolo 3.16, mentre sarà maggiore di 1 se uso 3.17, κ > 0.
32
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
• simulazione del sistema retroazionato.
Per svolgere il primo passo è stato usato un solver per programmazione lineare, mentre per la
simulazione è stato usato un modello simulink.
Scelta del solver per l’ottimizzazione
La scelta iniziale per il nostro solver è stata quella di usare il “linprog” e cioè il solver standard
fornito con Matlab. La cosa che si nota immediatamente usando questo solver è che le soluzioni
che trova hanno una scarsa accuratezza (dell’ordine di 10−5 ), oltre a questo si vede che variando
i vincoli le soluzioni che fornisce non sono sempre valide, capita infatti di fornire dei vincoli
e di vederli non rispettati (pur sapendo che esiste una soluzione valida che sfrutta tutti i
vincoli). Il linprog è inoltre molto dipendente dalla funzione obiettivo del problema. Per
ovviare a queste tematiche è stato impiegato il solver “Nag”. Oltre a presentare una maggiore
robustezza computazionale ed una migliore accuratezza nei calcoli (dell’ordine di 10−15 ) in
questo solver i vincoli vengono soddisfatti sempre dove possibile. È stato grazie a queste sue
qualità che è stato possibile raggiungere la formulazione finale del problema. Ovviamente la
maggiore accuratezza richiede un alto numero di iterazioni per poter raggiungere la precisione
richiesta sui risultati finali. Prendendo ad esempio un sistema con 3 veicoli, impiegando per il
solver un numero di iterazioni pari a 3500 si raggiunge una precisione di 10−13 , in questo caso
il solver conclude i suoi calcoli fornendo una soluzione “parziale”, in realtà questa soluzione
rispetta tutti i vincoli ed è quindi utilizzabile per le simulazioni.
Adesso che è stata definita l’accuratezza del solver è facile giustificare il fattore moltiplicativo κ nell’equazione 3.16, se questo mancasse il solver potrebbe tranquillamente maggiorare
nella sua soluzione quella cifra di un infinitesimo. Così facendo il vincolo sarebbe ugualmente
soddisfatto ma nella mia matrice di retroazione avrei elementi che non si annullano per via di
un infinitesimo e quindi la dinamica sarebbe comunque assente.
Passaggio del problema alla forma canonica.
Come noto dalla teoria la forma standard per un problema di ottimizzazione lineare è la
seguente:
funz.obiettivo : {J(x)
Aeq x = beq ,
vincoli :
Ax ≤ b,
Osservando il problema di ottimizzazione su cui noi lavoriamo si vede immediatamente che
deve essere trasformato per raggiungere una forma standard. La prima cosa da fare per
passare in una forma canonica è definire quali sono le variabili che devono essere calcolate
come soluzione del problema e che quindi andranno a costituire il vettore x. Una delle nostre
incognite è indubbiamente F, mentre l’altra è H. Il vettore x sarà quindi costituito da tutti gli
elementi di H ed F ordinati opportunamente. Prendiamo adesso alcuni vincoli del problema
originario ed adattiamoli per passare in forma standard. Iniziamo con il vincolo Hw ≤ 0,
questo si trasforma in:
wm · x ≤ 0,
33
3.3. RISULTATI DELLE SIMULAZIONI E SCELTE PROGETTUALI
dove wm 9 è una matrice fatta come segue

0 ··· 0 w
 0 ··· 0 0

wm =  .
.
 .. · · · .. · · ·
0 ··· 0 ···
0 ···
w ···
.. . .
.
.
0 ···
0
0
..
.
w



.

Prendiamo adesso un altro vincolo ovvero GA + GBF = HG, questo
tm2 )x = mga , dove10 tm1 e tm2 sono matrici fatte come segue:



G1,:
0 ···
0 ···
0 G:,1
0
 0 G1,: · · ·


0 ··· 
0

 0 G:,2



..
.
.
 0
 0
.
0 ··· 
···
.
0



 0


···
0 G1,: · · · 
0

 0 G:,n
 G2,:
 0

0
·
·
·
0
·
·
·
0
G
:,1



 0 G2,: · · ·
 0
0 ··· 
0
G:,2






..
..

.
0 ··· 
···
tm1 = 
0
.
 0
 , tm2 =  0
 0


···
0 G2,: · · · 
0
G:,n

 0
 ..


.
.
.
.
..
..
..
.. 
 .
 0
0
···



 Gq,:
 0

0
·
·
·
0
·
·
·
0
0



 0


G
·
·
·
0
·
·
·
0
0
0
q,:






.
..
 0
 0
0 ··· 
···
0
0
0
···
0 Gq,: · · ·
0
0
0
si trasforma in (tm1 +
···
···
0
0
···
···
···
···
0
0
0
0
···
···
..
.
0
0
···
···
···
···
0
G:,1
G:,2
..
.
G:,n


























Riguardo a tm2 si deve precisare che la prima colonna di zeri sono in realtà dei vettori di zeri,
di lunghezza pari al numero di elementi corrispondenti ad F su x. Per quanto riguarda mga
questo è un vettore colonna che contiene tutti gli elementi della moltiplicazione matriciale
G × A ordinati per righe. Esplicitati questi due vincoli si può vedere chiaramente la struttura
di x, la prima parte è composta dagli elementi di F ordinati per colonna e a seguire ci sono gli
elementi di H ordinati per riga. Un’altro vincolo che richiede un minimo di attenzione è quello
sull’ampiezza del segnale di ingresso ovvero −ρL ≤ F xi ≤ ρU . Al fine di soddisfare questo
vincolo è richiesta la conoscenza di tutti i vertici del poliedro all’interno del quale si vuole
vincolare lo stato. Per ottenere le coordinate dei vertici noi utilizziamo il multi-parametric
toolbox11 . I comandi per creare il poliedro definito dall’equazione Gx ≤ w ed estrarne i vertici
sono i seguenti:
P = polytope(G,w); %crea il politopo a partire dalla matrice G e dal vettore w.
V = extreme(P).
Tralascio la trasformazione in forma standard degli altri vincoli perchè banale a questo
punto. È adesso possibile risolvere il problema di ottimizzazione, l’unica cosa rimasta da fare
è estrarre gli elementi di F dal vettore soluzione e rimetterli in forma matriciale in modo da
poterli usare per le simulazioni. Questo è facilmente fattibile con un ciclo for annidato oppure
usando la seguente istruzione:
9
In questa matrice gli elementi w corrispondono al vettore precedentemente definito.
Con G1,: si intende la prima riga di G che è la matrice dei vincoli sullo stato iniziale.
11
vedi [1]
10
34
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
F = reshape(X(1:m*n),n,m)’;
Nell’istruzione matlab sopra m ed n sono banalmente il numero di righe e colonne su cui
vogliamo ridisporre F.
Simulazione del sistema
Come già precedentemente citato per la simulazione viene usato Simulink, questo tool ci
permette di ricostruire tramite appositi blocchi il nostro sistema. Le dinamiche che utilizziamo
sono a singolo integratore quindi si sfrutteranno gli appositi blocchi “integrator”, in totale
saranno usati n · 2 blocchi dove n è sempre il numero di veicoli.
Ad ogni blocco corrisponde quindi una coordinata di un veicolo, una cosa importante
di cui è necessario assicurarsi prima di cominciare la simulazione è settare come condizione
iniziale per il blocco integratore la variabile legata alla coordinata del blocco. Si procede poi a
collegare opportunamente i blocchi ottenendo così il modello in figura 3.7. Le coordinate dei
Figura 3.7: Modello Simulink del sistema con 4 agenti
veicoli vengono salvate ad ogni passo di campionamento in apposite variabili per poter essere
successivamente sfruttate nelle rappresentazioni grafiche.
3.3. RISULTATI DELLE SIMULAZIONI E SCELTE PROGETTUALI
3.3.2
35
Risultati del problema di ottimizzazione
Andando a risolvere il problema di ottimizzazione sopra definito si vede che a prescindere dal
numero di veicoli impiegati, la soluzione che trovo ha una struttura “canonica”. Mostriamo
con un esempio questa struttura (prendo il caso a 3 veicoli)


−a
sin(α)
a
−sin(α)
0
0
 −sin(α)

−b
sin(α)
b
0
0




0
0
−a
sin(α)
a
−sin(α)
,
F =
(3.19)


0
0
−sin(α)
−b
sin(α)
b



a
−sin(α)
0
0
−a
sin(α) 
sin(α)
b
0
0
−sin(α)
−b
dove i parametri a e b sono i due termini definiti con i vincoli di maggiorazione al fine di non
ottenere un sistema retroazionato nullo. Osservando l’equazione 3.19 si vede chiaramente che
la F è schematizzabile come circ [−F1 , F1 , 0, · · · , 0] , in cui il blocco F1 è visibile nell’equazione
appena citata. Ovviamente questa struttura si ottiene se il problema ha dei valori di r ed M
opportuni e tali per cui esista una matrice H che soddisfa i vincoli. Nel caso questi parametri
non avessero valori opportuni la nostra soluzione diventerebbe banalmente F = −A, che
fornendoci un sistema globalmente nullo mantiene sempre i veicoli connessi12 . Tralasciando
il caso banale che risulta essere di scarso interesse, abbiamo visto come è fatta la matrice di
retroazione nel caso in cui il problema di ottimizzazione viene risolto completamente. Possiamo
quindi osservare che il sistema ad anello chiuso ha una struttura a doppia diagonale come quella
che segue:


−c 0
c
0
0 ··· 0
 0 −d 0
d
0 ··· 0 


 0
0 −c 0
c ··· 0 




..
..
..
..
A + BF =  ...
,
. 0 
.
. 0
.


 0
0 ···
0 −d 0
d 


 c
0
0 ···
0 −c 0 
0
d
0
0 ···
0 −d
dove c = a + cos(α) e d = b + cos(α).
Essendo noto il sistema possiamo iniziare a fare speculazioni sul comportamento o più in
particolare sugli autovalori del sistema. Al fine di poter definire gli autovalori in forma chiusa
si sfruttano le proprietà della matrice ciclica. Una volta definito
ϕi = e
2(i−1)πj
n
sfruttiamo l’equazione 3.2 per trasformare la matrice A + BF in una matrice diagonale a
blocchi. Per poter meglio esplicare questa trasformazione sfrutterò un esempio, scelgo il caso
a 3 veicoli13 per brevità e semplicità.
D1 = −F1 + ϕ1 F1 , ϕ1 = e0 = 1
D2 = −F1 + ϕ2 F1 , ϕ2 = e
D3 = −F1 + ϕ3 F1 , ϕ3 = e
12
13
2πj
n
4πj
n
.
La nostra condizione iniziale è per ipotesi con tutti i veicoli in comunicazione.
Con F1 indico il blocco che si ripete all’interno della matrice A + BF .
36
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
Ora che i blocchi diagonali sono definiti posso riscrivere A+BF come segue:


(ϕ1 − 1)c
0
0
0
0
0


0
(ϕ1 − 1)d
0
0
0
0




0
0
(ϕ2 − 1)c
0
0
0


(A + BF )diag = 

0
0
0
(ϕ
−
1)d
0
0
2




0
0
0
0
(ϕ3 − 1)c
0
0
0
0
0
0
(ϕ3 − 1)d
(3.20)
Una volta terminata la trasformazione, la teoria sulle matrici circolanti mi dice che gli autovalori del sistema sono la collezione degli autovalori di ogni singolo blocco Di . Osservando
la matrice in 3.20 e vedendo la struttura diagonale dei miei blocchi risulta banale definire gli
autovalori del sistema, che sono i seguenti:
λ+
i = (ϕi − 1)c,
λ−
i
= (ϕi − 1)d.
(3.21)
(3.22)
A questo punto non ci resta che richiamare una piccola notazione di analisi complessa per
poter trasformare gli autovalori.
Noto che:
ejx = cos(x) + j sin(x),
posso riscrivere gli autovalori in funzione di seni e coseni ottenendo ad esempio:
2π
2π
+
− 1 + j sin
c.
λ2 = cos
n
n
Sono stati così esplicitati gli autovalori del mio sistema in retroazione. Ovviamente il mio
scopo è quello di avere un sistema stabile14 , andiamo quindi a vedere quali sono le condizioni
tali per cui questo si verifica. Guardando la forma degli autovalori nelle equazioni 3.21 e 3.22
si vede che per ottenere autovalori stabili si devono utilizzare dei parametri c e d positivi.
Questo è dovuto al fatto che la parte reale di ϕ − 1 è sempre negativa. Vediamo di giustificare
questa affermazione, partiamo riprendendo la forma generale del parametro ϕi :
2(i−1)πj
2(i − 1)π
2(i − 1)π
+ j sin
.
= cos
ϕi = e n
n
n
La parte reale di questo parametro è banalmente quella definita dal coseno. Questa funzione
può assumere valori appartenenti al range [−1, 1]. Visto che nei nostri autovalori compare il
termine (ϕi − 1) questo ci dice che la parte reale di questo termine potrà appartenere al range
[-2,0] per via dello shift dato dal termine unitario. Andiamo a vedere quando la parte reale
di questo termine potrebbe assumere valore nullo. Questo avviene ogni qualvolta l’argomento
del coseno assume un valore pari a 2π o un suo multiplo. Nel nostro caso si vorrebbe quindi
che:
2(i − 1)π
= 2π,
n
ovvero che:
i−1
= 1 → (i − 1) = n.
n
14
Il sistema è più precisamente marginalmente stabile in quanto ci sono due autovalori nulli per costruzione.
3.3. RISULTATI DELLE SIMULAZIONI E SCELTE PROGETTUALI
37
Ma il parametro i può assumere solo valori tra 1 ed n, quindi questa evenienza non si verificherà
mai ed abbiamo così dimostrato quanto precedentemente affermato.
Ora che sappiamo che i termini c e d devono essere positivi per la stabilità possiamo capire
per quale motivo il vincolo dell’equazione 3.17 non portava ad avere soluzione. Questo andava
ad imporre dei valori positivi per la parte reale degli autovalori e quindi c’erano dei problemi
con la matrice H che moltiplicata per il vettore w non poteva assumere valori positivi in quanto
vincolata.
3.3.3
Definizione dei parametri liberi c e d
Osservando la matrice ad anello chiuso del sistema abbiamo studiato le caratteristiche degli
autovalori e visto che i parametri c e d devono essere positivi. Non abbiamo però definito dei
valori specifici per questi termini, che possono essere posizionati arbitrariamente nell’intervallo
tra 0 e +∞ (a patto di scegliere degli opportuni valori per r ed M ). Essendo il nostro scopo
originario quello di fare in modo che il sistema in anello chiuso fosse il più verosimigliante
possibile a quello ad anello aperto potremmo sfruttare questi due parametri liberi per fare
in modo che gli autovalori del sistema retroazionato siano il più vicino possibile a quelli
originali. Per poter svolgere questo compito possiamo sfruttare svariati metodi, qui di seguito
ne verranno presentati due:
1. il metodo analitico;
2. il metodo empirico.
In entrambi i metodi si cercherà di minimizzare la distanza tra gli autovalori dei due sistemi.
Al fine di minimizzare questa differenza noi sfruttiamo una norma, tra le opzioni disponibili
potremmo scegliere:
• norma 1, ||(x1 , x2 )||1 = |x1 | + |x2 |
p
• norma 2, ||(x1 , x2 )||2 = |x1 |2 + |x2 |2
• norma ∞, max(|x1 |, |x2 |).
Dovendo noi lavorare su valori complessi la norma 2 ci resta particolarmente comoda in quando
ci restituisce un valore per la distanza che è interamente reale. Ora che abbiamo definito la
norma da utilizzare vediamo i due metodi.
Metodo analitico
Per poter usare questo metodo definisco due vettori z e w, il primo contenente tutti gli autovalori ad anello aperto ed il secondo contenente quelli ad anello chiuso. Quello che vogliamo
fare è effettuare la seguente minimizzazione:
min ||z − w||2 ,
c,d
dove ||z − w||2 =
p
|z1 − w1 |2 + |z2 − w2 |2 + · · · + |zn − wn |2 .
38
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
Sviluppo i calcoli nel caso di 3 robot per poter fornire una soluzione (trascuro i due
autovalori nulli perchè tanto sono sempre presenti in entrambi i vettori). I miei vettori
diventano:




(ϕ2 − 1)c
(ϕ2 − 1)e+jα

 (ϕ2 − 1)e−jα 

 , w =  (ϕ2 − 1)d 
z=
 (ϕ3 − 1)c 
 (ϕ3 − 1)e+jα 
−jα
(ϕ3 − 1)e
(ϕ3 − 1)d
Andandoli a sostituire nella norma ottengo:
v
2 2
u
√ !
√ !
u
3
3
3
3
(cos(α + j sin(α) − c)) + − −
(cos(α − j sin(α) − d)) +
||z − w||2 = t − −
2
2
2
2
2 2
√ !
√ !
3
3
3
3
(cos(α + j sin(α) − c)) + − +
(cos(α − j sin(α) − c)) .
+ − +
2
2
2
2
Svolgendo i conti sotto radice posso semplificare fino ad ottenere:
p
||z − w||2 = 12 + 6c2 + 6d2 − 12c · cos(α) − 12d · cos(α)
A questo punto posso usare il metodo del gradiente per trovare il minimo, in pratica cerco il
valore del parametro che mi annulla la sua derivata parziale (in quel punto il problema ha un
minimo).
12c − 12 cos(α)
∂||z − w||2
= p
,
2
∂c
2 12 + 6c + 6d2 − 12c · cos(α) − 12d · cos(α)
12d − 12 cos(α)
∂||z − w||2
= p
.
2
∂d
2 12 + 6c + 6d2 − 12c · cos(α) − 12d · cos(α)
(3.23)
(3.24)
Osservando le equazioni 3.23 e 3.24, si vede che queste si annullano in un valore comune
che è c = d = cos(α). Con passaggi analoghi è possibile verificare il risultato per qualunque
dimensione dei vettori z e w.
Metodo empirico
Nel metodo sopra descritto si vanno a minimizzare le distanze tra autovalori correlati, vengono
associati insieme i λ+
i corrispondenti all’anello chiuso e all’anello aperto (lo stesso vale per i
λ−
).
Si
crea
così
un
metodo
di assegnamento fisso, ma questa non è la scelta migliore possibile.
i
Infatti osservando la proprietà degli autovalori riportata nell’equazione 3.3, si vede che gli
autovalori sono a coppie complesse coniugate. Visto che questa proprietà si riscontra anche
nel sistema ad anello chiuso conviene assegnare le coppie di autovalori complesse coniugate a
termini nel sistema retroazionato con lo stesso parametro. Una volta associate le coppie di
autovalori si definisce la seguente notazione, ℜ(aa) e ℑ(aa) sono rispettivamente la parte reale
e la parte immaginaria dell’autovalore ad anello aperto, mentre ℜ(ac) e ℑ(ac) sono parte reale
e immaginaria dell’autovalore ad anello chiuso. Si vanno adesso ad analizzare una ad una le
coppie di autovalori e si controlla in quale dei seguenti casi siamo:
• ℜ(aa) > ℜ(ac) e ℑ(aa) > ℑ(ac),
3.3. RISULTATI DELLE SIMULAZIONI E SCELTE PROGETTUALI
39
• ℜ(aa) < ℜ(ac) e ℑ(aa) > ℑ(ac),
• ℜ(aa) > ℜ(ac) e ℑ(aa) < ℑ(ac),
• ℜ(aa) < ℜ(ac) e ℑ(aa) < ℑ(ac).
Da test sperimentali si trova che il parametro con il quale ottengo la
tra una
h minore distanza
i
ℜ(aa) ℑ(aa)
coppia di autovalori appartiene all’intervallo di valori compreso tra ℜ(ac) , ℑ(ac) , ovviamente
con gli estremi che si possono scambiare a seconda del
h caso in cui
i sono (per esempio nel
ℑ(aa) ℜ(aa)
terzo dei casi sopra riportati il mio intervallo sarà tra ℑ(ac) , ℜ(ac) ). Si provvede quindi a
calcolare gli intervalli sopra descritti per ogni coppia di autovalori, una volta fatto questo si
effettua un’unione tra gli intervalli riguardanti uno stesso parametro facendo attenzione al
fatto che l’intervallo di valori ottenuto sia continuo, se ad esempio ci sono delle discontinuità
nell’intervallo che si ottiene queste vengono rimosse ([−1, −0.5] ∪ [−0.4, 0.2] diventa [-1,0.2]).
Ottenuti questi insiemi di valori non resta che trovare il valore che mi minimizza la distanza
tra tutti gli autovalori.
Riassumendo, i passaggi da svolgere sono i seguenti:
1. creare associazioni tra le coppie di autovalori,
2. trovare intervalli con soluzione migliore per ogni coppia di autovalori,
3. fare unione degli intervalli trovati al passo precedente,
4. trovare i valori migliori per c e d all’interno dei due intervalli trovati al passo 3.
Commenti sui 2 metodi: Il primo metodo presenta una soluzione univoca che però risulta
essere banale, il secondo metodo d’altra parte dà soluzioni migliori ma che variano ad ogni
angolo α. In più questo secondo metodo ha un problema notevole che è quello dell’assegnamento, infatti seppur facilmente fattibile da un operatore umano questa tematica risulta essere
complessa da implementare algoritmicamente per un numero di autovalori non piccolo.
Grafici delle simulazioni
Ora che sono stati definiti anche dei valori per c e d possiamo mostrare i risultati delle nostre
modifiche sul controllo in anello chiuso. Le simulazioni riportate in figura 3.8 sono state fatte
con 3 veicoli, un angolo α = π4 , un valore del raggio r = 2 ed M = 10. I valori di c e d sono
nel caso del metodo analitico c = d = 0.7070, mentre nel caso del metodo empirico c = 0.7070
e d = 0.9660. Come si può vedere nonostante il comportamento differisca ancora da quello
in anello aperto questa volta il risultato è decisamente migliore, da notare che tutti i vincoli
sulle distanze sono rispettati.
Per quanto riguarda i valori del metodo empirico si vede che i robot convergono più rapidamente lungo un asse, questo è dovuto al fatto di avere valori per c e d diversi (la convergenza
più rapida è sull’asse cui corrisponde il parametro maggiore). Facendo delle simulazioni al
variare dei parametri c e d è inoltre possibile notare che più questi crescono e maggiore è la
velocità di convergenza, esempi di quanto detto sono visibili in figura 3.9
40
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
4
4
R1
3.8
R1
3.8
R
R
2
3.6
R3
3.6
3.4
3.4
3.2
3.2
y [m]
y [m]
2
R3
3
3
2.8
2.8
2.6
2.6
2.4
2.4
2.2
2.2
2
3
3.5
4
x [m]
4.5
2
5
3
3.5
4
x [m]
(a)
5
(b)
Distanza lungo asse x
Distanza lungo asse x
2
2
1.8
1.8
1.6
1.6
1.4
1.4
|x − x |
1
2
|x2 − x3|
1.2
1
1
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
1
2
3
4
5
t [s]
6
7
8
|x − x |
1
9
10
0
|x1 − x3|
0
1
2
3
4
(c)
5
t [s]
6
7
8
9
10
(d)
Distanza lungo asse y
Distanza lungo asse y
2
2
1.8
1.8
1.6
1.6
1.4
1.4
|y1 − y2|
|y1 − y2|
|y − y |
2
1.2
|y − y |
3
2
1.2
|y − y |
1
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
3
|y − y |
3
1
0
2
|x2 − x3|
1.2
|x1 − x3|
0.8
0
4.5
3
0.2
0
1
2
3
4
5
t [s]
6
7
8
9
10
(e)
0
0
1
2
3
4
5
t [s]
6
7
8
9
10
(f)
4
3.5
R1
R2
R
3
y [m]
3
2.5
2
3
3.5
4
x [m]
4.5
5
5.5
(g)
Figura 3.8: (a) Traiettorie con i parametri del metodo analitico; (b) Traiettorie con i parametri
del metodo empirico; (c),(d) Distanze dei veicoli lungo l’asse x con i parametri rispettivamente
del metodo analitico e di quello empirico; (e),(f) Distanze dei veicoli lungo l’asse y con i
parametri rispettivamente del metodo analitico e di quello empirico; (g) Traiettorie originali
ad anello aperto.
41
4
4
3.5
3.5
3
3
y [m]
y [m]
3.4. CONCLUSIONI
2.5
2
2.5
2
R1
1.5
R1
1.5
R2
R2
R3
1
R4
0
0.5
1
1.5
2
2.5
x [m]
3
3.5
4
4.5
R3
1
5
(a)
R4
0
0.5
1
1.5
2
2.5
x [m]
3
3.5
4
4.5
5
(b)
4
3.5
y [m]
3
2.5
2
R1
1.5
R2
R3
1
R
4
0
0.5
1
1.5
2
2.5
x [m]
3
3.5
4
4.5
5
(c)
Figura 3.9: (a) Comportamento con parametri pari a c = d = 5; (b)Comportamento con
parametri pari a c = d = 40; (c)Comportamento con parametri pari a c = d = 300.
3.4
Conclusioni
Il metodo di controllo retroazionato nasce come un’alternativa alla condizione di connettività
sul raggio del disco di comunicazione. L’obiettivo che ci eravamo prefissati era di trovare una
soluzione in cui il valore del raggio dei dischi di comunicazione potesse essere fissato da un
utente, rimuovendo così la sua dipendenza dagli autovalori e dagli autovettori del sistema.
Quello che abbiamo ottenuto è un sistema funzionante che garantisce la situazione di
comunicazione tra i veicoli mentre questi svolgono le loro funzioni. Andiamo adesso a valutare
pregi e difetti del risultato trovato.
Tra i pregi troviamo sicuramente l’utilizzo di un problema di programmazione lineare
per la risoluzione del problema, questo automatizza il calcolo dei risultati per qualunque
scenario, lasciando all’utente come unica richiesta il settaggio dei parametri dello spazio di
lavoro iniziale15 . Inoltre la libertà di scelta dei due parametri c e d ci permette di variare la
velocità con cui i robot si muovono permettendoci di adattare il task alle nostre necessità.
15
Tra questi parametri troviamo ad esempio l’angolo di vista α, il raggio r ed il parametro fittizio M .
42
CAPITOLO 3. CONTROLLO DECENTRALIZZATO IN RETROAZIONE
Purtroppo l’utilizzo di un problema di ottimizzazione per risolvere la computazione della
matrice di retroazione porta anche delle problematiche. Infatti la soluzione al problema non è
univoca ma solo una delle ∞ possibili, basta osservare il vincolo 3.18 per capirlo. Chi sceglie
la soluzione all’interno dell’insieme di quelle possibili è ovviamente il solver che và quindi
ad assumere una discreta importanza, solver diversi potrebbero portare a soluzioni diverse
in quanto sfruttano algoritmi con filosofie non identiche. Si deve poi anche citare il fatto
che il raggio di comunicazione dei veicoli non influisce direttamente sulla soluzione, ma solo
sulla matrice H, assume quindi una notevole importanza il rapporto tra i parametri r ed M ,
inoltre la H viene modificata a seconda dell’evenienza dal solver che talvolta sceglie soluzioni
limitanti per i nostri scopi. Questa eccessiva importanza del solver nella scelta delle soluzioni
ci ha impedito di replicare il comportamento che il sistema originale ottiene con un angolo di
vista α = nπ .
Pensando a sviluppi futuri per questo metodo possiamo immaginare di rimuovere il vincolo
che limita le dipendenze di ogni veicolo a se stesso ed al suo vicino che lo precede in formazione.
Un possibile scenario potrebbe essere quello di avere dipendenze da tre veicoli, ci sarebbe
quindi “dialogo” tra un veicolo ed i suoi due vicini, quello che lo precede in formazione e quello
che lo segue. Questo andrebbe a cambiare la struttura stessa del sistema retroazionato, non
saremmo più nel caso del cyclic pursuit. Questa modifica sostanziale ci permetterebbe però
una possibile svolta in quanto i parametri da controllare non sarebbero più 8, ma ben 1216 .
Con questi nuovi parametri il mio sistema si complicherebbe ed anche le semplici dinamiche
di anello chiuso di cui disponiamo al momento ne uscirebbero modificate. Non sarebbe quindi
da escludere che il solver avendo a disposizione una maggiore libertà facesse scelte più proficue
permettendoci quindi un ulteriore avvicinamento verso le dinamiche originali di anello aperto.
16
Questo vale nel caso con dipendenza da 3 veicoli, se il numero di dipendenze aumentasse ancora i parametri
liberi sarebbero 4 · num.dipendenze.
Bibliografia
[1] M. Kvasnica, P. Grieder, and M. Baotić. Multi-Parametric Toolbox (MPT), 2004.
[2] Giovanni Marro. Teoria dei Sistemi e del Controllo. Zanichelli, 1993.
[3] F. Morbidi, G. Ripaccioli, and D. Prattichizzo. On Connectivity Maintenance in Linear
Cyclic Pursuit. In Proc. IEEE Int. Conf. Robot. Automat, 2009, to appear.
[4] M. Pavone and E. Frazzoli. Decentralized Policies for Geometric Pattern Formation and
Path Coverage. ASME J. Dyn. Syst. Meas. Contr., 129(5):633–643, 2007.
[5] Marina Vassilaki and George Bitsoris. Constrained regulation of linear continuous-time
dynamical systems. Syst. Control Lett., 13(3):247–252, 1989.
[6] Nathan Sorensen Wei Ren. Distributed coordination architecture for multi-robot formation
control, 2004.
43

Compensazione di una rete planimetrica

Corso di Geometria Prof. F.Podest`a Foglio 7 1] Calcolare le matrici

aa 2004-2005 Geometria 2. Matrici simmetriche e forme quadratiche

MATLAB teoria - Il-Cubo

testi

Calcola la circonferenza dell`anello ritagliando il misuratore e

compito 15-09-14

Prova scritta del 3 febbraio 2014, testo e soluzioni

Scarica il modulo per richiedere un preventivo - Assiprev