Controle Preditivo (MPC) de um Sistema Massa-Mola-Amortecedor

Controle
MPC
Modelagem em espaço de estados, formulação matricial do MPC de horizonte finito e simulação em Python de um sistema massa-mola-amortecedor.
Data de Publicação

7 de outubro de 2026

Introdução

O Controle Preditivo Baseado em Modelo (Model Predictive Control — MPC) otimiza, a cada instante de amostragem, uma sequência de ações de controle sobre um horizonte finito de predição, aplica apenas a primeira ação e repete o processo no instante seguinte (estratégia de horizonte deslizante, ou receding horizon). Esta nota deduz a formulação matricial do MPC irrestrito para sistemas lineares invariantes no tempo (LIT) e a aplica a um sistema clássico: massa-mola-amortecedor.

Modelo do sistema

A equação diferencial de um sistema massa-mola-amortecedor com força de entrada \(u(t)\) e posição \(x_1(t)\) é

\[ m\,\ddot{x}_1 + c\,\dot{x}_1 + k\,x_1 = u(t), \]

com \(m\) a massa, \(c\) o coeficiente de amortecimento e \(k\) a constante elástica. Definindo o estado \(x = [x_1,\; x_2]^\top = [\text{posição},\; \text{velocidade}]^\top\), o modelo em espaço de estados contínuo é

\[ \dot{x} = A_c x + B_c u, \qquad y = C x, \]

\[ A_c = \begin{bmatrix} 0 & 1 \\ -k/m & -c/m \end{bmatrix}, \qquad B_c = \begin{bmatrix} 0 \\ 1/m \end{bmatrix}, \qquad C = \begin{bmatrix} 1 & 0 \end{bmatrix}. \]

Para uso digital, discretiza-se por retentor de ordem zero (ZOH) com período de amostragem \(T_s\), obtendo o modelo discreto

\[ x_{k+1} = A_d x_k + B_d u_k, \qquad y_k = C x_k . \]

Formulação matricial do MPC

Seja \(N\) o horizonte de predição. Empilhando as predições de estado a partir do instante atual \(x_0\) e da sequência de controles futuros \(U = [u_0, u_1, \dots, u_{N-1}]^\top\):

\[ \underbrace{\begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{bmatrix}}_{X} = \underbrace{\begin{bmatrix} A_d \\ A_d^2 \\ \vdots \\ A_d^N \end{bmatrix}}_{S} x_0 + \underbrace{\begin{bmatrix} B_d & 0 & \cdots & 0 \\ A_d B_d & B_d & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ A_d^{N-1} B_d & A_d^{N-2} B_d & \cdots & B_d \end{bmatrix}}_{\Gamma} U . \]

Aplicando a saída (\(\bar{C} = I_N \otimes C\)), as saídas preditas são \(Y = \bar C X = S_y x_0 + G_y U\), com \(S_y = \bar C S\) e \(G_y = \bar C \Gamma\).

A função custo pondera o erro de rastreamento em relação a uma referência \(Y_{ref}\) e o esforço de controle:

\[ J(U) = (Y - Y_{ref})^\top Q (Y - Y_{ref}) + U^\top R\, U, \]

com \(Q \succeq 0\) (peso do rastreamento) e \(R \succ 0\) (peso do esforço de controle). Substituindo \(Y = S_y x_0 + G_y U\) e resolvendo \(\partial J/\partial U = 0\), obtém-se a solução ótima em malha aberta (caso irrestrito):

\[ U^{\star} = \big(G_y^\top Q\, G_y + R\big)^{-1} G_y^\top Q\,(Y_{ref} - S_y x_0). \]

Na estratégia de horizonte deslizante, aplica-se apenas \(u_0^\star\) (primeiro elemento de \(U^\star\)) ao processo, mede-se o novo estado \(x_0\) no instante seguinte e o problema é resolvido novamente. Quando há restrições em \(u\) ou \(x\), o mesmo problema é resolvido por programação quadrática (QP) em vez da forma fechada acima — não tratado nesta nota para manter o foco na formulação central.

Implementação computacional

import numpy as np
from scipy.signal import cont2discrete
import matplotlib.pyplot as plt

# Parâmetros físicos do sistema massa-mola-amortecedor
m, c, k = 1.0, 0.6, 4.0

Ac = np.array([[0, 1],
               [-k / m, -c / m]])
Bc = np.array([[0],
               [1 / m]])
