Stabilita zpětnovazební soustavy

©  Jiří Hospodka

Autorská práva: Uživatel může tento text používat jako studijní materiál bez omezení. Distribbuce a tisk je možný pouze se svolením autora.

>    restart;

Nejdříve načteme potřebné knihovny.

>    with(Syrup):with(plots):

Warning, the name changecoords has been redefined

Stabilitu budeme zkoumat na soustavě, obsahující zesilovač se třemi póly a odporovou zápornou zpětnou vazbu (ZZV) podle následujícího obrázku. Pro další analýzy budeme potřebovat jednak přenos celkového zapojení (se ZV) - Pzv  a jednak přenos otevřené zpětnovazevbní smyčky - Pos .

[Maple OLE 2.0 Object]

Zadání soustav pro analýzu

 První pól zesilovače omega1  je zadán obecně, další dva jsou na 2MHz a 4MHz. Zpětnovazební odpor Rp  je také zadán obecně (bez numerické hodnoty). Ostatní prvky jsou zadány s numerickými hodnotami.  

>    ZZ:= "zesilovač se ZZV
V 1 0
R01 1 2 10E3
Rp1 2 3
Rz1 3 0 10E3
Rp2 4 5
Rz2 4 0 10E3
R02 5 0 10E3
Ra2 5 0 100E3
X1 2 3 0 zesilovac
X2 1 4 0 zesilovac
 .subckt zesilovac a b c
 Ra a c 100E3
 Ri d b 100
 E1 d c a c -2E4/(1+s/omega1)/(1+s/(2*Pi*2E6))/(1+s/(2*Pi*4E6))
 .ends
.end":

Analýza systému pro omega1 = 2*Pi*3000  a Rp1 = 10000

Analýza a výpočet napěťového přenosu zesilovače se ZV a otevřené ZV smyčky.

>    napeti:=syrup(ZZ, ac):

>    Pzv:=simplify(subs(napeti,v[3]/v[1])):

>    Pos:=simplify(subs(napeti,v[5]/v[1])):

parsedeck:   Analyzing SPICE deck "zesilovač se ZZV" (ignoring this line)

Modulová charakteristika zesilovače pro omega1 = 2*Pi*.3e4  a Rp1 = .10e5 (zesilení přibližně 1). Ale je to opravdu modulová charakteristika, resp. je to stabilní soustava? Takto to nepoznám.  

>    modul_zv:=simplify(evalc(abs(subs({omega1=2*Pi*3E3,Rp1=10E3,s=I*omega},Pzv)))):

>    semilogplot(20*log10(modul_zv),omega=1E5..1E8,title=`MODUL`);

[Maple Plot]

Podle polohy pólů přenosu Pzv  to kmitá!!!

>    poly:=solve(denom(subs({omega1=2*Pi*3E3,Rp1=10E3},Pzv)));

poly := -52170908.01, 7226473.304-31911430.30*I, 7226473.304+31911430.30*I

>    plot

Vykreslení polohy pólů.

>    complexplot([poly],style=POINT,symbol=CROSS,symbolsize=30);

[Maple Plot]

Podle vykreslení přenosu Pos  to musí samozřejmě kmitat také.  

>    modul_os:=simplify(evalc(abs(subs({omega1=2*Pi*3E3,Rp2=10E3,s=I*omega},Pos)))):

>    faze_os:=simplify(evalc(argument(subs({omega1=2*Pi*3E3,Rp2=10E3,s=I*omega},Pos)))):

>    semilogplot(20*log10(modul_os),omega=1E7..1E8,title=`MODUL otevřené smyčky`);

>    semilogplot(faze_os*180/Pi,omega=1E7..1E8,title=`FAZE otevřené smyčky ve stupnich`);

[Maple Plot]

[Maple Plot]

Tady to je vidět lépe.

>    complexplot(subs({omega1=2*Pi*3E3,Rp2=10E3,s=I*omega},Pos),omega=1E7..1E8);

[Maple Plot]

Takto je potom přenos otevřené smyčky vykreslen ještě pro větší rozsah kmitočtů.

>    complexplot(subs({omega1=2*Pi*3E3,Rp2=10E3,s=I*omega},Pos),omega=1E2..1E7,numpoints=200);

