Dal modello Simulink al codice della centralina: un controllo digitale con la trasformata Zeta, anche in Python

Questa pagina mostra come un regolatore digitale passa dal progetto al codice che gira in un microcontrollore o in una centralina (ECU) di un’automobile, seguendo il flusso di lavoro usato nell’industria: il modello in Simulink, i blocchi a tempo discreto basati sulla trasformata Zeta, la simulazione, la generazione automatica del codice C e le prove in più fasi. Lo stesso controllo è poi rifatto in Python, con librerie gratuite, e il codice è verificato. L’esempio è il controllo di velocità di un motore elettrico, ma la struttura è identica per il minimo di un motore termico, la pressione del turbo, il cruise control o la coppia di un motore di trazione.

Perché la trasformata Zeta è il linguaggio delle centraline

Una centralina non risolve equazioni differenziali: a intervalli fissi (per esempio ogni 10 ms) legge i sensori, esegue un pezzo di codice e aggiorna le uscite. Il regolatore è quindi un sistema a tempo discreto, descritto da un’equazione alle differenze, e la sua funzione di trasferimento è in (lezione e approfondimento sulla trasformata Zeta). In Simulink questo si vede direttamente: i blocchi della libreria *Discrete* (Unit Delay, che è esattamente , Discrete Transfer Fcn, Discrete-Time Integrator, Discrete PID Controller) hanno tutti un parametro *Sample time* e corrispondono uno a uno a righe di codice. Il processo fisico, invece, resta continuo (blocchi *Transfer Fcn*, *Integrator*), e il simulatore lo integra con un passo fine.

Il modello

Il motore, comandato in tensione , ha velocità con funzione di trasferimento

e il driver di potenza può dare al massimo V. Il regolatore è un PI che gira a ms.

In Simulink lo schema si costruisce così:

1. Step per il riferimento (0 → 200 rad/s a s), Sum con i segni ;
2. Discrete PID Controller con forma *Parallel*, , , *Sample time* , *Integrator method* Forward Euler (il valore predefinito);
3. nella scheda *Output Saturation* i limiti e, nella scheda *Anti-windup*, il metodo *clamping*;
4. Transfer Fcn con numeratore 20 e denominatore [0.5 1] per il motore;
5. il mantenitore non va disegnato: un blocco discreto che pilota un blocco continuo mantiene l’uscita costante tra i campioni, che è proprio lo ZOH;
6. nelle impostazioni del modello: solutore a passo fisso (per esempio ode4) con passo 1 ms, oppure a passo variabile per la sola simulazione. Per generare codice serve il passo fisso, e il regolatore deve essere un multiplo del passo base.

Dentro il blocco PI: equazioni alle differenze

Con l’integrazione di Eulero in avanti il PI è

cioè, in , . Lo stato è memorizzato in un Unit Delay: è la variabile che “ricorda” il valore del passo precedente.

L’anti-windup. Quando l’attuatore è in saturazione, l’errore resta grande e l’integratore continua a crescere: alla fine, quando l’uscita si avvicina al riferimento, l’integrale accumulato porta a una forte sovraelongazione. Il metodo *clamping* smette di integrare finché l’uscita è saturata. Nella simulazione sotto, senza anti-windup la velocità arriva a 235 rad/s (17,5 % di sovraelongazione) e si assesta in 3,2 s; con l’anti-windup non supera 200 rad/s e si assesta in 1,7 s.

Dal modello al codice

Con Embedded Coder (il generatore di codice di MathWorks per sistemi embedded) il sottosistema discreto diventa codice C. La struttura tipica, semplificata, è questa:

/* stati del modello (DWork): uno per ogni Unit Delay o integratore */
typedef struct { float I; } DW_motor_ctrl;
static DW_motor_ctrl dw;

void motor_ctrl_initialize(void) { dw.I = 0.0f; }

/* chiamata dallo scheduler ogni 10 ms */
float motor_ctrl_step(float ref, float w)
{
    float e = ref - w;
    float u_lib = 0.25f * e + dw.I;
    float u = u_lib > 12.0f ? 12.0f : (u_lib < -12.0f ? -12.0f : u_lib);
    if (u == u_lib) {                    /* anti-windup: clamping */
        dw.I += 0.5f * 0.01f * e;        /* I[k+1] = I[k] + Ki*T*e[k] */
    }
    return u;
}

Il generatore produce anche un esempio di programma principale (ert_main.c) in cui una funzione chiamata da un timer, rt_OneStep(), esegue model_step() a ogni periodo. In una centralina automobilistica la funzione di passo diventa un runnable AUTOSAR assegnato a un task periodico del sistema operativo (per esempio il task a 10 ms), che legge i segnali dal bus CAN o dai convertitori e scrive le uscite.

Dettagli importanti nella pratica:

