Zatvaranje povratne sprege, primer

za početak, da učitamo toolbox

In [1]:
pkg load control

da prebacimo slike na vektorsku grafiku, ovo nije Octave komanda, jupyter je!

In [2]:
%plot inline -f "svg"

ovo je objekt upravljanja, njegova funkcija prenosa

In [3]:
T0 = tf([1],[1, 3, 3, 1])

Transfer function 'T0' from input 'u1' to output ...

                1          
 y1:  ---------------------
      s^3 + 3 s^2 + 3 s + 1

Continuous-time model.
In [4]:
bode(T0)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -120 -100 -80 -60 -40 -20 0 10-2 10-1 100 101 102 Magnitude [dB] Bode Diagram gnuplot_plot_1a -250 -200 -150 -100 -50 0 10-2 10-1 100 101 102 Phase [deg] Frequency [rad/s] T0 T0

funkcija prenosa ima malo pojačanje, nema čak ni crossover frequency (ako izuzmemo DC), stalno je manja od 1 po amplitudi; moramo da podignemo pojačanje i postavimo crossover frequency; zvuči jednostavno, čak i jeste, ali je problem gde postaviti crossover frequency? možda se zbog toga vratimo u drugom prolazu na ovo mesto, ako velik propusni opseg suviše zahteva od sistema (napona, struja, duty ratio van opsega $0 \le d \le 1$, . . . ); postavimo crossover frequency na $\omega_C = 10$ radijana u sekundi, $k = 1000$

In [5]:
T1 = 1000 * T0

Transfer function 'T1' from input 'u1' to output ...

              1000         
 y1:  ---------------------
      s^3 + 3 s^2 + 3 s + 1

Continuous-time model.
In [6]:
bode(T1)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -60 -40 -20 0 20 40 60 10-2 10-1 100 101 102 Magnitude [dB] Bode Diagram gnuplot_plot_1a -250 -200 -150 -100 -50 0 10-2 10-1 100 101 102 Phase [deg] Frequency [rad/s] T1 T1
In [7]:
nyquist(T1)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -500 0 500 -200 0 200 400 600 800 1000 Imaginary Axis Real Axis Nyquist Diagram T1 T1 gnuplot_plot_2a gnuplot_plot_3a

liči na nestabilno, ali čik razaznajte šta je oko (-1, 0) sa ove slike? detaljnije je dole prikazano, zoom je česta komanda kada se radi nyquist()

In [8]:
nyquist(T1)
axis([-1, 1, -1, 1])
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Imaginary Axis Real Axis Nyquist Diagram T1 T1 gnuplot_plot_2a gnuplot_plot_3a

vrlo nestabilno, ali treba da se tumači; da probamo sa mojom kompresiojom?

In [9]:
compressed(T1)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

krajnje očigledno nestabilno, obuhvata kritičnu tačku; da napravimo lead kompenzator da popravimo faznu marginu? ovde je glavno da kompenzator ne pokvari $\omega_C$, tako se pravi funkcija prenosa lead kompenzatora, da na $\omega_C$ ima pojačanje 1, odnosno 0 dB

In [10]:
Hlead = 0.1 * tf([1, 1], [0.01, 1])

Transfer function 'Hlead' from input 'u1' to output ...

      0.1 s + 0.1
 y1:  -----------
      0.01 s + 1 

Continuous-time model.
In [11]:
bode(Hlead)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -20 -10 0 10 20 10-1 100 101 102 103 Magnitude [dB] Bode Diagram gnuplot_plot_1a 20 40 60 80 10-1 100 101 102 103 Phase [deg] Frequency [rad/s] Hlead Hlead

na $\omega_C$ je amplituda lead kompenzatora 0 dB; neće uticati na $\omega_C$

In [12]:
T2 = T1 * Hlead

Transfer function 'T2' from input 'u1' to output ...

                      100 s + 100                
 y1:  -------------------------------------------
      0.01 s^4 + 1.03 s^3 + 3.03 s^2 + 3.01 s + 1

Continuous-time model.

