Sistema Wiener–Hammerstein¶

Nota: O exemplo mostrado neste notebook é retirado do livro Nonlinear System Identification and Forecasting: Theory and Practice with SysIdentPy.

O conteúdo da descrição deriva principalmente do site do benchmark - Nonlinear Benchmark e do artigo associado - Wiener-Hammerstein benchmark with process noise. Para uma descrição detalhada, os leitores são encaminhados às referências vinculadas.

O site de benchmarks não lineares representa uma contribuição significativa para a comunidade de identificação de sistemas e aprendizado de máquina. Os usuários são encorajados a explorar todos os artigos referenciados no site.

Este benchmark foca em um circuito eletrônico Wiener-Hammerstein onde o ruído de processo desempenha um papel significativo na distorção do sinal de saída.

A estrutura Wiener-Hammerstein é um sistema orientado a blocos bem conhecido que contém uma não linearidade estática intercalada entre dois blocos Lineares Invariantes no Tempo (LTI) (Figura 2). Este arranjo apresenta um problema de identificação desafiador devido à presença desses blocos LTI.

Figura 2: o sistema Wiener-Hammerstein

Na Figura 2, o sistema Wiener-Hammerstein é ilustrado com ruído de processo \(e_x(t)\) entrando antes da não linearidade estática \(f(x)\), intercalado entre blocos LTI representados por \(R(s)\) e \(S(s)\) na entrada e saída, respectivamente. Além disso, pequenas fontes de ruído desprezíveis \(e_u(t)\) e \(e_y(t)\) afetam os canais de medição. Os sinais de entrada e saída medidos são denotados como \(u_m(t)\) e \(y_m(t)\).

O primeiro bloco LTI \(R(s)\) é efetivamente modelado como um filtro passa-baixa de terceira ordem. O segundo subsistema LTI \(S(s)\) é configurado como um filtro Chebyshev inverso com atenuação de banda de parada de \(40 dB\) e frequência de corte de \(5 kHz\). Notavelmente, \(S(s)\) inclui um zero de transmissão dentro da faixa de frequência operacional, complicando sua inversão.

A não linearidade estática \(f(x)\) é implementada usando uma rede de diodo-resistor, resultando em não linearidade de saturação. O ruído de processo \(e_x(t)\) é introduzido como ruído gaussiano branco filtrado, gerado a partir de um filtro Butterworth passa-baixa de terceira ordem em tempo discreto seguido por zero-order hold e filtragem de reconstrução passa-baixa analógica com corte de \(20 kHz\).

As fontes de ruído de medição \(e_u(t)\) e \(e_y(t)\) são mínimas comparadas a \(e_x(t)\). As entradas do sistema e o ruído de processo são gerados usando um Gerador de Forma de Onda Arbitrária (AWG), especificamente o Agilent/HP E1445A, amostrando a \(78125 Hz\), sincronizado com um sistema de aquisição (Agilent/HP E1430A) para garantir coerência de fase e prevenir erros de vazamento. O buffering entre as placas de aquisição e as entradas e saídas do sistema minimiza a distorção do equipamento de medição.

O benchmark fornece dois sinais de teste padrão através do site de benchmarking: um multisine de fase aleatória e um sinal de varredura senoidal. Ambos os sinais têm um valor \(rms\) de \(0.71 Vrms\) e cobrem frequências de DC a \(15 kHz\) (excluindo DC). A varredura senoidal abrange esta faixa de frequência a uma taxa de \(4.29 MHz/min\). Estes conjuntos de teste servem como alvos para avaliar o desempenho do modelo, enfatizando representação precisa sob condições variadas.

O benchmark Wiener-Hammerstein destaca três desafios principais de identificação de sistemas não lineares:

Ruído de Processo: Significativo no sistema, influenciando a fidelidade da saída.
Não Linearidade Estática: Indiretamente acessível a partir de dados medidos, apresentando desafios de identificação.
Dinâmicas de Saída: Inversão complexa devido à presença de zero de transmissão em \(S(s)\).

O objetivo deste benchmark é desenvolver e validar modelos robustos usando dados de estimação separados, garantindo caracterização precisa do comportamento do sistema Wiener-Hammerstein.

Pacotes Necessários e Versões¶

Este estudo de caso foi verificado com o SysIdentPy 0.9.0 no Python 3.12.12 e nonlinear-benchmarks==1.0.1. Instale explicitamente o checkout do repositório e o carregador oficial do benchmark:

python -m pip install -e .
python -m pip install nonlinear-benchmarks==1.0.1

Use um ambiente virtual para isolar o carregador opcional. Os resultados numéricos devem ser recalculados se o ambiente ou a configuração do modelo forem alterados.

Configuração do SysIdentPy¶

Nesta seção, demonstraremos a aplicação do SysIdentPy ao dataset do sistema Wiener-Hammerstein. O código a seguir guiará você através do processo de carregamento do dataset, configuração dos parâmetros do SysIdentPy e construção de um modelo para o sistema Wiener-Hammerstein.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sysidentpy.model_structure_selection import FROLS, AOLS, MetaMSS
from sysidentpy.basis_function import Polynomial, Fourier
from sysidentpy.utils.display_results import results
from sysidentpy.parameter_estimation import (
    LeastSquares,
    BoundedVariableLeastSquares,
    NonNegativeLeastSquares,
    LeastSquaresMinimalResidual,
)

from sysidentpy.metrics import root_mean_squared_error
from sysidentpy.utils.plotting import plot_results

