Ligação de 4 cidades de comprimento mínimo: análise de uma família de soluções

Prof. Adriano Azevedo Filho (azevedofilho@usp.br)

Conteúdo

 1 - Descrição do problema
 
 2 - Análise de uma família de soluções
 
 3 - Solução por força bruta (Simulação Monte Carlo)
 
 4 - Solução por cálculo resolvendo condição de primeira ordem
 
 5 - Solução por Newton-Raphson (ponto fixo e mapa de contração) 
 
 6 - Solução pelo uso direto do Mapa de Contração
 
 7 - Conclusão
 
 8 - Referência recomendada

1 - Descrição do problema

Suponha que deseja construir estradas para ligação de 4 cidades situadas nos vértices de um quadrado, ao menor custo possível. Suponha que a distância entre os vértices do quadrado é 1 (a unidade não importa) e o custo para construção é sempre o mesmo sobre a área ao redor das cidades. O custo mínimo será dado pela ligação que exigir o menor comprimento total de estradas, de forma que a partir de cada uma das 4 cidade seja possível chegar a qualquer uma das outras três cidades.

A figura a seguir ilustra essa situação, mostrando 3 possíveis maneiras de ligar as cidades (há infinitas maneiras).

VA

Para as alternativas \(A\), \(B\) e \(C\) podemos calcular o comprimento com facilidade. Os comprimentos das estradas nas alternativas \(A\) e \(B\) exige somente a soma dos comprimentos de arestas do quadrado. A avaliação da distância na alternativa \(C\) exige o cálculo do comprimento da estrada diagonal, algo que envolve algum conhecimento de trigonometria, em particular o teorema de Pitágoras:

Teorema de Pitágoras: Para um triângulo retângulo qualquer com \(h\) representando o comprimento da hipotenusa e \(c_1, c_2\) o comprimento dos catetos, é verdade que \(h^2=c_1^2 + c_2^2.\)

Como o comprimento dos lados do quadrado é 1, o valor da diagonal, correspondente à hipotenusa será dado por * \(\sqrt{1^2+1^2}=\sqrt{2}\)

Se representarmos por \(d(z)\) o comprimento total de uma alternativa \(z\) qualquer temos, para as alternativas \(A\), \(B\) e \(C\):

Dessas 3, a alternativa \(C\) parece a melhor pois temos

Essa notação diz que dentre todos os possíveis “valores” de \(z\) dentro do domínio \(\{A,B,C\}\), que inclui essas 3 alternativas analisadas, \(C\) leva \(D(\centerdot)\) ao valor mínimo (arg min representa o argumento que leva a função ao mínimo global dentro das possibilidades analisadas).

Do ponto de vista técnico a função \(d\), nesse caso restrito, é um mapa dos pontos no domínio (infinitas alternativas possíveis) a valores no contra-domínio (comprimentos) genericamente representado pelos \(\mathbb{R}^+\), ou seja os numeros reais positivos, ou seja,

A solução obtida é muito restritiva, dado que só considerou 3 alternativas de um conjunto “infinito” de possibilidades, que é difícil de ser especificado. Uma solução melhor é discutida a seguir.

2 - Análise de uma família de soluções

Um procedimento formal que possa levar ao “desenho” das estradas que leve ao comprimento mínimo não é nada óbvio e pode exigir técnicas matemáticas avançadas. A própria caracterização algébrica das possíveis alternativas é algo complicado.

Podemos, contudo, pesquisar uma “família de alternativas”, representada por \(F_x\), que inclui \(B\) e \(C\) como casos particulares, ilustrada na figura a seguir:

VA

Os (infinitos) membros dessa família \(F_x\) de alternativas podem ser particularizados de acordo com o valor de \(x\) definido no intervalo fechado \([0, 1/2]\). A alternativa \(B\) corresponde ao caso particular em que \(x=0\) e a alternativa \(C\) corresponde ao caso em que \(x=1/2\).

Para um dado valor de \(x\), o comprimento total será dado por \(b+4a\). Os valores de \(a\) e \(b\) podem ser calculados a partir do valor de \(x\), da caracterização geométrica de \(F_x\) e aplicação do Teorema de Pitágoras:

para \(x\in [0, 1/2]\).

Se usarmos \(c(x)\) para representar o comprimento das estradas na alternativa \(F_x\) temos que

