function u_opt = SimulatedMultipleShooting(x0,u0,Qin,Betain,h,Q,R,xSP,P2,P3)

import casadi.*

%% Model variables

x1 = SX.sym('x1');    %fraction of oil in overflow
x2 = SX.sym('x2');    %fraction of oil in underflow

x = [x1;x2];
nx = 2;

%% System inputs

u = SX.sym('D'); %[h^-1] - Dilution Rate: F/V
du = SX.sym('du'); %input variation
s = SX.sym('s');

%% System constraints

u_min = 0;      
u_max = 1;      

du_max = 0.05;     
du0 = 0;

s_min = 0;
s_max = inf;
s0 = 0;

% control horizon
mh = 10;
% prediction horizon
ph = 15;

%% Parameters
par = initHC();
Vo = par.VRF();
Vf = par.Vol_HC-par.VRF;
Zu = 0.5;       % Constant underflow valve opening

%% Model variables

x1 = SX.sym('x1');
x2 = SX.sym('x2');
x = [x1;x2];

%% Model inputs
u = SX.sym('u'); % Overflow valve opening

cv1 = pi*0.005^2;
cv2 = pi*0.001^2;


%% Model equations

Qu = cv1*Zu*sqrt(2/par.Rho_w*(P3-par.Patm));
Qo = cv2*u*sqrt(2/par.Rho_o*(P2-par.Patm));

% Qsep polynomial approximation
p2 =  -9.447e+07;
p1 =        9024;
p0 =      0.7648;

% Oil inflow
Qino = Betain.*Qin;
Qsep = Qino.*(p2*Qo.^2+p1*Qo+p0);

xdot1 = (1/Vo)*(Qsep-x1*Qo); 
xdot2 = (1/Vf)*(Qino-Qsep-x2*Qu);

xdot = [xdot1;xdot2];
   
% Objective term
L = Q*s^2 + R*du^2;     % Soft Constraints

Uu=[u;du;s];

% Formulate discrete time dynamics
% CVODES from the SUNDIALS suite

ode = struct('x',x,'p',Uu,'ode',xdot,'quad',L);

deltaX = 10E-6;

%% Integrator
 
opts = struct('tf',h);
F = integrator('F', 'cvodes', ode, opts);

% preparing symbolic variable

w={};

% preparing numeric variables and bounds
w0 = [];
lbw = [];
ubw = [];

g={};
lbg = [];
ubg = [];
J = 0;

% Formulate the NLP
% initial state as variable

Xk = MX.sym('X0', nx);
w = {w{:}, Xk};
lbw = [lbw; x0(1); x0(2)];
ubw = [ubw; x0(1); x0(2)];
w0 = [w0; x0(1); x0(2)];

%loop for discretized dynamic
for k=0:ph-1
    
    % variable for the control in this control interval
    Uk = MX.sym(['U_' num2str(k)]);
    DUk = MX.sym(['DU_' num2str(k)]);
    Sk = MX.sym(['S_' num2str(k)]);
    Uk1=[Uk;DUk;Sk];
    w = {w{:}, Uk1};
    
    % bounds
    lbw = [lbw; u_min; -du_max; s_min];
    ubw = [ubw; u_max; du_max; s_max];
    
    %initial
    w0 = [w0;u0;du0;s0];
    
    
     if k <= mh
            if k == 0
                DUk = Uk - u0;
            else
                DUk = Uk(1) - Uk_1;
            end
            
            % Add input movement constraint
            g = {g{:}, DUk};
            lbg = [lbg; -du_max];
            ubg = [ubg; du_max];
            
      else 
            Uk(1) = Uk_1;
            DUk = 0;
      end
    
    % Integrate till the end of the interval
    
    Fk = F('x0',Xk,'p',Uk1);

    J = J + Fk.qf; 
    
    Pk = Fk.xf;
    

    
    
    Xk = MX.sym(['X_' num2str(k+1)], nx);
    
    w = {w{:}, Xk};
    lbw = [lbw; 0; 0];
    ubw = [ubw; 1; 1];
    w0 = [w0; x0(1); x0(2)];
    
 %calculating delta U
       
   
    if k ~= ph - 1
            % for computing du
            Uk_1 = Uk1(1);
    end
    % Add continuinity constraint (closing shooting gap)
    g = {g{:}, Pk-Xk};
    lbg = [lbg; zeros(nx,1)];
    ubg = [ubg; zeros(nx,1)];
    
    % Soft constraint
    g = {g{:}, (Xk(2)-Sk-xSP)};
    lbg = [lbg; -deltaX];
    ubg = [ubg; 0]; 
    
end

% Create an NLP solver
prob = struct('f', J, 'x', vertcat(w{:}), 'g', vertcat(g{:}));
solver = nlpsol('solver', 'ipopt', prob);

% Solve the NLP
sol = solver('x0', w0, 'lbx', lbw, 'ubx', ubw,'lbg', lbg, 'ubg', ubg);
w_opt = full(sol.x);

%first input
u_opt = w_opt(nx + 1);

end

