Capítulo 3

Processos de Otimização Numérica

“La vie n’est bonne qu’à deux choses: à faire des mathématiques et à les professer.”

– François Arago in Notices Biographiques, 1854

3.1. Introdução

O ponto de partida para a análise de um sistema consiste na construção de um modelo que permita realizar previsões sobre seu comportamento. Um modelo pode ser compreendido como uma abstração de um sistema real, concebida com o propósito de facilitar a obtenção de previsões e a formulação de estratégias de controle. Tais modelos são, com frequência, de natureza matemática — representados por equações, funções ou relações formais — ou gráfica. Neste sentido, para que um modelo seja efetivamente útil, é necessário que ele concilie duas características frequentemente conflitantes: realismo e simplicidade. Modelos excessivamente detalhados podem se tornar intransponíveis do ponto de vista analítico ou computacional, enquanto modelos demasiadamente simplificados podem deixar de capturar comportamentos relevantes do sistema.

Uma vez construído e validado o modelo, a obtenção de soluções pode ser realizada por métodos analíticos ou numéricos. Soluções analíticas decorrem diretamente da estrutura matemática do modelo e são, muitas vezes, desejáveis por sua generalidade e exatidão. Todavia, em muitos casos práticos, tais soluções não estão disponíveis, o que impõe a necessidade do emprego de métodos numéricos, capazes de fornecer aproximações computacionalmente viáveis. Este capítulo dedica-se ao estudo dos conceitos e técnicas de otimização numérica, abordando tanto os fundamentos teóricos quanto aspectos práticos de implementação, com ênfase em problemas estatísticos que envolvem, por exemplo, técnicas de máxima verossimilhança e procedimentos baseados em simulação estocástica.

3.2. Otimização Numérica

A otimização constitui, em termos gerais, um processo fundamentado em modelos matemáticos, cujo objetivo é determinar a solução mais eficiente para um dado problema, conforme critérios previamente estabelecidos. A etapa inicial desse processo envolve a definição rigorosa da função objetivo, a qual representa uma medida quantitativa de desempenho — como lucro, tempo, custo, entre outras — a ser maximizada ou minimizada. De acordo com a natureza dessa função e das restrições envolvidas, os problemas de otimização podem ser classificados, de forma abrangente, em três categorias principais:

• Programação Linear (PL): tanto a função objetivo quanto as restrições são expressas por relações lineares.
• Programação Quadrática (PQ): a função objetivo é uma forma quadrática, enquanto as restrições permanecem lineares.
• Programação Não-Linear (PNL): ao menos uma das componentes — seja a função objetivo ou alguma restrição — apresenta comportamento não-linear.

A resolução de problemas de otimização é realizada por meio de um otimizador, isto é, um algoritmo ou procedimento computacional desenvolvido para identificar, de forma eficiente, o ponto ótimo da função objetivo, respeitando as restrições impostas. Em ambientes computacionais, a interface de chamada de um otimizador é geralmente definida por um conjunto básico de argumentos, conforme ilustrado no Código 3.1.

Código 3.1. Estrutura computacional geral de otimizador.

# ---------------------------------
# Estrutura Geral de um Otimizador
# ---------------------------------

optimizer(objective, 
          start,
          constraints, 
          bounds = NULL, 
          method,
          tol = 1E-6,
          maximum = FALSE)

em que:

• objective: função objetivo a ser minimizada ou maximizada.
• start: vetor de valores iniciais a partir do qual se inicia o processo iterativo de otimização. A escolha de bons valores iniciais pode influenciar significativamente a convergência do algoritmo.
• constraints: restrições lineares ou não lineares que a solução deve satisfazer. Quando NULL, assume-se ausência de restrições explícitas.
• bounds: limites inferiores e superiores admissíveis para as variáveis, geralmente definidos como intervalos fechados ou semiabertos.
• method: método numérico a ser empregado na otimização. A escolha depende do tipo de problema (existência ou não de restrições, convexidade, diferenciabilidade da função objetivo, entre outros aspectos).
• tol: tolerância empregada no critério de parada, como, por exemplo, a diferença entre valores sucessivos da função objetivo ou entre estimativas dos parâmetros.
• maximum:valor lógico que indica se a função objetivo deve ser maximizada (TRUE) ou minimizada (FALSE, valor padrão).

Naturalmente, a escolha adequada do otimizador e a formulação rigorosa do problema são elementos essenciais para garantir soluções robustas. Isso implica que, dada a heterogeneidade e complexidade intrínsecas desses problemas, a eficiência dos métodos de otimização depende do equilíbrio entre a fidelidade do modelo, as limitações computacionais e as características dos algoritmos empregados. É importante ressaltar que, antes de se trabalhar com esses otimizadores, é indispensável um conhecimento aprofundado dos métodos numéricos que os fundamentam. Por essa razão, alguns destes métodos serão discutidos nas próximas seções.

3.3. Métodos Numéricos

Os métodos numéricos são técnicas matemáticas empregadas na resolução de problemas de otimização, especialmente quando soluções analíticas são inviáveis ou, até mesmo, inexistentes. Esses métodos empregam algoritmos e procedimentos computacionais para obter aproximações da solução ótima, controlando rigorosamente os níveis de precisão. Nesta seção, serão explorados quatro métodos básicos: o método da bissecção (Seção 3.3.1), o método de Newton-Raphson (Seção 3.3.2), o recozimento simulado (em inglês, simulated annealing - Seção 3.3.3) e o algoritmo genético (Seção 3.3.4). Cada um desses métodos apresenta particularidades quanto à convergência, às hipóteses exigidas sobre a função objetivo e às estratégias adotadas para a exploração do espaço de busca. A eficácia das soluções obtidas está diretamente relacionada à calibração adequada dos parâmetros iterativos e à definição precisa dos critérios de parada.

3.3.1. Método da Bissecção

---

O problema de localizar raízes de funções ou pontos críticos de funções deriváveis é recorrente em diversas áreas da matemática aplicada, estatística, engenharia e ciências computacionais. Quando não há solução analítica disponível, torna-se imprescindível recorrer a métodos numéricos que assegurem a aproximação do valor procurado com precisão arbitrária. Nesse contexto, o método da bissecção destaca-se pela simplicidade conceitual, robustez e garantia de convergência, desde que a função considerada seja contínua e que exista mudança de sinal em um intervalo inicial conhecido. Embora apresente taxa de convergência relativamente lenta em comparação a outros métodos, como Newton-Raphson ou métodos quasi-Newtonianos, sua estabilidade numérica e segurança o tornam uma ferramenta essencial, sobretudo em estágios iniciais da análise ou quando há incerteza sobre a regularidade da função derivada.

Definição 3.1 (Método da Bissecção - Ruggiero e Rocha Lopes 1998; Srivastava e Dixit 2012; Ehiwario e Aghamie 2014; Solanki, Thapliyal, e Tomar 2014; Milne 2015). Seja \(f'\) uma função contínua no intervalo fechado \([a_0, b_0]\), tal que \(f'(a_0), f'(b_0) \leqslant 0\). Pelo Teorema do Valor Intermediário (Guidorizzi 2002; Andreescu et al. 2017), existe ao menos um ponto \(x_k \in [a_0, b_0]\) tal que \(f'(x_k) = 0\). Nessas condições, o ponto \(x_k\) constitui um candidato a ponto de máximo (ou mínimo) da função \(f\). O método da bissecção baseia-se na redução iterativa do intervalo de busca \([a_0, b_0]\), mediante a construção de uma sequência de subintervalos \([a_t, b_t]\), para \(t = 0, 1, 2, \ldots\), de forma que:

\[\begin{align}\\ [a_0, b_0] \supset [a_1, b_1] \supset \cdots \supset [a_t, b_t]\\\\ \end{align}\]

até que um critério de parada previamente estabelecido seja satisfeito. Para iniciar o processo, define-se o ponto médio do intervalo como:

\[\begin{align}\\ x^{(0)} = \dfrac{a_0 + b_0}{2}\\\\ \end{align}\]

A cada iteração, o intervalo é atualizado com base no sinal do produto \(f'(a_t)f'(x^{(t)})\), conforme a regra:

\[\begin{align}\\ [a_{t+1}, b_{t+1}] = \begin{cases} [a_t, x^{(t)}], & \text{se } f'(a_t)f'(x^{(t)}) \leqslant 0, \\\\ [x^{(t)}, b_t], & \text{se } f'(a_t)f'(x^{(t)}) > 0, \end{cases}\\\\ \end{align}\]

e o novo ponto médio é então calculado por:

\[\begin{align}\\ x^{(t+1)} = \dfrac{1}{2}(a_{t+1} + b_{t+1})\\\\ \end{align}\]

Este procedimento é repetido até que a largura do intervalo \([a_t, b_t]\) se torne suficientemente pequena. Para tanto, é fundamental estabelecer uma regra de parada que determine o término do processo iterativo. Em cada iteração, verifica-se se as condições de convergência foram atingidas. Quando isso ocorre, o ponto \(x^{(t+1)}\) é adotado como aproximação final da solução. Os critérios de parada mais comuns baseiam-se nas seguintes condições:

• Critério de Convergência Absoluta: \[\begin{align}\\ \left|x^{(t+1)} - x^{(t)}\right| < \epsilon\\\\ \end{align}\]

• Critério de Convergência Relativo: \[\begin{align}\\ \dfrac{\left|x^{(t+1)} - x^{(t)}\right|}{\left|x^{(t)}\right|} < \epsilon\\\\ \end{align}\]

em que \(\epsilon > 0\) denota uma constante de tolerância previamente estipulada. Cabe ressaltar que, na prática, limitações decorrentes da aritmética de ponto flutuante podem comprometer a estabilidade e a precisão do método. Para mitigar tais efeitos numéricos, recomenda-se, sempre que possível, reescrever o cálculo do ponto médio em uma forma numericamente mais estável, tal como:

\[\begin{align}\\ a_{t+1} = a_t + \frac{b_t - a_t}{2}\\\\ \end{align}\]

em vez da forma simétrica \(a_{t+1} = (a_t + b_t)/2\), com o intuito de evitar cancelamentos numéricos ou perda significativa de precisão durante o processo iterativo.

Exemplo 3.1. Considere a função:

\[\begin{align}\\ f(x) = \dfrac{\log(x)}{1 + x}\\\\ \end{align}\]

Suponha que nosso objetivo seja determinar o ponto de máximo em relação à variável \(x\). Note que, esta função, em particular, é uma função cuja derivada não admite solução analítica explícita para a equação \(f'(x) = 0\), impossibilitando, portanto, a determinação exata do ponto de máximo por métodos analíticos. Nesse cenário, torna-se necessário recorrer a procedimentos numéricos de otimização. Contudo, antes da aplicação de qualquer algoritmo iterativo, recomenda-se uma análise gráfica preliminar (Figura 3.1) da função \(f(x)\), com o intuito de compreender seu domínio, identificar possíveis descontinuidades ou assimetrias, e localizar regiões de interesse para a busca do ótimo.

Figura 3.1.. Ilustração do comportamento da função \(f(x)\) no intervalo \([1,5]\).

Nesta figura, observa-se que o valor máximo de \(f(x)\) ocorre em um ponto próximo de \(x = 3.5\). Com base nessa evidência empírica, é razoável escolher \(x^{(0)} = 3.0\) como ponto inicial para um processo iterativo de maximização. O objetivo é ajustar progressivamente este valor até atingir a raiz de \(f'(x)\), isto é, o ponto em que a derivada da função se anula. A derivada de \(f(x)\) é dada por:

\[\begin{align}\\ f'(x) = \dfrac{1 + \dfrac{1}{x} - \log(x)}{(1 + x)^2}\\\\ \end{align}\]

Para a aplicação do método da bissecção, estabelece-se inicialmente o intervalo de busca \([a_0, b_0] = [1, 5]\). É importante garantir que \(f'(a_0) \cdot f'(b_0) < 0\), pois tal condição assegura a existência de, pelo menos, uma raiz de \(f'\) no interior do intervalo, em virtude do Teorema do Valor Intermediário. O método da bissecção, por sua vez, procede iterativamente, dividindo o intervalo corrente ao meio em cada iteração, sendo que, em cada passo \(k\), calcula-se o ponto médio do intervalo \([a_k, b_k]\). Em seguida, avalia-se o sinal do produto \(f'(a_k) \cdot f'(m_k)\). Se o produto for negativo, então a raiz encontra-se em \([a_k, m_k]\), e define-se \(b_{k+1} = m_k\) e \(a_{k+1} = a_k\). Caso contrário, a raiz está em \([m_k, b_k]\), e define-se \(a_{k+1} = m_k\) e \(b_{k+1} = b_k\). Este processo se repete até que o comprimento do intervalo \(|b_k - a_k|\) seja inferior a uma tolerância pré-fixada ou até que se atinja um número máximo de iterações.

O Código 3.2 apresenta uma implementação do método da bissecção, na qual se fixa o número de iterações em itr = 40. Ressalta-se, contudo, que tal implementação não está otimizada, pois não contempla critérios adicionais de parada baseados, por exemplo, na magnitude de \(|f'(m_k)|\), que poderiam acelerar a detecção de convergência e evitar cálculos desnecessários. Por fim, a Figura 3.2 ilustra o gráfico da função objetivo juntamente com o ponto de máximo estimado por meio do método da bissecção.

Código 3.2. Implementação do método da bissecção para a maximização da função objetivo, \(f(x)\).

# -------------------------------------
# Implementação do Método da Bissecção
# -------------------------------------


# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Valores iniciais ---

a    <- 1              # Extremo esquerdo do intervalo inicial
b    <- 5              # Extremo direito do intervalo inicial
x    <- a + (b - a)/2  # Ponto médio inicial
itr  <- 40             # Número de iterações

# --- 2. Definição da função objetivo ---

f    <- function(x)
{
  log(x) / (1 + x)
}

# --- 3. Derivada da função objetivo ---

df   <- function(x) 
{
  (1 + 1/x - log(x)) / (1 + x)^2
}

# --- 4. Implementação do método da bissecção ---

for (i in 1:itr) 
{
  if (df(a) * df(x) < 0) {
    b <- x
  } else {
    a <- x
  }
  x <- a + (b - a)/2
}

# --- 5. Resultados ---

list('Ponto Ótimo' = x,
     'Função Objetivo no Ponto Ótimo' = f(x),
     'Derivada no Ponto Ótimo' = round(df(x), 0)
     )

$`Ponto Ótimo`
[1] 3.591121

$`Função Objetivo no Ponto Ótimo`
[1] 0.2784645

$`Derivada no Ponto Ótimo`
[1] 0

# --- 6. Visualização gráfica ---

x_vals            <- seq(1, 5, by = 0.0001)
df_plot           <- data.frame(x = x_vals,
                                fx = f(x_vals))

x_k               <- x
df_plot$highlight <- ifelse(abs(df_plot$x - x_k) < 1e-4, "Ótimo", "Não")

ggplot(df_plot, aes(x = x, y = fx)) +
  geom_line(size = 0.5, color = "steelblue") +
  geom_point(data = subset(df_plot, highlight == "Ótimo"),
             aes(x = x, y = fx), color = "red", size = 3) +
  geom_vline(xintercept = x_k, linetype = "dashed") +
  annotate("text", x = x_k + 0.05, y = 0.29,
           label = "x[k] %~~% 3.591", angle = 0, parse = TRUE) +
  labs(x = "x", y = "f(x)", title = "") +
  theme_minimal() +
  theme(axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10))) +
  coord_cartesian(ylim = c(0, 0.3))

Figura 3.2.. Ponto de máximo estimado via método da bissecção no intervalo \([1,5]\).

O resultado obtido pelo método da bisseção indica que o valor aproximado de \(x_k\) que maximiza a função \(f(x)\) é aproximadamente \(3{,}591\). Nesse ponto, observa-se que a derivada da função, \(f'(x)\), se anula, confirmando a condição necessária de otimalidade para máximos locais em funções diferenciáveis. Tal resultado está em conformidade com a teoria, na qual pontos críticos — definidos por \(f'(x) = 0\) — correspondem a candidatos a máximos ou mínimos locais.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.3.2. Método de Newton-Rapshon

---

Diferentemente do método da bisseção apresentado na seção anterior, que apresenta convergência linear ao subdividir sucessivamente o intervalo onde ocorre a mudança de sinal de \(f'\), o método de Newton-Raphson utiliza informações locais da função por meio da razão entre \(f'(x)\) e \(f''(x)\), alcançando convergência — geralmente quadrática — quando a função é suficientemente diferenciável e a aproximação inicial é adequada. Essa maior eficiência, contudo, requer o cálculo das derivadas e pode ser prejudicada pela presença de pontos em que \(f''(x)\) se anula ou pela escolha inadequada do ponto inicial, fatores que podem provocar instabilidade ou falta de convergência. O conceito do método está formalizado na Definição 3.2 a seguir.

Definição 3.2 (Método de Newton-Raphson - Wallis 1685; Raphson 1697; Newton 1711; Ruggiero e Rocha Lopes 1998; Ehiwario e Aghamie 2014; Milne 2015). Seja \(g: \mathbb{R} \rightarrow \mathbb{R}\) uma função cuja derivada de primeira ordem, \(g'\), seja continuamente diferenciável e tal que \(g''(x^*) \neq 0\), onde \(x^*\) é o ponto crítico desejado. O objetivo é encontrar um valor de \(x\) que satisfaça a condição \(g'(x) = 0\). No entanto, quando a solução analítica é inviável, é necessário empregae uma abordagem numérica baseada em aproximações locais. Neste caso, na iteração \(t\), considera-se uma aproximação de \(g'(x^*)\) por meio de uma expansão de Taylor de primeira ordem centrada em \(x^{(t)}\):

\[\begin{align}\\ 0 = g'(x^*) \approx g'(x^{(t)}) + (x^* - x^{(t)})g''(x^{(t)})\\\\ \end{align}\]

Essa expressão representa a substituição da derivada por sua reta tangente no ponto \(x^{(t)}\), o que sugere que a raiz de \(g'\) pode ser aproximada pela raiz de sua linearização local. Resolvendo a equação para \(x^*\), obtém-se:

\[\begin{align}\\ x^* \approx x^{(t)} - \frac{g'(x^{(t)})}{g''(x^{(t)})} = x^{(t)} + h^{(t)}\\\\ \end{align}\]

em que \(h^{(t)} = -g'(x^{(t)})/g''(x^{(t)})\) é denominado incremento de Newton. Assim, a regra de atualização do método é dada por:

\[\begin{align}\\ x^{(t+1)} = x^{(t)} + h^{(t)}\\\\ \end{align}\]

Este processo de otimização é conhecido por método de Newton-Raphson, que, sob condições regulares, apresenta convergência local quadrática. Especificamente, se \(g'''\) é contínua e \(x^*\) é uma raiz simples de \(g'\), então existe uma vizinhança de \(x^*\) tal que a sequência \(\{x^{(t)}\}\), iniciada em qualquer ponto \(x^{(0)}\) dessa vizinhança, converge para \(x^*\). Além disso, o método converge para uma raiz de \(g'\) em qualquer ponto \(x^{(0)} \in [a, b]\), desde que sejam satisfeitas as seguintes condições:

• \(g''(x) \neq 0\) para todo \(x \in [a,b]\);
• \(g'''(x)\) não muda de sinal em \([a,b]\);
• \(g'(a)g'(b) < 0\);
• Para todo \(x\in [a,b]\), vale que \(\left|\dfrac{g'(x)}{g''(x)}\right| < b-a\).

Exemplo 3.2. Considere a função:

\[\begin{align}\\ f(x) = xe^{-x}\\\\ \end{align}\]

Suponha que o objetivo, neste caso, seja determinar o ponto de máximo da função, isto é, o valor de \(x\) tal que \(f'(x) = 0\), condição necessária para ocorrência de máximos ou mínimos locais. Note que esta função é diferenciável em todo o domínio real, mas a equação \(f'(x)=0\) não possui solução explícita simples, impossibilitando, portanto, a determinação exata do ponto de máximo por métodos puramente algébricos. Nesse cenário, torna-se necessário recorrer a procedimentos numéricos para localizar o ponto em que a derivada primeira se anula. Contudo, antes da aplicação de qualquer algoritmo iterativo, recomenda-se uma análise gráfica preliminar (Figura 3.3) da função \(f(x)\) e também de sua derivada \(f'(x)\), com o intuito de compreender o comportamento da função, identificar possíveis descontinuidades ou assimetrias, e delimitar regiões de interesse para a busca da solução.

Figura 3.3.. Ilustração do comportamento da função \(f(x)\) (painel esquerdo) e de \(f'(x)\) (painel direito) no intervalo \([0,2]\).

Nesta figura,observa-se que o máximo de \(f(x)\) ocorre em um ponto próximo de \(x = 1\). Com base nessa evidência empírica, é razoável escolher \(x^{(0)} = 0.5\) como ponto inicial para um processo iterativo de localização do ponto de máximo. O objetivo é ajustar progressivamente este valor até que \(f'(x)\) se anule ou se torne suficientemente pequeno em valor absoluto. As derivadas de \(f(x)\) são dadas por:

\[\begin{align}\\ f'(x) &= e^{-x}(1 - x)\\\\ f''(x) &= e^{-x}(x - 2)\\\\ \end{align}\]

Para a aplicação do método de Newton-Raphson, utiliza-se a regra iterativa:

\[\begin{align}\\ x^{(k+1)} = x^{(k)} - \dfrac{f'\left(x^{(k)}\right)}{f''\left(x^{(k)}\right)}\\\\ \end{align}\]

Assim, as iterações sucessivas, considerando \(x^{(0)} = 0.5\) como ponto inicial, são:

\[\begin{align}\\ x^{(1)} = x^{(0)} - \dfrac{f'\left(x^{(0)}\right)}{f''\left(x^{(0)}\right)} = 0.5000 + \dfrac{0.3033}{0.9098} = 0.8333\\\\ x^{(2)} = x^{(1)} - \dfrac{f'\left(x^{(1)}\right)}{f''\left(x^{(1)}\right)} = 0.8333 + \dfrac{0.0724}{0.5070} = 0.9762\\\\ x^{(3)} = x^{(2)} - \dfrac{f'\left(x^{(2)}\right)}{f''\left(x^{(2)}\right)} = 0.9762 + \dfrac{0.0090}{0.3857} = 0.9994\\\\ x^{(4)} = x^{(3)} - \dfrac{f'\left(x^{(3)}\right)}{f''\left(x^{(3)}\right)} = 0.9994 + \dfrac{0.0002}{0.3683} = 1.0000\\\\ x^{(5)} = x^{(4)} - \dfrac{f'\left(x^{(4)}\right)}{f''\left(x^{(4)}\right)} = 1.0000 + \dfrac{0.0000}{0.3679} = 1.0000\\\\ \end{align}\]

O resultado obtido indica que o valor aproximado de \(x^{(k)}\) que maximiza a função \(f(x)\) é aproximadamente \(1\). O Código 3.3 apresenta uma implementação do método da Newton-Raphson, na qual se fixa o número de iterações em itr = 40. Por fim, a Figura 3.4 ilustra o gráfico da função \(f(x)\), juntamente com o ponto de máximo estimado via método de Newton-Raphson, evidenciando a precisão da solução numérica obtida.

Código 3.3. Implementação do método de Newton-Raphson para encontrar o ponto de máximo da função \(f(x) = x e^{-x}\).

# ------------------------------------------
# Implementação do Método de Newton-Raphson
# ------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Valores iniciais ---

x    <- 0.5   # Ponto inicial
itr  <- 5     # Número de iterações

# --- 2. Definição da função objetivo ---

f    <- function(x)
{
  x * exp(-x)
}

# --- 3. Derivadas da função objetivo ---

df   <- function(x) 
{
  exp(-x) * (1 - x)
}

ddf  <- function(x)
{
  exp(-x) * (x - 2)
}

# --- 4. Implementação do método de Newton-Raphson ---

for (i in 1:itr) {
  cat(sprintf("Iteração %d:\n", i))
  cat(sprintf("  x^(k)       = %.8f\n", x))
  cat(sprintf("  f'(x^(k))   = %.8f\n", df(x)))
  cat(sprintf("  f''(x^(k))  = %.8f\n", ddf(x)))
  
  x <- x - df(x) / ddf(x)
  
  cat(sprintf("  Atualizado para x^(k+1) = %.8f\n\n", x))
}

Iteração 1:
  x^(k)       = 0.50000000
  f'(x^(k))   = 0.30326533
  f''(x^(k))  = -0.90979599
  Atualizado para x^(k+1) = 0.83333333

Iteração 2:
  x^(k)       = 0.83333333
  f'(x^(k))   = 0.07243303
  f''(x^(k))  = -0.50703124
  Atualizado para x^(k+1) = 0.97619048

Iteração 3:
  x^(k)       = 0.97619048
  f'(x^(k))   = 0.00897009
  f''(x^(k))  = -0.38571367
  Atualizado para x^(k+1) = 0.99944629

Iteração 4:
  x^(k)       = 0.99944629
  f'(x^(k))   = 0.00020381
  f''(x^(k))  = -0.36828701
  Atualizado para x^(k+1) = 0.99999969

Iteração 5:
  x^(k)       = 0.99999969
  f'(x^(k))   = 0.00000011
  f''(x^(k))  = -0.36787967
  Atualizado para x^(k+1) = 1.00000000

# --- 6. Visualização gráfica ---

x_vals            <- seq(0, 2, by = 0.0001)
df_plot           <- data.frame(x = x_vals,
                                fx = f(x_vals))

x_k               <- x
df_plot$highlight <- ifelse(abs(df_plot$x - x_k) < 1e-4, "Ótimo", "Não")

ggplot(df_plot, aes(x = x, y = fx)) +
  geom_line(size = 0.5, color = "steelblue") +
  geom_point(data = subset(df_plot, highlight == "Ótimo"),
             aes(x = x, y = fx), color = "red", size = 3) +
  geom_vline(xintercept = x_k, linetype = "dashed") +
  annotate("text", x = x_k + 0.05, y = max(df_plot$fx) * 0.95,
           label = "x[k] %~~% 1.000", angle = 0, parse = TRUE) +
  labs(x = "x", y = "f(x)", title = "") +
  theme_minimal() +
  theme(axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.4.. Ponto de máximo estimado via método de Newton-Raphson no intervalo \([0,2]\).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 3.3 (Distribuição Weibull - Weibull 1951). Considere a distribuição Weibull com função densidade de probabilidade dada por:

\[\begin{align}\\ f(x; \alpha, \beta) = \frac{\alpha}{\beta} \left( \frac{x}{\beta} \right)^{\alpha - 1} \exp\left\{ -\left( \frac{x}{\beta} \right)^\alpha \right\}, \quad x > 0, \ \alpha > 0, \ \beta > 0\\\\ \end{align}\]

em que \(\alpha\) representa o parâmetro de forma (shape) e \(\beta\) o parâmetro de escala (scale). Neste exemplo, o objetivo consiste em estimar o vetor de parâmetros \(\boldsymbol{\theta} = (\alpha, \beta)^\top\) por meio do método da máxima verossimilhança. Para tal, considere uma amostra \(x_1, x_2, \ldots, x_n\) independente e identicamente distribuída segundo a distribuição Weibull. A função de log-verossimilhança associada é dada por:

\[\begin{align}\\ \ell(\alpha, \beta; \mathbf{x}) = n \ln(\alpha) - n \alpha \ln(\beta) + (\alpha - 1) \sum_{i=1}^n \ln(x_i) - \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha\\\\ \end{align}\]

Para obter os estimadores de máxima verossimilhança, é necessário determinar, primeiramente, as derivadas de \(\ell(\alpha, \beta; \mathbf{x})\) com respeito a \(\alpha\) e \(\beta\), que compõe o vetor escore de Fisher, \(\mathbf{U}(\alpha, \beta; \mathbf{x})\), e resolver o sistema \(\mathbf{U}(\alpha, \beta; \mathbf{x}) = 0\). Para a distribuição Weibull, o vetor escore de Fisher é dado por:

\[\begin{align}\\ \mathbf{U}(\alpha, \beta; \mathbf{x}) = \begin{bmatrix} \displaystyle \frac{n}{\alpha} - n\ln (\beta) + \sum_{i=1}^n \ln(x_i) - \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha \ln\left( \frac{x_i}{\beta} \right) \\[10pt] \displaystyle -\frac{n\alpha}{\beta} + \frac{\alpha}{\beta} \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha \end{bmatrix}\\\\ \end{align}\]

Note que a equação \(\mathbf{U}(\alpha, \beta; \mathbf{x}) = 0\), não admite solução analítica, o que implica que os estimadores de máxima verossimilhança de \(\alpha\) e \(\beta\) devem ser determinados por métodos numéricos. Sendo assim, será empregado o método de Newton-Raphson que, no contexto de verossimilhança, atualiza iterativamente os parâmetros utilizando a seguinte regra:

\[\begin{align}\\ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - \left[ \nabla^2_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}^{(t)}; \mathbf{x}) \right]^{-1} \nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}^{(t)}; \mathbf{x})\\\\ \end{align}\]

em que \(\nabla_{\boldsymbol{\theta}} \ell(\cdot)\) representa o vetor escore e \(\nabla^2_{\boldsymbol{\theta}} \ell(\cdot)\), a matriz Hessiana (matriz das derivadas parciais de segunda ordem) da log-verossimilhança. No caso da distribuição Weibull, a Hessiana observada é dada por:

\[\begin{align}\\ \nabla^2_{\boldsymbol{\theta}} \ell(\alpha, \beta) = \begin{bmatrix} \displaystyle -\frac{n}{\alpha^2} - \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha \left[ \ln\left( \frac{x_i}{\beta} \right) \right]^2 & - \dfrac{n}{\beta} + \dfrac{1}{\beta} \displaystyle \sum_{i=1}^n \ln\left( \frac{x_i}{\beta} \right) \left( \frac{x_i}{\beta} \right)^{\alpha} \\[20pt] - \dfrac{n}{\beta} + \dfrac{1}{\beta} \displaystyle \sum_{i=1}^n \ln\left( \frac{x_i}{\beta} \right) \left( \frac{x_i}{\beta} \right)^{\alpha} & \displaystyle \frac{n\alpha}{\beta^2} - \frac{\alpha(\alpha+1)}{\beta^2} \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha \end{bmatrix}\\\\ \end{align}\]

Esta forma de atualização do método de Newton-Raphson, que emprega a matriz Hessiana observada da log-verossimilhança, é comumente denominada método de Newton-Raphson-Fisher (também referido, em alguns textos, como método do escore de Fisher — Fisher 1922; Longford 1987). Para fins ilustrativos, consideraremos uma amostra simulada proveniente de uma distribuição Weibull com parâmetros \(\alpha = 2\) e \(\beta = 1.5\), sobre a qual será implementado o procedimento iterativo de maximização da log-verossimilhança utilizando este método. O Código 3.4 apresenta a implementação detalhada desse processo, enquanto a Figura 3.5 ilustra o histograma da amostra simulada juntamente com a curva da função densidade Weibull ajustada.

Código 3.4. Implementação do método de Newton-Raphson-Fisher para a maximização da log-verossimilhança da distribuição Weibull.

# -------------------------------------------------------------------------------------------------------------------
# Implementação do Método da Newton-Raphson-Fisher: Estimadores de Máxima de Verossimilhança da Distribuição Weibull
# -------------------------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)
library(knitr)

# --- 1. Implementação da log-verossimilhança ---

logvero     <- function(n, x, theta)
{
  alpha     <- theta[1]   
  beta      <- theta[2]   
  l         <- n * log(alpha) - n * alpha * log(beta) + (alpha - 1) * sum(log(x)) - sum((x / beta)^alpha)
  return(l)
}

# --- 2. Implementação do vetor escore ---

G           <- function(n, x, theta)
{
  alpha     <- theta[1]
  beta      <- theta[2]
  A         <- x / beta
  U         <- numeric(2)
  U[1]      <- n / alpha - n * log(beta) + sum(log(x)) - sum(A^alpha * log(A))
  U[2]      <- -n * alpha / beta + alpha / beta * sum(A^alpha)
  return(U)
}

# --- 3. Implementação da hessiana ---

J           <- function(n, x, theta)
{
  alpha     <- theta[1]
  beta      <- theta[2]
  A         <- x / beta
  H         <- matrix(nrow = 2, ncol = 2)
  H[1,1]    <- -n / alpha^2 - sum(A^alpha * (log(A))^2)
  H[1,2]    <- n / beta - 1 / beta * sum(log(A) * A^alpha)
  H[2,1]    <- H[1,2]
  H[2,2]    <- n * alpha / beta^2 - (alpha * (alpha + 1) / beta^2) * sum(A^alpha)
  return(H)
}

# --- 4. Implementação do método de Newton-Raphson-Fisher ---

set.seed(12345)
n           <- 10000
alpha_0     <- 1.9
beta_0      <- 1.4
theta_0     <- matrix(c(alpha_0, beta_0), nrow = 2)
x           <- rweibull(n, shape = 2, scale = 1.5) # Dados simulados
erro        <- 0.005
itr         <- 1

while(itr < 100)
{
  # Atualização de Newton-Raphson-Fisher
  
  theta_1   <- theta_0 - solve(J(n, x, theta_0)) %*% G(n, x, theta_0)
  
  # Critério de parada
  
  if(max(abs(G(n, x, theta_1) - G(n, x, theta_0))) < erro) break
  
  # Impressão da iteração no formato padronizado
  
  cat(sprintf("Iteração %d:\n", itr))
  cat(sprintf("  alpha^(k)            = %.8f\n", theta_0[1]))
  cat(sprintf("  beta^(k)             = %.8f\n", theta_0[2]))
  cat(sprintf("  Log-verossimilhança  = %.8f\n", logvero(n, x, theta_0)))
  cat(sprintf("  Atualizado para:\n"))
  cat(sprintf("    alpha^(k+1)        = %.8f\n", theta_1[1]))
  cat(sprintf("    beta^(k+1)         = %.8f\n\n", theta_1[2]))
  
  theta_0   <- theta_1
  itr       <- itr + 1
}

Iteração 1:
  alpha^(k)            = 1.90000000
  beta^(k)             = 1.40000000
  Log-verossimilhança  = -10021.55223712
  Atualizado para:
    alpha^(k+1)        = 2.03280209
    beta^(k+1)         = 1.50113984

Iteração 2:
  alpha^(k)            = 2.03280209
  beta^(k)             = 1.50113984
  Log-verossimilhança  = -9937.98621331
  Atualizado para:
    alpha^(k+1)        = 2.00655111
    beta^(k+1)         = 1.49288942

Iteração 3:
  alpha^(k)            = 2.00655111
  beta^(k)             = 1.49288942
  Log-verossimilhança  = -9937.49448546
  Atualizado para:
    alpha^(k+1)        = 2.01942175
    beta^(k+1)         = 1.50003163

Iteração 4:
  alpha^(k)            = 2.01942175
  beta^(k)             = 1.50003163
  Log-verossimilhança  = -9937.35644527
  Atualizado para:
    alpha^(k+1)        = 2.01064808
    beta^(k+1)         = 1.49577095

Iteração 5:
  alpha^(k)            = 2.01064808
  beta^(k)             = 1.49577095
  Log-verossimilhança  = -9937.30942837
  Atualizado para:
    alpha^(k+1)        = 2.01629322
    beta^(k+1)         = 1.49860503

Iteração 6:
  alpha^(k)            = 2.01629322
  beta^(k)             = 1.49860503
  Log-verossimilhança  = -9937.28767176
  Atualizado para:
    alpha^(k+1)        = 2.01263267
    beta^(k+1)         = 1.49677875

Iteração 7:
  alpha^(k)            = 2.01263267
  beta^(k)             = 1.49677875
  Log-verossimilhança  = -9937.27906269
  Atualizado para:
    alpha^(k+1)        = 2.01501916
    beta^(k+1)         = 1.49796964

Iteração 8:
  alpha^(k)            = 2.01501916
  beta^(k)             = 1.49796964
  Log-verossimilhança  = -9937.27527167
  Atualizado para:
    alpha^(k+1)        = 2.01347224
    beta^(k+1)         = 1.49719703

Iteração 9:
  alpha^(k)            = 2.01347224
  beta^(k)             = 1.49719703
  Log-verossimilhança  = -9937.27371340
  Atualizado para:
    alpha^(k+1)        = 2.01447940
    beta^(k+1)         = 1.49769963

Iteração 10:
  alpha^(k)            = 2.01447940
  beta^(k)             = 1.49769963
  Log-verossimilhança  = -9937.27304352
  Atualizado para:
    alpha^(k+1)        = 2.01382566
    beta^(k+1)         = 1.49737319

Iteração 11:
  alpha^(k)            = 2.01382566
  beta^(k)             = 1.49737319
  Log-verossimilhança  = -9937.27276381
  Atualizado para:
    alpha^(k+1)        = 2.01425086
    beta^(k+1)         = 1.49758542

Iteração 12:
  alpha^(k)            = 2.01425086
  beta^(k)             = 1.49758542
  Log-verossimilhança  = -9937.27264479
  Atualizado para:
    alpha^(k+1)        = 2.01397467
    beta^(k+1)         = 1.49744753

Iteração 13:
  alpha^(k)            = 2.01397467
  beta^(k)             = 1.49744753
  Log-verossimilhança  = -9937.27259477
  Atualizado para:
    alpha^(k+1)        = 2.01415423
    beta^(k+1)         = 1.49753716

Iteração 14:
  alpha^(k)            = 2.01415423
  beta^(k)             = 1.49753716
  Log-verossimilhança  = -9937.27257357
  Atualizado para:
    alpha^(k+1)        = 2.01403756
    beta^(k+1)         = 1.49747891

Iteração 15:
  alpha^(k)            = 2.01403756
  beta^(k)             = 1.49747891
  Log-verossimilhança  = -9937.27256464
  Atualizado para:
    alpha^(k+1)        = 2.01411339
    beta^(k+1)         = 1.49751677

Iteração 16:
  alpha^(k)            = 2.01411339
  beta^(k)             = 1.49751677
  Log-verossimilhança  = -9937.27256086
  Atualizado para:
    alpha^(k+1)        = 2.01406412
    beta^(k+1)         = 1.49749217

Iteração 17:
  alpha^(k)            = 2.01406412
  beta^(k)             = 1.49749217
  Log-verossimilhança  = -9937.27255926
  Atualizado para:
    alpha^(k+1)        = 2.01409614
    beta^(k+1)         = 1.49750815

Iteração 18:
  alpha^(k)            = 2.01409614
  beta^(k)             = 1.49750815
  Log-verossimilhança  = -9937.27255859
  Atualizado para:
    alpha^(k+1)        = 2.01407533
    beta^(k+1)         = 1.49749776

Iteração 19:
  alpha^(k)            = 2.01407533
  beta^(k)             = 1.49749776
  Log-verossimilhança  = -9937.27255831
  Atualizado para:
    alpha^(k+1)        = 2.01408886
    beta^(k+1)         = 1.49750452

Iteração 20:
  alpha^(k)            = 2.01408886
  beta^(k)             = 1.49750452
  Log-verossimilhança  = -9937.27255819
  Atualizado para:
    alpha^(k+1)        = 2.01408007
    beta^(k+1)         = 1.49750013

Iteração 21:
  alpha^(k)            = 2.01408007
  beta^(k)             = 1.49750013
  Log-verossimilhança  = -9937.27255814
  Atualizado para:
    alpha^(k+1)        = 2.01408578
    beta^(k+1)         = 1.49750298

Iteração 22:
  alpha^(k)            = 2.01408578
  beta^(k)             = 1.49750298
  Log-verossimilhança  = -9937.27255811
  Atualizado para:
    alpha^(k+1)        = 2.01408207
    beta^(k+1)         = 1.49750113

Iteração 23:
  alpha^(k)            = 2.01408207
  beta^(k)             = 1.49750113
  Log-verossimilhança  = -9937.27255810
  Atualizado para:
    alpha^(k+1)        = 2.01408448
    beta^(k+1)         = 1.49750233

Iteração 24:
  alpha^(k)            = 2.01408448
  beta^(k)             = 1.49750233
  Log-verossimilhança  = -9937.27255810
  Atualizado para:
    alpha^(k+1)        = 2.01408291
    beta^(k+1)         = 1.49750155

Iteração 25:
  alpha^(k)            = 2.01408291
  beta^(k)             = 1.49750155
  Log-verossimilhança  = -9937.27255810
  Atualizado para:
    alpha^(k+1)        = 2.01408393
    beta^(k+1)         = 1.49750206

# --- 5. Estimativas finais ---

valores     <- c(theta_1[1], theta_1[2])
nomes       <- c("α", "β")

res_f       <- data.frame(Par = nomes,
                          Estimativa = valores,
                          Valor_Nominal = c(2.0, 1.5))

kable(res_f, 
      digits = 7,
      align = "c",
      format = "html", 
      escape = FALSE,
      col.names = c("Parâmetros <br> (α, β)", "Estimativa <br> (Newton-Raphson-Fisher)", "Valor Nominal <br> (Verdadeiro)"))

Parâmetros (α, β)	Estimativa (Newton-Raphson-Fisher)	Valor Nominal (Verdadeiro)
α	2.014083	2.0
β	1.497502	1.5

# --- 6. Histograma da amostra gerada + curva Weibull ajustada ---

# Amostra gerada

df          <- data.frame(x = x)

# Densidade Weibull ajustada

weibull_den <- function(x)
{
  dweibull(x, shape = theta_1[1], scale = theta_1[2])
}

# Histograma + Curva ajustada

ggplot(df, aes(x = x)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "lightgray", color = "black") +
  stat_function(fun = weibull_den, aes(color = "Curva Weibull"), size = 0.5) +
  scale_color_manual(name = NULL, values = c("Curva Weibull" = "blue")) +
  labs(x = "x", y = "Densidade", title = "") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.5.. Histograma da amostra gerada junto a curva ajustada da distribuição Weibull, via método de Newton-Raphson-Fisher.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.3.3. Recozimento Simulado (Simulated Annealing)

---

Diferentemente de métodos determinísticos como o algoritmo de Newton-Raphson, cuja trajetória de busca é rigidamente definida pelas propriedades locais da função objetivo e que, por conseguinte, pode ficar restrito a mínimos locais, o recozimento simulado (em ingles, simulated annealing - SA) baseia-se em princípios probabilísticos inspirados no processo físico de recozimento de metais, em que a temperatura do sistema é gradualmente reduzida para favorecer a organização dos átomos em estados de menor energia. No contexto de otimização, essa analogia traduz-se na aceitação controlada de soluções piores ao longo das iterações, o que confere ao algoritmo a capacidade de escapar de ótimos locais e, assim, aumentar a probabilidade de alcançar o ótimo global, sobretudo em funções multimodais ou de alta complexidade. Contudo, a eficiência do método está fortemente ligada à escolha da função de temperatura e à taxa de resfriamento, parâmetros que, se mal calibrados, podem comprometer tanto a qualidade da solução quanto o tempo computacional necessário. O fundamento do algoritmo encontra-se formalizado na Definição 3.3 a seguir.

Definição 3.3 (Recozimento Simulado - Simulated Annealing - Kirkpatrick, Gelatt Jr, e Vecchi 1983; Van Laarhoven et al. 1987; Bertsimas e Tsitsiklis 1993). Seja um problema de otimização cujo objetivo é encontrar:

\[\begin{align}\\ x^* = \arg\min_{x \ \in \ \mathcal{X}} f(x)\\\\ \end{align}\]

onde \(f: \mathcal{X} \rightarrow \mathbb{R}\) é uma função objetivo, possivelmente não-convexa e multimodal, e \(\mathcal{X}\) é o espaço de soluções admissíveis. O método de recozimento simulado (em ingles, simulated annealing - SA) é um algoritmo estocástico baseado em um processo de busca probabilística, inspirado no fenômeno físico do recozimento metálico (annealing), no qual um material é gradualmente resfriado para minimizar sua energia interna, buscando seu estado de menor energia (estado fundamental). No contexto da otimização, o SA gera uma cadeia de soluções \(\left\{x^{(t)}\right\}_{t=0}^\infty\), com atualização a partir da solução atual \(x^{(t)}\) por meio dos seguintes passos:

• 1º Passo (Ponto Inicial). Escolha um ponto inicial \(x^{(0)}\), e determine a função objetivo neste ponto. Caso não seja possível a escolha de um ponto inicial, a escolha deve ser feita de maneira aleatória.

• 2º Passo (Geração da Solução Candidata). A partir da solução atual \(x^{(t)}\), é amostrada uma solução candidata \(x' \in \mathcal{N}\left(x^{(t)}\right)\), onde \(\mathcal{N}\left(x^{(t)}\right)\) denota a vizinhança de \(x^{(t)}\). Esta geração é feita segundo uma distribuição de proposta \(q(x^{(t)}, \cdot)\), que deve garantir irredutibilidade e a possibilidade de explorar todo o espaço \(\mathcal{X}\).

• 3º Passo (Avaliação da Variação do Custo). \[\begin{align}\\ \Delta f = f(x') - f\left(x^{(t)}\right)\\\\ \end{align}\]

• 4º Passo (Regra de Aceitação). A solução candidata é aceita com probabilidade: \[\begin{align}\\ P_{\text{aceitação}}(x^{(t)}, x') = \begin{cases} 1, & \text{se } \Delta f \leqslant 0, \\ \exp\left\{-\Delta f/T^{(t)}\right\}, & \text{se } \Delta f > 0, \end{cases}\\\\ \end{align}\] onde \(T^{(t)} > 0\) é a temperatura no passo \(t\), que decresce monotonamente segundo um cronograma de resfriamento (em inglês, cooling schedule), por exemplo, \[\begin{align}\\ T^{(t)} = \frac{T_0}{\log(t + t_0)}\\\\ \end{align}\] com constantes \(T_0 > 0\) e \(t_0 \geqslant 1\), e com a propriedade \(T^{(t)} \downarrow 0\) quando \(t \to \infty\).

• 5º Passo (Atualização da Solução). \[\begin{align}\\ x^{(t+1)} = \begin{cases} x', & \text{com probabilidade } P_{\text{aceitação}}(x^{(t)}, x'), \\ x^{(t)}, & \text{caso contrário} \end{cases}\\\\ \end{align}\]

Este procedimento, em particular, define uma cadeia de Markov homogênea para cada temperatura fixa \(T^{(t)}\), cuja distribuição estacionária é a distribuição de Gibbs dada por:

\[\begin{align}\\ \pi_T(x) = \frac{\exp\left\{-f(x)/T\right\}}{Z(T)}\\\\ \end{align}\]

onde \(Z(T) = \sum_{x \ \in \ \mathcal{X}} \exp\left\{-f(x)/T\right\}\) é a constante de normalização (função partição). O algoritmo progressivamente reduz a temperatura, aproximando a distribuição \(\pi_T\) para concentrações cada vez maiores em torno dos mínimos globais de \(f\). É importante destacar que a eficácia do simulated annealing depende criticamente da escolha do cronograma de resfriamento, da estrutura da vizinhança \(\mathcal{N}(\cdot)\) e da forma de geração das soluções candidatas, que devem permitir um equilíbrio adequado entre exploração global — permitindo escapar de mínimos locais — e refinamento local — garantindo a convergência para soluções de qualidade.

Exemplo 3.4. Considere a função objetivo definida por:

\[\begin{align}\\ \ell(\alpha) = \sum_{i=1}^{n} \ln\left[\beta^2 + (x_i - \alpha)^2\right]\\\\ \end{align}\]

onde \(\beta = 0{.}1\) é um parâmetro fixado e \(x_1, \ldots, x_n \in \mathbb{R}\) representam observações de uma amostra. O objetivo consiste em estimar o valor do parâmetro \(\alpha \in \mathbb{R}\) que minimiza a função \(\ell(\alpha)\), ou seja,

\[\begin{align}\\ \alpha^* = \arg\min_{\alpha \ \in \ \mathbb{R}} \sum_{i=1}^{n} \ln\left[\beta^2 + (x_i - \alpha)^2\right]\\\\ \end{align}\]

Como \(\ell(\alpha)\) apresenta uma superfície de otimização não-convexa e com múltiplos mínimos locais, a aplicação de métodos determinísticos de otimização local podem ser ineficientes neste caso. Então, para lidar com essa dificuldade, pode-se recorrer a métodos estocásticos globais, como o algoritmo simulated annealing. Para tal, considere a amostra observada:

\[\begin{align}\\ x = (-4{.}20,\, -2{.}85,\, -2{.}30,\, 1{.}02,\, 0{.}70,\, 0{.}98,\, 2{.}72,\, 3{.}50)\\\\ \end{align}\]

com \(n = 8\) observações. O algoritmo simulated annealing, neste caso, será implementado com os seguintes parâmetros de controle:

• Temperatura Inicial: \(T^{(0)} = 20\),
• Fator de Resfriamento: \(\rho = 0{.}9\),
• Valor Inicial: \(\alpha^{(0)} = 4\),
• Número de Iterações: \(100\).

A cada iteração \(t\), uma nova proposta \(\alpha'\) é gerada a partir da atual \(\alpha^{(t)}\), utilizando uma distribuição de proposta simétrica (por exemplo, a distribuição normal centrada em \(\alpha^{(t)}\)). A aceitação desta proposta ocorre com uma probabilidade que depende tanto da variação na função objetivo quanto da temperatura vigente, mecanismo que possibilita, sobretudo nas fases iniciais do algoritmo, a aceitação ocasional de soluções com maior valor de custo. Tal característica favorece a exploração do espaço de busca e reduz a probabilidade de aprisionamento em mínimos locais. À medida que a temperatura diminui progressivamente, o algoritmo passa a priorizar soluções com menores valores de \(\ell(\alpha)\), conduzindo a um processo de refinamento local em torno de regiões promissoras do espaço paramétrico. O Código 3.5 apresenta a implementação deste procedimento no ambiente R, enquanto a Figura 3.6 ilustra a trajetória percorrida pelo algoritmo e as estimativas obtidas para a função objetivo.

Código 3.5. Implementação do algoritmo simulated annealing para a minimização de \(\ell(\alpha)\).

# --------------------------------------------------------------------------------------------------
# Implementação do Simulated Annealing: Estimação do Parâmetro α via Minimização da Função Objetivo
# --------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)
library(patchwork)

# --- 1. Definição da função objetivo ---

ell             <- function(alpha, x, beta = 0.1) 
{
  sum(log(beta^2 + (x - alpha)^2))
}

# --- 2. Definição dos parâmetros iniciais ---

x               <- c(-4.20, -2.85, -2.30, 1.02, 0.70, 0.98, 2.72, 3.50) # Amostra
beta            <- 0.1
alpha           <- 4      # Valor inicial
T0              <- 20     # Temperatura inicial
rho             <- 0.9    # Fator de resfriamento
n_iter          <- 100    # Número de iterações
sd_prop         <- 1      # Desvio-padrão da distribuição proposta

# --- 3. Algoritmo simulated annealing ---

set.seed(2025)
alphas          <- costs <- numeric(n_iter + 1)
alphas[1]       <- alpha
costs[1]        <- ell(alpha, x, beta)

for (t in 1:n_iter) 
{
  T             <- T0 * rho^(t - 1)
  prop          <- alpha + rnorm(1, 0, sd_prop)
  delta         <- ell(prop, x, beta) - ell(alpha, x, beta)
  if (delta <= 0 || runif(1) < exp(-delta / T)) {
    alpha       <- prop
  }
  
  alphas[t + 1] <- alpha
  costs[t + 1]  <- ell(alpha, x, beta)
}

# --- 4. Determinação do Mínimo Global da Função Objetivo ---

alpha_grid      <- seq(min(x) - 2, max(x) + 2, length.out = 1000)
ll_grid         <- sapply(alpha_grid, ell, x = x, beta = beta)
global_idx      <- which.min(ll_grid)


# --- 5. Valor Ótimo de α Obtido pelo Simulated Annealing ---

best_idx        <- which.min(costs)

res_f           <- data.frame(Estimativa = alphas[best_idx],
                              Funcao_Obj = ell(alpha = alphas[best_idx], x = x, beta = beta),
                              Minimo = alpha_grid[global_idx])

kable(res_f, 
      digits = 7,
      align = "c",
      format = "html", 
      escape = FALSE,
      col.names = c("Valor Ótimo <br> de α", "Função Objetivo no <br> Valor Ótimo de α", "Mínimo <br> Global"))

Valor Ótimo de α	Função Objetivo no Valor Ótimo de α	Mínimo Global
1.000492	-0.1262409	0.9792793

# --- 6. Visualização Gráfica: Trajetória e Valor Ótimo de α ---

df_traj         <- data.frame(Iter = 0:n_iter, Alpha = alphas)
df_obj          <- data.frame(Alpha = alpha_grid, LL = ll_grid)
df_best         <- data.frame(Alpha = alphas[best_idx], LL = costs[best_idx])
df_global       <- data.frame(Alpha = alpha_grid[global_idx], LL = ll_grid[global_idx])
df_sa           <- data.frame(Alpha = alphas, LL = costs)

p1              <- ggplot(df_traj, aes(Iter, Alpha)) +
                    geom_line(color = "blue", size = 0.5) +
                    geom_hline(yintercept = df_global$Alpha, 
                               linetype = "dashed", 
                               color = "darkgreen", 
                               size = 0.5) +
                    annotate("point", x = best_idx - 1, 
                             y = df_best$Alpha, 
                             color = "red", 
                             size = 1.5) +
                    labs(title = "",
                         x = "Nº de Iterações",
                         y = expression(alpha)) +
                    theme_minimal() + 
                    theme(axis.title.x = element_text(margin = margin(t = 10)),
                          axis.title.y = element_text(margin = margin(r = 10)))

p2              <- ggplot() +
                    geom_line(data = df_obj, 
                              aes(Alpha, LL), 
                              color = "gray50", 
                              size = 0.5) +
                    geom_point(data = df_sa, 
                               aes(Alpha, LL, color = "Trajetória"), 
                               size = 0.5, 
                               alpha = 0.6) +
                    geom_point(data = df_best, 
                               aes(Alpha, LL, color = "Melhor Estimativa"), 
                               size = 1.5) +
                    geom_point(data = df_global, 
                               aes(Alpha, LL, color = "Mínimo Global"), 
                               size = 1.5,
                               shape = 4,
                               stroke = 1.2) +
                    scale_color_manual(values = c("Trajetória" = "blue",
                                                  "Melhor Estimativa" = "red",
                                                  "Mínimo Global" = "darkgreen")) +
                    labs(title = "",
                         x = expression(alpha),
                         y = expression(l(alpha)),
                         color = "Legenda") +
                    theme_minimal() + 
                    theme(axis.title.x = element_text(margin = margin(t = 10)),
                          axis.title.y = element_text(margin = margin(r = 10)))

p1 + p2

Figura 3.6.. Trajetória de α percorrida no simulated annealing (painel esquerdo) e as estimativas obtidas de α para a função objetivo no intervalo \([-6, 8]\) (painel direito).

Após a execução do procedimento, o algoritmo produziu a solução \(\widehat{\alpha} = 1{.}000492\), à qual corresponde o valor da função objetivo \(\ell(\widehat{\alpha}) = -0{.}1262409\). Embora a função objetivo apresente múltiplos mínimos locais, a estimativa obtida encontra-se bastante próxima do mínimo global, localizado em \(\alpha^{*} = 0{.}9792793\). Ademais, o gráfico resultante evidencia a trajetória da sequência \(\left\{\alpha^{(t)}\right\}_{t=1}^{100}\), evidenciando, em um primeiro momento, uma ampla exploração do espaço paramétrico, seguida por uma progressiva estabilização em torno do valor ótimo estimado.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.3.4. Algoritmo Genético

---

O algoritmo genético (em inglês, genetic algorithm – GA), assim como o recozimento simulado (em inglês, simulated annealing – SA), fundamenta-se em princípios probabilísticos inspirados em processos naturais, permitindo uma exploração mais ampla do espaço de busca. No caso do GA, tais princípios derivam dos mecanismos evolutivos observados na natureza, em particular da seleção natural postulada por Darwin. Nesse paradigma, soluções candidatas ao problema de otimização são representadas como indivíduos de uma população, os quais evoluem ao longo das gerações mediante operadores como seleção, cruzamento (ou recombinação) e mutação, visando à preservação e combinação de características favoráveis. No contexto de otimização, essa analogia biológica traduz-se na capacidade do algoritmo de escapar de ótimos locais e de explorar regiões extensas e complexas do espaço de busca, assim como o SA. Essses algoritmos foram estudados, inicialmente, por J. H. Holland em 1975 (Holland 1975), o qual formalizou a noção de esquemas (schemata) e estabeleceu bases teóricas para o entendimento do comportamento dinâmico dessas populações de soluções.

3.3.4.1. Terminologia Básica

---

O entendimento adequado dos algoritmos genéticos exige familiaridade com a terminologia biológica que fundamenta sua estrutura conceitual. Nesta seção, serão apresentados os principais elementos que compõem a representação das soluções, estabelecendo a analogia entre conceitos da biologia molecular e seus equivalentes computacionais. Serão definidos, em particular, os termos cromossomo (Definição 3.4), que corresponde à estrutura codificadora de uma solução; gene (Definição 3.5), que corresponde a unidade elementar de informação.

Definição 3.4 (Cromossomo - Sivanandam et al. 2008; Givens e Hoeting 2012; Kramer e Kramer 2017). Seja \(\mathcal{S} \subseteq \mathbb{R}^d\) o espaço de busca de um problema de otimização. Um cromossomo é uma estrutura de dados que codifica uma solução candidata \(x \in \mathcal{S}\) por meio de uma cadeia finita de símbolos extraídos de um alfabeto \(\Sigma\), usualmente \(\Sigma = {0,1}\) (representação binária), mas podendo também ser \(\Sigma = \mathbb{Z}\) ou \(\Sigma = \mathbb{R}\) em codificações inteiras ou reais, respectivamente. Formalmente, um cromossomo é um vetor \(c = (c_1, c_2, \dots, c_L) \in \Sigma^L\), onde \(L\) é o comprimento fixo da codificação, e existe uma função de decodificação \(\phi: \Sigma^L \to \mathcal{S}\) tal que \(\phi(c) = x\), com \(x\) sendo a solução representada. Os cromossomos são os elementos sobre os quais atuam os operadores genéticos — avaliação de aptidão, seleção, recombinação (crossover) e mutação — no processo evolucionário. Ao longo de sucessivas gerações, espera-se que a população de cromossomos evolua, produzindo representações de soluções progressivamente mais aptas no espaço \(\mathcal{S}\).

Exemplo 3.5. Considere o problema de otimização univariada:

\[\begin{align}\\ \text{Maximizar } f(x) = x \cdot \sin(x), \quad \text{para } \quad x \in [0, 10]\\\\ \end{align}\]

O comportamento dessa função objetivo no intervalo \([0, 10]\) é exibido na Figura 3.7, onde podem ser observadas suas oscilações e picos ao longo do domínio considerado. Nota-se a existência de múltiplos pontos de máximo, incluindo um máximo global e pelo menos um máximo local, o que caracteriza a função como multimodal neste intervalo.

Figura 3.7.. Ilustração do comportamento da função \(f(x)\) no intervalo \([0,10]\).

Nosso objetivo, neste exemplo, é ilustrar o comportamento dos cromossomos em um algoritmo genético com representação binária. Para tanto, adotamos uma codificação de \(x\) em uma cadeia binária de comprimento \(L = 10\), de modo que cada cromossomo \(c \in {0,1}^{10}\) representa um número inteiro no intervalo de \(0\) até \(2^{10} - 1 = 1023\) (pois, com 10 bits, existem \(2^{10} = 1024\) combinações binárias distintas, que representam números inteiros no intervalo \([0, 1023]\)). A função de decodificação, denotada por \(\phi: {0,1}^{10} \to [0,10]\), é definida como:

\[\begin{align}\\ \phi(c) = \frac{\mathrm{int}(c)}{1023} \cdot 10\\\\ \end{align}\]

onde \(\mathrm{int}(c)\) é o valor decimal correspondente à cadeia binária \(c\). Por exemplo, considere o cromossomo, \(c^{(1)}\), definido por:

\[\begin{align}\\ c^{(1)} = (0, 1, 1, 0, 0, 1, 0, 1, 1, 0)\\\\ \end{align}\]

cujo valor decimal é calculado como:

\[\begin{align}\\ \mathrm{int}(c^{(1)}) &= 0 \cdot 2^9 + 1 \cdot 2^8 + 1 \cdot 2^7 + 0 \cdot 2^6 + 0 \cdot 2^5 + 1 \cdot 2^4 + 0 \cdot 2^3 + 1 \cdot 2^2 + 1 \cdot 2^1 + 0 \cdot 2^0 \\\\ &= 256 + 128 + 0 + 0 + 16 + 0 + 4 + 2 + 0 \\\\ &= 294\\\\ \end{align}\]

Para esse cromossomo, a solução real correspondente é então:

\[\begin{align}\\ x^{(1)} = \phi\left(c^{(1)}\right) = \frac{294}{1023} \cdot 10 \approx 2{.}874\\\\ \end{align}\]

O valor da função objetivo para esse indivíduo é:

\[\begin{align}\\ f\bigl(x^{(1)}\bigr) &= x^{(1)} \cdot \sin\bigl(x^{(1)}\bigr) \\\\ &\approx 2{.}874 \times \sin(2{.}874) \\\\ &\approx 2{.}874 \times 0{.}267 \\\\ &\approx 0{.}767\\\\ \end{align}\]

Considere, agora, um cromossomo \(c^{(2)}\) definido por:

\[\begin{align}\\ c^{(2)} = (1, 1, 0, 0, 1, 0, 0, 1, 1, 1)\\\\ \end{align}\]

cujo valor decimal é calculado como:

\[\begin{align}\\ \mathrm{int}(c^{(2)}) &= 1 \cdot 2^9 + 1 \cdot 2^8 + 0 \cdot 2^7 + 0 \cdot 2^6 + 1 \cdot 2^5 + 0 \cdot 2^4 + 0 \cdot 2^3 + 1 \cdot 2^2 + 1 \cdot 2^1 + 1 \cdot 2^0 \\\\ &= 512 + 256 + 0 + 0 + 32 + 0 + 0 + 4 + 2 + 1\\\\ &= 807\\\\ \end{align}\]

Para esse cromossomo, a solução real correspondente é então:

\[\begin{align}\\ x^{(2)} = \phi(c^{(2)}) = \frac{807}{1023} \cdot 10 \approx 7{.}8885\\\\ \end{align}\]

O valor da função objetivo para esse indivíduo é:

\[\begin{align}\\ f\bigl(x^{(2)}\bigr) &= x^{(2)} \cdot \sin\bigl(x^{(2)}\bigr) \\\\ &\approx 7{.}8885 \times \sin(7{.}8885) \\\\ &\approx 7{.}8885 \times 0{.}9994 \\\\ &\approx 7{.}8838\\\\ \end{align}\]

Observando os cromossomos analisados, nota-se que o primeiro cromossomo, cujo valor decodificado é aproximadamente \(x^{(1)} \approx 2{.}874\), encontra-se próximo a um máximo local da função \(f(x) = x \cdot \sin(x)\) no intervalo \([0,10]\). Já o segundo cromossomo, com valor decodificado em torno de \(x^{(2)} \approx 7{.}8838\), situa-se próximo do máximo global da função dentro do mesmo intervalo. Essa distinção ilustra como diferentes representações cromossômicas podem mapear regiões distintas do espaço de busca, ressaltando a capacidade do algoritmo genético de explorar múltiplos ótimos locais e potencialmente aproximar-se do ótimo global. \[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Definição 3.5 (Gene - Sivanandam et al. 2008; Givens e Hoeting 2012; Kramer e Kramer 2017). Seja \(\Sigma\) um alfabeto finito ou infinito, por exemplo, \(\Sigma = {0,1}\) (no caso binário), \(\Sigma = \mathbb{Z}\) (no caso inteiro) ou \(\Sigma = \mathbb{R}\) (no caso real). Um gene é definido como um elemento individual \(c_i \in \Sigma\) que ocupa a \(i\)-ésima posição de um cromossomo \(c = (c_1, c_2, \dots, c_L) \in \Sigma^L\), representando a unidade básica de informação codificada utilizada na descrição de uma solução candidata.

Exemplo 3.6. Considere o cromossomo binário:

\[\begin{align}\\ c = (0, 1, 1, 0, 0, 1, 0, 1, 1, 0) \in \{0,1\}^{10}\\\ \end{align}\]

Cada posição \(i\) do vetor \(c\) corresponde a um gene \(c_i\). Por exemplo:

\[\begin{align}\\ c_1 &= 0, \\\\ c_2 &= 1, \\\\ c_3 &= 1, \\\\ &\hspace{0.3cm}\vdots \\\\ c_{10} &= 0\\\\ \end{align}\]

No contexto da representação binária, cada gene \(c_i\) assume um valor em \(\Sigma = \{0,1\}\), refletindo a presença ou ausência de uma característica codificada naquele locus específico. Suponha que, em um ciclo particular do algoritmo genético, ocorra uma mutação no quinto gene, alterando seu valor original de \(0\) para \(1\). Essa alteração resulta em um novo cromossomo,

\[\begin{align}\\ c' = (0, 1, 1, 0, 1, 1, 0, 1, 1, 0)\\\\ \end{align}\]

Neste cenário, o gene \(c_5\) é a unidade mínima de informação que foi modificada pela operação de mutação, demonstrando seu papel central como elemento fundamental de variação genética. Essa modificação pontual pode gerar novas características na solução codificada, potencialmente impactando a aptidão do indivíduo na população e contribuindo para a exploração do espaço de busca pelo algoritmo genético.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Definição 3.6 (Alelo - Sivanandam et al. 2008; Givens e Hoeting 2012; Kramer e Kramer 2017). Seja \(\Sigma\) um alfabeto finito ou infinito, cujos elementos são os possíveis valores que um gene pode assumir em uma dada posição do cromossomo. O alelo é, portanto, um valor específico \(a \in \Sigma\) que representa a variante do gene naquela posição, correspondendo ao elemento do alfabeto utilizado na codificação.

Exemplo 3.7. Considere um cromossomo binário:

\[\begin{align}\\ c = (0, 1, 1, 0, 0, 1, 0, 1, 1, 0) \in \{0,1\}^{10}\\\ \end{align}\]

O alfabeto utilizado é \(\Sigma = \{0,1\}\). Cada gene \(c_i\) assume um valor que é um alelo, ou seja, um elemento de \(\Sigma\). Por exemplo, os alelos presentes neste cromossomo são apenas \(0\) e \(1\), que são os possíveis valores do gene em cada posição. Se, por outro lado, considerarmos uma codificação inteira, onde \(\Sigma = \{0,1,2,3,4,5\}\), um gene pode assumir qualquer um desses seis valores, e cada valor seria um alelo distinto. Este exemplo evidencia, em particular, que os alelos definem as possíveis variantes para cada gene, enquanto os genes representam posições fixas no cromossomo que armazenam essas variantes. A Figura 3.8 ilustra de forma mais clara a distinção entre esses conceitos.

Figura 3.8.. Representação de um cromossomo biológico e sua analogia na computação, ilustrando a correspondência entre cromossomo, gene e alelo.

Fonte: Genetic Algorithm | Urso, A., Fiannaca, A., La Rosa, M., Ravì, V., & Rizzo, R. (2019). Data Mining: Prediction Methods. Encyclopedia of Bioinformatics and Computational Biology.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.3.4.2. Estrutura do Algoritmo Genético

---

A partir das definições anteriores de cromossomo, gene e alelo, pode-se agora formalizar os conceitos de genótipo e fenótipo, essenciais para a compreensão da estrutura dos algoritmos genéticos. O genótipo, denotado por \(\vartheta\), representa a informação codificada em um indivíduo, isto é, corresponde à sequência ordenada dos alelos que compõem o cromossomo. Já o fenótipo, denotado por \(\theta\), é a manifestação observável dessa informação codificada, ou seja, a interpretação do genótipo no espaço do problema, representando a solução real associada.

Com base nesses conceitos, a estrutura do algoritmo genético pode ser entendida como um procedimento iterativo, cujas iterações são indexadas por um parâmetro discreto \(t\), que representa a geração atual. Suponha que a \(t\)-ésima geração seja composta por uma população de tamanho \(P\), formada pelos genótipos \(\vartheta_1^{(t)}, \dots, \vartheta_P^{(t)}\), correspondentes aos fenótipos \(\theta_1^{(t)}, \dots, \theta_P^{(t)}\). Assim, o algoritmo genético evolui iterativamente essa população ao longo das gerações indicadas por \(t\).

A seleção natural Darwiniana (Darwin 1859) orienta o processo evolutivo, favorecendo os indivíduos com maior aptidão (em inglês, fitness). Em termos biológicos, essa aptidão representa a pressão ambiental, classificando os organismos conforme sua adaptação ao ambiente. Já no contexto algoritmo genético, a aptidão de cada genótipo \(\vartheta_i^{(t)}\) é determinada pelo valor da função objetivo avaliada no respectivo fenótipo \(\theta_i^{(t)}\). A função fitness, que quantifica essa medida de aptidão, é apresentada na Definição 3.7.

Definição 3.7 (Função Fitness - Sivanandam et al. 2008; Givens e Hoeting 2012; Kramer e Kramer 2017). Seja \(\mathcal{F}: \mathcal{S} \to \mathbb{R}\\\) a função fitness, onde \(\mathcal{S}\) denota o espaço das soluções candidatas (fenótipos). Para cada fenótipo \(\theta \in \mathcal{S}\), a função fitness é dada por:

\[\begin{align}\\ \mathcal{F}(\theta) = f(\theta)\\\\ \end{align}\]

onde \(f\) é a função objetivo do problema de otimização. No caso de problemas de maximização, a função fitness é definida diretamente pelo valor da função objetivo. Já para problemas de minimização, a função fitness pode ser obtida por meio de transformações algébricas, tais como, por exemplo,

\[\begin{align}\\ \mathcal{F}(\theta) = \frac{1}{1 + f(\theta)}\\\\ \end{align}\]

assumindo que \(f(\theta) \geqslant 0\), para garantir que maiores valores de fitness indiquem soluções melhores.

Exemplo 3.8. Considere o problema de otimização univariada:

\[\begin{align}\\ \max_{x \ \in \ [0,10]} f(x) = x \cdot \sin(x)\\\\ \end{align}\]

Note que a função objetivo \(f\) apresenta múltiplos máximos locais devido à sua natureza oscilatória, caracterizando um problema multimodal. Para este caso, a função fitness pode ser definida diretamente como:

\[\begin{align}\\ \mathcal{F}(x) = f(x) = x \cdot \sin(x)\\\\ \end{align}\]

Assim, soluções cujo valor de \(f(x)\) seja maior indicam indivíduos com maior aptidão (fitness), que têm maior probabilidade de serem selecionados nas operações do algoritmo genético.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 3.9. Considere o problema de minimização:

\[\begin{align}\\ \min_{x \ \in \ [0,10]} g(x) = (x - 5)^2\\\\ \end{align}\]

Para este exemplo, para que a função fitness reflita corretamente a qualidade das soluções — isto é, soluções com menor custo tenham maior aptidão — é necessário realizar uma transformação da função objetivo. Assim, assumindo que \(g(x) \geqslant 0\) para todo \(x \in [0,10]\), uma transformação comum é:

\[\begin{align}\\ \mathcal{F}(x) = \frac{1}{1 + g(x)} = \frac{1}{1 + (x - 5)^2}\\\\ \end{align}\]

Neste caso, valores menores de \(g(x)\) geram valores maiores de \(\mathcal{F}(x)\), indicando maior aptidão e favorecendo soluções próximas a \(x = 5\).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Os algoritmos genéticos utilizam mecanismos de seleção para identificar e privilegiar os indivíduos de maior aptidão dentro da população. Esses mecanismos, baseados na avaliação da função fitness, selecionam preferencialmente as soluções que apresentam maior qualidade, assegurando que apenas os indivíduos mais adaptados prossigam para o processo de reprodução. A reprodução, por sua vez, consiste na transmissão total ou parcial do material genético desses indivíduos por meio da aplicação dos operadores genéticos, como o crossover e a mutação, os quais promovem diversidade genética e ampliam a capacidade do algoritmo em explorar eficientemente o espaço de busca.

Definição 3.8 (Crossover - Sivanandam et al. 2008; Givens e Hoeting 2012; Kramer e Kramer 2017). O crossover é um procedimento de recombinação genética que gera novos indivíduos (descendentes) a partir da combinação do material genético de dois ou mais indivíduos preexistentes (progenitores). Formalmente, dado um par de cromossomos progenitores, \(c^{(p_1)}\) e \(c^{(p_2)}\), o operador crossover produz um ou mais descendentes \(c^{(d)}\) por meio da troca e reorganização dos genes entre esses progenitores, preservando características genéticas e possibilitando a exploração de novas regiões do espaço de busca. Dentre os esquemas mais comuns de crossover destacam-se:

• Crossover de um Ponto: Consiste na escolha aleatória de um ponto de corte ao longo da sequência de genes dos progenitores. A partir desse ponto, os segmentos dos cromossomos são trocados entre os indivíduos para formar os descendentes. Por exemplo, em cromossomos binários, o primeiro descendente recebe os genes do primeiro progenitor até o ponto de corte e, a partir daí, os genes do segundo progenitor; o segundo descendente é formado da maneira inversa.

• Crossover Uniforme: Cada gene dos descendentes é escolhido aleatoriamente, com igual probabilidade, entre os genes correspondentes dos progenitores. Esse esquema promove uma maior mistura genética e diversidade nas soluções geradas.

• Crossover Aritmético: Apropriado para problemas de otimização numérica, este operador gera descendentes por meio de combinações lineares dos vetores genéticos dos progenitores. Isto é, dados dois indivíduos \(\vartheta_1\) e \(\vartheta_2\), os descendentes \(\vartheta_1'\) e \(\vartheta_2'\) são calculados como: \[\begin{align}\\ \vartheta_1' &= a \vartheta_1 + (1 - a) \vartheta_2 \\\\ \vartheta_2' &= (1 - a) \vartheta_1 + a \vartheta_2 \\\\ \end{align}\] onde \(a \in [0,1]\) é um número aleatório, geralmente amostrado de uma distribuição uniforme.

Exemplo 3.9 (Crossover de um Ponto em Cromossomos Binários). Considere dois cromossomos progenitores, cada um representado por uma cadeia binária de comprimento 8:

\[\begin{align}\\ c^{(p_1)} &= (1, 0, 1, 1, 0, 0, 1, 0) \\\\ c^{(p_2)} &= (0, 1, 0, 0, 1, 1, 0, 1) \\\\ \end{align}\]

Suponha que seja escolhido aleatoriamente o ponto de corte após o terceiro gene. O operador crossover de um ponto troca os segmentos após esse ponto entre os progenitores para formar os descendentes, isto é,

\[\begin{align}\\ c^{(d_1)} &= \big(1, 0, 1,\; {\color{red}{0, 1, 1, 0, 1}}\big) \\\\ c^{(d_2)} &= \big(0, 1, 0,\; {\color{blue}{1, 0, 0, 1, 0}}\big) \\\\ \end{align}\]

ou seja, o primeiro descendente mantém os três primeiros genes do primeiro progenitor e incorpora os cinco últimos genes do segundo progenitor; o segundo descendente faz o inverso.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 3.10 (Crossover Uniforme em Cromossomos Binários). Considere os mesmos progenitores do exemplo anterior:

\[\begin{align}\\ c^{(p_1)} &= (1, 0, 1, 1, 0, 0, 1, 0) \\\\ c^{(p_2)} &= (0, 1, 0, 0, 1, 1, 0, 1) \\\\ \end{align}\]

No crossover uniforme, cada gene dos descendentes é selecionado aleatoriamente, com igual probabilidade, entre os genes correspondentes dos progenitores. Um possível resultado para os descendentes é:

\[\begin{align}\\ c^{(d_1)} &= (1, 1, 1, 0, 0, 1, 1, 1) \\\\ c^{(d_2)} &= (0, 0, 0, 1, 1, 0, 0, 0) \\\\ \end{align}\]

Essa seleção gene a gene permite uma mistura genética mais detalhada, aumentando a diversidade da população e a capacidade de exploração do espaço de soluções.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 3.11 (Crossover Aritmético em Cromossomos Numéricos). Considere dois indivíduos representados por vetores numéricos, cujos cromossomos correspondem às coordenadas no espaço da solução:

\[\begin{align}\\ \vartheta_1 &= (0.8, 1.2, 0.5)\\\\ \vartheta_2 &= (1.4, 0.6, 1.0)\\\\ \end{align}\]

Para realizar o crossover aritmético, escolhe-se um parâmetro \(a \in [0,1]\), aqui fixado como \(a = 0.3\). Os descendentes \(\vartheta_1'\) e \(\vartheta_2'\) são calculados como combinações lineares dos vetores progenitores:

\[\begin{align}\\ \vartheta_1' &= a \, \vartheta_1 + (1 - a) \, \vartheta_2 = 0.3 \times (0.8, 1.2, 0.5) + 0.7 \times (1.4, 0.6, 1.0) = (1.22, 0.84, 0.85) \\\\ \vartheta_2' &= (1 - a) \, \vartheta_1 + a \, \vartheta_2 = 0.7 \times (0.8, 1.2, 0.5) + 0.3 \times (1.4, 0.6, 1.0) = (1.0, 1.08, 0.65) \\\\ \end{align}\]

Esse operador é especialmente adequado para problemas de otimização numérica, pois permite a geração de soluções intermediárias entre os progenitores, enriquecendo a exploração do espaço de busca e aumentando a diversidade genética da população.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Definição 3.9 (Mutação - Sivanandam et al. 2008; Givens e Hoeting 2012; Kramer e Kramer 2017). A mutação é um procedimento que introduz variações aleatórias nos indivíduos da população, modificando o valor de um ou mais genes em um cromossomo. A mutação é essencial para manter a diversidade genética e evitar a convergência prematura a ótimos locais. Formalmente, seja \(c = (c_1, c_2, \dots, c_L)\) um cromossomo com genes \(c_i \in \Sigma\), onde \(\Sigma\) é o alfabeto utilizado (binário, inteiro, real etc.). A mutação gera um cromossomo mutado \(c' = \mathcal{M}(c)\) tal que:

\[\begin{align}\\ c_i' = \begin{cases} \mathcal{M}(c_i), & \text{com probabilidade } p_m, \\\\ c_i, & \text{com probabilidade } 1 - p_m, \end{cases}\\\\ \end{align}\]

onde \(p_m \in [0,1]\) é a taxa de mutação, um parâmetro que controla a frequência das alterações genéticas. Os tipos mais comuns de operadores de mutação incluem:

• Mutação em Codificação Binária: Neste caso, a mutação inverte o valor de um gene. Isto é, se \(c_i = 1\), então \(\mathcal{M}(c_i) = 0\) após a mutação, e vice-versa.

• Mutação Gaussiana (Codificação em Ponto Flutuante): Aplica perturbações gaussianas aos genes. Para um cromossomo \(\mathbf{x} = [x_1, \dots, x_n]\), a versão mutada é: \[\begin{align}\\ \mathbf{x}' = \mathbf{x} + \mathcal{N}(0,\sigma)\\\\ \end{align}\] onde \(\mathcal{N}(0, \sigma)\) é um vetor de variáveis aleatórias gaussianas independentes, cada uma com média zero e desvio padrão \(\sigma\).

• Mutação Uniforme: Seleciona aleatoriamente um gene \(k \in {1, 2, \dots, n}\) do cromossomo \(\mathbf{x} = [x_1, \dots, x_n]\) e o substitui por um novo valor \(x_k'\) amostrado uniformemente no domínio permitido, produzindo \(\mathbf{x}' = [x_1, \dots, x_{k}^{'}, \ldots, x_n]\).

• Mutação Não-Uniforme: É um operador dinâmico, no qual a magnitude da mutação decresce ao longo das gerações. Seja \(\mathbf{x}^{(t)} = [x_1, \dots, x_k, \dots, x_n]\) o cromossomo na geração \(t\). O gene \(x_k\) é alterado para: \[\begin{align}\\ x_k' = \begin{cases} x_k + \Delta(t, L_{\text{sup}} - x_k), & \text{com probabilidade } 0.5 \\\\ x_k - \Delta(t, x_k - L_{\text{inf}}), & \text{com probabilidade } 0.5 \end{cases}\\\\ \end{align}\] onde \(L_{\text{sup}}\) e \(L_{\text{inf}}\) são os limites superior e inferior do domínio do gene, e \(\Delta(t, y)\) é uma função decrescente em \(t\) que define a amplitude máxima de perturbação permitida na geração \(t\).

Exemplo 3.12 (Mutação Binária). Considere um problema no qual os cromossomos sejam codificados em binário. Seja o cromossomo:

\[\begin{align}\\ c = (0, 1, 1, 0, 0, 1, 0, 1, 1, 0) \in \{0,1\}^{10}\\\\ \end{align}\]

Suponha uma taxa de mutação \(p_m = 0,1\). Ao percorrer cada gene, apenas o quinto gene \(c_5\) seja escolhido para sofrer mutação. No operador binário, isso significa inverter o bit. Assim:

\[\begin{align}\\ c_5 = 0 \quad \longrightarrow \quad \mathcal{M}(c_5) = 1\\\\ \end{align}\]

O cromossomo mutado passa a ser:

\[\begin{align}\\ c' = (0, 1, 1, 0, 1, 1, 0, 1, 1, 0)\\\\ \end{align}\]

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 3.13 (Mutação Não-Uniforme). Considere um cromossomo real:

\[\begin{align}\\ \mathbf{x}^{(t)} = [3.0,\; 6.5,\; 2.2]\\\\ \end{align}\]

Suponha que o algoritmo esteja na geração \(t\), e que o domínio de cada gene seja \([0, 10]\). Escolhe-se o segundo gene \(x_2\) para sofrer mutação. O operador não-uniforme modifica esse gene conforme:

\[\begin{align}\\ x_2' = \begin{cases} x_2 + \Delta(t, 10 - x_2), & \text{com probabilidade } 0.5\\\\ x_2 - \Delta(t, x_2 - 0), & \text{com probabilidade } 0.5 \end{cases}\\\\ \end{align}\]

onde \(\Delta(t, y)\) é uma função que diminui com o tempo, por exemplo:

\[\begin{align}\\ \Delta(t, y) = y \cdot (1 - r^{(1 - t / T)})\\\\ \end{align}\]

com \(r\) sorteado uniformemente em \([0,1]\) e \(T\) sendo o número total de gerações. Seja \(r = 0.3\), \(t = 50\), \(T = 100\), então:

\[\begin{align}\\ \Delta(50, 10 - 6.5) = 3.5 \times \left(1 - 0.3^{0.5}\right) \approx 3.5 \times (1 - 0.5477) \approx 1.58\\\\ \end{align}\]

Admita que, neste caso, ocorra um incremento. Então:

\[\begin{align}\\ x_2' = 6.5 + 1.58 = 8.08\\\\ \end{align}\]

Logo, o cromossomo após mutação torna-se:

\[\begin{align}\\ \mathbf{x}' = [3.0,\; 8.08,\; 2.2]\\\\ \end{align}\]

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Para dar continuidade ao ciclo evolutivo, o algoritmo genético precisa definir quais indivíduos serão preservados ou propagados para a próxima geração. Essa etapa é denominada seleção e tem por objetivo favorecer a reprodução dos indivíduos mais adaptados, guiando o processo de busca para regiões promissoras do espaço de soluções. O principal método utilizado para este propósito é o chamado método da roleta (roulette wheel selection). Nesse esquema, cada indivíduo da população recebe uma probabilidade de seleção proporcional ao seu valor de fitness em relação à soma total dos valores de fitness de todos os indivíduos. Formalmente, se a população na geração \(t\) é composta pelos indivíduos \(\vartheta_1^{(t)}, \ldots, \vartheta_P^{(t)}\), e se \(\mathcal{F}(\theta_i^{(t)})\) representa o fitness do fenótipo correspondente \(\theta_i^{(t)}\), então a probabilidade \(p_i^{(t)}\) de seleção do indivíduo \(\vartheta_i^{(t)}\) é dada por:

\[\begin{align}\\ p_i^{(t)} = \frac{\mathcal{F}\bigl(\theta_i^{(t)}\bigr)}{\displaystyle\sum_{j=1}^{P} \mathcal{F}\bigl(\theta_j^{(t)}\bigr)}\\\\ \end{align}\]

Dessa forma, indivíduos com maior aptidão possuem maiores chances de serem escolhidos para participar das operações de crossover e mutação na formação da próxima geração. No entanto, indivíduos menos aptos ainda mantêm alguma probabilidade de seleção, o que preserva a diversidade genética e evita uma convergência prematura do algoritmo. A Figura 3.9 ilustra, de modo esquemático, um resumo da estrutura do algoritmo genético, evidenciando o fluxo iterativo das etapas de seleção, crossover e mutação que regem a evolução da população ao longo das gerações.

Figura 3.9.. Estrutura geral de um algoritmo genético, destacando o ciclo evolutivo composto pelas etapas de avaliação da função fitness, seleção dos indivíduos mais aptos, aplicação dos operadores genéticos (crossover e mutação) e formação da nova população.

Fonte: Genetic Algorithm | Givens, G. H. & Hoeting, J. A. (2013). Computational Statistics. 2\(^{\text{nd}}\) Edition. John Wiley & Sons, Inc.

Outros métodos de seleção também podem ser empregados, como a seleção por torneio (tournament selection), em que grupos aleatórios de indivíduos competem entre si e o melhor de cada grupo é selecionado; a seleção estocástica universal (stochastic universal sampling), que distribui proporcionalmente a seleção entre os indivíduos garantindo menor variância amostral; ou a seleção por classificação (ranking selection), na qual os indivíduos são ordenados segundo seu fitness e recebem probabilidades de seleção baseadas apenas nessa ordem, reduzindo o risco de convergência prematura.

Exemplo 3.14. Suponha uma população na geração \(t\) formada por quatro indivíduos, cujos fenótipos apresentam os seguintes valores de fitness:

\[\begin{align}\\ \mathcal{F}(\theta_1^{(t)}) &= 12 \\\\ \mathcal{F}(\theta_2^{(t)}) &= 8 \\\\ \mathcal{F}(\theta_3^{(t)}) &= 5 \\\\ \mathcal{F}(\theta_4^{(t)}) &= 15 \\\\ \end{align}\]

A soma total dos valores de fitness é:

\[\begin{align}\\ \sum_{i=1}^{4} \mathcal{F}\bigl(\theta_i^{(t)}\bigr) = 12 + 8 + 5 + 15 = 40\\\\ \end{align}\]

Assim, as probabilidades de seleção dos indivíduos serão:

\[\begin{align}\\ p_1^{(t)} &= \frac{12}{40} = 0.30 \\\\ p_2^{(t)} &= \frac{8}{40} = 0.20 \\\\ p_3^{(t)} &= \frac{5}{40} = 0.125 \\\\ p_4^{(t)} &= \frac{15}{40} = 0.375\\\\ \end{align}\]

Logo, o indivíduo \(\vartheta_4^{(t)}\), possuidor do maior valor de fitness, possui a maior chance (37.5%) de ser selecionado para compor a próxima geração. Entretanto, todos os demais indivíduos mantêm alguma probabilidade de serem escolhidos, o que preserva a diversidade genética e possibilita a exploração de regiões distintas do espaço de busca.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.3.4.3. Implementação Computacional

---

Nesta seção, abordaremos a implementação computacional completa de um algoritmo genético, contemplando todos os seus componentes fundamentais: inicialização da população, avaliação da função fitness, seleção, aplicação dos operadores genéticos de crossover e mutação, e formação das gerações subsequentes. Para ilustrar esse processo de forma prática e didática, utilizaremos como base o problema apresentado no Exemplo 3.15, aplicando passo a passo cada etapa do algoritmo.

Exemplo 3.15 (Pedras Preciosas em um Universo de RPG - Blørstad 2019; Pêcheux 2020). Em um universo de RPG, considere seis pedras preciosas distintas (Figura 3.10), cada uma com peso \(w_i\) (em gramas) e valor \(v_i\) (em moedas de ouro) que representam, respectivamente, sua raridade e seu poder ou preço no mercado de aventureiros. O desafio consiste em escolher um conjunto dessas pedras preciosas de modo que o peso total transportado pelo personagem não ultrapasse 16 gramas, enquanto se busca maximizar o valor total obtido, seja para uso em magias poderosas ou para troca por ouro e suprimentos. Este problema constitui, em particular, uma adaptação do clássico problema da mochila (Knapsack Problem), amplamente estudado em otimização combinatória (Pisinger e Toth 1998; Martello e Toth 1990; Bellman 1957; Dantzig 1957).

Figura 3.10.. Representação das seis pedras preciosas do problema Knapsack. Cada pedra preciosa é identificada por uma cor específica, apresentando valores distintos de peso (W - em gramas) e valor (V - em moedas de ouro), compondo o cenário do problema adaptado da mochila.

Fonte: A Peek at Genetic Algorithms | Medium - Mina Pêcheux. Acesso: Junho, 2025.

Para modelar esse problema usando o algoritmo genético, considere que cada pedra preciosa é representada por um gene binário, onde o valor 1 indica que a pedra foi incluída na mochila (inventário) do personagem, e o valor 0 indica que ela foi deixada para trás na masmorra. Assim, o cromossomo corresponde a uma cadeia binária que codifica uma solução candidata para o problema de otimização proposto. Os dados relativos às seis pedras preciosas deste cenário — incluindo seus pesos e valores correspondentes — estão apresentados no Código 3.6.

Código 3.6. Dados do problema de seleção de pedras preciosas, contendo os vetores de pesos e valores associados a cada pedra, que serão utilizados na implementação do algoritmo genético.

# -----------------------------------------------
# Dados do Problema: Seleção de Pedras Preciosas 
# -----------------------------------------------

# --- 0. Pacotes necessários ---

library(rlist)
library(knitr)

# --- 1. Dados ---

set.seed(2025)

Gems <- data.frame(Color = c("Green Triangle", "Yellow Square", "Gray Pear", 
                             "Blue Pentagon", "Red Hexagon", "Purple Octogon"),
                   Weight = c(1, 4, 5, 6, 8, 25),
                   Value = c(3, 4, 6, 8, 15, 20))

# --- 2. Exibição dos dados ---

kable(Gems, 
      format = 'html', 
      align = 'c', 
      escape = FALSE,
      col.names = c('Cor da Pedra <br> Preciosa', 'Peso (w) <br> (em g)', 'Valor (v) <br> (em moedas de ouro)'))

Cor da Pedra Preciosa	Peso (w) (em g)	Valor (v) (em moedas de ouro)
Green Triangle	1	3
Yellow Square	4	4
Gray Pear	5	6
Blue Pentagon	6	8
Red Hexagon	8	15
Purple Octogon	25	20

No processo de implementação do algoritmo genético, a etapa inicial consiste na geração de uma população inicial de soluções candidatas. Cada indivíduo dessa população é representado por um vetor binário, onde cada posição indica a decisão estratégica de incluir (1) ou não incluir (0) determinada pedra preciosa no inventário do personagem. O procedimento de inicialização da população, apresentado no Código 3.7, consiste em gerar aleatoriamente uma matriz binária, em que as linhas correspondem aos diferentes indivíduos (possíveis combinações de pedras, isto é, inventários) e as colunas representam as pedras preciosas disponíveis para seleção.

Código 3.7. Função de inicialização da população para o problema de seleção de pedras preciosas.

# ---------------------------------------
# Função Para Inicialização da População 
# ---------------------------------------

InitializePop   <- function(size) 
{
  pop           <- t(sapply(1:size, function(i) round(runif(nrow(Gems), 0, 1), 0)))
  colnames(pop) <- Gems$Color
  return(pop)
}

Após a inicialização da população, torna-se necessária a avaliação da aptidão (fitness) de cada indivíduo gerado. No nosso problema, a função fitness calcula o valor total acumulado das pedras preciosas selecionadas por cada indivíduo, garantindo que o peso total não ultrapasse o limite máximo permitido para o inventário do personagem. Caso a soma dos pesos exceda esse limite, a solução é penalizada com fitness igual a zero, sendo considerada inviável para utilização no jogo. O procedimento dessa avaliação encontra-se detalhado no Código 3.8.

Código 3.8. Função fitness para o problema de seleção de pedras preciosas.

# ---------------
# Função Fitness
# ---------------

fitness_function  <- function(population, weightRestriction)
{
  score           <- numeric(nrow(population))
  
  for (i in 1:nrow(population))
  {
    temp          <- apply(Gems[which(population[i, 1:nrow(Gems)] == 1), 2:3], 2, sum)
    score[i]      <- ifelse(temp[1] > weightRestriction, 0, temp[2])
  }
  
  pop             <- cbind(population, score)
  pop             <- pop[order(score, decreasing = TRUE), ]
  
  return(pop)
}

Após a avaliação da fitness dos indivíduos, aplica-se o operador de crossover, responsável pela recombinação genética entre pares selecionados da população. Para o nosso caso, essa operação combina as estratégias genéticas de dois indivíduos para gerar um novo indivíduos (inventário de pedras preciosas), potencialmente mais valioso. O crossover adotado é do tipo um ponto, no qual a primeira metade da seleção de pedras do primeiro indivíduo é combinada com a segunda metade da seleção do segundo indivíduo, originando um novo conjunto de pedras. O procedimento correspondente encontra-se detalhado no Código 3.9.

Código 3.9. Operação de crossover para o problema de seleção de pedras preciosas.

# --------------------
# Função de Crossover
# --------------------

crossover     <- function(dna1, dna2) 
{
  len         <- 1:ceiling(length(dna1) / 2)
  offspring   <- c(dna1[len], dna2[-len])
  return(offspring)
}

Após a aplicação do crossover, realiza-se a operação de mutação, a qual introduz variações aleatórias no material genético dos descendentes. No nosso contexto, essa mutação corresponde a alterar a decisão de incluir ou não uma determinada pedra preciosa no inventário do personagem, invertendo o valor binário na posição selecionada — ou seja, uma pedra que estava incluída (1) pode ser removida (0), ou vice-versa. O procedimento detalhado dessa operação está implementado no Código 3.10.

Código 3.10. Operação de mutação para o problema de seleção de pedras preciosas.

# ------------------
# Função de Mutação
# ------------------

mutate      <- function(dna)
{
  ind       <- sample(1:length(dna), 1)
  dna[ind]  <- ifelse(dna[ind] == 1, 0, 1)
  return(dna)
}

Chegamos agora ao procedimento de combinação genética ou breeding, responsável por gerar novos indivíduos a partir de estratégias genéticas existentes na população. A ideia é produzir novos inventários de pedras preciosas por meio da recombinação genética, preservando características vantajosas. Para manter a diversidade genética e evitar a convergência prematura, a combinação genética incorpora a possibilidade de mutação, ainda que com baixa probabilidade. No nosso caso, estabelece-se uma probabilidade de 95% para aplicação do crossover e 5% para ocorrência de mutação direta em um inventário. O procedimento está implementado no Código 3.11.

Código 3.11. Combinação genética (breeding) para o problema de seleção de pedras preciosas.

# --------------------
# Combinação Genética
# --------------------

Breed         <- function(parent1, parent2, population) 
{
  dna1        <- population[parent1, -ncol(population)]
  dna2        <- population[parent2, -ncol(population)]
  
  if (runif(1, 0, 100) > 5) {
    offspring <- crossover(dna1, dna2) # 95% de chance de crossover
  } else {
    offspring <- mutate(dna1)          # 5% de chance de mutação
  }
  
  return(offspring)
}

Por fim, chegamos à etapa de criação da nova geração, responsável por formar o conjunto de inventários que irão compor a próxima iteração do algoritmo genético. A função apresentada a seguir implementa essa lógica combinando três estratégias fundamentais:

• Seleção: a população atual é ordenada com base na aptidão (fitness), de modo que as combinações de pedras mais valiosas e viáveis tenham maior probabilidade de serem mantidas para a próxima geração.

• Cruzamentos: pares de inventários bem classificados são combinados sistematicamente para gerar novos conjuntos de pedras preciosas, assegurando a transmissão de características vantajosas entre as soluções.

• Introdução de Diversidade: além dos descendentes gerados, novos inventários aleatórios são introduzidos na população, garantindo variabilidade genética e evitando a convergência prematura em soluções subótimas.

O procedimento está implementado no Código 3.12.

Código 3.12. Criação da nova geração para o problema de seleção de pedras preciosas.

# ------------------------
# Criação da Nova Geração
# ------------------------

BreedPopulation   <- function(population) 
{
  # Avaliação da população atual

  population      <- fitness_function(population, 16)
  
  NewPopulation   <- list()
  len             <- floor(nrow(population) / 4)
  
  for (i in 1:(len - 1)) 
  {
    NewPopulation <- list.append(NewPopulation, Breed(i, i + 1, population))
    NewPopulation <- list.append(NewPopulation, Breed(i + 1, i, population))
    NewPopulation <- list.append(NewPopulation, Breed((len + 1 - i), i, population))
    NewPopulation <- list.append(NewPopulation, Breed(i, (len + 1 - i), population))
  }
  
  # Inclusão de novos indivíduos aleatórios para manter diversidade genética
  
  NewPopulation   <- list.append(NewPopulation, InitializePop(1))
  NewPopulation   <- list.append(NewPopulation, InitializePop(1))
  NewPopulation   <- list.append(NewPopulation, InitializePop(1))
  NewPopulation   <- list.append(NewPopulation, InitializePop(1))
  
  # Concatenação da nova população em uma matriz
  
  pop             <- Reduce(rbind, NewPopulation)
  
  return(pop)
}

Agora que todas as funções necessárias do algoritmo genético foram construídas, podemos implementar o procedimento final. A função Genetic_Algorithm executa todo o fluxo do algoritmo, desde a inicialização da população de inventários até a evolução das combinações de pedras preciosas ao longo das gerações, acompanhando a performance das melhores seleções obtidas. Essa função opera da seguinte forma:

• Inicialização: é criada uma população inicial de inventários, com tamanho definido por popsize.

• Iterações Evolutivas: o algoritmo realiza, no máximo, Gen gerações, produzindo a cada ciclo uma nova população por meio dos operadores genéticos de crossover, mutação e seleção, que combinam e modificam os inventários para buscar combinações mais valiosas.

• Monitoramento de Desempenho: em cada geração, são registrados o valor médio e o valor máximo do fitness da população, que refletem, respectivamente, a qualidade média dos inventários e a melhor combinação de pedras preciosas encontrada até o momento. Esses valores são armazenados na matriz Performance e podem ser visualizados dinamicamente em um gráfico que acompanha o progresso do algoritmo.

• Critério de Parada: além do limite máximo de gerações, o algoritmo pode ser interrompido caso a melhor solução não apresente melhoria por um número consecutivo de gerações especificado por BreakCond.

• Retorno dos Resultados: ao final, a função retorna o histórico de desempenho por geração e os detalhes da melhor solução encontrada, incluindo os genes (pedras) selecionados, o peso total e o valor total acumulado das pedras preciosas.

O procedimento está implementado no Código 3.13.

Código 3.13. Função final do algoritmo genético para o problema de seleção de pedras preciosas.

# -----------------------------------
# Função Final do Algoritmo Genético
# -----------------------------------

Genetic_Algorithm   <- function(popsize, Gen, BreakCond)
{
  # Inicialização da população
  
  population        <- InitializePop(popsize)
  
  # Data.frame para armazenar o desempenho por geração
  
  Performance       <- data.frame(Generation = integer(),
                                  Average = numeric(),
                                  Best = numeric(),
                                  stringsAsFactors = FALSE)
  
  breakCond         <- 0
  
  for (j in 1:Gen) 
  {
    # Geração de nova população
    
    population      <- BreedPopulation(population)
    
    # Avaliação de fitness
    
    score           <- fitness_function(population, 16)[, nrow(Gems) + 1]
    
    # Armazena valores da geração
    
    Performance     <- rbind(Performance,
                             data.frame(Generation = j,
                                        Average = mean(score),
                                        Best = max(score)))
    
    # Verificação do critério de parada
    
    if (Performance$Best[j] > max(Performance$Best[1:(j - 1)], na.rm = TRUE)) {
      breakCond     <- 0
    } else {
      breakCond     <- breakCond + 1
    }
    
    if (breakCond >= BreakCond) {
      break
    }
  }
  
  # Visualização gráfica
  
  Performance_long  <- reshape2::melt(Performance,
                                      id.vars = "Generation",
                                      variable.name = "Metric",
                                      value.name = "Value")
  
  p                 <- ggplot(Performance_long, aes(x = Generation, y = Value, color = Metric)) +
                        geom_line(size = 1) +
                        geom_point(size = 2) +
                        labs(title = "",
                             y = "Fitness Value",
                             color = "Metric") +
                        theme_minimal() + 
                        theme(axis.title.x = element_text(margin = margin(t = 10)),
                              axis.title.y = element_text(margin = margin(r = 10)))
  print(p)
  
  # Cálculo do resultado final
  
  result            <- list(GenerationPerformance = Performance,
                            Solution = list(Genes = population[1, ],
                                            Result = apply(Gems[which(population[1, ] == 1), 2:3], 2, sum)))
  return(result)
}

Para finalizar, o algoritmo genético foi aplicado ao problema das pedras preciosas, configurado com uma população inicial de 20 inventários, um máximo de 200 gerações e um critério de parada que interrompe a execução caso a melhor solução não se modifique após 150 gerações consecutivas. O procedimento encontra-se implementado no Código 3.14.

Código 3.14. Aplicação do algoritmo genético para o problema de seleção de pedras preciosas.

# --------------------------------
# Aplicação do Algoritmo Genético
# --------------------------------

GA_Result <- Genetic_Algorithm(popsize = 20, Gen = 200, BreakCond = 150)

A Figura 3.11 ilustra o gráfico gerado pelo algoritmo genético, mostrando a evolução da média de fitness da população (Average) e o valor do melhor indivíduo (Best) ao longo das gerações. Observa-se, assim, o comportamento do processo evolutivo, permitindo identificar a convergência do algoritmo para regiões promissoras do espaço de busca.

Figura 3.11.. Evolução do algoritmo genético para o problema das pedras preciosas. As linhas mostram o valor médio de fitness da população (em vermelho) e o valor do melhor indivíduo (em azul) ao longo das gerações.

A melhor solução encontrada, exibida pelo objeto GA_Result$Solution, apresenta a combinação de pedras preciosas selecionadas, juntamente com o peso total e o valor total associado a essa seleção, conforme descrito no Código 3.15.

Código 3.15. Resultados da aplicação do algoritmo genético para o problema da seleção ótima de pedras preciosas.

# ---------------------------------------------------------------------
# Resultados do Algoritmo Genético Para o Problema de Pedras Preciosas
# ---------------------------------------------------------------------

GA_Result$Solution

$Genes
Green Triangle  Yellow Square      Gray Pear  Blue Pentagon    Red Hexagon 
             1              0              1              0              1 
Purple Octogon 
             0 

$Result
Weight  Value 
    14     24 

3.4. Otimizadores no Ambiente R

Agora que conhecemos alguns dos métodos numéricos, podemos passar ao próximo conceito de interesse: os otimizadores. Um otimizador consiste em um algoritmo ou procedimento computacional projetado para localizar, de maneira eficiente, o ponto ótimo da função objetivo, observando as restrições impostas ao problema. Na prática, existe uma variedade de otimizadores empregados na resolução de problemas de otimização, abrangendo desde métodos determinísticos clássicos — como algoritmos baseados em gradientes (BFGS, Newton-Raphson, gradiente conjugado), métodos de programação quadrática e programação sequencial — até abordagens estocásticas e metaheurísticas voltadas para otimização global, tais como algoritmos genéticos, recozimento simulado (simulated annealing), otimização por enxame de partículas (particle swarm optimization) e evolução diferencial (differential evolution). Nas seções seguintes, apresentam-se os principais otimizadores disponíveis no ambiente R, destacando-se pela aplicabilidade em distintas classes de problemas e por oferecerem interfaces robustas para a formulação e solução de problemas de otimização.

3.4.1. optim

---

O optim é um dos otimizadores mais tradicionais e amplamente utilizados no ambiente R. Este otimizador oferece uma interface flexível e eficiente para a minimização (ou maximização) de funções objetivo, permitindo a inclusão de restrições simples, como limites inferiores e superiores para os parâmetros. Além disso, o optim disponibiliza diversos métodos numéricos para a otimização, abrangendo tanto algoritmos baseados em gradiente — como BFGS (Broyden-Fletcher-Goldfarb-Shanno) e gradiente conjugado — quanto técnicas que não requerem informações derivadas, como o método Nelder-Mead e o SANN (Simulated Annealing), também conhecido como recozimento simulado.

Otimizador 3.1 (optim). O optim é um otimizador nativo da linguagem R que oferece um conjunto abrangente de métodos numéricos para a resolução de problemas de otimização geral, com ênfase em funções diferenciáveis e não-lineares. Os principais argumentos desse otimizador são:

• par: é o vetor que contém os valores iniciais (ou chute inicial) do processo de otimização.
• fn: é a função que deve ser maximizada (ou minimizada). A função deve ter a forma fn(x, ...), onde x é um vetor que contém os valores das variáveis que a função fn tenta maximizar (ou minimizar) e ... são os argumentos adicionais da função fn.
• gr: é a função gradiente, isto é, a derivada da função fn em relação a x. Se gr não for fornecido, a função optim() usará uma estimativa numérica do gradiente.
• method: é o método de otimização a ser usado. Existem vários métodos disponíveis, como BFGS, L-BFGS-B, CG, Nelder-Mead, SANN, entre outros. O método padrão é Nelder-Mead.
• control: é uma lista de argumentos de controle para o método de otimização. Neste caso, para maximização, utiliza-se control = list(fnscale = -1), e, para minimização, utiliza-se control = list(fnscale = 1).

Exemplo 3.16. Considere a função escalar definida por:

\[\begin{align}\\ f(x_1, x_2) = 2(x_1 - 1)^2 + 5(x_2 - 3)^2 + 10\\\\ \end{align}\]

cujo domínio é \(\mathbb{R}^2\) e cujos coeficientes foram definidos de forma a simular diferentes curvaturas nas direções das variáveis \(x_1\) e \(x_2\). Nosso objetivo consiste em encontrar o ponto de mínimo desta função, ou seja, o vetor \(\mathbf{x} = (x_1, x_2)^\top \in \mathbb{R}^2\) que minimiza \(f(\mathbf{x})\). Para tal, será considerado o otimizador optim, em que será implementada uma função f que define a função objetivo a ser minimizada. O ponto inicial para o processo iterativo de otimização será estabelecido pelo vetor par = c(1, 1), a partir do método numérico padrão do optim (isto é, Nelder-Mead). Além disso, será especificado fnscale = 1 no argumento control para indicar explicitamente que se trata de um problema de minimização. Por fim, o objeto resultado irá armazenar a saída do otimizador, contendo informações relevantes sobre o processo, tais como:

• par: o ponto ótimo estimado (isto é, a solução numérica do problema);
• value: o valor da função objetivo no ponto ótimo encontrado;
• convergence: um indicador da convergência do algoritmo. No caso do optim, convergence = 0 denota que a otimização foi bem-sucedida.

O Código 3.16 apresenta uma rotina em ambiente R que exemplifica a aplicação do otimizador optim para a resolução do problema de minimização da função \(f(x_1, x_2) = 2(x_1 - 1)^2 + 5(x_2 - 3)^2 + 10\). Já a Figura 3.12 exibe o mapa de contorno da função com a marcação da solução ótima obtida.

Código 3.16. Implementação em R da solução de um problema de minimização da função \(f(x_1, x_2) = 2(x_1 - 1)^2 + 5(x_2 - 3)^2 + 10\) utilizando o otimizador optim.

# ------------------------------------------------------------
# Otimizador `optim`: Minimização da Função Objetivo f(x1,x2)
# ------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Definição da função objetivo ---

f         <- function(x) 
{
  2 * (x[1] - 1)^2 + 5 * (x[2] - 3)^2 + 10
}

# --- 2. Aplicação do otimizador ---

res       <- optim(par = c(1, 1), fn = f, control = list(fnscale = 1))
print(res)

$par
[1] 1.000168 3.000232

$value
[1] 10

$counts
function gradient 
      75       NA 

$convergence
[1] 0

$message
NULL

# --- 3. Visualização gráfica da solução ótima ---

# Definição do ponto ótimo

ponto_otimo   <- data.frame(x = res$par[1], y = res$par[2], grupo = "Solução Ótima")

# Construção da grade para o mapa de contorno

x_seq         <- seq(-2, 4, length.out = 150)
y_seq         <- seq(0, 6, length.out = 150)
grid          <- expand.grid(x = x_seq, y = y_seq)
grid$z        <- apply(grid, 1, function(row) f(c(row[1], row[2])))

# Preparação dos dados para o gráfico

df_surface    <- data.frame(x = grid$x, y = grid$y, z = grid$z, grupo = "Superfície da Função")
df_otimo      <- data.frame(x = res$par[1], y = res$par[2], grupo = "Solução Ótima")

# Gráfico

ggplot() +
  geom_raster(data = df_surface, aes(x = x, y = y, fill = z),
              interpolate = TRUE, alpha = 0.85) +
  geom_contour(data = df_surface, aes(x = x, y = y, z = z),
               color = "gray30", size = 0.3) +
  geom_point(data = df_otimo, aes(x = x, y = y, color = grupo), size = 3) +
  scale_fill_gradientn(name = "Valor da\nFunção Objetivo",
                       colours = rainbow(7, s = 0.5, v = 0.9)) +
  scale_color_manual(name = "Legenda",
                     values = c("Solução Ótima" = "red")) +
  labs(title = "",
       x = expression(x[1]),
       y = expression(x[2])) +
  theme_minimal() +
  theme(legend.position = "right",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.12.. Mapa de contorno da função objetivo \(f\), com indicação do ponto da solução ótima obtida pelo otimizador optim (marcado em vermelho).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.2. lp

---

O lp, disponível no ambiente R por meio do pacote lpSolve é um otimizador projetado para lidar com funções objetivo lineares, sujeitas a restrições também lineares, abrangendo tanto problemas de minimização quanto de maximização. Esse otimizador, em particular, permite especificar coeficientes da função objetivo, matriz de restrições, sentidos das desigualdades e limites das variáveis, oferecendo grande flexibilidade para a modelagem de problemas em diversas áreas, como logística, finanças, engenharia e pesquisa operacional. Além disso, possibilita a resolução de problemas de programação inteira, nos quais algumas ou todas as variáveis de decisão são restritas a valores inteiros, ampliando ainda mais seu campo de aplicação.

Otimizador 3.2 (lp). O lp é um otimizador pertencente ao pacote lpSolve, voltado à resolução de problemas de programação linear e inteira. Proporciona uma interface estruturada para a formulação de modelos lineares de maximização ou minimização sujeitos a restrições, utilizando algoritmos clássicos, tais como o método Simplex e o Branch-and-Bound. Os principais argumentos desse otimizador são:

• direction: string de caracteres que indica a direção da otimização: min (padrão) ou max.
• objective.in: vetor numérico de coeficientes da função objetivo.
• const.mat: matriz de coeficientes de restrição numéricos, uma linha por restrição e uma coluna por variável (a menos que transpose.constraints = FALSE).
• const.dir: vetor de strings de caracteres que indicam a direção da restrição: cada valor deve ser um dos <, <=, =, ==, >, ou >=.
• const.rhs: vetor de valores numéricos para os lados direitos das restrições.
• transpose.constraints: por padrão, cada restrição ocupa uma linha de const.mat, e essa matriz precisa ser transposta antes de ser passada para o código de otimização. Para matrizes de restrição muito grandes, pode ser mais sábio construir as restrições em uma matriz coluna por coluna. Nesse caso, defina transpose.constraintscomo FALSE.

Exemplo 3.17. A otimização linear é uma técnica de otimização matemática que busca maximizar ou minimizar uma função linear sujeita a um conjunto de restrições lineares. A função objetivo e as restrições, neste caso, são expressas em termos de equações e desigualdades lineares. Por exemplo, considere um sistema linear, cujo objetivo é minimização, descrito na seguinte forma:

\[\begin{align}\\ \min_\mathbf{x}(c^T \mathbf{x}) = \min_\mathbf{x}(c_1 x_1 + ... + c_n x_n)\\\\ \end{align}\]

sujeito as restrições: \(Ax \geqslant b\), \(x \geqslant 0\). Por questões de praticidade, podemos reescrever esse sistema linear, em forma matricial, como:

\[\begin{align}\\ \min_\mathbf{x} \left[\begin{matrix} c_1 \\ c_2 \\ \ldots \\ c_n\end{matrix}\right]^T\left[\begin{matrix} x_1 \\ x_2 \\ \ldots \\ x_n\end{matrix}\right]\\\\ \end{align}\]

sujeito a:

\[\begin{align}\\ \left[\begin{matrix} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n}\\ \ldots & \ldots & \ldots & \ldots \\ a_{m1} & a_{m2} & \ldots a_{mn}\end{matrix}\right]\left[\begin{matrix} x_1 \\ x_2 \\ \ldots \\ x_n\end{matrix}\right]\geqslant \left[\begin{matrix} b_1 \\ b_2 \\ \ldots \\ b_n\end{matrix}\right]; \left[\begin{matrix} x_1 \\ x_2 \\ \ldots \\ x_n\end{matrix}\right] \geqslant 0\\\\ \end{align}\]

Um sistema linear escrita dessa forma é chamado de problema de programação linear (PL). Para compreender o funcionamento do otimizador linear na resolução de um problema de PL, considere a seguinte situação-problema:

• Função Objetivo:
  • O objetivo é maximizar o lucro total.
  • Os produtos A e B são vendidos \(R\$25\) e \(R\$20\), respectivamente.
• Restrições de Recursos:
  • O produto A requer 20 unidades de recursos, o produto B precisa de 12.
  • Apenas 1800 unidades de recursos estão disponíveis por dia.
• Restrições de Tempo:
  • Ambos os produtos requerem um tempo de produção de 1/15 hora.
  • Um dia de trabalho tem um total de 8 horas.

Considerando o objetivo do problema de PL apresentado, a formulação matemática da função objetivo, que busca maximizar o lucro total, pode ser expressa da seguinte forma:

\[\begin{align}\\ \text{Vendas}_{\max} = \max_{x_1,x_2} \ [25 x_1 + 20 x_2] = \max_{x_1,x_2} \left[\begin{matrix}25 \\ 20\end{matrix}\right]^T\left[\begin{matrix}x_1 \\ x_2\end{matrix}\right]\\\\ \end{align}\]

sujeita as restrições:

\[\begin{align}\\ 20 x_1 + 12 x_2 &\leqslant 1800\\\\ \dfrac{1}{15} x_1 + \dfrac{1}{15} x_2 &\leqslant 8\\\\ \end{align}\]

em que \(x_1\) representa a quantidade produzida do produto A e \(x_2\) representa a quantidade produzida do produto B. O Código 3.17 apresenta uma rotina em ambiente R que exemplifica a aplicação do otimizador lp do pacote lpSolve para a resolução do problema de PL proposto. Já a Figura 3.13 exibe o mapa de contorno da função com a marcação da solução ótima obtida.

Código 3.17. Implementação em R da solução de um problema de programação linear utilizando o otimizador lp do pacote lpSolve.

# --------------------------------------------------------
# Otimizador `lp` do pacote `lpSolve`: Programação Linear
# --------------------------------------------------------

# --- 0. Pacotes necessários ---

library(lpSolve)
library(ggplot2)
library(dplyr)
library(knitr)

# --- 1. Definição da função objetivo (coeficientes) ---

coef_objetivo       <- c(25, 20)

# --- 2. Definição da matriz de restrições ---

matriz_restricoes   <- matrix(c(20, 12,
                                1/15, 1/15),
                              nrow = 2, byrow = TRUE)

vetor_restricoes    <- c(1800, 8)
direcao_restricoes  <- c("<=", "<=")

# --- 3. Aplicação do otimizador `lp` ---

solucao_lp          <- lp(direction = "max",
                          objective.in = coef_objetivo,
                          const.mat = matriz_restricoes,
                          const.dir = direcao_restricoes,
                          const.rhs = vetor_restricoes)

# --- 4. Exibição dos resultados ---

list("Solução Ótima" = solucao_lp$solution,
     "Valor da Função Objetivo na Solução" = solucao_lp$objval,
     "Status de Convergência" = solucao_lp$status)

$`Solução Ótima`
[1] 45 75

$`Valor da Função Objetivo na Solução`
[1] 2625

$`Status de Convergência`
[1] 0

# --- 5. Preparação dos dados para o gráfico ---

# Solução ótima

x1_opt              <- solucao_lp$solution[1]
x2_opt              <- solucao_lp$solution[2]

# Função objetivo

f_obj               <- function(x1, x2) 
{
  25 * x1 + 20 * x2
}

# Geração do grid para visualização

x1_vals             <- seq(0, 100, length.out = 300)
x2_vals             <- seq(0, 100, length.out = 300)

grid                <- expand.grid(x1 = x1_vals, x2 = x2_vals)
grid$restr1         <- 20 * grid$x1 + 12 * grid$x2
grid$restr2         <- (1/15) * (grid$x1 + grid$x2)
grid$viavel         <- with(grid, restr1 <= 1800 & restr2 <= 8)
grid$obj            <- with(grid, f_obj(x1, x2))
grid_viavel         <- grid %>% filter(viavel)

# --- 6. Geração do gráfico ---

ggplot() +
  geom_raster(data = grid_viavel, aes(x = x1, y = x2, fill = obj),
              interpolate = TRUE, alpha = 0.7) +
  geom_contour(data = grid_viavel, aes(x = x1, y = x2, z = obj),
               bins = 20,
               color = "black",
               size = 0.3) +
  geom_point(aes(x = x1_opt, y = x2_opt), color = "red", size = 3) +
  annotate("text", x = x1_opt + 3, y = x2_opt + 3,
           label = paste0("Solução Ótima: (",
                          round(x1_opt, 2), ", ",
                          round(x2_opt, 2), ")"),
           hjust = 0, size = 4) +
  scale_fill_gradientn(name = "Valor da\nFunção Objetivo",
                       colours = rainbow(7, s = 0.5, v = 0.9)) +
  labs(title = "",
       x = expression(x[1]),
       y = expression(x[2])) +
  theme_minimal() +
  theme(legend.position = "right",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.13.. Mapa de contorno da função objetivo com a região viável determinada pelas restrições lineares do problema. O ponto em vermelho indica a solução ótima obtida pelo otimizador lp.

Portanto, a solução ótima do problema de PL determina que, para maximizar o lucro, devem ser produzidas 45 unidades do produto A e 75 unidades do produto B, ou seja, \(x_1 = 45\) e \(x_2 = 75\). Essa alocação de produção, respeitando integralmente as restrições de recursos e tempo estabelecidas, resulta em um lucro total máximo de \(R\$2625\). Vale destacar que o status de retorno igual a 0 fornecido pelo otimizador lp indica que o algoritmo convergiu para uma solução ótima viável.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.3. solve.QP

---

O otimizador solve.QP, disponível no ambiente R através do pacote quadprog, é uma ferramenta especializada na resolução de problemas de programação quadrática. Esse otimizador, em particular, requer a especificação da matriz de coeficientes da parte quadrática da função objetivo, o vetor de coeficientes lineares, bem como a matriz e o vetor que definem as restrições do problema. Além de resolver problemas com restrições, o algoritmo é capaz de tratar casos em que a matriz de covariância (ou matriz Hessiana) é positiva definida, condição necessária para garantir a convexidade do problema e, consequentemente, a unicidade da solução ótima.

Otimizador 3.3 (solve.QP). Otimizador do pacote quadprog dedicado à resolução de problemas de programação quadrática. Este otimizador, em particular, emprega algoritmos baseados em métodos de pontos interiores e técnicas de ativação de restrições ativas para minimizar uma função objetivo sujeita a restrições lineares de igualdade e desigualdade. Os principais argumentos desse otimizador são:

• Dmat: matriz \(D\), que deve ser simétrica e definida positiva. Define os termos quadráticos da função objetivo.
• dvec: vetor \(\mathbf{d}\), com os coeficientes lineares da função objetivo.
• Amat: matriz de restrições, cujas colunas representam os vetores de coeficientes das desigualdades. Isto é, cada coluna define uma restrição da forma \(\mathbf{a}_j^\top \mathbf{x} \geqslant b_j\).
• bvec: vetor com os lados direitos das restrições, ou seja, o vetor \(\mathbf{b}\).
• meq: número de restrições de igualdade. Definir meq = 0 indica que todas as restrições são do tipo “maior ou igual”.

Exemplo 3.18. Considere o problema de otimização quadrática formulado como:

\[\begin{align}\\ \min_{\mathbf{x} \in \mathbb{R}^2} \quad & f(\mathbf{x}) = \frac{1}{2} \mathbf{x}^\top D \mathbf{x} - \mathbf{d}^\top \mathbf{x}, \\\\ \text{sujeito a} \quad & A^\top \mathbf{x} \geqslant \mathbf{b}\\\\ \end{align}\]

em que \(D\) é uma matriz simétrica definida positiva que representa os coeficientes quadráticos da função objetivo; \(\mathbf{d}\) é um vetor que define os coeficientes lineares da função objetivo; \(A^\top \mathbf{x} \geqslant \mathbf{b}\) impõe as restrições lineares do problema. Este tipo de estrutura é clássica em problemas de programação quadrática (PQ), com aplicações frequentes em otimização com regularização, economia e aprendizado de máquina. Para compreender o funcionamento do otimizador quadrático na resolução de um problema de PQ, considere a seguinte situação problema:

• Função Objetivo:

  • O objetivo é minimizar o custo total de operação, expresso em reais, de uma unidade produtiva composta por dois processos industriais interdependentes, cujos níveis de operação são representados pelas variáveis de decisão \(x_1\) e \(x_2\). A função de custo incorpora termos quadráticos e de interação entre os processos, bem como termos lineares associados a cada unidade de produção. A função objetivo a ser minimizada, representando o custo operacional total em reais, é dada por: \[\begin{align}\\ f(x_1, x_2) = 20x_1^2 + 15x_1x_2 + 20x_2^2 - 3350x_1 - 4350x_2\\\\ \end{align}\]

• Restrições de Capacidade:

  • A capacidade operacional da planta é limitada por uma combinação de recursos utilizados por ambos os processos. A primeira restrição impõe um nível mínimo de atividade global, sendo definida como \(10x_1 + 15x_2 \geqslant 1850\). Já a segunda restrição impõe um limite superior à soma dos esforços produtivos, sendo definida por \(x_1 + x_2 \leqslant 140\). Além disso, as variáveis de decisão devem ser não negativas, isto é, \(x_1 \geqslant 0\) e \(x_2 \geqslant 0\).

Considerando o objetivo do problema de PQ apresentado, a formulação matemática da função objetivo, que busca minimizar o custo operacional, pode ser expressa da seguinte forma:

\[\begin{align}\\ \min_{x_1, x_2} \quad f(x_1, x_2) = 20x_1^2 + 15x_1x_2 + 20x_2^2 - 3350x_1 - 4350x_2\\\\ \end{align}\]

sujeita as restrições:

\[\begin{align}\\ 10x_1 + 15x_2 &\geqslant 1850 \\\\ x_1 + x_2 &\leqslant 140 \\\\ x_1, x_2 &\geqslant 0 \\\\ \end{align}\]

em que \(x_1\) representa o nível de produção do primeiro processo industrial e \(x_2\) o nível de produção do segundo processo industrial. O Código 3.18 apresenta uma rotina em ambiente R que exemplifica a aplicação do otimizador solve.QP do pacote quadprog para a resolução do problema de PQ proposto. Já a Figura 3.14 exibe o mapa de contorno da função com a marcação da solução ótima obtida.

Código 3.18. Implementação em R da solução de um problema de programação quadrática utilizando o otimizador solve.QP do pacote quadprog.

# -------------------------------------------------------------------
# Otimizador `solve.QP` do pacote `quadprog`: Programação Quadrática
# -------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(quadprog)
library(ggplot2)
library(dplyr)

# --- 1. Definição da função objetivo ---

# Matriz D da parte quadrática (simétrica e positiva definida)

Dmat <- matrix(c(40, 15,
                 15, 40),
               nrow = 2, byrow = TRUE)

# Vetor d (note que solve.QP resolve 1/2 x' D x - d' x)

dvec <- c(-3350, -4350)

# --- 2. Definição das restrições ---

# Matriz de restrições Amat (cada coluna é uma restrição)

Amat <- matrix(c(10,    15,     # Primeira coluna (restrição 1)
                 -1,    -1,     # Segunda coluna (restrição 2)
                  1,     0,     # Terceira coluna (restrição 3)
                  0,     1),    # Quarta coluna (restrição 4)
               nrow = 2, byrow = FALSE)

# Vetor b das restrições

bvec <- c(1850, -140, 0, 0)

# --- 3. Aplicação do otimizador `solve.QP` ---

solucao_qp    <- solve.QP(Dmat, dvec, Amat, bvec, meq = 0)

# --- 4. Exibição dos resultados ---

list('Solução Ótima' = c(round(solucao_qp$solution[1]), round(solucao_qp$solution[2])),
     "Valor da Função Objetivo na Solução" = solucao_qp$value)

$`Solução Ótima`
[1]  26 106

$`Valor da Função Objetivo na Solução`
[1] 827779.4

# --- 5. Preparação dos dados para o gráfico ---

# Solução ótima

x1_opt <- solucao_qp$solution[1]
x2_opt <- solucao_qp$solution[2]

# Função objetivo

f_obj <- function(x1, x2) 
{
  0.5 * c(x1, x2) %*% Dmat %*% c(x1, x2) - sum(dvec * c(x1, x2))
}

# Geração do grid para visualização

x1_vals <- seq(0, 110, length.out = 300)
x2_vals <- seq(0, 110, length.out = 300)

grid <- expand.grid(x1 = x1_vals, x2 = x2_vals)

grid$restr1 <- 10 * grid$x1 + 15 * grid$x2
grid$restr2 <- -grid$x1 - grid$x2
grid$restr3 <- grid$x1
grid$restr4 <- grid$x2

grid$viavel <- with(grid, restr1 >= 1850 &
                          restr2 >= -140 &
                          restr3 >= 0 &
                          restr4 >= 0)

grid$obj <- NA
grid$obj[grid$viavel] <- apply(grid[grid$viavel, c("x1", "x2")], 1,
                               function(row) f_obj(row[1], row[2]))
grid_viavel <- filter(grid, viavel)

# --- 6. Geração do gráfico ---

ggplot() +
  geom_raster(data = grid_viavel, aes(x = x1, y = x2, fill = obj),
              interpolate = TRUE, alpha = 0.8) +
  geom_contour(data = grid_viavel, aes(x = x1, y = x2, z = obj),
               bins = 10,
               color = "black",
               size = 0.3) +
  geom_point(aes(x = x1_opt, y = x2_opt), color = "red", size = 3) +
  geom_text(aes(x = x1_opt, y = x2_opt,
                label = paste0("Solução Ótima:\n(",
                               round(x1_opt, 0), ", ",
                               round(x2_opt, 0), ")")),
            hjust = 0.5, vjust = -0.5, size = 4, color = "black") +
  scale_fill_gradientn(name = "Valor da\nFunção Objetivo",
                       colours = rainbow(7, s = 0.5, v = 0.9)) +
  labs(title = "",
       x = expression(x[1]),
       y = expression(x[2])) +
  theme_minimal() +
  theme(legend.position = "right",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.14.. Mapa de contorno da função objetivo associada ao problema de programação quadrática, com indicação da região viável definida pelas restrições lineares. O ponto em vermelho representa a solução ótima obtida pelo otimizador solve.QP.

Portanto, a solução ótima do problema de programação quadrática determina que, para minimizar o custo operacional, devem ser adotados níveis de produção de aproximadamente 26 unidades para o processo \(x_1\) e 106 unidades para o processo \(x_2\). Essa alocação de produção, que respeita as restrições de capacidade e não negatividade, resulta em um custo operacional total mínimo de aproximadamente R$ 827779,40.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.4. lm

---

O otimizador lm (Linear Models), disponível no ambiente R, constitui o procedimento padrão para a estimação de parâmetros em modelos lineares por meio do método dos mínimos quadrados ordinários. Em modelos lineares, a relação entre a variável resposta e as variáveis explicativas é expressa como uma combinação linear dos parâmetros, ainda que as variáveis preditoras possam ser submetidas a transformações não lineares, desde que a relação permaneça linear em relação aos coeficientes do modelo. Essa flexibilidade confere aos modelos lineares ampla aplicabilidade em diversas áreas do conhecimento, destacando-se particularmente na ciência de dados, onde são empregados tanto em análises exploratórias quanto em processos preditivos e inferenciais.

Otimizador 3.4 (lm). O lm é um otimizador nativo do R destinado ao ajuste de modelos lineares por mínimos quadrados ordinários. Este otimizador, em particular, resolve de forma analítica o sistema normal de equações oriundo do método dos mínimos quadrados, sendo amplamente empregado em regressão linear simples, múltipla e em modelos lineares que envolvam transformações das variáveis preditoras, desde que a relação seja linear nos parâmetros. Os principais argumentos desse otimizador são:

• formula: especificação do modelo linear na forma de uma fórmula do tipo Y ~ X1 + X2 + ....
• data: conjunto de dados utilizado para avaliar os termos da fórmula.
• weights: vetor de pesos para estimação ponderada dos parâmetros, opcional.
• subset: especifica subconjuntos do conjunto de dados a serem utilizados no ajuste.
• na.action: define o tratamento de valores ausentes no conjunto de dados.

Exemplo 3.19. Em Estatística, um modelo de regressão linear possui, em sua forma escalar, a seguinte estrutura: \[\begin{align}\\ y_i = \beta_0 + \sum_{j=1}^p \beta_j x_{ij} + \epsilon_i\\\\ \end{align}\] em que \(y_i\), para \(i = 1, \dots, n\), são as observações da variável resposta; \(x_{ij}\) representam as observações das variáveis preditoras; \(\beta_k\), com \(k = 0, \dots, p\), são os parâmetros de regressão; e os erros \(\epsilon_i \sim N(0, \sigma^2)\) são assumidos como independentes e homocedásticos. A partir desse modelo, obtém-se que a média condicional de \(y_i\) dado \(\mathbf{x}_i = (x_{i1}, \dots, x_{ip})^\top\) é: \[\begin{align}\\ E[y_i \mid \mathbf{x}_i] = \beta_0 + \sum_{j=1}^p \beta_j x_{ij}\\\\ \end{align}\] e a variância condicional de \(y_i\) dado \(\mathbf{x}_i\) é \(\text{Var}(y_i \mid \mathbf{x}_i) = \sigma^2\). Por simplicidade, o modelo pode ser reescrito na forma matricial: \[\begin{align}\\ \mathbf{y} = \mathbf{X} \beta + \epsilon \quad \rightarrow \quad \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & \dots & x_{1p} \\ 1 & x_{21} & \dots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \dots & x_{np} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix} + \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix}\\\\ \end{align}\] de modo que \(E(\mathbf{y} \mid \mathbf{X}) = \mathbf{X} \beta\) e \(\text{Cov}(\mathbf{y} \mid \mathbf{X}) = \sigma^2 \mathbf{I}_n\), onde \(\mathbf{I}_n\) é a matriz identidade de ordem \(n\). Os parâmetros são estimados pelo método dos mínimos quadrados ordinários, o qual consiste em minimizar: \[\begin{align}\\ Z^2 = \lVert \mathbf{y} - \mathbf{X} \beta \rVert^2\\\\ \end{align}\] em que \(\lVert \cdot \rVert\) representa a norma Euclidiana. A solução deste problema leva ao estimador dos mínimos quadrados ordinários para o vetor de parâmetros \(\beta\): \[\begin{align}\\ \hat{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}\\\\ \end{align}\] Se os pressupostos do modelo forem válidos, o valor esperado do vetor de coeficientes é: \[\begin{align}\\ E[\hat{\beta}] &= E[(\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}] \\\\ &= (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{X} E[\mathbf{y}] = \beta\\\\ \end{align}\] uma vez que \(E[\mathbf{y}] = \mathbf{X} \beta\). Assim, \(\hat{\beta}\) é um estimador não-viciado para o vetor de parâmetros \(\beta\). Por sua vez, a matriz de covariância de \(\hat{\beta}\), dado que \(\text{Cov}(\mathbf{y}) = \sigma^2 \mathbf{I}_n\), é descrita por: \[\begin{align}\\ \text{Cov}(\hat{\beta}) &= \text{Cov}((\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}) \\\\ &= (\mathbf{X}^\top \mathbf{X})^{-1} \text{Cov}(\mathbf{y}) \mathbf{X} (\mathbf{X}^\top \mathbf{X})^{-1} \\\\ &= \sigma^2 (\mathbf{X}^\top \mathbf{X})^{-1}\\\\ \end{align}\]

Para ilustrar a aplicação do otimizador lm em um contexto prático, considera-se o conjunto de dados mtcars, disponível no ambiente R, o qual reúne informações técnicas de diversos modelos de automóveis fabricados na década de 1970. Esses dados incluem, entre outras variáveis, o consumo de combustível medido em milhas por galão (mpg), o peso do veículo em milhares de libras (wt) e a potência do motor em cavalos-vapor (hp). O estudo da relação entre essas variáveis reveste-se de grande relevância em engenharia automotiva, economia de combustíveis e estudos ambientais, dada a busca constante por maior eficiência energética e redução de emissões de gases poluentes. Neste sentido, nos interessa aqui compreender como o consumo de combustível (\(Y_i\)), representado pela variável mpg, é influenciado pelo peso do automóvel (\(X_{1i}\), variável wt) e pela potência do motor (\(X_{2i}\), variável hp) na base mtcars. O modelo linear proposto para esta análise é expresso por:

\[\begin{align}\\ y_i = \beta_0 + \beta_1 \, x_{1i} + \beta_2 \, x_{2i} + \varepsilon_i, \quad i = 1, \dots, n\\\\ \end{align}\]

em que \(y_i\) denota a variável resposta correspondente ao consumo em milhas por galão, \(\beta_0\) representa o intercepto, \(\beta_1\) quantifica o efeito do peso do veículo sobre o consumo, e \(\beta_2\) expressa a variação média no consumo associada a alterações na potência do motor. O termo de erro \(\varepsilon_i\) captura fontes de variação não explicadas pelo modelo, tais como diferenças na aerodinâmica, tipo de transmissão, ou variabilidade inerente às medições. Para o ajuste do modelo linear proposto, será considerado o otimizador lm. O Código 3.19 apresenta uma rotina, em ambiente R, que descreve essa estimação, enquanto a Figura 3.15 ilustra os gráficos de diagnóstico do modelo, oferecendo subsídios para avaliar a adequação do mesmo aos dados do conjunto mtcars.

Código 3.19. Implementação em R da estimação dos parâmetros de um modelo de regressão linear múltipla para descrever a relação entre o consumo de combustível (mpg), o peso do veículo (wt) e a potência do motor (hp), utilizando o otimizador lm.

# --------------------------------------------
# Otimizador `lm`: Ajuste de Modelos Lineares
# --------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)
library(performance)
library(see)

# --- 1. Dados ---

data(mtcars)

# --- 2. Ajuste do modelo linear múltiplo ---

modelo_lm     <- lm(mpg ~ wt + hp, data = mtcars)

# --- 3. Exibição do sumário do ajuste ---

summary(modelo_lm)

Call:
lm(formula = mpg ~ wt + hp, data = mtcars)

Residuals:
   Min     1Q Median     3Q    Max 
-3.941 -1.600 -0.182  1.050  5.854 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
wt          -3.87783    0.63273  -6.129 1.12e-06 ***
hp          -0.03177    0.00903  -3.519  0.00145 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared:  0.8268,    Adjusted R-squared:  0.8148 
F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

# --- 4. Execução do diagnóstico ---

check_model(modelo_lm)

Figura 3.15.. Painel de diagnóstico gerado pela função check_model, do pacote performance, aplicado ao modelo linear múltiplo ajustado aos dados mtcars. Os gráficos permitem avaliar suposições como normalidade dos resíduos, linearidade, homocedasticidade e presença de valores influentes.

A partir dos resultados obtidos, observa-se que tanto o peso do veículo quanto a potência do motor exercem efeito negativo sobre o consumo de combustível, indicando que veículos mais pesados ou mais potentes tendem a apresentar menores valores de mpg. O parâmetro \(\hat{\beta}_1\) associado a wt quantifica a redução média do consumo a cada unidade adicional de peso, enquanto \(\hat{\beta}_2\) associado a hp revela o impacto do aumento da potência na eficiência energética. Por outro lado, a análise do painel de diagnóstico indica que, para o modelo ajustado, os resíduos se distribuem aproximadamente de forma simétrica e sem padrões sistemáticos quando plotados contra os valores ajustados, sugerindo que os pressupostos de linearidade e homocedasticidade estão, em grande parte, satisfeitos. Além disso, a ausência de observações altamente influentes corrobora a estabilidade dos coeficientes estimados.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.5. nls

---

O otimizador nls (Nonlinear Least Squares), disponível no ambiente R, é uma ferramenta essencial para a estimação de parâmetros em modelos não-lineares por mínimos quadrados. Diferentemente dos métodos lineares, em que a relação entre variáveis explicativas e a variável resposta é expressa como combinação linear dos parâmetros, o nls permite modelar situações em que essa relação é intrinsecamente não-linear, sendo particularmente útil em áreas como farmacocinética, crescimento populacional, química, engenharia e muitas outras aplicações científicas.

O nls, em particular, busca encontrar os valores dos parâmetros que minimizam a soma dos quadrados dos resíduos entre os valores observados e os valores previstos pelo modelo. Para tanto, requer a definição explícita da expressão matemática do modelo, valores iniciais plausíveis para os parâmetros e, em alguns casos, opções específicas de algoritmos de otimização, como Gauss-Newton ou algoritmo de Levenberg-Marquardt, para lidar com a eventual complexidade da superfície de mínimos.

Otimizador 3.5 (nls). O nls é um otimizador nativo do R para ajuste de modelos não-lineares por mínimos quadrados. Este otimizador, em particular, utiliza algoritmos iterativos de otimização numérica para estimar parâmetros em modelos com estrutura funcional não-linear, sendo empregado, em geral, para ajuste de curvas e modelagem não-linear. Os principais argumentos desse otimizador são:

• formula: especificação do modelo não-linear na forma de uma fórmula do tipo Y ~ f(x).
• data: refere-se ao conjunto de dados utilizado para avaliar os termos da fórmula.
• start: lista com os valores iniciais para os parâmetros do modelo a serem estimados.
• control: lista de opções de controle, incluindo tolerância, número máximo de iterações, etc.
• algorithm: define o algoritmo utilizado na estimação, podendo ser default (Gauss-Newton), plinear ou port.

Exemplo 3.20. Em Estatística, modelos de regressão não-linear têm por objetivo descrever a relação entre uma variável resposta \(Y\) e um conjunto de variáveis explicativas por meio de uma função determinística não-linear nos parâmetros, incorporando um componente estocástico que captura a variabilidade residual observada nos dados. A formulação geral desse tipo de modelo é dada por:

\[\begin{align}\\ y_i = f(x_i; \theta) + \varepsilon_i, \quad i = 1, \dots, n\\\\ \end{align}\]

em que \(y_i\) representa as observações da variável resposta, \(f(\cdot)\) é uma função contínua e não-linear em relação ao vetor de parâmetros \(\theta\), definida sobre o vetor de covariáveis \(x_i\), e os termos de erro \(\varepsilon_i\) são variáveis aleatórias independentes e identicamente distribuídas. Em muitas aplicações, assume-se que \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\), embora essa suposição não seja estritamente necessária para a estimação dos parâmetros via mínimos quadrados não-lineares.

Para ilustrar a aplicação do otimizador nls em um contexto prático, considera-se um problema relacionado à modelagem de ondas oceânicas. Tais ondas são perturbações na superfície dos oceanos que desempenham papel essencial na dinâmica costeira, afetando a erosão, a estabilidade de obras marítimas e a segurança de operações em mar aberto. Em condições de mar totalmente desenvolvido, com vento constante ao longo de uma extensão oceânica (Fetch), estudos empíricos indicam que a altura significativa das ondas, adimensionalizada por \(gH/U^2\), denotada por \(y_i\), tende a crescer de forma exponencial em função do Fetch, também adimensionalizado por \(gF/U^2\), denotado por \(x_i\), sobretudo nas fases iniciais de desenvolvimento das ondas. Neste contexto, considere o seguinte modelo não-linear:

\[\begin{align}\\ y_i = \theta_1 \, x^{\theta_2} + \varepsilon_i, \quad i = 1, \dots, n\\\\ \end{align}\]

em que os parâmetros \(\theta_1\) e \(\theta_2\) representam, respectivamente, um fator de escala associado à altura basal das ondas e a taxa de crescimento da altura em função do Fetch. O termo de erro \(\varepsilon_i\) incorpora fontes de variação não observadas, como flutuações na intensidade do vento, correntes marinhas e possíveis imprecisões instrumentais. Para fins ilustrativos, considera-se uma situação em que a intensidade do vento é constante em \(2{.}5 \, \text{m/s}\). Os valores da variável preditora \(x_i\) serão gerados a partir de uma distribuição exponencial com parâmetro \(\lambda = 3\), ao passo que os valores da variável resposta \(y_i\) serão simulados segundo o modelo acima, utilizando \(\theta_1 = 2\), \(\theta_2 = 0{.}75\), e erros \(\varepsilon_i \sim \mathcal{N}(0, 0{.}1^2)\). O ajuste do modelo será realizado utilizando o otimizador nls, que implementa um procedimento iterativo baseado em mínimos quadrados não-lineares, a partir de valores iniciais especificados para os parâmetros do modelo. O Código 3.20 ilustra esse processo, e a Figura 3.16 exibe o a curva de ajuste do modelo proposto para estudar a relação entre a altura de ondas oceânicas e a medida de Fetch.

Código 3.20. Implementação em R da estimação dos parâmetros de um modelo de regressão não-linear para descrever a relação entre a altura das ondas oceânicas e a medida de Fetch, em condições de mar totalmente desenvolvido, utilizando o otimizador nls.

# -------------------------------------------------
# Otimizador `nls`: Ajuste de Modelos Não-Lineares
# -------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Geração dos dados simulados ---

set.seed(123)

n       <- 5000
x       <- rexp(n, rate = 3)          # Medida de Fetch (gF/U)
theta1  <- 2
theta2  <- 0.75

y       <- theta1 * x^theta2 + rnorm(n, mean = 0, sd = 0.1)  # Altura (gH/U)

# --- 2. Ajuste do modelo não linear via `nls` ---

modelo_nls <- nls(y ~ a * x^b,
                  start = list(a = 1, b = 0.5))

# --- 3. Exibição do sumário do ajuste ---

summary(modelo_nls)

Formula: y ~ a * x^b

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a 2.001513   0.003201   625.3   <2e-16 ***
b 0.751231   0.002004   374.8   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1008 on 4998 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 2.154e-06

# --- 4. Preparação dos dados para o gráfico ---

# Coeficientes estimados

param <- coef(modelo_nls)

# Valores ajustados para a curva

x_seq <- seq(min(x), max(x), length.out = 200)
y_fit <- param["a"] * x_seq^param["b"]

# Data frames para plotagem

df_dados  <- data.frame(x = x, y = y, grupo = "Dados Observados")
df_ajuste <- data.frame(x = x_seq, y = y_fit, grupo = "Curva Ajustada")

# --- 5. Geração do gráfico ---

ggplot() +
  geom_point(data = df_dados, aes(x = x, y = y, color = grupo),
             alpha = 0.4, size = 1) +
  geom_line(data = df_ajuste, aes(x = x, y = y, color = grupo),
            size = 1) +
  scale_color_manual(name = "Legenda",
                     values = c("Dados Observados" = "blue",
                                "Curva Ajustada" = "red")) +
  labs(title = "",
       x = "Medida de Fetch (gF/U)",
       y = "Altura das Ondas (gH/U)") +
  theme_minimal() +
  theme(legend.position = "right",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.16.. Ajuste do modelo não-linear aos dados simulados, relacionando a altura das ondas (gH/U) à medida de fetch (gF/U). A curva em vermelho representa a função estimada pelo otimizador nls, enquanto os pontos azuis correspondem aos dados observados.

A partir dos resultados obtidos, conclui-se que o modelo ajustado descreve adequadamente como a altura das ondas cresce com a medida de Fetch em condições de mar totalmente desenvolvido. O parâmetro \(\hat{\theta}_1 \approx 2\) indica que, para valores iniciais de Fetch, a altura base das ondas gira em torno de 2 unidades (em escala adimensional, gH/U). Já o parâmetro \(\hat{\theta}_2 \approx 0.75\) mostra que esse crescimento segue uma taxa exponencial relativamente intensa à medida que o Fetch aumenta. Na prática, isso significa que, à medida que o vento sopra sobre uma distância maior (maior Fetch), a energia transferida para a superfície do mar também aumenta, fazendo com que as ondas cresçam mais rapidamente. O bom ajuste do modelo e a baixa variabilidade dos resíduos indicam que essa relação é consistente, podendo ser usada para prever alturas de onda com base em medidas de Fetch em cenários reais de engenharia oceânica e planejamento costeiro.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.6. ga

---

O otimizador ga, disponível no ambiente R por meio do pacote GA, é uma implementação robusta dos algoritmos genéticos (Genetic Algorithms), técnica metaheurística inspirada nos processos evolutivos naturais, como seleção, cruzamento e mutação. Esse otimizador, em particular, permite a otimização de funções objetivo definidas pelo usuário, sendo capaz de lidar com variáveis contínuas, inteiras ou binárias, e possibilita a inclusão de restrições através de funções penalidade. O pacote oferece também ferramentas para monitoramento visual do processo evolutivo, análise de convergência e ajustes finos dos parâmetros do algoritmo, como tamanho da população, probabilidade de mutação e critérios de parada.

Otimizador 3.6 (ga). O ga é um otimizador disponível no pacote GA que implementa algoritmos genéticos para resolver problemas de otimização global. Este otimizador, em particular, adota uma abordagem heurística baseada em princípios da evolução biológica, sendo especialmente eficaz em espaços de busca complexos, não-lineares e não-convexos. A otimização ocorre por meio da evolução iterativa de populações de soluções, utilizando operadores genéticos como seleção, cruzamento e mutação para explorar e refinar o espaço de soluções em direção ao ótimo global. Os principais argumentos desse otimizador são:

• type: especifica o tipo de problema a ser resolvido, como "real-valued" para variáveis contínuas.
• fitness: função objetivo a ser maximizada, fornecida pelo usuário.
• lower, upper: vetores que definem os limites inferior e superior do espaço de busca para cada variável.
• popSize: número de indivíduos na população inicial.
• maxiter: número máximo de gerações (iterações).
• run: número de gerações sem melhoria na melhor solução antes da parada automática.
• pmutation, pcrossover: probabilidades de mutação e cruzamento aplicadas nas gerações sucessivas.

Exemplo 3.21 (Otimização do Crescimento Populacional de Bactérias). Considere um experimento de microbiologia no qual se observa o crescimento de uma colônia de bactérias ao longo do tempo. A variável resposta \(y_i\) representa a densidade populacional (em unidades arbitrárias) no instante de tempo \(x_i\) (em horas). Em muitos sistemas biológicos, o crescimento de populações segue um padrão sigmoidal, frequentemente modelado pela função logística:

\[\begin{align}\\ y_i = \frac{\theta_1}{1 + \exp[-\theta_2(x_i - \theta_3)]} + \varepsilon_i\\\\ \end{align}\]

em que \(\theta_1\) é o valor assintótico máximo da população (capacidade de suporte), \(\theta_2\) controla a taxa de crescimento, \(\theta_3\) representa o tempo correspondente ao ponto de inflexão da curva (crescimento máximo), e \(\varepsilon_i\) representa o erro aleatório associado a variações experimentais. Para o ajuste deste modelo, será considerado o otimizador ga do pacote GA, que implementa algoritmos genéticos para otimização global. O ajuste visa encontrar os valores de \(\theta_1\), \(\theta_2\) e \(\theta_3\) que minimizam a soma dos quadrados dos resíduos entre os valores observados \(y_i\) e os valores preditos pelo modelo logístico. O Código 3.21, em ambiente R, implementa esse processo para \(\theta_1 = 100\), \(\theta_2 = 0.5\) e \(\theta_3 = 12\), e uma população de tamanho \(N=1000\). Já a Figura 3.17 exibe o a curva de ajuste do modelo proposto obtida pela aplicação do algoritmo genético.

Código 3.21. Implementação em R para estimação dos parâmetros de um modelo de crescimento populacional de bactérias utilizando o otimizador ga do pacote GA.

# -----------------------------------------------------------
# Otimizador `ga` do pacote `GA`: Ajuste de Modelo Logístico
# -----------------------------------------------------------

# --- 0. Pacotes necessários ---

library(GA)
library(ggplot2)

# --- 1. Geração dos dados simulados ---

set.seed(123)

n       <- 1000
x       <- seq(0, 24, length.out = n)

theta1  <- 100   # Capacidade de Suporte
theta2  <- 0.5   # Taxa de Crescimento
theta3  <- 12    # Ponto de Inflexão

y       <- theta1 / (1 + exp(-theta2 * (x - theta3))) + rnorm(n, 0, 3)

# --- 2. Definição da função objetivo ---

func_obj  <- function(par) 
{
  y_hat   <- par[1] / (1 + exp(-par[2] * (x - par[3])))
  sum((y - y_hat)^2)
}

# --- 3. Ajuste via Algoritmo Genético ---

ajuste_ga <- ga(type = "real-valued",
                fitness = function(par) -func_obj(par),
                lower = c(50, 0.1, 5),
                upper = c(150, 1.5, 20),
                popSize = 100,
                maxiter = 200,
                run = 50,
                pmutation = 0.2,
                crossover = gareal_blxCrossover,
                selection = gareal_lrSelection,
                mutation = gareal_raMutation,
                seed = 42)

# --- 4. Parâmetros estimados ---

parametros_estimados        <- ajuste_ga@solution[1, ]
names(parametros_estimados) <- c("theta1", "theta2", "theta3")
print(parametros_estimados)

     theta1      theta2      theta3 
100.0004150   0.5060415  12.0061799 

# --- 5. Predições do modelo ajustado ---

y_pred <- parametros_estimados["theta1"] / 
          (1 + exp(-parametros_estimados["theta2"] * 
                   (x - parametros_estimados["theta3"])))
head(y_pred, n = 50)

 [1] 0.2292950 0.2320931 0.2349252 0.2377918 0.2406933 0.2436302 0.2466027
 [8] 0.2496115 0.2526569 0.2557393 0.2588592 0.2620171 0.2652135 0.2684487
[15] 0.2717232 0.2750377 0.2783924 0.2817879 0.2852247 0.2887033 0.2922242
[22] 0.2957880 0.2993950 0.3030459 0.3067412 0.3104814 0.3142671 0.3180988
[29] 0.3219770 0.3259024 0.3298755 0.3338968 0.3379670 0.3420867 0.3462564
[36] 0.3504767 0.3547483 0.3590718 0.3634478 0.3678769 0.3723597 0.3768970
[43] 0.3814894 0.3861375 0.3908420 0.3956037 0.4004231 0.4053009 0.4102380
[50] 0.4152349

# --- 6. Preparação dos dados para o gráfico ---

df_plot <- data.frame(Tempo = rep(x, 2),
                      Valor = c(y, y_pred),
                      Tipo = factor(rep(c("Observado", "Predito"), each = length(x))))

# --- 7. Geração do gráfico ---

ggplot(df_plot, aes(x = Tempo, y = Valor, color = Tipo)) +
  geom_point(data = subset(df_plot, Tipo == "Observado"),
             size = 2, alpha = 0.6) +
  geom_line(data = subset(df_plot, Tipo == "Predito"), linewidth = 0.5) +
  scale_color_manual(values = c("Observado" = "blue", "Predito" = "red")) +
  labs(title = "",
       x = "Tempo (horas)",
       y = "Densidade Populacional",
       color = "Legenda") +
  theme_minimal() +
  theme(legend.position = "right",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.17.. Ajuste do modelo logístico ao crescimento bacteriano utilizando o algoritmo genético do pacote GA. Os pontos azuis representam os dados observados, enquanto a curva vermelha corresponde aos valores preditos pelo modelo ajustado.

O modelo logístico ajustado descreve de forma precisa o crescimento bacteriano ao longo do tempo, capturando as fases características do processo biológico, que incluem a fase exponencial inicial, o ponto de inflexão onde ocorre a desaceleração do crescimento e a estabilização da população próxima à capacidade máxima do ambiente. Os parâmetros estimados pelo algoritmo genético indicam que a densidade máxima da população bacteriana, representada por \(\hat{\theta}_1 \approx 100\), corresponde à capacidade de suporte do meio; a taxa intrínseca de crescimento, dada por \(\hat{\theta}_2 \approx 0{.}5\), reflete a velocidade do aumento populacional durante a fase exponencial; e o tempo do ponto de inflexão, \(\hat{\theta}_3 \approx 12\), indica o momento em que o crescimento atinge sua aceleração máxima antes de desacelerar. A partir desses valores, conclui-se que a população bacteriana cresce rapidamente até atingir cerca de 100 unidades, atingindo o ponto de desaceleração do crescimento próximo ao tempo de 12 horas, quando fatores limitantes começam a influenciar o processo.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.7. maxLik

---

O otimizador maxLik, disponível no ambiente R através do pacote maxLik, é uma ferramenta projetada para a estimação de parâmetros via métodos de máxima verossimilhança. Este pacote oferece uma infraestrutura versátil que permite ao usuário especificar funções de log-verossimilhança de forma flexível, sendo capaz de lidar tanto com funções diferenciáveis, para as quais gradientes e Hessianas podem ser fornecidos ou calculados numericamente, quanto com funções não diferenciáveis.

O maxLik também disponibiliza diversos métodos numéricos para a otimização, incluindo os baseados em gradiente, como BFGS e Newton-Raphson, bem como métodos robustos, como SANN, para casos em que a função objetivo apresenta comportamento complexo ou superfícies irregulares. Além do procedimento de otimização propriamente dito, o pacote também fornece facilidades para obtenção das matrizes de variância-covariância dos estimadores, cálculo de erros padrão e testes de hipótese.

Otimizador 3.7 (maxLik). O maxLik é um otimizador disponível no pacote maxLik, desenvolvido especificamente para estimação de parâmetros via máxima verossimilhança. Este otimizador, em particular, fornece uma estrutura flexível para lidar com funções de log-verossimilhança de forma geral, permitindo especificação diretamente a função objetivo (log-verossimilhança), bem como (opcionalmente) suas derivadas de primeira e segunda ordem. A otimização ocorre por meio de algoritmos numéricos bem estabelecidos. Os principais argumentos desse otimizador são:

• logLik: especifica a função de log-verossimilhança. Deve receber como primeiro argumento o vetor de parâmetros e retornar um único valor escalar (log-verossimilhança total) ou um vetor numérico com os valores da log-verossimilhança por observação.
• hess: especifica a matriz Hessiana da log-verossimilhança. Opcional. Deve retornar uma matriz quadrada com dimensões iguais ao número de parâmetros. Se ausente, a Hessiana será estimada numericamente.
• start: vetor numérico contendo os valores iniciais dos parâmetros a serem otimizados. Pode conter nomes, os quais são herdados pelos resultados, facilitando a interpretação dos coeficientes estimados.
• method: especifica o algoritmo de otimização a ser utilizado. Os métodos disponíveis são: NR (Newton-Raphson), BFGS (Broyden-Fletcher-Goldfarb-Shanno), BFGSR (algoritmo BFGS implementado em R), BHHH (Berndt-Hall-Hall-Hausman), SANN (Simulated ANNealing), CG (Gradiente Conjugado) ou NM (Nelder-Mead).

Exemplo 3.22 (Distribuição Exponencial). Considere uma distribuição exponencial com função densidade de probabilidade dada por:

\[\begin{align}\\ f(x;\lambda) = \lambda e^{-\lambda x}, \quad x > 0, \ \lambda > 0\\\\ \end{align}\]

onde \(\lambda\) é o parâmetro de taxa. Neste exemplo, busca-se estimar o parâmetro \(\lambda\) por meio do método da máxima verossimilhança. Suponha-se uma amostra \(x_1, x_2, \ldots, x_n\) independente e identicamente distribuída segundo a distribuição exponencial. A função de log-verossimilhança associada é dada por:

\[\begin{align}\\ \ell(\lambda; \mathbf{x}) = n \log(\lambda) - \lambda \sum_{i=1}^n x_i\\\\ \end{align}\]

a qual será maximizada numericamente utilizando o otimizador maxLik, disponível no pacote maxLik. O Código 3.22, em ambiente R, implementa esse processo considerando uma amostra de tamanho \(n = 1000\) gerada a partir de uma distribuição exponencial com \(\lambda = 5\), e o otimizador maxLik() do pacote maxLik. A Figura 3.18 exibe o a curva de ajuste da distribuição exponencial obtida via maxLik().

Código 3.22. Implementação em R para estimação do parâmetro de uma distribuição exponencial via máxima verossimilhança, utilizando o otimizador maxLik do pacote homônimo.

# -----------------------------------------------------------------------------
# Otimizador `maxLik` do pacote `maxLik`: Estimação por Máxima Verossimilhança
# -----------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(maxLik)
library(ggplot2)

# --- 1. Geração dos dados simulados ---

set.seed(123)
x       <- rexp(1000, rate = 5)

# --- 2. Definição da função log-verossimilhança ---

ll_exp  <- function(lambda) 
{
  n     <- length(x)
  ll    <- n * log(lambda) - lambda * sum(x)
  return(ll)
}

# --- 3. Estimação via `maxLik` ---

lambda_inicial  <- 4.5
ajuste          <- maxLik(logLik = ll_exp, start = c(lambda = lambda_inicial))

# --- 4. Sumário dos resultados ---

summary(ajuste)

--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 3 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-Likelihood: 579.8992 
1  free parameters
Estimates:
       Estimate Std. error t value Pr(> t)    
lambda   4.8545     0.1532    31.7  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------

# --- 5. Preparação dos dados para o gráfico ---

lambda_hat  <- coef(ajuste)["lambda"]
df          <- data.frame(x = x)

# --- 6. Geração do gráfico ---

ggplot(df, aes(x = x)) +
  geom_histogram(aes(y = ..density..),
                 bins = 30,
                 fill = "lightblue",
                 color = "black") +
  stat_function(aes(color = "Curva Ajustada"),
                fun = dexp,
                args = list(rate = lambda_hat),
                size = 0.5) +
  scale_color_manual(name = "", values = c("Curva Ajustada" = "red")) +
  labs(title = "",
       x = "x",
       y = "Densidade") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.18.. Histograma dos dados simulados a partir da distribuição exponencial, com sobreposição da curva ajustada obtida por máxima verossimilhança utilizando o otimizador maxLik. A linha vermelha representa a densidade teórica estimada.

A partir dos resultados obtidos, conclui-se que o parâmetro da taxa da distribuição exponencial foi estimado com precisão adequada. A estimativa \(\hat{\lambda} \approx 4.85\) indica que a taxa de ocorrência dos eventos modelados está próxima do valor verdadeiro da simulação, que era 5. O erro padrão relativamente baixo (\(\approx 0.15\)) e o valor-p significativo confirmam a confiabilidade dessa estimativa. Na prática, isso significa que o modelo ajustado representa bem o comportamento dos dados observados. A rápida convergência do método de Newton-Raphson e o gradiente próximo de zero evidenciam a eficiência e adequação do procedimento de máxima verossimilhança aplicado.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

3.4.8. fitdist

---

O otimizador fitdist, disponível em ambiente R através do pacote fitdistrplus, é uma ferramenta para o ajuste paramétrico de distribuições de probabilidade a conjuntos de dados observados. Esse otimizador permite estimar os parâmetros de diferentes distribuições contínuas ou discretas por métodos como máxima verossimilhança, método dos momentos ou outros critérios ajustáveis. Além disso, por meio de sua integração com outras funções do pacote, como plotdist e gofstat, o fitdist facilita uma análise completa do ajuste, promovendo uma compreensão mais profunda do comportamento dos dados e suporte à tomada de decisão baseada em modelos probabilísticos adequados.

Otimizador 3.8 (fitdist). O fitdist é um otimizador disponível no pacote fitdistrplus que realiza o ajuste de modelos de probabilidade a partir de dados observados. Diferentemente de otimizadores que recebem diretamente uma função de log-verossimilhança, o fitdist requer a especificação das funções de densidade de probabilidade e da função de distribuição acumulada do modelo, além da função quantil, opcionalmente, quando o método de otimização selecionado exige esse componente (como no caso do método qme). Esse otimizador permite diversas estratégias de estimação de parâmetros, incluindo máxima verossimilhança, método dos momentos, correspondência de quantis, máxima bondade de ajuste e espaçamento máximo. Os principais argumentos desse otimizador são:

• data: vetor numérico contendo os dados amostrais que serão usados para ajustar o modelo probabilístico.
• distr: string que especifica o nome do modelo de probabilidade a ser ajustado, considerando que as funções do modelo devem seguir o padrão dname para densidade (ou probabilidade), pname para distribuição acumulada e qname para quantis.
• method: string que define o método de estimação dos parâmetros, podendo assumir os valores mle (máxima verossimilhança), mme (método dos momentos), qme (correspondência de quantis), mge (máxima bondade de ajuste) e mse (espaçamento máximo).
• start: vetor numérico contendo os valores iniciais para os parâmetros do modelo, que servem como ponto de partida para os algoritmos de otimização utilizados.
• discrete: argumento lógico que especifica se o modelo de probabilidade é discreto. O padrão é FALSE.

Exemplo 3.23 (Distribuição Exponencial). Considere, novamente, uma distribuição exponencial com função densidade de probabilidade dada por:

\[\begin{align}\\ f(x;\lambda) = \lambda e^{-\lambda x}, \quad x > 0, \ \lambda > 0\\\\ \end{align}\]

onde \(\lambda\) é o parâmetro de taxa. Neste exemplo, objetiva-se estimar o parâmetro \(\lambda\) por meio do método da máxima verossimilhança, utilizando o otimizador fitdist, disponível no pacote fitdistrplus. Como a distribuição exponencial é nativamente suportada pelo ambiente R, não há necessidade de implementar explicitamente as funções densidade, distribuição acumulada ou quantil, necessárias para o fitdist.O Código 3.23, em ambiente R, ilustra esse processo considerando uma amostra de tamanho \(n = 1000\) gerada a partir de uma distribuição exponencial com \(\lambda = 5\), e o otimizador fitdist do pacote fitdistrplus. A Figura 3.19 exibe o a curva de ajuste da distribuição exponencial obtida via fitdist().

Código 3.23. Implementação em R para estimação do parâmetro de uma distribuição exponencial via máxima verossimilhança, utilizando o otimizador fitdist do pacote fitdistrplus.

# ----------------------------------------------------------------------------------
# Otimizador `fitdist` do pacote `fitdistrplus`: Ajuste de Distribuição Paramétrica
# ----------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(fitdistrplus)
library(ggplot2)

# --- 1. Geração dos dados simulados ---

set.seed(123)
x       <- rexp(1000, rate = 5)

# --- 2. Ajuste da distribuição via `fitdist` ---

ajuste  <- fitdist(data = x, distr = "exp", method = "mle")

# --- 3. Sumário dos resultados ---

summary(ajuste)

Fitting of the distribution ' exp ' by maximum likelihood 
Parameters : 
     estimate Std. Error
rate 4.854467  0.1535117
Loglikelihood:  579.8992   AIC:  -1157.798   BIC:  -1152.891 

# --- 4. Preparação dos dados para o gráfico ---

df        <- data.frame(x = x)
rate_est  <- ajuste$estimate["rate"]

# --- 5. Geração do gráfico ---

ggplot(df, aes(x = x)) +
  geom_histogram(aes(y = ..density..),
                 bins = 30,
                 fill = "lightblue",
                 color = "black") +
  stat_function(aes(color = "Curva Ajustada"),
                fun = dexp,
                args = list(rate = rate_est),
                size = 0.5) +
  scale_color_manual(name = "", values = c("Curva Ajustada" = "red")) +
  labs(title = "",
       x = "x",
       y = "Densidade") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 3.19.. Histograma dos dados simulados provenientes da distribuição exponencial, com a curva ajustada (em vermelho) obtida por máxima verossimilhança via função fitdist do pacote fitdistrplus. A sobreposição evidencia a adequação do modelo ajustado aos dados observados.

A partir dos resultados obtidos, observa-se que a estimativa do parâmetro taxa \(\hat{\lambda} \approx 4.85\) (valor obtido em fit$estimate) está próxima do valor verdadeiro utilizado na simulação (\(\lambda = 5\)), demonstrando boa capacidade do método de máxima verossimilhança via fitdist para recuperar parâmetros verdadeiros. O ajuste gráfico da curva densidade sobre o histograma evidencia a adequação do modelo aos dados simulados.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Andreescu, Titu, Cristinel Mortici, Marian Tetiva, Titu Andreescu, Cristinel Mortici, e Marian Tetiva. 2017. «The Intermediate Value Theorem». Mathematical Bridges, 189–200.

Bellman, Richard. 1957. Dynamic Programming. Princeton, NJ: Princeton University Press.

Bertsimas, Dimitris, e John Tsitsiklis. 1993. «Simulated annealing». Statistical Science 8 (1): 10–15.

Blørstad, Morten. 2019. «Genetic Algorithm From Scratch Using R». Online article. https://medium.com/@mortenblorstad/genetic-algorithm-from-scratch-using-r-a36dc664c46e.

Cinlar, Erhan. 2013. Introduction to Stochastic Processes. Courier Corporation.

Dantzig, George B. 1957. «Discrete-variable Extremum Problems». Em Operations Research, editado por A. J. Hoffman, R. L. Graves, e H. W. Kuhn, 331–43. New York: John Wiley & Sons.

Darwin, Charles. 1859. On the Origin of Species by Means of Natural Selection. London: John Murray.

Dobrow, Robert P. 2016. Introduction to Stochastic Processes with R. John Wiley & Sons.

Ehiwario, JC, e SO Aghamie. 2014. «Comparative study of bisection, Newton-Raphson and secant methods of root–finding problems». IOSR J. of Engineering 4 (4): 1–7.

Fisher, Ronald A. 1922. «On the mathematical foundations of theoretical statistics». Philosophical transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character 222 (594-604): 309–68.

Givens, Geof H, e Jennifer A Hoeting. 2012. Computational Statistics. John Wiley & Sons.

Guidorizzi, Hamilton Luiz. 2002. Um Curso de Cálculo. Rio de Janeiro: LTC.

Hoel, Paul G, Sidney C Port, e Charles J Stone. 1986. Introduction to Stochastic Processes. Waveland Press.

Holland, John H. 1975. Adaptation in Natural and Artificial Systems. MIT Press.

Kirkpatrick, Scott, C Daniel Gelatt Jr, e Mario P Vecchi. 1983. «Optimization by simulated annealing». Science 220 (4598): 671–80.

Kramer, Oliver, e Oliver Kramer. 2017. Genetic Algorithms. Springer.

Lawler, Gregory F. 2018. Introduction to Stochastic Processes. Chapman; Hall/CRC.

Longford, Nicholas T. 1987. «A fast scoring algorithm for maximum likelihood estimation in unbalanced mixed models with nested random effects». Biometrika 74 (4): 817–27.

Markov, Andrei Andreevich. 1906. «Rasprostranenie zakona bol’shih chisel na velichiny, zavisyaschie drug ot druga». Izvestiya Fiziko-matematicheskogo obschestva pri Kazanskom universitete 15 (135-156): 18.

Martello, Silvano, e Paolo Toth. 1990. Knapsack Problems: Algorithms and Computer Implementations. Chichester: John Wiley & Sons.

Meyn, Sean P, e Richard L Tweedie. 2012. Markov Chains and Stochastic Stability. Springer Science & Business Media.

Milne, William Edmund. 2015. Numerical Calculus. Princeton University Press.

Newton, Isaac. 1711. De analysi per aequationes numero terminorum infinitas. London: Typis Johannis Streater.

Pêcheux, Mina. 2020. «A Peek at Genetic Algorithms». Online article. https://medium.com/@mina.pecheux/a-peek-at-genetic-algorithms-123456789abc.

Pisinger, David, e Paolo Toth. 1998. «Knapsack Problems». Em Handbook of Combinatorial Optimization, editado por Ding-Zhu Du e Panos M. Pardalos, 299–428. Boston, MA: Springer.

Raphson, Joseph. 1697. Analysis Æquationum Universalis. 2nd ed. London: Thomas Bradyll.

Ruggiero, Márcia Aparecida Gomes, e Vera Lúcia da Rocha Lopes. 1998. Cálculo Numérico: Aspectos Teóricos e Computacionais. Pearson/Makron.

Sivanandam, SN, SN Deepa, SN Sivanandam, e SN Deepa. 2008. Genetic Algorithms. Springer.

Solanki, Chitra, Pragati Thapliyal, e Komal Tomar. 2014. «Role of bisection method». International Journal of Computer Applications Technology and Research 3 (8): 533–35.

Srivastava, RB, e Narendra Deo Dixit. 2012. «Numerical accuracy of bisection method». International Multidisciplinary Research Journal 2 (1).

Van Laarhoven, Peter JM, Emile HL Aarts, Peter JM van Laarhoven, e Emile HL Aarts. 1987. Simulated annealing. Springer.

Wallis, John. 1685. A Treatise of Algebra, Both Historical and Practical. Oxford: Richard Davis.

Weibull, Waloddi. 1951. «A statistical distribution function of wide applicability». Journal of applied mechanics.