[Maple Plot]

Kmitočet pro který je Im(Pos) = 0  lze spočítat jednoduše např. takto. Reálná složka přenosu pak vyjde větší než 4.  

>    om:=solve(evalc(Im(subs({omega1=2*Pi*3E3,Rp2=10E3,s=I*omega},Pos))));
evalf(evalc(Re(subs({omega1=2*Pi*3E3,Rp2=10E3,s=I*om[3]},Pos))));

om := 0., -17791513.49, 17791513.49

4.679896491

Soustava tedy kmitá a kmitočet oscilací lze určit z pólů přenosu Pzv . Pozor! Kmitočet oscilací se NEROVNÁ kmitočtu om , pro který je Im(Pos(om)) = 0  !! (Pozor! Výsledky nemusí být vždy stejně uspořádány.)  

>    om[3];
uhlovy_kmitocet_oscilaci:=Im(poly[3]);

17791513.49

uhlovy_kmitocet_oscilaci := 31911430.30

Skutečnost, že je systém nestabilní lze také jednoduše určit Laplaceovou transformací přenosu Pzv  (impulzní charakteristikou). Takto lze samozřejmě zjistit i kmitočet oscilací.  

>    with(inttrans);

[addtable, fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, invhilbert, invlaplace, invmellin, laplace, mellin, savetable]

Impulzní charakteristika obsahuje i Diracův impulz, protože výstupní odpor zesilovače není nulový (není nulvový přenos K3). Úhlový kmitočet oscilací se správně rovná imag. hodnotě kořenů jmenovatele přenosu Pzv  (vyšlo).  

>    odezva:=simplify(invlaplace(subs({omega1=2*Pi*3E3,Rp1=10E3},Pzv),s,t));

odezva := .4690431520e-2*Dirac(t)-12341256.64*exp(-52170908.01*t)+12341256.64*exp(7226473.304*t)*cos(31911430.30*t)-22971026.98*exp(7226473.304*t)*sin(31911430.30*t)
odezva := .4690431520e-2*Dirac(t)-12341256.64*exp(-52170908.01*t)+12341256.64*exp(7226473.304*t)*cos(31911430.30*t)-22971026.98*exp(7226473.304*t)*sin(31911430.30*t)

Že násobná konstanta je skutečně koeficientem K3 se lze přesvědčit jednoduchým výpočtem.

>    K3:=Ra/(R0+Ra)*Ri*Rz/(Ri+Rz)/(R0*Ra/(R0+Ra)+Rp+Ri*Rz/(Ri+Rz));
evalf(subs({Ra=100E3,R0=10E3,Rp=10E3,Ri=100,Rz=10E3},%));

K3 := Ra/(R0+Ra)*Ri*Rz/(Ri+Rz)/(R0*Ra/(R0+Ra)+Rp+Ri*Rz/(Ri+Rz))

.4690431520e-2

Takto pak vypadá impulzní charakteristika. Díky poloze pólů se amplituda kmitů musí stále zvětšovat (pro lineární systém).  

>    plot(odezva,t=1E-10..1E-6, title=`Impulzni charakteristika`);

[Maple Plot]

Soustava na mezi stability

Nyní vyřešíme minimální (kritickou) hodnotu zpětnovazebního odporu Rp_krit  (pro minimální zesílení systému), pro který systém ještě nekmitá, resp. je na mezi stability.  

>    rov1:=evalc(Re(subs({omega1=2*Pi*3E3,s=I*omega},Pos)))=1:
rov2:=evalc(Im(subs({omega1=2*Pi*3E3,s=I*omega},Pos)))=0:

>    reseni1:=solve({rov1,rov2},{Rp2,omega});

