Es folgen Programme des Beispielkapitels.
web901
CODESYS-Programm PI-Regler
Listing 9.1 und Listing 9.2: CODESYS Hauptprogramm und Funktionsbaustein Beispiel 9.1.2
PROGRAM PLC_PRG
VAR
INIT:BYTE:=0;
Regler:PI_REGLER; (*Instanz des Funktions-
blocks PI_REGLER*)
A_IN AT %IW1:INT; (*Analog-Kanal 2*)
A_OUT AT %QW1:WORD; (*Analog-Kanal 2 Ausgang*)
Drehzahl:REAL; (*in U/min*)
Fuehrungsgroesse: REAL := 0; (*in U/min*)
Stellgroesse: REAL := 0; (*in Volt*)
Kp:REAL := 0;
Tn:REAL := 0;
T0:REAL := 0.05;
d0:REAL := 0;
d1:REAL := 0;
END_VAR
(**Initialisierung Regler und Zähler (falls gesteckt)******)
CASE INIT OF
0: (*Initialisierung Zähler (falls gesteckt)*)
INIT:=1;
RETURN;
1:
(*********Initialisierung des Reglers********************)
(*Die folgenden Variablen werden über Winfact gesetzt *)
Kp:= (INT_TO_REAL(%IW256))*0.001;
Tn:= (INT_TO_REAL(%IW257))*0.001;
Fuehrungsgroesse:=(INT_TO_REAL(%IW259))*0.1;
IF Kp = 0 THEN (* ist WinFACT nicht in Betrieb! *)
Kp:=6; 217
Tn:=0.9;
Fuehrungsgroesse := 1;
END_IF
(*Diskretisierung des Reglers*)
d0 := Kp *T0 /(2*Tn) + Kp;
d1 := Kp *T0 /(2*Tn) - Kp;
INIT:=2;
RETURN;
2: (* CASE 2: Ausführung wenn INIT = 2*)
(******** Einlesen Daten Analogeingang***************)
Drehzahl:=INT_TO_REAL(A_IN)*0.00030518509;
(******************Aufruf Funktionsblock PI-REGLER ********)
Regler(xk:=Drehzahl, wk:=Fuehrungsgroesse, rd0:=d0, rd1:=d1);
Stellgroesse:=Regler.yk;
%QW257:=REAL_TO_INT(Drehzahl* 1000);
(*Speichere Regelgröße in %QW257 für Boris*)
%QW258:=REAL_TO_INT(Stellgroesse*1000);
(*Speichere Stellgröße in %QW258 für Boris*)
(**********analoge Ausgabe*************************)
A_OUT:=REAL_TO_INT( Stellgroesse * 3276.7);
(*Umrechnung Real in Zahlenformat für die analoge Ausgabe*)
(* Zur Sicherheit keine negative Spannungsausgabe *)
IF A_OUT > 32767 THEN
A_OUT:= 0; (* 0V = Motor stopp *)
END_IF
END_CASE
web902
Algorithmus PI-Regler
PI-Regler ohne AWR-Maßnahme, nur mit Stellgrößenbeschränkung
FUNCTION_BLOCK PI_REGLER
VAR_INPUT
wk: REAL := 0;
xk: REAL := 0;
rd0:REAL := 0;
rd1:REAL := 0;
END_VAR
VAR_OUTPUT
yk: REAL;
END_VAR
VAR
Min_Stell: REAL := 0.0; (*kleinster Wert von y k *)
Max_Stell: REAL := 10.0; (*größter Wert von y k *) 218
ek: REAL := 0.0; (*Regelfehler *)
ek1: REAL := 0.0; (*Regelfehler k-1*)
yk1: REAL:=0.0; (*Stellgröße k-1*)
END_VAR
ek := wk - xk;
yk := yk1 + rd0*ek +rd1*ek1;
(*Stellgrößenbegrenzung ohne AWR*)
IF yk < Min_Stell THEN
yk := Min_Stell;
ELSIF yk > Max_Stell THEN
yk := Max_Stell;
END_IF
yk1:=yk;
ek1:=ek;
%QW260:=REAL_TO_INT(xk*100);
%QW261:=REAL_TO_INT(yk*100);
web903
DC-Motor Beispiel 9.1.4
Listing 9.3: Der DC-Motor wird inklusive Regelung im Zustandsraum dargestellt und simuliert.
# -*- coding: utf-8 -*-
"""
DC-Motor mit P-PI-Lageregler im Zustandsraum und G(s)
Kaskade gemäß Faulhaber MC5005
Created on 12.4. 2019, mod. 20.8.24
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
Tc = 1e-3 # Ersatzzeitkonstante Stromregelkreis
zpiJ = 47.119e-6 # Meine Berechnung 2pi7,5e-6
KL = 1 # Integrierbeiwert Lagestrecke
cPsi= 0.0318 # Datenblatt Faulhaber
Ks = cPsi/zpiJ # Verstärkung Teilstrecke
a= 16
KR= 1/(a*Ks*Tc) # Verstärkung Drehzahlregler gemäß SymOpt
Tnv = a*a*Tc # Nachstellzeit Drehzahlregler
Tnv = 2*a*Tc # Nachstellzeit Drehzahlregler modifiziert
Kreib=1e-5
a1 = KR*Ks; a2 = a1*Tnv; a3 = Tnv; a4 = a3*Tc
null = np.roots([a4*a2, 0, (a3*a1 - a2*a2 ), 0, -a1*a1]);
wi = np.max(null)
Kkrit = (1/KL)*(a3*a2*wi**4 - a2*a1*wi**2 -a4*a1*wi**4 +a2*a1*wi**2)/( (a2**2)*(wi**2) + a1**2 )
Kv= Kkrit/40 # Verstärkung Lageregler
Ar = np.array([[0, KL, 0, 0],
[0, -Kreib/zpiJ, cPsi/zpiJ, 0],
[-KR*Kv/Tc, -KR/Tc, -1/Tc, 1/Tc],
[-KR*Kv/Tnv, -KR/Tnv, 0, 0]])
br = np.array([[ 0. ],[ 0], [KR*Kv/Tc],[KR*Kv/Tnv]])
cr = np.array([1, 0, 0, 0])
C = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
D = np.array([[ 0. ],[ 0], [0],[0]])
DCkreis = ct.ss(Ar, br, cr, 0)
t, y, x =ct.step_response(DCkreis,return_x=True)
plt.subplot(3,1,1)
plt.plot(t,x[0,:])
plt.title('Lageregelkreis DC-Motor'); plt.ylabel('Lage')
plt.xlabel('t [s]')
plt.grid(True)
plt.subplot(3,1,2)
plt.plot(t,x[1,:],'g')
plt.xlabel('t [s]'); plt.ylabel('n(t)')
plt.grid(True)
plt.subplot(3,1,3)
plt.plot(t,x[2,:],'r')
plt.xlabel('t [s]'); plt.ylabel('Strom')
plt.grid(True)
plt.figure() # neue Plot-Funktion
DCkreis = ct.ss(Ar, br, C, D)
ct.step_response(DCkreis).plot(
title='Sprungantwort Lageregelkreis')
web904
Symbolische Berechnung RC3
Listing 9.4:
# -*- coding: utf-8 -*-
"""
Created on 12.11. 2024, mod 13.7.25
RC3 ohne Re Buchbeispiel
@author: philippsen
"""
import sympy as sy
import sympy.physics.control as co
s = sy.symbols("s", complex=True)
T, T2 = sy.symbols("T, T2", real=True, positive=True)
# RC3 Zustandsraumbeschreibung
AA = sy.Matrix([[-2/T, 1/T, 0], [1/T, -2/T, 1/T],[0, 1/T, -(1/T +1/T2)]])
bb = sy.Matrix([[1/T], [0],[0]])
cc = sy.Matrix([0,0,1])
dd = sy.Matrix([0])
RC3ss = co.StateSpace(AA, bb, cc.T, dd)
# Wandlung in Übertragungsfuktion
Gvons = RC3ss.rewrite(co.TransferFunction)[0][0]
# Alternative Rechnung für G(s)
sI = sy.Matrix([[s , 0, 0],[0, s, 0],[0, 0, s]])
g = cc.T*(sI - AA).inv()*bb
print('g = '); sy.pprint(g); print()
# Eingabe Bauteildaten
R4 = 100e3 # Widerstand 4 ohne Störung
Con = 100e-6 # Kondensator
R = 10e3 # Widerstand 1, 2 und 3
TT = R*Con
TT2 = Con*R4
GgT2 = Gvons.subs({T: TT})
Gg = Gvons.subs({T: TT, T2: TT2})
# Mit Python Control Lib
import numpy as np
import control as ct
import matplotlib.pyplot as plt
C = 100e-6 # Kond
R = 10e3 # Widerstand 1, 2 und 3
R4 = 100e3 # Widerstand 4
T = R*C
T2 = C*R4
Am = np.array([[-2/T, 1/T, 0],
[1/T, -2/T, 1/T],
[0, 1/T, -(1/T + 1/T2)]])
bm = np.array([[ 1/T],[0],[0]])
cm = np.array([ 0, 0, 1])
dm = 0
RC3 = ct.ss(Am, bm, cm, 0)
RC3g = ct.ss2tf(RC3)
# Reglerentwurf Betragsoptimum
a3 = 1; a2= RC3g.den[0][0][1]; a1 = RC3g.den[0][0][2] ; a0 = RC3g.den[0][0][3]
kp = (a1/2) * (a1**2 - a0*a2) /(a1*a2 - a0*a3) - a0/2 # Formel für Strecke 3. Ordnung
tn = a1/a0 - (a1*a2 - a0*a3)/(a1*a1 - a0*a2)
print('Entwurf Betragsoptimum Kp = ', kp,' Tn = ', tn)
#Reglerentwurf Tsumme
Ver = ct.dcgain(RC3g)
Tsum = sum(-1/ct.poles(RC3g))
print('Entwurf Tsum: Ver = ', Ver,' Tsum = ', Tsum)
kkp = 1/Ver; ttn = 0.7 * Tsum
kkp = 0.5/Ver; ttn = 0.5 * Tsum
print('Entwurf Tsum Kp = ', kkp,' Tn = ', ttn)
# Evo liefert bessere Ergebnisse
regler = ct.tf([kp*tn, kp], [tn, 0])
G0 = regler*RC3g
web905
Raspberry pico
Listing 9.5:
# Aufnahme Sprungantwort Regelkreis mit RC3 und pico
# Ohne besondere Massnahmen fuer exakte Echtzeit
# programmiert in microPython MIT-License, H.-W. Philippsen
import os
import machine
import time
schalt = machine.Pin(19,machine.Pin.OUT) # Schalten Lastwiderstand
schalt.value(False) # Rlast aus
pin = machine.Pin(25,machine.Pin.OUT) # LED-Pin
ana1 = machine.ADC(machine.Pin(26)) # Analoger Eingang 1
ana2 = machine.ADC(machine.Pin(27)) # Analoger Eingang 2
ana3 = machine.ADC(machine.Pin(28)) # Analoger Eingang 3
pwm0 = machine.PWM(machine.Pin(16)) # erzeuge PWM Objekt für Pin 16
pwm0.freq(10000); pwm0.duty_u16(0) # Vorgabe Frequenz duty cycle auf 0
T0 = 0.05
delay = T0 - 0.0028 # Rechenzeit von 2.8ms beruecksichtigen
um = 3.3/65535 # kalibriert 15.11.24
logfile = open("RC3picoReglog2Z3.txt", "w")
pin.value(True) # LED an
Kp = 2.1; Tn = 4.09 # Reglerparameter
ymax = 3.4 # Begrenzung Stellg in Volt
Kid= Kp*T0/(2*Tn) # Ki diskretisiert
Tikorr = 0.96*Tn/T0 # Korrekturfaktor I-Anteil
ykk =0; yk =0.0; yi1 = 0.0; ek1 = 0.0 # Anfangswerte
w1 = const( 1.0 ); w2 = const( 2.0 ) # Sollwerte
start = time.ticks_ms() # get millisecond counter
for i in range(600):
x = ana3.read_u16()*um # x3 = x
if i < 2000:
ek = w1 - x
else:
ek = w2 - x # eventuell Sollwertsprung
# Anfang PI_AWR_REGLER, yk in Volt
Intkorr = (ykk - yk)/Tikorr # Korrektur I-Anteil falls Begrenzung
yi = yi1 + Kid*ek + Kid*ek1 - Intkorr # Berechnung I-Anteil
ek1 = ek
yi1=yi
ykk = Kp*ek + yi # Stellgroesse = P + I
if ykk < 0 : # Begrenzung Stellgroesse
yk = 0
elif ykk > ymax :
yk = ymax
else:
yk = ykk # oder yk ohne Begrenzung
# Ende PI-Regler
dutyint = int(65535*yk/3.4) # Umrechnunung in Wertebereich 0-65535
pwm0.duty_u16(dutyint) # set duty cycle
s = str(x) + '\t' + str(yk) + '\n'
logfile.write(s)
print(x,yk)
if i==300:
schalt.value(True) # Rlast ein
time.sleep(delay) # Warten
delta = time.ticks_diff(time.ticks_ms(), start) # compute time difference
print(delta)
logfile.close() # Datei schliessen
pin.value(False) # Led ausschalten
pwm0.duty_u16(0) # AUS
schalt.value(False) # Rlast aus
web906
Beobachter für RC3
Listing 9.6:
# -*- coding: utf-8 -*-
"""
RC3 im Zustandsraum mit Beobachter
Wird gestört versagt Beobachter
Created on 12.11.24 für Buch ohne Re
@author: philippsen
"""
import numpy as np
import control as ct
import matplotlib.pyplot as plt
# Regelstrecke in Zustandsraumdarstellung
C = 100e-6 # Kond
R = 10e3 # Widerstand 1, 2 und 3
R4 = 100e3 # Widerstand 4
T = R*C
T2 = C*R4
Am = np.array([[-(2/T), 1/T, 0],
[1/T, -(2/T), 1/T],
[0, 1/T, -(1/T + 1/(T2))]])
bm = np.array([[ 1/T],[0],[0]])
cmt = np.array([ [0], [0], [1]]) # transponierter Ausgangsvektor
dm = 0
CC = np.identity(3)
D = np.zeros((3,1))
RC3 = ct.ss(Am, bm, CC, D) #, name='rc3',inputs='u', outputs=['y1','y2','y3'])
t, y = ct.step_response(RC3)
def RC3_update(t, x, u, params={}):
# Parameter übergeben, falls erforderlich
R4 = params.get('R4', 100e3)
T = params.get('T', 1.0)
C = params.get('C', 100e-6)
T2 = C*R4
x1 = x[0]
x2 = x[1]
x3 = x[2]
# Zuweisung von Variablen und Zuständen
uein = u[0] # Eingang
zR = u[1]
x1p = -(2/T)*x1 + x2/T + uein/T
x2p = x1/T -(2/T)*x2 +x3/T
x3p = x2/T -(1/T + 1/(T2-zR*C))*x3 # mit Störterm
return [x1p, x2p, x3p]
io_RC3 = ct.NonlinearIOSystem(
RC3_update, None, name='RC3',
inputs = ['u','zR'], outputs = ['x1','x2','x3'], states = ['x1','x2','x3'],
params = {'R4':100e3, 'T':1.0, 'C':100.0e-6})
# Beobachterentwurf
p = -1/0.7 # Vorgabe Dynamik/Pole
# Rückkoppelvektor mit Hilfe der Formel von Ackermann berechnen
h = ct.acker(Am,cmt,[p,p,p])
# Bildung Eingangsmatrix Beobachter
bo = np.zeros((3,2))
bo[:,0:]=bm
#bo[:,1:]=h.T
bo[:,1:]=h.reshape(3,1) # mod für 0.10.2
DD = np.zeros((3,2))
#Beobachtermodell
bmod = ct.ss(Am, bo, CC, DD, name='beob',inputs=['u','e'], outputs=['xb1','xb2','xb3'])
sumblk = ct.summing_junction(inputs=['x3', '-xb3'], output='e')
# Strecke mit Beobachter
Sb = ct.interconnect([io_RC3,sumblk, bmod
], inplist=['u','zR'], outlist=['x1','x2','x3','xb1','xb2','xb3','e'])
t=np.arange(0,30,0.01)
ww = np.zeros([2,len(t)])
for i in range(len(t)):
if t[i] > 0.0:
ww[0,i] = 1.0
else:
ww[0,i] = 0
if t[i] > 15.0: # Störung
ww[1,i] = 50e3
else:
ww[1,i] = 0
x0 = [0.4,0.4,0.4,0,0,0] # Vorgabe Anfangswerte
t, yy = ct.input_output_response(Sb, t, ww,x0,solve_ivp_method='BDF')
plt.figure()
plt.plot(t,yy[0:6,:].T)
plt.title('Beobachter für RC3')
plt.xlabel('t [s]'); plt.ylabel('alle x')
plt.grid()
modd = ct.c2d(bmod,0.1)
print(modd)
web907
Raspberry pico Beobachter RC3
Listing 9.7:
# Aufnahme Sprungantwort auf dem pico Entwicklerbord
# H.-W. Philippsen für Beobachter RC3
import machine, time
pin = machine.Pin(25,machine.Pin.OUT) # LED-Pin
schalt = machine.Pin(19,machine.Pin.OUT) # Schalten Lastwiderstand
ana1 = machine.ADC(machine.Pin(26)) # Analoger Eingang 1
ana2 = machine.ADC(machine.Pin(27)) # Analoger Eingang 2
ana3 = machine.ADC(machine.Pin(28)) # Analoger Eingang 3
um = 3.3/65535 # Umrechnung in Volt
pwm0 = machine.PWM(machine.Pin(16)) # create PWM object from a pin
pwm0.freq(10000) # set frequency
T0 = 0.1
delay = T0 - 0.009 # Rechenzeit von 9ms beruecksichtigen
logfile = open("picoBeo4.txt", "w")
xk = [0.0, 0.0, 0.0]
xk1 = [0.0, 0.0, 0.0]
# Modell mit Re
bd = [9.07739346e-02, 4.38766072e-03, 1.46952622e-04]
Ad = [[0.82274902, 0.08214118, 0.00422601],
[0.08214118, 0.82697503, 0.08594882],
[0.00422601, 0.08594882, 0.9001889]]
h = [-1.59942729e-01, 2.59171609e-01, -6.42183175e-02 ] # Rückkoppelvektor disk
e=0.0
pin.value(True) # LED an
duty = 30.30303/2 # 30.30303 in Prozent für 1Volt
print('duty = ',duty)
dutyint = int(65535*duty/100) # Umrechnunung in range 0-65535
start = time.ticks_ms() # get millisecond counter
pwm0.duty_u16(dutyint) # set duty cycle für u = 0,5
for i in range(100):
pin.value(not pin.value()) # LED blinken Strecke auf Anfangswert bringen
time.sleep(delay)
duty = 30.30303 # 30.30303 in Prozent für 1Volt
dutyint = int(65535*duty/100) # Umrechnunung in range 0-65535
pwm0.duty_u16(dutyint) # set duty cycle für u = 1
u = 1
start = time.ticks_ms() # get millisecond counter
for i in range(400):
pin.value(not pin.value()) # LED blinken
x1 = ana1.read_u16()*um
x2 = ana2.read_u16()*um
x3 = ana3.read_u16()*um
s = str(T0*i) + '\t' + str(x1) + '\t' + str(x2) + '\t' + str(x3) + '\t' + str(xk[0]) + \
'\t' + str(xk[1]) + '\t' + str(xk[2]) + '\t' + str(u) + '\n'
xk[0] = Ad[0][0]*xk1[0] + Ad[0][1]*xk1[1] + Ad[0][2]*xk1[2] + bd[0]*u + h[0]*e
xk[1] = Ad[1][0]*xk1[0] + Ad[1][1]*xk1[1] + Ad[1][2]*xk1[2] + bd[1]*u + h[1]*e
xk[2] = Ad[2][0]*xk1[0] + Ad[2][1]*xk1[1] + Ad[2][2]*xk1[2] + bd[2]*u + h[2]*e
e= x3 - xk[2]
xk1 = xk
print('x1 = ',x1 ,'xk1 = ',xk[0] ,'x2 = ',x2 ,'xk2 = ',xk[1] , 'x3 = ',x3 ,'xk3 = ',xk[2])
logfile.write(s)
if i==200:
schalt.value(True) # Rlast ein
time.sleep(delay) # Warten
logfile.close() # Datei schliessen
delta = time.ticks_diff(time.ticks_ms(), start) # compute time difference
print(delta)
pin.value(False) # Led ausschalten
pwm0.duty_u16(0) # AUS
