NASTAWIANIE REGULATORA PID
12.1. PID dla układu bez szumu
12.2. PID o czasie regulacji jak P
12.3. PID jako PI+PD
12.4. PID dla układu z szumem
Rozważymy kilka wariantów doboru nastaw regulatora PID w zależności od tego, czy w układzie występują szumy pomiarowe oraz czy PID ma zapewniać taki sam, czy krótszy czas regulacji niż regulator P.
12.1.
PID dla układu bez szumu
1
T s
d
PID:
k p 1+
+ Td
T s
i
s + 1
D
D ≥ 5 :
2
T
k
s
p
( i +
T
2
)1
dla Td = i , ∆ =
1
0
2
k p T T s
T s
T
s
i
d
+ i +1
k
1
T s
p
+
+ d
=
=
4
i
T s
T
s
i
i
k p ( T s
T s
T
1
+ )
1 ( 2 + )
1
dla Td < i , ∆ > 0
T
s
4
i
• Charakterystyka czę stotliwoś ciowa
Uwaga. Matlab wymaga, aby stopień licznika nie przekraczał stopnia mianownika. Dlatego pozostawiamy D, ale przyjmiemy je dostatecznie duże, np. D = 10.
Matlab
kp=1; Ti=2; Td=0.5; D=10;
–
∆ ≅ 0
15
w=logspace(-1,1,100);
dB
dB
−
+
20
20
10
dek
dek
l=kp*[Ti*Td Ti 1];
– w przybliżeniu
m=Ti*[Td/D 1 0];
5
∆ = 0
∆ > 0
[M, F]=bode(l,m,w);
0
dB=20*log10(M);
-5 -1
0
1
Ti=4; Td=0.25;
–
∆ > 0
10
10
10
l=kp*[Ti*Td Ti 1];
100
m=Ti*[Td/D 1 0];
50
[M1, F1]=bode(l,m,w);
∆ = 0
dB1=20*log10(M1);
0
subplot(211)
∆ > 0
-50
semilogx(w,dB,w,dB1),grid
subplot(212)
-100 -1
0
1
10
10
10
semilogx(w,F,w,F1),grid
∆ =
∆ > 0
dB
0
dB
-20
+20
-20
+20
ω
2
1
ω
1
T 1
T
i
T
2
F
F
+90
+90
ω
ω
-90
-90
Dla niskich częstotliwości zachowanie regulatora PID jest podobne do PI, a dla wysokich do PD.
Projektowanie
Problem
–
dane: p%, tr
p
→ ξ →
%
PM
Dane
− szukane: k , T , T
p
i
d
t
→ ω
r
1
pośrednie
1
k p 1
k
PID:
k
1+
+ T s = k +
+ k T s = k
i
+
+ k s
p
d
p
p
d
p
d
{
T s
T s
s
i
i
{
k
k
d
i
Wzory:
j(−180+ PM )
j P
∠ ID
j G
∠
G
( ω )
ω
ω
1 =
(
)
1 ⋅
(
)
1 = 1 ⋅
= |
|⋅
⋅
⋅
otw j
PID j
G j
e
PID e
M e
2
k s + k + k s
p
i
d
j G
∠
j (−180+ PM )
⋅ M ⋅ e
= 1⋅ e
s
jω1
ω
2
1
j (−180+ PM −∠ G) k −ω k + jω k = j e
i
1
d
1
p
M
θ ≡ −180 + PM − ∠ G = P
∠ ID( jω )
– kąt regulatora
1
θ
e j = cosθ + j sin θ
– wzór Eulera
ω
2
1
ki − ω kd +
sin
1
θ + jω kp − j cosθ = 0
1
1
M
M
1
4
4
4
2
4
4
4
3
1
4
4
4
2
4
4
4
3
Re=0
Im=0
cosθ
sin θ
ki
k =
,
k =
+
p
d
2
M
ω M ω
1
1
Wniosek. Otrzymaliśmy dwa warunki projektowe z dwóch danych – PM, ω1. Chcąc wyznaczyć trzy nastawy ( k , T , T ) trzeba wprowadzić dodatkowy warunek.
p
i
d
Rozważymy dwa przypadki:
ki = 0
–
regulator PD
i
T
T =
d
–
regulator PID o podwójnym zerze
4
1)
PD
– ki = 0
cosθ
sin θ
k
tg
d
θ
k
k
T
,
dla
θ > 0 .
p =
,
d =
→
d =
=
M
ω M
k p
ω
1
1
Przykład 12.1.
Regulator PD dla obiektu całkującego o dwóch stałych czasowych.
2
G( s) =
, p
≤
% = 16.3,
tr
4
s( s + )
1 ( .
0 5 s + )
1
o
p
= 16 3
.
→ ξ = 5
.
0
→ PM = 50
%
8
8
t
–
przyjęto ω = .
1 7
r =
→ ω =
→ ω ≥ .
1 67
ω tg
1
PM
t tg
1
PM
1
1
r
Matlab – jako kalkulator
s=1.7*sqrt(-1);
–
jω1
G=2/(s*(s+1)*(0.5*s+1))
-0.44773 + 0.078133i
M=abs(G)
0.45449
fi=angle(G)*180/pi
–
∠ G = 9
− 0 − arctgω − arctg( 5
.
0 ω)
170.1
teta=-180+50-170.1
– Matlab podał fazę powiększoną o 3600. Faktycznie faza
-300.1
wynosi 170.1 - 360 = -189.90
– uwzględniono korektę Matlaba (doprowadzenie do kąta ostrego) 59.9
kp=cos(teta*pi/180)/M
1.1035
kd=sin(teta*pi/180)/(M*1.7)
PD :
.
1 1 1
( + 0
.
1 1 s)
1.1197
Td=kd/kp
1.0148
Uwaga. Z matematycznego punktu widzenia nie było potrzeby zmniejszania fazy o 3600, ponieważ ze względu na okresowość sin θ
( ± 360o ) = sin θ ,
cos θ
( ± 360o ) = cosθ
1.4
Symulacja
1.2
l=2*kp*[Td 1];
1
m=conv([1 1 0],[0.5 1]);
l=[0 0 l];
0.8
t=0:0.1:8;
y1=step(l, l+m,t);
0.6
plot(t,y1), grid
max(y1)
0.4
0.2
18.06%
0
0
1
2
3
4
5
6
7
8
T
2) PID dla
i
d
T =
– podwójne zero (Ziegler-Nichols)
4
k p
k p
k
k
1 k p
k =
d
→
d =
i
,
k =
→
d
k pTd
i
T =
,
T =
d
T
i
T
ki
k p
i
k p
4 ki
d
T = 4
cosθ
sin θ
k
Uwzględniamy poprzednie wzory, tzn.: k =
i
=
p
,
k
+
d
M
2
ω1 M ω1
M sin θ
k
1 cosθ 1
i
+
=
cosθ ω M
2
ω 4 M k
1
1
i
θ
θ
2
ω
2
sin
cos2
2
k
ω
ω
→
1
∆ =
i + ki
−
= 0
1
1
M
4
2
M
2
M
ω
k p
2
co θ
s
1
k
→ Ti =
=
i =
(1− sinθ )
2 M
ki
ω 1− sin
1
θ
2
cosθ
T
k
=
,
T =
,
i
=
o
θ ∈
p
i
d
T
dla
(0, 90 )
M
ω 1− sinθ
4
1
• c.d. – przykład 12.1
2
cos 59 9
.
4 3
. 7
k
Td =
=
i
T =
=
p = 1 1
. ,
4.37 ,
1 0
. 9
1 7
. 1 − sin 59.9
4
Matlab
1.4
l=2*1.1*[4.37*1.09 4.37 1];
m=4.37*conv([1 1 0 0],[0.5 1]);
1.2
y2
l=[0 0 l];
y2=step(l,l+m,t);
1
max(y2)
plot(t,y1,t,y2), grid
0.8
y1
0.6
23.49%
0.4
0.2
0
0
1
2
3
4
5
6
7
8
12.2.
PID o czasie regulacji jak P
Jest to podstawowy sposób doboru nastaw regulatora PID.
Problem
Dodatkowo tutaj przyjmujemy:
–
dane: p%, tr,PID = tr,P
D ≥ 5
– niski poziom szumu (j.p.)
− szukane: k , T , T
p
i
d
i
T
d
T =
– podwójne zero
4
Projektowanie
• p
→
ξ →
%
PM
•
8
t =
→
ω
= ω
r
ω
,
1 P
,
1 PID
– taka sama częstotliwość dla układu tg PM
1
z regulatorem P i PID
j∠ G( jω )
j(−180+
)
G
(
1
PM
ω =
ω =
ω ⋅
= ⋅
otw j
1)
k G(
p
j 1) k G(
p
j 1) e
1 e
ω :
∠ ( ω )
ω
1 = −180 +
1
G j
PM
– równanie dla
1
ω
→
1
8
( ω )
→
=
=
1
≡
1
G j
M
k p
,
tr, P
M
ω tg PM
1
• Układ z regulatorem PID
j P
∠ ID
j G
∠
j(−180+
)
G
(
PM
ω =
ω ⋅
ω =
ω ⋅
ω ⋅
⋅
= ⋅
otw j 1)
PID( j 1) G( j 1) PID( j 1) G( j 1) e e
1 e
P
∠ ID( jω )
ω
1 + ∠ G( j
)
1
= −180 + PM
Ponieważ ω1 ma pozostać niezmienione, to musi zachodzić warunek: o
∠ PID( jω ) = θ = 0
1
Dla
o
θ = 0 poprzednie wzory przyjmują postać: 1
2
T
k
=
,
T =
,
i
=
p
i
d
T
M
ω
4
1
Przykład 12.2. Obiekt dwuinercyjny z opóźnieniem, regulator PID taki, że tr,PID = tr,P
−
e s
P/PID
p% = 20.5
(10 s + )2
1
• p = 20 5
.
→ ξ =
→
o
PM =
%
.
0 45
45
• Regulator P
G
∠ ( jω )
1 = −180 + PM
180
− 2arctg 1
( 0ω ) − ω
− 45 = 1
− 80
ω
1
1 π
– równianie dla
1
w=0.1:0.01:1;
fi=(-2*atan(10*w)-w)*180/pi-45;
[w' fi']
...
8
0.19
-180.37
– ω = .
0 19
→
tr P =
= 4 .
2 1
1
,
.
0 19 ⋅ tg 45
s=0.19*sqrt(-1);
–
jω1
M=abs(1/((10*s+1)^2))
0.21692
kp=1/M
4.61
3
2
s
s
s
Obiekt
−
+
− +1
−
120
10
2
e s ≅
lp=[-1/120, 1/10, -1/2, 1];
– Pade
3
2
s
s
s
mp=[1/120, 1/10, 1/2, 1];
+
+ +1
lo=lp;
– licznik obiektu
120
10
2
mo=conv([100, 20, 1],mp);
l=kp*[0 0 lo];
m=mo;
t=0:1:100;
1.4
yP=step(l,l+m,t);
max(yP/yP(100))
1.2
PID
34.26%
1
•
0.8
Regulator PID
P
0.6
kp=1/M; Ti=2/0.19; Td=Ti/4;
0.4
l=kp*conv([Ti*Td, Ti, 1], lo);
m=Ti*[mo 0];
0.2
l=[0 l];
yPID=step(l,l+m,t);
0
plot(t,yP,t,yPID), grid
-0.2
max(yPID)
0
10
20
30
40
50
60
70
80
90
100
26.19%
• Regulator PI – porównanie z P, PID
Kroki projektowania (pkt. 11.3):
1.
α = 0.2
– wybór z przedziału 0.1...1.0 (0.2 oznacza dość silne działanie P) 1
2.
∠ PI(α = 0.2) = 9
− 0 + arctg
= 1
− 1.3
0.2
3.
ω
G
∠
jω
= −
+ PM − P
∠ I α
1 :
(
)
180
( )
1
– powiększenie zapasu fazy
180
o
− 2arctg 1
( 0ω ) − ω
− 45 −11.3 = −180
1
1 π
ω = .
0 156
1
– wyznaczone, jak w przykładzie powyżej
z = αω = 2
.
0 ⋅ 1
.
0 56 =
0
.
0 31
i
T =
1
,
32 2
.
1
5.
M = G( jω ) =
=
1
(10⋅ 1.
0 56)
2
.
0 91
2 +1
1
s + 0
.
0 31
1
k =
= 3
.
3 6
PI:
.
3 36
= 3
.
3 6 1 +
s
3 .
2 2 s
M 1
2
+α
8
6.
tr PI =
= 5 .
1 3
,
1
.
0 56 ⋅ tg 45
1.4
Matlab
PID
1.2
l=3.36*conv([1 0.031], lo);
1
m=[mo 0];
PI
l=[0 0 l];
0.8
P
yPI=step(l,l+m,t);
0.6
plot(t,yP,t,yPI,t,yPID), grid
0.4
0.2
0
-0.2
0
10
20
30
40
50
60
70
80
90
100
12.3. PID jako PI+PD
1
T s
s + z
s + z
PID:
d
1
2
k 1+
+
=
⋅
p
1
k
k 2
T s
Td
∆>0
+
+
i
s 1
s
s
p 2
4
1
4
2 3
4
1
4
2
3
D
PI
PD
Projektowanie
• PI
s + z
k
1
G( s)
1
s
α
1
, ω
∠ = −180 +
− ∠
z = αω
k =
1 :
G
PM
PI ,
1 , M,
– jak powyżej
2
M 1 + α
G΄
s + z
+
2
s
z
k
1
2
k
G( )
1
s
s + p 2
s
s + z
2
2
ω
4
z +
1
ω + p
2
ω tg
1
θ
=
=
1
2
=
1 ,
G
∠ ' , z 2
,
p
,
k
– pkt. 11.1
t
2
2
1
z 2
−
θ
2
2
ω +
r
M
z
ω tg
1
2
1
s + z
gdzie θ , M dotyczą G'= k 1 G
1
s
Przykład 12.3. Obiekt i przeregulowanie jak poprzednio (przykład 12.2), ale tr = 20 , czyli 2-krotnie krócej (było tr = 4 .
2 1).
1)
PI
–
jeżeli α = 0.2 , to wynik jak poprzednio, tzn.
s + 0 0
. 31
PI: 3.36
s
2)
PD
G΄
s + z
s
2
s +
−
0 0
. 31
e
k 2
3 3
. 6
s + p 2
s
(10 s + )2
1
s + z
p
= 20 5
.
→ ξ =
→
o
PM =
%
.
0 45
45
8
t
→ ω =
= 4
.
0
r = 20
1
20 tg 45
ω
180
1
o
∠ G'( jω ) = 9
− 0 + arctg
− 2arctg 1
( 0ω ) − ω
= −179 3
.
1
1
1
0.031
π ω = .04
1
o
θ = 1
− 80 + PM − ∠ G'= 1
− 80 + 45 − ( 1
− 79. )
3 = 44.3
4
4
z =
=
= 2
.
0
2
t
20
r
.
0 2 + 4
.
0 tg 4 .
4 3
2
1
p =
=
= 1
.
1 53
2
z
0.2
1
2
−
tgθ
1 −
tg 44 3
.
ω
0.4
1
2
ω + 0.0312
1
1
M = 3 3
. 6
=
ω1
(10ω +
1 )
0 1
. 98
2
1 ω =0 4.
1
1
2
2
ω +
1
p 2
k =
= 13 7
.
2
2
2
M
ω +
1
z 2
s +
−
0.2
s + 0.031
e s
G
( s)
otw
= 13 7
.
3 3
. 6
s +1.153
s
1
4
4
4
4
4
2
4
4
4
4
4
3 (10 s +
)2
1
PID = PI + PD
Matlab
1.4
lr=13.7*3.36*conv([1 0.2],[1 0.031]); 1.2
l=conv(lr,lo);
1
l=[0 0 l];
0.8
m=conv([1 1.153 0], mo);
t=0:0.6:60;
0.6
y=step(l,l+m,t);
0.4
plot(t,y), grid
max(y)
0.2
22.5%
0
-0.20
10
20
30
40
50
60
12.4. PID dla układu z szumem
• Założenia
1) Wartość D jest wybierana, np. jako 0.5, 1, 2, zależnie od poziomu szumu 2) t
=
r , PID
tr, P – czas regulacji jak dla regulatora P
3) PID ma mieć podwójne zero ( ∆ = 0 w liczniku)
• Transmitancja
1 2
T
T T 1
d
i d
+ s + Ti +
s +
1
T s
d
D
D
1
– PID:
k p 1+
+
= k p
T s
Td
T
i
s + 1
T s
d s
i
+
D
D
1
– W liczniku ∆ = 0 występuje dla T = T ⋅ D 2
+1− 2
(
+ )
1
=
d
i
( D
D
) T
D
i
1
4
4
4
4
2
4
4
4
4
3
d
1
d
D
0.5
1
2
5
8
d
7.46
5.828
4.949
4.391
4.246
D = 5 – Siemens
ν
2.36
3.414
5.448
11.47
17.48
D = 8 – Honeywell
ω '1
2.552
1.55
1.254
1.102
1.062
T
Wniosek. d dąży do 4 wraz ze wzrostem D (czyli i
T →
d
)
4
T
–
Pierwiastek podwójny z dla
i
T =
d
4
d
1
d +
i
T + T
1
z =
D
=
D
1
T
1
2
i
2 1 +
i
T d
T
1 +
D
D
–
Licznik
1
T
T T 1 +
⋅ +
=
+ ⋅ +
i d
( s z)
2
2
1
i
1
( s z)2
D
T
d
D
T
i
=
d
d
–
Mianownik
T
T
2
D
T
dD
T s d
i
s +
1 =
d
i
T
s s +
=
i
s s +
D
D
d
T
T
dD
i
T
Td =
i
d
ν z
ν = 1 dD = D +1
2
– ν rośnie wraz z D (tabela) z
1
i
T
1 +
ν stanowi daną pośrednią
dD
( s + z)2
( s + z)2
–
PID:
k ( D
k
p
+ )
1
=
– niewiadome k, z
s( s + z
ν )
s( s + z
ν )
(ν jest dane na podstawie D)
– dla P
k p
G( s)
p
→ ξ →
%
PM
j(−180+ PM )
j G
∠ ( ω )
1
G
(
= ⋅
=
otw jω
j
1)
1 e
k pMe
ω
∠ ( ω )
ω
1 = −180 +
1 :
G j
PM
– równanie dla
1
1
8
k
=
=
p
,
tr
M
ω t P
g M
1
• Zero z
– z warunku t
≡
ω
r , PID
tr, P
( 1 j.w.)
PID
G( s)
∠
+∠
−
+
–
j( PID
G)
j( 180
)
G
(
PM
ω =
⋅ ⋅
= ⋅
otw j 1)
| PID | M e
1 e
– Aby ω
∠ PID jω =
1 nie uległo zmianie, to musi zachodzić (
)
0
1
. Z równania tego
oblicza się z. Równanie ma postać: ω
ω
∠ PID( jω ) = − 90 + 2arctg 1 − arctg 1 = 0
1
z
ν z
ω
Nowa niewiadoma
1
ω =
1'
– częstotliwość względem zera z
z
∠ PID( jω ') = 0
1
ω
ω '
'
1
−
+
ω −
=
ω
1 :
90
2 arctg
' arctg
0
1
ν
– równanie dla
'
1
Mając ω '
1 obliczamy
ω1
z = ω '1
| PID( jω ) | ⋅ M = 1
1
– M otrzymano dla regulatora P
2
ω + 2
z
ω '2 +1
| PID( jω
1
1) =
1
| k
= k
2
ω ω
ν
2
2
1
1 +
2
( z)
ω ' ω ' ν
+
1
1
2
1 4
4 2
4
4 3
( s'+ )
1
M
- moduł
r
s' ( s' ν
+ ) s= j 1ω
1 ω ' ω '2 ν
1
1
+ 2
1
k =
=
M
ω '2
1
+1
M ⋅ M r
Przykład 12.4. Obiekt i przeregulowanie jak poprzednio (przykład 12.2), t
≡
r , PID
tr, P ,
założone D = 2.
−
( s + z)2
e s
k
2
p
= 20 5
.
%
s( s +ν z)
1
( 0 s + )
1
• ω1 i M jak dla P – zob. przykład 12.2
ω = .
0 19
t
M =
r =
1
,
4 .
2 1,
0 2
. 169
• d, ν , ω '1, Mr D = 2 →
d = 4.949,
ν = 5.448 (tabela)
Matlab
– równanie ∠ PID( jω ' ) = 0
ω
1
dla
'
1
ni=5.448;
lr=[1 2 1];
2
mr=[1 ni 0];
( s'+ )
1
PID:
bez k
wp=logspace(-1,1,500);
s'( s' ν
+ )
[Mr, Fr]=bode(lr, mr, wp);
[wp' Mr Fr]
...
1.2537
0.3669
-0.1135
ω '
≅
ω =
1
Mr
0
' 1.2537
1
Uwaga. Dotąd wszystko wynika z przyjętego D = 2.
ω
.
0 19
1
z =
=
= .
0 1515
ω '
.
1 2537
1
1
1
d +
4 9
. 49 +
1
1
D
2
i
T =
=
= 12 0
.
z
1
0.1515
1
2 1 +
2 1 +
D
2
T
12
T = i
d
=
= 2 4
. 2
d
4 9
. 49
• k, kp
1
1
k =
=
= 12 5
. 6
M ⋅ M
.
0 2169
r
⋅ .
0 3669
= k
k p
= 4 1
. 9
D + 1
( s + 1
.
0 51 )
5 2
1
2 42
.
s
PID:
1 .
2 56
= 4 19
.
1+
+
s( s + .
0 8 )
2
12
2 4
. 2
s
s + 1
{
2
z
ν
Matlab
1.4
1.2
lr=12.56*conv([1 0.1515],[1 0.1515]); 1
l=conv(lr,lo);
l=[0 0 l];
0.8
m=conv([1 0.82 0], mo);
0.6
y=step(l,l+m,t);
plot(t,y), grid
0.4
max(y)
0.2
1.2554
25.5%
0
-0.2
0
10
20
30
40
50
60