(!*********************************************************************
Mosel NL examples
=================
file moonshot.mos
`````````````````
Minimise fuel consumption of a rocket from earth to moon.
Nonlinear objective and constraints
Based on AMPL model moonshot.mod
Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/rocket/
(c) 2013 Fair Issac Corporation
author: S. Heipcke, Mar. 2013
*********************************************************************!)
model "moonshot"
uses "mmxnlp"
parameters
N = 100 ! Number of discrete times
end-parameters
declarations
D = 1..2 ! Dimensions
Nx: set of real ! Discrete times for positions
Nv: set of real ! Discrete times for velocity
Na: set of real ! Discrete times for acceleration
T: real ! Total time
G: real ! Constant factor
MASS_Earth: real ! Mass of Earth
RADIUS_Earth: real ! Radius of Earth
POS_Earth: array(D) of real ! Position(x,y) of Earth
MASS_Moon : real ! Mass of Moon
RADIUS_Moon : real ! Radius of Moon
POS_Moon : array(D) of real ! Position(x,y) of Moon
X0: array(D) of real ! Initial position
Xn: array(D) of real ! Final position
V0: array(D) of real ! Initial velocity
Vn: array(D) of real ! Final velocity
ALPHA0: real ! Initial degree
ALPHAn: real ! Final degree
end-declarations
Nx:= union(i in 0..N) {real(i)}
Na:= union(i in 1..N-1) {real(i)}
Nv:= union(i in 1..N) {i-0.5}
initializations from "moonshot.dat"
T G MASS_Earth MASS_Moon RADIUS_Earth RADIUS_Moon POS_Earth POS_Moon
end-initializations
! Initial and final velocity
ALPHA0 := 1.5*M_PI
ALPHAn := M_PI/2.0
X0(1) := POS_Earth(1) + RADIUS_Earth*cos(ALPHA0)
X0(2) := POS_Earth(2) + RADIUS_Earth*sin(ALPHA0)
Xn(1) := POS_Moon(1) + RADIUS_Moon*cos(ALPHAn)
Xn(2) := POS_Moon(2) + RADIUS_Moon*sin(ALPHAn)
V0(1) := 5*cos(ALPHA0)
V0(2) := 5*sin(ALPHA0)
Vn(1) := -5*cos(ALPHAn)
Vn(2) := -5*sin(ALPHAn)
declarations
x: array(D, Nx) of mpvar ! Position
v: array(D, Nv) of mpvar ! Velocity
a: array(D, Na) of mpvar ! Acceleration
theta: array(D, Na) of mpvar ! Thrust
fuelcost: nlctr
end-declarations
! Objective function: total fuel consumption
Fuelcost:= sum(d in D, i in Na) theta(d,i)^2
! Setting start values for positions
forall(d in D, j in Nx) do
setinitval(x(d,j), (1-j/N)*X0(d)+(j/N)*Xn(d) )
x(d,j) is_free
end-do
! Fixing start and final position and velocity
forall(d in D) do
x(d, 0) = X0(d) ! Initial position
x(d, N) = Xn(d) ! Final position
v(d, 0.5) = V0(d) ! Initial velocity
v(d, N-0.5)^1 = Vn(d) ! Final velocity
end-do
! Velocity constraint
forall(d in D, j in Nv) do
VelDef(d,j) := N*( x(d, j+0.5) - x(d, j-0.5) ) = T*v(d,j)
v(d,j) is_free
end-do
! Acceleration constraint
forall(d in D, j in Na) do
AccDef(d,j) := N*( v(d, j+0.5) - v(d, j-0.5) ) = T*a(d,j)
a(d,j) is_free
theta(d,j) is_free
end-do
! Force balance constraint
forall(d in D, j in Na)
Force(d,j) := a(d,j) =
-G*MASS_Earth*(x(d,j)-POS_Earth(d))/
(sum(dd in D) (x(dd,j)-POS_Earth(dd))^2)^(3/2) -
G*MASS_Moon*(x(d,j)-POS_Moon(d))/
(sum(dd in D) (x(dd,j)-POS_Moon(dd))^2)^(3/2) +
theta(d,j)
! Setting up and solving the problem
setparam("xnlp_verbose",true)
minimize(Fuelcost)
! Solution printing
forall(i in Na)
writeln("i=", i, ": position=(", x(1,i).sol, ",", x(2,i).sol, ")",
" thrust=(",theta(1,i).sol, ",", theta(2,i).sol,"), ",
"acceleration=(",getsol(a(1,i)), ",", getsol(a(2,i)),")" )
forall(i in Nv) writeln("velocity=(",getsol(v(1,i)), ",",getsol(v(2,i)),")" )
writeln("Fuelcost = ", getobjval, "\n")
end-model
|
(!*********************************************************************
Mosel NL examples
=================
file moonshot_graph.mos
```````````````````````
Minimise fuel consumption of a rocket from earth to moon.
Nonlinear objective and constraints
Based on AMPL model moonshot.mod
Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/rocket/
- Graphical representation of results -
(c) 2013 Fair Issac Corporation
author: S. Heipcke, Mar. 2013, rev. Sep. 2017
*********************************************************************!)
model "moonshot"
uses "mmxnlp", "mmsvg"
parameters
N = 100 ! Number of discrete times
end-parameters
declarations
D = 1..2 ! Dimensions
Nx: set of real ! Discrete times for positions
Nv: set of real ! Discrete times for velocity
Na: set of real ! Discrete times for acceleration
T: real ! Total time
G: real ! Constant factor
MASS_Earth: real ! Mass of Earth
RADIUS_Earth: real ! Radius of Earth
POS_Earth: array(D) of real ! Position(x,y) of Earth
MASS_Moon : real ! Mass of Moon
RADIUS_Moon : real ! Radius of Moon
POS_Moon : array(D) of real ! Position(x,y) of Moon
X0: array(D) of real ! Initial position
Xn: array(D) of real ! Final position
V0: array(D) of real ! Initial velocity
Vn: array(D) of real ! Final velocity
ALPHA0: real ! Initial degree
ALPHAn: real ! Final degree
end-declarations
Nx:= union(i in 0..N) {real(i)}
Na:= union(i in 1..N-1) {real(i)}
Nv:= union(i in 1..N) {i-0.5}
initializations from "moonshot.dat"
T G MASS_Earth MASS_Moon RADIUS_Earth RADIUS_Moon POS_Earth POS_Moon
end-initializations
! Initial and final velocity
ALPHA0 := 1.5*M_PI
ALPHAn := M_PI/2.0
X0(1) := POS_Earth(1) + RADIUS_Earth*cos(ALPHA0)
X0(2) := POS_Earth(2) + RADIUS_Earth*sin(ALPHA0)
Xn(1) := POS_Moon(1) + RADIUS_Moon*cos(ALPHAn)
Xn(2) := POS_Moon(2) + RADIUS_Moon*sin(ALPHAn)
V0(1) := 5*cos(ALPHA0)
V0(2) := 5*sin(ALPHA0)
Vn(1) := -5*cos(ALPHAn)
Vn(2) := -5*sin(ALPHAn)
declarations
x: array(D, Nx) of mpvar ! Position
v: array(D, Nv) of mpvar ! Velocity
a: array(D, Na) of mpvar ! Acceleration
theta: array(D, Na) of mpvar ! Thrust
fuelcost: nlctr
end-declarations
! Objective function: total fuel consumption
Fuelcost:= sum(d in D, i in Na) theta(d,i)^2
! Setting start values for positions
forall(d in D, j in Nx) do
setinitval(x(d,j), (1-j/N)*X0(d)+(j/N)*Xn(d) )
x(d,j) is_free
end-do
! Fixing start and final position and velocity
forall(d in D) do
x(d, 0) = X0(d) ! Initial position
x(d, N) = Xn(d) ! Final position
v(d, 0.5) = V0(d) ! Initial velocity
v(d, N-0.5)^1 = Vn(d) ! Final velocity
end-do
! Velocity constraint
forall(d in D, j in Nv) do
VelDef(d,j) := N*( x(d, j+0.5) - x(d, j-0.5) ) = T*v(d,j)
v(d,j) is_free
end-do
! Acceleration constraint
forall(d in D, j in Na) do
AccDef(d,j) := N*( v(d, j+0.5) - v(d, j-0.5) ) = T*a(d,j)
a(d,j) is_free
theta(d,j) is_free
end-do
! Force balance constraint
forall(d in D, j in Na)
Force(d,j) := a(d,j) =
-G*MASS_Earth*(x(d,j)-POS_Earth(d))/
(sum(dd in D) (x(dd,j)-POS_Earth(dd))^2)^(3/2) -
G*MASS_Moon*(x(d,j)-POS_Moon(d))/
(sum(dd in D) (x(dd,j)-POS_Moon(dd))^2)^(3/2) +
theta(d,j)
! Setting up and solving the problem
setparam("xnlp_verbose",true)
minimize(Fuelcost)
! Solution printing
forall(i in Na)
writeln("i=", i, ": position=(", x(1,i).sol, ",", x(2,i).sol, ")",
" thrust=(",theta(1,i).sol, ",", theta(2,i).sol,"), ",
"acceleration=(",getsol(a(1,i)), ",", getsol(a(2,i)),")" )
forall(i in Nv) writeln("velocity=(",getsol(v(1,i)), ",",getsol(v(2,i)),")" )
writeln("Fuelcost = ", getobjval, "\n")
! **** Solution display as IVE user graph ****
declarations
plot1,plot2,earth,moon: integer
end-declarations
SCALE:=2 ! Scaling factor between x and y axes
! Draw Earth and Moon
! Overlay circles of different size+colour to obtain 3D effect
svgaddgroup("E", "Earth", SVG_BLUE)
svgaddgroup("M", "Moon", SVG_GREY)
forall(sz in union(i in 0..50) {i*1/50}) do
offset:=0.2*sz
svgaddcircle("E", POS_Earth(1)+sz+offset, POS_Earth(2)+sz+offset/2, RADIUS_Earth*SCALE*(1.02-sz) )
col:=svgcolor(75+round(50*sz),100+round(100*sz),225+round(25*sz))
svgsetstyle(svggetlastobj, SVG_FILL, col)
svgsetstyle(svggetlastobj, SVG_STROKE, col)
svgaddcircle("M", POS_Moon(1)+sz-offset, POS_Moon(2)+sz-offset/2, RADIUS_Moon*SCALE*(1.02-sz))
col:=svgcolor(100+round(125*sz),100+round(125*sz),100+round(125*sz))
svgsetstyle(svggetlastobj, SVG_STROKE,col)
svgsetstyle(svggetlastobj, SVG_FILL,col)
end-do
! Draw acceleration at selected points
svgaddgroup("A", "Acceleration", SVG_RED)
svgsetstyle(SVG_STROKEWIDTH, 2)
forall (i in Na | i mod 4 = 0 and
abs(getsol(a(1,i)))+abs(getsol(a(2,i)))>1 ) do
col:=svgcolor(150+round(0.75*getsol(a(1,i))), 100-minlist(round(0.75*getsol(a(1,i))),100),0)
svgaddarrow(x(1,i).sol, x(2,i).sol*SCALE, x(1,i).sol+0.02*getsol(a(1,i)), x(2,i).sol+0.02*getsol(a(2,i))*SCALE)
svgsetstyle(svggetlastobj, SVG_STROKE,col)
svgsetstyle(svggetlastobj, SVG_FILL,col)
end-do
! Draw trajectory from Earth to Moon
svgaddgroup("T", "MoonShot Orbit", SVG_BLACK)
svgaddline(sum(i in Nx) [x(1,i).sol, x(2,i).sol*SCALE])
! Scale the size of the displayed graph
svgsetgraphscale(20)
svgsetgraphviewbox(-5,-10,30,30)
svgsave("moonshot.svg")
svgrefresh
svgwaitclose("Close browser window to terminate model execution.", 1)
end-model
|
