Relazione Sistemi Ambientali

UNIVERSITÁ DEGLI STUDI DI CATANIA Facoltà di Ingegneria Dipartimento di Ingegneria Elettrica Elettronica e dei Sistemi Corso di Laurea Specialistica in Ingegneria dell'Automazione e del Controllo dei Sistemi Complessi

_______________________________________________________

Modellistica e controllo dei sistemi ambientali Implementazione di un sistema di warning: l’indice ALICE Docente: Prof. G.Nunnari Studenti: L’Episcopo Gaetano, Nicolosi Giuseppe, Privitera Vanessa.

________________________________________________________

Anno Accademico 2007/2008

1. INTRODUZIONE Nel presente lavoro si è provveduto alla realizzazione di un sistema di warning basato sull’analisi di serie temporali. La serie temporale di cui ci si è occupati riguarda il rilevamento eseguito a mezzo di sensori di posizione (nord-est-quota), posizionati sulle pendici del vulcano Stromboli. Lo scopo di tali sensori è quello di rilevare le eventuali variazioni subite dalla superficie vulcanica dovute a fenomeni legati all’attività vulcanico-sismica. Date le particolari condizioni i dati rilevati non costituiscono una serie temporale regolare, non permettendo quindi di definire un tempo di campionamento costante, queste particolari condizioni impongono di dover definire un intervallo di campionamento che rappresenti il giusto compromesso in funzione della serie temporale. Lo scopo principale del software è quello di evidenziare a partire dall’analisi della serie temporale l’insorgenza di fenomeni particolari evidenziabili da variazioni repentine delle misurazioni, in quanto sono previste nella normale attività vulcanica movimenti della superficie ma se queste variazioni subiscono delle accelerazioni risulta utile generare l’evento di warning. Il sistema proposto provvede all’analisi della serie, evidenziando la presenza di eventuali “picchi” senza però gestire la segnalazione diretta del warning. Il sistema di basa sulla valutazione dell’indice “ALICE” (Absolute Local Index of Change Environment).

2. INDICE “ALICE” L’indice “ALICE” rappresenta un metodo di valutazione utilizzato per l’analisi dei warning. Questo metodo provvede alla definizione di una soglia calcolata dinamicamente, quindi non a priori, ma sulla base delle misurazioni effettuate. Di fatto l’indice “ALICE” evidenzia una variazione particolarmente rilevante del parametro monitorizzato al tempo t , comparandolo con il valore medio delle misurazioni effettuate, ma soprattutto tenendo conto della “storia” delle osservazioni eseguite. L’indice “ALICE” è definibile secondo la relazione:

A∆TX (T ) =

∆T x (t ) − ∆T x (t )

σ ∆Tx (t )

Dove ΔTx(t)= | Tx(t)- Tx(t-tc) | rappresenta la differenza in valore assoluto tra la misura effettuata all’istante corrente t e la misura acquisita all’istante precedente di campionamento (tc ), mentre <ΔTx(t)> e σΔT(t) rappresentano rispettivamente la media e la deviazione standard di ΔTx(t) agli istanti precedenti.

3. ORGANIZZAZIONE DEI DATI I dati della serie temporale sono stati forniti in formato Access, si è resa necessaria quindi una fase di pre-elaborazione dei dati. All’interno del database non era possibile distinguere le misure effettuate dai vari sensori (in totale nove) e i relativi istanti di rilevamento. Si è provveduto alla suddivisione, tramite Query, dei dati disponibili per sensore organizzati in differenti tabelle Access, convertite in formato Excel adattando la precisione delle misurazioni e generando un file di testo (codifica ANSI) per ogni sensore per renderne possibile l’elaborazione a mezzo MatLab. Il file di testo realizzato ha una struttura tabellare; è composto da 8 colonne e da un numero di righe variabile, ma pari al numero di rilevazioni effettuate dal sensore considerato. Le colonne sono organizzate secondo il seguente ordine: I. II.

Valore rilevamento Nord Valore rilevamento Est

III.

Valore rilevamento Quota

IV.

Giorno del rilevamento

V.

Mese del rilevamento

VI.

Anno del rilevamento

VII.

Ora del rilevamento

VIII.

Minuti del rilevamento

Un esempio di file ottenuto è il seguente: 4295126.341634540000000 4295126.317036850000000 4295126.324089920000000 4295126.309820860000000 4295126.298361020000000 4295126.298295440000000

518042.310167631000000 518042.325392020000000 518042.321469928000000 518042.331903958000000 518042.338954070000000 518042.338241171000000

25.854967496315200 25.856474917576500 25.853802051868600 25.851031732036700 25.852013330176200 25.853596825290000

28 6 07 17 49 28 6 07 18 48 28 6 07 19 47 28 6 07 20 26 28 6 07 20 36 28 6 07 20 49

4. DESCRIZIONE E CARATTERISTICHE DELLO SCRIPT MATLAB

Il codice è suddiviso in diverse funzioni; prevede un’interfaccia grafica per l’inserimento dei parametri di esecuzione (file di testo da caricare, nome del sensore, ampiezza della finestra temporale in minuti, ore o giorni, offset iniziale e percentuale di campioni scartate). La finestra d’inserimento è la seguente:

Come già detto, i dati vengono acquisiti da file di testo e organizzati in fasce temporali; a seconda se tale fascia è superiore o inferiore a un giorno, vengono richiamate due diverse funzioni: “alice_giorni” o “alice_min_ore”: nel primo caso si hanno dei gruppi di giorni tra i quali verrà calcolata la differenza (fra due intervalli consecutivi); nel secondo caso l’intera giornata viene

suddivisa in fasce orarie, e la differenza viene calcolata fra due fasce orarie uguali ma per giorni consecutivi. Per ogni fascia temporale viene determinato un valore rappresentativo che sarà nullo se non vi sono misurazioni, sarà pari al valore dato se esso è l’unico per quella fascia e sarà pari alla media fra il primo e l’ultimo valore della misurazione se si hanno più rilevazioni. Ottenuti tali valori è possibile calcolare il valore dell’indice ALICE, come illustrato nella formula precedentemente commentata, facendo attenzione ad aggiornare ad ogni passo le medie e le deviazioni standard dell’insieme delle differenze determinate nelle iterazioni precedenti. Inoltre è prevista, per evitare fluttuazioni dell’indice stesso dovute all’esiguità dei dati, soprattutto nelle iterazioni iniziali, l’impostazione di una soglia relativa al numero di iterazioni, al di sotto della quale l’indice non viene calcolato, ma vengono solo aggiornate la media e la deviazione standard; viceversa viene calcolato anche l’indice tenendo conto anche dei dati al di sotto della soglia. Il valore di tale soglia è impostabile attraverso la voce “percentuali di valori scartati” del prompt grafico iniziale. Calcolato l’indice ALICE (per i dati relativi ai tre riferimenti-nord, est, quota), si procede al plottaggio grafico degli indici ottenuti, in seguito per evidenziare i movimenti anomali del suolo, dovuti ad eventuali attività sismico-vulcaniche, si estrapolano dall’insieme degli indici ALICE ottenuti quelli che costituiscono un’anomalia. L’individuazione della anomalia viene eseguita se si verifica la seguente relazione:

A∆TX (T ) ≥ 〈 A∆TX 〉 + 2σ ∆Tx

Dove A∆TX (T ) rappresenta il valore dell’indice ALICE alla data iterazione, mentre 〈 A∆TX 〉 e σ ∆Tx rappresentano rispettivamente il valore medio e la deviazione standard di tutti gli indici ALICE relativi alla fascia temporale corrispondente. Dai valori così ottenuti si calcola il valore percentuale di anomalie presenti per le singole fasce temporali con relativi grafici ad istogramma.

