Einstieg in die Regelungstechnik mit Python

Fachbuch von Hans-Werner Philippsen

Kapitel 2

In den folgenden Beispielen des Kapitels 2 wird eine Einführung in die Programmierung mit Python gegeben und es werden wichtige Pakete vorgestellt.


web201

Einführung in Python

Listing 2.1: Gezeigt wird der typischer Aufbau eines Python-Programms, wobei das Programm mit Hilfe der Entwicklungsumgebung Spyder erstellt wurde. Es beginnt mit einem knappen Kommentar, Datum der Erstellung swie Autor. Gezeigt werden Datentypen, Typumwandlung, Arrays und Plotfunktionen.

# -*- coding: utf-8 -*-
"""

Einführung in Python: Typischer Aufbau eines Programms 
Datentypen, Typumwandlung, Arrays, Plotfunktionen

Created on 16.4. 2019
@author: philippsen
"""
# Zuerst werden die Pakete vollständig oder nur einzelne 
# Funktionen importiert!

import cmath
import numpy as np
import matplotlib.pyplot as plt

# Float-Variablen werden durch Zuweisung deklariert
Kp = 2.6        # z.B. Verstärkung
T1 = 0.2        # z.B. Zeitkonstante
print('Kp= 2.6 Datentyp Kp: ',type(Kp))

# Typumwandlung float nach int -> stets abrunden
K = int(Kp) 
print('K = int(Kp) = ',K,' Datentyp K: ',type(K))
Kr = round(Kp)  # Typumwandlung float nach int -> runden
print('Kr = round(Kp) = ',Kr,' Datentyp Kr: ',type(Kr))

# Komplexe Variablen ergeben sich oft im Verlauf einer Berechnung
a = cmath.sqrt(-2)
print('a = cmath.sqrt(-2) = ',a,' Datentyp a: ',type(a))

# Eine Deklaration eines Arrays, Vektors oder Matrix erfordert 
# das NumPy-Paket
N1  = np.array([T1, 1])
# Ein Array kann vorbelegt werden, hier mit aufsteigenden Zahlen
# von 0 bis 2pi mit dem Inkrement 0,01
time = np.arange(0.0, 2*np.pi, 0.01)
print('Datentyp time: ',type(time))
kk=np.size(time)
print('Anzahl Elemente time: ',np.size(time))

# Berechnung der Sinusfunktion
sinus = np.sin(time)

# Die Plotfunktionen stellt das Matplotlib-Paket zur Verfügung
# Erzeugung eines neuen Ausgabefensters mit plt.figure()
plt.figure()
# Plot der Sinus-Funktion über die Zeit:
plt.plot(time, sinus, "r") 
# Beschriftung der Achsen und des Titels:
plt.ylabel('sin t', color="red", fontsize=14) 
plt.xlabel('t [s]', fontsize=12) 
plt.title('Sinus', fontsize=12)
# Erzeugung eines Gitternetzes:
plt.grid(b=True)

web202

Berechnung der Sprungantwort

Listing 2.2: Die Standardregelstrecke des Buches wird im Zustandsraum dargestellt und die Berechnung der Sprungantwort mit dem Python-Paket SciPy durchgeführt.

# -*- coding: utf-8 -*-
"""

Standardregelstrecke im Zustandsraum Berechnung der Sprungantwort mit SciPy

Created on 14.12. 2018

@author: philippsen
"""

import scipy.signal as sig
import numpy as np
import matplotlib.pyplot as plt

a0 = 1/30.
a1 = 10/30.
a2 = 31/30.
b0 = 1/30.
A  = np.array([[0, 1.,0 ],[ 0, 0, 1. ],[ -a0, -a1, -a2]])
b  = np.array([[0.],[ 0. ],[ 1.]])
c  = np.array([b0, 0,0 ])
stand = sig.lti(A, b, c, 0)
t = np.linspace(0, 40,100)
u = np.ones_like(t)
tout, y, x = sig.lsim(stand,u,t)

plt.plot(t,y)
plt.title('Regelstrecke')
plt.xlabel('t [s]'); plt.ylabel('h(t)')
plt.grid(True)

web203

Plot mit zwei y-Achsen

Listing 2.3: Zweite y-Achse.

# -*- coding: utf-8 -*-
"""
Zweite y-Achse für die Plotfunktion in Matplotlib

Created on Fri Mar 15 17:08:01 2019
@author: philippsen
"""

import numpy as np 
import matplotlib.pyplot as plt

fig, ax1 = plt.subplots()
x = np.arange(0,10,0.1) 
ax1.plot(x, 2 * np.sin(x), "--b",lw=2) 
ax1.set_ylabel('2*Sinus (- -)', color="blue") 
plt.grid()
ax2 = ax1.twinx() 
ax2.plot(x, np.cos(x), "-.r") 
ax2.set_ylabel('Cosinus (-.)', color="red") 

ax1.set_xlabel('t [s]', fontsize=15) 

web204

3D-Plot mit Python

Listing 2.4: Die Standardregelstrecke wird mit einem PI-Regler geregelt und die Regler-Parameter werden in einem definierten Wertebereich variert (Batch-Betrieb). Die Güte wird gemäß ITAE-Kriterium berechnet und in einem 3D-Plot dargestellt. Die optimalen Parameter sind im Minimum zu finden.

# -*- coding: utf-8 -*-
"""

Standardregelstrecke mit PI-Regler, die Regler-Parameter werden in einem
Wertebereich variert (Batch-Betrieb) und die Güte gemäß ITAE-Kriterium 
in einem 3D-Plot dargestellt. Die optimalen Parameter sind im Minimum zu finden

Created on 19.7. 2022
@author: philippsen
""" 