ovde je popravljena faza; uočite kako red sistema raste! nažalost, to će nas pratiti do kraja

In [13]:
bode(T2)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -100 -80 -60 -40 -20 0 20 40 10-2 10-1 100 101 102 103 Magnitude [dB] Bode Diagram gnuplot_plot_1a -250 -200 -150 -100 -50 0 10-2 10-1 100 101 102 103 Phase [deg] Frequency [rad/s] T2 T2

i, šta kaže Nyquist? čik razaznajte odavde!

In [14]:
nyquist(T2)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -60 -40 -20 0 20 40 60 0 20 40 60 80 100 Imaginary Axis Real Axis Nyquist Diagram T2 T2 gnuplot_plot_2a gnuplot_plot_3a

jedan zoom . . .

In [15]:
nyquist(T2)
axis([-5 5 -5 5])
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -4 -2 0 2 4 -4 -2 0 2 4 Imaginary Axis Real Axis Nyquist Diagram T2 T2 gnuplot_plot_2a gnuplot_plot_3a

i još jedan . . .

In [16]:
nyquist(T2)
axis([-1 1 -1 1])
axis square
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Imaginary Axis Real Axis Nyquist Diagram T2 T2 gnuplot_plot_2a gnuplot_plot_3a

uz zoom se vidi da je jedva stabilno; da probamo kompresiju?

In [17]:
compressed(T2)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

lako i jasno se vidi da je stabilno, ali i mala margina faze; da probamo step response u zatvorenoj sprezi?

In [18]:
H2 = minreal(T2 / (1 + T2))

Transfer function 'H2' from input 'u1' to output ...

                            1e+04 s^4 + 4e+04 s^3 + 6e+04 s^2 + 4e+04 s + 1e+04                      
 y1:  -----------------------------------------------------------------------------------------------
      s^7 + 106 s^6 + 615 s^5 + 1.152e+04 s^4 + 4.201e+04 s^3 + 6.151e+04 s^2 + 4.06e+04 s + 1.01e+04

Continuous-time model.

sistem od četvrtog reda posta sistem sedmog reda zbog konačne dužine reči, tzv. numeričke greške! kasnije ćemo se time više baviti, sada samo da najavim:

In [19]:
H2a = unitnegativefeedback(T2)

Transfer function 'H2a' from input 'u1' to output ...

                      100 s + 100                 
 y1:  --------------------------------------------
      0.01 s^4 + 1.03 s^3 + 3.03 s^2 + 103 s + 101

Continuous-time model.

ovo je sistem četvrtog reda, uz malu pomoć symbolic precomputation, biće više reči kasnije

In [20]:
step(H2)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.5 1 1.5 0 5 10 15 20 y1 Time [s] Step Response H2 H2

uh, suvise oscilatorno! nema smisla sa ovim ići dalje, bolje da još popravimo marginu faze, još jedan lead kompenzator dodajemo na red!

a da vidimo onaj drugi sistem, četvrtog reda?

In [21]:
step(H2a, 20)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.5 1 1.5 0 5 10 15 20 y1 Time [s] Step Response H2a H2a

odziv je isti! naravoučenije: zaboravite na minreal ako baš ne mora; biće još reči o ovome!

In [22]:
T22 = T2 * Hlead

Transfer function 'T22' from input 'u1' to output ...

                          10 s^2 + 20 s + 10                    
 y1:  ----------------------------------------------------------
      0.0001 s^5 + 0.0203 s^4 + 1.06 s^3 + 3.06 s^2 + 3.02 s + 1

Continuous-time model.
In [23]:
bode(T22)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -80 -60 -40 -20 0 20 10-2 10-1 100 101 102 103 Magnitude [dB] Bode Diagram gnuplot_plot_1a -250 -200 -150 -100 -50 0 10-2 10-1 100 101 102 103 Phase [deg] Frequency [rad/s] T22 T22

dođosmo do sistema petog reda; šta kaže Nyquist?

In [24]:
nyquist(T22)
axis([-1 1 -1 1])
axis square
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Imaginary Axis Real Axis Nyquist Diagram T22 T22 gnuplot_plot_2a gnuplot_plot_3a

