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 9.1. 2026
@author: philippsen
"""
import numpy as np
import control as ct
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
stand = ct.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 = ct.tf([Kp*Tn, Kp], [Tn, 0])
MKp[ii,jj] = Kp; MTn[ii,jj] = Tn
G0 = reg*stand
Gw = ct.feedback(G0, 1) # % Regelkreis
tt, xk = ct.step_response(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 = ct.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')
plt.title('Q Batch Optimierung')
Gw = ct.feedback(R*stand,1)
stell = ct.feedback(R, stand)
t, x = ct.step_response(Gw)
t, y = ct.step_response(stell,t)
fig2, 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 Optimierung')
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')
web207
SymPy
Listing 2.7: Verwendung von SymPy
# -*- coding: utf-8 -*-
"""
Created on 8.10.24 Modellbildung Reihenschwingkreis mit R2 parallel C
Tiefsetzersatzmodell nicht lückender Strom
@author: philippsen
"""
import sympy as sy
import sympy.physics.control as co
from sympy.physics.control.control_plots import bode_plot,step_response_plot
import matplotlib.pyplot as plt
s = sy.symbols("s", complex=True)
R1, R2,C, L = sy.symbols("R1, R2, C, L", real=True, positive=True)
Xc, Xl = 1/(s*C), s*L
def par(a, b): # Funktion Paralellschaltung
return a*b/(a+b)
# RLC - Schaltung Tiefsetzersatz
G = ((par(R2, Xc)) / ( R1 + Xl + par(R2, Xc))).ratsimp()
# Bauteildaten einfügen
Gg = G.subs({R1: 0.8, R2: 20, L: 400e-6, C: 1000e-6})
# wandeln in sympy control Objekt
rlc = co.TransferFunction(*Gg.as_numer_denom(), s)
bode_plot(rlc, 1, 4,freq_unit='Hz',phase_unit='deg')
plt.figure()
step_response_plot(rlc,prec=16, upper_limit=0.01)
# step besser mit numerik Paket
Gr2 = G.subs({R1: 0.8, L: 400e-6, C: 1000e-6})
