Catenary: Determine chain shape

(!*********************************************************************
   Mosel NL examples
   =================
   file catenary.mos
   `````````````````
   QCQP problem (linear objective, convex quadratic constraints)
   Based on AMPL model catenary.mod
   (Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/ )

   This model finds the shape of a hanging chain by
   minimizing its potential energy.

   (c) 2008 Fair Issac Corporation
       author: S. Heipcke, May 2008
*********************************************************************!)

model "catenary"
 uses "mmxprs", "mmnl"

 parameters
  N = 100                       ! Number of chainlinks
  L = 1                         ! Difference in x-coordinates of endlinks
  H = 2*L/N                     ! Length of each link
 end-parameters

 declarations
  RN = 0..N
  x: array(RN) of mpvar         ! x-coordinates of endpoints of chainlinks
  y: array(RN) of mpvar         ! y-coordinates of endpoints of chainlinks
 end-declarations

 forall(i in RN) x(i) is_free
 forall(i in RN) y(i) is_free

! Objective: minimise the potential energy
 potential_energy:= sum(j in 1..N) (y(j-1)+y(j))/2

! Bounds: positions of endpoints
! Left anchor
  x(0) = 0; y(0) = 0
! Right anchor
  x(N) = L; y(N) = 0

! Constraints: positions of chainlinks
 forall(j in 1..N) 
  Link(j):= (x(j)-x(j-1))^2+(y(j)-y(j-1))^2 <= H^2

! Setting start values
 forall(j in RN) setinitval(x(j), j*L/N)
 forall(j in RN) setinitval(y(j), 0)

(! For exporting the matrix file
 setparam("XPRS_loadnames", true)
 loadprob(potential_energy)
 writeprob("catenary.mat","l")
!)

 setparam("XPRS_verbose", true)
 minimise(potential_energy)

 writeln("Solution: ", getobjval)
 forall(j in RN) 
  writeln(strfmt(getsol(x(j)),10,5), " ", strfmt(getsol(y(j)),10,5))

end-model