ok, stabilno, čak prevelika margina faze, možda bude presporo, ali rešićemo kasnije, bolje je za sada da bude stabilnije; da probamo kompresiju?

In [25]:
compressed(T22)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

vrlo jasno sta je, margina faze skoro $90^\circ$, nalik na sistem prvog reda; da probamo odziv na step pobudu u zatvorenoj sprezi?

In [26]:
H22 = minreal(T22 / (1 + T22))

Transfer function 'H22' from input 'u1' to output ...

                                1e+05 s^5 + 5e+05 s^4 + 1e+06 s^3 + 1e+06 s^2 + 5e+05 s + 1e+05                          
 y1:  -------------------------------------------------------------------------------------------------------------------
      s^8 + 206 s^7 + 1.122e+04 s^6 + 1.63e+05 s^5 + 6.54e+05 s^4 + 1.203e+06 s^3 + 1.151e+06 s^2 + 5.602e+05 s + 1.1e+05

Continuous-time model.

peti red sistema posta osmi, optimizacija ili skracivanje (minreal) nije baš sjajno odradila posao; ima i tu rešenja za tipične slučajeve, preko polinoma i delimičnog simboličkog računanja, to je bilo planirano za jedan master rad, ali je zavrseno za pola sata, pa je tema potrošena . . .

In [27]:
step(H22)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.2 0.4 0.6 0.8 0 1 2 3 4 5 6 y1 Time [s] Step Response H22 H22

očekivano, smirene oscilacije, usporen sistem; samo, vidite li da je pojačanje na niskim frekvencijama malo, imate grešku ustaljenog stanja! Da bi to popravili dodajemo lag kompenzator; on fazu kvari, pa ga postavljamo tako da je ne pokvari puno; $\omega_C$ isto ne sme da pokvari!

a da opet probamo alternativu?

In [28]:
H22a = unitnegativefeedback(T22)

Transfer function 'H22a' from input 'u1' to output ...

                           10 s^2 + 20 s + 10                      
 y1:  -------------------------------------------------------------
      0.0001 s^5 + 0.0203 s^4 + 1.06 s^3 + 13.06 s^2 + 23.02 s + 11

Continuous-time model.

ovo je sistem očekivanog petog reda; da mu probamo step?

In [29]:
step(H22a)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.2 0.4 0.6 0.8 0 1 2 3 4 5 6 y1 Time [s] Step Response H22a H22a

dakle, alternativni metod opet radi; da se vratimo na lag kompenzator

In [30]:
Hlag = tf([1, 1], [1 0])

Transfer function 'Hlag' from input 'u1' to output ...

      s + 1
 y1:  -----
        s  

Continuous-time model.
In [31]:
bode(Hlag)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 10 20 30 40 10-2 10-1 100 101 102 Magnitude [dB] Bode Diagram gnuplot_plot_1a -80 -60 -40 -20 0 10-2 10-1 100 101 102 Phase [deg] Frequency [rad/s] Hlag Hlag

na $\omega_C$ pojačanje je blisko 1, ne kvari se crossover frequency

In [32]:
T3 = T22 * Hlag

Transfer function 'T3' from input 'u1' to output ...

                      10 s^3 + 30 s^2 + 30 s + 10                 
 y1:  ------------------------------------------------------------
      0.0001 s^6 + 0.0203 s^5 + 1.06 s^4 + 3.06 s^3 + 3.02 s^2 + s

Continuous-time model.

eto sistema šestog reda!

In [33]:
bode(T3)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -80 -60 -40 -20 0 20 40 60 10-2 10-1 100 101 102 103 Magnitude [dB] Bode Diagram gnuplot_plot_1a -250 -200 -150 -100 10-2 10-1 100 101 102 103 Phase [deg] Frequency [rad/s] T3 T3

šta kaže Nyquist? razaznajte iz slike, ako možete . . .

