optimization - How to conduct optimisation in R when there are constraints -

i need optimise following object function. qgpd package fextremes.

var.asym <- function(alpha1, alpha2, xi, beta, n){   term11 <- alpha1*(1-alpha1)^(2*xi-1)   term12 <- alpha1*(1-alpha1)^(xi-1)*(1-alpha2)^xi   term22 <- alpha2*(1-alpha2)^(2*xi-1)   sigma <- matrix(c(term11, term12, term12, term22), nrow=2, byrow=true)   sigma*beta^2/n }  mop.jacob.inv <- function(alpha1, alpha2, xi, beta){   term11 <- -qgpd(alpha1, xi, beta)/xi - beta*(1-alpha1)^xi*log(1-alpha1)/xi   term12 <- qgpd(alpha1, xi, beta)/beta   term21 <- -qgpd(alpha2, xi, beta)/xi - beta*(1-alpha2)^xi*log(1-alpha2)/xi   term22 <- qgpd(alpha2, xi, beta)/beta   jacob <- matrix(c(term11, term12, term21, term22), nrow=2, byrow=true)   jacob.inv <- solve(jacob)   jacob.inv }  var.asym2 <- function(alpha1, alpha2) var.asym(alpha1, alpha2, 0.2, 1, 1000)  mop.jacob.inv2 <- function(alpha1, alpha2) mop.jacob.inv(alpha1, alpha2, 0.2, 1)  # function optimised:  object <- function(alpha1, alpha2){   term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))   sum(diag(term1)) } 

to minimise object, have additional constraint of 0 < alpha1 < alpha2 < 1. question whether can using generic optim function in r. if so, syntax, i.e., how set problem in r? , there other and/or better way? if not, there way in r? thanks.

update:

with in comment, have following:

object <- function(alpha1, alpha2){   term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))   1/sum(diag(term1))*(alpha1>0)*(alpha2>alpha1)*(alpha2<1) }  optim(c(0.01, 0.75), object) 

then got error error in stopifnot(min(p, na.rm = true) >= 0) : argument "alpha2" missing, no default. went wrong?

you can use constroptim(...) this. need change definition of object bit.

object <- function(alpha){   alpha1 <- alpha[1]   alpha2 <- alpha[2]   term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))   sum(diag(term1)) } ui <- matrix(c(1,0,-1,0,-1,1),nc=2) ci <- c(0,-1,0) result <- constroptim(th=c(0.4,0.6),object, grad=null, ui=ui, ci=ci) result$par # [1] 1.962097e-10 7.962686e-01 

the constraints applied using ui=... , ci=... arguments. ui k x p matrix p number of parameters (2 in case) , k number of constraints (3 in case), , ci vector of length k. constraints must specified that:

ui × alpha - ci ≥ 0

so in case constraints are:

α₁ ≥ 0

-α₂ + 1 ≥ 0

-α₁ + α₂ ≥ 0

the definitions of ui , ci in code above enforce constraints.

we can check result using grid search.

# check using grid search x <- seq(0,1,by=0.1) m <- expand.grid(a1=x,a2=x) m <- m[m$a1<m$a2,] grid <- apply(m,1,object) m[which.min(grid),] #    a1  a2 # 89  0 0.8 

which yields result close optimized solution.