Nachfolgend ist ein Python-Programm gelistet, das ein Mehrgrößensystem in Form einer Übertragungsmatrix darstellt und eine statische Entkopplung berechnet. Darüber hinaus finden Berechnungen in der Zustandsraumdarstellung statt, z.B. Transformation auf RNF.
web701
Mehrgrößensystem G(s)
Listing 7.1: Es erfolgt die Erzeugung eines Mehrgrößensystems mit der Python Control Library und Berechnung der statischen Entkopplungsmatrix. Für die Berechnung von Mehrgrößensystem ist das Paket slycot erforderlich.
# -*- coding: utf-8 -*-
"""
Erzeugung eines Mehrgrößensystems mit Sympy und der
Python Control Systems Library für die Berechnung der statischen
sowie dyn. Entkopplung, mt PI-Regelkreis, slycot erforderlich
Entkopplung funktioniert nicht im Fall integralwirkende Strecke!
Created on 10.1.25
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
n11= [2.0]; d11= [4., 9., 6., 1.]
n22= [2.0]; d22= [2., 5., 4., 1.]
n12= [-1.5]; d12= [0.4, 1.4, 1.]
n21= [1.5]; d21= [0.8, 2.4, 1.]
g11 = ct.tf(n11,d11); g110 = ct.dcgain(g11)
g12 = ct.tf(n12,d12); g120 = ct.dcgain(g12)
g22 = ct.tf(n22,d22); g220 = ct.dcgain(g22)
g21 = ct.tf(n21,d21); g210 = ct.dcgain(g21)
num = [[n11, n12], [n21, n22]]
den = [[d11, d12], [d21, d22]]
strecke = ct.tf(num, den,name='streck',
inputs = ['yf1','yf2'], outputs = ['x1','x2'])
num = [[n11, [0]], [[0], n22]]
den = [[d11, [1]], [[1], d22]]
streckeohne = ct.tf(num, den,name='streckohnekopp',
inputs = ['yf1','yf2'], outputs = ['x1','x2'])
# ohne Entkopplung, DiagonalPGlied
ne = [[[1], [0]], [[0], [1]]]
de = [[[1.], [1.]], [[1.], [1.]]]
noent = ct.tf(ne, de,name='noent',
inputs = ['y1','y2'], outputs = ['yf1','yf2'])
# stat Entkopplung
k= g110*g220/(g110*g220 - g210*g120)
ne = [[[k], [-k*g120/g110]], [[-k*g210/g220], [k]]]
de = [[[1.], [1.]], [[1.], [1.]]]
ent = ct.tf(ne, de,name='ent',
inputs = ['y1','y2'], outputs = ['yf1','yf2'])
# dyn Entkopplung lead
k= g110*g220/(g110*g220 - g210*g120)
leadnum = [-k*g120/g110, -k*g120/g110]; leadnum2 = [-2*k*g210/g220, -k*g210/g220]
leadden = [0.2, 1]; leadden2 = [0.4, 1]
ne = [[[k], leadnum], [leadnum2, [k]]]
de = [[[1.], leadden], [leadden2, [1.]]]
entlead = ct.tf(ne, de,name='entl',
inputs = ['y1','y2'], outputs = ['yf1','yf2'])
# PI-Regler
kr = 0.9; tn = 5
krtn=kr*tn
reg1 = ct.tf([krtn, kr],[tn,0],name='reg1',
inputs = ['e1'], outputs = ['y1'])
kr = 0.8; tn = 3
krtn=kr*tn
reg2 = ct.tf([krtn, kr],[tn,0],name='reg2',
inputs = ['e2'], outputs = ['y2'])
sum1 = ct.summing_junction(inputs=['w1', '-x1'], output='e1')
sum2 = ct.summing_junction(inputs=['w2', '-x2'], output='e2')
#Regelkreise aufbauen
RkohneKopp = ct.interconnect([streckeohne, noent, reg1, reg2,
sum1, sum2], inplist=['w1', 'w2'], outlist=['x1','x2'])
Rk0 = ct.interconnect([strecke, noent, reg1, reg2,
sum1, sum2], inplist=['w1', 'w2'], outlist=['x1','x2'])
Rk = ct.interconnect([strecke, ent, reg1, reg2,
sum1, sum2], inplist=['w1', 'w2'], outlist=['x1','x2'])
Rk2 = ct.interconnect([strecke, entlead, reg1, reg2,
sum1, sum2], inplist=['w1', 'w2'], outlist=['x1','x2'])
# Simulationen
T0 = 0.05
T = np.arange(0.0, 40.0, T0) # Simulationsdauer
ww = np.ones([2,len(T)]) # Eingangsgröße
ww[0,:] = ww[0,:]*4
ww[1,:] = ww[1,:]*3
t, yy = ct.input_output_response(RkohneKopp, T, ww)
# Plot Regelkreis-Sprungantwort ohne Verkopplung und Entkopplung
plt.figure()
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel(' x1, x2')
plt.title('Regelkreise ohne Verkopplung in der Strecke')
plt.grid()
t, yy = ct.input_output_response(Rk0, T, ww)
# Plot Regelkreis-Sprungantwort ohne Entkopplung
plt.figure()
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel('x1, x2')
plt.title('Regelkreise mit Verkopplung in der Strecke')
plt.grid()
t, yy = ct.input_output_response(Rk, T, ww)
# Plot Regelkreis-Sprungantwort
plt.figure()
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel('x1, x2')
plt.title('Regelkreise mit stat. Entkopplung')
plt.grid()
t, yy = ct.input_output_response(Rk2, T, ww)
# Plot Regelkreis-Sprungantwort
plt.figure()
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel('x1, x2')
plt.title('Regelkreise mit dyn. Lead-Entkopplung')
plt.grid()
# SymPy vollständige dynamische Entkopplung
from sympy.abc import s
from sympy import pprint
import sympy as sy
import sympy.physics.control as co
#gtt = sy.exp(-s)
gg11 = 2 / ( 4*s**3 + 9*s**2 + 6*s + 1)
gg22 = 2 / ( 2*s**3 + 5*s**2 + 4*s + 1)
gg12 = -1.5 / (0.4*s**2 + 1.4*s + 1)
gg21 = 1.5 / ( 0.8*s**2 + 2.4*s + 1)
gg = sy.Matrix([
[gg11, gg12],
[gg21, gg22]])
gdiag = sy.Matrix([
[gg11, 0],
[0, gg22]])
ff = sy.simplify(gg.inv()*gdiag,rational=True )
print('Entkopplungsmatrix = ')
pprint(sy.simplify(ff,rational=True ))
# Probe mit sy.simplify(gg*sy.simplify(ff),rational=True)
# Sympy Matrix in control Mehrgrößensystem wandeln
def zzz(fff,ii,jj):
print(ii,jj)
fe=sy.simplify(ff[ii,jj],rational=True)
nenn = sy.Poly.from_expr(fe.as_numer_denom()[1]) # Nenner holen
#sy.nroots(nenn) # wenn Pole gewünscht
nenk=sy.Poly.coeffs(nenn) # Koeffizienten als Liste mit Floats bilden
for i in range(len(nenk)) : # sympyfloats in echte Floats wandeln
nenk[i] = float(nenk[i])
if type(fe.as_numer_denom()[0]) == sy.core.add.Add:
zae = sy.Poly.from_expr(fe.as_numer_denom()[0]) # Zähler holen
zaek=sy.Poly.coeffs(zae) # Koeffizienten als Liste mit Floats bilden
for i in range(len(zaek)) :
zaek[i] = float(zaek[i])
print('zaehler poly')
else:
zaek = float(fe.as_numer_denom()[0])
#print('zaehler float')
return zaek, nenk
nn,dd = zzz(ff,0,0)
ge11 = ct.tf(nn,dd)
nn,dd = zzz(ff,0,1)
ge12 = ct.tf(nn,dd)
nn,dd = zzz(ff,1,0)
ge21 = ct.tf(nn,dd)
nn,dd = zzz(ff,1,1)
ge22 = ct.tf(nn,dd)
nume = [[ge11.num[0][0], ge12.num[0][0]], [ge21.num[0][0], ge22.num[0][0]]]
dene = [[ge11.den[0][0], ge12.den[0][0]], [ge21.den[0][0], ge22.den[0][0]]]
entfilter = ct.tf(nume,dene,name='entfilter',
inputs = ['y1','y2'], outputs = ['yf1','yf2'])
# ct.poles(entfilter[0,0]) # wenn Pole gewünscht
Rk3 = ct.interconnect([strecke, entfilter, reg1, reg2,
sum1, sum2], inplist=['w1', 'w2'], outlist=['x1','x2','yf1','yf2','y1','y2'])
# Plot Regelkreis-Sprungantwort vollständige dynamische Entkopplung
# Vergleiche mit Figure 2
t, yy = ct.input_output_response(Rk3, T, ww,solve_ivp_method='RK23')
plt.figure()
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel(' x1, x2')
plt.title('Sprungantwort Regelkreise mit vollst. dyn Ent')
plt.grid()
plt.figure()
plt.plot(t,yy[2,:],t,yy[3,:],t,yy[4,:],t,yy[5,:])
plt.xlabel('t [s]');plt.ylabel(' y1, y2')
plt.title('Sprungantwort Stellgrößen mit vollst. dyn Ent')
plt.grid()
# Modellredunktion für den Entkopplungsfilter
ge11rz = ct.balred(ct.tf2ss(ge11),2, method='matchdc')
ge12rz = ct.balred(ct.tf2ss(ge12),2, method='matchdc')
ge21rz = ct.balred(ct.tf2ss(ge21),2, method='matchdc')
ge22rz = ct.balred(ct.tf2ss(ge22),2, method='matchdc')
ge11r = ct.ss2tf(ge11rz)
ge12r = ct.ss2tf(ge12rz)
ge21r = ct.ss2tf(ge21rz)
ge22r = ct.ss2tf(ge22rz)
numer = [[ge11r.num[0][0], ge12r.num[0][0]], [ge21r.num[0][0], ge22r.num[0][0]]]
dener = [[ge11r.den[0][0], ge12r.den[0][0]], [ge21r.den[0][0], ge22r.den[0][0]]]
entfiltred = ct.tf(numer,dener,name='entfiltred',
inputs = ['y1','y2'], outputs = ['yf1','yf2'])
Rk4 = ct.interconnect([strecke, entfiltred, reg1, reg2,
sum1, sum2], inplist=['w1', 'w2'], outlist=['x1','x2','yf1','yf2','y1','y2'])
# Plot Regelkreis-Sprungantwort reduzierte dynamische Entkopplung
# Vergleiche mit Figure 6
t, yy = ct.input_output_response(Rk4, T, ww,solve_ivp_method='RK23')
plt.figure()
plt.plot(t,yy[0,:],t,yy[1,:])
plt.xlabel('t [s]');plt.ylabel(' x1, x2')
plt.title('Sprungantwort Regelkreise mit reduzierter dyn Ent')
plt.grid()
web702
Standardregelstrecke im Zustandsraum
Listing 7.2: Die Standardregelstrecke wird im Zustandsraum mit der Regelungsnormalform dargestellt und simuliert, gemäß Beispiel 7.7.1.
# -*- coding: utf-8 -*-
"""
Standardregelstrecke im Zustandsraum in RNF
Created on 12.4. 2019, mod. 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
a0 = 1./30.
a1 = 10./30.
a2 = 31./30.
b0 = 1./30.
A = np.array([[0, 1.,0 ],[ 0, 0, 1. ],[ -a0, -a1, -a2]])
b = np.array([[0.],[ 0. ],[ 1.]])
c = np.array([b0, 0,0 ])
stand = ct.ss(A, b, c, 0)
#T = np.linspace(0, 40,num=200)
#t, y, x = ct.forced_response(stand,T,U=1,X0=0,return_x=True)
t, y, x = ct.step_response(stand,return_x=True) # mit Zuständen
plt.plot(t,y)
plt.title('Standardregelstrecke im Zustandsraum, RNF')
plt.xlabel('t [s]'); plt.ylabel('h(t)')
plt.grid()
web703
DC-Motor im Zustandsraum
Listing 7.3: Der DC-Motor 2237 wird im Zustandsraum dargestellt, wobei die Simulation der Drehzahl erfolgt.
# -*- coding: utf-8 -*-
"""
DC-Motor im Zustandsraum Faulhaber Motor 2237
Created on 16.10. 2018, 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
Am = np.array([[-R/L, -2*pi*cPsi/L],
[cPsi/zpiJ, 0]])
bm = np.array([[ 1/L],[0]])
cm = np.array([ 0, 1])
C = np.array([[1, 0 ],
[0, 1 ]]) # Falls alle Zustände ausgegeben werden sollen
D = np.array([[ 0. ],[0]])
DCmot = ct.ss(Am, bm, cm, 0)
t, y, x = ct.step_response(DCmot, return_x=True)
plt.gcf().clear() # löscht den alten Plot
plt.plot(t,y)
#plt.plot(t,x[1,:],t,x[0,:]) # mit Stromverlauf
plt.title('DC-Motor im Zustandsraum')
plt.xlabel('t [s]'); plt.ylabel('n(t)')
plt.grid()
web704
Transformation auf Regelungsnormalform
Listing 7.4: Eine Transformation auf RNF wird mit Hilfe des Algorithmus aus Kapitel 7.7.3 durchgeführt und verglichen mit dem Ergebnis der Control Toolbox.
# -*- coding: utf-8 -*-
"""
Standardregelstrecke im Zustandsraum und Transformation auf RNF
Steuerbarkeitsmatrix S, Vergleich mit Control Lib
Created on 25.4. 2019, mod. 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
T1 = 2; T2 = 3; T3 = 5
A = np.array([[-1/T1, 0.,0 ],[ 1/T2, -1/T2, 0 ],[ 0, 1/T3, -1/T3]])
b = np.array([[1/T1],[ 0. ],[ 0]])
c = np.array([0, 0,1 ])
stand = ct.ss(A, b, c, 0)
# Berechnung Transformationsmatrix auf RNF
n = np.size(b)
AA = A
S = np.zeros((n,n)); T = np.zeros((n,n))
S[:,0] = b.T
S[:,1] = np.dot(A,b).T # Aufbau Steuerbarkeitsmatrix
for i in range(2,n):
AA = np.dot(AA,A)
S[:,i] = np.dot(AA,b).T
S1 = np.linalg.inv(S)
T[0,:] = S1[n-1,:] #letzte Zeile Übernehmen
for i in range(1,n):
T[i,:] = np.dot(T[i-1,:],A).T
print('T = ')
print(T,'\n')
Arnf = np.mat(T)*np.mat(A)*np.linalg.inv(T)
print('Arnf = ')
print(Arnf,'\n')
zu = ct.ss(A,b,c,0.0)
Arnf2, T2 = ct.canonical_form(zu)
print('System mit Control Lib auf RNF transformiert:')
print(Arnf2)
web705
Pendeldämpfung mit Zustandsregler
Listing 7.5:
# -*- coding: utf-8 -*-
"""
Kran im Zustandsraum Reglerentwurf mit Ackermann
Created on 24.8. 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
T = 0.1 ; l = 0.5 ; g = 9.81
Am = np.array([[0, 1, 0, 0],
[0, -1/T, 0, 0],
[0, 0, 0, 1],
[0, 1/(l*T), -g/l, 0]])
bm = np.array([[0],[ 1/T],[0],[-1/(l*T)]])
#cm = np.array([ 0, 1])
C = np.eye(4) # Falls alle Zustände ausgegeben werden sollen
D = np.array([[ 0. ],[0],[0],[0]])
linKran = ct.ss(Am, bm, C, D)
print('Pole : ',ct.poles(linKran))
# Zustandsreglerentwurf
pol = np.array([-1,-1,-1,-1])
rr = ct.acker(Am,bm,pol*2.0)
# nicht lineares Modell Kran
def kran_update(t, x, u, params={}):
# Reglerparameter übergeben, falls erforderlich
kp = params.get('kp', 0.05)
T1 = params.get('T1', 0.1)
l = params.get('l', 0.5)
# Zuweisung von Variablen und Zuständen
vsoll = u[0] # Vorgabe v
lage = x[0] # Zustand Lage
v = x[1] # Zustand v = Lage punkt
alpha = x[2] # Zustand alpha
alphap = x[3] # Zustand erste Ableitung alpha
# Berechnung Ableitungen
lagep = v
vp = -v/T1 + vsoll*kp/T1
alphapp = -np.sin(alpha)*9.81/l - vp*np.cos(alpha)/l # zweite Abl. alpha
return [lagep, vp, alphap , alphapp]
io_kran = ct.nlsys(
kran_update, None, inputs=['u'],
outputs=['lage','v','alpha','alphap'],
states=['lage','v', 'alpha', 'alphap'], name='kran',
params = { 'kp':1.0})
def x_rueck(t, x, u, params={}):
# parameter übergeben, falls erforderlich
vor = params.get('vor', 1)
# Zuweisung von Variablen und Zuständen
x0 = u[0]; x1 = u[1]; x2 = u[2]; x3 = u[3]
# Rückführung berechnen später Skalarprodukt
return rr[0,0]*x0 + rr[0,1]*x1 + rr[0,2]*x2 + rr[0,3]*x3
rueck = ct.nlsys(
None, x_rueck, name='rueck',
inputs = ['lage','v','alpha','alphap'], outputs = ['rue'],
params = {'vor':1})
sumblk = ct.summing_junction(inputs=['w', '-rue'],
output=['u'])
#Regelkreis aufbauen
Rk = ct.interconnect([io_kran, rueck, sumblk],
inplist=['w'], outlist=['lage','alpha','u'])
T = np.arange(0.0, 10.0, 0.01) # Simulationsdauer
ww = np.ones([len(T)])
t, yy = ct.input_output_response(Rk, T, rr[0,0]*ww)
#resp = control.input_output_response(Rk, T, ww) #rr[0,0]*ww)
# Umrechnung rad Grad
ur = 180/np.pi
plt.figure(1)
plt.plot(t,yy[0,:],t,ur*yy[1,:],t,yy[2,:])
plt.title('Kran mit Zustandsregler')
plt.xlabel('t [s]'); plt.ylabel('lage, alpha, y')
plt.grid(True)
web706
Modellreduktion mit balred()
Listing 7.6:
# -*- coding: utf-8 -*-
"""
Modellreduktion Standardregelstrecke
Created on 4.1. 2022, 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
g3 = ct.tf(1,[30, 31, 10, 1]) # Standardregelstrecke
s3 = ct.tf2ss(g3)
t, x = ct.step_response(g3)
plt.figure(1)
plt.plot(t,x)
b2 = ct.balred(s3,2, method='matchdc')
#b2 = balred(s3,2, method='truncate') #
gb2 = ct.ss2tf(b2)
print('G(s) = ',gb2)
print(ct.poles(gb2))
t, xx = ct.step_response(b2,t)
plt.plot(t,xx)
plt.title('Standardregelstrecke und balred Ersatzmodell')
plt.xlabel('t [s]')
plt.ylabel('x')
plt.grid()
g2 = ct.tf(1,[25,10,1]) #Ersatzmodell 2. Ordnung T1=T2=5
t, xxx = ct.step_response(g2,t)
plt.figure(2)
plt.plot(t,x,t,xxx)
plt.title('Standardregelstrecke und P-T2 Ersatzmodell')
plt.xlabel('t [s]')
plt.ylabel('x')
plt.grid()
