diff --git a/server/callbacks.py b/server/callbacks.py
index 508a7de..2f34a16 100644
--- a/server/callbacks.py
+++ b/server/callbacks.py
@@ -2,18 +2,30 @@ import casadi as cs
 import numpy as np
 import tensorflow as tf
 
+from helpers import get_combined_evaluator
+
 # Package the resulting regression model in a CasADi callback
 class GPR(cs.Callback):
-    def __init__(self, name, model, opts={}):
+    def __init__(self, name, model, n_in, opts={}):
         cs.Callback.__init__(self)
 
         self.model = model
-        self.n_in = model.data[0].shape[1]
+        self.combined_evaluator = get_combined_evaluator(model)
+        self.n_in = n_in
+        self.tf_var = tf.Variable(np.ones((1, self.n_in)), dtype = tf.float64)
+        self.grads = None
         # Create a variable to keep all the gradient callback references
         self.refs = []
 
         self.construct(name, opts)
     
+    # Update tf_evaluator
+    def update_model(self, model):
+        self.model = model
+        self.combined_evaluator = get_combined_evaluator(model)
+
+    def uses_output(self): return True
+    
     # Number of inputs/outputs
     def get_n_in(self): return 1
     def get_n_out(self): return 1
@@ -27,14 +39,15 @@ class GPR(cs.Callback):
 
 
     def eval(self, arg):
-        inp = np.array(arg[0])
-        inp = tf.Variable(inp, dtype=tf.float64)
-        [mean, _] = self.model.predict_f(inp)
+        self.tf_var.assign(arg[0])
+        preds, grads = self.combined_evaluator(self.tf_var)
+        [mean, _] = preds
+        self.grads = grads
         return [mean.numpy()]
     
     def has_reverse(self, nadj): return nadj==1
     def get_reverse(self, nadj, name, inames, onames, opts):
-        grad_callback = GPR_grad(name, self.model)
+        grad_callback = GPR_grad(name, self.n_in, self.combined_evaluator)
         self.refs.append(grad_callback)
         
         nominal_in = self.mx_in()
@@ -43,10 +56,12 @@ class GPR(cs.Callback):
         return cs.Function(name, nominal_in+nominal_out+adj_seed, grad_callback.call(nominal_in), inames, onames)
         
 class GPR_grad(cs.Callback):
-    def __init__(self, name, model, opts={}):
-        cs.Callback.__init__(self)  
-        self.model = model
-        self.n_in = model.data[0].shape[1]
+    def __init__(self, name, n_in, combined_evaluator, opts={}):
+        cs.Callback.__init__(self)
+        
+        self.combined_evaluator = combined_evaluator
+        self.n_in = n_in
+        self.tf_var = tf.Variable(np.ones((1, self.n_in)), dtype = tf.float64)
 
         self.construct(name, opts)
 
@@ -58,14 +73,10 @@ class GPR_grad(cs.Callback):
         return cs.Sparsity.dense(1,self.n_in)
     def get_sparsity_out(self,i):
         return cs.Sparsity.dense(1,self.n_in)
-
+    
 
     def eval(self, arg):
-        inp = np.array(arg[0])
-        inp = tf.Variable(inp, dtype=tf.float64)
-        
-        with tf.GradientTape() as tape:
-            preds = self.model.predict_f(inp)
+        self.tf_var.assign(arg[0])
+        _, grads = self.combined_evaluator(self.tf_var)
 
-        grads = tape.gradient(preds, inp)
         return [grads.numpy()]
diff --git a/server/controllers.py b/server/controllers.py
index 4fa91cf..2c57038 100644
--- a/server/controllers.py
+++ b/server/controllers.py
@@ -8,6 +8,10 @@ import pandas as pd
 import gpflow
 import tensorflow as tf
 
+from gpflow.ci_utils import ci_niter
+
+import matplotlib.pyplot as plt
+
 from sklearn.preprocessing import MinMaxScaler
 
 import callbacks
@@ -43,44 +47,14 @@ class PIDcontroller:
         return sig_P + sig_I + sig_D
 
 
-
-class GP_MPCcontroller:
+class Base_MPCcontroller(object):
     def __init__(self, dict_cols, model = None, scaler = None, N_horizon = 10, recover_from_crash = False):
 
