Transfer orbital optimal folosind propulsie electrică [628668]

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

1

Transfer orbital optimal folosind propulsie electrică

Absolvent: [anonimizat]. Simioana Mădălin Vior el

Conducător : Prof. Dr. Ing. Teodor Viorel Chelaru

2016-2018

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

2

CUPRINS
1 Introducere ………………………….. ………………………….. ………………… 3
1.1 Structura lucrării ………………………….. ………………………….. ……………………. 4
1.2 Unelte folosite ………………………….. ………………………….. ……………………….. 4
2 Baze teo retice ………………………….. ………………………….. ……………… 5
2.1 Introducere ………………………….. ………………………….. ………………………….. .. 5
2.2 Coordonatele echinocțiale ………………………….. ………………………….. ……….. 5
2.3 Formularea problemei ………………………….. ………………………….. …………….. 5
2.4 Problema de control optimal. Metoda directă contra metodei indirecte. … 8
2.5 Aplicarea calculului variațional într -o problemă de co ntrol optimal ……. 10
2.5.1 Ecuațiile Euler -Lagrange ………………………….. ………………………….. ………………………….. ……………… 10
2.5.2 Timpul final neconstrâns ………………………….. ………………………….. ………………………….. ……………… 12
2.5.3 Obținerea problemei cu condiții la limită ………………………….. ………………………….. ……………………. 13
2.5.4 Normalizarea ………………………….. ………………………….. ………………………….. ………………………….. …. 15
2.5.5 Schimbarea variabilelor independente ………………………….. ………………………….. ……………………….. 16
2.5.6 Hamiltonianul ………………………….. ………………………….. ………………………….. ………………………….. … 16
2.5.7 Legea de control optimal ………………………….. ………………………….. ………………………….. ……………… 17
2.5.8 Ecuațiile diferențiale adjuncte ………………………….. ………………………….. ………………………….. ………. 18
3 Metode numerice ………………………….. ………………………….. ………. 20
3.1 Introducere ………………………….. ………………………….. ………………………….. 20
3.2 Algoritm ul Evoluției Diferențiale ………………………….. ……………………….. 20
3.3 Shooting Method ………………………….. ………………………….. ………………….. 22
3.4 Propagarea numerică ………………………….. ………………………….. …………….. 24
4 Bibliografie ………………………….. ………………………….. ……………….. 29

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

3

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 mecanici i 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 sa telitului. Datorită
acestor probleme, o activitate importantă o repr ezintă 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 ch imică, 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 u nui 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, d in 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 ch imic 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 prop ulsiei 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 tr ansfer 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 propulsi ei electrice în acest domeniu, și aplicarea acesteia în misiuni spațiale.
Cu toate acestea, revenim la problema anterioară și anume : mode lul de orbitare pentru propulsi e
chimic ă nu poate fi aplicat, din moment ce acele delta -V-uri instantanee aplicate sate litului 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ă prod use 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, me todele 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 INGINE RIE AEROSPAȚIALĂ

4

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 c u 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 implmen tarea 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 INGINE RIE AEROSPAȚIALĂ

5

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 co rp 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. Pen tru 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 cos2aa
he
ke
ip
iq
Lv

