1 Introduction

I was thinking about a good test case to evaluate different versions of auglag, say with optim, nlminb, or ucminf. I found the catenary curve to be a good one. The aim is to solve the “hanging chain” problem not as a differential equation, but as an optimization problem.

1.1 The theoretical solution

First we reproduce the exact solution for a chain of length \(L = 2\) in the \((x,y)\)-plane, from \((0,0)\) to \((1,0)\). The theoretical solution is \[ y = a\,\cosh(\frac{x}{a}) \] where the distance \(w\) between the fixed points, the length \(L\) of the chain, and the curvature radius \(a\) are related to each other through \[ L = 2 a\,\sinh(\frac{w}{2a}) \]

We calculate the parameter \(a\) in R as a solution of this as a function in \(a\):

w <- 1; L <- 2
fun <- function(x) L - 2*x*sinh(w/(2*x))
aa <- uniroot(fun, c(0.1, 0.5))$root        # 0.22965

The chain curve is symmetric around \(0.5\) and the fixed points shall have a height of 1.0 above the ground.

bb <- aa * cosh(0.5/aa)                        # 0.02597928
cy <- function(x) aa * cosh((x - 0.5)/aa) - bb + 1.0

1.2 The catenary as an optimization problem

So we are assuming a chain of length 2 is bound to the points \(P1 = (0,1)\) and \(P2 = (1,1)\) with gravitational force acting in negative y-direction.

N <- 101        # no. of points
L <- 2          # total length of chain
h <- L / (N-1)  # maximal length of each chain link

The decision variables are the x- and y-coordinates of the beginning and end of the chain elements. We will look at the cases of 101 points, that is 100 chain links.

The parameter vector c(x, y) has dimension \(2 N\), the concatenated vector of x- and y-coordinates. The objective function is simply to minimize the potential energy.

fnobj <- function(p)  sum(p[(N+1):(2*N)])  # sum(y)

grobj <- function(p)  c(rep(0, N), rep(1, N))

We will use the exaxt gradient, as it is so easy to write down.

We have two kinds of constraints, equality and inequality constraints. The equality constraints fix the beginning and end of the chain to the points \(P1\) an¡“¡“d \(P2\). Therefore:

heq <- function(p) {
    c(p[1], p[N] - 1, p[N+1] - 1, p[2*N] - 1)
}

The inequality constraints fix the length of the individual chain links to be \(h\) maximally, so the total length is \(L = (N-1) h = 2\).

hin <- function(p) {
    x <- p[1:N]; y <- p[(N+1):(2*N)]
    h^2 - diff(x)^2 - diff(y)^2
}

The starting configuration will be with all links on the zero level.

x0 <- seq(0, 1, length.out=N)
p0 <- c(x0, rep(0, N))

2 The ‘augmented Lagrangian’ approach

2.1 auglag() with ‘optim’

First we will solve it with alabama::auglag. There are two choices for the inner solver, optim with method BFGS, or nlminb (both in base R).

system.time(
sol1 <- alabama::auglag(p0, fnobj, grobj, hin=hin, heq=heq,
                        control.outer=list(trace=FALSE, eps=1e-08),
                        control.optim=list(maxit=500, reltol=1e-08))
)
##    user  system elapsed 
##  29.933   4.971  35.364

This is not a good approximation of the true catenary line, as we will see in a minute. Changing the tolerances will not help, but I had to increase the maximum number of iterations, the default value 100 is much too small.

optim uses its own calculation of the ‘Jacobians’. Instead we will apply the finite-difference Jacobians of the pracma package here (using the Jacobians from the numDeriv package will slow down the calculations too much).

library(pracma)
heq.jac <- function(p) pracma::jacobian(heq, p)
hin.jac <- function(p) pracma::jacobian(hin, p)

system.time(
sol2 <- alabama::auglag(p0, fnobj, grobj, hin=hin, heq=heq,
                        hin.jac = hin.jac, heq.jac = heq.jac,
                        control.outer=list(trace=FALSE, eps=1e-08),
                        control.optim=list(maxit=500, reltol=1e-08))
)
##    user  system elapsed 
##  45.937   7.288  53.352

Let’s plot the generated curves to check plausibility for the moment.

The red curve is the solution generated by optim with method ‘BFGS’ and Jacobians calculated by the inner solver. What is striking is that the solution is not symmetric as it should be. optim still returns sol1$convergence=0, that is successful convergence.

The blue curve uses Jacobians calculated with higher precision and lies almost exactly on top of the exact solution (in gray). The maximal distance between this curve and the true catenary line is \(1.5\cdot 10^{-4}\) and the RMS error is \(7.3\cdot10^{-5}\).

xx <- x2[1:N]; yy <- x2[(N+1):(2*N)]
sqrt(sum((yy - cy(xx))^2) / (N+1))
## [1] 7.300808e-05
max(abs(yy - cy(xx)))
## [1] 0.0001544463

2.2 auglag() with ‘nlminb’

Calling auglag with nlminb as inner solver is a viable alternative. The control parameters for the inner solver are different, though, especially the step.min has been found by trial and error.

system.time(
sol3 <- alabama::auglag(p0, fnobj, grobj, hin=hin, heq=heq,
                    hin.jac = hin.jac, heq.jac = heq.jac,
                    control.outer=list(trace=FALSE, method="nlminb"),
                    control.optim=list(iter.max=500, rel.tol=1e-10, 
                                       step.min=0.5))
)
##    user  system elapsed 
##  56.086   8.659  64.965

The generated curve looks like this: