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") ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.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
## Testthat Results: No output is goodAdditional 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.
## Testthat Results: No output is goodSession Info
sessionInfo()## R version 4.4.1 (2024-06-14)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sonoma 14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] ggplot2_3.5.1 CVXR_1.0-15 testthat_3.2.1.1 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] gmp_0.7-4 clarabel_0.9.0 sass_0.4.9 utf8_1.2.4
## [5] generics_0.1.3 slam_0.1-50 blogdown_1.19 lattice_0.22-6
## [9] digest_0.6.36 magrittr_2.0.3 evaluate_0.24.0 grid_4.4.1
## [13] bookdown_0.40 pkgload_1.4.0 fastmap_1.2.0 rprojroot_2.0.4
## [17] jsonlite_1.8.8 Matrix_1.7-0 ECOSolveR_0.5.5 brio_1.1.5
## [21] Rmosek_10.2.0 fansi_1.0.6 scales_1.3.0 codetools_0.2-20
## [25] jquerylib_0.1.4 cli_3.6.3 Rmpfr_0.9-5 crayon_1.5.3
## [29] rlang_1.1.4 Rglpk_0.6-5.1 bit64_4.0.5 munsell_0.5.1
## [33] withr_3.0.0 cachem_1.1.0 yaml_2.3.9 tools_4.4.1
## [37] Rcplex_0.3-6 rcbc_0.1.0.9001 dplyr_1.1.4 colorspace_2.1-0
## [41] gurobi_11.0-0 assertthat_0.2.1 vctrs_0.6.5 R6_2.5.1
## [45] lifecycle_1.0.4 bit_4.0.5 desc_1.4.3 cccp_0.3-1
## [49] pkgconfig_2.0.3 bslib_0.7.0 pillar_1.9.0 gtable_0.3.5
## [53] glue_1.7.0 Rcpp_1.0.12 highr_0.11 xfun_0.45
## [57] tibble_3.2.1 tidyselect_1.2.1 knitr_1.48 farver_2.1.2
## [61] htmltools_0.5.8.1 labeling_0.4.3 rmarkdown_2.27 compiler_4.4.1