Capitolo 23 ) Identificazione di un sistema

Vista l’importanza della conoscenza della funzione di trasferimento del sistema da controllare per una buon progetto del regolatore, è naturale l’interesse ai metodi di identificazione, cioè alle procedure di determinazione dei parametri che definiscono il comportamento del sistema stesso.

Dato un blocco incognito, occorre innanzitutto determinarne la stabilità in anello aperto, cioè l’applicazione di un segnale impulsivo all’ingresso deve provocare una uscita che tende ad uno stato stabile (anche se con eventuali oscillazioni, che però devono essere ‘smorzate’).

Supponendo poi che, come nella maggioranza dei casi pratici, il blocco sia definibile da un ‘guadagno’ K a da una sola costante di tempo T, si tratta di ricavare i valori di questi due parametri.

Ovviamente la cosa più semplice è l’applicazione di un ingresso a gradino, con registrazione dell’andamento dell’uscita: il rapporto fra il valore a regime dell’uscita rispetto all’ampiezza del gradino d’ingresso determina K, mentre la retta tangente alla curva d’uscita nella sua origine consente di determinare T (analogamente a quanto già visto nell’esempio del forno all’inizio del capitolo 21).

Tale procedimento è semplice ma può risultare troppo grossolano per la l’intrinseca imprecisione nella valutazione dell’asintoto e della tangente.

Un metodo più preciso è la ricerca delle curva più approssimata passante per N punti: determinate sperimentalmente le coordinate di N punti della funzione di uscita si cerca la curva esponenziale che si discosta il meno possibile da questi.

La Fig. 23.1 riporta questo metodo che si basa sulla minimizzazione dell’errore quadratico medio.

Il principio è semplice: si scrive l’espressione dello scostamento fra la curva cercata e ciascuno dei punti considerati e se ne minimizza il quadrato, cioè si rende minima la somma dei quadrati di tutti gli scostamenti (questo elimina il segno e amplifica l’effetto degli errori più grandi).

La figura mostra i dati di partenza (4 punti), l’espressione del quadrato degli errori, le derivate parziali rispetto a Y (massimo asintotico) e T (costante di tempo) e il sistema di equazioni che determina questi due parametri (ottenuto uguagliando a zero le rispettive derivate parziali).

Il Mathcad^â permette la soluzione automatica del sistema, quindi consente di trovare i parametri incogniti. L’esempio è completato dal grafico che mostra i punti originali (indicati con piccoli quadrati) e la curva ottenuta con i parametri calcolati.

Fig. 23.1 - Ricerca dei parametri di un blocco dati N punti.

Quando vi è la possibilità di acquisire i dati direttamente in forma digitale (con scheda A / D direttamente collegata ad un calcolatore) e può essere quindi disponibile una grande quantità di dati da elaborare è però preferibile ricorrere ad altri metodi per l’identificazione dei parametri del sistema.

Un metodo abbastanza semplice è la valutazione delle aree nelle registrazioni in funzione del tempo dei segnali d’ingresso e di uscita di un blocco (semprechè si supponga rappresentabile con una sola costante di tempo): applicando un gradino di ampiezza u all’ingresso si avrà in uscita il segnale y ad andamento esponenziale (teoricamente y = u·(1-e^-t/T) ).

Poichè l’integrale dell’errore, per t abbastanza lungo (almeno t > 5·T ) tende al valore u·T [1]⁾, si può ricavare T semplicemente dividendo l’integrale di (u-y)·dt per u.

Fig. 23.2 - Determinazione di T col metodo dell’area.

La Fig. 23.2 mostra un esempio di applicazione di tale metodo, in cui l’integrale è stato sostituito dalla sommatoria dell’errore nei vari campionamenti. Ovviamente la precisione nella valutazione di T dipende anche dal numero di campionamenti considerato.

Ma la raccolta di centinaia di campionamenti sotto forma di file di dati consente anche di ricavare direttamente il diagramma di Bode del sistema.

Tale risultato è possibile con il calcolo degli spettri di frequenza dei due segnali (rispettivamente di ingresso e di uscita), facendo poi il rapporto dell’ampiezza di ogni frequenza d’uscita con il rispettivo valore della frequenza d’ingresso, secondo la definizione stessa di diagramma di Bode (vedi capitolo 7).

