Predisposizione dei controllori nei sistemi multivariabili: relè, Niederlinski, metodo sequenziale e BLT

Negli impianti di processo (colonne di distillazione, reattori, caldaie) la maggior parte dei sistemi ha più ingressi e più uscite, eppure quasi sempre si controllano con PID separati, uno per anello: è il controllo decentralizzato. È semplice, ogni operatore capisce cosa fa ogni regolatore e, se un anello va in manuale, gli altri continuano a lavorare. Il prezzo è che gli anelli si disturbano a vicenda: un PID tarato bene da solo può rendere instabile l’impianto quando anche gli altri anelli sono chiusi. Questo approfondimento mostra come si predispongono (si tarano) i PID in questo caso: l’identificazione con il relè, la verifica dell’accoppiamento con RGA e indice di Niederlinski, il metodo sequenziale, il metodo BLT con la dimostrazione del criterio e i metodi di detuning, tutto sull’esempio classico della colonna di Wood e Berry, con il codice Python che riproduce i numeri. La lezione collegata spiega i concetti di accoppiamento, RGA e disaccoppiamento.

Il tema è stato oggetto della tesi di laurea di Eliana Li Santi all’Università di Pisa (relatori Aldo Balestrino e Alberto Landi, 2004), che confronta i metodi in simulazione: i capitoli sono in fondo alla pagina.

L’esempio: la colonna di Wood e Berry

È il processo di prova più usato nella letteratura sul controllo multivariabile: una colonna di distillazione metanolo-acqua in cui si controllano la purezza del distillato e quella del fondo agendo sulla portata di riflusso e sul vapore al ribollitore . Il modello identificato (tempi in minuti) è

Ogni elemento è un primo ordine con ritardo: ogni ingresso agisce su entrambe le uscite, con ritardi diversi.

Accoppiamento. Il guadagno relativo dell’abbinamento , è

Positivo, quindi l’abbinamento è quello giusto, ma lontano da 1: chiudere l’anello 2 dimezza il guadagno che l’anello 1 vede (). Un PID tarato sull’anello 1 da solo risulterà circa il doppio più aggressivo del previsto, se l’altro è chiuso.

Indice di Niederlinski. Per un sistema con regolatori a azione integrale su ogni anello,

Teorema (Niederlinski). Se , con quell’abbinamento il controllo decentralizzato con azione integrale non può funzionare in sicurezza: o il sistema è instabile con tutti gli anelli chiusi, o diventa instabile quando uno degli anelli viene aperto (messo in manuale, o saturato). I regolatori hanno guadagni dello stesso segno dei rispettivi , come sempre nella pratica.

Dimostrazione (caso dei regolatori lenti). Si considerino regolatori puramente integrali molto lenti, con piccolo: è la taratura più prudente possibile. I poli lenti dell’anello chiuso sono vicini all’origine, dove ; l’equazione caratteristica diventa , con , quindi i poli lenti sono , con gli autovalori di . Il prodotto degli autovalori è . Poiché ha il segno di , il prodotto ha il segno di e quindi ha il segno di . Se il determinante è negativo: una matrice reale con determinante negativo ha un numero dispari di autovalori reali negativi, quindi almeno un e un polo . Il sistema è instabile anche con la taratura più lenta. Il risultato generale, per qualunque taratura, è di Grosdidier, Morari e Holt (1985). ∎

Per la colonna : l’abbinamento non è escluso. Per un sistema 2×2 vale sempre , e qui infatti .

Identificazione con il relè

Per tarare un PID con le regole di Ziegler-Nichols servono il guadagno critico e il periodo critico . Portare un impianto al limite di stabilità aumentando il guadagno è pericoloso; Åström e Hägglund (1984) hanno proposto di sostituire il regolatore con un relè: l’uscita del regolatore commuta tra e a seconda del segno dell’errore. L’anello entra in un’oscillazione stabile di ampiezza controllata (si sceglie ), da cui si ricavano e . È la base della funzione di autotaratura presente in quasi tutti i regolatori industriali.

Teorema (funzione descrittiva del relè). Se l’uscita oscilla in modo quasi sinusoidale con ampiezza , allora

Dimostrazione. L’uscita del relè è un’onda quadra di ampiezza . La sua serie di Fourier è . Il processo è un filtro passa-basso: le armoniche superiori sono attenuate e conta solo la prima, di ampiezza . Il relè si comporta quindi come un guadagno equivalente (rapporto tra prima armonica dell’uscita e ampiezza dell’ingresso), in fase. Un’oscillazione si autosostiene quando : la fase di deve essere ° e . Ma questa è proprio la definizione della frequenza critica e del guadagno critico: . ∎

Quanto è precisa? Simulando la prova sui due anelli della colonna (codice più sotto) e confrontandola con il calcolo esatto dalla risposta in frequenza:

	esatto	relè	esatto	relè
