R/essR.R
In jaegea/MEIGOR: MEIGOR - MEtaheuristics for bIoinformatics Global Optimization
Defines functions
essR
Documented in
essR
essR <-
function(problem,opts=list(maxeval=NULL,maxtime=NULL),...){
# Function : essR
# R code of the eSS optimization code from: Process Engineering Group IIM-CSIC
# Created on : 09/06/2011
# Last Update: 25/07/2012
# Email : josea.egea@gmail.com
#
# (c) CSIC, Spanish Council for Scientific Research, Spain
#
# Global optimization algorithm for MINLPs based on Scatter Search
#
# essR attempts to solve problems of the form:
# min F(x) subject to: ceq(x) = 0 (equality constraints)
# x c_L <= c(x) <= c_U (inequality constraints)
# x_L <= x <= x_U (bounds on the decision variables)
#
# Constraint functions, if applicable, must be declared in the same script as the
# objective function as a second output argument:
# e.g., myfunction <- function(x){
# calculate fx - scalar containing the objective function value
# calculate gx - vector (or empty) containing the constraints values
# return(list(fx,gx))
# }
#
# Please have a look at the manual before using essR
#
# USAGE: Results <- essR(problem,opts,p1,p2,....,pn);
#
# INPUT PARAMETERS:
# ****************
# problem - List containing problem settings
# problem$f = Name of the file containing the objective
# function (String)
# problem$x_L = Lower bounds of decision variables (vector)
# problem$x_U = Upper bounds of decision variables (vector)
# problem$x_0 = Initial point(s) (optional; vector or matrix)
# problem.$f_0 = Function values of initial point(s) (optional) These values MUST correspond to feasible points.
#
# NOTE:The dimension of f_0 and x_0 may be different. For example, if
# we want to introduce 5 initial points but we only know the values
# for 3 of them, x_0 would have 5 rows whereas f_0 would have only 3
# elements. In this example, it is mandatory that the first 3 rows of x_0 correspond
# to the values of f_0
#
# Additionally, fill the following fields if your problem has
# non-linear constraints
# problem$neq = Number of equality constraints (Integer; do not define it
# if there are no equality constraints)
# problem$c_L = Lower bounds of nonlinear inequality constraints
# (vector)
# problem$c_U = Upper bounds of nonlinear inequality constraints
# (vector)
# problem$int_var = Number of integer variables (Integer)
# problem$bin_var = Number of binary variables (Integer)
# problem$vtr = Objective function value to be reached (optional)
#
#
# NOTE: The order of decision variables is x <- c(cont, int, bin)
#
# opts - List containing options (if set as opts <- numeric(0) default options
# will be loaded). See the script of default options to know their values
#
# User options
# opts$maxeval = Maximum number of function evaluations
# (Default 1000)
# opts$maxtime = Maximum CPU time in seconds (Default 60)
# opts$iterprint = Print each iteration on screen: 0-Deactivated
# 1-Activated (Default 1)
# opts$weight = Weight that multiplies the penalty term added
# to the objective function in constrained
# problems (Default 1e6)
# opts$log_var = Indexes of the variables which will be used
# to generate diverse solutions in different
# orders of magnitude (vector)
# opts$tolc = Maximum absolute violation of the constraints
# (Default 1e-5)
# opts$prob_bound = Probability (0-1) of biasing the search towards
# the bounds (Default 0.5)
# opts$save_results = Saves final results in eSSR_report.RData. (Binary; Default = 1)
# opts$inter_save = Saves results of intermediate iterations in eSSR_report.RData
# Useful for very long runs. (Binary; Default = 0)
#
# Global options
# opts$dim_refset = Number of elements in Refset
# (Integer; automatically calculated)
# opts$ndiverse = Number of solutions generated by the
# diversificator (Default 10*nvar)
# opts$combination = Type of combination of Refset
# elements (Default 1)
# 1: hyper-rectangles
# 2: linear combinations
#
# Local options
# opts$local_solver = Choose local solver
# 0: Local search deactivated (Default),
# "NM", "BFGS", "CG", "LBFGSB",
# "SA","SOLNP", "DHC", "NLS2SOL"
# opts$local_tol = Level of tolerance in local
# search
# opts$local_iterprint = Print each iteration of local
# solver on screen (Binary;
# default = 0)
# opts$local_n1 = Number of iterations
# before applying
# local search for the 1st time
# (Default 1)
# opts$local_n2 = Minimum number of iterations
# in the global
# phase between 2 local calls
# (Default 10)
# opts$local_balance = Balances between quality (=0)
# and diversity (=1) for
# choosing initial points for
# the local search (default
# 0.5)
# opts$local_finish = Applies local search to the
# best solution found once the
# optimization if finished
# (same values as
# opts.local.solver)
# opts$local_bestx = When activated (i.e. =1) only
# applies local search to the
# best solution found to
# date,ignoring filters
# (Default=0)
#
#
# p1,p2... : optional input parameters to be passed to the objective
# function
#
#
# OUTPUT PARAMETERS:
# *****************
# A file called "eSSR_report.Rdata" is generated containing.
#
# Results - Data frame containing results
# problem - Data frame containing problem settings
# opts - Data frame containing all options
#
#
# Fields in Results
# Results$fbest = Best objective function value
# found after the optimization
# Results$xbest = Vector providing the best
# function value
# Results$cpu_time = Time in seconds consumed in the
# optimization
# Results$f = Vector containing the best
# objective function value after each
# iteration
# Results$x = Matrix containing the best vector
# after each iteration
# Results$time = Vector containing the cpu time
# consumed after each iteration
# Results$neval = Vector containing the number of
# function evaluations after each
# iteration
# Results$numeval = Number of function evaluations
# Results$local_solutions = Local solutions found by the
# local solver (in rows)
# Results$local_solutions_values = Function values of the local
# solutions
# Results$end_crit = Criterion to finish the
# optimization
# 1: Maximal number of function
# evaluations achieved
# 2: Maximum allowed CPU Time
# achieved
# 3: Value to reach achieved
#
#
# REFERENCES
#
# If you use essR and publish the results, please cite the following papers:
#
# Egea, J.A., Marti, R., Banga, J.R. (2010) An evolutionary method for
# complex-process optimization. Computers & Operations Research 37(2):
# 315-324.
#
# Egea, J.A., Balsa-Canto, E., Garcia, M.S.G., Banga, J.R. (2009) Dynamic
# optimization of nonlinear processes with an enhanced scatter search
# method. Industrial & Engineering Chemistry Research 49(9): 4388-4401.
# Function : eSS Release 2014A
print('eSS R2014A - Enhanced Scatter Search');
stopOptimization<-0;
# Initialize time
cpu_time=proc.time()[3];
#Take the 3rd component of cpu_time for the total elapsed time
#Initialize the number of function evaluations and a stop variable
nfuneval <- 0;
fin <- 0;
if (!is.numeric(opts$maxeval) && !is.numeric(opts$maxtime)){
cat("WARNING:Either opts.maxeval or opts.maxtime must be defined as a stop criterion \n")
cat("Define at least one of these two options and rerun \n")
Results <- numeric(0);
return(Results)
stop()
} else{
if (!is.numeric(opts$maxeval)){
maxeval <- 1e12;
} else{
maxeval <- opts$maxeval;
}
if (!is.numeric(opts$maxtime)){
maxtime <- 1e12;
} else{
maxtime <- opts$maxtime;
}
if(!(maxtime>0)){
Results <- numeric(0);
return(Results);
}
if(!(maxeval>0)){
Results <- numeric(0);
return(Results);
}
}
cat("\n ------------------------------------------------------------------------------ \n")
cat(" essR - Enhanced Scatter Search in R \n");
cat("<c> IIM-CSIC, Vigo, Spain - email: gingproc@iim.csic.es \n");
cat("------------------------------------------------------------------------------ \n\n");
#Load default values for the options
default <- ssm_defaults();
nargin <- length(as.list(match.call())) -1;
if (nargin<2){opts <- numeric(0)};
#Set all optiions
opts<-ssm_optset(default,opts);
x_U<-problem$x_U;
x_L<-problem$x_L;
#Check if the bounds have the same dimension
if (length(x_U)!= length(x_L)){
cat("Upper and lower bounds have different dimensions!!! \n")
cat("EXITING")
Results<-numeric(0);
return(Results)
stop()
} else{
#Number of decision variables
nvar<-length(x_L);
}
prob_names<-names(problem);
temp<-match("x_0",prob_names);
if (is.na(temp)){x_0<-numeric(0)} else{x_0<-problem$x_0}
temp<-match("f_0",prob_names);
if (is.na(temp)){
f_0<-numeric(0);
fbest<-Inf;
xbest<-rep(Inf,nvar)
} else{
f_0<-problem$f_0;
fbest<-min(f_0);
iii<-which.min(f_0);
xbest<-x_0[iii,];
}
temp<-match("neq",prob_names);
if (is.na(temp)){neq<-0} else{neq<-problem$neq}
temp<-match("c_U",prob_names);
if (is.na(temp)){c_L<-numeric(0); c_U<-numeric(0)} else{c_U<-problem$c_U; c_L<-problem$c_L}
temp<-match("int_var",prob_names);
if (is.na(temp)){int_var<-0} else{int_var<-problem$int_var}
temp<-match("bin_var",prob_names);
if (is.na(temp)){bin_var<-0} else{bin_var<-problem$bin_var}
temp<-match("vtr",prob_names);
if (is.na(temp)){vtr<-numeric(0)} else{vtr<-problem$vtr}
#User options
iterprint=opts$iterprint;
prob_bound=opts$prob_bound;
save_results=opts$save_results;
inter_save=opts$inter_save;
plot_results=opts$plot;
weight=opts$weight;
tolc=opts$tolc;
log_var=opts$log_var;
#Global options
dim_refset=opts$dim_refset;
if (is.character(dim_refset)){
z<-c(-10*nvar/1,-1,1);
nnn<-polyroot(z);
iii<-which(Re(nnn)>0);
dim_refset<-ceiling(Re(nnn[iii]));
#%% is the modulus of a division e.g. 7%%4 provides 3 (like mod(7,4) in Matlab)
if (dim_refset%%2){
dim_refset<-dim_refset+1
if (iterprint){
cat("Refset size automatically calculated:",dim_refset, "\n");
}
}
}
ndiverse<-opts$ndiverse;
if (is.character(ndiverse)){
ndiverse=10*nvar;
if (iterprint){
cat("Number of diverse solutions automatically calculated:",ndiverse,"\n");
}
}
if (ndiverse<dim_refset){ndiverse<-dim_refset}
combin=opts$combination;
#Local options
local_solver=opts$local_solver;
if(is.character(local_solver)){local_solver<-toupper(local_solver)}
local_tol=opts$local_tol;
local_iterprint=opts$local_iterprint;
local_n1<-opts$local_n1;
local_n2<-opts$local_n2;
if (is.character(local_n1)){local_n1<-1}
if (is.character(local_n2)){local_n2<-10}
local_balance=opts$local_balance;
if (!length(opts$local_finish)){
local_finish<-local_solver;
} else{
local_finish=opts$local_finish;
}
local_bestx=opts$local_bestx;
#If there are equality constraints
if (neq){
#Set the new bounds for all the constraints
c_L<-c(rep(0,neq),c_L);
c_U<-c(rep(0,neq),c_U)
}
nconst<-length(c_U);
fobj<-problem$f;
#Check if the objective function has 2 output arguments in constrained problems
if (nconst){
n_out_f<-do.call(fobj,list(x_L,...));
if (length(n_out_f)<2){
cat("For constrained problems the objective function must have at least 2 output arguments \n");
cat("EXITING \n");
Results<-numeric(0);
return(Results)
stop()
}
}
#Change to log
xl_log<-x_L;
xu_log<-x_U;
if (length(log_var)){
aaa<-which(!xl_log[log_var]);
if (length(aaa)){
xl_log[log_var[aaa]]<-1e-10;
}
xl_log[log_var]<-log(xl_log[log_var]);
xu_log[log_var]<-log(xu_log[log_var]);
}
#Check if the objective function has 3 output arguments when NL2SOL is invoked
if(match(local_solver,"NL2SOL", nomatch=0) | match(local_finish,"NL2SOL", nomatch=0)){
n_out_f<-do.call(fobj,list(x_L,...));
nfuneval<-nfuneval+1;
if (length(n_out_f)<3){
cat("NL2SOL requires 3 output arguments: f, g, and the vector of residuals \n");
cat("EXITING \n");
Results<-numeric(0);
return(Results)
stop()
}
}
# TRANSFORM THIS INTO R CODE WHEN THE LOCAL SEARCHES ARE IMPLEMENTED!!!!!!
# %Build auxiliary files in case local search is activated
# if local_solver
# ssm_aux_local(problem.f,x_L,x_U,c_L,c_U,neq,local_solver,nvar,varargin{:});
# end
#Initialize variables
stage_1<-1;
n_minimo<-0;
n_critico<-local_n1;
local_solutions<-numeric(0);
local_solutions_values<-numeric(0);
initial_points<-numeric(0);
Results<-list(f=numeric(0), x=numeric(0), time=numeric(0), neval=numeric(0));
#Possible combinations among the dim_refset elements in Refset, taken in pairs
ncomb<-t(combn(1:dim_refset,2));
MaxSubSet<-(dim_refset^2-dim_refset)/2;
MaxSubSet2<-2*MaxSubSet;
#We generate 5 solutions
solutions<-matrix(0,ndiverse+5,nvar);
solutions[1:5,]<-(matrix(runif(nvar*5),5,nvar)+1:5-1)/5;
solutions[6:(ndiverse+5),]<-matrix(runif(ndiverse*nvar),ndiverse,nvar);
solutions<-solutions*kronecker(matrix(1,ndiverse+5,1),t((xu_log-xl_log)))+kronecker(matrix(1,ndiverse+5,1),t(xl_log));
for (i in 1:(ndiverse+5)){
solutions[i,log_var]=exp(solutions[i,log_var]);
}
l_f_0<-length(f_0);
if (is.vector(x_0) && length(x_0)){l_x_0<-1} else if(is.matrix(x_0)){l_x_0<-nrow(x_0)} else{l_x_0<-0};
if (!length(l_x_0)){l_x_0<-0}
#Put x_0 WITHOUT their corresponding f_0 values at the beginning of solutions
solutions<-rbind(x_0, solutions)
if (l_f_0){solutions<-solutions[-seq(1,l_f_0),]}
temp<-l_x_0-l_f_0;
sol<-matrix(0,ndiverse+5+temp,nvar);
sol_val<-rep(0,ndiverse+5+temp);
sol_val_pen<-rep(0,ndiverse+5+temp);
sol_pen<-rep(0,ndiverse+5+temp);
sol_nlc<-matrix(0,ndiverse+5+temp,nconst);
counter<-1;
for (i in 1:(ndiverse+5+temp)){
res<-ssm_evalfc(solutions[i,],x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
sol[counter,]<-x;
sol_val[counter]<-val;
sol_val_pen[counter]<-val_penalty;
sol_pen[counter]<-pena;
sol_nlc[counter,]<-nlc;
counter<-counter+1;
}
}
nsol<-nrow(sol);
if ((counter-1)<nsol){
sol<-sol[-seq(counter-1,nsol),];
sol_val<-sol_val[-seq(counter-1,nsol)];
sol_val_pen<-sol_val_pen[-seq(counter-1,nsol)];
sol_pen<-sol_pen[-seq(counter-1,nsol)];
sol_nlc<-sol_nlc[-seq(counter-1,nsol),];
}
#Add points with f_0 values. We assume feasiblity
if (l_f_0){
sol<-rbind(x_0[1:l_f_0,], sol);
sol_val<-c(f_0, sol_val);
#We assume feasible points in f_0
sol_val_pen<-c(f_0, sol_val_pen);
sol_pen<-c(rep(0,l_f_0),sol_pen);
#We do not know these values, but we do not care. We assume feasibility
sol_nlc<-rbind(NaN*matrix(0,l_f_0,nconst),sol_nlc);
}
#Initialize Refset
Refset<-matrix(0,dim_refset,nvar);
Refset_values<-rep(0,dim_refset);
Refset_values_penalty<-rep(0,dim_refset);
Refset_nlc<-numeric(0);
penalty<-numeric(0);
temp<-sort(sol_val_pen,index.return=TRUE);
iii<-temp[[2]];
first_members<-ceiling(dim_refset/2);
last_members<-dim_refset-first_members;
#Create the initial sorted Refset
Refset<-sol[iii[1:first_members],];
Refset_values<-sol_val[iii[1:first_members]];
Refset_values_penalty<-sol_val_pen[iii[1:first_members]];
penalty<-sol_pen[iii[1:first_members]];
Refset_nlc<-sol_nlc[iii[1:first_members],];
#The rest of Refset members are chosen randomly
iii<-iii[-(1:first_members)];
temp<-1:length(iii);
ooo<-sample(temp);
#The comment below comes from matlab code
#ooo=randperm(ndiverse+2+size(x_0,1)-first_members);
Refset<-rbind(Refset, sol[iii[ooo[1:last_members]],]);
Refset_values<-c(Refset_values, sol_val[iii[ooo[1:last_members]]]);
Refset_values_penalty<-c(Refset_values_penalty, sol_val_pen[iii[ooo[1:last_members]]]);
penalty<-c(penalty, sol_pen[iii[ooo[1:last_members]]]);
Refset_nlc<-rbind(Refset_nlc, sol_nlc[iii[ooo[1:last_members]],]);
#Check best value
ggg<-which(penalty<=tolc);
if (length(ggg)){
fbest<-min(Refset_values_penalty[ggg]);
iii<-which.min(Refset_values_penalty[ggg]);
xbest<-Refset[ggg[iii],];
}
iter<-0;
if (iterprint){
cat("Initial Pop: NFunEvals:", nfuneval, "Bestf:",fbest,"CPUTime:",proc.time()[3]-cpu_time,"Mean:", mean(Refset_values_penalty),"\n");
}
Results$f=c(Results$f,fbest);
Results$x=rbind(Results$x,xbest);
Results$time=c(Results$time, proc.time()[3]-cpu_time);
Results$neval=c(Results$neval, nfuneval);
#Sort the Refset
temp<-sort(Refset_values_penalty,index.return=TRUE);
Refset_values_penalty<-temp[[1]];
I<-temp[[2]];
Refset_values<-Refset_values[I];
penalty<-penalty[I];
Refset<-Refset[I,];
Refset_nlc<-Refset_nlc[I,];
if (combin==1){nrand=nvar} else if (combin==2) {nrand=1}
use_bestx<-0;
index1<-ncomb[,1];
index2=ncomb[,2];
index<-c(index1,index2);
diff_index<-index2-index1;
#Minimum ppp
lb_p<-1;
#Defines maximum ppp. max(ppp)= lb_p+ st_p
st_p<-0.75;
#Example: lb_p=1.25 and st_p=0.5 varies ppp from 1.25 to 1.75
ppp<-st_p*(diff_index-1)/(dim_refset-2)+lb_p;
hyper_x_L<-kronecker(matrix(1,MaxSubSet,1),t(x_L));
hyper_x_U<-kronecker(matrix(1,MaxSubSet,1),t(x_U));
refset_change<-rep(0,dim_refset);
while (!fin){
child<-matrix(0,MaxSubSet2,nvar);
child_values<-rep(0,MaxSubSet2);
child_values_penalty<-rep(0,MaxSubSet2);
child_penalty<-rep(0,MaxSubSet2);
child_nlc<-matrix(0,MaxSubSet2,nconst);
child_parent_index<-rep(0,MaxSubSet2);
members_update<-rep(0,dim_refset);
counter<-1;
continuar<-0;
while (!continuar){
#Sort Refset
#This is mandatory because the variable index is also sorted
counter<-1;
temp<-sort(Refset_values_penalty,index.return=TRUE);
Refset_values_penalty<-temp[[1]];
I<-temp[[2]];
Refset_values<-Refset_values[I];
penalty<-penalty[I];
Refset<-Refset[I,];
Refset_nlc<-Refset_nlc[I,];
refset_change<-refset_change[I];
parents_index1<-Refset[index1,];
parents_index2<-Refset[index2,];
parents_index1_values_penalty<-Refset_values_penalty[index1];
parents_index2_values_penalty<-Refset_values_penalty[index2];
denom22<-pmax(abs(parents_index1),abs(parents_index2));
denom22[which(!denom22)]<-1;
AAA<-abs((parents_index1-parents_index2)/denom22);
#Compute the maximum value of AAA by rows
temp<-apply(AAA,1,max);
#if any of these maximum values are lower than a threshold, then keep the indexes
BBB<-which(temp<1e-3);
if (length(BBB)){
#Check the index of worst element among the twins of index2
index_refset_out<-max(index2[BBB]);
includ=0;
while (!includ){
#Generate a new random solution
new_refset_member<-runif(nvar)*(xu_log-xl_log)+xl_log;
new_refset_member[log_var]=exp(new_refset_member[log_var]);
res<-ssm_evalfc(new_refset_member,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
Refset[index_refset_out,]<-x;
Refset_values[index_refset_out]<-val;
Refset_values_penalty[index_refset_out]<-val_penalty;
Refset_nlc[index_refset_out,]<-nlc;
penalty[index_refset_out]<-pena;
members_update[index_refset_out]<-0;
}
}
} else{
continuar=1;
}
candidate<-Refset;
candidate_values<-Refset_values;
candidate_values_penalty<-Refset_values_penalty;
candidate_nlc<-Refset_nlc;
candidate_penalty<-penalty;
candidate_update<-rep(0,dim_refset);
members_update<-rep(1,dim_refset);
factor<-kronecker(matrix(1,1,nvar),ppp)*(parents_index2-parents_index1)/1.5;
v1<-parents_index1-factor;
v2<-parents_index2-factor;
v3<-2*parents_index2-parents_index1-factor;
aaa<-which(v1<hyper_x_L);
if (length(aaa)){
rand_aaa<-runif(length(aaa));
AAA<-which(rand_aaa>prob_bound);
aaa<-aaa[AAA];
v1[aaa]=hyper_x_L[aaa];
}
bbb<-which(v1>hyper_x_U);
if (length(bbb)){
rand_bbb<-runif(length(bbb));
BBB<-which(rand_bbb>prob_bound);
bbb<-bbb[BBB];
v1[bbb]=hyper_x_U[bbb];
}
ccc<-which(v3<hyper_x_L);
if (length(ccc)){
rand_ccc<-runif(length(ccc));
CCC<-which(rand_ccc>prob_bound);
ccc<-ccc[CCC];
v3[ccc]=hyper_x_L[ccc];
}
ddd<-which(v3>hyper_x_U);
if (length(ddd)){
rand_ddd<-runif(length(ddd));
DDD<-which(rand_ddd>prob_bound);
ddd<-ddd[DDD];
v3[ddd]=hyper_x_U[ddd];
}
if (nrand>1){
new_comb1<-v1+(v2-v1)*matrix(runif(MaxSubSet*nrand),MaxSubSet,nrand);
new_comb2=v2+(v3-v2)*matrix(runif(MaxSubSet*nrand),MaxSubSet,nrand);
} else{
new_comb1=v1+(v2-v1)*runif(1);
new_comb2=v2+(v3-v2)*runif(1);
}
new_comb<-rbind(new_comb1,new_comb2);
for (i in 1:MaxSubSet2){
#Evaluate
res<-ssm_evalfc(new_comb[i,],x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
child[counter,]<-x;
child_values[counter]<-val;
child_values_penalty[counter]<-val_penalty;
child_penalty[counter]<-pena;
child_nlc[counter,]<-nlc;
child_parent_index[counter]<-index[i];
counter<-counter+1;
if (val_penalty<candidate_values_penalty[index[i]]){
#Update candidate
#This can be accelerated saving the index and replacing at the end with the best one
candidate[index[i],]<-x;
candidate_values_penalty[index[i]]<-val_penalty;
candidate_values[index[i]]<-val;
candidate_penalty[index[i]]<-pena;
candidate_nlc[index[i],]<-nlc;
candidate_update[index[i]]<-1;
members_update[index[i]]<-0;
}
}
}
#Check the stop criterion
if (nfuneval>= maxeval){
fin<-1;
} else if(proc.time()[3]-cpu_time>=maxtime){
fin<-2;
} else if(length(vtr)){
if (fbest<=vtr){
fin<-3;
}
}
child[-(1:(counter-1)),]<-numeric(0);
child_values[-(1:(counter-1))]<-numeric(0);
child_values_penalty[-(1:(counter-1))]<-numeric(0);
child_penalty[-(1:(counter-1))]<-numeric(0);
child_parent_index[-(1:(counter-1))]<-numeric(0);
members_to_update<-which(candidate_update==1);
if (length(members_to_update)){
for (i in 1:length(members_to_update)){
res <- ssm_beyond(Refset[members_to_update[i],],candidate[members_to_update[i],],candidate_values_penalty[members_to_update[i]],
fobj,nrand,tolc,weight,x_L,x_U,c_L,c_U,nconst,int_var,bin_var,nfuneval,prob_bound,nvar,...)
vector<-res[[1]];
vector_value<-res[[2]];
vector_penalty<-res[[3]];
vector_value_penalty<-res[[4]];
vector_nlc<-res[[5]];
new_child<-res[[6]];
new_child_value<-res[[7]];
new_child_penalty<-res[[8]];
new_child_value_penalty<-res[[9]];
new_child_nlc<-res[[10]];
nfuneval<-res[[11]];
if (length(vector)){
Refset[members_to_update[i],]<-vector;
Refset_values[members_to_update[i]]<-vector_value;
Refset_values_penalty[members_to_update[i]]<-vector_value_penalty;
penalty[members_to_update[i]]<-vector_penalty;
Refset_nlc[members_to_update[i],]<-vector_nlc;
} else{
Refset[members_to_update[i],]<-candidate[members_to_update[i],];
Refset_values[members_to_update[i]]<-candidate_values[members_to_update[i]];
Refset_values_penalty[members_to_update[i]]<-candidate_values_penalty[members_to_update[i]];
penalty[members_to_update[i]]<-candidate_penalty[members_to_update[i]];
Refset_nlc[members_to_update[i],]<-candidate_nlc[members_to_update[i],];
}
}
}
fbest_lastiter=fbest;
ggg<-which(penalty<=tolc);
if (length(ggg)){
fbest_new<-min(Refset_values_penalty[ggg]);
hhh<-which.min(Refset_values_penalty[ggg]);
if (fbest_new<fbest){
xbest<-Refset[ggg[hhh],];
fbest<-fbest_new;
use_bestx<-1;
if (inter_save){
Results$fbest<-fbest;
Results$xbest<-xbest;
Results$numeval<-nfuneval;
Results$end_crit<-fin;
Results$cpu_time<-cpu_time;
Results$Refset$x<-Refset;
Results$Refset$f<-Refset_values;
Results$Refset$fpen<-Refset_values_penalty;
Results$Refset$const<-Refset_nlc;
Results$Refset$penalty<-penalty;
if (save_results){save(opts, problem, Results, file = "eSSR_report.RData")}
}
}
}
iter<-iter+1;
n_minimo<-n_minimo+1;
if (iterprint){
cat("Iteration:",iter, "NFunEvals:", nfuneval, "Bestf:",fbest,"CPUTime:",proc.time()[3]-cpu_time,"Mean:", mean(Refset_values_penalty),"\n");
}
ddd<-which(members_update>0);
eee<-which(!members_update);
refset_change[ddd]<-refset_change[ddd]+1;
refset_change[eee]<-0;
fff<-which(refset_change>=20);
if (length(fff)){
to_replace<-length(fff);
replaced<-0;
while (replaced<to_replace){
newsol<-runif(nvar)*(xu_log-xl_log)+xl_log;
newsol[log_var]=exp(newsol[log_var]);
res<-ssm_evalfc(newsol,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
Refset[fff[replaced+1],]<-x;
Refset_values[fff[replaced+1]]<-val;
Refset_values_penalty[fff[replaced+1]]<-val_penalty;
penalty[fff[replaced+1]]<-pena;
Refset_nlc[fff[replaced+1],]<-nlc;
refset_change[fff[replaced+1]]<-0;
replaced<-replaced+1;
}
}
}
if (fbest<fbest_lastiter){
Results$f<-c(Results$f,fbest);
Results$x<-rbind(Results$x,xbest);
Results$time<-c(Results$time, proc.time()[3]-cpu_time);
Results$neval<-c(Results$neval,nfuneval);
#if plot_results==1
# stairs(Results.time,Results.f)
# xlabel('CPU Time (s)');
# ylabel('Objective Function Value (-)');
# title(strcat('Best f: ',num2str(fbest)))
# drawnow
# end
}
#******************************************************************************************
#************************* LOCAL FILTERS ROUTINE ******************************************
call_filters<-0;
if(is.character(local_solver) && n_minimo>=n_critico){
call_filters<-1;
if (local_bestx && !use_bestx){
call_filters<-0;
}
}
if (call_filters){
if (stage_1==1){
#In this first stage we apply the local search over the best solution found so far
minimo<-min(Refset_values_penalty);
I<-which.min(Refset_values_penalty);
x0<-Refset[I,];
f0<-minimo;
#Switch to stage 2
n_critico<-local_n2;
if (local_bestx){
#not to go directly to stage 2 the next time we pass over here unless we have really improved fbest
use_bestx<-0;
}
#use_bestx is a flag indicating that the best solution has been improved
}
else if (stage_2==1){
if (local_bestx && use_bestx){
x0<-xbest;
f0<-fbest;
use_bestx<-0;
} else{
#Calculate distance between our children and local_solutions+initial points
local_init<-rbind(local_solutions, initial_points);
xuxl<-x_U-x_L;
#child<-rbind(child,new_child);
# child_values<-c(child_values,new_child_value);
# child_values_penalty<-c(child_values_penalty,new_child_value_penalty);
# child_penalty<-c(child_penalty,new_child_penalty);
# child_nlc<-rbind(child_nlc,new_child_nlc);
child_norm<-child/kronecker(matrix(1,nrow(child),1),t(xuxl));
local_init_norm<-local_init/kronecker(matrix(1,nrow(local_init),1),t(xuxl));
distance<-eucl_dist(child_norm,local_init_norm);
#Now we need to sort the children according the their distances to local minima and intial points
#Here we compute the minimum euclidean distance between child i and all the vector in initial_points+local_solutions
temp<-apply(distance,1,min);
#Now we sort those minimum distances and favour the highest ones to favour diverse points
res<-sort(temp,decreasing=TRUE,index.return=TRUE);
ooo<-res[[2]];
res<-sort(child_values_penalty,index.return=TRUE);
uuu<-res[[2]];
n_child<-length(child_values_penalty)
child_points<-rep(0,n_child);
www<-local_balance;
for (i in 1:n_child){
child_points[ooo[i]]<-child_points[ooo[i]]+www*i;
child_points[uuu[i]]<-child_points[uuu[i]]+(1-www)*i;
}
I<-which.min(child_points);
x0<-child[I,];
f0<-child_values_penalty[I];
}
}
if (iterprint){
cat("Call local solver:", toupper(local_solver),"\n")
if (child_penalty[I]>tolc){
cat("Initial point function value:",child_values[I], "(INFEASIBLE)","\n");
} else{
cat("Initial point function value:",f0,"\n");
}
tic<-proc.time()[3];
}
#AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
res<-ssm_localsolver(x0,x_L,x_U,c_L,c_U,neq,int_var,bin_var,fobj,local_solver,local_iterprint,local_tol,weight,nconst,tolc,...);
x<-res[[1]];
fval<-res[[2]];
numeval<-res[[3]];
initial_points<-rbind(initial_points, x0);
nfuneval<-nfuneval+numeval;
#Evaluate the local solution
res<-ssm_evalfc(x,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ && pena<=tolc){
if (iterprint){
cat("Local solution function value:",val,"\n");
}
if (stage_1==1){
stage_1=0;
stage_2=1;
}
if (val_penalty<fbest){
fbest<-val_penalty;
xbest<-x;
if (fbest<fbest_lastiter){
Results$f<-c(Results$f,fbest);
Results$x<-rbind(Results$x,xbest);
Results$time<-c(Results$time,proc.time()[3]-cpu_time);
Results$neval<-c(Results$neval,nfuneval);
#if plot_results==1
# stairs(Results.time,Results.f)
# xlabel('CPU Time (s)');
# ylabel('Objective Function Value (-)');
# title(strcat('Best f: ',num2str(fbest)))
# drawnow
# end
}
}
} else{
cat("Local solution function value:","UNFEASIBLE","\n");
}
cat("Number of function evaluations in the local search:",numeval,"\n");
cat("CPU Time of the local search:",proc.time()[3]-tic,"seconds \n\n");
#Check the stopping criterion
if (nfuneval>= maxeval){
fin<-1;
} else if(proc.time()[3]-cpu_time>=maxtime){
fin<-2;
} else if(length(vtr)){
if (fbest<=vtr){
fin<-3;
}
}
if (includ){
if (length(local_solutions)){
#Add the found local solution to the list (if it is a new one)
adicionar_local<-1;
res<-ssm_isdif2(x,local_solutions,1e-2,1);
f<-res[[1]];
ind<-res[[2]];
#If it is too close to a previously found local solution, we do not add it
if (f){
adicionar_local<-0
}
} else {
#Add it if no local solutions have been found yet
adicionar_local<-1;
}
if (adicionar_local){
#choose either this (replace the worst element in Refset)
# [aaa,jjj]=max(Refset_values_penalty);
#
# if val_penalty<Refset_values_penalty(jjj)
# Refset(jjj,:)=x;
# Refset_values(jjj)=val;
# Refset_values_penalty(jjj)=val_penalty;
# Refset_nlc(jjj,:)=nlc;
# % penalty(jjj)=pena;
# % end
#
# %or this (Replace the initial point)
#
# % if val_penalty<f0
# % child(I,:)=x;
# % child_values(I)=val;
# % child_values_penalty(I)=val_penalty;
# % child_nlc(I,:)=nlc;
# % end
#
#or this (Replace the parent of the initial point)
if (val_penalty<Refset_values_penalty[child_parent_index[I]]){
Refset[child_parent_index[I],]<-x;
Refset_values[child_parent_index[I]]<-val;
Refset_values_penalty[child_parent_index[I]]<-val_penalty;
Refset_nlc[child_parent_index[I],]<-nlc;
penalty[child_parent_index[I]]<-pena;
refset_change[child_parent_index[I]]<-0;
}
local_solutions<-rbind(local_solutions,x);
local_solutions_values<-c(local_solutions_values,fval);
}
}
n_minimo<-0
}
#************************* END OF LOCAL FILTERS ROUTINE ***********************************
#******************************************************************************************
fbest_lastiter<-fbest;
#Check the stopping criterion
if (nfuneval>= maxeval){
fin<-1;
} else if(proc.time()[3]-cpu_time>=maxtime){
fin<-2;
} else if(length(vtr)){
if (fbest<=vtr){
fin<-3;
}
}
if (fin){
cat("\n");
if (is.character(local_finish) && fin<3){
local_tol=3;
if (is.numeric(xbest)){
x0=xbest;
f0=fbest;
}else{
x0=Refset[1,];
f0=Refset_values_penalty[1];
}
if (iterprint){
cat("Final local refinement: ", toupper(local_finish),"\n")
cat("Initial point function value: ", f0,"\n");
tic<-proc.time()[3];
}
res<-ssm_localsolver(x0,x_L,x_U,c_L,c_U,neq,int_var,bin_var,fobj,local_finish,local_iterprint,local_tol,weight,nconst,tolc,...);
x<-res[[1]];
fval<-res[[2]];
numeval<-res[[3]];
#Evaluate the local solution
res<-ssm_evalfc(x,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ && pena<=tolc){
if (iterprint){
cat("Local solution function value:",val,"\n");
}
if (val_penalty<fbest){
fbest<-val_penalty;
xbest<-x;
}
} else{
cat("Local solution function value:","UNFEASIBLE","\n");
}
cat("Number of function evaluations in the local search:",numeval,"\n");
cat("CPU Time of the local search:",proc.time()[3]-tic,"seconds \n\n");
}
cpu_time<-proc.time()[3]-cpu_time;
Results$f<-c(Results$f,fbest);
Results$x<-rbind(Results$x, xbest);
Results$neval<-c(Results$neval,nfuneval);
Results$time<-c(Results$time,cpu_time);
if (iterprint){
msg<-switch(fin, "1" = "Maximum number of function evaluations achieved",
"2" = "Maximum computation time achieved", "3"="Desired function value achieved")
cat(msg,"\n");
cat("Best solution value", fbest, "\n");
cat("Decision vector",xbest, "\n");
cat("CPU time", cpu_time,"\n");
cat("Number of function evaluations", nfuneval,"\n");
}
# fclose all;
Results$fbest<-fbest;
Results$xbest<-xbest;
Results$numeval<-nfuneval;
Results$end_crit<-fin;
Results$cpu_time<-cpu_time;
Results$Refset$x<-Refset;
Results$Refset$f<-Refset_values;
Results$Refset$fpen<-Refset_values_penalty;
Results$Refset$const<-Refset_nlc;
Results$Refset$penalty<-penalty;
if (save_results){save(opts, problem, Results, file = "eSSR_report.RData")}
#if plot_results==1 |plot_results==2
# stairs(Results.time,Results.f)
# xlabel('CPU Time (s)');
# ylabel('Objective Function Value (-)');
# title(strcat('Best f: ',num2str(fbest)))
# end
# save ess_report problem opts Results
# %Delete aux files
# if local_solver & isempty(strmatch(local_solver,local_finish));
# ssm_delete_files(local_solver,c_U);
# end
# % keyboard
return(Results)
# end
}
}
}
}
jaegea/MEIGOR documentation built on April 8, 2024, 9:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions essR
Documented in essR
essR <-
function(problem,opts=list(maxeval=NULL,maxtime=NULL),...){
# Function : essR
# R code of the eSS optimization code from: Process Engineering Group IIM-CSIC
# Created on : 09/06/2011
# Last Update: 25/07/2012
# Email : josea.egea@gmail.com
#
# (c) CSIC, Spanish Council for Scientific Research, Spain
#
# Global optimization algorithm for MINLPs based on Scatter Search
#
# essR attempts to solve problems of the form:
# min F(x) subject to: ceq(x) = 0 (equality constraints)
# x c_L <= c(x) <= c_U (inequality constraints)
# x_L <= x <= x_U (bounds on the decision variables)
#
# Constraint functions, if applicable, must be declared in the same script as the
# objective function as a second output argument:
# e.g., myfunction <- function(x){
# calculate fx - scalar containing the objective function value
# calculate gx - vector (or empty) containing the constraints values
# return(list(fx,gx))
# }
#
# Please have a look at the manual before using essR
#
# USAGE: Results <- essR(problem,opts,p1,p2,....,pn);
#
# INPUT PARAMETERS:
# ****************
# problem - List containing problem settings
# problem$f = Name of the file containing the objective
# function (String)
# problem$x_L = Lower bounds of decision variables (vector)
# problem$x_U = Upper bounds of decision variables (vector)
# problem$x_0 = Initial point(s) (optional; vector or matrix)
# problem.$f_0 = Function values of initial point(s) (optional) These values MUST correspond to feasible points.
#
# NOTE:The dimension of f_0 and x_0 may be different. For example, if
# we want to introduce 5 initial points but we only know the values
# for 3 of them, x_0 would have 5 rows whereas f_0 would have only 3
# elements. In this example, it is mandatory that the first 3 rows of x_0 correspond
# to the values of f_0
#
# Additionally, fill the following fields if your problem has
# non-linear constraints
# problem$neq = Number of equality constraints (Integer; do not define it
# if there are no equality constraints)
# problem$c_L = Lower bounds of nonlinear inequality constraints
# (vector)
# problem$c_U = Upper bounds of nonlinear inequality constraints
# (vector)
# problem$int_var = Number of integer variables (Integer)
# problem$bin_var = Number of binary variables (Integer)
# problem$vtr = Objective function value to be reached (optional)
#
#
# NOTE: The order of decision variables is x <- c(cont, int, bin)
#
# opts - List containing options (if set as opts <- numeric(0) default options
# will be loaded). See the script of default options to know their values
#
# User options
# opts$maxeval = Maximum number of function evaluations
# (Default 1000)
# opts$maxtime = Maximum CPU time in seconds (Default 60)
# opts$iterprint = Print each iteration on screen: 0-Deactivated
# 1-Activated (Default 1)
# opts$weight = Weight that multiplies the penalty term added
# to the objective function in constrained
# problems (Default 1e6)
# opts$log_var = Indexes of the variables which will be used
# to generate diverse solutions in different
# orders of magnitude (vector)
# opts$tolc = Maximum absolute violation of the constraints
# (Default 1e-5)
# opts$prob_bound = Probability (0-1) of biasing the search towards
# the bounds (Default 0.5)
# opts$save_results = Saves final results in eSSR_report.RData. (Binary; Default = 1)
# opts$inter_save = Saves results of intermediate iterations in eSSR_report.RData
# Useful for very long runs. (Binary; Default = 0)
#
# Global options
# opts$dim_refset = Number of elements in Refset
# (Integer; automatically calculated)
# opts$ndiverse = Number of solutions generated by the
# diversificator (Default 10*nvar)
# opts$combination = Type of combination of Refset
# elements (Default 1)
# 1: hyper-rectangles
# 2: linear combinations
#
# Local options
# opts$local_solver = Choose local solver
# 0: Local search deactivated (Default),
# "NM", "BFGS", "CG", "LBFGSB",
# "SA","SOLNP", "DHC", "NLS2SOL"
# opts$local_tol = Level of tolerance in local
# search
# opts$local_iterprint = Print each iteration of local
# solver on screen (Binary;
# default = 0)
# opts$local_n1 = Number of iterations
# before applying
# local search for the 1st time
# (Default 1)
# opts$local_n2 = Minimum number of iterations
# in the global
# phase between 2 local calls
# (Default 10)
# opts$local_balance = Balances between quality (=0)
# and diversity (=1) for
# choosing initial points for
# the local search (default
# 0.5)
# opts$local_finish = Applies local search to the
# best solution found once the
# optimization if finished
# (same values as
# opts.local.solver)
# opts$local_bestx = When activated (i.e. =1) only
# applies local search to the
# best solution found to
# date,ignoring filters
# (Default=0)
#
#
# p1,p2... : optional input parameters to be passed to the objective
# function
#
#
# OUTPUT PARAMETERS:
# *****************
# A file called "eSSR_report.Rdata" is generated containing.
#
# Results - Data frame containing results
# problem - Data frame containing problem settings
# opts - Data frame containing all options
#
#
# Fields in Results
# Results$fbest = Best objective function value
# found after the optimization
# Results$xbest = Vector providing the best
# function value
# Results$cpu_time = Time in seconds consumed in the
# optimization
# Results$f = Vector containing the best
# objective function value after each
# iteration
# Results$x = Matrix containing the best vector
# after each iteration
# Results$time = Vector containing the cpu time
# consumed after each iteration
# Results$neval = Vector containing the number of
# function evaluations after each
# iteration
# Results$numeval = Number of function evaluations
# Results$local_solutions = Local solutions found by the
# local solver (in rows)
# Results$local_solutions_values = Function values of the local
# solutions
# Results$end_crit = Criterion to finish the
# optimization
# 1: Maximal number of function
# evaluations achieved
# 2: Maximum allowed CPU Time
# achieved
# 3: Value to reach achieved
#
#
# REFERENCES
#
# If you use essR and publish the results, please cite the following papers:
#
# Egea, J.A., Marti, R., Banga, J.R. (2010) An evolutionary method for
# complex-process optimization. Computers & Operations Research 37(2):
# 315-324.
#
# Egea, J.A., Balsa-Canto, E., Garcia, M.S.G., Banga, J.R. (2009) Dynamic
# optimization of nonlinear processes with an enhanced scatter search
# method. Industrial & Engineering Chemistry Research 49(9): 4388-4401.
# Function : eSS Release 2014A
print('eSS R2014A - Enhanced Scatter Search');
stopOptimization<-0;
# Initialize time
cpu_time=proc.time()[3];
#Take the 3rd component of cpu_time for the total elapsed time
#Initialize the number of function evaluations and a stop variable
nfuneval <- 0;
fin <- 0;
if (!is.numeric(opts$maxeval) && !is.numeric(opts$maxtime)){
cat("WARNING:Either opts.maxeval or opts.maxtime must be defined as a stop criterion \n")
cat("Define at least one of these two options and rerun \n")
Results <- numeric(0);
return(Results)
stop()
} else{
if (!is.numeric(opts$maxeval)){
maxeval <- 1e12;
} else{
maxeval <- opts$maxeval;
}
if (!is.numeric(opts$maxtime)){
maxtime <- 1e12;
} else{
maxtime <- opts$maxtime;
}
if(!(maxtime>0)){
Results <- numeric(0);
return(Results);
}
if(!(maxeval>0)){
Results <- numeric(0);
return(Results);
}
}
cat("\n ------------------------------------------------------------------------------ \n")
cat(" essR - Enhanced Scatter Search in R \n");
cat("<c> IIM-CSIC, Vigo, Spain - email: gingproc@iim.csic.es \n");
cat("------------------------------------------------------------------------------ \n\n");
#Load default values for the options
default <- ssm_defaults();
nargin <- length(as.list(match.call())) -1;
if (nargin<2){opts <- numeric(0)};
#Set all optiions
opts<-ssm_optset(default,opts);
x_U<-problem$x_U;
x_L<-problem$x_L;
#Check if the bounds have the same dimension
if (length(x_U)!= length(x_L)){
cat("Upper and lower bounds have different dimensions!!! \n")
cat("EXITING")
Results<-numeric(0);
return(Results)
stop()
} else{
#Number of decision variables
nvar<-length(x_L);
}
prob_names<-names(problem);
temp<-match("x_0",prob_names);
if (is.na(temp)){x_0<-numeric(0)} else{x_0<-problem$x_0}
temp<-match("f_0",prob_names);
if (is.na(temp)){
f_0<-numeric(0);
fbest<-Inf;
xbest<-rep(Inf,nvar)
} else{
f_0<-problem$f_0;
fbest<-min(f_0);
iii<-which.min(f_0);
xbest<-x_0[iii,];
}
temp<-match("neq",prob_names);
if (is.na(temp)){neq<-0} else{neq<-problem$neq}
temp<-match("c_U",prob_names);
if (is.na(temp)){c_L<-numeric(0); c_U<-numeric(0)} else{c_U<-problem$c_U; c_L<-problem$c_L}
temp<-match("int_var",prob_names);
if (is.na(temp)){int_var<-0} else{int_var<-problem$int_var}
temp<-match("bin_var",prob_names);
if (is.na(temp)){bin_var<-0} else{bin_var<-problem$bin_var}
temp<-match("vtr",prob_names);
if (is.na(temp)){vtr<-numeric(0)} else{vtr<-problem$vtr}
#User options
iterprint=opts$iterprint;
prob_bound=opts$prob_bound;
save_results=opts$save_results;
inter_save=opts$inter_save;
plot_results=opts$plot;
weight=opts$weight;
tolc=opts$tolc;
log_var=opts$log_var;
#Global options
dim_refset=opts$dim_refset;
if (is.character(dim_refset)){
z<-c(-10*nvar/1,-1,1);
nnn<-polyroot(z);
iii<-which(Re(nnn)>0);
dim_refset<-ceiling(Re(nnn[iii]));
#%% is the modulus of a division e.g. 7%%4 provides 3 (like mod(7,4) in Matlab)
if (dim_refset%%2){
dim_refset<-dim_refset+1
if (iterprint){
cat("Refset size automatically calculated:",dim_refset, "\n");
}
}
}
ndiverse<-opts$ndiverse;
if (is.character(ndiverse)){
ndiverse=10*nvar;
if (iterprint){
cat("Number of diverse solutions automatically calculated:",ndiverse,"\n");
}
}
if (ndiverse<dim_refset){ndiverse<-dim_refset}
combin=opts$combination;
#Local options
local_solver=opts$local_solver;
if(is.character(local_solver)){local_solver<-toupper(local_solver)}
local_tol=opts$local_tol;
local_iterprint=opts$local_iterprint;
local_n1<-opts$local_n1;
local_n2<-opts$local_n2;
if (is.character(local_n1)){local_n1<-1}
if (is.character(local_n2)){local_n2<-10}
local_balance=opts$local_balance;
if (!length(opts$local_finish)){
local_finish<-local_solver;
} else{
local_finish=opts$local_finish;
}
local_bestx=opts$local_bestx;
#If there are equality constraints
if (neq){
#Set the new bounds for all the constraints
c_L<-c(rep(0,neq),c_L);
c_U<-c(rep(0,neq),c_U)
}
nconst<-length(c_U);
fobj<-problem$f;
#Check if the objective function has 2 output arguments in constrained problems
if (nconst){
n_out_f<-do.call(fobj,list(x_L,...));
if (length(n_out_f)<2){
cat("For constrained problems the objective function must have at least 2 output arguments \n");
cat("EXITING \n");
Results<-numeric(0);
return(Results)
stop()
}
}
#Change to log
xl_log<-x_L;
xu_log<-x_U;
if (length(log_var)){
aaa<-which(!xl_log[log_var]);
if (length(aaa)){
xl_log[log_var[aaa]]<-1e-10;
}
xl_log[log_var]<-log(xl_log[log_var]);
xu_log[log_var]<-log(xu_log[log_var]);
}
#Check if the objective function has 3 output arguments when NL2SOL is invoked
if(match(local_solver,"NL2SOL", nomatch=0) | match(local_finish,"NL2SOL", nomatch=0)){
n_out_f<-do.call(fobj,list(x_L,...));
nfuneval<-nfuneval+1;
if (length(n_out_f)<3){
cat("NL2SOL requires 3 output arguments: f, g, and the vector of residuals \n");
cat("EXITING \n");
Results<-numeric(0);
return(Results)
stop()
}
}
# TRANSFORM THIS INTO R CODE WHEN THE LOCAL SEARCHES ARE IMPLEMENTED!!!!!!
# %Build auxiliary files in case local search is activated
# if local_solver
# ssm_aux_local(problem.f,x_L,x_U,c_L,c_U,neq,local_solver,nvar,varargin{:});
# end
#Initialize variables
stage_1<-1;
n_minimo<-0;
n_critico<-local_n1;
local_solutions<-numeric(0);
local_solutions_values<-numeric(0);
initial_points<-numeric(0);
Results<-list(f=numeric(0), x=numeric(0), time=numeric(0), neval=numeric(0));
#Possible combinations among the dim_refset elements in Refset, taken in pairs
ncomb<-t(combn(1:dim_refset,2));
MaxSubSet<-(dim_refset^2-dim_refset)/2;
MaxSubSet2<-2*MaxSubSet;
#We generate 5 solutions
solutions<-matrix(0,ndiverse+5,nvar);
solutions[1:5,]<-(matrix(runif(nvar*5),5,nvar)+1:5-1)/5;
solutions[6:(ndiverse+5),]<-matrix(runif(ndiverse*nvar),ndiverse,nvar);
solutions<-solutions*kronecker(matrix(1,ndiverse+5,1),t((xu_log-xl_log)))+kronecker(matrix(1,ndiverse+5,1),t(xl_log));
for (i in 1:(ndiverse+5)){
solutions[i,log_var]=exp(solutions[i,log_var]);
}
l_f_0<-length(f_0);
if (is.vector(x_0) && length(x_0)){l_x_0<-1} else if(is.matrix(x_0)){l_x_0<-nrow(x_0)} else{l_x_0<-0};
if (!length(l_x_0)){l_x_0<-0}
#Put x_0 WITHOUT their corresponding f_0 values at the beginning of solutions
solutions<-rbind(x_0, solutions)
if (l_f_0){solutions<-solutions[-seq(1,l_f_0),]}
temp<-l_x_0-l_f_0;
sol<-matrix(0,ndiverse+5+temp,nvar);
sol_val<-rep(0,ndiverse+5+temp);
sol_val_pen<-rep(0,ndiverse+5+temp);
sol_pen<-rep(0,ndiverse+5+temp);
sol_nlc<-matrix(0,ndiverse+5+temp,nconst);
counter<-1;
for (i in 1:(ndiverse+5+temp)){
res<-ssm_evalfc(solutions[i,],x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
sol[counter,]<-x;
sol_val[counter]<-val;
sol_val_pen[counter]<-val_penalty;
sol_pen[counter]<-pena;
sol_nlc[counter,]<-nlc;
counter<-counter+1;
}
}
nsol<-nrow(sol);
if ((counter-1)<nsol){
sol<-sol[-seq(counter-1,nsol),];
sol_val<-sol_val[-seq(counter-1,nsol)];
sol_val_pen<-sol_val_pen[-seq(counter-1,nsol)];
sol_pen<-sol_pen[-seq(counter-1,nsol)];
sol_nlc<-sol_nlc[-seq(counter-1,nsol),];
}
#Add points with f_0 values. We assume feasiblity
if (l_f_0){
sol<-rbind(x_0[1:l_f_0,], sol);
sol_val<-c(f_0, sol_val);
#We assume feasible points in f_0
sol_val_pen<-c(f_0, sol_val_pen);
sol_pen<-c(rep(0,l_f_0),sol_pen);
#We do not know these values, but we do not care. We assume feasibility
sol_nlc<-rbind(NaN*matrix(0,l_f_0,nconst),sol_nlc);
}
#Initialize Refset
Refset<-matrix(0,dim_refset,nvar);
Refset_values<-rep(0,dim_refset);
Refset_values_penalty<-rep(0,dim_refset);
Refset_nlc<-numeric(0);
penalty<-numeric(0);
temp<-sort(sol_val_pen,index.return=TRUE);
iii<-temp[[2]];
first_members<-ceiling(dim_refset/2);
last_members<-dim_refset-first_members;
#Create the initial sorted Refset
Refset<-sol[iii[1:first_members],];
Refset_values<-sol_val[iii[1:first_members]];
Refset_values_penalty<-sol_val_pen[iii[1:first_members]];
penalty<-sol_pen[iii[1:first_members]];
Refset_nlc<-sol_nlc[iii[1:first_members],];
#The rest of Refset members are chosen randomly
iii<-iii[-(1:first_members)];
temp<-1:length(iii);
ooo<-sample(temp);
#The comment below comes from matlab code
#ooo=randperm(ndiverse+2+size(x_0,1)-first_members);
Refset<-rbind(Refset, sol[iii[ooo[1:last_members]],]);
Refset_values<-c(Refset_values, sol_val[iii[ooo[1:last_members]]]);
Refset_values_penalty<-c(Refset_values_penalty, sol_val_pen[iii[ooo[1:last_members]]]);
penalty<-c(penalty, sol_pen[iii[ooo[1:last_members]]]);
Refset_nlc<-rbind(Refset_nlc, sol_nlc[iii[ooo[1:last_members]],]);
#Check best value
ggg<-which(penalty<=tolc);
if (length(ggg)){
fbest<-min(Refset_values_penalty[ggg]);
iii<-which.min(Refset_values_penalty[ggg]);
xbest<-Refset[ggg[iii],];
}
iter<-0;
if (iterprint){
cat("Initial Pop: NFunEvals:", nfuneval, "Bestf:",fbest,"CPUTime:",proc.time()[3]-cpu_time,"Mean:", mean(Refset_values_penalty),"\n");
}
Results$f=c(Results$f,fbest);
Results$x=rbind(Results$x,xbest);
Results$time=c(Results$time, proc.time()[3]-cpu_time);
Results$neval=c(Results$neval, nfuneval);
#Sort the Refset
temp<-sort(Refset_values_penalty,index.return=TRUE);
Refset_values_penalty<-temp[[1]];
I<-temp[[2]];
Refset_values<-Refset_values[I];
penalty<-penalty[I];
Refset<-Refset[I,];
Refset_nlc<-Refset_nlc[I,];
if (combin==1){nrand=nvar} else if (combin==2) {nrand=1}
use_bestx<-0;
index1<-ncomb[,1];
index2=ncomb[,2];
index<-c(index1,index2);
diff_index<-index2-index1;
#Minimum ppp
lb_p<-1;
#Defines maximum ppp. max(ppp)= lb_p+ st_p
st_p<-0.75;
#Example: lb_p=1.25 and st_p=0.5 varies ppp from 1.25 to 1.75
ppp<-st_p*(diff_index-1)/(dim_refset-2)+lb_p;
hyper_x_L<-kronecker(matrix(1,MaxSubSet,1),t(x_L));
hyper_x_U<-kronecker(matrix(1,MaxSubSet,1),t(x_U));
refset_change<-rep(0,dim_refset);
while (!fin){
child<-matrix(0,MaxSubSet2,nvar);
child_values<-rep(0,MaxSubSet2);
child_values_penalty<-rep(0,MaxSubSet2);
child_penalty<-rep(0,MaxSubSet2);
child_nlc<-matrix(0,MaxSubSet2,nconst);
child_parent_index<-rep(0,MaxSubSet2);
members_update<-rep(0,dim_refset);
counter<-1;
continuar<-0;
while (!continuar){
#Sort Refset
#This is mandatory because the variable index is also sorted
counter<-1;
temp<-sort(Refset_values_penalty,index.return=TRUE);
Refset_values_penalty<-temp[[1]];
I<-temp[[2]];
Refset_values<-Refset_values[I];
penalty<-penalty[I];
Refset<-Refset[I,];
Refset_nlc<-Refset_nlc[I,];
refset_change<-refset_change[I];
parents_index1<-Refset[index1,];
parents_index2<-Refset[index2,];
parents_index1_values_penalty<-Refset_values_penalty[index1];
parents_index2_values_penalty<-Refset_values_penalty[index2];
denom22<-pmax(abs(parents_index1),abs(parents_index2));
denom22[which(!denom22)]<-1;
AAA<-abs((parents_index1-parents_index2)/denom22);
#Compute the maximum value of AAA by rows
temp<-apply(AAA,1,max);
#if any of these maximum values are lower than a threshold, then keep the indexes
BBB<-which(temp<1e-3);
if (length(BBB)){
#Check the index of worst element among the twins of index2
index_refset_out<-max(index2[BBB]);
includ=0;
while (!includ){
#Generate a new random solution
new_refset_member<-runif(nvar)*(xu_log-xl_log)+xl_log;
new_refset_member[log_var]=exp(new_refset_member[log_var]);
res<-ssm_evalfc(new_refset_member,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
Refset[index_refset_out,]<-x;
Refset_values[index_refset_out]<-val;
Refset_values_penalty[index_refset_out]<-val_penalty;
Refset_nlc[index_refset_out,]<-nlc;
penalty[index_refset_out]<-pena;
members_update[index_refset_out]<-0;
}
}
} else{
continuar=1;
}
candidate<-Refset;
candidate_values<-Refset_values;
candidate_values_penalty<-Refset_values_penalty;
candidate_nlc<-Refset_nlc;
candidate_penalty<-penalty;
candidate_update<-rep(0,dim_refset);
members_update<-rep(1,dim_refset);
factor<-kronecker(matrix(1,1,nvar),ppp)*(parents_index2-parents_index1)/1.5;
v1<-parents_index1-factor;
v2<-parents_index2-factor;
v3<-2*parents_index2-parents_index1-factor;
aaa<-which(v1<hyper_x_L);
if (length(aaa)){
rand_aaa<-runif(length(aaa));
AAA<-which(rand_aaa>prob_bound);
aaa<-aaa[AAA];
v1[aaa]=hyper_x_L[aaa];
}
bbb<-which(v1>hyper_x_U);
if (length(bbb)){
rand_bbb<-runif(length(bbb));
BBB<-which(rand_bbb>prob_bound);
bbb<-bbb[BBB];
v1[bbb]=hyper_x_U[bbb];
}
ccc<-which(v3<hyper_x_L);
if (length(ccc)){
rand_ccc<-runif(length(ccc));
CCC<-which(rand_ccc>prob_bound);
ccc<-ccc[CCC];
v3[ccc]=hyper_x_L[ccc];
}
ddd<-which(v3>hyper_x_U);
if (length(ddd)){
rand_ddd<-runif(length(ddd));
DDD<-which(rand_ddd>prob_bound);
ddd<-ddd[DDD];
v3[ddd]=hyper_x_U[ddd];
}
if (nrand>1){
new_comb1<-v1+(v2-v1)*matrix(runif(MaxSubSet*nrand),MaxSubSet,nrand);
new_comb2=v2+(v3-v2)*matrix(runif(MaxSubSet*nrand),MaxSubSet,nrand);
} else{
new_comb1=v1+(v2-v1)*runif(1);
new_comb2=v2+(v3-v2)*runif(1);
}
new_comb<-rbind(new_comb1,new_comb2);
for (i in 1:MaxSubSet2){
#Evaluate
res<-ssm_evalfc(new_comb[i,],x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
child[counter,]<-x;
child_values[counter]<-val;
child_values_penalty[counter]<-val_penalty;
child_penalty[counter]<-pena;
child_nlc[counter,]<-nlc;
child_parent_index[counter]<-index[i];
counter<-counter+1;
if (val_penalty<candidate_values_penalty[index[i]]){
#Update candidate
#This can be accelerated saving the index and replacing at the end with the best one
candidate[index[i],]<-x;
candidate_values_penalty[index[i]]<-val_penalty;
candidate_values[index[i]]<-val;
candidate_penalty[index[i]]<-pena;
candidate_nlc[index[i],]<-nlc;
candidate_update[index[i]]<-1;
members_update[index[i]]<-0;
}
}
}
#Check the stop criterion
if (nfuneval>= maxeval){
fin<-1;
} else if(proc.time()[3]-cpu_time>=maxtime){
fin<-2;
} else if(length(vtr)){
if (fbest<=vtr){
fin<-3;
}
}
child[-(1:(counter-1)),]<-numeric(0);
child_values[-(1:(counter-1))]<-numeric(0);
child_values_penalty[-(1:(counter-1))]<-numeric(0);
child_penalty[-(1:(counter-1))]<-numeric(0);
child_parent_index[-(1:(counter-1))]<-numeric(0);
members_to_update<-which(candidate_update==1);
if (length(members_to_update)){
for (i in 1:length(members_to_update)){
res <- ssm_beyond(Refset[members_to_update[i],],candidate[members_to_update[i],],candidate_values_penalty[members_to_update[i]],
fobj,nrand,tolc,weight,x_L,x_U,c_L,c_U,nconst,int_var,bin_var,nfuneval,prob_bound,nvar,...)
vector<-res[[1]];
vector_value<-res[[2]];
vector_penalty<-res[[3]];
vector_value_penalty<-res[[4]];
vector_nlc<-res[[5]];
new_child<-res[[6]];
new_child_value<-res[[7]];
new_child_penalty<-res[[8]];
new_child_value_penalty<-res[[9]];
new_child_nlc<-res[[10]];
nfuneval<-res[[11]];
if (length(vector)){
Refset[members_to_update[i],]<-vector;
Refset_values[members_to_update[i]]<-vector_value;
Refset_values_penalty[members_to_update[i]]<-vector_value_penalty;
penalty[members_to_update[i]]<-vector_penalty;
Refset_nlc[members_to_update[i],]<-vector_nlc;
} else{
Refset[members_to_update[i],]<-candidate[members_to_update[i],];
Refset_values[members_to_update[i]]<-candidate_values[members_to_update[i]];
Refset_values_penalty[members_to_update[i]]<-candidate_values_penalty[members_to_update[i]];
penalty[members_to_update[i]]<-candidate_penalty[members_to_update[i]];
Refset_nlc[members_to_update[i],]<-candidate_nlc[members_to_update[i],];
}
}
}
fbest_lastiter=fbest;
ggg<-which(penalty<=tolc);
if (length(ggg)){
fbest_new<-min(Refset_values_penalty[ggg]);
hhh<-which.min(Refset_values_penalty[ggg]);
if (fbest_new<fbest){
xbest<-Refset[ggg[hhh],];
fbest<-fbest_new;
use_bestx<-1;
if (inter_save){
Results$fbest<-fbest;
Results$xbest<-xbest;
Results$numeval<-nfuneval;
Results$end_crit<-fin;
Results$cpu_time<-cpu_time;
Results$Refset$x<-Refset;
Results$Refset$f<-Refset_values;
Results$Refset$fpen<-Refset_values_penalty;
Results$Refset$const<-Refset_nlc;
Results$Refset$penalty<-penalty;
if (save_results){save(opts, problem, Results, file = "eSSR_report.RData")}
}
}
}
iter<-iter+1;
n_minimo<-n_minimo+1;
if (iterprint){
cat("Iteration:",iter, "NFunEvals:", nfuneval, "Bestf:",fbest,"CPUTime:",proc.time()[3]-cpu_time,"Mean:", mean(Refset_values_penalty),"\n");
}
ddd<-which(members_update>0);
eee<-which(!members_update);
refset_change[ddd]<-refset_change[ddd]+1;
refset_change[eee]<-0;
fff<-which(refset_change>=20);
if (length(fff)){
to_replace<-length(fff);
replaced<-0;
while (replaced<to_replace){
newsol<-runif(nvar)*(xu_log-xl_log)+xl_log;
newsol[log_var]=exp(newsol[log_var]);
res<-ssm_evalfc(newsol,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ){
Refset[fff[replaced+1],]<-x;
Refset_values[fff[replaced+1]]<-val;
Refset_values_penalty[fff[replaced+1]]<-val_penalty;
penalty[fff[replaced+1]]<-pena;
Refset_nlc[fff[replaced+1],]<-nlc;
refset_change[fff[replaced+1]]<-0;
replaced<-replaced+1;
}
}
}
if (fbest<fbest_lastiter){
Results$f<-c(Results$f,fbest);
Results$x<-rbind(Results$x,xbest);
Results$time<-c(Results$time, proc.time()[3]-cpu_time);
Results$neval<-c(Results$neval,nfuneval);
#if plot_results==1
# stairs(Results.time,Results.f)
# xlabel('CPU Time (s)');
# ylabel('Objective Function Value (-)');
# title(strcat('Best f: ',num2str(fbest)))
# drawnow
# end
}
#******************************************************************************************
#************************* LOCAL FILTERS ROUTINE ******************************************
call_filters<-0;
if(is.character(local_solver) && n_minimo>=n_critico){
call_filters<-1;
if (local_bestx && !use_bestx){
call_filters<-0;
}
}
if (call_filters){
if (stage_1==1){
#In this first stage we apply the local search over the best solution found so far
minimo<-min(Refset_values_penalty);
I<-which.min(Refset_values_penalty);
x0<-Refset[I,];
f0<-minimo;
#Switch to stage 2
n_critico<-local_n2;
if (local_bestx){
#not to go directly to stage 2 the next time we pass over here unless we have really improved fbest
use_bestx<-0;
}
#use_bestx is a flag indicating that the best solution has been improved
}
else if (stage_2==1){
if (local_bestx && use_bestx){
x0<-xbest;
f0<-fbest;
use_bestx<-0;
} else{
#Calculate distance between our children and local_solutions+initial points
local_init<-rbind(local_solutions, initial_points);
xuxl<-x_U-x_L;
#child<-rbind(child,new_child);
# child_values<-c(child_values,new_child_value);
# child_values_penalty<-c(child_values_penalty,new_child_value_penalty);
# child_penalty<-c(child_penalty,new_child_penalty);
# child_nlc<-rbind(child_nlc,new_child_nlc);
child_norm<-child/kronecker(matrix(1,nrow(child),1),t(xuxl));
local_init_norm<-local_init/kronecker(matrix(1,nrow(local_init),1),t(xuxl));
distance<-eucl_dist(child_norm,local_init_norm);
#Now we need to sort the children according the their distances to local minima and intial points
#Here we compute the minimum euclidean distance between child i and all the vector in initial_points+local_solutions
temp<-apply(distance,1,min);
#Now we sort those minimum distances and favour the highest ones to favour diverse points
res<-sort(temp,decreasing=TRUE,index.return=TRUE);
ooo<-res[[2]];
res<-sort(child_values_penalty,index.return=TRUE);
uuu<-res[[2]];
n_child<-length(child_values_penalty)
child_points<-rep(0,n_child);
www<-local_balance;
for (i in 1:n_child){
child_points[ooo[i]]<-child_points[ooo[i]]+www*i;
child_points[uuu[i]]<-child_points[uuu[i]]+(1-www)*i;
}
I<-which.min(child_points);
x0<-child[I,];
f0<-child_values_penalty[I];
}
}
if (iterprint){
cat("Call local solver:", toupper(local_solver),"\n")
if (child_penalty[I]>tolc){
cat("Initial point function value:",child_values[I], "(INFEASIBLE)","\n");
} else{
cat("Initial point function value:",f0,"\n");
}
tic<-proc.time()[3];
}
#AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
res<-ssm_localsolver(x0,x_L,x_U,c_L,c_U,neq,int_var,bin_var,fobj,local_solver,local_iterprint,local_tol,weight,nconst,tolc,...);
x<-res[[1]];
fval<-res[[2]];
numeval<-res[[3]];
initial_points<-rbind(initial_points, x0);
nfuneval<-nfuneval+numeval;
#Evaluate the local solution
res<-ssm_evalfc(x,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ && pena<=tolc){
if (iterprint){
cat("Local solution function value:",val,"\n");
}
if (stage_1==1){
stage_1=0;
stage_2=1;
}
if (val_penalty<fbest){
fbest<-val_penalty;
xbest<-x;
if (fbest<fbest_lastiter){
Results$f<-c(Results$f,fbest);
Results$x<-rbind(Results$x,xbest);
Results$time<-c(Results$time,proc.time()[3]-cpu_time);
Results$neval<-c(Results$neval,nfuneval);
#if plot_results==1
# stairs(Results.time,Results.f)
# xlabel('CPU Time (s)');
# ylabel('Objective Function Value (-)');
# title(strcat('Best f: ',num2str(fbest)))
# drawnow
# end
}
}
} else{
cat("Local solution function value:","UNFEASIBLE","\n");
}
cat("Number of function evaluations in the local search:",numeval,"\n");
cat("CPU Time of the local search:",proc.time()[3]-tic,"seconds \n\n");
#Check the stopping criterion
if (nfuneval>= maxeval){
fin<-1;
} else if(proc.time()[3]-cpu_time>=maxtime){
fin<-2;
} else if(length(vtr)){
if (fbest<=vtr){
fin<-3;
}
}
if (includ){
if (length(local_solutions)){
#Add the found local solution to the list (if it is a new one)
adicionar_local<-1;
res<-ssm_isdif2(x,local_solutions,1e-2,1);
f<-res[[1]];
ind<-res[[2]];
#If it is too close to a previously found local solution, we do not add it
if (f){
adicionar_local<-0
}
} else {
#Add it if no local solutions have been found yet
adicionar_local<-1;
}
if (adicionar_local){
#choose either this (replace the worst element in Refset)
# [aaa,jjj]=max(Refset_values_penalty);
#
# if val_penalty<Refset_values_penalty(jjj)
# Refset(jjj,:)=x;
# Refset_values(jjj)=val;
# Refset_values_penalty(jjj)=val_penalty;
# Refset_nlc(jjj,:)=nlc;
# % penalty(jjj)=pena;
# % end
#
# %or this (Replace the initial point)
#
# % if val_penalty<f0
# % child(I,:)=x;
# % child_values(I)=val;
# % child_values_penalty(I)=val_penalty;
# % child_nlc(I,:)=nlc;
# % end
#
#or this (Replace the parent of the initial point)
if (val_penalty<Refset_values_penalty[child_parent_index[I]]){
Refset[child_parent_index[I],]<-x;
Refset_values[child_parent_index[I]]<-val;
Refset_values_penalty[child_parent_index[I]]<-val_penalty;
Refset_nlc[child_parent_index[I],]<-nlc;
penalty[child_parent_index[I]]<-pena;
refset_change[child_parent_index[I]]<-0;
}
local_solutions<-rbind(local_solutions,x);
local_solutions_values<-c(local_solutions_values,fval);
}
}
n_minimo<-0
}
#************************* END OF LOCAL FILTERS ROUTINE ***********************************
#******************************************************************************************
fbest_lastiter<-fbest;
#Check the stopping criterion
if (nfuneval>= maxeval){
fin<-1;
} else if(proc.time()[3]-cpu_time>=maxtime){
fin<-2;
} else if(length(vtr)){
if (fbest<=vtr){
fin<-3;
}
}
if (fin){
cat("\n");
if (is.character(local_finish) && fin<3){
local_tol=3;
if (is.numeric(xbest)){
x0=xbest;
f0=fbest;
}else{
x0=Refset[1,];
f0=Refset_values_penalty[1];
}
if (iterprint){
cat("Final local refinement: ", toupper(local_finish),"\n")
cat("Initial point function value: ", f0,"\n");
tic<-proc.time()[3];
}
res<-ssm_localsolver(x0,x_L,x_U,c_L,c_U,neq,int_var,bin_var,fobj,local_finish,local_iterprint,local_tol,weight,nconst,tolc,...);
x<-res[[1]];
fval<-res[[2]];
numeval<-res[[3]];
#Evaluate the local solution
res<-ssm_evalfc(x,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
val<-res[[1]];
val_penalty<-res[[2]];
pena<-res[[3]];
nlc<-res[[4]];
includ<-res[[5]];
x<-res[[6]];
nfuneval<-nfuneval+1;
if (includ && pena<=tolc){
if (iterprint){
cat("Local solution function value:",val,"\n");
}
if (val_penalty<fbest){
fbest<-val_penalty;
xbest<-x;
}
} else{
cat("Local solution function value:","UNFEASIBLE","\n");
}
cat("Number of function evaluations in the local search:",numeval,"\n");
cat("CPU Time of the local search:",proc.time()[3]-tic,"seconds \n\n");
}
cpu_time<-proc.time()[3]-cpu_time;
Results$f<-c(Results$f,fbest);
Results$x<-rbind(Results$x, xbest);
Results$neval<-c(Results$neval,nfuneval);
Results$time<-c(Results$time,cpu_time);
if (iterprint){
msg<-switch(fin, "1" = "Maximum number of function evaluations achieved",
"2" = "Maximum computation time achieved", "3"="Desired function value achieved")
cat(msg,"\n");
cat("Best solution value", fbest, "\n");
cat("Decision vector",xbest, "\n");
cat("CPU time", cpu_time,"\n");
cat("Number of function evaluations", nfuneval,"\n");
}
# fclose all;
Results$fbest<-fbest;
Results$xbest<-xbest;
Results$numeval<-nfuneval;
Results$end_crit<-fin;
Results$cpu_time<-cpu_time;
Results$Refset$x<-Refset;
Results$Refset$f<-Refset_values;
Results$Refset$fpen<-Refset_values_penalty;
Results$Refset$const<-Refset_nlc;
Results$Refset$penalty<-penalty;
if (save_results){save(opts, problem, Results, file = "eSSR_report.RData")}
#if plot_results==1 |plot_results==2
# stairs(Results.time,Results.f)
# xlabel('CPU Time (s)');
# ylabel('Objective Function Value (-)');
# title(strcat('Best f: ',num2str(fbest)))
# end
# save ess_report problem opts Results
# %Delete aux files
# if local_solver & isempty(strmatch(local_solver,local_finish));
# ssm_delete_files(local_solver,c_U);
# end
# % keyboard
return(Results)
# end
}
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.