5. SIMULAZIONI E ANALISI DEI RISULTATI OTTENUTI

Disponendo delle serie temporali relative a numero cinque sensori (da SDF18 a SDF22), si è provveduto ad effettuare l’analisi delle stesse utilizzando il procedimento sopra descritto, automatizzato a mezzo dello script MatLab di cui all’Appendice. Per ogni sensore l’analisi dei dati per la ricerca delle anomalie, cioè dei possibili casi di warning, è stata eseguita utilizzando diverse finestre temporali, in successione: 15 min, 30 min, 60 min, 120 min, 180 min, 240 min, 360 min, 720 min, 1 giorno, 7 giorni. Per ogni sessione sono disponibili i grafici relativi al valore di Est, di Nord e di Quota il tutto per ogni singola fascia oraria, di conseguenza il numero di plot risulta essere considerevole. Lo script realizzato provvede a salvare in automatico i plot generati nominandoli a seconda del sensore considerato, dal parametro monitorizzato e dalla fascia oraria (ad es: SDF18_2100-2115-quota.png, questo file contiene il plot relativo al sensore SDF18, alla fascia oraria dalle 21:00 alle 21:15 e relative al valore di quota). Si ritiene utile presentare la valutazione dei risultati ottenuti relativi solamente ad un sensore (SDF18), considerando tre finestre temporali (15 minuti, 120 minuti e 7 giorni) e tutti e tre i parametri:

-

Finestra temporale di 15 minuti, fascia oraria dalle 12:15 alle 12:30

-

Finestra temporale di 120 minuti, fascia oraria dalle 12:00 alle 14:00

-

Fascia di 7 giorni:

Come si evince dai grafici sopra riportati la presenza di picchi evidenzia l’avvenuta anomalia che rappresenta un probabile caso di warning. Ovviamente man mano che la finestra temporale si allarga la dinamica del segnale si regolarizza e il numero di picchi presenti si riduce. Risulta quindi più conveniente valutare gli eventi con una finestra temporale dell’ordine di quindici minuti poiché con valori maggiori si perde un certo numero di campioni causando una minore accuratezza nella valutazione di eventuali reali anomalie. Si è provveduto inoltre alla realizzazione di alcuni istogrammi volti a rappresentare la percentuale di anomalie rilevate nelle varie finestre temporali, da esse è possibile osservare in quali fasce si sono osservate più anomalie si riportano di seguito gli istogrammi: -

Per intervalli di 15 minuti

-

Per intervalli di 2 ore:

-

Per intervalli di 7 giorni

6. CONCLUSIONI

Prescindendo dalla qualità dei dati in ingresso, di cui non si conosce la relativa attendibilità, è stato comunque possibile osservare una certa efficacia del sistema di warning realizzato. La caratteristica saliente è costituita dalla possibilità di tenere conto della “storia” delle misurazioni durante il calcolo dell’indice ALICE, limitando così l’insorgere di errate segnalazioni, dovute ad esempio alle incertezze delle misurazioni e/o al rumore inevitabilmente presente in tutte le misure. Si è osservata una progressiva riduzione della media delle percentuali di anomalie rivelate all’interno di ogni fascia al crescere dell’intervallo di osservazione, come mostra la seguente tabella, relativa al sensore SDF18: Ampiezza Intervallo 15 minuti 30 minuti 1 ora 2 ore 3 ore 4 ore 6 ore 12 ore 1 giorno 7 giorni

Anomalie Nord 5,48% 11,53% 4,76% 3,1% 3,72% 3,58% 3,48% 2,9% 5,77% 11,11%

Anomalie Est 4,86% 11,53% 4,75% 3,26% 3,2% 3,58% 3,49% 2,9% 5,77% 5,55%

Anomalie Quota 4,63% 11,53% 4,75% 3,51% 2,32% 3,03% 3,23% 2,9% 2,9% 5,55%

In particolare, si osserva che per un intervallo di 15 minuti le percentuale media di anomalie rilevate sia minore rispetto a 30 minuti, ampiezza per la quale si ha un picco nelle anomalie (ciò può essere spiegato col fatto che le misure della serie temporale sono state acquisite in media ogni 20-30 minuti), mentre per ampiezze crescenti si ha un trend decrescente, con qualche piccola fluttuazione. Tale comportamento è spiegabile col fatto che aumentando l’intervallo di osservazione la media fra la prima e l’ultima acquisizione nella data fascia tende a smorzare eventuali picchi nelle misure. Per un intervallo di una settimana, invece, la percentuale media di anomalie rilevate è più alta. Tale trend è confermato se confrontiamo più sensori, come è evidente dai seguenti grafici, in cui si mostra la variazione delle percentuali medie di anomalie rilevate, rispettivamente per il nord, l’est e la quota:

Rimane comunque da realizzare un sistema che a partire dai segnali di warning generati dal presente sistema attui le opportune segnalazioni volte a rendere efficace il sistema.

7. APPENDICE – codice MatLab. -

Interfaccia_alice.m

%% Interfaccia grafica da eseguire per il calcolo dell'idice ALICE. % Tale script richiama le funzioni alice_min_ore.m e alice_giorni.m che a loro volta richiedono gli script calcoloGiorno.m, pulisciZeri2.m e pulisciZeri.m %% AUTORI: L'Episcopo Gaetano, Mangano Claudia, Nicolosi Giuseppe, Privitera Vanessa, Sansone Giuseppe %% Data: Dicembre 2007 clear all close all clc %% CREAZIONE DI UNA FINESTRA DI DIALOGO DI INPUT prompt={'Nome target ','Percorso file da caricare','Valore della finestra temporale','minuti,ore o giorni?','offset (filtraggio del rumore)','percentuale di campioni scartati'}; dlg_title='Settaggio dei Parametri'; num_lines=1; def={'SDF18','C:\Tesina\sdf18.txt','3','giorni','650','30'}; %parametri di default options.Resize='on'; answer=inputdlg(prompt,dlg_title,num_lines,def,options); %% CARICAMENTO DATI DAL FILE DI TESTO DESIDERATO target=answer{1,1}; % Nome del sensore, utile per i plot ed i file da salvare percorso=answer{2,1}; % Percorso da cui leggere il file di testo in ingresso datiCaricati=load(percorso); % Caricamento dati da file di testo dal percorso indicato offset=str2double(answer{5,1}); % lettura offset percentuale=str2double(answer{6,1}); % lettura percentuale dati da scartare %% VERIFICA DATI IMMESSI e DEFINIZIONE DELLE FINESTRE TEMPORALI ampiezza=str2double(answer(3,1)); % Conversione della stringa in numero switch(answer{4,1}) % Scelta dello script da caricare a seconda dell'ampiezza della finestra scelta case {'minuti','minuto'} % Finestra temporale espressa in minuti if((ampiezza<=0)||(ampiezza>1440)) % Verifica sulla correttezza dell'inserimento dei dati della finestra temporale msgbox('Errore nel valore impostato della finestra temporale!','Errore nei dati inseriti!'); else passoMinuti=ampiezza; [dati,numeroCampioniFascia,anomalieNord,anomalieEst,anomalieQuota, percentualeAnomalie,percentualeMediaAnomalie]=alice_min_ore(target,datiCaricat i,passoMinuti,offset,percentuale); % Script richiamato: alice_min_ore.m