Poichè questo rapporto deve essere espresso in dB, occorre ovviamente calcolarne il logaritmo decimale e moltiplicarlo per 20 per tracciare il diagramma del modulo (spesso è sufficiente questo, trascurando il diagramma delle fasi).

Può quindi essere previsto nel calcolatore di acquisizione un particolare programma che applichi l’ FFT (o anche il DFT, se il numero di campioni non è una potenza di 2 e se si accettano tempi di calcolo più lunghi).

La Fig. 23.3 illustra il metodo, simulando la risposta ad un’onda quadra di un blocco con costante di tempo nota ( T = 40 ms) e guadagno unitario.

Come segnale d’ingresso si è scelto un’onda quadra di ampiezza unitaria e di periodo P = 4 sec, il che fornisce uno spettro con fondamentale di 0.25 Hz e con armoniche quindi multiple di questa. Si noti che data la simmetria del segnale, l’ampiezza delle armoniche pari risulta zero.

Il numero di campionamenti è stato fissato in 1024 (uno ogni circa 4 ms), il che permetterebbe uno spettro di 512 frequenze, limitato però in pratica a 60. Per la trasformazione si è utilizzata la funzione fft del Mathcad^â.

Il segnale d’uscita è stato ricavato come funzione esponenzale di T e pure trasformato con fft.

I grafici mostrano appunto l’andamento dei rispettivi segnali nel tempo ed i loro spettri.

Si è poi calcolato il modulo (G) per le armoniche dispari, interpolandolo per quelle pari, tracciando infine il diagramma di Bode.

Poichè la scala delle frequenze è in Hz e per definizione la banda passante corrisponde al valore per cui l’attenuazione è 3 dB, risulta che questa è di circa 4 Hz.

Infatti la frequenza di taglio del sistema è pari a 1/(2·p·T) = 3.97887 Hz.

La Fig. 23.4 riporta il risultato di acquisizioni reali su un impianto industriale, applicando lo stesso metodo (e la stessa forma d’onda): è evidente la frequenza di taglio a circa 1.8 Hz (T@ 88 ms), ma è anche evidente il cambiamento di pendenza a circa 4 Hz , il che comporta una seconda costante di tempo di circa 40 ms.

Come si può vedere anche in questo caso le valutazioni risultano alquanto approssimate, ma spesso sufficienti nelle applicazioni pratiche.

Fig. 23.3 - Tracciamento del diagramma di Bode dai dati acquisiti, mediante rapporto degli spettri di frequenza dei segnali d’ingresso e d’uscita.

Fig. 23.4 - Rilievo sperimentale della funzione di trasferimento di un sistema e sua interpretazione.

L’espressione della trasformata z di un sistema consente valutazioni matematicamente più esatte rispetto ai metodi visti.

Data la struttura del blocco, possono essere valutati i singoli parametri scrivendo un sistema di equazioni in cui sono sperimentalmente noti i valori di u e di y in vari istanti. Il numero di equazioni è ovviamente legato al numero di parametri incogniti.

Nel caso più semplice di una sola costante di tempo T e di un guadagno K, l’espressione risulta:

H(z) =

dove a₀ = K · DT / (T+ DT) e b₁ = T / (T+ DT)

Poichè l’equivalente espressione nel tempo è

y_t = b₁· y_t-1 + a₀ · u_t

sarebbero sufficienti due campionamenti per determinare le incognite a₀ e b₁ .

Si avrebbe cioè il sistema:

y₁ = b₁· y₀ + a₀ · u₁

y₂ = b₁· y₁ + a₀ · u₂

con y₀ , y₁ , y₂ , u₁ , u₂ noti.

Generalizzando il calcolo, si può ricorrere al metodo matriciale per la soluzione di sistemi di equazioni lineari.

Il riquadro seguente illustra l’uso delle matrici applicato al caso in esame:

In altri termini la soluzione si riduce al calcolo della ‘matrice inversa’ (M ^-1 ) degli stati del sistema per ricavare i coefficienti incogniti.

Per meglio illustrare il metodo, la Fig. 23.5 riporta un caso completo con 4 coefficienti:

Fig.23.5 - Determinazione dei coefficienti con metodo matriciale.

