[R] constraining functions for optim()

[Ravi Varadhan]

It makes no sense to set your constraints as the solution that you want.  I
also don’t like your non-smooth objective function: abs(y1 + y2).  It is
also wrong.  It should be abs(y1) + abs(y2), but even that is not preferable
since it is non-smooth.

What was wrong with the approach that I already proposed?  I think it makes
sense to minimize y1^2 + y2^2 for various reasons (e.g. roughly same scales,
smoothness of objective function). 

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml



----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Guillaume Théroux Rancourt
Sent: Monday, February 08, 2010 2:40 PM
To: r-help at r-project.org
Subject: [R] constraining functions for optim()

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.8
9189)/gi+164.30765+272.38*(1+2*10/165.82))^2-4*(-1/gi)*(0.89189*(164.30765+2
72.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

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.