MPCmatrices_Example_1.m

S.nx=3;
S.nu=6;
S.nd=4;

S.A=rand(S.nx,S.nx);
S.B=rand(S.nx,S.nu);
S.Gd=rand(S.nx,S.nd);
S.E=[ 1    -1    -1     0     0    -1
      0     1     0     0    -1     0 ];
S.Ed=[ 0     0     0     0
       0    -1     0     0 ];
S.umax = (2*abs(rand)+1)*ones(S.nu,1);
S.umin = -(2*abs(rand)+1)*ones(S.nu,1);
S.xmax = (2*abs(rand)+10)*ones(S.nx,1);
S.xmin = -(2*abs(rand)+10)*ones(S.nx,1);

% The following variables are loaded from the file dwn.mat:
% S (system), P (problem), D (water demand)

% Initial state:
f0 = 0.3; % (repletion factor)
x0 = S.xmin + f0*(S.xmax-S.xmin);
k=450; % (initial time instant)

P.beta=0.15;
mpc_ops = struct('use_soft_constraints',false,'use_beta',true);

% Prediction and Control Horizons   
P.Hp=24;
P.Hu=P.Hp-1;
P.alpha1= [rand;0;0;rand;0;0];
P.alpha2 = [zeros(8760,2) rand(8760,3)/100 zeros(8760,1)];
P.Wx=ones(S.nx,S.nx);
P.Wu=ones(S.nu,S.nu)*3;
P.xs=0.5*(S.xmax+S.xmin);
P.beta=0.15;
% Tolerance (for the assertions):
eps = 1e-7;

% DemandData : T-by-nd, where T=8760 (1 year, hourly)
% D : The sequence of predicted demands (here we make use of the actual
% demand data).
DemandData = rand(8760, S.nd)/100;
D=vec(DemandData(k:k+P.Hp-1,:)');

% MPC Matrices:
tic;
MPC = mpc_matrices(S, P, x0, D, k, mpc_ops);
time_mpc_matrices=toc;
disp(['MPC Optimisation problem (Formulated in ' num2str(time_mpc_matrices) ' s)...']);
disp([' -  ' num2str(length(MPC.f)) ' decision variables']);
disp([' -  ' num2str(size(MPC.F,1)) ' inequalities (generic)']);
disp([' -  ' num2str(length(MPC.ymin)) ' inequalities (box)']);
disp([' -  ' num2str(size(MPC.G,1)) ' equalities']);

% System response obtained by x+=A*x+B*u+Gd*d and X=Sx*x+Su*U+Sd*D
U=.2*rand(S.nu*(P.Hu+1),1);
X=MPC.Sx*x0+MPC.Su*U+MPC.Sd*D;
X2=zeros(P.Hp*S.nx,1);
for i=1:P.Hp
    if(i==1)
        X2(1:S.nx,1)=S.A*x0+S.B*U(1:S.nu,1)+S.Gd*D(1:S.nd,1);
    else
        if(i<P.Hu+1)
            X2((i-1)*S.nx+1:i*S.nx,1)=...
                S.A*X2((i-2)*S.nx+1:(i-1)*S.nx,1) + ...
                S.B*U((i-1)*S.nu+1:i*S.nu,1) + ...
                S.Gd*D((i-1)*S.nd+1:i*S.nd,1);
        else
            X2((i-1)*S.nx+1:i*S.nx,1)=...
                S.A*X2((i-2)*S.nx+1:(i-1)*S.nx,1) + ...
                S.B*U(P.Hu*S.nu+1:(P.Hu+1)*S.nu,1) + ...
                S.Gd*D((i-1)*S.nd+1:i*S.nd,1);
        end
    end
end


% Assertions:
assert(sum(size(MPC.Sx)==[S.nx*P.Hp, S.nx])==2, 'MPC.Sx is of wrong size.');
assert(sum(size(MPC.Sd)==[S.nx*P.Hp, S.nd*P.Hp])==2, 'MPC.Sd is of wrong size.');
assert(sum(size(MPC.Su)==[S.nx*P.Hp, S.nu*(P.Hu+1)])==2, 'MPC.Su is of wrong size.');
ne=size(S.E,1);
assert(sum(size(MPC.G)==[ne*(P.Hu+1), S.nu*(P.Hu+1)+S.nx*P.Hp])==2, 'MPC.G is of wrong size.');
assert(sum(size(MPC.gamma)==[ne*(P.Hu+1), 1])==2, 'MPC.gamma is of wrong size.');
if (mpc_ops.use_soft_constraints)
    assert(sum(size(MPC.phi)==[3*S.nx*P.Hp, 1])==2, 'MPC.phi is of wrong size.');
    assert(sum(size(MPC.F)==[3*S.nx*P.Hp, S.nx*P.Hp+S.nu*(P.Hu+1)])==2, ...
    'MPC.phi is of wrong size.');
else
     assert(sum(size(MPC.phi)==[2*S.nx*P.Hp, 1])==2, 'MPC.phi is of wrong size.');
     assert(sum(size(MPC.F)==[2*S.nx*P.Hp, S.nx*P.Hp+S.nu*(P.Hu+1)])==2, ...
    'MPC.phi is of wrong size.');
end

assert(length(MPC.Wat)==S.nu*(P.Hu+1),'MPC.Wat has incorrect size.');
assert(size(MPC.Wat,2)==1,'MPC.Wat is not a vector.');
assert(sum(size(MPC.Wxt)==[S.nx*P.Hp, S.nx*P.Hp])==2, 'MPC.Wxt is of wrond size.');
assert(sum(size(MPC.H)==[S.nu*(P.Hu+1)+S.nx*P.Hp, S.nu*(P.Hu+1)+S.nx*P.Hp])==2, ...
    'MPC.H is of wrond size.');
assert(sum(size(MPC.f)==[S.nu*(P.Hu+1)+S.nx*P.Hp, 1])==2, 'MPC.f is of wrond size.');
assert(norm(X-X2)<eps,'The response was not interpolated properly');