Wykład 9

CHARAKTERYSTYKI CZĘSTOTLIWOŚCIOWE

9.1.

Generowanie charakterystyk

9.2.

Zapas fazy i zapas modułu

9.3.

Przeregulowanie a zapas fazy

Czas regulacji a częstotliwość

9.4.

Logarytmiczne charakterystyki Bodego

9.1.

Generowanie charakterystyk częstotliwościowych Wprowadzenie

• Historia

1932 - Nyquist - badanie stabilności na podstawie charakterystyk częstotliwościowych 1942 - Bode - charakterystyki logarytmiczne i dobór wzmocnienia 1943 - Ziegler i Nichols – eksperymentalny dobór nastaw regulatorów PID po doprowadzeniu do oscylacji granicznych (granica stabilności)

• Cechy metod czę stotliwoś ciowych 1. Do projektowania służą charakterystyki częstotliwościowe a nie transmitancja, nie ma więc ograniczeń na rząd obiektu (inaczej niż w liniach pierwiastkowych) 2. Charakterystyki częstotliwościowe otrzymuje się eksperymentalnie za pomocą szybkiej transformaty Fouriera (FFT) pobudzając obiekt sygnałem o rosnącej częstotliwości

3. Uzupełnienie układu sterowania o filtry eliminujące zakłócenia nie komplikuje projektowania (inaczej niż w liniach pierwiastkowych) 4. Automatyczne strojenie metodą częstotliwościową daje lepsze rezultaty niż strojenie na podstawie odpowiedzi skokowej – gwarantowana stabilna praca 5. Metody częstotliwościowe są rozpowszechnione w pokrewnych dziedzinach – teoria sygnałów, telekomunikacja, elektronika.

- chirp

FFT – Fast Fourier Transform

moduł – M(ω) F(ω) - faza

Wady

-

Analiza i projektowanie odbywa się w dziedzinie częstotliwości a nie czasu. Nie występują więc wprost tak naturalne pojęcia jak stała czasowa, przeregulowanie, czas regulacji itp.

-

Nie można stosować najprostszego sposobu doboru nastaw jakim jest eliminacja stałej czasowej.

-

Opracowanie wyników eksperymentu częstotliwościowego i projektowanie wymaga pakietu CAD (w przypadku linii pierwiastkowych prostsze układy udało się projektować ręcznie – Siemens, Honeywell).

Instrukcje Matlaba

• w = logspace(d1,d2,n) - generowanie n punktów częstotliwości ω rozmieszczonych równomiernie w skali logarytmicznej w przedziale 10d1...10d2

w =logspace(d1,d2)

-

standardowo 50 punktów; np. logspace(-1,1) wygeneruje 50

punktów w przedziale 0.1...10

•

l

[M,F] = bode(l,m,w)

-

wyznaczenie modułu M i fazy F transmitancji dla

m

częstotliwości ω, gdzie F jest w stopniach (wyjątkowo w Matlabie)

• subplot(211)

-

wybór górnej połowy ekranu do umieszczenia wykresu semilogx(w,M), grid

-

wykres modułu w skali półlogarytmicznej

subplot(212)

-

wybór dolnej połowy ekranu

semilogx(w,F), grid

-

wykres fazy

clg

-

ekran standardowy (następny wykres pojedynczy).

• Zalecenia

1. W typowych problemach wystarczy 50 punktów na dekadę.

2. Przedział 10d1...10d2 powinien objąć częstotliwości charakterystyczne lub graniczne, którymi są odwrotności najmniejszej i największej stałej czasowej.

Transmitancja II rzędu

2

ω

G( s)

n

= 2

2

s + 2ξω s + ω

n

n

2ξω nω

j(− arctg

)

2

2

2

2

ω n ω

−

1

4

4

2

4

4

3

ω n

ω

G( s = jω) =

=

n

e

F(ω) - faza

2

2

2

2 2

2

