The Catenary Problem

Introduction

A chain with uniformly distributed mass hangs from the endpoints \((0,1)\) and \((1,1)\) on a 2-D plane. Gravitational force acts in the negative \(y\) direction. Our goal is to find the shape of the chain in equilibrium, which is equivalent to determining the \((x,y)\) coordinates of every point along its curve when its potential energy is minimized.

This is the famous catenary problem.

A Discrete Version

To formulate as an optimization problem, we parameterize the chain by its arc length and divide it into \(m\) discrete links. The length of each link must be no more than \(h > 0\). Since mass is uniform, the total potential energy is simply the sum of the \(y\)-coordinates. Therefore, our (discretized) problem is

\[ \begin{array}{ll} \underset{x,y}{\mbox{minimize}} & \sum_{i=1}^m y_i \\ \mbox{subject to} & x_1 = 0, \quad y_1 = 1, \quad x_m = 1, \quad y_m = 1 \\ & (x_{i+1} - x_i)^2 + (y_{i+1} - y_i)^2 \leq h^2, \quad i = 1,\ldots,m-1 \end{array} \]

with variables \(x \in {\mathbf R}^m\) and \(y \in {\mathbf R}^m\). This discretized version which has been studied by Griva and Vanderbei (2005) was suggested to us by Hans Werner Borchers.

The basic catenary problem has a well-known analytical solution (see Gelfand and Fomin (1963)) which we can easily verify with CVXR.

## Problem data
m <- 101
L <- 2
h <- L / (m - 1)

## Form objective
x <- Variable(m)
y <- Variable(m)
objective <- Minimize(sum(y))

## Form constraints
constraints <- list(x[1] == 0, y[1] == 1,
                    x[m] == 1, y[m] == 1,
                    diff(x)^2 + diff(y)^2 <= h^2)

## Solve the catenary problem
prob <- Problem(objective, constraints)
result <- solve(prob)

We can now plot it and compare it with the ideal solution. Below we use alpha blending and differing line thickness to show the ideal in red and the computed solution in blue.

xs <- result$getValue(x)
ys <- result$getValue(y)

catenary <- ggplot(data.frame(x = xs, y = ys)) +
    geom_line(mapping = aes(x = x, y = y), color = "blue", size = 1) +
    geom_point(data = data.frame(x = c(xs[1], ys[1]), y = c(xs[m], ys[m])),
               mapping = aes(x = x, y = y), color = "red") 

ideal <- function(x) { 0.22964 *cosh((x -0.5) / 0.22964) - 0.02603 }

catenary + stat_function(fun = ideal , colour = "brown", alpha = 0.5, size = 3)

Figure 1: Analytic (red) and computed solution (blue) to the catenary problem

Additional Ground Constraints

A more interesting situation arises when the ground is not flat. Let \(g \in {\mathbf R}^m\) be the elevation vector (relative to the \(x\)-axis), and suppose the right endpoint of our chain has been lowered by \(\Delta y_m = 0.5\). The analytical solution in this case would be difficult to calculate. However, we need only add two lines to our constraint definition,

constr[[4]] <- (y[m] == 0.5)
constr <- c(constr, y >= g)

to obtain the new result.

Below, we define \(g\) as a staircase function and solve the problem.

## Lower right endpoint and add staircase structure
ground <- sapply(seq(0, 1, length.out = m), function(x) {
    if(x < 0.2)
        return(0.6)
    else if(x >= 0.2 && x < 0.4)
        return(0.4)
    else if(x >= 0.4 && x < 0.6)
        return(0.2)
    else
        return(0)
})
constraints <- c(constraints, y >= ground)
constraints[[4]] <- (y[m] == 0.5)
prob <- Problem(objective, constraints)
result <- solve(prob)

to obtain the new result.

The figure below shows the solution of this modified catenary problem for \(m = 101\) and \(h = 0.04\). The chain is shown hanging in blue, bounded below by the red staircase structure, which represents the ground.

xs <- result$getValue(x)
ys <- result$getValue(y)

ggplot(data.frame(x = xs, y = ys)) +
    geom_line(mapping = aes(x = x, y = y), color = "blue") +
    geom_point(data = data.frame(x = c(xs[1], ys[1]), y = c(xs[m], ys[m])),
               mapping = aes(x = x, y = y), color = "red") +
    geom_line(data.frame(x = xs, y = ground),
              mapping = aes(x = x, y = y), color = "brown")

Figure 2: Asymmetric catenary problem with ground constraints.

Session Info

sessionInfo()

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin21.6.0 (64-bit)
## Running under: macOS Ventura 13.0
## 
## Matrix products: default
## BLAS:   /usr/local/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
## LAPACK: /usr/local/Cellar/r/4.2.1_4/lib/R/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
## [1] ggplot2_3.3.6 CVXR_1.0-11  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0 xfun_0.34        bslib_0.4.0      slam_0.1-50     
##  [5] lattice_0.20-45  Rmosek_10.0.25   colorspace_2.0-3 vctrs_0.5.0     
##  [9] generics_0.1.3   htmltools_0.5.3  yaml_2.3.6       gmp_0.6-6       
## [13] utf8_1.2.2       rlang_1.0.6      jquerylib_0.1.4  pillar_1.8.1    
## [17] glue_1.6.2       Rmpfr_0.8-9      withr_2.5.0      DBI_1.1.3       
## [21] Rcplex_0.3-5     bit64_4.0.5      lifecycle_1.0.3  stringr_1.4.1   
## [25] munsell_0.5.0    blogdown_1.13    gtable_0.3.1     gurobi_9.5-2    
## [29] codetools_0.2-18 evaluate_0.17    labeling_0.4.2   knitr_1.40      
## [33] fastmap_1.1.0    fansi_1.0.3      cccp_0.2-9       highr_0.9       
## [37] Rcpp_1.0.9       scales_1.2.1     cachem_1.0.6     jsonlite_1.8.3  
## [41] farver_2.1.1     bit_4.0.4        digest_0.6.30    stringi_1.7.8   
## [45] bookdown_0.29    dplyr_1.0.10     Rglpk_0.6-4      grid_4.2.1      
## [49] ECOSolveR_0.5.4  cli_3.4.1        tools_4.2.1      magrittr_2.0.3  
## [53] sass_0.4.2       tibble_3.1.8     pkgconfig_2.0.3  rcbc_0.1.0.9001 
## [57] Matrix_1.5-1     assertthat_0.2.1 rmarkdown_2.17   R6_2.5.1        
## [61] compiler_4.2.1

Source

R Markdown