In [34]:
nyquist(T3)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -300 -200 -100 0 100 200 300 -1 -0.8 -0.6 -0.4 -0.2 0 Imaginary Axis Real Axis Nyquist Diagram T3 T3 gnuplot_plot_2a gnuplot_plot_3a

da probamo kompresiju? bolje, ali ne sjajno, nema obilaznice, ovo je samo komprimovan dijagram koji daje nyquist() funkcija iz paketa control za octave; naša funkcija prenosa u otvorenoj sprezi sada ima pol u koordinatnom početku (nuli, (0, 0))

In [35]:
compressed(T3)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

a da probamo nešto bolje, pyquist funkciju koja ima i obilaznicu?

In [36]:
pyquist(T3)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

closed loop funkcija prenosa i njen step odziv?

In [37]:
H3 = minreal(T3 / (1 + T3))

Transfer function 'H3' from input 'u1' to output ...

                             1e+05 s^6 + 6e+05 s^5 + 1.5e+06 s^4 + 2e+06 s^3 + 1.5e+06 s^2 + 6e+05 s + 1e+05                        
 y1:  ------------------------------------------------------------------------------------------------------------------------------
      s^9 + 206 s^8 + 1.122e+04 s^7 + 1.63e+05 s^6 + 7.54e+05 s^5 + 1.703e+06 s^4 + 2.151e+06 s^3 + 1.56e+06 s^2 + 6.1e+05 s + 1e+05

Continuous-time model.

šesti red postao deveti, minreal() nije baš odradio posao, numerička greška radi svoje; sugestija za "developers": napraviti funkciju koja za T pravi H, funkciju prenosa u jediničnoj "negativnoj" povratnoj sprezi, da zaobiđemo numeriku; hm, da uradim to ja, kao što već jesam? pogledajte funkciju unitnegativefeedback.m, u direktorijumu vam je

In [38]:
H3p = unitnegativefeedback(T3)

Transfer function 'H3p' from input 'u1' to output ...

                           10 s^3 + 30 s^2 + 30 s + 10                      
 y1:  ----------------------------------------------------------------------
      0.0001 s^6 + 0.0203 s^5 + 1.06 s^4 + 13.06 s^3 + 33.02 s^2 + 31 s + 10

Continuous-time model.

sad šesti red osta šesti red! da im vidimo pzmaps?

In [39]:
[p, z] = pzmap(H3)

p =

  -127.95569 +   0.00000i
   -58.73944 +   0.00000i
   -13.30487 +   0.00000i
    -1.00434 +   0.00000i
    -1.00216 +   0.00376i
    -1.00216 -   0.00376i
    -0.99783 +   0.00373i
    -0.99783 -   0.00373i
    -0.99569 +   0.00000i

z =

  -1.00408 + 0.00235i
  -1.00408 - 0.00235i
  -1.00001 + 0.00471i
  -1.00001 - 0.00471i
  -0.99592 + 0.00236i
  -0.99592 - 0.00236i
In [40]:
[p, z] = pzmap(H3p)

p =

  -127.95569 +   0.00000i
   -58.73944 +   0.00000i
   -13.30487 +   0.00000i
    -1.00001 +   0.00001i
    -1.00001 -   0.00001i
    -0.99999 +   0.00000i

z =

  -1.00000 + 0.00001i
  -1.00000 - 0.00001i
  -1.00001 + 0.00000i

uporedite liste vidite šta radi numerička greška; zašto je treba izbeći? bypass je urađen preko analitičke pripreme sa simboličkim računom; samo par linija koda rešava naš partikularan problem jedinične povratne sprege; možda problem jeste partikularan, ali je naš, nas muči, to mu je značaj

In [41]:
step(H3)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 y1 Time [s] Step Response H3 H3

da uporedimo rezultat sa mojom funkcijom?

In [42]:
step(H3, H3p)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 y1 Time [s] Step Response H3 H3 H3p H3p

ok, to je to, odzivi su isti, funkcija radi; kao što vidite, sistem može lako da se unapredi, prilagodite ga svojim potrebama

