1

Metodi statistici multivariati per il monitoraggio delle prestazioni dei processi produttivi

Prof. Massimiliano BaroloCAPE-Lab – Computer-Aided Process Engineering LaboratoryDipartimento di Ingegneria IndustrialeUniversità di Padova

max.barolo@unipd.ithttp://www.capelab.dii.unipd.it

Due tipiche problematiche nell’industria di processo

Come monitorare le prestazioni di un processo di produzione? capire se qualcosa “sta andando male” rilevare anomalie di funzionamento, anticipando le azioni correttive assegnare una causa a eventuali deviazioni dalle specifiche di produzione (diagnostica)

È possibile stimare in tempo reale la qualità di un prodotto senza poterla misurare?

Mon

itora

ggio

Sens

orist

ica

CWR

VAPORE

CONDENSAREATTORE

SEPARATORE

STRIPPER

CONDENSATORE

COMPRESSORE

CWR

CWS

ANALIZZATORE

ANALIZZATORE

PI

FI

FI

FI

FI FI

FI

FIFI

FITI

TI

TI

JI

TI

TI

LI

LI

LI

PI

PI

FI

XA

XH

XG

XF

XE

XD

XC

XB

XA

XF

XE

XD

XC

XB

ANALIZZATORE

XD

XH

XG

XF

XE

SC

CWS

2

3

612

7

4 11

10

98

5

1

13

A

D

E

C PRODOTTO

SPURGO

2

Variabili di prodotto e di processoImpiego esteso di calcolatori in linea con il processo disponibilità di basi di dati molto estese centinaia/migliaia di variabili di processo raccolte ogni minuto

temperature; portate; livelli; pressioni; … decine di variabili di prodotto/qualità ogni ora/giorno

concentrazioni; viscosità; peso molecolare medio; indice di rifrazione; … trend storici disponibili rapidamente e a costo nullo su parecchi mesi/anni

Le variabili di processo contengono tutte le informazionisull’andamento del processo tutti gli eventi che accadono lasciano una “traccia” sulle variabili di processo prima

ancora che sulla qualità del prodotto una volta individuata dai dati di processo una situazione anomala, è possibile diagnosticarne la

causa i dati di qualità spesso non sono disponibili con frequenza sufficiente

spesso la qualità può essere determinata solo in uno stadio successivo esempio: determinazione del titolo virale in processi di produzione di vaccini

le misure di qualità disponibii possono non essere sufficienti per caratterizzare il prodotto: la qualità è una proprietà multivariabile esempio: se si misurano e si tengono in specifica solo proprietà reologiche di un polimero (p.e.

viscosità), le proprietà del polimero che dipendono dalla struttura chimica (p.e. ramificazione) non possono essere desunte dalle sole misure disponibili

Analizzare i dati di processo (SPC)

Il controllo statistico di processo (SPC) intende monitorare le prestazioni del processo produttivo usando tutte le misuredisponibili (non solo quelle di qualità)

Esistono alcune difficoltà nell’impiego di queste “masse” di dati la dimensione del problema

centinaia/migliaia di dati da osservare il “sovraccarico” di informazione rende nascosta l’informazione stessa

spesso ci si riduce a osservare le serie temporali di un numero ridotto di variabili

le misure disponibili sono spesso molto correlate l’una con l’altra (multicollinearità) il numero di variabili indipendenti di un processo (“gradi di libertà”) è molto inferiore al

numero delle misure disponibili molte delle variabili che si misurano rappresentano, anche se in modo diverso, la risposta del

processo agli stessi (pochi) eventi fondamentali cui è sottoposto (cioè: sono rappresentazioni diverse degli stessi eventi)

c’è difficoltà a interpretare gli effetti simultanei già di sole 2-3 variabili correlate tra loro

le misure sono affette da rumore strumentale difficile distinguere il rumore dalla variabilità dovuta a cause esterne

non tutte sono disponibili alla stessa frequenza

3

Natura “numerica” delle basi di dati industrialiDimensionalmente enormi centinaia di variabili di processo misurate ogni pochi secondi per mesi/anni misure ridondanti

Non di rango pieno le misure sono altamente correlate l’una con l’altra

la correlazione è dovuta ai numerosi vincoli attivi: bilanci materia/energia, specifiche sulla qualità, vincoli operativi, vincoli di sicurezza, regolazioni in retroazione

il rango statistico è basso è indipendente dal numero delle misure fatte dipende dal numero di “sorgenti di variabilità” indipendenti alle quali il processo è

sottoposto

Dati mancanti “buchi” del 10÷20 % di dati sono più norma che eccezione

Basso rapporto segnale-su-rumore è difficile distinguere la variabilità “naturale” da quella indotta da cause esterne

Non causali non sono il risultato di esperimenti progettati ad hoc (experiment design)

in generale, non è possibile individuare l’effetto causale di una variabile su un’altra

Correlazione e causalità

La correlazione è fondamentalmente diversa dalla causalità

4

Correlazione e causalità

La correlazione è fondamentalmente diversa dalla causalità la causalità implica che un cambiamento in una variabile (fattore) determina il

cambiamento in un’altra variabile (risposta) la correlazione indica che il cambiamento in una variabile è associato al cambiamento

in un’altra variabile

Il modo migliore per stabilire relazioni di causalità è condurre una serie di esperimenti pianificati (designed experiments) queste relazioni di solito sono “nascoste” nelle applicazioni multivariate a causa

dell’enorme quantità di dati disponibili

120 140 160 180 200 220 240 26040

45

50

55

60

65

70

75

80

n° a

bita

nti

n° cicogne

Oldenburg (Germania)1930-1936

Che due variabili siano tra loro correlate non significa che l’una sia causa dell’altra!

anno n° cicogne n° abitanti

1930 132 471931 141 471932 169 591933 185 631934 240 651935 248 691936 250 73

Oldenburg (Germania)

Problematiche e tecnicheConsideriamo due tipologie di problema connesse con l’assicurazione di qualità1. analisi di dati di processo per rilevare se il processo è “sotto

controllo” sono disponibili dati di un unico “blocco” ( X ), e si vogliono evidenziare le

informazioni “nascoste” all’interno del blocco di dati analisi delle componenti principali (PCA; principal component analysis)

2. stime di variabili di qualità di prodotto non misurate sono disponibili dati di due “blocchi” ( X di processo e Y di qualità), e si vogliono

stimare in linea le Y disponendo solo delle X senza un modello del processo proiezione su strutture latenti (PLS; projection to latent structures)

Il problema principale deriva dalla ridondanza delle misure disponibili è un approccio opposto rispetto a quello della progettazione degli

esperimenti (experiment design), dove si sceglie a priori di far variare un numero limitato di variabili tra loro indipendenti

Come affronteremo la questione in questa presentazione1. esempio applicativo (motivazione)2. sviluppi teorici3. altre applicazioni

5

0 5 10 15 20 25 301200

1300

1400

1500

1600

LCL

UCL

varia

bile

2

n° lotto0 5 10 15 20 25 30

1200

1300

1400

1500

1600

LCL

UCL

varia

bile

2

n° lotto

Il controllo statistico (convenzionale) di qualitàEsempio: qualità di una fibra sintetica perché la fibra abbia la qualità desiderata, il produttore fa il monitoraggio di (ad

esempio) due variabili di qualità valore medio e limiti di controllo (superiore UCL, e inferiore LCL) derivano da dati pregressi

di produzioni “sotto controllo” un cliente si è lamentato delle prestazioni della fibra fornita col lotto n°8

il n°8 è l’unico lotto (su 30) col quale ci sono stati problemi

0 5 10 15 20 25 300

5

10

15

20

25

30

35

40

LCLLCL

varia

bile

1

n° di lotto

UCL

0 5 10 15 20 25 300

5

10

15

20

25

30

35

40

LCLLCL

varia

bile

1

n° di lotto

UCL

Carte di controllo di Shewart

Un limite delle carte di controllo tradizionaliLe due variabili sono molto correlate l’analisi non può essere di tipo

univariato (in una sola direzione) esiste un’area all’interno della quale

c’è una data probabilità (ad ex., 99%) che il processo sia in-control

al di fuori dell’intervallo di fiducia si hanno situazioni in cui il processo è statisticamente “fuori controllo” (out-of-control)

Il campione n°8 sta ben al di là dell’intervallo di fiducia è una situazione “fuori controllo”,

anche se le misure di qualitàsingolarmente sarebbero entro i limiti

intervallo di fiducia del 99%

il problema con quel lotto avrebbe potuto essere identificato anche prima di venderlo al cliente!

0 5 10 15 20 25 30 35 401200

1250

1300

1350

1400

1450

1500

1550

1600

vari

abile

2

variabile 1

6

Carte convenzionali di controllo – Caso multivariato

Per saggiare la qualità si effettua un test di ipotesi Si hanno q variabili di qualità aventi media (non nota) (q 1) e matrice di

covarianza (q q ) Dalle osservazioni delle variabili, valutare se la media della popolazione assume il

valore ipotizzato ("specifica") 0 H0 (ipotesi nulla): = 0 H1 (ipotesi alternativa): ≠ 0

Test statistico si suppone che le q misure di qualità q abbiano una distribuzione q-normale si preleva un campione qk della popolazione e se ne costruisce una statistica

(della quale si conosce la distribuzione di probabilità) la statistica considerata è lo stimatore 2 :

questa statistica ha una distribuzione chi-quadro con q gradi di libertà si può calcolare a priori un limite di controllo superiore per la statistica (UCL; upper

control limit, che dipende dal limite di fiducia desiderato) risulta invece: LCL=0

se, per il nuovo campione, la statistica supera UCL, esiste una certa probabilità (ex.: 99%) che la media della popolazione non sia 0

la matrice di covarianza (non nota) è stimata con la matrice di covarianza campionaria S di un set di riferimento di dati “sotto controllo” dalla statistica 2

alla statistica T 2 di Hotelling (che ha una distribuzione F )

01

02 μqΣμq

kkk

Molto spesso non si hanno (sufficienti) misure di qualità disponibiliVolendo usare le variabili di processo invece di quelle di qualità il numero di variabili da monitorare aumenta (anche enormemente) le variabili di processo sono molto correlate tra loro non è detto che la loro distribuzione sia multinormale

La matrice di covarianza campionaria di processo S diventa molto mal condizionata è quasi singolare la sua inversione è problematica l’approccio T 2 diventa inapplicabile

Carte convenzionali di controllo – Caso multivariato cont.

01

02 μqSμq

kkkT statistica di Hotelling

Le carte di controllo rimangono a 2 dimensioni• tutti i punti su UCL stanno sul perimetro dell’ellisse• generalizzazione al caso M ≥ 2 dell’ellisse di controllo

Dove sta il problema?0 4 8 12 16 20 24 28

0

5

10

15

20

25

30

campione n°8

Hot

ellin

g T

2

numero di osservazione (k)

UCL (interv. confidenza 99%)

7

Il controllo statistico di processo (SPC)

Sono necessari metodi che “estraggano” dai dati l’informazione utile per la descrizione del processo

Idea di base “comprimere” i dati in uno spazio di dimensione ridotta, all’interno del quale si

possano rappresentare tutti i fenomeni significativi, non appena essi avvengono

le coordinate di questo nuovo spazio possono essere interpretate come delle pseudo-misure (variabili latenti; componenti principali) “apparentate” con tutte le grandezze misurate, ma in numero molto inferiore ad esse in grado di ricostruire la struttura della correlazione esistente tra le variabili di processo

Proiezione su sottospazi latenti

È possibile modellare la variabilità naturalepresente nelle misure industriali di processo?

Esempio: come rilevare anomalie di marciaI dati disponibili 4 diverse variabili di processo, misurate in 20 situazioni diverse di marcia

“normale” successivi aggiornamenti in linea delle stesse misure non disponibile la misura diretta della “qualità” del prodotto

Il problema dall’analisi dei set di misure di processo future, rilevare in linea se il processo è

“sotto controllo” (variabilità naturale, dovuta a cause comuni) oppure “fuori controllo” (variabilità non naturale, indotta da cause esterne) attribuire cause alle eventuali anomalie di marcia, e rimuovere le cause anticipare le azioni correttive

La difficoltà le misure disponibili sono correlate tra loro e rumorose quali tra le variazioni nelle misure di processo sono realmente rappresentative

di variazioni della qualità di produzione?

8

Esempio: come rilevare anomalie di marcia /cont.

Sono davvero necessarie tutte e quattro le misure per capire in quale stato si trova il processo?Quanto correlate tra loro (ridondanti) sono le misure?La variabilità delle misure è “naturale” (common cause variability) oppure dovuta a qualche evento “nuovo” (assignable cause variability)?Da questo tipo di misure, si può capire in linea anche quando il processo è fuori controllo? la variabile primaria di qualità non è misurata in linea con

frequenza sufficiente

In definitiva: è possibile modellare la variabilità “naturale” di questi dati e distinguerla da quella indotta da eventi “esterni”?

x1 x2 x3 x41 10 20.7 13.6 15.5

2 10.5 19.9 18.1 14.8

3 9.7 20 16.1 16.5

4 9.8 20.2 19.1 17.1

5 11.7 21.5 19.8 18.3

6 11 20.9 10.3 13.8

7 8.7 18.8 16.9 16.8

8 9.5 19.3 15.3 12.2

9 10.1 19.4 16.2 15.8

10 9.5 19.6 13.6 14.5

11 10.5 20.3 17 16.5

12 9.2 19 11.5 16.3

13 11.3 21.6 14 18.7

14 10 19.8 14 15.9

15 8.5 19.2 17.4 15.8

16 9.7 20.1 10 16.6

17 8.3 18.4 12.5 14.2

18 11.9 21.8 14.1 16.2

19 10.3 20.5 15.6 15.1

20 8.9 19 8.5 14.7

misure di processo

osse

rvaz

ioni

Esempio: come rilevare anomalie di marcia /cont.

La “variabilità” del processo descritta dalle 4 misure può essere rappresentata mediante 2 sole nuove coordinate (t1 e t2) le nuove coordinate si ottengono “elaborando” le misure disponibili ogni set di misure viene rappresentato con 2 coordinate invece che 4 è l’estensione al caso multivariabile delle tradizionali carte di Shewhart per variabili di

qualità

-4 -3 -2 -1 0 1 2 3 4 5 6 7-4.0-3.5-3.0-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.03.54.0

t 2

t1

Si assume che i dati rappresentino assetti di processo “sotto controllo” esiste un’area all’interno della quale c’è

una data probabilità (ad ex., 95%) che il processo sia in-control

al di fuori dell’intervallo di fiducia si hanno situazioni in cui il processo è “fuori controllo” (out-of-control )

esempio: è cambiata un’alimentazione; si è disattivato un catalizzatore; uno scambiatore è sporco; …

9

-4 -3 -2 -1 0 1 2 3 4 5 6 7-4.0-3.5-3.0-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.03.54.0

t 2

t1

Esempio: come rilevare anomalie di marcia /cont.

Quando arriva un nuovo set di misure di processo… le misure vengono trasformate nelle nuove

coordinate t1 e t2 si verifica se le nuove misure cadono dentro

oppure fuori dell’intervallo di fiducia

x1 x2 x3 x421 9.9 20 15.4 15.9

22 8.7 19 9.9 16.8

23 11.5 21.8 19.3 12.1

24 15.9 24.6 14.7 15.3

25 12.6 23.9 17.1 14.2

26 14.9 25 16.3 16.6

27 9.9 23.7 11.9 18.1

28 12.8 26.3 13.5 13.7

29 13.1 26.1 10.9 16.8

30 9.8 25.8 14.8 15

A partire dal set n°24 c’è stato qualche cambiamento nel processo che ne ha alterato le prestazioni il prodotto sta andando fuori

controllo non è necessario misurare

in linea la qualità!

Il concetto di variabile latente

t1 da (x1, x3, x4)

t2 da (x2, x5)

x1

x2

x3

x4

x5

Variabili di processo(misure reali)

Variabili latenti(misure fittizie)

10

Analisi di dati: PCA (Principal Component Analysis )

PCA è una tecnica per rappresentare la varianza di una matrice X di dati attraverso un numero ridotto di nuove variabili (“latenti”)dette componenti principali (PC)Esempio 2 misure di natura diversa (x1 ; x2) campionate a 9 istanti diversi matrice X di dati [9×2] ogni campione di 2 misure si può rappresentare come un punto in uno spazio di

dimensione 2 (piano), nelle coordinate x1 e x2

x1

x2

12

3

4

5

6

78

9

La prima PC è il vettore che individua la direzione di massima variabilità dei datiOgni set di dati può essere rappresentato mediante una solacoordinata ti (score), invece che da due xi ogni score è la coordinata nel nuovo

sistema di riferimento gli scores informano su come i campioni

sono in relazione l’uno con l’altro

I loadings pi sono invece le coordinate del nuovo riferimento rispetto all’originale p11 = cos(1) ; p21 = cos(2) i loadings informano su come le misure

sono in relazione l’una con l’altra

set 1set 2

set i

set 9

......

......

x1 x2

matrice X

p2p1

EptX 1 T1

si osserva che le misure non sono linearmente indipendenti

PCA (Principal Component Analysis)

Se si considerassero entrambe le componenti principali risulterebbe, esattamente (cioè senza errore residuo):

x1

x2

12

3

4

5

6

78

9

p1

EptX 1 T1 X

2

9

= t1

1

9

pT11

2

+ E

9

2

TT221 ptptX 1

