[R] constraining functions for optim()

[Guillaume Théroux Rancourt]

Hi list!

I am optimizing one function with two data sets where the unknown values (Vcmax and gi) are the same in each dataset (see script below). The script works, but I would like to add other constraints to this function.

In this function, the optimization is carried out like:

### code and functions before - see script below
y.1 = (12.7313) - f.1
y.2 = (10.52935) - f.2
y = abs(y.1 + y.2)
return(y)
}
res <- optim(par=c(0.1,130), fn=f, lower=c(0,0), upper=c(1,300), method="L-BFGS-B", control = list(maxit=10000))
# I need lower and upper boundaries, so I use the L-BFGS-B method.

I would like to also optimize where y.1 and y.2 are equal or near zero. So the three conditions here are:
y.1 = 0
y.2 = 0
y = 0

In the data set below, we get:

> res$value
[1] 0.0004200628

but: 
> y.1
[1] 4.772095
> y.2
[1] -2.379089

Thank you for your help.

guillaume


################
#### Complete script

f=function(x) { 
	
Vcmax = x[2]
gi = x[1]
		
# First data set
f.1=(-((Vcmax-0.89189)/gi+164.30765+272.38*(1+2*10/165.82))+sqrt(((Vcmax-0.89189)/gi+164.30765+272.38*(1+2*10/165.82))^2-4*(-1/gi)*(0.89189*(164.30765+272.38*(1+2*10/165.82))-Vcmax*(164.30765-(5*2/2.60546)))))/(-2/gi)
		
if (is.nan(f.1)) f.1 = 1e30
		
# Second data set
f.2=(-((Vcmax-0.89189)/gi+161.61032+272.38*(1+10*10/165.82))+sqrt(((Vcmax-0.89189)/gi+161.61032+272.38*(1+10*10/165.82))^2-4*(-1/gi)*(0.89189*(161.61032+272.38*(1+10*10/165.82))-Vcmax*(161.61032-(5*10/2.60546)))))/(-2/gi)
		
if (is.nan(f.2)) f.2 = 1e30
		
		
# Verification with measured values
y.1 = (12.7313) - f.1
y.2 = (10.52935) - f.2
		
y = abs(y.1 + y.2)
		
return(y)
}
	
res <- optim(par=c(0.1,130), fn=f, lower=c(0,0), upper=c(1,300), method="L-BFGS-B", control = list(maxit=10000))


## VERIFICATION 
# First 
f.1=function(gi, Vcmax){
(-((Vcmax-Rd)/gi+Ci.a+272.38*(1+O.a*10/165.82))+sqrt(((Vcmax-Rd)/gi+Ci.a+272.38*(1+O.a*10/165.82))^2-4*(-1/gi)*(Rd*(Ci.a+272.38*(1+O.a*10/165.82))-Vcmax*(Ci.a-(5*O.a/Sco)))))/(-2/gi)}
	
# Second 
f.2=function(gi, Vcmax){
(-((Vcmax-Rd)/gi+Ci.b+272.38*(1+O.b*10/165.82))+sqrt(((Vcmax-Rd)/gi+Ci.b+272.38*(1+O.b*10/165.82))^2-4*(-1/gi)*(Rd*(Ci.b+272.38*(1+O.b*10/165.82))-Vcmax*(Ci.b-(5*O.b/Sco)))))/(-2/gi)}
	

### y.1 and y.2 should be equal or near zero
	y.1 = 12.7313 - f.1a(res$par[1],res$par[2])
	y.2 = 10.52935 - f.2a(res$par[1],res$par[2])



Guillaume Théroux Rancourt
Ph.D. candidate --- Plant Biology
Université Laval, Québec, QC, Canada
guillaume.theroux-rancourt.1 at ulaval.ca