end case {'ore','ora'} % Finestra temporale espressa in ore, viene poi convertita in minuti if((ampiezza<=0)||(ampiezza>24)) % Verifica sulla correttezza dell'inserimento dei dati della finestra temporale msgbox('Errore nel valore impostato della finestra temporale!','Errore nei dati inseriti!'); else passoMinuti=ampiezza*60; [dati,numeroCampioniFascia,anomalieNord,anomalieEst,anomalieQuota, percentualeAnomalie,percentualeMediaAnomalie]=alice_min_ore(target,datiCaricat i,passoMinuti,offset,percentuale); % Script richiamato: alice_min_ore.m end case {'giorni','giorno'} % Finestra temporale espressa in giorni if((ampiezza<=0)||(ampiezza>366)) % Verifica sulla correttezza dell'inserimento dei dati della finestra temporale msgbox('Errore nel valore impostato della finestra temporale!','Errore nei dati inseriti!'); else passoGiorni=ampiezza; [dati,percentualeAnomalie]=alice_giorni(target,datiCaricati,passoG iorni,offset,percentuale); % Script richiamato: alice_giorni.m end otherwise disp('formato dell/''ora o dei minuti non corretto'); end

-

alice_min_ore.m

function [dati,numeroCampioniFascia,anomalieNord,anomalieEst,anomalieQuota,percentualeA nomalie,percentualeMediaAnomalie]=alice_min_ore(target,datiCaricati,passoMinut i,offset,percentuale) %% Script che calcola l'indice ALICE, ne rivela le anomalie e salva i plot, per intervalli inferiori o pari ad un giorno. % Richiede gli script calcoloGiorno.m, pulisciZeri2.m e pulisciZeri.m %% AUTORI: L'Episcopo Gaetano, Mangano Claudia, Nicolosi Giuseppe, Privitera Vanessa, Sansone Giuseppe %% Data: Dicembre 2007 %% Dati in ingresso % target=Nome del sensore % datiCaricati=dati letti dal file di testo e passati allo script. Il file dev'essere strutturato in 8 colonne: la 1a: misure nord, la 2a: misure est, la 3a: misure quota, la 4a: giorno dell'acquisizione, la 5a: mese dell'acquisizione, la 6a: anno dell'acquisizione, la 7a: ora dell'acquisizione, la 8a: minuto dell'acquisizione % passoMinuti=ampiezza dellintervallo temporale espresso in minuti % offset=Parametro per eliminare eventuali misure iniziali affette da rumore % percentuale=numero in percentuale dei dati iniziali su cui non calcolare l'indice ALICE %% Dati in uscita % dati=dati caricati in memoria % numeroCampioniFascia=numero di campioni contati per fascia oraria % anomalieNord=anomalie rilevate nei dati del Nord

% anomalieEst=anomalie rilevate nei dati dell'Est % anomalieQuota=anomalie rilevate nei dati della Quota % percentualeAnomalie=percentale di anomalie rilevate per fascia oraria e per Nord, Est e Quota % percentualeMediaAnomalie=media delle percentali di anomalie rilevate nelle varie fascie orarie per Nord, Est e Quota %% Impostazione degli intervalli temporali intervalliGiorno=fix(1440/passoMinuti);% Numero di intervalli in un giorno (1440=numero di minuti in una giornata) resto=rem(1440,passoMinuti); % Calcola eventuali minuti residui giornalieri % Controllo se i minuti residui sono abbastanza da creare un altro intervallo o devono essere inglobati nell'ultimo intervallo della giornata if(resto>(passoMinuti/2)) intervalliGiorno=intervalliGiorno+1; index=1; elseif(resto<=(passoMinuti/2)) index=0; end %% Caricamento ed organizzazione dei dati [rig,col]=size(datiCaricati); % calcolo numero di righe e colonne della matrice da cui leggere i dati for i=offset:1:rig % ciclo sulla riga ora=datiCaricati(i,7); % lettura dati min=datiCaricati(i,8); % lettura dati pointMin=ora*60+min; % calcolo del minuto della giornata in cui avviene l'acquisizione % ricerca dell'intervallo temporale (giornaliero) in cui cade la misura acquisita for h=1:1:intervalliGiorno % Ciclo sul numero degli intervalli del giorno StartMin=(h-1)*passoMinuti; % Estremo inferiore dell'intervallo EndMin=h*passoMinuti; % Estremo superiore dell'intervallo if((h==intervalliGiorno)&&(index==1)) % ampliamento dell'intervallo se i minuti residui sono minori al 50% del passo EndMin=StartMin+resto; elseif((h==intervalliGiorno)&&(index==0)) % creazione di un intervallo se i minuti residui sono maggiori al 50% del passo EndMin=h*passoMinuti+resto; end if((pointMin>=StartMin)&&(pointMin<EndMin)) % Verifica se la misura ricade nell'intervallo in questione codiceIntervallo=h; % Calcolo del codice dell'intervallo temporale codiceGiorno=calcoloGiorno(datiCaricati(i,4),datiCaricati(i,5),dat iCaricati(i,6)); % calcolo del giorno dell'anno % Assemblaggio della matrice dei dati dati(i-offset+1,1)=codiceIntervallo; % Prima Colonna=Codice dell'intervallo temporale dati(i-offset+1,2)=datiCaricati(i,1); % Seconda Colonna=Dati del Nord dati(i-offset+1,3)=datiCaricati(i,2); % Terza Colonna=Dati dell'Est dati(i-offset+1,4)=datiCaricati(i,3); % Quarta Colonna=Dati della Quota dati(i-offset+1,5)=codiceGiorno; % Quinta Colonna=Numero del giorno dell'anno break; end end end

%% Calcolo dell'indice ALICE per Nord, Est e Quota [rig2,col2]=size(dati); % calcolo numero di righe e colonne della matrice dati g=0; % contatore dei campioni per fascia oraria al giorno % ordinamento dei campioni a seconda della fascia di appartenenza for j=1:1:intervalliGiorno % ciclo lungo gli intervalli k=1; % contatore del numero totale di campioni per fascia oraria for i=1:1:rig2 % scorro lungo le righe della matrice dati codInt=dati(i,1); % memorizzo il codice intervallo if(j==codInt) % verifico se il dato letto fa parte del codice % creazione di tabelle per dati del Nord, Est, Quota e giorni in cui k è il numero totale di misure per fascia oraria j datiNord(k,j)=dati(i,2); % Matrice dei dati del Nord datiEst(k,j)=dati(i,3); % Matrice dei dati dell'Est datiQuota(k,j)=dati(i,4); % Matrice dei dati della Quota datiGiorno(k,j)=dati(i,5); % Matrice dei dati del giorno giorno=datiGiorno(k,j); % memorizzo il giorno in cui cade la misura % Conteggio del numero di campioni all'interno di ogni fascia per ogni giorno if((k==1)||((k>1)&&(giorno~=datiGiorno(k-1,j)))) g=1; elseif((k>1)&&(giorno==datiGiorno(k-1,j))) g=g+1; end numeroCampioniFascia(giorno,j)=g; % Assemblaggio della matrice che contiene il numero di campioni per fascia k=k+1; end end end % creazione delle matrici con le medie per ogni intervallo, naturalmente si fa la media tra il primo e l'ultimo campione della fascia... [rigDati,colDati]=size(datiNord); % calcolo il numero max di campioni per fascia oraria che è dato dalla dimensione delle righe delle matrici dati... [ncfrig,ncfcol]=size(numeroCampioniFascia); % calcolo le dimensioni della tabella numeroCampioniFascia for j=1:1:ncfcol % calcolo della media dei campioni di misura per fascia y=1; % indice di lettura dalla matrice dati for ind=1:1:ncfrig x=numeroCampioniFascia(ind,j); % Calcolo del numero di misure per fascia oraria if(x==0) % Se non ci sono misure in una fascia, la media è posta a zero mediaNord(ind,j)=0; mediaEst(ind,j)=0; mediaQuota(ind,j)=0; end if((x==1)&&(y<=rigDati)) % Se c'è una sola misura in una fascia, la media è pari a quel campione mediaNord(ind,j)=datiNord(y,j); mediaEst(ind,j)=datiEst(y,j); mediaQuota(ind,j)=datiQuota(y,j); y=y+x; % aggiornamento dell'indice di lettura dalla matrice dati end if((x>1)&&(y<=rigDati)) % Se ci sono 2 o più misure in una fascia, la media è pari alla media fra il primo e l'ultimo campione della fascia mediaNord(ind,j)=(datiNord(y+x-1,j)-datiNord(y,j))/2; mediaEst(ind,j)=(datiEst(y+x-1,j)-datiEst(y,j))/2; mediaQuota(ind,j)=(datiQuota(y+x-1,j)-datiQuota(y,j))/2;

