Minimize Energy of Piecewise Linear Mass-Spring System Using Cone Programming, Solver-Based

Mathematical Formulation

The location of mass $$i$$ is $$x(i)$$, with horizontal coordinate $$x_1(i)$$ and vertical coordinate $$x_2(i)$$. Mass $$i$$ has potential energy due to gravity of $$gm(i)x_2(i)$$. The potential energy in spring $$i$$ is $$k(i)(d(i)-l(i){)}^2/2$$, where $$d(i)$$ is the length of the spring between mass $$i$$ and mass $$i-1$$. Take anchor point 1 as the position of mass 0, and anchor point 2 as the position of mass $$n+1$$. The preceding energy calculation shows that the potential energy of spring $$i$$ is

$$Energy(i)=\frac{k(i)(\Vert x(i)-x(i-1)\Vert -l(i){)}^2}{2}$$.

Reformulating this potential energy problem as a second-order cone problem requires the introduction of some new variables, as described in Lobo [1]. Create variables $$t(i)$$ equal to the square root of the term $$Energy(i)$$.

Let $$e$$ be the unit column vector $$\left[\genfrac{}{}{0ex}{}{0}{1}\right]$$. Then $$x_2(i)=e^{T}x(i)$$. The problem becomes

$$\underset{x,t}{\min }(\sum _{i}gm(i)e^{T}x(i)+\Vert t{\Vert }^2).$$               (1)

Now consider $$t$$ as a free vector variable, not given by the previous equation for $$t(i)$$. Incorporate the relationship between $$x(i)$$ and $$t(i)$$ in the new set of cone constraints

$$\Vert x(i)-x(i-1)\Vert -l(i)\le \sqrt{\frac{2}{k(i)}}t(i).$$   (2)

The objective function is not yet linear in its variables, as required for coneprog. Introduce a new scalar variable $$y$$. Notice that the inequality $$\Vert t{\Vert }^2\le y$$ is equivalent to the inequality

$$\Vert \left[\genfrac{}{}{0ex}{}{2t}{1-y}\right]\Vert \le 1+y$$.                               (3)

Now the problem is to minimize

$$\underset{x,t,y}{\min }(\sum _{i}gm(i)e^{T}x(i)+y)$$                     (4)

subject to the cone constraints on $$x(i)$$ and $$t(i)$$ listed in (2) and the additional cone constraint (3). Cone constraint (3) ensures that $$\Vert t{\Vert }^2\le y$$. Therefore, problem (4) is equivalent to problem (1).

The objective function and cone constraints in problem (4) are suitable for solution with coneprog.

MATLAB® Formulation

Define six spring constants $$k$$, six length constants $$l$$, and five masses $$m$$.

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$$.

g = 9.807;

The variables for optimization are the ten components of the $$x$$ vectors, the six components of the $$t$$ vector, and the $$y$$ variable. Let v be the vector containing all these variables.

Using these variables, create the corresponding objective function vector f.

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)

$$\Vert x(i)-x(i-1)\Vert -l(i)\le \sqrt{\frac{2}{k(i)}}t(i)$$.

The coneprog solver uses cone constraints for a variable vector $$v$$ in the form

$$\Vert A_{sc}\cdot v-b_{sc}\Vert \le d_{sc}^{T}v-\gamma$$.

In the following code, the Asc matrix represents the term $$\Vert x(i)-x(i-1)\Vert$$, with bsc = [0;0]. The cone variable dsc = $$\sqrt{2/k(i)}$$ and the corresponding gamma = $$-l(i).$$

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 $$y$$ variable

by creating the matrix Asc which, when multiplied by the v vector, gives the vector $$\left[\genfrac{}{}{0ex}{}{2t}{-y}\right]$$. The bsc vector corresponds to the constant 1 in the term $$1-y$$. The dsc vector, when multiplied by v, returns $$y$$. And gamma = $$-1$$.

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 $$t$$ and $$y$$ 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 coneprog.

[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 m, l, and 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.

References

[1] 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. https://doi.org/10.1016/S0024-3795(98)10032-0.