Control Optimal
92
UNIVERSITATEA POLITEHNICA DIN BUCUREŞTI Facultatea de Ştiinţe Aplicate Specializarea: Matematică şi informatică aplicate în inginerie CONTROL OPTIMAL Profesor coordonator: Student: Prepeliţă Valeriu Petcu Andreea Cristina 1
-
Upload
andreea-cristina-petcu -
Category
Documents
-
view
159 -
download
12
Transcript of Control Optimal
UNIVERSITATEA POLITEHNICA DIN BUCUREŞTI
Facultatea de Ştiinţe Aplicate
Specializarea: Matematică şi informatică aplicate în inginerie
CONTROL OPTIMAL
Profesor coordonator: Student:
Prepeliţă Valeriu Petcu Andreea Cristina
Bucureşti 2009
1
CUPRINS
1. Introducere
2. Optimizare Statică
2.1 Optimizarea statică în circuit deschis
2.2 Optimizarea statică în circuit închis
3. Optimizare Dinamică
3.1 Optimizare prin metode variaţionale
3.2 Programarea dinamică
3.3 Principiul maximului al lui Pontriaghin
4. Aplicaţii
4.1 Aplicaţie optimizare prin metode variaţionale utilizând teorema de optimalitate
4.2 Aplicaţie optimizare prin metode variaţionale utilizând ecuaţia matricială Riccati
4.3 Aplicaţie optimizare prin metoda programării dinamice
5. Concluzii
6. Anexe
7. Bibliografie
2
1. Introducere
O anumită problemă inginerească nu are soluţie unică; în cele mai multe cazuri, ea posedă o infinitate de “soluţii posibile”. Dacă aceste soluţii variate sunt evaluate pe baza unei “funcţii de cost” specifice (numită şi funcţie de pierderi sau criteriu de calitate), fiecare soluţie va fi caracterizată printr-o anumită valoare a acestei funcţii. Soluţia care oferă o valoare optimă se numeste soluţie optimală. Căutarea unei astfel de soluţii optimale pentru o problemă dată se justifică din ce in ce mai mult in contextul actual al dezvoltării tehnologiei, al restricţiilor energetice si de resurse naturale, al pretenţiilor tot mai ridicate privind calitatea produselor.[1]
În această lucrare vom studia în profunzime optimizarea dinamică a proceselor. Vom prezenta trei metode de optimizare dinamică: optimizarea prin metode variaţionale (metoda multiplicatorilor Lagrange şi metoda ecuaţiei matriciale Riccati) ([2], [3]), metoda programării dinamice a lui Bellman [1] şi metoda principiului maximului al lui Pontriaghin. [1]
Studiul sistemelor optimale se face pornind de la următoarele elemente fundamentale:
1. Dinamica procesului este cunoscută, fiind exprimată sub forma unei ecuaţii de stare
(1.1)
presupusă în cele mai multe cazuri a fi liniară.
2. Se limitează atenţia la cazul determinist, adică se admite că sistemul nu este supus unor perturbaţii stohastice (iar măsurarea vectorului de stare nu este alterată de zgomote).
3. Vectorul de stare x şi/sau comanda u sunt supuse unor restricţii diverse – de exemplu, componentele lui u trebuie să se încadreze între anumite limite (restricţii de tip inegalitate), sau se pot impune anumite valori iniţiale şi finale pentru x (restricţii de tip egalitate), sau se impun chiar restricţii privind durate tf a conducerii procesului respectiv.
4. Se alege un indice de performanţă J, dependent în general de x şi de u şi crescător cumulativ cu timpul, adică de forma:
(1.2)
în raport cu care se realizează optimizarea.
5. Nu se face nici o presupunere cu privire la structura dispozitivului de conducere optimală (spre diferenta de cazul proiectării prin încercări).
Cu acestea, problema comenzii optimale se formulează concis astfel:
Să se găsească comanda optimă uopt (numită şi “strategia de conducere optimală”) care optimizează indicele de performanţă J, în condiţiile unor restricţii impuse.
În cadrul problemelor de optimizare se disting două categorii:
3
1. Probleme de optimizare statică.
2. Probleme de optimizare dinamică.
Optimizarea statică se întalneşte atunci când se cere conducerea optimală a unui proces a cărui stare normală este regimul staţionar. Fie astfel un proces descris de ecuaţia (1.1), pe care dorim să-l menţinem în starea staţionară constantă
(1.3)
Într-un astfel de regim, ecuaţia diferenţială vectorială (1.1) se reduce la o ecuaţie algebrică vectorială (staţionară):
(1.4)
ecuaţie a cărei rezolvare ne dă valoarea comenzii u0 constante, care poate menţine sistemul în regimul staţionar dorit. În acest tip de probleme suntem interesaţi mai mult de viteza de variaţie a funcţiei de cost (indicelui de performanţă) J decât de valorea însăşi a lui J. Din (1.2) se obţtine direct:
(1.5)
Dacă x şi u se reduc la nişte constante, rezultă că viteza de variaţie a indicelui de performanţă capată valoarea constantă
(1.6)
ceea ce indică posibilitatea optimizării procesului prin optimizarea vitezei de variaţtie a lui J.
Exemplu:
Să consideram un sistem electric compus din n generatoare G1, G2… Gn care furnizează
unor consumatori respectiv puterile (în kW). Întrucât energia electrică nu poate fi
stocată, va trebui ca puterea totală, în fiecare moment, să fie egală cu cererea P a tuturor consumatorilor, cerere care nu poate fi influenţată de către producătorii de energie electrică. Aşadar, valorile puterilor din sistem vor trebui să satisfacă restricţia de tip egalitate
(1.7)
Această condiţie este satisfăcută automat prin reglarea frecvenţei (la valoarea de 50 Hz). O crestere a frecvenţei denotă un surplus de putere produsă în sistem, surplusul de energie din sistem fiind utilizat la creşterea energiei cinetice a acestuia (creşterea vitezei generatoarelor).
4
Problema care se pune este aceea a repartiţiei producţiei pe fiecare generator. Se va remarca faptul că sarcina P cerută de consumatori reprezintă o marime lent variabilă în raport cu
vitezele cu care se pot varia ieşirile ale generatoarelor, ceea ce permite considerarea relaţiei
(1.7) drept o ecuaţie staţionară. În rezolvarea acestei probleme, un indice de calitate judicious de ales apare ca fiind costul total al energiei produse în sistem, în condiţiile în care preţul de producţie al energiei electrice variază de la centrală la centrală (după tipul energiei primare
utilizate, puterea generatoarelor, costul exploatării etc.). Fie rata (viteza) de variaţie a
costului (exprimată de exemplu în lei/h) pentru a produce puterea (în kW) la centrala i. Se
obţine deci rata de cost totală – staţionară – in sistem:
(1.8)
(în lei/h). Rata de cost poate fi considerată ca jucând rolul unei “comenzi” pentru generatorul
Gi şi la care el răspunde prin marimea de ieşire . Aşadar, problema se reduce la minimizarea
ratei de cost totale (1.8) în condiţiile respectării restricţiei de tip egalitate (1.7).
Pentru ilustrarea problemei optimizării dinamice, ne vom referi la un alt exemplu.
Exemplu:
Fie problema poziţionării pe o orbită circumterestră a unui vehicul spaţial, în scopul cuplării cu un alt vehicul spaţial. Se va considera cazul simplificat al poziţionării doar în planul orizontal local. La o rotire a vehiculului spaţial cu unghiul φ, cauzată de plidă, de trecerea pe o altă orbită, nişte reactoare laterale vor acţiona cu cuplul C, care să readucă vehiculul în poziţia corectă. În absenţa frecărilor şi pentru perturbaţii neglijabile, considerând totodată consumul de combustibil al reactoarelor de poziţionare mult mai mic decât masa vehiculului, ecuaţia mişcării vehiculului va fi dată de:
(1.9)
unde J este momentul de inerţie al vehiculului în jurul axei de rotaţie. Relaţia (1.9) se mai poate scrie:
(1.10)
şi, adoptând ca variabile de stare si , avem ecuaţtiile de stare:
(1.11)
unde
5
(1.12)
reprezintţ mărimea de comandă. Dar cuplul C este proporţional cu viteza v de expansiune a gazelor (din reactoarele de poziţionare) faţă de vehicul şi cu viteza de ardere a combustibilului, după relaţia:
(1.13)
unde m reprezintă masa de combustibil. Aşadar, comanda poate fi influenţată prin viteza de ardere . Se va propune drept criteriu de performanţă o valoare cât mai mică de combustibil
consumat pentru poziţionare, adică
(1.14)
sau, ţinând cont de (1.12) şi (1.13)
(1.15)
Dar extremul valorii lui J[u] nu depinde de constantele J şi v, astfel că în loc de criteriul (1.15) putem lua
(1.16)
Comparând cu (1.2), rezultă că funcţia L (numită şi funcţie de pierdere, căci ea reprezintă o măsură a îndepărtării instantanee faţă de valoarea ideală a indicelui de performanţă) va fi îin acest caz:
(1.17)
Problema de optimizare va consta în a găsi variaţia optimă a lui u(t) – numită şi strategie optimă – care să minimizeze integrala (1.16), satisfăcând totodată diversele restricţii impuse de cazul fizic concret. Fireşte, astfel de restricţii – de tip egalitate – se vor impune variabilelor de stare:
(1.18)
Se observă, din acest punct de vedere, că problema de optimizare trebuie să rezolve transferarea sistemului dintr-o stare iniţială dată într-o stare finală dată (astfel de probleme sunt cunoscute sub numele de probleme “cu condiţii pe frontiere”). În plus, s-ar putea ca să nu avem disponibil decât
6
un reactor de poziţtionare; de asemenea, este clar că valoarea lui u(t) este limitată (fizic) la o valoare maximă umax. Aşadar şi comanda poate fi supusş unor restricţii (de tip inegalitate) având
forma: .
În fine, trebuie remarcat că nu a fost pusă în discuţie până acum nici o limitare asupra duratei totale tf a procesului de optimizare. Fireşte, dacă comanda u ar putea fi oricât de mare, durata tf ar putea fi facută arbitrar de mică. Dar, după cum s-a menţionat, limitările fizice asupra
valorii comenzii ridică chestiunea încadrării într-un interval de timp specificat .
În exemplul analizat starea sistemului se modifică neîncetat, iar realizarea strategiei optime uopt(t) trebuie realizată pe măsura şi în concordanţă cu aceste modificări. Se vede astfel că avem de-a face cu o optimizare dinamică.
2. OPTIMIZARE STATICĂ
Optimizarea statică a proceselor poate fi realizată în două moduri:
- ca optimizare în circuit deschis, utilizabilă atunci când modelul matematic al unui proces este bine cunoscut (inclusiv valorile parametrilor), iar efectele perturbaţiilor sunt neglijabile.
- ca optimizare îin circuit închis, pentru procesele în care parametrii sunt insuficient cunoscuţi şi/sau variază cu timpul, iar perturbaţiile – de tip aleator – au o puternică influenţă asupra stării sistemului; această metodă se bazează pe experimentare.
2.1 OPTIMIZARE STATICĂ ÎN CIRCUIT DESCHIS
Întrucât, după cum am arătat, în cazul optimizării statice, mărimile u0 si x0 sunt legate printr-o ecuaţie algebrică staţionară (1.4), este posibil ca – cel puţin în principiu – să se rezolve această ecuaţie în raport cu u0 şi, după înlocuirea acestei valori în (1.6), să se exprime viteza de variaţie a indicelui de performanţă J ca o funcţie numai de starea staţionară x0:
(2.1.19)
funcţia H luând diverse valori în funcţie de vectorul particular x0 ales. Obţinerea valorii lui x0
care să optimizeze funcţia H depinde de: natura particulară a funcţiei H, dimensiunea lui x0, restricţiile impuse.
Dacă se reperzintă funcţia într-un spaţiu (n+1) – dimensional, se obţin nişte
contururi de valori Hi=const.
7
Figura 1 Contururi H=const la sisteme bidimensionale
Această situaţie se poate vizualiza uşor în cazul n=2, adică pentru . Se observă
că unele suprafeţte prezintă un maxim de tip “neted” (fig. 1 b), altele au un vârf
ascuţit (fig. 1 d), altele prezintă mai multe maxime de valori diferite (fig. 1 c) sau nu prezintă nici
un maxim dacă coordonatele nu sunt supuse la restricţii (cazul din fig. 1 a, corespunzător la
cu α1, α2 =const)
O situaţie cu totul diferită se întâlneşte dacă se impun unele restricţii asupra variabilelor de stare – de tip egalitate:
(2.1.20)
sau, sub forma vectorială:
(2.1.21)
8
ori de tip inegalitate:
(2.1.22)
Astfel, de exemplu, suprafata H din fig. 1 b are – în lipsa restricţiilor – un maximum în punctul
P1. Dacă adoptăm însă restricţia care impune ca şi să ia valori numai pe o
suprafaţă ce intersectează suprafaţa H conform conturului zonei haşurate din fig. 1 e, atunci rezultă un maxim în condiţii de restricţii, care este punctul P2 de tangentă între curba H0 de pe
suprafaţta H şi suprafaţa. .
a) O primă situaţie care poate apare este aceea a unei funcţii de performanţă H de tip liniar, adică de forma:
(2.1.23)
iar restricţiile – de tip egalitate sau tip inegalitate – sunt de asemenea liniare. O astfel de problema reprezintă o problema de “programare liniară”.
b) Un caz important în practică este acela când funcţia de performanţă H posedă un maxim de tip “neted” (numit şi “calculabil”), ca de exemplu în fig. 1 b.
În caz că nu există restricţii, punctul de optim se poate identifica prin faptul că valoarea
lui H nu se modifică la o deplasare arbitrar de mică pe suprafaţa adică
este staţionară. Din (2.1.19) rezultă:
(2.1.24)
Întrucât sunt libere (deplasarea din punctul de optim putându-se face în orice direcţie cu
orice valoare mică), condiţia:
(2.1.25)
implică
(2.1.26)
Pentru a cunoaşte dacă optimul obţinut reprezintă un maxim, un minim sau un punct “în şa”, trebuie utilizate alte mijloace; deci condiţia (2.1.26) este doar necesară (nu şi suficientă).
Mult mai complicat este însă cazul determinării unui optim în prezenţa restricţiilor, în cazul n-dimensional, cu p restricţii de tip (2.1.20). În acest caz avem nevoie de p multiplicatori Lagrange
9
(2.1.27)
Se defineşte o “funcţie majorată” H*:
(2.1.28)
Maximul supus la restricţii de tip (2.1.20) va fi caracterizat prin
(2.1.29)
Exemplu:
Reluăm exemplul iniţial de optimizare statică a repartiţiei puterilor între generatoarele unui sistem. În acest caz funcţia H este dată de (1.8), iar restricţia de forma (1.7). Figura 2 a)
reprezintă o curba tipică [lei/h](x0i[kW]). Pentru cazul a două generatoare, suprafaţa H este
reprezentată în figura 2 b).
Figura 2 Curbe tipice de cost în raport cu puterea electrică produsă în centrale termoelectrice
Funcţia de performanţă (cost) majorată va fi
şi, conform (2.1.29):
10
respectiv
Derivatele parţiale reprezintă aşa-numitele “costuri incrementale de producţie”
ale generatoarelor.
Rezultatul obţinut în acest exemplu se traduce astfel: pentru realizarea unui beneficiu maxim, sistemul va trebui condus astfel încât costurile incrementale să fie egale pentru toate generatoarele din sistem (ceea ce nu implică neapărat ca şi costul total al producţiei de energie este egal pentru toate generatoarele) şi, fireşte, suma puterilor generate în sistem să fie egală cu puterea P a consumatorilor.
În concluzie, o astfel de optimizare se spune că este “în circuit deschis”, întrucât optimul rezultat este corect atâta timp cât modelul matematic al procesului este corect şi acest model este complet precalculat şi fix, fără a ţine cont de perturbaţiile dinamice (ci doar de cele staţionare, modelul însuşi fiind staţionar = fix).
2.2 OPTIMIZARE STATICĂ ÎN CIRCUIT ÎNCHIS
Dacă nu dispunem de un model matematic acceptabil pentru procesul ce trebuie optimizat static, este posibil a-l face pe acesta din urmă să funcţioneze optim prin experimentare în timp real. În literatura de specialitate au fost propuse numeroase metode de căutare experimentală a optimului static, metode încadrate în denumirea de “tehnici de urcarea pantei”. Toate aceste metode procedează mai întâi la căutarea direcţiei către vârful pantei şi apoi produc deplasarea în acea direcţie. Unele metode sofisticate caută nu numai direcţia către vârful pantei, ci direcţia cea mai bună în sensul atingerii optimului în timp minim; astfel de metode reprezintă “metode de gradient sau metode de urcare (respectiv coborâre, în cazul că optimul reprezintă un minim) pe linia de cea mai mare pantă”.
Împărţind (2.1.24) în ambii membri cu dt, se obţine:
(2.2.30)
care se mai poate scrie compact sub forma produsului scalar
(2.2.31)
unde
11
(2.2.32)
reprezintă gradientul funcţiei scalare H(x0) în raport cu vectorul x0.
Pentru o valoare dată a “vitezei” x0, viteza de variaţie a lui H va fi maximă atunci când x0
şi sunt paralele, adică atunci când mişcarea se face în direcţia gradientului.
Figura 3 Structura generală a unui sistem de căutare a optimului static
Un sistem automat cu optimizare statică prin căutarea experimentală va trebui să folosească semnale speciale de testare (perturbaţii de valoare mică şi lent variabile), injectate în sistem; trebuie însă menţionat că, dacă sistemul conţine în mod normal zgomote interne de valoare suficient de mare, se pot chiar ele folosi drept semnale de test. Schema din figura 3 reprezintă structura generală a unui sistem automat de căutare a optimului, în care semnalele perturbatoare de testare sunt introduse aditiv peste cele m intrări staţionare u0
i. Se presupune că toate variabilele de stare care afectează valoarea funcţiei H sunt accesibile. Se foloseşte un generator al funcţiei H pe baza valorilor rezultate din proces pentru variabilele de stare. Din cauza variaţiilor – lente – ale perturbaţiilor, variabilele de stare îşi modifică valorile şi H îşi va modifica valoarea. Aceste modificări vor fi interpretate de dispozitivul de optimizare static, care va efectua modificări ale valorilor comenzilor u0
i astfel încât valoarea lui H să se apropie de optim. Optimizarea se face în circuit închis şi, dacă dispozitivul de optimizare static este conceput corespunzător, sistemul va evolua în sensul atingerii optimului.
12
3. OPTIMIZAREA DINAMICĂ
Problema optimizării dinamice este caracterizată de vectori de stare şi de comandă variabili în timp. Altfel spus, în cazul optimizării statice avem de optimizat o funcţie (2.1.19); în cazul optimizării dinamice trebuie căutate condiţiile în care se extremizează o funţională de tip (1.2), adică
(3.33)
În cele ce urmează vom prezenta metode de optimizare dinamică.
3.1 OPTIMIZAREA PRIN METODE VARIAŢIONALE
Tot aşa cum “metoda multiplicatorilor lui Lagrange” permitea – în cadrul optimizării statice – determinarea condiţiei necesare pentru ca funcţia H – posedând un extrem de tip “neted” – să fie optimizată, în optimizarea dinamică calculul variaţional ne oferă posibilitatea determinării condiţiilor necesare (dar nu suficiente) de obţinere a unui optim pentru o funcţională (integrală) prezentând un extrem de tip “neted”.
Calculul variaţional în contextul teoriei sistemelor prezintă următoarea problemă generală (problema Bolza):
a) Procesul este descris de vectorul de stare x(t) satisfăcând ecuaţia de stare (1.1), cu
, unde funcţiile (componentele lui f) sunt continue şi
admit derivate parţiale continue într-un domeniu al spaţiului - ce va fi precizat prin
restricţii. Se acceptă, de asemenea, că sistemul (1.1) admite o soluţie unică pentru condiţii iniţiale
date şi o comandă u(t) dată.
b) Stările şi valorile comenzilor sunt supuse unor restricţii de tip
(3.1.37)
unde funcţiile ri (componentelele lui r) sunt presupuse continue şi continuu derivabile, restricţii ce definesc un domeniu (în general închis) al aşa-numitelor stări şi comenzi “admisibile” (cazul
este de acest tip, scriindu-se sub forma: ).
Un alt tip de restricţii au forma unei integrale
(3.1.38)
13
(un exemplu îl constituie limitarea energiei disponibile în sistem, exprimată prin si
care se mai poate scrie sub forma:
.
c) Sistemul evoluează între momentele t0 şi tf; acestea, ca şi x(t0), x(tf), pot fi fixate aprioric, pot fi libere, ori pot fi obligate să raspundă unor condiţii la limită (pe frontieră) de tip
(3.1.39)
unde ki şi li sunt presupuse continue şi continuu derivabile.
d) Se impune un criteriu tip Bolza – mai general ca (1.2):
(3.1.40)
unde primul termen este o integrală evaluată de-a lungul curbei integrale ce rezultă pentru x(t0) şi u(t) aleşi, iar cel de-al doilea termen este dependent de extremităţi în cazul când acestea nu sunt impuse; funcţiile L şi M sunt presupuse continue şi continuu derivabile.
În cele ce urmează vom prezenta condiţiile necesare de optimalitate şi vom enunţa teorema de optimaliate.
Fie un sistem cu ecuaţia de stare de forma ecuaţiei (1.1) cu functionala de cost dată de ecuaţia (3.1.40) unde:
1. f este o funcţie măsurabilă în raport cu t, de clasă C1 în raport cu x şi de clasă C1 în raport cu u;
2. L funcţie de clasă C1 în raport cu toate argumentele şi pozitivă (Lagrangean)
3. M funcţie de clasă C1, convexă (Meyerian)
Definim mulţimile: mulţimea comenzilor admisibile şi
mulţimea starilor admisibile.
Definim mulţimea traiectoriilor (perechilor) admisibile ca fiind
.
14
Problema este să se determine o traiectorie optimală, adicăo traiectorie astfel
încât .
Fie o traiectorie optimală.
Considerăm o variaţie mică a traiectoriei optimale . Atunci:
Definim funcţia:
cu
proprietăţile:
1. F este derivabilă;
2. . Rezultă că este un punct de minim
pentru funcţia F. Conform teoremei lui Fermat rezultă că (3.1.41).
Derivăm funcţia :
,
unde derivatele parţiale sunt calculate pentru .
Conform relaţiei (3.1.41) rezultă că:
. (3.1.42)
Definim Hamiltonianul acestei probleme:
, unde vectorul p este de forma
, iar p1…pn sunt multiplicatorii lui Lagrange.
Atunci .
Derivăm această ecuaţie după x şi după u pentru a afla şi .
Calculăm mai întâi:
15
Rezultă că:
.
Analog obţinem:
.
Astfel am obţinut următoarele formule pentru si :
.
Ecuaţia (3.1.42) devine:
(3.1.43)
Pentru a calcula integrala notăm şi
. Înlocuim în integrala I şi obţinem:
16
Conform ecuaţiei (1.1) . Astfel . Cum obţinem că:
.
Înlocuim în integrala I: . Integrăm prin părţi şi obţinem
.
Înlocuind în ecuaţia (3.1.43) cele obţinute mai sus aceasta devine:
Ţinând cont de rezultatele obţinute mai sus enunţăm teorema de optimalitate:
Dacă este o pereche optimală a sistemului dat prin ecuaţia de stare
cu funcţională de cost
atunci există multiplicatorii Lagrange astfel încât verifică următoarele
condiţii:
În cele ce urmează considerăm un caz particular de problemă de optimizare, respectiv problema liniar patratică. Astfel:
Considerăm un sistem liniar :
(3.1.44)
cu o funcţională de cost patratică
17
(3.1.45)
unde si este o matrice nenegativ
definităţ simetricăţ reală , , este o matrice
simetrică nenegativ definită şi este o matrice simetrică pozitiv definită
.
Aplicăm teorema de optimalitate:
1. Scriem Hamiltonianul
(3.1.46)
2. Scriem sistemul canonic
(3.1.47)
3. Determinăm comanda optimală
(3.1.48)
de aici
(3.1.49)
( există din faptul că R>0).
Verificăm dacă realizează minimum: Hessiana este , deci realizează
minimul problemei.
4. Găsim traiectoria optimă, introducând comanda optimală (3.1.49) în (3.1.47)
Obţinem sistemul:
(3.1.50)
18
(3.1.51)
sau
unde matricea Hamilton este de forma
(3.1.52)
5. Scriem condiţia de transversalitate (cu starea iniţială fixată şi cu starea finală liberă)
(3.1.53)
Rezolvăm sistemul
(3.1.54)
pentru a obţine perechea .
O metodă de a rezolva acest sistem este cea în care ne folosi de ecuaţia matricială Riccati.
Astfel căutăm o matrice P(t) astfel încât vectorul p să fie de forma cu condiţia
finală , pentru a obţine comanda optimală .
În acest caz formula (3.1.50) devine
(3.1.55)
şi vom avea
(3.1.56)
De aici - prin înlocuirea lui - obţinem
(3.1.57)
Rezultă că P trebuie să satisfacă ecuaţia Riccati (RDE) de forma
19
(3.1.58)
cu condiţia finală
(3.1.59)
Astfel am demostrat că dacă este comanda optimală atunci matricea P(t)
este soluţia ecuaţiei Riccati (3.1.58) cu condiţia finală (3.1.59).
Acum facem demostratia în sens invers. Presupunem că P(t) este soluţie a ecuaţiei
(3.1.58) cu condiţia finală (3.1.59). Arătăm că realizează minimul funcţionalei.
Calculăm în două moduri integrala .
Metoda 1:
(3.1.60)
Metoda 2
(3.1.61)
După ce facem toate reducerile posibile obţinem
(3.1.62)
Înmulţim cu ½ ecuaţiile (3.1.60) şi (3.1.62) şi le egalăm.
Astfel funcţionala de cost (3.1.45) se rescrie:
(3.1.63)
20
Deoarece rezultă că integrantul este minim când este 0. Acest lucru se obţine
pentru . Astfel comanda optimală este de forma .
Am demostrat astfel urmatoarea teoremă cunoscută şi sub numele de sinteza Kalman-Letov :
Problema liniar patratică (3.1.44)-(3.1.45) are soluţia optimală unică ,
unde matricea reală simetrică P este soluţia ecuaţiei diferenţiale Riccati cu condiţia finală
.
În cazul în care timpul final tf tinde spre infinit ecuaţia matricială Riccati devine:
În cele ce urmează voi prezenta soluţtia ecuaţiei matriciale Riccati.
Fie matricea fundamentală a sistemului
cu condiţia finală
Soluţia ecuaţiei diferenţiale Riccati
este
(3.1.64)
Demonstratie:
Φ şi Ψ verifică ecuaţiile
(3.1.65)
Prin derivarea egalităţii (3.1.66) obţinem
(3.1.67)
Ţinând cont de ecuaţiile (3.1.65) şi (3.1.66), ecuaţia (3.1.67) devine:
21
(3.1.68)
Dacă înmulţim la dreapta cu , vom obţine
. De aici unica soluţie a ecuaţiei matriciale Riccati este
.
Demonstraţia fiind încheiată trecem la următorul pas.
Împărţim ecuaţia fundamentală corespunzătoare matricei Hamilton
astfel
(3.1.69)
Astfel soluţia problemei Cauchy
(3.1.70)
(3.1.71)
este
(3.1.72)
Soluţia (3.1.72) poate fi scrisă desfăsurat astfel:
(3.1.73)
Ţinând cont de cele obţinute mai sus soluţia ecuaţiei matriciale Riccati poate fi calculată astfel
(3.1.74)
iar comanda optimală are forma
(3.1.75)
22
unde
(3.1.76)
Pentru a calcula soluţia ecuaţiei Riccati putem folosi algoritmul Potter [4] care constă în următoarele:
1. Se formează matricea
(3.1.77)
2. Se calculează valorile proprii matricei M. Se demonstrează că aceste valori proprii au
următoarea proprietate de simetrie: dacă λ este o valoare proprie, atunci (unde λ*
este conjugatul complex) sunt şi ele valori proprii.
Ca urmare, n valori proprii sunt situate în semiplanul stâng şi n în semiplanul drept.
Notăm cu valorile proprii din semiplanul stâng, astfel încât:
(3.1.78)
3. Se calculează vectorii prorpii asociaţi valorilor proprii
Aşadar:
(3.1.79)
4. Se partiţinonează vectorii proprii în doi vectori de dimensiune n, notaţi astfel:
(3.1.80)
5. Se formează matricea
(3.1.81)
6. Se formează matricea
(3.1.82)
7. Soluţia ecuaţiei algebrice Riccati, unde R este simetrică şi pozitiv definită, rezultă:
(3.1.83)
Pentru a înţelege mai bine acest algoritm am luat următorul exemplu numeric:
23
Formăm matricea M:
Ecuaţia caracteristică este , astfel că obţinem următoarele valori proprii:
Cele două valori proprii din semiplanul stâng sunt
Calculăm vectorii proprii asociaţi acestor valori proprii:
Un minor caracteristic nenul de ordin maxim este
Atunci este necunoscută secundară iar sistemul devine:
Rezolvând sistemul obţinem vectorul:
24
Un minor caracteristic nenul de ordin maxim este:
Atunci este necunoscută secundară iar sistemul devine:
Rezolvând sistemul obţinem vectorul:
Formăm matricile S si N:
Calculăm
Soluţia ecuaţiei Riccati este:
Matricea de reglaj este
25
În raport cu metodele iterative, acest algoritm prezintă avantajul unui calcul direct, cu soluţie unică.
3.2 PROGRAMAREA DINAMICĂ
După cum este denumită, această metodă reprezintă o metodă ingenioasă de programare pe calculator a unor probleme de optimizare, fiind introdusă de Bellman. Întrucât un calculator numeric acceptă doar date discrete, este necesar ca, pentru aplicarea metodei, să se discretizeze procesul, în caz că acesta reprezintă un proces continuu. Astfel discretizat, el devine un aşa-numit proces secvenţial de decizie, în care, la momente discrete de timp k – numite şi etape -se elaborează câte o decizie în sensul alegerii unei stări xv (valori discrete alcătuind mulţimea
stărilor posibile), conform cu utilitatea (numită şi increment de cost) asociată pasului v.
Numărul total de etape N este numit orizont. Numărul cardinal al mulţimii stărilor posibile este notat cu r şi defineşte puterea de rezoluţie în procesul secvenţial de decizie.
În procesul continuu, trecerea de la starea iniţială x0 la o stare finală xf se putea face printr-un număr infinit de traiectorii (fig. 4 a). Prin optimizare se alegea numai acea traiectorie care minimizează criteriul de performanţă.
În cadrul procesului de decizie secvenţială, la fiecare dintre paşii v=0,1,…,N-2 se elaborează câte o decizie; o vom nota cu qv (la pasul N-1 nu mai avem posibilitatea de opţiune,
căci din starea în care ne aflăm trebuie sa trecem neapărat în starea finală impusă - fig.
4 b). Deciziile qv se pot aranja într-un vector de decizie care determină în mod unic traiectoria urmată între starea iniţială şi cea finală.
(3.2.84)
26
Figura 4 Proces secvenţial de decizie
Echivalent criteriului de performanţă de tip integral (1.2), utilizat la sistemele continue, îl constituie în cazul sistemelor discrete:
(3.2.85)
cu:
(3.2.86)
unde xv şi uv sunt valorile medii ale lui x şi u în pasul respectiv. În acest caz funcţia J se numeşte funcţie obiectiv. Traiectoria care optimizează funcţia obiectiv se numeşte strategie optimă.
O metodă “simplă” de a alege strategia optimă ar consta în calcularea funcţiei obiectiv J de-a lungul tututror traiectoriilor posibile. Din punct de vedere al preciziei apare firească alegerea unor valori r şi N cât mai mari. Efectuând pasul 1, avem r decizii posibile pentru un proces unidimensional (o singură variabilă de stare). În pasul al doilea vor exista r2 decizii posibile. Întregul proces de decizie secvenţială va implica rN-1 traiectorii posibile. Dacă adoptăm valorile r=10 şi N=10 (ceea ce nu oferă o prea mare precizie), vom avea 109 traiectorii posibile. Presupunând că avem un calculator care poate efectua calcularea a unui million de valori J pe secundă, calculul valorilor J pentru toate traiectoriile posibile va dura cam 1000 secunde (acceptabil). Dar dacă avem, în aceleaşi condiţii, un proces bidimensional numărul de traiectorii
posibile va fi , adică 1018 şi vor necesita un timp de calcul al tuturor valorilor lui
27
J de ordinal a 1012secunde, adică aproximativ 32000 ani. Iar pentru un proces n dimensional, unde rezultă
(3.2.87)
traiectorii, nici numai poate fi vorba de găsirea unei traiectorii optime.
Avantajul esenţial al metodei programării dinamice este adus de “principiul optimalităţii” şi “metoda scufundării”, care permit reducerea procesului de decizie în N paşi la un proces secvenţial de N procese de decizie cu un singur pas.
Principiul optimalităţii spune că: o strategie care este optimă pe intervalul este
în mod necesar optimă si pe orice subinterval cu .
Metoda “scufundării” constă într-o procedură de descompunere a problemei de rezolvat într-o familie de probleme de aceeşi natură, dar mai simple şi deci mai uşor de rezolvat în paşi succesivi (prin iteraţii); soluţia obtinută pentru intreaga clasă de probleme oferă totodată soluţia problemei originale.
La fiecare pas trebuie efectuate comparaţii – în număr de r la primul şi la ultimul pas, iar la paşii intermediari câte r2 comparaţii (în cazul unui proces bidimensional). Aşadar, numărul
total de comparaţii va fi egal cu . Dacă am fi avut de-a face cu un proces n-
dimensional, ar fi rezultat
(3.2.88)
comparaţii, adică un număr de ordinul . Astfel, pentru n=2, r=10 şi N=10 trebuie
efectuate circa operaţii, necesitând un timp de numai 0.1 secunde pentru calcul
(faţă de 32000 de ani). Această reducere este datorată aplicării principiului optimalităţii, care elimină la fiecare pas deciziile ce nu pot candida în continuare pentru strategia optimă.
Dacă în rezolvarea prin metoda calculului variaţional a problemelor de optimizare dinamică, existenţa restricţiilor îngreuna obţinerea soluţiei, în programarea dinamică aceste restricţii sunt avantajoase căci reduc numărul soluţiilor premise în procesul de căutare a strategiei optime.
Fie comanda optimă care, pentru condiţiile pe frontieră x0 şi xf date, dă o valoare
minimă
(3.2.89)
unde J este indicele de performanţă (1.2). Vom remarca că J0 este unic, deşi u0 şi x0 pot să nu fie unice; totuşi, în practică, u0 şi x0 sunt de obicei unice, iar atunci când nu e aşa, u0 poate fi ales în mod arbitrar dintre soluţiile disponibile, pe baza unor considerente legate de aplicaţia respectivă.
28
Să presupunem că a fost obţinută. Rezultă aşadar că J0 nu depinde de u(t), ci doar de
x0, t0, xf şi tf, adică
(3.2.90)
J0 fiind valoarea minimă a lui J, avem
(3.2.91)
Aplicând “principiul scufundării”, relaţia (3.2.91) se mai poate scrie:
(3.2.92)
Aplicarea principiului optimalităţii asupra ultimei porţiuni de traiectorie optimală (de la
la tf) ne permite să scriem (3.2.92) sub forma:
(3.2.93)
Dezvoltând integrala din (3.2.93) în serie Taylor în jurul lui t=t0, obţinem:
(3.2.94)
unde reprezintă termenii conţinând pe cu . Substituind (3.2.94) in (3.2.93)
obţinem:
(3.2.95)
Întrucât procesul de minimizare nu depinde de Δt, iar J0 nu depinde de u, relaţia (3.2.95) mai poate fi scrisă sub forma:
(3.2.96)
unde:
(3.2.97)
29
Întrucât
(3.2.98)
rezultă că expresia minimizată din (3.2.90) trebuie să fie egală cu 0. Trecând la limită această expresie pentru şi presupunând că funcţia L este astfel derivabilă încât
(3.2.99)
rezultă
(3.2.100)
care este cunoscută sub numele de “forma compactă a ecuaţiei lui Bellman” şi care reprezintă o condiţie necesară pentru existenţa unei soluţii optimale.
Dacă tf şi xf sunt constante, atunci derivata totală a lui J0 în raport cu timpul se poate exprima sub forma
(3.2.101)
Substituind acest rezultat în (3.2.100), obţinem:
(3.2.102)
unde am ţinut cont că t0 este un moment arbitrar şi că, deci, relaţia (3.2.102) trebuie să fie
satisfăcută pentru orice .
Scriind ecuaţia de stare (1.1) pe componente:
(3.2.103)
şi introducând-o în (3.2.102) rezultă:
(3.2.104)
care reprezintă o ecuaţie de tip Hamilton-Jacobi.
Dacă u nu este supus restricţiilor, atunci este necesar, în scopul obţinerii unui minim, ca
(3.2.105)
sau, ţinând cont că J0 este independent de u:
30
(3.2.106)
ecuaţie ce trebuie rezolvată simultan cu (3.2.102).
Pentru procese liniare descries de ecuaţii de stare de tip , avem
(3.2.107)
(unde bij este elementul din rândul i şi coloana j a matricei B), ceea ce face ca ecuaţiile (3.2.106) să se simplifice întrucâtva.
Observăm însă că ecuaţiile (3.2.104) şi (3.2.106) sunt ecuaţii neliniare cu derivate parţiale, rezolvarea lor reclamând utilizarea calculatoarelor numerice.
3.3 PRINCIPIUL MAXIMULUI
Există o clasă largă de probleme de tip variaţional în optimizarea sistemelor dinamice, probleme care nu pot fi soluţionate prin metoda calculului variaţional – ecuaţii de tip Euler-Lagrange.
Ele ar putea fi abordate sub forma discretă a programării dinamice, dar, această metodă reprezintă în fond o tehnică numerică care furnizează o informaţie săracă asupra naturii fizice a problemei de optimizare în ansamblu. Iar aplicarea variantei continue a principiului optimalităţii se izbeşte de dificultăţi de calcul de îndată ce sistemul nu este liniar, criteriul nu este de tip cuadratic, iar comanda suferă restricţii de tip inegalitate.
Principiul optimului, elaborat de Pontriaghin, Boltianskii şi Gamkrelidze, este similar calculului variaţional şi este strâns legat de programarea dinamică – fiind dealtfel posibilă obţinerea principiului optimalităţii din principiul maximului printr-o schimbare de variabile.
Principiul maximului se bazează pe următoarele elemente:
1. Sistemul este descris de ecuaţia de stare (3.2.103) – sau forma compactă (1.1).
2. Timpul iniţial t0 şi cel final tf sunt date.
3. Unele componente ale vectorului de stare au valori date la t0 şi tf.
4. Indicele de performanţă este de forma funcţionalei
(3.3.108)
5. Comanda u(t) este supusă restricţiilor de tip inegalitate.
31
6. Se caută strategia optimă (comanda optimală) care va transfera sistemul din
starea initială în starea finală de asemenea manieră încât să optimizeze
criteriul J dat de (3.3.108), respectând restricţiile impuse comenzii.
Coeficienţii de pondere αi din (3.3.108) pot fi priviţi ca fiind componentele unui vector de ponderare – numit şi vector obiectiv:
(3.3.109)
astfel încât criteriul (3.3.108) se poate exprima sub forma unui produs scalar
(3.3.110)
(Fară a impieta asupra generalităţii, putem considera că α este normalizat, adică
).
Expresia (3.3.110) are semnificaţia geometrică a proiecţiei vectorului pe vectorul α;
aşadar optimizarea criteriului J (3.3.108) se traduce ca încercarea de a optimiza proiecţia lui
pe direcţia vectorului obiectiv.
Vom considera mai întâi cazul:
a) tf este fixat, dar este liber.
Vom presupune că strategia optimă a fost găsită şi că traiectoria optimală
corespunzătoare se cunoaşte. Vom studia modificarea dJ a criteriului de performanţă cauzată de
o variaţie de la strategia optimă. O astfel de situaţie, pentru un caz bidimensional, este
ilustrată în figura 5.
32
Figura 5. Abaterea de la strategia optimă şi abaterea rezultantă de la traiectoria optimală
Întrucât comanda optimă este de tip “bang-bang”, adică având valori extreme între care se face comutaţia (teoretic) instantanee (în fig. 5 b se arată o astfel de comutaţie), este firesc a
presupune că mica abatere de la strategia optimă, adică va reprezenta doar o
modificare redusă a momentului comutaţiei valorii comenzii. Aşadar, variaţia va consta
din unul sau mai multe impulsuri de lăţime finite dar foarte mică (un astfel de impuls este ilustrat în fig. 5 c). Se poate formula diferenţa dintre principiul maximului şi metoda calculului variaţional astfel: dacă vectorul m-dimensional u(t) este supus restricţiei
(3.3.111)
unde {U} este domeniul comenzilor admisibile în spaţiul m-dimensional, atunci metoda
calculului variaţional caută soluţia optimă în interiorul lui {U}, în timp ce principiul
maximului caută soluţia optimă pe frontiera lui {U}.
Introducem conceptul de vector adjunct de stare p(t) ale cărui componente sunt definite astfel:
(3.3.112)1
(3.3.112)2
33
În cazul sistemelor liniare, având ecuaţia de stare , relaţia (3.3.112)1
ia forma:
(3.3.113)
care defineşte aşa-numitul sistem adjunct al sistemului cu ecuaţia de stare de tipul
.
Vectorul de stare x şi vectorul adjunct de stare p introduce un număr total de 2n variabile, care trebuie să satisfacă ecuaţiile diferenţiale ordinare de ordinul unu (3.3.103) şi (3.3.112)1; în plus, ele trebuie să satisfacă cele 2n condiţii terminale specificate, în maniera “sparte în două puncte”:
şi .
Analizăm variaţia componentei xi a vectorului x cauzată de variaţia mai sus menţionată Δu a lui u. Conform cu (3.3.103) avem:
(3.3.113)
Multiplicând ambii membri ai relaţiei (3.3.113) cu pi, însumând relaţia astfel obţinută pentru toate cele n componente ale lui x şi, în final, integrând între t0 şi tf se obţine:
(3.3.114)
Membrul stâng al relaţiei (3.3.114) poate fi integrat prin părţi:
(3.3.115)
Întrucât, fireşte, avem:
(3.3.116)
şi ţinând cont de (3.3.112)2, vom avea:
(3.3.117)
Dar memebrul drept din (3.3.117) reprezintă variaţia – cu semn schimbat – a criteriului J, adică:
(3.3.118)
Din (3.3.112)1, (3.3.114), (3.3.115), (3.3.117) şi (3.3.118) rezultă:
34
(3.3.119)
Factorul reprezintă o diferenţă mică (de ordinul unu)
astfel că el poate fi dezvoltat în serie Taylor:
(3.3.120)
unde R este suma termenilor de ordin superior lui 1. Introducând (3.3.120) în (3.3.119) se obţine:
(3.3.121)
Se consideră doar funcţiile fi care sunt liniare în raport cu x şi depind de u în mod aditiv, adică cele pentru care (1.1) ia forma:
(3.3.122)
Pentru această clasă de funcţii vom avea:
(3.3.123)
astfel că cea de-a doua integrală din (3.3.119) va fi nulă şi relaţia (3.3.119) se reduce la:
(3.3.124)
unde
(3.3.125)
reprezintă aşa-numita funcţie de stare a lui Pontriaghin – denumire justificată prin aceea că H depinde exclusiv de vectorul de stare x(t) – şi întâlnită uneori şi sub denumirea de “Hamiltonian” (dată fiind similaritatea sa cu expresia hamiltonianului din mecanica clasică).
Relaţia (3.3.124) se interpretează astfel: Dacă se poate dovedi că o variaţie Δu arbitrară face ca integrala din (3.3.124) să fie pozitivă (negativă), aceasta va implica faptul că dJ să fie negativă (pozitivă) pentru aceste deviatii şi u0 va corespunde atunci la un maxim (minim) al lui J.
35
Dar integrala din (3.3.124) va fi pozitivă (negativă) dacă integrandul va fi pozitiv
(negativ) pe întregul interval , adică dacă
(3.3.126)
sau
(3.3.127)
Aşadar, principiul maximului se formulează astfel:
O condiţie necesară şi suficientă pentru a avea un maxim (minim) al criteriului de performanţă J – dat de (3.3.108) – este aceea ca: comanda optimă u0 să minimizeze
(maximizeze) funcţia de stare H a lui Pontrianghin pe tot intervalul .
În expunerea metodelor de optimizare dinamică 3.1 şi 3.2 am întâlnit criterii de performanţă de tip integral (forma 1.2). Principiul maximului este aplicabil şi în aceste cazuri dacă introducem o variabilă de stare “suplimentară” definită astfel:
(3.3.128)1
unde, evident:
(3.3.128)2
Atunci, criteriul (3.3.108) devine:
(3.3.129)
şsi H ia forma:
(3.3.130)
care reprezintă forma generalizată a funcţiei de stare a lui Pontriaghin.
O observaţie importantă este aceea că: din ecuaţia reprezentând condiţia necesară de existentă a unui optim pentru indicele de performanţă J dat de (1.2) rezultă direct funcţia de stare generalizată a lui Pontriaghin ca fiind
36
(3.3.131)
dacă facem schimbarea de variabilă
(3.3.132)
(relaţie ce se mai scrie compact sub forma:
(3.3.133)
unde reprezintă gradientul în raport cu x).
Observăm, de asemenea, din (3.3.131) si (3.3.132), că
(3.3.134)
iar din expresia (3.3.130) a lui H şi ţinând cont de (1.1), avem:
(3.3.135)
Relaţtiile (3.3.134) şi (3.3.135), care se mai pot scrie sub forma compactă
(3.3.136)
(3.3.137)
constituie ecuaţiile Hamilton-Pontriaghin. Soluţionarea lor, pentru condiţii la limită convenabile, furnizează comanda optimă.
Dacă u nu este limitat (supus restricţiilor de tip inegalitate), atunci condiţia necesară de extremizare a lui H prin u este:
(3.3.138)
Dacă u este însă supus unei restricţii, exprimabilă sub forma (3.3.111), atunci condiţia necesară pentru extremizarea lui H se exprimă prin relaţia generală
sau (3.3.139)
unde
(3.3.140)
37
b) Dacă x(tf) este parţial impus, adică doar q componente
q<n (3.3.141)
ale vectorului de stare sunt supuse unor restricţii de tip egalitate la tf – cu tf fixat – adică
(3.3.142)
Celelalte n-q componente sunt “libere” şi atunci criteriul (3.3.108) va suferi o variaţie
(3.3.143)
dată doar de variabilele de stare care nu sunt supuse restricţiilor (se spune că “am pierdut q grade de libertate”). Cu acestea, cele expuse la punctual a) rămân valabile, cu condiţia ca cele q variabile de stare adjuncte pv(t) ce corespund variabilelor de stare xv(t) supuse restricţiilor să nu mai fie obligate a avea valorile finale αv, ci să fie lăsate libere (compensăm astfel cele q grade de libertate pierdute prin ridicarea a q restricţii).
4 Aplicaţii
4.1 Aplicaţie optimizare prin metode variaţionale utilizând teorema de optimalitate [3]
Ne punem problema să găsim traiectoria ce trebuie urmată de un satelit pentru a ajunge dintr-o stare iniţială dată într-o stare finală dată cu consum minim de energie. Adică vom afla comanda
optimală uopt(t) ce minimizează indicele de performanţă .
Legea de mişcare a satelitului este dată prin sistem de ecuaţii diferenţiale de oridin 2 de mai jos:
unde u1 este forţa ce acţionează asupra satelitului pe direcţie radială, iar u2 este foarţa ce acţionează pe direcţie tangenţială.
Sistemul dinamic este neliniar şi poate fi tratat ca atare, dar poate fi liniarizat şi studiat ca sistem liniar, dupã cum vom vedea în cele ce urmeazã.
Observãm cã pentru sistemul de mai sus admite soluţia
, unde r0 şi ω sunt constante, iar .
38
Dacã alegem r0 = 1 deducem: , , şi . Alegând stãrile
x1, x2, x3 şi x4 ca fiind: , , şi
; mai precis , , şi
partea neliniarã a sistemului devine:
iar partea liniarã:
cu soluţia iniţialã x10 = x20 = x30 = x40 = 0 .
Prin urmare
Calculând coeficienţii de mai sus corespunzători soluţiei iniţiale obţinem ecuaţiile de stare ale sistemului liniarizat:
Astfel am adus sistemul sub forma:
unde:
39
Avand sistemul care descrie legea de mişcare a satelitului scris sub forma de mai sus facem următoarele notaţii:
1. - vectorul multiplicatorilor lui Lagrange;
2. ,
3. - Lagrangeanul
Folosindu-ne de aceste notaţii scriem Hamiltonianul:
.
Pentru a determina vectorul p(t) al multiplicatorilor lui Lagrange rezolvăm următorul sistem de ecuaţii diferenţiale:
având condiţiile iniţiale: .
După ce aflăm vectorul p(t), pentru a determina comanda optimală uop(t) rezolvăm sistemul:
40
Pentru a afla traiectoria optimală nu ne mai rămâne decât să înlocuim u(t) cu uop(t) în ecuaţia de stare a sistemului. Astfel obţinem xop(t) rezolvând:
având condiţiile iniţiale . Astfel obţinem coordonatele polare ce
descriu traiectoria optimală a satelitului de forma:
Pentru a trece la coordonate carteziene nu ne rămâne decât să facem următoarele transformări :
Rezultate şi concluzii
Pentru a rezolva numeric această problemă am folosit atât programul Maple cât şi programul Matlab. Programele propriu-zise se găsec în Anexa 1, respectiv Anexa 2.
În calculele efectuate am considerat .
În urma calculelor am obţinut comanda optimală de forma:
Traiectoria ce trebuie urmată de satelit este dată prin coordonatele polare de mai jos:
41
a) b)
Figura 6 a) traiectoria optimală reprezentată folosindu-ne de coordonatele carteziene x şi y;
b) traiectoria optimală reprezentată folosindu-ne de coordonatele polare r şi θ
În figura de mai sus am reprezentat grafic traiectoria ce trebuie urmată de satelit, folosindu-ne atât de coordonatele carteziene ( a) ) cât şi de coordonatele polare ( b) ).
4.2. Aplicaţie optimizare prin metode variaţionale utilizând ecuaţia matricială Riccati [5]
Ştim că marea majoritate a rachetelor ghidate au nevoie ca unghiul de ruliu să fie menţinut constant pe toată durata zborului pentru ca sistemul de ghidare să funcţioneze corect. Din acest motiv este necesar un autopilot de ruliu care să menţină constantă poziţia rachetei în timpul zborului.
Astfel noi vom construi un sistem de comandă cu feedback care va ţine unghiul de ruliul aproape de zero. În acelaşi timp mişcarea eleronului nu va depăşi un anumit prag dat pentru unghiul de deviere.
Ecuaţiile de mişcare pentru mişcarea de ruliu a rachetei sunt:
42
unde este momentul de ruliu datorat mişcării eleroanelor; iar este momentul de
amortizare la mişcarea de ruliu.
Făcând notaţiile: , putem rescrie ecuaţiile de mişcare astfel:
Pentru a aduce ecuaţiile de mai sus sub forma le scriem astfel:
unde
Funcţionala de minimizat în acest caz are forma:
unde este unghiul maxim de ruliu dorit, pmax este valoarea maximă admisă a coeficientului de
ruliu , iar este valoarea maximă admisă a deviaţia eleroanelor.
Comparând această expresie cu forma generală a funcţionalei de cost:
, observăm că în cazul nostru avem:
Determinăm comanda optimă rezolvând ecuaţia matricială Riccati
43
Deoarece în cazul nostru timpul final este infinit matricea S devine o matrice constant (nu
depinde de timp) şi astfel termenul dispare. Astfel noi trebuie să rezolvăm ecuaţia, având ca
necunoscute elementele matricii S:
Înlocuind în ecuaţia de mai sus matricile A, B, Q şi R obţinem următorul sistem de ecuaţii, cu necunoscutele S11, S12, S21 şi S22. Dar cum S este o matrice simetrică avem egalitatea S12,=S21.
Comanda optimală este dată de ecuaţia:
unde .
Rezultate si concluzii
Am folosit următoarele valori numerice pentru caracteristicile aerodinamice ale rachetei:
Am rezolvat numeric această aplicaţie folosind atât Maple cât si Matlab. Programele propriu-zise se găsec în Anexa 3, respectiv Anexa 4.
În urma calculelor am obţinut comanda optimală de forma:
Valorile optimale pentru şi p(t) sunt:
Funcţionala de cost corespunzatoare comenzii optimale are valoarea:
44
Pentru a scoate în evidenţă că această valoare a funcţionalei este cea minima considerăm o altă valoare pentru comandă şi comparăm noul cost cu costul obţinut pentru comanda optimală.
Fie noua comandă . Costul corespunzător acestei comenzi are valoarea . Observăm că noul cost este mai mare decât cel obţinut pentru comanda optimală.
În figura de mai jos am reprezentat şi funcţie de timp. Din grafice observăm că valorile acestora tind la 0 lucru pe care îl şi doream. Scopul nostru a fost să găsim comanda optimala ce va menţine ughiul de rului cât mai aproape de valoarea 0.
a) b)
Figura 7 a) funcţie de t; b) funcţie de t
Acelaşi lucru l-am verificat calculând şi . Aceste limite ne-au dat 0, lucru care
ne confirmă că am obţinut comanda care menţine unghiul de ruliu cât mai aproape de valoarea 0.
Făcând o comparaţie între programul realiyat în Maple şi cel realiyat în Matlab putem spune că rezolvarea problemei utilizând programul Matlab este mai avantajoasă deoarece există funcţia lqr (linear-quadratic regulator) care calculează direct soluţia S a ecuaţiei matriciale Riccatri şi matricea K. Sintaxa acestei funcţii este [K,S,e]=lqr(A,B,Q,R) unde matricile A şi B sunt cele din ecuaţia de stare a sistemului, iar matricile Q şi R sunt cele din funcţionala de cost. Funcţia lqr întoarce ca rezultat matricile K şi S. Astfel ne este uşurată foarte mult munca.
4.3. Aplicaţie optimizare folosind metoda programării dinamice [6]
Fie un avion care poate zbura de la stânga la dreapta de-a lungul rutelor reprezentate în figura 8. Intersecţiile a, b, c,… reprezintă oraşe iar numerele reprezintă combustibilul necesar pentru a parcurge distanţa dintre oraşe. Vom folosi principiul programării dinamice al lui Bellman pentru a determina ruta ce trebuie urmată de avion pentru a ajunge din oraşul a în orasul i cu consum minim de combustibil.
45
Figura 8
Numerotăm nivelurile de decizie ale procesului de la până la , aşa cum se
vede în figura 9. La fiecare nivel se ia o decizie, acest lucru ne mai fiind necesar
la ultimul nivel deoarece starea finală este impusă.
Figura 9
Starea curentă este nodul in care luăm decizia curentă. Astfel, starea initială este .
La nivelul 1, starea poate fi sau . Similar ori g; sau h; iar starea
finală este constransă să fie .
Pentru a controla uk considerăm că la fiecare nivel k uk poate luam valorile +1 şi -1 după
regula: dacă ne mişcăm în sus pe diagrama rutelor iar dacă ne mişcăm în jos pe
diagrama rutelor (nivelul de referinţă fiind nivelul nodului în care se ia decizia).
Pentru a rezolva problema vom parcurge invers diagrama: de la nodul final i la nodul de start a.
46
Pornim cu . La acest nivel nu este necesară nicio decizie aşa că trecem la
nivelul . Dacă atunci controlul optimal este iar costul (adică combustibilul
necesar pentru a ajunge din nodul i in nodul f) este 4. Dacă atunci controlul optimal este
iar costul este 2.
Acum trecem la nivelul următor . Aici avem 3 stări posibile. Dacă atunci
controlul este iar costul este 3+4=7. Dacă trebuie să luăm o decizie. Dacă luam
o luăm pe ruta şi apoi ajungem în nodul i, având un cost de 3+4=7. Pe de altă
parte, dacă luăm o luăm pe ruta şi apoi ajungem în nodul i, vom avea un cost de
2+2=4. Astfel că decizia optimală este cu un cost asociat de 4. Dacă avem o
singură opţiune: cu un cost de 2+4=6.
Trecem la nivelul . Avem 2 stări posibile. Dacă trebuie să luăm o decizie.
Dacă luăm o luăm pe ruta şi apoi ajungem in nodul i, având un cost de
2+3+4=9. Dacă luăm o luăm pe ruta şi apoi continuăm pe ruta optimală aflată
mai sus , având un cost de 1+2+2=5. Astfel decizia optimală este cu un cost
asociat de 5. Şi în cazul în care trebuie luată o decizie. Dacă luăm o luăm pe ruta
şi apoi continuăm pe ruta optimală , având un cost de 3+2+2=7. Dacă luăm
o luăm pe ruta şi apoi ajungem în nodul i, având un cost de 2+4+2=8. Astfel
decizia optimală este cu un cost asociat de 7.
Ajunge astfel la nivelul . Avem o singură stare posibilă (starea iniţială). Dacă
luăm o luăm pe ruta şi apoi continuăm pe ruta optimală , având
un cost de 3+1+2+2=8. Dacă luăm o luăm pe ruta şi apoi continuăm pe ruta
optimală , având un cost de 1+3+2+2=8. Observăm că indiferent de valoarea lui
u0 obţinem acelaşi cost. Astfel controlul optimal pentru nu este unic.
Rezultate si concluzii
Conform rezultatelor noastre există două rute de la nodul a la nodul i cu acelaşi cost 8 (consumul minim de combustibil): şi . Astfel soluţia
problemei de cost minim nu este unică.
Am evidenţiat rezultatele obţinute mai sus în figura de mai jos (figura 9). Numerele dintre paranteze reprezintă costul minim necesar pentru a ajunge din nodul aosicat numarului din paranteză în nodul final i.
47
Figura 10
Observăm că principiul de optimalitate al lui Bellman a redus numărul calculelor necesare pentru a rezolva problema prin reducerea numărului deciziilor ce trebuie luate.
5. Concluzii
În această lucrare am prezentat trei metode de optimizare dinamică a proceselor : optimizarea prin metode variaţionale (metoda multiplicatorilor lui Lagrange şi metoda ecuaţiei matriciale Riccati), metoda programării dinamice a lui Bellman şi metoda principiului maximului
al lui Pontriaghin.
Pentru a întelege mai bine aceste metode am considerat trei aplicaţii.
Prima aplicaţie am rezolvat-o utilizând metoda multiplicatorilor lui Lagrange. Ne-am pus problema să găsim traiectoria ce trebuie urmată de un satelit pentru a ajunge dintr-o stare iniţială dată într-o stare finală dată cu consum minim de energie. Adică sa aflam comanda optimală
uopt(t) ce minimizează indicele de performanţă .
În urma calculelor am obţinut comanda optimală de forma:
iar traiectoria optimală de forma :
48
a) b)
Figura 11 a) traiectoria optimală reprezentată folosindu-ne de coordonatele carteziene x şi y;
b) traiectoria optimală reprezentată folosindu-ne de coordonatele polare r şi θ
În a doua aplicaţie am exemplificat metoda ecuaţiei matriciale Riccati. Noi am construit un sistem de comandă cu feedback care va ţine unghiul de ruliu al rachetei aproape de zero. Rachetele ghidate au nevoie ca unghiul de ruliu să fie menţinut constant pe toată durata zborului pentru ca sistemul de ghidare să funcţioneze corect
În urma calculelor am obţinut comanda optimală (unghiul eleroanelor) de forma:
Valorile optimale pentru (unghiul de ruliu) şi p(t) (coeficientul de ruliu) sunt:
Reprezentând grafic şi funcţie de timp observăm că valorile acestora
tind la 0 lucru pe care îl şi doream.
49
a) b)
Figura 12 a) funcţie de t; b) funcţie de t
În a treia aplicaţie am folosit principiul programării dinamice al lui Bellman pentru a determina ruta ce trebuie urmată de un avion pentru a ajunge din oraşul a în orasul i cu consum minim de combustibil.
Am observat faptul ca problema nu are soluţie unică (lucru conform cu teoria). Conform rezultatelor noastre există două rute de la nodul a la nodul i cu acelaşi cost 8 (consumul minim de combustibil).
Am evidenţiat rezultatele obţinute în figura de mai jos. Numerele dintre paranteze reprezintă costul minim necesar pentru a ajunge din nodul asociat numărului din paranteză în nodul final i.
Figura 13
50
Prin utilizarea principiul de optimalitate al lui Bellman am redus numărul calculelor necesare pentru a rezolva problema prin micşorarea numărului deciziilor ce trebuie luate.
Metodele variaţionale de optimizare se utilizează în cazul porceselor continue, având un indice de performanta de tip cuadratic (probleme de tip Bolza).
Metoda programării dinamice a lui Bellman se aplică în cazul proceselor discrete. Avantajul esenţial al metodei programării dinamice este adus de “principiul optimalităţii” şi “metoda scufundării”, care permit reducerea procesului de decizie în N paşi la un proces secvenţial de N procese de decizie cu un singur pas.
Principiul lui Pontriaghin se utilizează în cazul problemelor care nu pot fi soluţionate prin metoda calculului variaţional – de exemplu ecuaţii de tip Euler-Lagrange. Ele ar putea fi abordate sub forma discretă a programării dinamice, dar, această metodă reprezintă în fond o tehnică numerică care furnizează o informaţie săracă asupra naturii fizice a problemei de optimizare în ansamblu. Principiul optimului, elaborat de Pontriaghin, Boltianskii şi Gamkrelidze, este similar calculului variaţional şi este strâns legat de programarea dinamică.
Diferenţa dintre principiul maximului şi metoda calculului variaţional constă în faptul că :
în timp ce metoda calculului variaţional caută soluţia optimă în interiorul lui {U},
principiul maximului caută soluţia optimă pe frontiera lui {U}, unde {U} este domeniul comenzilor admisibile
În concluzie, în această lucrare am prezentat metode de rezolvare a problemei comenzii
optimale: Să se găsească comanda optimă uopt (numită şi “strategia de conducere optimală”) care optimizează indicele de performanţă J, în condiţiile unor restricţii impuse; în cazul proceselor dinamice în timp.
51
6. Anexe
Anexa 1
Codul programului în Maple pentru rezolvarea aplicaţiei 4.1
> restart;
> with(linalg):
> with(DEtools):
> with(LinearAlgebra):
> tf:=2;omega:=1;
:= tf 2
:= 1
> r0:=1;k:=r0^3*omega^2;
:= r0 1
:= k 1
> A:=Matrix([[0,1,0,0],[3*omega^2,0,0,2*omega],[0,0,0,1],[0,-2*omega,0,0]]);
:= A
0 1 0 0
3 0 0 2
0 0 0 1
0 -2 0 0
> B:=Matrix([[0,0],[1,0],[0,0],[0,1]]);
:= B
0 0
1 0
0 0
0 1
> x:=Matrix([[x1],[x2],[x3],[x4]]);
:= x
x1
x2
x3
x4
52
> u:=Matrix([[u1],[u2]]);
:= u
u1
u2
> M:=A.x+B.u;
:= M
x2
3 x1 2 x4 u1
x4
2 x2 u2
> f1:=M[1,1];
:= f1 x2
> f2:=M[2,1];
:= f2 3 x1 2 x4 u1
> f3:=M[3,1];
:= f3 x4
> f4:=M[4,1];
:= f4 2 x2 u2
> L:=Transpose(u).u;
:= L [ ]u12 u22
> p:=Matrix([[p1(t)],[p2(t)],[p3(t)],[p4(t)]]);
:= p
( )p1 t
( )p2 t
( )p3 t
( )p4 t
> f:=Matrix([[f1],[f2],[f3],[f4]]);
:= f
x2
3 x1 2 x4 u1
x4
2 x2 u2
> H:=L+transpose(p).f;
:= H [ ] u12 u22 ( )p1 t x2 ( )p2 t ( ) 3 x1 2 x4 u1 ( )p3 t x4 ( )p4 t ( ) 2 x2 u2
53
Determinăm vectorul p(t):
> sis_ode:=diff(p1(t),t)=-diff(H[1,1],x1),diff(p2(t),t)=-diff(H[1,1],x2),diff(p3(t),t)=-diff(H[1,1],x3),diff(p4(t),t)=-diff(H[1,1],x4);
sis_ode ddt
( )p1 t 3 ( )p2 t ddt
( )p2 t ( )p1 t 2 ( )p4 t ddt
( )p3 t 0, , , :=
ddt
( )p4 t 2 ( )p2 t ( )p3 t
> ics:=p1(0)=c1,p2(0)=c2,p3(0)=c3,p4(0)=c4;
:= ics , , ,( )p1 0 c1 ( )p2 0 c2 ( )p3 0 c3 ( )p4 0 c4
> s:=dsolve([sis_ode,ics]);
s ( )p4 t 2 ( )cos t ( ) c1 2 c4 2 ( )sin t ( )2 c3 c2 3 c3 t 2 c1 3 c4,{ :=
( )p1 t 3 ( )cos t ( ) c1 2 c4 3 ( )sin t ( )2 c3 c2 6 c3 t 4 c1 6 c4,
( )p3 t c3 ( )p2 t ( )sin t ( ) c1 2 c4 ( )cos t ( )2 c3 c2 2 c3, }
> p1(t):=-3*cos(omega*t)*(c1-2*c4*omega)-3*sin(omega*t)*(2*c3+c2*omega)+6*c3*omega*t+4*c1-6*c4*omega;
:= ( )p1 t 3 ( )cos t ( )c1 2 c4 3 ( )sin t ( )2 c3 c2 6 c3 t 4 c1 6 c4
> p2(t):=(-sin(omega*t)*(c1-2*c4*omega)+cos(omega*t)*(2*c3+c2*omega)-2*c3)/omega;
:= ( )p2 t ( )sin t ( )c1 2 c4 ( )cos t ( )2 c3 c2 2 c3
> p3(t):=c3;
:= ( )p3 t c3
> p4(t):= -2*cos(omega*t)/omega*(c1-2*c4*omega)-2*sin(omega*t)*(2*c3+c2*omega)/omega+3*c3*t+1/2/omega*(4*c1-6*c4*omega);
:= ( )p4 t 2 ( )cos t ( )c1 2 c4 2 ( )sin t ( )2 c3 c2 3 c3 t 2 c1 3 c4
Determinăm comanda optimală:
> ec5:=diff(H[1,1],u1)=0;
:= ec5 2 u1 ( )p2 t 0
> ec6:=diff(H[1,1],u2)=0;
54
:= ec6 2 u2 ( )p4 t 0
> s1:=solve({ec5,ec6},[u1,u2]);
s1 u1 12
( )sin t c1 ( )sin t c4 ( )cos t c312
( )cos t c2 c3,
:=
u2 ( )cos t c1 2 ( )cos t c4 2 ( )sin t c3 ( )sin t c23 c3 t
2c1
3 c42
> u1:=-1/2*(-sin(omega*t)*(c1-2*c4*omega)+cos(omega*t)*(2*c3+c2*omega)-2*c3)/omega;
:= u1 12
( )sin t ( )c1 2 c412
( )cos t ( )2 c3 c2 c3
> u2:=-1/2*(-2*cos(omega*t)*c1+4*cos(omega*t)*c4*omega-4*sin(omega*t)*c3-2*sin(omega*t)*c2*omega+3*c3*omega*t+2*c1-3*c4*omega)/omega;
:= u2 ( )cos t c1 2 ( )cos t c4 2 ( )sin t c3 ( )sin t c23 c3 t
2c1
3 c42
Comanda optimală:
> u;
12
( )sin t ( )c1 2 c412
( )cos t ( )2 c3 c2 c3
( )cos t c1 2 ( )cos t c4 2 ( )sin t c3 ( )sin t c23 c3 t
2c1
3 c42
Determinăm constantele p10
, p20
, p30
şi p40
:
> ec7:=diff(p1(t),t)=-diff(H[1,1],x1);
ec7 3 ( )sin t ( )c1 2 c4 3 ( )cos t ( )2 c3 c2 6 c3 :=
3 ( )sin t ( )c1 2 c4 3 ( )cos t ( )2 c3 c2 6 c3
> ec8:=diff(p2(t),t)=-diff(H[1,1],x2);
ec8 :=
( )cos t ( )c1 2 c4 ( )sin t ( )2 c3 c2 ( )cos t ( )c1 2 c4 ( )sin t ( )2 c3 c2
> ec9:=diff(p3(t),t)=-diff(H[1,1],x3);
:= ec9 0 0
55
> ec10:=diff(p4(t),t)=-diff(H[1,1],x4);
ec10 2 ( )sin t ( )c1 2 c4 2 ( )cos t ( )2 c3 c2 3 c3 :=
2 ( )sin t ( )c1 2 c4 2 ( )cos t ( )2 c3 c2 3 c3
> s2:=solve({ec7,ec8,ec9,ec10},[c1,c2,c3,c4]);
:= s2 [ ][ ], , ,c1 c1 c2 c2 c3 c3 c4 c4
> c1:=1;c2:=1;c3:=1;c4:=1;
:= c1 1
:= c2 1
:= c3 1
:= c4 1
Comanda optimală devine:
> u;
12
( )sin t32
( )cos t 1
( )cos t 3 ( )sin t3 t2
12
Condiţiile iniţiale:
>x10:=0;x20:=0;x30:=0;x40:=0;x0:=Matrix([[x10],[x20],[x30],[x40]]);
:= x10 0
:= x20 0
:= x30 0
:= x40 0
:= x0
0
0
0
0
56
Determinăm traiectoria optimală:
> q:=t-g;
:= q t g
> C:=exponential(A,t).x0;
:= C
0
0
0
0
> E:=exponential(A,q).B.u;
E :=
( )sin t g
12
( )sin t32
( )cos t 1
( ) 2 ( )cos t g 2
( )cos t 3 ( )sin t
3 t2
12
( )cos t g
12
( )sin t32
( )cos t 1
2 ( )sin t g
( )cos t 3 ( )sin t
3 t2
12
( )2 ( )cos t g 2
12
( )sin t32
( )cos t 1
( ) 3 t 3 g 4 ( )sin t g
( )cos t 3 ( )sin t
3 t2
12
2 ( )sin t g
12
( )sin t32
( )cos t 1
( ) 3 4 ( )cos t g
( )cos t 3 ( )sin t
3 t2
12
> x1(t):=C[1,1]+int(E[1,1],g=0..t);
( )x1 t52
( )cos t ( )sin t32
( )cos t 2 52
( )cos t 6 ( )sin t 2 9 t ( )sin t32
( )sin t 1 :=
2 ( )cos t t 3 t2 t
> x2(t):=C[2,1]+int(E[2,1],g=0..t);
57
( )x2 t :=
12
( )sin t 2 152
( )cos t ( )sin t 7 ( )sin t 2 ( )cos t 2 3 ( )cos t t 3 ( )cos t 1 3 t
> x3(t):=C[3,1]+int(E[3,1],g=0..t);
( )x3 t ( )sin t 2 15 ( )cos t ( )sin t 14 ( )sin t 4 ( )cos t 2 9 ( )cos t t 6 ( )cos t 2 :=
t ( )sin t 8 t32
( )cos t t2 92
( )sin t t2 9 t3
43 t2
4
> x4(t):=C[4,1]+int(E[4,1],g=0..t);
( )x4 t 5 ( )cos t ( )sin t 3 ( )cos t 2 5 ( )cos t 12 ( )sin t 2 15 t ( )sin t 3 ( )sin t 2 :=
3 ( )cos t t9 t2
23 t2
Traiectoria optimală (dată în coordonate polare):
> r(t):=x1(t)+1;
( )r t52
( )cos t ( )sin t32
( )cos t 2 52
( )cos t 6 ( )sin t 2 9 t ( )sin t32
( )sin t 2 :=
2 ( )cos t t 3 t2 t
> theta(t):=x3(t)+omega*t;
( ) t ( )sin t 2 15 ( )cos t ( )sin t 14 ( )sin t 4 ( )cos t 2 9 ( )cos t t 6 ( )cos t 2 :=
t ( )sin t 7 t32
( )cos t t2 92
( )sin t t2 9 t3
43 t2
4
Traiectoria optimală (dată prin coordinate carteziene):
> xc(t):=r(t)*cos(theta(t));
( )xc t52
( )cos t ( )sin t32
( )cos t 2 52
( )cos t 6 ( )sin t 2 9 t ( )sin t32
( )sin t 2 :=
2 ( )cos t t 3 t2 t cos ( )sin t 2 15 ( )cos t ( )sin t 14 ( )sin t 4 ( )cos t 2
9 ( )cos t t 6 ( )cos t 2 t ( )sin t 7 t32
( )cos t t2 92
( )sin t t2 9 t3
43 t2
4
58
> yc(t):=r(t)*sin(theta(t));
( )yc t52
( )cos t ( )sin t32
( )cos t 2 52
( )cos t 6 ( )sin t 2 9 t ( )sin t32
( )sin t 2 :=
2 ( )cos t t 3 t2 t sin ( )sin t 2 15 ( )cos t ( )sin t 14 ( )sin t 4 ( )cos t 2
9 ( )cos t t 6 ( )cos t 2 t ( )sin t 7 t32
( )cos t t2 92
( )sin t t2 9 t3
43 t2
4
Costul corespunzător comenzii optimale:
> J:=evalf(simplify(1/2*int(L[1,1],t=0..tf)));
:= J 1.001491393
Reprezentarea grafică a traiectoriei:
> with(plots):
> spacecurve([xc(t),yc(t),t],t=0..3,axes=FRAME);
> polarplot([r(t),theta(t),t=0..3]);
59
Anexa 2
Codul programului în Matlab pentru rezolvarea aplicaţtiei 4.1
clear all;
r0=1
omega=1
syms omega;syms x1;syms x2;syms x3;syms x4;
syms u1;syms u2;syms t;syms g;syms q;
syms p1;syms p2;syms p3;syms p4;
syms p10;syms p20;syms p30;syms p40;
k=r0^3*omega^2
A=[0 1 0 0; 3*omega^2 0 0 2*omega; 0 0 0 1; 0 -2*omega 0 0]
A = 0 1 0 0
3 0 0 2
0 0 0 1
0 -2 0 0
B=[0 0; 1 0; 0 0; 0 1]
B = 0 0
1 0
60
0 0
0 1
x=[x1;x2;x3;x4]
x = x1
x2
x3
x4
Condiţiile iniţiale:
x10=0;x20=0;x30=0;x40=0;
x0=[x10;x20;x30;x40]
u=[u1;u2]
u = u1
u2
M=A*x+B*u
M = x2
3*x1+2*x4+u1
x4
-2*x2+u2
f1=M(1,1)
f1 = x2
f2=M(2,1)
f2 = 3*x1+2*x4+u1
f3=M(3,1)
f3 = x4
f4=M(4,1)
f4 = -2*x2+u2
f=[f1;f2;f3;f4]
61
p=[p1;p2;p3;p4]
p = p1
p2
p3
p4
p0=[p10;p20;p30;p40]
L=u.'*u
L = u1^2+u2^2
H=L+p.'*f
H = u1^2+u2^2+p1*x2+p2*(3*x1+2*x4+u1)+p3*x4+p4*(-2*x2+u2)
Determinăm vectorul p(t):
diff(H,x1)
3*p2
diff(H,x2)
p1-2*p4
diff(H,x3)
0
diff(H,x4)
2*p2 +p3
s=dsolve ('Dp1=-3*p2*omega^2','Dp2=-(p1-2*p4*omega)','Dp3=0','Dp4=-(2*p2*omega+p3)',’p1 (0) =p10',’p2 (0) =p20',’p3 (0) =p30',’p4 (0) =p40')
p1=eval(s.p1)
p1 = -3*cos(t)*(p10-2*p40)-3*sin(t)*(2*p30+p20)+6*p30*t+4*p10-6*p40
p2=eval(s.p2)
p2 = (sin (t)*(-p10+2*p40)+cos (t)*(2*p30+p20)-2*p30)
p3=eval(s.p3)
p3 = p30
p4=eval(s.p4)
p4 = (-2*p10+4*p40)*cos(t)-(4*p30+2*p20)*sin(t)+3*p30*t+2*p10-3*p40
62
Determinăm comanda optimală:
diff (H,u1)
2*u1+p2
diff (H,u2)
2*u2+p4
u1=eval (solve ('2*u1+p2=0','u1'))
u1 = 1/2*sin(t)*(p10-2*p40)-1/2*cos(t)*(2*p30+p20)+p30
u2=eval (solve ('2*u2+p4=0','u2'))
u2 = -1/2*(-2*p10+4*p40)*cos(t)+1/2*(4*p30+2*p20)*sin(t)-3/2*p30*t-p10+3/2*p40
diff(s.p1,t)
diff(s.p2,t)
diff(s.p3,t)
diff(s.p4,t)
eval(diff(H,x1))
eval(diff(H,x2))
eval(diff(H,x3))
eval(diff(H,x4))
s1=solve('-3*sin(t) *(-p10+2*p40)-3*cos(t) *(2*p30+p20)+6*p30 =-(3*(sin(t)*(- p10+2*p40)+cos(t)*(2*p30+p20)-2*p30','(cos(t) *(-p10+2*p40)-sin(t) *(2*p30+p20)) =-(3*cos(t)*(-p10+2*p40)-3*sin(t)*(2*p30+p20)+6*p30* t+4*p10-6*p40 -(4*(-p10+2*p40) *cos(t)-4*(2*p30+p20) *sin(t)+6*p30*t+1/ 4*p10-6*p40)))','0=0','-2*sin(t)*(-p10+2*p40)-2*cos(t)*(2*p30+p20)+3*p30=-(2*sin(t)*(-p10+2*p40)+2*cos(t)*(2*p30+p20)-3*p30)','p10','p20','p30','p40')
s1.p10
s1.p20
s1.p30
s1.p40
p10=1 p20=1 p30=1 p40=1
p10=1 p20=1 p30=1 p40=1
63
Comanda optimală:
u1opt=eval (u1)
u1opt = -1/2*sin(t)-3/2*cos(t)+1
u2opt=eval (u2)
u2opt = -cos(t)+3*sin(t)-3/2*t+1/2
Determinăm traiectoria optimală:
q=t-g;
C=expm(t*A)*x0
D=expm(q*A)*B*uopt
x1_opt=C(1,1)+int(D(1,1),g,0,t)
x1_opt =-3/2*sin(t)-5/2*cos(t)+1-2*cos(t)*t+9*sin(t)*t-3*t^2+t+5/2*cos(t)*sin(t)+3/2*cos(t)^2-6*sin(t)^2
x2_opt=C(2,1)+int(D(2,1),g,0,t)
x2_opt =1-3*cos(t)+7*sin(t)-3*t-1/2*sin(t)^2-15/2*cos(t)*sin(t)+2*cos(t)^2+3*cos(t)*t
x3_opt=C(3,1)+int(D(3,1),g,0,t)
x3_opt = 2+sin(t)*t+9*cos(t)*t-8*t+3/2*cos(t)*t^2-6*cos(t)-9/2*sin(t)*t^2+14*sin(t)+9/4*t^3-3/4*t^2-sin(t)^2-15*cos(t)*sin(t)+4*cos(t)^2
x4_opt=C(4,1)+int(D(4,1),g,0,t)
x4_opt = 3*sin(t)+5*cos(t)-2+3*cos(t)*t-15*sin(t)*t+9/2*t^2-3/2*t-5*cos(t)*sin(t)-3*cos(t)^2+12*sin(t)^2
Traiectoria optimală (în coordonate polare):
r_opt=x1_opt+1
r_opt = -3/2*sin(t)-5/2*cos(t)+2-2*cos(t)*t+9*sin(t)*t-3*t^2+t+5/2*cos(t)*sin(t)+3/2*cos(t)^2-6*sin(t)^2
theta_opt=x3_opt+omega*t
theta_opt = 2+sin(t)*t+9*cos(t)*t-7*t+3/2*cos(t)*t^2-6*cos(t)-9/2*sin(t)*t^2+14*sin(t)+9/4*t^3-3/4*t^2-sin(t)^2-15*cos(t)*sin(t)+4*cos(t)^2
Traiectoria optimală (în coordonate carteziene):
64
x=r_opt.*cos(theta_opt)
x = (-3/2*sin(t)-5/2*cos(t)+2-2*cos(t)*t+9*sin(t)*t-3*t^2+t+5/2*cos(t)*sin(t)+3/2*cos(t)^2-6*sin(t)^2)*cos(2+sin(t)*t+9*cos(t)*t-7*t+3/2*cos(t)*t^2-6*cos(t)-9/2*sin(t)*t^2+14*sin(t)
+9/4*t^3-3/4*t^2-sin(t)^2-15*cos(t)*sin(t)+4*cos(t)^2)
y=r_opt.*sin(theta_opt)
y = (-3/2*sin(t)-5/2*cos(t)+2-2*cos(t)*t+9*sin(t)*t-3*t^2+t+5/2*cos(t)*sin(t)+3/2*cos(t)^2-6*sin(t)^2)*sin(2+sin(t)*t+9*cos(t)*t-7*t+3/2*cos(t)*t^2-6*cos(t)-9/2*sin(t)*t^2+14*sin(t)
+9/4*t^3-3/4*t^2-sin(t)^2-15*cos(t)*sin(t)+4*cos(t)^2)
tf=2;
J=eval(1/2*int(uopt.'*uopt,0,tf))
J = 1.0015
Reprezentăm grafic traiectoria optimală:
clear all
t=0:0.1:3;
r_opt=-3/2.*sin(t)-5/2.*cos(t)+2-2.*cos(t).*t+9.*sin(t).*t-3.*t.^2+t+5/2.*cos(t).*sin(t)+3/2.*cos(t).^2-6*sin(t).^2;
theta_opt=2+sin(t).*t+9.*cos(t).*t-7.*t+3/2.*cos(t).*t.^2-6.*cos(t)-9/2.*sin(t).*t.^2+14.*sin(t)+9/4.*t.^3-3/4.*t.^2-sin(t).^2-15.*cos(t).*sin(t)+4.*cos(t).^2;
figure, polar (theta_opt,r_opt,'--r')
x=r_opt.*cos (theta_opt);
y=r_opt.*sin (theta_opt);
figure,plot3(x,y,t),xlabel('x'),ylabel('y'),zlabel('t'),grid on
65
Anexa 3
Codul programului în Maple pentru rezolvarea aplicaţiei 4.2
> restart;
Încărcăm pachetul LinearAlgebra:
> with(LinearAlgebra):
Definim constantele ce caracterizează racheta
> Lp:=-2;Lda:=9000;phi_max:=0.174;pmax:=5.23;damax:=0.522; := Lp -2
:= Lda 9000 := phi_max 0.174
:= pmax 5.23 := damax 0.522
Definim matricile ce descriu sistemul nostru:
> A:=<<0,0>|<1,Lp>>;
:= A
0 1
0 -2
> B:=<<0,Lda>>;
:= B
0
9000
> Q:=<<1/phimax^2,0>|<0,1/pmax^2>>;
66
:= Q
33.02946228 0
0 0.03655919482
> R:=<<1/damax^2>>;
:= R [ ]3.669940253
Matricea S este necunoscuta din ecuaţia Riccati:
> S:=<<s11,s12>|<s12,s22>>;
:= S
s11 s12
s12 s22
Rezolvăm ecuaţia Riccati
> M:=S.A+Transpose(A).S+Q-S.B.MatrixInverse(R).Transpose(B).S;
M :=
[ , ]33.02946228 0.2207120400108 s122 s11 2 s12 0.2207120400108 s12 s22
s11 2 s12 0.2207120400108 s12 s22[ ,
2 s12 4 s22 0.03655919482 0.2207120400108 s222 ]
> ec1:=M[1,1]=0;
:= ec1 33.02946228 0.2207120400108 s122 0
> ec2:=M[1,2]=0;
:= ec2 s11 2 s12 0.2207120400108 s12 s22 0
> ec3:=M[2,1]=0;
:= ec3 s11 2 s12 0.2207120400108 s12 s22 0
> ec4:=M[2,2]=0;
:= ec4 2 s12 4 s22 0.03655919482 0.2207120400108 s222 0
> s1:=solve(ec1,s12);
:= s1 ,-0.001223313418 0.001223313418
> s12:=s1[2];
:= s12 0.001223313418
> s2:=solve(ec4,s22);
:= s2 ,-0.00004212964825 0.00004194841663
67
> s22:=s2[2];
:= s22 0.00004194841663
> s3:=solve(ec2,s11);
:= s3 1.135053876
> s11:=s3;
:= s11 1.135053876
Soluţia ecuaţiei Riccati este:
> S;
1.135053876 0.001223313418
0.001223313418 0.00004194841663
> x:=<<phi(t),p(t)>>;
:= x
( ) t
( )p tAstfel obţinem matricea K:
> K:=MatrixInverse(R).Transpose(B).S;
:= K [ ]3.000000001 0.1028724512
Comanda optimală este:
> uop:=-MatrixInverse(R).Transpose(B).S.x; := uop [ ] 3.000000001 ( ) t 0.1028724512 ( )p t
Înlocuim comanda optimală în ecuaţia de stare a sistemului şi rezolvăm sitemul liniar astfel obţinut
> with(DEtools):
> M:=A-B.MatrixInverse(R).Transpose(B).S;
:= M
0. 1.
-27000. -927.852060899999970
Condiţii iniţiale
> t0:=0;
:= t0 0
> phi(t0):=0.087;p(t0):=2.616; := ( ) 0 0.087
68
:= ( )p 0 2.616
> with(linalg):
> sol:=exponential(M,t).<<phi(t0),p(t0)>>;
:= sol
0.006030239587e
( ) 897.7778044t0.09303023959e
( ) 30.07425654t
2.797815291e( ) 30.07425654t
5.413815290e( ) 897.7778044t
> phi(t):=sol[1,1];
:= ( ) t 0.006030239587e( ) 897.7778044t
0.09303023959e( ) 30.07425654t
> limit(phi(t),t=infinity);0.
> p(t):=sol[2,1];
:= ( )p t 2.797815291e( ) 30.07425654t
5.413815290e( ) 897.7778044t
> limit(p(t),t=infinity);
0.
Înlocuind şi p(t) obţinute mai sus comanda optimală devine:
> uop;
[ ] 0.5388417304e( ) 897.7778044t
0.0087273981e( ) 30.07425654t
Calculăm funcţionala de cost corespunzătoare comenzii optimale
> F:=(phi(t)/phi_max)^2+(p(t)/pmax)^2+(uop[1,1]/damax)^2;
F 33.02946228( ) 0.006030239587e( ) 897.7778044t
0.09303023959e( ) 30.07425654t 2
:=
0.03655919482( ) 2.797815291e( ) 30.07425654t
5.413815290e( ) 897.7778044t 2
3.669940253( ) 0.5388417304e( ) 897.7778044t
0.0087273981e( ) 30.07425654t 2
> J:=1/2*int(F,t=0..infinity);
:= J 0.004717563808
Considerăm o altă comandă şi comparăm noul cost cu costul obţinut pentru comanda optimală
> u2:=-5*csi(t)-10*p(t);
:= u2 54.10800170e( ) 897.7778044t
27.51300171e( ) 30.07425654t
> F2:=(phi(t)/phi_max)^2+(p(t)/pmax)^2+(u2/damax)^2;
69
F2 33.02946228( ) 0.006030239587e( ) 897.7778044t
0.09303023959e( ) 30.07425654t 2
:=
0.03655919482( ) 2.797815291e( ) 30.07425654t
5.413815290e( ) 897.7778044t 2
3.669940253( ) 54.10800170e( ) 897.7778044t
27.51300171e( ) 30.07425654t 2
> J2:=1/2*int(F2,t=0..infinity);
:= J2 20.20120104
Observăm astfel că noul cost este mai mare decât cel obţinut pentru comanda optimală.
Anexa 4
Codul programului în Matlab pentru rezolvarea aplicaţiei 4.2
clear all
Definim constantele ce caracterizează racheta:
Lp=2;Lda=9000;
phi_max=0.174;pmax=5.23;dmax=0.522;
Definim matricile ce descriu sistemul nostru:
A=[0,1;0,Lp]
B=[0;Lda]
Q=[1/phi_max^2,0;0,1/pmax^2]
R=[1/dmax^2]
A = 0 1
0 2
B = 0
9000
Q = 33.0295 0
0 0.0366
R = 3.6699
Rezolvăm ecuaţia Riccati asociată:
[K,S,e]=lqr(A,B,Q,R)
K = 3.0000 0.1033
70
S = 1.1351 0.0012
0.0012 0.0000
e = -897.7778
-30.0743
phi=sym('phi');
p=sym('p');
x=[phi;p];
Determinăm comanda optimală folosindune de matricea K obţinută mai sus:
uop=[];
uop=-K*x
uop =-6755399441055787/2251799813685248*phi-1861191731142797/18014398509481984*p
Condiţii iniţiale
x0=[0.087;2.616]
x0 = 0.0870
2.6160
Înlocuim comanda optimală în ecuaţia de stare a sistemului şi determinăm x(t):
t=sym('t');
M=A-B*K
M = 1.0e+004 * 0 0.0001
-2.7000 -0.0928
x=eval(expm(t*M)*x0)
x =
418970952517043201/4503599627370496000*exp(-264535958050349/8796093022208*t)-27157784935810049/4503599627370496000*exp(-7896937080653855/8796093022208*t)
761926766614425899/140737488355328000*exp(-7896937080653855/8796093022208*t)-393757497076887851/140737488355328000*exp(-264535958050349/8796093022208*t)
phi=x(1)
phi =
71
418970952517043201/4503599627370496000*exp(-264535958050349/8796093022208*t)-27157784935810049/4503599627370496000*exp(-7896937080653855/8796093022208*t)
p=x(2)
p =
761926766614425899/140737488355328000*exp(-7896937080653855/8796093022208*t)-393757497076887851/140737488355328000*exp(-264535958050349/8796093022208*t)
Înlocuind şi p(t) comanda optimală devine de forma:
uop=-K*x
uop =
101116652087706599686274811782801/10141204801825835211973625643008000*exp(-264535958050349/8796093022208*t)-
5488905505860864342848128732194449/10141204801825835211973625643008000*exp(-7896937080653855/8796093022208*t)
Calculăm funcţionala de cost corespunzătoare comenzii optimale
F=(phi/phi_max)^2+(p/pmax)^2+(uop/dmax)^2;
J=eval(1/2*int(F,t,0,inf))
J = 0.0047
6. BILIOGRAFIE
[1] D. Mihoc, M. Ceaparu, S. St. Iliescu, I. Borangiu, “Teoria si elementele sistemelor de reglare automata”, Editura Didactica si Pedagogica, Bucuresti, 1980
[2] Prof. Dr. V. Prepeliţă, T. Vasilache, M. Doroftei,” Control Theory”
[3] Cursuri Teoria sistemelor si control optimal, anul III, Facultatea de Stiinte Aplicate
[4] Nicolae Racoveanu, “Automatica”, Editura Militară, Buc, 1980
[5] Dr. Robert C. Nelson, “Flight Stability and Automatic Control”, McGraw-Hill, Inc., 1989
[6] Frank L. Lewis, Vassilis L. Syrmos, “Optimal Control”, Second Edition, John Wiley Sons, New York, 1995
72
73