import numpy as np
import control
import control.matlab
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

n = 300 # Anzahl Schritte
T0 = 0.1 # Abtastzeit
kk = 30  # Anzahl Parametersätze
xx = np.zeros(n)    # Regelgrößenvektor
u = np.zeros(n)     # Stellgrößenvektor
tk= np.arange(0,T0*n + T0,T0)    # k-Achse
MKp = np.zeros((kk,kk))
dKp_kk = 3/kk; dTn_kk = 6/kk
MTn = np.zeros((kk,kk))
Q = np.zeros((kk,kk)) # Güte-Matrix

# reg = control.tf(Kp*[Tn, 1], [Tn, 0])
stand = control.tf([1],[30,31, 10, 1]) # Standardregelstrecke
w = 1.0    #Sollwert
for ii in range(0, kk): # Liste von 0 bis kk-1
    for jj in range(0, kk):
        Kp = dKp_kk*(ii + 1); Tn = 5.0 + dTn_kk*(jj +1)      # Reglerparameter
        reg = control.tf([Kp*Tn, Kp], [Tn, 0])
        MKp[ii,jj] = Kp; MTn[ii,jj] = Tn
        G0 = reg*stand
        Gw = control.feedback(G0, 1) # % Regelkreis
        xk , tt = control.matlab.step(Gw,tk)
        e = w - xk
        e = e*tk
        Q[ii,jj] = sum( abs(e))*T0 #  Nachbildung Gütefunktion ITAE

print('Minimum = ',Q.min())
pa = np.where(Q == Q.min())
print('Kp = ',MKp[pa[0],pa[1]],'   Tn = ',MTn[pa[0],pa[1]] )
Kp = MKp[pa[0],pa[1]] ;    Tn = MTn[pa[0],pa[1]] 
Kp=Kp[0]; Tn=Tn[0]
R = control.tf([Kp*Tn, Kp], [Tn, 0])

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

ax.plot_surface(MKp, MTn, Q, rstride=1, cstride=1, cmap=cm.viridis)
plt.xlabel('Kp'); plt.ylabel('Tn'); ax.set_zlabel('Q')
#fig2 = plt.figure()
Gw = control.feedback(R*stand,1)
stell = control.feedback(R, stand)
x, t = control.matlab.step(Gw)
y, t = control.matlab.step(stell,t)


fig, ax1 = plt.subplots()
ax1.plot(t, x, "b") 
ax1.set_ylabel('x', color="blue", fontsize=14) 
plt.grid()
ax2 = ax1.twinx() 
ax2.plot(t, y, "r") 
ax2.set_ylabel('y', color="red", fontsize=14) 
ax1.set_xlabel('t [s]', fontsize=12) 
plt.title('Entwurf gemäß Batch Opt')

web205

Interpolation mit Python

Listing 2.5: Das Beispiel zeigt die Möglichkeit zwischen Messwerten zu interpolieren, wobei verschiedene Methoden Anwendung finden. Die Interpolation erfolgt linear, Kubisch und mit einem Polynom 5. Ordnung.

# -*- coding: utf-8 -*-
"""

Interpolation von Messwerten

Created on 2.10. 2018
@author: philippsen
"""


from scipy.interpolate import interp1d as intpo
import numpy as np
import matplotlib.pyplot as plt

x  = np.array([0,1,2,3,4,5,6,7,8,9,10])
y = np.array([0,0,0.3,0.85,0.77,1.2,1.26,1.9,2.0,2.25,2.2])
li = intpo(x,y,kind='linear')
cu = intpo(x,y,kind='cubic')
p5 = intpo(x,y,kind=5)
xx = np.linspace(0,10,100)
plt.plot(x,y,'ro')
plt.plot(xx,cu(xx))
plt.plot(xx,li(xx))
plt.plot(xx,p5(xx))
plt.xlabel('Messpunkte x')
plt.ylabel('Interpolation Messwerte')
plt.grid('true')

web206

Fast Fourier Transformation

Listing 2.6: Eine Fast Fourier Transformation erfolgt mit dem Python-Paket SciPy. Es werden drei Sinussignale mit 25, 50 und 50Hz sowie normalverteiltem Rauschen überlagert.

# -*- coding: utf-8 -*-
"""
Fast Fourier Transformation mit Python SciPy
Überlagerung von drei Sinussignalen 25, 50 und 50Hz 
sowie normalverteiltem Rauschen

Created on 1.10. 2018

@author: philippsen
"""

from scipy.fftpack import fft
import numpy as np
import matplotlib.pyplot as plt

k = 50000
TT = 10.
t  = np.linspace(0.0, TT,k)
y = 1.2*np.sin(100*np.pi*t)+ 0.8*np.sin(200*np.pi*t)+ \
    0.6*np.sin(50*np.pi*t) + np.random.normal(size=k)

ff = fft(y)
f = np.linspace(0.0,k/(2*TT),k//2)
z = 2.0/k*np.abs(ff[0:k//2])
plt.subplot(2,1,1)
plt.plot(t,y)
plt.title('25 +50 + 100 Hz + Rauschen')
plt.ylabel('y(t)')
plt.xlabel('t [s]')
plt.grid('true')
plt.subplot(2,1,2)
plt.plot(f,z)
plt.xlim(0,120)
plt.xlabel('f in Hertz')
plt.ylabel('Betrag Spektrum')
plt.grid('true')

Theme von Anders Norén