=
= + 
= + 
= 
= 
= + +  (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= (
22limih pq− + =  ), 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ția le directe. [1]
2.3 Formularea problemei
Pe scurt, problema considerată luată în considerare este : obținerea vectorului optim de tracțiune în
timp,
()ut
, care va permite satelitului sa căl ătorească de la o poziție inițială
0()xt
, la o poziție finală
diferită
()fxt
în cel mai scurt timp posibil, adică minimizarea
ft . Tracțiunea disponibilă a satelitului este
foarte mică în fi ecare moment de timp al transferului. Tracțiunea este considerată ca fiind mică atunci când
raportată cu atracția gr avitaț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 perturbații în afară de tracțiune, cu atât problema devine mai dificil de rezolvat. Pentru cel mai

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

6

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ă fo rma Gauss a variației element elor echinocțiale în timp , când subiectul de interes
este tracțiunea, este următoarea [2]:

( , )dx uF x Ldt m=
(2.2)

01( , ) ( , ).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
22
22
22
22200
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 LBB
aA D k L hk L L h k h LF x L h p L q LD BB
Lpq
Lpq


+ − − + + − − −
+ + − + +  = −

 ++

 ++
(2.4)

2
0 33( , )Dg x LaA=
(2.5)

1( , ) 0 0 ( sin cos )aAg x L p L q LD=−

(2.6)
Unde, A, B și D sunt definiți astfel :

22
221
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
32/km s . [2]

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.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

7

Longitudinea reala, L, este separată de celel alte element e 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 perpendic ulară 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
zuu
uu
uu

=
=
=
(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 /ms ș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 c omplet 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
traiector ie, 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 :

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

8

01
00
00
00min
( , )
( , ) ( , )
||
( ) ( )
( ) ( )
( ) ( )P SP
ff
f
f
fmizeazatf
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 liber=
=+
=−
==
=
=
(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âng erile. În cazul nostru,
f Jt= ș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 u neltelor matematice capabile să rezolve problema de
control optimal. În următorul subcapitol vom det alia principalele abordări asupra problemei.

2.4 Probl ema de control optimal. Metoda directă contra metodei indirecte.
Din punct de vedere matematic, problemele de control optimal sunt formulate astfel [ 3]: determina ți
starea,
()xt
, vector real de dimensiune n, controlul,
()ut
, vector real de dimensiune m, unde valoarea
inițială
0t și valoarea finală
ft a unei variabile reale independente
0,f t t t care optimizează indexul de
performanță :

 
000( ), , ( ), ( ), ( ),ft
ff
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,

  min max ( ), ( ), C C x t u t t C 
(2.13)

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

9

ș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
ft= pentru o problema de timp de transfer minim.
 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 /vr   = − .
Constrângerile traiectoriei C, reprezintă condiții pl asate î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 vizibilitat e.
Î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 înce rcă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 ind irecte 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 toa tă 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 pertu batoare, 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ă acest e 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 to ate 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
detali ată a teoriei calculul variațiilor.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

10

2.5 Aplicarea calculului variațional î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 a sociem sistemul dinamic indexului de performanță J, obținând :

     
000( ), , ( ), ( ), ( ), ( ) ( ), ( ),ft
ff
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
()ut
pentru timpi stabiliți
ft și
0t :

000( ) ( ) .f
ft
t ttHHJ 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 com ponentă 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 :

.Hf
x x x    = − = − −
(2.19)

unde
/fx
este o matrice de dimensiune
nn , 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
0tt= , 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()xt
produce o extra condiție pentru co -stările n -i:

0( ) 0t =
(2.20)
Î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

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

11

(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
()fxt
trebuie să îndeplinească condiția:

0
ftt x=−=
(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
u
pentru minimizarea lui J:

0 0, ,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 fo losind 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
HHHtu    
    

= +  +  + 
= +  +   
=+
(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
0 xf−=
. Atfel ecuația ce trebuie îndeplinită ajunge la forma:

 
022
2 2
2
2 22
211022f
ft
t ttHH
x x x uJ x x x u dtu x HH
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 INGINE RIE AEROSPAȚIALĂ

12

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:

()()()()
000f
ft
f f f
t ttffJ dt d x 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
()fxt

dispar, urmărind teoria detaliată în secțiunea anterioară. Astfel derivata indicelui de performanță devine:

0*f
ft
f
tt 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( ) 0xt =
, 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
() ifxt ;
0,…,iq= sunt definiți, iar
() jfxt ;
1 ,…, j q n=+ sunt liberi. Atunci, este
posibil să aflăm influența
u
asupra
() ifxt utilizând conceptul funcțiilor adjunct , obținând:

()()()
0*f
ft
i
i f i f tt
tfx t f dt t u dtu  ==+  (2.29)
unde
()() 1,
0,i
jfijtij==
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
tt tfdJ dx t f f dt udtt u u          = + = +  + + + + +   
(2.30)
Alegând:

()
1 *
fj
f i i
ttdt k f ft== − +  + +
(2.31)

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

13

și,

( )( ) ( )
2ji
ifukuu   = − + + 
(2.32)
unde k1 si k2 sunt constante pozitive, și introducând aceste expreii in 2.30 rezultă:

() ( )
022
( ) ( ) ( )
12 *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
i care satisfac ecuația 2.28, ecuațiile 2.31, 2.32 sunt
introduse in 2.28, și după câteva calcule rezultă:

() ( )
001
( ) ( ) ( ) ( ) ( ) 11
22ff
fftt
j j i j j
i j i itt ttttkk 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
ii
ttfft=+  + + =
(2.35)

( ) ()( ) ( )
0 0, ,ji
ifft t tuu  + +  =  
(2.36)
Astfel, impunând ecuațiile 2.35 si 2.36 ne permit să obținem o soluție staționară ce va satisfac e ș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 INGINE RIE AEROSPAȚIALĂ

14

01
00
00
00min
( , )
( , ) ( , )
||
( ) ( )
( ) ( )
( ) ( )P SP
ff
f
f
fmizeazatf
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 liber=
=+
=−
==
=
=
(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 Hamilton iene 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
x si
.
5. Rescrierea H introducând rezultatul anterior, și scriindu -l dor in funcție de
x si
.
6. Derivarea parțială a lui H în funcție de
x
și obținerea ecuațiilor diferențiale pentru

ca în e cuația
2.19.
7. Medierea ecuațiilor diferențiale.
8. Rezolvarea sistemului pentru
x si
, 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 med ierea 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 altfe l
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 INGINE RIE AEROSPAȚIALĂ

15

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 ffx mux 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
fa este folosită ca lungime de referință. Masa inițială
0m adimensionalizeaza 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
00
32
max max
0 1 0 1 3
00
max max
33
001
1ff
ff
f f f
f
f P sp f P spa a u uu dx uFFdt a a m m m m
a a a u uu 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
0fau
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:

00
1
0 max
23
max 01
f f P sp f P spdx u uFFd m m
gg dL ugdm
mu dmuud a u a m g I a g I

  
 ==
= + =
= − = − = −
(2.41)
Chiar dacă variabilele sunt adimesionale, barele nu au fost folosite pentru a păstra ecu ațiile curate. Din
acest moment acest criteriu va fi folosit. Al doilea termen al ecuației longitudinii este mult mai mic decat

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

16

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
00
00d
dL g
dx d u dx uFFd dL m g dL g m
dm d dmuud 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
10L dx uFd g m
L dmudg
L d
dg
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ță
f Jt= , și
0= :

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

17

11
1 1 1
0 0 0
1
0x m L
xm
xmdL dx dm dHd d d d
L L L uFug m g g
L uFugm

      
  
   =  + + +
=   − +
=   − +
(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
Conform calculului variațional, 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.

11
000
1xm
x
mLL HuFu g m g u
uFum 
=  − =
=
(2.46)
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ă:
( )
( )( )( )222
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 12
2 2 2 2 2 2 2 2 2
1 2 3 1 2 3 1 2 3,,
2 22 111,,222uuuu
uu
u u u u u u u u u
u u u
u uu
u u u u u u u u u
u
u


  ++
=
+ + + + + +=


=
 + + + + + +
=

cum apare in ecuația 2.46.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

18

Deoarece
u
este 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.47)
În ce privește modulul tracțiunii, este ușor de înțeles că, în cazul unui timp minim de transfer, tracțiunea
maximă este necesară în fiecare punct al traiectoriei. Din acest motiv
1optu= .
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
xm
x
x
m
x
xmF LHFgm F
F L
gm F
LFgm


   


  

     =   − +
 
= − + 
=  − +
(2.48)
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.

( )( )//
10 1
2
00, 1,…,5
1ix
i
xx
xm
xd Hidx
FF Lg LFg m g m F 


   
= − =
= −  − + − 
(2.49)

1
2
0m
xd L HFd m g m = − = 
(2.50)

1
1 0 111 L
xmd HHFd L g m L    = − = −  − + = −
(2.51)

0d H
d 
 = − = (2.52)
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ă.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

19

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
( )0 ms = . Adițional, orbita dorită este de asemenea definită,
( )1 , 1,…,5ix s i== . Cele
3 con diții ramase sunt obținute pentru adjuncte. Masa este liberă la extrema finală, deci, conform
calculului variațiilor, valorile adjunctei
( )1 / 0msm   = = = deoarece pentru problema noastră
ft= .
În final,
1L este liberă la ambele exteme ale traiectoriei, deci:
11(0) (1) 0LL == .
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 INGINE RIE AEROSPAȚIALĂ

20

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 traiec torii, 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 dec at 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 valoril e 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, am estecarea 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 Rung e-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ă nu merica de integrare Euler -Lagrange este
implementată pentru a efectua medierea.
În ultimul rând, un instrument este necesar pentru a rezolva pro bleme 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 ev oluției diferențiale efectuează câteva

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

21

evaluări în 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 m enț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 cel e 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.53)
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.55)
Cu
( )
1
,1
,
,,…, 1
, 1,G
ijG
ijG
ijv pentru j n n L
u
x j D+
+ = + −= 
(2.56)
n și L fiind scalari generați aleatoriu, urmând întotdeauna regula 0 < n < L < D.
Asta înseamnă ca, pentru fiecare parametru al vectorului
ix
, 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 a lese 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:
11GG
iixu++=
. În caz contrar, vectorul anterior este reținut:
1GG
iixx+=
.
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:
( )15222
1f obj
i i m L
iF x x 
== − + +
(2.57)
Această funcție obiectiv este validă ca orice altă funcție continua a cărui minim este obținut doar atunci
când orbita fina lă este egală cu orbita obiectiv, și variază ușo r în apropierea acestui minim.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

22

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 folositoa re 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.58)
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.59)
Și
21n n n=− condiții la limită pentru punctul final
2x
( )2 2 2 , 0, 1,…,kB x y k n ==
(2.60)
Procesul de găsire a soluției prin această metodă este următorul:
1. Estimarea tuturor variabilelor de stare la starea inițială
1()yx
. 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ă impl ementată 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
21n 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
()1yx
.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

23

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
()1yx
sunt obținute
variabilele finale
()2yx
.
În continuare , un nou vector de dimensiune
2,nd
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 destu l de
simplă și anume, folosirea condițiilor la limită la punctul final,
( )2 2 2 , , 1,…,kkd 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
v
astfel încât
d
sa fie zero. Acest lucru
este obținut prin calculul soluției:

A v d = −
(2.61)

Ș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
ij
jdAv
=
(2.62)
Valori ce nu pot fi obținute în mod analitic. În loc, la fiecare iterație aceste deriva te parțiale sunt estimate
folosind diferențe finite.
( )( )22 11,…, ,…, ,…, ,…,i j j n i j ni
jjd v v v v d v v v d
vv
+  −
=
(2.63)
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, c onvergenț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

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

24

scalară
dd
descrește strict monoton în timpul iterațiilor. Asta înseamnă că după fiecare pas al iterației
valoarea
dd
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
propa garea 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 i ntegrării. Sistemul de ecuații diferențiale ordinare,
ramâne deci același:
()(),i
idy xf x ydx=
(2.64)
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 m etode 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 expansiun e 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ăr at 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 INGINE RIE AEROSPAȚIALĂ

25

( )()2
1 ,n n n ny y hf x y O h+= + +
(2.65)
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 ord inul
4:
( )
( )
()1
1
2
2
3
43
5 3 1 2 4
1,
,22
,22
,
6 3 3 6nn
nn
nn
nn
nnk hf x y
k hk hf x y
k hk hf x y
k hf x h y k
k k k ky y O h+=
= + + 

= + + 

= + +
= + + + + +
(2.66)
Formule Runge -Kutta de ordin superior sunt rareori folosi te, de obicei deoarece pentru m ordin sunt
necesare mai mult de m eval uă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
2225!
225!yxy x h y h O h
yxy x h y h O h
+ − = + 

+ − = + 
(2.67)
Erori de trunchiere pot fi estimate prin diferența celor două propagări:

21yy = −
(2.68)
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 INGINE RIE AEROSPAȚIALĂ

26

( ) ()6
2 215y x h y O h+ = + +
(2.69)
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
ihinh ==
(2.70)
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.71)
Î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.72)
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 ma i 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 for mule de ordin superior cu același număr de evaluări.

UNIVERSITATEA POLITEHNICĂ BUCURE ȘTI
FACULTATEA DE INGINE RIE AEROSPAȚIALĂ

27

Formulele create din această filozofie sunt denumite formulele Gaussian patratice. Forma lor generală fiind
umătoarea:
()() ()
1b N
ii
i aW x f x dx w f x
= 
(2.73)
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 calc ul 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.74)
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
11
11
11
1
111
i i i i
i i i
i i i ipx
xp p
p x x p x
pp
xp p xp p
p x x p x x p x
p p p p−
−−=

=−

   
= − − −   
      
(2.75)

Poate fi dovedit [8] că polinomialul
()npx 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 util izatorul 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
11
110
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.76)
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 INGINE RIE AEROSPAȚIALĂ

28

()()1
, 1,…,ni
i n idp xw x i ndx−== 
 (2.77)
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.78)
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 secula re. 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 INGINE RIE AEROSPAȚIALĂ

29

4 Bibliografie
1. „Semianalytic satellite theory”, D. A. Danielson, C. P. S agovac, 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 Socie ty, 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 mathema tics 12. ” J. Stoer, R. Bulirsch, R.
Bartels, W. Gautschi, and C. Witzgall. Springer, 3rd ed edition, 2002