Einstieg in die Regelungstechnik mit Python

Fachbuch von Hans-Werner Philippsen

Kapitel 8

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
@author: philippsen
"""

import numpy as np
import control
# from   control.matlab import *
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(b=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
@author: philippsen
"""

import numpy as np 
import control 
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 = control.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 = control.forced_response(DCmo,[tk[i], tk[i]+T0], 
[u[i], u[i]], X0=Xvor) 
     xx[i] = xxk[1] 
     Xvor = xk[1] 
 
plt.plot(tk,xx, tk, u, tk,yyi) 
plt.xlabel('t [s]'); plt.ylabel('x, y, yi') 
plt.grid(b=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
@author: philippsen
"""

import numpy as np
import control
# from   control.matlab import *
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 = control.tf(dz,dn,t0)
t=np.arange(0,2,t0)
t1, y = control.step_response(diskSys,t)
plt.step(t1,y[0],where='post')
plt.xlabel('k*T0'); plt.ylabel('x(k)')
plt.grid(b=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
@author: philipp
"""

import numpy as np
import control
from control.matlab import *
import matplotlib.pyplot as plt

den  = np.array([3, 1])
pt1 = tf(1,den)
y, t = control.matlab.step(pt1)

plt.plot(t,y)

t0 = 1 # Abtastzeit 
diskPT1=control.c2d(pt1,t0,'tustin')

tt=np.arange(0,t.max(),t0)
t1, y = control.step_response(diskPT1,tt)
plt.step(t1,y[0],where='post')
plt.xlabel('t [s] und k'); plt.ylabel('h(t) und h(k)')
plt.grid(b=True)

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 Toolbox.

# -*- 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
@author: philippsen
"""

import numpy as np
import control
# from   control.matlab import *
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 = control.tf(dz,dn,t0)
#t=np.arange(0,4,t0)
#t1, y = control.step_response(dnotch,t)
#plt.step(t1,y[0],where='post')

notch2 = control.tf([b2, b1, b0],[a2, a1, a0])
dnotch2 = control.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(b=True)
plt.subplot(3,1,2) 
plt.plot(tt,v)
plt.ylabel('Geschwindigkeit'); plt.grid(b=True)
plt.subplot(3,1,3) 
plt.plot(tt,a)
plt.ylabel('Beschleunigung'); plt.grid(b=True)
plt.xlabel('Zeit in Sek'); plt.grid(b=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(b=True)
plt.subplot(4,1,2) 
plt.plot(tt,v,'y')
plt.ylabel('Geschwindigkeit'); plt.grid(b=True)
plt.subplot(4,1,3) 
plt.plot(tt,a,'g')
plt.ylabel('Beschleunigung'); plt.grid(b=True)
plt.subplot(4,1,4) 
plt.plot(tt,r,'r')
plt.xlabel('Zeit in Sek')
plt.ylabel('Ruck'); plt.grid(b=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)


Theme von Anders Norén