In this notebook we show how to solve the following problem: Find the equlibrium of a system of masses connected by a system of strings, with some masses being assigned fixed coordinates (attached to the wall, say). See the next picture.

Suppose we have $n$ masses with weights $w_1,\ldots,w_n$, and the length of the string between $i$ and $j$ is $\ell_{ij}$ for some set $L$ of pairs of indices $(i,j)$ (we assume $\ell_{ij}$ is not defined if there is no connection). The strings themselves have no mass. We also have a set $F$ of indices such that the $i$-th point is fixed to have coordinates $f_i$ if $i\in F$. The equilibrium of the system is a configuration which minimizes potential energy. With this setup we can write our problem as:

\begin{equation} \begin{array}{ll} minimize & g\cdot \sum_i w_ix_i^{(2)} \\ s.t. & \|x_i-x_j\|\leq \ell_{ij},\ ij\in L \\ & x_i = f_i,\ i\in F \end{array} \end{equation}

where $x\in (\mathbf{R}^n)^2$, $x_i^{(2)}$ denotes the second (vertical) coordinate of $x_i$ and $g$ is the gravitational constant.

Here is a sample problem description.

In [1]:

```
w = [0.0, 1.1, 2.2, 0.0, 2.1, 2.2, 0.2]
l = {(0,1): 1.0, (1,2): 1.0, (2,3): 1.0, (1,4): 1.0, (4,5): 0.3, (5,2): 1.0, (5,6): 0.5, (1,3): 8.0}
f = {0: (0.0,1.0), 3: (2.0,1.0)}
g = 9.81
```

Now we can formulate the problem using Mosek Fusion:

In [2]:

```
from mosek.fusion import *
# w - masses of points
# l - lengths of strings
# f - coordinates of fixed points
# g - gravitational constant
def stringModel(w, l, f, g):
n, m = len(w), len(l)
starts = [ lKey[0] for lKey in l.keys() ]
ends = [ lKey[1] for lKey in l.keys() ]
M = Model("strings")
# Coordinates of points
x = M.variable("x", [n, 2])
# A is the signed incidence matrix of points and strings
A = Matrix.sparse(m, n, range(m)+range(m), starts+ends, [1.0]*m+[-1.0]*m)
# ||x_i-x_j|| <= l_{i,j}
c = M.constraint("c", Expr.hstack(Expr.constTerm(l.values()), Expr.mul(A, x)),
Domain.inQCone() )
# x_i = f_i for fixed points
for i in f:
M.constraint(x.slice([i,0], [i+1,2]), Domain.equalsTo(list(f[i])))
# sum (g w_i x_i_2)
M.objective(ObjectiveSense.Minimize,
Expr.mul(g, Expr.dot(w, x.slice([0,1], [n,2]))))
# Solve
M.solve()
if M.getProblemStatus(SolutionType.Interior) == ProblemStatus.PrimalAndDualFeasible:
return x.level().reshape([n,2]), c.dual().reshape([m,3])
else:
return None, None
```

Here is a quick description of how we use vectorization to deal with all the conic constraints in one go. The matrix $A$ is the incidence matrix between the masses and the strings, with coefficients $+1, -1$ for the two endpoints of each string. It is chosen so that the product $Ax$ has rows of the form

$$ (x_i^{(1)} - x_j^{(1)}, x_i^{(2)} - x_j^{(2)}) $$

for all pairs $i,j$ for which $\ell_{ij}$ is bounded. Stacking the values of $\ell$ in the left column produces a matrix with each row of the form

$$ (\ell_{ij}, x_i^{(1)} - x_j^{(1)}, x_i^{(2)} - x_j^{(2)}) $$

and a conic constraint is imposed on all the rows, as required.

The objective and linear constraints show examples of slicing the variable $x$.

The function returns the coordinates of the masses and the values of the dual conic variables. A zero dual value indicates that a particular string is hanging loose, and a nonzero value means it is fully stretched.

All we need now is to define a display function and we can look at some plots.

In [3]:

