Identificação de Modelos ARX por Mínimos Quadrados

Da Reta ao ARX

Na aula 01, ajustamos \(y = \theta_1 x + \theta_0 + e\) a dados estáticos, via mínimos quadrados:

\[\boldsymbol{\theta} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{y}\]

Na aula 02, simulamos um sistema dinâmico ARX e coletamos dados \((u[k], y[k])\) excitando-o com um sinal PRBS.

Hoje: mostrar que o mesmo aparato de mínimos quadrados — só que com uma matriz de regressores diferente — permite estimar os parâmetros de um modelo dinâmico.

Modelo ARX

O modelo AutoRegressive with eXogenous input, com \(n_a\) termos autorregressivos e \(n_b\) termos de entrada:

\[y[k] + a_1 y[k-1] + \dots + a_{n_a} y[k-n_a] = b_1 u[k-1] + \dots + b_{n_b} u[k-n_b] + e[k]\]

Ou, na forma polinomial em \(q^{-1}\):

\[A(q)\, y[k] = B(q)\, u[k] + e[k]\]

\[A(q) = 1 + a_1 q^{-1} + \dots + a_{n_a} q^{-n_a}, \qquad B(q) = b_1 q^{-1} + \dots + b_{n_b} q^{-n_b}\]

No nosso exemplo (aula 02): \(n_a = n_b = 1\), com \(a_1=-0{,}7\) e \(b_1=1{,}5\).

Diagrama de Blocos (caso geral)

\(B(q)\) realimenta a entrada; o bloco \(-[A(q)-1]\) realimenta a própria saída — a estrutura recursiva que caracteriza um modelo dinâmico.

Reescrevendo como Regressão Linear

Isolando \(y[k]\):

\[y[k] = \underbrace{-a_1 y[k-1] - \dots - a_{n_a} y[k-n_a] + b_1 u[k-1] + \dots + b_{n_b} u[k-n_b]}_{\boldsymbol{\varphi}[k]^T \boldsymbol{\theta}} + e[k]\]

Definindo o vetor de regressores e o vetor de parâmetros:

\[\boldsymbol{\varphi}[k] = \begin{bmatrix} -y[k-1] \\ \vdots \\ -y[k-n_a] \\ u[k-1] \\ \vdots \\ u[k-n_b] \end{bmatrix}, \qquad \boldsymbol{\theta} = \begin{bmatrix} a_1 \\ \vdots \\ a_{n_a} \\ b_1 \\ \vdots \\ b_{n_b} \end{bmatrix}\]

chega-se à mesma forma da aula 01:

\[y[k] = \boldsymbol{\varphi}[k]^T \boldsymbol{\theta} + e[k]\]

Montando a Matriz de Regressores \(\boldsymbol{\Phi}\)

Empilhando \(\boldsymbol{\varphi}[k]^T\) para \(k = n+1, \dots, N\) (onde \(n=\max(n_a,n_b)\)), obtém-se o mesmo formato matricial da aula 01:

\[\begin{bmatrix} y[n+1] \\ y[n+2] \\ \vdots \\ y[N] \end{bmatrix} = \begin{bmatrix} \boldsymbol{\varphi}[n+1]^T \\ \boldsymbol{\varphi}[n+2]^T \\ \vdots \\ \boldsymbol{\varphi}[N]^T \end{bmatrix} \boldsymbol{\theta} + \begin{bmatrix} e[n+1] \\ e[n+2] \\ \vdots \\ e[N] \end{bmatrix}\]

\[\mathbf{Y} = \boldsymbol{\Phi}\, \boldsymbol{\theta} + \mathbf{E}\]

A diferença para a aula 01 é só o conteúdo de \(\boldsymbol{\Phi}\): antes eram colunas \([x_i, 1]\); agora são saídas e entradas passadas.

Reconstruindo os Dados da Aula 02

import numpy as np
from scipy.signal import max_len_seq

a1_true, b1_true = -0.7, 1.5
sigma_e = 0.3

def simula_arx(u, a1, b1, sigma_e, seed=None):
    rng = np.random.default_rng(seed)
    N = len(u)
    y = np.zeros(N)
    e = rng.normal(0, sigma_e, N)
    for k in range(1, N):
        y[k] = -a1 * y[k-1] + b1 * u[k-1] + e[k]
    return y

n_bits, Ts_bit, A = 7, 2, 1.5
bits, _ = max_len_seq(n_bits)
u_prbs = A * np.repeat(2*bits - 1, Ts_bit)

