######## UNCONSTRAINED OPTIMIZATION ##########################
################### One dimensional example ###################

f <- function(x) -exp(-( (x-2)^2 ))
optim(1, f)$par
df <- function(x) -2*(x-2)*f(x)
optim(1, f, df, method="CG")$par
#Notice the derivative free method appears more sensitive to the starting value and gives a warning message.
#But, the converged point looks about right when the starting value is reasonable.


################### Two dimensional examples ###################
a)
f <- function(x1,y1) (1-x1)^2 + (y1 - x1^2)^2
x <- seq(-2,2,by=.1)
y <- seq(-1,3,by=.1)
z <- outer(x,y,f)
persp(x,y,z,phi=45,theta=-45,col="yellow",shade=.065,ticktype="detailed")
f <- function(x) (1-x[1])^2 + 100*(x[2]-x[1]^2)^2
optim( c(0,0), f )                         # starting values must be a vector now

b)
f <- function(x1,y1) (x1^2 + y1 - 11)^2 + (x1 + y1^2 - 7)^2      #this is called Himmelblau’s function
x <- seq(-4.5,4.5,by=.2)
y <- seq(-4.5,4.5,by=.2)
z <- outer(x,y,f)
persp(x,y,z,phi=-45,theta=45,col="yellow",shade=.65 ,ticktype="detailed")
f <- function(x) (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
optim(c(-4,-4), f)$par
[1] -3.779347 -3.283172
optim(c(2,-2), f)$par
[1]  3.584370 -1.848105
optim(c(2,2), f)$par
[1] 3.000014 2.000032
optim(c(-4,4),f)$par
[1] -2.805129  3.131435

c)
f <- function(x1,y1) 2*x1*y1 + 2*x1 - x1^2 - 2*y1^2
x <- seq(-2,2,by=.1)
y <- seq(-1,3,by=.1)
z <- outer(x,y,f)
persp(x,y,z,phi=45,theta=-45,col="yellow",shade=.065,ticktype="detailed")
f <- function(x) 2*x[1]*x[2] + 2*x[1] - x[1]^2 - 2*x[2]^2
optim( c(2,1), f,control=list(fnscale=-1) )         # starting values must be a vector now

######## OPTIMIZATION WITH GRADIENT ##########################

f1<-function(x)
{
     return(100 * (x[1] - 15)^2 + 20 * (28 - x[1])^2 + 100 * (x[2] - x[1])^2 + 20 * (38 - x[1] - x[2])^2)
}
g1=function(x){
     return(c(200 * (x[1] - 15) - 40 * (28 - x[1]) - 200 * (x[2] - x[1])-40 * (38 - x[1] - x[2]), 200*(x[2]-x[1])-40 * (38 - x[1] - x[2])))
}
optim(c(10.0,14.0),f1,g1,method="BFGS",control=list(trace=1))
optim(c(10.0,14.0),f1,g1,method="BFGS",control=list(trace=2, REPORT=1))
optim(c(10.0,14.0),f1,g1,method= "Nelder-Mead", control=list(trace=2,REPORT=1))

######## CONSTRAINED OPTIMIZATION ##########################
f2 <- function(x) {
				x1 <- x[1]
    				x2 <- x[2]
 			   -(log(x1) + x1^2/x2^2)  # need negative since constrOptim is a minimization routine
}
constrOptim(c(0.99,0.001), f2, NULL, ui=rbind(c(-1,-1),   # the -x-y > -1
                                              c(1,0),     # the x > 0
                                              c(0,1) ),   # the y > 0
                                           ci=c(-1,0, 0)) # the thresholds :-)
gr <- function(x) { ## Gradient of 'fr'
    x1 <- x[1]
    x2 <- x[2]
    c(-(1/x[1] + 2 * x[1]/x[2]^2),
       2 * x[1]^2 /x[2]^3 )
}
constrOptim(c(0.99,0.001), f2, gr, ui=rbind(c(-1,-1),  # the -x-y > -1
                                               c(1,0),    # the x > 0
                                               c(0,1) ),  # the y > 0
                                            ci=c(-1,0, 0) )
constrOptim(c(0.25,0.25), f2, NULL, #We succeeded in getting the corner values by setting the constraint slightly away from the edge
              ui=rbind( c(-1,-1), c(1,0),   c(0,1) ),  
              ci=c(-1, 0.0001, 0.0001))