Minimize Energy of Piecewise Linear Mass-Spring System Using Cone Programming, Solver-Based
This example shows how to find the equilibrium position of a mass-spring system hanging from two anchor points. The springs have piecewise linear tensile forces. The system consists of masses in two dimensions. Mass is connected to springs and . Springs and are also connected to separate anchor points. In this case, the zero-force length of spring is a positive length , and the spring generates force when stretched to length . The problem is to find the minimum potential energy configuration of the masses, where potential energy comes from the force of gravity and from the stretching of the nonlinear springs. The equilibrium occurs at the minimum energy configuration.
This illustration shows five springs and four masses suspended from two anchor points.
The potential energy of a mass at height is , where is the gravitational constant on Earth. Also, the potential energy of an ideal linear spring with spring constant stretched to length is . The current model is that the spring is not ideal, but has a nonzero resting length .
The mathematical basis of this example comes from Lobo, Vandenberg, Boyd, and Lebret . For a problem-based version of this example, see Minimize Energy of Piecewise Linear Mass-Spring System Using Cone Programming, Problem-Based.
The location of mass is , with horizontal coordinate and vertical coordinate . Mass has potential energy due to gravity of . The potential energy in spring is , where is the length of the spring between mass and mass . Take anchor point 1 as the position of mass 0, and anchor point 2 as the position of mass . The preceding energy calculation shows that the potential energy of spring is
Reformulating this potential energy problem as a second-order cone problem requires the introduction of some new variables, as described in Lobo . Create variables equal to the square root of the term .
Let be the unit column vector . Then . The problem becomes
Now consider as a free vector variable, not given by the previous equation for . Incorporate the relationship between and in the new set of cone constraints
The objective function is not yet linear in its variables, as required for
coneprog. Introduce a new scalar variable . Notice that the inequality is equivalent to the inequality
Now the problem is to minimize
subject to the cone constraints on and listed in (2) and the additional cone constraint (3). Cone constraint (3) ensures that . Therefore, problem (4) is equivalent to problem (1).
The objective function and cone constraints in problem (4) are suitable for solution with
Define six spring constants , six length constants , and five masses .
k = 40*(1:6); l = [1 1/2 1 2 1 1/2]; m = [2 1 3 2 1];
Define the approximate gravitational constant on Earth .
g = 9.807;
The variables for optimization are the ten components of the vectors, the six components of the vector, and the variable. Let
v be the vector containing all these variables.
[v(1),v(2)]corresponds to the 2-D variable .
[v(3),v(4)]corresponds to the 2-D variable .
[v(5),v(6)]corresponds to the 2-D variable .
[v(7),v(8)]corresponds to the 2-D variable .
[v(9),v(10)]corresponds to the 2-D variable .
[v(11):v(16)]corresponds to the 6-D vector .
v(17)corresponds to the scalar variable .
Using these variables, create the corresponding objective function vector
f = zeros(size(m)); f = [f;g*m]; f = f(:); f = [f;zeros(length(k)+1,1)]; f(end) = 1;
Create the cone constraints corresponding to the springs between the masses (2)
coneprog solver uses cone constraints for a variable vector in the form
In the following code, the
Asc matrix represents the term , with
[0;0]. The cone variable
dsc = and the corresponding
d = zeros(1,length(f)); Asc = d; Asc([1 3]) = [1 -1]; A2 = circshift(Asc,1); Asc = [Asc;A2]; ml = length(m); dbase = 2*ml; bsc = [0;0]; for i = 2:ml gamma = -l(i); dsc = d; dsc(dbase + i) = sqrt(2/k(i)); conecons(i) = secondordercone(Asc,bsc,dsc,gamma); Asc = circshift(Asc,2,2); end
Create the cone constraints corresponding to the springs between the end masses and the anchor points by using the anchor points for the locations of the end masses, as in the preceding code.
x0 = [0;5]; xn = [5;4]; Asc = zeros(size(Asc)); Asc(1,(dbase-1)) = 1; Asc(2,dbase) = 1; bsc = xn; gamma = -l(ml); dsc = d; dsc(dbase + ml) = sqrt(2/k(ml)); conecons(ml + 1) = secondordercone(Asc,bsc,dsc,gamma); Asc = zeros(size(Asc)); Asc(1,1) = 1; Asc(2,2) = 1; bsc = x0; gamma = -l(1); dsc = d; dsc(dbase + 1) = sqrt(2/k(1)); conecons(1) = secondordercone(Asc,bsc,dsc,gamma);
Create the cone constraint (3) corresponding to the variable
by creating the matrix
Asc which, when multiplied by the
v vector, gives the vector . The
bsc vector corresponds to the constant 1 in the term . The
dsc vector, when multiplied by
v, returns . And
gamma = .
Asc = 2*eye(length(f)); Asc(1:dbase,:) = ; Asc(end,end) = -1; bsc = zeros(size(Asc,1),1); bsc(end) = -1; dsc = d; dsc(end) = 1; gamma = -1; conecons(ml+2) = secondordercone(Asc,bsc,dsc,gamma);
Finally, create lower bounds corresponding to the and variables.
lb = -inf(size(f)); lb(dbase+1:end) = 0;
Solve Problem and Plot Solution
The problem formulation is complete. Solve the problem by calling
[v,fval,exitflag,output] = coneprog(f,conecons,,,,,lb);
Optimal solution found.
Plot the solution points and the anchors.
pp = v(1:2*ml); pp = reshape(pp,2,); pp = pp'; plot(pp(:,1),pp(:,2),'ro') hold on xx = [x0,xn]'; plot(xx(:,1),xx(:,2),'ks') xlim([x0(1)-0.5,xn(1)+0.5]) ylim([min(pp(:,2))-0.5,max(x0(2),xn(2))+0.5]) xxx = [x0';pp;xn']; plot(xxx(:,1),xxx(:,2),'b--') legend('Calculated points','Anchor points','Springs','Location',"best") hold off
You can change the values of the parameters
k to see how they affect the solution. You can also change the number of masses; the code takes the number of masses from the data you supply.
 Lobo, Miguel Sousa, Lieven Vandenberghe, Stephen Boyd, and Hervé Lebret. “Applications of Second-Order Cone Programming.” Linear Algebra and Its Applications 284, no. 1–3 (November 1998): 193–228.