```
%matplotlib inline
# x - coordinates of the points
# c - dual values of string length constraints
# d - pairs of points to connect
def display(x, c, d):
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
# Plot points
ax.scatter(x[:,0], x[:,1], color="r")
# Plot fully stretched strings (nonzero dual value) as solid lines, else dotted lines
for i in range(len(c)):
col = "b" if c[i][0] > 1e-4 else "b--"
ax.plot([x[d[i][0]][0], x[d[i][1]][0]], [x[d[i][0]][1], x[d[i][1]][1]], col)
ax.axis("equal")
plt.show()
```

In [4]:

```
x,c = stringModel(w, l, f, g)
if x is not None:
display(x, c, l.keys())
```

How about we find a discrete approximation to the catenary:

In [5]:

```
n = 1000
w = [1.0]*n
l = {(i,i+1): 1.0/n for i in range(n-1)}
f = {0: (0.0,1.0), n-1: (0.7,1.0)}
g = 9.81
x,c = stringModel(w, l, f, g)
if x is not None:
display(x, c, l.keys())
```

We can also have more suspension points and more complicated shapes:

In [6]:

```
n = 20
w = [1.0]*n
l = {(i,i+1): 0.09 for i in range(n-1)}
l.update({(5,14): 0.3})
f = {0: (0.0,1.0), 13: (0.5,0.9), 17: (0.7,1.1)}
g = 9.81
x,c = stringModel(w, l, f, g)
if x is not None:
display(x, c, l.keys())
```

The dual problem is as follows:

\begin{equation} \begin{array}{ll} maximize & -\sum_{ij\in L}\ell_{ij}y_{ij} - \sum_{i\in F}f_i\circ z_i\\ s.t. & y_{ij}\geq \|v_{ij}\|,\ ij\in L \\ & \sum_{j~:~ij\in L} v_{ij}\mathrm{sgn}_{ij} + \left(\begin{array}{c}0\\ gw_i\end{array}\right) +z_i = 0, \ i=1,\ldots,n \end{array} \end{equation}

where $\mathrm{sgn}_{ij}=+1$ if $i>j$ and $-1$ otherwise and $\circ$ is the dot product. The variables are $(y_{ij},v_{ij})\in \mathbf{R}\times\mathbf{R}^2$ for $ij\in L$ and $z_i\in\mathbf{R}^2$ for $i\in F$ (we assume $z_i=0$ for $i\not\in F$).

Obviously (!) the linear constraints describe the equilibrium of forces at every mass. The ingredients are: the vectors of forces applied through adjacent strings ($v_{ij}$), gravity, and the attaching force holding a fixed point in its position. By proper use of vectorization this is much easier to express in Fusion than it looks:

In [7]:

```
def dualStringModel(w, l, f, g):
n, m = len(w), len(l)
starts = [ lKey[0] for lKey in l.keys() ]
ends = [ lKey[1] for lKey in l.keys() ]
M = Model("dual strings")
x = M.variable(Domain.inQCone(m,3)) #(y,v)
y = x.slice([0,0],[m,1])
v = x.slice([0,1],[m,3])
z = M.variable([n,2])
# z_i = 0 if i is not fixed
for i in range(n):
if i not in f:
M.constraint(z.slice([i,0], [i+1,2]), Domain.equalsTo(0.0))
B = Matrix.sparse(m, n, range(m)+range(m), starts+ends, [1.0]*m+[-1.0]*m).transpose()
w2 = Matrix.sparse(n, 2, range(n), [1]*n, [-wT*g for wT in w])
# sum(v_ij *sgn(ij)) + z_i = -(0, gw_i) for all vertices i
M.constraint(Expr.add( Expr.mul(B, v), z ), Domain.equalsTo(w2))
# Objective -l*y -fM*z
fM = Matrix.sparse(n, 2, f.keys()+f.keys(), [0]*len(f)+[1]*len(f),
[pt[0] for pt in f.values()] + [pt[1] for pt in f.values()])
M.objective(ObjectiveSense.Maximize, Expr.neg(Expr.add(Expr.dot(l.values(), y),Expr.dot(fM, z))))
M.solve()
```