− ω + j 2ξω nω + ω n

(ω n − ω ) + (2ξω nω)

1

4

4

4

4

2

4

4

4

4

3

M(ω) - moduł

Matlab

-

ω n = 1

l = 1

w = logspace(-1,1,100);

ksi = 1

m = [1 2*ksi 1]

[M, F] = bode(l,m,w);

subplot(211)

semilogx(w,M), grid

subplot(212)

semilogx(w,F), grid

ksi = 0.6

..............

ksi = 0.4

..............

ksi = 1.5

..............

Uwaga. Częstotliwością charakterystyczną jest ω n = 1. Przedział 0.1...10 w logspace (-1,1,100) obejmuje ją po dekadzie z lewej i prawej strony.

Wnioski

1. Wzrost tłumienia ξ powoduje

-

zmniejszenie modułu szczytowego Mp ( peak) i odpowiadającej mu częstotliwości ω p

-

zmniejszenie pasma przenoszenia BW ( band-width) 2. Częstotliwość ω p jest nieco niższa niż częstotliwość naturalna ω n = 1

3. Faza F ustala się na wartości -180°, ponieważ stopnie licznika i mianownika różnią się o 2 (2 ⋅ (−90o ) = 1

− 80o ) .

Wzrost rzędu a charakterystyki

1

1

1

G( s) =

,

,

2

3

4

( s + )

1

( s + )

1

( s + )

1

Matlab

l = 1

m = [1 2 1]

w = logspace(-1,1,100);

[M, F]=bode(l,m,w);

m1=[1 3 3 1]

[M1, F1]=bode(l,m1,w);

m2=conv(m,m)

[M2, F2]=bode(l,m2,w);

subplot(211)

semilogx(w,M,w,M1,w,M2),grid

subplot(212)

semilogx(w,F,w,F1,w,F2),grid

Wniosek. Ze wzrostem rzędu, a w zasadzie różnicy między stopniem licznika i mianownika zachodzi:

-

charakterystyka modułowa maleje coraz gwałtowniej

-

charakterystyka fazowa spada poniżej -180° (do o

( n

n

)

l −

)

m

⋅90

9.2.

Zapas fazy i zapas modułu

Serwomechanizm napięciowy

k

2

k

ω

G

( s)

T

n

=

=

=

,

otw

s( Ts + )

1

1

s( s + 2ξω )

s( s +

)

n

T

2

k

1

gdzie ω =

,

2ξω =

n

n

T

T

2

ω n

2

s( s + 2ξω )

ω

G

( s)

n

n

=

=

- standardowy układ II rzędu

zam

2

2

2

ω

s + 2ξω s + ω

1

n

n

n

+ s( s + 2ξω ) n

• Charakterystyki układu otwartego

2

2

ω

ω

ω

o

G

( s = j

n

n

(−90

arctan

)

ω

−

=

=

otw

)

e j

ξω

jω( jω + 2ξω )

2

2

2

n

n

ω ω + (2ξω n

1

4

4

4

2

4

4

4

3

)

1

4

4

4

2

4

4

4

3

F

M

Dane: ξ = ,

1

ω n = 1 (czyli p% = 0, tr = 4 ) Matlab

l = 1

m = [1 2 0]

w = logspace(-1,1,100);

[M, F] = bode(l,m,w);

subplot(211)

semilogx(w,M), grid

subplot(212)

semilogx(w,F), grid

[ w’

M

F ]

...........................

0.4863

0.9991

-103.6

ω1

≈1

F(ω1)

Definicja - zapas fazy

Niech M(ω), F(ω) będą charakterystykami częstotliwościowymi układu otwartego, a ω 1

częstotliwością, taką że M(ω1) =1. Zapasem fazy PM ( phase margin) nazywamy kąt określony wzorem

PM = 180o + F(ω )

1

Tok obliczeń:

G

( j

=

→

→ ∠

→

=

o + ∠

otw

ω ) 1

ω

G

( j

otw

ω )

PM

180

G

( j

otw

ω )

1

1

1

1

°

°

°

PM=180 -103.6 = 76.4

• Obliczenia „rę czne”

ω

( 9

− 0o

j

−arctg )

1

1

1

2

G

( s)

otw

=

,

G

(

otw jω) =

=

e

s( s + )

2

jω( jω + )

2

2

ω 4 + ω

ω : G (

otw

ω

j ) =1 → ω 4

2

+ω =1

4

→ ω + 4 2

ω −1= 0 → ω = 5 − 2 ≅ 4

.

0 858

1

1

0.4858

o

o

o

o

o

∠ G( jω ) = 9

− 0 − atan

= −103 6

.

→

PM = 180 − 103 6

.

= 76 4

.

1

2

Układ III rzędu

• Okreś lenie kkr

k

( s + )

1 3

k

G

( s)

zam

=

= 3

k

s + 3 2

s + 3 s + k + 1

1 + ( s + )13

Hurwitz: 3 ⋅ 3 > k + 1 → k < 8, kkr = 8

Wybieramy np. k = 4, aby układ był stabilny.

• Charakterystyki układu otwartego

Matlab

l = 4; m = [1 3 3 1];

w = logspace(-1,1,200);

[M, F] = bode(l,m,w);

subplot(211)

semilogx(w,M), grid

subplot(212)

semilogx(w,F), grid

[ w’

M

F ]

...........................

1.722

0.5062

-179.59

ω2

≈1/2

-180°

Definicja - zapas modułu

Niech ω 2 będzie częstotliwością, taką że F(ω2) = -180°. Zapasem modułu GM nazywamy odwrotność modułu M(ω2), tzn.

1

GM = M(ω )2

1

Powyżej M ≅ , zatem GM = 2.

2

• Obliczenia rę czne

4

4

j( 3 arctgω)

G

(

otw jω

−

) =

=

e

3

3

( jω + )

1

2 2

1

( + ω )

ω −

ω = −

o

→ ω =

2:

3arctg( )

180

3

2

G

( j

GM = 2

otw

ω2 =

1

)

= 4 = 4 = 1

→

3

3

8

2

1

( + )

3 2

4 2

Uwaga. Zapas modułu mówi, ile razy należy zwiększyć wzmocnienie, aby osiągnąć granicę stabilności, tzn. k ⋅ GM = k .

kr

Z Hurwitza: kkr = 8 → 4 ⋅ 2 = 8

W typowych układach sterowania zapas fazy PM wynosi 40°...75°, a zapas modułu GM

od 2 do 4. Im którykolwiek zapas mniejszy, tym większa skłonność do oscylacji.

Zapas modułu rozważa się tylko wtedy, gdy różnica stopni licznika i mianownika transmitancji układu otwartego wynosi przynajmniej 3, albo gdy w układzie występuje opóźnienie. Tylko wtedy charakterystyka fazowa schodzi poniżej -180°.

9.3.

Przeregulowanie a zapas fazy. Czas regulacji a częstotliwość

• Układ II rzę du - serwomechanizm 2

ω

G

( s)

n

=

zam

2

2

s + 2ξω s + ω

n

n

ω2

ω

n

G

( ω) =

,

∠

( ω = −

o −

otw j

o

G tw j )

90

arctg

2

2

ξ

2 ω n

ω ω + ( ξ

2 ω n )

• Czę stotliwość ω 1

4

ω n

4

2

G

( j

otw

ω ) = 1

→

= 1

→

ω = ω

4

n

ξ +1 − 2ξ

1

1

2

2

2

ω [ω + (2ξω ) ]

1

1

n

Matlab

ksi = 0:0.01:1;

KW=ksi. * ksi;

w1=sqrt(sqrt(4*KW. *KW+1)-2*KW);

clg

- ekran z jednym wykresem

