Here's my code :

import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras
from keras.models import Sequential
from keras.layers import Dense
from firebase import firebase
import firebase_admin
from firebase_admin import credentials
from firebase_admin import db
import numpy as np
from sklearn.model_selection import train_test_split
from keras.layers import Dropout, Dense, LSTM, SimpleRNN
import pandas
from sklearn.preprocessing import MinMaxScaler


def generate_future_timesteps(data1, data2, data3, data4, num_timesteps):
    # Use the last `num_timesteps` data points as input

    max_length = max(len(data1), len(data2), len(data3), len(data4))
    future_input = np.zeros((4, max_length, 1))
    future_input[0, :len(data1), :] = np.array(data1).reshape(-1, 1)
    future_input[1, :len(data2), :] = np.array(data2).reshape(-1, 1)
    future_input[2, :len(data3), :] = np.array(data3).reshape(-1, 1)
    future_input[3, :len(data4), :] = np.array(data4).reshape(-1, 1)

    # Reshape the input data to match the model's input shape
    future_input = future_input[:, :num_timesteps, :]
    future_input = np.swapaxes(future_input, 0, 1)

    return future_input


def create_dataset(dataset, look_back):
    dataX, dataY = [], []
    for i in range(len(dataset) - look_back - 1):
        a = dataset[i:(i + look_back), 0]
        dataX.append(a)
        dataY.append(dataset[i + look_back, :])
    return np.array(dataX), np.array(dataY)


def create_rnn_model(trainX, trainY, testX, testY, rnn_units, look_back, input_shape = 0, name='RNN Model'):
    # Define the RNN model

    batch = 16
    if input_shape == 0:
        input_shape = (1, look_back)
    model = Sequential()
    model.add(SimpleRNN(rnn_units[0], activation='tanh', input_shape=input_shape, return_sequences=True))
    model.add(SimpleRNN(rnn_units[1], activation='tanh'))
    model.add(Dense(rnn_units[2]))
    model.add(Dense(4))

    # Compile the model
    model.compile(optimizer='adam', loss='mae')

    # Train the model on the training data
    history = model.fit(trainX, trainY, epochs=25, batch_size=batch, verbose=1, validation_split=0.2)

    # Evaluate the model on the testing data
    loss = model.evaluate(testX, testY, batch_size=batch)
    print(f'{name} loss: {loss}')

    return model


def create_lstm_model(trainX, trainY, testX, testY, lstm_units, look_back, input_shape = 0, name='LSTM Model'):
    # Define the LSTM autoencoder model

    batch = 16
    if input_shape == 0:
        input_shape = (1, look_back)
    model = Sequential()
    model.add(LSTM(lstm_units[0], activation='relu', input_shape = input_shape, return_sequences = True))
    model.add(LSTM(lstm_units[1], activation='relu', return_sequences = True))
    model.add(LSTM(lstm_units[2], activation='relu'))
    model.add(Dense(lstm_units[3]))
    model.add(Dense(4))

    # Compile the model
    model.compile(optimizer='adam', loss ='mae')

    # Train the model on the training data
    history = model.fit(trainX, trainY, epochs=25, batch_size=batch, verbose=1, validation_split=0.2)

    # Evaluate the model on the testing data
    loss = model.evaluate(testX, testY, batch_size=batch)
    print(f'{name} loss: {loss}')

    return model