reseni1 := {omega = 0., Rp2 = -180027191.7}, {Rp2 = 80616.91549, omega = -17791513.49}, {omega = -.5501824491e13*I, Rp2 = -9189.918992}, {Rp2 = -9189.918992, omega = .5501824491e13*I}, {Rp2 = 80616.915...
reseni1 := {omega = 0., Rp2 = -180027191.7}, {Rp2 = 80616.91549, omega = -17791513.49}, {omega = -.5501824491e13*I, Rp2 = -9189.918992}, {Rp2 = -9189.918992, omega = .5501824491e13*I}, {Rp2 = 80616.915...
reseni1 := {omega = 0., Rp2 = -180027191.7}, {Rp2 = 80616.91549, omega = -17791513.49}, {omega = -.5501824491e13*I, Rp2 = -9189.918992}, {Rp2 = -9189.918992, omega = .5501824491e13*I}, {Rp2 = 80616.915...

Reálná je tato hodnota. (Pozor! Výsledky nemusí být opět vždy stejně uspořádány.)

>    Rp_krit:=subs(reseni1[2],Rp2);

Rp_krit := 80616.91549

Takto vyjdou póly přenosu Pzv , pro kritickou hodnotu zpětnovazebního odporu (systém na mezi stability), samozřejmě je zde konečná přesnost.

>    reseni2:=solve(denom(subs({omega1=2*Pi*3E3,Rp1=Rp_krit},Pzv)));

reseni2 := -37717961.40, .1010270550e-3-17791513.49*I, .1010270550e-3+17791513.49*I

Právě a pouze v tomto případě (mez stability) se rovná kmitočet, pro který platí Im(Pos) = 0   kmitočtu oscilací!! (porovnej s imaginární hodnotou pólu přenosu - reseni2 ). (Pozor! Výsledky nemusí být opět vždy stejně uspořádány.)  

>    reseni1[5][2];

>    Im(reseni2[3]);

omega = 17791513.49

17791513.49

Přenos pro Rp = Rp_krit . Na modulové charakteristice to je teď znát. Patrné to musí být i z přenosu otevřené smyčky (nyní vynesené jen v komplexní rovině).

>    modul_krit:=simplify(evalc(abs(subs({omega1=2*Pi*3E3,Rp1=Rp_krit,s=I*omega},Pzv)))):

>    semilogplot(20*log10(modul_krit),omega=1E5..1E8,title=`MODUL`);
complexplot(subs({omega1=2*Pi*3E3,Rp2=Rp_krit,s=I*omega},Pos),omega=1E7..1E8);

[Maple Plot]

[Maple Plot]

Odvodili jsme, hodnotu Rp1 = Rp_krit  pro systém, který je právě na mezi stability. "Zesílení" takovéhoto systému pro stejnosměrný signál je potom:

>    Pzv0:=evalf(subs({Rp1=Rp_krit,s=0},Pzv));

Pzv0 := -8.057671260

Z předchozích výpočtů je zřejmé, že pro stabilní systém musí platit Rp_krit < Rp1 , resp. Pzv0 < Pzv . Pokusme se nyní kompenzovat zesilovač tak aby byl stabilní i pro zesílení Pzv = -1  (ve SKP).   

Závislost polohy pólů zesilovače na vlivu ZV  

Výpočet pólů přenosu ZV soustavy pro různé hodnoty odporu Rp1  

>    poly1:=seq([solve(denom(subs({omega1=2*Pi*3E3,Rp1=R0},Pzv)),s)],R0=[100,Rp_krit,1e6,1.5e6,10e6]);

poly1 := [-62332415.05, 12307226.83-41223154.43*I, 12307226.83+41223154.43*I], [-37717961.40, .1010270550e-3-17791513.49*I, .1010270550e-3+17791513.49*I], [-27674337.34, -5021812.029-3656325.475*I, -50...
poly1 := [-62332415.05, 12307226.83-41223154.43*I, 12307226.83+41223154.43*I], [-37717961.40, .1010270550e-3-17791513.49*I, .1010270550e-3+17791513.49*I], [-27674337.34, -5021812.029-3656325.475*I, -50...
poly1 := [-62332415.05, 12307226.83-41223154.43*I, 12307226.83+41223154.43*I], [-37717961.40, .1010270550e-3-17791513.49*I, .1010270550e-3+17791513.49*I], [-27674337.34, -5021812.029-3656325.475*I, -50...
poly1 := [-62332415.05, 12307226.83-41223154.43*I, 12307226.83+41223154.43*I], [-37717961.40, .1010270550e-3-17791513.49*I, .1010270550e-3+17791513.49*I], [-27674337.34, -5021812.029-3656325.475*I, -50...

>    barva:=[violet,red,magenta,green,orange];

barva := [violet, red, magenta, green, orange]

Je nutné sledovat póly podle barev. Stejné barvy = stejná hodnota Rp , resp. vratného rozdílu F .

>    for i from 1 to 5 do
graf[i]:=complexplot([poly1[i][]], style=point, thickness=20, title=`Poloha pólů přenosu ZV soustavy v závislosti na Rp`,xtickmarks=0,ytickmarks=0,color=barva[i],symbolsize=20);
od:
display(graf[1],graf[2],graf[3],graf[4],graf[5]);

[Maple Plot]

Lépe to bude vidět na následující animaci.

Nejprve si připravíme sekvenci čísel - jednotlivé hodnoty odporu Rp  (to je nutné udělat nelineárně, abyby potom póly byly rozloženy téměř rovnoměrně v clém rozsahu).

>    Rp_seq:=evalf(seq(8^((R0+18)/4),R0=1..15)):

Výpočet pólů a vykreslení jednotlivých obrázků (zde obrátíme pořadí, abychom postupovali od reálných pólů k pólům komplexně združeným).

>    vv:=seq([fsolve(denom(subs({omega1=2*Pi*3E3,Rp1=Rp_seq[nops([Rp_seq])+1-i]},Pzv)),s,complex)],i=1..nops([Rp_seq])):

>    obrazky:=seq([op(1,complexplot(vv[i]))],i=1..nops([Rp_seq])):

Vykreslení (v Maple klikněte na obrázek a potom vyberte volbu "přehrávat" sekvenci dokola).

>    PLOT(ANIMATE(obrazky), AXESTICKS(0,0), TITLE(`Poloha pólů přenosu ZV soustavy v závislosti na Rp`), THICKNESS(20), AXESLABELS("",""), SYMBOL(DEFAULT,20), STYLE(POINT), VIEW(DEFAULT,DEFAULT));

[Maple Plot]

Kompenzace zesilovače pro stabilni provoz i při zesílení -1

Aby měl zesilovač zesílení rovno -1, je nutné zapojit do ZV odpor Rp1_1 . (Pozor! Výsledky nemusí být opět vždy stejně uspořádány.)

>    Rp1_1:=solve(simplify(subs(s=0,Pzv))=-1);

Rp1_1 := 10001.07106

Aby zesilovač při tomto zesílení nekmital, a měl naopak fázovou jistotu např. 45^o , je nutné změnit kmitočet prvního zlomu omega1  na hodnotu omega1_1  (tj. výše uvedená kompenzace). (Pozor! Výsledky nemusí být opět vždy stejně uspořádány.)

>    rov1:=evalc(abs(subs({Rp2=Rp1_1,s=I*omega},Pos)))=1:
rov2:=evalc(argument(subs({Rp2=Rp1_1,s=I*omega},Pos)))=Pi/4:
fsolve({rov1,rov2},{omega1,omega},omega=1E6..1E9);omega1_1:=subs(%,omega1);

{omega1 = 896.3634708, omega = 7057984.966}

omega1_1 := 896.3634708

Že to Maple spočítal dobře se můžeme ještě přesvědčit z průběhu modulové a fázové charakteristiky přenosu otevřené ZV smyčky.

>    modul_os_45:=simplify(evalc(abs(subs({omega1=omega1_1,Rp2=Rp1_1,s=I*omega},Pos)))):

>    faze_os_45:=simplify(evalc(argument(subs({omega1=omega1_1,Rp2=Rp1_1,s=I*omega},Pos)))):

>    semilogplot(20*log10(modul_os_45),omega=1E6..1E8,title=`MODUL otevřené smyčky`);

>    semilogplot(faze_os_45*180/Pi,omega=1E6..1E8,title=`FAZE otevřené smyčky ve stupnich`);

[Maple Plot]

[Maple Plot]

Mužeme si to ještě ověřit numerickým výpočtem. (Pozor! Výsledky nemusí být opět vždy stejně uspořádány.)

>    om_45:=solve(modul_os_45=1,omega);

>    evalf(subs(omega=om_45[6],faze_os_45)*180/Pi);

om_45 := -24609944.80*I, 24609944.80*I, -15288365.56*I, 15288365.56*I, -7057984.966, 7057984.966

45.00000000

Není to elegantní a přitom jednoduché?

>