I am trying to apply the predictive control model using an already trained neural network model and SciPy for optimization. My model takes a window of past command + noise + system output to predict the system output. Then, I use these predictions to calculate the objective function. My question is: I have written the code, and it works, but the optimizer does not provide the necessary commands for my output to stay within the defined range. Here is the code and the results provided by the solver

this is my code:

#%% optimal control 
from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint

def mpc_cost(U_optimal,noise_window, u_window, y_window, window,P, M ,sp,ny, nu,a, b,  s1, s2, model_multi): 
        
        """
    
        Args:
        U_optima : optimal control obtained by MPC
        noise_window : noise 
        y_window : output window
        M and P : are the horizon of control and prediction
        ny : nombre of output, nu nombre of coontrol input
        model_lstm : trained model Returns:Predicted output.y(k+1),...,y(k+p)
        
        return : cost function
        """
        U_optimal = np.reshape(U_optimal,(M,1))
        noise_window = np.reshape(noise_window, (window, nu))
        u_window = np.reshape(u_window, (window, nu))
        input_window = np.concatenate((noise_window,u_window ),axis = 1)
        
        # Prepare the input sequence for prediction
        y_window = np.reshape(y_window,(window,ny))
        X_seq = np.concatenate((input_window, y_window), axis=1)
        
        X_seq = s1.transform(X_seq)
        X_seq = np.expand_dims(X_seq, axis=0)
        
        # Predict with the model
        y_pred = model_multi.predict(X_seq)
        y_pred = np.reshape(y_pred , (-1, 1))
        
        # denormalize the predicted output
        y_pred = s2.inverse_transform(y_pred)
        
        
        # calcul incremental control moves delta_u
        
        U_optimal = np.reshape(U_optimal,(M,1))
        
        # add a last values of u for calculate delta_u
        d_u = np.append(u_window[-1], U_optimal) 
        d_u = d_u.reshape((-1,1))
        SP = np.ones((P, ny)) *sp
        
        # cost function 
        W_CV = np.array([a]) 
        W_MV = np.array([b]) 
         
        pred_nn = {}
        pred_nn["y_hat_multi"] = y_pred

    
        Obj = np.sum(((pred_nn["y_hat_multi"] - SP ) ** 2).dot(W_CV)) + np.sum(
            ((d_u[1:] - d_u[0:-1]) ** 2).dot(W_MV))
        
        
        return Obj.item()
    
# MPC calculation
                                                                   
def mpc_solver(U_optimal, noise_window, u_window, y_window, window, P, M, sp,ny, nu,a,b, s1, s2, model_multi):
    
    U_optimal = U_optimal.flatten()
    bnds = [(0, 1)] * len(U_optimal)
    
    solution = minimize(mpc_cost, U_optimal, args=(noise_window, u_window, y_window, window,P, M ,sp,ny, nu,a, b, s1, s2, model_multi),
                        method='SLSQP', 
                        bounds=bnds, 
                        options={'eps': 1e-04, 'maxiter': 100,'ftol': 1e-6, 'disp': True})

    u = solution.x
    print("Message retourné par minimize :", solution.message)

    u = np.reshape(u, (M, 1))

    return u

def pred_output(noise_window, u_window, y_window, window,P, M ,sp,ny, nu,a, b,  s1, s2, model_multi): 
        
        """
        Prédit la sortie du modèle LSTM.
        
        Args:
            
        noise_window : noise window
        u_window : control window
        y_window : output window
        M and P : are the horizon of control and prediction
        ny : nombre of output, nu nombre of coontrol input
        model_lstm : trained model Returns:Predicted output.y(k+1),...,y(k+p)
        
        """
        noise_window = np.reshape(noise_window, (window, 1))
        u_window = np.reshape(u_window, (M, 1))
        input_window = np.concatenate((noise_window,u_window),axis = 1)
        
        # Prepare the input sequence for prediction
        y_window = np.reshape(y_window,(window,ny))
        X_seq = np.concatenate((input_window, y_window), axis=1)
        
        X_seq = s1.transform(X_seq)
        X_seq = np.expand_dims(X_seq, axis=0)
        
        # Predict with the model
        y_pred = model_multi.predict(X_seq)
        y_pred = np.reshape(y_pred , (-1, 1))
        
        # Inverse transform the predicted output
        y_pred = s2.inverse_transform(y_pred)
        
        return y_pred


    


#%% execution 
# nbr of control,output,weight for cost function
nu = 1
ny = 1
a = 1
b = 1

# control, prediction horizon
M = window
P = 10

# setpoint
sp = 17

#initial guess
U_optimal = np.zeros((M,1))

#define the : noise and the control,output window
noise = np.reshape(jours200_Tex.iloc[:],(-1,1))                                
u_window = np.reshape(jours200_beta.iloc[:][0:5],(-1,1))
y_window  = np.reshape(jours200_Tint.iloc[:][0:5],(-1,1))

# --------Simulatiopn -------------------
for i in range(0,30):
    
    noise_w = noise[i:i+window]
    u_w = u_window[i:i+window]
    y_w = y_window[i:i+window]
    
    U_mpc = mpc_solver(U_optimal, noise_w, u_w, y_w, window, P, M, sp,ny, nu,a,b, s1, s2, model_lstm)
    
    first_u = U_mpc[0]
    u_window = np.append(u_window,first_u)
    
    y_pred = pred_output(noise_w, u_w, y_w, window,P, M ,sp,ny, nu,a, b,  s1, s2, model_lstm)
    
    first_y = y_pred[0]
    y_window = np.append(y_window,first_y)
    
    U_optimal = U_mpc

#%% plot the result 
y1 = np.ones(len(y_window)) * 16
y2 = np.ones(len(y_window)) * 18

# Tracer y0 avec les lignes horizontales y1 et y2
plt.figure(1)
plt.plot(y_window, label='ypred')
plt.plot(y1, 'r--', label='lower bnd')  
plt.plot(y2, 'r--', label='upper bnd') 
plt.xlabel('Time')
plt.ylabel('T(c°)')
plt.legend()


# Tracer jours200_beta
plt.figure(2)
plt.plot(u_window, label='control optimal')
plt.xlabel('Time')
plt.ylabel('control optimal')
plt.legend()
# Afficher les graphiques
plt.show()