La Mosca dell`olivo - Dipartimento di Ingegneria dell`Informazione

Università degli Studi di Firenze
Facoltà di Ingegneria
C.di L. in Ambiente e Territorio
Anno Accademico 2001-2002
Corso di Teoria dei Sistemi
Elaborato di Modellazione matematica
La Mosca dell’olivo
Prof. Alessandro Casavola
Studenti: Gabriele De Rulli, Francesco Innocenti
INDICE
• Capitolo 1
• Capitolo 2
pag. 1
1.1
1.2
Introduzione 1
Dati Raccolti 7
pag.12
2.1 Costruzione del modello matematico 12
2.1.1 Modello a sei componenti 13
2.1.2 Conclusione modello sei componenti 18
2.2 Costruzione schema Simulink 19
2.3 Verifica modello 20
2.4 Analisi modale 21
2.5 Stabilità, raggiungibilità e osservabilità 25
2.6 Sintesi del sistema 26
2.6.1 Conclusioni sull’allocazione degli
autovalori 29
• Capitolo 3 Programmi MATLAB
pag.30
3.1 Calcolo errore ai minimi quadrati del modello a
2 componenti 30
3.2 Calcolo errore ai minimi quadrati del modello a
4 componenti 31
3.3 Calcolo errore ai minimi quadrati del modello a
6 componenti 32
3.4 Modello a sei componenti 32
3.5 Modello con temperatura costante 34
3.6 Simulazione ottenuta utilizzando i dati reali agli
istanti precedenti 35
3.7 Modello retroazionato con
autovalori 0.7 0.5 0.2 37
3.8 Calcolo autovalori e autovettori 38
3.9 Calcolo della risposta libera del modo
dominante 39
3.10
Calcolo della risposta libera 40
• Grafici Simulink
• Bibliografia
pag.41
pag.43
Capitolo 1
1.1 Introduzione
La mosca dell'olivo, nome scientifico Bactrocera Oleæ (Gmelin), un
tempo classificata come Dacus Oleæ, è un dittero (cioè un insetto con due sole ali)
diffuso principalmente in tutta l'area del Mediterraneo; solo recentemente ne è
stata
rilevata
la
presenza
nell'area
olivicola
della
California.
Il limite settentrionale dell'espansione nell'area mediterranea è il Lago di Garda.
La mosca (Bactrocera oleae) rappresenta il parassita più pericoloso per l’olivo,
anche se l'intensità degli attacchi è assai variabile in relazione alla posizione
geografica,
all'orientamento
del
singolo
appezzamento
e
all'andamento
meteorologico annuale.
Il danno è arrecato dalle larve dell’insetto che scavano gallerie nel frutto per
nutrirsi della polpa. Una parte, anche consistente, delle olive attaccate cade a terra,
mentre su quelle che restano sull’albero si sviluppano marciumi fungini e batterici
che conferiscono all’olio un sapore sgradevole, un'acidità elevata e un alto numero
di perossidi. Tra le varietà particolarmente sensibili ci sono quelle a frutto grosso
come l’Ascolana Tenera, l’Intosso, la Cucco, il Leccino, mentre Dritta e Gentile di
Chieti risultano mediamente meno infestate.
Gli stadi evolutivi della mosca sono:
Svernamento
La mosca dell'olivo sverna prevalentemente allo stadio di pupa nelle drupe
rimaste
appese
sulla
pianta
o
nel
terreno
nelle
vicinanze
dell’olivo.
In zone o in anni con inverni particolarmente miti, con drupe che rimangono sulla
pianta anche in inverno (per esempio incolto produttivo, olivi selvatici), l'insetto può
sopravvivere anche allo stadio larvale o addirittura di adulto.
Sfarfallamento e ovideposizione
Indipendentemente dal comportamento dell'insetto nel periodo invernale primaverile (gennaio-giugno) la femnmina della mosca non ovidepone finché le
olive non abbiano raggiunto la fase fenologica dell'indurimento del nocciolo.
In questa fase iniziano le ovideposizioni, spesso precedute dalle punture di prova,
che
in
generale
possono
risultare
sterili
con
delle
rare
eccezioni.
Il punto perforato appare come una macchia triangolare brunastra lunga 1-1.5 mm,
che presenta una fessura a mezza luna alla base, la ferita dopo una ventina di
giorni suberifica. Una femmina può infestare diverse olive (in condizioni di
laboratorio
che
risultano
ottimali,
arriva
anche
a
qualche
centinaia).
Ciascuna oliva normalmente subisce una sola ovideposizione, anche se, in annate di
scarica delle olive con forte infestazione non è raro trovare olive infestate da due o
più larve.
Sviluppo larvale
Dalle uova dopo qualche giorno di incubazione (soglia termica minima di 910 gradi) fuoriesce la larva.
Le larve presentano 3 stadi di sviluppo, che crescono
alimentandosi della polpa dell'oliva formando una galleria
filiforme e superficiale (larva di I età).
La galleria poi aumenta di spessore avvicinandosi al
nocciolo senza penetrarlo o intaccarlo (larva di II età).
La velocità di sviluppo dipende dalle temperature, dalla
coltivazione e dal grado di maturazione delle drupe.
La larva di III età scava nella drupa una camera a
ridosso del nocciolo realizzando un foro verso l'esterno di
1,5-2 mm che rimane chiuso da l'epidermide.
Impupamento
In estate (luglio-agosto) la larva preferibilmente si inpupa nell'oliva.
A partire dalla seconda generazione (da settembre) la larva matura fuoriesce dalla
camera aprendo il foro, lasciandosi cadere a terra dove diventa pupa, di solito non
molto in profondità (3-5 cm).
Nel periodo che va da luglio a novembre la mosca dell'olivo compie
generalmente tre generazioni. Naturalmente nei vari anni e nelle varie zone ci sono
delle notevoli variazioni, ma possono essere evidenziate delle caratteristiche
generali.
Prima generazione
Definiamo prima generazione quella che generalmente si ha nel periodo di
Luglio
-
Agosto.
Temperature invernali e primaverili basse tendono a posticipare l'inizio della
generazione; è oramai provato che le prime punture coincidono con la fase
fenologia dell'indurimento del nocciolo della drupa che si ha con temperature più
elevate.
Frequentemente questa generazione subisce una forte mortalità a carico delle
forme larvali di I e II età causata dalle alte temperature e favorita da valori bassi di
umidità relativa.
Seconda e terza generazione
Le generazioni successive tendono ad accavallarsi, e ciò le rende a volte
difficilmente distinguibili i vari stadi e conseguentemente anche gli interventi.
In ogni caso esse rappresentano quelle che possono provocare maggiori danni sia
quantitativi (cascola) che qualitativi (ossidazione della polpa in presenza di larve di
III età,che insieme hanno un numero maggiore di pupe e fori di uscita).
In genere da settembre a novembre si hanno quindi due generazioni e le
ovideposizioni si hanno finché le temperature medie rimangono sopra i 14 ºC; in
Liguria normalmente si scende sotto questo valore dopo la prima decade di
novembre.
In condizioni favorevoli possono verificarsi le ovideposizioni della quarta
generazione, generazione che termina il ciclo con gli sfarfallamenti primaverili.
Per la difesa delle colture individuiamo una fase fenologica che va dalla
invaiatura fino alla maturazione. Per ogni provincia viene riportata la soglia di
infestazione attiva oltre la quale è consigliato un trattamento, dipendente dallo
stadio della mosca e da altri fattori, non da ultimo quello economico spesso in
disaccordo con quello ecologico. La soglia di intervento si ha con il 10% di
infestazione attiva, ricordandosi che temperature minime di circa 7°-10° non sono
in grado di bloccare lo sviluppo della mosca, ma lo rallentano in modo consistente.
I trattamenti possono essere di natura chimica o biologica e possono essere fatti
durante i vari stadi, quelli che fino ad oggi hanno dato maggiore esito sono quelli
che si attuano durante lo sviluppo larvale.
Un metodo per il controllo della mosca dell’olivo è quello mediante la cattura
massale degli adulti; il programma prevede l’individuazione di zone olivate
omogenee, il contatto con le aziende per concordare le modalità operative,
sistemazione in loco delle aziende campione di sacchetti esca (fig.sottostante) per il
controllo della popolazione della mosca (giugno-novembre).
Tuttavia il metodo più usato per il controllo della mosca dell’olivo è la stima
dell’infestazione attiva: tale metodo prevede la raccolta e l’analisi di 100 olive per
ogni azienda. Lo studio della singola oliva è volto a determinare se si è avuta
ovideposizione e se c’è stato sviluppo della larva fino alla seconda età; in tutti
questi casi si considera che l’oliva è stata infestata. A seconda delle percentuali di
infestazione attiva registrate si decide o meno di trattare; l’effetto di tale azione è
quello di uccidere le larve presenti nell’oliva. Nel caso in cui le larve siano di terza
età l’efficacia del trattamento si riduce notevolmente perché i danni arrecati al
frutto sono già di notevoli dimensioni.
1.2 Dati raccolti
Abbiamo deciso di studiare le percentuali di infestazione attiva nell’area
agricola della provincia di Firenze. I dati sono stati raccolti dall’ARSIA (Agenzia
Regionale per lo sviluppo e l’Innovazione nel settore Agricolo-forestale).
28
zona1
zona1
zona1
zona1
zona1
zona2
zona2
zona2
zona2
zona2
zona3
zona3
zona3
zona3
zona3
zona4
zona4
zona4
zona4
zona4
zona5
zona5
zona5
zona5
zona5
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
!!
!!
!!
!!
!!
!!
!!
!!
!!
/
!!
!!
!!
/
/
29
0
/
/
/
/
/
/
/
/
/
/
/
/
/
0
/
/
/
/
/
/
/
0 /
0
2
0
0
0
/
0 /
0
0 /
/
0,78
0
/
/
/
2
/
0 /
/
0
0
0
/
31
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
3,85
2,5
0
/
/
4
1
3
4
0
1
0
0
7,55
1
0
0
0
1
0
0
1
1
0,94
1,49
0
0,5
0
0
0 /
0
0 /
0
/
30
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
4
1,75
/
2,88
3,85
9.09
5
0
9
2,86
0,93
6,8
3
0
2
1
1,98
3
2
2
0
tab.1
32
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
3,7
0
0
0
2
1,82
2
15,84
3,88
1
3
2
0,75
12
2
2
4
1
2
3
5
1
1
33
0
0
0
0
0
0
0
0
0
0
0
0
0 !!
0
0
3,77
1
0
2
1
1
0
16.36
9
2
2
1
1
3,64 /
1
5
3
3
3
5
2
2
0
34
0
0
0
1
0
0
0
0
1
0
3
0
/
0
1
5,45
0
2
1
3
1
2 /
22 /
5,77
0 /
3
0
0
2
5
1
1
4
3
1
1
0
35
1
1
1
3
0
0
2
0
1
0
2
0
1
0
10,09
2,59
4
1
2
1,54
26
1,82
0
2
0
3
3
1
2
2
2
6
0
0
zona1
zona1
zona1
zona1
zona1
zona2
zona2
zona2
zona2
zona2
zona3
zona3
zona3
zona3
zona3
zona4
zona4
zona4
zona4
zona4
zona5
zona5
zona5
zona5
zona5
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
zona6
36
3
0
3
2
2
2
0
1
3
1
3
3
1
3
2
15,89
5/
15
3
6
13,21 /
10,83 /
/
/
6
12
13
12,5
4
2/
9,17
4
24 /
6
5
9
13
10
4
5
13
8
/
37
3
3
4
2
3
2
2
0
3
3
4
3
1
2
3
27,5
16,04
10
35
22
4
19
19
14
15
36
23
0
!!
38
1
4
1
1
1 /
2
2
2
0
4
/
3
0 /
!!
!!
/
1
4
10
12
0
2 /
0,94
8,33
1
5
2
5
4
11
6
8
21
7
35
1
39
1
2
3
2
0
1
0
0
1
20
1
0
0
0
2
6
10
2
4
21 /
1
6
1
4
0
0
7
4
1
7
28
2
33 /
5
40
3
0
1 /
0
0
0
2
2
1
1
0
0
3
0
0
4
7
2
0
3
2
1,82
7,55
8,18
10,71 /
1
15 /
7
7
12
11
5
9
23
0
6
41
2
1
1
1
1
0
1
1
2
1
0
4/
0
0
11,82
21
4
2
8
12
5
10
1
14
4
6
9
18
16
10
12
12
2
2
6
42
6
0
0
0
0
3
4
1
0
0
0
0
/
0
0
10 /
18
4
4
2
16
8
7
3 /
17
0
0
3
9
15
18
14
17
11
5
10
2
8
tab.1
Per il nostro studio abbiamo dovuto considerare le temperature medie
settimanali dal 13 Maggio al 17 Novembre del 2002; il periodo in cui la mosca
43
2
0
1
0
0
0
1
0
0
0
0
0
0
0
12
1
3
3
18
16
14
5
4
2
7
13
13
17
15
16
10
7
10
3
9
dell’ulivo ha il suo ciclo vitale. La stazione meteorologica di riferimento è stato
l’Osservatorio Ximeniano di Firenze.
Settimane
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
T.Media
19,21
19,83
19,53
19,26
22,5
28,67
25,78
24,04
25,14
22,08
24,41
24,94
21,73
22,98
24,27
21,8
20,53
20,18
20,58
12,48
15,22
15,4
16,23
16,63
15,16
tab.2
Temperature
Temperatura
30
25
20
15
10
5
0
Temperatura
28
33
38
43
Settimane
graf.1
I dati presi in considerazione per lo sviluppo del modello sono quelli delle
aziende in cui è stato necessario effettuare un trattamento, escludendo quelli
posteriori agli interventi di disinfestazione, dal momento che questi ultimi non
possono essere studiati perché si riferiscono ad una popolazione il cui andamento è
stato modificato dall’uomo.
Sett.
28
29
0
0
0
2
0
0
0
0
0
0,78
IA senzaT 0,222222
IA con T
0,222222
1,118333
1,118333
Az. 1
Az. 2
Az. 3
Az. 4
Az. 5
Az. 6
Az. 7
Az. 8
Az. 9
Az.10
1
0,94
1,49
0,5
2
30
31
32
33
34
35
0
4
4
1
3
4
1
0
7,55
1
0
3,85
2,88
3,85
9.09
5
9
0,93
6,8
2
0
3,7
1,82
2
15,84
3,88
3
0,75
12
1
0
3,77
1
0
16.36
9
2
1
3,64
2
0
5,45
1
2
22
5,77
3
0
0
10,09
1,54
2,555 3,469125
2,55
3,47
1
0
26
1,82
2
0
0
4,399 2,310167 4,468889
4,4
2,31
4,47
5,18125
5,18
tab.3
Sett.
Az. 1
Az. 2
Az. 3
Az. 4
Az. 5
Az. 6
Az. 7
Az. 8
Az. 9
Az.10
IA senzaT
IA con T
36
37
38
39
40
41
42
43
3
15,89
13,21
10,83
3
0
0
2
0,94
8,33
0
2
0
35
20
0
2
4
21
6
4
0
33
0
4
2
1,82
7,55
0
10,71
15
7
0
0
11,82
12
5
10
1
0
0
6
2
0
10
16
8
7
3
0
3
9
2
0
0
9,17
24
6
8
3
27,5
0
0
12,5
4
16,04
0
35
23
11,2625
11,26
20,908
17,29
19
7,32
26,5
10
6,01
5,98
5,8
tab.3
18
16
14
4
7
13
3
7,5
Infestazione Attiva
Andamento IA considerando o non considerando
i dati dopo il trattamento
30
25
20
15
10
5
0
Infestazione Attiva (si
escludono i dati dopo il
trattamento)
28
30
32
34
36
38
40
42
44
Infestazione Attiva (non
si escludono i dati dopo
il trattamento)
Settimane
graf.2
Correlazione IA e Temperature
Temperature-IA
30
25
Temperature
20
15
Infestazione Attiva (si
escludono i dati dopo il
trattamento)
10
5
0
28
30
32
34
36
38
40
42
Settimane
graf.3
Capitolo 2
44
2.1 Costruzione del modello matematico
Abbiamo deciso di trattare il nostro problema utilizzando un modello ARMA
perché si ritiene che la popolazione della mosca dell’ulivo ad un certo istante k
dipenda dalla popolazione stessa e dalle temperature valutate negli istanti
precedenti. La parte iniziale del nostro lavoro è volta a stimare i coefficienti arma
di tale relazione utilizzando il Metodo dei Minimi Quadrati che minimizza la somma
dei quadrati degli scostamenti esistenti tra previsioni e misure.
Per determinare la dimensione del sistema n avendo una serie temporale di N
(N=12) valori successivi di ingresso e uscita abbiamo effettuato quattro tentativi
dovendo verificare che N
3n :
1) n=1 N 3n 12 3 —> relazione verificata
errore stimato con metodo dei minimi quadrati E=132,4045
2) n=2 N 3n 12 6 —> relazione verificata
errore stimato con metodo dei minimi quadrati E=116,3824
3) n=3 N 3n 12 9 —> relazione verificata
errore stimato con metodo dei minimi quadrati E=97,1923
4) n=4 N 3n 12 12 —> relazione verificata con uguaglianza
in questo caso, dato che gli ingressi sono sufficientemente informativi,
la soluzione esiste ed è unica quindi non è possibile stimare
ai minimi quadrati l’errore minimo.
Errore del modello in funzione del numero dei coefficienti
140
120
100
80
60
40
20
0
2 coefficienti
Errore
4 coefficienti
6 coefficienti
numero dei coefficienti
graf.4
E’ stato scelto il modello a sei componenti (n=3) perché è quello che
minimizza l’errore tra quelli che è possibile utilizzare.
2.1.1 Modello a sei componenti
L’equazione del modello ARMA è la seguente :
y(t) = 0.7803 * y(t-1) – 0.0281 * y(t-2) + 0.9799 * y(t-3) -
(1)
- 0.5563 * u(t-1) + 0.0284 * u(t-2) + 0.5465 * u(t-3)
t = Settimane
y = Infestazione Attiva
u = Temperatura
La funzione di trasferimento W(z) è :
- 0.5563 ⋅ z2 + 0.0284 ⋅ z + 0.5465
W(z) = 3
z - 0.7803 ⋅ z 2 + 0.0281 ⋅ z - 0.9799
z = Operatore anticipo unitario
(2)
Nota la funzione di trasferimento è possibile determinare il modello di stato
utilizzando i comandi della Control System Toolbox di Matlab
sys = tf ( right, left )
[A, B, C, D] =ssdata (sys)
dove il vettore right è il vettore dei coefficienti alla destra dell’equazione del
modello ARMA mentre il vettore left è quello alla sinistra.
La rappresentazione di stato quindi ottenuta è :
 0.7803