Let us quickly discuss the possible situations regarding feasibility:

- The system has an equilibrium --- the problem is
**primal feasible**and**dual feasible**. - The strings are too short and it is impossible to stretch the required distance between fixed points --- the problem is
**primal infeasible**. - The system has a component that is not connected to any fixed point, hence some masses can keep falling down indefinitely, causing the problem
**primal unbounded**. Clearly the forces within such component cannot be balanced, so the problem is**dual infeasible**.

We can extend this to consider infinitely strechable springs instead of fixed-length strings connecting the masses. The next model appears in Applications of SOCP by Lobo, Boyd, Vandenberghe, Lebret. We will now interpret $\ell_{ij}$ as the base length of the spring and assume that the elastic potential energy stored in the spring at length $x$ is

$$ E_{ij}=\left\{\begin{array}{ll}0 & x\leq \ell_{ij}\\ \frac{k}{2}(x-\ell_{ij})^2 & x>\ell_{ij}\end{array}\right. $$

That leads us to consider the following second order cone program minimizing the total potential energy:

\begin{equation} \begin{array}{ll} minimize & g\cdot \sum_i w_ix_i^{(2)} + \frac{k}{2}\sum_{ij\in L} t_{ij}^2 \\ s.t. & \|x_i-x_j\|\leq \ell_{ij}+t_{ij},\ ij\in L \\ & 0\leq t_{ij},\ ij\in L \\ & x_i = f_i,\ i\in F \end{array} \end{equation}

If $t$ denotes the vector of $t_{ij}$ then using a rotated quadratic cone for $(1,T,t)$:

$$ 2\cdot 1\cdot T\geq \|t\|^2 $$

will place a bound on $\frac12\sum t_{ij}^2$. We now have a simple extension of the first model.

In [8]:

```
# w - masses of points
# l - lengths of strings
# f - coordinates of fixed points
# g - gravitational constant
# k - stiffness coefficient
def elasticModel(w, l, f, g, k):
n, m = len(w), len(l)
starts = [ lKey[0] for lKey in l.keys() ]
ends = [ lKey[1] for lKey in l.keys() ]
M = Model("strings")
x = M.variable("x", [n, 2]) # Coordinates
t = M.variable(m, Domain.greaterThan(0.0)) # Streching
T = M.variable(1) # Upper bound
M.constraint(Expr.vstack(T, Expr.constTerm(1.0), t), Domain.inRotatedQCone())
# A is the signed incidence matrix of points and strings
A = Matrix.sparse(m, n, range(m)+range(m), starts+ends, [1.0]*m+[-1.0]*m)
# ||x_i-x_j|| <= l_{i,j} + t_{i,j}
c = M.constraint("c", Expr.hstack(Expr.add(t, Expr.constTerm(l.values())), Expr.mul(A, x)),
Domain.inQCone() )
# x_i = f_i for fixed points
for i in f:
M.constraint(x.slice([i,0], [i+1,2]), Domain.equalsTo(list(f[i])))
# sum (g w_i x_i_2) + k*T
M.objective(ObjectiveSense.Minimize,
Expr.add(Expr.mul(k,T), Expr.mul(g, Expr.dot(w, x.slice([0,1], [n,2])))))
# Solve
M.solve()
if M.getProblemStatus(SolutionType.Interior) == ProblemStatus.PrimalAndDualFeasible:
return x.level().reshape([n,2]), c.dual().reshape([m,3])
else:
return None, None
```

In [9]:

```
n = 20
w = [1.0]*n
l = {(i,i+1): 0.09 for i in range(n-1)}
l.update({(5,14): 0.3})
f = {0: (0.0,1.0), 13: (0.5,0.9), 17: (0.7,1.1)}
g = 9.81
k = 800
x, c = elasticModel(w, l, f, g, k)
if x is not None:
display(x, c, l.keys())
```

This work is licensed under a Creative Commons Attribution 4.0 International License. The **MOSEK** logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API`

are not guaranteed. For more information contact our support.