L’errore o residuo E è la distanza ortogonale di ciascun punto dalla direzione della PC selezionata è quindi la quantità di informazione che si “perde” riducendo la dimensione dello spazio da 2 a 1 (cioè

considerando una sola PC invece che due) al residuo si può associare il rumore dei dati (che “corrompe” l’informazione contenuta in X), cioè tutta la

variazione non deterministica nei dati

informazione originaria

cambio coordinate

informazione “compressa”

informazione “persa”

11

PCA su uno spazio 3-D

Anche quando lo spazio di origine ha una dimensione molto elevata (es. 15 e oltre), di solito bastano solo 2 o 3 componenti principali per catturare tutta la variabilità originaria

• ••

•

•

•

•

• ••

•

•

•

•

••

•

•

•••

•

•

x1

x2

x3

p1

p2

p3

•

•

•

Supponiamo ora di avere terne di misure campionate possono essere sufficienti (ad esempio)

due sole componenti principali per descrivere tutta la variabilità delle misure

si passa da uno spazio a 3 dimensioni a uno a 2 dimensioni, ortogonali tra loro, che contiene la stessa informazione

EptptX 1 TT221

Come definire formalmente la procedura?Come valutare la quantità di informazione “persa”?

PCA: visualizzazione del caso multi-dimensionale

Matrice delle misure di processo:

X

M

K

= T

A

K

PTAM

+

KMKK

M

xxx

xx

xxx

21

2221

11211

X

misure

osservazioni

X KM

E

M

K

ETPX T

pseudo-misure (variabili latenti)

(A << M)

(M misure distinte; K osservazioni di ciascuna)

A

i

Tii

T

1

ˆ ptTPX

(approssimazione della matrice X) informazione

originaria

cambio coordinate (loadings)

informazione “compressa”

(scores)

informazione “persa”

(residuals)

12

PCA: formalizzazioneSono date K osservazioni di M variabili:

Obiettivo: approssimare X nel prodotto di due matrici, con A << M (A è il n° di componenti principali):

M variabili

K rilevazioni

KMKK

M

xxx

xx

xxx

21

2221

11211

X

X KM

T

MAM

A

KAK

AT

pp

pp

tt

tt

...

...

...

...ˆ

1

111

1

111

TPX

Matrice degli scores(colonne di T)

Matrice dei loadings(colonne di P)

(K × A ) (M × A )T

EXX ˆ

Formalmente, data la matrice X delle misure, la prima PC è quella combinazione lineare t1 = Xp1 delle M misure che ne spiega la massima varianza, col vincolo |p1|=1

N.B. I vettori colonna di P :

• sono ortogonali• sono le direzioni di

massima variabilità dei dati

KMMKKk

MM

MM

xpxpxpt

xpxpxpt

xpxpxpt

12211111

212221211121

111221111111

• I coefficienti (pi 1) sono le incognite, e rappresentano le componenti del primo vettore di loading (direzione di massima variabilità dei dati): p1 = [p11 p21 … pM1]T

• I coefficienti (tj 1) sono le proiezioni (scores) di ciascun set di dati lungo questa prima direzione principale

PCA: formalizzazione /cont.

021 M

10021

M

i

• Sommando i contributi delle varie PC, si ottiene la totale % di varianza espressa nel nuovo sistema di coordinate

• Spesso, bastano 2‒3 PC per spiegare oltre l’80% della totale variabilità iniziale, indipendentemente da M

La variabilità originaria che non è descritta dalla prima PC è: E1 = X t1p1T

cioè è la matrice residua ottenuta rimuovendo le proiezioni sulla prima PCLa seconda PC è quella combinazione lineare t2 = E1p2 della restante quantità di informazione E1 che ne spiega la massima varianza, coi vincoli di essere ortogonale (= non correlata) alla prima PC e che sia |p2|=1Si procede in questo modo fino alle A <<M componenti principali desiderate; si assume che le (M – A ) PC non utilizzate rappresentino il rumore (noise) nei dati

Come trovare i vettori pi (cioè le direzioni principali)? se S è la stima della matrice di covarianza di X , cioè S = (XT X)/(M – 1) , allora le

pi sono gli autovettori associati a ciascun autovalore i di S gli autovalori di S sono positivi, pari al n° delle misure (M ) e possono essere ordinati in

ordine decrescente (definiscono la varianza spiegata da ciascuna PC):

La porzione di variabilità dei dati originari espressa dalla i -esima PC è:

13

PCA: formalizzazione /cont.

Matrice di covarianza di X :

Decomposizione PCA di X :

I vettori di loading pi sono gli autovalori i della matrice di covarianza, cioè:

I ti sono ortogonali, mentre i pi sono ortonormali

Il vettore score ti è la combinazione lineare dei dati X definita da pi

per ogni coppia (ti ; pi ) :

i i sono una misura della frazione di variabilità descritta dalla coppia (ti ; pi )