x(t + 1) =  1.0000
 0

-0.0281
0
1.0000
0.9799 
1 

 0  u (t )
0
x
(
t
)
+

 

 0
0

 
(3)
y(t ) = ( -0.5563 0.0284 0.5465 ) x(t )
Ponendo :
 0.7803

A =  1.0000
0

-0.0281
0
1.0000
0.9799

0


0

1 
b =  0 
0
 
c=( -0.5563 0.0284 0.5465)
Si ottiene la seguente notazione più compatta del sistema ( 3 )
x(t + 1) = A ⋅ x(t ) + b ⋅ u (t )
y (t ) = c ⋅ x (t )
Lo stato iniziale è stato ricavato imponendo che le prime tre risposte del
modello fossero uguali a quelle reali.
 -38.7565 
x(0) =  -37.0932 
 -37.1213 


graf.5
Per valutare se il nostro modello è utilizzabile per temperature diverse da quelle con
cui è stato tarato si sono fatte varie prove utilizzando come ingressi temperature
costanti.
1) Temperatura di ingresso pari a 0°C. Nella realtà la mosca dell’ulivo non
sopravvive a temperature così basse e quindi il nostro modello dovrebbe
fornire un andamento dell’infestazione attiva tendente a 0, invece la
simulazione a T = 0°C presenta un andamento all’infinito (come mostrato
dal grafico 6 nella pagina seguente).
2) Temperatura di ingresso pari a 45°C: anche in questo caso il nostro modello
dovrebbe fornire un andamento dell’infestazione attiva che tende a 0 perché
le uova della mosca muoiono quando la temperatura supera i 40°C, mentre,
come
si vede dal grafico 7 nella
pagina seguente, l’evoluzione
dell’infestazione attiva simulato tende a meno infinito (tale andamento non
ha senso fisico perché l’ infestazione attiva è una variabile positiva).
3) Temperatura di ingresso è 20°: in questo caso la simulazione risulta
efficiente dato che rispecchia abbastanza bene l’andamento reale (vedere
grafico 8 ).
graf.6
graf.7
graf.8
2.1.2 Conclusione modello a sei componenti
Dallo studio dei grafici si osserva che questo modello ARMA con sei parametri
può essere utilizzato per la simulazione solo se le temperature di ingresso sono
sufficientemente vicine a quelle con cui è stato tarato. In caso contrario la
simulazione fornisce risultati che si discostano completamente da quelli reali.
2.2 Costruzione schema simulink
Lo schema simulink
fornisce un’interfaccia grafica che usa vari tipi di
elementi chiamati blocks per creare una simulazione di un sistema dinamico che
può essere modellato con equazioni differenziali o alle differenze la cui variabile
indipendente è il tempo.
schema 1
I ritardi della variabile w(t) definita nel seguente modo
w(t) = u(t) – a1*w(t-1) – a2*w(t-2) – a3*w(t-3)
individuano le componenti del vettore di stato del sistema (3)
x1(t)=w(t-1)
x2(t)=w(t-2)
x3(t)=w(t-3)
E’ tuttavia difficile determinare il significato fisico di tali componenti a
differenza dell’ingresso e dell’uscita del nostro sistema che sono rispettivamente
temperature medie settimanali [°C] ed infestazioni attive [% di olive danneggiate].
2.3 Verifica modello
Un’ulteriore verifica del
modello Arma
è il confronto dei dati reali
dell’infestazione attiva con quelli calcolati utilizzando ad ogni passaggio i dati reali
delle infestazioni attive precedenti.
y(t) = -a1*IA(t-1) –a2*IA(t-2) -a3*IA(t-3) + b1*u(t-1)+ b2*u(t-2) +b3*u(t-3)
t = settimane
y = dati infestazione attiva calcolati
IA = dati infestazione attiva reali
u = temperature
graf.9
Nelle prime settimane la differenza tra dati reali e quelli calcolati è piuttosto
elevata perché all’inizio della simulazione non sono completamente disponibili le
informazioni relative agli istanti precedenti riguardanti le temperature e le
infestazioni attive.
2.4 Analisi Modale
L’analisi modale nel nostro caso si riconduce allo studio delle soluzioni del
sistema
x(t+1)= A x(t)
(4)
dove A rappresenta la matrice del sistema ( 3 )
 0.7803
