(!*********************************************************************
   Mosel NL examples
   =================
   file catenary.mos
   `````````````````
   Find the shape of a hanging chain by minimising its potential energy
   QCQP problem (linear objective, convex quadratic constraints)

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

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

model "catenary"
 uses "mmxnlp"

 parameters
  N = 100                       ! Number of chainlinks
  L = 1                         ! Difference in x-coordinates of endlinks
  H = 2*L/N                     ! Length of each link
  A = 0.5                       ! Height of start point, value in [-1,1]
  B = 0.1                       ! Height of end point, value in [-1,1]
 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 RN | j>0) (y(j-1)+y(j))/2

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

! 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)

 setparam("XNLP_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