function [UU,YY,XP,XMPC,STATUS,REGION]=mpcloop_engine(isexplicit,xp0,nxp,lastx,lastu,nx,nxQP,nv,ndp,nym,nu,ny,...
    ypoff,yoff,xpoff,xoff,dxpoff,uoff,voff,vpoff,upoff,...
    myindex,ref_signal,md_signal,ud_signal,mn_signal,un_signal,mv_signal,openloop,...
    Kx,Kr,Ku1,Kv,Kut,utarget,degrees,isQP,isPreviewREF,isPreviewMD,...
    CustomSolver,CustomSolverCodeGen,UseSuboptimalSolution,...
    QPAlgorithm,ASOptions,IPOptions,ADOptions,...
    MIQPOptions,isHyb,Mdis,Ndis,Vdis,numdis,maxdis,Mrows,...
    Hinv,Linv,Mlim,Mx,Mu1,Mv,Ac,p,A,Bu,Bv,Cm,Dvm,L,M,...
    Ap,Bup,Bvp,Bdp,Cp,Dvp,Ddp,Tf,simBar,...
    F,G,H,K,i1,i2,nr,np,xunfindex,refunfindex,uunfindex,mdunfindex,thmin,thmax,...
    uscale,yscale,vscale,Jm,SuJm,Su1,Sx,Hv)     %#ok<INUSL>
% Internal use (by "sim" command)

% Author:  A. Bemporad
% Copyright 1994-2020 The MathWorks, Inc.
% $Revision: 1.1.8.5.8.4 $

%#codegen

%% DEfine waitbar as an extrinsic function
coder.extrinsic('waitbar');
coder.extrinsic('mpcCustomOptimization'); % used by economic MPC

%% Initialize variables
XP=coder.nullcopy(zeros(Tf,nxp));
XMPC=coder.nullcopy(zeros(Tf,nx));
YY=coder.nullcopy(zeros(Tf,ny));
UU=coder.nullcopy(zeros(Tf,nu));
STATUS=coder.nullcopy(zeros(Tf,1));
REGION=coder.nullcopy(zeros(Tf,1)); % Only used by Explicit MPC

v=coder.nullcopy(zeros(nv-1,1));
d=coder.nullcopy(zeros(ndp,1));
yn=coder.nullcopy(zeros(nym,1));
un=coder.nullcopy(zeros(nu,1));
u=coder.nullcopy(zeros(nu,1));
ym=coder.nullcopy(zeros(nym,1));
y_innov = coder.nullcopy(zeros(nym,1));
xQP=coder.nullcopy(zeros(nxQP,1));

if ~isexplicit
    vseq=coder.nullcopy(zeros(nv*(p+1),1));
    rseq=coder.nullcopy(zeros(ny*p,1));
else
    vseq=coder.nullcopy(zeros(nv-1,1));
    rseq=coder.nullcopy(zeros(ny,1));
end

%% Initialize states
old_u=lastu;
xmpc=lastx;

%% Simulate in loop 
deltavoff=vpoff-voff; % delta in v_off between plant and predictive model
deltayoff=ypoff-yoff; % delta in y_off between plant and predictive model
xp=xp0-xpoff;
nCon = size(Ac,1);
iA=false(nCon,1); % active constraint vector, length = #constr
bar_time=0; % initialize status bar