plot(ksi,w1), grid

Wniosek. Im większe tłumienie ξ tym niższa częstotliwość ω 1, dla której układ osiąga zapas fazy.

• Zapas fazy PM

ω

ξ

4 4 + 1 − ξ 2

PM = 180o + ∠ G

(

o

o

otw jω

2

) = 180 − 90 − arctan

1

= ° −

1

ξ

2 ω

90

arctan

n

ξ

2

Przeregulowanie a zapas fazy

Matlab

p%

ξ

PM

PM=90-atan(w1. /(2*ksi))*180/pi;

37.2

0.3

33.3

plot(ksi, PM, ksi, 100*ksi), grid

25.4

0.4

43.1

16.3

0.5

51.8

9.5

0.6

59.2

4.6

0.7

65.1

4.3

1/ 2

65.5

0

1

76.3

Wniosek. „Inżynierski” wzór PM ≅ 100 ⋅ξ jest dobrym przybliżeniem pełnego wzoru na zapas fazy PM, ale tylko dla przeregulowań nie mniejszych niż 10% ( ξ ≤ 0.6 ) .

ξ

4 4 + 1 − ξ

2 2

PM = 90° − atan

≅ 100 ⋅ξ

ξ

2

• Czas regulacji tr a ω 1

t α

g

= x = x

→

arctg x = α

1

1

1

tg(90o − α ) =

→

90o − α = arctg

x

x

1

90o − arctg x = arctg x

ω

2ξω

Ponieważ

PM = 90o − arctg

1

, więc również PM

arctg

n

=

albo

2ξω

ω

n

1

2ξω

1

tg

n

PM =

ξω = ω

ω

oraz

PM

n

tg

1

.

1

2

4

4

Czas regulacji:

t =

=

r

ξω

1

n

ω tg PM

1

2

8

8

t =

ω =

r

ω

lub

1

tg PM

1

t

PM

r tg

• Idea projektowania

Mając dane p% i tr za pomocą powyższych wzorów przechodzimy na PM i ω1

w dziedzinie częstotliwości.

Podobnie mając charakterystyki częstotliwościowe układu otwartego w dziedzinie częstotliwości wyznaczamy z nich PM i ω1. Na podstawie powyższych wzorów PM

określamy ξ ≅

, przeregulowanie p% i czas regulacji tr.

100

-

p% a ξ :

p%

−πξ

ln

2

100

1− ξ

p

= e

100 %,

ξ =

%

2

p% 2

π + | ln

|

100

9.4.

Logarytmiczne charakterystyki Bodego

Idea metody

- Bode – 1943, USA

jF ω

-

( )

G( jω) = M (ω) e

,

M

= 20log M (

- wartości w decybelach [dB]

dB

ω)

- Wykresy MdB(ω) dla ω w skali logarytmicznej są liniami prostymi lub łamanymi o dB

dB

standardowych nachyleniach, tzn. ,

0 ± 20

, ± 40

itd.

dek

dek

- Transmitancja złożona

jF

j( F + F + F )

G = 1

G G 2 G 3 = Me

= M 1 M 2 M 3

1

2

3

e

20 log M = 20log M + 20log M + 20 log M

1

2

3

M

= M

+ M

+ M

dB

,

1 dB

2, dB

,

3 dB

F = F + F + F

1

2

3

Wniosek. Charakterystyki logarytmiczne transmitancji złożonej tworzy się jako sumę charakterystyk transmitancji składowych (prostych).

Charakterystyki elementarne

• Wzmacniacz

o

G( s) = k ,

j 0

G( jω) = k = k ⋅ e o

M dB = 20log k,

F (ω) = 0

• Całkowanie i róż niczkowanie

1

1

G( s) = , s,

, L

2

s

s

1

1

j( 90o

−

)

G( jω) =

= e

jω

ω

dB

M

= −

ω

→

−

db

20 log

20 dek

• Inercja i blok PD

