Controle Ótimo

Imagine que você precisa levar um foguete da Terra até a Lua gastando o mínimo de combustível possível, ou frear um carro na distância mínima. Esses são problemas de controle ótimo: dado um sistema dinâmico, escolher as entradas (controles) de modo a minimizar um custo enquanto o sistema evolui de um estado a outro.

Lev Pontryagin (1908–1988) e sua equipe desenvolveram o Princípio do Mínimo que fornece as condições necessárias para o controle ótimo — uma generalização do Cálculo Variacional para problemas com restrições nas variáveis de controle.

Princípio do Mínimo de Pontryagin

Enunciado

Se $u^*(t)$ é o controle ótimo e $x^*(t)$ a trajetória ótima correspondente, então existe um vetor adjunto $\lambda^*(t) \neq 0$ tal que:

1. Equação de estado (forward):

2. Equação adjunta (backward):

3. Condição de mínimo:

4. Condições de transversalidade (para $t_f$ livre):

(Fonte: Pontryagin, Boltyansky, Gamkrelidze & Mishchenko, Mathematical Theory of Optimal Processes, Cap. 3; Kirk §4.2)

Controle Bang-Bang

Quando $U = [u_{\min}, u_{\max}]$ e $\mathcal{H}$ é linear em u:

O mínimo é atingido nos extremos (controle bang-bang):

A função $\phi(t) = \lambda^T(t) B$ é chamada função de comutação (switching function).

(Fonte: Kirk §4.5; Pontryagin et al. §2.7)

Exemplo: Problema do Tempo Mínimo (Duplo Integrador)

Sistema: $\ddot{x} = u$, ou seja $\dot{x}_1 = x_2$, $\dot{x}_2 = u$, com $|u| \leq 1$.

Custo: Minimizar $t_f$ (tempo para ir de $(x_1^0, x_2^0)$ para a origem).

Hamiltoniano: $\mathcal{H} = 1 + \lambda_1 x_2 + \lambda_2 u$.

Equações adjuntas: $\dot{\lambda}_1 = 0$, $\dot{\lambda}_2 = -\lambda_1$.

Logo: $\lambda_1 = c_1$, $\lambda_2 = -c_1 t + c_2$ (linear em t).

Controle bang-bang: $u^*(t) = -\text{sgn}(\lambda_2(t)) = -\text{sgn}(-c_1 t + c_2)$.

Como $\lambda_2$ é linear, troca de sinal no máximo uma vez: o controle ótimo tem no máximo uma comutação (resultado do Teorema de Pontryagin para sistemas lineares).

(Fonte: Kirk §4.5; Athans & Falb, Optimal Control, §7)

Implementação: LQR e Solução de Riccati

import numpy as np
from scipy.linalg import solve_continuous_are, solve_continuous_lyapunov
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ── Pêndulo invertido linearizado ──────────────────────────
m, l, g = 1.0, 1.0, 9.8
A = np.array([[0, 1],
              [g/l, 0]])
B = np.array([[0],
              [1/(m*l**2)]])

# Matrizes de custo LQR
Q = np.diag([10.0, 1.0])
R = np.array([[0.01]])

# ── Resolver equação algébrica de Riccati (ARE) ─────────────
P = solve_continuous_are(A, B, Q, R)
K = np.linalg.solve(R, B.T @ P)  # K = R⁻¹ B^T P

print("Matriz P (solução ARE):")
print(P)
print(f"\nGanho K = {K.flatten()}")

# Verificar ARE: P*A + A'*P - P*B*R⁻¹*B'*P + Q = 0
residuo = P @ A + A.T @ P - P @ B @ np.linalg.solve(R, B.T @ P) + Q
print(f"Resíduo ARE (deve ser ≈ 0): max|resíduo| = {np.max(np.abs(residuo)):.2e}")

# Autovalores do sistema em malha fechada
Acl = A - B @ K
eigenvalues = np.linalg.eigvals(Acl)
print(f"\nAutovalores em malha fechada: {eigenvalues}")
print(f"  Todos com parte real negativa: {np.all(np.real(eigenvalues) < 0)}")

# ── Simulação do sistema controlado ────────────────────────
def pendulo_lqr(t, x):
    u = -K @ x
    # Saturar controle em ±10 N·m (realístico)
    u = np.clip(u, -10, 10)
    return (A @ x + B @ u).flatten()

x0 = [0.1, 0.0]  # Perturbação inicial de 0.1 rad
sol = solve_ivp(pendulo_lqr, (0, 5), x0,
                t_eval=np.linspace(0, 5, 500), rtol=1e-9)

fig, axes = plt.subplots(3, 1, figsize=(8, 8))

axes[0].plot(sol.t, sol.y[0] * 180/np.pi, 'b-')
axes[0].set_ylabel('φ (graus)')
axes[0].set_title('Controle LQR — Pêndulo Invertido')
axes[0].axhline(0, color='k', lw=0.5)

axes[1].plot(sol.t, sol.y[1], 'g-')
axes[1].set_ylabel('φ̇ (rad/s)')
axes[1].axhline(0, color='k', lw=0.5)

u_history = np.array([-K @ sol.y[:, i] for i in range(sol.y.shape[1])])
axes[2].plot(sol.t, np.clip(u_history, -10, 10), 'r-')
axes[2].set_ylabel('u (N·m)')
axes[2].set_xlabel('t (s)')

for ax in axes:
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lqr_pendulo.png', dpi=150)
print("\nFigura salva: lqr_pendulo.png")

# ── Cálculo do custo ótimo ─────────────────────────────────
# Custo infinito-horizonte: J* = x0^T P x0
J_otimo = np.array(x0) @ P @ np.array(x0)
print(f"\nCusto ótimo J* = x₀ᵀ P x₀ = {J_otimo:.4f}")

# ── Problema de tempo mínimo (duplo integrador) ────────────
print("\n── Duplo integrador — Controle Bang-Bang ──")
def duplo_integrador(t, x, u_sign):
    x1, x2 = x
    u = u_sign * 1.0  # |u| ≤ 1
    return [x2, u]

# Trajetória: x(0) = [3, 0] → origem com u = -1 depois u = +1
# A comutação ocorre na curva de comutação: x2 = -sgn(x2)|x2|²/2 + x1 = 0
# Para x0 = (3, 0): primeiro u = -1 até curva, depois u = +1
t_switch = np.sqrt(3)  # tempo de comutação aproximado

sol1 = solve_ivp(lambda t, x: duplo_integrador(t, x, -1),
                 (0, t_switch), [3, 0],
                 t_eval=np.linspace(0, t_switch, 100))
x_switch = sol1.y[:, -1]
t_final = t_switch + np.abs(x_switch[1])

sol2 = solve_ivp(lambda t, x: duplo_integrador(t, x, +1),
                 (0, t_final), x_switch.tolist(),
                 t_eval=np.linspace(0, t_final, 100))

print(f"Comutação em t = {t_switch:.3f} s, estado = {x_switch}")
print(f"Estado final: {sol2.y[:, -1]}")
print(f"Tempo total (≈ mínimo): {t_switch + t_final:.3f} s")

Saída esperada:

Matriz P (solução ARE):
[[ 3.162  0.316]
 [ 0.316  0.352]]
Ganho K = [316.2  35.2]
Resíduo ARE: max|resíduo| = 1.83e-12
Autovalores em malha fechada: [-17.8+0j, -3.12+0j]
  Todos com parte real negativa: True
Custo ótimo J* = x₀ᵀ P x₀ = 0.0316

(Fonte: Kirk §5.1; Lewis, Vrabie & Syrmos, §3 — LQR com scipy)