Matematică și Informatică Aplicate în Inginerie [606011]
1
UNIVERSITATEA POLITEHNICA DIN BUCUREȘTI
FACULTATEA DE ȘTIINȚE APLICATE
Matematică și Informatică Aplicate în Inginerie
Proiect de Licență
Coordonator științific:
Lect.dr. Tiberiu VASILACHE
Absolvent: [anonimizat]
2017
2
UNIVERSITATEA POLITEHNICA DIN BUCUREȘTI
FACULTATEA DE ȘTIINȚE APLICATE
Matematică și Informatică Aplicate în Inginerie
Aprobat
Decan:
Prof. dr. Emil Petrescu
Sisteme optimale și metode rezolvare a
ecuației Riccati
Coordonator științific: Absolvent: [anonimizat].dr. Tiberiu VASILACHE Ramona -Maria DUMITRESCU
București
2017
3
CUPRINS
1. INTRODUCERE ………………………….. ………………………….. ………………………….. ………………….. 4
2. CONTROL OPTIMAL ………………………….. ………………………….. ………………………….. ………… 6
2.1 Noțiuni de calcul variațional ………………………….. ………………………….. ………………………….. . 6
2.2 Formularea problemei fundamentale în teoria controlului optimal ………………………….. …. 12
3. OPTIMIZAREA CU CRITE RIU PĂTRATIC A SISTE MELOR LINIARE ……………… 16
3.1 Problema de optimizare liniar pătratică ………………………….. ………………………….. ………….. 16
3.2 Soluția problemei liniar pătratice cu timp final finit (EDMR) ………………………….. ……….. 20
3.3 Problema reglării optimale a mărimii de ie șire ………………………….. ………………………….. … 25
3.4 Regulatorul liniar pătratic ………………………….. ………………………….. ………………………….. … 27
4. ECUAȚIA MARICEALĂ RI CCATI CU METODE DE R EZOLVARE ……………………. 29
4.1 Ecuația matriceală algebrică Riccati (EMAR) ………………………….. ………………………….. …. 29
4.2 Metode de rezolvare a ecuației matriceale algebrice Riccati ………………………….. ………….. 37
4.2.1 Metoda Potter ………………………….. ………………………….. ………………………….. …………. 37
4.2.2 Metoda Schur ………………………….. ………………………….. ………………………….. …………. 41
4.2.3 Metoda matricei de semn ………………………….. ………………………….. ……………………… 46
5. APLICAȚII ………………………….. ………………………….. ………………………….. ……………………….. 48
5.1 Aplicații teoretice ale controlului optimal ………………………….. ………………………….. ………. 48
5.2 Sistemul de suspensie la autobuz ………………………….. ………………………….. …………………… 59
5.2.1 Noțiuni introductive pentru reprezentarea în spațiu l stărilor ………………………….. …… 59
5.2.2 Compararea metodelor de rezolvare a ecuației Riccati ………………………….. ………….. 62
6. CONCLUZII ………………………….. ………………………….. ………………………….. ……………………… 67
7. ANEXE ………………………….. ………………………….. ………………………….. ………………………….. ….. 68
8. BIBLIOGRAFIE ………………………….. ………………………….. ………………………….. ………………… 70
4
1.Introducere
Controlul optimal devine cunoscut după al doilea război mondial odată cu apariția
războiului rece, când matematicienii au fost implicați în noua confruntare dintre SUA și URSS.
O provocare matematică a celor vremuri cu care s -au confruntat a fost: Care est e traiectoria
optimală a unui avion ce urmează a fi direcționat dintr -o poziție de croazieră dată, într -o poziție
favorabilă necesară răspunderii unui posibil atac?
Răspunsul la acestă problemă a fost ușor de sesizat. Matematicienii au recunoscut imediat
asemănarea cu întrebarea lui John Bernoulli din 1969: Care este curba celei mai rapide coborâri
între două puncte date într -un plan vertical?, care și -a imaginat în acest fel problema
brachistochronei (în grecește, brachistos=cel mai scurt, chrones=timp).
În ultimele decade domeniile teoriei controlului și optimizării au fost la interfața dintre
creativitatea matematică, inginerie si informatică. Scopul teoriei controlului este să ajute la
îınțelegerea principiilor fundamentale ale controlului și să le carac terizeze matematic într -un
mod ce poate fi folosit pentru obținerea exactă a unor controale ce ne ajută să îndeplinim un scop
precizat. De fapt, scopul optimizării performanțelor ne duce către controlul optimal ce ilustrează
afinitatea dintre teoria contro lului, calculul variațional și optimizare.
Optimizarea poate fi de finită ca știința determinării celei mai bune soluții la anumite
probleme definite matematic, care sunt adesea modele ale realității fizice.
A face cel mai bine posibil este sensul oricăr ei atitudini naturale în viața de zi cu zi.
Pentru un inginer a face „ cel mai bine posibil ” ar trebui să fie un obiectiv permanent
atunci când are de conceput o clădire, de dimensionat o instalație etc.
În aproape orice activitate în care sunt prelucrate i nformațiile numerice: știință, inginerie,
matematică, economie etc. metodele de optimizare au o largă aplicabilitate. O alegere a
domeniilor în care apar probleme de optimizare ar cuprinde: proiectarea reactoarelor chimice, a
aparatelor aerospațiale, a clă dirilor, a podurilor, în diferite ramuri ale analizei numerice, principii
în sisteme de ecuații diferențiale și cu derivate parțiale etc.
Lucrarea de față își propune să prezinte controlul optimal aflat în strânsă legatură cu
teoria sistemelor și modul în care ecuația matriceală Riccati poate ajuta în găsirea unei comenzi
optimale ce minimizează o funcțională asociată unei reprezentari în spațiul stărilor a unui sistem
fizic. Conținutul lucrării se axează în principal pe rezolvarea problemei liniar pătratic e și
5
aparițiile acesteia în diverse domenii, mai multă atenție a fost acordată ecuației matriceale
algebrice Riccati.
În acest sens al 2 -lea capitol abordează aspectele teoretice cu privire la conceptul general
de control optimal. De asemenea, sunt ilustar te și noțiunile de calcul variațional care au o temelie
esențială asupra controlului optimal și condițiile necesare de optimalitate.
În capitolul al 3 -lea este realizată o prezentare generală a o ptimizării cu criteriu pătratic a
sistemelor liniare, făcându -se referie la problema de optimizare pătratică, la soluția problemei
liniar pătratice cu timp final finit, p roblema reglării optimale a mărimii de ie șire si a r egulatorului
liniar pătratic.
Capitolul al 4 -lea intitulat “ Ecuația matriceală Riccati cu metode de rezolvare ” se
centrează pe ecuația matriceală algebrică Riccati și prezintă trei metode de rezolvare ale acesteia.
Metodele de rezolvare folosite pe p arcursul acestei lucrări sunt: metoda Potter, metoda Schur și
metoda maticei de semn .
Capitolul al 5 -lea conține prezentarea în întregime a aplicațiilor. Sunt prezentate exemple
simple ce pot fi tratate ușor prin metode analitice, dar și cu ajutorul programului MATLAB
utilizând comenzile specific e din Control System Toolbox . Această lucrare prezintă ca și aplicație
centrală sistemul de suspensie la autobu z, unde sunt aplicate metodele ecuației Riccati pentru a
verfica faptul că, acest sistem este unul cu aplicabilitate în realitate, și faptul că se poate ilustra
oscilația pe care oamenii care stau în autobuz o resimt sau nu în timpul deplasării pe un drum cu
perturbații.
6
2. Control Optimal
2.1 Noțiuni de calcul variațional
Se dă sistemul Σ descris de ecuația de stare:
),(xtfx
,
mR x tx 0 0)( (2.1)
Aici,
)()(
)(1
txtx
tx
n este vectorul de stare și
f o funcție vectorială de clasă .
Căutăm construirea unui sistem Hamilton canonic, pentru găsirea extremalelor
funcționalei
2
1),,,,,,( ],,[][1 1 1t
tn n n dtx xx xtL x xIxI
(2.2)
în mulțimea
},1, )(, )(|]),([ )),(( ),,({2 1 210
212
1 n ib txa txttC ttC x xx Mi i i i n , unde
abi nii,, ,1
sunt numere reale date. Integrând prin părți termenii ținând seama de
proprietățile funcțiilor și aplicând lema fundamentală, deducem sistemul de ecuații diferențiale
ordinare de ordinul al doilea
0. ………. ………. ……….00
2 21 1
i i xL
dtd
xLxL
dtd
xLxL
dtd
xL ,
n i,1 . (2.3)
care este un sistem de n ecuații diferențiale de ordinul al doilea cu n funcții necunoscut
x xn 1,,
și care mai poate fi scris și sub forma
7
2
12
12
0 1L
xxxL
xxxL
xtL
xi n
i jj
jn
i jj
jn
i i
, ,. (2.4)
Ecuațiile sistemului (2.3) se numesc ecuații Euler –Lagrange, iar orice soluții ale
acestui sistem se numesc extremale ale funcționalei (2.2) .
Acest sistem poate fi rezolvat în raport cu
x xn 1, dacă Hessianul este nenul, adică
0 det
,12
njij ixxL
. (2.5)
Luând ca variabile canonice (hamiltoniene)
pL
xi
i
și
xi ni, ,1 , iar
x xn 1,, ca
variabile, putem scrie sistemul
),,,,,,(1 1 n n
ii x xx xtxLp
(2.6)
Întrucât (2.5) este îndeplinită, teorema funcțiilor implicite afirm ă că soluția (locală) a
sistemului (2.6) există și are forma
x ftx xp p j j n n (,,,,,,) 1 1
(2.7)
Cu alte cuvinte, Hamiltonianul funcționalei (2.2) este funcția
Htx xp p xL
xLn n j
j jn
(,,,,,,)1 1
1
, (2.8)
adică
H xp Ljj
jn
1
(2.9)
unde
xj sunt date de ( 2.6).
8
Reconsiderând Lagrangianul ca fiind
L xp Htx xp pjj
jn
n n
11 1 (,,,,,,)
(2.10)
rescriem sistemul Euler –Lagrange corespunzător:
0 ) (0 ) (
i ii i
pL
pL
dtdxL
xL
dtd
d
dtpH
x
d
dtxH
pi
i
i
i()
()
0
0 0 .
Astfel, se obț ine sistemul Hamilton canonic
iiii
pHppHx
,
care înlocuiește sistemul Euler –Lagrange (2.3) și care este un sistem de 2 n ecuații diferențiale de
ordinul I cu 2 n funcții necunoscute
xp i nii,, ,1 . Acesta mai poate fi scris și în forma
pxHJpx
dtd
, unde
n nn n
O II OJ .
Exemplul 2.1.
Considerăm funcționala:
02
22
1 2 1 2 1 ]))(( ))((2)()(2[(],[ dttx tx txtx xxI
,
pentru care dorim să obținem sistemul Euler -Lagrange și sistemul Hamilton canonic.
Sistemul Euler -Lagrange devine:
9
00
2 21 1
xL
dtd
xLxL
dtd
xL
0 2 20 2 4 2
0 2 20 2 4 2
2 11 1 2
2 11 1 2
x xx x x
xdtdxxdtdx x
Introducem momentul generalizat
),(2 1pp p
1 2 1 2 1 1 2),,,,( x xxxxtp
2 2 1 2 1 2 2),,,,( x xxxxtp
și Lagran gianul care trebuie să fie regulat, adică
0 det
,12
njij ixxL
De aici rezultă că:
042 002det
L este regul at
Atunci
2 21 1
2
221
11
2121
22
p xp x
xxLpxxLp
,
deci Hamiltonianul corespunzător este:
2
22
12
1 212
12
12
1 212
22
1 22 1141
412 2 )41
412 2(21
21p p x xx p p x xx p p L pxpx H
Obținem sistemul Hamilton canonic:
10
01 22 4 r r
1
221 2
112
221
11
24 22121
xxHpx xxHpppHxppHx
În continuare căutăm soluția generală pentru
1x și pentru
2x . Pornim de la sistemul Euler –
Lagrange obținut deja
0 2 20 2 4 2
2 11 1 2
x xx x x
Împărțim prin 2 și rescriem sistemul:
00 2
2 11 1 2
xxxx x
Din a doua ecuație rezultă următoarele:
ivx x x x xx2 1 2 1 2 1 0 și astfel
prima ecuație dev ine:
0 22 2 2 ivxx x .
Pornind de la acestă ecuație încercăm să obținem un sistem fundamental de soluți
0 22 2 2 ivxx x
4
2 r xiv
,
2
2r x ,
12r x
Notăm
t r2 .
Avem solu țiile:
t tet ei rt tet eir
tttt
coscossinsin
21
2 , 1 010)1 (0)1( 012
2,12 22 22 2
ami r r rrt t t
11
Pentru
0 și
1 obținem sistemul fundamental de soluții :
Scriem soluțiile generale ale funcționalei:
ttCttCt Ct Ctx cos sin cos sin )(4 3 2 1 2
ttCt CttCt Ct Ct Ctx sin cos cos sin sin cos )(4 4 3 3 2 1 2
ttCttCt C Ct C C tx cos sin cos)2 ( sin)2 ()(4 3 3 2 4 1 1
unde,
. ,,,4 3 2 1 R CCCC
În programul MATLAB pentru a rezolva acest exemplu trebuie să avem în vedere
Lagrangianul problemei și dimensiunea acestuia astfel , folosi m comenzi le specifice pachetul ui
Control System Toolbox , a programului . În cazul nostru Lagrangianul este de dimensiune 2×2 .
ttCttCt C C t C Ct Ct Ct CttCt Ct Ct Ct C tx
cos sin cos)2 ( sin)2 (cos sin sin sin cos cos cos sin )(
4 3 3 2 4 14 4 4 3 3 3 2 1 2
t ttrtttr
coscossinsin
21
12
2.2 Formularea problemei fundamentale în teoria controlului optimal
Considerăm dat un proces evolutiv, a cărui stare este descrisă printr -un numă r de
parametri de stare
ix , ansamblul lor constituind un punct figurativ
x (sau variabilă de stare), al
unei varietăți
X , spațiul stărilor (în particular domeniu din
nR ), ca și printr -un număr de
parametrii de comandă
au , ansamblul lor constituind un punct figurativ numit variabilă de
control sau de comandă și aparținând unei varietăți
U (domeniu din
mR ). O dată cu v arietățile
X
și
Uvom considera și varietățile lor tangente și cotangente.
Considerăm sistemul:
(,,), () xftxu xt x 1 1
(2.11)
pentru
tItt[,]12 și funcția
nUXIf R : continuă și cu derivate continue pe
porțiuni
))(,),(),(,),(),(,( :1 2 1 tu tutx txtxtf tfl n i i ,
n i ,,2,1 .
Ea se identifică cu evoluția unui sistem dinamic, descrisă prin ecuația:
)(),(,)(tutxtfdttdx
,
10, t tt (2.12)
numită sistem de control.
Numind în general o funcție
UIu: funcție de control sau de comandă , care
introdusă în sistemul (2.12) îl transformă pe acesta într -un sistem dinamic obișnuit pentru care o
soluție a
XIx: este o funcție de stare. Vom spune că sistemul (2.12) guvernează un proces de
evoluție.
Sistemul (2.12) se numește sistem controlat , iar perechea de funcț ii
)( : txtx ,
)( : tutu
(definte pe
I ) care verifică sistemul (2.11), poartă respectiv numele de traiectorie
admisibilă și de control admisibil . Vom spune că sistemul (2.11) guvernează procesul evolutiv,
t
reprezintă timpul și
dtdxx este viteza de variație a parametrului de stare
x .
Dacă funcția
f nu depinde explicit de timp, atunci procesul evolutiv descris de sistemul
(2.11) se numeș te staționar sau autonom , în caz contrar, este nestaționar sau neautonom .
13
Notăm cu
T spațiul traiectorilor admisibile, U spațiul comenzilor admisibile și X spațiul
stărilor.
Alegem funcționala de cost dată de ecuația :
Ju LtxtutdtMtxttxt
tt
() (,(),()) (,(),,())
12
1 12 2
, (2.13)
unde primul termen dat de integrală este criteriul de tip Lagrange, iar cel de -al doilea termen
este criteriul de tip Mayer. O astfel de funcțională de cost, conținând ambele criterii este nu mită
criteriu de tip Bolza (L este Lagrangianul, iar M este Mayerianul).
Vom presupune că,
R UXIL: este o funcție de clasă C1 pozitivă, iar
R R R 2 1 : X X M
o funcție convexă de clasă C1, unde
X XX X 1 2 , .
În contextul următo r, problema controlului optimal este determinarea unei comenzi
~()u
U care să realizeze minimul funcționalei de cost J(u), astfel încât
Ju Ju () () ~ ,
)(u
U , unde,
Ju Ltxtutdt Mtxttxt
tt
(~) (,~(),~()) (,~(),,~())
12
1 12 2
și
)(~tx este soluția sistemului (2.11) ce corespunde comenzii optimale
)(~tu .
Având o traiectoria optimală
(~(),~()) xu și o variație mică a sa
((),()) xu T de forma:
)( )(~)()( )(~)(
tk tututh txtx
(2.14)
cu R și h(t), k(t) funcții vectoriale arbitrare fixate de dimensiuni respectiv n și m. Costul
corespunzător poate fi privit ca funcție de :
.))( )(~,),( )(~,( ))( )(~),( )(~,( )~( )( )(
2 2 2 1 1 12
1
th txtth txtMdttk tuth txtL k uJ uJ Ft
t
14
Astfel,
(~(),~()) xu este o traiector ie optimală,
Ju F (~) () 0 , deci = 0 este un punct de
minim pentru F(). Prin urmare, conform teoremei lui Fermat rezultă că,
F()0 0 . (2.15)
Derivăm funcția
()F :
)( )( )( )( )(2 2 1 12
1th
xMth
xMdttkuLthxLFT T t
tT T
, (2.16)
unde derivatele parțiale sunt calculate pentru
(~() (),~() ()) xt htut kt .
Conform relației (2.15) rezultă că, traiectoria optimală verifică ecuația :
0)( )( )( )(2 2 1 12
1
th
xMth
xMdttkuLthxLT T t
tT T
(2.17)
Formăm Hamilt onianul problemei de control optimal (2.11), (2.13) prin
R RRRR n m nH:
:
),,()( ),,( )()( ),,( ),,,( uxtftp uxtLtxtp uxtL puxtHT T
,
unde
] [1 nTp p p este vectorul multiplicatorilor Lagrange; se numește și vectorul adjunct al
lui x.
Atunci
Ltxu HtxupptftxuT(,,) (,,,)()(,,) și ecuația (2.16) devine:
FH
xhtH
uktdt ptf
xhtf
uktdt
M
xhtM
xhtT T
tt
T
tt
T T() () () () () ()
() ().0
12
12
1122
Dacă notăm
htdxktdu () ,() , avem
f
xhtf
uktdf dxht () () ()
.
Integrând prin părți cea de -a doua integrală din (2.16) obținem
15
pthtdt ptht pthtdtT
tt
T
tt T
tt
()() ()()()()
12
12
12
deci, pentru
(~,~)xu avem
0)( )( )( )( )( )( )0(2 221 112
1
th tp
xMth tp
xMdttkuHthpxHFT T t
tT T
pentru orice h(t), k(t) (suficient de mici).
Având la bază rezultatele obținute mai sus, enunțăm teorema de optimalitate:
Teorema 2.2.1 Dacă
(~(),~()) xu este o traiectorie optimală a sistemului (2.11) cu funcționala de
cost (2.13), atunci există multiplicatorii Lagrange
p pn 1 , astfel încât
(~(),~()) xu verifică
următoarele condiții necesare:
i. sistemul canonic:
),~,~,(~puxtpHx
(2.18)
),~,~,( puxtxHp
(2.19)
ii. ecuația de control optimal:
H
utxup (,~,~,)0 ; (2.20)
iii. condițiile de transversalitate:
0)()( )()(2 2 2 1 1 1
th tp
xMth tp
xMT T
. (2.21)
16
3. Optimizarea cu criteriu pătratic a sistemelor liniare
3.1 Problema de optimizare liniar pătratică
Considerăm sistemul liniar (variabil în timp) :
() ()()()() xtAtxtBtut
(2.22)
cu condiția ințială fixată
xt x ()1 1 și funcționala de cost pătratică
)()(21)()()( )()()(21)(2 22
1tMxtx dttutRtutxtQtx uJTt
tT T
, (2.23)
unde
o
ACInn(, ) R ,
BCInm(, ) R matrice continue pe
Itt[,] 12 .
o
nnMR este o matrice reală, nenegativ definită (
0 , = M M MT );
o
QCInn(, ) R ,
RCImm(, ) R și pentru orice
],[21ttt ,
0 )(tQ
(matrice reală, simetrică, semipozitiv definită
n TT
Rx QxxtQ tQ
,0)( )( , iar
0 )(tR (matrice
reală, simetrică, pozitiv definită
0 ,0)( )(
u RxutR tR
TT .
Dorim determinarea unei comenzii optimale
)(u în circuit închis pe intervalul
],[21tt
care să minimizeze funcționala de cost (2.19).
Ținem cont și de urmatoarele observații, referitoare la matricele funcțion alei de cost și
anume:
– Q,R și M sunt dependente de
t, atunci când se dorește o ponderare în timp a
influențelor componentelor respective în diverse perioade ale intervalului de optimizare
],[21tt
;
– R nu poate fi nulă, c ăci indicele de calitate nu ar mai depinde de comandă
17
– una din matricele M sau Q poate fi nulă, dar nu ambele, căci problema de
optimizare ar admite în acest caz soluția trivială
0)(tu .
Aplicăm teorema de optimalitate prin etapele următoare :
Etapa 1. Scriem Hamiltonianul:
Htxup xQx uRu pAx BuT T T(,,,) [ ] ( ) 1
2
(2.24)
Etapa 2. Scriem sistemul canonic (2.18), (2.19):
Bu AxpHx
(2.25)
pA Qx pA tQx xtQxHpTTQQ
T T T
) ))( (21)(21(
(2.26)
Etapa 3. Determinăm ecuația comenzii optimale
)(~tu din (2.20) :
0 )( 0 pButRuHT
(2.27)
înmulțind la stânga cu
R1 (există, deoarece R >0), rezultă comanda optimală:
pBR tuT1)(~
(2.28)
Prin intermediul matricei Hessiane verificăm dacă
)(~tu realizează minimizarea:
2
20H
uRt() ,
cu alte cuvinte,
)(~tu realizează minimul problemei și, conform (2.28) este unică.
Etapa 4. Înlocuim (2.28)
)(~tu în ecuația (2.25),(2.26) și obținem tripletul
))(~),(~),(~( tptutx
, care verifică sistemul de ecuații diferențiale:
ptBtRtBxtAxT)()()( )(1
(2.29)
ptAxtQ pT)( )(
(2.30)
18
sau
x
pHx
p
, unde matricea hamiltoniană
H are forma
)( )()()( )(1
tA tQtBt BR tAHTT (2.31)
Etapa 5. Scriem condițiile de transversalitate, cu stare inițiala fixată și stare finală
liberă:
xt x ()1 1
)( )(2
22 tMxxMtp
.
Astfel, sistemul (9.29), (9.30 ) poate fi rezolvat pentru a obține
))(~),(~( tptx .
Exemplul 3 .1
Considerând indicele de performanță :
1
2)(212t
ttu J
cu sistemul dublu integrator
)( )()( )(
22 1
tutxtutx
găsiți controlul optimal și stările optimale, având condițiile:
1)0(1 x
,
2)0(2 x ,
1)2(1 x ,
0)2(2 x .
Soluție:
Formăm Hamiltonianul:
)()( )()( )(21
2 2 12tut txt tu H
;
19
t tu t tuuH);(~)(~0)(~)(~02 2
)(~
21)(~)(~)(~)(~)(~)(~
21))(~),(~),(~),(~(
2
2 2 12
2 2 12
2 2 1 2 1
t txtt txt t t t txtxH
Obținem ecuațiile:
)(~)(~)(~)(~
2 22 1
t txtxtx
)(~)(~0)(~
1 21
t tt
.
Rezolvăm ecuațiile de mai sus și obținem:
2 423
21 224 33
1
2)(~2 6)(~
CtCtCtxCtCtCtCtx
4 3 23 1
)(~)(~
CtC tCt
.
Obținem controlul optimal de forma :
4 3 2)(~)(~CtCt tu
Folosindu -ne de comenzi le programului MATLAB specifice pachetului Control System
Toolbox încercăm să găsim soluția problemei cu un grafic sugestiv care să ne indice controlul
optimal și stările.
20
3.2 Soluția problemei liniar pătratice cu timp final finit (EDMR)
Considerăm legea de comandă optimală
)()( )( txtFtu , cu
)(tF o m× n matrice continuă
pe
].,[21tt Astfel, căutăm o matrice P(t) unde vectorul adjunct p să fie de forma
ptPtxt () ()() ,
cu condiția finală
M tP)(2 .
Sistemul de ecuații diferențiale (2.29) și (2.30) devin:
xPB BRAxT1
.
PxA Qx xPxPpT
Înlocuim
x și obținem:
PAPPA QPBR BPxT T 10
, pentru orice
xX .
Rezultă că P trebuie să fie soluție a ecuației diferențiale matriceale Riccati (EDMR)
PAPPA QPBR BPT T 10
(2.32)
cu condiția finală
M tP)(2
, deoarece
)( )()( )(2 2 2 2 tMx txtP tp . (2.33)
Considerând i mplementată legea de comandă optimală (Fig. 1), se constată imediat că
aceasta furnizează, corespunzător inițializării arbitrare
1x , o soluție în circuit deschis a problemei
liniar pătratice. Altfel spus, în timp ce determinarea unei so luții în circuit deschis a problemei
liniar pătratice trebuie făcută ori de câte ori se schimbă inițializarea
1x , cunoașterea unei soluții
în circuit închis, explicitată prin matricea F(t), permite auto generarea, prin schema din Fig. 1, a
unei soluții în circuit deschis, oricare ar fi inițializarea
1x .
(Fig.1)
1 1)( t tx
)()( )()( )( tutBtxtAtx
)(tx
)(tu
)(tF
21
În consecință, am demonstrat că
1Tu R B Px
este comanda optimală, iar în acest caz
matricea P(t) este soluția ecuației Riccati (2.32) cu condiția finală (2.33).
Presupunem în continuare că P(t) este soluție a ecuației (2.32) cu condiția finală (2.33). și
arătăm că
1Tu R B Px
realizează minimul funcționalei. Astfel că încercăm demonstrația în
sens invers.
Remarcăm că (2.32) este o ecuație simetrică, deci soluția sa îndeplinește
PPT .
Calculăm:
d
dtxtPtxt xPx xPxxPxT T T T()()()
) (1Bu AxPxxPB PBRQPA PA x PxAx BuT T T T T T T T
uBPx xPBu xQx xPBR BPxTT T T T T1
.
Integrăm pe
[,] tt12 și înmulțim cu ½ ecuațiile după care le egalăm și obținem:
)()()(21)()()(21|)()()(21)()()(21
1 1 1 2 2 2
12
12
txtPtx txtPtx txtPtx
tdttxtPtxdtdT T t
tTt
T
)()()(21)()(21)33.2(
1 1 1 2 2 txtPtx tMxtxT T
. ) (212
11 t
tT T T T Tdt PBux PxBuxQPB PBRx
Astfel, funcționala de cost se rescrie:
Ju uRBPx RuRBPxdt xtPtxtTTT
tt
T() ()()() 1
21
21 1
1 1 1
12
.
Deoarec e, R(t) >0 intergrandul este o formă pătratică strict pozitiv definită, rezultă că
J(u) este minim dacă
u RBPxT1 . Altfel, comanda optimală are forma
PxBR uT1 ~ , iar
costul minim devine:
Ju Ju xtPtxtT(~)min () ()()() 1
21 1 1
. (2.34)
Cu cele enunțate mai sus, am demonstrat sinteza Kalman –Letov enunțată prin următorea
teoremă:
22
Teorema 3.2.1 Problema liniar pătratică (2.22), (2.23) are o unică soluție optimală de
tip feedback
~()ut
RtBtPtxtT() ()()()1 (2.35) , unde matricea reală și simetrică P(t) este
soluția ecuației diferențiale Riccati
PAPPA QPBR BPT T 10 , cu condiția finală
MtP)(
.
Pentru cazul în care timpul final
2t tinde spre infinit ecuația matriceală Riccati devine:
01 PB PBRQPA PAT T
.
În continuare căutam soluția EDMR considerând astfel
),(1tt și
(,)tt1 matricele
fundamentale ale sistemelor (2.29) și (2.30)
xAx BR Bp
p Qx ApT
T
1
,
pt Mxt () () 2 2 .
Atunci soluția ecuației diferențiale maticeale Riccati (2.32) este
),( ),( )(21
2 tt tt tP . (2.36)
Demonstrație:
și verifică egalitățile de mai jos:
M tt tttA tttQ ttI tt tttBtRtB tttA tt
TT
),(,),()( ),()( ),(),(,),()()()( ),()( ),(
22 2 2 222 21
2 2
Prin derivarea egalității
),()( ),(2 2 tttP tt obținem
),()( ),( ),()(2 2 2 tttP tt tttP
(înlocuim
),()( ),(2 2 tttP tt )
),()()()()()( )()( )()( )(21tt tPtBtRtBtPtAtPtPtAtQT T
,
iar înmulțind la dreapta cu
),( ),(2 21tt tt obținem (2.32). Dec i unica soluție a ecuației
Riccati este (2.36).
Partiționăm matricea fundamentală
Utt(,)2 corespunzătoare matricei hamiltoniene
H
(2.31) sub forma
23
UttU tt U tt
U tt U tt(,)(,) (,)
(,) (,)211 2 12 2
21 2 22 2
;
prin urmare, soluția problemei Cauchy
),(),()(
),(),(
22
22
tttttH
tttt
,
MI
tttt
),(),(
2 22 2
este
MIttUtttt),(),(),(
2
22 , deci
Mtt U ttU ttMttU ttU tt
),( ),( ),(),( ),( ),(
2 22 2 21 22 12 2 11 2
.
Așadar soluția (2.36) a ecuației Riccati poate fi scrisă și ca
MttUttUtP ),( ),( )(2 22 2 21
1
2 12 2 11 ),( ),( MttU ttU (2.37)
iar sinteza Kalm an–Letov (2.35) capătă forma
)()( )(~txtK tu , unde
)()()( )(1tPtBtRtKT .
Obesrvație:
Dacă matricele A, B, Q și R sunt matrice constante, există matricea hamiltoniană
H
constantă, deci matricea fundamentală este
) (
22),(ttHe ttU .
Din (2.37) reiese că
)(tP (și
)(tK ) sunt matrice variable.
Exemplu 3.2
Fie sistemul
uxx cu
ux, și criteriul
1
02 2 2) ()1( dtux xJ .
Soluție:
Ecuația R iccati este:
1 22 PP P ,
1)1( P .
Atunci matricea hamiltoniană este:
1 11 1H
și rezultă
24
4)2 2( )2 2(
42 242 2
4)2 2( )2 2(
2 2 2 22 2 2 2
t t t tt t t t
tH
e e e ee e e e
e.
Aplicând (2.37) obținem soluția
)1(2 )1(2)1(2 )1(2)2 1( )2 1()(
t tt t
e ee etP
comanda optimală fiind
)()( )( txtP tu .
25
3.3 Problema reglării optimale a mărimii de ie șire
Considerăm problema determinării unei sinteze pentru sistemul complet observabil :
utBxtAx )( )(
,
1 1)( x tx (2.33)
)()(txtCy
(2.34)
care să minimizeze funcționala pătratică de cost referitor la ieșirea
y .
Funcționala pătratică de cost asociată sistemului este:
Ju ytQtytutRtutdt ytMytT
yT
tt
T
y () ()()() ()()() () () 1
21
2
12
2 2
, (2.35)
unde
Qty ()0 , R(t) > 0,
My0 .
Teorema 3.3.1 Comanda optimală (de tip feedback) a sistemului e ste
~() ()()()() ut R tBtPtxtT1
(2.36)
unde P(t) este soluția ecuației diferentiale Riccati
PAPPA CQCPBR BPT T
yT 10
(2.37)
cu condiția finală
Pt CMCT
y ()2
. (2.38)
Demonstrație:
Perechea ( C, A) este observabilă
ieșirea y determină în mod unic starea x .
Înlocuind în (2.35) ieșirea
yCtx() vom obține funcționala (2.23) cu matricele
.CMC MCQCQ
yTyT
Aici, matricele sunt simetrice pozitiv definite pentru că
0yM și
0yQ .
Astfel că , rezultatul este o cosecință a teoremei 3.2.1.
Observații:
26
– Desi problema este formulată în legătură cu mărimea de iesire y
,regulatorul optimal este realizat cu reacție după stare, presupusă accesibilă si măsurabilă.
– Dacă stările nu sunt accesibile si măsura bile, pentru rezolvarea problemei
trebuie introdus un estimator.
27
3.4 Regulatorul liniar pătratic
Considerăm sistemul complet observabil :
() ()() ()(),() xtAtxtBtut xt x 1 1
(2.39)
)()( )( txtCty
(2.40)
si z(t) ieșirea dorită.
Dorim să construim o comandă optimală cu feedback astfel încât aceasta să minimizeze
funcționala pătratică de cost J(u). În acest scop considerăm eroarea
etztyt ()()() și
introducem funcționala dată de
Ju etQtetutRtutdt etMeteT
eT
tt
T
e () ()()() ()()() () () 1
21
2
12
2 2
.
Eroarea se poate exprima si în funcție de stare
)()( )( )( txtCtzte și atunci funcționala
se rescrie corespunzăto r.
Considerăm Hamiltonianul acestei probleme:
Htxup zCx QzCx uRu pAx BuT
eT T(,,,) [( ) ( ) ] ( ) 1
2
.
Condiția de optimalitate (2.20) implică
~u RBpT1 , deci sistemul canonic (2.18),
(2.19) devine
xH
pAx BR BpT
1
( ) pH
xCQzCx Ap CQCx ApCQzT
eT T
eT T
e
.
Deoarece ultima ecuație este n eomogenă, conținând un termen în z căutăm o soluție de
forma
ptPtxtrt () ()()()
(2.41)
care să verifice condiția de transversalitate (2.21) cu stare iniți ală fixată
28
ptxetMetxztCtxt MztCtxtT
e
ttT
ett() () () () ()() () ()()21
21
2
22
Ct MCtxt Ct Mzt Ptxt rtT
eT
e () ()() () () ()()() 2 2 2 2 2 2 2 2
,
unde
xt()2 și
zt()2 sunt vectori arbitrari. De aici rezultă
Pt Ct MCtT
e () () ()2 2 2
,
rt Ct MztT
e () () ()2 2 2 . (2.42)
În consecință, din (2.41) obținem
zQCr PxA CxQC rr PxB PBR PAxxPrxPxPpeT T
eT T ) ( ) (1
.
O soluție de forma (2.41) este dată de soluțiile ecuației diferențiale Riccati și ecuației
diferențiale liniare
PAPPA CQCPBR BPT T
eT 10
(2.43)
zQCrA B PBRreT T T 1
(2.44)
cu condițiile finale (2.42).
Astfel, comanda optimală va fi
~()ut
RtBtPtxtRtBtrtT T() ()()() () ()()1 1
, (2.45)
unde prima componentă este obținută prin reacția de stare și este independentă de z(t), iar a doua
componentă reprezintă o intrare r(t) impusă de valoarea dorită z(t).
O reprezentare a sistemului de mai sus este dat în urmatoarea figură:
)(tz
QCT
)(t MT
e
)(0tr
)(tr
)(1tx
)(tMe
)(tx
)(ty
TB BR1
C
29
4. Ecuația matriceală Riccati cu metode de rezolvare
4.1 Ecuația marticeală algebrică Riccati (EMAR)
Considerăm ecuația marticeală algebrică Riccati (EMAR)
0 XRX XAXAQT
(4.1)
unde maticele
XRQA ,,, și
XQ, sunt matrice simerice (
TQQ ,
TRR ).
În investigarea existenței și unicități i soluțiilor pentru ecuația Riccati, studiul este strâns
legat de subspațiile invariante ale matricei hamiltoniene
nn22 .
Definim matricea hamiltoniană
TA QR AH
(4.2)
și enunțăm următoarea lemă în scopul de a verifica dacă valorile proprii ale lui H sunt
distribuite simetric în raport cu axa imaginară .
Lema 4.2.1 Rădăcinile polinomului
) det()(2HI p sunt distribuite simetric în raport
cu axa imaginară
iR .
Demonstrație: Verificăm
TH JHJ1 , unde J este de forma:
O II OJ
nn . Atunci:
J J1
și
T
nnT
T
nnHO II O
R AAQJA QR A
O II OJHJ
1 1 . Astfel,
Lema 4.2.1 este o conse cință a identității verificate.
Dacă , atunci se admite o forma Jordan a matricei H:
1
21
00
UJJUH
, (4.3)
unde
)(1J și
)(2J și ținem seama că
și
sunt respectiv semiplanul
stang si semiplanul drept din planul complex.
30
Fiind
21
XX
OIUn și
2 1,JJ ,
1J , atunci avem
21
21
XX
XXH
(4.4)
și
)( .
Lema 4.2.2 Dacă și (4.4) este bine definită pentru o matrice
(nu neapărat în formă Jordan), atunci:
1 2 2 1 XX XXT T
. (4.5)
Demonstrație: Având în vedere matricea
1 2 2 1 XX XXZT T calculăm
Z :
T TX X Z2 1
21
XXJ
T TX X2 1
21
XXJH
T TX X2 1
21
XXJHT
T T TX X2 1
21
XXJ
ZT
Prin urmare, matricea Z este o soluție a ecuației
OZ ZT .C
Cu toate acestea, (4.5) rezultă din faptul că
OZ este singura soluție a ecuației de mai
sus.
Teoremă 4.2.1 Dacă și matricea
1X din (4.4) este inversabilă,
atunci:
a) Matricea
1
1 2 XX X este o soluție a ecuației maticeale algebrice
Riccati (4.1).
b)
TXX .
c)
) ( RXA .
Demonstrație:
1X este inversabilă și
1
1 2 XX X , atunci:
1
1 1
XXXI
XIHn n
.
31
Deci,
1
1 1 XX RXA
) (1
1 1 XXX XAQT
.
Prin urmare,
) ( RXAX XAQT
.
Astfel, se verifică a). Următoarea afirmație b) este evidentă prin lema 4.2.2, pe când c)
rezultă din formula
1
1 1 XX RXA și
)( .
Cu toate acestea teorema 4.2.1 stabilește condițiile ce garantează existența unei soluții X a
ecuației matriceale Riccati (4.1) astfel încât
) ( RXA . O teoremă reciprocă a acesteia
este:
Teoremă 4.2.2 Dacă
TXX este o soluție a ecuației EMAR astfel încât
) ( RXA
, atunci matricea Hamilton nu are valori proprii pe axa imaginară (
) și matricea
1X din (4.4) este inversabilă.
Demonstrație: Dacă
X este o soluție a EMAR cu proprietățile aferente atunci:
Tn
A QR A
XIH
RXAXI
XAQRXA
XIn
Tn
.
În plus, după următoarea descompunere Jordan a matricei H obținem:
1
1 PPJ RXA
.
Observăm cu ușurință că
1PJXIPXIHn n
,
cu ajutorul căruia se identifică spațiul n-dimensional generat de coloanele matricei:
XPPPXIn
.
32
Acesta fiind spațiul generat de vectorii proprii ai matricei H corespunzători valorilor
proprii din
.
Deci,( ) și prin urmare,
P X1 este o matrice inversabilă.
Teoremă 4.2.3 Ecuația Riccati are cel mult o soluție
X astfel încât
TXX și
) ( RXA
.
Demonstrație: Considerăm X și Y soluții Hermitice ale ecuației Riccati, astfel încât
) ( RXA
și
) ( RYA . Atunci avem:
00
YRY YAYAQXRX XAXAQ
TT
0 ) () ( YRY XRXAYX YXAT
.
Cu toate astea, deoarece
TYY ecuația rezultată mai sus poate fi rescrisă astfel:
0) )( () () ( RXAYX YX RYAT
.
unde,
YX este soluția unei ecuații de forma
0ZC BZ
.
Aici,
)(B și
)(C . Prin urmare, această ecuație are cel mult o soluție.
Evident
nnOZ , dec i
0Zeste unica soluție, de unde rezultă
YX . Astfel, ecuația Riccati
nu are decât o soluție.
În continuare, următoarea teoremă oferă condițiile în care constrângerile impuse matricei
H din teorema 4.2.1 sunt îndepli nite atunci când:
HBB B și
CCQH .
Teoremă 4.2.4 Considerăm matricele
A ,
B ,
C și presupunem că
nCI Arangn
, pentru
)(
. (4.6)
și
n BI A rangn
, pentru
)(
. (4.7)
Atunci există exact o singură soluție hermitică X pentru EMAR având forma:
33
0 CCX XBB XAXAT T T (4.8)
astfel încât
) ( XBBAT . În plus, acesată so luție X este pozitiv semidefinită în
, iar dacă A, B și C sunt matrice reale
X .
Demonstrație: Fie
T TT
A CCBB AH
și presupunem condițiile din teoremă îndeplinite:
T TT
A CCBB A
yx
yx
,
pentru o anumită alegere
yx, și
. Atunci
yBBxI AT
n ) ( și
CxC yI AT
nT ) (
.
Prin urmare,
2
2, ,) ( yB yyBB yxI AT T
n
și
2
2, ) (, ,) ( Cx CxCx yI Ax yxI AT
nT
n
.
Astfel,
yxI A yxI A yxn n ,) ( ,) ( ,) (
2
22
2Cx yBT
de unde reiese că, dacă
, atunci:
0 0 yBT
și
0Cx .
Asta implică:
0
xCI An
și
0 BI AynT ,
34
când
.. Însă, ținând seama de condiți ile impuse de teoremă acest lucru este valabil
doar dacă
0x și
0y , deci .
În continuare arătăm dacă (4.6) și (4.7) îndeplinesc condițiile din teorema și verificăm
dacă
T TT
A CCBB A
21
XX
21
XX
și
nXXrang
21 ,
unde
,,2 1XX și
)( . Atunci
1X este inversabilă.
Presupunem că,
)(1X Keru . Atunci
uXuXBBT 1 2
și, deoarece
2 1 1 2 XX XXT T conform lemei 4.2.2 avem:
02 1 1 2 2 22
2 uXXuuXXuuXBBXu uXBT T T T T T T
.
Obținem,
02uXBT și
02 1 UXBB UXT
,
care înseamnă faptul că nucleul lui
1X este invariant față de
și
0)(1X Ker , sau
v v
, pentru un punct din și un vector nenul
)(1X Kerv .
În ultimul caz avem
01 1 2 vX vXvXBBT și
vX vXvXAT
2 2 2
echivalent cu
T
nT TBI AXv 02 , pentru
)(
.
Deci,
02vX implică:
21
XX
1 1 0)( 0 0 X X Ker v v
este inversabilă.
Prin urmare, luând în vedere teoremele enunțate anterior, există o soluție hermitică a
ecuației Riccati, astfel încât
) ( XBBAT . Dacă matricele A,B și C sunt reale, atunci
X
este o soluție hermitică a ecuației Riccati, astfel încât
) ( XBBAT .
XX ,
X .
Ceea ce rămâne de verificat este dacă X este o soluție pozitiv semidefinită în raport cu , iar
pentru acest lucru este mai convenabil să rescriem ecuația Riccati
0 CCX XBB XAXAT T T
.
Astfel,
X XBBCC XBBAXXXBBAT T T T T ) ( ) ( care este de forma
Q XAXAT1 1 ,
35
unde
)(1A și
Q 0.
Teoremă 4.2.5 Fie
A ,
B ,
Q ,
R și presupunem că
Q 0,
R
0.
nQI Arangn
, pentru
)(
(4.10)
și
n BI A rangn , pentru
)(
. Atunci există exact o soluție hermitică pentru
ecuația Riccati
01 T TB XBR XAXAQ
, (4.11)
astfel încât
) (1XB BRAT . Prin urmare, această soluție X este pozitiv semidefinită
peste , și dacă A, B și C sunt matice reale, atunci
X .
Demonstrație: Dacă
Q 0, atunci există o matice
C astfel încât
QCCT și
r rangQ rangC
. Astfel, la stabilirea
21
1
BR B putem observa că:
TT
A QB BR A1
T TT
A CCBB A1
.
În plus, ajutându -ne de condiția
CI An
0u
0
uQI An
ne rezultă faptul că:
nCI Arangn
, pentru
)(
,
Din condiția
BI A rangn
21
BRI A rangn
,
ajungem la concluzia că:
n BI A rangn 1
, pentru
)(
.
36
Pentru cel e enunțate mai sus , luăm ca exemplu următorul cod MATLAB cu comenzile
specifice pachetului Control System Toolbox și calculăm stabilizarea soluției ecuației matriceale
algebrice Riccati (EMAR) pentru anumite matrice A,B,R,Q date.
37
4.2 Metode de r ezolvare a ecuației matriceale algebrice
Riccati
Fie ecuația algebrică matriceală Riccati cu forma generală:
01 QPB PBRPA PAT T
Ecuația matriceală Riccati degenerată este neliniară, motiv pentru care în majoritatea
cazurilor nu s e obțin soluții anal itice.
Prezentăm în continuare câteva metode de rezolvare numerică:
4.2.1 Metoda Potter
În scopul enunțării acestei metode prezetăm următorii pași:
o formăm matricea M:
TT
A QB BR AM1
o calculă m valorile proprii ale matricei M astfel încât să aibă proprietate a de
simetrie ; deci
n valori proprii sunt situate în semiplanul stâng și în semiplanul
drept. Dacă notăm cu
n ,…,,2 1 valorile proprii din semiplanul stâng,
avem
0) Re(i ,
n i,…,1 .
o Calculăm vectorii propii
nv vv ,…,,2 1 asociați valorilor proprii
n ,…,,2 1 . Prin
urmare,
ii i v Mv ,
n i,…,1 .
o Partiționăm acești vectorii proprii
iv în doi vectori de dimensiune
n :
ii
insv .
o Formăm matricele S și N :
ns ssS …2 1 ,
nn nn N …2 1 . Astfel,
rezultă s oluția EMAR :
1NSP .
Luăm următorul exemplu spre rezolvare pe ntru o mai bună înțelege a metodei.
Exemplul 4.2 .1
Se dau matricele:
38
0010A,
10B ,
0004Q ,
1R .
Formăm matricea M :
01 00000400000010
M
.
Determinăm v alorile proprii din caracteristicile ecuației:
01 0 00 0 0 41 0 0 00 0 1 0
Ecuația caract eristică este
440 . Se obține valori le proprii:
i11
,
i12 ,
i13 ,
i14 .
Deci, valorile proprii din semiplanul stâng sunt :
i11 și
i12 .
Calculăm vectorii proprii asociați pentru:
i11
0 ) (1 1 vMI
0000
1 1 0 00 1 0 41 0 1 00 0 1 1
4321
xxxx
iiii
.
Un minor caracteristic nenul de ordin maxim este
0 1 01 0 10 0 1
ii
.
39
Avem
1x (necunoscută secundară ) , deci rezolvăm sistem ul :
0 )1( 40 )1(0 )1(
34 22
xix xix i
Obținem
))1(,14,)1(,(2
1 iii v
,
iar pentru
1 găsim :
)2,22,1,1(1 ii i v
.
i12
0 ) (1 1 vMI
0000
1 1 0 00 1 0 41 0 1 00 0 1 1
4321
xxxx
iiii
Un minor caracter istic n enul de ordin maxim este
0 1 01 0 10 0 1
ii
.
Avem
1x (necunoscută secundară) , deci rezolvăm sistemul :
0 )1( 40 )1(0 )1(
34 22
xix xix i
Obținem :
))1(,14,)1(,(2
1 iii v
40
iar pentru
1 găsim:
)2,22,1,1(1 ii i v
.
Formăm matricele S și N:
i iS1 11 1
;
i ii iN2 2)1(2 )1(2
Calculăm
1 11 1
211
ii
iS .
Găsim soluția ecuației Riccati :
2224
1 11 1
1 11 11
ii i iNSP .
Matricea de reglaj este :
2222241011
PBRKT .
În continuare, folosim programul MATLAB cu comenz ile specifice ale pachetului
Control System Toolbox încercăm rezolvarea acestuia după cum ilustrează în print screen -ul de
mai jos :
41
4.2.2 Metoda Schur
În problema determinării valorilor pr oprii ale unei matrice esențială este așa numita
formă Schur care este de structură triunghiulară și, se poat e obține, pentru orice matrice numai cu
transformări unitare (ortogonale în cazul real) ; din această structură particulară, valorile proprii
sunt chiar elementele diagonale.
Teoremă 4.2.2 .1 (FS -Forma Schur ) Pentru orice matrice
A există o matrice
unitară
Q~ astfel încât
AQQSH
este superior triunghiulară și se numește forma Schur a matricei A, iar elementele
diagonale ale matricei S sunt valorile proprii ale matricei A; coloanele matricei Q se
numesc vectorii Schur a matricei A.
Teoremă 4.2.2.2 (FSR -Forma Schur reală ) Pentru orice matrice
A există o
matrice ortogonală
Q astfel încât
p
iii
pppp
TS S A
SS SS S S
AQQS
12 221 12 11
)( )( )(,
…0 0… … … …… 0…
,
unde
iiS sau
iiS după cum sunt asociate valorilor proprii reale , respectiv
complex conjugate , iar valorile proprii ale matricei A sunt evident, elementele diagonale ale
matricei S.
Matricea superior triunghiulară pe blocuri S se numește forma Schur reală a matricei A,
iar coloanele matricei de transformare Q formează vectorii Schur ai matricei A asociați formei
Schur reale S. În multe aplicații însă, vectorii Schur pot înl ocui cu succes vectorii proprii.
Forma Schur reală precum și matricea de transformare se obțin cu comanda schur
specifică pachetului Control System Toolbox a programulul MATLAB :
H schur TU,
,
42
astfel încât
TUTUH și
))(( A sizeeyeUUT ; iar forma Schur complexă, se obține
printr -un apel suplimentar al funcției rsf2csf :
TUcsfrsf TU , 2 ,
.
Alte comenzi specifice sunt:
)(H schurT returnează matricea T , reprezentând forma Schur H.
'',, , lhpHU ordschurTU reordoneză factorizarea Schur dată de comanda
anterioară, astfel încât toate valorile proprii aflate în să fie incluse în factorizare.
Pentru o exemplificare mai bună a metodei Sch ur avem print screen -urile ca exemplele :
Exemplu 4.2.2
Metoda Schur aplicată unei marice cu valori proprii reale
Metoda Schur aplicată unei marice cu valori proprii complexe
43
)3 )(1(3 42421 424) 3)( 7(3 122 7
22
Exemplul 4.2. 3 Aflați decompunerea Schur a matricei A
3 122 7A
.
Soluție:
Găsim valorile proprii:
I A det
Deci, valorile proprii sunt
3,1 .
Pentru fiecare λ vom găsi vectorii lui proprii :
1
4 122 6IAI A
.
Atunci vom avea un sistem omogen de ecuații liniare pe care îl rezolvăm prin metoda
Gauss și obține m:
031
2 1 x x .
Din ecuația sistemului de mai sus găsim variabila
1x :
2 131x x , de unde rezultă că
2231
xxX
.
44
Luăm
32x și obținem
31X , normalizat la
31
101 .
3
6 122 43I AI A
.
Atunci vom avea un sistem omogen de ecuații liniare pe care îl rezolvăm prin metoda
Gauss și obținem:
021
2 1 x x .
Din ecuația sistemului de mai sus găsim variabila
1x :
2 121x x , de unde rezultă că
2221
xxX
.
Luăm
22x și obținem
21X , normalizat la
21
51 .
Alegem bază ortonormată
2 1,vu =
12
51,21
51 .
Astfel,
122 1
51
2,1vu U , iar matricea A în baza ortonormată
2 1,vu este
AUU AUU AT
vu
1 014 35 070 15
51122 1
128 31
51122 1
51
3 122 7
1221
511
,21
45
Deoarece A este
22 putem alege
2 2v u și algoritmul se oprește . Astfel, găsim
descompunere schur
TUTUA , unde
1 014 3T și
122 1
51U .
Prin comenzi le specifice p achetului Control System Toolbox din programului MATLAB
punem următorul print screen cu modul în care putem să rezolvăm exemplul:
46
4.2.3 Metoda matr icei de semn
Considerăm matricea
A cu forma canonică Jordan:
ND TATJ 1
,
unde M este o matrice de vectori proprii,
) ,…,,(2 1 nd dd diagD este o matrice diagonală,
iar N este nilpotentă.
Funcția sign(A) este definită de:
1)(XYX A sign
,
unde
) ,…,,(2 1 ny yy diagY este o matrice diagonală dată de:
0) Re( ,10) Re( ,1
ii
id dacăd dacăy
.
Consideră m
1
22 2112 11)(
n nn
n nn
n nn
O II X
I OZ I
O II X
W WW WK sign
.
Aici, Z satisface ecuația Lyapunov
F FXAZZFXAT2 ) ( ) ( și astfel determină m:
2111
WI WMn
și
2212
WIWN
n .
Deci,
2212
2111
WIWXWI W
nn
(4.12)
Prin urmare, soluția ecuației (4.12), adică
N MX este chiar s oluția aproximativă a
ecuaț iei al gebrice Riccati.
47
Calculul funcției sign(K) îl putem realiza fie prin aproxi marea acesteia utilizând relațiile:
,2,
1
10
k k
k kW WW WA W
iar cât timp matricea
kW este inversabilă șirul de matrice va tinde c ătre valoarea exactă a funcției
sign(K) . Deci, calculând
NM Xk1 , șirul
kX va tinde către soluția exactă a EMAR.
Exemplificăm mai jos printr -un print screen un program ul al matricei de semn construit
cu ajutor MATLAB prin comenzile specifice p achetului Control System Toolbox ale acestiua:
48
u(t) l m
θ(t)
x(t) 5.Aplicații
5.1 Aplicații teoretice ale controlului optimal
Exercițiul 5.1.1 (Stabilizarea pendulului invers)
Considerăm pendul invers descris prin figura de mai jos
unde:
l = lungimea pendulului;
m = masa pendulului;
M = masa căruciorului;
θ = unghiul de deviație al pendulului față de poziția verticală;
u = comanda motorului asupra căruciorului.
și sistemul de ecuații:
lgu
nn
22
.
Soluți e:
Definim variablele de stare
x1 ,
2x și obținem ecuația
uxx
xxB AX X
n
10
010
21
2
21
M
49
asociată funcționalei de cost
02
221
21dtu
cJ
, unde
0001Q și
21
cR .
Luăm ecuația marticeală algebrică Riccati (EMAR)
01 Q B PBRPA PAT T
cu maticea simetrică
3 22 1
p pp pP și găsim:
0000
0001
2 12
32
2
22
312
2
p pp p
p pp pn n
nn
.
Rezolvând operațiile de mai sus găsim ecuțiile:
2 32
32
2322 2
3 12 4 2
2 22
22 2
2
210 20101 2
pcp pc pppc p pc
cp pc p
nn n n
Cu toate acestea,
3p este un termen di agonală, care trebuie să fie real și pozitiv.
Prin urmare,
2p trebuie să fie pozitiv. Așadar,
p ppc pppc p p
nn
2
3 322
1322 2
3 1 0
Găsim matricea K și comanda optimală u de forma:
pc pc PBRKT
32
22 1
Trecând în spațiul stărilor , avem variabilele de stare:
P P PP PP PP
x x xx xx xx
3 431 21
uMgxMmxx xuMlgxMlm Mxx x
P PP PP PP P
11
1 44 31 22 1
Utilizăm forma matricială vectorială standard putem scrie:
50
u
MgMl
xxxx
gMmgMlm M
xxxx
PPPP
PPPP
1010
000100 0000001 0
4321
4321
Considerăm
și
Px ieșirile sistemului:
31
21
P Pxx
x yy
yP
PP
P
4321
21
01000001
PPPP
PP
xxxx
yy
ce reprezintă modelul de stare.
Funcția de transfer a sistemului în buclă deschisă este obținută utilizând programul
MATLAB prin utilizarea comenzilor specifice p achetului Control System Toolbox , după cum
reiese din print screen -ul de mai jos :
Din răspunsul sistemului în buclă deschisă se poate observa că acesta nu este stabil.
Exercițiul 5.1. 2
Considerăm sistemul dubu integrator
51
)()( )(2 )()( )(
2 1 22 1
tutxtx txtxtx
cu condițiile ințiale ,
2)0(1 x ,
3)0(2 x și funcționala de cost asociată
dttu tx txtx tx x x x x J 5
02 2
2 2 12
12
2 2 12
1 )( 25,0)(5)()(6)(221)5(2)5()5( )5(21
pentru care dorim să aflăm comanda optimală și un grafic reprezentativ al ei.
Soluție:
Avem matricele asociate sistemului:
1210A
,
10B ,
25.05.01)(ftF ,
5332Q ,
41rR ,
00t ,
5ft .
Luăm
22 1212 11
p pp pP o matrice simetrică a ecuației EMAR:
01 QPB PBRPA PAT T
Astfel, controlul optimal este dat de:
*
2 22*
1 12 *
2*
1
22 1212 11 *4 104 )( xp xp
xx
p pp ptu
, unde,
533210410
1210
1210
22 1212 11
22 1212 11
22 1212 11
22 1212 11
22 1212 11
p pp p
p pp p
p pp p
p pp p
p pp p
satisface condiția finală
25.05.01)(ftF astfel încât , obținem sistemul de ecuații:
, 5 4 2 2, 3 4 2, 2 4 4
222
22 22 1212 22 12 22 12 1111 122
12
p p p pp pp p p pp p p
2)5(5.0)2(1)5(
221211
ppp
În programul MATLAB folosim comanda lqr și initial specific e pachetului Control
System Toolbox pentru a rezolva acest exemplu . Prin urmare, pentru a pute a efectua acest lucru
este necesar să precizăm anumite caracteristici:
52
– comanda optimală va fi
)( )(*tKxtu , unde
PBR KT1 ;
– e = eig(A-B*K) matricea valorilor p ropii și matricea P soluția ecuației matriceale Riccati;
– definim matricele C,D și intervalul de timp
5,0t pentru realizarea graficului
traiectoriei optimale;
– definim G un vector fictiv necesar pentru comanda matlab "initial";
Ilustrăm p rint screen -ul realizat cu ajutor MATLAB pentru o mai bună observare :
Exercițiul 5.1.3
Considerăm sistemul
ux x
și funcționala de cost
dtu x JT
02 2) (21
.
pentru care se cere să găsim soluția ecuației diferențiale Ri ccati corespu nzătoare.
Soluție:
Avem matricele:
53
1A,
1B ,
1Q ,
1R
Necunoscuta P a ecuației Riccati se va nota cu p care este o funcție scalară. Ecuația
diferențială matriceală Riccati devine o e cuație diferențială ordinară de forma:
12 ppp p
adică,
01 22 pp p cu condiția la limită
0)(Tp .
Întrucît termenul
2)(21
STx nu există în expresi a indicelui de performanță, avem S = 0 .
Prin urmare ,
yu xy u x uyxtH ) (21),,,(2 2 rezultând din
0uH și
y yxtcu ),,(
.
Hamiltonianul problemei va fi:
xy y x yxtH 2 2
021
21),,(
Sistemul canonic rezultă ca fiind
yx yyx x
care coincide cu sistemul canonic matriceal, având condițiil e
1)(Tx ,
0)(Ty
corespunzătoare respectiv cu
I TX)( ,
S TY)( .
Matricea, M a sistemului canonic este:
111 1M
care are valorile proprii
21 ,
22
De aici, rezultă:
54
) ( ) () ( ) (22)2 1( )2 1(
2222 22)2 1( )2 1(
22 2112 11) (2 ) (2 ) (2 ) (2) (2 ) (2 ) (2 ) (2
) (
Tt TtTt Tte e e ee e e e
eTt Tt Tt TtTt Tt Tt Tt
TtM
Ținem cont că S=0, atunci avem soluția ecuației Riccati:
) (2 ) (2) (2 ) (2) (2 ) (2) (2 ) (21
11 21
)12( )12()2 1( )2 1(22
22) () ( )0,,(
Tt TtTt TtTt TtTt Tt
e ee ee ee eTt Tt Ttp
Comanda optimală este dată de:
xTtp u )0,,( și observăm că
2 1
121)0,,( lim
Ttp
T
este chiar soluția ecuațiri Riccati
01 22p p .
Exercițiul 5.1. 4
Considerăm sistemul
u x x
10
0010
și funcționala de cost
dtu x J
02 2) 4(21
.
Soluție:
Matricele asociate sistemului:
0010A
,
10bB ,
0004Q ,
0002C ,
1R
55
Obervăm că perechea (A,B)este controlabilă,
20110
rang ABB rang
,
iar perechea (A,C) este complet observabilă deci, matricea
02000002T T TCAC
are rangul maximal 2.
Prin urmare, ecuația algebrică Riccati
) (T are deci o soluție unică pozitiv definită:
01 QPB PBR PAPAT T
cu
22 1212 11
p pp pP . Astfel, obținem:
000041010
0010
0100
22 1212 11
22 1212 11
22 1212 11
22 1212 11
p pp p
p pp p
p pp p
p pp p
Prin rezolvare operațiior de mai sus găsim sistemul de ecuații algebrice:
0 2004
2
22 1222 12 112
12
p ppp pp
a cărei sol uție pozitiv definită este
411p
,
212p ,
222p
adică,
2224P .
Atunci comanda optimală va fi:
56
) (2222410 ),(2 1
21 1xxxxPxb PxBR xtkuT T
.
În programul MATLAB folosim comenzile specifice pachetului Control System Toolbox
putem rezolva ușor acest exemplu. Mai jos arătăm print screen -ul relalizat în acest program:
Exercițiul 5.1. 5
Considerăm sistemul descris de:
uxx
xx
10
2 11 0
21
21
,
x y 01
și indicele de performanță
dtux x JT
02
1002
.
Dorim să determinăm matricea Riccati P , matricea feedback K și valorile proprii ale
buclei închise .
Soluție:
57
Avem matricele asociate:
2 11 0A ,
10B ,
1002Q ,
1R .
Calculăm ecuația Riccati
01 PB PBRQPA PAT T .
22 21 2212 11 12
22 1212 11
22
2 11 0
p p pp p p
p pp pPA
22 12 21 1122 21
22 1212 11
2 2 2 11 0
p p p pp p
p pp pPAT
2
22 21 2222 12 21 12
22 21
2212
22 1212 11
22 1212 11 110110
p pppp ppp ppp
p pp p
p pp pPB PBRT
Combinăm ecuațiile de mai sus și avem:
01002
2 2 22
2
22 21 2222 12 21 12
22 12 21 1122 21
22 21 2212 11 12
p pppp pp
p p p pp p
p p pp p p
.
Deoarece, P este simetric
12 21p p obținem sistemul de ecuații:
0 1 2 20 20 20 2
2
22 22 12 22 1222 12 12 11 2222 12 22 12 112
12 12 12
p p p p ppp p p ppp p p pp p p
Rezolvăm sistemul și obținem
732.021 12p p ,
542.022p ,
403.211p .
Astfel, ecuația Riccati este
542.0 732.0732.0 403.2P .
Găsim apoi și feedback -ul matricei:
542.0 732.0732.0 403.21011PBRKT
Prin urmare,
542.0 732.0K .
Valorile prop rii ale buclei închise sunt:
58
0 542.0 732.010
2 11 0
00
ss
0542.0 732.00 0
2 11
ss
0542.2 732.11
ss
0 732.1 542.22 s s ,
i ss 340.0 271.1 ,21 .
În programul MATLAB folosim comenzile specifice pachetului Control System Toolbox
și rezolvăm acest exemplu ilustrăndu -l printr -un print screen :
59
2K
1K5.2 Sistemul de suspensie la autobuz
5.2.1 Noțiuni introductive pentru reprezentarea în spațiul
stărilor
Controlul optimal este cea mai răspândită tehnică de contro l pentru sistemul de suspensie
deci, vehiculele sunt compuse din mai multe sisteme , iar unul dintre acestea fiind chiar sistemul
de suspensie. Proiectarea unui astfel de sistem de suspensii pentru autovehicule este o problemă
de control interesantă ș i provocatoare , deoarece p rincipalele funcții ale sis temului de suspensii
pentru vehicule sunt de a asigura suportul vehiculului, stabilitatea și controlul direcțional în
timpul manevrelor . Atunci când sistemul de suspensie este proiectat, pentru un anumit model , la
una din cele patru roți, este utilizat pentru simplificarea problemei un sistem 1D de amortizare cu
arc.
Prin urmare , prezentăm un model ce este destinat unui sistem de suspensie activă , în care
este inclus un servomotor care poate genera forța de comandă u pentru a controla mișcarea
autobuzului. O diagramă a acestui sistem este prezentată în imaginea de mai jos:
Aici, parametrii sistemului sunt:
1M
= masa autobuzului
1M
2M
u
1b
u
2b
W
1X
2X
60
2M = masa suspensiei
1K
= constanta de elasticitate a sistemului de suspensie
2K
= constanta de elasticitate a roții și a cauciucului
1b
= constanta de amortizare a sistemului de su spensie
2b
= constanta de amortizare a roții și a cauciucului
W = perturbația drumului
Considerăm ecuațiile de mișcare :
u xxK xxb xM ) ( ) (2 11 2 11 11
(1)
u x WK x Wb xxK xxb xM ) ( ) ( ) ( ) (2 2 2 2 2 1 1 2 1 1 22 (2)
În continuare efectuăm calcule asupara celor două ecuații de mai sus pentru a găsi spațiul
stărilor și vom începe prin a împărțiți (1) și (2) cu M1 , și respectiv M2 :
22
22
2
22
2 1
21
2 1
21
212 1
11
2 1
11
1
) ( ) ( ) ( ) () ( ) (
MuxWMKxWMbxxMKxxMbxMuxxMKxxMbx
Înlocuim în ecuațiile de mai sus
2 1xx =
1y :
22
22
2
22
1
21
1
21
211
11
1
11
1
) ( ) (MuxWMKxWMbyMKyMbxMuyMKyMbx
Găsim ecuația
uM MxWMKxWMbyMK
MKyMb
Mbxx y
2 12
22
2
22
1
21
11
1
21
11
2 1 11 1) ( ) (
care prin integare obținem:
dtuM MxWMKyMK
MKxWMbyMb
Mby
2 12
22
1
21
11
2
22
1
21
11
11 1) ( ) (
Căutăm mai departe ecuația de stare pentru
2y :
61
uM MxWMKyMK
MKy
2 12
22
1
21
11
21 1) ( .
Ținem cont de faptul că,
1 1yx =
2x și obținem:
uM MyxWMKyMK
MKy
2 11 1
22
1
21
11
21 1) (
.
Înlocuim
1 1 2 yx x în
1y pentru a o bține ecuația de stare pentru
1y :
2 1 1
22
1
21
11
2 ) ( y yxWMbyMb
Mby
.
Mai departe , vom înlocui
1y în ecuația
1x :
WMMbb
MuyMbyMK
Mb
Mb
Mb
MbxMMbbx
2 121
12
11
1
11
22
21
11
11
1
2 121
1
Prin urmare, ne rezultă reprezentarea în spațiul stă rilor a sistemului de suspensie la
autobuz
BW Axx
DW Cxy
care va fi pusă în următoarea formă în spațiul stărilo r:
2111
yyxx
=
Wu
MK
M MMbMMbb
M
yyxx
MK
MK
MK
MKMb
Mb
Mb
MbMb
MK
Mb
Mb
Mb
Mb
MMbb
22
2 1222 121
1
1111
22
21
11
2222
21
11
2211
11
22
21
11
11
2 121
1 1010 0
0 01 000 0 1 0
Wu
yyxx
y 00 0100
1111
.
62
5.2.2 Compararea metodelor de rezolvare a ecuației
Riccati și aplicarea lor
În continuare, pentru a putea compara met odele de rezolvare a ecuației Riccati este
necesar să definim ur matoarele condiții referitoare la suspensia autobuzului:
1M =2.500 kg ;
2M
=320 kg;
1K =80.000 N/m;
2K =500.000 N/m;
1b= 350
smN ;
2b = 15 .020
smN ;
Astfel, obținem matricele sistemului de suspensie :
0 5. 1844 0 5. 15621 17125.48 0 9375.4614.0 256025.25 0 57125.60 0 1 0
A
,
5. 1562 003525.09375.46 057125.6 0004.00 0
B
0100C
,
00D .
De asemenea, prin calcule găsim și matricele:
0000010000000000
Q ,
1001103R .
Ne propunem să găsim o comandă optimală fosindu -ne de metodele ecuației Riccati și să
cercetăm dacă acest sistem de suspensie al autobuzului este unul stabil și se poate încadra în
tiparul practic întâlnit în realita te.
Pornim ca și prim pas prin verificarea matricelor, astfel încât acestea să respecte
condițiile pentru existenț a unei soluții a ecuației matriceale algebrice Riccati, precum ș i
stabilitatea matricei A folosin du-ne de MATLAB prin comenzile specific acest uia din pachetului
Control System Toolbox realizând pentru observarea mai detaliată un print screen a datelor
introduse în veder ea obținerii rezultatelor .
63
În urma rulării acestor date vom observa că toate valorile proprii ale matricei A au partea
reală negativă , iar controlabilitatea și observabilitatea este îndeplinită. Așadar, matricea a
rezultat a fi stabilizabilă și sistemul este încadrat în tiparul practic din realitate . Cu alte cuvinte
putem spune că există o unică soluție pozitiv semidefinită a e cuației matriceale algebrice Riccati.
Mai departe, tot cu ajutor MATLAB folosind comanda lqr specifică Control System
Toolbox ilustrăm căteva grafice specifice modul în care f uncționează sistemul (fără n ici un
control al feedback -ului) la introducerea u nui pas de perturbare , cum ar fi graficul răspunsului
sistemului în buclă deschisă și răspunsul sistemului în buclă închisă. Din graficul buclă deschisă
se poate observa că sistemul este sub amortizat , deci c orpul autobuzului va oscila pentru o
perioadă ina cceptabil de lung. Spre exemplu, o amenii care stau în autobuz nu se vor simți
confortabil cu o astfel de oscilație. Pe când la graficul buclă închisă oamenii care stau în autobuz
vor simți o foarte mică oscilație într-un timp mai scurt.
64
Aplicăm metodele:
Metoda Potter
Realizăm un print screen din programul MA TLAB pentru a ilusta acestă metodă și
observăm faptul că în comparație cu cu metodele iterative aceasta prezintă avant ajul unui calcul
direct având soluție unică .
Metoda Schur
Realizăm un print screen în care vom afla matricea hamiltoniană H, precum și matricele
U și T așa cum sunt descrise de metodă . Apoi , continuăm cu găsirea unei grafic pentru a arăta
65
prezența frecve nței la suspensia unui autobuz și găsirea unei aproximări a ecuației matriceale
Riccati pe care o vom verifica prin calculul unei norme a diferenței celor două matrice obținute
anterior.
Din rularea programului observăm că aproximarea este convenabilă, ordinul normei are o
valore negativă de
1210 1013.2 .
Metoda matricei de semn
66
În print screen -ul de mai sus f olosim func ția sign( K) și calculăm matricea de semn pe
care o vom nota cu Ws și ne reiese din rularea programului prin încercarea calcul ării aproximării
soluției Riccati cu această metodă că există o diferență în normă a soluției, atunci când se
folosește comanda lqr și sign. Însă cu toate acestea, metoda nu este adecvată pentru acest s istem
de suspensie. De asemenea, tot pentru funcția sign(K) intervine în acest caz și calculul direct
folosindu -ne de funcția Jordan dată în programul MATLAB , dar această metodă are un
dezavantaj, deoarece în cazul când matricea K are valori întregi obținem un rezultat corect pentru
că valorile se pot rotunji , dar în schimb când avem valori reale cu mai multe zecimale rezultatul
nu va fi corect. Cu alte cuvinte aceasta metodă nu poate fi folosită în cazuri reale.
Comparația dintre aceste trei metode ar fi aceea că metoda Potter se bazează pe calculul
valorilor proprii, iar aproximarea nu e cea mai bună; metoda matricei de semn poate avea un
rezultat mult mai apropiat de soluția exactă pentru sistemul și pare a fi cea mai bună aproximare
a soluției ecuației Riccati, însă metoda Schur este cea mai potrivită ca fi ind aproximarea a
soluției ecuației Riccati , deoarece prezența ordinului negative al normei face acest lucru posibil
pentru acest sistem de suspensie la autobuz , dar nu numai, pentru că de orice tip ar fi sistemul
rezultatul va fi unul bun .
67
6. CONCLUZII
68
7. ANEXE
clear all; clc;
A=[0 1 0 0; -6.57125 0 -25.256025 -0.14;46.9375 0 -48.17125 1;
1562.5 0 -1844.5 0]
B=[0 0;0.0004 6.57125; 0 -46.9375; 0.003525 -1562.5]
C=[0 0 1 0]
D=[0 0]
Q=[0 0 0 0; 0 0 0 0; 0 0 1 0;0 0 0 0]
R=10^(-3)*[1 0;0 1]
%matricea Kamilton
K = [A -B*R^(-1)*B'; -Q -A']
%calculul comenzii optimale si a rasp unsului sistemului in bucla
inchisa folosind metoda Sch ur
[U,T]=schur(K)
[U,T]=ordschur(U,T, 'lhp')
%aproximarea so lutiei Riccati prin metoda ScKur
[Ks,S,e]=lqr(A,B,Q,R)
n=length(A)
aprox_sol = U(n+1:end,1:n)*U(1:n,1:n)^( -1) %solutia de aprox a
ecuatiei Riccati
K1s= inv(R)*B'*aprox_sol
K2s=inv(R)*B'*inv((A' -aprox_sol*B*inv(R)*B'))*aprox_sol*e
%calculul diferentei in norma intre matricea S calculata cu comenzi
implicite matlab si cea calc ulata folosind descompunerea Sch ur
norma_schur=norm(S -aprox_sol,1)
M1=320;M2=2500;b1=350;b2=15020;K1=80000;K2=500000
% definire numitorului pentru introducerea fortei actionate
numitor=[(M1+M2) b2 K2];
% definirea numaratorului pentru introducerea fortei actionate
numarator=[(M1*M2) (M1*(b1+b2))+(M2*b1) (M1*(K1+K2))+(M2*K1)+(b1*b2)
(b1*K2)+(b2*K1) K1*K2];
figure(2)
w = logspace( -1,2) %genereaza un vector de rând de logaritmic
bode(numi tor,numarator,w)
numitor=100000*numitor
grid on
69
%calcularea solutiei ecuatiei Riccati folosind functia creata, Sign
%foloseste descompunerea Jordan proprie matlab
n=length(A)
I=eye(n)
Ws=sign(K)
W11 = Ws(1:n,1:n);
W12 = Ws(1:n,(n+1):end);
W21 = Ws((n +1):end,1:n);
W22 = Ws((n+1):end,(n+1):end);
M1 = [W12; W22+I];
N1 = [W11+I; W21];
Ssemn_1= -M1\N1;
norma_semn1=norm(S -Ssemn_1)
%aproximarea functiei de semn
W_1 = K
Ws = W_1+ones(length(W_1))*10^( -13);
W_1 = Ws;
Ws=W_1-(W_1-inv(W_1))/2;
W11_2 = Ws(1:n,1:n) ;
W12_2 = Ws(1:n,(n+1):end);
W21_2 = Ws((n+1):end,1:n);
W22_2 = Ws((n+1):end,(n+1):end);
M_2 = [W12_2; W22_2+I];
N_2 = [W11_2+I; W21_2];
Ssemn=-M_2\N_2
disp('norma aproximarii folosind matricea de semn' )
norma_semn2=norm(S -Ssemn)
Ksemn=inv(R)*B'*Ssemn;
Kffsemn=inv(R)*B'*inv((A' -Ssemn*B*inv(R)*B'))*Ssemn*e;
figure(3)
syssemn=ss(A -B*Ksemn,B,C,D)
step(A-B*Ksemn,B,C,D,2, 'r')
axis([0 0.35 -10^2 4])
grid on
70
8. BIBLIOGRAFIE
Copyright Notice
© Licențiada.org respectă drepturile de proprietate intelectuală și așteaptă ca toți utilizatorii să facă același lucru. Dacă consideri că un conținut de pe site încalcă drepturile tale de autor, te rugăm să trimiți o notificare DMCA.
Acest articol: Matematică și Informatică Aplicate în Inginerie [606011] (ID: 606011)
Dacă considerați că acest conținut vă încalcă drepturile de autor, vă rugăm să depuneți o cerere pe pagina noastră Copyright Takedown.