[R] None-linear equality constrained optimisation problems

Berend Hasselman bhh at xs4all.nl
Mon Dec 19 12:27:10 CET 2011


Michael Griffiths wrote
> 
> Dear R users,
> 
> I have a problem. I would like to solve the following:
> I have
> 
> pL = 1/(1+e^(-b0+b1))
> pM = 1/(1+e^(-b0))
> pH = 1/(1+e^(-b0-b1))
> 
> My target function is
> 
> TF= mean(pL,pM,pH) which must equal 0.5%
> 
> My non-linear constraint is
> 
> nl.Const = 1-(pM/pH), which must equal 20%, and would like the values of
> both b0 and b1 where these conditions are met.
> 
> I have searched widely for an answer, and did think that Rdonlp2 would
> suffice, only to find it no longer supported. I have solved this using
> Excel's solver function, however, because of the non-linear constraint I
> am
> having problems finding the solution in R.
> 
> Can the community suggest another method by which this might be solved?
> 

>From your description I gather that this is actually a problem of solving a
system of non-linear equations.
So you could for example use package nleqslv. As follows.

library(nleqslv)

f <- function(x) {
    b0 <- x[1]
    b1 <- x[2]
    
    pL <- 1/(1+exp(-b0+b1)) 
    pM <- 1/(1+exp(-b0)) 
    pH <- 1/(1+exp(-b0-b1)) 
    
    CR <-  (1 - pM/pH) - 0.2
    TF <-  mean(pL,pM,pH) - .005## which must equal 0.5% 
    
    c(CR,TF)
}
 
bstart <- c(.5,.5)
nleqslv(bstart ,f, control=list(trace=0))

with output (excerpt)

$x
[1] -5.0685870  0.2247176
$fvec
[1] 1.337569e-09 8.879125e-10

bstart <- c(.85,.85)
nleqslv(bstart ,f, control=list(trace=0))

with output (excerpt)

$x
[1] 1.384720 6.678025
$fvec
[1]  1.402461e-11 -3.595026e-11

So your solution for b0 and b1  depends heavily on the starting values.

BTW: your non-linear constraint is not really non-linear.
1-pM/pH=.2 is equivalent to .8 * pH = pM which is linear.


Berend


--
View this message in context: http://r.789695.n4.nabble.com/None-linear-equality-constrained-optimisation-problems-tp4213442p4213731.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list