step odziv je ok, nema premašenje, sve je očekivano i dobro; a da probamo da cedimo suvu drenovinu i još malo proširimo frekvencijski opseg velikih kružnih pojačanja, tj. da proširimo opseg frekvencija u kome je kružno pojačanje veliko?

In [43]:
Hlag2 = tf([1, 10], [1 0])

Transfer function 'Hlag2' from input 'u1' to output ...

      s + 10
 y1:  ------
        s   

Continuous-time model.
In [44]:
bode(Hlag2)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 10 20 30 40 10-1 100 101 102 103 Magnitude [dB] Bode Diagram gnuplot_plot_1a -80 -60 -40 -20 0 10-1 100 101 102 103 Phase [deg] Frequency [rad/s] Hlag2 Hlag2

na $\omega = 10$ nije baš pojačanje 0 db

In [45]:
T32 = T22 * Hlag2

Transfer function 'T32' from input 'u1' to output ...

                     10 s^3 + 120 s^2 + 210 s + 100               
 y1:  ------------------------------------------------------------
      0.0001 s^6 + 0.0203 s^5 + 1.06 s^4 + 3.06 s^3 + 3.02 s^2 + s

Continuous-time model.
In [46]:
bode(T32)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -50 0 50 10-2 10-1 100 101 102 103 Magnitude [dB] Bode Diagram gnuplot_plot_1a -250 -200 -150 -100 10-2 10-1 100 101 102 103 Phase [deg] Frequency [rad/s] T32 T32

crossover frequency je malo podignut; posledica preterivanja u $\omega_L$; šta kaže Nyquist? razaznajte, ovde to još i može

In [47]:
nyquist(T32)
axis([-1, 1, -1, 1])
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Imaginary Axis Real Axis Nyquist Diagram T32 T32 gnuplot_plot_2a gnuplot_plot_3a

da probamo kompresiju?

In [48]:
compressed(T32)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

opet problem sa polom u nuli; da uradimo obilaznicu?

In [49]:
pyquist(T32)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 -2 -1 0 1 2 -2 -1 0 1 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a

potpuno jasno, margina faze oko $45^\circ$; closed loop transfer function i jedinični odziv?

In [50]:
H32 = minreal(T32 / (1 + T32))

Transfer function 'H32' from input 'u1' to output ...

                            1e+05 s^6 + 1.5e+06 s^5 + 6e+06 s^4 + 1.1e+07 s^3 + 1.05e+07 s^2 + 5.1e+06 s + 1e+06                       
 y1:  ---------------------------------------------------------------------------------------------------------------------------------
      s^9 + 206 s^8 + 1.122e+04 s^7 + 1.63e+05 s^6 + 1.654e+06 s^5 + 6.203e+06 s^4 + 1.115e+07 s^3 + 1.056e+07 s^2 + 5.11e+06 s + 1e+06

Continuous-time model.

opet isti komentari za funkciju minreal() i red sistema; dižem ruke od toga zauvek, imam alternativu, svoj program!

In [51]:
H32 = unitnegativefeedback(T32)

Transfer function 'H32' from input 'u1' to output ...

                          10 s^3 + 120 s^2 + 210 s + 100                    
 y1:  ----------------------------------------------------------------------
      0.0001 s^6 + 0.0203 s^5 + 1.06 s^4 + 13.06 s^3 + 123 s^2 + 211 s + 100

Continuous-time model.

sada je sistem šestog reda, a isti je, tri reda manje! uz to, ovo rešenje je tačno, bez numeričke greške; tačnije i jednostavnije, toliko o izboru algoritma!

In [52]:
step(H32)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 1 2 3 4 5 6 y1 Time [s] Step Response H32 H32

očekivano, imamo brzinu, ali i ringing! da uporedimo odzive . . .

In [53]:
step(H3,H32,1.5)
Gnuplot Produced by GNUPLOT 5.2 patchlevel 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 y1 Time [s] Step Response H3 H3 H32 H32

"o odzivima ne valja raspravljati"; očigledno je oscilatoran nešto brži, ali da li se to isplati imajući u vidu ringing?