Package 'NlcOptim'

Solve Optimization problem with Nonlinear Objective and Constraints

Description

Sequential Quatratic Programming (SQP) method is implemented to find solution for general nonlinear optimization problem (with nonlinear objective and constraint functions). The SQP method can be find in detail in Chapter 18 of Jorge Nocedal and Stephen J. Wright's book. Linear or nonlinear equality and inequality constraints are allowed. It accepts the input parameters as a constrained matrix. The function solnl is to solve generalized nonlinear optimization problem:

$min f(x)$

$s.t. ceq(x)=0$

$c(x)\le 0$

$Ax\le B$

$Aeq x \le Beq$

$lb\le x \le ub$

Usage

solnl(X = NULL, objfun = NULL, confun = NULL, A = NULL, B = NULL,
  Aeq = NULL, Beq = NULL, lb = NULL, ub = NULL, tolX = 1e-05,
  tolFun = 1e-06, tolCon = 1e-06, maxnFun = 1e+07, maxIter = 4000)

Arguments

X	Starting vector of parameter values.
objfun	Nonlinear objective function that is to be optimized.
confun	Nonlinear constraint function. Return a ceq vector and a c vector as nonlinear equality constraints and an inequality constraints.
A	A in the linear inequality constraints.
B	B in the linear inequality constraints.
Aeq	Aeq in the linear equality constraints.
Beq	Beq in the linear equality constraints.
lb	Lower bounds of parameters.
ub	Upper bounds of parameters.
tolX	The tolerance in X.
tolFun	The tolerance in the objective function.
tolCon	The tolenrance in the constraint function.
maxnFun	Maximum updates in the objective function.
maxIter	Maximum iteration.

Value

Return a list with the following components:

par	The optimum solution.
fn	The value of the objective function at the optimal point.
counts	Number of function evaluations, and number of gradient evaluations.
lambda	Lagrangian multiplier.
grad	The gradient of the objective function at the optimal point.
hessian	Hessian of the objective function at the optimal point.

Author(s)

Xianyan Chen, Xiangrong Yin

References

Nocedal, Jorge, and Stephen Wright. Numerical optimization. Springer Science & Business Media, 2006.

Examples

library(MASS)
###ex1
objfun=function(x){
 return(exp(x[1]*x[2]*x[3]*x[4]*x[5]))  
}
#constraint function
confun=function(x){
 f=NULL
 f=rbind(f,x[1]^2+x[2]^2+x[3]^2+x[4]^2+x[5]^2-10)
 f=rbind(f,x[2]*x[3]-5*x[4]*x[5])
 f=rbind(f,x[1]^3+x[2]^3+1)
 return(list(ceq=f,c=NULL))
}

x0=c(-2,2,2,-1,-1)
solnl(x0,objfun=objfun,confun=confun)

####ex2
obj=function(x){
 return((x[1]-1)^2+(x[1]-x[2])^2+(x[2]-x[3])^3+(x[3]-x[4])^4+(x[4]-x[5])^4)  
}
#constraint function
con=function(x){
 f=NULL
 f=rbind(f,x[1]+x[2]^2+x[3]^3-2-3*sqrt(2))
 f=rbind(f,x[2]-x[3]^2+x[4]+2-2*sqrt(2))
 f=rbind(f,x[1]*x[5]-2)
 return(list(ceq=f,c=NULL))
}

x0=c(1,1,1,1,1)
solnl(x0,objfun=obj,confun=con)


##########ex3
obj=function(x){
 return((1-x[1])^2+(x[2]-x[1]^2)^2)  
}
#constraint function
con=function(x){
 f=NULL
 f=rbind(f,x[1]^2+x[2]^2-1.5)
 return(list(ceq=NULL,c=f))
}

x0=as.matrix(c(-1.9,2))
obj(x0)
con(x0)
solnl(x0,objfun=obj,confun=con)