import nonlinear_benchmarks

train_val, test = nonlinear_benchmarks.WienerHammerBenchMark(atleast_2d=True)
x_train, y_train = train_val
x_test, y_test = test

Usamos o pacote nonlinear_benchmarks para carregar os dados. O usuário é encaminhado à documentação do pacote para verificar os detalhes de como usá-lo.

O gráfico a seguir detalha os dados de treinamento e teste do experimento.

plot_n = 800

plt.figure(figsize=(15, 4))
plt.plot(x_train[:plot_n])
plt.plot(y_train[:plot_n])
plt.title("Experiment: training data")
plt.legend(["x_train", "y_train"])
plt.show()

plt.figure(figsize=(15, 4))
plt.plot(x_test[:plot_n])
plt.plot(y_test[:plot_n])
plt.title("Experiment: testing data")
plt.legend(["x_test", "y_test"])
plt.show()

O benchmark externo fornece um contexto útil para o experimento. Uma comparação direta, porém, exige a mesma divisão dos dados, janela de inicialização e normalização; essas condições são explicitadas abaixo antes de qualquer comparação.

Resultados estado-da-arte apresentados no artigo de benchmarking. Nesta seção estamos trabalhando apenas com os resultados Wiener-Hammerstein, que são apresentados na coluna \(W-H\).

Resultados¶

Começaremos com uma configuração básica do FROLS usando uma função de base polinomial com grau igual a 2. O xlag e ylag são definidos como \(7\) neste primeiro exemplo. Como o dataset é consideravelmente grande, começaremos com n_info_values=50. Isso significa que o algoritmo FROLS não incluirá todos os regressores ao calcular os critérios de informação usados para determinar a ordem do modelo. Embora esta abordagem possa resultar em um modelo sub-ótimo, é um ponto de partida razoável para nossa primeira tentativa.

n = test.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=7,
    ylag=7,
    basis_function=basis_function,
    estimator=LeastSquares(unbiased=False),
    n_info_values=50,
)

model.fit(X=x_train, y=y_train)
if model.max_lag > n:
    raise ValueError("The model lag exceeds the benchmark initialization window.")
start = n - model.max_lag
yhat = model.predict(X=x_test[start:], y=y_test[start:n])
yhat = yhat[model.max_lag :]
rmse = root_mean_squared_error(y_test[n:], yhat)
nrmse = rmse / np.std(y_test[n:])
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=1000,
    title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(nrmse, 4)}",
)

A primeira configuração produz RMSE \(0.020007\) e NRMSE \(0.082029\). Começamos com xlag=ylag=7 para estabelecer uma referência compacta. O artigo de benchmarking usa memórias mais longas em alguns modelos; por isso, a próxima configuração define xlag=ylag=10.

x_train, y_train = train_val
x_test, y_test = test

n = test.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=10,
    ylag=10,
    basis_function=basis_function,
    estimator=LeastSquares(unbiased=False),
    n_info_values=50,
)

model.fit(X=x_train, y=y_train)
if model.max_lag > n:
    raise ValueError("The model lag exceeds the benchmark initialization window.")
start = n - model.max_lag
yhat = model.predict(X=x_test[start:], y=y_test[start:n])
yhat = yhat[model.max_lag :]
rmse = root_mean_squared_error(y_test[n:], yhat)
nrmse = rmse / np.std(y_test[n:])
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=1000,
    title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(nrmse, 4)}",
)

A configuração com 10 lags melhora o resultado para RMSE \(0.015202\) e NRMSE \(0.062328\). Por enquanto, não estamos otimizando a complexidade do modelo. Ainda assim, o traçado do critério de informação mostra que o modelo com 50 regressores muda pouco após várias das últimas inclusões.

plt.plot(model.info_values)

Então, o que acontece se definirmos um modelo com metade dos regressores?

x_train, y_train = train_val
x_test, y_test = test

n = test.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=10,
    ylag=10,
    basis_function=basis_function,
    estimator=LeastSquares(unbiased=False),
    n_info_values=50,
    n_terms=25,
    order_selection=False,
)

model.fit(X=x_train, y=y_train)
if model.max_lag > n:
    raise ValueError("The model lag exceeds the benchmark initialization window.")
start = n - model.max_lag
yhat = model.predict(X=x_test[start:], y=y_test[start:n])
yhat = yhat[model.max_lag :]
rmse = root_mean_squared_error(y_test[n:], yhat)
nrmse = rmse / np.std(y_test[n:])
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=1000,
    title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(nrmse, 4)}",
)

O modelo fixado em 25 termos produz RMSE \(0.018809\) e NRMSE \(0.077117\). Ele é mais compacto que o modelo de 50 termos selecionado automaticamente, ao custo de um erro maior. Os três resultados atuais são resumidos a seguir.

Configuração	RMSE	NRMSE
7 lags, ordem automática	0.020007	0.082029
10 lags, ordem automática	0.015202	0.062328
10 lags, 25 termos fixos	0.018809	0.077117

As figuras históricas armazenavam quatro casas decimais, e os valores reproduzidos permanecem iguais nessa precisão. Tabelas publicadas podem usar outra normalização ou divisão dos dados, portanto não são tratadas como uma comparação direta. Quem quiser investigar alternativas baseadas em modelos de estado profundos pode explorar o pacote deepsysid.

Esta configuração básica pode servir como ponto de partida para os usuários desenvolverem modelos ainda melhores usando o SysIdentPy. Experimente!