function [objvalsim,yfs,yrs] = analyze_dro(Nsim)
%%  simulation

% clear all

load trajectory_results

% load initial worst-case disturbance. Note that we use subscript 's'
% for simulation
d0s = double(d0);
xfs(:,1) = Bd*d0s; 
xrs(:,1) = xfs(:,1);	% initial states

% Do a closed-loop simulation

for i = 2:Nsim
    x0if = xfs(:,i-1);
    x0ir = xrs(:,i-1);
    % set up the controllers
    [Matricesf Matricesr] = matrices_analyze_dro(nfull,nred,N,x0if,x0ir)
    
    % calculate inputs
    ufs(:,i-1) = MPC_solv(Matricesf);
    urs(:,i-1) = MPC_solv(Matricesr);
    uf(:,i-1)  = double(Uf{i-1}(1:2));
    ur(:,i-1)  = double(Ur{i-1}(1:2));
    % update states
    % results for full-order controller:
    xfs(:,i)     = Af*xfs(:,i-1) + Bf*ufs(:,i-1);
    % results for low-order controller:
    xrs(:,i)   = Af*xrs(:,i-1) + Bf*urs(:,i-1);
    
end

% calculate outputs
for i = 1:Nsim
% full-order controller:
yfs(:,i)	 = Cf*xfs(:,i);
yf(:,i)	 = Cf*double(xf{i});

% low-oerder controller:
yrs(:,i)	 = Cf*xrs(:,i);
yr(:,i)	 = Cf*double(xr{i});
end


% plot closed loop simulation

subplot(2,2,[1:2])	% outputs
    plot(yfs(1,:),'b')
    hold
    plot(yrs(1,:),'r--')
    plot(yfs(2,:),'b')
    plot(yrs(2,:),'r--')
    legend('full','reduced')
    ylabel('Outputs')
    xlabel('Time')
subplot(2,2,[3:4])	% inputs
    plot(ufs(1,:),'b')
    hold
    plot(urs(1,:),'r--')
    plot(ufs(2,:),'b')
    plot(urs(2,:),'r--')
    ylabel('Inputs')
    xlabel('Time')
    XLIM([1 Nsim]);

objvalsim = norm((yfs(:,Nsim)-yrs(:,Nsim)),inf);

