Modelarea Si Proiectarea Unui Regulator Digital Ptr Sistemele de Suspensie ale Automobilelor

Cap.I Introducere

I.1 Introducere

I.2 Modelarea sistemului de suspensie activă utilizând spațiul stărilor

I.2.1 Cerințe de proiectare:

I.2.2 Reprezentarea sistemului utilizând spațiul stărilor

I.2.3 Reprezentarea sistemului în spațiul stărilor utilizând Matlab

I.2.4 Răspunsul sistemului în buclă deschisă

Cap.II Modelarea sistemului utilizând funcția de transfer

II.1 Conversia modelului matematic din spațiul stărilor în funcție de transfer

II.2 Modelarea sistemului de suspensie activă utilizând funcția de transfer

II.2.1 Funcția de Transfer:

II.2.2 Reprezentarea sistemului prin funcții de transfer utilizând Matlab:

II.2.3 Răspunsul sistemului în buclă deschisă (fără control în buclă)

Cap.III Sinteza în domeniul frecvența a unui regulator PID utilizând metoda locului rădăcinilor

III.1 Considerații teoretice

III.2 Trasarea locului rădăcinilor

III.3 Adăugarea unui filtru notch ( filtru stop bandă, filtru crestătură)

III.4 Trasarea răspunsului sistemului închis

Cap.IV Sinteza în frecvență a regulatorului utilizând caracteristicile Bode

IV.1 Trasarea răspunsului în frecvență utilizând comanda bode

IV.2 Utilizarea unui regulator cu dublu avans de fază

IV.2.1 Considerații teoretice

IV.2.2 Compensatorul cu avans de fază proiectat utilizând locul rădăcinilor

IV.2.3 Compensatorul cu avans de fază proiectat prin utilizarea răspunsului în frecvență

Cap.V Sinteza unui regulator optimal (LQR) pentru SSA

V.1 Considerații teoretice

V.2 Rezolvarea problemei de control optimal în Matlab

V.3 Regulator LQ pentru sistemul de suspensie activă

V.4 Trasarea răspunsului sistemului în buclă închisă

Cap.VI Proiectarea unui regulator digital cu reacție de la stare

VI.1 Metoda de reproiectare în timp discret

VI.1.1 Definirea cerințelor de proiectare

VI.1.2 Selectarea perioadei de eșantionare

VI.1.3 Conversia continuu-discret

VI.1.4 Adăugarea unui integrator(pentru eroare staționară nulă)

VI.1.5 Proiectarea controlerului

VI.1.6 Simularea răspunsului sistemului închis 

Cap.VII Controlul robust al suspensiei active

VII.1 Modelul de stare al suspensiei active pentru un sfert de mașină

VII.2 Modelul suspensiei active, pentru proiectarea controlerului robust liniar H∞

VII.3 Proiectarea controlerului H∞ liniar cu funcția hinfsyn, pentru a minimiza funcția de transfer a abaterii corpului mașinii(x1)

VII.4 Proiectarea controlerului H∞ liniar cu funcția hinfsyn pentru a minimiza funcția de transfer a abaterii suspensiei(x1-x3)

VII.5 Analiza în domeniul timp

VII.6 Proiectarea controlerului robust Kdk(K3) aferent suspensiei active a autovehiculelor utilizând toolbox-ul μ(iterația D-K)

VII.7 Analiza în domeniul timp utilizând K1 și Kdk

Cap.VIII Concluzii

Bibliografie

ANEXĂ

77 pagini

=== Proiect ===

Cap.I Introducere

I.1 Introducere

Performanța sistemului de suspensie auto este tipic stabilită de capacitatea acestuia de a asigura conducerea îmbunătățită pe șosea și de confortul sporit al pasagerilor .

Sistemele de suspensie ale automobilelor actuale utilizând componente pasive pot oferi doar un compromis al acestor două criterii contradictorii prin prevederea coeficienților de elasticitate și de amortizare cu măsuri stabilite.

De obicei, automobilele de curse au suspensii solide care conferă un confort redus călătorilor, în vreme ce berlinele luxoase prezintă suspensii estompate dar cu capacități nesatisfăcătoare de conducere pe șosea. Randamentul redus de conducere pe șosea și confortul scăzut al pasagerilor sunt datorate fără măsură vibrațiilor corpului automobilului.

Sistemele de suspensie activă urmăresc să amelioreze aceste efecte nedorite prin izolarea corpului vehiculului de vibrațiile roții determinate prin traversarea unui teren cu asperități (gropi, ridicături, etc).

Numeroasele structuri și implementări ale sistemelor de suspensie a automobilelor pot fi clasificate ca suspensii pasive, suspensii semi-active, sau ca suspensii active.

Încercările de îmbunătățire a performanțelor utilizând suspensii pasive necesită adăugarea unor arcuri suplimentare în afara roților, sau reproiectarea completă a sistemului de suspensie mecanic.

Obiectivul principal al sistemelor de suspensie activă este acela de a reduce mișcările masei corpului autovehiculului. Este de asemenea cunoscut faptul că, mișcările corpului autovehiculului, la caracteristicile de frecvență ale roții, nu pot fi reduse dacă singura intrare de control este doar o forță aplicată între masa corpului vehiculului și masa fără arcuri (acesta este cazul pentru sistemele de suspensie auto).

I.2 Modelarea sistemului de suspensie activă utilizând spațiul stărilor

Proiectarea unui sistem de suspensie activă pentru autovehicule se transformă într-o interesantă problemă de control .

Sistemul de suspensie activă este utilizat pentru simplificarea problemei sistemului de amortizare a vehiculului. Schema simplificată a unui astfel de sistem este prezentat în continuare în figura 1,

Figura I.1 – 1/4 din modelul suspensiei active a unui automobil (o roată)

unde, de exemplu, in cazul unui microbuz, avem:

M = m1 = 290 kg – masa corpului autovehiculului,

m = m2 = 59 kg – masa suspensiei(dar fără arcuri și amortizor),

k1 = 16182N/m – constanta elastică a arcului sistemului de suspensie,

k2 =190000 N/m – constanta elastică a roții pneumatice,

d1 = 1000 Ns/m – constanta de amortizarea a amortizorului sistemului de suspensie,

d2 = 0 Ns/m – constanta de amortizarea a roții și a jantei,

u – forța[kN] (mărimea) de comandă/control = forța de suspensie activă (mărime de xxxxxxxxxxintrare),

x1 – variația corpului vehiculului (mărime de ieșire = mărime reglată/controlată),

x2 – variația roții de mers pe drum (mărime de intrare = mărime necontrolată),

w – perturbația (înălțimea căii de rulare, deci tot o mărime de intrare necontrolată).

I.2.1 Cerințe de proiectare:

Un sistem de suspensie bun trebuie să aibă calități de confort, nu numai când automobilul circulă pe o șosea dreaptă ci și atunci când autovehiculul circulă peste ridicăturile și gropile șoselei. Când autovehiculul trece peste orice denivelare (gaură, ridicătură) a șoselei, corpul acestuia nu trebuie să aibă oscilații largi, iar oscilațiile trebuie să dispară imediat.

Întrucât distanța (X1-W) este dificil de măsurat, iar deformația roții și a jantei (X2-W) este neglijabilă, vom utiliza distanța (X1-X2) în loc de (X1-W), ca mărime de ieșire (care este doar o estimare) în problema noastră.

Perturbația șoselei notată cu W o vom simula cu o treaptă pe intrare. Cerințele de proiectare se referă la dorința de a proiecta un regulator (controller) într-o buclă cu reacție negativă, astfel încât ieșirea (X1-X2) să aibă o suprareglare σ cat mai mică și un timp tranzitoriu 3 secunde.

De exemplu, când autovehiculul va rula peste o treaptă înaltă de 10 cm, corpul vehiculului va trebui să oscileze în banda de 5 mm și să revină la o rulare normală fără oscilații, într-un timp mai mic de 3.5 secunde(adică σ=5% iar tt=3.5 sec.).

I.2.2 Reprezentarea sistemului utilizând spațiul stărilor

Reprezentarea sistemului în spațiul stărilor se poate obține din ecuațiile diferențiale ce caracterizează mișcarea automobilului obținute ca urmare a aplicării legilor de conservare a energiei, ale lui Newton și a lui Hooke.

Ecuațiile dinamicii autobuzului și masa suspensiei sunt:

(1.2)

Pentru a fi o reprezentare valabilă a spațiului stărilor, derivatele tuturor stărilor trebuie să fie funcții de intrări si de stări insele. Pentru început se împart prima și a doua ecuație prin M, respectiv m. De notat că apare în ecuație pentru .

(1.3)

(1.4)

Prima stare va fi X1. Din moment ce derivatele intrărilor nu apar în ecuația pentru, se alege a doua stare ca fiind . A treia stare se alege ca fiind diferența dintre X1 și X2 și apoi se determină a patra stare. Astfel se înlocuiește cea de-a treia stare cu Y1=X1-X2 în ecuația de mai jos:

(1.5)

(1.6)

Scăzând a doua ecuație din prima, se obține expresia pentru :

(1.7)

Deoarece nu putem utiliza derivatele de ordinul doi pentru reprezentarea în spațiul stărilor, se va integra această ecuație pentru a obține :

(1.8)

Notăm integrala cu Y2. Ecuația de stare pentru Y2 este:

(1.9)

(1.10)

Se înlocuiește X2 = X1-Y1 în pentru a obține ecuația de stare pentru Y1:

(1.11)

Apoi se înlocuiește derivata lui Y1 în ecuația derivatei lui X1 :

(1.12)

Variabilele de stare sunt X1,, Y1 și Y2. Forma matriceală a ecuațiilor prezentate mai sus este:

(1.13)

(1.14)

Pentru marea majoritate a problemelor, metoda ce utilizează modelul în spațiului stărilor este mult mai ușoară decât metoda transformatei Laplace.

I.2.3 Reprezentarea sistemului în spațiul stărilor utilizând Matlab

Se pot pune ecuațiile de mai sus într-un fișier .m în Matlab prin definirea a patru matrice A,B,C, D ale ecuațiilor de stare standard:

(1.15)

Se creează fișierul .m cu următoarele linii de comandă care se află în Anexă :

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

I.2.4 Răspunsul sistemului în buclă deschisă

Poate fi utilizat Matlab pentru a vedea cum evoluează sistemul original în buclă deschisă. Se adaugă următoarele comenzi în fișierul .m și apoi se rulează întreg programul în fereastra de lucru a Matlab, pentru a observa răspunsul sistemului când la intrare este aplicată o treaptă unitate:

step(A, B, C, D, 1)

Figura I.2 – Răspunsul sistemului în buclă deschisă la intrare treaptă unitate

Din acest grafic, se poate observa că răspunsul sistemului în buclă deschisă la intrare treaptă este neamortizat deci călătorii autovehiculului vor simți oscilațiile . Mai mult, automobilul necesită un timp tranzitoriu () destul de mare pentru a atinge eroarea de regim staționar (). Soluția pentru această problemă este aceea de a adăuga un controller într-o buclă cu reacție, în schema bloc a sistemului.

