Einstieg in die Regelungstechnik mit Python

Fachbuch von Hans-Werner Philippsen

Kapitel 7

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 Toolbox und Berechnung der statischen Entkopplungsmatrix. Für die Berechnung von Mehrgrößensystem ist das Paket slycot erforderlich. Die Installation erfolgt im Anaconda-Prompt:

  1. conda config –append channels conda-forge
  2. conda install -c repa slycot
# -*- coding: utf-8 -*-
"""
Erzeugung eines Mehrgrößensystems mit der
Python Control Toolbox und Berechnung der statischen Entkopplung.
Mehrgrößensystem und PI-Regler slycot erforderlich
1. conda config --append channels conda-forge
2. conda install -c repa slycot

statische Entkopplung funktioniert nicht im Fall
integralwirkende Strecke!

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

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

#Gii(s) = 2/2s^2 + 3s +1    Gij = 0.5/2s+1
n11= [2.0]; d11= [2., 3., 1.]
n22= [2.0]; d22= [2., 3., 1.]
n12= [0.5]; d12= [2., 1.]
n21= [0.5]; d21= [2., 1.]
g11 = tf(n11,d11); g110 = dcgain(g11)
g12 = tf(n12,d12); g120 = dcgain(g12)
g22 = tf(n22,d22); g220 = dcgain(g22)
g21 = tf(n21,d21); g210 = dcgain(g21)
        
num = [[n11, n12], [n21, n22]]
den = [[d11, d12], [d21, d22]]
strecke = tf(num, den)
k= g110*g220/(g110*g220 - g210*g120)

ne = [[[k], [-k*g120/g110]], [[-k*g210/g220], [k]]]
de = [[[1.], [1.]], [[1.], [1.]]]
ent = tf(ne, de)
y, t = control.matlab.step(strecke)
# Plot zeigt Wirkung von Eingang 1 auf beide Ausgänge
# neu ab 0.9
plt.plot(t,y[:,0],t,y[:,1])
# Version < 0.9 plt.plot(t,y)
plt.title('Mehrgrößensystem')
plt.xlabel('t [s]'); plt.ylabel('h(t)')
plt.grid()
plt.figure()
yy, tt = control.matlab.step(ent*strecke)
# Plot zeigt Wirkung von Eingang 1 auf beide Ausgänge
plt.plot(tt,yy[:,0],tt,yy[:,1])
# für Version < 0.9 plt.plot(tt,yy)
plt.title('statisch entkoppeltes Mehrgrößensystem')
plt.xlabel('t [s]'); plt.ylabel('h(t)')
plt.grid()


web702

Standardregelstrecke im Zustandsraum

Listing 7.2: Die Standardregelstrecke wird im Zustandsraum mit der Regelungsnormalform dargestellt und simuliert. Modifizieren sie die Darstellung gemäß Beispiel 7.7.1.

# -*- coding: utf-8 -*-
"""
Standardregelstrecke im Zustandsraum in RNF
Created on 12.4. 2019
@author: philippsen
"""

import numpy as np
import control
from  control.matlab import *
from matplotlib.pyplot import *
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 = ss(A, b, c, 0)

y, t = control.matlab.step(stand)

plt.plot(t,y)
plt.title('Standardregelstrecke im Zustandsraum')
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
@author: philippsen
"""

import numpy as np; from math import pi
import control
from   control.matlab import *
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 = ss(Am, bm, cm, 0)

y, t = control.matlab.step(DCmot)

plt.gcf().clear() # löscht den alten Plot
plt.plot(t,y)
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 Toolbox

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

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

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 = 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 = ss(A,b,c,0.0)
Arnf2, T2 = control.canonical_form(zu)
print('System mit Control Toolbox 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
@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

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, 0,0,1])
C  = np.eye(4)   # Falls alle Zustände ausgegeben werden sollen
D  = np.array([[ 0 ],[0],[0],[0]])
linKran = ss(Am, bm, C, D)

print(pole(linKran))
# Zustandsreglerentwurf
pol  = np.array([-1,-1,-1,-1])
rr = control.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          #  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 = control.NonlinearIOSystem(
    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 = control.NonlinearIOSystem(
    None, x_rueck, name='rueck',
    inputs = ('lage','v','alpha','alphap'), outputs = ('rue'),
    params = {'vor':1}) 
 
sumblk = control.summing_junction(inputs=['w', '-rue'], output='u')

#Regelkreis aufbauen
Rk = control.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 = control.input_output_response(Rk, T, 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
@author: philippsen
"""

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

g3 = tf(1,[30, 31, 10, 1])  # Standardregelstrecke
s3 = tf2ss(g3)

x,t = step(g3)
plt.figure(1)
plt.plot(t,x)

b2 = balred(s3,2, method='matchdc') 
#b2 = balred(s3,2, method='truncate') #
gb2 = ss2tf(b2)
print('G(s) = ',gb2)
print(pole(gb2))
xx,t = step(b2,t)
plt.plot(t,xx)

plt.title('Standardregelstrecke und balred Ersatzmodell')
plt.xlabel('t [s]')
plt.ylabel('x')
plt.grid()

g2 = tf(1,[25,10,1]) #Ersatzmodell 2. Ordnung T1=T2=5
xxx,t = step(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()

Theme von Anders Norén