Digitale Systeme werden mit Python und der Python Control Toolbox berechnet und simuliert.
web801
Diskretisierung System 2.Ordnung
Listing 8.1: Die Sprungantwort eines Systems 2. Ordnung wird in eine Matrix überführt und die Koeffizienten der DZG berechnet
# -*- coding: utf-8 -*-
"""
Diskretisierung System 2. Ordnung
Created on 5.11. 2018
@author: philippsen
"""
import numpy as np
x1 = 0.0185330; x2 = 0.0645967
x3 = 0.1270349; x4 = 0.19799
y1 = 1.0; y2=y1; y3 = y2; y4=y3
M = np.array([[0, 0, y1, 0],
[-x1, 0, y2, y1],
[-x2, -x1, y3, y2],
[-x3, -x2, y4, y3]])
h = np.array([[ x1 ],[ x2], [x3],[x4]])
par = np.linalg.solve(M,h)
print(par)
web802
PI-Regler-Algorithmus
Listing 8.2: Vorgestellt wird ein einfacher Algorithmus für den PI-Regler in Python.
# -*- coding: utf-8 -*-
"""
# Algorithmus PI-Regler und Simulation Sprungantwort
Created on 29.10. 2018, mod 20.8.24
@author: philippsen
"""
import numpy as np
#import control as ct
import matplotlib.pyplot as plt
# Parameter zeitkontinuierlicher Regler
kp = 2; tn = 1.5
t0 = 0.1 # Abtastzeit
# Diskretisierung
d0 = kp*(t0/(tn*2) + 1); d1 = kp*(t0/(tn*2) - 1)
# Eingangsgröße
e = 1; y = np.zeros(21)
e1 = 0 # Setzen der Anfangswerte
# Berechnung der Stellgröße
for i in range(1, 21):
y[i] = y[i-1] + d0*e + d1*e1
y1 = y[i]
e1 = e
t=np.arange(0,2+t0,t0)
plt.step(t,y,where='pre')
plt.xlabel('k*T0'); plt.ylabel('y(k)')
plt.grid(True)
web803
PI-Algorithmus mit ARW-Maßnahme.
Listing 8.3: PI-Algorithmus gemäß Bild 5.6. Mit der Funktion forced_reponse wird die Reaktion des DC-Motors für einen Abtastschritt berechnet.
# -*- coding: utf-8 -*-
"""
PI-Regler mit verbesserter AWR I-Rückführung
Created on Thu Apr 25 15:18:58 2019, mod 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
def PI_AWR_REGLER(ek,yminmax,Kp):
""" PI-Regler mit verbesserter AWR I-Rückführung """
global ykk , yk, Kid, Tikorr
global yi1, ek1
Intkorr = (ykk - yk)/Tikorr # Korrektur I-Anteil falls Begrenzung
yi = yi1 + Kid*ek + Kid*ek1 - Intkorr # Berechnung I-Anteil
ek1 = ek
yi1=yi
ykk = Kp*ek + yi # Stellgroesse = P + I
if ykk < -yminmax : # Begrenzung Stellgroesse
yk = -yminmax
elif ykk > yminmax :
yk = yminmax
else:
yk = ykk # y ohne Begrenzung
return yk
# Sim mit forced response und digitalem PI-Regler
n = 500 # Anzahl Schritte
T0 = 0.01 # Abtastzeit
xx = np.zeros(n) # Regelgrößenvektor
u = np.zeros(n) # Stellgrößenvektor
yyi = np.zeros(n) # I-Anteil Stellgröße
tk= np.arange(0,T0*n,T0) # k-Achse
DCmo = ct.tf([5],[1.2, 1]) # Motormodell 1. Ordnung
Kp = 2; Tn = 1.4 # Reglerparameter
Kid=Kp*T0/(2*Tn) # Ki diskretisiert
Tikorr = 0.96*Tn/T0 # Korrekturfaktor I-Anteil
ykk =0; yk =0.0; yi1 = 0.0; ek1 = 0.0; Xvor = 0 # Anfangswerte
w = 20.0 #Sollwert
for i in range(1, n):
e = w - xx[i-1]
u[i] = PI_AWR_REGLER(ek=e,yminmax=10,Kp=Kp)
yyi[i]=yi1
tt, xxk , xk = ct.forced_response(DCmo,[tk[i], tk[i]+T0],
[u[i], u[i]], X0=Xvor,return_x=True)
xx[i] = xxk[1]
Xvor = xk[0,1]
plt.plot(tk,xx, tk, u, tk,yyi)
plt.xlabel('t [s]'); plt.ylabel('x, y, yi')
plt.grid(True)
MicroPython Listing 8.3: PI-Algorithmus gemäß Bild 5.6 für das Pyboard D
# PI-Regler auf dem pyboard d
import os
import machine
import time
from pyb import Pin, ADC, DAC
def PI_AWR_REGLER(ek,yminmax,Kp):
""" PI-Regler mit verbesserter AWR I-R鐪塩kf鐪塰rung """
global ykk , yk, Kid, Tikorr
global yi1, ek1
Intkorr = (ykk - yk)/Tikorr # Korrektur I-Anteil falls Begrenzung
yi = yi1 + Kid*ek + Kid*ek1 - Intkorr # Berechnung I-Anteil
ek1 = ek
yi1=yi
ykk = Kp*ek + yi # Stellgroesse = P + I
if ykk < 0 : # Begrenzung Stellgroesse
yk = 0
elif ykk > yminmax :
yk = yminmax
else:
yk = ykk # y ohne Begrenzung
return yk
T0 = 0.2
Kp = 0.9; Tn = 2.4 # Reglerparameter
Kid=Kp*T0/(2*Tn) # Ki diskretisiert
Tikorr = 0.96*Tn/T0 # Korrekturfaktor I-Anteil
ykk =0; yk =0.0; yi1 = 0.0; ek1 = 0.0; Xvor = 0 # Anfangswerte
adc = ADC(Pin('X1')) # read value, 0-4095
d = DAC(Pin('X5'), bits=12)
pin = machine.Pin('LED_GREEN',Pin.OUT)
os.remove('plog1.txt')
logfile = open('plog1.txt', 'w')
for i in range(100):
start = time.ticks_ms() # get millisecond counter
pin.value(not pin.value())
x = adc.read()*3.29/4095 # noch kalibrieren
e = 1.0 - x
yk = PI_AWR_REGLER(ek=e,yminmax=3.0,Kp=Kp)
d.write(int(yk*4095/3.29))
time.sleep(T0)
delta = time.ticks_diff(time.ticks_ms(), start) # compute time difference
print('T0 [ms] = ',delta,'x = ',x ,' y = ', yk )
s = str(T0*i) + '\t' + str(x) + '\t' + str(yk) + '\n'
logfile.write(s)
logfile.close()
pin.value(False)
d.write(0)
MicroPython Listing 8.3: PI-Algorithmus gemäß Bild 5.6 für den ESP32
# PI-Regler auf dem ESP32 WROOM Entwicklerbord
# programmiert in microPython
import os
import machine
import time
def PI_AWR_REGLER(ek,yminmax,Kp):
""" PI-Regler mit verbesserter AWR I-Rückführung """
global ykk , yk, Kid, Tikorr
global yi1, ek1
Intkorr = (ykk - yk)/Tikorr # Korrektur I-Anteil falls Begrenzung
yi = yi1 + Kid*ek + Kid*ek1 - Intkorr # Berechnung I-Anteil
ek1 = ek
yi1=yi
ykk = Kp*ek + yi # Stellgroesse = P + I
if ykk < 0 : # Begrenzung Stellgroesse
yk = 0
elif ykk > yminmax :
yk = yminmax
else:
yk = ykk # y ohne Begrenzung
return yk
T0 = 0.2
Kp = 1.5; Tn = 5.2 # Reglerparameter
Kid=Kp*T0/(2*Tn) # Ki diskretisiert
Tikorr = 0.96*Tn/T0 # Korrekturfaktor I-Anteil
ykk =0; yk =0.0; yi1 = 0.0; ek1 = 0.0; Xvor = 0 # Anfangswerte
pin = machine.Pin(2,machine.Pin.OUT)
ana = machine.ADC(machine.Pin(32))
ana.atten(machine.ADC.ATTN_11DB)
ana.width(machine.ADC.WIDTH_12BIT)
d=machine.DAC(machine.Pin(25))
d.write(0) # analogen Ausgang auf null setzen
#os.remove('elog3.txt')
logfile = open("elog3.txt", "w")
for i in range(100):
start = time.ticks_ms() # get millisecond counter
pin.value(not pin.value())
x = ana.read()*3.6/4095 # noch kalibrieren
e = 1.0 - x
yk = PI_AWR_REGLER(ek=e,yminmax=3.0,Kp=Kp)
d.write(int(yk*255/3.13))
time.sleep(T0)
delta = time.ticks_diff(time.ticks_ms(), start) # compute time difference
print('T0 [ms] = ',delta,'x = ',x ,' y = ', yk )
s = str(T0*i) + '\t' + str(x) + '\t' + str(yk) + '\n'
logfile.write(s)
logfile.close() # Datei schließen
pin.value(False) # Led ausschalten
d.write(0) # analogen Ausgang auf null setzen
# mit ampy -p com3 -d 2 get elog1.txt > elog1.txt hochladen
# disconnect nicht vergessen
# mit f1 = np.loadtxt('elog1.txt') unter Python in Numpy-Array wandeln
# plt.plot(f1[:,0],f1[:,1],f1[:,0],f1[:,2]) und plotten
web804
Hidden Oscilations
Listing 8.4: Python-Code zur Erzeugung von verdeckten Schwingungen.
# -*- coding: utf-8 -*-
"""
versteckte Schwingungen
Created on 29.10. 2018, mod 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
# Zeitdiskretes System 1. Ordnung
db1 = 0.0; db0 = 0.2; da1 = 1; da0 = 0.8
t0 = 0.1 # Abtastzeit
dz = np.array([db1, db0])
dn = np.array([da1, da0])
diskSys = ct.tf(dz,dn,t0)
t=np.arange(0,2,t0)
t1, y = ct.step_response(diskSys,t)
plt.step(t1,y,where='post')
# für Version < 0.9 plt.step(t1,y[0],where='post')
plt.xlabel('k*T0'); plt.ylabel('x(k)')
plt.grid(True)
web805
Diskretisierung gemäß Tustin
Listing 8.5: Ein P-T1 wird mit der Python Control Toolbox gemäß Tustin diskretisiert und die Sprungantwort des diskreten und analogen Systems geplottet.
# -*- coding: utf-8 -*-
"""
Diskretisierung P-T1 mit Toolbox und Tustin
Created on 2.4. 2019, mod. 20.8.24
@author: philipp
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
den = np.array([3, 1])
pt1 = ct.tf(1,den)
t, y = ct.step_response(pt1)
plt.plot(t,y)
t0 = 1 # Abtastzeit
diskPT1=ct.c2d(pt1,t0,'tustin')
tt=np.arange(0,t.max(),t0)
t1, y = ct.step_response(diskPT1,tt)
plt.step(t1,y,where='post')
plt.xlabel('t [s] und k'); plt.ylabel('h(t) und h(k)')
plt.grid()
web806
Diskretisierung Notch-Filter
Listing 8.6: Die Diskretisierung erfolgt mit der in Kapitel 8.9 angegebenen Formel undim Vergleich mit der c2d-Funktion der Python Control Lib.
# -*- coding: utf-8 -*-
"""
Diskretisierung mittels bilinearer Transformation
System 2. Ordnung mit der Übertragungsfunktion
b2s^2 + b1s + b0 db2 + db1z + db0z^2
G(s) = ------------------- G(z) = ---------------------
a2s^2 + a1s + a0 da2 + da1z + z^2
Created on 29.10. 2018, mod. 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
b2 = 1; b1 = 0.01; b0 = 1; a2 = 1; a1 = 1; a0 = 1
t0 = 0.01 # Abtastzeit
db2 = (b0 - 2*b1/t0 + 4*b2/(t0*t0))/(a0 + 2*a1/t0 + 4*a2/(t0*t0))
db1 = (2*b0 - 8*b2/(t0*t0))/(a0 + 2*a1/t0 + 4*a2/(t0*t0))
db0 = (b0 + 2*b1/t0 + 4*b2/(t0*t0))/(a0 + 2*a1/t0 + 4*a2/(t0*t0))
da2 = (a0 - 2*a1/t0 + 4*a2/(t0*t0))/(a0 + 2*a1/t0 + 4*a2/(t0*t0))
da1 = (2*a0 - 8*a2/(t0*t0))/(a0 + 2*a1/t0 + 4*a2/(t0*t0))
dz = np.array([db0, db1, db2])
dn = np.array([1, da1, da2])
dnotch = ct.tf(dz,dn,t0)
#t=np.arange(0,4,t0)
#t1, y = ct.step_response(dnotch,t)
#plt.step(t1,y[0],where='post')
notch2 = ct.tf([b2, b1, b0],[a2, a1, a0])
dnotch2 = ct.sample_system(notch2,t0,method='tustin') # Diskretisierung
# Vergleich Formel mit Python-Funktion
print('Diskretisierung mit Formel = ',dnotch)
print('Diskretisierung mit Python Funktion = ',dnotch2)
web807
Trapez-Profil Motion Control
Listing 8.7: Der Python Algorithmus berechnet das Trapez-Profil der Geschwindigkeit und die Fahrkurve für die Position.
# -*- coding: utf-8 -*-
"""
Trapez Profil v, Fahrkurve Weg
Created on 16.10. 2018
@author: philippsen
"""
import numpy as np
import matplotlib.pyplot as plt
n = 400; s0 = 10.0; TH = 10; TVS = 2.5
v0 = s0/(TH - TVS)
T0 = TH/n
x = np.zeros(n+1); v = np.zeros(n+1) ; a = np.zeros(n+1); r = np.zeros(n+1)
tt = np.linspace(0,TH,n+1)
# (* Anfangswerte setzen*)
s1 = 0; s2 = 0; s3 = 0; str = 0
for k in range(1,n+1):
tk = k*T0
if tk <= TVS :
s1 = tk*tk*s0/(2*TVS*(TH - TVS))
elif tk <= (TH - TVS) :
s2 = (tk - TVS)*v0
elif tk <= TH :
s3 = (TVS*TVS - (tk - TH)*(tk - TH))*s0/(2*TVS*(TH - TVS))
x[k]= s1 + s2 + s3
for k in range(1,n):
v[k] = (x[k+1]-x[k])/T0 # ableiten
for k in range(1,n):
a[k] = (v[k+1]-v[k])/T0
plt.subplot(3,1,1)
plt.plot(tt,x)
plt.ylabel('Weg x'); plt.grid(True)
plt.subplot(3,1,2)
plt.plot(tt,v)
plt.ylabel('Geschwindigkeit'); plt.grid(True)
plt.subplot(3,1,3)
plt.plot(tt,a)
plt.ylabel('Beschleunigung'); plt.grid(True)
plt.xlabel('Zeit in Sek'); plt.grid(True)
#p= a*v #Leistung
#plt.plot(tt,p)
#plt.plot(a,v)
#from scipy.fftpack import fft
#ff = fft(a)
#z = 2.0/n*np.abs(ff[0:n//2])
#plt.plot(z)
web808
Sin^2-Profil Motion Control
Listing 8.8: Berechnet wird das Sinusquadrat-Profil der Geschwindigkeit. Der zurückgelegte Weg, die Beschleunigung und der Ruck werden außerdem in einem Diagramm geplottet. In einem weiteren Fenster wird der Leistungsverlauf dargestellt.
# -*- coding: utf-8 -*-
"""
sin^2 Profil
Created on 12.4. 2019
@author: philippsen
"""
from math import sin, pi
import numpy as np
import matplotlib.pyplot as plt
n = 200; s0 = 10.0; TH = 10; TA = 3.5
v0 = s0/(TH - TA)
w = pi/(2*TA)
T0 = TH/n
x = np.zeros(n+1); v = np.zeros(n+1) ; a = np.zeros(n+1); r = np.zeros(n+1)
tt = np.linspace(0,TH,n+1)
# (* Anfangswerte setzen*)
s1 = 0; s2 = 0; s3 = 0; str = 0
for k in range(1,n+1):
tk = k*T0
if tk <= TA :
s1 = (tk/2 - (1/(4*w))*sin(2*w*tk))*v0
elif tk <= (TH - TA) :
s2 = (tk - TA)*v0
elif tk <= TH :
s3 = ((tk - TH + 2*TA)/2 - (1/(4*w))*sin(2*w*(tk - TH + 2*TA)))*v0 - s1
x[k]= s1 + s2 + s3
for k in range(1,n):
v[k] = (x[k+1]-x[k])/T0 # ableiten
for k in range(1,n):
a[k] = (v[k+1]-v[k])/T0
for k in range(1,n):
r[k] = (a[k+1]-a[k])/T0
plt.subplot(4,1,1)
plt.plot(tt,x)
plt.ylabel('Weg x'); plt.grid(True)
plt.subplot(4,1,2)
plt.plot(tt,v,'y')
plt.ylabel('Geschwindigkeit'); plt.grid(True)
plt.subplot(4,1,3)
plt.plot(tt,a,'g')
plt.ylabel('Beschleunigung'); plt.grid(True)
plt.subplot(4,1,4)
plt.plot(tt,r,'r')
plt.xlabel('Zeit in Sek')
plt.ylabel('Ruck'); plt.grid(True)
p= a*v #Leistung
plt.figure()
plt.plot(tt,p)
plt.xlabel('Zeit in Sek')
plt.ylabel('Leistung')
#from scipy.fftpack import fft
#ff = fft(a)
#z = 2.0/n*np.abs(ff[0:n//2])
#plt.plot(z)