1)cov(

M

TXXX

EptptptX TAA

TT 2211

KMKK

M

xxx

xx

xxx

21

2221

11211

X

variabili

osservazioni

(A << M)

iii ppX )cov(

ii pXt

Tornando al primo esempio1° passo: si calcola la matrice di covarianza S dei dati originari (processo

“sotto controllo”): S = (XT X)/(M – 1) si calcolano tutti gli autovettori pi e autovalori i di S si decide quanti autovettori (cioè PC) trattenere, in base alla

totale varianza spiegata

x1 x2 x3 x4

1 10 20.7 13.6 15.5

2 10.5 19.9 18.1 14.8

3 9.7 20 16.1 16.5

4 9.8 20.2 19.1 17.1

5 11.7 21.5 19.8 18.3

6 11 20.9 10.3 13.8

7 8.7 18.8 16.9 16.8

8 9.5 19.3 15.3 12.2

9 10.1 19.4 16.2 15.8

10 9.5 19.6 13.6 14.5

11 10.5 20.3 17 16.5

12 9.2 19 11.5 16.3

13 11.3 21.6 14 18.7

14 10 19.8 14 15.9

15 8.5 19.2 17.4 15.8

16 9.7 20.1 10 16.6

17 8.3 18.4 12.5 14.2

18 11.9 21.8 14.1 16.2

19 10.3 20.5 15.6 15.1

20 8.9 19 8.5 14.7

M = 4 misureK = 30 set (osservazioni)

1 2 3 4Autovalori 2.3181 1.0118 0.6088 0.0613

% varianza spiegata 58.0 25.3 15.2 1.5

% cumulativa 58.0 83.3 98.5 100

Autovettori

p1 p2 p3 p4

x1 0.59410 -0.33393 0.25699 0.68519

x2 0.60704 -0.32960 0.08341 -0.71826

x3 0.28553 0.79369 0.53368 -0.06092

x4 0.44386 0.38717 -0.80137 0.10440

Bastano 2 sole PCper catturare l’83% della varianza

Coordinate del nuovo riferimento (2-D) rispetto allo spazio (4-D) originale100

21

M

i

Varianza spiegata =

14

-4 -3 -2 -1 0 1 2 3 4 5 6 7-4.0-3.5-3.0-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.03.54.0

t 2

t1

coordinata lungo la prima PC

coor

dina

ta lu

ngo

la s

econ

da P

CTornando al primo esempio /cont.

2° passo: nota P, si calcolano le coordinate

(scores) di ciascun set di dati rispetto al nuovo riferimento: T = XP

si diagrammano i set di dati nel nuovo riferimento (diagramma degli score)

si traccia l’intervallo di fiducia

x1 x2 x3 x41 10 20.7 13.6 15.5

2 10.5 19.9 18.1 14.8

3 9.7 20 16.1 16.5

4 9.8 20.2 19.1 17.1

5 11.7 21.5 19.8 18.3

6 11 20.9 10.3 13.8

7 8.7 18.8 16.9 16.8

8 9.5 19.3 15.3 12.2

9 10.1 19.4 16.2 15.8

10 9.5 19.6 13.6 14.5

11 10.5 20.3 17 16.5

12 9.2 19 11.5 16.3

13 11.3 21.6 14 18.7

14 10 19.8 14 15.9

15 8.5 19.2 17.4 15.8

16 9.7 20.1 10 16.6

17 8.3 18.4 12.5 14.2

18 11.9 21.8 14.1 16.2

19 10.3 20.5 15.6 15.1

20 8.9 19 8.5 14.7

In pratica, per ragioni numeriche spesso si preferisce calcolare le PC per una via diversa, che non richiede la valutazione degli autovalori

t1 t2

0.2917 -0.6034

0.2943 0.4915

0.1973 0.6409

0.8390 1.4696

3.2049 0.8791

0.2033 -2.2951

-0.9921 1.6704

-1.7024 -0.3609

-0.1425 0.5608

-0.9950 -0.3149

0.9447 0.5047

-1.2195 -0.0913

2.6087 -0.4218

-0.1238 -0.0877

-1.1042 1.4726

-0.2783 -0.9476

-2.6561 0.1353

2.3653 -1.3049

0.4113 -0.2189

-2.1466 -1.1785

Tornando al primo esempio /cont.

3° passo: quando arriva un nuovo set Xnew di misure lo si proietta nel nuovo spazio,

calcolandone gli scores: Tnew = XnewP si diagrammano i punti nello score plot

e si verifica se cadono dentro o fuori l’intervallo di fiducia

x1 x2 x3 x4 t1 t2

21 9.9 20 15.4 15.9 0.0742 0.2394

22 8.7 19 9.9 16.8 -1.5176 -0.2112

23 11.5 21.8 19.3 12.1 1.4085 -0.8759

24 15.9 24.6 14.7 15.3 6.2980 -3.6740

25 12.6 23.9 17.1 14.2 3.8020 -1.9958

26 14.9 25 16.3 16.6 6.4907 -2.7314

27 9.9 23.7 11.9 18.1 2.7388 -1.3762

28 12.8 26.3 13.5 13.7 4.9587 -3.9485

29 13.1 26.1 10.9 16.8 5.6780 -3.8584

30 9.8 25.8 14.8 15 3.3697 -2.1088

-4 -3 -2 -1 0 1 2 3 4 5 6 7-4.0-3.5-3.0-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.03.54.0

t 2

t1

I punti al di fuori dell’intervallo di fiducia indicano che è avvenuto qualche cambiamento nel processo

• la variabilità delle misure non è più “naturale”, ma è dovuta a eventi anomali che stanno perturbando la marcia del processo (out-of-control)

• questi eventi vengono rilevati appena insorgono, senza aspettare che abbiano effetto sul prodotto i diagrammi dei loadings possono invece dare

indicazioni su quali variabili di processo siano maggiormente responsabili del cambiamento in atto

-4 -3 -2 -1 0 1 2 3 4 5 6 7-4.0-3.5-3.0-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.03.54.0

t 2

t1

15

“Dove” si analizza la variabilità? E di cosa?

Quindi nella PCA si cerca di descrivere come ciascuna misura (“oggetto”), tra diverse rilevazioni (“osservazioni”), varia attorno al proprio valore medio tutte le misure vengono centrate rispetto alla media e scalate a varianza unitaria per attribuire loro lo

stesso peso (effetto delle unità di misura)

La PCA è in effetti un modello di questa variabilità “naturale” (dovuta a cause “comuni”)

X KM

KMKK

M

xxx

xx

xxx

21

2221

11211

X

variabili

osservazioni

d e s c r i t t o r i

Direzione lungo la quale si dispongono i descrittori del processo

Direzione lungo la quale si vuole analizzare la variabilità

v a

r i a

b i

l i t

à

Indici diagnostici del modello PCA

• ••

•

•

•

•

• ••

•

••

••

•

•

•

••

•

•

• T 2Q

x3

p1

•p2x2

x1

•

•

Hotelling T 2 rappresenta la distanza pesata di ciascun set di

misure dalla media attesa per quel set

il test viene fatto sulla media degli scores, e non più sulle misure “reali” di processo gli scores sono distribuiti “più normalmente” delle

variabili originarie, e hanno media nulla

un valore elevato di T 2 indica che la variabilità delle misure nell’iperpiano delle prime A PC è superiore rispetto a quanto previsto in condizioni “normali” se Tk

2 > UCL, allora la variabilità nel set k di misure è maggiore di quanto possa essere spiegato da cause comuni allarme

T 2 è una misura della variabilità all’interno del modello PCA

M

Ai t

iA

i t

i

iis

t

s

tT

12

2

12

22 utile specie quando le

PC sono più di 2

nuova osservazione

16

Indici diagnostici del modello PCA /cont.

Residuo Q (squared prediction error ; SPE) misura la distanza ortogonale fra un nuovo set di

misure e la sua previsione in base al (cioè la sua proiezione sul) modello già costruito

in sostanza, corrisponde a un test di ipotesi sulla media della popolazione dei residui E ci si attende che questa media sia pari a zero

una valore elevato di Q indica che è accaduto un evento “nuovo” che non è descritto nel piano delle 2 (A ) PC usate nel modello di riferimento

Q misura la quantità di variabilità noncatturata dal modello PCA (variabilità “esterna” al modello)

M

iinewinewx xxSPEQ

1

2,, ˆ

• ••

•

•

•

•

• ••

•

••

••

•

•

•

••

•

•

• T 2Q

x3

p1

•p2x2

x1

•

•

Quindi, tipicamente, le carte di controllo statistico multivariato sono 2: una carta di T 2 sulle A PC dominanti vs. il numero di osservazione una carta di Q vs. il numero di osservazione

nuova osservazione

Il ruolo “asimmetrico” di Q e T 2L’indice Q misura la variabilità che rompe la correlazione che è “normalmente”

presente nel processo questo di solito è indice dell’esistenza di una situazione anomala (fault ) la regione “normale” definita dal limite di controllo su Q include le componenti

residue, che sono dovute principalmente a rumore anomalie di entità moderata o anche piccola possono facilmente eccedere il limite queste anomalie possono essere quindi evidenziate con molta facilità

L’indice T 2

misura la distanza dall’origine (media) nel sottospazio delle PC questo sottospazio contiene variazioni “normali” delle variabili di processo

aventi varianza elevata quindi, la regione “normale” definita da T 2 è di solito molto più ampia di quella

definita da Q l’anomalia deve essere molto significativa prima che possa essere rilevata da T 2

Un campione eccede il limite su T 2 ma non quello su Q ? allora non segnala una rottura della struttura della correlazione, ma solo un

allontanamento dall’origine nel sottospazio delle PC potrebbe essere un fault, ma anche semplicemente un cambio nella regione

operativa

17

Il ruolo “asimmetrico” di Q e T 2 /cont.

Modello PCA con una PC cerchi vuoti

“normale” esercizio stazionario triangolo

anomalia (fault ) segnalata da un valore eccessivo di Q

l’anomalia rompe il bilancio di materia ed è identificata chiaramente da Q

allo stesso tempo, T 2 è nei limiti per questa anomalia

il valore di portata è nella media attesa, ma la correlazione è variata

quadrato T 2 è oltre ai limiti “normali”,

mentre Q è nei limiti è segnalato un cambiamento

consistente col modello transizione a una regione

operativa diversa

unità

FI FI

x1 x2

Fase liquida, no reazione chimicaT = cost.

anomalia

x 2

x1

sotto

spaz

io

delle

PC

regione di T 2

cambio di

produzione

Q=SPE

45°

Altro esempio: reattore per polietilene (LDPE)

Produzione di polietilene a bassa densità reattore tubolare ad alta pressione; 3-4 sezioni (~1.5m)

etilene e solvente preriscaldati + iniziatore alla prima sezione etilene “freddo”, solvente e iniziatore alle sezioni successive ogni sezione è raffreddata con camicia

Qualità del polimero dovuta a: proprietà misurate fuori linea molto

raramente: peso molecolare; distribuzione dei

pesi molecolari; frequenza di ramificazione di catenecorte e lunghe

proprietà misurate fuori linea raramente melt index; densità

Variabili di processo misurate in linea: profilo temperatura nel reattore temperatura fluido servizio portate iniziatore e solvente … in tutto 14 misure

MacGregor et al., AIChE J., 40, 826-838 (1994)

18

Altro esempio: reattore per polietilene (LDPE) /cont.

MonitoraggioRilevazione dell’ingresso di impurezze nell’etilene alimentato

Diagramma degli scores Diagramma dei residui

MacGregor et al., AIChE J., 40, 826-838 (1994)

Altro esempio: reattore per polietilene (LDPE) /cont.

DiagnosticaRilevazione di anomalie di marcia

Contributi delle variabili di processoCause delle anomalie

sporcamento del reattore

surriscald.fluido servizio

impurezze in alimentazione

MacGregor et al., AIChE J., 40, 826-838 (1994)

19

Controllo statistico di processi batchProcesso batch: sequenza di operazioni “elementari” condotte in un intervallo di tempo assegnato

carico mescolamento riscaldamento reazione separazione raffreddamento scarico

il coordinamento tra le operazioni è ottenuto per applicazione di una ricetta

Ciascuna sequenza completa di operazioni coordinate è detta essa stessa “batch” (o “lotto”)

La qualità del prodotto è definita dal profilo nel tempo delle variabili operative, e non dai singoli valori a un certo istante di tempoLa qualità spesso non è saggiata durante un batch, ma solo alla fine end-point properties

Se il prodotto finale non è in specifica, va rilavorato, o venduto a un grado inferiore, o eliminato perdita economica

Riproducibilità dei batchSpesso da batch a batch c’è mancanza di riproducibilità effetto di disturbi condizioni iniziali variabili (materie prime; stato delle apparecchiature) mancanza di misure in linea di qualità scarsa automazione (effetto degli operatori)

È difficile legare “visivamente” le variazioni nei profili, tra un batch e l’altro, a possibili variazioni nella qualità finale del prodotto

Come distinguere un batch “normale” da uno anomalo?

Tipici profili rilevati per una temperatura in una

serie di batch diversi

García-Muñoz et al., Ind. Eng. Chem. Res., 42, 3592-3601 (2003)

20

Il problema del tempoIn un processo batch, le misure disponibili sono caratterizzate da una dimensione in più: il tempo il valore medio che, per una certa variabile, viene considerato “normale” varia ad ogni

istante di tempo

Le matrici dei dati assumono quindi una forma tri-dimensionale

J variabili

K osservazioni nel tempoIba

tch

Matrice X

È possibile studiare il problema ancora con PCA, che è un metodo bi-dimensionale?

“Srotolare” la matrice dei dati: Multiway PCA

Ciascuna variabile in ciascun istante di tempo viene riguardata come un “oggetto” (descrittore) distintoPer ciascuna variabile, si hanno tanti valori medi quanti sono gli istanti di campionamento di questa variabileNella matrice bidimensionale risultante, ciascuna riga corrisponde alla “storia” completa di un batchA questa matrice si applica una PCA “classica”

…..

…..

2

Set 1

3

200

…Si studia qual è la variabilità attraverso i batch di una certa variabile, ad un certo istante di tempo v

a r i

a b

i l i

t à

La procedura richiede che tutti i batch abbiano la stessa durata

Srotolamento nella direzione dei batch (batch-wise unfolding)

21

“Srotolare” la matrice dei dati: batch non sincroni

I profili nel tempo di ciascuna variabile in diversi batch vengono accatastati l’uno sull’altroLa media rispetto alla quale si valuta la variabilità non ha molto significato fisicoIl monitoraggio di processo tramite PCA è poco convenientePuò invece essere molto utile per effettuare stime mediante PLS

K i or H

i

cam

pion

i

Srotolamento nella direzione delle variabili (variable-wise unfolding)

Esempio: monitoraggio di un processo batchReattore industriale batch per produzione materia plastica

9 variabili di processo misurate in linea ogni 5 mindurata del batch: 1000 min (16.7 h)disponibili dati per 50 batch “normali”

FI

FI

TI

TI

TI

CI

I

TI

iMac

Acqua di raffreddamento

Stirene

Butadiene

Il problemaÈ possibile sapere in anticipo se da un batch si ottiene un prodotto in

specifica, ancor prima di misurarne la qualità finale?

Nomikos and MacGregor, AIChE J., 40, 1361-1375 (1994)

22

Monitoraggio: anomalie e diagnostica /cont.

Non è possibile riconoscere le situazioni anomale con una semplice ispezione visuale dei profili delle variabili misurate!

0 50 100 150 20044

45

46

47

48

49

50

51

Tem

pera

tura

acq

ua r

affr

edda

men

to (

°C)

Campione

0 50 100 150 2000.0080

0.0085

0.0090

0.0095

Por

tata

en

tran

te d

i but

adie

ne

Campione

temperatura H2O camicia portata butadiene

Come distinguere un batch normale da uno anomalo?

Monitoraggio: anomalie e diagnostica /cont.

Nel batch 99 è individuata una diversa struttura della correlazione tra i dati il batch no. 99 è anomalo!

Causa più probabile: variabilità “anormale” nel butadiene (F_buta)

alimentato rispetto alla variabilità media

Cambiato fornitore? Aumentato livello impurezze? Misuratore portata guasto?

0 50 100 150 2000.0080

0.0085

0.0090

0.0095

Por

tata

ent

rant

e d

i but

adie

ne

Campione

Fstir Fbuta Tf Tr Tcw Trj ro cTOT Er0

50

100

150

200

Con

trib

uto

rela

tivo

a Q

99

Variabile

F_buta

-60 -40 -20 0 20 40 60-60

-40

-20

0

20

40

60

varia

bile

late

nte

t 2

variabile latente t1

0 10 20 30 40 500

1000

2000

3000

4000

5000

6000

Res

idui

Qi

Serie di batch di riferimento e Batch 99

batch n° 99

23

PCA: riassuntoL’analisi di dati di processo col metodo PCA passa attraverso questi stadi si collezionano dati di riferimento quando dal processo esce un prodotto in

specifica la variabilità intrinseca delle condizioni di processo, quando la marcia è

nominale, viene “catturata” dal modello PCA il modello PCA fornisce quindi una “impronta” del processo e definisce un

riferimento rispetto al quale verranno valutate le condizioni future di esercizio quando nuovi dati di processo sono disponibili, vengono confrontati col

riferimento per verificare se sono consistenti con l’esercizio “normale” l’analisi degli indici di prestazione (squared prediction error Q ; Hotelling T2), unita alla

conoscenza ingegneristica del processo, permettono di rilevare eventuali anomalie e di risalire alle cause fisiche che le hanno determinate

Le conoscenze statistiche non sono sufficienti è sempre essenziale la conoscenza “a priori” del processo

quali sono le variabili di esercizio come esse influiscono sulla qualità del prodotto

Un passo indietro

24

Stime di qualità: sensori “virtuali”Colonna di distillazione ternaria batchI dati disponibili misure di temperatura sui piatti della colonna;

portate; livelli; … composizioni note solo su test run ad-hoc, e

misurate off-line (oppure da modello dettagliato)

Il problema sviluppare un sensore “virtuale” per stimare in

linea le composizioni dei prodotti e inviarle al sistema di controllo

Le difficoltà non è noto nulla sulla termodinamica del sistema o

sull’impianto

P1 S1 P2 S2

P3

Feed

TI

TI

TI

Steam

Water

.V

R.

D.

LC AC

P1 S1 P2 S2

P3

Feed

softsensor

softsensor

TI

TI

TI

Steam

Water

.V

R.

D.

LC

Esempio: sensore virtuale di composizione /cont.

I calcoli da eseguire in linea sono puramente algebrici, e quindi “rapidi” le stime possono sostituire le misure

quando queste non sono disponibili

0 200 400 600 800 1000 1200 1400 1600

50

60

70

80

90

100

110

120

Time sample

Tra

y te

mp

erat

ure

[°C

]

0 400 800 1200 16000.00.20.40.60.81.0

0.00.20.40.60.81.0

0.00.20.40.60.81.0

Validation data

Time sample

He

avy

Inte

rmed

iate

actual estimated

Ligh

t

Dalle temperature misurate in linea …

… alle composizioni stimate in linea

25

La tecnica: PLS (Partial Least-Squares regression o Projection onto Latent Structures )

Obiettivo note le misure di processo X, stimare in linea le variabili di qualità Y (di importanza

primaria, ma misurate non frequentemente) la PCA descrive la variabilità delle variabili di processo X qui non si vuole fare questo, ma piuttosto usare X per descrivere Y

Metodo “standard” immaginare una relazione lineare del tipo: Y = X B + E determinare la matrice B dei parametri del modello mediante minimizzazione dell’errore E

regressione ai minimi quadrati: B = (XTX)-1 XTY soluzione esplicita

Problema quando le misure di processo sono correlate tra loro, l’inversione di (XTX) è molto

difficile la stima dei parametri B è imprecisa e quindi la capacità di predire Y tramite la regressione è

molto modesta la soluzione è mal condizionata, cioè molto sensibile a piccole variazioni nei dati (ad esempio, in

presenza di rumore di misura)

La tecnica: PLS (projection to latent structures) /cont.

La tecnica PLS è concettualmente simile alla PCA vengono ridotte simultaneamente le dimensioni di X e di Y i vettori latenti dello spazio X sono quelli più predittivi dello spazio Y (massim.ne

della covarianza tra ti e ui)

FUQY

ETPXT

T

• I vettori score ti e ui sono ortogonali e legati tra loro da una relazione interna di linearità:

• I parametri bi del modello si determinano per regressione ai minimi quadrati

iiii b εtu

N.B. Non si tratta di due trasformazioni PCA separate la PLS cerca di trovare una diversa rappresentazione “ridotta” dei dati X in modo da dare una

migliore predizione dei dati Y cerca cioè quelle direzioni per le quali la variabilità dei dati X è meglio correlata con la

variabilità dei dati Y questo è vantaggioso specie quando non tutte le principali sorgenti di variabilità in X sono

correlate con la variabilità in Y il singolo set di dati X è rappresentato con la PLS peggio che con la PCA, ma il set di scores di

X dà una migliore rappresentazione dei dati Y (è il vantaggio rispetto alla PCR)

Per valutare le prestazioni del modello PLS si usano diagrammi e indici simili a quelli impiegati per PCA anche in questo caso si riescono a rilevare eventuali anomalie di funzionamento

prima che la qualità del prodotto venga alterata

26

Correlazione e causalità

La correlazione è fondamentalmente diversa dalla causalità la causalità implica che un cambiamento in una variabile (fattore) determina il

cambiamento in un’altra variabile (risposta) la correlazione indica che il cambiamento in una variabile è associato al cambiamento

in un’altra variabile

Il modo migliore per stabilire relazioni di causalità è condurre una serie di esperimenti pianificati (designed experiments) queste relazioni di solito sono “nascoste” nelle applicazioni multivariate a causa

dell’enorme quantità di dati disponibili

120 140 160 180 200 220 240 26040

45

50

55

60

65

70

75

80

n° a

bita

nti

n° cicogne

Oldenburg (Germania)1930-1936

Che due variabili siano tra loro correlate non significa che l’una sia causa dell’altra!

anno n° cicogne n° abitanti

1930 132 471931 141 471932 169 591933 185 631934 240 651935 248 691936 250 73

Oldenburg (Germania)

PLS: altro esempioReattore per l’ossidazione di benzene ad anidride maleica

T0 T1 Tn T31

benzene+ariabenzene+aria+anidridemaleica+sottoprodotti

z z

0 2 4 6 8 10 12 14 16 18 203.103.153.203.253.303.353.403.453.503.553.60

profilo reale profilo stimato

Ani

drid

e M

alei

ca (

mol

/min

) x1

03

tempo (min)

Stima della portata di prodotto

27

PLS: altro esempio ̶ Industria alimentare

Previsione del livello di ricoprimento di uno snack food

Modello PLS su analisi d’immagine

Yu and MacGregor, Chemom. Intell. Lab. Sys., 67, 125-144 (2003)

Sommario su PCA/PLSLa PCA proietta i dati originali di processo (X) su uno spazio di dimensione molto inferiore lo spazio è descritto dai primi A autovettori della matrice di covarianza di X viene descritta comunque la maggior parte della variabilità dei dati di partenza

La PLS trova in più le correlazioni esistenti tra le grandezze misurate in linea (X) e quelle non misurate in linea (Y)

Applicazioni principali: monitoraggio di processo rilevazione anticipata di anomalie di marcia (early detection) attribuzione di cause alle anomalie (diagnostica) “compressione” dei dati da memorizzare sensoristica “virtuale” caratterizzazione di superficie manutenzione preventiva progettazione di prodotto trasferimento di scala …

28

Riferimenti bibliograficiGeladi, P. and R. Kowalski (1986). Partial least squares regression: a tutorial. Anal. Chim. Acta, 185, 1-17. Kourti, T. and J.F. MacGregor (1995). Process analysis, monitoring and diagnosis, using multivariate projection methods. Chemom. Intell. Lab. Sys., 28, 3-21.Nomikos, P. and J.F. MacGregor (1994). Monitoring batch processes using multiway principal component analysis. AIChE J., 40, 1361-1375.Wise, B.M. and N.B. Gallagher (1996). The process chemometrics approach to process monitoring and fault detection. J. Process Control, 6, 329-348.