Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I often use constrOptim to quickly solve nonlinear optimization problems. constrOptim works well as a general tool to tackle constrained problems like
$\min_{Ax -b \geq0} f(x)$
There are many other options and packages for specific problems but constrOptim is likely to be the first choice when little is known on the problem at hand or for exploratory/quick-and-dirty analysis.

The main problem I have with it is the nature of the constraints: they must be linear (written in matrix form as $$Ax-b≥ 0$$). Often I would like to compute (local) solutions of more general problems of the form
$\min_{g(x) -b \geq0} f(x)$
where g(x) takes vector values, one component for each scalar constraint, and is in general nonlinear, i.e., $$g(x)$$ cannot be written as $$Ax$$.

Hence, I tweaked the code of constrOptim replacing all the occurrences of the linear constraints with my “new” non linear constraints. (In order to do that, I had to discard the chance to use “BFGS” as a solving method, for now…)  The result is a function, called constrOptimNL, see the code below, which can handle non linear constraints if a vector function g(x) and b are provided. All the drawbacks of the old function are inherited, that’s life, but now non linear problems can be solved as shown by the following simple but non entirely trivial example in two dimensions.

 f <- function(x,y) x^2*y+x-x*y   fb <- function(x) f(x[1],x[2]) #vectorial version of f   x <- seq(0,3,length=31)   y <- seq(0,3,length=31)   z <- outer(x,y,f)   #draws a graph   image(x,y,z);contour(x,y,z,add=T)   polygon(c(0,0,3),c(0,3,0),col="white",density=20)  
 Heatmap and contour lines of the objective function, the feasible domain is filled in white.

Consider the constraints $$x\geq0,y\geq 1, x+y\leq 3$$, the feasible domain is a triangle with vertices (0,0), (0,3), (3,0): due to linearity, the problem can be solved with

 A <- matrix(c(1,0,-1,0,1,-1),3,2)   b <- c(0,1,-3)   res <- constrOptim(c(0.5,1.5),fb,NULL,A,b)   res
$par[1] 0.2792393 2.7207607$value[1] -0.2683538$countsfunction gradient246 NA  ... The same problems can be solved with constrOptimNL, which trivially can handle linear constraints as well:  g <- function(x,y) c(x,y,-x-y) gb <- function(x) g(x[1],x[2]) source("constrOptimNL.R") resNL <- constrOptimNL(c(0.5,1.5),fb,NULL,gb,b) resNL $par[1] 0.2792393 2.7207607$value[1] -0.2683538$countsfunction gradient246 NA
... 

Assume now that some non linear constraint is involved in the optimization problem: $$x\geq0,y\geq1,(x-1)^2+(y-1)^2\leq1$$, the feasible domain is now the upper half-circle centered at (1,1) with unit radius. This problem has nonlinearities in the constraints and cannot be solved with the standard constrOptim.