(!*********************************************************************
   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 -   

   *** This model cannot be run with a Community Licence 
       for the provided data instance ***

   (c) 2013 Fair Issac Corporation
       author: S. Heipcke, Mar. 2013, rev. Jun. 2023
*********************************************************************!)
    
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)  
        
! In this example we will use a local solver, since it can be time consuming to solve it to global optimality
  setparam("xprs_nlpsolver", 1)
   
! 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