(!*********************************************************************
Mosel NL examples
=================
file springs.mos
````````````````
Find the shape of a hanging chain by minimising its potential energy
where each chain link is a spring that buckles under compression and
each node (connection between links) has a weight hanging from it.
The springs are assumed weightless.
SOCP problem (convex quadratic objective, convex second-order cone
constraints)
Based on AMPL model springs.mod
Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/noncute/springs/
Reference: "Applications of Second-Order Cone Programming",
M.S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, 1998
(c) 2013 Fair Issac Corporation
author: S. Heipcke, Sep. 2013
*********************************************************************!)
model "springs"
uses "mmxnlp"
parameters
N = 15 ! Number of chainlinks
! A = (AX, AY) and B = (BX, BY) are the positions of end nodes
AX = 0
AY = 0
BX = 2
BY = -1
G = 9.8 ! Acceleration due to gravity
K = 100 ! Stiffness of springs
end-parameters
declarations
RN = 0..N
RN1 = 1..N
x: array(RN) of mpvar ! x-coordinates of endpoints of chainlinks
y: array(RN) of mpvar ! y-coordinates of endpoints of chainlinks
MASS: array(RN) of real ! Mass of each hanging node
t: array(RN1) of mpvar ! Extension of each spring
LENGTH0: real ! Rest length of each spring
end-declarations
! Initializing data
LENGTH0 := 2*sqrt((AX-BX)^2+(AY-BY)^2)/N
! Try different settings for node weights:
forall(j in RN | j>0 and j0 and j0 and j0) t(j)>=0
! Objective: minimise the potential energy
potential_energy:= sum(j in RN) MASS(j)*G*y(j) + (K/2)*sum(j in RN1) t(j)^2
! Bounds: positions of endpoints
! Left anchor
x(0) = AX; y(0) = AY
! Right anchor
x(N) = BX; y(N) =BY
! Constraints: positions of chainlinks
forall(j in RN| j>0)
Link(j):= sqrt((x(j)-x(j-1))^2+(y(j)-y(j-1))^2) <= LENGTH0+t(j);
! Setting start values
forall(j in RN) setinitval(x(j), (j/N)*BX+(1-j/N)*AX);
forall(j in RN) setinitval(y(j), (j/N)*BY+(1-j/N)*AY);
! Solving
setparam("XPRS_verbose", true)
minimise(potential_energy)
! Solution reporting
writeln("Solution: ", getobjval)
forall(j in RN)
writeln(j, ": (", strfmt(getsol(x(j)),10,5), ", ",
strfmt(getsol(y(j)),10,5), ")")
end-model