Transfer orbital optimal folosind propulsie electrică [628667]
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
1
Transfer orbital optimal folosind propulsie electrică
Absolvent: [anonimizat]. Simioana Mădălin Viorel
Conducător: Prof. Dr. Ing. Teodor Viorel Chelaru
2019
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
2
CUPRINS
1 Introducere …………………………………………………………………………. 4
1.1 Structura lucrării …………………………………………………………………………….. 5
1.2 Unelte folosite ………………………………………………………………………………… 5
2 Baze teoretice ………………………………………………………………………. 6
2.1 Introducere …………………………………………………………………………………….. 6
2.2 Coordonatele echinocțiale ………………………………………………………………… 6
2.3 Formularea problemei ……………………………………………………………………… 6
2.4 Problema de control optimal. Metoda directă contra metodei indirecte. … 9
2.5 Aplicarea calculului variațiilor într-o problemă de control optimal …….. 11
2.5.1 Ecuațiile Euler-Lagrange ……………………………………………………………………………………….. 11
2.5.2 Timpul final neconstrâns ……………………………………………………………………………………….. 13
2.5.3 Obținerea problemei cu condiții la limită …………………………………………………………………. 14
2.5.4 Normalizarea ………………………………………………………………………………………………………… 16
2.5.5 Schimbarea variabilelor independente ……………………………………………………………………… 17
2.5.6 Hamiltonianul ………………………………………………………………………………………………………. 17
2.5.7 Legea de control optimal în problema masei minime ………………………………………………… 18
2.5.8 Ecuațiile diferențiale adjuncte ………………………………………………………………………………… 22
3 Metode numerice ……………………………………………………………….. 23
3.1 Introducere …………………………………………………………………………………… 23
3.2 Algoritmul Evoluției Diferențiale ……………………………………………………. 23
3.3 Shooting Method …………………………………………………………………………… 25
3.4 Propagarea numerică ……………………………………………………………………… 27
3.5 Medierea ecuațiilor ……………………………………………………………………….. 29
4 SOFTWARE ……………………………………………………………………… 32
4.1 Prezentare generală ……………………………………………………………………….. 32
4.2 Programul principal ………………………………………………………………………. 32
4.2.1 WRITE PANEL ……………………………………………………………………………………………………. 33
4.2.2 PLOT DATA AND VISUALIZE RESULTS …………………………………………………………… 34
4.2.3 READ INPUT ………………………………………………………………………………………………………. 35
4.2.3.1 READ PANELS ………………………………………………………………………………………………… 35
4.2.3.2 READ INITIAL GUESS …………………………………………………………………………………….. 36
4.2.4 PROCESS INITIAL GUESS ………………………………………………………………………………….. 37
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
3
4.2.5 DEFINE CONSTRAINTS……………………………………………………………………………………… 38
4.2.5.1 CONSTRÂNGERI LINIARE DE EGALITATE ……………………………………………………. 38
4.2.5.2 CONSTRÂNGERI ALE STĂRII INIȚIALE ORIGINALE …………………………………….. 38
4.2.5.3 CONSTRÂNGERI LEGATE DE DURATA ARCURILOR …………………………………… 38
4.2.5.4 CONSTRÂNGERI NELINIARE DE EGALITATE ………………………………………………. 39
4.2.6 EVALUATE CONSTRAINTS AND COST FUNCTION …………………………………………. 39
4.2.6.1 EVALUAREA FUNCȚIEI COST ……………………………………………………………………….. 39
4.2.6.2 EVALUAREA CONSTRÂNGERILOR MODULULUI VITEZELOR ÎNVECINATE ȘI
GRADIENTUL ACESTORA …………………………………………………………………………………………….. 39
4.2.6.3 EVALUAREA CONSTRÂNGERILOR DE DEFECT DEFINITE DE SFÂRȘITUL
UNUI ARC ………………………………………………………………………………………………………………………. 40
4.2.6.4 EVALUAREA CONSTRÂNGERILOR DEFINITE DE STAREA FINALĂ ……………. 40
4.2.7 RUNGE-KUTTA ………………………………………………………………………………………………….. 41
4.2.8 CHECK-STOP CRITERIA ……………………………………………………………………………………. 41
4.2.9 PERFORM AUXILIARY COMPUTATIONS …………………………………………………………. 41
5 SOFTWARE : Interfața grafică ………………………………………….. 42
5.1 Introducere …………………………………………………………………………………… 42
5.2 Meniul Principal ……………………………………………………………………………. 43
5.2.1 Date de intrare ………………………………………………………………………………………………………. 43
5.2.2 Controlul procesului ……………………………………………………………………………………………… 44
5.3 Consola de Execuție ………………………………………………………………………. 44
5.4 Afișarea rezultatelor grafice ……………………………………………………………. 46
6 Rezultate ……………………………………………………………………………. 48
6.1 Introducere …………………………………………………………………………………… 48
6.2 Transfer GTO-GEO ………………………………………………………………………. 48
6.3 Rezultatele transferului ………………………………………………………………….. 48
7 Concluzii ……………………………………………………………………………. 56
7.1 Obiective atinse …………………………………………………………………………….. 56
7.2 Implementări viitoare …………………………………………………………………….. 56
8 Bibliografie ………………………………………………………………………… 57
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
4
1 Introducere
Mișcarea corpurilor cerești sub acțiunea fortelor gravitaționale este una din cele mai discutate probleme
a fizicii moderne. Imediat după descoperirea lui sir Isaac Newton, mulți matematicieni înzestrați și-au dedicat
timp și au dezvoltat domeniul mecanicii orbitale, încercând să rezolve, de exemplu, problema n-body. Acest
domeniu al științei a devenit mult mai important în ultimul deceniu, datorită dezvoltării rachetelor și a
tehnologiilor satelitare.
O operație critică în domeniul ingineriei aerospațiale este plasarea unui vehicul spațial într-o orbită cât
mai apropiată de cea dorită inițial. Un satelit plasat în afara orbitei dorite va fi foarte neproductiv, dacă va
produce în vreun fel, totuși, un transfer orbital ineficient poate reduce drastic viața satelitului. Datorită acestor
probleme, o activitate importantă o reprezintă calculul unei traiectorii de transfer optimale pentru un corp
mic ce orbitează o planetă.
În prezent, mare parte din transferuri sunt aplicabile doar pentru rachete cu propulsie chimică, acestea
furnizând o creștere masivă în viteză, într-o perioadă foarte scurtă de timp, fiind modelate foarte precis ca
fiind impulsuri instantanee. Un astfel de model nu poate fi aplicat unei rachete cu propulsie electrică, din
acest motiv, găsirea unui transfer optimal folosind un motor cu tracțiune redusă este una din cele mai
complexe probleme de rezolvat pentru un inginer din domeniul aerospațial.
Încă de la începutul erei satelitare, s-a remarcat o orbită deosebit de atractivă în jurul Pămantului: inelul
geostaționar. Perioada acestei orbite este egală cu cea a Pămantului, astfel păstrând satelitul în același punct
deasupra Ecuatorului. Această orbită este cea mai favorabilă pentru o multitudine de aplicații, cum ar fii
sistemele de comunicații, din moment ce sunt necesari doar 3-4 sateliți pentru a acoperi aproape întreaga
suprafață a Pămantului în orice moment de timp (regiunile cu laitudine înaltă sunt mai slab acoperite de
această orbită, din fericire populația în aceste regiuni este foarte scăzută), totuși, dacă sistemul de comunicații
ar orbita într-o orbită joasă, ar fi nevoie de aproximativ 70 de sateliți pentru o acoperire totală a Pămantului.
Totuși există un mare dezavantaj al orbitei geostaționare: altitudinea de minim 35786 km, deci un satelit
cu o rachetă cu combustibil chimic are nevoie de o cantitate de combustibil de același ordin de magnitudine
ca și masa întregului ansamblu fără combustibil, astfel crescând enorm costul total al misiunii. Acest fapt
stârnește interes asupra propulsiei electrice, care, comparativ cu propulsia chimică, scade tracțiunea și crește
timpul de transfer în schimbul combustibilului. Cu alte cuvinte, un motor electric consumă cu mult mai puțin
combustibil, dar orbitează și mult mai mult timp față de un transfer chimic. Acest timp de transfer poate părea
inacceptabil la prima vedere, dar reducerea masei la lansare împreună cu viața aproape nelimitată a unui
satelit geostaționar mai mult decât compensează, acesta fiind motivul principal al creșterii propulsiei electrice
în acest domeniu, și aplicarea acesteia în misiuni spațiale.
Cu toate acestea, revenim la problema anterioară și anume: modelul de orbitare pentru propulsie chimică
nu poate fi aplicat, din moment ce acele delta-V-uri instantanee aplicate satelitului sunt o ipoteză departe de
a fi adevarată în cazul unui motor electric cu tracțiune redusă. Din acest motiv, este necesar să dezvoltăm o
noua teorie pentru a folosi în mod eficient propulsia electrică.
Contrastând cu schimbări bruște de viteză produse de un motor chimic, un sistem propulsat electric poate
fi modelat ca fiind o mică perturbație în timp, nu foarte diferită de frânarea atmosferică sau oblativitatea
Pământului. Acest lucru făcând problema de optimizare să fie cu atât mai dificilă, găsirea unei legi ale
tracțiunii care va furniza ce mai eficient transfer, în termeni de consum al masei sau al timpului de transfer.
În această lucrare, procesul găsirii unui transfer optimal între două orbite, folosind propulsie electrică
este descris în amănunt, incluzând modelarea matematică a problemei, metodele considerate pentru a o
rezolva, algoritmii numerici implementați într-o interfată grafică și câteva rezultate obținute.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
5
1.1 Structura lucrării
Această lucrare este imparțită în 7 capitole, fiecare conținând descrierea unei părți de lucru.
Capitolul 2 – Reprezintă teoria aplicată pentru rezolvarea problemei de optimizare, furnizând date despre
pricipalele abordări folosite în ziua de azi cu avantajele și dezavantajele acestora. De asemenea include și
descrierea metodei alese și toate funcțiile analitice dezvoltate.
Capitolul 3 – Descrie procesele numerice folosite pentru a rezolva problema.
Capitolul 4 – Descrie programul numeric și toate modulele sale, care implementează toate procesele
descrise în capitolele anterioare.
Capitolul 5 – Descrie interfața cu programul numeric din capitolul anterior, incluzând metoda de folosire
a acestuia.
Capitolul 6 – Reprezintă rezultatele obținute cu acest program pentru diferite teste.
Capitolul 7 – Concluzii
1.2 Unelte folosite
Limbajul de programare folosit pentru implmentarea solverului este Fortran 90. Acesta este unul dintre
cele mai folosite limbaje în domeniul aerospațial datorită vitezei de procesare, a usurinței de folosire și a
librăriilor freesource deja existente.
Solverul va fi integrat într-o interfață grafică aceasta fiind dezvoltată în Visual Basic și C#.
Pentru plotare este folosit GNUPLOT în legatură cu interfața grafică pentu a oferi utilizatorului o
experiență cât mai placută.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
6
2 Baze teoretice
2.1 Introducere
În acest capitol vor fi prezentate în amănunt conceptele și ecuațiile folosite pentru a rezolva problema de
transfer optimal, de la setul de coordonate ales până la descrierea metodei ce efectuează optimizarea.
2.2 Coordonatele echinocțiale
Există diferite seturi de coordonate ce pot defini starea unui corp ceresc. Unul din cele mai folosite seturi
și de asemenea unul dintre cele mai cunoscute sunt elementele clasice ,,,,,ieä , care permit aflarea
orbitei descrise de corp cu usurință, având totuși câteva dezavantaje importante, cum ar fi nedeterminarea
ascensiunii drepte a nodului ascendent, Ω, și a argumentului la perigeu, ω, când orbita este equatorială, de
exemplu i = 0. Pentru a preveni această nedeterminare, sunt folosite coordonate echinocțiale, care pot fi
definite în funcție de elementele clasice astfel:
cos( )
sin( )
tan sin2
tan cos2a a
h e
k e
ip
iq
L v
(2.1)
Unde a ramane neschimbat ca semi-axa mare, h și k reprezinta componente ale vectorului eccentricității
e, p și q sunt componente ale vectorului nodului ascendent iar L reprezintă longitudinea medie sau
argumentul longitudinii reale. De fapt, aceasta reprezintă definiția elementelor echinocțiale directe. Setul
direct are o singularitate când i ( 2 2limi hp q ), pe când pentru elementele retrograde aceasta
singularitate apare cand i = 0. Deoarece i este foarte rar folosit, sau deloc pentru misiunile geostaționare,
elementele alese au fost echinocțiale directe. [1]
2.3 Formularea problemei
Pe scurt, problema considerată luată în considerare este: obținerea vectorului optim de tracțiune în timp,
( )u t, care va permite satelitului sa călătorească de la o poziție inițială 0( )x t , la o poziție finală diferită ( )fx t
folosind cât mai puțin combustibil, adică maximizarea fm. Tracțiunea disponibilă a satelitului este foarte
mică în fiecare moment de timp al transferului. Tracțiunea este considerată ca fiind mică atunci când raportată
cu atracția gravitațională a corpului orbitat valoarea atinge 310. Datorită acestui fapt, traiectoria nu poate fi
calculată ca o schimbare instantanee a orbitei făcută de către impulsuri instantanee a motorului, calculul fiind
mult mai precis când tracțiunea este considerată ca o mică perturbație ce acțiunează asupra navetei spațiale,
aceasta “plutind” ușor spre orbita obiectiv.
Pentru a rezolva această problemă, este necesar să știm cum variază vectorul de stare x, construit din
elementele echinocțiale descrise anterior (2.1), sub acțiunea tracțiunii. Desigur, cu cât considerăm mai multe
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
7
perturbații în afară de tracțiune, cu atât problema devine mai dificil de rezolvat. Pentru cel mai simplu caz,
considerând doar perturbațiile datorate tracțiunii, ecuațiile diferențiale pot fi tratate analitic pentru majoritatea
pașilor, totuși pentru a oferi o perspectivă mai bună a întregului proces de optimizare, și pentru a obține
rezultate mai interesante este necesară implementarea unor modele mai complexe.
Având în vedere că forma Gauss a variației elementelor echinocțiale în timp, când subiectul de interes
este tracțiunea, este următoarea [2]:
( , )dx uF x Ldt m (2.2)
0 1( , ) ( , ).dL ug x L g x Ldt m (2.3)
unde,
m – masa satelitului;
x – este compus din toate elementele echinocțiale mai puțin L, [ , , , , ]x a h k p q;
iar, 0( , ), ( , )F x L g x L și 1( , )g x L sunt descrise de următoarele ecuații [2]:
2
2 2
2 2
2 2
2 220 0
2 ( cos ) 2 cos sin ( ) 2 sin( sin cos )
2 ( sin ) 2 sin cos ( ) 2 cos( , ) ( sin cos )
cos0 0 (1 )2
sin0 0 (1 )2aBD
A
D h L hk L L h k h Lk p L q LB B
a AD k L hk L L h k h LF x L h p L q LD B B
Lp q
Lp q
(2.4)
2
0 3 3( , )Dg x La A (2.5)
1( , ) 0 0 ( sin cos )a Ag x L p L q LD
(2.6)
Unde, A, B și D sunt definiți astfel:
2 2
2 21
1 2 cos 2 sin
1 cos sinA h k
B h L k L h k
D h L k L
(2.7)
iar este constanta gravitațională, pentru pământ având valoarea 398600.47 3 2/km s. [2]
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
8
Fiecare element al matricii F descrie influența unei componente a tracțiunii asupra unui element
echinocțial. Scalarul 0g reprezintă schimbarea longitudinii reale în timp datorită miscării satelitului în jurul
propriei axe, iar 1g arată influența tracțiunii asupra variației longitudinii reale.
Longitudinea reala, L, este separată de celelalte elemente echinocțiale din x , pentru că variază mult
mai rapid decât celelalte element, variind de la 0 la 2 la fiecare rotație.
Vectorul tracțiunii u este construit din 3 componente, și este referențiat la un cadru propriu, în care,
prima axă este îndreptată către direcția tangențială la orbită, a doua formează un sistem ortogonal și este
îndreptata catre concavitatea orbitei iar a treia este perpendiculară pe planul orbitei, obținută prin produs
vectorial dintre vectorul poziției și vectorul vitezei. În afară de asta, se consideră că nu există restricții asupra
direcției vectorului u, dar modulul acestuia va avea un maxim prevăzut de producătorul motorului.
Fig. 1 – Orientarea sistemului de referintă local
Unghiurile și care apar în figura de mai sus sunt unghiurile de direcție, care descriu direcția
vectorului tracțiunii cu privire la sistemul de referință descris. Folosind aceste unghiuri, componentele
tracțiunii pot fi definite astfel:
| |cos cos
| |cos sin
| |sint
n
zu u
u u
u u
(2.8)
O a șaptea variabilă de stare este necesară pentru a putea descrie complet mișcarea satelitului: masa
acestuia. Pentru asta, o a șaptea ecuație privind variația masei este adăugată în sistem:
| |
P SPdm u
dt g I
(2.9)
unde Pg este accelerația gravitațională a planetei Pământ, având valoarea 29.81 /m s și SPI este impulsul
specific al motorului electric, în secunde. SPI este de asemenea dat de către producătorul motorului.
Acum că am construit toate ecuațiile necesare pentru a descrie problema, tot ce avem de făcut este
setarea unor condiții la limită pentru închiderea ei. Starea inițială este complet cunoscută, atât orbita cât și
masa, de asemenea și orbita dorită. Longitudinea reală finală este necunoscută deoarece variază prea rapid
pentru a putea oferi un rezultat satisfăcător pe parcursul propagării. Deoarece depinde foarte mult de
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
9
traiectorie, masa finală este de asemenea necunoscută. În final, timpul final de transfer este tot o valoare
liberă datorită faptului că la sfârșitul propagării, este valoarea ce va fi optimizată. Astfel problema ce va fi
rezolvată va arăta astfel:
0 1
0 0
0 0
0 0min
( , )
( , ) ( , )
| |
( ) ( )
( ) ( )
( ) ( )P SP
f f
f
f
fimize mf
dx uF x Ldt m
dL ug x L g x Ldt m
dm u
dt g I
x t x x t x
L t L L t liber
m t m m t liber
t fix
(2.10)
Acesta este un exemplu de problemă de control optimal și se formulează astfel: se dă un sistem
dinamic de x variabile de stare, afectate de u variable de control și supusă la k constrangeri, găsiți starea și
valorile variabilelor de control pentru a minimiza o funcție cost ( , , )J x u t respectând dinamica și
constrângerile. În cazul nostru, fJ m și nu se vor aplica constrângeri. Totuși, pentru a modela o problemă
de control mai realistică pentru o misiune spațială, constrângerile devin foarte relevante pentru problemă,
permițând să limiteze variația tracțiunii, de exemplu pentru a evita ca instrumente sensibile sa fie îndreptate
direct spre soare, sau pentru a asigura vizibilitate asupra satelitului din diferite stații terestre. Cu toate acestea,
astfel de constrângeri ar crește semnificativ complexitatea problemei, nu doar în timpul aflării soluției
numerice, dar și în implementarea funcțiilor matematice.
Problema principală acum este găsirea uneltelor matematice capabile să rezolve problema de control
optimal. În următorul subcapitol vom detalia principalele abordări asupra problemei.
2.4 Problema de control optimal. Metoda directă contra metodei indirecte.
Din punct de vedere matematic, problemele de control optimal sunt formulate astfel [3]: determinați
starea, ( )x t, vector real de dimensiune n, controlul, ( )u t, vector real de dimensiune m, unde valoarea inițială
0t și valoarea finală ft a unei variabile reale independente 0,ft t t care optimizează indexul de
performanță:
00 0( ), , ( ), ( ), ( ),ft
f f
tJ x t t x t t x t u t t dt (2.11)
supus dinamicii (construită din ecuații diferențiale) constrânse,
[ ( ), ( ), ]x f x t u t t (2.12)
constrângerile rutei sau constrângerile traiectoriei,
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
10
min max( ), ( ), C C x t u t t C (2.13)
și de condițiile la limită,
min max( ), ( ),x t u t t (2.14)
Indexul de performanță, J, este format din doi termeni, primul, , este funcția care va fi minimizată.
De exemplu fm pentru o problema de masă minimă de combustibil. este funcția care va fi minimizată
pe tot parcursul traiectoriei. Poate fi folosită de exemplu când avem ca obiectiv minimizarea energiei în
fiecare punct 2/2 /v r .
Constrângerile traiectoriei C, reprezintă condiții plasate într-un punct intermediar al orbitei sau într-o
serie de puncte, cum ar fi setarea unui minim asupra semi axei mari sau constrângerea satelitului să treaca
printr-un anumit punct pentru vizibilitate.
În cele din urmă, condițiile la limită reprezintă constrângerile ce vor fi îndeplinite la extremitățile
traiectoriei. De exemplu, elementele orbitei inițiale și finale sunt stabilite aici.
Există doua abordări generale în încercările de a rezolva o problemă de control optimal, și pot fi descrise
pe scurt astfel: fie aplicarea calcului variațiilor pentru a găsi condițiile de ordinul întâi urmând să constrângă
starea și variabilele de control, obținând astfel o problemă cu două puncte la limită, sau a încerca iterativ
găsirea unui set de variabile care satisfac ecuațiile știute, transformând problema într-o problema de
programare neliniară (NLP). A doua abordare este mai directă, pentru că se dedică găsirii variabilelor
necesare rezolvării problemei, din cauza aceasta se numește și metoda directă. Prima abordare mai degrabă
rezolvă problema crescând ordinul sistemului, și din cauza asta se numește metoda indirectă.
Ambele metode prezintă avantaje dar și dezavantaje. Metodele indirecte sunt mai ușor de înțeles și
furnizează aproximări excelente atunci când converg. Din păcate, această metodă este foarte sensibilă la
condițiile inițiale ale problemei, necesitând valori inițiale foarte precise pentru mai multe variabile. Mai mult
de atât, unele dintre variabilele de intrare nu au niciun înțeles fizic prin urmare lipsind orice experiență
anterioară care ar putea ajuta la găsirea stării inițiale a unei noi probleme. Cea mai bună cale de a rezolva
această problemă este de a străbate toată gama de valori inițiale pentru a găsi una destul de buna pentru a fi
rezolvată, astfel crescând foarte mult timpul de calcul. Altă problemă cu metodele indirecte ar fi că orice
schimbare în definirea problemei, de exemplu, adăugarea unei noi forțe pertubatoare, schimbă toate condițiile
producând o nouă problemă.
Pe de altă parte, metodele directe sunt robuste și o condiție inițială, ca o linie dreaptă între punctul de
start și punctul final, este destul de bun pentru a asigura o soluție, pentru că aceste metode sunt făcute sa
conveargă. Principalul dezavantaj al metodelor directe este timpul de calcul foarte îndelungat față de
metodele indirecte, totuși o metoda ingenioasă și eficientă de manipulare a matricilor rare pot rezolva această
problemă.
Cu toate acestea, problema principală a metodelor indirecte, găsirea unei stări inițiale precise, poate
fi minimizată daca un algoritm care caută soluții inițiale într-o regiune fezabilă este implementat, permițând
obținerea soluției controlului optimal într-un timp mult mai scurt față de metoda directă. Mai mult de atât,
modele simple de perturbații pot fi transformate analitic relativ ușor într-o problema cu două puncte de
extremă urmând a fi implementate, acestea pot oferi niște rezultate foarte fezabile, de asemenea și o
aproximare destul de realistă.
Datorită motivelor enunțate anterior, metoda indirectă este folosită pentru a rezolva problema de
control optimal. În subcapitolul următor, pentru o mai bună întelegere, va fi prezentată o descriere mai
detaliată a teoriei calculul variațiilor.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
11
2.5 Aplicarea calculului variațiilor într-o problemă de control optimal
2.5.1 Ecuațiile Euler-Lagrange
Pentru a rezolva problema de control optimal (2.1.1) – (2.1.4) folosind calculul variațiilor, este în primul
rând necesar sa asociem sistemul dinamic indexului de performanță J, obținând :
00 0( ), , ( ), ( ), ( ), ( ) ( ), ( ),ft
f f
tJ x t t x t t x t u t t t f x t u t t x dt (2.15)
Este important să observăm că pentru a executa acest pas n, noi variabile
au fost adaugate în sistem. Acum
este folositor să definim o funcție scalară, hamiltonianul H, astfel :
( ), ( ), ( ), ( ), ( ), ( ) ( ), ( ),H x t u t t t x t u t t t f x t u t t (2.16)
Acum integrând ecuația (2.15) prin părți obținem :
00 0 0 0( ), , ( ), ( ) ( ) ( ) ( ) ( ), ( ), ( ), ( ) ( )ft
f f f f
tJ x t t x t t t x t t x t H x t u t t t t x t dt
(2.17)
O diferențiere poate fi făcută acum a indexului J în funcție de ( )u t pentru timpi stabiliți ft și 0t :
00 0( ) ( ) .f
ft
t t tH HJ x t x t x u dtx x x (2.18)
De remarcat este faptul că o derivată a unei funcții scalare în funcție de un vector, cum ar fi H
x
, este un
vector ce conține derivatele parțiale a funcției, în funcție de fiecare componentă a vectorului.
Pentru a evita calculul variațiilor lui x produse de un vector u dat, funcțiile de multiplicare
sunt alese
astfel încât coeficientii lui x să dispară. Astfel, derivata lui
în timp devine :
.H f
x x x (2.19)
unde /f x este o matrice de dimensiune n n , ce conține dericate parțiale a fiecărei ecuații diferențiale
în funcție de toate variabilele de stare x .
Avem acum un set de 2n ecuații diferențiale, (2.12) și (2.19), ce vor descrie n variabile de stare x și
n variabile de co-stare,
. Este necesar să obținem viariația controlului u în timp, și condițiile la limită
pentru variabilele de stare și co-stare pentru a închide problema.
În general presupunem un număr i < n de variabile de stare definite pe 0t t, astfel încât variațiile lor
ix la momentul inițial sa fie zero, anulând componentele celul de-al doilea termen din partea dreaptă a
ecuației (2.18). Urmând obiectivul de a anul cel de-al doilea termen, și lăsând libere componente n-i rămase
în vectorul inițial 0( )x t produce o extra condiție pentru co-stările n-i:
0( ) 0t
(2.20)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
12
În ceea ce privește condițiile la limită pentru timpul final, este considerată de asemenea o mărime j<n
de variabile de stare definite, făcând variația acestora zero și anulând parțial primul termen al ecuației (2.18).
Pentru a elimina complet acest termen, astfel micșorând J, variabilele de co-stare n-j ce corespund
componentelor nedefinite n-j a lui ( )fx t trebuie să îndeplinească condiția:
0
ft tx
(2.21)
Pe scurt, setând o componentă liberă a lui x , o condiție la limită este impusă unei componente
asociate a lui
.
După cum este bine știut, la o extremă, funcția J nu suferă nici o variație când sunt făcute schimbări
foarte mici asupra variabilelor de control u asta înseamnă ca J nu depinde de u când acesta este la
minim. Prinvind ecuația (2.18), putem extrage o condiție necesară despre upentru minimizarea lui J:
00, ,fHt t tu
(2.22)
În calculul variațiilor, ecuațiile (2.19), (2.21) și (2.22) sunt numite ecuațiile Euler-Lagrange.
În acest moment, pentru a obține vectorul de control optimal u, ecuațiile diferențiale 2n definite în
(2.12) și (2.19) trebuie rezolvate obținând u prin ecuația (2.22) și folosind ca și condițiile la limită valori
știute ale lui x și ecuațiile (2.20) și (2.21) ca și corespondenți. Acestea definesc o problemă cu două condiții
la limită.
Dacă f
și nu sunt funcții explicite de timp, poate fi dovedit că prima integrală a problemei valorii
la limită există:
H H HH x u ft x u
H H HH u ft u x
H HHt u
(2.23)
Dacă nici una dintre funcțiile f
sau nu sunt funcții explicite de timp, iar u satisface condiția (2.22), atunci
0H iar prima integrală a problemei există. Dar pentru a minimiza J, nu este necesar doar să obținem
0J dar și derivata de ordinul doi trebuie să fie pozitivă pentru toate valorile lui u , în timp ce 0x f
. Atfel ecuația ce trebuie îndeplinită ajunge la forma:
02 2
2 2
2
2 2 2
21 102 2f
ft
t t tH H
x x x uJ x x x u dtu x H H
u x u
(2.24)
Ajungându-se la relația următoare între u și x :
d f fx x udt x u
(2.25)
Ce ne va permite sa obținem, printr-o metodă destul de dificilă x în funcție de u.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
13
2.5.2 Timpul final neconstrâns
Până acum, un factor critic al problemei de calcul variațional nu a fost luat în calcul: lăsarea liberă a variabilei
independente, timpul, pentru punctul final. Acest pas este fundamental pentru a putea minimiza timpul de
transfer din problemă.
Derivata funcției de performanță, J, considerând schimbările în timpul final, devine:
00 0f
ft
f f f
t t tf fJ dt dx t x t t x t x u dtt x x x u u
(2.26)
Știind că variațiile lui x
sunt făcute în timp ce timpul este fix, derivata sa poate fi definită astfel:
f f f fx t x t x t dt
(2.27)
În plus, vectorul cu variabile de co-stare
¨, poate fi ales astfel încât coeficienții lui x și fx t
dispar, urmărind teoria detaliată în secțiunea anterioară. Astfel derivata indicelui de performanță devine:
0*f
ft
f
t t tfJ f dt udtt u u (2.28)
unde doar componentele lui care au eliberat variabilele de stare corespondente sunt diferite de 0.
De asemenea, știind că 0( ) 0x t, acesta dispare.
Considerând ca sunt un număr de variabile de stare, ce sunt definite la timpul final, definite de către
indicele i, astfel i fx t;0,…,i q sunt definiți, iar j fx t;1,…,j q n sunt liberi. Atunci, este posibil
să aflăm influența u asupra i fx t utilizând conceptul funcțiilor adjunct, obținând:
0*f
ft
i
i f i f t t
tfx t f dt t u dtu (2.29)
unde ( ) 1,
0,i
j fi jti j
Ecuația 2.29 poate fi obținută din 2.28 egalând 0 și înlocuind icu x.
Următorul pas este obținerea lui u și găsirea unei valori pentru fdt care va minimiza J. Pentru acest
lucru ecuația 2.29 este înmulțită cu o constantă, , ce va fi determinată ulterior, și adăugarea acesteia la
ecuația 2.28, rezultând:
0( ) ( ) ( )*f
ft
j j i
i i f i i f i
t t tfdJ dx t f f dt udtt u u (2.30)
Alegând:
( )
1 *
fj
f i i
t tdt k f ft
(2.31)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
14
și,
( ) ( )
2j i
ifu ku u (2.32)
unde k1 si k2 sunt constante pozitive, și introducând aceste expreii in 2.30 rezultă:
022
( ) ( ) ( )
1 2 * 0ft
j j i
i i f i i i
tfdJ dx t k f f k dtt u u
(2.33)
Pentru a determina valorile tuturor constantelor icare satisfac ecuația 2.28, ecuațiile 2.31, 2.32 sunt
introduse in 2.28, și după câteva calcule rezultă:
0 01
( ) ( ) ( ) ( ) ( ) 1 1
2 2f f
fft t
j j i j j
i j i it t t tt tk k f f f fdt f f dt f fu u k u u u k
(2.34)
examinând ecuația 2.33, putem trage concluzia că indexul de performanță, J, nu poate fi scăzut doar în cazul
în care:
( )* 0
fj
i i
t tf ft
(2.35)
( ) ( )
00, ,j i
i fft t tu u
(2.36)
Astfel, impunând ecuațiile 2.35 si 2.36 ne permit să obținem o soluție staționară ce va satisface și
constrângerile finale.
În concluzie, problema de timp minim de transfer are aceleași conditii necesare ca și problema cu timp
final fixat, exceptând ecuația 2.35 necesară pentru determinarea unei valori optime pentru ft.
Știind că aceste condiții necesare rămân neschimbate după eliberarea timpului final de transfer, altă cale poate
fi folosită pentru obținerea timpului minim, evitând folosirea ecuației 2.35. Printr-o schimbare isteață de
variabile, timpul poate deveni o variabilă de stare. Aceasta este metoda preferată, și va fi descrisă în secțiunea
următoare.
2.5.3 Obținerea problemei cu condiții la limită
Din ecuațiile obținute în secțiunea anterioară, problema de optimizare poate fi transformată într-o problemă
cu condiții la limită, care poate fi rezolvată numeric relativ ușor. În această secțiune, acest proces este aplicat
problemei de interes:
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
15
0 1
0 0
0 0
0 0min
( , )
( , ) ( , )
| |
( ) ( )
( ) ( )
( ) ( )P SP
f f
f
f
fimize mf
dx uF x Ldt m
dL ug x L g x Ldt m
dm u
dt g I
x t x x t x
L t L L t liber
m t m m t liber
t fix
(2.37)
Pe scurt, operațiile necesare sunt:
1. Adimesionalizarea variabilelor de stare.
2. Efectuarea unei schimbare de variabile, astfel încât timpul să devină o variabilă de stare.
3. Producerea funcției Hamiltoniene H folosind ecuația 2.16
4. Derivarea parțială a lui H, în funcție de vectorului de control, u, si egalarea rezultatului cu 0,
urmărind ecuația 2.22, și obținând legea de control optim în funcție de variabilele de stare și
adjuncte xsi.
5. Rescrierea H introducând rezultatul anterior, și scriindu-l dor in funcție de xsi.
6. Derivarea parțială a lui H în funcție de x și obținerea ecuațiilor diferențiale pentru
ca în ecuația
2.19.
7. Medierea ecuațiilor diferențiale.
8. Rezolvarea sistemului pentru xsi, aplicând condițiile la limită.
Este comodă adimensionalizarea sistemului înainte de a face orice operație. Mai mult, dacă o schimbare de
variabilă este folosită pentru longitudinea adevarată, L, ca variabilă independentă, timpul poate fi considerat
o variabilă de stare urmând a fi optimizat pentru a obține timpul final cu ușurință.
Mai multe manipulări a longitudinii adevărate pot fi făcute, definind-o în termeni de variabilă normalizată,
având valori cuprinse între 0 și 1, și o variabilă independentă de timp, pentru a folosi variabila normalizată
ca o variabilă independentă, astfel permițând propagarea de la 0 la 1, în loc de un timp inițial și final. Această
îmbunătățire, face posibilă rezolvarea problemei valorilor inițiale. După cum vom vedea în continuare este
crucială găsirea unor valori inițiale destul de bune pentru rezolvarea problemei cu condiții la limită.
Un alt pas convenabil de făcut este medierea ecuațiilor diferențiale. Medierea înseamnă integrarea în funcție
de toate variabilele ce se schimbă rapid pe parcursul unei revoluții, minimizând erorile care altfel vor înrăutăți
soluția. Datorită faptului că singura variabilă ce se schimbă rapid în sistemul considerat este L, integrarea
este necesară doar în funcție de aceasta variabilă. Acest pas trebuie parcurs chiar înainte de rezolvarea
problemei cu condiții la limită.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
16
2.5.4 Normalizarea
Primul pas este adimensionalizarea tuturor variabilelor. Știind că sunt folosite coordonate echinocțiale, doar
semiaxa mare are dimensiuni, pe lângă masă, tracțiune și timp. Următoarele schimbări sunt făcute:
1
1 3
0 max f fx m ux m u t ta m u a (2.38)
Variabilele adimensionalizate au fost notate cu o bară deasupra, iar modulul tracțiunii a fost scris ca u,
Valoarea semiaxei mari faeste folosită ca lungime de referință. Masa inițială 0madimensionalizeaza masa,
iar tracțiunea maximă este folosită pentru adimensionalizarea tracțiunii. În final referința timpului folosită
este parametrul gravitațional .
Cu această alegere a parametrilor dimensionali, semiaxa mare, masa, și tracțiunea sunt normalizate. Semiaxa
mare va varia de la o fracțiune a lui 1 la 1, atâta timp cât orbit inițială este mai joasă decât orbita obiectiv.
Masa va fi 1 la timpul inițial și va scădea pe măsura ce propagarea avansează, iar modului tracțiunii nu poate
avea valori mai mari de 1.
Aplicând aceste schimbări de variabilă la ecuațiile diferențiale a problemei noastre, acestea devin:
2
max max
3
0 0
3 2
max max
0 1 0 13
0 0
max max
3 3
0 01
1f f
f f
f f f
f
f P sp f P spa a uu u dx uF Fdt a a m m m m
a a a u u u dL ug g g gdt a m m m m
u u u dmudt a m g I a m g I
(2.39)
Observăm că un parametru adimensional important apare în ecuațiile anterioare:
2
max
0fa u
m (2.40)
Acest parametru, , reprezintă raportul dintre accelerația produsă de tracțiune și accelerația gravitațională.
Prin urmare, aceasta are o valoare foarte mică pentru un transfer GEO a unui satelit de 2000 kg cu o tracțiune
de 0.2 N, 410.
Este folositor să scalăm timpul cu acest nou parametru, astfel creând o nouă variabila independentă
reprezentând timpul care variază mult mai puțin și mai incet, t . Introducând aceste noi variabile in 2.39
acestea devin:
0 0
1
0 max
2 3
max 01
f f P sp f P spdx u uF Fd m m
g gdL ugd m
m u dmu ud a u a m g I a g I
(2.41)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
17
Chiar dacă variabilele sunt adimesionale, barele nu au fost folosite pentru a păstra ecuațiile curate. Din acest
moment acest criteriu va fi folosit. Al doilea termen al ecuației longitudinii este mult mai mic decat primul
și a fost eliminat pentru simplificarea operațiilor, permițând rezultate analitice pentru ecuațiile diferențiale
adjuncte.
2.5.5 Schimbarea variabilelor independente
Este timpul să efectuăm schimbarea variabilelor independente, de la la L. Inversând a doua ecuație din 2.41
și înmulțind-o cu primele 2 ecuații obținem:
0
0 0
0 0d
dL g
dx d u dx uF Fd dL m g dL g m
dm d dmu ud dL g dL g
(2.42)
Schimbarea finală a variabilelor este să introducem o noua variabilă independentă 0,1, care este 0 la
timpul inițial și 1 la timpul final. Dar, datorită faptului că numărul de revoluții este necunoscut, o nouă,
necunoscută, variabilă de stare este adăugată în sistem: 1L, o variabilă independenta de timp care ne da
numarul de revoluțtii cand este înmulțit cu . Cu alte cuvinte schimbarea de variabile făcuta este:
1LL (2.43)
din care va rezulta următoarea schimbare în ecuațiile diferențiale:
1
0
1
0
1
0
10Ldx uFd g m
Ldmud g
Ld
d g
dL
d
(2.44)
2.5.6 Hamiltonianul
Ecuațiile anterioare oferă un sistem solid pentru construirea Hamiltonianului. Pentru asta vom introduce
ecuațiile ecuațiile 2.44 în ecuația hamiltonianului, considerând că în problema noastră indexul de performanță
fJ t, și 0:
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
18
11
1 1 1
0 0 0
1
0x m L
x m
x mdL dx dm dHd d d d
L L L uF ug m g g
L uF ug m
(2.45)
Unde vectorul adjunct ste folosit pentru prima dată. În acest caz este un vector cu 8 randuri, conținând câte o
variabilă pentru fiecare variabilă de stare
1 2 3 4 5 1, , , , , , ,x x x x x m L
. Pentru o nomenclatură mai
compacta folosim
1 2 3 4 5, , , ,x x x x x x
.
Derivând ecuația 2.45 în funcție de fiecare variabilă de stare vom avea derivatele în funcție de a fiecărei
variabile adjuncte, dar înainte de asta este convenabil să obținem legea de control optimal, optu.
2.5.7 Legea de control optimal în problema masei minime
Problema abordată în această lucrare este problema masei minime, și anume, găsirea unei legi de control
optime, astfel încât satelitul sa fie transferat de la o orbită la alta folosind cât mai puțin combustibil. Ecuațiile
diferențiale rămân aceleași ca și în cazul unei probleme de timp minim de transfer. Diferența apare la indexul
de performanță, care va deveni fJ m . Ca și consecință, condițiile la limită vor suferi de asemenea
schimbări. Problema rezultată poate fi descrisă astfel:
0 1
0 0
0 0
0 0min
( , )
( , ) ( , )
| |
( ) ( )
( ) ( )
( ) ( )P SP
f f
f
f
fimize mf
dx uF x Ldt m
dL ug x L g x Ldt m
dm u
dt g I
x t x x t x
L t L L t liber
m t m m t liber
t fix
(2.46)
Deoarece sistemul de ecuații diferențiale rămâne același, iar indexul de performanță nu include nicio
integrală, hamiltonianul rămâne același ca în cazul timpului de transfer minim.
Următorul pas este derivarea funcției Hamiltoniene în raport cu thrustul, u¨, pentru a obține legea optimă de
control. Diferența principală față de problema de timp de transfer minim, este ca modulul thrustului nu mai
este constant și egal cu 1, deoarece acea reducere a masei trebuie să fie direct legată de folosirea a mai puțin
combustibil.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
19
Conform calculului variațiilor, legea de control optimal a unei probleme de optimizare este obținută
minimizând Hamiltonianul, H, în funcție de vectorul de control, u. Cu alte cuvinte este necesară aplicarea
ecuației 2.22.
1 1
0 00
1x m
x
mL LH uFu g m g u
uFu m
(2.47)
Observați că, în ecuațiile anterioare, derivând în funcție de fiecare component al vectorului, rezultă un vector;
în cazul nostru 1 2 3, ,u u u u însemnând că:
1 2 3, ,H H H H
u u u u
Având asta în minte este simplu de realizat că:
2 2 2
1 2 3
2 2 2 2 2 2 2 2 2
1 2 3 1 2 3 1 2 3
1 2 3
3 1 2
2 2 2 2 2 2 2 2 2
1 2 3 1 2 3 1 2 3, ,
2 2 21 1 1, ,2 2 2u u uu
u u
u u u u u u u u u
u u u
u u u
u u u u u u u u u
u
u
cum apare in ecuația 2.46.
Deoarece ueste un vector normalizat, singura informație de care avem nevoie este direcția, deci din ecuația
2.46 vom extrage:
x
opt
xFu
F
(2.48)
Cu următoarele condiții la limită:
1
0x
m
x
mF
u dacăm
F
u dacăm
(2.49)
Asta înseamnă că propulsia trebuie să fie pornită și oprită instantaneu pentru a obține legea de control pentru
a rezolva problema masei minime. Acest comportament mai este numit și bang-bang.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
20
Pentru a rezolva problema apariției comportamentului bang-bang va trebui să schimbăm ecuația fluxului de
masă (2). Astfel ecuația fluxului de masă va avea forma:
2 1
0(1 )Ldmu uds g (2.50)
Modificând funcția Hamiltoniana astfel:
2 1
0(1 )¨x mL uH F u ug m
(2.51)
Derivând ecuația anterioară în raport cu u , obținem legea de control optimal:
12 (1 )x mHF uu m
(2.52)
Care, după cum am menționat anterior, legea de control optimal trebuie sa fie paralelă cu produsul
variabilelor de stare adjuncte și matricea F:
x
xFu
uF
(2.53)
Introducând acest rezultat în Hamiltonian, obținem:
2 1
0(1 )x mL uH F u ug m
(2.54)
Expresie ce poate fi derivată în raport cu u, pentru a obține condiția pentru thrust optim pentru modulul
thrustului:
2 (1 ) 0
2 (1 )x
m
x
mFHuu m
F
muu
(2.55)
Această expresie scoate în evidență singularitatea ce apare pentru =1. Totuși, știind că 0,1u, condiția
anterioară poate fi modificată într-o formă mai versatilă:
0 0
2 (1 ) 0 1x
m
x
mF
Dacă atunci um
F
Dacă u atunci um
(2.56)
Aceste condiții apar impunând o valoare pozitivă mai mică ca 1 pentru u. Dacă niciuna din condiții nu este
adevărată, atunci ecuația (2.55) va fi folosită. Astfel, pentru cazul crucial în care =1, ambele condiții sunt
legate în una singură, fiind imposibil să nu fie indeplinită cel puțin una.
Funcția Hamiltoniană diferă de cazul problemei masei minime. Asta înseamnă că ecuațiile obținute anterior
pentru derivarea ecuațiilor adjuncte vor trebui modificate, derivând ecuația (2.54) în raport cu cele opt
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
21
variabile de stare 1, , ,x m L, luând în calcul rezultatul condițiilor (2.56): ( , , , )x m u u x m . Rezultând
ecuațiile următoare:
0 1
0 0( ) 2 (1 )ix x
x m x xF dg LH u u u uuds g g m m
(2.57)
2 (1 )x x
x
i m xF F uuxm F
(2.58)
1
2
02 (1 )m m
x m m md u L uF u uuds g m m
(2.59)
22 (1 )x
m
mF uum m
(2.60)
1 2
0 11(1 )L
x md u HF u uds g m L
(2.61)
0d
ds (2.62)
Ecuațiile de mai sus completează sistemul de șaisprezece ecuații diferențiale ce definesc problema masei
minime.
Cât despre condițiile la limită, șaisprezece sunt necesare pentru a închide problema. Orbita inițială și finală
oferă zece, și timpul final una. 1L este liber la ambele capete, și ca rezultat adjunctele sunt 0 la începutul și
sfârșitul transferului. Condițiile finale la limită vin din masa finală care este liberă, astfel condiția la limită
pentru adjunctă devine: ( ) ( )/ ( ) 1m f f ft J t m t . Problema este astfel închisă, și o soluție poate fi
găsită.
Se poate observa că, la fel ca în cazul problemei timpului minim de transfer, variabila de timp , nu apare
explicit în funcția Hamilton, rezultând ca valoarea adjunctei este constantă pe parcursul transferului.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
22
2.5.8 Ecuațiile diferențiale adjuncte
Introducând legea de control optimal în hamiltonian acesta ia următoarea formă:
1
0
2
1
0
1
01
1
1x
x m
x
x
m
x
x mF LH Fg m F
F L
g m F
LFg m
(2.63)
Acum este posibil să derivăm Hamiltonianul în funcție de fiecare variabilă de stare, obținând astfel ecuațiile
diferențiale pentru variabilele adjuncte.
//
1 0 1
2
0 0, 1,…,5
1ix
i
x x
x m
xd Hid x
F F L g LFg m g m F
(2.64)
1
2
0m
xd LHFd m g m
(2.65)
1
1 0 11 1 L
x md H HFd L g m L
(2.66)
0d H
d
(2.67)
Aceste ecuații formează un sistem cu 16 ecuații diferențiale și 16 variabile, opt variabile de stare 1, , ,x m L
și adjunctele lor
1, , ,x m L
. Totuși dacă examinăm mai atent sistemul putem observa ca variabila
timpului, , nu apare explicit în nici una din ecuații, însemnând că și adjuncta sa este constantă pe întreaga
traiectorie.
Din aceste informații putem extrage ca adjuncta timpului este -1. Combinând aceasta cu decuplarea ecuației
diferențiale a lui , rezultă că avem nevoie doar de 14 ecuații diferențiale pentru a rezolva problema de
control optimal. Desigur, când timpul de transfer este dorit, și ecuația sa diferențială trebuie propagată.
Pentru a completa problema, avem nevoie de cele 14 condiții la limită. Orbita inițială si masa sunt știute,
0 , 1,…,5ix s i și 0m s. Adițional, orbita dorită este de asemenea definită, 1 , 1,…,5ix s i . Cele
3 condiții ramase sunt obținute pentru adjuncte. Masa este liberă la extrema finală, deci, conform calculului
variațiilor, valorile adjunctei 1 / 0ms m deoarece pentru problema noastră ft. În final, 1L
este liberă la ambele exteme ale traiectoriei, deci: 1 1(0) (1) 0L L .
Această problemă este prea complexă pentru a fi rezolvată analitic. În urmatorul capitol vor fi descrise
metodele numerice folosite pentru a o rezolva rapid și precis.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
23
3 Metode numerice
3.1 Introducere
Problema bilocală enunțată anterior este unul din cele mai complexe sisteme de rezolvat. Prin urmare, pentru
rezolvarea acesteia avem nevoie de un mediu de calcul puternic, și anume, un computer.
Pentru a rezolva o problemă bilocală, există două abordări posibile: estimarea valorilor inițiale, propagarea
până la sfarșit și ajustarea valorilor inițiale pentru a satisface constrângerile la limită și itera, metodă numită
popular „shooting method”, sau discretizarea întregii traiectorii, estimând variabilele de stare în fiecare punct
al discretizării, și corectarea iterativă a estimărilor pentru a satisface condițiile la limită, această metodă
numindu-se „relaxation method”.
Metodele prin relaxare necesita mult mai multe resurse decat metoda prin tragere, și în general sunt folosite
doar atunci când metoda prin tragere nu poate găsi o soluție. Problema principală cu metoda prin tragere este
găsirea unui (initial guess) valori inițiale atât de bun încât să permită convergența.
Pentru a rezolva această problemă, o metodă în doi pași este folosită: se face o căutare pentru a găsi toate
soluțiile posibile pentru valorile inițiale și imediat ce o valoare relativ bună a fost găsită, o nouă reiterație mai
precisă este facută până ce valorile inițiale ating o acuratețe de aproximativ 1210.
Metoda aleasă pentru a efectua o cautare rapidă în aria valorilor inițiale optime este Evoluția Differențială
[7]. Aceasta este o metodă simplă și robustă pentru a rezolva ecuații cu variabile multiple evaluând funcția
cu variabile aleatorii si, depinzând de cât de aproape suntem de rezultatul dorit, amestecarea acestora.
Odată ce au fost obținute rezultate inițiale bune, este timpul sa aplicăm metoda tragerii, de obicei prin metoda
Newton-Raphson. Dacă valorile obținute inițial sunt destul de bune, această metodă va ajunge la acuratețea
dorită în doar câteva iterații.
Desigur, în spatele ambelor metode, care, atunci când sunt combinate ne permit rezolvarea problemei, stă un
alt proces numeric, și anume propagarea sistemului de ecuații diferențiale de la starea inițială la starea finală.
Un algoritm Runge-Kutta de gradul 4 ne permite o acuratețe sporită fără o solicitare intensă a computerului.
Este de asemenea important să nu uităm că, la fel cum a fost menționat anterior, o metodă buna pentru a
reduce erorile seculare în timpul propagării este să efectuăm o mediere a sistemului diferențial de ecuații.
Din moment ce aceste ecuații sunt neliniare, o metodă numerica de integrare Euler-Lagrange este
implementată pentru a efectua medierea.
În ultimul rând, un instrument este necesar pentru a rezolva probleme mult mai complexe cum ar fi
minimizarea masei finale. Pentru aceasta, metode de omotopie au fost folosite.
3.2 Algoritmul Evoluției Diferențiale
Acest algoritm este o schemă pentru a minimiza o funcție obiectiv, continuă, cu variabile multiple. Conform
[7], ideea de bază din spatele acestui algoritm este de a genera câteva seturi de vectori cu valori aleatorii,
evaluând funcția obiectiv cu valorile obținute și apoi a crea noi vectori pentru a adauga o pondere dintre doi
vectori intr-un al treilea. Daca acest nou vector da un randament mai bun și funcția obiectiv scade, acesta este
salvat, în caz contrar, este respins.
Acest proces este supranumit procesul de selecție lacom, deoarece păstrează doar vectorul care minimizează
funcția obiectiv scufundând astfel funcția într-un minim local, nereușind să ajungă la minimul global dorit.
În încercarea de a evita acest comportament schema evoluției diferențiale efectuează câteva evaluări în
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
24
paralel, astfel încât vectorul ce ne oferă minimul local sa poată fi înlocuit de un alt vector care oferă încet
convergență către un minim global.
Printre avantajele Evoluției Diferențiale, sunt demne de menționat, viteza de convergență și robustețea
datorită paralelismului și ușurința de folosire. Pe lângă acestea, codul sursă este gratuit și poate fi găsit în
numeroase limbaje de programare (cum ar fi Java, Matlab, Fortran, C++).
În mod oficial, Evoluția Diferențială este o metodă de căutare paralelă directă ce folosește un vector de n
parametri , 1,…,G
ix i n pentru orice generație G. Parametrul n rămâne constant pe durata întregului proces.
Populația inițială este aleasă aleatoriu, de obicei asumându-se o distribuție uniformă. Vectorii sunt apoi
combinați folosind o pondere, și vectorul obținut este păstrat doar dacă este mai aproape de soluție decât cel
anterior.
Mai multe formule au fost încercate pentru a genera un nou vector. Una din cele mai simple și totuși de succes
formule este următoarea:
1; 1,G G G G
i r s iv x F x x r s t n (2.68)
Unde F este o constantă reală, pozitivă, diferită de 0 ce controleaza amplificarea diferenței. Pentru fiecare
vector al populației, trei vectori sunt aleși aleatoriu printre componentele rămase ale populației și combinate
urmând expresia anterioară.
Pentru a crește diversitatea vectorilor următori, un concept numit crossover sau mutație este introdus. Asta
înseamnă că un nou vector, u, este creat prin amestecarea la întâmplare a componentelor lui xși v.
Matematic, această operație poate fi scrisă astfel:
1 1 1
,1 ,,…,G G G
i i i Du u u (2.70)
Cu
1
,1
,
,,…, 1
, 1,G
i jG
i jG
i jv pentru j n n L
u
x j D
(2.71)
n și L fiind scalari generați aleatoriu, urmând întotdeauna regula 0 < n < L < D.
Asta înseamnă ca, pentru fiecare parametru al vectoruluiix, un nou vector iu este creat folosind câteva
componente din vectorul iv și restul din ix.Pentru fiecare membru al populației, numărul și locul
componentelor alese din iv sunt diferite.
Odată ce 1G
iu a fost creat, este comparat cu G
ix și dacă oferă o valoare mai mică a funcței obiectiv, este
folosită pentru urmatoarea generație astfel: 1 1G G
i ix u . În caz contrar, vectorul anterior este reținut:
1G G
i ix x .
Funcția obiectiv este unul din factorii ce depind în mod special de problema pe care o rezolvăm. În cazul de
față, funcția este minimă atunci când orbita de la punctul final al propagării este egală cu orbita obiectiv. Prin
urmare, cea mai simplă și totuși de efect funcție obiectiv poate fi definită astfel:
1522 2
1f obj
i i m L
iF x x
(2.72)
Această funcție obiectiv este validă ca orice altă funcție continua a cărui minim este obținut doar atunci când
orbita finală este egală cu orbita obiectiv, și variază ușor în apropierea acestui minim.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
25
Observați că aceasta funcție obiectiv nu este directă, deoarece valorile de intrare sunt valori găsite anterior și
folosește variabile de stare finale. Asta înseamnă că propagarea trebuie reluată de fiecare dată când funcția
este evaluată.
Această metodă implică generarea a mai multor vectori pentru a ajunge la soluția dorită, având fiecare câteva
chemări către funcția obiectiv și prin urmare câteva propagări. Cu toate acestea, aceasta este o metoda foarte
rapidă și folositoare chiar și pentru această problemă complexă.
3.3 Shooting Method
Metoda tragerii este un proces numeric pentru rezolvarea problemelor bilocale. După cum ne putem aminti
din capitolele anterioare aceasta poate fi enunțată astfel: găsiți soluția pentru n ecuții diferențiale cuplate:
, , 1,…,i
idy xf x y i ndx (2.73)
Ce satisfac 1n condiții la limită la punctul de pornire 1x de forma:
1 1 1, 0, 1,…,jB x y j n (2.74)
Și 2 1n n n condiții la limită pentru punctul final 2x
2 2 2, 0, 1,…,kB x y k n (2.75)
Procesul de găsire a soluției prin această metodă este următorul:
1. Estimarea tuturor variabilelor de stare la starea inițială 1( )y x. Aceste valori trebuie să respecte
condițiile la limita în acel punct, 1B
.
2. Integrarea sistemului de ecuații diferențiale ordinare urmând Metoda Valorilor Inițiale (ce va fi
enunțată în următorul subcapitol) până la punctul final, 2x.
3. Evaluarea conformității cu condițiile la limită la timpul final, 2B
.
4. Ajustarea iterativă a valorilor absolute, neconstrânse de condițiile la limită 1B
, astfel încât
condițiile la limită finale 2B
pot fi respectate.
Problema este redusă astfel la o problemă de găsire a răcinii ce poate fi enunțată de exemplu, ajustați
parametrii liberi inițiali astfel încât condițiile la limită finale să fie respectate.
Cea mai comună metodă numerică implementată pentru rezolvarea problemei găsirii rădăcinii, este Newton-
Raphson, care încearcă să aducă la 0, 2n funcții cu 2n variabile.
Pentru o problemă bilocală, știm ca sunt 1n condiții la momentul inițial și 2 1n n n variabile libere la
momentul inițial. Un vector v de dimensiune 2n poate fi format cu aceste variabile libere, astfel încât având
v împreună cu condițiile la limită inițiale 1B
, ne permit să construim un set de ecuații diferențiale complet
1y x.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
26
Având toate valorile inițiale, este posibil să integrăm ecuațiile diferențiale ordinare pentru momentul de timp
inițial și final, practic rezolvând problema valorilor inițiale. Conceptual, din 1y x sunt obținute variabilele
finale 2y x.
În continuare, un nou vector de dimensiune 2,n d
este definit, numit vectorul de discrepanță, folosit pentru a
măsura cât de departe sunt valorile obținute de soluția dorită. Definirea acestui vector este destul de simplă
și anume, folosirea condițiilor la limită la punctul final, 2 2 2, , 1,…,k kd B x y k n , dar sunt și alte definiții
acceptabile, atâta timp cât d
este zero dacă și numai dacă toate condițiile la limită la momentul 2x sunt
satisfăcute.
Problema în rezolvarea Newton-Raphson este găsirea valorilor lui vastfel încât d
sa fie zero. Acest lucru
este obținut prin calculul soluției:
A v d (2.76)
Și adăugând corectia v valorilor lui v. Procesul este repetat iterativ scăzând gradual valoarea lui v (în
cazul în care avem convergență) în timp ce condițiile la limită sunt mai aproape de a fi împlinite.
Matricea pătrată A este construită din:
,i
i j
jdAv
(2.77)
Valori ce nu pot fi obținute în mod analitic. În loc, la fiecare iterație aceste derivate parțiale sunt estimate
folosind diferențe finite.
2 2 1 1,…, ,…, ,…, ,…,i j j n i j ni
j jd v v v v d v v vd
v v
(2.78)
Prin urmare, o iterație a metodei prin tragere necesită 21n integrări ale sistemului de ecuații diferențiale: o
integrare pentru a evalua discrepanțele și 2n pentru a obține matricea A.
Este crucial să selectăm o valoare adecvată pentru creșterea jv. Dacă este prea mică, derivata nu va avea
niciun sens datorită pierderii acurateții, fie prin rotunjire, fie prin acuratețea limitată a procesului de
propagare. În caz contrar, dacă jv este prea mare, aproximarea făcută de ecuația anterioară nu va arăta
comportamentul real al derivatei.
Din fericire, deși valoarea lui jv trebuie aleasă, manual și aleatoriu, nu este foarte greu să găsim o valoare
destul de mică cu ajutorul căreia să obținem convergența, doar în cazul în care problema nu este potrivită
pentru această metodă, caz în care jv nu mai are oricum niciun sens.
Metoda tragerii va converge daca ecuațiile diferențiale ordinare vor varia foarte puțin în jurul soluției și dacă
valorile inițiale sunt destul de bune. Ca orice altă aplicație Newton-Raphson, această metoda va converge
pătratic în jurul soluției.
După cum enunțam anterior, convergența globală slabă este principala problemă a acestei metode. Totuși,
există un număr de mecanisme ce pot rezolva această problemă [8]. Una dintre cele mai cunoscute modificări
aduse metodei Newton este adăugarea unui coeficient creșterii jv, calculat astfel încât valoarea scalară
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
27
d d
descrește strict monoton în timpul iterațiilor. Asta înseamnă că după fiecare pas al iterației valoarea
d d
va scădea.
Fără a mai intra in detalii matematice greu de digerat, următorul proces este folosit pentru a obține o pondere
k pentru fiecare kv la fiecare pas al iterației k:
Calculăm k k k kh v d d , kd
fiind vectorul de discrepanță obținut cu ajutorul lui kv.
Găsirea celui mai mic număr întreg 0m, astfel încât 2m
k k k k kh v h v v
2m
k, și în final 1k k k kv v v
3.4 Propagarea numerică
Secțiunile anterioare, 3.2 și 3.3 ne oferă o descriere la nivel superior a proceselor necesare pentru rezolvarea
problemelor bilocale. Totuși, ambele procese se bazează pe un al treilea mecanism pentru a obține starea
finală când este dată starea inițială. Prin urmare, rezolvarea problemei valorilor inițiale, sau propagarea
variabilelor de stare de la momentul inițial la momentul de timp final, este necesară.
Problema valorilor inițiale poate fi văzută ca o problemă bilocală, în care condițiile la limită sunt impuse într-
un singur punct – de obicei momentul inițial al integrării. Sistemul de ecuații diferențiale ordinare, ramâne
deci același:
,i
idy xf x ydx (2.79)
singurele schimbări aparând doar în condițiile la limită.
Ideea de bază a acestei metode de rezolvare este rescrierea operatorilor diferențiali dy si dx ca diferențe finite,
y si x , și apoi înmulțind funcțiile if cu x. În acest fel, formule algebrice sunt obținute pentru a măsura
schimbarea în derivate când este adăugată în x o mică creștere sau un mic pas. Dacă aceste creșteri sunt destul
de mici, atunci o bună aproximare a derivatelor este obținută.
Sunt mai multe metode numerice funcționează după acest principiu. Cele trei metode numerice de propagare
principale sunt:
Runge-Kutta: câțiva pași din metoda Euler sunt combinați pentru a obține infomații despre
intervalul propagării, pentru a potrivi mai apoi o expansiune a seriei Taylor.
Extrapolarea Richardson: rezultatul calculat este extrapolat către o valoare care ar fi putut fi
obținută daca pasul ar fi fost mult mai mic, de preferat zero.
Predictor-corector: soluția unui pas este folosită pentru a extrapola către urmatorul; extrapolarea
este corectată cu informația noii derivate.
Din cele trei metode, metoda predictor-corector este cel mai greu de implementat, datorită dificultăților de
pornire. Mai mult de atât, doar cele mai sofisticate metode predictor-corector sunt competitive cu celelalte
două metode și doar atunci când funcțiile sunt foarte omogene. Din cele două rămase, metoda Runge-Kutta
este mult mai robustă, totuși nu este neapărat mai eficientă. Pentru a asigura pe cât posibil convergența
problemei valorilor inițiale, metoda Runge-Kutta a fost aleasă.
Miezul metodei Runge-Kutta, este, după cum am enunțat anterior, metoda Euler. Matematic, aceasta este
enunțată astfel:
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
28
2
1 ,n n n ny y hf x y O h (2.80)
Impementarea acestei formule nu este cea mai bună abordare, dat fiind ca pentru acelasi pas h alte formule
oferă o acuratețe și o stabilitate mai bună.
De departe, cea mai folosită și fără îndoială cea mai folositoare metodă este metoda Runge-Kutta de ordinul
4:
1
1
2
2
3
4 3
5 3 1 2 4
1,
,2 2
,2 2
,
6 3 3 6n n
n n
n n
n n
n nk hf x y
khk hf x y
khk hf x y
k hf x h y k
k k k ky y O h
(2.81)
Formule Runge-Kutta de ordin superior sunt rareori folosite, de obicei deoarece pentru m ordin sunt necesare
mai mult de m evaluări ale derivatelor.
Această formulă poate fi implementată cu ușurință cu un pas constant h de la momentul inițial la momentul
final. Totuși, având în vedere că folosim computerul, acest lucru va fi departe de a fi optim. Funcțiile pot
avea mici zone unde se schimba rapid, separate de zone foarte mari unde variază foarte puțin. Ar fi mult mai
interesant să evaluăm eroarea după fiecare pas și apoi sa schimbăm pasul pentru a ne menține pe lângă
valoarea dorită. Cu alte cuvinte, un pas adaptiv ar putea scădea timpul de găsire a soluției fără a sacrifica
acuratețe, chiar dacă asta înseamnă să evaluăm de mai multe ori derivatele.
Cu Runge-Kutta de ordinul 4, una din cele mai directe metode de impelementare a pasului adaptiv este
dublarea: folosirea pasului de 2h sau a doi pași de mărime h in același punct. Notam cu 1y rezultatul cu un
singur pas și 2y rezultatul propagării în doi pași, diferența respectivă în raport cu soluția exactă y fiind:
5
5 6
1
5
5 6
22 25!
2 25!y xy x h y h O h
y xy x h y h O h
(2.82)
Erori de trunchiere pot fi estimate prin diferența celor două propagări:
2 1y y (2.83)
Această valoare trebuie menținută sub o anumită valoare pentru a asigura acuratețea și viteza procesului.
Acest lucru poate fi obținut făcând legătura dintre
și h. Deoarece există o expresie pentru eroarea de
trunchiere pentru formula Runge-Kutta de ordinul 4, aceasta poate fi folosită pentru a îmbunătăți propagarea.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
29
6
2215y x h y O h
(2.84)
Obținându-se astfel o expresie de ordinul cinci. Dar, având în vedere ca nu este nicio cale de a monitoriza
expresia erorii de trunchiere, pentru a verifica aceasta poate fi folosită direct în locul lui 2y și a verifica dacă
sunt îmbunătățiri în rezultatele propagării. În orice caz, nu poate degrada precizia.
Din ecuațiile anterioare putem observa că
poate fi scalat ca 5h. O relație poate fi stabilită între eroarea 1
produsă de un anumit pas, 1h, si pasul 0h care ar fi trebuit să fie folosit pentru a obține eroarea dorită 0
:
0.2
0, 0
1 1,min , 1,…,i
ihi nh (2.85)
Luând în considerare ecuația de mai sus, este ușor să obținem un pas mai mic dacă eroarea obținută este prea
mare, și respectiv să sugerăm un pas mai mare pentru următorul pas dacă eroarea este prea mică. Putem
observa că, din moment ce
are la fel de multe componente ca și variabilele de stare din problemă, cea mai
restrictivă variabilă va ajusta pasul.
Singura problemă ramasă este sa obținem o valoare funcțională pentru 0
, acuratețea dorită. Una din cele
mai directe opțiuni este folosirea valorilor variabilelor de stare la pasul curent, y, scalat cu un număr mic
. Adăugând estimarea variației lui y, o valoare mai constantă și sensibilă a erorii este obținută. O buna
definiție pentru vectorul erorilor este:
0,i i iy h f (2.86)
În acest moment avem toate elementele necesare pentru a începe propagarea. Pașii sunt repetați, crescând
variabila independentă, pană se ajunge la momentul final.
3.5 Medierea ecuațiilor
Sistemul de ecuații diferențiale ordinare obținut în capitolul anterior, deși precis este precis pentru modelul
ales, este posibil sa nu se comporte corect în urma popagarii numerice. Erori seculare pot fi acumulate,
datorită prezenței unei variabile ce variază foarte rapid în raport cu restul, și anume L. Acest comportament
poate fi atenuat prin medierea efectului lui L asupra sistemului pe întreaga revolutie. Matematic, acest lucru
poate fi descris astfel:
21, , ,2L
LF y x F y L x dL
(2.87)
unde x este variabila independentă, y este format din toate variabilele de stare, F
este format din toate
ecuațiile diferențiale, iar F
reprezintă ecuațiile diferențiale mediate.
Privind către ecuațiile de interes, poate fi dedus rapid acestea nu pot fi mediate analitic. În schimb, acestea
pot fi integrate.
Cea mai simplă abordare este să presupunem că avem două puncte la distanțe egale pe parcursul integrării.
Mai multe metode știute, de exemplu Simpson, urmează această abordare. De asemenea, setând spațiul ca
valoare liberă sau presetată, ne permit să obținem formule de ordin superior cu același număr de evaluări.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
30
Formulele create din această filozofie sunt denumite formulele Gaussian patratice. Forma lor generală fiind
umătoarea:
1b N
i i
i aW x f x dx w f x
(2.88)
Funcția W(x) poate fi folosită pentru a îndepărta singularitățile integralei funcției ce trebuie mediată., dar în
general este mai simplu să ne asumăm W(x) = 1. Aceasta alegere este știută ca integrarea Gauss-Legendre,
când a = -1 și b = 1.
Metoda de calcul a absciselor ix și ponderile iw a ecuației anterioare necesită polinomiale ortogonale.
Un produs scalar a doua funcții f și g față de funcția ponderilor W, este definită astfel:
b
af g W x f x g x dx (2.89)
un număr independent de x. Două funcții sunt ortogonale dacă funcția lor scalară este zero.
Este posibil să gasim un set de funcții polinomiale ortogonale ce includ un singur polinomial de ordin j pentru
fiecare j = 0,1,2,… în umătoarele relații:
0
1 1
1 1
1 1
1
1 11
i i i i
i i i
i i i ip x
xp p
p x x p x
p p
xp p xp p
p x x p x x p x
p p p p
(2.90)
Poate fi dovedit [8] că polinomialul np x are n rădăcini diferite în intervalul (a,b). Mai important este faptul
că, aceaste rădăcini vor fi abscisele a n puncte în cuadratura Gausiană, cu funcția ponderilor W în intervalul
(a,b). Acest lucru este crucial pentru cuadraturile Gausiene, lăsând utilizatorul să găsească abscisele în funcție
de caz.
Având abscisele, tot ce ne rămâne este să găsim echivalentul acestora în n ponderi. Pentru acestea o nouă
secvență de polinomiale poate fi definită, jx , de forma:
0
1
1
1 1
1 10
b
a
i i i i
i i i
i i i ix
dpx W x dxdx
xp p xp p
x x x x
p p p p
(2.91)
polinomialele jx urmează aceeași recurență ca și jp, dar cu valori de început diferite, de asemenea având
și un ordin inferior. Cea mai importantă folosință a acestor polinomiale este următoarea relație:
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
31
1
, 1,…,n i
i n idp xw x i ndx
(2.92)
ce ne permite să calculăm cu ușurință ponderea pentru fiecare cuadratură Gausiană. Pentru fiecare caz de
interes din cuadraturile Gauss-Legendre, totuși, o expresie mai simplă poate fi folosită pentru aceste ponderi:
22 /2
1i
i n iw
x p x
(2.93)
unde /
np este o contracție a /ndp dx.
În acest moment avem toate elementele necesare pentru mediere. Acesta este un pas ce trebuie făcut la fiecare
pas al propagarii pentru fiecare funcție propagată. Singura întrebare rămasă este, care este numarul minim de
abscise ce micșorează erorile seculare. Conform [1], numărul de abscise poate varia de 12 până la 96,
depinzând de frecvența relativă a variabilelor de stare.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
32
4 SOFTWARE
Capitolele anterioare au fost dedicate pentru a prezenta problema unui transfer orbital optimal folosind
propulsie scăzută și pentru a detalia metodele necesare rezolvării unei astfel de probleme. Acest capitol
descrie software-ul dezvoltat pentru rezolvarea acestei probleme.
Programul a fost numit LTO, acronym pentru Low-Thrust Optimizer. Limbajele de programare folosite
pentru acest program sunt Matlab, Vb.NET and C#. Acest program are rol demonstrativ și nu are
implementate toate metodele din capitolele anterioare.
4.1 Prezentare generală
Principalele sarcini ale programului sunt următoarele:
Citirea informatiilor inițiale despre masa satelitului, masa de combustibil folosită, detalii despre
sistemul de propulsie (thrust și ISP), numărul maxim de iteratii, eroare maximă a constrângerilor și
pasul datelor de ieșire.
Citirea fișierului traiectoriei inițiale ce conține un transfer aproximativ și procesarea fișierului
pentru a defini arcurile traiectoriei.
Folosirea unui propagator bazat pe Runge-Kutta, folosit de Newton-Raphson pentru a obtine
rezultatele necesare, dar și pentru obținerea masei finale și a tuturor rezultatelor intermediare. Ce
este de remarcat este calculul legii de control care are loc în aceste iterații.
Convertirea informației obținute într-un format ușor de citit de către utilizator (redimensionarea
variabilelor de stare, reprezentarea grafică a rezultatelor).
4.2 Programul principal
Programul principal “LTO” se ocupă cu inițializarea modulelor auxiliare necesare pentru pregătirea datelor
inițiale folosite de către programul de calcul.
Prima acțiune a programului este reprezentat de interogarea eventualelor erori în datele de intrare. Apoi, dacă
în urma interogarii nu sunt descoperite erori, acesta va scrie într-un fișier extern toate datele de intrare
adăugate de către utilizator ce vor fi folosite în programul de calcul.
Următorul pas este pornirea programului de calcul ce va începe procesul prin citirea fișierului scris de către
programul principal prin intermediul modulului config_file.
După aceea, modul Optimizer este initializat, modul ce va efectua toate procesele necesare rezolvării
problemei de control optimal.
În pasul următor, toate datele obținute ce au legătură cu traiectoria optimă sunt trimise în modulul Output
pentru a fi procesate în vederea reprezentării grafice ale rezultatelor.
La finalul procesului, programul de calcul este închis cu mesajul OPTIMIZATION COMPLETE, iar
rezultatele vor putea fi vizualizate în programul principal.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
33
Fig. 2 – Prezentarea generală a procesului
4.2.1 WRITE PANEL
Modulul WRITE PANEL este folosit pentru scrierea datelor de intrare necesare programului de calcul.
Acesta va scrie un fisier „CONFIG” în interiorul folderului simulării, bazat pe inputul definit de către
utilizator, ce va fi mai apoi folosit de către programul de calcul pentru a începe procesul de optimizare.
Formatul fișierului scris de către modulul WRITE PANEL este cel din următorul tabel.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
34
Rând Definiție
1 ISP (s)
2 Thrust (N)
3 Masa satelitului (kg)
4 Masa de combustibil (kg)
5 Violarea maximă a constrângerilor
6 Număr de iterații
7 Output spacing (s)
8 Degradarea motorului (-)
9 Timpul de simulare (zile)
4.2.2 PLOT DATA AND VISUALIZE RESULTS
Acest modul este folosit pentru prelucrarea fisierelor de iesire pentru a rescrie fișierele într-un format ușor
de interpretat pentru utilizator.
Acest modul este folosit în special pentru reprezentarea grafică a rezultatelor făcând conexiunea dintre
fișierul de ieșire al programului de calcul și interfața grafică.
De asemenea, acest modul oferă utilizatorului pe lângă posibilitatea de a vizualiza rezultatele în interiorul
interfeței grafice, și vizulaizarea acestora în afara ei generând următoarele fișiere în format PNG:
Fișier Definiție
Eccentricity.png Eccentricitatea vs Timp (zile)
Inclination.png Înclinația vs Timp (zile)
Mass.png Masa (kg) vs Timp (zile)
Hodograph.png Vizualizare 3D: Vx (km/s), Vy (km/s), Vz
(km/s)
SMA.png SMA, Altitudinea la Perigeu și Apogeu (km)
vs Timp (zile)
Trajectory.png Vizualizare 3D: X (km), Y (km), Z (km)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
35
Thrust.png Thrust XYZ și Modulul (kN)
Thrust_Direction.png Direcția Thrustului vs Timp (zile)
xyzR.png Poziția XYZ și Raza (km) vs Timp (zile)
uvwV.png Viteza XYZ și Viteza Parcursului de zbor
(km/s) vs Timp (zile)
4.2.3 READ INPUT
Acest modul are următoarele sarcini:
Citirea infomațiilor din panul principal al interfeței grafice
Citirea traiectoriei de referință
Citirea informațiilor hard codate despre corpul central și alți paramentri folosiți în metodele numerice
4.2.3.1 READ PANELS
Acest sub-modul este folosit pentru citirea informațiilor din panoul principal al interfeței grafice enunțate în
descrierea modulului WRITE PANEL (4.2.1), dar și pentru citirea inputurilor hard codate în interiorul
programului de calcul.
Acest modul permite folosirea programului de calcul în afara intefeței grafice, permițându-i utilizatorului să
folosească programul în CMD. Acest element constitue un ajutor imens, permițând utilizatorului integrarea
executabilului într-un proces automatizat mult mai complex.
Pentru a folosi programul în afara interfeței grafice utilizatorul trebuie să respecte cu strictețe formatul
fișierului CONFIG, dar și denumirea traiectoriei de referință și formatul acesteia.
Inputurile următoare sunt hard-codate pentru a facilita experiența utilizatorului cu programul:
Parametru Definitie Valoare
MU Constanta gravitațională a Pământului conform
Goddard Earth Models 5&6 398601.3 (km^3/s^2)
RE Raza Pământului la Ecuator 6378.144 (km)
F Earth flatenning 1 / 298.257
C Viteza luminii 299792.458 (km/s)
RAD Factor de conversie din grade în radian atan2( 1 , 1 ) / 45
DEG Factor de conversie din radian în grade 1/RAD
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
36
PI Numărul Pi 4 * atan2( 1 , 1 )
EPSILON Număr mic 1e-16
dt0 Pasul maxim de integrare 90 s
JD Julian Date 2451544.5
MUS Constanta gravitațională a soarelui 132712440018 (km^3/s^2)
MUM Constanta gravitațională a lunii 4902.8000
MCFG Schimbarea minimă a funcției cost 1e-12
atm_dens_fcn Calculul densității atmosferei la altitudini mari
considerând activitatea solară Conform Orbit &
Constellation Design and
Management, JR Wertz
p68-69
revperarc Număr folosit pentru crearea numărului de arcuri 13/139
4.2.3.2 READ INITIAL GUESS
Acest modul este folosit pentru citirea traiectoriei de referință provenită din orice alt software capabil de
producearea unui astfel de transfer orbital. Foarte important este ca formatul acestui fișier să fie cel specificat
în următorul tabel. Acest format a fost folosit deoarece este suportat de software-ul cu care am generat
traiectoria de referință și anume LOTOS, creat de ASTOS SOLUTIONS. Pentru a putea fi citită de către
program traiectoria de referință trebuie să fie plasată în folderul scenariului, să fie în format text cu denumirea
„initial_guess.txt”, și să aibe formatul din tabelul următor:
Coloana Valoare în sistem de referință inerțial
1 MJD2000
2 Coordonate X (km)
3 Coordonate Y (km)
4 Coordonate Z (km)
5 Componentele vitezei pe X (km/s)
6 Componentele vitezei pe Y (km/s)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
37
7 Componentele vitezei pe Z (km/s)
8 Componentele thrustului pe X (kN)
9 Componentele thrustului pe Y (kN)
10 Componentele thrustului pe Z (kN)
11 Masa (kg)
4.2.4 PROCESS INITIAL GUESS
Acest modul este folosit pentru definirea arcurilor traiectoriei folosite în procesul de optimizare prin definirea
matematică a vectorului ce conține variabilele de optimizare (X_opt). Acest modul v-a crea un fișier numit
„arc_initialisation.txt” ce va conține toate arcurile procesate. Numărul de arcuri create este hard codat și este
13 revoluții la 139 de arcuri.
Formatul acestui fișier este următorul:
Coloana Valoare în sistem de referință inerțial
1 Timp (s)
2 Coordonate X (km)
3 Coordonate Y (km)
4 Coordonate Z (km)
5 Componentele vitezei pe X (km/s)
6 Componentele vitezei pe Y (km/s)
7 Componentele vitezei pe Z (km/s)
8 Multiplicatorul poziției X la începutul arcului (1/s)
9 Multiplicatorul poziției Y la începutul arcului (1/s)
10 Multiplicatorul poziției Z la începutul arcului (1/s)
11 Multiplicatorul vitezei X la începutul arcului (-)
12 Multiplicatorul vitezei Y la începutul arcului (-)
13 Multiplicatorul vitezei Z la începutul arcului (-)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
38
14 Masa (kg)
În acest modul se calculează și se afișează starea inițială și starea finală.
4.2.5 DEFINE CONSTRAINTS
4.2.5.1 CONSTRÂNGERI LINIARE DE EGALITATE
Constrângerile liniare sunt construite sunt formulate simplificat astfel:
* _A X opt b (2.94)
Aceste constrângeri sunt de două tipuri:
Constrângeri pe starea inițială
Constrângeri pe durata arcurilor
4.2.5.2 CONSTRÂNGERI ALE STĂRII INIȚIALE ORIGINALE
Aceste constrângeri sunt în număr de opt, sunt constrângeri create automat pe baza stărilor inițiale și sunt
cele de mai jos:
Poziție (X,Y,Z)
Viteză (X,Y,Z)
Masă
Timp
4.2.5.3 CONSTRÂNGERI LEGATE DE DURATA ARCURILOR
Aceste constrângeri sunt bazate pe durata arcurilor, mai exact pe raportul dintre durata unui arc și durata
totală.
Constrângerile sunt formulate astfel:
Prima constrângere:
(1 (1))* (1) (1)*( (2) (3) … ( _ ))) 0k k n arcs (2.95)
Celelalte constrângeri:
(1)* (1) ( )* ( ) 0, 2,…, _k k j j j n arcs (2.96)
Astfel, numărul total de constrângeri de egalitate liniară va fi numărul de constrangeri ale stării inițiale plus
numărul de constrângeri legate de durata arcurilor și anume:
_ _ _ _ _ 8 _ 1N LEC STATES N LEC ARC DURATION n arcs (2.97)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
39
4.2.5.4 CONSTRÂNGERI NELINIARE DE EGALITATE
Constrângerile neliniare luate în considerare sunt bazate pe valorile finale ale stărilor din traiectoria de
referință și pe defectele apărute la sfârșitul și începutul arcurilor sunt definite astfel:
Constrângeri de egalitate apărute la modulul inițial al vectorului vitezelor învecinate. Fiecare
constrângere depinzând doar de multiplicatorii vitezei (3 variabile).
Constrângeri de defect la sfârșitul tuturor arcurilor (mai puțin ultimul). Fiecare constrângere
depinzând de numarul de variabile ce definesc două arcuri învecinate (8 variabile).
Constrângeri definite de starea finală a ultimului arc. Fiecare constrângere depinzând de toate
variabilele ce definesc ultimul arc (15 variabile), durata arcului cu zbor liber (1 variabilă),
schimbarea nodului ascendent (1 variabilă) și argumentului la perigeu (1 variabilă).
4.2.6 EVALUATE CONSTRAINTS AND COST FUNCTION
4.2.6.1 EVALUAREA FUNCȚIEI COST
Acest modul este folosit pentru a evalua funcția cost, în acest caz fiind minus masa minimă, fiind calculată
astfel:
Date de intrare:
o 0m- masa inițială, în kg
o T – modulul propulsiei citit din meniul principal, în kN
o 0SPg I
o i – durata arcurilor, în secunde
o COSTSCALE = masa combustibil / 2 ;
o COSTOFFSET = -( masa initială – masa combustibil / 2 ) ;
Date de ieșire:
o Funcția cost:
0 0( ( / )* ( ))
( )/kg SP i
SCALED kg COST COSTCOST m T g I sum
COST COST OFFSET SCALE
(2.98)
4.2.6.2 EVALUAREA CONSTRÂNGERILOR MODULULUI VITEZELOR ÎNVECINATE ȘI
GRADIENTUL ACESTORA
Acest modul calculează valoarea constrângerii și gradientul modulului unitar a vectorului multiplicatorilor
vitezei, astfel:
Date de intrare:
o X_arc – Variabilele stocate în memoria RAM, dintr-un arc al traiectoriei. Dintre aceste
variabile se folosesc doar multiplicatorii vitezei X, Y, Z.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
40
Date de ieșire:
o Valoarea constrângerii:
2 2 21vmx Multiplicatorulcomponentei X avitezei
vmy MultiplicatorulcomponenteiY avitezei
vmz MultiplicatorulcomponenteiZ avitezei
CONSTRAINT vmx vmy vmz
(2.99)
o Gradientul funcției constrângerii respectând formatul componentelor datelor de intrare:
000000000…
2* …
2* …
2* …
000vmx
GRAD vmy
vmz
(2.100)
4.2.6.3 EVALUAREA CONSTRÂNGERILOR DE DEFECT DEFINITE DE SFÂRȘITUL UNUI
ARC
Acest modul calculează valoarea constrângerilor și gradientul pentru diferența dintre două stări la nodul
dintre două arcuri ale traiectoriei, astfel:
Date de intrare:
o X_arc_1 – Variabilele stocate în memoria RAM, din primul arc al traiectoriei.
o X_arc_2 – Variabilele stocate în memoria RAM, din cel de-al doilea arc al traiectoriei.
o Tf – Timpul final al primului arc
Date de ieșire:
o Constrângeri pe diferențele de poziție (X, Y , Z)
o Constrângeri pe diferența de viteză (X, Y , Z)
o Constrângere pe diferența de masă
o Constrângere pe diferența de timp
4.2.6.4 EVALUAREA CONSTRÂNGERILOR DEFINITE DE STAREA FINALĂ
Acest modul calculează valoarea constrângerilor de poziție și viteză definite de starea de la timpul final al
traiectoriei.
Date de intrare:
o X_prop_arc – Variabilele de stare stocate în memoria RAM, dintr-un arc al traiectoriei.
o Xf_prop_arc – Variabilele de stare stocate în memoria RAM, din arcul final al traiectoriei.
o Tf – Timpul de la sfârșitul ultimului arc
Date de ieșire:
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
41
o Constrângeri pe diferențele de poziție (X, Y , Z)
o Constrângeri pe diferența de viteză (X, Y , Z)
4.2.7 RUNGE-KUTTA
Propagarea se face folosind o metoda Runge-Kutta de ordinul 8. Metoda implementată este cea din referința
[10].
Date de intrare:
o Valoarea inițială a variabilei indepente
o Pasul variabilei indepente
o Funcția de derivare în timp a stărilor
o Parametrii globali ce vor fi folosiți de funcția de derivare în timp a stărilor
Date de ieșire:
o Valoarea finală a variabilei independente
o Valoarea finală a vectorului viariabilelor dependente
4.2.8 CHECK-STOP CRITERIA
Programul are implementate două funcții de oprire a procesului de optimizare, și acestea sunt definite după
cum urmează:
Valoarea absolută a constrângerilor violate este mai mică decât cea specificată în panoul
programului, în secțiunea „Max Violation”
Schimbarea funcției cost de la o iterație la alta este mai mică decat 1e-12. Acestă valoare este
hard-codată pentru a minimiza datele de intrare oferite de utilizator
4.2.9 PERFORM AUXILIARY COMPUTATIONS
Acest modul este folosit pentru transformarea stărilor în elemente Kepleriene, afișarea acestora în consola de
lucru, și pentru scrierea traiectoriei optimizate într-un fișier text. Denumirea traiectoriei optimizate este
„trajectory_optimized.txt” iar formatul acesteia este similar cu formatul specificat pentru traiectoria de
referință în sub-capitolul 4.2.3.2.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
42
5 SOFTWARE : Interfața grafică
5.1 Introducere
Acest capitol este dedicat părții de software ce este vizibil către utilizator, unde informațiile și opțiunile sunt
introduse în codul executabilului. De asemenea, fiecare parametru de intrare este explicat, incluzând valori
recomandate și ponturi oricând este posibil.
Programul este divizat în trei tab-uri :
Meniu principal (date de intrare, executare)
Consola de execuție
Afișarea rezultatelor grafice
Figure 3 – Interfața grafică
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
43
5.2 Meniul Principal
Acest panou conține un număr de nouă date de intrare și trei butoane pentru a controla procesul de optimizare,
respectiv afișare a rezultatelor.
Figure 4 – Meniul Principal
5.2.1 Date de intrare
Datele de intrare sunt datele ce sunt inserate de către utilizator și sunt necesare procesului de optimizare.
Pentru o funcționare corectă a programului, toate datele de intrare sunt necesare.
Impulsul specific, ISP : valoarea impulsului specific a motorului satelitului folosit în transfer, în
secunde.
Propulsia maximă, Thrust : valoarea maximă a modulului propulsiei motorului satelitului folosit în
transfer, în Newton
Rata de degradare, Degrading rate : valoarea ratei de degradare a motorului satelitului folosit în
transfer, în procente. Rata de degradare folosită în mod frecvent este de 5%/an.
Masa, Mass : valoarea masei inițiale a satelitului folosit în transfer, în kilograme.
Timpul de transfer, Transfer Time : valoarea estimativă a duratei totale de timp în care are loc
transferul, în zile.
Masa de combustibil, Propellant : valoarea estimativă a masei de combustibil folosită pentru
îndeplinirea cu succes a transferului, în kilograme.
Violarea maximă a constrângerii, Max Violation : valoarea absolută maximă a violării
constrângerilor, folosită ca și criteriu de oprire a procesului de optimizare, în cazul în care această
valoare a fost atinsă.
Numarul maxim de iterații, Iterations : valoarea maximă a numărului de iterații în procesul de
optimizare. Acest număr este folosit pentru a nu lăsa programul sa ruleze la nesfârșit în cazul în care
niciun criteriu de oprire nu este atins.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
44
Pasul de ieșire, Output step : valoarea pasului de timp folosit în fișierul de ieșire al traiectoriei
optimizate, în secunde.
5.2.2 Controlul procesului
În meniul principal, pe lângă datele de intrare, găsim de asemenea și cele trei butoane de control ale procesului
de optimizare.
Optimize : Pornește procesul de optimizare folosind toate datele de intrare provenite de la
utilizator.
Stop : Opreste întreg procesul de optimizare, fără a salva datele interemediare obținute în timpul
procesului.
Plot : Pornește afișarea grafică a rezultatelor obținute în urma procesului de optimizare. Atenție :
pentru ca acest buton să meargă, cel putin un proces de optimizare trebuie sa fie finalizat anterior cu
mesajul „OPTIMIZATION COMPLETE!”
5.3 Consola de Execuție
Consola de execuție va afișa pe parcursul rulării rezultatele intermediare ale optimizării cât și rezultatele
finale la sfârșitul execuției procesului.
Figure 5 – Consola de execuție
Mesajele afișate în consola de execuție în urma unei rulări cu succes al procesului de optimizare, sunt cele
de mai jos și pot fi definite astfel:
READING USER INPUT: se citesc datele provenite de la utilizator.
PROCESSING INITIAL GUESS: se citește traiectoria de referință și se creează fișierul ce contine
arcurile traiectoriei
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
45
INITIAL GUESS VALUES: se afișează starea inițială și finală citite din traiectoria de referință
Iteration i : se afișează numărul iterației din procesul de optimizare
o Cost function: valoarea scalată a funcției cost
o Final time: durata transferului, în zile
o Final position: poziția finală, în km
o Final velocity: viteza finală, în km/s
o Final mass: masa finală
o Max constraint violation: valoarea absolută a violării constrângerilor
o Cost function change : schimbarea funcției cost față de iterația anterioară
o Number of constraints: numărul de constrângeri
o Number of affected arcs: numărul total de arcuri afectate
o Number of optimised variables : numărul de variabile optimizate
Propagating and writing the final trajectory: se scrie fișierul ce conține traiectoria finală
Auxiliary calculations have been performed: se calculează orbita finală în elemente kepleriene
INITIAL ORBIT: se afișeaza orbita inițială optimizată
o Initial SMA[km]
o Initial e [-]
o Initial i [deg]
o Initial AOP [deg]
o Initial RAAN [deg]
o Initial Geo longitude [deg]
o Initial mass [kg]
FINAL ORBIT: se afișează orbita finală optimizată
o Final SMA [km]
o Final ecc [-]
o Final i [deg]
o Final AOP [deg]
o Final raan[deg]
o Final Geo longitude [deg]
o Final mass [kg]
OPTIMAL SOLUTION FOUND : se afișează valorile găsite ca fiind optime
o Total propellant used [kg] : masa de combustibil folosită, în kilograme
o Transfer duration [days]: durata transferului, în zile
OPTIMIZATION COMPLETE! : Acest mesaj reprezintă finalizarea cu succes al procesului de
optimizare. Dacă vor apărea erori pe parcursul procesului de optimizare, acestea vor fi afișate în
consolă, iar procesul va fi finalizat cu mesajul „OPTIMIZATION FAILED!”.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
46
5.4 Afișarea rezultatelor grafice
Cel de-al treilea tab al interfeței grafice este folosit pentru afișarea rezultatelor obținute în urma procesului
de optimizare. Acest tab v-a afișa rezultatele doar după apăsarea butonului „Plot” după ce procesul de
optimizare este complet. În cazul în care niciun proces de optimizare nu a fost realizat anterior cu succes,
rezultatele grafice vor afișa mesajul „NO DATA”.
Figure 6 – Afișarea rezultatelor grafice
Rezultatele disponibile pentru afișarea grafică sunt cele din tabelul următor:
TAB Rezultat afișat
Eccentricity Eccentricitatea vs Timp (zile)
Inclination Înclinația vs Timp (zile)
Mass Masa (kg) vs Timp (zile)
Hodograph Vizualizare 3D: Vx (km/s), Vy (km/s), Vz (km/s)
SMA SMA, Altitudinea la Perigeu și Apogeu (km) vs Timp (zile)
Trajectory Vizualizare 3D: X (km), Y (km), Z (km)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
47
Thrust Thrust XYZ și Modulul (kN)
Thrust Direction Direcția Thrustului vs Timp (zile)
Thrust Angles Unghiurile Thrustului vs Timp (zile)
Position Poziția XYZ și Raza (km) vs Timp (zile)
Velocity Viteza XYZ și Viteza Parcursului de zbor (km/s) vs Timp (zile)
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
48
6 Rezultate
6.1 Introducere
Având problema definită și drumul către rezolvarea acesteia bine pavat, este timpul să selectăm un caz
relevant pentru a vedea performanța întregului proces. Cazul selectat este derivat din referința 5, iar traiectoria
de referință a fost calculată cu ajutorul programului LOTOS, dezvoltat de Astos Solutions.
6.2 Transfer GTO-GEO
Cazul ales este următorul :
Transfer GTO-GEO : acesta este un caz crucial, datorită interesului tot mai crescut al
dezvoltatorilor și investitorilor din acest domeniu pentru a realiza astfel de transferuri folosind
propulsie electrică. Satelitul considerat în acest caz are o masă de 2126 kg, un thrust de 0.2 N și un
impuls specific de 2000 s.
Orbiă a (km) e i (°) (°) (°)
Inițială 24474 0.729 6 360 0
Obiectiv 42114 0.0001 0.01 180 180
6.3 Rezultatele transferului
După ce o traiectorie de referință bună a fost obținută de către programul folosit, în acest caz LOTOS,
optimizarea are nevoie doar de câteva iterații pentru a indeplini condițiile la limită cu o precizie foarte bună.
Transferul optim durează 255 de zile, iar combustibilul folosit este de 225 kg. Aceste rezultate corespund cu
cele exportate de LOTOS, după optimizarea traiectoriei de referință.
INITIAL ORBIT FINAL ORBIT
Initial SMA[km]: 2.447433e+04
Initial e [-]: 7.291625e-01
Initial i [deg]: 5.999550e+00
Initial AOP [deg]: 0
Initial RAAN [deg]: 3.600000e+02
Initial Geo longitude [deg]: 3.588866e+02
Initial mass [kg]: 2.125995e+03 Final SMA [km]: 4.211306e+04
Final ecc [-]: 9.665853e-04
Final i [deg]: 2.334825e-02
Final AOP [deg]: 2.805681e+02
Final raan[deg]: 9.781970e-02
Final Geo longitude [deg]: 1.800010e+02
Final mass [kg]: 1.900677e+03
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
49
OPTIMAL SOLUTION FOUND:
Total propellant used [kg]: 2.253175e+02
Transfer duration [days]: 2.557411e+02
Figure 7 Profilul optim al masei vs Timp
Figure 8 – Traiectoria satelitului, arată traiectoria urmată de către satelit în planul X-Y-Z. Această traiectorie
nu are nicio schimbare bruscă de plan, mai degrabă o scădere constantă a înclinației pentru a atinge orbita
ecuatorială dorită, deci în acest plan sunt văzute cele mai mari schimbări.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
50
Figure 8 – Traiectoria satelitului XYZ
Figure 9 – Poziție și Rază vs Timp
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
51
Figure 10 Hodograph XYZ
Figure 11 Viteză vs Timp
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
52
În Figure 12 – SMA, Raza la apogeu și perigeu, se poate observa cum raza apogeului este crescută peste
valoarea obiectiv, ajutând la reducerea combustibilului ce va fi folosit pentru schimbările de înclinație și
eccentricitate.
Figure 12 – SMA, Raza la apogeu și perigeu
Setul de reprezentări grafice de mai jos, arată mai multe detalii cu privire la evoluția elementelor orbitale.
Există două jumătăți ale traiectoriei ce pot fi diferențiate : prima jumătate în care raza apogeului crește, și a
doua jumătate în care scade până la valoarea inițială. Raza perigeului crește pe toată durata transferului.
Combinația dintre cele două produc oscilația din semi-axa mare în prima parte a transferului.
Inclinația și eccentricitatea scad monoton pe tot parcursul transferului. Totuși se poate observa cum înclinația
este redusă mai mult în prima parte a transferului, fiind în contrast cu eccentricitatea.
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
53
Figure 13 -Eccentricitate
Figure 14 – Înclinație
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
54
Cât despre legea de control ce produce transferul optim, componentele pe X, Y, Z sunt reprezentate în figurile
de mai jos. Aceste componente arată doar direcția vectorului propulsiei, și nu modulul acestuia, fiind
normalizate. Toate trei componentele arată schimbări rapide în comportament ale vectorului, coordonate cu
ciclul true anomaly, și un trend de schimbare mai mic în valorile extreme ale fiecărei revoluții. Analizând
valorile fiecărui component este ușor să realizăm că vectorul thrust formează un unghi foarte mic cu planul
orbital pe durata întregului transfer ceea ce este de așteptat din moment ce variațiile necesare schimbarii de
altitudine și eccentricitate sunt mai greu de obținut decât cele din schimbarea înclinației.
Toate rezultatele prezentate au fost obținute propagând doar cu integratorul Runge-Kutta, fără a include
procesul de mediere a ecuațiilor descris în capitolul 3.5.
Figure 15 Componenta X a legii de control în TOD
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
55
Figure 16 Componenta Y a legii de control în TOD
Figure 17 Componenta Z a legii de control TOD
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
56
7 Concluzii
7.1 Obiective atinse
Problema unui transfer orbital optimal a fost intens studiată. Până în acest moment, s-a căutat întelegerea
diferitelor metode pentru a rezolva astfel de probleme.
O metodă robustă, rapidă, și un program ușor de folosit a fost creat, care, odată ce primește o traiectorie
inițială bună, poate rezolva problema de masă minimă dintre două orbite, calculând legea de control ce
ghidează satelitul pe durata întregului transfer.
Acest obiectiv a fost atins cu o formulare analitică, folosind ecuațiile de mișcare Gauss pentru a obține
condițiile la limită pentru problema de control optimal, aplicând calculul variațiilor, și rezolvând problemă
cu condiții la limită aplicând metode numerice, un integrator Runge-Kutta pentru propagarea sistemului de
ecuții și un solver Newton-Rapshon pentru a obține soluția exactă a problemei.
Procesul a fost implementat într-un program de calcul scris în Matlab, capabil să rezolve problema în cateva
zeci de minute. De asemenea, include metode de control a acurateții și generează reprezentări grafice ale
rezultatelor.
Rezultatele nu doar că includ reprezentări ale rezultatelor ușor de citit pentru utilizatori, dar și o lege de
control care îi permite acestuia să producă un profil al atitudinii fezabil pentru satelit.
Cu ajutorul acestui software, odată ce problema de timp minim de transfer a fost rezolvată, este posibil să
studiem problema masei minime, pentru un transfer fixat cu un timp de transfer mai mare. Acest lucru aduce
o valoare foarte mare în faza de analiză a misiunii, pentru a calcula cu precizie trade-offul dintre timpul de
transfer si masa satelitului, și poate reduce costul misiunii spațiale.
Interesul pentru acest tip de transferuri a crescut exponențial în ultimii ani, și cu siguranță panta va păstra un
ritm ascendent și în următorii ani.
7.2 Implementări viitoare
Posibilitățile de îmbunătățire a aplicației sunt nenumărate, acestea nefiind incluse datorită limitării legate de
timp.
Câteva din viitoarele posibile implementări sunt următoarele :
Adăugarea unui modul de initializare pentru generarea unei traiectorii de referință
Adăugarea unor constrângeri de inegalitate pe apogeu si perigeu
Adăugarea unei constrângeri pentru longitudinea geografică finală
Adăugarea unei constrângeri pentru perioadele de eclipsă
UNIVERSITATEA POLITEHNICĂ BUCUREȘTI
FACULTATEA DE INGINERIE AEROSPAȚIALĂ
57
8 Bibliografie
1. „Semianalytic satellite theory”, D. A. Danielson, C. P. Sagovac, B. Neta and L. W. Early, Naval
Postgraduate School, 1994
2. „Automation of multi-revolution low-thrust transfer optimization via diffrential evolution”, Jose M.
Sanchez Perez, Andrea Campa, American Astronautical Society, 2014
3. „A survey of numerical methods for optimal control, V. Anil Rao, Advances in the Astronautical
Sciences, 2009
4. „Optimal Control: An Introduction to the Theory and Its Applications”, M. A. Athans and P. L. Falb.,
Dover Publications, New York, 2006.
5. „Optimal low-thrust transfer with constraints – generalization of averaging techniques”, Sophie Geroy
and Richard Epenoy, Acta Astronautica, 1997.
6. „Applied Optimal Control”, A. E. Bryson and Y.-C. Ho., Hemisphere Publishing, New York, 1975
7. “Minimizing the real functions of the icec’96 contest by differential evolution. Proceedings of IEEE
International Conference on Evolutionary Computation”, Rainer Storm and Kenneth Price, pages
842–844, 1996a.
8. „Introduction to numerical analysis. Texts in applied mathematics 12.” J. Stoer, R. Bulirsch, R.
Bartels, W. Gautschi, and C. Witzgall. Springer, 3rd ed edition, 2002
9. Orbit & Constellation Design and Management, JR Wertz p68-69
10. E.B. Shanks, Math.Comp.20 (1966), 21-38
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: Transfer orbital optimal folosind propulsie electrică [628667] (ID: 628667)
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.