1

G( s) =

,

Ts + 1

Ts + 1

1

j( arctg T

ω )

G( jω

−

) =

e

2

1 + ( T

ω )

1

ω < : M ≅ ,1

M dB = 0

T

1

1

ω > : M ≅

− jak integrator

T

T

ω

1

1

1

ω = : M =

→ 20log

≅ − d

3 B

T

2

2

o

F = − arctg 1 = 4

− 5

• Opóź nienie

s

τ

−

G( s) = e

j( ωτ )

G( jω

−

) = 1⋅ e

M dB = ,

0

F = ω

− τ [rad]

1

o

ω = :

F = 1

− rad = 5

− 7.3

τ

2

o

ω = :

F = −2 rad = 1

− 14.6

τ

Faza silnie maleje ze wzrostem częstotliwości.

• Blok II rzę du

2

ω

G( s)

n

=

- pkt. 9.1.

2

2

s + 2ξω s + ω

n

n

Skala logarytmiczna – „ręcznie”

• Oś ω

ω = 1..10 - dekada jako D jednostek (kratek, centymetrów itp.), gdzie D = 3,6,... –

wielokrotność 3. Odległość od początku dekady na osi w skali logarytmicznej wyraża wzór d = D logω.

Np. D = 3

ω

1

1.5

2

3

5

7

1

3logω

0

0,528

0.903

1.431

2.097

2.535

3

Odległość około:

0.5

1

1.5

2

2.5

Ręczne skalowanie

• Tangensoida

ω

0.1

0.2

1

5

10

arctgω

5.7°

11.3°

45°

78.7°

84.3°

Około

5°

10°

45°

80°

85°

Wniosek

Tangensoidę można ręcznie

wykreślać korzystając z

odcinków prostych.

Bode i Matlab

• Przykład 9.1. Transmitancja złoż ona 4 s + 2

G( s) =

0

.

0 02 s 3 + .

0 12 s 2 + s

Dobór zakresu czę stotliwoś ci następuje z zapasem nie przekraczającym dekady w prawo i w lewo od skrajnych częstotliwości granicznych (charakterystycznych). Częstotliwościami

tymi są rzeczywiste pierwiastki licznika i mianownika, a w przypadku pierwiastków zespolonych częstotliwość naturalna ω n.

Mianownik:

roots( [ 0.002 0.12 1 0 ] ) →

10,

50,

0

(

4 s + .

0 )

5

(

2 2 s + )

1

G( s) =

=

0

.

0 02 s( s + 1 )

0 ( s + 5 )

0

s( .

0 1 s + )

1 ( .

0 02 s + )

1

0.5, 10, 50 - częstotliwości charakterystyczne (odwrotności stałych czasowych) Zakres 0.1...100

- 3 dekady (150 punktów)

Matlab

l = [ 4 2 ]

m = [ 0.002 0.12 1 0 ]

w = logspace(-1,2,150);

[M, F] = bode(l,m,w);

subplot(211)

semilogx(w,M), grid

subplot(212)

semilogx(w,F), grid

decybele

dB = 20*log10(M)

semilogx(w,dB),grid

• Przykład 9.2. Obiekt z opóź nieniem 5

.

0 6 −0 5

. 2 s

G( s) =

e

2 s +1

3

2

1

G (

′ s)

∠ G( jω) = ∠ G (′ jω) − .

0

ω 180

52

π

1

1

Częstotliwości charakterystyczne:

= 0 5

. ,

= 1.92 → zakres 0.1 ... 10

2

0.52

Matlab

w = logspace(-1,1,100);

[M, Fprim] = bode(0.56,[2 1],w);

F=Fprim-0.52*w’*180/pi;

subplot(211)

semilogx(w,20*log10(M)), grid

subplot(212)

semilogx(w,F), grid

Uwaga. Charakterystyka fazowa obiektu z opóźnieniem osiąga silnie ujemne wartości.