Using nonlinear programming for solving mechanical engineering problems: Designing a two-bar bracket

In my previous blog posts I demonstrated how to use nonlinear programming in R for solving an optimization problem. To be more specific, I used the Augmented Lagrange Multiplier (ALM) method for solving a specific engineering optimization problem.

In this blog post I will again employ the ALM method for solving an engineering optimization problem. This time I will analyze a two-bar bracket design subjected to a load. The goal of this optimization problem is to find the minimal mass of the design so that the bracket will support the load without structural failure.

This two-bar bracket design problem is described and formulated in Arora’s 2012 book Introduction to optimum design (3rd. Ed., pp. 30-36).

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

 

Finding the minimum design weight of a two-bar bracket
The goal of this optimization problem is to find the minimum mass of a two-bar bracket subjected to a load W. The load is applied under angle θ (between 0 and 90 degrees). The two-bar bracket and the applied load are depicted in the following figure (with on the right-hand side a free-body diagram of the bracket at joint c).

nonlinear programming bracket

From the above figure we see that sin(α)=(S/2)/L, cos(α)=H/L, and L=\sqrt{H^{2} + (\frac{S}{2})^2}

Using the free-body diagram we can state the following equilibrium conditions:

  1. The sum of the forces in the horizontal direction:
    -F1*sin(α)+F2*sin(α)=W*cos(θ)
  2. The sum of the forces is the vertical direction:
    -F1*cos(α)-F2*cos(α)=W*sin(θ)

Solving these two equilibrium equations simultaneously yields:

F1=-0.5*W*L*( sin(θ)/H + 2*cos(θ)/S )
F2=-0.5*W*L*( sin(θ)/H – 2*cos(θ)/S )

From the free-body diagram we can see that vectors F1 and F2 represent tensile forces. In case the bars are subjected to tension then we assume that these force vectors take positive values. However, should a bar be in compression then the force vector will take a negative value (that is, in reality the force vector points in the opposite direction).

For this design problem we further assume that the two bars are made from the same material. Additionally, both bars are circular tubes (see figure below) having identical outer and inner diameters.

nonlinear programming tube

Note that the cross-sectional area of this circular tube is equal to π/4*(Do2-Di2). Based on this cross-sectional area the stress in each bar is:

  1. σ1 = F1/A1
    where σ1 is the stress in bar 1 (in Pa), F1 is the force in bar 1 (in N), and A1 the cross-sectional area (in m2)
  2. σ2 = F2/A2
    where σ2 is the stress in bar 2 (in Pa), F2 is the force in bar 2 (in N), and A2 the cross-sectional area (in m2)

We are now ready to state the objective function and the constraints for this optimization problem.

Minimize the nonlinear objective function:

mass = ρ*(L*(A1+A2)) = \rho*\sqrt{H^2 + (\frac{S}{2})^2}*\frac{\pi}{4}*(D_{o1}^2-D_{i1}^2 + D_{o2}^2-D_{i2}^2)
where ρ=material mass density (in kg/m3), H= height of the bracket (in m), S=span of the bracket (in m), A=cross sectional area of a bar (in m2), Do=outer diameter of a bar (in m), and Di=inner diameter of a bar (in m).

Subject to the following nonlinear constraint:

The stress in each bar must not exceed the allowable stress, that is:
F/A ≤ σallowable
where F=area force (in N), A=cross sectional area of a bar (in m2), σallowable=material allowable stress (in N/m2).

Remember that the force in a bar is either tensile (positive value) or compressive (negative value).
Therefore, when a bar is under tension (i.e., F is positive) we employ:
F/A ≤ σallowable
On the other hand, when a bar is under compression we substitute a negative value for F in:
-F/A ≤ σallowable

Finally, the design variables are bounded:

  1. Outer diameter values: (Do)min ≤ Do ≤ (Do)max
  2. Inner diameter values: (Di)min ≤ Di ≤ (Di)max

R code

#Minimum design weight of two-bar bracket
#Arora, pp.30-36, and exercise 3.26, p.86

library(nloptr)

#initialize variables
rho <- 7850 #density material (kg/m^3)
sigmaAllow <- 250e+6 #allowable stress (Pa)
E <- 210e+9 #Young's modulus (Pa)
h <- 1 #height of bracket (m)
s <- 1.5 #base of bracket (m)
W <- 10e+3 #load (N)
theta <- 30 #angle of load (degrees)
thetaRad <- 30*pi/180 #angle of load (rad)


#objective function (minimize mass of two-bar bracket)
objFun <- function(x) {
  
  #initialize variables
  #(assume that the two bars are identical circular tubes)
  Do <- x[1] #outer diameter
  Di <- x[2] #inner diameter
  
  #objective function
  ((pi*rho)/4)*sqrt(h^2+(0.5*s)^2)*(Do^2-Di^2+Do^2-Di^2)}


