I am implementing a Model Predictive Control (MPC) to maintain the temperature of a room within a defined range while minimizing the control input. I started by training an LSTM neural network with Keras that takes a window of past data to predict the next output, but I couldn’t implement it with the existing optimization libraries. Recently, I read that MPC can be implemented with L4CasADi, but the model needs to be in PyTorch. I tried to convert my Keras model to a PyTorch model, but I couldn’t succeed, so I decided to retrain another model with PyTorch. For this, I read that the model needs to be traceable and differentiable L4CasADi. Can someone please explain this to me? Here is my model that I trained with Keras:

#%% Import library 
import os 
os.chdir(r'C:algoritmeMPC')

import warnings
import matplotlib.pyplot as plt

# with warnings.catch_warnings():
warnings.simplefilter("ignore")
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import time

import matplotlib.pyplot as plt
from pickle import dump, load
from sklearn.preprocessing import MinMaxScaler

# For LSTM model
import tensorflow as tf
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from keras.layers import Reshape
from keras.callbacks import EarlyStopping
from keras.models import load_model

# construire dataset
data = pd.DataFrame({'T_external': T_ext, 'beta': control,'T_int': T_int})
#%% Scale a feature 
window = 10
P = 10

# Scale features
s1 = MinMaxScaler(feature_range=(0,1))
Xs = s1.fit_transform(data[['T_external', 'beta', 'T_int']])
s2 = MinMaxScaler(feature_range=(0,1))
Ys = s2.fit_transform(data[['T_int']])

# # Save model parameters
model_params = dict()
model_params['Xscale'] = s1
model_params['yscale'] = s2
model_params['window'] = window
model_params['P'] = P # prediction horizon

path = os.getcwd()
print(path)
dump(model_params, open(path + 'model_param_MIMO_W5_palie=10.pkl', 'wb'))

#%% Prepa data
nstep = Xs.shape[0]
val_ratio = 0.8
cut_index = int(nstep*val_ratio) 
print(cut_index)
Xs_train = Xs[0:cut_index]
Ys_train = Ys[0:cut_index]
Xs_val = Xs[cut_index:]
Ys_val = Ys[cut_index:]
#%% preparer les donnees du test et entrainement

X_train, Y_train, X_val, Y_val = [], [], [], []
# Construct Y_train and Y_val correctly

for i in range(len(Xs_train)):
        end_ix = i + window
        out_end_ix = end_ix + P
        if out_end_ix > len(Xs_train):
            break
        X_train.append(Xs_train[i:end_ix,:])
        Y_train.append(Ys_train[end_ix:out_end_ix])

for i in range(len(Xs_val)):
        end_ix = i + window
        out_end_ix = end_ix + P
        if out_end_ix > len(Xs_val):
            break
        X_val.append(Xs_val[i:end_ix,:])
        Y_val.append(Ys_val[end_ix:out_end_ix])
        
        
    # Reshape data to format accepted by LSTM
     X_train, Y_train = np.array(X_train), np.array(Y_train)
     X_val, Y_val = np.array(X_val), np.array(Y_val)
    
    return (X_train,Y_train,X_val,Y_val)
        
    
(X_train,Y_train,X_val,Y_val) = diviser_donnees(Xs_train,Ys_train,Xs_val,Ys_val)


#%% Initiliser LSTM network
model_lstm  = Sequential()
model_lstm .add(LSTM(128, input_shape=(X_train.shape[1],X_train.shape[2]), return_sequences=True))
model_lstm .add(tf.keras.layers.LeakyReLU(alpha=0.01))
model_lstm .add(LSTM(128,return_sequences=True))
model_lstm .add(tf.keras.layers.LeakyReLU(alpha=0.01))
model_lstm .add(Dropout(0.2))
model_lstm .add(LSTM(64,return_sequences=False))
model_lstm .add(Dropout(0.2))
model_lstm.add(tf.keras.layers.Flatten())
model_lstm .add(Dense(units = Y_train.shape[1]*Y_train.shape[2])) 
model_lstm.add(Reshape((Y_train.shape[1], Y_train.shape[2])))

model_lstm .summary()

#%% compile the modele
from tensorflow.keras.optimizers import Adam
early_stopping = EarlyStopping(monitor='loss', patience=10, mode='min')
model_lstm .compile(optimizer=Adam(learning_rate=0.0001), loss='mse', metrics=['mae'])

#%%  Fit (and time) LSTM model pour window = 10
t0 = time.time()
history  = model_lstm .fit(X_train, Y_train, epochs = 300, batch_size = 32,
                             callbacks=[early_stopping],
                             validation_data=(X_val, Y_val))
t1 = time.time()
print('Runtime: %.2f s' %(t1-t0))