anello 1	2,099	1,709	3,91 min	3,89 min
anello 2	0,422	0,349	11,13 min	10,97 min

Il periodo è quasi esatto; il guadagno è sottostimato del 17-19 %, perché con un ritardo piccolo rispetto alla costante di tempo l’uscita è più triangolare che sinusoidale (lo si vede nel grafico) e l’ipotesi della prima armonica è approssimata. L’errore va nella direzione prudente: si ottiene un regolatore un po’ meno aggressivo. I regolatori commerciali correggono la stima o aggiungono un’isteresi al relè per renderla robusta al rumore.

Ziegler-Nichols anello per anello

Con e esatti, le regole di Ziegler-Nichols per il PI (, ) danno

• anello 1: , min;
• anello 2: , min (negativo perché ).

Ciascuno, da solo, funziona. Ma chiudendo entrambi gli anelli il sistema diventa instabile: nella simulazione l’oscillazione dopo il gradino di riferimento non si spegne e la sua ampiezza cresce lentamente (0,248, 0,257, 0,266 in tre intervalli successivi di 100 minuti). È la conseguenza di : ciascun anello vede un processo diverso da quello su cui è stato tarato.

Il metodo sequenziale

Si tarano gli anelli uno alla volta, chiudendoli in ordine: prima l’anello più veloce (qui l’anello 1, min), con la prova a relè e le regole di taratura; poi, con l’anello 1 chiuso, si fa la prova a relè sull’anello 2, che vede così il processo “efficace” che comprende l’effetto dell’altro anello. Se ci sono più anelli si continua, e si può ripetere il giro finché i parametri non cambiano più. Vantaggio: ogni prova a relè tiene conto delle interazioni reali dell’impianto. Svantaggi: il risultato dipende dall’ordine, e il primo anello, tarato quando gli altri erano aperti, può dover essere ritarato.

Il metodo BLT

Luyben (1986) ha proposto un metodo semplice e molto usato, chiamato BLT (biggest log-modulus tuning):

1. Si calcolano le impostazioni di Ziegler-Nichols di ciascun anello, e .

2. Si sceglie un unico fattore di detuning e si pongono , : tutti gli anelli diventano più lenti nella stessa proporzione.

3. Si calcola la funzione

4. Si sceglie in modo che il massimo di sia dB, con il numero di anelli.

Perché funziona. Per il teorema di Nyquist multivariabile, il sistema ad anelli chiusi è stabile se il diagramma di non circonda l’origine (con processo e regolatori stabili ad anello aperto). è definita in modo che : la condizione diventa “ non circonda “, come nel Nyquist di un solo anello con funzione d’anello . Allora gioca il ruolo della funzione di sensitività complementare, e il suo picco (il modulo di risonanza in anello chiuso) misura quanto il diagramma si avvicina a . Per un solo anello la regola pratica è un picco di +2 dB; per anelli Luyben usa dB perché contiene i prodotti di funzioni d’anello e cresce più in fretta. ∎

Per la colonna, con si cerca un picco di 4 dB:

	1	2	2,5	3
picco di	50,1 dB	6,5 dB	4,2 dB	2,6 dB

Con (Ziegler-Nichols senza correzione) il picco è altissimo: il diagramma passa quasi su , coerente con l’oscillazione che non si smorza. Risolvendo per 4 dB si trova (Luyben lo arrotonda a 2,55), cioè

• anello 1: , min;
• anello 2: , min,

gli stessi valori pubblicati da Luyben.

Nella simulazione del gradino di riferimento su , con le tarature BLT ha una sovraelongazione del 10 % e un disturbo massimo di 0,67; l’integrale dell’errore assoluto (IAE, somma sulle due uscite in 120 minuti) è 19,7, contro 76,3 di Ziegler-Nichols, che però continuerebbe a crescere perché il sistema è instabile.

Altri metodi di detuning

Il fattore unico del BLT è prudente ma non ottimale: rallenta allo stesso modo anche l’anello che non ne avrebbe bisogno. I metodi successivi, confrontati anche nella tesi citata, cercano compromessi migliori:

• detuning con criteri integrali: si sceglie minimizzando un indice come IAE o ISE calcolato in simulazione sul sistema completo, invece che con il picco di ;
• detuning con fattori diversi: un fattore per ogni anello, , scelti con lo stesso criterio; gli anelli più accoppiati vengono rallentati di più;
• detuning iterativo: si parte da e lo si aumenta a passi sull’impianto vero, finché la risposta non è accettabile; non richiede il modello ma solo prove;
• integrazione condizionale: è l’anti-windup. Quando l’uscita di un regolatore satura (una valvola tutta aperta), l’integrale smette di accumulare; senza questa precauzione, nei sistemi multivariabili un anello in saturazione trascina gli altri in lunghe oscillazioni.

