Einstieg in die Regelungstechnik mit Python

Fachbuch von Hans-Werner Philippsen

Kapitel 6

Auf dieser Seite finden sich Python-Programme, die Berechnungen zur Stabilität durchführen und Entwurfsverfahren zeigen.


web601

Stabilität

Listing 6.1: Berechnet werden die Pole und Nullstellen eines Systems und ein Pol-Nullstellen-Plan wird geplottet.

# -*- coding: utf-8 -*-
"""
Stabilitätsanalyse: Berechnung der Pole bzw. Plot
der Pole und Nullstellen

Created on 17.4. 2019
@author: philippsen
"""

import numpy as np
#import control
import control.matlab as ma
import matplotlib.pyplot as plt

num = np.array([1,0.4, 1])
den  = np.array([1,2, 1, 4])
print(den)
print("Pole = ", str(np.roots(den)) )

sys = ma.tf(num,den)
print(sys)
print("Pole = ", ma.pole(sys) )
print("Nullstellen = ", ma.zero(sys) )
ma.pzmap(sys)
plt.grid()

web602

Wurzelortskurve

Listing 6.2: Für die Standardregelstrecke erfolgt die Berechnung der Wurzelortskurve.

# -*- coding: utf-8 -*-
"""
Standardregelstrecke G(s), Pole und rlocus

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

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

a0 = 1.; a1 = 10. ; a2 = 31. ; a3 = 30.
b0 = 1.
num  = np.array([b0])
den  = np.array([a3, a2, a1, a0])

stand = control.tf(num,den)
po=control.pole(stand)
print(stand,'Die Pole sind : ',po)
control.root_locus(stand)

plt.title('Wurzelortskurve')
plt.grid()

web603

Ortskurven offene Kreis

Listing 6.3: Es werden drei Ortskurven und Sprungantworten für unterschiedliche Verstärkungen berechnet.

# -*- coding: utf-8 -*-
"""
Ortskurven, h(t) Standardregelstrecke und PI-Regler

Created on 21.3. 2019
@author: philippsen
"""

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

a0 = 1.; a1 = 10. ; a2 = 31. ; a3 = 30.
b0 = 1.
num  = np.array([b0])
den  = np.array([a3, a2, a1, a0])

strecke = control.matlab.tf(num,den)
regler = 1.0
G0 = regler*strecke
real, imag, ww = control.matlab.nyquist(G0,Plot=False)
plt.plot(real, imag, "b")
real, imag, ww = control.matlab.nyquist(2*G0,Plot=False)
plt.plot(real, imag, "r")
real, imag, ww = control.matlab.nyquist(4*G0,Plot=False)
plt.plot(real, imag, "g")
plt.title('Ortskurven')
plt.grid()
plt.show()
plt.figure()
Gw = feedback(G0, 1) #   Kreis
xx, tt = control.matlab.step(Gw)
plt.plot(tt,xx,"b")
Gw = feedback(2*G0, 1) #   Kreis
xx, tt = control.matlab.step(Gw)
plt.plot(tt,xx,"r")
Gw = feedback(4*G0, 1) #   Kreis
xx, tt = control.matlab.step(Gw)
plt.plot(tt,xx,"g")
plt.title('Regelkreis')
plt.xlabel('t [s]'); plt.ylabel('x(t)')
plt.grid()
plt.show()
#plt.gcf().clear() # So wird Plot gelöscht

web604

Betragsoptimum

Listing 6.4: Gezeigt wird der PID-Entwurf gemäß Betragsotimum für die Standardregelstrecke.

# -*- coding: utf-8 -*-
"""
Standardregelstrecke  G(s) Betragsoptimum