-
-        self.recover_from_crash = recover_from_crash
-
-        if self.recover_from_crash:
-            self.model = pickle.load(open("controller_model.pkl", 'rb'))
-            self.scaler = pickle.load(open("controller_scaler.pkl", 'rb'))
-            self.X_log = pickle.load(open("controller_X_log.pkl", 'rb'))
-            df = pd.read_pickle("controller_df.pkl")
-            self.recovery_signal = iter(df['SimulatedHeat'])
-
-
-        if model is not None:
-            # Model is already trained. Using as is.
-            if scaler is None: raise ValueError("Not allowed to pass a model without a scaler")
-            self.model = model
-            self.cs_model = callbacks.GPR("gpr", self.model)
-            self.scaler = scaler
-            self.scaler_helper = ScalerHelper(self.scaler)
-        else:
-            # No model has been passed. Setting up model initialization
-            self.model = None
-            self.nb_data = 500 + 1
-            self.ident_signal = get_random_signal(self.nb_data, signal_type = 'analog')
-
-
-            self.Pel = 2 * 6300
-            self.COP = 5.0
-
-            self.ident_signal = iter(self.Pel * self.COP * self.ident_signal)
-
-
         self.dict_cols = dict_cols
         self.max_lag = max([lag for lag,_ in self.dict_cols.values()])
         self.N_horizon = N_horizon
+        self.n_states = np.sum([len(cols) * lags for lags,cols in
+            self.dict_cols.values()])
         self.X_log = []
 
 
@@ -91,119 +65,50 @@ class GP_MPCcontroller:
             self.data_cols += cols
         self.data = np.empty((0, len(self.data_cols)))
 
+        # Dataset used for training
+        self.dataset_train_minsize = 500
+        self.dataset_train_maxsize = 500
+        self.dataset_train = np.empty((0, self.n_states))
+
         # The current weather forcast
         self.weather_forecast = None
 
         # Current measurements
         self.w, self.u, self.y = None, None, None
 
+        # Recover from a previous crash with precomputed values and continue
+        self.recover_from_crash = recover_from_crash
 
+        if self.recover_from_crash:
+            self.model = pickle.load(open("controller_model.pkl", 'rb'))
+            self.scaler = pickle.load(open("controller_scaler.pkl", 'rb'))
+            self.scaler_helper = ScalerHelper(self.scaler)
+            self.X_log = pickle.load(open("controller_X_log.pkl", 'rb'))
+            df = pd.read_pickle("controller_df.pkl")
+            self.recovery_signal = iter(df['SimulatedHeat'])
+            return
 
-    ###
-    # GPflow model training and update
-    ###
-    def _train_model(self):
-        """ Identify model from gathered data """
+        # Pre-existing model passed. Load all the necessary objects
+        if model is not None:
+            # Model is already trained. Using as is.
+            if scaler is None: raise ValueError("Not allowed to pass a model without a scaler")
+            self.model = model
+            self.cs_model = callbacks.GPR("gpr", self.model, self.n_states)
+            self.scaler = scaler
+            self.scaler_helper = ScalerHelper(self.scaler)
+        # No pre-existing model. Set up data acquisition and model training
+        else:
+            # No model has been passed. Setting up model initialization
+            self.model = None
 
-        nb_train_pts = self.nb_data - 1
-        ###
-        # Dataset
-        ###
-        df = pd.DataFrame(self.data[:nb_train_pts], columns = self.data_cols)
-        self.scaler = MinMaxScaler(feature_range = (-1, 1)) 
-        self.scaler_helper = ScalerHelper(self.scaler)
-        df_sc = get_scaled_df(df, self.dict_cols, self.scaler)
-        df_gpr_train = data_to_gpr(df_sc, self.dict_cols)
+            # Define an identification signal to be used first
+            self.Pel = 2 * 6300
+            self.COP = 5.0
 
-        df_input_train = df_gpr_train.drop(columns = self.dict_cols['w'][1] + self.dict_cols['u'][1] + self.dict_cols['y'][1])
-        df_output_train = df_gpr_train[self.dict_cols['y'][1]]
+            self.ident_signal = get_identification_signal(size = self.dataset_train_minsize)
+            self.ident_signal = iter(self.COP * self.Pel * self.ident_signal)
 