A =  1.0000
 0

0.9799 

0


0

-0.0281
0
1.0000
Gli autovalori della matrice A sono
ë1= 1.3208
| ë1 | = 1.3208
ë2= - 0.2702 + 0.8179i
| ë2 | = 0.8613
ë3= - 0.2702 - 0.8179i
| ë3 | = 0.8613
arg ë1 = 0
arg ë2 =
1.8899
arg ë3 = - 1.8899
si pone
ñ2 = 0.8613;
ù2 =
1.8899;
Gli autovettori della matrice A sono
-0.7251
v1 = -0.5490

-0.4157




 0.3936 + 0.2919i
v2 =  0.1785 - 0.5402i
-0.6605




 0.3936 - 0.2919i
v3 =  0.1785 + 0.5402i
-0.6605




Determiniamo la matrice cambiamento di coordinate T = [v1 , á 2 , â 2 ] dove
il vettore v1 è l’autovettore associato all’autovalore reale, mentre vettori á2 e â2
rappresentano rispettivamente la parte reale ed immaginaria dell’autovettore v2.
-0.7251 
v1 = -0.5490 
-0.4157 
 0.3936 
α2 =  0.1785 
 -0.6605
 0.2919 
β2 = -0.5402 
 0

ÄR = T-1 A T
0
0
1.3208