Created on 27.1. 2020
@author: philippsen
"""

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

a0 = 1.
a1 = 10.
a2 = 31.
a3 = 30. ; a4 =0; a5 = 0
b0 = 1.
num  = np.array([b0])
den  = np.array([a3, a2, a1, a0])
strecke = tf(num,den)
D =np.linalg.det([[a1, -a0, 0],
 	    [a3, -a2, a1],
	    [a5, -a4, a3]])

r0=np.linalg.det([[a0**2,   -a0, 0],
 	 [(-a1**2 + 2*a0*a2),    -a2, a1],
	 [(a2**2  + 2*a0*a4 -2*a1*a3),-a4, a3]])
r0 = r0/D

r1=np.linalg.det([[a1,   a0**2,		0],
 	  [a3, 	(-a1**2 + 2*a0*a2),	 	a1],
	  [a5,	(a2**2  + 2*a0*a4 -2*a1*a3),	 a3]])
r1 = r1/D

r2=np.linalg.det([[a1,    -a0,	 a0**2],
 	  [a3, 	-a2,	(-a1**2 + 2*a0*a2) ],
	  [a5,	-a4,	(a2**2  + 2*a0*a4 -2*a1*a3)]	])
r2 = r2/D

kr = r1/2
tn = r1/r0
tv = r2/r1 ; t1 = tv/5.0
#regler = tf([kr*tn*tv, kr*tn, kr],[tn, 0]) # ohne D-T1
regler = tf([kr*(tv+t1)*tn, kr*(tn+t1), kr],[tn*t1, tn, 0])
regelkreis = feedback(regler*strecke,1)
x, t = step(regelkreis)

stell = feedback(regler, strecke)
y, t = step(stell)

fig, ax1 = plt.subplots()

ax1.plot(t, x, "b") 
ax1.set_ylabel('x', color="blue") 
plt.grid()
plt.plot([0, 36], [0.97, 0.97], 'k--',color="blue")
plt.plot([0, 36], [1.03, 1.03], 'k--',color="blue")
ax2 = ax1.twinx() 
ax2.plot(t, y, "r") 
ax2.set_ylabel('y', color="red") 
ax1.set_xlabel('t [s]', fontsize=12) 
plt.title('Entwurf gemäß Betragsoptimum, Tol. +-3%')

web605

I-T1 -Strecke mit PI-Regler

Listing 6.5: Symmetrisches Optimum

# -*- coding: utf-8 -*-
"""
Symmetrisches Optimum PI-Reg
P-T1 in Reihenschaltung mit I-Glied
Created on 17.8. 2021
@author: philippsen
"""

import numpy as np
import control
from   control.matlab import *
import matplotlib.pyplot as plt
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42

Kp= 1.0    #    Verstärkung Strecke
T1 = 0.2   #    Zeitkonstante
G1 = tf(Kp,[T1, 1])  # P-T1
Gi = tf(1.0,[1.0, 0])  # Integrator
G = G1*Gi          # I-T1
#K = 1.0; Tn = 0.4   # Regler
a = 4;K = 1.0/(a*Kp*T1); Tn = T1*a**2   # Regler gemäß Sym. Optimum
R = tf([K*Tn, K],[Tn, 0])
G0=R*G
bode(G0)
gm, pm, wg, wp = control.matlab.margin(G0)
gmdb = 20*np.log10(gm); wphz = wp/(2*np.pi)
print('Reglerverstärkung = ',K,' Tn = ',Tn)
print('Amplitudenreserve = ',gm,' in dB = ',gmdb)
print('Phasenreserve = ',pm)
print('Durchtrittsfrequenz = ',wp,' in Hertz = ',wphz)

Gw = control.feedback(G0,1)
x, t = control.matlab.step(Gw)

Gy = control.feedback(R,G) # Stellgröße
y, tt = control.matlab.step(Gy)

fig, ax1 = plt.subplots()

ax1.plot(t, x, "b") 
ax1.set_ylabel('x', fontsize=14, color="blue") 
plt.grid()
ax2 = ax1.twinx() 
ax2.plot(tt, y, "r") 
ax2.set_ylabel('y', fontsize=14,color="red") 
ax1.set_xlabel('t [s]', fontsize=12) 
plt.title('Entwurf gemäß symm. Optimum')    
#Entwurf Vorfilter
vt1=0.5 ; vt2=0.5
Gvg = tf(1,[vt1*vt2, vt1+vt2, 1]) # Vorgabe
F = Gvg/Gw
print(F)

web606

I-T1-Strecke mit PI-Regler und Lead

Listing 6.6:

# -*- coding: utf-8 -*-
"""
Created on 4.1.22     IT1-Strecke
PI-Regler und PDT1 als extra Rückführung
@author: philippsen
"""
import control as ct
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42

K=1; T1 = .2
strecke = ct.tf2io(ct.tf(K, [T1, 1.0, 0.0]), inputs='y', outputs='xs')
Kr = 0.7; Tn = 0.6261
piregler = ct.tf2io(ct.tf([Kr*Tn, Kr], [Tn, 0.0]), inputs='e', outputs='ypi')
sumblk = ct.summing_junction(inputs=['w', '-xs'], output='e')
r1 =  2*np.sqrt(T1/K)-1/K
Tf = abs(r1)/10
pdt1 = ct.tf2io(ct.tf([Tf+r1, 1.0], [Tf, 1.0]), inputs='xs', outputs='yd')
sumblk2 = ct.summing_junction(inputs=['ypi', '-yd'], output='y')

#Regelkreis aufbauen
Rk = ct.interconnect([strecke, pdt1,
    piregler, sumblk, sumblk2], inplist='w', outlist=['xs','y'])

T = np.arange(0.0, 10.0, 0.01) # Simulationsdauer und T0 

ww = np.ones([len(T)]) # Eingangsgröße

t, yy = ct.input_output_response(Rk, T, ww)

# Plot Regelkreis-Sprungantwort
plt.figure(1)
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel('y (rot) , x (blau)')
plt.title('Sprungantwort I-T1 mit PI und Lead')
plt.grid()


web607

Totzeit Strecke mit Smith-Prädiktor

Listing 6.7:

# -*- coding: utf-8 -*-
"""
Smith-Prädiktor
Created on 6.1.22
@author: philippsen
"""

import numpy as np
import control as ct
# from   control.matlab import *
import matplotlib.pyplot as plt
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42
  
pt = ct.tf(1,[1,1])
#pt = ct.tf(1,[30,31,10,1])
ptm = ct.tf(1,[1,1])
K=2; Tn = 1.1
pi = ct.tf([K*Tn,K],[Tn,0])     #PI-Regler
t0 = 0.1  # Abtastzeit
totzeit = 1
nn = int(totzeit/t0)

z=ct.tf('z')
tot=(z**(-nn) )/1
totM = tot
print(tot)
ptdisk = ct.sample_system(pt,t0,method='zoh')	# Diskretisierung
ptdiskM = ct.sample_system(pt,t0,method='zoh')	# Diskretisierung
pidisk = ct.sample_system(pi,t0,method='zoh')	# Diskretisierung
#Strecke in ioSys aufbauen
iotot = ct.tf2io(tot, inputs='ys', outputs='xs')
iopt = ct.tf2io(ptdisk, inputs='y', outputs='ys')
#Streckenmodell in ioSys aufbauen
iototM = ct.tf2io(totM, inputs='ym', outputs='xm')
ioptM = ct.tf2io(ptdiskM, inputs='y', outputs='ym')

iopi = ct.tf2io(pidisk, inputs='e', outputs='y')

sumblk = ct.summing_junction(inputs=['r', '-xx'], output='e')
sumblk2 = ct.summing_junction(inputs=['xs', '-xm'], output='em')
sumblk3 = ct.summing_junction(inputs=['ym', 'em'], output='xx')

Rk = ct.interconnect([iotot, iopt,iototM, ioptM,iopi,
                      sumblk,sumblk2,sumblk3], inplist='r', outlist=['xs','y'])
t=np.arange(0,20,t0)
ww = np.ones(len(t))
t, yy = ct.input_output_response(Rk, t, ww)

t1, y = ct.step_response(tot*ptdisk,t)
#plt.step(t1,y)  #,where='post')
#plt.plot(t1,y+0.1)
plt.plot(t,yy[0,:],t,yy[1,:])
plt.grid(True)

web608

Totzeit Strecke mit Smith-Prädiktor in Julia

Listing 6.8: In der Julia-Konsole zunächst die ControlSystems Toolbox installieren mit

using Pkg; Pkg.add("ControlSystems")

und gegebenenfalls noch weitere Pakete wie z.B. „Plots“.

# Simulation Smith-Prädiktor-Regelung für Strecke mit Totzeit
# Verwendung: "ControlSystems Toolbox for Julia" 
# Philippsen, Einstieg in die RT 4. Auflage
# siehe dort auch [Bag21]
#

using ControlSystems
using OrdinaryDiffEq, DelayDiffEq
using Plots

## Zwei P-T1 mit Totzeit
g1=tf(1.0, [1, 1])
g2=tf(0.8, [2.0, 1])
g=g1*g2
gt = delay(5)
ggt = g*gt
#print(ggt)

# PI - Regler
K=1.0; Tn = 2
r=tf([K*Tn, K],[Tn,0])

# Formel Smith-Prädiktor und Aufbau Regelkreis
prae = feedback(r, g*(1 - gt))
rekreis = feedback(prae*ggt, 1)

# Simulation
t = 0:0.1:30
y, = step(rekreis, t)
plot(t,y',w=2,grid=true)

web609

Numerische Opt. mit evolutionärem Algorithmus

Listing 6.9:

# -*- coding: utf-8 -*-
"""
Created on Sun Sep 26 13:45:39 2021

