Fourier Basis Function¶
Example created by Wilson Rocha Lacerda Junior
This example shows how changing or adding a new basis function could improve the model
In [1]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
from sysidentpy.model_structure_selection import FROLS
from sysidentpy.basis_function import Polynomial, Fourier
from sysidentpy.utils.plotting import plot_results
from sysidentpy.metrics import root_relative_squared_error
np.seterr(all="ignore")
np.random.seed(1)
%matplotlib inline
import numpy as np import matplotlib.pyplot as plt from sysidentpy.model_structure_selection import FROLS from sysidentpy.basis_function import Polynomial, Fourier from sysidentpy.utils.plotting import plot_results from sysidentpy.metrics import root_relative_squared_error np.seterr(all="ignore") np.random.seed(1) %matplotlib inline
Defining the system¶
In [2]:
Copied!
# Simulated system
def system_equation(y, u):
yk = (
(0.2 - 0.75 * np.cos(-y[0] ** 2)) * np.cos(y[0])
- (0.15 + 0.45 * np.cos(-y[0] ** 2)) * np.cos(y[1])
+ np.cos(u[0])
+ 0.2 * u[1]
+ 0.7 * u[0] * u[1]
)
return yk
repetition = 5
random_samples = 200
total_time = repetition * random_samples
n = np.arange(0, total_time)
# Generating input
x = np.random.normal(size=(random_samples,)).repeat(repetition)
_, ax = plt.subplots(figsize=(12, 6))
ax.step(n, x)
ax.set_xlabel("$n$", fontsize=18)
ax.set_ylabel("$x[n]$", fontsize=18)
plt.show()
# Simulated system def system_equation(y, u): yk = ( (0.2 - 0.75 * np.cos(-y[0] ** 2)) * np.cos(y[0]) - (0.15 + 0.45 * np.cos(-y[0] ** 2)) * np.cos(y[1]) + np.cos(u[0]) + 0.2 * u[1] + 0.7 * u[0] * u[1] ) return yk repetition = 5 random_samples = 200 total_time = repetition * random_samples n = np.arange(0, total_time) # Generating input x = np.random.normal(size=(random_samples,)).repeat(repetition) _, ax = plt.subplots(figsize=(12, 6)) ax.step(n, x) ax.set_xlabel("$n$", fontsize=18) ax.set_ylabel("$x[n]$", fontsize=18) plt.show()
Simulate the system¶
In [3]:
Copied!
y = np.empty_like(x)
# Initial Conditions
y0 = [0, 0]
# Simulate it
y[0:2] = y0
for i in range(2, len(y)):
y[i] = system_equation(
[y[i - 1], y[i - 2]], [x[i - 1], x[i - 2]]
) + np.random.normal(scale=0.1)
# Plot
_, ax = plt.subplots(figsize=(12, 6))
ax.plot(n, y)
ax.set_xlabel("$n$", fontsize=18)
ax.set_ylabel("$y[n]$", fontsize=18)
ax.grid()
plt.show()
y = np.empty_like(x) # Initial Conditions y0 = [0, 0] # Simulate it y[0:2] = y0 for i in range(2, len(y)): y[i] = system_equation( [y[i - 1], y[i - 2]], [x[i - 1], x[i - 2]] ) + np.random.normal(scale=0.1) # Plot _, ax = plt.subplots(figsize=(12, 6)) ax.plot(n, y) ax.set_xlabel("$n$", fontsize=18) ax.set_ylabel("$y[n]$", fontsize=18) ax.grid() plt.show()
Adding noise to the system¶
In [4]:
Copied!
# Noise free data
ynoise_free = y.copy()
# Generate noise
v = np.random.normal(scale=0.5, size=y.shape)
# Data corrupted with noise
ynoisy = ynoise_free + v
# Plot
_, ax = plt.subplots(figsize=(14, 8))
ax.plot(n, ynoise_free, label="Noise-free data")
ax.plot(n, ynoisy, label="Corrupted data")
ax.set_xlabel("$n$", fontsize=18)
ax.set_ylabel("$y[n]$", fontsize=18)
ax.legend(fontsize=18)
plt.show()
# Noise free data ynoise_free = y.copy() # Generate noise v = np.random.normal(scale=0.5, size=y.shape) # Data corrupted with noise ynoisy = ynoise_free + v # Plot _, ax = plt.subplots(figsize=(14, 8)) ax.plot(n, ynoise_free, label="Noise-free data") ax.plot(n, ynoisy, label="Corrupted data") ax.set_xlabel("$n$", fontsize=18) ax.set_ylabel("$y[n]$", fontsize=18) ax.legend(fontsize=18) plt.show()
Generating training and test data¶
In [5]:
Copied!
n_train = 700
# Identification data
y_train = ynoisy[:n_train].reshape(-1, 1)
x_train = x[:n_train].reshape(-1, 1)
# Validation data
y_test = ynoise_free[n_train:].reshape(-1, 1)
x_test = x[n_train:].reshape(-1, 1)
n_train = 700 # Identification data y_train = ynoisy[:n_train].reshape(-1, 1) x_train = x[:n_train].reshape(-1, 1) # Validation data y_test = ynoise_free[n_train:].reshape(-1, 1) x_test = x[n_train:].reshape(-1, 1)
Polynomial Basis Function¶
As you can see bellow, using only the polynomial basis function with the following parameters results in a divergent model
In [6]:
Copied!
basis_function = Polynomial(degree=2)
sysidentpy = FROLS(
order_selection=True,
n_info_values=70,
xlag=2,
ylag=2,
basis_function=basis_function,
model_type="NARMAX",
)
sysidentpy.fit(X=x_train, y=y_train)
yhat = sysidentpy.predict(X=x_test, y=y_test)
frols_loss = root_relative_squared_error(
y_test[sysidentpy.max_lag :], yhat[sysidentpy.max_lag :]
)
print(frols_loss)
plot_results(y=y_test[sysidentpy.max_lag :], yhat=yhat[sysidentpy.max_lag :])
basis_function = Polynomial(degree=2) sysidentpy = FROLS( order_selection=True, n_info_values=70, xlag=2, ylag=2, basis_function=basis_function, model_type="NARMAX", ) sysidentpy.fit(X=x_train, y=y_train) yhat = sysidentpy.predict(X=x_test, y=y_test) frols_loss = root_relative_squared_error( y_test[sysidentpy.max_lag :], yhat[sysidentpy.max_lag :] ) print(frols_loss) plot_results(y=y_test[sysidentpy.max_lag :], yhat=yhat[sysidentpy.max_lag :])
c:\Users\wilso\Desktop\projects\GitHub\sysidentpy\sysidentpy\utils\deprecation.py:37: FutureWarning: Passing a string to define the estimator will rise an error in v0.4.0. You'll have to use FROLS(estimator=LeastSquares()) instead. The only change is that you'll have to define the estimator first instead of passing a string like 'least_squares'. This change will make easier to implement new estimators and it'll improve code readability. warnings.warn(message, FutureWarning) c:\Users\wilso\Desktop\projects\GitHub\sysidentpy\sysidentpy\model_structure_selection\forward_regression_orthogonal_least_squares.py:569: UserWarning: n_info_values is greater than the maximum number of all regressors space considering the chosen y_lag, u_lag, and non_degree. We set as 15 self.info_values = self.information_criterion(reg_matrix, y)
nan
Ensembling a Fourier Basis Function¶
In this case, adding the Fourier Basis Function solves the problem and returns a model capable to predict the defined system
In [7]:
Copied!
basis_function = Fourier(degree=2, n=2, p=2 * np.pi, ensemble=True)
sysidentpy = FROLS(
order_selection=True,
n_info_values=70,
xlag=2,
ylag=2, # the lags for all models will be 13
basis_function=basis_function,
model_type="NARMAX",
)
sysidentpy.fit(X=x_train, y=y_train)
yhat = sysidentpy.predict(X=x_test, y=y_test)
frols_loss = root_relative_squared_error(
y_test[sysidentpy.max_lag :], yhat[sysidentpy.max_lag :]
)
print(frols_loss)
plot_results(y=y_test[sysidentpy.max_lag :], yhat=yhat[sysidentpy.max_lag :])
basis_function = Fourier(degree=2, n=2, p=2 * np.pi, ensemble=True) sysidentpy = FROLS( order_selection=True, n_info_values=70, xlag=2, ylag=2, # the lags for all models will be 13 basis_function=basis_function, model_type="NARMAX", ) sysidentpy.fit(X=x_train, y=y_train) yhat = sysidentpy.predict(X=x_test, y=y_test) frols_loss = root_relative_squared_error( y_test[sysidentpy.max_lag :], yhat[sysidentpy.max_lag :] ) print(frols_loss) plot_results(y=y_test[sysidentpy.max_lag :], yhat=yhat[sysidentpy.max_lag :])
c:\Users\wilso\Desktop\projects\GitHub\sysidentpy\sysidentpy\utils\deprecation.py:37: FutureWarning: Passing a string to define the estimator will rise an error in v0.4.0. You'll have to use FROLS(estimator=LeastSquares()) instead. The only change is that you'll have to define the estimator first instead of passing a string like 'least_squares'. This change will make easier to implement new estimators and it'll improve code readability. warnings.warn(message, FutureWarning)
0.3742244715879481
Important¶
Currently you can't get the model representation using sysidentpy.regressor_code for Fourier NARX models. Actually, you can use the method, but the representation is not accurate because we don't make clear what are the regressors related to the polynomial or related to the fourier basis function. This is a improvement to be done in future updates!