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, modifiziert für Version 0.10.1 am 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
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 = ct.tf(num,den, name='sys')
print(sys)
print("Pole = ", ct.poles(sys) )
print("Nullstellen = ", ct.zeros(sys) )
response = ct.pole_zero_map(sys)
ct.pole_zero_plot(response)
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, mod 15.4.24
@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.poles(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, mod 15.4.24, modifiziert für Version 0.10.1 am 19.8.24
@author: philippsen
"""
import numpy as np
import control as ct
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 = ct.matlab.tf(num,den)
regler = 1.0
G0 = regler*strecke
# Hinweis: mit print(G0.name) kann interne Systemnummer angezeigt werden
res = ct.nyquist_response(G0)
plt.plot(res.response.real,res.response.imag, "b")
res = ct.nyquist_response(2*G0)
plt.plot(res.response.real,res.response.imag, "r")
res = ct.nyquist_response(4*G0)
plt.plot(res.response.real,res.response.imag, "g")
plt.title('Ortskurven G0 Verstärkung 1, 2, 4')
plt.xlabel('Real'); plt.ylabel('Imag')
plt.grid()
plt.show()
plt.figure()
Gw = ct.matlab.feedback(G0, 1) # Kreis
xx, tt = ct.matlab.step(Gw)
plt.plot(tt,xx,"b")
Gw = ct.matlab.feedback(2*G0, 1) # Kreis
xx, tt = ct.matlab.step(Gw)
plt.plot(tt,xx,"r")
Gw = ct.matlab.feedback(4*G0, 1) # Kreis
xx, tt = ct.matlab.step(Gw)
plt.plot(tt,xx,"g")
plt.title('Sprungantwort Regelkreis für Kp=1, 2, 4')
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, mod. 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
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 = ct.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 = ct.tf([kr*tn*tv, kr*tn, kr],[tn, 0]) # ohne D-T1
regler = ct.tf([kr*(tv+t1)*tn, kr*(tn+t1), kr],[tn*t1, tn, 0])
regelkreis = ct.feedback(regler*strecke,1)
t, x = ct.step_response(regelkreis)
stell = ct.feedback(regler, strecke)
t, y = ct.step_response(stell,t)
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], '--',color="blue")
plt.plot([0, 36], [1.03, 1.03], '--',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, modifiziert für Version 0.10.1 am 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
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 = ct.tf(Kp,[T1, 1]) # P-T1
Gi = ct.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 = ct.tf([K*Tn, K],[Tn, 0])
G0=R*G
ct.bode_plot(G0,title='Bode-Plot G0',label=[''])
#bode(G0)
gm, pm, wg, wp = ct.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 = ct.feedback(G0,1)
tt, x = ct.step_response(Gw)
Gy = ct.feedback(R,G) # Stellgröße
tt, y = ct.step_response(Gy,tt)
fig, ax1 = plt.subplots()
ax1.plot(tt, 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 = ct.tf(1,[vt1*vt2, vt1+vt2, 1]) # Vorgabe
F = Gvg/Gw
print('Vorfilter = ')
print(F)
web606
I-T1-Strecke mit PI-Regler und Lead
Listing 6.6:
# -*- coding: utf-8 -*-
"""
Created on 4.1.22 IT1-Strecke, modifiziert für Vers 0.10.1 am 19.8.24
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.tf(K, [T1, 1.0, 0.0], inputs='y', outputs='xs')
Kr = 0.7; Tn = 0.6261
piregler = 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.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 13.6.23
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
pt = ct.tf(1,[1,1]) # Strecke
ptm = ct.tf(1,[1,1]) # Modell der Strecke
K=2; Tn = 1.1
pi = ct.tf([K*Tn,K],[Tn,0]) #PI-Regler
t0 = 0.1 # Abtastzeit
totzeit = 1 # Sekunden
# Totzeit als zeitdiskretes System
nn = int(totzeit/t0)
z=ct.tf('z')
tot=(z**(-nn) )/1
tot.input_labels = 'ys' ; tot.output_labels = 'xs'
totM = (z**(-nn) )/1
totM.input_labels = 'ym' ; totM.output_labels = 'xm'
# Diskretisierung Strecke Regler
ptdisk = ct.sample_system(pt,t0,method='zoh', name='ptd',inputs='y', outputs='ys')
ptdiskM = ct.sample_system(pt,t0,method='zoh', name='ptdM',inputs='y', outputs='ym') # Diskretisierung
pidisk = ct.sample_system(pi,t0,method='zoh', 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([tot, ptdisk,totM, ptdiskM,pidisk,
sumblk,sumblk2,sumblk3], inplist='r', outlist=['xs','y'])
t=np.arange(0,10,t0)
ww = np.ones(len(t))
t, yy = ct.input_output_response(Rk, t, ww)
plt.figure(1)
plt.plot(t,yy[0,:],t,yy[1,:])
plt.ylabel('y (rot) , x (blau)')
plt.xlabel('t [s]')
plt.title('Sprungantwort Regelkreis mit Smith Prädiktor')
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, mod. 20.8.24
@author: philippsen
"""
# Evo Opt Rglerentwurf
from scipy.optimize import differential_evolution
import numpy as np
import control as ct
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 = ct.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 = ct.tf([K*Tn, K],[Tn, 0.0])
G0 = G*R
Gw = ct.feedback(G0,1)
t,y = ct.step_response(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 = ct.tf([K*Tn, K],[Tn, 0.0])
G0 = G*R
Gw = ct.feedback(G0,1)
t, x = ct.step_response(Gw,T)
stell = ct.feedback(R, G)
t, y = ct.step_response(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 , mod. 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
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 = ct.tf(1.0, [0.2, 1.0])
strecke1 = ct.tf([1],[1, 1]) *pt1 # Strecke 1
strecke2 = ct.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 = ct.tf([Kp*Tn, Kp], [Tn, 0])
MKp[ii,jj] = Kp; MTn[ii,jj] = Tn
G0 = reg*strecke1
Gw = ct.feedback(G0, 1) # % Regelkreis
tt, xk = ct.step_response(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 = ct.feedback(reg, strecke1) # % Stellgröße
tt, yk = ct.step_response(Gy,tk)
g = np.dot(yk,yk)
Q1[ii,jj] = g # Stellenergie
G0 = reg*strecke2
Gw = ct.feedback(G0, 1) # % Regelkreis
tt, xk = ct.step_response(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 = ct.feedback(reg, strecke2) # % Stellgröße
tt, yk = ct.step_response(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 = ct.tf([K*Tn, K],[Tn, 0.0])
G0 = R*strecke1
Gw = ct.feedback(G0,1)
t, x = ct.step_response(Gw,tk)
stell = ct.feedback(R, strecke1)
t, y = ct.step_response(stell,tk)
G0 = R*strecke2
Gw = ct.feedback(G0,1)
t, x2 = ct.step_response(Gw,tk)
stell2 = ct.feedback(R, strecke2)
t, y2 = ct.step_response(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)
