$TITLE Equilibrium of System with Piecewise Linear Springs - Conic Formulation (SPRINGS)

$ONTEXT

 This model finds the shape of a hanging chain consisitng of
 N springs and N-1 nodes. Each spring buckles under compression and each
 node has a weight hanging from it.  The springs are assumed 
 weightless. The goal is minimize the potential energy of the
 system.

 We use rotated quadratic cone constraints to model the extension 
 of each spring.

 M. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret ,  
 "Applications of second-order cone programming", Linear Algebra and its 
 Applications, 284:193-228, November 1998, Special Issue on Linear Algebra 
 in Control, Signals and Image Processing. 

 Vanderbei, R, online at 
 http://www.princeton.edu/~rvdb/ampl/nlmodels/facloc/springs.mod 

$OFFTEXT

*=== Number of chainlinks
$if not set N $set N 10

Set n "spring index"  /n0*n%N%/;


*=== (a_x, a_y) and (b_x, b_y) are coordinates of end nodes.
Scalar a_x "x coordinate of beginning node" /0/
       a_y "y coordinate of beginning node" /0/      
       b_x "x coordinate of end node"       /2/ 
       b_y "y coordinate of end node"       /-1/
;       

Scalar L0 "rest length of each spring";
L0 = 2*sqrt[sqr(a_x-b_x) + sqr(a_y-b_y)]/%N%;  

Scalar g    "acceleration due to gravity" /9.8/
       k    "stiffness of springs"        /100/; 
Parameter 
       m{n} "mass of each hanging node";
m(n)$[ord(n)>1 and ord(n)<card(n)] = 1;


Variable 
       obj
       x(n)       "x-coordinates of nodes"
       y(n)       "y-coordinates of nodes"
       t_L0(n)
       unit
       delta_x(n)
       delta_y(n) 
;        
Positive variable 
       t(n)  "extension of each spring"
       v    
;


Equation
       pot_energy
       delta_x_eq(n)
       delta_y_eq(n)
       unit_eq
       link_L0(n)
       link_up(n)
       cone_eq
;

pot_energy..       obj =E=    sum(n$[ord(n)>1 and ord(n)<card(n)], m[n]*g*y[n]) + k*v;

delta_x_eq(n)..    delta_x(n)     =E= x[n] - x[n-1];
delta_y_eq(n)..    delta_y(n)     =E= y[n] - y[n-1];

unit_eq..          unit =E= 1;

link_L0(n)..       t_L0[n]        =E= L0 + t[n];
link_up(n)$[ord(n)>1]..       
                   t_L0[n]        =C= delta_x[n] + delta_y[n];

cone_eq..          v + unit       =C= sum(n$[ord(n)>1], t[n]);

Model spring /all/;

x.L(n) = ( (ord(n)-1)/%N% )*b_x + (ord(n)/%N%)*a_x;
y.L(n) = ( (ord(n)-1)/%N% )*b_y + (ord(n)/%N%)*a_y;

x.FX['n0']   = a_x;
y.FX['n0']   = a_y;
x.FX['n%N%'] = b_x;
y.FX['n%N%'] = b_y;

Solve spring using lp minimizing obj;

display x.L,y.L;