function [Mu, Mv, Mu1, Mx, Mlim] = mpc_customconstraintcoef(SuJm, Sx, Su1, Hv, C, Dv, Jm, nCC, E, F, S, G, mvoff, mdoff, yoff)
% MPC_CONSTRAINTCOEF compute output prediction model coefficients and
% constraint coefficients that depend on the prediction model. 
%
%   Inputs:
%       A, Bu, Bv, C, Dv: prediction model parameters, always 3-D signals
%       Jm: constant matrix
%       nCC, E, F, S, G: number of custom constraints and parameters
%       mvoff, mdoff, yoff: scaled offsets of MV, MD and OV signals
%       isLTV:; whether LTV mode is on (model varying at prediction steps)
%
%   Note: "S" and "mdoff" are ignored when MD does not exist (i.e. nv==1).

% Author: A. Bemporad, N.L. Ricker, R. Chen
% Copyright 1986-2016 The MathWorks, Inc.

%#codegen

% Notation:
%   Prediction model:
%       y = Su*du + Hv*v + Sx*x0 + Su1*u1
%   where:
%       y     p*ny, future outputs at intervals k+1:k+p
%       du    p*nu, MV increments, k:k+p-1
%       v     (p+1)*nv, current and future MD values, k:k+p
%       x0    nx, current states
%       u1    nmv, MV values used over the previous sampling instant.
%       Su, Hv, Sx, Su1 are constant model coefficients.
%
%   Relationship between MV increments and optimization decisions:
%       du = Jm*z
%   where:
%       z     blocked MV increments (actual decision variables for QP)
%       Jm    constant matrix mapping dum to du.
%
%   Relationship between MV increments and MVs:
%       u = I2*du + I1*u1 = I2*Jm*z + I1*u1
%   where:
%       u     p*nu MV values, k:k+p-1
%       I1, I2 are constant matrices
%
%   Ultimately, we will solve (at each control interval)
%       Min f'*z + 0.5*z'*H*z
%           s.t. Ac*z >= b
%   where:
%       z1 = [z; e]  includes the slack variable (for constraint softening)
%       Vectors f and b depend on the controller and plant state at each interval (i.e., v, x0, u1)
%       Matrices H and Ac depend on the prediction model and are constant unless the model changes
%
%   Standard bound constraints:
%       -y + Vymax*e  >= -ymax
%        y + Vymin*e  >=  ymin
%       -u + Vumax*e  >= -umax
%        u + Vumin*e  >=  umin
%      -du + Vdumax*e >= -dumax
%       du + Vdumin*e >=  dumin
%   Custom constraints:
%       -E*u(k) - F*y(k)- S*v(k) + V*e >= -G
%   where k = 0:p with the stipulation u(p) = u(p-1).
%
% Some or all constraints may be irrelevant, e.g., redundant or with bounds
% at -Inf.  A key assumption is that the QP constraints retain their
% identity but their coefficients may change if the prediction model is
% updated.

% Prediction model coefficients Su*Jm, Su1, Sx, Hv.  Special handling for
% the case in which v is constant over the prediction horizon.
ny = size(C,1);     % number of OV
nx = size(C,2);     % number of augmented states (nxp+nxdist+nxnoise)
nu = length(mvoff); % number of MV
nv = size(Dv,2);    % number of MD + 1 (last column is for DX offset)
p = size(Jm,1)/nu;  % prediction horizon

% If there are mixed constraints, construct the QP coefficient matrix
% contributions.
if nCC > 0
    % Extend matrices for k=0:p (p+1 steps)
    SuJmp1 = [zeros(ny,size(SuJm,2),'like',C); SuJm];  % du doesn't affect y(0)
    Su1p1 =  [zeros(ny,nu,'like',C); Su1];
    Sxp1 =   [C(:,:,1); Sx];
    Hvp1 = [Dv(:,:,1), zeros(ny,p*nv,'like',C); Hv];
    I1     = kron(ones(p+1,1,'like',C),eye(nu,'like',C));
    I2     = kron(tril(ones(p+1,'like',C)),eye(nu,'like',C));
    I2Jmp1 = I2*[Jm; zeros(nu,size(Jm,2),'like',C)];
    
    Ep1    = kron(eye(p+1,'like',C), E);
    Fp1    = kron(eye(p+1,'like',C), F);
    % Add a column of zeros to each block in Sp1 to account for the DX
    % offset term, which is the last v and doesn't contribute to the
    % mixed constraints.
    if nv>1
        Sp1 = kron(eye(p+1,'like',C), [S, zeros(nCC,1,'like',C)]);
    else
        Sp1 = kron(eye(p+1,'like',C), zeros(nCC,1,'like',C));
    end
    % Calculate Mu, Mu1, Mv, and Mx for k = 0:p.
    Mu  =  Ep1*I2Jmp1 + Fp1*SuJmp1;
    Mv  = -Fp1*Hvp1 - Sp1;
    Mx  = -Fp1*Sxp1;
    Mu1 = -Ep1*I1 - Fp1*Su1p1;
    % Mlim accounts for G and effect of offsets
    Mlim1 = E*mvoff(:) + F*yoff(:) - G;
    if nv > 1
        Mlim1 = Mlim1 + S*mdoff(:);
    end
    Mlim = -reshape(Mlim1(:,ones(1,p+1,'like',C)), nCC*(p+1), 1);
else
    Mu = zeros(1,size(SuJm,2),'like',C);
    Mv = zeros(1,nv,'like',C);
    Mx = zeros(1,nx,'like',C);
    Mu1 = zeros(1,nu,'like',C);
    Mlim = zeros(1,1,'like',C);
end
