In den Beispielprogrammen des vierten Kapitels werden lineare und nicht lineare Regelstrecken simuliert, Kennwerte ermittelt und Frequenzgänge berechnet.
web401
P-Tn – Strecke
Listing 4.1: Reihenschaltung von drei P-T1-Gliedern, Simulation der Sprungantwort der Teilstrecken und Gesamtstrecke sowie Plot der Ortskurve.
# -*- coding: utf-8 -*-
"""
P-Tn Reihenschaltung 3 P-T1
Created on 15.4. 2019, mod 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
Kp= 3.0 # Verstärkung
T1 = 2.0 # Zeitkonstanten
T2 = 3.0
T3 = 5.0
N1 = np.array([T1, 1])
N2 = np.array([T2, 1])
N3 = np.array([T3, 1])
G1 = ct.tf(Kp,N1)
G2 = ct.tf(1.0,N2)
G3 = ct.tf(1.0,N3)
Gges = G1*G2*G3
t, x3 = ct.step_response(Gges)
t, x2 = ct.step_response(G1*G2,t)
t, x1 = ct.step_response(G1,t)
plt.plot(t, x1, t, x2,t, x3, "b")
plt.ylabel('x1, x2, x3', color="blue", fontsize=14)
plt.xlabel('t [s]', fontsize=12)
plt.title('Reihenschaltung drei P-T1', fontsize=12)
plt.grid()
plt.figure()
response = ct.nyquist_response(Gges)
plt.plot(response.response.real,response.response.imag)
plt.title('Ortskurve P-T3, pos. Frequenzen')
plt.xlabel('Real'); plt.ylabel('Imag')
plt.grid()
web402
P-T2s – Strecke
Listing 4.2: Eine P-T2s Strecke wird simuliert, wobei die Dämpfung variert wird. Darüber hinaus erfolgt die Berechnung des Plots des Pol-Nullstellen-Plans.
# -*- coding: utf-8 -*-
"""
P-T2s Modifikation Dämpfung Batchbetrieb
Created on 2.4. 2019, 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
k=11
wurz = np.zeros((k,2),dtype=np.complex_)
for ii in range(0, k): # Liste von 0 bis kk-1
d = 1.2 - 0.1*ii
den = np.array([1,2*d, 1])
print("d = %3.2f, Pole = %s " %(d,str(np.roots(den))) )
wurz[ii,:] = np.roots(den)
pt2 = ct.tf(1,den)
t, y = ct.step_response(pt2)
plt.plot(t,y)
plt.title('P-T2s Batch')
plt.xlabel('t [s]'); plt.ylabel('h(t)')
plt.grid()
plt.figure()
response = ct.pole_zero_map(pt2)
ct.pole_zero_plot(response)
plt.grid()
web403
Totzeit mit P-T1 mit Pade – Approximation
Listing 4.3: Die Simulation und Berechnung des Frequenzgangs einer Strecke mit Totzeit kann mit Hilfe der Pade-Approximation angenähert werden.
# -*- coding: utf-8 -*-
"""
Totzeit mit P-T1 Pade-Approximation
Created on 10.4. 2019, mod 15.4.24 und 20.8.24
Neues Plot-Paradigma ab Version 0.10.0 für bode und nyquist
@author: philippsen
"""
import control as ct
import numpy as np
import matplotlib.pyplot as plt
T = 1 # Zeitkonstante P-T1
TT = 0.5 # Totzeit
n = 500 # Anzahl Schritte
T0 = 0.01 # Abtastzeit
u = np.zeros(n) # Vektor Eingangsgröße
tk= np.arange(0,T0*n,T0) # Zeit- bzw. k-Achse
ab = round(TT/T0)
for i in range(ab, n): # Verschiebung um Totzeit
u[i] = 1.0
G1 = ct.tf([ 1],[T, 1]) # P-T1
tk, yout = ct.forced_response(G1, tk, u, X0=0)
num, den = ct.pade(TT,3) # Pade Appr. 3.Ordnung
Gpade = ct.tf(num,den)
Gers = G1*Gpade
tk, x = ct.step_response(Gers, tk)#, X0=0)
plt.plot(tk,yout, tk, x)
plt.title('Totzeit-P-T1 und Pade-Approx.')
plt.xlabel('t [s]'); plt.ylabel('x')
plt.grid(True)
plt.figure()
ct.bode_plot(Gers,title='Bode-Plot G Ersatz',label=[''])
web404
Integrale Strecke mit Verzögerung (I-T1)
Listing 4.4: Eine I-T1-Strecke wird simuliert, das Bode-Diagramm berechnet und die Orstkurve mit Markierung der Frequenz geplottet.
# -*- coding: utf-8 -*-
"""
P-T1 in Reihenschaltung mit I-Glied
Created on 15.4. 2019 and 7.4.2024 and 20.8.24
Neues Plot-Paradigma ab Version 0.10 für bode und nyquist
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
Kp= 2.0 # Verstärkung
T1 = 0.2 # Zeitkonstante
N1 = np.array([T1, 1])
G1 = ct.tf(Kp,N1) # P-T1
Gi = ct.tf(1.0,[1.0, 0]) # Integrator
G4 = G1*Gi # I-T1
time = np.arange(0.0, 1.0, 0.001)
tt, x4 = ct.step_response(G4,time)
plt.figure()
plt.plot(tt, x4, "r")
plt.ylabel('x4', color="red", fontsize=14)
plt.xlabel('t [s]', fontsize=12)
plt.title('I-T1 Sprungantwort', fontsize=12)
plt.grid()
plt.figure()
ct.bode_plot(G4)
plt.figure() # Ortskurve mit neuer Funktion
response = ct.nyquist_response(G4)
response.plot()
plt.title('Control Nyquist_response')
plt.figure() # Nur pos. Frequenzen
response = ct.nyquist_response(G4,omega_limits=[0.1 ,100])
plt.plot(response.response.real,response.response.imag)
plt.title('Ortskurve I-T1 mit limits, pos. Frequenzen')
plt.xlabel('Real'); plt.ylabel('Imag')
plt.grid()
# Eine andere Lösung mit dem SciPy - Paket
from scipy import signal
num = G4.num[0][0][:] # Zugriff auf Zähler
den = G4.den[0][0][:] # Zugriff auf Nenner
gg=signal.TransferFunction(num,den)
w, H = signal.freqresp(gg)
plt.figure()
plt.plot(H.real, H.imag, "b")
plt.title('Ortskurve I-T1 scipy')
#plt.plot(H.real, -H.imag, "r")
plt.grid()
web405
Lead-Lag Bode-Diagramm
Listing 4.5: Mit wenigen Programmierzeilen kann das Bode-Diagramm eines Lag-Glieds geplottet werden. Verändern Sie a auf 10, dann erfolgt die Berechnung eines Leads.
# -*- coding: utf-8 -*- """ Lead Lag Created on 10.4. 2019, mod. 15.4.24 Neues Plot-Paradigma ab Version 0.10.0 für bode und nyquist @author: philippsen """ import control import numpy as np import matplotlib.pyplot as plt T = 1 a = 0.10 GLL = control.tf([a*T, 1],[T, 1]) # Lead Lag #w, mag, phase = control.bode(GLL) out=control.bode_plot(GLL) w, mag, phase = control.frequency_response(GLL)
web406
Fahrzeug Querdynamik/Lenkmodell
Listing 4.6: Das nicht lineare Modell der Fahrzeuglenkung wird mit dem SciPy DGL-Löser berechnet, wobei die Fahrt in der Ebene geplottet wird.
# -*- coding: utf-8 -*-
"""
Fahrzeug Querdynamik Lenkmodell
Created on 27.1. 2020
@author: philippsen
nach 10 Sek wird Lenkwinkel auf null gestellt
"""
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
def deriv(xa, t, alpha):
beta=xa[0]
#x = xa[1]; z = xa[2]
v = 0.79
Rs = 0.5
bp = np.sin(alpha)*v/Rs
xp = v*np.sin(alpha + beta)
zp = v*np.cos(alpha + beta)
xapunkt = [bp, xp, zp]
return xapunkt
time = np.arange(0.0, 10.0, 0.01)
xaAnfang = [0,0,0]
u = 0.2
xa = odeint(deriv, xaAnfang, time, args=(u,))
xb = odeint(deriv, xa[999,:], time, args=(0,))
plt.figure()
plt.plot(xa[:,2],xa[:,1] ) #
plt.plot(xb[:,2],xb[:,1] )
plt.axis('scaled')
plt.ylabel('x', fontsize=14)
plt.xlabel('z', fontsize=14)
plt.title('Lösung Lenkung-DGL mit Python')
plt.grid()
web407
DC-Motor
Listing 4.7: Berechnet werden Strom und Drehzahl im Fall einer sprungförmigen Anpassung der Ankerspannung für den Faulhaber-Motor 2237.
# -*- coding: utf-8 -*-
"""
DC-Motor G(s) Faulhaber Motor 2237
Created on 2.4. 2019, mod 20.8.24
@author: philippsen
"""
import numpy as np; from math import pi
import control as ct
import matplotlib.pyplot as plt
cPsi= 0.0318 # Datenblatt Faulhaber 2237
R = 15.7 # Widerstand Anker
L = 590e-6 # Induktivität Anker
zpiJ = 47.119e-6 # Schwungmasse Berechnung 2pi7,5e-6
#zpiJ = 2*pi*75e-7
Kr = 35.82e-6
zN = np.array([cPsi/R ])
nN = np.array([ zpiJ, Kr+(2*pi*cPsi**2)/R])
zI = np.array([zpiJ/R , Kr/R])
DCmot = ct.tf(zN,nN)
t, n = ct.step_response(DCmot)
DCmotI = ct.tf(zI,nN)
tt, ii = ct.step_response(DCmotI)
fig, ax1 = plt.subplots()
ax1.plot(tt, n*60*5, "b")
ax1.set_ylabel('n [1/min]', color="blue", fontsize=14)
plt.grid()
ax2 = ax1.twinx()
ax2.plot(tt, ii*5*1000, "r")
ax2.set_ylabel('i [mA]', color="red", fontsize=14)
ax1.set_xlabel('t [s]', fontsize=12)
ax1.set_title('DC-Motor 2237', fontsize=12)
web408
RLC-Ersatzschaltung Tiefsetzsteller
Listing 4.8: Esatzmodell für den Tiefsetzsteller
# -*- coding: utf-8 -*-
"""
Created on 7.10.24 Tiefsetz RLC
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
# RLC - Schaltung Tiefsetzersatz
R1= 0.8; R2 = 20; L = 400e-6; C = 1000e-6
G = ct.tf(1,[L*C, R1*C + L/R2, 1 + R1/R2])
po = ct.poles(G)
print(po)
print('d = ',ct.damp(G,doprint=False)[1][0])
print(ct.step_info(G))
t, yy = ct.step_response(G)
plt.plot(t,yy)
plt.title(' Sprungantwort RLC Tiefsetz')
plt.xlabel('t [s]'); plt.ylabel('u2')
plt.grid()
k = 40
re = np.zeros(k)
im = np.zeros(k)
daem = np.zeros(k)
R2i = np.zeros(k)
for i in range(k):
R2 = R2i[i] = i+1
G = ct.tf(1,[L*C, R1*C + L/R2, 1 + R1/R2])
daem[i] = ct.damp(G,doprint=False)[1][0]
# po = G.damp()[2][0]
po = ct.poles(G)
re[i] = po[0].real
im[i] = po[0].imag
plt.figure()
plt.plot(re,im,re,-im)
plt.title('Eigenwerte RLC Tiefsetz R2 variabel 1-40')
plt.xlabel('Realteil'); plt.ylabel('Imag')
plt.grid()
plt.figure()
plt.plot(R2i,daem)
plt.title('Dämpfung RLC Tiefsetz')
plt.xlabel('R2 [Ohm]'); plt.ylabel('d')
plt.grid()
web409
Identifikation RC3
Listing 4.9:
# -*- coding: utf-8 -*-
"""
RC3 Identifikation mit leastsquare
Created on 25.11.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
from scipy.optimize import leastsq, least_squares
# Einlesen Messdaten
a = np.genfromtxt('d:\\phi\\Buch5teAuflage\\micropython\\RC3picolog4.txt',skip_header=0)
# Aufbereitung Daten für Ident. und Plot
t = np.arange(0,25,0.05)
N= len(t)
x= a[:,2]
# Definition Funktion für leastsquare mit Vorgabe Typ und Ordnung
def fungs(p,y,t): # U0 = 1; für Sprungantwort
[b0, a2, a1, a0]=p
g = ct.tf(b0,[1,a2,a1,a0])
t,yo = ct.step_response(g,t)
err = y - yo
return err
p0 =[1,2,1,1] # Startparameter setzen und Parameter schätzen
parls = leastsq(fungs, p0, args=(x,t))
b0 = parls[0][0]
a2 = parls[0][1]
a1 = parls[0][2]
a0 = parls[0][3]
gg = ct.tf(b0,[1,a2,a1,a0])
print('G(s) = ',gg); print('Verstärkung = ',gg.dcgain())
print('Tsumme = ',sum(-1/gg.poles()))
t,xx = ct.step_response(gg,t)
plt.plot(t,x,t,xx)
plt.title('RC3 Messung und Simulation ERsatzmodell')
plt.xlabel('t [s]'); plt.ylabel('ua(t)')
plt.grid()
plt.figure()
# Cosinusantwort
T0 = 0.05
T = np.arange(0.0, 50.0, T0) # Simulationsdauer
ww = np.zeros([len(T)])
for i in range(len(T)):
ww[i] = 1.0- np.cos(i*T0)
tt, yy = ct.forced_response(gg, T, ww)
b = np.genfromtxt('d:\\phi\\Buch5teAuflage\\micropython\\RC3picosinus.txt',skip_header=0)
xcos= b[:,2]
plt.plot(tt,ww,tt,yy*1.32,tt,xcos)
plt.title('RC3 Messung und Simulation Cosinus')
plt.xlabel('t [s]'); plt.ylabel('y, x')
plt.grid('True')
#plt.figure()
fig, ax1 = plt.subplots()
ax1.plot(tt, yy*1.32,tt,xcos)
ax1.set_ylabel('xsim, xmess')
plt.grid()
ax2 = ax1.twinx()
ax2.plot(tt, ww)
ax2.set_ylabel('y', color="blue")
ax1.set_title('RC3 Messung und Simulation Sinusantwort')
ax1.set_xlabel('t [s]', fontsize=15)