• tempi di esecuzione: lettura dell’A/D, calcolo e scrittura del D/A occupano una parte del periodo; il comando viene applicato con un ritardo che si somma al mezzo campione del mantenitore e toglie margine di fase (approfondimento sulla stabilità nel piano z);
• aritmetica: molte centraline economiche non hanno l’unità in virgola mobile, e il codice si genera in virgola fissa, con le scale scelte per non perdere risoluzione né andare in overflow;
• più frequenze: regolatori veloci (corrente, 100 µs) e lenti (temperatura, 100 ms) convivono; i blocchi *Rate Transition* gestiscono il passaggio dei dati tra task diversi;
• sicurezza funzionale: nell’automotive il processo di sviluppo segue la norma ISO 26262, che richiede tracciabilità tra requisiti, modello, codice e prove.

Le prove: MIL, SIL, PIL, HIL

Il modello è anche il riferimento per verificare il codice, in fasi successive:

• MIL (*model in the loop*): si simula il modello completo, regolatore e processo, sul PC;
• SIL (*software in the loop*): il codice C generato, compilato per il PC, sostituisce il blocco del regolatore e i risultati devono coincidere con il MIL;
• PIL (*processor in the loop*): lo stesso codice gira sul processore vero (o su un suo simulatore), collegato al modello del processo sul PC: si misurano tempi di esecuzione e memoria;
• HIL (*hardware in the loop*): la centralina completa è collegata a un simulatore in tempo reale che riproduce il veicolo e i sensori; si provano anche i guasti, senza rischi;
• infine le prove sul veicolo.

Lo stesso controllo in Python

Python non ha un equivalente completo di Simulink, con la generazione automatica di codice certificabile per centraline. Per il progetto e la simulazione, però, bastano librerie gratuite: numpy e scipy per il calcolo, python-control per le funzioni di trasferimento, matplotlib per i grafici. Il ciclo del regolatore si scrive esplicitamente, campione per campione, ed è proprio lo stesso codice che finirà nel microcontrollore. Il processo continuo si simula con il suo modello esatto con mantenitore, con (esercizio 8 della trasformata Zeta):

import numpy as np

# processo: motore, velocità w [rad/s] comandata in tensione u [V]
Km, tau = 20.0, 0.5          # guadagno [rad/s/V] e costante di tempo [s]
T = 0.01                     # periodo del regolatore: 10 ms
a = np.exp(-T / tau)         # modello esatto con mantenitore (ZOH)

# regolatore PI discreto (integratore con Eulero in avanti, come il blocco
# Discrete PID Controller di Simulink nelle impostazioni predefinite)
Kp, Ki = 0.25, 0.5           # [V/(rad/s)], [V/(rad/s)/s]
Umax = 12.0                  # tensione massima del driver

def simula(antiwindup, t_fine=4.0):
    n = int(t_fine / T)
    w = 0.0                  # stato del processo
    I = 0.0                  # stato dell'integratore (Unit Delay)
    out = []
    for k in range(n):
        t = k * T
        rif = 200.0 if t >= 0.1 else 0.0
        e = rif - w                       # lettura del sensore e errore
        u_lib = Kp * e + I                # uscita del PI
        u = min(max(u_lib, -Umax), Umax)  # saturazione
        if not antiwindup or u == u_lib:  # anti-windup per bloccaggio
            I = I + Ki * T * e
        w = a * w + (1 - a) * Km * u      # il motore evolve per un periodo
        out.append((t, rif, w, u))
    return np.array(out)

for aw in (False, True):
    r = simula(aw)
    t, w = r[:, 0], r[:, 2]
    fuori = np.where(abs(w - 200) > 4)[0]          # campioni fuori dalla fascia del 2 %
    print(f"anti-windup={aw}: massimo {w.max():.1f} rad/s, "
          f"sovraelongazione {max(0.0, 100 * (w.max() - 200) / 200):.1f} %, "
          f"assestamento al 2 % a t = {t[fuori[-1] + 1]:.2f} s")

L’esecuzione stampa:

anti-windup=False: massimo 235.1 rad/s, sovraelongazione 17.5 %, assestamento al 2 % a t = 3.22 s
anti-windup=True: massimo 200.0 rad/s, sovraelongazione 0.0 %, assestamento al 2 % a t = 1.72 s

e i grafici di questa pagina sono stati disegnati proprio con questi dati. Per sistemi più complessi, python-control permette di collegare blocchi continui e discreti con ct.interconnect e di simularli con ct.input_output_response. Alternative gratuite a Simulink con un editor grafico a blocchi sono Xcos (incluso in Scilab) e OpenModelica; la libreria Python bdsim permette di descrivere schemi a blocchi in codice.

Come sono stati scelti i parametri. Lo zero del PI, in rad/s, cancella il polo del motore in rad/s; resta un anello del primo ordine , con costante di tempo s. Con ms ci sono dieci campioni per costante di tempo: il ritardo di campionamento costa circa 3° di margine di fase, trascurabile (approfondimento sulla stabilità nel piano z).

Gli esercizi su stabilità, luogo delle radici e margini dei sistemi digitali sono nella pagina degli esercizi collegata in cima.