Algoritmo para encontrar máximo ou mínimo de funções

Computação Completa

•

UNB

Erique Pereira Neto

13/01/2018

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

##########################
### ###
### TRABALHO PRATICO ###
### 4 ###
##########################
#######################################################
##########################
### Metodo Gradiente ###
##########################
a=Sys.time()
### Calculo Numerico - Gradiente
g <- function(f,x) {
 deltha <- 10^(-4)
 e <- diag(length(x))
 g <- c((f(x + deltha*e[,1]) - f(x))/deltha,(f(x + deltha*e[,2]) - f(x))/deltha)
 return(g)
}
### Secao Aurea 
aurea <- function(f, x, d) {
 a <- 0
 b <- 100
 ep <- 0.0001
 golden.ratio <- 2/(sqrt(5) + 1)
 
 while((b - a) > ep) {
 xa <- b - golden.ratio*(b - a)
 xb <- a + golden.ratio*(b - a)
 fa <- f(x + xa*d)
 fb <- f(x + xb*d)
 if(fa < fb){
 b <- xb
 }else{
 a <- xa
 }
 }
 est.min <- (a + b)/2
 return(list(est.min = est.min))
}
### Método Gradiente
gradiente <- function(f,x,erro) {
 k <- 0
 grad <- g(f,x)
 ep <- erro
 
 criterio.Parada <- ifelse(sqrt(grad[1]^2+grad[2]^2)<ep, T, F)
 
 
 while(sum(criterio.Parada)!=1 ){
 grad <- g(f,x)
 d <- -grad
 alpha <- aurea(f,x, d)$est.min
 xk=x
 x <- xk+alpha*d
 
 
 #df=(1/erro)*abs((f(x)-f(xk))/f(xk))
 #criterio.Parada <- ifelse(df<1, T, F)
 criterio.Parada <- ifelse(sqrt(grad[1]^2+grad[2]^2)<ep, T, F)
 y=f(x)
 k <- k+1
 if(k>999){
 break()
 }
 }
 return(list(X_Otimo=x,F_Otimo=y,Iteracoes=k))
}
##########################
### Metodo BFGS ###
##########################
### Calculo Numerico - Hessiana
#gh <- function(x) {
# deltha <- 10^(-4)
# h <- diag(length(x))
# gh <- c((f(x + deltha*e[,1]) - f(x))/deltha,(f(x + deltha*e[,2]) - f(x))/deltha)
# return(gh)
#}
### Calculo Numerico - Hessiana
#gh <- function(x) {
# deltha <- 10^(-4)
# e <- diag(length(x))
### Método BFGS
bfgs <- function(f, x,erro){
 H <- diag(length(x)) # Identidade
 grad <- g(f, x)
 k <- 0
 ep <- erro
 criterio.Parada <- ifelse(sqrt(grad[1]^2+grad[2]^2)<ep, T, F)
 while(sum(criterio.Parada)!=1){
 d <- -H%*%grad 
 golden <- aurea(f, x, d)
 xk=x
 x <- xk + (golden$est.min)*d
 gradk=grad
 grad <- g(f, x)
 V <- xk - x
 R <- matrix(gradk - grad)
 part1 <- as.numeric(1 + ((t(R)%*%H%*%R)/(t(R)%*%V)))*((V%*%t(V)))/as.numeric((t(V)%*%R))
 part2 <- ((V%*%(t(R))%*%H)+(H%*%R%*%(t(V))))/as.numeric(t(R)%*%V)
 Hk=H
 H <- Hk + part1 - part2
 #df=(1/erro)*abs((f(x)-f(xk))/f(xk))
 #criterio.Parada <- ifelse(df<1, T, F)
 criterio.Parada <- ifelse(sqrt(grad[1]^2+grad[2]^2)<ep, T, F)
 k <- k + 1 
 
 if(k>999){
 break()
 }
 
 y=f(x)
 }
 return(list(X_Otimo=x,F_Otimo=y,Iteracoes=k))
}
#PSO method
pso=function(x,ITER,erro,NPART){
 
 X=x
 Z=list(rep(0,NPART))
 Y=rep(Z,10000)
 j=1
 W=c()
 
 
 g=X[which.min(Y[[j]]),]
 criterio.Parada=ifelse(pi<0,T,F)
 
 while(sum(criterio.Parada)!=1 ){
 
 for(i in 1:NPART){
 
 Y[[j]][i]=f(as.vector(X[i,1:2]))
 
 }
 
 if(j==1){
 p=X
 g=X[which.min(Y[[j]]),]
 V=matrix(0,NPART,2)
 Vi=V
 gb=matrix(g,NPART,2)
 for(k in 1:NPART){
 gb[k,1:2]=g
 }
 V=0.73*(0.3*Vi+runif(1)*2.05*(p-X)+runif(1)*2.05*(gb-X))
 Xi=X
 X=Xi+V
 W[j]=f(g)
 }else{
 px=matrix(0,NPART,2)
 for(i in 1:NPART){
 if((f(p[i,])-f(X[i,]))<0){px[i,]=p[i,]}else{px[i,]=X[i,]}
 }
 p=px
 gk=g
 gx=X[which.min(Y[[j]]),]
 if(f(gx)-f(g)<0){g=gx}else{g=g}
 
 gb=matrix(g,NPART,2)
 for(k in 1:NPART){
 gb[k,1:2]=g
 }
 Vi=V
 V=0.73*(0.3*Vi+runif(1)*2.05*(p-X)+runif(1)*2.05*(gb-X))
 Xi=X
 X=Xi+V
 
 }
 W[j]=f(g)
 
 
 if(j>ITER){
 maxj=max(W)
 minj=min(W)
 maxj5=max(W[(j-5):j])
 minj5=min(W[(j-5):j])
 
 var.global=maxj-minj
 varj5=maxj5-minj5
 vard=varj5/var.global
 }
 
 criterio.Parada=ifelse(j>ITER,ifelse((1/erro)*vard<1,T,F),F)
 j=j+1
 
 }
 saida=list(X_Otimo=g,F_otimo=W[j-1],Iteracoes=j-1)
 return(saida)
}
# gerador de vetores aleatorios
gera_vet_aleatorio=function(a,b,n){
 
 x1=(b-a)*runif(n)+a
 
 x2=(b-a)*runif(n)+a
 
 matrix(c(x1,x2),n,2)
}
# ERRO
erro=10^(-4)
NPART=1000
### Definindo funçoes para otimizar
# Funcao Quadratica
f <- function(x) {
 x[1]^2 + x[2]^2
}
X01_grad=gera_vet_aleatorio(-100,100,100)
Aplic_Gradiente=list()
for(j in 1:100){
 Aplic_Gradiente[[j]]=gradiente(f,X01_grad[j,],erro)
}
Aplic_BFGS=list()
for(j in 1:100){
 Aplic_BFGS[[j]]=bfgs(f,X01_grad[j,],erro)
}
Aplic_PSO=list()
for(j in 1:100){
 
 Aplic_PSO[[j]]=pso(gera_vet_aleatorio(-100,100,NPART),19,erro,NPART)
}
#Funcao Quadratica Malcondicionada
f <- function(x) {
 x[1]^2 + 100*x[2]^2
}
X02_grad=gera_vet_aleatorio(-100,100,100)
Aplic_Gradiente2=list()
for(j in 1:100){
 Aplic_Gradiente2[[j]]=gradiente(f,X02_grad[j,],erro)
}
Aplic_BFGS2=list()
for(j in 1:100){
 Aplic_BFGS2[[j]]=bfgs(f,X02_grad[j,],erro)
}
Aplic_PSO2=list()
for(j in 1:100){
 
 Aplic_PSO2[[j]]=pso(gera_vet_aleatorio(-100,100,NPART),19,erro,NPART)
}
# Funcao Nao-Diferenciavel ### cuidado com a f3
f <- function(x) {
 max(c((x[1]^2 - 1)^2 + x[2]^2,(x[1]^2 + 1)^2 + x[2]^2))
}
X03_grad=gera_vet_aleatorio(-100,100,100)
Aplic_Gradiente3=list()
for(j in 1:100){
 Aplic_Gradiente3[[j]]=gradiente(f,X03_grad[j,],erro)
}
Aplic_BFGS3=list()
for(j in 1:100){
 Aplic_BFGS3[[j]]=bfgs(f,X03_grad[j,],erro)
}
Aplic_PSO3=list()
for(j in 1:100){
 
 Aplic_PSO3[[j]]=pso(gera_vet_aleatorio(-100,100,NPART),19,erro,NPART)
}
# Funcao Rosenbrock
f <- function(x) {
 (1- x[1])^2 + 100*(x[2] - x[1]^2)^2
}
X04_grad=gera_vet_aleatorio(-5,5,100)
Aplic_Gradiente4=list()
for(j in 1:100){
 Aplic_Gradiente4[[j]]=gradiente(f,X04_grad[j,],erro)
}
Aplic_BFGS4=list()
for(j in 1:100){
 Aplic_BFGS4[[j]]=bfgs(f,X04_grad[j,],erro)
}
Aplic_PSO4=list()
for(j in 1:100){
 
 Aplic_PSO4[[j]]=pso(gera_vet_aleatorio(-5,5,NPART),19,erro,NPART)
}
# Funcao Rastrigin
f <- function(x) {
 20 + (x[1]^2 - 10*cos(2*pi*x[1])) + (x[2]^2 - 10*cos(2*pi*x[2])) 
}
X05_grad=gera_vet_aleatorio(-5.12,5.12,100)
Aplic_Gradiente5=list()
for(j in 1:100){
 Aplic_Gradiente5[[j]]=gradiente(f,X05_grad[j,],erro)
}
Aplic_BFGS5=list()
for(j in 1:100){
 Aplic_BFGS5[[j]]=bfgs(f,X05_grad[j,],erro)
}
Aplic_PSO5=list()
for(j in 1:100){
 
 Aplic_PSO5[[j]]=pso(gera_vet_aleatorio(-5.12,5.12,NPART),19,erro,NPART)
}
distancia=function(v1,v2){
 
 D=sqrt((v2[1]-v1[1])^2+(v2[2]-v1[2])^2)
 D
}
v1=c(0,0)
v2=Aplic_Gradiente[[1]]$X_Otimo
distancia(v1,v2)
D1Gf1=c()
for(i in 1:100){
 v1=c(0,0)
 D1Gf1[i]=distancia(v1,Aplic_Gradiente[[i]]$X_Otimo)
}
D1Bf1=c()
for(i in 1:100){
 v1=c(0,0)
 D1Bf1[i]=distancia(v1,Aplic_BFGS[[i]]$X_Otimo)
}
D1Pf1=c()
for(i in 1:100){
 v1=c(0,0)
 D1Pf1[i]=distancia(v1,Aplic_PSO[[i]]$X_Otimo)
}
D2Gf2=c()
for(i in 1:100){
 v1=c(0,0)
 D2Gf2[i]=distancia(v1,Aplic_Gradiente2[[i]]$X_Otimo)
}
D2Bf2=c()
for(i in 1:100){
 v1=c(0,0)
 D2Bf2[i]=distancia(v1,Aplic_BFGS2[[i]]$X_Otimo)
}
D2Pf2=c()
for(i in 1:100){
 v1=c(0,0)
 D2Pf2[i]=distancia(v1,Aplic_PSO2[[i]]$X_Otimo)
}
D3Gf3=c()
for(i in 1:100){
 v1=c(0,0)
 D3Gf3[i]=distancia(v1,Aplic_Gradiente3[[i]]$X_Otimo)
}
D3Bf3=c()
for(i in 1:100){
 v1=c(0,0)
 D3Bf3[i]=distancia(v1,Aplic_BFGS3[[i]]$X_Otimo)
}
D3Pf3=c()
for(i in 1:100){
 v1=c(0,0)
 D3Pf3[i]=distancia(v1,Aplic_PSO3[[i]]$X_Otimo)
}
D4Gf4=c()
for(i in 1:100){
 v1=c(1,1)
 D4Gf4[i]=distancia(v1,Aplic_Gradiente4[[i]]$X_Otimo)
}
D4Bf4=c()
for(i
in 1:100){
 v1=c(1,1)
 D4Bf4[i]=distancia(v1,Aplic_BFGS4[[i]]$X_Otimo)
}
D4Pf4=c()
for(i in 1:100){
 v1=c(1,1)
 D4Pf4[i]=distancia(v1,Aplic_PSO4[[i]]$X_Otimo)
}
D5Gf5=c()
for(i in 1:100){
 v1=c(0,0)
 D5Gf5[i]=distancia(v1,Aplic_Gradiente5[[i]]$X_Otimo)
}
D5Bf5=c()
for(i in 1:100){
 v1=c(0,0)
 D5Bf5[i]=distancia(v1,Aplic_BFGS5[[i]]$X_Otimo)
}
D5Pf5=c()
for(i in 1:100){
 v1=c(0,0)
 D5Pf5[i]=distancia(v1,Aplic_PSO5[[i]]$X_Otimo)
}
convergiu_Gf1=which(D1Gf1<=10^(-3))
convergiu_Gf2=which(D2Gf2<=10^(-3))
convergiu_Gf3=which(D3Gf3<=10^(-3))
convergiu_Gf4=which(D4Gf4<=10^(-3))
convergiu_Gf5=which(D5Gf5<=10^(-3))
convergiu_Bf1=which(D1Bf1<=10^(-3))
convergiu_Bf2=which(D2Bf2<=10^(-3))
convergiu_Bf3=which(D3Bf3<=10^(-3))
convergiu_Bf4=which(D4Bf4<=10^(-3))
convergiu_Bf5=which(D5Bf5<=10^(-3))
convergiu_Pf1=which(D1Pf1<=10^(-3))
convergiu_Pf2=which(D2Pf2<=10^(-3))
convergiu_Pf3=which(D3Pf3<=10^(-3))
convergiu_Pf4=which(D4Pf4<=10^(-3))
convergiu_Pf5=which(D5Pf5<=10^(-3))
#Numero de iteracoes
Iter_GF1=c()
for (i in convergiu_Gf1){
 Iter_GF1[i]=Aplic_Gradiente[[i]]$Iteracoes
}
Iter_GF2=c()
for (i in as.integer(convergiu_Gf2)){
 Iter_GF2[i]=Aplic_Gradiente2[[i]]$Iteracoes
}
Iter_GF3=c()
for (i in convergiu_Gf3){
 Iter_GF3[i]=Aplic_Gradiente3[[i]]$Iteracoes
}
Iter_GF4=c()
for (i in convergiu_Gf4){
 Iter_GF4[i]=Aplic_Gradiente4[[i]]$Iteracoes
}
Iter_GF5=c()
for (i in convergiu_Gf5){
 Iter_GF5[i]=Aplic_Gradiente5[[i]]$Iteracoes
}
Iter_BF1=c()
for (i in convergiu_Bf1){
 Iter_BF1[i]=Aplic_BFGS[[i]]$Iteracoes
}
Iter_BF2=c()
for (i in convergiu_Bf2){
 Iter_BF2[i]=Aplic_BFGS2[[i]]$Iteracoes
}
Iter_BF3=c()
for (i in convergiu_Bf3){
 Iter_BF3[i]=Aplic_BFGS3[[i]]$Iteracoes
}
Iter_BF4=c()
for (i in convergiu_Bf4){
 Iter_BF4[i]=Aplic_BFGS4[[i]]$Iteracoes
}
Iter_BF5=c()
for (i in convergiu_Bf5){
 Iter_BF5[i]=Aplic_BFGS5[[i]]$Iteracoes
}
Iter_PF1=c()
for (i in convergiu_Pf1){
 Iter_PF1[i]=Aplic_PSO[[i]]$Iteracoes
}
Iter_PF2=c()
for (i in convergiu_Pf2){
 Iter_PF2[i]=Aplic_PSO2[[i]]$Iteracoes
}
Iter_PF3=c()
for (i in convergiu_Pf3){
 Iter_PF3[i]=Aplic_PSO3[[i]]$Iteracoes
}
Iter_PF4=c()
for (i in convergiu_Pf4){
 Iter_PF4[i]=Aplic_PSO4[[i]]$Iteracoes
}
Iter_PF5=c()
for (i in convergiu_Pf5){
 Iter_PF5[i]=Aplic_PSO5[[i]]$Iteracoes
}
# tAXa de convergência
T_Gf1=sum(Iter_GF1>0 ,na.rm=T)/100
T_Gf2=sum(Iter_GF2>0,na.rm = T)/100
T_Gf3=sum(Iter_GF3>0,na.rm = T)/100
T_Gf4=sum(Iter_GF4>0,na.rm = T)/100
T_Gf5=sum(Iter_GF5>0,na.rm = T)/100
T_Bf1=sum(Iter_BF1>0,na.rm=T)/100
T_Bf2=sum(Iter_BF2>0,na.rm = T)/100
T_Bf3=sum(Iter_BF3>0,na.rm = T)/100
T_Bf4=sum(Iter_BF4>0,na.rm = T)/100
T_Bf5=sum(Iter_BF5>0,na.rm = T)/100
T_Pf1=sum(Iter_PF1>0,na.rm=T)/100
T_Pf2=sum(Iter_PF2>0,na.rm = T)/100
T_Pf3=sum(Iter_PF3>0,na.rm = T)/100
T_Pf4=sum(Iter_PF4>0,na.rm = T)/100
T_Pf5=sum(Iter_PF5>0,na.rm = T)/100
# Número de iteracoes
M_Gf1=mean(Iter_GF1[Iter_GF1>0 & Iter_GF1<5],na.rm=T)
M_Gf2=mean(Iter_GF2,na.rm = T)
M_Gf3=mean(Iter_GF3,na.rm = T)
M_Gf4=mean(Iter_GF4,na.rm = T)
M_Gf5=mean(Iter_GF5,na.rm = T)
M_Bf1=mean(Iter_BF1[Iter_BF1>0 & Iter_BF1<5],na.rm=T)
M_Bf2=mean(Iter_BF2,na.rm = T)
M_Bf3=mean(Iter_BF3,na.rm = T)
M_Bf4=mean(Iter_BF4,na.rm = T)
M_Bf5=mean(Iter_BF5,na.rm = T)
M_Pf1=mean(Iter_PF1,na.rm=T)
M_Pf2=mean(Iter_PF2,na.rm = T)
M_Pf3=mean(Iter_PF3,na.rm = T)
M_Pf4=mean(Iter_PF4,na.rm = T)
M_Pf5=mean(Iter_PF5,na.rm = T)
metodos=c("GRADIENTE","BFGS","PSO")
funcoes=c("Quadratica","Quadratica_Malcond.","Nao_Diferenciavel","Rosenbrock","Rastrigin")
G=100*c(T_Gf1,T_Gf2,T_Gf3,T_Gf4,T_Gf5 )
B=100*c(T_Bf1,T_Bf2,T_Bf3,T_Bf4,T_Bf5 )
P=100*c(T_Pf1,T_Pf2,T_Pf3,T_Pf4,T_Pf5 )
Taxa_CONVERGENCIA=data.frame(G,B,P)
colnames(Taxa_CONVERGENCIA)=metodos
row.names(Taxa_CONVERGENCIA)=funcoes
#Media de iteracoes
GM=c(M_Gf1,M_Gf2,M_Gf3,M_Gf4,M_Gf5 )
BM=c(M_Bf1,M_Bf2,M_Bf3,M_Bf4,M_Bf5 )
PM=c(M_Pf1,M_Pf2,M_Pf3,M_Pf4,M_Pf5 )
NUMERO_ITERACOES=data.frame(GM,BM,PM)
colnames(NUMERO_ITERACOES)=metodos
row.names(NUMERO_ITERACOES)=funcoes
#Avaliacoes da funcao objetivo
GMf=(2+1)*c(M_Gf1,M_Gf2,M_Gf3,M_Gf4,M_Gf5 )
BMf=(2+1)^2*c(M_Bf1,M_Bf2,M_Bf3,M_Bf4,M_Bf5 )
PMf=NPART*c(M_Pf1,M_Pf2,M_Pf3,M_Pf4,M_Pf5 )
NUMERO_AVAL_FO=data.frame(GMf,BMf,PMf)
colnames(NUMERO_AVAL_FO)=metodos
row.names(NUMERO_AVAL_FO)=funcoes
Taxa_CONVERGENCIA
NUMERO_ITERACOES
NUMERO_AVAL_FO
b=Sys.time()
b-a
#save(Taxa_CONVERGENCIA,NUMERO_ITERACOES,NUMERO_AVAL_FO, file = "tabelas.RData")