end

y=y+x; % aggiornamento dell'indice di lettura dalla matrice dati

end

end

%% Costruzione di stringhe ed impostazioni per il plot dei risultati (titolo e nome dei file salvati) fascia=' - Fascia Oraria:'; duepunti=':'; minutiInizio=0; oreInizio=0; minutiFine=0; oreFine=0; spazio=' '; nord='-nord'; est='-est'; quota='-quota'; zero='0'; indexN=1; % Indice per il conteggio delle anomalie per il Nord indexE=1; % Indice per il conteggio delle anomalie per l'Est indexQ=1; % Indice per il conteggio delle anomalie per la Quota for c=1:1:ncfcol % ciclo lungo le colonne per calcolare l'indice ALICE per ogni fascia - c=indice di fascia oraria [N2,G2]=pulisciZeri2(mediaNord(:,c));% pulisco gli zeri e memorizzo i giorni corrispondenti al campione E=pulisciZeri(mediaEst(:,c))'; %pulisco solo gli zeri Q=pulisciZeri(mediaQuota(:,c))'; %pulisco solo gli zeri N=N2'; G=G2'; righe=length(N); % calcolo del numero di campioni solo sui valori Nord perchè è uguale per i 3 tipi di dato xmax=righe; % Calcolo del numero max di righe sogliaPlot=fix(xmax*percentuale/100); % Calcolo della soglia di dati a partire dai quali si calcola l'indice Alice % inizializzo gli indici a zero per aggiornarli ad ogni fascia c ALICENord=0; ALICEEst=0; ALICEQuota=0; for r=1:1:righe % scorro i campioni per ogni fascia oraria if(r==1) % Verifica se si distNord(r)=0; primo elemento distMediaNord(r)=0; distanze per il Nord distSommaNord(r)=0; Nord elemento

distEst(r)=0;

tratta del primo elemento della serie % calcolo della distanza per il Nord del % inizializzazione della media delle % inizializzo la somma delle distanze per il % calcolo della distanza per l'Est del primo

distMediaEst(r)=0; distanze per l'Est distSommaEst(r)=0; distanze per l'Est

% inizializzazione della media delle

distQuota(r)=0; primo elemento distMediaQuota(r)=0; distanze per la Quota distSommaQuota(r)=0; distanze per la Quota else % Nel caso in cui

% calcolo della distanza per la Quota del

% inizializzazione della media delle

% inizializzazione della media delle % inizializzazione della somma delle non si tratta del primo elemento della serie

Nord

distNord(r)=(N(r)-N(r-1)); % aggiornamento della distanza per il

distSommaNord(r)=distSommaNord(r-1)+distNord(r-1); % aggiornamento della somma delle distanze per il Nord distMediaNord(r)=distSommaNord(r)/(r-1); % aggiornamento della media delle distanze per il Nord distEst(r)=(E(r)-E(r-1)); % aggiornamento della distanza per l'Est distSommaEst(r)=distSommaEst(r-1)+distEst(r-1); % aggiornamento della somma delle distanze per l'Est distMediaEst(r)=distSommaEst(r)/(r-1); % aggiornamento della media delle distanze per l'Est distQuota(r)=(Q(r)-Q(r-1)); % aggiornamento della distanza per la Quota

distSommaQuota(r)=distSommaQuota(r-1)+distQuota(r-1); % aggiornamento della somma delle distanze per la Quota distMediaQuota(r)=distSommaQuota(r)/(r-1); % aggiornamento della media delle distanze per la Quota end if(r>=sogliaPlot) % Da questo punto in poi si calcola anche l'indice ALICE sigmaMediaNord(r)=std(distNord(1:r-1),1); % calcolo la deviazione standard dei campioni precedenti del Nord sigmaMediaEst(r)=std(distEst(1:r-1),1); % calcolo la deviazione standard dei campioni precedenti dell'Est sigmaMediaQuota(r)=std(distQuota(1:r-1),1); % calcolo la deviazione standard dei campioni precedenti della Quota ALICENord(r-sogliaPlot+1,1)=(distNord(r)distMediaNord(r))/sigmaMediaNord(r); % Calcolo effettivo dell'indice ALICE per il Nord ALICEEst(r-sogliaPlot+1,1)=(distEst(r)distMediaEst(r))/sigmaMediaEst(r); % Calcolo effettivo dell'indice ALICE per l'Est ALICEQuota(r-sogliaPlot+1,1)=(distQuota(r)distMediaQuota(r))/sigmaMediaQuota(r); % Calcolo effettivo dell'indice ALICE per la Quota end end %% Calcolo della media e della deviazione standard degli indici ALICE mediaALICENord(c)=mean(ALICENord(:,1),1); mediaALICEEst(c)=mean(ALICEEst(:,1),1); mediaALICEQuota(c)=mean(ALICEQuota(:,1),1); stdALICENord(c)=std(ALICENord(:,1),1); stdALICEEst(c)=std(ALICEEst(:,1),1); stdALICEQuota(c)=std(ALICEQuota(:,1),1); %% Verifica delle anomalie anomalieFasciaN=0; % Numero di anomalie per fascia oraria - Nord (inizializzazione) anomalieFasciaE=0; % Numero di anomalie per fascia oraria - Est (inizializzazione) anomalieFasciaQ=0; % Numero di anomalie per fascia oraria - Quota (inizializzazione) for r=1:1:length(ALICENord) % Ciclo sulle colonne dell'indice Alice calcolato if((ALICENord(r))>(mediaALICENord(c)+2*stdALICENord(c)))|| ((ALICENord(r))<(mediaALICENord(c)-2*stdALICENord(c))) % Verifica delle anomalie per il Nord % Assemblaggio della matrice delle anomalie per il Nord anomalieNord(indexN,1)=G(sogliaPlot+r-1); % Prima colonna=Codice del giorno dell'anomalia

