(!*********************************************************************
   Mosel NL examples
   =================
   file steiner_graph.mos
   ``````````````````````
   Steiner tree problem.
   SOCP formulation. 

   Based on AMPL model steiner_socp.mod
   Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/steiner/ 

   - Graphical representation of results -   

   (c) 2013 Fair Issac Corporation
       author: S. Heipcke, Sep. 2013, rev. Jun. 2023
*********************************************************************!)
model "steiner"
 uses "mmxnlp", "mmsvg"

 parameters
  DATAFILE = "steiner.dat"
 end-parameters

 declarations
  NNODES: integer           ! Number of nodes (original plus steiner)
  NSTEINER: integer         ! Number of steiner nodes
 end-declarations

 initialisations from DATAFILE
  NNODES
  NSTEINER
 end-initialisations

 declarations
  SNODES = 1..NSTEINER                    ! Steiner nodes
  ONODES = NSTEINER+1..NNODES             ! Original nodes
  NODES = 1..NNODES                       ! All nodes
  DIM: range                              ! Dimensions of coordinates
  POS: array(ONODES,DIM) of real          ! Coordinates of original nodes
  ARCS: dynamic array(NODES,NODES) of boolean    ! Arcs between nodes
  x: array(NODES,DIM) of mpvar            ! Node positions
  t: dynamic array(NODES,NODES) of mpvar  ! Distance between nodes
 end-declarations

 initialisations from DATAFILE
  ARCS
  POS
 end-initialisations

! Decision variables
 forall(i in NODES,k in DIM) x(i,k) is_free

 forall(i,j in NODES | exists(ARCS(i,j))) do
  create(t(i,j))
  t(i,j)>=0
 end-do

! Objective: total distance
 Dist:= sum(i,j in NODES | exists(ARCS(i,j))) t(i,j)

! Constraints
 forall(i,j in NODES | exists(ARCS(i,j))) 
   Bnd_t(i,j):= sqrt(sum(k in DIM) (x(i,k)-x(j,k))^2) <= t(i,j)

! Fix position for original nodes
 forall(j in ONODES,k in DIM) x(j,k)= POS(j,k)

! Start values
 setrandseed(3)
 forall(j in SNODES,k in DIM) setinitval(x(j,k),random)
 
! Since this is a convex problem, it is sufficient to call a local solver
 setparam("xprs_nlpsolver", 1)

! Solving
 setparam("xnlp_verbose", true)
 minimize(Dist)

! Solution reporting
 writeln("Solution: ", Dist.sol)
 forall(i in NODES) do
   write(i, ": ")
   forall(k in DIM) write(x(i,k).sol, " ")
   writeln
 end-do  

! **** Display the solution as user graph ****
 
! Draw the nodes and arcs (2-dim)
 if getsize(DIM)<2 then exit(0); end-if
 svgaddgroup("SGr", "Steiner nodes", SVG_BLUE)
 svgaddgroup("OGr", "Original nodes", SVG_SILVER)
 forall(i in ONODES) do
   svgaddpoint("OGr", POS(i,1), POS(i,2))
   svgaddtext("OGr", POS(i,1)+0.2, POS(i,2)+0.2, text(i))
 end-do
 forall(i in SNODES) do
   svgaddpoint("SGr", x(i,1).sol, x(i,2).sol)
   svgaddtext("SGr",  x(i,1).sol-0.4, x(i,2).sol+0.2, text(i))
 end-do
 forall(i,j in NODES | exists(ARCS(i,j)))
   svgaddline("SGr", x(i,1).sol, x(i,2).sol, x(j,1).sol, x(j,2).sol)

! Scale the size of the displayed graph
 svgsetgraphscale(30)
 svgsetgraphpointsize(3)
 
 svgsave("sphere.svg")
 svgrefresh
 svgwaitclose("Close browser window to terminate model execution.", 1)
end-model