(!********************************************************************* 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, rev. Jun. 2023 *********************************************************************!) 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); ! 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) ! 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