Universidade Estadual de Campinas
FACULDADE DE ENGENHARIA ELÉTRICA E COMPUTAÇÃO
Otimização Não-Linear
BUSCA LINEAR: MÉTODOS PARA BUSCA UNIDIMENSIONAL
Autor: Tiago Agostinho de Almeida
INTRODUÇÃO Considere o seguinte problema irrestrito:
minimizar ƒ(x), x Є R n Muitos algoritmos de otimização irrestrita assumem a seguinte estrutura geral: Escolha ε > 0, xo e faça k = 0
while || ∇ f ( x k ) ||≥ ε 1: Encontre dk Є R n tal que ∇ f ( x k )T d k < 0 2: Determine α k := arg min 0 f ( x k + α d k ) 3: Faça x k 1 := x k + α k d k e k:= k + 1 end x* = xk α >
+
Notas 1. No passo 1 determina-se uma direção de descida. A condição ∇ f ( x k )T d k < 0 garante que f decresce na direção dk a partir de xk para α > 0 2. No passo 2 é realizada um busca linear para encontrar o tamanho do passo α k que minimiza f na direção dk a partir de xk 3. No passo 3 um novo ponto é calculado. O critério de parada é || ∇ f ( x ) ||< ε (idealmente, || ∇ f ( x ) ||= 0 ) k
k
Exemplo No método do gradiente, se || ∇ f ( x k ) ||≥ ε , então d := −∇ f ( x ) é tal que T ∇ f ( x ) d < 0 e a busca linear assume a forma k
k
k
k
α
k
:= arg min
α k >
(
k
( k ))
0 f x − α ∇ f x
Existem vários métodos para encontrar o valor de k k α := arg min k 0 f ( x − α ∇ f ( x )) , dois deles foram abordados e testados neste trabalho, eles são: o Método da Falsa Posição e o Método da Seção Áurea. α k >
MÉTODO DA FALSA POSIÇÃO Dada uma função f ( x ) , o método da falsa posição procura encontrar um tamanho de passo α tal que α k := arg min 0 f ( x k + α d k ) em uma determinada direção d . Em primeiro lugar, devemos encontrar a direção d para “caminharmos” de tal forma que melhore o valor da função objetivo. Então, como sabemos que o gradiente da função indica a direção de máximo crescimento, se queremos minimizar f ( x ) devemos caminhar na direção contrária ao gradiente, sendo assim d := −∇ f ( x ) indica a direção de máximo decrescimento da função na direção d . Partindo de um ponto inicial x0 damos um passo α na direção d := −∇ f ( x ) de tal forma que α k := arg min 0 f ( x k + α d k ) e encontraremos um novo ponto x1 := x 0 + α d que melhorará o valor da função f ( x ) . Analisando o gráfico da função f ( x ) formado pela variação de α obtemos: α >
k
k
k
k
α >
A inclinação de f ( x ) em um determinado ponto relação: a (α ) := ∇ f ( x 0 + α .d )T d
α
pode ser obtido pela seguinte
A esta inclinação chamamos de derivada direcional de f ( x ) em relação à α . Assim sendo, o método de Falsa Posição pode ser resumido nas seguintes etapas: Dado x0 e um critério de parada ε , fazer:
1. d = −∇ f ( x) ; // direção de descida 2. α n = 0; // passo inicial igual a zero 3. dd n = ∇ f ( x 0 + α n d ) T d // derivada direcional inicial 4. α n 1 = 0.1 ; // dando um pequeno passo na direção d 5. dd n 1 = ∇ f ( x 0 + α n 1 d ) T d ; // encontrando valor da derivada direcional para α 1 6. Enquanto |dd n+1 | > ε faça: dd n 1 (α n 1 − α n ) ; // calculando valor da variação de α 7. ∆α = n n 1 +
+
+
+
+
dd
−
dd +
8. α 1 = α + ∆α ; // encontrando novo valor de α 9. dd n 1 = ∇ f ( x 0 + α n 1 d ) T d ; // derivada direcional para α n+1 10. Fim n+
n
+
+
Se f ( x ) for uma função quadrática então existe um derivada direcional dd = 0 de tal forma que a convergência do método ocorrerá em apenas duas iterações. Alguns resultados podem ser visualizados na próxima seção.
RESULTADOS Nesta seção, analisamos duas funções e obtemos vários resultados interessantes. Para a implementação dos algoritmos, bem como para a interpretação gráfica foi utilizada a ferramenta MatLab 6.0, em um PC com processador Atlon 2.0Ghz e 256Mb de memória RAM.
Função 1: f ( x1 , x2 ) = x12 + x2 2 O gráfico formado pela função 1 no intervalo [-1,1] é:
Como vemos, a função 1 é quadrática e portanto, independente do ponto x0 que partirmos o método convergirá em duas únicas iterações. Analisando a variação de f ( x ) em relação à α obtemos:
Graficamente, percebemos que para α = -1 f ( x ) =18, e que quando aumentamos o valor de α o valor de f ( x ) vai decrescendo até α =0.5 onde f ( x ) =0, e que a partir daí, conforme aumentamos α , o valor de f ( x ) também aumenta. Portanto, baseado nesta análise, deduzimos que α k := arg min 0 f ( x k + α d k ) = 0.5. Utilizando o método de falsa posição obtemos os seguintes resultados. α >
0
Teste
x
1 2 3 4 5 6
[1 1] [1 1] [1 1] [-1 -1] [-0.5 0.5] [1 -0.5]
n
x
10e-1 10e-3 10e-6 10e-3 10e-3 10e-3
0.500000 0.500000 0.500000 0.500000 0.500000 0.500000
[0 0] [0 0] [0 0] [0 0] [0 0] [0 0]
f ( x n ) Iterações
0 0 0 0 0 0
2 2 2 2 2 2
Como vemos, independente do ponto de partida x0 o método converge em duas iterações e, como prevíamos por meio da análise gráfica, o valor de k k α := arg min k 0 f ( x + α d ) = 0.5. α >
Função 2: f ( x1 , x2 ) = 100( x2 − x12 ) 2 + (1 − x1 ) 2 {Função de Rosenbrock) A Função 2 apresenta uma complexidade bem maior que a função 1 pois deixa de ser quadrática e isso faz com que a previsibilidade sobre ela seja menor. O gráfico formado pela função 2 no intervalo [-1 1] é:
Como vemos, a função 2 não é quadrática e, portanto, não podemos prever a quantidade de iterações bem como o valor de α k := arg min 0 f ( x k + α d k ) , no entanto, analisando o gráfico, percebemos que no ponto x=(1,1) o valor da função é igual a zero. Fica claro, calculando os pontos estacionários onde ∇ f ( x*) = 0 que o ponto x=(1,1) é mínimo da função, portanto, podemos prever que partindo de qualquer ponto x o método de falsa posição nos retornará um novo ponto xn+1 na direção do ponto x*=(1,1). Analisando a variação de f ( x ) em relação à α obtemos: α >
Como podemos ver, partindo de x=(0,0) para nenhum valor de α , f ( x ) = 0. O melhor valor de f ( x ) é aproximadamente igual à 0.777, quando α é aproximadamente igual à 0.090. Utilizando o método de falsa posição obtemos os seguintes resultados. Teste
1 2 3 4 5 6 7 8 9 10
0
n
x
[0 0] 10e-1 [0 0] 10e-3 [0 0] 10e-6 [-1 -1] 10e-3 [-0.5 0.5] 10e-3 [1 -0.5] 10e-3 [1 1] 10e-3 [0.99 0.99] 10e-3 [0.99 0.99] 10e-9 [0.99 0.99] 10e-14
x
0.072968 0.080627 0.080631 0.001564 0.002376 0.001613 0.010000 0.001013 0.001013 0.001013
f ( x n )
[0.145936 0.000000] 0.774783 [0.161254 0.000000] 0.771110 [0.161262 0.000000] 0.771110 [0.257184 -0.374535] 19.971568 [-0.611677 0.381195] 2.602467 [0.032216 -0.016108] 0.966004 [1.000000 1.000000] 0.000000 [0.993991 0.987994] 0.000036 [0.993991 0.987994] 0.000036 [0.993991 0.987994] 0.000036
Iterações
13 16 17 11 6 13 1 4 5 6
Como vemos, partindo do ponto x0 = (0,0) o valor de α tende à 0.900 tornando-se mais claro conforme aumentamos o critério de convergência, pois quando utilizamos ε =10e-1, α =0.072968 e f ( x ) =0.774783 e quando ε =10e-6, α =0.080631 e f ( x ) =0.771110 melhorando o valor da função. Analisando agora a função f ( x1 , x2 ) = 100( x2 − x12 ) 2 + (1 − x1 ) 2 notamos que quanto maior o valor de x2 em relação à x1 maior será o resultado e que quanto mais próximos estiverem x1 e x2 melhor será o valor da função. Fica evidente na tabela acima que o método busca sempre que possível reduzir o valor de x2 mesmo que para isso seja necessário aumentar x 1, porém sempre minimizando o valor de f ( x ) . Partindo de pontos próximos do mínimo da função o método melhora o valor da função lentamente pois o tamanho do passo calculado é muito pequeno e portanto f ( x ) reduzirá muito pouco. Isso explica a razão de que o método do gradiente quando aplicado a funções não-quadráticas com certas características quando chega próximo ao mínimo da função tem um rendimento baixo e uma convergência lenta, pois os passos são muito curtos.
ALGORITMO O algor itmo para Busca Linear do Método da Falsa Posição programado em MatLab 6.0 é: % Universidade Estadual de Campinas % Faculdade de Engenharia Eletrica e de Computacao % Nome: Tiago Agostinho de Almeida RA: 025625 % % ------------------------------------------------------------------------% METODO -> FALSA POSIÇÃO %-------------------------------------------------------------------------% LEGENDA: % ------------------------------------------------------------------------% BU() = Nome do Metodo = Busca Unidimensional % x1 = coordenada x do ponto x % x2 = coordenada y do ponto x % eps = epsilon, criterio de convergencia da derivada direcional % alfa(k) = tamanho do passo na direcao informada % x(k) = vetor-coluna x formada pelas coordenadas x1 e x2 informadas % G = gradiente da funcao no ponto x(k) % dd(k) = derivada direcional da funcao no ponto x(k) % d = direcao que minimiza a funcao % deltaalfa = alfa(k+1) - alfa(k) % F = valor da funcao no ponto x(k)
% ------------------------------------------------------------------------% Inicio function [x,y,z] = BU(); % Lim pando tela e inicializando as variaveis clear ; clc; % Fornecendo o ponto inicial e a funcao fprintf('\n'); fprintf('\n*** DADOS INICIAIS ***'); fprintf('\n'); fprintf('\nEscolha uma das funcoes abaixo:'); fprintf('\n[1] - f1 = x1^2 + x2^2 ......................'); fprintf('\n[2] - f2 = 100*(x2-x1^2)^2 + (1-x1)^2 .......'); op = input('\nOpcao numero: '); fprintf('\n'); x1 = input('Informe a coordenada inicial para x1: '); x2 = input('Informe a coordenada inicial para x2: '); e ps = input('Informe criterio de parada (epsilon): '); % Dados iniciais x0 = [x1; x2]; % Montando o vetor x com as coordenas informadas alfa0 = 0; % Tamanho do passo inicial alfa1 = 0.01; % Tamanho do passo seguinte % Calculando o valor do gradiente para a funcao escolhida switch op case {1} G = [2*x0(1); 2*x0(2)]; case {2} G = [-400*x0(1)*(x0(2)-(x0(1)^2)) - 2*(1-x0(1)); 200*(x0(2)-(x0(1)^2))]; end d = -G; %Direcao d % Derivada Direcional dd0 = G'*d; % Calculando novo x para alfa = 0.01 x1 = x0 + alfa1*d; % Grandiente para o novo x switch op case {1} G = [2*x1(1); 2*x1(2)]; case {2} G = [-400*x1(1)*(x1(2)-(x1(1)^2)) - 2*(1-x1(1)); 200*(x1(2)-(x1(1)^2))];
end % Derivada direcional para o alfa1 dd1 = G'*d; cont = 1; % Contador de iteracoes % Loop para encontrar o tamanho do passo (alfa1) while abs(dd1) > eps if cont >= 1000 % Solucionando o problema do loop infinito br eak end deltaalfa = (dd1*(alfa1-alfa0))/(dd0-dd1); % Calculando deltaalfa = alfa2 - alfa1 alfa0 = alfa1; % Atribuindo alfa1 em alfa0 alfa1 = alfa1 + deltaalfa; % Calculando novo tamanho do passo x1 = x0 + alfa1*d; % Calculando novo x % Calculando o Gradiente para o novo x switch op case {1} G = [2*x1(1); 2*x1(2)]; case {2} G = [-400*x1(1)*(x1(2)-(x1(1)^2)) - 2*(1-x1(1)); 200*(x1(2)-(x1(1)^2))]; end dd0 = dd1; % Atribuindo derivada direcional nova em derivada direcional anterior dd1 = G'*d; % Calculando nova derivada direcional % Incrementando contador de iteracoes cont = cont + 1; end % Valor da Funcao switch op case {1} F = x1(1)^2 + x1(2)^2; case {2} F = 100*(x1(2)-x1(1)^2)^2 + (1-x1(1))^2; end % Exibindo os resultados fprintf(1,' \n'); fprintf(1,'*** R ESULTADO FINAL ***'); fprintf(1,' \n'); fprintf(1,'%d iteracoes \n',cont); fprintf(1,'Tamanho do passo (alfa1) = %f \n',alfa1); fprintf(1,'Ponto final X(x1,x2) = X(%f,%f) \n',x1(1),x1(2)); fprintf(1,'Valor da Funcao f(x1,x2): %f \n', F); fprintf(1,' \n');
MÉTODO DA SEÇÃO ÁUREA
O Método da Seção Áurea utiliza um esquema de redução do intervalo de incerteza −1± 5 baseado na razão áurea. A Razão Áurea (r = ≅ 0.618 ) era considerada na 2 Antigüidade como a mais estética para os lados de um retângulo. Dados um in tervalo de incerteza I =[a,b], um critério de convergência ε e uma função do tipo f (α ) , o método da Seção Áurea realiza a redução do intervalo I a uma taxa de redução igual a razão áurea (r ) até que o critério de convergência seja satisfeito. Do b'+ a ' intervalo restante I’ =[a’,b’] o valor de α é obtido pela relação α = . 2 Quando utilizamos um critério de convergência pequeno, como 10 6 por exemplo, provavelmente o valor de α obtido pelo método de Seção Áurea será ig ual (ou próximo) do valor obtido pelo método de Falsa Posição. −
R ESULTADOS O Método da Seção Áurea foi testado para as duas funções utilizadas para o Método da Falsa Posição: 1 – Função 1 : f ( x1 , x2 ) = x12 + x2 2 2 – Função 2 : f ( x1 , x2 ) = 100( x2 − x12 ) 2 + (1 − x1 ) 2 {Função de Rosenbrock) Os resultados obtidos para a Função 1 utilizando o Intervalo I foram: 0
Teste
x
I
1 2 3 4 5 6
[1 1] [1 1] [1 1] [-1 -1] [-0.5 0.5] [1 -0.5]
[0 1] [0 1] [0 1] [0 1] [0 1] [0 1]
x
10e-1 10e-3 10e-6 10e-6 10e-6 10e-6
0.500000 0.497488 0.500000 0.500000 0.500000 0.500000
[0.000000 [0.005025 [0.000000 [0.000000 [0.000000 [0.000000
n
f ( x n )
0.000000] 0 0.005025] 0.000051 0.000000] 0 0.000000] 0 0.000000] 0 0.000000] 0
Iterações
1 11 25 25 25 25
No teste nº 1, utilizamos um critério de convergência relativamente grande e por causa disso o intervalo nem chegou a ser reduzido, e, como sabemos que o método utiliza a b'+ a' metade do intervalo como α ( α = ), o resultado por mero acaso deu que α =0.5, que 2 corresponde ao α k := arg min 0 f ( x k + α d k ) em apenas uma iteração, porém é válido ressaltar que este resultado só foi obtido por mero acaso, pois foi coincidência α estar localizado exatamente na metade do intervalo inicial. No teste nº 2 aumentamos o valor de ε e notamos que o valor de α chegou próximo do valor obtido pelo método da falsa posição. “Apertando” um pouco mais o critério de convergência, como no teste nº 3 o método da Seção Áurea obteve o mesmo valor do α obtido pela Falsa Posição. Como a análise gráfica foi realizada antes da implementação, podemos perceber que o método retornou o valor esperado e teve um comportamento adequado. Os testes seguintes provam que partindo de um mesmo ponto inicial e utilizando um critério de convergência adequado, o Método de Seção Áurea encontra um k k α := arg min k 0 f ( x + α d ) igual ao encontrado pelo Método de Falsa Posição. Fica claro que para funções quadráticas o Método da Falsa Posição converge mais rápido que o Método da Seção Áurea, pois o número de iterações utilizadas por esta é bem maior que as duas iterações necessárias para a convergência pelo cálculo da derivada direcional. α >
α >
Os resultados obtidos pela Função 2 utilizando o Intervalo I foram: 0
Teste
x
I
1 2 3 4 5 6 7 8 9 10
[0 0] [0 0] [0 0] [-1 -1] [-0.5 0.5] [1 -0.5] [1 1] [0.99 0.99] [0.99 0.99] [0.99 0.99]
[0 1] [0 1] [0 1] [0 1] [0 1] [0 1] [0 1] [0 1] [0 1] [0 1]
x
10e-1 10e-3 10e-6 10e-6 10e-6 10e-6 10e-6 10e-3 10e-6 10e-9
0.500000 0.081080 0.080628 0.001562 0.002377 0.001613 0.000005 0.004065 0.001012 0.001013
n
[1.000000 0.000000] [0.162159 0.000000] [0.161257 0.000000] [0.256214 -0.375018] [-0.611733 0.381135] [0.032226 -0.016113] [1.000000 1.000000] [1.006019 0.981951] [0.993988 0.987996] [0.993991 0.987994]
f ( x n )
Iterações
100.00000 0.771123 0.771110 19.971652 2.602469 0.966004 0.000000 0.090778 0.000036 0.000036
Como havíamos analisado graficamente, sabemos que partindo do ponto x = [0 0] k k α := arg min 0 f ( x + α d ) é aproximadamente igual à 0.09 e f ( x ) é aproximadamente k igual à 0.77. No teste nº 1 o resultado obtido foi bastante insatisfatório pois o valor de ε utilizado foi muito inadequado, isso fica claro analisando-se os dois testes subseqüentes em que reduzimos ε e conseguimos obter resultados muito próximos daqueles que prevíamos e que obtivemos pelo método da falsa posição. α >
1 11 25 25 25 25 25 11 25 40
Observamos no teste nº 7 que quando partimos do ponto x* = [1 1] o resultado obtido foi satisfatório, porém o número de iterações utilizadas (i = 25) foi muito além do número utilizado pelo método da Falsa Posição (i = 1). Isso mostra que quando estamos próximo do ponto ótimo o método da Seção Áurea tem um desempenho menor que o método da Falsa Posição se utilizarmos o mesmo critério de convergência. Isso fica claro nos 3 últimos testes, no qual partimos de um ponto próximo do ótimo e aumentamos o valor de ε gradativamente. Como podemos ver, os dois métodos apresentaram resultados satisfatórios e muito coincidentes com os quais estamos esperando.
ALGORITMO O algoritmo para Busca Linear do Método da Seção Áurea programado em MatLab 6.0 é: % Universidade Estadual de Campinas % Faculdade de Engenharia Eletrica e de Computacao % Nome: Tiago Agostinho de Almeida RA: 025625 % % ------------------------------------------------------------------------------% METODO -> BUSCA UNIDIMENSIONAL - SECAO AUREA %-------------------------------------------------------------------------------% LEGENDA: % ------------------------------------------------------------------------------% SAurea() = Nome do Metodo = Secao Aurea % r = razao aurea =~ 0,618 % x1, x2 = ponto inicial x (x1,x2) % I1 = ponto inicial do intervalo I % I2 = ponto final do intervalo I % eps = epsilon, criterio de convergencia da derivada direcional % alfa = razao aurea dada por x1 + (1-r)*(x2-x1) % beta = razao aurea dada por x1 + r*(x2-x1) % xna = valor do novo x para um deslocamento alfa % xnb = valor do novo x para um deslocamento beta % y1 = valor da funcao para xna % y2 = valor da funcao para xnb % F = valor da funcao final obtida pelo metodo % -------------------------------------------------------------------------
% Inicio function result = SAurea(); % Limpando tela e inicializando as variaveis clear; clc; % Fornecendo o ponto inicial e a funcao fprintf('\n'); fprintf('\n*** DADOS INICIAIS ***'); fprintf('\n'); fprintf('\nEscolha uma das funcoes abaixo:'); fprintf('\n[1] - f1 = x1^2 + x2^2 ......................'); fprintf('\n[2] - f2 = 100*(x2-x1^2)^2 + (1-x1)^2 .......'); op = input('\nOpcao numero: '); fprintf('\n'); x1 = input('Informe o ponto x1: '); x2 = input('Informe o ponto x2: '); I1 = input('Informe o ponto inicial do intervalo (a): '); I2 = input('Informe o ponto final do intervalo (b): '); eps = input('Informe criterio de parada (epsilon): '); % Dados iniciais x = [x1;x2]; r = (sqrt(5)-1)/2; I = [I1; I2]; % Montando o vetor x com as coordenas informadas alfa = I1 + (1-r)*(I2-I1); beta = I1 + r*(I2-I1); % Calculando o valor da funcao y1 para alfa e y2 para beta switch op case {1} G1 = [2*x(1);2*x(2)]; %Valor do gradiente no ponto x d1 = -G1; %Direcao de descida xna = x + alfa*d1; %Valor do novo x para deslocamento alfa xnb = x + beta*d1; %Valor do novo x para deslocamento beta y1 = xna(1)^2 + xna(2)^2; %Valor da funcao para x com deslocamento alfa y2 = xnb(1)^2 + xnb(2)^2; %Valor da funcao para x com deslocamento beta case {2} G2 = [-400*x(1)*(x(2)-(x(1)^2)) - 2*(1-x(1)); 200*(x(2)-(x(1)^2))]; %Valor do gradiente no ponto x d2 = -G2; %Direcao de descida xna = x + alfa*d2; %Valor do novo x para deslocamento alfa xnb = x + beta*d2; %Valor do novo x para deslocamento beta y1 = 100*(xna(2)-xna(1)^2)^2 + (1-xna(1))^2; y2 = 100*(xnb(2)-xnb(1)^2)^2 + (1-xnb(1))^2; end
cont = 1; % Contador de iteracoes % Loop para encontrar o tamanho do passo (alfa) while (I2-I1) > eps if cont >= 1000 % Solucionando o problema do loop infinito break end if y1>y2 I1 = alfa; alfa = beta; y1 = y2; beta = I1 + r*(I2-I1); switch op case {1} xnb = x + beta*d1; y2 = xnb(1)^2 + xnb(2)^2; case {2} xnb = x + beta*d2; y2 = 100*(xnb(2)-xnb(1)^2)^2 + (1-xnb(1))^2; end else I2 = beta; beta = alfa; y2 = y1; alfa = I1 + (1-r)*(I2-I1); switch op case {1} xna = x + alfa*d1; y1 = xna(1)^2 + xna(2)^2; case {2} xna = x + alfa*d2; y1 = 100*(xna(2)-xna(1)^2)^2 + (1-xna(1))^2; end end % Incrementando contador de iteracoes cont = cont + 1; end alfa = (I2+I1)/2; % Valor da Funcao switch op case {1} xna = x + alfa*d1; F = xna(1)^2 + xna(2)^2;
case {2} xna = x + alfa*d2; F = 100*(xna(2)-xna(1)^2)^2 + (1-xna(1))^2; end % Exibindo os resultados fprintf(1,' \n'); fprintf(1,'*** RESULTADO FINAL ***'); fprintf(1,' \n'); fprintf(1,'%d iteracoes \n',cont); fprintf(1,'Tamanho do passo (alfa1) = %f \n',alfa); fprintf(1,'Ponto Final X(x1,x2) = X(%f,%f) \n',xna(1),xna(2)); fprintf(1,'Intervalo Final I(I1,I2) = I(%f,%f) \n',I1,I2); fprintf(1,'Valor da Funcao f(x1,x2): %f \n', F); fprintf(1,' \n');
