#Copyright (C) 2013-2015 Abdalkarim Awad
#if the package is not yet installed, uncomment the following 
#line to install the nonlinear optimization package
#install.packages("nloptr")
library('nloptr') #use the library
 


#objective function
#It is the sum of all costs of all generators 
#let PA=x1, PB=x2
#f(x)=399.8+11.69*x[1]+0.00334*x[1]^2+616.9+11.83*x[2]+0.00149*x[2]^2
function_f <- function(x){ 
return(400+9*100*x[1]+0.0015*10000*x[1]^2+100+8*100*x[2]+0.0048*10000*x[2]^2 )
}
# constraint function g(x),
#x1+x2=600 ==> x1+x2-600=0
#the optimization method COBYLA supports equality constraints 
#by transforming them into two inequality constraints.
#x1+x2-600=0 ==> x1+x2-600<=0 and x1+x2-600>=0
#==> x1+x2-600<=0  and -x1-x2+600<=0 
const_g<- function( x) {
    return( rbind(
	c(x[1]+x[2]-4.5), -c(+x[1]+x[2]-4.5), 
	c(0.5-x[2]-20*x[3]+10*x[4]), -c(0.5-x[2]-20*x[3]+10*x[4]),
	c(4+10*x[3]-20*x[4]), -c(4+10*x[3]-20*x[4])
	
	) )
}

 
# Solve using NLOPT_LN_COBYLA 
results <- nloptr( x0=c(1,0.5,0,0), #initial values
	eval_f=function_f,  #min (f(x))	
	lb = c(1,0.5,-Inf,-0.17), #lower bound
	ub = c(6,4,Inf,0.17), #upper bound
	eval_g_ineq= const_g, #equality and inequality constraints
	opts = list("algorithm"="NLOPT_LN_COBYLA",xtol_rel = 1e-4)) #options
print( results )
x=results$solution
PA=x[1]
PB=x[2]
theta1=0
theta2=x[3]
theta3=x[4]
PA
PB 
theta1
theta2
theta3