I have a complex function that I want to maximize (I have multiplied this
function by -1 so that it becomes a minimization problem in the code below).
This function has two equality constraints.
I get the programs to run but the answer isn't correct because, when it
does converge, at least one of the constraints is violated.
Any suggestions?
Code below Violated constraint (an easy check): x[3]+x[4]+x[5]=L (global
parameter); L=12600; you may have to change the max_eval, and the
algorithms.
I included the constraint gradients (I'm pretty sure they are correct) in
case other algorithms are needed.
Thanks!
############## Load Packages (mosaic not used here, but was used to get
gradients/derivatives) #################
install.packages("nloptr")
library(nloptr)
install.packages ("mosaic")
library(mosaic)
###### Global Parameters ############
beeta=0.8
pq=10000
pf=10000
F=20
L=12600
theta=0.6
psale=0.6
mu=psale*(1-theta)
alphah=0.15
Cg=6240
Cs=2820
A= 100
D=0.0001
greekp=0.43
K=100000
##### Species Parameters##########
b1=0.38
p1=16654
v1 = 0.28
N1=6000
g1=1
delta1=1
b2=0.4
p2=2797
v2 = 0.31
N2=10000
g2=1
delta2=1
####################################### Objective function
############################################
eval_f = function (x) {
return(-1 * ( ( (1-alphah) *log(F) ) + ( alphah* (
(x[1]*(((N1/A)*(g1^greekp)*(x[3]^b1))+((1-exp(-2*D*v1*N1))*x[4])))
+ (x[2]*(((N2/A)*(g2^greekp)*(x[3]^b2))+((1-exp(-2*D*v2*N2))*x[4])
)) ) ) ) )
}
############################### Objective function gradient
########################################
eval_grad_f = function (x) {
return( c(-(alphah * (((N1/A) * (g1^greekp) * (x[3]^b1)) + ((1 - exp(-2
*
D * v1 * N1)) * x[4]))),
-(alphah * (((N2/A) * (g2^greekp) * (x[3]^b2)) + ((1 - exp(-2 *
D * v2 * N2)) * x[4]))),
-(alphah * (x[1] * ((N1/A) * (g1^greekp) * (x[3]^(b1 - 1) * b1)) +
x[2] * ((N2/A) * (g2^greekp) * (x[3]^(b2 - 1) * b2)))),
-(alphah * (x[1] * (1 - exp(-2 * D * v1 * N1)) + x[2] * (1 -
exp(-2 * D * v2 * N2)))),
0))
}
############################## 2 Equality Constraint Functions
##############################
eval_g_eq = function (x) {
return( c( ( ( ((pq*(x[5]^beeta))+((1-x[1]) * (
(p1*((N1/A)*(g1^greekp)*(x[3]^b1)) ) -
(delta1*((theta+mu)*(((N1/A)*g1^greekp*x[3]^b1)+K))) ))
+((1-x[2]) * ( (p2*((N2/A)*(g2^greekp)*(x[3]^b2))) -
(delta2*((theta+mu)*(((N2/A)*g2^greekp*x[3]^b2)+K))) ))
+((1-x[1]) * ( (p1*((1-exp(-2*D*v1*N1))*x[4])) -
(delta1*((theta+mu)*(((1-exp(-2*D*v1*N1))*x[4])+K))) ))
+((1-x[2]) * ( (p2*((1-exp(-2*D*v2*N2))*x[4]) ) -
(delta2*((theta+mu)*(((1-exp(-2*D*v2*N2))*x[4])+K))) )))
- ((pf*F)
+(x[1] * delta1 * theta * (((N1/A)*(g1^greekp)*(x[3]^b1)) +
((1-exp(-2*D*v1*N1))*x[4]) + K))
+(x[2] * delta2 * theta * (((N2/A)*(g2^greekp)*(x[3]^b2)) +
((1-exp(-2*D*v2*N2))*x[4]) + K))
+(Cg * ((N1/A)*(g1^greekp)*(x[3]^b1)) )
+(Cg * ((N2/A)*(g2^greekp)*(x[3]^b2)))
+(Cs * ((1-exp(-2*D*v1*N1))*x[4]))
+(Cs * ((1-exp(-2*D*v2*N2))*x[4]) )) ) ),
((L-x[3]-x[4]-x[5]))
) )
}
############################## 2 Equality Constraint Gradients
##############################
eval_jac_g_eq = function (x) {
return(rbind (c( -(((p1 * ((1 - exp(-2 * D * v1 * N1)) * x[4])) -
(delta1 * ((theta +
mu) * (((1 - exp(-2 * D * v1 * N1)) * x[4]) + K)))) + ((p1 *
((N1/A) * (g1^greekp) * (x[3]^b1))) - (delta1 * ((theta + mu) *
(((N1/A) * g1^greekp * x[3]^b1) + K)))) + delta1 * theta *
(((N1/A) * (g1^greekp) * (x[3]^b1)) + ((1 - exp(-2 * D * v1 *
N1)) * x[4]) + K)),
-(((p2 * ((1 - exp(-2 * D * v2 * N2)) * x[4])) - (delta2 * ((theta +
mu) * (((1 - exp(-2 * D * v2 * N2)) * x[4]) + K)))) + ((p2 *
((N2/A) * (g2^greekp) * (x[3]^b2))) - (delta2 * ((theta + mu) *
(((N2/A) * g2^greekp * x[3]^b2) + K)))) + delta2 * theta *
(((N2/A) * (g2^greekp) * (x[3]^b2)) + ((1 - exp(-2 * D * v2 *
N2)) * x[4]) + K)),
(1 - x[1]) * (p1 * ((N1/A) * (g1^greekp) * (x[3]^(b1 - 1) * b1)) -
delta1 * ((theta + mu) * ((N1/A) * g1^greekp * (x[3]^(b1 -
1) * b1)))) + (1 - x[2]) * (p2 * ((N2/A) * (g2^greekp) *
(x[3]^(b2 - 1) * b2)) - delta2 * ((theta + mu) * ((N2/A) *
g2^greekp * (x[3]^(b2 - 1) * b2)))) - (x[1] * delta1 * theta *
((N1/A) * (g1^greekp) * (x[3]^(b1 - 1) * b1)) + x[2] * delta2 *
theta * ((N2/A) * (g2^greekp) * (x[3]^(b2 - 1) * b2)) + Cg *
((N1/A) * (g1^greekp) * (x[3]^(b1 - 1) * b1)) + Cg * ((N2/A) *
(g2^greekp) * (x[3]^(b2 - 1) * b2))),
(1 - x[1]) * (p1 * (1 - exp(-2 * D * v1 * N1)) - delta1 * ((theta +
mu) * (1 - exp(-2 * D * v1 * N1)))) + (1 - x[2]) * (p2 *
(1 - exp(-2 * D * v2 * N2)) - delta2 * ((theta + mu) * (1 -
exp(-2 * D * v2 * N2)))) - (x[1] * delta1 * theta * (1 -
exp(-2 * D * v1 * N1)) + x[2] * delta2 * theta * (1 - exp(-2 *
D * v2 * N2)) + Cs * (1 - exp(-2 * D * v1 * N1)) + Cs * (1 -
exp(-2 * D * v2 * N2))),
pq * (x[5]^(beeta - 1) * beeta) ),
c(0,0, -1, -1,-1)))
}
############################## Optimization ##############################
x0 = c(0.5, 0.5, 4200, 4200, 4200) ## start values
lb=c(0, 0, 0,0,0) ## lower bound for each variable in vector x for
objective function, Inf for infinity
ub=c(1, 1, L, L, L) ## upper bound for each variable in vector x for
objective function
local_opts=list("algorithm" = "NLOPT_LN_AUGLAG",
"xtol_rel"=0.01) ## some
algorithms in nloptr require a local algorithm too
opts=list("algorithm" = "NLOPT_GN_ISRES",
"xtol_rel" = 0.01,
"maxeval"=100000 ,"local_opts" = local_opts,
"print_level"=2)
nloptr(x0=x0, eval_f = eval_f, lb=lb, ub=ub, eval_g_eq = eval_g_eq,
opts=opts)
--
Alicia Ellis
Postdoc
Gund Institute for Ecological Economics
University of Vermont
617 Main Street
Burlington, VT 05405
(802) 656-1046
http://www.uvm.edu/~aellis5/
<http://entomology.ucdavis.edu/faculty/scott/aellis/>
[[alternative HTML version deleted]]