Răspunsul sistemului în buclă deschisă la perturbația “w” poate fi vizualizat prin comanda:

step(A, 0.1*B, C, D, 2),

adăugată fișierului .m de la paginile 5 și 6.

Figura I.3 – Răspunsul sistemului în buclă deschisă la perturbație treaptă ”w”

Matricea B a fost multiplicată cu valoarea 0.1 pentru a simula treapta de 10 cm, care reprezintă perturbația “w” aplicată sistemului.

Din acest grafic se observă că, atunci când automobilul trece peste o denivelare (ridicătură, groapă) de pe șosea de 10 cm, corpul automobilului va oscila pentru un timp tranzitoriu () lung. Suprareglarea mare (de la impactul inițial) și timpul tranzitoriu mare, vor provoca avarierea sistemului de suspensie. Soluția acestei probleme este utilizarea unui sistem automat care să crească performanțele sistemului.

Schema sistemului în buclă închisă este următoarea:

w

Figura I.4 – Schema sistemului în buclă închisă

Cap.II Modelarea sistemului utilizând funcția de transfer

II.1 Conversia modelului matematic din spațiul stărilor în funcție de transfer

Înainte de a începe proiectarea regulatorului, este necesară găsirea funcțiilor de transfer pentru proces, respectiv pentru perturbație. Funcția de transfer a procesului poate fi găsită (de la intrarea u, la ieșirea (X1-X2) prin utilizarea în Matlab a comenzii ss2tf.

Se adaugă următoarele comenzi în fișierul .m creat mai sus:

[nump,denp] = ss2tf (A, 0.1*B, C, D, 1)

Valoarea 1 arată care intrare trebuie să fie utilizată pentru calculul funcției de transfer. Prin rularea fișierului .m în fereastra de lucru Matlab, se obțin următoarele:

nump = [0 0.0000 0.0004 0.0019 0.0625]

denp = 1.0e+004 *[0.0001 0.0048 0.1851 0.1719 5.0000]

De notat faptul că numărătorul este un polinom de gradul 5 cu două zerouri. Acest lucru ar putea cauza câteva probleme în viitor, când se va adăuga un regulator sistemului. Când numărătorul este combinat cu numărătorul regulatorului, cele două zerouri dominante pot crește ordinul polinomului numărătorului al buclei închise astfel încât acesta este de ordin mai mare față de polinomul numitorului. Această situație nu este permisă în Matlab, el returnând un mesaj de eroare. Pentru a elimina aceste două cifre, se redefinește num astfel încât să fie format doar din ultimele trei cifre prin adăugarea următoarei linii în fișierul .m :

nump=[nump(3) nump(4) nump(5)]

Prin urmare funcția de transfer a procesului devine:

P(s) =

Acum este cunoscută funcția de transfer a procesului numpr/denpr, dar nu cunoaștem care este funcția de transfer a lui F. Mai întâi, este necesară aflarea funcției de transfer de la intrarea W la ieșirea X1-X2 prin utilizarea comenzii ss2tf:

[ num1, den1] = ss2tf ( A, 0.1*B, C, D, 2)

De notat faptul că valoarea 2 este utilizată pentru aflarea funcției de transfer de la cea de-a doua intrare (perturbația), iar valoarea 0.1 corespunde intrării treaptă de 10 cm.

Matlab va returna:

num1 =[ 0 -55100000 0 0]

den1 = 1.0e+009 *[ 0.0000 0.0003 0.0607 0.1900 3.0746]

Se adaugă următoare linie în fișierul .m:

num1 = [ num1(2) num1(3) num1(4) num1(5)]

Din schema prezentată mai sus se observă că:

F*Proces = ,

astfel încât

F =,

sau

.

Funcția de transfer pentru controllerul buclei cu reacție numc/denc, este singurul lucru care a mai rămas de aflat pentru a satisface specificațiile cerute sistemului automobilului.

Deci, schema bloc a sistemului automobilului devine:

w

Figura II.1 – Schema bloc a sistemului automobilului în buclă închisă

II.2 Modelarea sistemului de suspensie activă utilizând funcția de transfer

II.2.1 Funcția de Transfer:

Ecuațiile dinamicii vehiculului sunt:

(2.1)

(2.2)

Presupunem că toate condițiile inițiale sunt nule, așa că aceste ecuații reprezintă situația când roata vehiculului întâlnește o denivelare. Aplicăm transformata Laplace în cazul celor două ecuații de mai sus și găsim funcțiile de transfer G1(s) și G2(s):

(2.3)

(2.4)

= (2.5)

A = (2.6)

Δ = det (2.7)

sau

Δ = . (2.8)

Trebuie determinată inversa matricei A și apoi înmulțită cu intrările U(s) și W(s) astfel:

(2.9)

(2.10)

Când se consideră doar intrarea U(s), se setează W(s)=0. Astfel se obține funcția de transfer G1(s):

(2.11)

Când se consideră doar intrarea perturbatoare W(s), se setează U(s) = 0. Funcția de transfer G2(s) este:

(2.12)

II.2.2 Reprezentarea sistemului prin funcții de transfer utilizând Matlab:

Cele două funcții de transfer, G1(s) și G2(s) mai sus determinate, pot fi introduse în Matlab prin definirea numitorului și a numărătorului funcțiilor de transfer în forma, numpr/denpr pentru forța de intrare de la elementul de execuție (control) și num1/den1 pentru intrarea perturbație :

G1(s) = numpr/denpr

G2(s) = num1/den1

Se creează un nou fișier .m care se află în Anexă :

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0; -(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) – (d1/m1); d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1;k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0]

B = [0;1/m1;0;(1/m1)+(1/m2)]

C = [0 0 1 0];

D = [0 0];

%Determinarea functiei de transfer a procesului

nump=[(m1+m2) d2 k2];

denp=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2];

'G1(s)'

printsys (nump, denp)

%Răspunsul sistemului în buclă deschisă la intrare treaptă

step (nump, denp)

title ('Răspunsul sistemului în buclă deschisă la intrare treaptă')

%Determinarea funcției de transfer pentru perturbație

num1=[-(m1*d2) -(m1*k2) 0 0];

den1=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2];

'G2(s)'

printsys (0.1*num1,den1).

II.2.3 Răspunsul sistemului în buclă deschisă (fără control în buclă)

Se adaugă următoarele comenzi în fișierul .m și se rulează programul în fereastra de comandă a Matlab, pentru a vedea răspunsul sistemului la intrare treaptă și la intrare perturbație “w”:

step (nump, denp)

Figura II.2 – Răspunsul sistemului în buclă închisă la intrare treaptă

Din acest grafic se observă că sistemul deschis are un răspuns neamortizat.

Timpul tranzitoriu () este mai mare decat cel dorit pentru a atinge eroarea de regim staționar (). Soluția pentru această problemă este de a adăuga un regulator în bucla cu reacție.

Răspunsul sistemului în buclă deschisă la perturbația “w” poate fi vizualizat prin comanda:

step (0.1*num1,den1)

Figura II.3 – Răspunsul sistemului în buclă deschisă la perturbația treaptă ”w”

Pentru a vedea câteva detalii, se pot schimba axele prin comanda axis:

axis ([0 10 –0.1 0.1])

Figura II.4 – Răspunsul sistemului în buclă deschisă la perturbația treaptă ”w”

Din acest grafic al răspunsului sistemului deschis pentru o treaptă de 0.1m pe perturbație, se observă că, atunci când automobilul trece peste o denivelare(ridicătură) pe șosea de 10 cm(0.1 m), corpul automobilului va oscila un timp tranzitoriu () având valoarea de 4 sec. Suprareglarea mare și timpul tranzitoriu mare, pot provoca avarierea sistemului de suspensie. Soluția acestei probleme este utilizarea unui sistem automat care să ridice performanțele sistemului.

Schema sistemului în buclă închisă este următoarea:

w

Figura II.5 – Schema sistemului în buclă închisă

Din funcțiile de transfer și schema de mai sus, putem desena diagrama bloc a sistemului autovehiculului astfel:

w

Figura II.6 – Schema sistemului în buclă închisă

Din schema de mai sus observăm că:

Proces =

F * Proces =

astfel încât:

F =

F = * = =

Cap.III Sinteza în domeniul frecvența a unui regulator PID utilizând metoda locului rădăcinilor

III.1 Considerații teoretice

Metoda locului rădăcinilor permite ca, pornind de la cunoașterea repartiției polilor și zerourilor funcției de transfer a sistemului în circuit deschis în planul s complex, să tragem concluzii asupra rădăcinilor ecuației caracteristice a sistemului de reglare automată în circuit închis . Se variază, de exemplu, un parametru al sistemului deschis, astfel încât se modifică poziția în planul s a rădăcinilor ecuației caracteristice a sistemului de reglare în circuit închis. Rădăcinile descriu în acest caz în planul s, curbe ce definesc locul rădăcinilor ecuației caracteristice a sistemului de reglare în circuit închis.

Cunoașterea locului rădăcinilor, care este reprezentat de cele mai multa ori în funcție de un parametru, ne permite să facem aprecieri asupra stabilității sistemului închis și asupra criteriului de stabilitate, de exemplu asupra distanței polilor de axa reală.

Locul rădăcinilor este potrivit nu numai în analiza, ci și în sinteza sistemelor de reglare automată. Pentru determinarea locului rădăcinilor se pornește de la funcția de transfer a sistemului de reglare automată în circuit deschis:

, (3.1)

unde >0, mn și . Expresia (3.1) poate fi reprezentată și sub forma:

. (3.2)

Ecuația caracteristică a sistemului de reglare în circuit închis este:

1+G(s)=0. (3.3)

Din aceasta rezultă:

G(s)= -1/. (3.4)

Totalitatea curbelor descrise de numărul complex , cu condiția ca 0, reprezintă locul rădăcinilor căutat .Prin separarea lui (3.2) în modul și fază se obține condiția de amplitudine:

(3.5)

și condiția de fază

(3.6)

cu k=0,1,2,……. Aici și reprezintă unghiul corespunzător numărului complex (), respectiv (). În mod obișnuit condiția de fază este dependentă de . Toate punctele planului complex s ce îndeplinesc condiția de fază reprezintă locul geometric al tuturor polilor sistemului de reglare în circuit închis la variația factorului de amplificare . Parametrizarea, adică atribuirea la fiecare punct de pe curbă a unei valori a lui se face prin utilizarea condiției de amplitudine (3.5).

În cazul sistemului de suspensie activă la autovehicule, ecuațiile dinamicii sistemului sub forma funcțiilor de transfer sunt cele prezentate în capitolul anterior:

(3.7)

(3.8)

unde

(3.9)

Schema sistemului este următoarea:

w

Figura III.1 – Schema sistemului în buclă închisă

Cerințele de proiectare rămân aceleași, astfel încât atunci când perturbația șoselei (W) este simulată cu o intrare treaptă, ieșirea (X1-X2) să aibă un timp tranzitoriu mai mic de 5 secunde și o suprareglare de maxim 5%.

Se creează un fișier .m, care se găsește în Anexă, unde se introduc următoarele linii de comandă:

m1=290;

m2=59;

k1=16182;

k2=190000;

d1=1000;