@author: philippsen
"""

# Evo Opt Rglerentwurf
from scipy.optimize import differential_evolution
import numpy as np
from   control.matlab import *
import matplotlib.pyplot as plt
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42

Kp= 1.0    #    Verstärkung Strecke
T1 = 1.0   #    Zeitkonstante
G = tf(Kp,[30,31,10, 1])  # Standardregelstrecke
T = np.arange(0.0, 40.0, 0.05) # Simulationsdauer
 
# Zielfunktion
def objective(v):
    K, Tn = v
    R = tf([K*Tn, K],[Tn, 0.0])    
    G0 = G*R
    Gw = feedback(G0,1)
    y, t = step(Gw,T) # Sprungantwort
    return np.sum(t*np.abs(1-y))*0.05 # ITAE

# Parameterbereich
p_min, p_max = 0.0, 10.0
bounds = [[p_min, p_max], [p_min, p_max]]
# Durchführung der evolutionären Optimierung
result = differential_evolution(objective, bounds)
# Ergebnisse
print('Status : %s' % result['message'])
print('Anzahl Iterationen: %d' % result['nfev'])
# Überprüfung der Optimierung
solution = result['x']
evaluation = objective(solution)
print('Lösung: f(%s) = %.5f' % (solution, evaluation))
K, Tn = solution
R = tf([K*Tn, K],[Tn, 0.0])    
G0 = G*R
Gw = feedback(G0,1)
x, t = step(Gw,T)
stell = feedback(R, G)
y, t = 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) 
s = str("Entwurf gemäß Evo Opt Kr = %1.2f "% K)
s2 = str("Tn = %1.2f" % Tn)
plt.title(s+s2)

web610

Numerische Optimierung im Batch-Betrieb

Listing 6.10

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

Regelstrecke mit einer wechselnden Zeitkonstanten und PI-Regler, die 
Regler-Parameter werden in einem Wertebereich variert (Batch-Betrieb) 
und die Güte in einem 3D-Plot dargestellt. 
Die optimalen Parameter sind im Minimum zu finden

Created on 2.12.  2021 
@author: philippsen
""" 