Nella pratica dell’industria la sequenza tipica è: RGA e Niederlinski per scegliere gli abbinamenti, prova a relè per ogni anello, BLT per una taratura sicura di partenza, ritocchi sull’impianto. Se l’accoppiamento è troppo forte ( molto diverso da 1) si passa ai disaccoppiatori o al controllo predittivo (MPC), descritti nella lezione.

Il codice Python

Il programma calcola RGA e Niederlinski, il guadagno e il periodo critici esatti, simula la prova a relè, trova il fattore BLT:

import numpy as np
from scipy.optimize import brentq

# colonna di Wood-Berry: G_ij = K e^{-theta s} / (tau s + 1), tempi in minuti
K = np.array([[12.8, -18.9], [6.6, -19.4]])
TAU = np.array([[16.7, 21.0], [10.9, 14.4]])
TH = np.array([[1.0, 3.0], [7.0, 3.0]])

def G(s):
    return K * np.exp(-TH * s) / (TAU * s + 1)

l11 = 1 / (1 - K[0, 1] * K[1, 0] / (K[0, 0] * K[1, 1]))
NI = np.linalg.det(K) / (K[0, 0] * K[1, 1])
print(f"RGA lambda11 = {l11:.3f}, Niederlinski = {NI:.3f}")

# guadagno e periodo critici esatti (fase = -180 gradi) e regole di Ziegler-Nichols per il PI
ZN = []
for i in range(2):
    fase = lambda w: -TH[i, i] * w - np.arctan(TAU[i, i] * w) + np.pi
    wu = brentq(fase, 1e-3, 10)
    Ku = np.sqrt(1 + (TAU[i, i] * wu) ** 2) / abs(K[i, i])
    Pu = 2 * np.pi / wu
    ZN.append((np.sign(K[i, i]) * Ku / 2.2, Pu / 1.2))
    print(f"anello {i+1}: Ku = {Ku:.3f}, Pu = {Pu:.2f} min, PI: Kc = {ZN[i][0]:.3f}, TI = {ZN[i][1]:.2f}")

# prova a relè simulata (d = 1, Eulero con passo 2 ms, ritardo come coda di campioni)
def rele(i, d=1.0, dt=0.002, T=150):
    n, nd = int(T / dt), int(round(TH[i, i] / dt))
    u_coda, y, ys = np.zeros(n + nd), 0.0, np.zeros(n)
    for k in range(n):
        u_coda[k + nd] = -d * np.sign(K[i, i]) if y > 0 else d * np.sign(K[i, i])
        y += dt * (-y + K[i, i] * u_coda[k]) / TAU[i, i]
        ys[k] = y
    coda = ys[n // 2:]
    a = (coda.max() - coda.min()) / 2
    salite = np.where(np.diff(np.sign(coda)) > 0)[0]
    return 4 * d / (np.pi * a), np.mean(np.diff(salite)) * dt

for i in range(2):
    Ku_r, Pu_r = rele(i)
    print(f"relè anello {i+1}: Ku = {Ku_r:.3f}, Pu = {Pu_r:.2f} min")

# BLT: picco del modulo logaritmico di Luyben uguale a 2N = 4 dB
def picco_Lc(F):
    picco = -np.inf
    for w in np.logspace(-3, 1, 4000):
        s = 1j * w
        Gc = np.diag([ZN[j][0] / F * (1 + 1 / (F * ZN[j][1] * s)) for j in range(2)])
        W = -1 + np.linalg.det(np.eye(2) + G(s) @ Gc)
        picco = max(picco, 20 * np.log10(abs(W / (1 + W))))
    return picco

F = brentq(lambda F: picco_Lc(F) - 4.0, 1.2, 5)
print(f"BLT: F = {F:.3f}")
for j in range(2):
    print(f"anello {j+1}: Kc = {ZN[j][0] / F:.3f}, TI = {ZN[j][1] * F:.2f} min")

L’esecuzione stampa , , i valori critici e di Ziegler-Nichols della tabella, e le tarature BLT. In MATLAB lo stesso calcolo si fa con tf con la proprietà InputDelay, allmargin per i valori critici e freqresp per .

La tesi originale

I capitoli della tesi di Eliana Li Santi, *Predisposizione dei controllori nei sistemi multivariabili*, Università di Pisa, 2004, con le implementazioni e i confronti in simulazione:

• Frontespizio e sommario
• Capitolo 1: modello del PID, identificazione con il relè
• Capitolo 2: metodo sequenziale
• Capitolo 3: BLT
• Capitolo 4: detuning con criteri integrali, robustezza al rumore, Nyquist multivariabile
• Capitolo 5: detuning con integrazione condizionale, saturazione
• Capitolo 6: detuning con fattori diversi
• Capitolo 7: detuning iterativo
• Conclusioni e bibliografia