for t=0:Tf-1
    
    %% save states
    XP(t+1,:)=(xp+xpoff)';          % save current Plant state (including offset)
    XMPC(t+1,:)=(xmpc+xoff)';       % save current MPC state (including offset)
    
    %% get actual disturbance and noise signals
    v=mpc_getrv(v,md_signal,t,t,0); % Actual MD signal minus offset
    d=mpc_getrv(d,ud_signal,t,t,0); % Actual UMD signal minus offset
    yn=mpc_getrv(yn,mn_signal,t,t,0); % Actual y noise signal minus offset
    un=mpc_getrv(un,un_signal,t,t,0); % Actual u noise signal minus offset
    
    %% compute actual plant output
    if nv>1
        vp=v+deltavoff; % Plant MD signal minus vpoff
    else
        vp=zeros(0,1);
    end
    y=Cp*xp+Dvp*[vp;1]+Ddp*d;   % Plant output minus ypoff

    YY(t+1,:)=(y+ypoff)';   % save current Plant output (including offset)
    
    %% compute u
    if openloop
        u=mpc_getrv(u,mv_signal,t,t,0); % External MV signal (including offset)
        status=1;
    else
        % get measured plant output (with noise)
        for i=1:nym
            j=myindex(i);
            ym(i)=y(j)+deltayoff(j)+yn(i); % Meas. output signal minus yoff
        end
        
        % update of state observer to x(k|k) from x(k|k-1) = xmpc
        y_innov = (ym-Cm*xmpc-Dvm*[v;1]);
        xkk=xmpc+M*y_innov;
        
        % Compute xQP
        for j=1:nxQP % xQP=xkk(1:nxQP)  only these first nx states are fed back to the QP problem (i.e., multiplied by the Kx gain)
            xQP(j)=xkk(j);
        end

        if ~isexplicit
            % process reference signal
            if (isPreviewREF==0)
                % Get current reference (minus yoff)
                rseq=mpc_getrv(rseq,ref_signal,t,t,0);
                % Repeat over prediction horizon
                for i=1:p-1
                    for j=1:ny
                        rseq(j+i*ny) = rseq(j);
                    end
                end
            else
                rseq=mpc_getrv(rseq,ref_signal,t+1,t+p,0);
            end
            
            % process measured disturbance signal
            if (isPreviewMD==0)
                vseq=mpc_getrv(vseq,md_signal,t,t,1);
                % Repeat over prediction horizon
                for i=1:p
                    for j=0:nv-1
                        vseq(j+i*nv+1) = vseq(j+1);
                    end
                end
            else
                vseq=mpc_getrv(vseq,md_signal,t,t+p,1);
            end
            
            % compute Bc
            if isQP
                Bc = -Mlim-Mx*xQP-Mu1*old_u-Mv*vseq;
            else
                Bc = zeros(size(Mlim));
            end
            
            % compute optimal control action zopt and update iA
            % Option to disable warm start
            if ~ASOptions.UseWarmStart
                iA=false(nCon,1);
            end
            % zopt is analytical solution: when isQP is false or when isQP is true and max iterations is 0.  status is 1.
            % zopt is 0 (freezing): when status<0 (infeasible QP) or when (status=0 && UseSuboptimalSolution is false).  isQP is true.
            [zopt, ~, status, iA] = mpc_solveQP(xQP,nCon,degrees,isQP,Kx,Kr,rseq,Ku1,...
                old_u,Kv,vseq,Kut,utarget,Linv,Hinv,Ac,Bc,iA,CustomSolver,CustomSolverCodeGen,H,UseSuboptimalSolution,QPAlgorithm,...
                ASOptions,IPOptions,ADOptions,MIQPOptions,isHyb,Mdis,Ndis,Vdis,uoff,numdis,maxdis);
        else
            %% Explicit MPC
            
            % Get current reference (minus yoff)
            rseq=mpc_getrv(rseq,ref_signal,t,t,0);
            if nv>1
                vseq=mpc_getrv(vseq,md_signal,t,t,0);
            else
                vseq=zeros(0,1);
            end
            
            openloopflag=false;
            no_mv=true;
            CustomEstimation=true; % This makes MPCBLOCK_PWA_EVALUATE use xQP as the current state estimate instead of the current measured output
            xoff=zeros(nxQP,1);   % dummy, xQP is offset-independent
            RMYscale=zeros(nym,1); % dummy
            myoff=zeros(nym,1);   % dummy
            
            % single controller
            switch_in = 1;
            switch_inport = 0;
            no_switch = 1;
            enable_value = 1;
            
            [zopt, status, region] = mpcblock_pwa_evaluate(xQP, rseq, vseq, xQP, old_u, old_u, no_mv, nu, nx, nv, uoff, ...
                openloopflag, F, G, H, K, i1, i2, nr, np, xunfindex, refunfindex, uunfindex, mdunfindex,...
                thmin,thmax, CustomEstimation, RMYscale, myoff, xoff, A, Bu, Bv, Cm, Dvm, M, L, switch_in, switch_inport, no_switch, enable_value);
            
            zopt=zopt-old_u; % To get the du signal
            
            % save current region
            REGION(t+1,:)=region;

        end
        % compute u (including offset)
        for i=1:nu
            u(i)=old_u(i)+zopt(i)+uoff(i);
        end
        % Clip the solution to satisfy MV constraints if necessary
        Suboptimal = ((status == 0 || (status == -3 && isHyb)) && UseSuboptimalSolution);
        if ~isexplicit && (isQP || isHyb) && (Suboptimal || (QPAlgorithm==int32(1) && ASOptions.MaxIterations == 0) || (QPAlgorithm==int32(3) && IPOptions.MaxIterations == 0))
            u = mpc_clipping(u, old_u+uoff, uoff, Mlim, Mrows, p, ny, degrees);
        end
        old_u = u - uoff;
    end
    
    % compute actual u used in the plant (with u noise, including offset)
    up=u+un; 
    % save current input to plant
    UU(t+1,:)=up'; 

    % save current status code
    STATUS(t+1,:)=status; 

    % update observer state x(k+1|k) from x(k|k-1)
    xmpc=A*xmpc+Bu*(u-uoff)+Bv*[v;1]+L*y_innov; % We might feed up instead of u here, depending if un is measured or not !
    
    % update plant state
    xp=Ap*xp+Bup*(up-upoff)+Bvp*[vp;1]+Bdp*d;
    
    % Update WAITBAR if required
    if simBar
        adder=(t+1)/Tf;
        if adder-bar_time>=0.4  % only allows updating waitbar up to 25 times
            bar_time=adder;
            waitbar(bar_time);
        end
    end

end