import numpy as np
import control
import control.matlab
import matplotlib.pyplot as plt
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

n = 400 # Anzahl Schritte
T0 = 0.02 # 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 = 3/kk # Änderungsstufen Parameter
MTn = np.zeros((kk,kk))
Q1 = np.zeros((kk,kk)); Q2 = np.zeros((kk,kk)) # Güte-Matrizen
Qs = np.zeros((kk,kk)) # Summe Güte-Matrizen
#Erzeugung untere und ober Grenzlinie
nn = len(tk); glT1 = 2.4; glT2 = 3.1; tolband = 0.03
gl = np.zeros(nn) # Grenzlinie untere
# i Werte für glT1 und glT2
w = 1.0    #Sollwert
glT1i = int(glT1/T0); glT2i = int(glT2/T0)
sh = (w - tolband)/(glT2i - glT1i) #Stufenhöhe

for i in range(glT1i, nn): # Einträge Werte in gl
    gl[i] = (i - glT1i)*sh
    if i > glT2i:
        gl[i] = w - tolband


pt1 = control.tf(1.0, [0.2, 1.0])
strecke1 = control.tf([1],[1, 1]) *pt1  # Strecke 1
strecke2 = control.tf([1],[2, 1]) *pt1  # Strecke 2

for ii in range(0, kk): # Liste von 0 bis kk-1
    for jj in range(0, kk):
        Kp = 1.0 + dKp_kk*(ii + 1); Tn = 0.0 + dTn_kk*(jj +1)      # Reglerparameter
        reg = control.tf([Kp*Tn, Kp], [Tn, 0])
        MKp[ii,jj] = Kp; MTn[ii,jj] = Tn
        G0 = reg*strecke1
        Gw = control.feedback(G0, 1) # % Regelkreis
        xk , tt = control.matlab.step(Gw,tk)
        drin = True
        for i in range(glT1i, nn): # Überschwingen nicht zulassen ab 0, sonst ab T2
            if xk[i] > (w+tolband) or xk[i] < gl[i]:
                drin = False
                g = 1000
                # Die untere Grenzlinie könnte durch w-tolband ersetzt werden!
        
        if drin:
            Gy = control.feedback(reg, strecke1) # % Stellgröße
            yk , tt = control.matlab.step(Gy,tk)
            g = np.dot(yk,yk)
        Q1[ii,jj] = g # Stellenergie
        
        G0 = reg*strecke2
        Gw = control.feedback(G0, 1) # % Regelkreis
        xk , tt = control.matlab.step(Gw,tk)
        drin = True
        for i in range(glT1i, nn): # Überschwingen nicht zulassen ab 0, sonst ab T2
            if xk[i] > (w+tolband) or xk[i] < gl[i]:
                drin = False
                g = 1000
                # Die untere Grenzlinie könnte durch w-tolband ersetzt werden!
        
        if drin:
            Gy = control.feedback(reg, strecke2) # % Stellgröße
            yk , tt = control.matlab.step(Gy,tk)
            g = np.dot(yk,yk)*T0
        Q2[ii,jj] = g # Stellenergie


