Introduzione

Desidero mostrare l'implementazione del modello matematico d-q del motore asincrono in linguaggio Python, evitando di utilizzare l'ambiente di programmazione Simulink, adottato nella quasi totalità delle università italiane.

Molteplici sono gli scopi di questo scritto:

1. mostrare un uso alternativo del linguaggio di programmazione Python;
2. approfondire la programmazione scientifica in Python;
3. suscitare l'interesse degli studenti attorno all'argomento.

Il modello della macchina, i dati numerici e le osservazioni riportate in questo scritto sono stati estrapolati dal testo in [1]. Io ho cercato solo di reinterpretarli in linguaggio Python.

Il programma, che riporta anche i dati di ingresso della simulazione, e' presentato subito dopo la breve trattazione teorica del modello matematico.

Il modello matematico della macchina asincrona

Dal punto di vista della teoria unificata delle macchine elettriche, quella asincrona presenta le maggiori difficoltà di trattazione, in quanto richiede il passaggio da avvolgimento trifase a avvolgimento equivalente lungo gli assi d-q sia per lo statore che per il rotore.

Il modello teorico, comunque, e' di semplice comprensione e implementazione.

Con i pedici D,Q e d,q si indicano le grandezze elettriche trasformate secondo tali assi, rispettivamente di statore e di rotore.

Le tensioni v_d e v_q sono nulle dato che il rotore della macchina e' a gabbia.

La coppia elettromagnetica Ce si considera positiva nel funzionamento da motore e quella resistente del carico, Cr, e' positiva se agisce in senso opposto rispetto alla velocità di rotazione.

Con ω_r si indica la velocità di rotazione elettrica.

Per qualsiasi approfondimento, si rimanda al testo citato [1] reperibile gratuitamente in rete.

Il modello della macchina e' costituito da 4 equazioni elettriche e una meccanica, di seguito riportate:

nelle quali Δ = L₁L₂ − M² e p e' la derivata rispetto al tempo.

L'equazione meccanica necessaria per completare il modello e' la seguente:

con

Queste equazioni possono essere risolte solo numericamente.

I dati di ingresso della simulazione sono riportate nel programma.

Il programma Python

In questo paragrafo riporto il codice Python che ho scritto per la risoluzione del sistema precedentemente descritto e per la creazione dei grafici.

Il programma e' migliorabile, ma scritto in questo modo pone in risalto la semplicità con cui e' possibile implementare le equazioni sopra descritte.

Per risolvere il sistema di equazioni ho utilizzato la libreria esterna SciPy. In particolare, ho fatto uso di una delle prime routine di SciPy: la scipy.integrate.odeint. Queste prime implementazioni sono compilate principalmente in Fortran (per utenti "normali" questo aspetto e' del tutto irrilevante), pertanto hanno la caratteristica di essere abbastanza efficienti e affidabili dato che sono ormai in uso da parecchi decenni.

il file

Il programma e' il seguente:

#Importo le librerie necessarie alla esecuzione del programma
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Il modello del motore asincrono
def MotoreAsincrono(z, t, R1, R2, L1, L2, M, WJ, NP, Vn, f, a, b):
	D   = L1*L2 - M**2
	W0  = 2*np.pi*f	#pulsazione
	RAD = W0*t
	
	vD  = Vn*np.cos(-RAD)
	vQ  = Vn*np.cos(np.pi/2 - RAD)
	
	piD = (-R1*L2*z[0] + z[4]*M**2*z[1] + R2*M*z[2] + z[4]*M*L2*z[3] + L2*vD)/D
	piQ = (-z[4]*M**2*z[0] - R1*L2*z[1] - z[4]*M*L2*z[2] + R2*M*z[3] + L2*vQ)/D
	pid = ( R1*M*z[0]- z[4]*M*L1*z[1]- L1*R2*z[2]- z[4]*L1*L2*z[3]- M*vD)/D
	piq = ( z[4]*M*L1*z[0] + R1*M*z[1] + z[4]*L1*L2*z[2]- L1*R2*z[3]- M*vQ)/D
	Ce  = NP*M*(z[1]*z[2]  - z[0]*z[3])
	Cr  = a + b*z[4]**2
	pWr = (Ce - Cr)*(NP/WJ)
	pz  = [piD, piQ, pid, piq, pWr]
	return pz


#Parte principale

########################


#Dati di ingresso

R1 =   1.45     # ohm - resistenza statore
R2 =   1.18     # ohm - resistenza rotore
L1 =   0.15088  # H - autoinduttanza statore
L2 =   0.15088  # H - autoinduttanza rotore
M  =   0.14324  # H - mutua induttanza statore-rotore
WJ =   0.1     # kg m2 - momento di inerzia
NP =   2        # Numero di paia di poli
Vn = 380        # V - tensione nominale
f  =  50        # Hz - frequenza nominale

#Le costanti che definiscono la coppia resistente
a  = 0.7
b  = 0
 
#Condizioni iniziali	
z0 = [0, 0, 0, 0, 0]

#Tempo di simulazione in secondi
t  = np.linspace(0, 0.8 , 4000)


#Risoluzione del modello
z  = odeint(MotoreAsincrono, z0, t, args=(R1, R2, L1, L2, M, WJ, NP, Vn, f, a, b))

#Corrente di statore sulla fase A
iA = np.sqrt(2/3)*z[:, 0]

#Coppia elettromagnetica
Ce  = NP*M*(z[:, 1]*z[:, 2]  - z[:, 0]*z[:, 3])
 

Cr  = a + b*z[:, 4]**2

#Grafici
plt.figure(1)
plt.plot(t, iA)
plt.xlabel('Tempo (s)')
plt.ylabel('Corrente statorica (A)')
plt.grid()

plt.figure(2)
plt.plot(t, Ce)
plt.xlabel('Tempo (s)')
plt.ylabel('Coppia (Nm)')
plt.plot(t, Cr)
plt.grid()

plt.figure(3)
plt.plot((z[:,4]/2)*(30/np.pi), Ce)
plt.xlabel("Velocita' (giri/min)")
plt.ylabel('Coppia (Nm)')
plt.grid()

plt.figure(4)
plt.plot(t, (z[:,4]/2)*(30/np.pi))
plt.ylabel("Velocita' (giri/min)")
plt.xlabel('Tempo (s)')
plt.grid()

plt.show() 

I grafici della simulazione, con condizioni iniziali nulle (avviamento a vuoto), sono i seguenti:

Corrente statorica della fase A

Andamento della coppia elettromagnetica (in blu) e della coppia resistente (in arancione)

Coppia elettromagnetica in funzione del velocità

Velocità' di rotazione

Conclusioni

L'articolo mostra che e' possibile risolvere un sistema di equazioni differenziali non lineari con pochissime righe di codice.

Come anticipato all'inizio, comunque, il modello della macchina asincrona e' stato solo il pretesto per far conoscere ai lettori le potenzialità di un linguaggio di programmazione come Python.

Riferimenti

1. G. Petrecca, E. Bassi, F. Benzi: La teoria unificata delle macchine elettriche rotanti, clup, 1976.
2. Travis E. Oliphant: Guide to NumPy, 2006
3. NumPy Reference 1.15, 2018.
4. SciPy Reference 1.1.0, 2018.
5. Matplotlib User's Guide 2.2.3, 2018.