Using nonlinear programming for solving mechanical engineering problems: Designing a column subjected to a buckling load

Nonlinear programming is used for solving optimization problems. The goal of such optimization problems is to minimize (or maximize) some objective function. However, these optimization problems often specify a number of constraints. Nonlinear programming is applied when the objective function or some of the specified constraints are nonlinear.

The following R code demonstrates how to solve an optimization problem using nonlinear programming. More specifically, the R code uses the Augmented Lagrange Multiplier (ALM) method for finding the minimum design weight of a circular thin walled column subjected to a buckling load.

An introduction to nonlinear programming (and the ALM method) can be found in Rao’s 2009 book Engineering Optimization: theory and practice (4th. Ed., Chapter 7).
The minimum design weight problem of the circular thin walled column is described in Arora’s 2012 book Introduction to optimum design (3rd. Ed., pp. 40-42).

Do you have any suggestions and/or questions? Please contact Stefan Gelissen (email: info at
See this page for an overview of all of Stefan’s R code blog posts.


Finding the minimum design weight of a circular thin walled column
The goal of this optimization problem is to find the minimum mass of a circular thin walled column of length L. This tube has to support a specified axial load P without buckling, and without the stress in the column exceeding the allowable stress. The column and its axial load are depicted in the figure below.

nonlinear programming buckling

The optimization problem can be formulated in the following way:

Minimize the nonlinear objective function:

mass = ρ*(L*A) = ρ*L*2*R*t ,
where ρ=material mass density (in kg/m3), L=length of column (in m), A=cross sectional area (in m2), R=mean radius of the column (in m), and t=thickness of the wall of the column (in m).

Subject to the following nonlinear constraints:

  1. The stress in the column must not exceed the allowable stress, that is:
    P/A ≤ σallowable
    where P=area force (in N), A=cross sectional area of column (in m2), σallowable=material allowable stress (in N/m2)
  2. The applied axial load must not exceed the Euler (or critical) load, that is:
    P ≤ (π2*E*I)/(4*L2)
    where P=axial load (in N), E=Young’s modulus (in Pa), I=moment of inertia of the column (in m4), L=length of column (in m)
  3. Prevent local buckling of the thin walled column, that is:
    R/t ≤ k
    where R=mean radius of the column (in m), t=wall thickness of the column (in m), k=a constant.

Finally, the design variables are bounded:

  1. R values: Rmin ≤ R ≤ Rmax
  2. t values: tmin ≤ t ≤ tmax

R code

#Minimum design weight of circular thin walled column
#Arora, pp.40-42, and exercise 3.23, p.85


#initialize variables
rho <- 7850 #density material (kg/m^3)
sigmaAllow <- 250e+6 #allowable stress (Pa)
E <- 210e+9 #Young's modulus (Pa)
l <- 5 #length of column (m)
P <- 50e+3 #axial load (N)

#objective function: mass = rho*(l*A) = 2*rho*l*pi*R*t
objFun <- function(x) {
  R <- x[1] #mean radius
  t <- x[2] #wall thickness

constrFun <- function(x) {
  #initialize variables
  R <- x[1] #mean radius
  t <- x[2] #wall thickness
  A <- 2*pi*R*t #material cross-sectional area
  I <- pi*R^3*t #moment of inertia (assuming R>>t)
  #define constraints
  z1 <- sigmaAllow - P/A #allowable stress
  z2 <- (pi^2*E*I)/(4*l^2) - P #buckling load
  z3 <- 50 - R/t #local buckling
  return(c(z1, z2, z3))}

#Augmented Lagrange Multiplier (ALM) method
#note: for this ALM method to converge to the optimal design values
#the start values (starting design values) have to be in the infeasible region
#thus, first check infeasibility of the start values
constrFun(c(2e-2, 2e-2))
#apply ALM method
opt <- auglag(c(2e-2, 2e-2), fn=objFun, hin=constrFun, localsolver="LBFGS",
              lower=c(.01, 5e-3), upper=c(1,200e-3)) #.01<=R<=1 m; 5<=t<=200 mm
#check convergence: successful (integer>0) or failure (integer<0)
#optimal solution
#(this solution is identical to the solution computed by Arora: R=53.6mm; t=5.0mm)
#function value (is again identical to Arora's computed value of 66kg)

#graphical representation of the optimization problem and its solution

#contour plot
x1values <- seq(0,.06,length.out=75)
x2values <- seq(0,.03,length.out=75)
lbgrid <- expand.grid(x1=x1values, x2=x2values)
lbgrid$cl <- apply(lbgrid, 1, objFun)
cl.matrix <- matrix(lbgrid$cl, nrow=length(x1values))
#plot contours
contour(x1values, x2values, cl.matrix, xlab="R", ylab="t",

#include constraints:
#lower bound of R
abline(v=.01, lty=2, col="blue", lwd=2)
#lower bound of t
abline(h=5e-3, lty=2, col="blue", lwd=2)
#allowable stress
constr1 <- function (x) {P/(sigmaAllow*2*pi*x)}
lines(x1values,constr1(x1values), lty=2, col="blue")
#buckling load
constr2 <- function (x) {(P*4*l^2)/(pi^3*E*x^3)}
lines(x1values,constr2(x1values), lty=2, col="blue")
#local buckling
constr3 <- function (x) {50/x}
lines(x1values,constr3(x1values), lty=2, col="blue")

#mark feasible region (FR)
text(.05,.020, "FR")

#include start value
points(2e-2, 2e-2, col="blue")
#include optimal (final) solution
points(opt$par[1],opt$par[2], col="red", lwd=2)

Next blog post?
In my next blog post I will again employ the ALM method and solve another engineering optimization problem in R.