anomalieNord(indexN,2)=c; % Seconda Colonna=Codice della fascia oraria dell'anomalia indexN=indexN+1; anomalieFasciaN=anomalieFasciaN+1; end if((ALICEEst(r))>(mediaALICEEst(c)+2*stdALICEEst(c)))|| ((ALICEEst(r))<(mediaALICEEst(c)-2*stdALICEEst(c))) % Verifica delle anomalie per l'Est % Assemblaggio della matrice delle anomalie per l'Est anomalieEst(indexE,1)=G(sogliaPlot+r-1); % Prima colonna=Codice del giorno dell'anomalia anomalieEst(indexE,2)=c; % Seconda Colonna=Codice della fascia oraria dell'anomalia indexE=indexE+1; anomalieFasciaE=anomalieFasciaE+1; end if((ALICEQuota(r))>(mediaALICEQuota(c)+2*stdALICEQuota(c)))|| ((ALICEQuota(r))<(mediaALICEQuota(c)-2*stdALICEQuota(c))) % Verifica delle anomalie per la Quota % Assemblaggio della matrice delle anomalie per la Quota anomalieQuota(indexQ,1)=G(sogliaPlot+r-1); % Prima colonna=Codice del giorno dell'anomalia anomalieQuota(indexQ,2)=c; % Seconda Colonna=Codice della fascia oraria dell'anomalia indexQ=indexQ+1; anomalieFasciaQ=anomalieFasciaQ+1; end end %% Calcolo la percentuale delle anomalie e assemblaggio della matrice (percentualeAnomalie) percentualeAnomalie(c,1)=c; % Prima colonna=codice di fascia percentualeAnomalie(c,2)=(100*anomalieFasciaN)/length(ALICENord); % Seconda colonna=Percentuale Anomalie Nord percentualeAnomalie(c,3)=(100*anomalieFasciaE)/length(ALICEEst); % Terza colonna=Percentuale Anomalie Est percentualeAnomalie(c,4)=(100*anomalieFasciaQ)/length(ALICEQuota); % Quarta colonna=Percentuale Anomalie Quota %% Calcolo dell'intervallo delle fasce orarie e conversione in stringa per la formattazione e memorizzazione dei grafici if (c==1) % Intervallo di partenza minutiInizio=0; oreInizio=0; else % Incremento del passo per l'intervallo (per l'estremo inferiore) if(minutiInizio+passoMinuti<60) minutiInizio=minutiInizio+passoMinuti; else for ho=1:1:23 if((minutiInizio+passoMinuti)>=(ho*60)&&((minutiInizio+passoMi nuti)<((ho+1)*60))) minutiInizio=(minutiInizio+passoMinuti)-(ho*60); oreInizio=oreInizio+ho; break; end end end end if(minutiInizio+passoMinuti<60) % Calcolo dell'estremo superiore dell'intervallo (maniera incrementale) minutiFine=minutiInizio+passoMinuti; else for zh=1:1:23 if(index==0)

if((minutiInizio+passoMinuti)>=(zh*60)&&((minutiInizio+passoMi nuti)<((zh+1)*60))) minutiFine=(minutiInizio+passoMinuti)-(zh*60); oreFine=oreInizio+zh; if(c==ncfcol) oreFine=24; minutiFine=0; end break; end end if(index==1) if(c==ncfcol) oreFine=24; minutiFine=0; break; end if((c=(zh*60))&&((minuti Inizio+passoMinuti)<((zh+1)*60))) minutiFine=(minutiInizio+passoMinuti)-(zh*60); oreFine=oreInizio+zh; break; end end end end % Dettagli stilistici per il titolo e il nome del file if((oreFine==0)&&(minutiFine==0)) oreFine=24; end if(minutiInizio<10) minIniz=num2str(minutiInizio); minStart=strcat(zero,minIniz); else minStart=num2str(minutiInizio); end if(minutiFine<10) minFin=num2str(minutiFine); minEnd=strcat(zero,minFin); else minEnd=num2str(minutiFine); end if(oreInizio<10) oreIniz=num2str(oreInizio); oreStart=strcat(zero,oreIniz); else oreStart=num2str(oreInizio); end if(oreFine<10) oreFin=num2str(oreFine); oreEnd=strcat(zero,oreFin); else oreEnd=num2str(oreFine); end fasciaStart=strcat(oreStart,duepunti,minStart); fasciaEnd=strcat(oreEnd,duepunti,minEnd); titolo=strcat(target,fascia,spazio,fasciaStart,'-',fasciaEnd); %% Impostazione dei limiti sull'asse delle ascisse del plot xmin=G(sogliaPlot);% xmin=limite inferiore del plot (sempre asse x) if(xmin==G(xmax)) % Aggiunto per evitare che il programma si blocchi quando non ci sono misure in una fascia, evita che punto iniziale e punto finale del plot coincidano G(xmax)=xmin+1;

end %% Plot vero e proprio dell'indice Alice figure(1) % Figura 1=Indice ALICE per il Nord hold on; % Abilitazione del mantenimento dei plot sulla stessa figura box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico plot(G(sogliaPlot:xmax),ALICENord(:,1),'r-'); % plot dell'indice ALICE rispetto ai giorni plot([G(sogliaPlot) G(xmax)],[mediaALICENord(c)+stdALICENord(c) mediaALICENord(c)+stdALICENord(c)],'g--'); % plot della deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord(c)-stdALICENord(c) mediaALICENord(c)-stdALICENord(c)],'g--'); % plot della deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord(c)+2*stdALICENord(c) mediaALICENord(c)+2*stdALICENord(c)],'b--'); % plot della doppia deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord(c)-2*stdALICENord(c) mediaALICENord(c)-2*stdALICENord(c)],'b--'); % plot della doppia deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord(c) mediaALICENord(c)],'y-.'); % plot del valore medio dell'indice hold off; % Disabilitazione del mantenimento dei plot sulla stessa figura xlim([xmin G(xmax)]); % Set dell'intervallo sull'asse x xlabel('Giorni'); % Set del nome sull'asse x ylabel('Alice Nord'); % Set del nome sull'asse y titoloNord=strcat(titolo,' - Alice Nord vs Giorni'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(titoloNord); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',oreStart,minStart,'-',oreEnd,minEnd,nord); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente figure(2) % Figura 2=Indice ALICE per l'Est hold on; % Abilitazione del mantenimento dei plot sulla stessa figura box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico plot(G(sogliaPlot:xmax),ALICEEst(:,1),'r-'); % plot dell'indice ALICE rispetto ai giorni plot([G(sogliaPlot) G(xmax)],[mediaALICEEst(c)+stdALICEEst(c) mediaALICEEst(c)+stdALICEEst(c)],'g--'); % plot della deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst(c)-stdALICEEst(c) mediaALICEEst(c)-stdALICEEst(c)],'g--'); % plot della deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst(c)+2*stdALICEEst(c) mediaALICEEst(c)+2*stdALICEEst(c)],'b--'); % plot della doppia deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst(c)-2*stdALICEEst(c) mediaALICEEst(c)-2*stdALICEEst(c)],'b--'); % plot della doppia deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst(c) mediaALICEEst(c)],'y-.'); % plot del valore medio dell'indice hold off; % Disabilitazione del mantenimento dei plot sulla stessa figura xlim([xmin G(xmax)]); % Set dell'intervallo sull'asse x xlabel('Giorni'); % Set del nome sull'asse x ylabel('Alice Est'); % Set del nome sull'asse y titoloEst=strcat(titolo,' - Alice Est vs Giorni'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(titoloEst); % Set del titolo del grafico

set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',oreStart,minStart,'-',oreEnd,minEnd,est); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente figure(3) % Figura 3=Indice ALICE per la Quota hold on; % Abilitazione del mantenimento dei plot sulla stessa figura box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico plot(G(sogliaPlot:xmax),ALICEQuota(:,1),'r-'); % plot dell'indice ALICE rispetto ai giorni plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota(c)+stdALICEQuota(c) mediaALICEQuota(c)+stdALICEQuota(c)],'g--'); % plot della deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota(c)-stdALICEQuota(c) mediaALICEQuota(c)-stdALICEQuota(c)],'g--'); % plot della deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota(c)+2*stdALICEQuota(c) mediaALICEQuota(c)+2*stdALICEQuota(c)],'b--'); % plot della doppia deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota(c)-2*stdALICEQuota(c) mediaALICEQuota(c)-2*stdALICEQuota(c)],'b--'); % plot del valore medio dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota(c) mediaALICEQuota(c)],'y-. '); hold off; % Disabilitazione del mantenimento dei plot sulla stessa figura xlim([xmin G(xmax)]); %Commentare per plottare gli zeri iniziali xlabel('Giorni'); % Set del nome sull'asse x ylabel('Alice Quota'); % Set del nome sull'asse y titoloQuota=strcat(titolo,' - Alice Quota vs Giorni'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(titoloQuota); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',oreStart,minStart,'-',oreEnd,minEnd,quota); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente close all; % Chiusura di tutte le finestre di figure aperte per evitare problemi all'interno del ciclo for end %% Calcolo della media delle percentuali di anomalie per intervalli, sul totale delle rilevazioni, per tutte le fascie orarie percentualeMediaAnomalie(1,1)=mean(percentualeAnomalie(:,2),1); % Prima Colonna=Percentuale di Anomalie per il Nord percentualeMediaAnomalie(1,2)=mean(percentualeAnomalie(:,3),1); % Seconda Colonna=Percentuale di Anomalie per l'Est percentualeMediaAnomalie(1,3)=mean(percentualeAnomalie(:,4),1); % Terza Colonna=Percentuale di Anomalie per la Quota %% Plot delle percentuali di anomalie rilevate nelle varie fasce orarie sul totale di misurazioni (Istogrammi) % Set delle stringhe da utilizzare nelle figure titolo2='Percentuale anomalie su totale misurazioni vs fasce orarie'; % Titolo generico per tutti gli istogrammi subtitoloN=' - Dati Nord'; % Sottotitolo per il Nord subtitoloE=' - Dati Est'; % Sottotitolo per l'Est subtitoloQ=' - Dati Quota'; % Sottotitolo per la Quota

asseX='Fasce orarie [codice progressivo]'; % Nome per l'etichetta dell'asse delle x asseY='Percentuale anomalie su totale misurazioni [%]'; % Nome per l'etichetta dell'asse delle y anomalie='anomalie'; % Stringa utile per il nome del file da salvare figure(1) % Figura 1=Percentuale delle anomalie vs intervalli orari per il Nord box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico bar(percentualeAnomalie(:,1),percentualeAnomalie(:,2)); % Istogramma delle percentuali di anomalie vs intervalli xlim([0 (intervalliGiorno+1)]); % Set dell'intervallo sull'asse x xlabel(asseX); % Set del nome sull'asse x ylabel(asseY); % Set del nome sull'asse y tit=strcat(target,' - ',titolo2,subtitoloN); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(tit); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',anomalie,nord); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente figure(2) % Figura 2=Percentuale delle anomalie vs intervalli orari per l'Est box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico bar(percentualeAnomalie(:,1),percentualeAnomalie(:,3)); xlim([0 (intervalliGiorno+1)]); % Set dell'intervallo sull'asse x xlabel(asseX); % Set del nome sull'asse x ylabel(asseY); % Set del nome sull'asse y tit=strcat(target,' - ',titolo2,subtitoloE); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(tit); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',anomalie,est); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente figure(3) % Figura 3=Percentuale delle anomalie vs intervalli orari per la Quota box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico bar(percentualeAnomalie(:,1),percentualeAnomalie(:,4)); % Istogramma delle percentuali di anomalie vs intervalli xlim([0 (intervalliGiorno+1)]); % Set dell'intervallo sull'asse x xlabel(asseX); % Set del nome sull'asse x ylabel(asseY); % Set del nome sull'asse y tit=strcat(target,' - ',titolo2,subtitoloQ); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(tit); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',anomalie,quota); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente close all % Chiusura di tutte le finestre di figure aperte per evitare problemi

-

alice_giorni.m

function [dati,percentualeAnomalie]=alice_giorni(target,datiCaricati,numGiorni,offset,p ercentuale) %% Script che calcola l'indice ALICE, ne rivela le anomalie e salva i plot, per intervalli superiori o pari ad un giorno. % Richiede gli script calcoloGiorno.m, pulisciZeri2.m e pulisciZeri.m %% Autori: L'Episcopo Gaetano, Mangano Claudia, Nicolosi Giuseppe, Privitera Vanessa, Sansone Giuseppe %% Data: Dicembre 2007 %% Dati in ingresso % target=Nome del sensore % datiCaricati=dati letti dal file di testo e passati allo script. Il file dev'essere strutturato in 8 colonne: la 1a: misure nord, la 2a: misure est, la 3a: misure quota, la 4a: giorno dell'acquisizione, la 5a: mese dell'acquisizione, la 6a: anno dell'acquisizione, la 7a: ora dell'acquisizione, la 8a: minuto dell'acquisizione % passoMinuti=ampiezza dellintervallo temporale espresso in minuti % offset=Parametro per eliminare eventuali misure iniziali affette da rumore % percentuale=numero in percentuale dei dati iniziali su cui non calcolare l'indice ALICE %% Dati in uscita % dati=dati caricati in memoria % percentualeAnomalie=percentale di anomalie rilevate per Nord, Est e Quota %% Impostazione degli intervalli temporali giorniAnno=365; % Numero di giorni dell'anno intervalliAnno=fix(giorniAnno/numGiorni);% Numero di intervalli in un anno resto=rem(giorniAnno,numGiorni); % Calcola eventuali giorni residui annuali % Controllo se i minuti residui sono abbastanza da creare un altro intervallo o devono essere inglobati nell'ultimo intervallo della giornata if(resto>(numGiorni/2)) intervalliAnno=intervalliAnno+1; index=1; elseif(resto<=(numGiorni/2)) index=0; end %% Caricamento ed organizzazione dei dati [rig,col]=size(datiCaricati); % calcolo numero di righe e colonne della matrice da cui leggere i dati for(i=offset:1:rig) % ciclo sulla riga Giorno=calcoloGiorno(datiCaricati(i,4),datiCaricati(i,5),datiCaricati( i,6)); % calcolo del giorno dell'anno % ricerca dell'intervallo temporale (annuale) in cui cade la misura acquisita for(h=1:1:intervalliAnno) % Ciclo sul numero degli intervalli dell'anno StartGior=(h-1)*numGiorni; % Estremo inferiore dell'intervallo EndGior=h*numGiorni; % Estremo superiore dell'intervallo if((h==intervalliAnno)&&(index==1)) % ampliamento dell'intervallo se i giorni residui sono minori al 50% del passo EndGior=StartGior+resto;

