Controlli automatici in Python - scipy.signal e sympy

Il programma comprende "Applicazioni dei controlli automatici: SISO tool", cioè l'uso di un programma di calcolo per analizzare e sintetizzare controllori. Con Python si ottengono gli stessi risultati. Il codice qui sotto è stato eseguito e i numeri riportati sono quelli ottenuti. (I grafici si ottengono con matplotlib, per esempio plt.semilogx(w, 20*np.log10(abs(H))) per il modulo di Bode; non sono riprodotti.) L'uso più utile è la verifica di un progetto prima di consegnarlo: dopo avere scelto $C(s)$, si controlla la $\omega_A$ e il margine di fase effettivi.

Il sistema di esempio

$G(s)=\dfrac{10}{s(1+s)(1+\frac s{10})}$ (Criterio di Nyquist - poli sull'asse immaginario e guadagno variabileSe $G$ ha poli sull'asse immaginario il diagramma va all'infinito e si chiude al finito con archi di raggio infinito percorsi in senso orario: un polo di molteplicità $\mu$ in $j\omega_0$ ($\omega_0\ge0$) aggiunge un arco di ampiezza $\mu\pi$ tra $G(j(\omega_0-\epsilon))$ e $G(j(\omega_0+\epsilon))$ (per $\omega_0=0$ gli estremi sono i rami $\omega\to0^\mp$). Poi vale ancora $N=n_{G+}-n_{W+}$, con $n_{G+}$ che non conta i poli immaginari. Per $W=\frac{KG}{1+KG}$, $K\in\mathbb R\setminus{0}$, il punto critico è $-\frac1K$: si conta $N$ in ciascun intervallo dell'asse reale tra le intersezioni del diagramma; quando $-\frac1K$ è sul diagramma $W$ ha poli immaginari (caso critico).Criterio di Nyquist - poli sull'asse immaginario e guadagno variabile →). I polinomi si rappresentano come liste di coefficienti (potenze decrescenti).

import numpy as np
import sympy as sp
from scipy import signal, optimize

num = [10.0]
den = np.polymul(np.polymul([1, 0], [1, 1]), [0.1, 1])      # s (s+1) (0.1 s + 1)

Poli dell'anello chiuso e luogo delle radici

I poli di $W_K=\frac{KG}{1+KG}$ sono le radici di $d(s)+Kn(s)$ (Luogo delle radici - definizione e regole di tracciamentoPer $G=\frac{n(s)}{d(s)}$ (monici, coprimi, $n=\deg d\ge m=\deg n$) e $W_K=\frac{KG}{1+KG}$, il luogo positivo è l'insieme dei $s$ con $d(s)+Kn(s)=0$ per qualche $K>0$ (i poli di $W_K$ al variare di $K$). Regole: $n$ rami simmetrici rispetto all'asse reale, partono dai poli ($K=0$) e finiscono negli $m$ zeri o negli $n-m$ asintoti; un punto reale è nel luogo positivo se alla sua destra ci sono un numero dispari di poli+zeri reali (con molteplicità); asintoti con direzioni $\frac{(2i+1)\pi}{n-m}$ e centro $x_B=\frac{\sum p_i-\sum z_i}{n-m}$; il guadagno in un punto è $K=-\frac{d(s)}{n(s)}$.Luogo delle radici - definizione e regole di tracciamento →): per ogni $K$ si calcolano con np.roots e il luogo è l'insieme di questi punti al variare di $K$. Il valore critico si trova dove la massima parte reale si annulla:

poli = lambda K: np.roots(np.polyadd(den, K * np.array(num)))
for K in (0.1, 0.3, 1.1, 1.5):
    print(f"K={K}: poli {np.round(poli(K), 3)}")
Kc = optimize.brentq(lambda K: poli(K).real.max(), 0.5, 2.0)
print("K critico =", round(Kc, 4))

Risultato: $K=0.1$: poli $-10{,}109$ e $-0{,}446\pm0{,}889j$; $K=0.3$: $-0{,}344\pm1{,}671j$; $K=1{,}1$: $-11$ e $\pm3{,}162j$; $K=1{,}5$: $0{,}145\pm3{,}642j$ (instabile); $K_c=1{,}1$, uguale al limite di stabilità trovato con Routh e con Nyquist.

Pulsazione di attraversamento e margine di fase

$|G(j\omega)|$ si ottiene con signal.freqs(num, den, worN=[w]). $\omega_A$ è lo zero di $\log_{10}|KG(j\omega)|$ (si cerca in scala logaritmica: il modulo varia di molti ordini di grandezza), e la fase di Bode si ricostruisce dai termini (qui $-90^\circ-\arctan\omega-\arctan\frac\omega{10}$) per evitare il salto di $360^\circ$ di np.angle:

K0 = 0.1
H = lambda w: signal.freqs(K0 * np.array(num), den, worN=[w])[1][0]
wA = 10 ** optimize.brentq(lambda lw: np.log10(abs(H(10 ** lw))), -3, 3)
fase = -90 - np.degrees(np.arctan(wA)) - np.degrees(np.arctan(wA / 10))
print(f"wA = {wA:.3f} rad/s, margine di fase = {180 + fase:.1f} gradi")

Risultato: $\omega_A=0{,}784$ rad/s e $m_\varphi=47{,}4^\circ$ (Pulsazione di attraversamento, margine di fase e criterio di BodeLa pulsazione di attraversamento $\omega_A$ è la $\omega>0$ in cui $|G(j\omega)|=1$ (0 dB); il margine di fase è $m_\varphi=180^\circ+\arg G(j\omega_A)$: quanta fase si può ancora perdere prima di arrivare a $-180^\circ$. Criterio di Bode: se $G$ è strettamente propria, ha guadagno di Bode $K_B>0$, nessun polo con $\mathrm{Re}>0$ e $\omega_A$ è ben definita (esiste ed è unica), allora $W=\frac G{1+G}$ è BIBO stabile $\iff m_\varphi>0$. Più $m_\varphi$ è grande, più il sistema è lontano dall'instabilità e meno oscilla.Pulsazione di attraversamento, margine di fase e criterio di Bode →).

Risposta al gradino e parametri

signal.step(W, T=t) simula $W$ con condizioni iniziali nulle; con la definizione al 10% del corso (Risposta al gradino e suoi parametriPer $W$ strettamente propria, BIBO stabile e con $W(0)\ne0$ la risposta al gradino $w_{-1}$ parte da $0$ e tende a $W(0)$. I suoi parametri (al 10%) sono il tempo di salita $t_r=\min{t:|w_{-1}(t)-W(0)|\le0{,}1|W(0)|}$, il tempo di assestamento $t_s$ (da cui in poi la risposta resta nella banda del 10%) e la sovraelongazione $S=\frac{\max w_{-1}-W(0)}{W(0)}\cdot100%$. Se la risposta è monotona $t_r=t_s$ e $S=0$. Per $W=\frac K{s+p}$: $t_r=t_s=\frac{\ln10}p$. Uno zero instabile produce sottoelongazione (undershoot).Risposta al gradino e suoi parametri →):

W = signal.TransferFunction(K0 * np.array(num), np.polyadd(den, K0 * np.array(num)))
t = np.linspace(0, 80, 80001)
t, y = signal.step(W, T=t)
yf = y[-1]
tr = t[np.argmax(np.abs(y - yf) <= 0.1 * abs(yf))]
ts = t[np.where(np.abs(y - yf) > 0.1 * abs(yf))[0][-1] + 1]
print(f"tr = {tr:.2f} s, ts = {ts:.2f} s, S = {(y.max() - yf) / yf * 100:.1f} %")

Risultato: $t_r=2{,}14$ s, $t_s=4{,}93$ s, $S=20{,}6\%$. È coerente con la tabella di Margine di fase, banda e picco di risonanzaPer $G$ strettamente propria, senza poli instabili, con $K_B>0$ e $\omega_A$ ben definita, la banda di $W=\frac G{1+G}$ è molto vicina a $\omega_A$ ($B_p\approx\omega_A$) e il tempo di salita si stima con $t_r\approx\frac{2{,}3}{\omega_A}$. Il margine di fase controlla le oscillazioni: $|W(j\omega_A)|=\frac1{2\sin(m_\varphi/2)}$, che approssima il picco di risonanza; per $\frac{\omega_n^2}{s(s+2\xi\omega_n)}$ si ha $m_\varphi\approx100,\xi$ gradi. Specifiche tipiche: $m_\varphi\ge45^\circ$ per poche oscillazioni, $m_\varphi\approx90^\circ$ per risposta tipo primo ordine.Margine di fase, banda e picco di risonanza →: $m_\varphi=47{,}4^\circ$ corrisponde a $\xi\approx0{,}45$ e $S\approx20\%$. Nota: l'orizzonte di simulazione deve essere molto più lungo della costante di tempo più lenta; con poli $-0{,}04\pm3j$ (cioè $K=1$) servirebbero oltre $100$ s.

Diagramma di Nyquist

Il diagramma $\{G(j\omega)\}$ si ottiene valutando signal.freqs(num, den, worN=w) su una griglia fitta di $w$ (parte reale e immaginaria dei risultati; con matplotlib: plt.plot(H.real, H.imag)). L'intersezione con l'asse reale è dove la parte immaginaria cambia segno:

w = np.logspace(-2, 3, 200001)
Hn = signal.freqs(num, den, worN=w)[1]
i = np.where(np.diff(np.sign(Hn.imag)))[0][0]
print(f"asse reale in w = {w[i]:.3f}, punto {Hn.real[i]:.4f}")

Risultato: $\omega=3{,}162=\sqrt{10}$, punto $-0{,}9092=-\frac{10}{11}$ ✓ (Diagrammi di Nyquist - tracciamento, asintoti e intersezioniIl diagramma di Nyquist di $G(s)$ è la curva ${G(j\omega):\omega\in\mathbb R}$ nel piano complesso, percorsa per $\omega$ crescente; per $\omega<0$ è il simmetrico rispetto all'asse reale di quella per $\omega>0$. Si traccia combinando modulo e fase di Bode: per $\omega\to0^+$ parte da $K_B$ (se $h=0$) o da infinito con angolo $\arg K_B-90^\circ h$ (asintoto verticale $\mathrm{Re}=K_B(\sum T_i'-\sum T_i)$ per $h=1$), per $\omega\to\infty$ arriva nell'origine con angolo $-90^\circ(n-m)$. Le intersezioni con gli assi si trovano da $\mathrm{Im},G(j\omega)=0$ e $\mathrm{Re},G(j\omega)=0$.Diagrammi di Nyquist - tracciamento, asintoti e intersezioni →).

Routh con parametro (sympy)

sympy permette di tenere $K$ simbolico e ottenere direttamente gli intervalli di stabilità (Criterio di Routh-HurwitzIl criterio di Routh stabilisce se un polinomio $P(s)=p_ns^n+\dots+p_0$ è di Hurwitz senza calcolarne le radici. Si costruisce la tabella (righe $n,n-1,\dots,0$) con $\text{nuovo}=\frac{\text{pivot}\cdot a-\text{prec.}\cdot b}{\text{pivot}}$; se la tabella va a compimento, il numero di variazioni di segno nella prima colonna è il numero di radici con $\mathrm{Re}>0$, le permanenze il numero con $\mathrm{Re}<0$. Elemento nullo in prima colonna: $P$ non è Hurwitz; riga nulla: radici simmetriche (spesso immaginarie) date dal polinomio ausiliario. Con parametro $K$ nei coefficienti dà gli intervalli di stabilità.Criterio di Routh-Hurwitz →):

s, K = sp.symbols('s K')
c = sp.Poly(s * (1 + s) * (1 + s / 10) + 10 * K, s).all_coeffs()      # [1/10, 11/10, 1, 10K]
prima_colonna = [c[0], c[1], sp.simplify((c[1] * c[2] - c[0] * c[3]) / c[1]), c[3]]
print(prima_colonna)
print(sp.solve_univariate_inequality(prima_colonna[2] > 0, K))

Risultato: prima colonna $\frac1{10},\ \frac{11}{10},\ 1-\frac{10K}{11},\ 10K$; condizione $1-\frac{10K}{11}>0$ e $K>0$, cioè $0<K<\frac{11}{10}$.

Che cosa si può fare con la sola libreria standard

Con numpy e scipy si possono: fattorizzare polinomi (np.roots), costruire funzioni di trasferimento (signal.TransferFunction), calcolare risposte in frequenza (signal.freqs, per sistemi a tempo continuo), simulare (signal.step, signal.impulse, signal.lsim), moltiplicare e sommare polinomi (np.polymul, np.polyadd). Non c'è una funzione per luoghi delle radici o margini: si costruiscono come sopra. Con sympy si ottengono fratti semplici (sp.apart), trasformate (sp.laplace_transform) e limiti simbolici (sp.limit), utili per controllare i calcoli a mano.

Errori comuni

• Usare np.angle per la fase senza srotolarla o ricostruirla: salta di $\pm360^\circ$ e dà margini sbagliati.
• Cercare $\omega_A$ con un intervallo lineare: il modulo varia di ordini di grandezza, serve la scala logaritmica.
• Simulare con un orizzonte troppo corto e dedurre un $t_s$ sbagliato.
• Fidarsi del numero invece del metodo: il risultato numerico serve a verificare il progetto, la giustificazione resta quella analitica richiesta nel compito.