-        np_input_train = df_input_train.to_numpy()
-        np_output_train = df_output_train.to_numpy().reshape(-1, 1)
-
-        data_train = (np_input_train, np_output_train)
-
-        df_test = pd.DataFrame(self.data[nb_train_pts:], columns = self.data_cols)
-        df_test_sc = get_scaled_df(df_test, self.dict_cols, self.scaler)
-        df_gpr_test = data_to_gpr(df_test_sc, self.dict_cols)
-        df_input_test = df_gpr_test.drop(columns = self.dict_cols['w'][1] + self.dict_cols['u'][1] + self.dict_cols['y'][1])
-        df_output_test = df_gpr_test[self.dict_cols['y'][1]]
-        np_input_test = df_input_test.to_numpy()
-        np_output_test = df_output_test.to_numpy()
-
-        ###
-        # Kernel
-        ###
-
-        nb_dims = np_input_train.shape[1]
-        rational_dims = np.arange(0, (self.dict_cols['t'][0] + 1) * len(self.dict_cols['t'][1]), 1)
-        nb_rational_dims = len(rational_dims)
-        squared_dims = np.arange(nb_rational_dims, nb_dims, 1)
-        nb_squared_dims = len(squared_dims)
-
-        default_lscale = 75
-        while True:
-            try:
-
-                squared_l = np.linspace(default_lscale, 2*default_lscale, nb_squared_dims)
-                rational_l = np.linspace(default_lscale, 2*default_lscale, nb_rational_dims)
-
-                variance = tf.math.reduce_variance(np_input_train)
-
-                k0 = gpflow.kernels.SquaredExponential(lengthscales = squared_l, active_dims = squared_dims, variance = variance)
-                k1 = gpflow.kernels.Constant(variance = variance)
-                k2 = gpflow.kernels.RationalQuadratic(lengthscales = rational_l, active_dims = rational_dims, variance = variance)
-                k3 = gpflow.kernels.Periodic(k2)
-
-                k = (k0 + k1) * k2
-                k = k0
-
-                ###
-                # Model
-                ###
-
-                m = gpflow.models.GPR(
-                    data = data_train,
-                    kernel = k,
-                    mean_function = None,
-                )
-
-                ###
-                # Training
-                ###
-                print(f"Training a model with lscale:{default_lscale}")
-                opt = gpflow.optimizers.Scipy()
-                opt.minimize(m.training_loss, m.trainable_variables)
-                break
-            except:
-                print(f"Failed.Increasing lengthscale")
-                default_lscale += 5
-
-        ###
-        # Save model
-        ###
-        self.model = m
-        self.n_states = self.model.data[0].shape[1]
-
-#        # Manual model validation
-#        import matplotlib.pyplot as plt
-#
-#        plt.figure()
-#        plt.plot(data_train[1], label = 'real')
-#        mean, var = self.model.predict_f(data_train[0])
-#        plt.plot(mean, label = 'model')
-#        plt.legend()
-#        plt.show()
-#
-#        plt.figure()
-#        plt.plot(np_output_test, label = 'real')
-#        mean, var = self.model.predict_f(np_input_test)
-#        plt.plot(mean, label = 'model')
-#        plt.legend()
-#        plt.show()
-#
-#        import pdb; pdb.set_trace()
-        pass
+        return
 
     ###
     # Update measurements
@@ -218,18 +123,16 @@ class GP_MPCcontroller:
     def _add_input_measurement(self, u):
         self.u = np.array(u).reshape(1, -1)
 
+    def _add_measurement_set(self):
+        new_data = np.hstack([self.w, self.u, self.y])
+        self.data = np.vstack([self.data, new_data])
+        print(f"{self.data.shape[0]} data points. Newest: {new_data}")
+        self.w, self.u, self.y = None, None, None
+
     def set_weather_forecast(self, W):
         assert (W.shape[0] == self.N_horizon)
         self.weather_forecast = W
 
-    def update_model(self):
-        new_data = np.hstack([self.w, self.u, self.y])
-        print(f"Adding new data: {new_data}")
-        self.data = np.vstack([self.data, new_data])
-        print(f"Data size: {self.data.shape[0]}")
-        self.w, self.u, self.y = None, None, None
-        print("---")
-
     ###
     # Set up optimal problem solver
     ###
