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.

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
```

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))
```

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`

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: