Exemplo de Recozimento Simulado (Simulated Annealing)

> restart;

> with(plots):
readlib(randomize);
randomize();

Vamos empregar o algoritmo de recozimento simulado para achar os mínimos da função F(x) abaixo.

> F:=k->(k+10)*(k+6)*(k+5)*(k+1)*(k-7)*(k-10)/9000;
plot(F(x),x=-11..11,-10..30,axes=boxed);

Funções auxiliares: a( ) calcula número real aleatório entre -1 e 1.
b( ) calcula número real aleatório entre 0 e 1.

> a:=rand(-1000..1000)/1000.;
b:=rand(0..1000)/1000.;

Esquema de resfriamento

> resfr:=T->alpha*T;

Esquema de Vizinhança: um vizinho é qualquer ponto que diste menos de dk do ponto atual

> dk:=1.0;
viz:=x->x+dk*a();

Definindo parâmetros ( sugestão: T=80 =0,75 n1=20 n2=50 )

> T0:=80; # temperatura inicial
alpha:=0.75; # parametro da função de resfriamento
k0:=-11; # chute inicial
n1:=15; # numero de reducoes de temperatura
n2:=60; #numero de iteracoes com a mesma temperatura dk:=1;
sols:=array(1..n1*n2);
bests:=array(1..n1*n2);

Algoritmo de Recozimento Simulado

> k:=k0:
T:=T0:
best:=evalf(k);

for i from 1 to n1 do
for j from 1 to n2 do

# computo de solucao vizinha
knew:=viz(k):
delta:=F(knew)-F(k):
if (evalf(delta)<0) then
k:=knew:
else
temp:=b()-exp(-delta/T):
if (evalf(temp)<0) then
k:=knew:
fi:
fi:

# armazenando o melhor ponto ja visitado
if (evalf(F(k))<evalf(F(best))) then
best:=evalf(k);
fi;
# armazenando solucoes intermediarias
# (apenas para fazer graficos depois!)
sols[(i-1)*n2+j]:=evalf(k):
bests[(i-1)*n2+j]:=evalf(best):

od:

# reducao da temperatura
T:=evalf(resfr(T)):

od:

RESULTADO FINAL:

Ultimo ponto visitado => azul
Melhor ponto visitado => verde

> plot([F(x),[[sols[n1*n2], F(sols[n1*n2])]],[[best, F(best)]]],x=-11..11,
style=[line,point,point],symbol=[circle,cross,box],axes=boxed,
color=[red,blue,green]);

Pontos visitados ao longo da execução do algoritmo.
verde: solução tentativa daquela iteração
vermelha: melhor já visitada

> plot([[seq([k,bests[k]],k=1..n1*n2)],[seq([k,sols[k]],k=1..n1*n2)]],x=1..n1*n2,-15..15,axes=boxed,color=[red,green],labels=[`iteracao`,``]);

O Gráfico abaixo mostra a evolução do valor da solução atual.
Observe que a amplitude de variação da função diminui com a redução da temperatura.

> plot([[seq([k,F(bests[k])],k=1..n1*n2)],[seq([k,F(sols[k])],k=1..n1*n2)]],axes=boxed,color=[red,green],labels=[`iteracao`,`F(k)`]);

A animação seguinte mostra o caminho percorrido ao longo da execução do algoritmo.
preto: chute inicial
azul: pontos visitados durante a execução do algoritmo

verde: solução final

> S1:=seq(plot([F(x),[[sols[1], F(sols[1])]]],x=-13..11,-10..30,
style=[line,point],symbol=circle,axes=boxed,
color=[red,black]), kk=1..n1):

> S2:=seq(plot([F(x),[[sols[kk], F(sols[kk])]]],x=-13..11,-10..30,
style=[line,point],symbol=circle,axes=boxed,
color=[red,blue]), kk=1..n1*n2):

> S3:=seq(plot([F(x),[[sols[n1*n2], F(sols[n1*n2])]]],x=-13..11,-10..30,
style=[line,point],symbol=circle,axes=boxed,
color=[red,green]), kk=1..n1):

> display([S1,S2,S3],insequence=true);

>

>

>