0;

nump=[(m1+m2) d2 k2]

denp=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2]

num1=[-(m1*d2) -(m1*k2) 0 0]

den1=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2]

numf=num1;

denf=nump;

Mai întâi trebuie determinați polii sistemului deschis utilizând comanda roots(denp), atașată fișierul .m de mai sus:

R=roots(denp)

Matlab va returna:

R =

-8.7197 +58.0573i

-8.7197 -58.0573i

-1.4790 + 7.0674i

-1.4790 – 7.0674i

De aici observăm polii dominanți, cei care sunt complex conjugați și mai apropiați de axa imaginară(și de origine) si au un factor de amortizare (ξ) mai mic.

III.2 Trasarea locului rădăcinilor

Ideea de bază a proiectării prin metoda locului rădăcinilor este estimarea răspunsului sistemului în buclă închisă, din trasarea locului rădăcinilor ecuației caracteristice a sistemului închis. Prin adăugarea zerourilor și/sau polilor la sistemul original deschis (adică adăugând un regulator), locul rădăcinilor și răspunsul sistemului închis vor fi modificate.

Se trasează locul rădăcinilor sistemului original, adăugând la fișierul .m de mai sus următoarea comandă :

rlocus( nump,denp)

După rularea programului în fereastra de lucru a Matlab, se obține următorul grafic:

Figura III.2 – Locul Rădăcinilor sistemului original

Din desenul de mai sus, se poate observa că, există două perechi de poli și zerouri care sunt foarte apropiate unele de altele. Aceste perechi de poli și zerouri sunt apropiate de axa imaginară, motiv pentru care s-ar putea ca sistemul cu suspensie activă al automobilului să fie marginal stabil (la limita de stabilitate), lucru care ar putea cauza o problemă de stabilitate.

Pentru a evita problema unui sistem instabil, trebuie mișcați toți polii și zerourile în semiplanul stâng al planului complex.

În acest scop, se vor pune două zerouri foarte apropiate de cei doi poli de lângă axa imaginară a sistemului necompensat, pentru a face anularea polilor prin zerouri. Mai mult, se vor pune alți doi poli, mai departe pe axa reală, pentru a se obține și un răspuns rapid.

III.3 Adăugarea unui filtru notch ( filtru stop bandă, filtru crestătură)

Se cunoaște faptul că, dacă există poli sau zerouri în semiplanul drept al planului complex, se poate utiliza tehnica compensării polilor prin zerouri sau invers, dar că această tehnică nu este atât de eficientă, deoarece, datorită naturii sistemelor reale, nu se poate obține niciodată o compensare exactă. Atunci când un pol este compensat printr-un zero (sau un zero printr-un pol), acesta nu dispare în mod real, ci pur și simplu devine unul necontrolabil sau neobservabil. Din acest motiv această tehnică este în general evitată.

Totuși, un remediu îl constituie utilizarea unui filtru (unui controller) “tăietură”, ”crestătură” sau “șanț” (notch), ce are o caracteristică de frecvență “permisivă” pentru toate frecvențele, mai puțin pentru o anumită valoare . De aceea, de câte ori o funcție de transfer a unui sistem controlat conține una sau mai multe perechi de poli complex conjugați ce sunt apropiați de axa imaginară din planul s, se poate spune că avem de-a face cu un sistem în buclă închisă care este aproape instabil, adică unul care este foarte puțin amortizat.

Va fi probabil necesar, în cazul sistemului nostru, de două zerouri lângă cei doi poli apropiați de axa imaginară pentru a trasa locul rădăcinilor, și să rămână ca dominanți polii compensatorului în locul zerourilor de pe axa imaginară.

De asemenea, va fi nevoie de doi poli plasați departe în stânga, pentru a împinge locul rădăcinilor în stânga semiplanului stâng al planului complex(SSPC).

Se va verifica dacă un filtru notch(care este un controller cu două anticipări, adică cu două zerouri), poate să realizeze acest lucru. De aceea, se vor pune(asigna) zerourile la valorile 3±3.5i și polii la 30 și 60.

În fișierul .m se adaugă următoarele linii:

z1=3+3.5i;

z2=3-3.5i;

p1=30;

p2=60;

numc=conv([1 z1],[1 z2]);

denc=conv([1 p1],[1 p2]);

rlocus (conv(nump,numc),conv(denp,denc))

Schimbând axele obținem:

axis ([-40 10 -30 30])

Figura III.3 – Locul rădăcinilor utilizând un filtru ”notch”

“Câștigul” K, poate fi determinat utilizând comanda rlocfind :

[K, poles]=rlocfind (conv (nump,numc) , conv (denp , denc))

Figura III.4 – Locul rădăcinilor cu filtru “notch”

Se adaugă “câștigul” K obținut sistemului:

numc=K*numc.

III.4 Trasarea răspunsului sistemului închis

Se utilizează o treaptă de 0.1m ca perturbație și pentru a simula răspunsul sistemului închis la treaptă cu acest compensator (notch), se multiplică numa cu 0.1. Se adaugă următoarele comenzi în fișierul .m și se comentează (%) toate comenzile rlocus și rlocfind:

step (0.1*numa,dena)

title (‘Răspunsul sistemului în buclă închisa la treapta "w" de 0.1m/ filtru “notch”’)

Din răspunsul sistemului la treaptă pe perturbație se poate observa, că atunci când automobilul întâlnește pe șosea o treaptă de 0.1m, deviația maximă a corpului automobilului față de roată (sau față de șosea) este mica și oscilațiile se stabilizează în 2 secunde. Deci, acest răspuns al sistemului este satisfăcător.

Cap.IV Sinteza în frecvență a regulatorului utilizând caracteristicile Bode

IV.1 Trasarea răspunsului în frecvență utilizând comanda bode

Ideea principală a proiectării în frecvență constă în utilizarea caracteristicilor Bode ale funcției de transfer a sistemului în buclă deschisă, pentru a estima răspunsul sistemului în buclă închisă. Adăugând un regulator sistemului, acesta schimbă caracteristicile Bode în buclă deschisă, astfel încât și răspunsul sistemului în buclă închisă va fi modificat.

