Alpha C. Chiang. Elements Of Dynamic Optimization. 1992: Respostas do Capítulo 7 em R

R

 

Exemplo 1, página 171.

maxV=∫T0−(1−u2)1/2dtSujeito a:˙y=uy(0)=Ay(T) livreA, 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.

maxV=∫20(2y−3u)dtSujeito a:˙y=y+uy(0)=4y(2) livreu(t)∈[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

maxV=∫403ydtSujeito a:˙y=y+uy(0)=5y(4) livreu(t)∈[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

maxV=∫20(y−u2)dtSujeito a:˙y=uy(0)=0u(t) 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∫10−12(y2+u2)dtSujeito a:˙y=u−yy(0)=1u(t) 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

maxV=∫10−u2dtSujeito a:˙y=y+uy(0)=1y(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

maxV=∫T0−1dtSujeito a:˙y=y+uy(0)=5y(T)=11T Livreu(t)∈[−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=∫T0−(t2+u2)dtSujeito a:˙y=uy(0)=4y(T)=5T 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=∫403ydtSujeito a:˙y=y+uy(0)=5y(4)>=300u∈[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)")