def ensemble_models(data1, data2, data3, data4, look_back):
    min_length = min(len(data1), len(data2), len(data3), len(data4))

    # Slice each array to the minimum length
    data1 = data1[len(data1) - min_length:]
    data2 = data2[len(data2) - min_length:]
    data3 = data3[len(data3) - min_length:]
    data4 = data4[len(data4) - min_length:]

    # Pad the sequences with zeros so that they all have the same length
    max_length = max(len(data1), len(data2), len(data3), len(data4))
    all_data = np.zeros((4, max_length, 1))
    for i, data in enumerate([data1, data2, data3, data4]):
        all_data[i, :len(data), :] = np.array(data).reshape(-1, 1)

    d = {'hum': data1, 'light': data2, 'temp': data3, 'mois': data4}
    dataset = pandas.DataFrame(d)
    dataset = dataset.values
    dataset = dataset.astype("float32")

    scaler = MinMaxScaler()

    train_size = int(len(dataset) * 0.7)
    test_size = len(dataset) - train_size
    train, test = dataset[0:train_size, :], dataset[train_size:len(dataset), :]

    scaler.fit(train)
    train = scaler.transform(train)
    test = scaler.transform(test)

    trainX, trainY = create_dataset(train, look_back=look_back)
    testX, testY = create_dataset(test, look_back=look_back)

    trainX = np.reshape(trainX, (trainX.shape[0], 1, trainX.shape[1]))
    testX = np.reshape(testX, (testX.shape[0], 1, testX.shape[1]))

    # Create the first LSTM model
    lstm1 = create_lstm_model(trainX, trainY, testX, testY, lstm_units=[128, 64, 32, 64], look_back=look_back)

    # Create the second LSTM model with a different architecture
    #lstm2 = create_lstm_model(trainX, trainY, testX, testY, lstm_units=[64, 32, 16, 32], name='LSTM Model 2')
    rnn = create_rnn_model(trainX, trainY, testX, testY, rnn_units=[128, 64, 32, 64], look_back=look_back)

    # Get predictions from both models
    predictions1 = lstm1.predict(testX)
    predictions2 = rnn.predict(testX)

    predictions1 = scaler.inverse_transform(predictions1)
    predictions2 = scaler.inverse_transform(predictions2)

    # Combine the predictions using a weighted average
    ensemble_weight_1 = 0.5
    ensemble_weight_2 = 1.0 - ensemble_weight_1
    ensemble_predictions = ensemble_weight_1 * predictions1 + ensemble_weight_2 * predictions2

    # Evaluate the ensemble model’s performance
    testY = scaler.inverse_transform(testY)

    print_accuracy(ensemble_predictions, testY)

    # Now make predictions using the reinitialized model
    future_input = generate_future_timesteps(data1, data2, data3, data4, look_back)
    future_input = np.reshape(future_input, (-1, 1, look_back))

    # Scale the input data
    future_input = future_input.reshape(-1, 4)
    future_input = scaler.transform(future_input)
    future_input = future_input.reshape(-1, 1, look_back)

    # Get predictions for future timesteps
    future_predictions1 = lstm1.predict(future_input)
    future_predictions2 = rnn.predict(future_input)


    # Revert the predictions back to original scaling
    future_predictions1 = scaler.inverse_transform(future_predictions1)
    future_predictions2 = scaler.inverse_transform(future_predictions2)

    # Combine the predictions using a weighted average
    future_ensemble_predictions = (ensemble_weight_1 * future_predictions1 + ensemble_weight_2 * future_predictions2)

    # Reshape array2 to match the first dimension of array1
    array2_reshaped = np.repeat(future_ensemble_predictions, ensemble_predictions.shape[0] // future_ensemble_predictions.shape[0], axis=0)

    # Concatenate the reshaped arrays along axis 0
    future_ensemble_predictions = np.concatenate((ensemble_predictions, array2_reshaped), axis=0)

    plot_data(ensemble_predictions, future_ensemble_predictions, testY)

    return lstm1, rnn


def print_accuracy(ensemble_predictions, testY):
    mae = np.mean(np.abs(ensemble_predictions - testY), axis=0)
    print(f'Ensemble MAE: {mae}')

    max_humidity = 100  # Set the maximum possible value for humidity
    max_light = 100  # Set the maximum possible value for light
    max_temperature = 50  # Set the maximum possible value for temperature
    max_moisture = 100  # Set the maximum possible value for moisture

    mae_humidity = np.mean(np.abs(ensemble_predictions[:, 0] - testY[:, 0]))
    mae_light = np.mean(np.abs(ensemble_predictions[:, 1] - testY[:, 1]))
    mae_temperature = np.mean(np.abs(ensemble_predictions[:, 2] - testY[:, 2]))
    mae_moisture = np.mean(np.abs(ensemble_predictions[:, 3] - testY[:, 3]))

    percentage_error_humidity = (mae_humidity / max_humidity) * 100
    percentage_error_light = (mae_light / max_light) * 100
    percentage_error_temperature = (mae_temperature / max_temperature) * 100
    percentage_error_moisture = (mae_moisture / max_moisture) * 100

    print(f'Humidity Percentage Error: {round(percentage_error_humidity, 2)}%')
    print(f'Light Percentage Error: {round(percentage_error_light, 2)}%')
    print(f'Temperature Percentage Error: {round(percentage_error_temperature, 2)}%')
    print(f'Moisture Percentage Error: {round(percentage_error_moisture, 2)}%')

    print(f'Mean error: {round(np.mean([percentage_error_light, percentage_error_humidity, percentage_error_moisture, percentage_error_temperature]), 2)}%')

    return 0


def plot_data(ensemble_predictions, future_ensemble_predictions, testY):
    # (Optional) Visualization code

    plt.figure(figsize=(16, 8))

    predic_input = input("Comparer les prédictions avec les mesures ou prédire les mesures futures ? c/f")

    if predic_input == 'c':
        plt.plot(ensemble_predictions[:, 0], label='Humidity Predictions', linestyle='-')
        plt.plot(ensemble_predictions[:, 1], label='Light Predictions', linestyle='-')
        plt.plot(ensemble_predictions[:, 2], label='Temperature Predictions', linestyle='-')
        plt.plot(ensemble_predictions[:, 3], label='Moisture Predictions', linestyle='-')

    elif predic_input == 'f':
        plt.plot(future_ensemble_predictions[:, 0], label='Humidity Predictions', linestyle='-')
        plt.plot(future_ensemble_predictions[:, 1], label='Light Predictions', linestyle='-')
        plt.plot(future_ensemble_predictions[:, 2], label='Temperature Predictions', linestyle='-')
        plt.plot(future_ensemble_predictions[:, 3], label='Moisture Predictions', linestyle='-')

    else:

        print("Incorrect input received, plotting measured predictions...")

        plt.plot(ensemble_predictions[:, 0], label='Humidity Predictions', linestyle='-')
        plt.plot(ensemble_predictions[:, 1], label='Light Predictions', linestyle='-')
        plt.plot(ensemble_predictions[:, 2], label='Temperature Predictions', linestyle='-')
        plt.plot(ensemble_predictions[:, 3], label='Moisture Predictions', linestyle='-')

    plt.plot(testY[:, 0], label='Humidity Real', alpha=0.7, linestyle='dashed')
    plt.plot(testY[:, 1], label='Light Real', alpha=0.7, linestyle='dashed')
    plt.plot(testY[:, 2], label='Temperature Real', alpha=0.7, linestyle='dashed')
    plt.plot(testY[:, 3], label='Moisture Real', alpha=0.7, linestyle='dashed')

    plt.xlabel('Time Steps')
    plt.ylabel('Sensor Values')

    plt.legend()
    plt.show()

    data = input("print predicted data ? y/n")
    if data == "y":
        print(future_ensemble_predictions)

    return 0


def start_firebase():
    # Fetch the service account key JSON file contents
    cred = credentials.Certificate('gunadarma-monitoring-firebase-adminsdk-t4zec-95d7f006d3.json')
    # Initialize the app with a service account, granting admin privileges
    firebase_admin.initialize_app(cred, {
        'databaseURL': "https://gunadarma-monitoring-default-rtdb.asia-southeast1.firebasedatabase.app"
    })

    humidity = db.reference('/live/Humidity').get()
    light = db.reference('/live/Luminosity').get()
    temp = db.reference('/live/Temperature').get()
    moisture = db.reference('/live/moisture').get()
    print("humidity : ", humidity)
    print("luminosity : ", light)
    print("temperature : ", temp)
    print("moisture : ", moisture)

    record_hum = db.reference('/previous_Humidity').get()
    humrecord_list = [i for i in record_hum.values()]

    record_light = db.reference('/previous_Luminosity').get()
    lightrecord_list = [i for i in record_light.values()]

    record_temp = db.reference('/previous_Temperature').get()
    temprecord_list = [i for i in record_temp.values()]

    record_mois = db.reference('/previous_moisture').get()
    moisrecord_list = [i for i in record_mois.values()]

    return humrecord_list, lightrecord_list, temprecord_list, moisrecord_list


def remove_outliers(data1, data2, data3, data4, threshold=1):

    humrecord_list = np.array(data1)
    lightrecord_list = np.array(data2)
    temprecord_list = np.array(data3)
    moisrecord_list = np.array(data4)

    z_humidity = np.abs((humrecord_list - np.mean(humrecord_list)) / np.std(humrecord_list))
    z_light = np.abs((lightrecord_list - np.mean(lightrecord_list)) / np.std(lightrecord_list))
    z_temperature = np.abs((temprecord_list - np.mean(temprecord_list)) / np.std(temprecord_list))
    z_moisture = np.abs((moisrecord_list - np.mean(moisrecord_list)) / np.std(moisrecord_list))

    # Remove any data point that has a z-score greater than the threshold
    humrecord_list = humrecord_list[z_humidity < threshold]
    lightrecord_list = lightrecord_list[z_light < threshold]
    temprecord_list = temprecord_list[z_temperature < threshold]
    moisrecord_list = moisrecord_list[z_moisture < threshold]

    return humrecord_list, lightrecord_list, temprecord_list, moisrecord_list


if __name__ == '__main__':

    humrecord_list, lightrecord_list, temprecord_list, moisrecord_list = start_firebase()

    humrecord_list, lightrecord_list, temprecord_list, moisrecord_list = remove_outliers(humrecord_list, lightrecord_list, temprecord_list, moisrecord_list, threshold=3)


    if input("plot les données ? ") == "o":

        plt.figure(figsize=(16, 8))
        plt.plot(humrecord_list, color='r', label='hum')
        plt.plot(lightrecord_list, color='g', label='light')
        plt.plot(temprecord_list, color='b', label='temp')
        plt.plot(moisrecord_list, color='y', label='mois')
        plt.legend()
        plt.show()



    print("humidity history : (size ", len(humrecord_list), ") ", humrecord_list)
    print("luminosity history : (size ", len(lightrecord_list), ") ", lightrecord_list)
    print("temperature history : (size ", len(temprecord_list), ") ", temprecord_list)
    print("moisture history : (size ", len(moisrecord_list), ") ", moisrecord_list)

    ensemble_models(humrecord_list, lightrecord_list, temprecord_list, moisrecord_list, look_back=200)


My problem is about future timesteps prediction. In fact predictions based on testY work fine and gives decent results BUT when I ask the model to predict future timesteps that have not been measured yet, it gives me really basic results like constant values through the time for each category. 

Here's some of the result, the first part being a part of the prediction that is already measured (testY) and the last part is the last bunch of predictions of the future. 
[[78.36468   79.17461   24.37408    6.0894547]
 [78.567215  78.888794  24.395845   5.9449997]
 [78.69241   78.38314   24.393354   5.8619213]
 ...
 [79.2414    83.293686  24.013659   4.751072 ]
 [79.2414    83.293686  24.013659   4.751072 ]
 [79.2414    83.293686  24.013659   4.751072 ]]