L’esempio simula un sistema con coefficienti assegnati per ricavarne il comportamento dinamico, cioè la risposta (y) ad un’eccitazione (u), che in questo caso è rappresentata da un impulso positivo seguito da uno negativo.

Determinati i valori nei vari istanti, così come potrebbero essere acquisiti sperimentalmente in una prova reale, questi vengono utilizzati nella matrice M per poi ricavarne la matrice inversa (il Mathcad^â la calcola automaticamente) ed infine la matrice q dei coefficienti.

Ovviamente, in assenza di disturbi o di errori di conversione (come inevitabilmente avverrebbe in una prova reale), i risultati devono coincidere con i coefficienti assegnati in partenza.

Fig. 23.6 - Ricerca parametri mediante minimizzazione degli scostamenti quadratici.

Risulta quindi chiaro che se si utilizzassero invece valori acquisiti sperimentalmente si otterrebbero i coefficienti relativi al blocco incognito.

Poichè però i valori sperimentali sono soggetti a disturbi, il calcolo dei coefficienti incogniti risulta affetto da errori e per ridurli si ricorre alla minimizzazione degli scostamenti quadratici (analogamente a quanto già visto in Fig. 23.1).

In Fig. 23.6 è riportato il metodo seguito: simulando un sistema con solo 2 parametri ma affetto da errori casuali, si ricava un numero di stati (quindi di equazioni) superiori al minimo necessario (3 anzichè 2) e , quadrando , derivando rispetto ai singoli parametri ed uguagliando a zero, si minimizzano gli errori di valutazione dei parametri stessi.

Questi vengono infatti ricavati da un sistema di 2 equazioni che contengono le sommatorie dei vari stati d’ingresso (u) e di uscita (y) del blocco .

Fig. 23.6a - Seguito della figura precedente, con l’uso del metodo matriciale per la ricerca dei parametri stimati.

La Fig. 23.6a mostra lo stesso caso, ma risolto con calcolo matriciale.

Quest’ultimo permette di pervenire ad una formula risolvente di grande semplicità formale, che può essere estesa a qualsiasi numero di parametri e, soprattutto, ad un numero di stati esaminati anche di molto superiore al minimo indispensabile, con evidenti vantaggi di pratica riduzione degli errori nella stima dei parametri, anche in presenza di errori nella rilevazione degli stati del sistema.

Fig. 23.7 - Applicazione del metodo matriciale ad un sistema con 4 parametri (e 8 equazioni degli stati).

La Fig. 23.7 ne illustra un’applicazione in cui si è eliminata l’influenza dei disturbi per mostrare l’esattezza del metodo.

Si noti che per l’effettivo svolgimento dei calcoli, il segnale di eccitazione u deve presentare ad ogni istante di campionamento una variazione rispetto al valore precedente (in questo esempio si è adottata una rampa a scalini), altrimenti la matrice degli stati presenta singolarità.

La semplicità formale dell’espressione risolvente non deve però far credere ad una semplicità di calcolo, e la figura riporta anche le matrici di passaggio per dare un’idea del grande numero di calcoli necessari ad arrivare ai risultati (soprattutto per l’inversione delle matrici) anche in casi con relativamente pochi dati.

Deve comunque risultare chiaro che tali elaborazioni sono possibili solo in ambienti di sviluppo opportuni (quali appunto il Mathcad^â) che permettono l’uso diretto del calcolo matriciale e non è quindi facilmente trasportabile in programmi applicativi.

La considerazione dei tempi di calcolo conseguenti fa capire inoltre che tali metodi sono difficilmente applicabili in ‘tempo reale’, cioè per l’identificazione ‘in linea’ dei parametri di un sistema.

Per scopi di ottimizzazione in linea, cioè per l’adattamento automatico dei parametri di regolazione in funzione della variazione degli effettivi parametri del sistema, si devono utilizzare forme ricorsive di valutazione che riducono l’ammontare dei calcoli.

Ad ogni modo, la possibilità di individuare i parametri di un sistema incognito attraverso semplici rilievi sperimentali ed utilizzando uno dei metodi illustrati in questo capitolo, consente già la possibilità di ottimizzare il progetto di regolazione: uniche assunzioni devono essere la struttura del sistema stesso (cioè il numero di parametri che si suppone lo caratterizzino) e la costanza nel tempo di tali parametri.