#constraints
constrFun <- function(x) {
  
  #initialize variables
  #(assume that the two bars are identical circular tubes)
  Do <- x[1] #outer diameter
  Di <- x[2] #inner diameter
  
  #stress in bar 1
  bar1stress <- (2*W*sqrt(h^2+(0.5*s)^2)/(pi*(Do^2-Di^2)))*
    (sin(thetaRad)/h + 2*cos(thetaRad)/s)
  #stress in bar 2
  bar2stress <- (2*W*sqrt(h^2+(0.5*s)^2)/(pi*(Do^2-Di^2)))*
    (sin(thetaRad)/h - 2*cos(thetaRad)/s)
  
  #define constraints
  #bar 1:
  if (bar1stress<0) {z1 <- sigmaAllow + bar1stress} #compression
  else if (bar1stress>=0) {z1 <- sigmaAllow - bar1stress} #tension
  #bar 2:
  if (bar2stress<0) {z2 <- sigmaAllow + bar2stress} #compression
  else if (bar2stress>=0) {z2 <- sigmaAllow - bar2stress} #tension
  
  #outer diameter > inner diameter
  z3 <- Do - Di
  
  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(.002,.001)) #infeasible
#apply ALM method
opt1 <- auglag(c(2e-3,1e-3), fn=objFun, hin=constrFun, localsolver="LBFGS",
               lower=c(0,0), upper=c(1,.9)) #0<=Do<=1 m; 0<=Di<=.9 m
#check convergence: successful (integer>0) or failure (integer<0)
opt1$convergence
#optimal solution
opt1$par
#function value (is identical to the minimum weight reported by Arora: 0.812kg)
opt1$value

constrFun(c(.02,0.001)) #feasible
opt2 <- auglag(c(2e-2,1e-3), fn=objFun, hin=constrFun, localsolver="LBFGS",
               lower=c(0,0), upper=c(1,.9))
opt2$convergence
opt2$value


#the "derivative free approach" COBYLA seems to be less sensitive to the
#feasibility/infeasility of the start values for this optimization problem
constrFun(c(.002,.001)) #infeasible
opt3 <- auglag(c(2e-2,1e-3), fn=objFun, hin=constrFun, localsolver="COBYLA",
               lower=c(0,0), upper=c(1,.9))
opt3$convergence
opt3$value

constrFun(c(.02,0.001)) #feasible
opt4 <- auglag(c(2e-2,1e-3), fn=objFun, hin=constrFun, localsolver="COBYLA",
               lower=c(0,0), upper=c(1,.9))
opt4$convergence
opt4$value



#graphical representation of the optimization problem and its solution

#contour plot objective function
x1values <- seq(0.001,0.02,length.out=75)
x2values <- seq(0,0.02,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="outer", ylab="inner",
        levels=c(.5, .82))


#include constraints:

#outer diameter > inner diameter
abline(a=0, b=1, lty=2, col="blue", lwd=2)

#allowable stress bar 1 (note: for theta=30 the force in bar 1 is compressive)
-0.5*W*sqrt(h^2+(0.5*s)^2)*(sin(thetaRad)/h + 2*cos(thetaRad)/s) #force in bar 1
constr1 <- function (x) {sqrt(x^2-
                                ((sin(thetaRad)/h + 2*cos(thetaRad)/s)*(2*W*sqrt(h^2+(0.5*s)^2)))/
                                (pi*sigmaAllow))}
lines(x1values,constr1(x1values), lty=2, lwd=2, col="red")

#allowable stress bar 2 (note: for theta=30 the force in bar 2 is tensile)
-0.5*W*sqrt(h^2+(0.5*s)^2)*(sin(thetaRad)/h - 2*cos(thetaRad)/s) #force in bar 2
constr2 <- function (x) {sqrt( (((sin(thetaRad)/h - 2*cos(thetaRad)/s)*(2*W*sqrt(h^2+(0.5*s)^2)))/
                                  (pi*sigmaAllow)) + x^2)}
lines(x1values,constr2(x1values), lty=2, lwd=2, col="red")

#mark feasible region (FR)
text(.015,.005, "FR")

#include start value
points(.002,.001, col="blue")

#include optimal (final) solution
points(opt1$par[1],opt1$par[2], col="red", lwd=2)
#finally, note that the line of the stress constraint for bar 1
#coincides exactly with the objective function having value 0.182kg
#this implies that there are actually an infinite number of
#optimum design values (or optimum solutions)
#a similar conclusion was reached by Arora