print('Minimum1 = ',Q1.min(),'Anzahl Iterationen = ',kk*kk )
pa = np.where(Q1 == Q1.min())
print('Kp = ',MKp[pa[0],pa[1]],'   Tn = ',MTn[pa[0],pa[1]] )
print('Minimum2 = ',Q2.min(),'Anzahl Iterationen = ',kk*kk )
pa = np.where(Q2 == Q2.min())
print('Kp = ',MKp[pa[0],pa[1]],'   Tn = ',MTn[pa[0],pa[1]] )
Qs = Q1 + Q2
print('MinimumSumme = ',Qs.min(),'Anzahl Iterationen = ',kk*kk )
pa = np.where(Qs == Qs.min())
print('Kp = ',MKp[pa[0],pa[1]],'   Tn = ',MTn[pa[0],pa[1]] )
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(MKp, MTn, Q1, rstride=1, cstride=1, cmap=cm.viridis)
plt.xlabel('Kp'); plt.ylabel('Tn'); ax.set_zlabel('Q1')
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(MKp, MTn, Q2, rstride=1, cstride=1, cmap=cm.viridis)
plt.xlabel('Kp'); plt.ylabel('Tn'); ax.set_zlabel('Q2')

K = MKp[pa[0],pa[1]][0]; Tn = MTn[pa[0],pa[1]][0]
R = control.matlab.tf([K*Tn, K],[Tn, 0.0])    
G0 = R*strecke1
Gw = control.matlab.feedback(G0,1)
x, t = control.matlab.step(Gw,tk)
stell = control.matlab.feedback(R, strecke1)
y, t = control.matlab.step(stell,tk)

G0 = R*strecke2
Gw = control.matlab.feedback(G0,1)
x2, t = control.matlab.step(Gw,tk)
stell2 = control.matlab.feedback(R, strecke2)
y2, t = control.matlab.step(stell2,tk)

fig, ax1 = plt.subplots()
ax1.plot(t, x,'b', t,x2,"b--") 
simEnde = t.max()
ax1.plot([0, simEnde], [w+tolband, w+tolband], 'g--') 
ax1.plot([glT2, simEnde], [w-tolband, w-tolband], 'g--') 

ax1.set_ylabel('x x2--', color="blue", fontsize=14) 
plt.grid()
ax2 = ax1.twinx() 
ax2.plot(t, y,'r', t,y2,"r--") 
ax2.set_ylabel('y y2--', color="red", fontsize=14) 
ax1.set_xlabel('t [s]', fontsize=12) 
s = str("Robuster Entwurf vollst. Parameterraum Kr = %1.2f "% K)
s2 = str("Tn = %1.2f" % Tn)
plt.title(s+s2)

Theme von Anders Norén