Applying SciPy.optimize to form finding
Previously, in our post “The force density method”, we described how to use the force density method to compute an equilibrium shape of the provided network data. In this post we will show how to apply SciPy’s nonlinear solvers to the same problem.
To solve the system of linear equations, Cholesky decomposition was used. The latter part of the code was:
from numpy import array, zeros, float64 p = zeros((n, 3), dtype=float64) q = array([1.0 for i in range(m)]).reshape((-1, 1)) Q = diags([q.flatten()], [0]) A = Cit.dot(Q).dot(Ci) b = p[free] - Cit.dot(Q.dot(Cf.dot(xyz[fixed]))) xyz[free] = cho_solve(cho_factor(A), b)
Imagine for a moment that the force densities Q are nonlinear, which happens when they need to be recalculated from the resulting geometry. In that case, we continuously update LHS matrix A and RHS vector b until the (norm of the) residual forces Ax-b are sufficiently small (within a chosen tolerance). This is equivalent to applying Newton’s or Newton-Rhapson’s method.
SciPy, a scientific library for Python, offers many solvers for nonlinear problems, typically more advanced than Newton’s method. When looking at the online documentation we find scipy.optimize.minimize for the “minimization of [a] scalar function of one or more variables” and scipy.optimize.root in order to “find a root of a vector function”.
We can actually apply both types. The latter can be used by formulating our problem, i.e. finding static equilibrium, as the following vector-valued function:
\begin{align}
f'(x_i) &= Ax_i-b \\
&= C_i^T Q C_i x_i – ( p – C_i^T Q C_f x_f )\\
&= C_i^T Q C x – p
\end{align}
With it’s derivative:
\begin{align}
f”(x_i) &= A \\
&= C_i^T Q C_i
\end{align}
Then we can solve for the root of this vector-valued function using SciPy. Note that for SciPy we have to reshape some of our variables (defined as a set of three vectors; one for each dimension) to be vectors (stacking the three vectors on top of each other).
from scipy import optimize
from scipy.sparse import block_diag
xi = xyz[free]
ni = len(free)
def fprime(xi):
A = Cit.dot(Q).dot(Ci)
b = p[free] - Cit.dot(Q.dot(Cf.dot(xyz[fixed])))
residuals = A.dot(xi.reshape(3,ni).T)-b
return residuals.T.reshape(3*ni,)
def fprime2(_):
A = Cit.dot(Q).dot(Ci)
return block_diag([A,A,A]).todense()
x0 = xi.T.reshape(3*ni,)
sol = optimize.root(fprime, x0, jac=fprime2, method='lm')
xyz[free] = sol.x.reshape(3,ni).T
For the former approach, minimization, we need to define a scalar-valued function to minimize. This is done by defining the energy within our cable net, which we also want to approach zero. Notice that the internal forces f'(x) are the derivative of energy f(x) with respect to x.
\begin{align}
f(x_i) &= \frac{1}{2}x_i^TAx_i-b^Tx_i \\
&= \frac{1}{2}x_i^TC_i^T Q C_i x_i – ( p – C_i^T Q C_f x_f )^Tx_i\\
\end{align}
def fun(xi):
A = block_diag([A,A,A])
b = (p[free] - Cit.dot(Q.dot(Cf.dot(xyz[fixed])))).T.reshape(3*ni,1)
return (0.5*xi.T*A+b.T).dot(xi)
sol = optimize.fmin_cg(fun,x0,fprime)
xyz[free] = sol.reshape(3,ni).T
Both approaches yield the same answer as our original code. However, neither approach is faster; in fact they are orders slower. Once the problem is truly nonlinear i.e. Q is no longer constant, will their use potentially pay off.