##########ex4
objfun=function(x){
 return(x[1]^2+x[2]^2)  
}
#constraint function
confun=function(x){
 f=NULL
 f=rbind(f,-x[1] - x[2] + 1)
 f=rbind(f,-x[1]^2 - x[2]^2 + 1)
 f=rbind(f,-9*x[1]^2 - x[2]^2 + 9)
 f=rbind(f,-x[1]^2 + x[2])
 f=rbind(f,-x[2]^2 + x[1])
 return(list(ceq=NULL,c=f))
}

x0=as.matrix(c(3,1))
solnl(x0,objfun=objfun,confun=confun)


##############ex5
rosbkext.f <- function(x){
   n <- length(x)
   sum (100*(x[1:(n-1)]^2 - x[2:n])^2 + (x[1:(n-1)] - 1)^2)
}
n <- 2
set.seed(54321)
p0 <- rnorm(n)
Aeq <- matrix(rep(1, n), nrow=1)
Beq <- 1
lb <- c(rep(-Inf, n-1), 0)
solnl(X=p0,objfun=rosbkext.f, lb=lb, Aeq=Aeq, Beq=Beq)
ub <- rep(1, n)
solnl(X=p0,objfun=rosbkext.f, lb=lb, ub=ub, Aeq=Aeq, Beq=Beq)


##############ex6
nh <- vector("numeric", length = 5)

Nh <- c(6221,11738,4333,22809,5467)
ch <- c(120, 80, 80, 90, 150)

mh.rev <- c(85, 11, 23, 17, 126)
Sh.rev <- c(170.0, 8.8, 23.0, 25.5, 315.0)

mh.emp <- c(511, 21, 70, 32, 157)
Sh.emp <- c(255.50, 5.25, 35.00, 32.00, 471.00)

ph.rsch <- c(0.8, 0.2, 0.5, 0.3, 0.9)

ph.offsh <- c(0.06, 0.03, 0.03, 0.21, 0.77)

budget = 300000
n.min <- 100
relvar.rev <- function(nh){
 rv <- sum(Nh * (Nh/nh - 1)*Sh.rev^2)
 tot <- sum(Nh * mh.rev)
 rv/tot^2
}

relvar.emp <- function(nh){
 rv <- sum(Nh * (Nh/nh - 1)*Sh.emp^2)
 tot <- sum(Nh * mh.emp)
 rv/tot^2
}

relvar.rsch <- function(nh){
 rv <- sum( Nh * (Nh/nh - 1)*ph.rsch*(1-ph.rsch)*Nh/(Nh-1) )
 tot <- sum(Nh * ph.rsch)
 rv/tot^2
}

relvar.offsh <- function(nh){
 rv <- sum( Nh * (Nh/nh - 1)*ph.offsh*(1-ph.offsh)*Nh/(Nh-1) )
 tot <- sum(Nh * ph.offsh)
 rv/tot^2
}

nlc.constraints <- function(nh){
 h <- rep(NA, 13)
 h[1:length(nh)] <- (Nh + 0.01) - nh
 h[(length(nh)+1) : (2*length(nh)) ] <- (nh + 0.01) - n.min
 h[2*length(nh) + 1] <- 0.05^2 - relvar.emp(nh)
 h[2*length(nh) + 2] <- 0.03^2 - relvar.rsch(nh)
 h[2*length(nh) + 3] <- 0.03^2 - relvar.offsh(nh)
 return(list(ceq=NULL, c=-h))
}

nlc <- function(nh){
 h <- rep(NA, 3)
 h[ 1] <- 0.05^2 - relvar.emp(nh)
 h[ 2] <- 0.03^2 - relvar.rsch(nh)
 h[3] <- 0.03^2 - relvar.offsh(nh)
 return(list(ceq=NULL, c=-h))
}

Aeq <- matrix(ch/budget, nrow=1)
Beq <- 1

A=rbind(diag(-1,5,5),diag(1,5,5))
B=c(-Nh-0.01,rep(n.min-0.01,5))

solnl(X=rep(100,5),objfun=relvar.rev,confun=nlc.constraints, Aeq=Aeq, Beq=Beq)

solnl(X=rep(100,5),objfun=relvar.rev,confun=nlc, Aeq=Aeq, Beq=Beq, A=-A, B=-B)

---