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 = dare(A,B,C'*Q*C,R);
else
    PN = dare(A,B,C'*Q*C,R,[],E);
end

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

% Calculate the residual if necessary:
if model.calc_residual
    tResidual = tic;
    [n,nn] = model.rb_residual_norm(reduced_data, rbsim);

    rbsim.time_residual = toc(tResidual);
    
    rbsim.residual = n;
    rbsim.nresidual = nn;
end

if model.enable_error_estimator
    tEstim = tic;
    
    if isempty(reduced_data.gamma_function)
        warning('Cannot perform error estimation as no suitable gamma function was specified')
        return;
    end
    
    % Check if the gamma function fits to the one which is set in the
    % model:
    if ~strcmp(model.RB_gamma_mode, reduced_data.gamma_mode)
        warning('The gamma mode does not fit: %s in model <-> %s in reduced data.', model.RB_gamma_mode, reduced_data.gamma_mode);
    end
    
    tGamma = tic;
    gamma = reduced_data.gamma_function(model, reduced_data, rbsim);
    tGamma = toc(tGamma);
    
    rbsim.gamma = gamma;
    rbsim.time_gamma = tGamma;
    
    alpha = 2*rbsim.gamma*rbsim.residual;
    estim = reduced_data.estim;
    
    
    A_coeff = model.A_coeff();
    B_coeff = model.B_coeff();
    E_coeff = model.E_coeff();
    [~,~,B,~,~,R] = model.assemble(reduced_data);
    
    BTB = 0;
    for i=1:length(B_coeff)
        for j=1:length(B_coeff)
            BTB = BTB + estim.BTB{i,j}*B_coeff(i)*B_coeff(j);
        end
    end
    R_X = R + B'*PN*B;
    traceATA = 0;
    for i=1:length(A_coeff)
        for j=1:length(A_coeff)
            traceATA = traceATA + estim.traceATA{i,j}*A_coeff(i)*A_coeff(j);
        end
    end
    
    BTAATW = 0;
    for i=1:length(B_coeff)
        for j=1:length(A_coeff)
            for k=1:length(A_coeff)
                BTAATW = BTAATW + estim.BTAATW{i,j,k}*B_coeff(i)*A_coeff(j)*A_coeff(k);
            end
        end
    end
    
    WTAATW = 0;
    for i=1:length(A_coeff)
        for j=1:length(A_coeff)
            WTAATW = WTAATW + estim.WTAATW{i,j}*A_coeff(i)*A_coeff(j);
        end
    end
    
    % Calculate closed loop norm:
    if strcmp(model.RB_closed_loop_norm, 'fro')
        tClosedLoopNorm = tic;
        normAXsquared = abs(traceATA-2*trace(B*inv(R_X)*BTAATW*PN) + ...
            trace(B*inv(R_X)*BTB*inv(R_X)*B'*PN*WTAATW*PN));
        rbsim.time_closed_loop_norm = toc(tClosedLoopNorm);
    else
        tClosedLoopNorm = tic;
        if ~isfield(estim, 'model_data')
            error('Cannot use 2-norm calculation')
        end
        rbrec = DARE.SimData;
        rbrec.Z = reduced_data.RB_W*real(sqrtm(rbsim.PN));
        normAXsquared = model.closed_loop_norm(estim.model_data, rbrec)^2;
        rbsim.time_closed_loop_norm = toc(tClosedLoopNorm);
    end
    
    % Validity Criterion:
    temp = alpha*max(abs(eig(inv(R_X - alpha*BTB)*BTB)));
    rbsim.validity_crit = temp*(2+temp)*normAXsquared;
    rbsim.valid = rbsim.validity_crit <= 1;
    
    % Stability Criterion
    rbsim.stab_crit = alpha/(2 - rbsim.validity_crit)*sqrt(normAXsquared)* ...
        max(abs(eig(inv(R_X - alpha/(2 - rbsim.validity_crit)*BTB)*BTB)));
    rbsim.stable = rbsim.stab_crit < sqrt(normAXsquared + 1/rbsim.gamma) - sqrt(normAXsquared);
    
    if model.RB_closed_loop_stable == true
        rbrec = DARE.SimData;
        rbrec.Z = reduced_data.RB_W*real(sqrtm(rbsim.PN));
        [clstable,cleig] = model.closed_loop_stable(estim.model_data, rbrec);
        rbsim.closed_loop_stable = clstable;
        rbsim.closed_loop_max_eig = cleig;
    end
    
    % Error estimator:
    rbsim.error_estimate = alpha/(2-rbsim.validity_crit);
    
    % Effectivity of the error estimator
    traceETE = 0;
    for i=1:length(E_coeff)
        for j=1:length(E_coeff)
            traceETE = traceETE + estim.traceETE{i,j}*E_coeff(i)*E_coeff(j);
        end
    end
    
    rbsim.time_error_estim = toc(tEstim);
end