y_prbs = simula_arx(u_prbs, a1_true, b1_true, sigma_e, seed=42)
N = len(u_prbs)
print(f"N = {N} amostras coletadas")

N = 254 amostras coletadas

Mesma semente (seed=42) e mesmos parâmetros da aula 02 — reprodutibilidade garantida.

Construindo \(\boldsymbol{\Phi}\) e \(\mathbf{Y}\) (caso \(n_a=n_b=1\))

# phi[k] = [-y[k-1], u[k-1]]  para k = 1, ..., N-1
Phi = np.column_stack((-y_prbs[:-1], u_prbs[:-1]))
Y = y_prbs[1:]

print("Phi (5 primeiras linhas):")
print(Phi[:5])
print("Shape de Phi:", Phi.shape)

Phi (5 primeiras linhas):
[[-0.          1.5       ]
 [-1.93800477  1.5       ]
 [-3.8317387   1.5       ]
 [-5.2143865   1.5       ]
 [-5.31476     1.5       ]]
Shape de Phi: (253, 2)

Cada linha de \(\boldsymbol{\Phi}\) contém \([-y[k-1],\ u[k-1]]\); cada elemento de \(\mathbf{Y}\) é o correspondente \(y[k]\).

Estimador de Mínimos Quadrados

Exatamente a mesma equação normal da aula 01, agora aplicada a \(\boldsymbol{\Phi}\):

\[\boldsymbol{\theta} = \arg\min_{\boldsymbol{\theta}} \|\mathbf{Y} - \boldsymbol{\Phi}\boldsymbol{\theta}\|^2 \quad\Longrightarrow\quad \boldsymbol{\theta} = (\boldsymbol{\Phi}^T \boldsymbol{\Phi})^{-1} \boldsymbol{\Phi}^T \mathbf{Y}\]

theta_hat = np.linalg.inv(Phi.T @ Phi) @ Phi.T @ Y
a1_hat, b1_hat = theta_hat

print(f"a1: verdadeiro={a1_true:.4f}  estimado={a1_hat:.4f}")
print(f"b1: verdadeiro={b1_true:.4f}  estimado={b1_hat:.4f}")

a1: verdadeiro=-0.7000  estimado=-0.6949
b1: verdadeiro=1.5000  estimado=1.4994

Com dados gerados por PRBS, a estimativa fica muito próxima dos valores verdadeiros usados na simulação.

Predição Um Passo à Frente

O modelo estimado prevê \(y[k]\) usando o valor real de \(y[k-1]\) (disponível no conjunto de dados):

y_pred_1passo = Phi @ theta_hat

Predição vs. Simulação Livre

Duas formas de usar o modelo estimado são bem diferentes:

• Predição um passo à frente: usa o \(y\) real do instante anterior — sempre disponível durante a estimação, mas “trapaceia” ao validar o modelo;
• Simulação livre: realimenta a própria saída predita, sem nunca olhar o \(y\) real — é o teste mais rigoroso, pois é assim que o modelo seria usado para prever o futuro sem medições.

Simulando em Malha Livre

def simula_livre(u, theta):
    a1, b1 = theta
    N = len(u)
    yhat = np.zeros(N)
    for k in range(1, N):
        yhat[k] = -a1 * yhat[k-1] + b1 * u[k-1]
    return yhat

y_livre = simula_livre(u_prbs, theta_hat)

Mesmo em malha livre — sem nenhuma realimentação de dados reais — o modelo acompanha bem a saída medida: sinal de que \(a_1, b_1\) foram bem estimados.

Qualidade do Ajuste: métrica FIT%

Uma métrica comum em identificação de sistemas normaliza o erro pela variação do próprio sinal:

\[\text{FIT}\% = 100 \left(1 - \frac{\|\mathbf{y} - \hat{\mathbf{y}}\|}{\|\mathbf{y} - \bar{y}\|}\right)\]

def fit_pct(y, y_hat):
    return 100 * (1 - np.linalg.norm(y - y_hat) / np.linalg.norm(y - np.mean(y)))

print(f"FIT% (predição 1 passo): {fit_pct(Y, y_pred_1passo):.2f}%")
print(f"FIT% (simulação livre):  {fit_pct(y_prbs, y_livre):.2f}%")

FIT% (predição 1 passo): 93.08%
FIT% (simulação livre):  90.41%

FIT% = 100% seria ajuste perfeito; FIT% = 0% equivale a prever apenas a média do sinal.

Resíduos

O resíduo é a parte de \(y[k]\) que o modelo não conseguiu explicar:

\[\hat{e}[k] = y[k] - \boldsymbol{\varphi}[k]^T \hat{\boldsymbol{\theta}}\]

residuos = Y - y_pred_1passo

Se o modelo captura bem a dinâmica, o resíduo deve parecer ruído branco — sem padrão temporal.

Checando a “Brancura” dos Resíduos

A autocorrelação mede se o resíduo em \(k\) tem relação com o resíduo em \(k-\ell\). Para ruído branco, ela deve ficar próxima de zero para \(\ell \neq 0\):

def acf(x, max_lag):
    x = x - np.mean(x)
    N = len(x)
    c0 = np.sum(x * x) / N
    return np.array([np.sum(x[:N-lag] * x[lag:]) / N / c0 for lag in range(max_lag + 1)])

r = acf(residuos, max_lag=20)
limite = 1.96 / np.sqrt(len(residuos))   # faixa de confiança ~95% para ruído branco

Os valores permanecem dentro da faixa de 95% (linhas vermelhas) — consistente com resíduo branco, indicando que a estrutura \(n_a=n_b=1\) foi suficiente.

Efeito da Ordem do Modelo

• Ordem baixa demais (\(n_a, n_b\) pequenos): o modelo não tem “graus de liberdade” para capturar a dinâmica → resíduos correlacionados (viés);
• Ordem alta demais: o modelo passa a ajustar o próprio ruído (overfitting) → parâmetros extras com pouco significado físico, sensíveis a novos dados;
• Critérios formais de escolha de ordem (ex.: AIC, critério de informação) penalizam modelos com parâmetros em excesso — tema de aulas futuras.

Na prática: aumente a ordem apenas enquanto o FIT% melhorar de forma expressiva e os resíduos permanecerem não-brancos.

Voltando ao PRBS: e se fosse um Degrau?

u_step = np.concatenate([np.zeros(10), A*np.ones(N-10)])
y_step = simula_arx(u_step, a1_true, b1_true, sigma_e, seed=42)

Phi_step = np.column_stack((-y_step[:-1], u_step[:-1]))
theta_step = np.linalg.inv(Phi_step.T @ Phi_step) @ Phi_step.T @ y_step[1:]

print(f"cond(Phi) com PRBS:  {np.linalg.cond(Phi):.1f}")
print(f"cond(Phi) com degrau: {np.linalg.cond(Phi_step):.1f}")
print(f"theta_hat (degrau): a1={theta_step[0]:.3f}, b1={theta_step[1]:.3f}  |  verdadeiro: a1={a1_true}, b1={b1_true}")

cond(Phi) com PRBS:  2.9
cond(Phi) com degrau: 50.4
theta_hat (degrau): a1=-0.674, b1=1.620  |  verdadeiro: a1=-0.7, b1=1.5

O número de condição de \(\boldsymbol{\Phi}\) é muito maior com o degrau, e a estimativa de \(b_1\) fica visivelmente pior — a falta de persistência de excitação compromete a identificação, mesmo usando o mesmo estimador de mínimos quadrados.

Resumo

• O modelo ARX vira uma regressão linear ao definir \(\boldsymbol{\varphi}[k]\) com saídas e entradas passadas;
• A mesma equação normal \(\boldsymbol{\theta}=(\boldsymbol{\Phi}^T\boldsymbol{\Phi})^{-1}\boldsymbol{\Phi}^T\mathbf{Y}\) da aula 01 se aplica diretamente;
• Predição 1 passo vs. simulação livre avaliam o modelo de formas diferentes — a simulação livre é o teste mais exigente;
• FIT% e autocorrelação dos resíduos ajudam a julgar a qualidade e a ordem do modelo;
• A qualidade da estimativa depende diretamente da persistência de excitação do sinal de entrada (aula 02).

Próximos temas: modelos com ruído colorido (ARMAX), mínimos quadrados recursivo, seleção de ordem.

Ver também

Este mesmo estimador (ARX por mínimos quadrados) aplicado a sistemas reais e generalizado a dinâmicas mais complexas:

• ControleCC2CC — ARX de uma bancada motor-gerador real.
• Projeto-Aeropendulo — mesmo estimador aplicado ao aeropêndulo.
• Controle MPC-FCS de gerador DFIG (ControleDFIG) — ARX multivariável num sistema trifásico, usado como modelo de predição num MPC.
• SINDy e Koopman (DataDrivenControl) — o mesmo estimador de mínimos quadrados, generalizado para dinâmicas não lineares.