Initial commit

Author: forgi86

README.md

@@ -1 +1,55 @@
-# sysid-transfer-functions-pytorch
+# Transfer functions and deep learning with dynoNet : new applications in system
+identification
+ 
+
+This repository contains the Python code to reproduce the results of the paper "Transfer functions and deep learning with dynoNet : new applications in system
+identification" by Marco Forgione and Dario Piga.
+
+We describe the linear dynamical operator as a differentiable layer compatible with back-propagation-based training. 
+The operator is parametrized as a rational transfer function and thus can represent an infinite impulse response (IIR)
+filtering operation, as opposed to the Convolutional layer of 1D-CNNs that is equivalent to finite impulse response (FIR) filtering.
+
+In the dynoNet architecture (already introduced [here](https://github.com/forgi86/dynonet)), linear dynamical operators are combined with static (i.e., memoryless) non-linearities which can be either elementary
+activation functions applied channel-wise; fully connected feed-forward neural networks; or other differentiable operators. 
+
+In this work, we show how to non-standard learning problems may be tackled using the differentiable 
+transfer function block, namely:
+
+* Learning with quantized measurements
+* Learning in the presence of colored noise
+
+# Folders:
+* [torchid](torchid_nb):  PyTorch implementation of the linear dynamical operator (aka G-block in the paper) used in dynoNet
+* [examples](examples): examples using dynoNet for system identification 
+* [util](util): definition of metrics R-square, RMSE, fit index 
+
+Two [examples](examples) discussed in the paper are:
+
+* [Parallel Wiener-Hammerstein](examples/ParWH): A circuit with Wiener-Hammerstein behavior. Experimental dataset from http://www.nonlinearbenchmark.org
+* [BW](examples/BW): Bouc-Wen. A nonlinear dynamical system describing hysteretic effects in mechanical engineering. Experimental dataset from http://www.nonlinearbenchmark.org
+
+
+For the [WH2009](examples/WH2009) example, the main scripts are:
+
+ *  ``WH2009_train.py``: Training of the dynoNet model
+ *  ``WH2009_test.py``: Evaluation of the dynoNet model on the test dataset,  computation of metrics.
+  
+Similar scripts are provided for the other examples.
+
+NOTE: the original data sets are not included in this project. They have to be manually downloaded from
+http://www.nonlinearbenchmark.org and copied in the data sub-folder of the example.
+# Software requirements:
+Simulations were performed on a Python 3.7 conda environment with
+
+ * numpy
+ * scipy
+ * matplotlib
+ * pandas
+ * pytorch (version 1.4)
+ 
+These dependencies may be installed through the commands:
+
+```
+conda install numpy scipy pandas matplotlib
+conda install pytorch torchvision cudatoolkit=10.2 -c pytorch
+```

examples/ParWH/__init__.py

@@ -0,0 +1,63 @@
+import torch
+from torchid_nb.module.lti import MimoLinearDynamicalOperator, SisoLinearDynamicalOperator
+from torchid_nb.module.static import MimoStaticNonLinearity, MimoChannelWiseNonLinearity
+
+
+class ParallelWHNet(torch.nn.Module):
+    def __init__(self):
+        super(ParallelWHNet, self).__init__()
+        self.nb_1 = 12
+        self.na_1 = 12
+        self.nb_2 = 13
+        self.na_2 = 12
+        self.G1 = MimoLinearDynamicalOperator(1, 2, n_b=self.nb_1, n_a=self.na_1, n_k=1)
+        self.F_nl = MimoChannelWiseNonLinearity(2, n_hidden=10)
+        #self.F_nl = MimoStaticNonLinearity(2, 2, n_hidden=10)
+        self.G2 = MimoLinearDynamicalOperator(2, 1, n_b=self.nb_2, n_a=self.na_2, n_k=0)
+        #self.G3 = SisoLinearDynamicalOperator(n_b=3, n_a=3, n_k=1)
+
+    def forward(self, u):
+        y1_lin = self.G1(u)
+        y1_nl = self.F_nl(y1_lin)  # B, T, C1
+        y2_lin = self.G2(y1_nl)  # B, T, C2
+
+        return y2_lin #+ self.G3(u)
+
+
+class ParallelWHNetVar(torch.nn.Module):
+    def __init__(self):
+        super(ParallelWHNetVar, self).__init__()
+        self.nb_1 = 3
+        self.na_1 = 3
+        self.nb_2 = 3
+        self.na_2 = 3
+        self.G1 = MimoLinearDynamicalOperator(1, 16, n_b=self.nb_1, n_a=self.na_1, n_k=1)
+        self.F_nl = MimoStaticNonLinearity(16, 16) #MimoChannelWiseNonLinearity(16, n_hidden=10)
+        self.G2 = MimoLinearDynamicalOperator(16, 1, n_b=self.nb_2, n_a=self.na_2, n_k=1)
+
+    def forward(self, u):
+        y1_lin = self.G1(u)
+        y1_nl = self.F_nl(y1_lin)  # B, T, C1
+        y2_lin = self.G2(y1_nl)  # B, T, C2
+
+        return y2_lin
+
+
+class ParallelWHResNet(torch.nn.Module):
+    def __init__(self):
+        super(ParallelWHResNet, self).__init__()
+        self.nb_1 = 4
+        self.na_1 = 4
+        self.nb_2 = 4
+        self.na_2 = 4
+        self.G1 = MimoLinearDynamicalOperator(1, 2, n_b=self.nb_1, n_a=self.na_1, n_k=1)
+        self.F_nl = MimoChannelWiseNonLinearity(2, n_hidden=10)
+        self.G2 = MimoLinearDynamicalOperator(2, 1, n_b=self.nb_2, n_a=self.na_2, n_k=1)
+        self.G3 = SisoLinearDynamicalOperator(n_b=6, n_a=6, n_k=1)
+
+    def forward(self, u):
+        y1_lin = self.G1(u)
+        y1_nl = self.F_nl(y1_lin)  # B, T, C1
+        y2_lin = self.G2(y1_nl)  # B, T, C2
+
+        return y2_lin + self.G3(u)

examples/ParWH/models.py

@@ -0,0 +1,41 @@
+import pandas as pd
+import numpy as np
+import os
+import matplotlib.pyplot as plt
+
+if __name__ == '__main__':
+
+    N = 16384  # number of samples per period
+    M = 20  # number of random phase multisine realizations
+    P = 2  # number of periods
+    nAmp = 5 #  number of different amplitudes
+
+    # Column names in the dataset
+    COL_F = ['fs']
+    TAG_U = 'u'
+    TAG_Y = 'y'
+
+    # Load dataset
+    #df_X = pd.read_csv(os.path.join("data", "WH_CombinedZeroMultisineSinesweep.csv"))
+    df_X = pd.read_csv(os.path.join("data", "ParWHData_Estimation_Level2.csv"))
+    df_X.columns = ['amplitude', 'fs', 'lines'] + [TAG_U + str(i) for i in range(M)] + [TAG_Y + str(i) for i in range(M)] + ['?']
+
+    # Extract data
+    y = np.array(df_X['y0'], dtype=np.float32)
+    u = np.array(df_X['u0'], dtype=np.float32)
+    fs = np.array(df_X[COL_F].iloc[0], dtype = np.float32)
+    N = y.size
+    ts = 1/fs
+    t = np.arange(N)*ts
+
+
+    # In[Plot]
+    fig, ax = plt.subplots(2, 1, sharex=True)
+    ax[0].plot(t, y, 'k', label="$y$")
+    ax[0].legend()
+    ax[0].grid()
+
+    ax[1].plot(t, u, 'k', label="$u$")
+    ax[1].legend()
+    ax[1].grid()
+

examples/ParWH/parWH_plot_signals.py

@@ -0,0 +1,125 @@
+import pandas as pd
+import numpy as np
+import os
+import matplotlib.pyplot as plt
+import torch
+import torch.nn as nn
+import control
+from torchid_nb.module.lti import MimoLinearDynamicalOperator
+from torchid_nb.module.static import MimoStaticNonLinearity
+import util.metrics
+from examples.ParWH.models import ParallelWHNet
+
+
+if __name__ == '__main__':
+
+    model_name = "PWH_quant"
+
+    # Dataset constants
+    amplitudes = 5  #  number of different amplitudes
+    realizations = 20  # number of random phase multisine realizations
+    samp_per_period = 16384  # number of samples per period
+    n_skip = 1000
+    periods = 1  # number of periods
+    seq_len = samp_per_period * periods  # data points per realization
+
+    # Column names in the dataset
+    TAG_U = 'u'
+    TAG_Y = 'y'
+
+    test_signal = "100mV" #"ramp" #ramp"#"320mV" #"1000mV"#"ramp"
+
+
+    # In[Load dataset]
+
+    dict_test = {"100mV": 0, "320mV": 1, "550mV": 2, "775mV": 3, "1000mV": 4, "ramp": 5}
+    dataset_list_level = ['ParWHData_Validation_Level' + str(i) for i in range(1, amplitudes + 1)]
+    dataset_list = dataset_list_level + ['ParWHData_ValidationArrow']
+
+    df_X_lst = []
+    for dataset_name in dataset_list:
+        dataset_filename = dataset_name + '.csv'
+        df_Xi = pd.read_csv(os.path.join("data", dataset_filename))
+        df_X_lst.append(df_Xi)
+
+
+    df_X = df_X_lst[dict_test[test_signal]]  # first
+
+    # Extract data
+    y_meas = np.array(df_X['y'], dtype=np.float32)
+    u = np.array(df_X['u'], dtype=np.float32)
+    fs = np.array(df_X['fs'].iloc[0], dtype=np.float32)
+    N = y_meas.size
+    ts = 1/fs
+    t = np.arange(N)*ts
+
+    # In[Set-up model]
+
+    net = ParallelWHNet()
+    model_folder = os.path.join("models", model_name)
+    net.load_state_dict(torch.load(os.path.join(model_folder, f"{model_name}.pt")))
+
+    # In[Predict]
+    u_torch = torch.tensor(u[None, :, None],  dtype=torch.float, requires_grad=False)
+
+    with torch.no_grad():
+        y_hat = net(u_torch)
+
+    # In[Detach]
+
+    y_hat = y_hat.detach().numpy()[0, :, 0]
+
+    # In[Plot]
+    fig, ax = plt.subplots(2, 1, sharex=True)
+    ax[0].plot(t, y_meas, 'k', label="$y$")
+    ax[0].plot(t, y_hat, 'r', label="$\hat y$")
+    ax[0].plot(t, y_meas - y_hat, 'g', label="$e$")
+    ax[0].legend()
+    ax[0].grid()
+
+    ax[1].plot(t, u, 'k', label="$u$")
+    ax[1].legend()
+    ax[1].grid()
+
+    # In[Inspect linear model]
+
+    # First linear block
+    # a_coeff_1 = net.G1.a_coeff.detach().numpy()
+    # b_coeff_1 = net.G1.b_coeff.detach().numpy()
+    # a_poly_1 = np.empty_like(a_coeff_1, shape=(2, 2, net.na_1 + 1))
+    # a_poly_1[:, :, 0] = 1
+    # a_poly_1[:, :, 1:] = a_coeff_1[:, :, :]
+    # b_poly_1 = np.array(b_coeff_1)
+    # G1_sys = control.TransferFunction(b_poly_1, a_poly_1, ts)
+    #
+    # plt.figure()
+    # mag_G1_1, phase_G1_1, omega_G1_1 = control.bode(G1_sys[0, 0])
+    # plt.figure()
+    # mag_G1_2, phase_G1_2, omega_G1_2 = control.bode(G1_sys[1, 0])
+    #
+    # # Second linear block
+    # a_coeff_2 = net.G2.a_coeff.detach().numpy()
+    # b_coeff_2 = net.G2.b_coeff.detach().numpy()
+    # a_poly_2 = np.empty_like(a_coeff_2, shape=(2, 1, net.na_2 + 1))
+    # a_poly_2[:, :, 0] = 1
+    # a_poly_2[:, :, 1:] = a_coeff_2[:, :, :]
+    # b_poly_2 = np.array(b_coeff_2)
+    # G2_sys = control.TransferFunction(b_poly_2, a_poly_2, ts)
+
+    # plt.figure()
+    # mag_G2_1, phase_G2_1, omega_G2_1 = control.bode(G2_sys[0, 0])
+    # plt.figure()
+    # mag_G2_2, phase_G2_2, omega_G2_2 = control.bode(G2_sys[0, 1])
+
+
+    # In[Metrics]
+
+    idx_test = range(n_skip, N)
+
+    e_rms = 1000*util.metrics.error_rmse(y_meas[idx_test], y_hat[idx_test])
+    mae = 1000 * util.metrics.error_mae(y_meas[idx_test], y_hat[idx_test])
+    fit_idx = util.metrics.fit_index(y_meas[idx_test], y_hat[idx_test])
+    r_sq = util.metrics.r_squared(y_meas[idx_test], y_hat[idx_test])
+    u_rms = 1000*util.metrics.error_rmse(u, 0)
+
+    print(f"RMSE: {e_rms:.2f}mV\nMAE: {mae:.2f}mV\nFIT:  {fit_idx:.1f}%\nR_sq: {r_sq:.1f}\nRMSU: {u_rms:.2f}mV")