@@ -239,10 +142,11 @@ class GP_MPCcontroller:
         ###
         # Initialization
         ###
-        self.cs_model = callbacks.GPR("gpr", self.model)
+        self.cs_model = callbacks.GPR("gpr", self.model, self.n_states)
 
         T_set = 21
         T_set_sc = self.scaler_helper.scale_output(T_set)
+        self.T_set_sc = T_set_sc
 
 
         X = cs.MX.sym("X", self.N_horizon + 1, self.n_states)
@@ -311,7 +215,7 @@ class GP_MPCcontroller:
 
         prob = {'f': J, 'x': cs.vec(X), 'g': cs.vertcat(*g), 'p': p}
         options = {"ipopt": {"hessian_approximation": "limited-memory", "max_iter": 100,
-                             #"acceptable_tol": 1e-6, "tol": 1e-6,
+                             "acceptable_tol": 1e-4, "tol": 1e-4,
                              #"linear_solver": "ma57",
                              #"acceptable_obj_change_tol": 1e-5,
                              #"mu_strategy": "adaptive",
@@ -341,6 +245,7 @@ class GP_MPCcontroller:
             except StopIteration:
                 # Finished passing in the pre-existing control
                 # Switching to normal operation
+                self._setup_solver()
                 self.recover_from_crash = False
 
         # Training pre-loop
@@ -360,6 +265,10 @@ class GP_MPCcontroller:
 
 
         data_scaled = self.scaler.transform(self.data)
+        # Append a dummy row to data, to align data_to_gpr cols (which include measureemnt
+        # at time t, with the current measurements, that stop at t-1)
+        dummy_row = np.nan * np.ones((1, self.data.shape[1]))
+        data_scaled = np.vstack([data_scaled, dummy_row])
         df = pd.DataFrame(data_scaled, columns = self.data_cols)
         
         x0 = data_to_gpr(df, self.dict_cols).drop(
@@ -384,6 +293,8 @@ class GP_MPCcontroller:
 
         X = np.array(res['x'].reshape((self.N_horizon + 1, self.n_states)))
         self.X_log.append(X)
+        df_X = pd.DataFrame(X)
+        #import pdb; pdb.set_trace()
         u_idx = self.dict_cols['w'][0] * len(self.dict_cols['w'][1])
         # Take the first control action and apply it
         u = X[1, u_idx]
@@ -403,6 +314,363 @@ class GP_MPCcontroller:
         return
 
 
+
+class GP_MPCcontroller(Base_MPCcontroller):
+    def __init__(self, dict_cols, model = None, scaler = None, N_horizon = 10, recover_from_crash = False):
+        super().__init__(dict_cols, model, scaler, N_horizon, recover_from_crash)
+
+    ###
+    # GPflow model training and update
+    ###
+    def _train_model(self):
+        """ Identify model from gathered data """
+
+        nb_train_pts = self.dataset_train_minsize - 1
+        nb_train_pts = 500
+        ###
+        # Dataset
+        ###
+        df = pd.DataFrame(self.data[:nb_train_pts], columns = self.data_cols)
+        self.scaler = MinMaxScaler(feature_range = (-1, 1)) 
+        self.scaler_helper = ScalerHelper(self.scaler)
+        df_sc = get_scaled_df(df, self.dict_cols, self.scaler)
+        df_gpr_train = data_to_gpr(df_sc, self.dict_cols)
+
+        df_input_train = df_gpr_train.drop(columns = self.dict_cols['w'][1] + self.dict_cols['u'][1] + self.dict_cols['y'][1])
+        df_output_train = df_gpr_train[self.dict_cols['y'][1]]
+
+        np_input_train = df_input_train.to_numpy()
+        np_output_train = df_output_train.to_numpy().reshape(-1, 1)
+
+        data_train = (np_input_train, np_output_train)
+
+        df_test = pd.DataFrame(self.data[nb_train_pts:], columns = self.data_cols)
+        df_test_sc = get_scaled_df(df_test, self.dict_cols, self.scaler)
+        df_gpr_test = data_to_gpr(df_test_sc, self.dict_cols)
+        df_input_test = df_gpr_test.drop(columns = self.dict_cols['w'][1] + self.dict_cols['u'][1] + self.dict_cols['y'][1])
+        df_output_test = df_gpr_test[self.dict_cols['y'][1]]
+        np_input_test = df_input_test.to_numpy()
+        np_output_test = df_output_test.to_numpy()
+
+        ###
+        # Kernel
+        ###
+
+        nb_dims = np_input_train.shape[1]
+        rational_dims = np.arange(0, (self.dict_cols['t'][0] + 1) * len(self.dict_cols['t'][1]), 1)
+        nb_rational_dims = len(rational_dims)
+        squared_dims = np.arange(nb_rational_dims, nb_dims, 1)
+        nb_squared_dims = len(squared_dims)
+
+        default_lscale = 150
+        while True:
+            try:
+
+                squared_l = np.linspace(default_lscale, 2*default_lscale, nb_squared_dims)
+                rational_l = np.linspace(default_lscale, 2*default_lscale, nb_rational_dims)
+
+                variance = tf.math.reduce_variance(np_input_train)
+
+                k0 = gpflow.kernels.SquaredExponential(lengthscales = squared_l, active_dims = squared_dims, variance = variance)
+                k1 = gpflow.kernels.Constant(variance = variance)
+                k2 = gpflow.kernels.RationalQuadratic(lengthscales = rational_l, active_dims = rational_dims, variance = variance)
+                k3 = gpflow.kernels.Periodic(k2)
+
+                k = (k0 + k1) * k2
+                k = k0
+
+                ###
+                # Model
+                ###
+
+                m = gpflow.models.GPR(
+                    data = data_train,
+                    kernel = k,
+                    mean_function = None,
+                )
+
+                ###
+                # Training
+                ###
+                print(f"Training a model with lscale:{default_lscale}")
+                opt = gpflow.optimizers.Scipy()
+                opt.minimize(m.training_loss, m.trainable_variables)
+                break
+            except:
+                print(f"Failed.Increasing lengthscale")
+                default_lscale += 10
+
+        ###
+        # Save model
+        ###
+        self.model = m
+
+#        plt.figure()
+#
+#        # Testing on training data
+#        mean, var = m.predict_f(np_input_train)
+#
+#        plt.plot(df_input_train.index, np_output_train[:, :], label = 'Measured data')
+#        plt.plot(df_input_train.index, mean[:, :], label = 'Gaussian Process Prediction')
+#        plt.fill_between(
+#            df_input_train.index,
+#            mean[:, 0] - 1.96 * np.sqrt(var[:, 0]),
+#            mean[:, 0] + 1.96 * np.sqrt(var[:, 0]),
+#            alpha = 0.2
+#        )
+#        plt.show()
+#
+#        plt.figure()
+#        # Testing on testing data
+#        mean, var = m.predict_f(np_input_test)
+#
+#        plt.plot(df_input_test.index, np_output_test[:, :], label = 'Measured data')
+#        plt.plot(df_input_test.index, mean[:, :], label = 'Gaussian Process Prediction')
+#        plt.fill_between(
+#            df_input_test.index,
+#            mean[:, 0] - 1.96 * np.sqrt(var[:, 0]),
+#            mean[:, 0] + 1.96 * np.sqrt(var[:, 0]),
+#            alpha = 0.2
+#        )
+#        plt.show()
+#
+#
+#        import pdb; pdb.set_trace()
+#        plt.figure()
+#
+#        start_idx = 25
+#        nb_predictions = 25
+#        N_pred = 20
+#
+#        plt.figure()
+#
+#        y_name = self.dict_cols['y'][1][0]
+#        for idx in range(start_idx, start_idx + nb_predictions):
+#            df_iter = df_input_test.iloc[idx:(idx + N_pred)].copy()
+#            for idxx in range(N_pred - 1):
+#                idx_old = df_iter.index[idxx]
+#                idx_new = df_iter.index[idxx+1]
+#                mean, var = m.predict_f(df_iter.loc[idx_old, :].to_numpy().reshape(1, -1))
+#                df_iter.loc[idx_new, f'{y_name}_1'] = mean.numpy().flatten()
+#                for lag in range(2, self.dict_cols['y'][0] + 1):
+#                    df_iter.loc[idx_new, f"{y_name}_{lag}"] = df_iter.loc[idx_old, f"{y_name}_{lag-1}"]
+#
+#            mean_iter, var_iter = m.predict_f(df_iter.to_numpy())
+#            plt.plot(df_iter.index, mean_iter.numpy(), '.-', label = 'predicted', color = 'orange')
+#        plt.plot(df_output_test.iloc[start_idx:start_idx + nb_predictions + N_pred], 'o-', label = 'measured', color = 'darkblue')
+#        plt.title(f"Prediction over {N_pred} steps")
+#
+#        plt.show()
+
+        return
+
+
+    def update_model(self):
+        self._add_measurement_set()
+
+
+class SVGP_MPCcontroller(Base_MPCcontroller):
+    def __init__(self, dict_cols, model = None, scaler = None, N_horizon = 10, recover_from_crash = False):
+        super().__init__(dict_cols, model, scaler, N_horizon, recover_from_crash)
+
+        # Adaptive models update parameters
+        self.model_update_frequency = 100
+        self.pts_since_update = 0
+        
+        # Model log
+        self.model_log = []
+
+    ###
+    # GPflow model training and update
+    ###
+    def _train_model(self):
+        """ Identify model from gathered data """
+
+        #nb_train_pts = self.dataset_train_minsize - 1
+        nb_train_pts = self.data.shape[0]
+        ###
+        # Dataset
+        ###
+        df = pd.DataFrame(self.data[:nb_train_pts], columns = self.data_cols)
+        self.scaler = MinMaxScaler(feature_range = (-1, 1)) 
+        self.scaler_helper = ScalerHelper(self.scaler)
+        df_sc = get_scaled_df(df, self.dict_cols, self.scaler)
+        df_gpr_train = data_to_gpr(df_sc, self.dict_cols)
+
+        df_input_train = df_gpr_train.drop(columns = self.dict_cols['w'][1] + self.dict_cols['u'][1] + self.dict_cols['y'][1])
+        df_output_train = df_gpr_train[self.dict_cols['y'][1]]
+
+        np_input_train = df_input_train.to_numpy()
+        np_output_train = df_output_train.to_numpy().reshape(-1, 1)
+
+        data_train = (np_input_train, np_output_train)
+
+        ###
+        # Kernel
+        ###
+
+        nb_dims = np_input_train.shape[1]
+        rational_dims = np.arange(0, (self.dict_cols['t'][0] + 1) * len(self.dict_cols['t'][1]), 1)
+        nb_rational_dims = len(rational_dims)
+        squared_dims = np.arange(nb_rational_dims, nb_dims, 1)
+        nb_squared_dims = len(squared_dims)
+
+        default_lscale = 1
+
+        squared_l = np.linspace(default_lscale, 2*default_lscale, nb_squared_dims)
+        rational_l = np.linspace(default_lscale, 2*default_lscale, nb_rational_dims)
+
+        variance = tf.math.reduce_variance(np_input_train)
+
+        k0 = gpflow.kernels.SquaredExponential(lengthscales = squared_l, active_dims = squared_dims, variance = variance)
+        k1 = gpflow.kernels.Constant(variance = variance)
+        k2 = gpflow.kernels.RationalQuadratic(lengthscales = rational_l, active_dims = rational_dims, variance = variance)
+        k3 = gpflow.kernels.Periodic(k2)
+
+        k = (k0 + k1) * k2
+        k = k0
+
+        ###
+        # Model
+        ###
+
+        N = data_train[0].shape[0]
+        M = 150 # Number of inducing locations
+        Z = data_train[0][:M, :].copy()
+
+        m = gpflow.models.SVGP(k, gpflow.likelihoods.Gaussian(), Z, num_data = N)
+
+        elbo = tf.function(m.elbo)
+
+        ###
+        # Training
+        ###
+
+        minibatch_size = 100
+        train_dataset = tf.data.Dataset.from_tensor_slices(data_train).repeat().shuffle(N)
+
+        # Turn off training for inducing point locations
+        gpflow.set_trainable(m.inducing_variable, False)
+
+        def run_adam(model, iterations):
+            """
+            Utility function running the Adam optimizer
+
+            :param model: GPflow model
+            :param interations: number of iterations
+            """
+            # Create an Adam Optimizer action
+            logf = []
+            train_iter = iter(train_dataset.batch(minibatch_size))
+            training_loss = model.training_loss_closure(train_iter, compile=True)
+            optimizer = tf.optimizers.Adam()
+
+            @tf.function
+            def optimization_step():
+                optimizer.minimize(training_loss, model.trainable_variables)
+
+            for step in range(iterations):
+                optimization_step()
+                if step % 10 == 0:
+                    elbo = -training_loss().numpy()
+                    logf.append(elbo)
+            return logf
+
+        
+        maxiter = ci_niter(10000)
+        logf = run_adam(m, maxiter)
+
+        ###
+        # Save model
+        ###
+        self.model = m
+
+#        plt.figure()
+#
+#        # Testing on training data
+#        mean, var = m.predict_f(np_input_train)
+#
+#        plt.plot(df_input_train.index, np_output_train[:, :], label = 'Measured data')
+#        plt.plot(df_input_train.index, mean[:, :], label = 'Gaussian Process Prediction')
+#        plt.fill_between(
+#            df_input_train.index,
+#            mean[:, 0] - 1.96 * np.sqrt(var[:, 0]),
+#            mean[:, 0] + 1.96 * np.sqrt(var[:, 0]),
+#            alpha = 0.2
+#        )
+#        plt.show()
+#
+#        plt.figure()
+#        # Testing on testing data
+#        mean, var = m.predict_f(np_input_test)
+#
+#        plt.plot(df_input_test.index, np_output_test[:, :], label = 'Measured data')
+#        plt.plot(df_input_test.index, mean[:, :], label = 'Gaussian Process Prediction')
+#        plt.fill_between(
+#            df_input_test.index,
+#            mean[:, 0] - 1.96 * np.sqrt(var[:, 0]),
+#            mean[:, 0] + 1.96 * np.sqrt(var[:, 0]),
+#            alpha = 0.2
+#        )
+#        plt.show()
+#
+#
+#        plt.figure()
+#
+#        start_idx = 25
+#        nb_predictions = 25
+#        N_pred = 20
+#
+#        plt.figure()
+#
+#        y_name = self.dict_cols['y'][1][0]
+#        for idx in range(start_idx, start_idx + nb_predictions):
+#            df_iter = df_input_test.iloc[idx:(idx + N_pred)].copy()
+#            for idxx in range(N_pred - 1):
+#                idx_old = df_iter.index[idxx]
+#                idx_new = df_iter.index[idxx+1]
+#                mean, var = m.predict_f(df_iter.loc[idx_old, :].to_numpy().reshape(1, -1))
+#                df_iter.loc[idx_new, f'{y_name}_1'] = mean.numpy().flatten()
+#                for lag in range(2, self.dict_cols['y'][0] + 1):
+#                    df_iter.loc[idx_new, f"{y_name}_{lag}"] = df_iter.loc[idx_old, f"{y_name}_{lag-1}"]
+#
+#            mean_iter, var_iter = m.predict_f(df_iter.to_numpy())
+#            plt.plot(df_iter.index, mean_iter.numpy(), '.-', label = 'predicted', color = 'orange')
+#        plt.plot(df_output_test.iloc[start_idx:start_idx + nb_predictions + N_pred], 'o-', label = 'measured', color = 'darkblue')
+#        plt.title(f"Prediction over {N_pred} steps")
+#
+#        plt.show()
+        return
+
+
+    def update_model(self):
+        self._add_measurement_set()
+        if self.model is not None and not self.recover_from_crash:
+            print(f"Points since last update: {self.pts_since_update}")
+            if self.pts_since_update < self.model_update_frequency:
+                self.pts_since_update += 1
+            else:
+                print(f"Updating model after {self.pts_since_update} measurements")
+                # Append old model to log
+                self.model_log.append(self.model)
+                # Train new model
+                self._train_model()
+                # Re-initialize CasADi solver
+                self._setup_solver()
+                self.pts_since_update = 0
+                
+    # Redefine save_data since now we're also saving the model log
+    def save_data(self):
+        df = pd.DataFrame(self.data, columns = self.data_cols)
+        df.to_pickle("controller_df.pkl")
+
+        pickle.dump(self.scaler, open(Path("controller_scaler.pkl"), 'wb'))
+        pickle.dump(self.model, open(Path("controller_model.pkl"), 'wb'))
+        pickle.dump(self.X_log, open(Path("controller_X_log.pkl"), 'wb'))
+        pickle.dump(self.model_log, open(Path("controller_model_log.pkl"), 'wb'))
+        
+        return
+
 class TestController:
     def __init__(self):
         return
diff --git a/server/helpers.py b/server/helpers.py
index 03af4f1..1a04e94 100644
--- a/server/helpers.py
+++ b/server/helpers.py
@@ -1,9 +1,36 @@
 import pandas as pd
 import numpy as np
+
 import gpflow
+import tensorflow as tf
 
 from sklearn.exceptions import NotFittedError
 
+# Generator for tensorflow functions for a given model
+def get_model_evaluator(model):
+    @tf.function
+    def model_evaluator(tf_input):
+        preds = model.predict_f(tf_input)
+        return preds
+    return model_evaluator
+
+def get_grad_evaluator(model):
+    @tf.function
+    def grad_evaluator(tf_input):
+        with tf.GradientTape() as tape:
+            preds = model.predict_f(tf_input)
+        grads = tape.gradient(preds, tf_input)
+        return grads
+    return grad_evaluator
+
+def get_combined_evaluator(model):
+    @tf.function
+    def combined_evaluator(tf_input):
+        with tf.GradientTape() as tape:
+            preds = model.predict_f(tf_input)
+        grads = tape.gradient(preds, tf_input)
+        return preds, grads
+    return combined_evaluator
 
 def get_random_signal(nstep, a_range = (-1, 1), b_range = (2, 10), signal_type = 'analog'):
 
@@ -47,6 +74,20 @@ def get_random_signal(nstep, a_range = (-1, 1), b_range = (2, 10), signal_type =
         raise ValueError(signal_type)
 
 
+def get_identification_signal(size):
+    # Base random signal
+    rand_signal = get_random_signal(size, signal_type = 'prbs')
+    # Integrator (cumulative sum)
+    cum_signal = 3/size * np.ones((1, size))
+    cum_signal = np.cumsum(cum_signal)
+    # Combine signals and clip signal to [-1, 1] range
+    ident_signal = rand_signal + cum_signal
+    ident_signal = np.where(ident_signal < -1, -1, ident_signal)
+    ident_signal = np.where(ident_signal > 1, 1, ident_signal)
+
+    return ident_signal
+
+
 def get_scaled_df(df, dict_cols, scaler):
 
     """
diff --git a/server/server.py b/server/server.py
index e51045b..4c5dc5d 100644
--- a/server/server.py
+++ b/server/server.py
@@ -6,11 +6,12 @@ import pickle
 
 import numpy as np
 
-from controllers import GP_MPCcontroller, TestController
+from controllers import *
 
 clients = ['measure', 'control', 'weather']
 
-N_horizon = 6
+N_horizon = 8
+print(f"[*] Controller Horizon {N_horizon}")
 
 HOST            = '127.0.0.1'
 PORT = {
@@ -55,7 +56,7 @@ dict_cols = {
     'y': (y_lags, y_cols)
 }
 
-controller = GP_MPCcontroller(dict_cols = dict_cols, N_horizon = N_horizon, recover_from_crash = False)
+controller = SVGP_MPCcontroller(dict_cols = dict_cols, N_horizon = N_horizon, recover_from_crash = False)
 
 # Enter TCP server loop
 while True:
diff --git a/server/test.py b/server/test.py
index 72892f4..a1b5c70 100644
--- a/server/test.py
+++ b/server/test.py
@@ -17,7 +17,7 @@ import casadi as cs
 
 import callbacks
 from helpers import *
-from controllers import GP_MPCcontroller
+from controllers import *
 
 from time import sleep
 
@@ -40,8 +40,7 @@ dict_cols = {
 
 N_horizon = 10
 
-
-controller = GP_MPCcontroller(dict_cols = dict_cols, N_horizon = N_horizon)
+controller = SVGP_MPCcontroller(dict_cols = dict_cols, N_horizon = N_horizon)
 
 
 
@@ -64,6 +63,4 @@ while True:
     
     idx += 1
 
-    if idx > 502:
-        import pdb; pdb.set_trace()
 pass