KTH - Institutionen för matematik
Kandidatexamensarbete
Implementering av en numerisk lösare
för optimal styrning
Mikael Andblom - mandblom@kth.se
Mai Nguyen - maing@kth.se
handledare Mattias Sandberg
Innehåll
1 Inledning 2
1.1 Rapportens struktur . . . 2
1.2 Formulering av ett exempel på optimalt styrproblem . . . 2
1.3 Raketlandning . . . 3
1.4 Mer generellt optimalt styrproblem . . . 3
1.5 Vårt arbete . . . 3
2 Ett testproblem 3 2.1 Formulering samt gissning . . . 3
2.2 Lagrangemultiplikator . . . 4
3 Pontryagins minimeringsprincip 6 4 Regularisering 7 4.1 Ett regulariserat testproblem . . . 7
4.2 Regulariserad Hamiltonian . . . 9
5 Diskretisering 10 5.1 Diskretisering av optimalt styrproblem i D dimensioner . . . 10
6 Användarens indata 11
7 Numerisk derivering 12
8 Vår implementation av den numeriska lösaren 13
9 Raketuppskjutning 14
1 Inledning
1.1 Rapportens struktur
I den här rapporten kommer vi introducera området optimal styrning samt visa hur man löser denna typ av problem både analytiskt och numeriskt. Vi kommer inleda med att visa hur ett optimalt styrproblem kan se ut i praktiken och sedan utvidga det till något mer allmänt. Varje nytt begrepp kommer motiveras med ett exempel, både av pedagogiska skäl men även för att senare kunna visa att lösarens svar stämmer överens med eventuella förväntningar.
1.2 Formulering av ett exempel på optimalt styrproblem
Tänk dig att en raket ska landa på jorden efter en resa i rymden med landningstiden T, raketens begynnelsehastighet v0 samt höjden h givet. Då är det rimligt att man vill
minimera bränsleförbrukningen för landningen. Om vi låter α : [0, T ] → [0, αmax] vara
den momentana bränsleförbrukningen [l/s] så fås den totala mängden förbrukat bränsle som
Z T
0
α(t) dt.
Om vi nu försöker minimera detta får vi lösningen α(t) ≡ 0. För visso resulterar detta i en väldigt låg förbrukning (ingen alls!) men däremot kommer raketen antagligen landa med en oönskad hög hastighet. Vi måste därmed straa landningshastigheten. Låt x1(t)
vara den momentana hastigheten. Då vill vi samtidigt minimera x2
1(T ), där kvadraten
säkerställer att hastigheter nära 0 belönas och inte −∞.
Hur vet vi att denna lösning medför att raketen landar på marknivå? För detta kan vi införa raketens altitud x2(t). På samma sätt som för hastigheten vill vi då minimera
x22(T ). Vi kan också tänka oss att minimera Cx22(T ) där konstanten C > 0 kan justeras så att landningsnivån blir acceptabel och på motsvarande sätt för hastigheten. Ett sätt att ta hänsyn till samtliga bidrag är att minimera summan, dvs
Z T 0
α(t) dt + Bx21(T ) + Cx22(T ).
Det sista vi måste göra nu innan vår modell är färdig är att koppla ihop styrsignalen α(t) med tillstånden x1(t) och x2(t). Detta gör vi med Newtons andra lag som säger
m ˙x1(t) = Ftot.
Kraften skulle kunna modelleras som
Ftot = mg − Dα(t),
där vi antagit att raketens lyftkraft är proportionell mot bränsleförbrukningen med pro-portionalitetskonstant D och dessutom försummat andra krafter så som luftmotståndet samt antagit konstant massa för att inte krångla till inledningen i onödan. Altituden x2
1.3 Raketlandning
Problemet som beskrevs i föregående sektion kan nu sammanfattas som att minimera Z T
0
α(t) dt + x21(T ) + Cx22(T )
med avseende på styrsignalen α : [0, T ] → [0, αmax] tillsammans med
begynnelsevärdes-problemen ( ˙ x1(t) = g − Dmα(t), x1(0) = v0 ˙ x2(t) = −x1(t), x2(0) = h,
där minustecknet i rad två kommer från att vi de nierat hastigheten positiv nedåt medan altituden som positiv uppåt. Ett optimalt styrproblem är alltså på formen att minimera en funktional där bivillkoren består av dierentialekvationer med begynnelsevärden.
1.4 Mer generellt optimalt styrproblem
De niera funktionalen F enligt F [x, α] =
Z T
0
h(α(t), x(t)) dt + g(x(T )). (1) Det optimala styrproblemet de nieras då som att minimera funktionalen F med avseende på styrsignalen α : [0, T ] → B där B ⊂ RDoch D är dimensionen för α, tillsammans med
begynnelsevärdesproblemet ˙
x(t) = m(x, α), x(0) = x0.
Notera att funktionerna h, g, m, konstanten T och vektorn x0är givet i problemet.
Uppgif-ten är alltså att nna α(t) och därmed blir även tillstånden x(t) via m och x0 bestämda.
1.5 Vårt arbete
Resultatet av vårt arbete är en numerisk lösare till optimala styrproblem. Hur koden fungerar och vilka indata som behövs kommer vi återkomma till senare i denna rapport. Först behövs en mycket användbar princip som kommer att motiveras med ett problem i nästa avsnitt. Principen kommer sedan formuleras mer allmänt i avsnitt 3.
2 Ett testproblem
2.1 Formulering samt gissning
För att ha ett relativt enkelt problem att bolla med samt testköra senare i vår kod har vi formulerat ett endimensionellt testproblem. De niera funktionalen F enligt
F [x, α] = Z 1
0
x2(t) + α2(t) dt + (x(T ) − 1)2. Vi vill nu minimera F givet
˙x(t) = α(t), x(0) = 1
1. slutvärde nära 1
2. funktionsvärden nära 0
3. derivata vars värden är nära 0
Om vi fokuserar på punkt 1 och 2 kan vi välja en parabel som når x-axeln på halva intervallets längd, dvs x(t) = 4(x − 0.5)2. För senare referens döper vi denna funktion till
x1(t). Sätter vi in detta i F fås
F [x1(t), ˙x1(t)] = F1 = 5.53.
Om vi istället prioriterar punkt 1 och 3 kan vi helt enkelt bara välja den konstanta funktionen x2(t) ≡ 1. Vi får då
F2 = 1.
Den optimala lösningen bör däremot ta hänsyn till alla tre punkter, men hur då?
2.2 Lagrangemultiplikator
Vill man minimera en funktion f(x, y) begränsad av bivillkoret g(x, y) = 0 kan man använda sig av metoden Lagrangemultiplikator. Denna metod använder sig av faktumet att i en extrempunkt måste gradienten av f och g vara parallella. Metoden sammanfattar detta som att extrempunkter antas då
∇x,y,λ L = 0.
Här är L Lagrangianen som de nieras genom
L(x, y, λ) = f (x, y) − λg(x, y). (2) I vårt fall vill vi inte minimera en funktion f utan en funktional F där bivillkoret istället blir en dierentialekvation. Det visar sig att denna metod blir väldigt användbar även för oss om vi gör följande översättningar
f (x, y) ⇒ F [x, α] x, y, λ ⇒ x(t), α(t), λ(t) g(x, y) ⇒ g(t) = ˙x(t) − m(x(t), α(t)).
Eftersom Lagrangianens argument nu är funktioner så måste vi se till att högerledet i (2) är en funktional. Genom (1) är redan F en funktional men däremot är λ(t)g(t) inte det. Multiplikationen väljs då att översättas till en L2 produkt och vi får
L[x(t), α(t), λ(t)] = F [x(t), α(t)] − Z T 0 λ(t)g(t) dt. Sätter vi in testproblemet fås L = Z 1 0 x2+ α2− λ( ˙x − α) dt + (x(1) − 1)2.
De skalära derivatorna på Lagrangianen översätts nu till riktningsderivator, även kallad Gâteaux derivata, längs v(t) enligt
L0x ⇒ lim h→0 L(x + hv, α, λ) − L(x, α, λ) h = lim h→0 1 h( Z 1 0 2hvx − λh ˙v dt + 2hv(1)(x(1) − 1) + O(h2)) = Z 1 0 2vx + ˙λv dt − λ(1)v(1) + 2v(1)(x(1) − 1). (3)
Där vi i sista steget har integrerat partiellt och använt oss av v(0) = 0 eftersom startpunk-ten för x är given och vi behöver därmed inte variera oss från denna. För den optimala lösningen x ska (3) vara noll. Detta innebär att integranden är identiskt noll samt att punktevalueringen är noll. Eftersom detta ska gälla oberoende av valet av v så har vi att
2x + ˙λ = 0, 2(x(1) − 1) − λ(1) = 0. Gör vi samma beräkning för α fås L0α ⇒ lim h→0 L(x, α + hv, λ) − L(x, α, λ) h = lim h→0 1 h( Z 1 0 2hvα + λhv dt + O(h2)) = Z 1 0 (2α + λ)v dt. (4)
Eftersom detta ska vara noll för alla v måste vi ha 2α + λ = 0. Slutligen gör vi samma beräkning för λ
L0λ ⇒ lim h→0 L(x, α, λ + hv) − L(x, α, λ) h = lim h→0 1 h( Z 1 0 −hv( ˙x − α) dt) = Z 1 0 (α − ˙x)v dt. (5)
Detta ger oss
α − ˙x = 0,
vilket är precis begynnelsevärdesproblemet. Sammanfattningsvis har vi ˙x = α, x(0) = 1 ˙λ = −2x, λ(1) = 2(x(1) − 1) λ = −2α. (6) Detta system är lösbart eftersom vi har tre okända funktioner men också tre linjärt obero-ende ekvationer, dessutom har varje första ordningens dierentialekvation ett tillhörande
randvillkor vilket gör att den i sig är bestämd. Detta går med papper och penna samt tålamod att lösa analytiskt. Lösningen x∗ fås då som
x∗(t) = e t−1 2 + (1 − 1 2e)e −t . Insättning i F ger nu F∗ = 1 2(3 + 1 e2 − 4 e) ≈ 0.832.
Figur 1: Gissningarna tillsammans med lösningen samt funktionalens värde Notera att ekvationerna (6) kan skrivas som
(
˙x = Hλ0 x(0) = 1 ˙λ = −H0
x λ(1) = ∂xT(xT − 1)2.
(7) Där xT = x(T ) = x(1) och H de nieras som
H(x, λ) = x2−λ
2
4 .
3 Pontryagins minimeringsprincip
Vi har precis observerat ett viktigt samband, nämligen Pontryagins minimeringsprincip. Principen säger att lösningen till ett optimalt styrproblem måste uppfylla
( ˙ x(t) = ∇λH(λ(t), x(t)) x(0) = x0 ˙ λ(t) = −∇xH(λ(t), x(t)) λ(T ) = ∇xg(x(T )). (8) Funktionen H kallas för Hamiltonianen och de nieras som
H(x, λ) = min
α∈B(h(x, α) + λ · m(x, α)).
Där B ⊂ RD. Kravet på den optimala styrsignalen är alltså att den ska lösa detta system
av första ordningens dierentialekvationer, vilket fungerar väldigt bra att lösas numeriskt. Systemet kallas ibland för det Hamiltonska systemet. Det är dock viktigt att notera att villkoret ovan enbart är ett nödvändigt villkor men ej tillräckligt. Det är alltså inte garan-terat att lösningen till det Hamiltonska systemet verkligen är minimum för motsvarande optimala styrproblem. Det kan även uppstå problem vid bestämningen av Hamiltonianen, vilket exempli eras i nästa avsnitt.
4 Regularisering
4.1 Ett regulariserat testproblem
För att illustrera ett vanligt problem vid bestämningen av Hamiltonianen inför vi ett nytt testproblem enligt följande. Låt
F [x(t), α(t)] = Z 1
0
x2+ α dt + (x(1) − 1)2.
Där α(t) ∈ [0, 2] för alla 0 ≤ t ≤ 1. Uppgiften är att minimera F med villkoret ˙x = α, x(0) = 0. Hamiltonianen fås nu som H = min α (x 2 + α + λα) = x2+ min α (α(1 + λ)) = x 2 + Q(λ). Där vi har infört Q(λ) enligt
Q(λ) = (
2(1 + λ) för λ < −1 0 annars.
Detta fås eftersom när 1 + λ är negativt då vill α anta sitt högsta värde medan när 1 + λ är positiv då vill α anta sitt minsta värde, dvs enbart randpunkterna kan antas för α. Att lösa det Hamiltonska systemet blir nu besvärligt på grund av att Q inte är slät. Systemet innehåller derivatan av Q som fås utav
Q0(λ) = (
2 för λ < −1 0 annars.
För att göra denna funktion slät (regularisera) inför vi (se Figur 2 nedan) Q0δ(λ) = 1 − tanhλ + 1
Figur 2: Q0 och två av dess regulariseringar.
Som syns i guren ovan blir regulariseringen bättre ju mindre vi gör δ. I den numeriska lösaren kommer man att iterativt lösa problemet med ett gradvis avtagande δ.
Om vi vill lösa testproblemet analytiskt gör vi antagandet att styrsignalen α(t) bara slås på en gång under hela intervallet, nämligen vid tidpunkten a ∈ [0, 1]. Eftersom vi har begynnelsevärdet x(0) = 0 och α(t) är derivatan som bara kan anta värdet 0 eller 2 får vi en graf enligt följande
Figur 3: Tillståndet x(t) där vi slår på styrningen α(t) vid tidpunkten a. Insättning i funktionalen F ger nu
F [x(t), α(t)] = f (a) = −4a
3
3 + 8a
2− 10a + 13
3 .
För att nna minimum måste vi testa randpunkterna samt extrempunkterna. Randpunk-terna ger
f (0) = 13
3 ≈ 4.33, f (1) = 1. De två extrempunkterna till denna funktion ges av
a = 2 ± r
3 2, men bara a = 2 −q3
2 ≈ 0.775ligger i det givna intervallet, insättning ger
f (2 − r
3
2) ≈ 0.768.
Tidpunkten då styrsignalen ska slås på är alltså ungefär för t = 0.775.
Figur 4: Tillståndet x(t) för några olika val av δ. Plotten visar att inte bara Hamiltonianen går mot den exakta motsvarigheten då δ går mot 0 utan även tillståndet x(t).
4.2 Regulariserad Hamiltonian
Vi ger i vår kod användaren möjligheten att mata in en regulariserad Hamiltonian enligt H → Hδ ; |Hδ(x, λ) − H(x, λ)| ≤ cδ ∀ x, λ.
Ett användbart trick när man försöker regularisera sin hamintonian är att använda sig utav tangens hyperbolicus tanh som har den trevliga egenskapen lim
x→±∞tanh x = ±1. Den
Hδ(x, λ) = x2+ (1 + λ)(1 − tanhλ + 1 δ ).
5 Diskretisering
I Matlab nns en mycket användbar funktion kallad fsolve som kan lösa ekvationssystem på formen F (y) = 0 där både F och y är vektorer. Metoden kan även ta emot en Jacobian för snabbare lösning. Det första vi måste göra för att kunna formulera optimala styrpro-blem på denna form är att gå över från kontinuerlig tid till diskret tid genom att dela upp det kontinuerliga tidsintervallet i ekvidistanta punkter enligt
t ∈ [0, T ] → t ∈ {t0, t1, t2, ..., tN}, tn= n∆t, n ∈ [0, N ], tN = T
x(tn) → xn
˙x(tn) →
xn+1− xn
∆t . Motsvarande samband gäller dessutom för λ.
5.1 Diskretisering av optimalt styrproblem i D dimensioner
Vi börjar med att införa kolonnvektorn xn ∈ RD vars element är samtliga tillstånd vid
tidpunkten tn. Motsvarande de nition görs även för vektorn λ.
Givet en Hamiltonian H(x, λ) så ska Pontryagins princip gälla för varje tidpunkt tn
(x n+1−xn ∆t = ∇λH(xn, λn) λn+1−λn ∆t = −∇xH(xn, λn). (9) Eftersom både x och λ har längd D så kommer alltså ekvation (8) bestå av 2D rader för varje n. Dessutom måste följande villkor uppfyllas
x0 = x00, λN = ∇xg(xN).
Här är x0
0 kolonnvektorn av alla initialvärden för x. Vi kan nu införa
F (y) = x0− x00 x1− x0− ∆t∇λH(x0, λ0) λ1 − λ0+ ∆t∇xH(x0, λ0) ... xN − xN −1− ∆t∇λH(xN −1, λN −1) λN − λN −1+ ∆t∇xH(xN −1, λN −1) λN − ∇xg(xN) , y = x0 λ0 ... xN λN . (10)
Observera att y är en vektor och inte en matris, ty xn och λn är kolonnvektorer. Rad n
av Jacobianen J beräknas nu som gradienten av rad n i F med avseende på y. Rad 1 i F ovan kommer då deriveras till enhetsmatrisen ID följt av nollor till höger. Rad 2 och 3
kommer deriveras till
följt av nollor till höger. Vi har infört en matris för andraderivatorna av Hamiltonianen enligt Hpp(xn, λn) = ∂λ1∂x1 · · · ∂λ1∂xD ∂λ1∂λ1 · · · ∂λ1∂λD ... ∂λD∂x1 ... −∂x1∂x1 ... ... −∂xD∂x1 −∂xD∂λD H(xn, λn). (11)
Minustecknet på den undre halvan av matrisen kommer från teckenskillnaden i Pontry-agins princip mellan derivatorna på x och λ. Nedan visas hur denna matris ser ut i två dimensioner Hpp(xn, λn) = ∂λ1x1 ∂λ1x2 ∂λ1λ1 ∂λ1λ2 ∂λ2x1 ∂λ2x2 ∂λ2λ1 ∂λ2λ2 −∂x1x1 −∂x1x2 −∂x1λ1 −∂x1λ2 −∂x2x1 −∂x2x2 −∂x2λ1 −∂x2λ2 H(xn, λn). (12)
Därefter kommer rad 4 och 5 deriveras på samma sätt förutom att Hpp ska evalueras i
(x1, λ1) samt så kommer motsvarande rader i J börja med en nollmatris av dimension
2D × 2D. För rad 6 och 7 kommer J börja med två stycken nollmatriser av samma storlek osv. Sista raden i F deriveras till
−gpp(xN), ID.
Där vi har infört en matris för andraderivatorna till g enligt
gpp(xn) = ∂x1∂x1 · · · ∂x1∂xD ... ∂xD∂x1 · · · ∂xD∂xD g(xn). (13) Sammanfattningsvis får vi J = ID 0 0 0 · · · 0 −I2D− ∆tHpp(x0, λ0) I2D 0 0 · · · 0 0 −I2D− ∆tHpp(x1, λ1) I2D 0 · · · 0 ... ... 0 · · · 0 0 −gpp(xN) ID . (14)
Matlabs metod fsolve har nu allt den behöver förutom en startgissning på y. Hur denna kan se ut beskrivs i avsnitt 6 Användarens indata.
6 Användarens indata
Användaren har en del valmöjligheter vid användningen av vår lösare. Följande inmatning behövs
• Hamiltonian, regulariserad eller ej. I det regulariserad fallet anger man ett önskat initialvärde på regulariseringsparametern δmax, förslagsvis något stort, ett undre
värde δmin samt ett startvärde på reduceringsfaktorn dq. Efter varje beräkning med
fsolve beräknas också ett nytt värde på regulariseringsparametern enligt δnew =
δold/dq. Beroende på om beräkningen är lyckad eller inte kommer dq automatiskt få
ett nytt värde. Har man en icke regulariserad Hamiltonian ska man sätta δmax = δmin
vilket då ser till att programmet enbart löser uppgiften en gång. • Sluttiden T
• Initialvärden för x(t), dvs x0
• Antalet tidsteg N. Större N motsvarar högre noggrannhet i svaret på bekostnad av beräkningstid
• Startgissning på λ(t) och x(t). Programmet är i nuläget kodat så att givet en start-gissning för x(t) beräknas en startstart-gissning för λ(t) baklänges via det Hamiltonska systemet. Detta går självklart att modi era om man hellre vill ange λ(t) eller både x(t) och λ(t). Många gånger krävs en bra startgissning för att få konvergens vilket kommer att illustreras med ett exempel senare i rapporten. Ibland har man turen att en trivial startgissning så som x(t) = 0 medför konvergens, så det kan vara värt att undersöka. Annars kan man behöva skaa sig en känsla för hur lösningen bör se ut och använda det som startgissning.
• Funktionen g(x) som motsvarar slutkostnaden i det optimala styrproblemet
• Om man kan och vill så nns möjligheten att mata in analytiska derivator för Ha-miltonaianen samt g annars kommer de approximeras med hjälp av centraldierens. Om de anges analytiskt får man en i många fall kortare beräkningstid. Anges deri-vatorna analytiskt behöver man för inte ange H eller g direkt.
7 Numerisk derivering
I många problem kan man få en Hamiltonian vars derivator man inte vill beräkna ana-lytiskt. Det kan vara så att den består av många termer eller av uttryck som helt enkelt är besvärliga att beräkna. Många funktionsuttryck blir också allt mer krångliga för var-je derivation. Då kan man spara tid och besvär genom att låta datorn beräkna dessa numeriskt.
Om man däremot har ett problem med en liknande Hamiltonian som ska lösas många gånger kan man spara tid på att beräkna derivatorna analytiskt. Ett exempel på ett sådant fall är om Hamiltonianen innehåller några obestämda konstanter man vill testa olika värden på mellan varje körning för att nna en önskad lösning.
Den numeriska metoden vi använder för numerisk derivering är centraldierens f0(x0) ≈
f (x0+ h) − f (x0− h)
2h .
Felet på centraldierensen är av ordning h2 medan avrundningsfelet från Matlab går som
h där ≈ 10
−15. Det totala felet är då proportionellt mot summan av dessa. Dvs
e ∝ h2+ h.
Derivering av detta uttryck ger att minimum fås för h = p3
2 ≈ 10
−5. Det värde på h som
vi använder oss utav är h = 10−5.
På motsvarande vis fås f00(x0) ≈ f (x0+ 2h) − 2f (x0) + f (x0 − 2h) 4h2 . Korsderivata i två dimensioner fås av fxy00 (x0, y0) ≈ f (x0+ h, y0+ h) − f (x0− h, y0+ h) − f (x0+ h, y0− h) + f (x0− h, y0− h) 4h2 .
8 Vår implementation av den numeriska lösaren
Lösaren vi har skrivit består av 8 ler enligt nedan.
• solver.m - Denna l initierar beräkningsprocessen. Det är även här man de nierar sluttiden T, initialvärdena x0, antalet tidssteg N, startgissning för x och λ, val av
numerisk eller analytisk derivering, regulariseringsparameterns maxvärde, minvärde och reduceringsfaktorns startvärde dqstart. I denna l nns ett metodanrop till len
regsolve.m som utför själva lösandet. Vill man plotta lösarens svar är det smidigt att skriva den koden i denna l efter att regsolve har bestämt ett y.
• regsolve.m - I denna l anropas Matlab metoden fsolve i en regulariseringsprocess. Processen går ut på att lösa uppgiften för en given regulariseringsparameter δ, om lösningen konvergerade använder man svaret som startgissning för beräkningen med ett lägre δ. Detta görs i en while-loop tills δ ≤ δmin. Om fsolve returnerar en exit ag
motsvarande en lyckad lösning så behålls det aktuella värdet på reduceringsfaktorn dq annars ökas den. Detta görs eftersom om δ sänks för fort kommer inte lösningen till det tidigare δ som startgissning kunna få lösaren att konvergera. Koden återgår till det senaste δ som medförde konvergens och sänker den med ett nytt dq som är lägre. Om dq sänks under en undre gräns dqmin kommer lösningsprocessen avslutas
och ett y motsvarande det lägsta δ som gav konvergens returneras.
• F.m - Detta är funktionen som bygger upp vektorn F som fsolve ska nna lösning till. Den returnerar även Jacobianen J. Funktionen använder sig även av två andra funktioner g.m och Hp.m.
• Hp.m - Denna l returnerar förstaderivatorna till H evaluerad i önskad punkt. Den är uppdelad i en del som kan beräkna dem numeriskt samt en del där användaren kan skriva in de analytiskt bestämda derivatorna om så önskas. Uppdelningen sker via en if-sats som tar emot argumentet numDer som antingen är true om man vill använda sig av numerisk derivering, annars false. Detta bestäms som tidigare nämnts i kör len.
• g.m - Här beräknas första och andra derivatorna till g. Har man beräknat deriva-torna analytiskt ska man skriva in dem i denna l annars nns kod för numerisk derivering i denna l.
• J.m - Filen som returnerar Jacobianen till F.m. Den behöver lerna g.m och Hpp.m för att bestämma matrisen.
• Hpp.m - Här bestäms samtliga andraderivator till Hamiltonianen H, med samma struktur som för Hp. De numeriska derivatorna beräknas från funktionen H men man hade lika gärna kunnat koda så de bestäms från funktionen Hp istället. • H.m - Som tidigare nämnts anropas denna funktion från funktionerna Hp och Hpp
om numerisk derivation är önskad. Här fyller man alltså in sin Hamiltonian om man vill ha numerisk derivering.
Observera att solver.m och g.m är problemberoende och alltid måste modi eras. Har man beräknat nödvändiga derivator analytiskt behöver man ange dessa i lerna g.m, hp.m och hpp.m annars modi erar man bara H.m. Filerna F.m och J.m behöver aldrig modi eras. Eftersom det bara nns en variabel, numDer, som har koll på vilken derivering som önskas medför detta att användaren måste antingen ange samtliga derivator analytiskt eller inga alls.
9 Raketuppskjutning
Vi ska nu betrakta ett problem som visar att optimala styrproblem inte alltid är så lätta att lösa även med den kod vi skrivit. Tänk dig följande problem. En raket väger 15 kg varav 5 kg består av bränslet. Denna raket ska skjutas upp så högt upp som möjligt. Vilken styrsignal α(t) : [0, T ] → [0, 1] kg/s för bränsleförbrukningen ska man välja?
Rörelseekvationen för raketen skrivs som
m ˙v = Fmotor− Fluf t− Fg, v(0) = 0 m/s.
Lyftkraften från motorn modelleras som Fmotor = cα där c = 1000 m/s. Luftmotståndet
modelleras som Fluf t = 12ρAC|v|v = D|v|v där tvärsnittsarean A = 0.5 m2,
luftmot-ståndskoe cienten C = 0.7 och luftens densitet ρ = 1.3 kg/m3. Detta ger oss värdet
D = 0.23 kg/m. Anledningen till varför vi har valt att skriva |v|v istället för v2 är för att
luftmotståndet är antiparallellt mot hastigheten. Vidare har vi Fg = mg där m är raketens
massa och gravitationskonstanten g = 9.82 m/s2. Raketens massa styrs av ˙m = −α.
Ett sätt att formulera detta som ett optimalt styrproblem är min
α(t)∈[0,1]
Z 25
0
−v+dt + K(m(25) − 10)2.
Där v+= max(0, v). Integralen av v+ ger då höjden medan den andra termen ser till att
allt bränsle är förbrukat (och inte mer!) vid sluttiden för tillräckligt stort värde på K. Argumentationen kring valet av sluttid på 25 s kommer senare i detta avsnitt.
Låt oss nu bestämma tillhörande Hamiltonian. De två tillståndsfunktionerna vi har är v(t) och m(t) och deras tillhörande dierentialekvationer ges av
(
˙v = m1(cα − Dv2) − g ˙
m = −α. (15)
För att nna den optimala styrsignalen behöver vi enbart korrekt modell då raketen är på väg uppåt. Det är då ok att använda oss av v2 istället för |v|v här. Däremot är denna
skillnad viktig vid den gra ska representationen av resultatet.
Låt nu dualfunktionerna λ1 och λ2 tillhöra tillstånden v respektive m. Per de nition
H(v, m, λ1, λ2) = min α (−v++ λ1 m(cα − Dv 2) − λ 1g − λ2α) = min α (( λ1 mc − λ2)α) − v+− λ1 mDv 2− λ 1g. (16) Vi inför nu q = cλ1
m − λ2. Precis som i det regulariserade testproblemet inser man då
att α vill anta sitt maximala värde då q < 0, dvs α = 1. I det andra fallet har vi α = 0. Vi vill alltså ha en slät funktion som är nära q för små q och annars 0. Vi kan då välja q
2(1 − tanh q
δ). För v+ vill vi istället göra tvärtom, dvs v δ + = v2(1 + tanh v δ). Sammanfattningsvis har vi nu Hδ = q 2(1 − tanh q δ) − v δ +− D mλ1v − gλ1.
Vi ska nu försöka nna en bra startgissning. Första tanken man får är kanske: hur högt kommer vi på att köra med full gas så länge bränslet räcker? Detta illustreras i Figur 5.
Figur 5: Resultatet för fallet då styrsignalen hålls konstant vid 100%.
Det visar sig tyvärr att vårt program har väldigt svårt att lösa detta och denna start-gissning duger inte. Problemet verkar vara väldigt känsligt för vilken startstart-gissning man ger och vi har aldrig lyckats få konvergens. Däremot kan vi pröva att exempelvis istället för att köra med full gas gå ned till hälften dvs α(t) = 0.5 kg/s. Detta illustreras i Figur 6.
Figur 6: Resultatet för fallet då styrsignalen hålls konstant vid 50%.
Det är tydligt att denna styrsignal är bättre. Den fysikaliska anledningen till det-ta är att vi håller en lägre hastighet under en längre tid och eftersom luftmotståndet är proportionellt mot hastigheten i kvadrat lönar detta sig. Men detta är inte den bästa styr-signalen. Det lönar sig nämligen att nå en konstant hastighet tidigt och behålla den tills bränslet tar slut och därmed minimera tiden raketen behöver accelerera. Detta illustreras i Figur 7.
Figur 7: Optimerad styrsignal med ett steg.
Men nns det en bättre styrsignal? Det visar sig att denna styrsignal ej löser det Hamiltonska systemet och därmed bör det nnas en bättre styrsignal än denna. Även denna startgissning duger tyvärr inte.
Figurerna ovan är plottade med ett sidoprogram vi skrev i Matlab med namnet hcalc.m som följer med lösaren. Vill man se hur högt man kommer med en styrsig-nal α kan man testa det med kommandot hcalc(α, dt, plots). Där dt är tidssteget mellan elementen i α och plots är antingen true eller false beroende på om man vill få resultatet plottat eller inte. Motsvarande maxhöjd returneras alltid.
10 Vidareutveckling av lösaren
Vi avslutar denna rapport med att nämna några förslag på vidareutveckling av koden. Vi är nöjda med vår numerisk lösare för optimala styrproblem men vi har samtidigt några idéer på vad som kan förbättras. Nedan kommer en lista på dessa idéer men på grund av tidsbrist har vi inte kunnat implementera dem.
• Newton direkt - Istället för att använda sig av Matlabs inbyggda lösare för ekva-tionssystem fsolve skulle man kunna skriva en lösare själv som utgår från Newtons metod. Eftersom valet av de nition på y medförde en bandad Jacobian verkar det vid första anblick vara en god idé då Matlabs 'backslash' löser bandade matrisekva-tioner eektivt. Fördelen med fsolve är förstås att den är lätt att implementera och har en hel del inbyggd funktionalitet.
• Valfrihet kring inmatning av analytiska derivator - I nuläget som nämnts tidigare har vi enbart en parameter som styr om numerisk eller analytisk derivata ska användas. I vissa fall kanske man bara beräknat första derivatan av Hamiltonianen eller något annat urval och vill då kunna använda sig av detta. Ett sätt är kanske att införa en array av parametrar som säger vilka delar som ska beräknas numeriskt respektive analytiskt.
• Metod som stegvis ökar sluttiden T - Ofta får man lättare konvergens om slut-tiden T väljs kort. Därmed skulle man kunna ha en metod som börjar med att lösa problemet för en kort tid T sedan använda lösningen som startgissning för ett större T. Detta skulle möjligen kunna hjälpa till att lösa raketuppskjutningsproblemet.
