Alpha C. Chiang. Elements Of Dynamic Optimization. 1992: Respostas do Capítulo 7 em R
\(\text{ }\)
Exemplo 1, página 171.
\[ \max \; V = \int_0^T - (1-u^2)^{1/2} dt\\ \text{Sujeito a:}\\ \dot{y} = u\\ y(0) = A\\y(T) \text{ livre}\\ \text{$A$, e $T$ dados} \]
O primeiro passo para a resolução do problema é instalar as bibliotecas necessárias e liberá-las para o uso. De início, considere usar a biblioteca deSolve, ideal para a resolução de equações diferenciais.
lapply(list("deSolve"), function(x){
if(x %in% installed.packages() == F){
install.packages(x)
require(x, character.only = T)
}else{
require(x, character.only = T)
}
})
FALSE Carregando pacotes exigidos: deSolve
FALSE [[1]]
FALSE [1] TRUE
Feito isso, o próximo passo é resolver o problema. Como \(A\) e \(T\) são dados, considere \(A=1\) e \(T=4\) para viabilizar a resolução.
# Parâmetros
T <- 4 # Tempo final
lambda <- 0.5 # Valor arbitrário para lambda
# Função de controle
u <- function(lambda) {
lambda / sqrt(1 - lambda^2)
}
# Equação diferencial
dy_dt <- function(t, y, params) {
list(u(params$lambda))
}
# Condição inicial
y0 <- 1
# Tempos
times <- seq(0, T, by = 0.01)
# Resolvendo a EDO
out <- ode(y = y0, times = times, func = dy_dt, parms = list(lambda = lambda))
# Plotando resultados
plot(out, xlab = "Tempo", col = "blue", lwd = 2, ylab = "y(t)", main = "Solução y(t)")
grid()

Exemplo 2, página 173.
\[ \max \; V = \int_0^2 (2y - 3u) dt\\ \text{Sujeito a:}\\ \dot{y} = y+u\\ y(0) = 4\\y(2) \text{ livre}\\ u(t) \in [0,2] \]
# Função do sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo
u <- max(0, min(2, (3/1) * lambda))
# Derivadas
dy <- y + u
dlambda <- -2 - lambda
list(c(dy, dlambda))
}
# Condições iniciais
state <- c(y = 4, lambda = 0) # Supondo lambda(2) = 0
time <- seq(0, 2, by = 0.01)
# Resolver o sistema
results <- data.frame(ode(y = state, times = time, func = dynamics, parms = NULL))
# Plotar os resultados
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória Ótima de y(t)")

Exercício 7.2: Questão 2, página 176
\[ \max \; V = \int_0^4 3y dt\\ \text{Sujeito a:}\\ \dot{y} = y+u\\ y(0) = 5\\y(4) \text{ livre}\\ u(t) \in [0,2] \]
# Função que define o sistema de equações diferenciais
# Função que define o sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo
u <- max(0, min(2, -lambda)) # Ajuste de u dentro do intervalo [0, 2]
# Derivadas
dy <- y + u
dlambda <- -3 - lambda
list(c(dy, dlambda))
}
# Condições iniciais
state <- c(y = 5, lambda = 0) # y(0) = 5
time <- seq(0, 4, by = 0.01) # Integração no tempo direto
# Resolver o sistema
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
# Plotar os resultados
results <- as.data.frame(results)
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória Ótima de y(t)")

Exercício 7.2: Questão 3, página 176
\[ \max \; V = \int_0^2(y - u^2) dt\\ \text{Sujeito a:}\\ \dot{y} = u\\ y(0) = 0\\ u(t) \text{ irrestrito} \]
# Função do sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo
u <- lambda / 2 # u não possui restrições
# Derivadas
dy <- u
dlambda <- -1
list(c(dy, dlambda))
}
# Condições iniciais
state <- c(y = 0, lambda = 0) # y(0) = 0, λ final será ajustado (integração reversa possível)
time <- seq(0, 2, by = 0.01)
# Resolver o sistema
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
# Plotar os resultados
results <- as.data.frame(results)
u_values <- results$lambda / 2
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória Ótima de y(t)")

Exercício 7.2: Questão 4, página 177
\[ \max \; \int_0^1 -\frac{1}{2}(y^2 + u^2) dt\\ \text{Sujeito a:}\\ \dot{y} = u - y\\ y(0) = 1\\ u(t) \text{ irrestrito} \]
library(deSolve)
# Função que define o sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo
u <- lambda # Condição de máximo: u = λ
# Derivadas
dy <- u - y
dlambda <- y - lambda
list(c(dy, dlambda))
}
# Condições iniciais
state <- c(y = 1, lambda = 0) # Condição inicial: y(0) = 1
time <- seq(0, 1, by = 0.01)
# Resolver o sistema
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
# Plotar os resultados
results <- as.data.frame(results)
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória Ótima de y(t)")

Exemplo 1: página 184
\[ \max \; V = \int_0^1 - u^2 dt\\ \text{Sujeito a:}\\ \dot{y} = y + u\\ y(0) = 1\\ y(1) = 0 \]
# Função para o sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo
u <- lambda / 2
# Derivadas
dy <- y + u
dlambda <- -lambda
list(c(dy, dlambda))
}
# Condições iniciais e finais
state <- c(y = 1, lambda = 0) # y(0) = 1; ajuste inicial para lambda
time <- seq(0, 1, by = 0.01)
# Resolver o sistema
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
# Processar resultados
results <- as.data.frame(results)
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória Ótima de y(t)")

Exemplo 3: página 187
\[ \max \; V = \int_0^T - 1 dt\\ \text{Sujeito a:}\\ \dot{y} = y + u\\ y(0) = 5\\ y(T) = 11\\ T \text{ Livre} \\ u(t) \in [-1,1] \]
# Função que define o sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo baseado no sinal de lambda
u <- ifelse(lambda > 0, -1, 1)
# Derivadas
dy <- y + u
dlambda <- -lambda
list(c(dy, dlambda))
}
# Parâmetros e condições iniciais
state <- c(y = 5, lambda = 0.1) # Ajuste inicial para lambda
time <- seq(0, 10, by = 0.01) # T estimado inicialmente
# Resolver o sistema
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
# Processar os resultados
results <- as.data.frame(results)
# Plotar as soluções
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória ótima de y(t)")

Exercício 7.4: Questão 1, página 191
\[ V = \int_0^T -(t^2+u^2) dt\\ \text{Sujeito a:}\\ \dot{y} = u\\ y( 0) = 4\\ y(T) = 5\\ T \text{ Livre} \]
# Função para o sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo
u <- lambda / 2
# Derivadas
dy <- u
dlambda <- 0 # λ é constante
list(c(dy, dlambda))
}
# Função objetivo para ajustar T
objective_function <- function(T_guess) {
# Tempo estimado
time <- seq(0, T_guess, by = 0.01)
# Condições iniciais
state <- c(y = 4, lambda = 1) # Suposição inicial para λ
# Resolver o sistema
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
results <- as.data.frame(results)
# Erro entre y(T) e o valor-alvo (5)
final_y <- tail(results$y, 1)
return(abs(final_y - 5)) # Queremos minimizar esse erro
}
# Otimizar o tempo T para satisfazer y(T) = 5
optimized_T <- optimize(objective_function, interval = c(0.1, 10))$minimum
# Resolver com o tempo otimizado
time <- seq(0, optimized_T, by = 0.01)
state <- c(y = 4, lambda = 1)
results <- ode(y = state, times = time, func = dynamics, parms = NULL)
results <- as.data.frame(results)
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória ótima de y(t)")

Exercício 7.4: Questão 2, página 191
\[ V = \int_0^4 3y dt\\ \text{Sujeito a:}\\ \dot{y} = y+u\\ y( 0) = 5\\ y(4) >= 300\\ u \in [0,2] \]
# Função que define o sistema de equações diferenciais
dynamics <- function(t, state, parameters) {
y <- state[1]
lambda <- state[2]
# Controle ótimo (ajustado com restrição u ∈ [0, 2])
u <- max(0, min(2, -lambda))
# Derivadas
dy <- y + u
dlambda <- -3 - lambda
list(c(dy, dlambda))
}
# Condições iniciais e finais
state <- c(y = 5, lambda = 0) # Estado inicial
time <- seq(0, 4, by = 0.01)
# Resolver o sistema usando 'ode'
results <- data.frame(ode(y = state, times = time, func = dynamics, parms = NULL))
# Garantir que y(T) ≥ 300
if (tail(results$y, 1) < 300) {
stop("Condição final y(4) ≥ 300 não foi atingida. Ajuste os parâmetros.")
}
# Plotar as soluções
plot(results$time, results$y, type = "l", col = "blue", lwd = 2, xlab = "Tempo", ylab = "y(t)", main = "Trajetória Ótima de y(t)")
