rb_simulation.m

function rbsim = rb_simulation(model, reduced_data)
% RB_SIMULATION Perform a reduced simulation
%
%  Andreas Schmidt, 2015
tsim = tic;
[E,A,B,C,Q,R] = model.assemble(reduced_data);

if norm(E - eye(size(E))) < 1e-8
    PN = care(A,B,C'*Q*C,R);
else
    PN = care(A,B,C'*Q*C,R,0*B,E);
end

rbsim = ARE.RBSimData;
rbsim.PN = PN;
rbsim.time_sim = toc(tsim);

% Calculate the residual if necessary:
if model.calc_residual
    testim = tic;
    tResidual = tic;
    [n,nn,nnn] = model.rb_residual_norm(reduced_data, rbsim);
    tResidual = toc(tResidual);
    
    rbsim.time_residual = tResidual;
    
    rbsim.residual = n;
    rbsim.nresidual = nn;
    rbsim.nnresidual = nnn;
end

% If the error estimator is enabled, calculate the stuff:
if model.enable_error_estimator
    if isempty(reduced_data.gamma_function)
        warning('Cannot perform error estimation as no suitable gamma function was specified')
        return;
    end
    
    if ~strcmp(reduced_data.gamma_mode, model.gamma_mode)
        warning('rbmatlab:are:gamma_mode_changed', 'The gamma mode in the ReducedData class is different from the gamma mode in the Model class!');
    end
    
    tGamma = tic;
    gamma = reduced_data.gamma_function(model, reduced_data, rbsim);
    tGamma = toc(tGamma);
    
    rbsim.gamma = gamma;
    rbsim.time_gamma = tGamma;

    normBRinvBt = 1;%norm(full(B))^2*norm(inv(R));
    L = 2*reduced_data.estim.normE^2*normBRinvBt;
    rbsim.L = L;
    
    rbsim.validity_crit = 4*gamma^2*rbsim.residual*L;
    rbsim.valid = rbsim.validity_crit <= 1;
    
    rbsim.error_estimate = 1/2/gamma/L*(1 - sqrt( 1 - rbsim.validity_crit ));
    rbsim.stab_crit = rbsim.error_estimate*2*gamma*normBRinvBt*reduced_data.estim.normE/reduced_data.estim.normEinv;
    rbsim.stable = rbsim.stab_crit < 1;
    
    rbsim.time_errorestim = toc(testim);
end