elseif((h==intervalliAnno)&&(index==0)) % creazione di un intervallo se i giorni residui sono maggiori al 50% del pass EndGior=h*numGiorni+resto; end if((Giorno>=StartGior)&&(Giorno<EndGior)) % Verifica se la misura ricade nell'intervallo in questione codiceIntervallo=h; % Calcolo del codice dell'intervallo temporale % Assemblaggio della matrice dei dati dati(i-offset+1,1)=codiceIntervallo; % Prima Colonna=Codice dell'intervallo dati(i-offset+1,2)=datiCaricati(i,1); % Seconda Colonna=Dati del Nord dati(i-offset+1,3)=datiCaricati(i,2); % Terza Colonna=Dati dell'Est dati(i-offset+1,4)=datiCaricati(i,3); % Quarta Colonna=Dati della quota dati(i-offset+1,5)=Giorno; % Quinta Colonna=Numero del giorno del'anno break; end end end %% Calcolo dell'indice ALICE per Nord, Est e Quota [rig2,col2]=size(dati); % calcolo numero di righe e colonne della matrice dati for(j=dati(1,1):1:intervalliAnno) % ordinamento dei campioni a seconda della fascia di appartenenza k=1; % contatore del numero totale di campioni per fascia temporale for(i=1:1:rig2) % scorro lungo le righe della matrice dati codInt=dati(i,1); % memorizzo il codice intervallo if(j==codInt) % verifico se il dato letto fa parte del codice % creazione di tabelle per dati del Nord, Est, Quota e giorni in cui k è il numero totale di misure per fascia temporal j datiNord(k,(j-dati(1,1)+1))=dati(i,2); % Matrice dei dati del Nord datiEst(k,(j-dati(1,1)+1))=dati(i,3); % Matrice dei dati dell'Est datiQuota(k,(j-dati(1,1)+1))=dati(i,4); % Matrice dei dati della Quota datiGiorno(k,(j-dati(1,1)+1))=dati(i,5); % Matrice dei dati del giorno gior=datiGiorno(k,(j-dati(1,1)+1)); % memorizzo il giorno in cui cade la misura k=k+1; end end end % creazione delle matrici con le medie per ogni intervallo, naturalmente si fa la media tra il primo e l'ultimo campione della fascia... [rigDati,colDati]=size(datiNord); % calcolo il numero max di campioni per fascia oraria che è dato dalla dimensione delle righe delle matrici dati... for(j=1:1:colDati) % Calcolo del numero di misure per fascia oraria if(datiNord(1,j)~=0) numeroCampioniFascia(j)=length(pulisciZeri(datiNord(:,j))); else numeroCampioniFascia(j)=0; end end y=1; % indice di lettura dalla matrice dati

for(j=1:1:colDati)

% calcolo della media dei campioni di misura per fascia

x=numeroCampioniFascia(j); % lettura del numero di campioni per fascia if(x==0) % Se non ci sono misure in una fascia, la media è posta a zero

mediaNord(j)=0; mediaEst(j)=0; mediaQuota(j)=0;

end if((x==1)) % Se c'è una sola misura in una fascia, la media è pari a quel campione mediaNord(j)=datiNord(1,j); mediaEst(j)=datiEst(1,j); mediaQuota(j)=datiQuota(1,j); end if((x>1)) % Se ci sono 2 o più misure in una fascia, la media è pari alla media fra il primo e l'ultimo campione della fascia mediaNord(j)=(datiNord(x,j)-datiNord(1,j))/2; mediaEst(j)=(datiEst(x,j)-datiEst(1,j))/2; mediaQuota(j)=(datiQuota(x,j)-datiQuota(1,j))/2; end end %% Costruzione di stringhe ed impostazioni per il plot dei risultati (titolo e nome dei file salvati) tit='Ampiezza Fascia= '; giorni=' giorni '; ampiezzaFascia=num2str(numGiorni); titolo=strcat(target,' - ',tit,ampiezzaFascia,' ',giorni); GiorniInizio=0; GiorniFine=0; oreFine=0; spazio=' '; nord='-nord'; est='-est'; quota='-quota'; zero='0'; indexN=1; % Indice per il conteggio delle anomalie per il Nord indexE=1; % Indice per il conteggio delle anomalie per l'Est indexQ=1; % Indice per il conteggio delle anomalie per la Quota [N2,G2]=pulisciZeri2(mediaNord(:));% pulizia degli zeri e memorizzazione dei giorni corrispondenti alla misura acquisita (per evitare la perdita d'informazione nella serie temporale) E=pulisciZeri(mediaEst(:))'; % pulizia degli zeri Q=pulisciZeri(mediaQuota(:))'; % pulizia degli zeri N=N2'; G=G2'; righe=length(N); % calcolo del numero di campioni solo sui valori Nord perchè è uguale per i 3 tipi di dato xmax=righe; % Calcolo del numero max di righe sogliaPlot=fix(xmax*percentuale/100); % Calcolo della soglia di dati a partire dai quali si calcola l'indice Alice % inizializzo gli indici a zero per aggiornarli ad ogni fascia c ALICENord=0; ALICEEst=0; ALICEQuota=0; for(r=1:1:righe) % scorro i campioni per ogni fascia temporale if(r==1) distNord(r)=0; % calcolo della distanza per il Nord del primo elemento

distMediaNord(r)=0;

% inizializzazione della media delle distanze

distSommaNord(r)=0;

% inizializzo la somma delle distanze per il

per il Nord Nord

distEst(r)=0; elemento

% calcolo della distanza per l'Est del primo

distMediaEst(r)=0; % inizializzazione della media delle distanze

per l'Est

distSommaEst(r)=0; % inizializzazione della media delle distanze

per l'Est distQuota(r)=0; % calcolo della distanza per la Quota del primo

elemento

distMediaQuota(r)=0; % inizializzazione della media delle distanze

per la Quota

distSommaQuota(r)=0; % inizializzazione della somma delle distanze per la Quota else % Nel caso in cui non si tratta del primo elemento della serie distNord(r)=(N(r)-N(r-1)); % aggiornamento della distanza per il Nord distSommaNord(r)=distSommaNord(r-1)+distNord(r-1); % aggiornamento della somma delle distanze per il Nord distMediaNord(r)=distSommaNord(r)/(r-1); % aggiornamento della media delle distanze per il Nord distEst(r)=(E(r)-E(r-1)); % aggiornamento della distanza per l'Est distSommaEst(r)=distSommaEst(r-1)+distEst(r-1); % aggiornamento della somma delle distanze per l'Est distMediaEst(r)=distSommaEst(r)/(r-1); % aggiornamento della media delle distanze per l'Est Quota

distQuota(r)=(Q(r)-Q(r-1)); % aggiornamento della distanza per la

distSommaQuota(r)=distSommaQuota(r-1)+distQuota(r-1); % aggiornamento della somma delle distanze per la Quota distMediaQuota(r)=distSommaQuota(r)/(r-1); % aggiornamento della media delle distanze per la Quota end if(r>=sogliaPlot) % Da questo punto in poi si calcola anche l'indice ALICE sigmaMediaNord(r)=std(distNord(1:r-1),1); % calcolo la deviazione standard dei campioni precedenti del Nord sigmaMediaEst(r)=std(distEst(1:r-1),1); % calcolo la deviazione standard dei campioni precedenti dell'Est sigmaMediaQuota(r)=std(distQuota(1:r-1),1); % calcolo la deviazione standard dei campioni precedenti della Quota ALICENord(r-sogliaPlot+1,1)=(distNord(r)distMediaNord(r))/sigmaMediaNord(r); % Calcolo effettivo dell'indice ALICE per il Nord ALICEEst(r-sogliaPlot+1,1)=(distEst(r)distMediaEst(r))/sigmaMediaEst(r); % Calcolo effettivo dell'indice ALICE per l'Est ALICEQuota(r-sogliaPlot+1,1)=(distQuota(r)distMediaQuota(r))/sigmaMediaQuota(r); % Calcolo effettivo dell'indice ALICE per la Quota end end %% Calcolo della media e della deviazione standard degli indici ALICE mediaALICENord=mean(ALICENord(:,1),1); mediaALICEEst=mean(ALICEEst(:,1),1); mediaALICEQuota=mean(ALICEQuota(:,1),1);