C = np.array([[1, 0]])
D = np.array([[0]])

Ts = 0.05  # período de amostragem [s]
Ad, Bd, Cd, Dd, _ = cont2discrete((Ac, Bc, C, D), Ts, method="zoh")
n, m_in = Ad.shape[0], Bd.shape[1]
def matrizes_predicao(Ad, Bd, C, N):
    n = Ad.shape[0]
    S = np.zeros((N * n, n))
    Gamma = np.zeros((N * n, N * Bd.shape[1]))
    Apow = np.eye(n)
    for i in range(N):
        Apow = Apow @ Ad
        S[i * n:(i + 1) * n, :] = Apow
        for j in range(i + 1):
            power = i - j
            Aij = np.linalg.matrix_power(Ad, power) @ Bd
            Gamma[i * n:(i + 1) * n, j * Bd.shape[1]:(j + 1) * Bd.shape[1]] = Aij

    Cbar = np.kron(np.eye(N), C)
    Sy = Cbar @ S
    Gy = Cbar @ Gamma
    return Sy, Gy

N = 15                      # horizonte de predição
Sy, Gy = matrizes_predicao(Ad, Bd, C, N)

Q = np.eye(N) * 50.0        # peso do rastreamento (saída)
R = np.eye(N) * 0.1         # peso do esforço de controle

# Ganho de MPC pré-calculado (sistema invariante -> reaproveitável a cada passo)
Kmpc = np.linalg.solve(Gy.T @ Q @ Gy + R, Gy.T @ Q)
T_sim = 6.0
n_steps = int(T_sim / Ts)

x = np.array([[0.0], [0.0]])   # estado inicial (repouso, posição 0)
ref = 1.0                       # referência de posição (degrau)

pos_hist = np.zeros(n_steps)
vel_hist = np.zeros(n_steps)
u_hist = np.zeros(n_steps)

for t in range(n_steps):
    Yref = np.ones((N, 1)) * ref
    U = Kmpc @ (Yref - Sy @ x)      # resolve o QP irrestrito (forma fechada)
    u0 = U[0, 0]                    # aplica apenas a primeira ação (receding horizon)

    pos_hist[t] = x[0, 0]
    vel_hist[t] = x[1, 0]
    u_hist[t] = u0

    x = Ad @ x + Bd * u0            # atualiza o estado real do processo

tempo = np.arange(n_steps) * Ts

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6.5, 5.5), sharex=True)
ax1.plot(tempo, pos_hist, label="Posição $x_1$")
ax1.axhline(ref, color="gray", linestyle="--", label="Referência")
ax1.set_ylabel("Posição [m]")
ax1.legend()

ax2.step(tempo, u_hist, where="post", color="tab:red")
ax2.set_xlabel("Tempo [s]")
ax2.set_ylabel("Força $u$ [N]")
plt.tight_layout()
plt.show()

Rastreamento em malha fechada de um degrau de referência via MPC.

Discussão

  • Horizonte \(N\): horizontes curtos tornam o controlador “míope” (pode gerar sobressinal ou instabilidade em sistemas com atraso); horizontes longos aproximam o comportamento do LQR de horizonte infinito, ao custo de mais computação.
  • Pesos \(Q\) e \(R\): aumentar \(Q\) relativamente a \(R\) produz rastreamento mais agressivo (resposta mais rápida, maior esforço de controle); aumentar \(R\) suaviza a ação de controle.
  • Reaproveitamento do ganho: como o sistema é LIT e o problema é irrestrito, o ganho \(K_{mpc} = (G_y^\top Q G_y + R)^{-1} G_y^\top Q\) é calculado uma única vez fora do laço de controle — equivalente, neste caso particular, a um controlador de realimentação de estados variante no tempo apenas pela referência.
  • Extensão a restrições: limites físicos (curso do atuador, força máxima) tornam \(U^\star\) a solução de uma QP restrita, tipicamente resolvida com bibliotecas como cvxpy, qpOASES ou OSQP a cada passo de amostragem.

Referências

  • Camacho, E. F., Bordons, C. (2013). Model Predictive Control. Springer.
  • Rawlings, J. B., Mayne, D. Q., Diehl, M. (2017). Model Predictive Control: Theory, Computation, and Design. Nob Hill Publishing.
  • Maciejowski, J. M. (2002). Predictive Control with Constraints. Prentice Hall.