∆R = 
0 -0.2702 0.8179 

0 -0.8179 -0.2702 
x(t)= T z(t)
z(t+1)= T-1 x(t+1) = T-1 A x (t) = ÄR z(t)
z(t+1)= ÄR
t
z(0)
z(0)= T-1 x(0)
64.9094 
z (0) = 15.3522 
 7.7699 
x(t)= T ÄR
t
T-1 x(0)
x(t)= ( ë1 )t z1 (0) v1 + ( ñ2 )t [ z2(0) cos(ù2 t) + z3 (0) sen( ù2 t )] á2 +
+ ( ñ2 )t [ - z2(0) sen(ù2 t) + z3 (0) cos(ù2 t)] â2
Come si può notare la soluzione in forma chiusa del sistema (4) è composta
da un modo reale e da due modi complessi che rappresentano la parte oscillatoria.
Dall’analisi dei moduli degli autovalori emerge che l’autovalore dominante è
ë1= 1.3208 .
Di conseguenza per t tendente all’infinito, il sistema evolverà secondo la relazione
x(t)= ( ë1 )t z1 (0) v1
La risposta libera del modo dominante è
y(t)= c x(t)
y(t)= [( ë1 )t z1 (0)] c v1
graf.10
La risposta libera del sistema ( 3 ) è secondo la formula di Lagrange
y(t)=c At x(0)
graf.11
2.5 Stabilità, raggiungibilità ed osservabilità
Il sistema è instabile dato che esiste un autovalore con modulo maggiore
di 1 (| ë1 |= 1.3208 ).
Per valutare la raggiungibilità e l’osservabilità del sistema occorre studiare le
matrici
R = [ b, A b, A2 b ]
matrice di raggiungibilità
O = [ c, c A, c A2 ]T
matrice di osservabilità
Dato che il rango di entrambe le matrici è massimo essendo pari all’ordine del
sistema n = 3 si può concludere che il sistema è sia raggiungibile che osservabile.
Le matrici R ed O si calcolano usando i comandi di Matlab
R = CTRB(A,B)
1.0000 0.7803 0.5808 
R = 
0 1.0000 0.7803