Caracteristicile Bode sunt reprezentări ale modulului și defazajul aferente lui G(j*w). Frecvența este în scară logaritmică, faza este dată în grade/modulul este dată în decibeli. Un decibel este definit astfel: 20*log10(|G(j*w|).

Pentru a trasa caracteristicile Bode aferente funcției de transfer a sistemului deschis original utilizăm fișierul .m creat în capitolele anterioare la care adăugăm următoarele linii și apoi rulăm:

w=logspace(-1,2);
bode(nump,denp,w);

Caracteristicile obținute sunt cele date în figura de mai jos:

Figura IV.1 – Caracteristicile Bode ale sistemului în buclă închisă

Pentru ușurință la reprezentarea sistemul cu diferite frecvențe naturale , normalizăm și scalăm rezultatele înainte de a trasa caracteristicile Bode, astfel încât asimptotele la infinit ale fiecărui termen să fie în 0 dB.

Această normalizare prin ajustarea amplificării K, face mai ușoară adăugarea componentelor pe caracteristicile Bode.

Efectul lui K este mișcarea curbei amplitudinii în sus (când K crește) sau în jos (când K descrește) printr-o valoare de 20logK, dar amplificarea K nu are nici un efect asupra curbei de fază. Deci, într-o caracteristică oarecare, K trebuie să fie egal cu 100 dB sau 100.000, pentru a deplasa caracteristica sus, deasupra lui 0 dB la 0.1 rad/s.

Pentru vizualizarea celor mai sus enunțate se adaugă următoarea linie înaintea comenzii bode și apoi se rulează programul care se află în Anexă:

nump=100000*nump

Se vor obține următoarele caracteristici Bode:

Figura IV.2 – Caracteristicile Bode ale sistemului în buclă deschisă

IV.2 Utilizarea unui regulator cu dublu avans de fază

IV.2.1 Considerații teoretice

Se știe că un compensator cu avans de fază (lead), poate crește stabilitatea sau viteza răspunsului sistemului; un compensator cu întârziere poate micșora (dar nu elimina) eroarea de reglare staționară. Rezultă că, funcție de efectul dorit, se pot utiliza unul sau mai multe compensatoare cu avans sau cu întârziere de fază, în diverse combinații.

IV.2.2 Compensatorul cu avans de fază proiectat utilizând locul rădăcinilor

Un compensator cu avans de ordinul I se poate proiecta utilizând locul rădăcinilor. Un compensator cu avans de fază sub forma locului rădăcinilor este dat de:

(4.1)

unde .

Un compensator cu avans de fază tinde să deplaseze locul rădăcinilor către semiplanul stâng al planului complex. Acest lucru duce la creștere stabilității sistemului și totodată la o creștere a vitezei răspunsului sistemului.

Ecuația care dă intersecția asimptotelor locului rădăcinilor cu axa reală este:

(4.2)

Când se adaugă un compensator cu avans la sistem, valoarea acestei intersecții va fi un număr negativ mai mare decât a fost cel anterior. Numărul net de zerouri și poli va fi același (un zero și un pol vor fi adăugați), dar polul adăugat este un număr negativ mai mare decât cel al zeroului adăugat(condiția).

Astfel, rezultatul adăugării unui pol compensatorului cu avans de fază este acela că intersecția asimptotelor se deplasează spre stânga, mai departe în semiplanul stâng al planului complex și, evident, întregul loc al rădăcinilor va fi deplasat către stânga.

Aceasta poate crește regiunea de stabilitate și, la fel, viteza de răspuns a sistemului.

În Matlab un compensator cu avans de fază este implementat utilizând funcția de transfer în forma:

numlead=Kc*[1 z];

denlead=[1 p];

și utilizând apoi funcția conv( ) pentru a-l implementa cu num-ul și den-ul procesului.

newnum=conv (num, numlead);

newden=conv (den, denlead)

IV.2.3 Compensatorul cu avans de fază proiectat prin utilizarea răspunsului în frecvență

Un compensator cu avans de fază de ordinul I se poate proiecta utilizând și răspunsul în frecvență. Un compensator cu avans de fază sub forma răspunsului în frecvență este dat de funcția de transfer:

(a>1) (4.3)

Se poate observa că, această formă este echivalentă cu forma locului rădăcinilor:

, (4.4)

în care p = 1/T, z = 1/aT și Kc = a.

La metoda proiectării compensatoarelor cu avans de fază utilizând răspunsul în frecvență, acest compensator adaugă un defazaj pozitiv sistemului în plaja de frecvențe de la 1/aT la 1/T.

Caracteristicile Bode ale unui compensator cu avans de fază sunt ilustrate mai jos: Figura IV.3 – Caracteristicile Bode ale unui compensator cu avans de fază

Cele două frecvențe de frângere sunt la 1/aT și 1/T. Reținem că defazajul pozitiv este adăugat sistemului între cele două frecvențe de frângere. Funcție de valoarea lui a, defazajul maxim adăugat poate fi de până la ; dacă însă este necesar un defazaj mai mare de , atunci pot fi utilizate două compensatoare cu avans de fază.

Valoarea maximă a defazajului este adăugată la frecvența naturală, care este plasată la:

(4.5)

Ecuația care determină defazajul maxim este:

(4.6)

Defazajul pozitiv suplimentar, conduce la creșterea marginii de fază și astfel crește stabilitatea sistemului. Acest tip de compensator este proiectat prin determinarea parametrului „a” din valoarea defazajului necesar satisfacerii cerințelor marginii de fază și prin determinarea parametrului „T”, pentru a-l plasa(adăuga) la noua pulsație de tăiere. Compensatorul cu avans crește amplificarea sistemului la frecvențe ridicate, valoarea acestei amplificări fiind egală cu „a”. Acest lucru duce la micșorarea timpului de creștere și a timpului tranzitoriu.

În Matlab, un compensator cu avans de fază în forma răspunsului la frecvență este implementat astfel:

numlead=[aT 1];

denlead=[T 1];

și apoi utilizând funcția conv ( ) pentru a multiplica num-ul și den-ul procesului:

newnum=conv(num,numlead);

newden=conv(den,denlead);

În cazul sistemului nostru, din caracteristicile Bode(figurile IV.1 și IV.2) se observă că, curba de fază este concavă la cca 5rad/sec. Mai întâi se încearcă adăugarea unei faze pozitive (defazaj pozitiv suplimentar) în jurul acestei pulsații de 5 rad/sec, astfel încât defazajul să rămână deasupra liniei de -180.

Din moment ce o marginea de fază mare conduce la o suprareglare (σ) mică, va trebui să adăugăm cel puțin +140(defazaj pozitiv) în zona apropiată de 5 rad/sec.

Pentru că un regulator cu avans de fază poate adăuga maximum +90, vom folosi un regulator cu dublu avans de fază, respectiv:

(4.7)

Pentru obținerea parametrilor a și T, se pot utiliza următorii pași:

Pas 1. Determinarea defazajului pozitiv necesar:

Deoarece se dorește o valoare totală de 140o, vor fi necesare două a câte 70o, pentru fiecare controller.

Pas 2. Determinarea frecvenței ”w” unde trebuie adăugat defazajul:

În cazul nostru, această frecvență ar trebui să fie în jurul frecvenței de 5 rad/sec.

Pas 3. Determinarea constantei” a” din ecuația de mai sus:

Această constantă determină de fapt spațiul necesar dintre zeroul și polul necesar pentru defazajul maxim ce trebuie adăugat.

(4.8)

Pas4. Determinarea parametrilor ”T” și ” aT” din ecuațiile următoare:

Acești parametri determină frecvențele de frângere, așa încât defazajul maxim va fi adăugat la frecvența dorită.

(4.9)

Pentru a vedea ce se întâmplă în continuare cu caracteristicile Bode când se introduce regulatorul cu dublu avans, se adaugă următoarele linii de cod în fișierul .m anterior:

numc = conv([1.13426 1], [1.13426 1]);

denc = conv([0.035265 1], [0.035265 1]); margin(conv(nump,numc),conv(denp,denc)

Prin rularea programului se obțin următoarele caracteristici Bode:

Figura IV.4 – Caracteristicile Bode când se adaugă regulator cu dublu avans de fază

Caracteristicile de mai sus, sunt echivalente cu caracteristicile obținută prin rularea programului, după ce se adaugă la acesta următoarele linii de cod:

w=logspace(-4,4);

[mag,phase,w]=bode(conv(nump,numc),conv(denp,denc),w);

subplot (2,1,1);

semilogx (w, 20*log10(mag));

grid

title (‘Caracteristicile Bode când se aduga un regulator cu dublu avans de fază’)

xlabel ('Frecventa (rad/s)')

ylabel ('20logM')

subplot (2,1,2);

semilogx (w,phase)

axis ([1e-4, 1e4, -180, 360])

grid

xlabel ('Frecventa (rad/s)')

ylabel ('Faza (grade)')

Figura IV.5 – Caracteristicile Bode când se adaugă un regulator cu dublu avans de fază

Din aceste caracteristici se observă că porțiunea concavă a caracteristicii fază-frecvență este acum deasupra liniei corespunzătoare lui -180 în jurul frecvenței de 5 rad/s și deci marginea de fază este suficient de mare pentru criteriul de proiectare.

Pentru a vedea cum ieșirea (X1-X2) răspunde acum unei denivelări (W) de pe șosea. Mai întâi trebuie să obținem funcția de transfer în buclă închisă de la (W) la (X1-X2), respectiv în programul principal vom adăuga următoarele linii de comandă:

numa=conv(numf,denc);

dena=polyadd (conv(denp,denc),conv(nump,numc));

IV.3 Trasarea răspunsului sistemului în buclă închisă

Pentru a simula treapta cu înălțimea de 0.1m, se înmulțește numa cu 0.1 și se adaugă următoarele linii de comandă în fișierul .m:

t=0:0.01:5;

step(0.1*numa,dena,t)

axis([0 5 -.01 .01])

După rularea programului se obține următorul răspuns:

Figura IV.6 – Răspunsul sistemului în buclă închisă

Amplitudinea răspunsului la treaptă este mult mai mică decât procentul suprareglării cerute și timpul tranzitoriu de asemenea este mai mic de 5 secunde. Se observă că amplitudinea ieșirii este mai mică de 0.0001m sau 1% din amplitudinea intrării după 4 secunde. Deci, putem vedea din desenul de mai sus că timpul tranzitoriu este de 4 secunde. Din caracteristicile Bode prezentate anterior se constată că prin creșterea amplificării buclei, va crește și frecvența de tăiere și acest lucru va face ca răspunsul să fie mai rapid.

Cap.V Sinteza unui regulator optimal (LQR) pentru SSA

V.1 Considerații teoretice

Termenul LQ este o prescurtare pentru regulatorul liniar optimal cu indice de performanță pătratic. Regulatorul LQ este cel mai important rezultat al controlului modern în abordare de stare.

La un regulator LQ, performanța este estimată cu un indicator de performanță scalar de forma:

(5.1)

unde R și Q sunt matrice simetrice, pozitiv semidefinită(R) respectiv pozitiv definită(Q).

Indicatorii de performanță utilizați în probleme de control optimal pătratic pot fi:

-indicatori de performanță complecși ( hermitieni – vezi integrala de mai sus) în care Q, R sunt matrici simetrice hermitiene semidefinită respectiv pozitiv definită.

-indicatori de performanță reali, în cazul sistemelor cu vectori reali și matrici reale.

Sistemul de control optimal trebuie să minimizeze indicatorul de performanță(J), în cazul în care un sistem este stabil. Există mai multe abordări pentru soluționarea acestui tip de problemă, dintre care cea mai întâlnită este cea bazată pe metoda a doua a lui Lyapunov (1892) de determinare a stabilității sistemelor dinamice.

Problema care se pune este una de reglare, în care obiectivul este de a conduce starea de la o valoare inițială oarecare către valoarea 0.

Termenul: penalizează abaterea de la origine (x(t)=0) a vectorului x(t). Penalizarea este pătratică, pentru că termenii și sunt forme pătratice.

(5.2)

Întrucât penalizarea este pătratică, rezultă că erorile mari sunt penalizate mai mult decât erorile mai mici.

Este de dorit ca termenul să aibă o valoare mică întrucât acest lucru implică faptul că starea x(t) este relativ lângă origine.

Încercările de a proiecta optimal sistemul, conduc la o proiectare ce utilizează valori mari la intrare. Rațiunea este faptul că intrările mari (impulsurile) cauzează schimbări rapide ale stării și astfel x(t) este dus la origine mult mai rapid decât dacă s-ar utiliza intrări de amplitudine moderată.

Termenul este inclus în ecuația lui J pentru a nu pierde din vedere amplitudinea intrării.

O proiectare care utilizează valori mari ale intrării va duce la o valoare mare a lui J, pentru că termenul este mare.

Pe de altă parte, utilizând intrări mici, soluția rezultată nu va fi una bună pentru conducerea lui x(t) către 0, astfel că primul termen din J va fi unul mare. O proiectare care minimizează pe J reprezintă un compromis între cele două deziderate conflictuale: adică ale unei bune reglări și a unor intrări rezonabile ca mărime. Privim problema de control ca și o problemă de optimizare nu din cauza faptului că se dorește minimizarea indicatorului de performanță (J), ci pentru că soluția optimală este ușor de rezolvat matematic. Cadrul de optimizare oferă o manipulare convenabilă a proprietăților soluției.

Regulatorul LQ optimal poate fi definit astfel:

-fiind dat un sistem liniar invariant în timp de forma:

, (5.3)

x(0) – stare inițială,

se cere să se calculeze o matrice K a câștigurilor pe reacție de la stare, astfel încât legea de control u(t) = -Kx(t) să minimizeze indicele de performanță J, unde Q și R sunt simetrice.

Sistemul închis și indicele de performanță sunt descrise de ecuațiile:

(5.4)

(5.5)

Observație. În scopul ca problema să fie de interes practic, trebuie să se obțină stabilizarea internă cu o reacție de la stare, caz în care sistemul este, adică trebuie să fie, stabilizabil. Oricare dintre modulele instabile trebuie să fie controlabile. Printr-un anume K putem obține un minim a lui J. Pentru ca un K să minimizeze pe J pentru toți trebuie satisfăcută condiția necesară :

(5.6)

unde P satisface ecuația matricială Riccati:

. (5.7)

Trebuie deci rezolvată ecuația de forma (5.7) pentru determinarea lui P. Se calculează apoi K cu relația (5.6) și se verifică A – B*K din punct de vedere al stabilității. Întrucât ecuația Riccati este neliniară, în general ea are mai multe soluții, deci va trebui găsită o soluție care să stabilizeze sistemul. Impunem condiția ca (A,) să fie detectabilă.

Dacă (A,) nu este detectabilă, rezultă că există un mod instabil și neobservabil cu valoarea proprie și vector propriu astfel încât:

și și (5.8)

Dacă (A,) este detectabilă, se arată ca o matrice pozitiv definită (semidefinită) P, poate conduce la găsirea unei matrice K stabilizabilă, care va conduce sistemul cel într-un minim local.

Se poate arăta, că ecuația Riccati are una și numai o singură soluție P>=0; rezultă că există întotdeauna un minim local și doar unul singur, care este și minim global.

ALGORITM:

Pas 1. Pentru a rezolva problema reglării, ecuațiile de stare sunt scrise în formă incrementală:

(5.9)

(5.10)

Indicele de performanță se definește deci cu privire la și .

Pas 2. Ecuația Riccati are o unică soluție pozitiv definită sau semidefinită.

Pas 3. Odată ce matricea P este cunoscută, matricea optimală K se calculează cu relația (5.6).

Pas 4. Dată fiind matricea ampificării K, rezultă răspunsul la o stare inițială dată, care se calculează rezolvând ecuația: pentru . (5.11)

Starea, comanda și ieșirea se obțin cu relațiile:

(5.12)

V.2 Rezolvarea problemei de control optimal în Matlab

Comanda K = lqr(A,B,O,R) – rezolvă problema regulatorului LQ liniar în timp continuu iar comanda “dlqr” în discret. Comanda “lqr” calculează matricea K a amplificărilor optimale de pe reacție de la stare, astfel încât legea de comandă u = -K*x să minimizeze J care este afectat ecuației cu constrângeri .

O altă comandă este: [K,P,E] = lqr(A,B,Q,R) care returnează matricea P, care este soluție pozitiv definită unică, a ecuației matriciale Riccati asociate. Dacă matricea A-B*K este stabilă, atunci o soluție P pozitiv definită există întotdeauna. Polii sistemului închis sau valorile proprii ale matricei A-B*K sunt obținute prin această comandă.

Observație. Uneori, la unele sisteme, matricea A-B*K nu poate fi făcută o matrice stabilă, indiferent de matricea K aleasă. În acest caz nu există o matricea P pozitiv definită pentru ecuația Riccati.

V.3 Regulator LQ pentru sistemul de suspensie activă

Fiind dat sistemul descris prin ecuațiile:

(5.13)

trebuie să se minimizeze indicele de performanță:

. (5.14)

Pentru determinarea unei legi de control optimal presupunem că referința este zero. Matricea k=[k1 k2 k3 k4] este aleasă astfel încât să fie minimizat indicele de performanță pentru care:

; R=1 și (5.15)

Vom alege pe q11 suficient de mare comparativ cu q12, q13, q14 conform celor spuse mai sus. Astfel, vom alege q11=1010, q12=q13=q14=1 și R=1.

Se creează fișierul .m cu următoarele linii de comandă, care se află și în Anexa A :

m1=290;

m2=59;

k1=16182;

k2=190000;

d1=1000;

d2=0;

A=[0 1 0 0; -(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) – (d1/m1); d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1;k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0]

B=[0;1/m1;0;(1/m1)+(1/m2)]

C=[0 0 1 0];

D=[0];

%Introducerea matricelor Q si R care pondereaza starea, respectiv comanda

Q=[10000000000 0 0 0 ;0 1000 0 0; 0 0 1000 0; 0 0 0 1000];

R=[1]

Pentru obținerea matricei K utilizăm comanda:

K=lqr(A,B,Q,R)

care este atașată fișierului .m de mai sus.

Prin rularea programului obținem următoarele valori pentru matricea K:

K=1.0e+004 *[6.5455 0.6962 2.3201 -0.0173]

V.4 Trasarea răspunsului sistemului în buclă închisă

Se poate obține răspunsul sistemului închis prin simpla adăugare în fișierul .m a următoarelor linii:

Aa=A-B*K;

Bb=B*K(1);

Cc=C;

Dd=0;

t=0:0.01:5;

[y,x,t]=step(Aa,0.1*Bb,Cc,Dd,1,t);

plot (t,y);

title ('Răspunsul indicial al a sistemului automat optimal pătratic');

xlabel ('Timpul(s)');

ylabel ('Ieșirea (X1-X2)');

grid;

Prin rularea fișierului .m se va obține următorul răspuns:

Figura V.1 – Răspunsul indicial al sistemului cu regulator optimal linear pătratic(LQR)

Din răspunsul de mai sus se observă că cerințele de proiectare sunt satisfăcute.

Cap.VI Proiectarea unui regulator digital cu reacție de la stare

Există două metode de proiectare a unui regulator digital cu reacție de la stare:

-metoda reproiectării în timp discret;

-metoda proiectării directe în domeniul timp.

VI.1 Metoda de reproiectare în timp discret

Pentru proiectarea unui regulator digital cu reacție de la stare folosim metoda de reproiectare in timp discret. În cadrul acestei metode indirecte avem șase pași importanți de parcurs:

VI.1.1 Definirea cerințelor de proiectare

-σ5% și tt<5[sec], adică suprareglarea să fie sub 5% la semnale treaptă de 0.1m, iar timpul tranzitoriu să nu dureze mai mult de 5 sec.

VI.1.2 Selectarea perioadei de eșantionare

Problema „indirectă” constă în convertirea modelului matematic în timp continuu, aferent procesului, în modelul matematic echivalent, dar în timp discret. Pentru aceasta se selecteaza o perioadă de eșantionare adecvată “T”.

În cazul suspensiei active, alegerea corectă a acestei perioade de eșantionare este foarte importantă întrucât o treaptă(o denivelare sau groapă) în suprafața șoselei pe care rulează mașina afectează foarte rapid mărimea de ieșire. Din punct de vedere fizic, ceea ce se intamplă este faptul că suprafața șoselei coboară(sau urcă) brusc roata automobilului, compresând(la ridicare) sau intinzând(la coborâre) arcul kt.

Întrucât masa suspensiei este relativ scazută și arcul kt este extrem de rigid, masa suspensiei mu se înalță sau coboară foarte rapid, crescând pe zu aproape imediat. Deoarece regulatorul poate “să vadă” efectul perturbației numai după o perioadă de eșantionare completă Te, aceasta trebuie aleasă suficient de mică, astfel încât ieșirea să nu depasească suprareglajul impus(5%), într-o singură perioadă de eșantionare. Pentru a selecta perioada de eșantionare Te, va trebui să examinăm cu atenție începutul răspunsului la treaptă. Se va simula doar începutul acestui răspuns, setând vectorul de timp pentru funcția treaptă la intrare, pentru o plajă de valori de la 0 la 0.0005 secunde(de la 0 la 0.5msec). Pentru o groapă de 0.1 metri, acesta este simulat prin înmulțirea matricei B cu factorul 0.1.

VI.1.3 Conversia continuu-discret

După ce se alege perioada de eșantionare, putem converti modelul matematic continuu într-unul de tip discret. Pentru a converti matricele A,B,C,D din continuu în discret(Ad, Bd, Cd, Dd), utilizăm Matlab(comenzile c2dm,c.da au 6 argumente: patru matrici de stare Ad, Bd, Cd, Dd, perioada de eșantionare Te și tipul circuitului de extrapolare).

Matlab va intoarce:

Ad =

1.0000 0.0005 0.0000 0.0000

-0.0116 0.9996 0.0113 0.0004

0.0001 0.0000 0.9995 0.0005

0.2951 0.0093 -1.8976 0.9902

Bd =

1.0e-005 *

0.0000

0.0332

-0.0002

-0.8432

Cd =

1 0 0 0

Dd =

0

VI.1.4 Adăugarea unui integrator(pentru eroare staționară nulă)

Integratorul are scopul de a conduce(aduce) răspunsul în regim staționar la 0. Integratorul se adăugă în serie cu procesul și are efectul adăugării unei alte stări, suplimentare, la proces.

Se adaugă integratorul la reprezentarea in spațiul stărilor și se utilizează comanda series. Această comandă ia matricile A, B, C, D ale celor două sisteme ce trebuie conectate în serie, ca argumente, și returnează un nou set de matrici A, B, C, D.

În timp discret și în spațiul stărilor, un integrator poate fi reprezentat ca o aproximare trapezoidală a integrării peste fiecare perioadă de eșantionare, adică în intervalul fiecărei perioade de eșantionare, după cum urmează:

x(k+1)=x(k)+Tu(k);

y(k)=x(k)+T/2u(k);

y=z;

u=f;

Ai =

1

Bi =

1

Ci =

5.0000e-004

Di =

2.5000e-004

Ada =

1.0000 1.0000 0 0 0

0 1.0000 0.0005 0.0000 0.0000

0 -0.0116 0.9996 0.0113 0.0004

0 0.0001 0.0000 0.9995 0.0005

0 0.2951 0.0093 -1.8976 0.9902

Bda =

1.0e-005 *

0

0.0000

0.0332

-0.0002

-0.8432

Cda =

1.0e-003 *

0.5000 0.2500 0 0 0

Dda =

0

Cda =

1 0 0 0 0

VI.1.5 Proiectarea controlerului

Structura acestuia este similară cu cea din timp continuu. Utilizăm comanda place pentru a calcula matricea K, a amplificărilor pe reacție de la stare, care va fi dată de polii doriți ai sistemului închis.

Mai întâi este necesar să decidem unde vor fi plasați polii sistemului închis. În particular, cei cinci poli ai sistemului închis pot fi plasati astfel încât să compenseze toate zerourile procesului și astfel să obținem răspunsul dorit.

Pentru început vom determina zerourile procesului făcând conversia ecuațiilor de stare generale într-o funcție de transfer echivalentă și apoi vom determina rădăcinile numărătorului acestei funcții de transfer. Vom utiliza comanda ss2tf care va avea ca argumente matricile de stare și intrarea selectată, iar ca ieșiri numărătorul și numitorul funcției de transfer.

Se adaugă comenzile următoare :

[num,den]=ss2tf(Ad,Bd,Cd,Dd,1);

zoros=roots(num)

Matlab va intoarce:

zeros =

-0.9968

0.9996 + 0.0285i

0.9996 – 0.0285i

Se vor selecta aceste trei zerouri ca fiind trei dintre cei cinci poli doriți ai sistemului închis. Unul dintre cei doi poli rămasi va fi selectat la locația 0.9992 întrucât acolo un pol stabilizează sistemul in aproximativ 10000 de eșantioane. Ultimul pol, cel de-al cincilea, va fi selectat la z=0.94, întrucât aceasta valoare este suficient de rapidă pentru a fi nesemnificativă

VI.1.6 Simularea răspunsului sistemului închis 

Simulăm pentru o treaptă de –0,1[m](o groapă) și o treaptă de +0,1[m](un dâmb):

ms=290;

mu=59;

kr=35000;

br=1100;

kt=190000;

g=9.8;

A=[0 1 0 0;-kr/ms -br/ms kr/ms br/ms;0 0 0 1;kr/mu br/mu -(kr/mu+kt/mu) -br/mu]

B=[0;1/ms;0;-1/mu]

C=[1 0 0 0]

D=[0]

p=[0;0;0;kt/mu]

comp_lib=[0;-g;0;-g]

T=0.0005;

[Ad,Bd,Cd,Dd]=c2dm(A,B,C,D,T,'zoh')

Ai=1

Bi=1

Ci=T

Di=T/2

[Ada,Bda,Cda,Dda]=series(Ad,Bd,Cd,Dd,Ai,Bi,Ci,Di)

Cda=[Cd 0]

[num,den]=ss2tf(Ad,Bd,Cd,Dd,1)

zeros=roots(num)

p1=zeros(1);

p2=zeros(2);

p3=zeros(3);

p4=0.9990;

p5=0.5;

K=place(Ada,Bda*[1 0],[p1,p2,p3,p4,p5])

yies=dstep(Ada-Bda*[1;0]'*K,-0.1*Bda,Cda,-0.1*Dda,1,10001)

t=0:0.0005:5

stairs(t,yies)

title('Răspunsul sistemului pentru semnal treaptă unitară de -0.1[m]')

xlabel(‘time(sec)’)

ylabel(‘amplitudinea[m]’)

Figura VI.1- Răspunsul sistemului pentru o treaptă(groapă) de -0.1[m]

Pentru situația unei trepte pozitive(a unui dâmb) de +0.1[m] rulăm următorul program:

ms=290;

mu=59;

kr=35000;

br=1100;

kt=190000;

g=9.8;

A=[0 1 0 0;-kr/ms -br/ms kr/ms br/ms;0 0 0 1;kr/mu br/mu -(kr/mu+kt/mu) -br/mu];

B=[0;1/ms;0;-1/mu];

C=[1 0 0 0];

D=[0];

p=[0;0;0;kt/mu];

comp_lib=[0;-g;0;-g];

T=0.0005;

[Ad,Bd,Cd,Dd]=c2dm(A,B,C,D,T,'zoh');

Ai=1;

Bi=1;

Ci=T;

Di=T/2;

[Ada,Bda,Cda,Dda]=series(Ad,Bd,Cd,Dd,Ai,Bi,Ci,Di);

Cda=[Cd 0];

[num,den]=ss2tf(Ad,Bd,Cd,Dd,1)

zeros=roots(num)

p1=zeros(1);

p2=zeros(2);

p3=zeros(3);

p4=0.9990;

p5=0.5;

K=place(Ada,Bda*[1 0],[p1,p2,p3,p4,p5])

yies=dstep(Ada-Bda*[1;0]’*K,-0.1*Bda,Cda,-0.1*Dda,1,10001)/0.00000016;

t=0:0.0005:5;

stairs(t,yies)

title('Răspunsul sistemului pentru semnal treaptă unitară de +0.1[m]')

xlabel(‘time(sec)’)

ylabel(‘amplitudinea[m]’)

Figura VI.2 – Răspunsul sistemului pentru o treaptă pozitivă(dâmb) de +0.1[m]

Cap.VII Controlul robust al suspensiei active

VII.1 Modelul de stare al suspensiei active pentru un sfert de mașină

Modelul de suspensie pentru o singură roată va fi folosit pentru a modela legile de control ale suspensiei active a întregului autovehicul.

Figura VII.1 – Modelul de suspensie activă pentru o singură roată

În figura de mai sus, pentru unui microbuz, notațiile reprezintă:

ms = M = m1 = 290 kg – masa corpului autovehiculului,

mus = m2 = 59 kg – masa suspensiei(dar fără arcuri și amortizor),

ks = k1 = 16182N/m – constanta elastică a arcului sistemului de suspensie,

kt = k2 =190000 N/m – constanta elastică a roții pneumatice,

bs = d1 = 1000 Ns/m – constanta de amortizarea a amortizorului sistemului de suspensie,

fs = u – forța[kN] (mărimea) de comandă/control = forța de suspensie activă (mărime de xxxxxxxx intrare),

xs = x1 – variația corpului vehiculului (mărime de ieșire = mărime reglată/controlată),

xus = x2 – variația roții de mers pe drum (mărime de intrare = mărime necontrolată),

r = w – perturbația (înălțimea căii de rulare, deci tot o mărime de intrare necontrolată).

Forța fs a suspensiei active , aplicată între masele corpului autovehiculului și axul(osia) șasiului este intrarea de control, control realizat prin feedback. Ea reprezintă componenta activă a sistemului de suspensie. În acest exemplu, dinamica actuatorului(elementului de execuție, adică a elementului hidraulic al suspensiei active), precum și forța fs sunt ignorate.

Definim:

X1 = xs,

X2 = s,

X3 = xus,

X4 = us,

și rezultă următoarea descriere(modelul matematic) pentru dinamica sfertului de autovehicul(1/4 car):

1 = X2,

2 =

3 = X4,

4 = .

Vom introduce următoarele valori în programul din Matlab:

ms=290; % kg

mus=59; % kg

bs=1000; % N/m/s

ks=16182 ; % N/m

kt=190000; % N/m

Un model liniar, invariabil în timp, al sfertului de mașină studiat, qcar, este construit apoi, mai jos, din ecuațiile de mișcare și valorile parametrilor. Intrările modelului, qcar, sunt denivelările drumului și forța actuatorului respectiv, iar ieșirile sunt abaterea corpului autovehiculului, accelerația și abaterea suspensiei(x1-x3).

A12=[ 0 1 0 0; [-ks -bs ks bs]/ms ];

A34=[ 0 0 0 1; [ks bs -ks-kt -bs]/mus];

B12=[0 0; 0 10000/ms];

B34=[0 0; [kt -10000]/mus];

C=[1 0 0 0; A12(2,:); 1 0 -1 0; 0 0 0 0];

D=[0 0; B12(2,:); 0 0; 0 1];

qcar=ss([A12; A34], [B12; B34], C, D);

Se știe că funcția de transfer a accelerației[Hedrick 90] are un punct invariabil la o frecvență de oscilație a pneului de 56.7 rad/sec. În mod similar, funcția de transfer a abaterii suspensiei are un punct invariabil la o frecvență de 23.3 rad/sec, „frecvența spațiului de accelerare”. Compromisul dintre confortul pasagerilor și abaterea suspensiei se datorează faptului că nu este posibil ca să fie ținute simultan amândouă mici și peste funcția de transfer din jurul frecvenței de oscilație a cauciucului, și în același timp, și la joasă frecvență.

VII.2 Modelul suspensiei active, pentru proiectarea controlerului robust liniar H∞

Modelul controlerelor liniare ale suspensiei trebuie să aleagă între confortul pasagerilor și abaterea suspensiei. Controlerele din acest subcapitol sunt proiectate folosind sinteza liniară H∞.

Ca un standard în cadrul H∞, obiectivele de performanță sunt atinse prin minimizarea normelor funcțiilor de transfer ponderate.

Funcțiile ponderate urmăresc atingerea a două scopuri in cadrul H∞:

-permit comparația directă a diferitelor obiective de performanță cu aceeași normă;

-permit informațiilor referitoare la frecvență să fie introduse în analiză.

O schemă bloc a legăturilor modelului proiectat pentru controlul H∞ în problema suspensiei active este prezentată in figura următoare:

Figura VII.2 – Schema bloc a modelului ce utilizează proiectarea unui controler H∞ liniar

Ieșirea măsurată sau semnalul de reacție, y, este abaterea suspensiei x1-x3. Controlerul acționează asupra acestui semnal pentru a produce intrarea de control, forța elementului de execuție hidraulic fs. Blocul Wn servește pentru a modela zgomotul senzorului pe canalul de măsură. Wn este setat la valoarea de zgomot a senzorului de 0.01 m.

Wn=0.01;

Într-o proiectare mai realistă, Wn ar trebui să fie dependent de frecvență și ar putea servi la modelarea zgomotului asociat cu senzorul de deplasare.

Ponderea Wref este utilizată pentru a grada magnitudinile perturbațiilor drumului. Presupunem că perturbația maximă a drumului este de 7 cm și deci alegem Wref=0.07.

Wref=0.07;

Magnitudinea si conținutul de frecvență al forței de control fs sunt limitate de către funcția de ponderare Wact. Alegem:

Wact

Modulul ponderii crește pentru valori mai mari de 50 rad/sec, în scopul de a limita variația lărgimii de bandă.

Wact=(100/13)*tf([1 50], [1 500]);

Deoarece se urmărește obținerea atât a unui confort ridicat pentru pasageri cât și eliminarea mișcărilor de ruliu și tangaj ale mașinii, se vor aborda trei proiectări de controlere H∞, urmând ca în final acestea să fie comparate pentru a se obține o soluție de compromis care să satisfacă întreaga dorință.

VII.3 Proiectarea controlerului H∞ liniar cu funcția hinfsyn, pentru a minimiza funcția de transfer a abaterii corpului mașinii(x1)

Scopul funcțiilor de ponderare Wx1 și Wx1-x3 este acela de a menține mici atenuarea copului mașinii, cât si a abaterii suspensiei, în plajele de frecvență dorite În această primă proiectare, vom proiecta controlerul în scopul creșterii confortului pasagerilor și deci, va fi penalizată abaterea corpului autovehiculului, x1, prin Wx1:

Wx1=8*tf(2*pi*5, [1 2*pi*5]);

Modulul ponderării este deasupra valorii de 5*2π rad/sec, pentru a respecta binecunoscuta regulă empirică H∞ care cere ca ponderile de performanță să fie deasupra unui zerou al buclei deschise(56.7 rad/sec în acest caz). Ponderea abaterii suspensiei, Wx2-x3, nu este inclusă în această formulare a problemei de control.

Putem construi modelul H∞ al procesului(automobilului) pentru proiectarea buclei de control, notat qcaric1, folosind comanda sysic. Semnalul de control corespunde intrării în qcaric1, adică lui fs. Accelerația corpului mașinii, care este zgomotoasă, este semnalul măsurat și corespunde ultimei ieșiri a lui qcaric1.

systemnames=’qcar Wn Wref Wact Wx1’;

inputvar=’[d1;d2;fs]’;

outputvar=’[Wact;Wx1;qcar(2)+Wn]’;

input_to_qcar=’[Wref;fs]’;

input_to_Wn=’[d2]’;

input_to_Wref=’[d1]’;

input_to_Wact=’[fs]’;

input_to_Wx1=’[qcar(1)]’;

qcaric1=sysic;

Un controler H∞ este sintetizat cu comanda hinfsyn. Această funcție are o intrare de control, ncont, forța elementului actuator hidraulic și o alta,un semnal de măsurat, nmeas, adică accelerația corpului mașinii.

ncont=1;

nmeas=1;

[K1,Scl1,gam1]=hinfsyn(qcaric1,nmeas,ncont);

CL1=lft(qcar([1:4 2], 1:2), k1);

sprintf(‘Controlerul H∞ K1 obținut are norma de %2.5g’, gam1)

ans=

Controlerul H∞ K1 obținut are norma de 0.56698

Putem analiza controlerul H∞ construind sistemul în buclă închisă, CL1. Diagramele Bode ale suspensiei pasive și ale celei active sunt prezentate in figura următoare:

Figura VII.3 – Diagramele Bode ale suspensiilor pasivă(-∙-∙) și activă(K1)

VII.4 Proiectarea controlerului H∞ liniar cu funcția hinfsyn pentru a minimiza funcția de transfer a abaterii suspensiei(x1-x3)

În continuare vom încerca să menținem funcția de transfer a abaterii suspensiei la o valoare cât mai mică. Efectele denivelărilor șoselei(perturbațiile șoselei) asupra abaterii suspensiei(x1-x3) sunt penalizate prin intermediul funcției de ponderare Wx1-x3. Magnitudinea(modulul) ponderii Wx1-x3 variază la valori mai mari(peste) decât valoarea de 10 rad/sec și în jurul(înainte și după) valorii unui zerou al funcției de transfer a sistemului deschis(23.3 rad/sec) în proiectare.

Wx1x3=25*tf(1, [1/10 1]);

Ponderea abaterii mașinii mașinii, Wx1, nu este inclusă in această formulare a problemei de control. Putem construi modelul procesului H∞ ponderat pentru proiectarea controlerului model, notat qcaric2, folosind comanda sysic. Sunt folosite aceleași mărimi și funcții ca în exemplul precedent.

M=iconnect;

d=icsignal(2);

fs=icsignal(1);

ycar=icsignal(size(qcar,1));

M.Equation{1}=equate(ycar,qcar*[Wref*d(1); fs]);

M.Input=[d;fs];

M.Output [Wact*fs;Wx1*ycar(1);ycar(2)+Wn*d(2)];

qcaric2=M.System;

Al doilea controler H∞ este sintetizat folosind aceeași comandă hinfsyn:

[K2,Scl2,gam2]=hinfsyn(qcaric2,nmeas,ncont);

CL2=lft(qcar([1:4 2],1:2),K2);

sprintf('Controlerul H∞ K2 obținut are norma de %2.5g', gam2)

ans =

Controlerul H∞ K2 obținut are norma de 0.89949

Trebuie reamintit că acest model de control ce utilizează H∞ nu pune pe primul plan minimizarea efectelor suspensiei asupra confortului pasagerilor așa cum făcea modelul de mai sus(cap.VII.3) care se axa numai pe confortul pasagerilor.

Putem analiza controlerul H∞ construind sistemul în buclă închisă(procesul qcaric2 și controlerul K2), adică CL2.Caracteristicile Bode ale funcției de transfer datorate denivelărilor șoselei asupra abaterii suspensiei pentru ambele tipuri de controlere cât și a sistemului de suspensie pasivă, sunt prezentate în figura următoare(fig. VII.4). Caracteristicile Bode care urmează se referă la cele două funcții de transfer de pe canalul intrării referinței r(perturbația șoselei) și cele două mărimi de ieșire, respectiv accelerația corpului mașinii și abaterea(deplasarea) suspensiei(fig. VII.4 și respectiv fig. VII.5).

Figura VII.4 – Diagramele bode ale suspensiei pasive și ale celor active utilizând controlerele K1 și K2

Liniile punctată și cea continuă din figură sunt răspunsurile în frecvență ale sistemelor deschise ce rezultă utilizând cele trei funcții de ponderare a performanței diferite(adică, sistemele deschise fără controler, cu controlerul K1 și cu controlerul K2, ultimele fiind controlerele rezultate din proiectările din cap.VII.3(K1) și respectiv cap.VII.4(K2). Observăm o reducere a abaterii suspensiei in jurul frecvenței ω1=56.7 rad/sec, și corespunzător, o creștere a răspunsului în frecvență a accelerației în această vecinătate. Mai observăm, în comparație cu prima proiectare(cap.VII.3), o reducere a abaterii suspensiei pentru frecvențe situate sub(mai mici decât) frecvența ω2=23.3 rad/sec.

Al doilea model de control atenuează ambele moduri de rezonanță, în timp ce primul model de controler își focalizează eforturile asupra primului mod, la frecvența de 23.3 rad/sec. Acestea se observă și in figura următoare:

Figura VII.5 – Diagramele Bode ale suspensiei pasive și ale celor active utilizând controlerele K1 și K2

VII.5 Analiza în domeniul timp

Toate analizele de până acum s-au făcut în domeniul frecvență. Caracteristicile de performanță în domeniul timp sunt hotărâtoare pentru succesul sistemului de suspensie activă la autovehicule. Răspunsurile în timp ale celor două controlere H∞ sunt arătate in figurile următoare. Liniile punctată, întreruptă și continuă corespund suspensiei pasive, controlerelor H∞ unu și doi(K1 și K2), respectiv.

Toate răspunsurile corespund perturbației șoselei r(t) ca funcție de timp, adică ele corespund pentru semnalul:

r(t)=a(1-cos8πt),

r(t)=0 în celelalte cazuri, unde coeficientul a=0.025 corespunde unei denivelări a șoselei de 5 cm.

Se observă că deși accelerația răspunsului în cazul primului model, la o denivelare de 5 cm este foarte bună, totuși abaterea(deplasarea) suspensiei este mai mare decât cea din cazul celui de-al doilea model. Aceasta se datorează faptului că abaterea suspensiei nu a fost penalizată în acest caz. Răspunsul abaterii suspensiei modelului doi(proiectarea 2, K2) la o denivelare bruscă de 5 cm este bun, dar răspunsul la denivelarea de 5 cm este mult inferior primului model. Aceasta se datorează(se observă încă o dată) faptului că, deplasarea și accelerația corpului mașinii nu au fost penalizate în proiectarea a doua, K2.

Figura VII.6 – Comparație între răspunsurile suspensiei pasive și a celor două modele de suspensie activă(K1 și K2)

Cele două proiectări, K1 și K2, reprezintă limitele extreme ale spectrului compromisului de performanță. Această secțiune a descris sinteza lui H∞ pentru a atinge obiectivele de performanță ale sistemului de suspensie activă. În mod egal, dacă nu chiar mai importantă, este proiectarea controlerelor robuste care modelează erorile sau incertitudinea.

În secțiunea următoare vom aborda un model de controler care atinge performanța robustă folosindu-se de metodologia sintezei μ.

VII.6 Proiectarea controlerului robust Kdk(K3) aferent suspensiei active a autovehiculelor utilizând toolbox-ul μ(iterația D-K)

Scopul oricărei proiectări a unui sistem automat de control/reglare este acela de a obține specificațiile de performanță dorite pe un model nominal, precum și pe alte procese care sunt apropiate de modelul nominal. Cu alte cuvinte, dorim să se obțină obiectivele de performanță în prezența erorilor de modelare , adică în prezența incertitudinilor. Aceasta se numește performanță robustă. După cum s-a mai spus și mai sus, în acest subcapitol se va proiecta un controler care obține performanță robustă, dar utilizând metodologia de proiectare a sistemelor automate de control din toolboxul „μ-Analysis & Synthesis Toolbox”(iterația D-K). Controlerul ce se va proiecta pentru suspensia activă a autovehiculului este aici notat cu Kdk(K3).

Controlerul H∞ al suspensiei active conceput în secțiunea precedentă nu a ținut cont de dinamica elementului hidraulic. În acest paragraf, în loc să fie introdus un element de execuție(actuator) perfect, se utilizează un model nominal de actuator care are erori de modulare, Hact(s), notat act(s) pentru ușurință.

Dacă în subcapitolele VII.3 și VII.4 dinamica elementului de execuție hidraulic(a actuatorului) a fost ignorată, la această a treia proiectare, dinamica elementului electronic hidraulic este un model de ordinul întâi, Hact(s). Acest model poate modela incertitudinea, adică poate ține seama da diferențele dintre modelul actuatorului și dinamica momentană a acestuia.

Hact(s)=act(s)=

Elementul hidraulic este modelat prin:

act=tf(1, [1/60 1]);

Modelul actuatorului însuși este incert. Putem descrie eroarea de modelare a elementului de execuție drept un set de modele posibile folosind o funcție de moderare. La frecvențe joase, sub 4 rad/sec, acesta poate varia cu până la 10% față de valoarea nominală. Începând de la 4 rad/sec procentul variației începe să crească și atinge 400%, la aproximativ 800 rad/sec. Incertitudinea de modelare este reprezentată de ponderea Wunc, care corespunde variației frecvenței modelului incert și obiectivului dinamic LTI, Δunc, definit ca unc.

Wunc=0.10*tf([1/4 1], [1/800 1]);

unc=ultidyn('unc', [1 1]);

actmod=act*(1+ Wunc*unc)

USS: 2 States, 1 Output, 1 Input, Continuous System

unc: 1×1 LTI, max. gain=1, 1 occurrence

Modelul actuatorului, numit aici actmod, este un sistem incert în spațiul stărilor. Următoarele diagrame Bode arată modelul nominal al actuatorului, notat cu semnul ”+”, și 50 de modele aleatoare descrise de actmod. Se utilizează comanda(funcția):

bode(actmod, 'b/r+', 50, logspace(-1, 3, 120))

Figura VII.7 – Diagramele Bode ale celor 50 de modele nominale aleatoare ale actuatorului

Modelul actuatorului incert, actmod, reprezintă modelul elementului hidraulic folosit pentru control. Diagrama modificată a interconexiunilor modelului(modelul generalizat P), necesară proiectării sistemelor automate de control robust, este:

Figura VII.8 – Schema bloc modificată a interconexiunilor modelului, necesară proiectării sitemelor automate de control robust(Kdk)

Dacă vom concepe controlerul pentru a oferi confort pasagerilor, ca și în prima proiectare de control H∞(K1), abaterea x1 a corpului autovehiculului va fi penalizată cu Wx1. Modelul incert ponderat al procesului necesar proiectării H∞, notat aici cu qcaricunc, utilizează aceeași comandă sysic. Așa cum s-a prezentat anterior, semnalul de control ce corespunde ultimei intrări la qcaric1, este fs. Accelerația corpului autovehiculului, care e zgomotoasă(conține zgomot), reprezintă semnalul măsurat și corespunde ultimei ieșiri a modelului qcaricunc.

systemnames='qcar Wn Wref Wact Wx1 actmod';

inputvar='[ d1; d2; fs ]';

outputvar='[ Wact; Wx1; qcar(2)+Wn ]';

input_to_actmod='[ fs ]';

input_to_qcar='[ Wref; fs]';

input_to_Wn='[ d2 ]';

input_to_Wref='[ d1 ]';

input_to_Wact='[ fs ]';

input_to_Wx1='[ qcar(1) ]';

qcaricunc=sysic;

Un controler care folosește sinteza μ(„μ-Analysis & Synthesis Toolbox”) este sintetizat folosind iterația D-K cu ajutorul comenzii dksyn. Procedeul iterației D-K este o aproximație a sintezei μ care încearcă să sintetizeze un controler care atinge performanțe robuste. Aici este o intrare de control, ncont, forța elementului hidraulic și un semnal de măsură, nmeas, accelerația corpului autovehiculului.

nmeas=1;

nant=1;

[Kdk, CLdk, gdk]=dksyn(qcaricunc, nmeas, ncont);

CLdkunc=lft(qcar([1:4 2], 1:2)*blkdiag(1, actmod), Kdk);

sprintf('sinteza mu(μ) a controlerului K3 a obținut o normă de

%2.5g', gdk)

ans=

sinteza mu(μ) a controlerului K3 a obținut o normă de 0.53946

Putem analiza controlerul obținut cu toolboxul de sinteza μ construind sistemul în buclă închisă, CLdkunc. Diagramele Bode ale sistemelor de suspensie pasivă și active atât ale modelului de actuator nominal, obținute cu controlerul din prima proiectare H∞(K1) cât și cu controlerul obținut prin sinteza μ(Kdk) sunt prezentate în figura VII.9. Se observă că, controlerul obținut prin sinteza μ atenuează mai bine primul mod rezonant în detrimentul scăderii performanței sub 3 rad/sec.

Figura VII.9 – Diagramele Bode ale sistemelor de suspensie pasivă și active(K1 și Kdk), cazul f.d.t. pe canalul perturbație-abaterea suspensiei

VII.7 Analiza în domeniul timp utilizând K1 și Kdk

Este important de înțeles cum ambele controlere robuste sunt în prezența erorii de modelare. Putem simula comportarea în timp a sistemului de suspensie activă utilizând primul controler H∞(K1) precum și controlerul obținut prin sinteză μ(Kdk). Sistemele în buclă închisă cu incertitudini, CL1unc și CLdkunc, sunt formate cu ajutorul controlerelor K1 și respectiv Kdk. Pentru fiecare sistem incert sunt simulate 40 de modele aleatoare din jurul acelui model nominal. După cum se observă, controlerele sunt robuste și au performanțe bune în prezența erorii de modelare a actuatorului(vezi răspunsurile sistemelor în buclă închisă, CL1unc și CLdkunc, în figura VII.10 și respectiv VII.11).

Controlerul de sinteză μ , Kdk, atinge performanțe ușor mai bune decât primul controler H∞(K1).

CL1unc=lft(qcar([1:4 2], 1:2)*blkdiag(1, actmod), K1);

[CLdkunc40, dksamples]=usample(CLdkunc, 40);

CL1unc40=usubs(CL1unc, dksamples);

Figura VII.10 – Răspunsurile celor 40 de modele aleatoare ale procesului, utilizând controlerul K1(sinteză hinfsyn)

Figura VII.11 – Răspunsurile celor 40 de modele aleatoare ale procesului, utilizând controlerul Kdk(sinteză μ-DK)

Cap.VIII Concluzii

Sistemul de suspensie pentru autovehicule este un sistem dificil de studiat și proiectat dar, totodată unul interesant, funcționarea acestuia realizându-se în timp real.

Sunt foarte greu de implementat forțe și frecvențe ale corpului autovehiculului fără a utiliza controlere complexe sau modificări mecanice asupra acestuia.

Din acest motiv, lucrarea de față se ocupă nu numai cu modelarea unui astfel de sistem de suspensie, dar și făcându-se sinteza unor regulatoare continue PI, PID, optimale LQR, robuste, precum și proiectarea unui regulator digital(prin reproiectarea unui regulator continuu).

Noua tehnologie oferă prin suspensiile active o soluție simplă și economică pentru îmbunătățirea controlul vehiculului cât și a performanțelor acestuia, fără a diminua confortul, așa cum se întâmplă în cazul sistemelor de suspensie convenționale. Aceasta se realizează în principal prin asigurarea unui contact permanent între pneu și drum.

Din lucrare reiese că sinteza robustă, deși este mai bună pentru obținerea performanțelor impuse inițial, totuși din cauza complexității acestuia, în practică, se preferă , probabil, utilizarea unui regulator PI sau PID mai simplu și mai ieftin.

Bibliografie

[1] BERNSTEIN, D. S.-HADDAD, W.M.: LQG Controll with an H∞ Performance Bound,IEEE Trans. on Automatic Control

[2] CROLLA, D. A.-ABDEL HADY, M. B. A.: Activa suspension Control; Performance Comparisons Using Control Laws Applied to a Full Vehicle Model, Vehicle System Dynamics

[3] DOYLE, J.-ZHOU, K.-GLOVER, K.-BODENHEIMER, B.: Mixed H2 and H∞ Performance Objectives. II. Optimal Control, IEEE Trans. on Automatic Contro,pp.831-847, 1989

[4] HROVAT, D.: Optimal Active Suspension Structures for Quarter-car Vehicle Models, Automatica

[5] MICHELBERGER, P.-BOKOR, J.-PALKOVICS, L.: Robust Design of Active Suspension System, International Journal of Vehicle Design

[6] YAMASHITA, M.-FUJIMORI, K.-HAYAKAWA, K.-KIMURA, H.: Applications of H∞ Control to Active Suspension System, Automatica

[7] LIN, J.-KANELLAKOPOULOS, I.: Road Adaptive Nonlinear Design of Active Suspension, Proc of the ACC., pp 714-718, 1997

[8] HEDRICK, J. T.-BATSNEN, T.: Invariant Properties of Automotive Suspensions, Proc. Of IME., pp.21-27, 1990

[9] FLORESCU, A.: Lucrare de licență, Universitatea „Dunarea de jos”, iunie 2004

[10] MATHWORKS, Mu-Analysis and Synthesis, Toolbox User’s Guide, Natick, MA, USA, pp. 734, 2001

ANEXĂ

PROGRAM 1

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

step(A, B, C, D, 1)

PROGRAM 2

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

[ num1, den1] = ss2tf ( A, 0.1*B, C, D, 2)

PROGRAM 3

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

nump=[nump(3) nump(4) nump(5)]

PROGRAM 4

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

[ num1, den1] = ss2tf ( A, 0.1*B, C, D, 2)

PROGRAM 5

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

num1 = [ num1(2) num1(3) num1(4) num1(5)]

PROGRAM 6

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

%Determinarea functiei de transfer a procesului

nump=[(m1+m2) d2 k2];

denp=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2];

'G1(s)'

printsys (nump, denp)

%Raspunsul sistemului în bucla deschisa la intrare treapta

step (nump, denp)

title ('Raspunsul sistemului in bucla deschisa la intrare treapta')

PROGRAM 7

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

%Determinarea functiei de transfer pentru perturbatie

num1=[-(m1*d2) -(m1*k2) 0 0];

den1=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2];

'G2(s)'

printsys (0.1*num1,den1).

PROGRAM 8

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

nump=[(m1+m2) d2 k2];

denp=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1)

step (nump, denp)

PROGRAM 9

m1=290;

m2=59;

k1 = 16182;

k2 = 190000;

d1 =1000;

d2 = 0;

A = [0 1 0 0

-(d1*d2)/(m1*m2) 0 ((d1/m1)*((d1/m1)+(d1/m2)+(d2/m2)))-(k1/m1) -(d1/m1)

d2/m2 0 -((d1/m1)+(d1/m2)+(d2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B = [0 0

1/m1 (d1*d2)/(m1*m2)

0 -(d2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

num1=[-(m1*d2) -(m1*k2) 0 0];

den1=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1)

step(0.1*num1,den1)

axis ([0 3.5 -0.1 0.1])

PROGRAM 10

m1=290;

m2=59;

k1=16182;

k2=190000;

d1=1000;

d2=0;

nump=[(m1+m2) d2 k2]

denp=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2]

num1=[-(m1*d2) -(m1*k2) 0 0]

den1=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1) k1*k2]

numf=num1;

denf=nump;

R = roots(denp)

rlocus( nump,denp)

z1 = 3+3.5i;

z2 = 3-3.5i;

p1 = 30;

p2 = 60;

numc = conv([1 z1],[1 z2]);

denc = conv([1 p1],[1 p2]);

rlocus (conv(nump,numc),conv(denp,denc))

axis ([-40 10 -30 30])

[K, poles] = rlocfind (conv (nump,numc) , conv (denp , denc))

PROGRAM 11

m1=290;

m2=59;

k1=16182;

k2=190000;

d1=1000;

d2=0;

nump=[(m1+m2) d2 k2]

denp=[(m1*m2) (m1*(d1+d2))+(m2*d1) (m1*(k1+k2))+(m2*k1)+(d1*d2) (d1*k2)+(d2*k1)

w=logspace(-1,2);

bode(nump,denp,w)

nump=100000*nump

numc = conv([1.13426 1], [1.13426 1]);

denc = conv([0.035265 1], [0.035265 1]);

margin(conv(nump,numc),conv(denp,denc))

w = logspace(-4,4);

[mag,phase,w] = bode(conv(nump,numc),conv(denp,denc),w);

subplot (2,1,1);

semilogx (w, 20*log10(mag));

grid

title ('Caracteristicile Bode când se aduga un regulator cu dublu avans de faza')

xlabel ('Frecventa (rad/s)')

ylabel ('20logM')

subplot (2,1,2);

semilogx (w,phase)

axis ([1e-4, 1e4, -180, 360])

grid

xlabel ('Frecventa (rad/s)')

ylabel ('Faza (grade)')

numa = conv(numf,denc);

e1=conv(denp,denc)

e2=conv(nump,numc)

e3=[e2 0 0]

dena = e1+e3

t=0:0.01:5;

step(0.1*numa,dena,t)

axis([0 5 -.01 .01])

PROGRAM 12

ms=290;

mu=59;

kr=16182;

br=1000;

kt=190000;

g=9.8;

A=[0 1 0 0;-kr/ms -br/ms kr/ms br/ms;0 0 0 1;kr/mu br/mu -(kr/mu+kt/mu) -br/mu]

B=[0;1/ms;0;-1/mu]

C=[1 0 0 0]

D=[0]

p=[0;0;0;kt/mu]

comp_lib=[0;-g;0;-g]

t=0:0.1:5

f=ones(size(t));

sys=ss(A,[p comp_lib],C,D)

zr=0.5*sin(6*t);

[z,t]=lsim(sys,[zr;f],t)

plot(z)

grid

title('Raspunsul sistemului pentru denivelari sinusoidale de forma zr=0.5*sin(w*t) unde w=6rad/sec ')

xlabel('timpul(sec)')

ylabel('amplitudinea[m]')

PROGRAM 13

ms=290;

mu=59;

kr=16182;

br=1000;

kt=190000;

g=9.8;

A=[0 1 0 0;-kr/ms -br/ms kr/ms br/ms;0 0 0 1;kr/mu br/mu -(kr/mu+kt/mu) -br/mu]

B=[0;1/ms;0;-1/mu]

C=[1 0 0 0]

D=[0]

p=[0;0;0;kt/mu]

comp_lib=[0;-g;0;-g]

t=0:0.1:5

f=ones(size(t));

sys=ss(A,[p comp_lib],C,D)

zr=0.5*sin(8*t);

[z,t]=lsim(sys,[zr;f],t)

plot(z)

grid

title('Raspunsul sistemului pentru denivelari sinusoidale de forma zr=0.5*sin(w*t) unde w=8rad/sec ')

xlabel('timpul(sec)')

ylabel('amplitudinea[m]')

PROGRAM 14

ms=290;

mu=59;

kr=16182;

br=1000;

kt=190000;

g=9.8;

A=[0 1 0 0;-kr/ms -br/ms kr/ms br/ms;0 0 0 1;kr/mu br/mu -(kr/mu+kt/mu) -br/mu]

B=[0;1/ms;0;-1/mu]

C=[1 0 0 0]

D=[0]

p=[0;0;0;kt/mu]

comp_lib=[0;-g;0;-g]

t=0:0.1:5

f=ones(size(t));

sys=ss(A,[p comp_lib],C,D)

zr=0.5*sin(58*t);

[z,t]=lsim(sys,[zr;f],t)

plot(z)

grid

title('Raspunsul sistemului pentru denivelari sinusoidale de forma zr=0.5*sin(w*t) unde w=58rad/sec ')

xlabel('timpul(sec)')

ylabel('amplitudinea[m]')

PROGRAM 15

ms=290

mus=59

bs=1000

ks=16182

kt=190000

A12=[0 1 0 0;[-ks -bs ks bs]/ms]

A34=[0 1 0 0;[ks bs -ks-kt -bs]/mus]

B12=[0 0;0 10000/ms]

B34=[0 0;[kt -10000]/ms]

C=[1 0 0 0;A12(2,:);1 0 -1 0;0 0 0 0]

D=[0 0;B12(2,:);0 0;0 1]

qcar=ss([A12;A34],[B12;B34],C,D)

wn=0.01

wref=0.07

wact=(100/13*tf([1 50],[1 500]))

systemnames='qcar wn wrefwact wx1'

wx1=8*tf(2*pi*5,[1 2*pi*5])

bode(qcar)