Essa função \(c:[0,1/2]\rightarrow \mathbb{R}\) é claramente contínua no domínio \([0,1/2]\subseteq \mathbb{R}\), que é um conjunto compacto, algo que possibilita estabelecer, pelo teorema de Weiestrass (máximo e mínimo) que essa função terá um ponto de mínimo global e um ponto de máximo global.

3 - Solução por força bruta (simulação Monte Carlo)

A estratégia mais geral possível para encontrar pontos próximos do ponto de mínimo global, poderia considerar a verificação de muitos valores \(x_i\in [0,1/2]\), \(i=1,\ldots,n\), extraídos aleatóriamente do intervalo, e verificação de \(c(x_i)\) para cada valor. Se

Ache uma estimativa do ponto de mínimo por esse processo. Funções uteis do R para essa empreitada: runif, which.min

Definindo a função \(c(x)\) no R (vamos chamá-la f_c):

f_c<-function(x){
  return(1-2*x+4*sqrt(0.25+x^2))
}

Implementando o procedimento, com \(n=100000\)

options(digits=15)
set.seed(100)
x<-runif(100000,0,0.5)
y<-f_c(x)
x[which.min(y)]
## [1] 0.288679967983626
x<-runif(1000000,0,0.5)
y<-f_c(x)
x[which.min(y)]
## [1] 0.288675682037137
y=rnorm(length(x),0,0.01)
plot(x,y,cex=0.2)
abline(h=0,col="red")

o gráfico acima mostra o grau de cobertura do domínio definido pelo conjunto \([0,1/2]\), considerado pelo procedimento de simulação. Os valores foram plotados com um certo deslocamento aleatório vertical, para facilitar a visualização da cobertura.

4 - Solução por cálculo, usando a continuidade e diferenciabilidade de \(c(x)\)

Como a função \(c: S\rightarrow S, c(x)\) é diferenciável (e portanto contínua) e \(S=[0,1/2]\subset \mathbb{R}\) é um conjunto compacto, sabemos pelo Teorema de Weirstrass (máximos e mínimos globais), que existirão pontos de máximo e mínimos globais nesse intervalo. Para encontrar esses pontos podemos usar o seguinte resultado:

Teorema: se \(x_*\) representa um ponto de máximo global ou ponto de mínimo global, de uma função diferenciável num intervalo fechado, \(x_*\) será: (a) um ponto na fronteira de \(S\), ou, (b) será um ponto interior de \(S\). Se \(x_*\) for um ponto interior, pela continuidade de \(c(x)\), deverá atender a \(c'(x)=0\) (veja Piskounov 1969, cap. 5 sec. 6 para detalhes).

Logo, para encontrar o máximo e mínimo global de \(c(x)\), devemos: (a) encontrar \(c(x)\) para \(x\in \text{fr}(S)\), com \(\text{fr}(S)=\{0,1/2\}\); (b) encontrar o(s) valor(es) de \(x\) que atende(m) a \(c'(x)\), para \(x\in (0,1/2)\) e obter o valor de \(c(x)\) nesses pontos (se existirem); (c) comparar os resultados de (a) e (b), algo que possibilitará o reconhecimento dos pontos de máximo global e mínimo global

Operacionalização dos procedimentos

Encontrando os valores de \(c(x)\) para \(x\in \text{fr}(S)\)

Como \(\text{fr}(S)=\{0,1/2\}\) temos os valores de \(c(0)\) e \(c(1/2)\) definidos respectivamente, por: (usando a função f_c no R, definida anteriormente)

## c(0)
f_c(0)
## [1] 3
## c(1/2)
f_c(1/2)
## [1] 2.82842712474619

Encontrando pontos \(x\) que atendem a \(c'(x)=0\) no interior de \([0,1/2]\)

Como já vimos, temos

  • \(c'(x)=\displaystyle \frac{d\, c(x)}{d\,x}\;\;=-\displaystyle 2+\frac{4\,x}{\sqrt{0{,}25 + x^2}}\)

Encontrar um ou mais valores de \(x\) que atendam a \(c'(x)=0\) pode ser um problema difícil em muitos casos. Neste caso é até fácil. Representando por \(x^*\) o valor de \(x\) que atende a última equação (derivada igual a zero), chegamos a

  • \(x^*=\dfrac{\sqrt{3}}{6}=0,288675134594813\)

e podemos encontrar \(c(x^*)\) para compararmos com os valores obtidos no teste com os pontos na fronteira, usando a função f_c no R, definida anteriormente)

## c(raiz(3)/6)
f_c(sqrt(3)/6)
## [1] 2.73205080756888

comparando esse resultado com as distâncias obtidas nos pontos de fronteira, concluimos que o mínimo é obtido no interior do domínio com

  • \(x^*=\dfrac{\sqrt{3}}{6}=0,288675134594813\), que leva a \(c(x^*)=2,73205080756888\) (mínima distância)

5 - Solução pelo método de Newton-Raphson

Em muitas situações, não é possível resolver algebricamente a condição \(c'(x)=0\), e nesses casos outros métodos são necessários. Vamos apresentar a seguir um processo de solução que é motivado pelo método de Newton-Raphson, aplicado à solução numérica de raízes de equações. Esse método é fundamentado na obtenção de uma sequência que, caso sejam atendidas condições de regularidade, converge para um ponto fixo, usando noções relacionadas a noção de Mapa de Contração. E esse ponto fixo é exatamente a raiz desejada da equação.

A sequência é estabelecida pela processo iterativo definido por

Em condições ideais (atendimento ao requisito para ser um mapa de contração), partindo-se de um ponto inicial \(x_0\), essa \(x_1,x_2,\ldots,x_n\) converge para uma raiz da \(f(x)\).

No caso,

inicialmente definindo funções no R que representem \(f(x)\) e sua derivada \(f'(x)\):

## definindo f(x)
m_f<-function(x){
  -2+4*x/sqrt(1/4+x^2)
}
## definido f'(x)
m_fderiv<-function(x,delta=0.0001){
  (m_f(x+delta)-m_f(x))/delta
}

Com essas funções no R podemos estabelecer o processo recursivo, através de

newton<-function(epsilon,x_0){
  x_velho<-x_0
  n<-0
  while (abs(m_f(x_velho))>=epsilon){
    x_novo<-x_velho-m_f(x_velho)/m_fderiv(x_velho)
    x_velho<-x_novo
    n<-n+1
    if (n>100) break
  }
  cat("x*=",x_velho,"f(x)",m_f(x_velho),"n=",n,"\n")
}

Vamos agora testar o processo, partindo de um ponto qualquer dentro do intervalo \([0,1/2]\), por exemplo, \(x_0=1/4\). A condição de término envolve chegar a uma solução em que \(|f(x)|<0.0000000001\) um valor suficientemente próximo de zero para fins práticos (pode ser alterado em caso de interesse).

newton(0.0000000001,1/4)
## x*= 0.28867513459475 f(x) -3.28403970684121e-13 n= 4

6 - Solução pelo uso direto do Mapa de Contração em Função Lipschitz Contínua

Uma outra possibilidade para solução do problema é aplicar diretamente a condição estabelecida pelo Teorema de Banach, descrito abaixo.

Trabalhando a condição \(c'(x)=0\), no tópico 4, chegamos à equação:

Claramente, solução desejada \(x_*\), se existir, deverá atender a \(f(x_*)=x_*\), ou seja, ser um ponto fixo de \(f(x_*)\).

Examinando essa função \(f(x)\), podemos concluir que ela é Lipschitz contínua no intervalo \([0,1/2]\) dado que

Note que imagem do domínio \([0,1/2]\) de \(f(x)\) é \([0,25;0,3536]\subseteq [0,1/2]\), logo, podemos afirmar que a função \(f\) é \(f:S\rightarrow S\), com \(S=[0,1/2]\). Considere agora o:

  • Teorema do ponto fixo de Banach: No contexto do espaço métrico \((X,d)\), uma função contínua \(g: S\rightarrow S\), em que \(S\subseteq X\) é um conjunto compacto, é um mapa de contração, é verdade que: (a) \(g(x)\) tem um único ponto fixo \(x_*\); e (b) o ponto fixo da função \(g\) pode ser encontrado pelo limite da sequência \(\{x_n\}\), com \(x_{i+1}=g(x_i)\), dado que \(x_n\rightarrow x_*\), em que \(x_0\in S\) é um ponto qualquer de \(S\). O teorema garante também que
  • \(d(x_n,x_*)\le \frac{K^n}{1-K} d(x_1,x_0)\), em que \(d\) é a métrica usada (nesse caso usaremos a métrica euclidiana, \(d(x,y)=|x-y|\)). Nas expressão, \(K\) é a melhor constante de Lipschitz associada à função \(g(x)\).

Os resultados obtidos quanto à função \(f(x)\) (ser mapa de contração e ser \(S\rightarrow S\) com \(S\) compacto) garantem a aplicabilidade do teorema de Banach para encontrar o ponto fixo \(x_*\) em que \(f(x_*)=x_*\). Pela definição de \(f(x)\) essa solução \(x_*\) levará a \(c'(x_*)=0\) e é única.

Implementação computacional da obtenção do ponto fixo \(x_*\)

Definição da função \(f(x)\) no R:

f_f <-function(x){
  return(1/4*sqrt(1+4*x^2))
}

Precisamos agora estabelecer a sequência \(\{x_n\}\) nas condições especificadas pelo Teorema de Banach, algo que é fácil do ponto de vista computacional.

f_recursao <-function(x0,n,K){
  x1<-f_f(x0)
  numdigits<-getOption("digits") ## salva número de digitos nos resultados
  options(digits=15)    ## define 15 digitos nos resultado
  cat("i=",0,"x_i=",x0,"\n")
  cat("i=",1,"x_i=",x1,"\n")
  erro_maximo=K^n/(1-K)*abs(x1-x0)
  x<-x1
  for(i in 2:(n-1)) {
    x<-f_f(x)
    cat("i=",i,"x_i=",x,"\n")
  } 
  xn<-f_f(x)
  cat("i=",n,"x_i=",xn,"\n")
  cat("d(x_n,x_(n-1)=",abs(xn-x),"\n")
  cat("d(x_n,x*) máxima:", erro_maximo,"\n") 
  options(digits=numdigits) ## retorna o valor inicial de digitos nos resultados
  return(xn)
}

Obtendo a sequência \(\{x_n\}\) considerando \(n=30\), observe a velocidade de convergência e o erro máximo que está sendo cometido ao se utilizar \(x_n\) como a estimativa de \(x_*\).

xn<-f_recursao(0.45,30,K=0.3536)
## i= 0 x_i= 0.45 
## i= 1 x_i= 0.336340601176843 
## i= 2 x_i= 0.301299269829849 
## i= 3 x_i= 0.29188236072089 
## i= 4 x_i= 0.289480272428019 
## i= 5 x_i= 0.288876629430714 
## i= 6 x_i= 0.288725521486779 
## i= 7 x_i= 0.288687732142281 
## i= 8 x_i= 0.288678284033218 
## i= 9 x_i= 0.288675921957635 
## i= 10 x_i= 0.28867533143572 
## i= 11 x_i= 0.288675183805052 
## i= 12 x_i= 0.288675146897374 
## i= 13 x_i= 0.288675137670453 
## i= 14 x_i= 0.288675135363723 
## i= 15 x_i= 0.28867513478704 
## i= 16 x_i= 0.28867513464287 
## i= 17 x_i= 0.288675134606827 
## i= 18 x_i= 0.288675134597816 
## i= 19 x_i= 0.288675134595564 
## i= 20 x_i= 0.288675134595001 
## i= 21 x_i= 0.28867513459486 
## i= 22 x_i= 0.288675134594825 
## i= 23 x_i= 0.288675134594816 
## i= 24 x_i= 0.288675134594814 
## i= 25 x_i= 0.288675134594813 
## i= 26 x_i= 0.288675134594813 
## i= 27 x_i= 0.288675134594813 
## i= 28 x_i= 0.288675134594813 
## i= 29 x_i= 0.288675134594813 
## i= 30 x_i= 0.288675134594813 
## d(x_n,x_(n-1)= 0 
## d(x_n,x*) máxima: 5.01731880155248e-15

6 - Conclusão

Os procedimentos apresentados na últimas seções mostram 4 maneiras distintas de resolver o problema original de otimização dentro da familia de soluções especificada.

Uma questão não resolvida aqui: seria essa a solução ótima geral, considerando todos os possíveis desenhos de caminhos entre as 4 cidades? o que resolvemos foi a solução ótima para a família de soluções proposta. Felizmente essa solução é a ótima dentre todos os possíveis desenhos, mas os argumentos necessários para demonstrar isso estão fora do escopo desta nota técnica.

8 - Referência Útil