0
0
1.0000
O = OBSV(A,C)
-0.5563 0.0284 0.5465 
O = -0.4057 0.5621 -0.5451 
 0.2456 -0.5337 -0.3975
2.6 Sintesi del sistema
Come si può osservare dal grafico della risposta libera il nostro sistema è
instabile, ciò significa che l’infestazione attiva tenderebbe ad infinito provocando la
distruzione del raccolto se l’uomo non intervenisse trattando con opportuni prodotti
l’oliveto. L’azione di tale trattamento può essere spiegata mediante una retroazione
lineare dello stato che consiste nell’applicare al sistema :
x(t+1) = A x(t) + b u(t)
un ingresso proporzionale allo stato
u(t) = - K x(t)
che rappresenta l’effetto del trattamento.
Il sistema diventa
x(t+1) = (A – b k) x(t)
dato che il sistema è completamente raggiungibile, scegliendo opportunamente il
vettore k, è possibile allocare arbitrariamente gli autovalori di (A – b k).
La risposta invece è
y(t+1) = c x(t+1)
e quindi andando a sostituire in x(t+1) si ottiene
y(t+1) = c (A – b k) x(t)
Lo schema simulink è rappresentato nello schema 2:
schema 2
graf.12
graf.13
graf.14
graf.15
2.6.1 Conclusioni sull’allocazione degli autovalori
Dallo studio dei grafici relativi ai sistemi retroazionati è possibile notare che
diminuendo il modulo dell’autovalore dominante si ottiene una stabilizzazione del
modello più rapida. Dato che il modello retroazionato indica il controllo dell’uomo,
mediante trattamento, della diffusione dell’infestazione attiva, si può ritenere che
gli autovalori della matrice (A+bk) rappresentino l’intensità del trattamento. Nella
realtà si preferisce utilizzare trattamenti che provocano l’abbattimento immediato
dell’infestazione attiva (comportamento simile al grafico di fig.15) perchè dal punto
di vista economico l’intensità del trattamento non incide sui costi se non in quantità
molto ridotta. Per l’aspetto ambientale sarebbero invece auspicabili interventi
multipli con dosaggi minori in modo da stabilizzare il sistema gradualmente
(comportamento simile al grafico di fig.12).
Capitolo 3
Programmi Matlab
3.1 Calcolo errore ai minimi quadrati del modello a 2 componenti
clear
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
Sett. 28
29
30
31
32
33
34
35
36
37
38
39
T=[ 25.1429 22.08571 24.4143 24.9429 21.7286 22.9857 24.2714 21.8000 20.5286 20.1857
20.5833,12.4857 ]
Y=[1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]'
F=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000
25.1429 22.0857 24.4143 24.9429 21.7286 22.9857 24.2714 21.8000 20.5286 20.1857 20.5833 ]'
rangoF = rank(F)
x=0;
x=(inv((F)'*(F)))*(F)'*(Y)
e=Y-(F)*(x)
E=norm(e)^2.
econf=e(3:11)
Econf=norm(econf)^2.
3.2 Calcolo errore ai minimi quadrati del modello a 4 componenti
clear
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
Sett. 28
29
30
31
32
33
34
35
36
37
38
39
T=[ 25.1429 22.08571 24.4143 24.9429 21.7286 22.9857 24.2714 21.8000 20.5286 20.1857
.
20.5833 12.4857]
Y=[ 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000 ]'
F=[ 1.1200
2.5500
3.4700 4.4000
2.3100
4.4700 5.1900 11.2600 20.9000 19.0000
0.2200
1.1200
2.5500 3.4700
4.4000
2.3100 4.4700 5.1900 11.2600 20.9000
22.08571 24.4143
25.1429
24.9429 21.7286 22.9857 24.2714 21.8000 20.5286 20.1857 20.5833
22.08571 24.4143 24.9429 21.7286 22.9857 24.2714 21.8000 20.5286 20.1857]'
rangoF=rank(F)
x=0;
x=(inv((F)'*(F)))*(F)'*(Y)
e=Y-(F)*(x)
E=norm(e)^2.
econf=e(2:10)
Econf=norm(econf)^2.
3.3 Calcolo errore ai minimi quadrati del modello a 6 componenti
clear
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
Sett. 28
29
30
31
32
33
34
35
36
37
38
39
T=[ 25.1429 22.08571 24.4143 24.9429 21.7286 22.9857 24.2714 21.8000 20.5286 20.1857
.
20.5833 12.4857]
Y=[ 3.4700 4.4000
2.3100
4.4700
5.1900 11.2600 20.9000 19.0000 26.5000 ]'
F=[ 2.5500
3.4700
4.4000
2.3100
4.4700
5.1900
11.2600 20.9000
19.0000
1.1200
2.5500
3.4700
4.4000
2.3100
4.4700
5.1900
11.2600
20.9000
0.2200
1.1200
2.5500
3.4700
4.4000
2.3100
4.4700
5.1900
11.2600
24.4143
24.9429
21.7286 22.9857 24.2714
21.8000
20.5286
20.1857
20.5833
22.08571 24.4143
24.9429 21.7286 22.9857
24.2714
21.8000 20.5286
20.1857
22.08571 24.4143 24.9429 21.7286
22.9857
24.2714 21.8000
20.5286]'
25.1429
x=0;
x=(inv((F)'*(F)))*(F)'*(Y)
e=Y-(F)*(x)
Econf=norm(e)^2.
3.4 Modello a sei componenti
clear
% [0.7803 -0.0281
left=[ 1
0.9799 -0.5563
-0.7803 +0.0281
right=[0 -0.5563
0.0284
sys=tf(right,left)
[A,B,C,D]=ssdata(sys)
-0.9799]
0.5465]
0.0284
0.5465] vettore coefficienti modello ARMA
u=[ 25.14285714
% 28° settimana
22.08571429
24.41428571
24.94285714
21.72857143
22.98571429
24.27142857
21.8
20.52857143
20.18571429
20.58333333
12.48571429
15.22
15.4
16.22857143
16.62857143
15.15714286];
%44° settimana
% Inizializzazione
O3(1,:)=C
O3(2,:)=C*A
O3(3,:)=C*(A^2)
O=OBSV(A,C)
P=[ 0
0
C*B
0
C*A*B
C*B]
U2=[25.14285714
22.08571429 ]
Y3=[ 0.2200
1.1200
2.5500]'
x=[inv(O3)]*(Y3-P*U2) % stato iniziale
y=0;
for i=1:17
y(i)=C*x;
x=A*x+B*u(i);
end
y
k=28:44
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000
0 0 0 0 0]
plot(k,y,'b*',k,IA,'go')
title('Confronto Dati Reali-Modello')
xlabel('Settimane')
ylabel('Infestazione Attiva')
legend('Dati forniti dal Modello','Dati Reali')
3.5 Modello con temperatura costante
clear
% [0.7803 -0.0281
left=[ 1
0.9799 -0.5563
-0.7803 +0.0281
right=[0 -0.5563
0.0284
0.0284
0.5465] vettore coefficienti modello ARMA
-0.9799]
0.5465]
sys=tf(right,left)
[A,B,C,D]=ssdata(sys)
u=0; % Temperatura costante pari a 0°C
%Inizializzazione
O3(1,:)=C
O3(2,:)=C*A
O3(3,:)=C*(A^2)
O=OBSV(A,C)
P=[0
C*B
0
0
C*A*B
C*B]
U2=[25.14285714
22.08571429 ]
Y3=[ 0.2200
1.1200
2.5500]'
x=[inv(O3)]*(Y3-P*U2) % stato iniziale del sistema
y=0;
for i=1:12
y(i)=C*x;
x=A*x+B*u;
end
y
k=28:39
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
plot(k,y,'b*',k,IA,'go')
title('Modello con ingresso T=0°')
xlabel('Settimane')
ylabel('Infestazione Attiva')
legend('Dati forniti dal Modello','Dati Reali')
3.6 Simulazione ottenuta utilizzando i dati reali agli istanti
precedenti
clear
% [0.7803 -0.0281
u=[ 25.14285714
0.9799 -0.5563
% 28° settimana
22.08571429
24.41428571
24.94285714
21.72857143
22.98571429
24.27142857
21.8
20.52857143
20.18571429
20.58333333
12.48571429
15.22
15.4
16.22857143
16.62857143
15.15714286];
% 44°
0.0284
0.5465] vettore coefficienti modello ARMA
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
y(1)=0.2200
% y(1) è nullo perché u(0) e y(0) ed i loro valori negativi si suppone che non
esistano!
y(2)=0.7803*IA(1)-0.5563 *u(1)
y(3)=0.7803*IA(2)-0.0281*IA(1)-0.5563 *u(2)+0.0284*u(1)
for i=4:1:12
y(i)=0.7803*IA(i-1)-0.0281*IA(i-2)+0.9799*IA(i-3)-0.5563*u(i-1)+0.0284*u(i-2)+0.5465*u(i-3);
end
y
k=28:39
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
plot(k,y,'b*',k,IA,'go')
title('Confronto Dati reali - Dati ottenuti utilizzando i dati reali agli istanti precedenti ')
xlabel('Settimane')
ylabel('Infestazione Attiva')
legend('Dati Calcolati','Dati Reali')
e1=y-IA
M1=norm(e1)^2.
a(1:9)=y(4:12)
b(1:9)=IA(4:12)
e2=a-b
M2=norm(e2)^2.
3.7 Modello retroazionato con autovalori 0.7 0.5 0.2
clear
% [0.7803 -0.0281
left=[ 1
0.9799 -0.5563
-0.7803 +0.0281
right=[0 -0.5563
0.0284
0.0284
-0.9799]
0.5465]
sys=tf(right,left)
[A,b,c,d]=ssdata(sys)
R=CTRB(A,b)
%Matrice di Raggiungibilità
rkR=rank(R)
O=OBSV(A,c)
%Matrice di Osservabilità
rkO=rank(O)
[k,PREC,MESSAGE]=place(A,b,[0.7,0.5,0.2])
%Inizializzazione
O3(1,:)=c
O3(2,:)=c*A
O3(3,:)=c*(A^2)
O=OBSV(A,c)
P=[0
c*b
0
0
c*A*b c*b]
U2=[25.14285714
22.08571429 ]
Y3=[ 0.2200
1.1200
2.5500]'
x=[inv(O3)]*(Y3-P*U2) % stato iniziale
y=0;
for i=1:30
y(i)=c*x;
ing(i)= -k*x;
x=A*x + b*ing(i);
end
0.5465] vettore coefficienti modello ARMA
y
ing
t=28:57
plot(t,y,'b*')
title('Modello Retroazionato Autovalori 0.7 0.5 0.2')
xlabel('Settimane')
ylabel('Infestazione Attiva')
legend('Dati forniti dal Modello','Dati Reali')
3.8 Calcolo autovalori e autovettori
clear
% [0.7803 -0.0281
left=[ 1
0.9799 -0.5563
-0.7803 +0.0281
right=[0 -0.5563
0.0284
0.0284
-0.9799]
0.5465]
sys=tf(right,left)
[A,B,C,D]=ssdata(sys)
[U,Lambda]=eig(A)
[T,D]=cdf2rdf(U,Lambda)
ad=norm(Lambda,1)
% Autovalore dominante
a1=Lambda(1,1)
a2=Lambda(2,2)
a3=Lambda(3,3)
M1=abs(a1)
% Moduli autovalori
M2=abs(a2)
M3=abs(a3)
f1=angle(a1)
% Fasi autovalori
f2=angle(a2)
f3=angle(a3)
v1=U(:,1)
alfa2=[U(:,2)+U(:,3)]/2
beta2=[U(:,2)-U(:,3)]/[2*i]
x0=[-38.7565 -37.0932 -37.1213]'
z0=inv(T)*x0
T
0.5465] vettore coefficienti modello ARMA
3.9 Calcolo della risposta libera del modo dominante
Clear
% v1 è l'autovettore dominante, 1.3208 è l'autovalore dominante
v1=[-0.7251
-0.5490
-0.4157]
c=[-0.5563
0.0284
0.5465]
a=c*v1
y=0;
for i=1:12
k(i)=i+27;
y(i)= [(64.9094)*(1.3208)^i]*a ;
end
y
plot(k,y,'*b');
title('Risposta Libera del Modo Dominante')
xlabel('Settimane')
ylabel('Infestazione Attiva')
3.10 Calcolo della risposta libera
clear
% [0.7803 -0.0281
left=[ 1
0.9799 -0.5563
-0.7803 +0.0281
right=[0 -0.5563
0.0284
0.0284
0.5465] vettore coefficienti modello ARMA
-0.9799]
0.5465]
sys=tf(right,left)
[A,B,C,D]=ssdata(sys)
%Inizializzazione
O3(1,:)=C
O3(2,:)=C*A
O3(3,:)=C*(A^2)
O=OBSV(A,C)
P=[0
C*B
% Matrice di Osservabilità
0
0
C*A*B
C*B]
U2=[25.14285714
22.08571429 ]
Y3=[ 0.2200
1.1200
2.5500]'
x=[inv(O3)]*(Y3-P*U2)
y=0;
for i=1:12
y(i)=C*x;
x=A*x;
end
y
k=28:39
IA=[0.2200 1.1200 2.5500 3.4700 4.4000 2.3100 4.4700 5.1900 11.2600 20.9000 19.0000 26.5000]
plot(k,y,'b*',k,IA,'go')
title('Confronto Dati Reali-Risposta Libera')
xlabel('Settimane')
ylabel('Infestazione Attiva')
legend('Dati forniti dal Modello','Dati Reali')
BIBLIOGRAFIA
Sergio Rinaldi Carlo Piccardi : I Sistemi Lineari, Città Studi
Edizioni , 1997
A. Cavallo R. Setola F. Vasca: Guida operativa a MATLAB
Simulink e Control toolbox , casa
editrice Liguori , 1996
William J. Palm III : Introduction to MATLAB 6 for Engineers,
Edizioni McGraw-Hill , 2001
A. Casavola : Appunti del corso di teori dei Sistemi, Anno
Accademico 2000-2001 , corso di laurea di ingegneria per
l’ambiente e il territorio Università degli studi di Firenze
SITI INTERNET:
-
www.arsia.it
www.arssa.it
www.sssup.it
www.taggiasca.com