stdALICENord=std(ALICENord(:,1),1); stdALICEEst=std(ALICEEst(:,1),1); stdALICEQuota=std(ALICEQuota(:,1),1); %% Verifica delle anomalie anomalieFasciaN=0; % Numero di anomalie per fascia - Nord (inizializzazione) anomalieFasciaE=0; % Numero di anomalie per fascia - Est (inizializzazione) anomalieFasciaQ=0; % Numero di anomalie per fascia - Quota (inizializzazione) for(r=1:1:length(ALICENord)) % Ciclo sulle colonne dell'indice Alice calcolato if((ALICENord(r))>(mediaALICENord+2*stdALICENord))|| ((ALICENord(r))<(mediaALICENord-2*stdALICENord)) % Verifica delle anomalie per il Nord % Assemblaggio della matrice delle anomalie per il Nord anomalieNord(indexN,1)=G(sogliaPlot+r-1); indexN=indexN+1; anomalieFasciaN=anomalieFasciaN+1; end if((ALICEEst(r))>(mediaALICEEst+2*stdALICEEst))|| ((ALICEEst(r))<(mediaALICEEst-2*stdALICEEst)) % Verifica delle anomalie per l'Est % Assemblaggio della matrice delle anomalie per l'Est anomalieEst(indexE,1)=G(sogliaPlot+r-1); indexE=indexE+1; anomalieFasciaE=anomalieFasciaE+1; end if((ALICEQuota(r))>(mediaALICEQuota+2*stdALICEQuota))|| ((ALICEQuota(r))<(mediaALICEQuota-2*stdALICEQuota)) anomalieQuota(indexQ,1)=G(sogliaPlot+r-1); indexQ=indexQ+1; anomalieFasciaQ=anomalieFasciaQ+1; end %% Calcolo la percentuale delle anomalie e assemblaggio della matrice (percentualeAnomalie) percentualeAnomalie(1)=(100*anomalieFasciaN)/length(ALICENord); % Prima colonna=Percentuale Anomalie Nord percentualeAnomalie(2)=(100*anomalieFasciaE)/length(ALICEEst); % Seconda colonna=Percentuale Anomalie Est percentualeAnomalie(3)=(100*anomalieFasciaQ)/length(ALICEQuota); % Terza colonna=Percentuale Anomalie Quota end xGior=G(sogliaPlot);% insieme a xmax rappresenta il limite sull'asse x %% Plot vero e proprio dell'indice Alice figure(1) % Figura 1=Indice ALICE per il Nord hold on; % Abilitazione del mantenimento dei plot sulla stessa figura box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico plot(G(sogliaPlot:xmax),ALICENord(:,1),'r-'); % plot dell'indice ALICE rispetto agli intervalli di giorni plot([G(sogliaPlot) G(xmax)],[mediaALICENord+stdALICENord mediaALICENord+stdALICENord],'g--'); % plot della deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord-stdALICENord mediaALICENordstdALICENord],'g--'); % plot della deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord+2*stdALICENord mediaALICENord+2*stdALICENord],'b--'); % plot della doppia deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICENord-2*stdALICENord mediaALICENord-2*stdALICENord],'b--'); % plot della doppia deviazione standard inferiore dell'indice

plot([G(sogliaPlot) G(xmax)],[mediaALICENord mediaALICENord],'y-.'); % plot del valore medio dell'indice hold off; % Disabilitazione del mantenimento dei plot sulla stessa figura xlim([xGior G(xmax)]); % Set dell'intervallo sull'asse x xlabel('Fascia di Giorni [codice]'); % Set del nome sull'asse x ylabel('Alice Nord'); % Set del nome sull'asse y titoloNord=strcat(titolo,' - Alice Nord vs Fascia di Giorni'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(titoloNord); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',ampiezzaFascia,'g',nord); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente figure(2) % Figura 2=Indice ALICE per l'Est hold on; % Abilitazione del mantenimento dei plot sulla stessa figura box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico plot(G(sogliaPlot:xmax),ALICEEst(:,1),'r-'); % plot dell'indice ALICE rispetto agli intervalli di giorni plot([G(sogliaPlot) G(xmax)],[mediaALICEEst+stdALICEEst mediaALICEEst+stdALICEEst],'g--'); % plot della deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst-stdALICEEst mediaALICEEststdALICEEst],'g--'); % plot della deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst+2*stdALICEEst mediaALICEEst+2*stdALICEEst],'b--'); % plot della doppia deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst-2*stdALICEEst mediaALICEEst2*stdALICEEst],'b--'); % plot della doppia deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEEst mediaALICEEst],'y-.'); % plot del valore medio dell'indice hold off; % Disabilitazione del mantenimento dei plot sulla stessa figura xlim([xGior G(xmax)]); % Set dell'intervallo sull'asse x xlabel('Fascia di Giorni [codice]'); % Set del nome sull'asse x ylabel('Alice Est'); % Set del nome sull'asse y titoloEst=strcat(titolo,' - Alice Est vs Fascia di Giorni'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(titoloEst); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',ampiezzaFascia,'g',est); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente figure(3) % Figura 3=Indice ALICE per la Quota hold on; % Abilitazione del mantenimento dei plot sulla stessa figura box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico plot(G(sogliaPlot:xmax),ALICEQuota(:,1),'r-'); % plot dell'indice ALICE rispetto agli intervalli di giorni plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota+stdALICEQuota mediaALICEQuota+stdALICEQuota],'g--'); % plot della deviazione standard superiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota-stdALICEQuota mediaALICEQuota-stdALICEQuota],'g--'); % plot della deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota+2*stdALICEQuota mediaALICEQuota+2*stdALICEQuota],'b--'); % plot della doppia deviazione standard superiore dell'indice

plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota-2*stdALICEQuota mediaALICEQuota-2*stdALICEQuota],'b--'); % plot della doppia deviazione standard inferiore dell'indice plot([G(sogliaPlot) G(xmax)],[mediaALICEQuota mediaALICEQuota],'y-.'); % plot del valore medio dell'indice hold off; % Disabilitazione del mantenimento dei plot sulla stessa figura xlim([xGior G(xmax)]); % Set dell'intervallo sull'asse x xlabel('Fascia di Giorni [codice]'); % Set del nome sull'asse x ylabel('Alice Quota'); % Set del nome sull'asse y titoloQuota=strcat(titolo,' - Alice Quota vs Fascia di Giorni'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(titoloQuota); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_',ampiezzaFascia,'g',quota); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente close all; %% Plot dell'istogramma della percentuale delle anomalie figure(1) % Figura 1=Percentuale delle anomalie vs intervalli orari per il Nord box on; % Abilitazione del box sul grafico zoom on; % Abilitazione dello zoom sul grafico bar(percentualeAnomalie); % Istogramma delle percentuali di anomalie vs dati xlim([0 4]); % Set dell'intervallo sull'asse x xlabel('Nord - Est - Quota'); % Set del nome sull'asse x ylabel('Percentuale anomalie su totale misurazioni [%]'); % Set del nome sull'asse y tit=strcat(target,' - ','Percentuale anomalie su totale misurazioni vs Dati'); % Uso della concatenazione di stringhe per costruire il titolo del grafico title(tit); % Set del titolo del grafico set(gca,'FontSize',11); % Set della dimensione dei caratteri della figura file=strcat(target,'_','percentuale-anomalie',ampiezzaFascia,'g'); % Uso della concatenazione di stringhe per costruire il nome del file per il grafico da salvare print('-dpng',file); % Salvataggio della figura su file .png col nome costruito nell'istruzione precedente close all