riccati_rb_simulation.m

Go to the documentation of this file.

function sim_data = riccati_rb_simulation(model, reduced_data)
%sim_data = riccati_rb_simulation(model, reduced_data)
%  This function performs a reduced basis simulation. The solution to the
%  reduced ARE is calculated by calling the function "care".
%
%  The resulting structure sim_data contains the following fields:
%    PN ...................... The reduced solution (N x N dimensional)
%    time.solve .............. Time for solving the reduced equation
%    time.estim .............. Time for the estimation procedure
%    time.residual ........... Time for the residual calculation
%    time.total .............. Total time for the whole rb simulation
%    P_residual .............. Frobenius norm of the residual
%    P_residual_normalized ... Frobenius norm of the residual, normalized
%    
%  Additionally, if model.RB_enable_estimator is true, the structure
%  estim is added to the structure, containing the following fields:
%    valid ................... 1/0 stating if the validity bound is valid
%    validity_bound .......... Value of the validity bound
%    bound ................... The error bound
%    normPN .................. Norm of the reduced solution matrix W PN W'
%    bound_relative .......... The relative error bound Delta_P(mu)/|| P ||
%
% Andreas Schmidt, 2015
sim_data = {};

% Assemble the reduced matrices
[E,A,B,C] = riccati_assemble_system_matrices(model, reduced_data);

time = {};
if norm(eye(size(A,1))-E) < 10e-6
    t = tic;
    sim_data.PN = care(A,B,C'*C);
    time.solve = toc(t);
else
    try
    R = eye(size(B,2));
    maxnm = 100;
    maxadi = 200;
    
    t = tic;
    sim_data.PN = real(care(full(A),full(B),full(C'*C),R,full(B*0),full(E)));
    time.solve = toc(t);
    
    catch e
        warning('Equation not solveable');
        display(e)
        sim_data.PN = zeros(size(A,1));
        sim_data.P_residual = -1;
        return
    end
end

t = tic;
[r,rn] = riccati_residual_norm(model, reduced_data, sim_data);
sim_data.P_residual = r;
sim_data.P_residual_normalized = rn;
time.residual = toc(t);

% Perform the estimation process
time.estim = 0;
if model.RB_enable_estimator == 1
    t = tic;
    
    gamma = riccati_rb_gamma(model, reduced_data, sim_data);
    sim_data.estim.gamma = gamma;
    
    [validity_bound, valid] = riccati_validity_bound(model, ...
        reduced_data, sim_data);
    sim_data.estim.validity_bound = validity_bound;
    sim_data.estim.valid = valid;
    
    % Calculate the norm:
    % Assemble the values for ||E||^2 and ||B||^2
    model.decomp_mode = 2;
    Ecoeff = model.E(model, reduced_data);
    if length(Ecoeff) > 1
        error('Q_E > 1 currently not implemented');
    end
    normE2 = (reduced_data.estim.normE*Ecoeff(1))^2;

    Bcoeff = model.B(model, reduced_data);
    if reduced_data.estim.normBBT_use2norm
        normBBT = (reduced_data.estim.normBBT*Bcoeff(1))^2;
    else
        normBBT = 0;
        for i = 1:length(Bcoeff)
            for j = 1:length(Bcoeff)
                for k = 1:length(Bcoeff)
                    for l = 1:length(Bcoeff)
                        normBBT = normBBT + Bcoeff(i)*Bcoeff(j)*Bcoeff(k)*Bcoeff(l)*reduced_data.estim.normBBT{i,j,k,l};
                    end
                end
            end
        end
    end
    
    epsilon = sim_data.P_residual;
    sim_data.estim.bound = gamma*epsilon/(1-4*gamma^2*epsilon*normE2*normBBT);
    sim_data.estim.bound2 = 2*gamma*epsilon;
    sim_data.estim.normPN = sqrt( trace(sim_data.PN'*sim_data.PN) );
    sim_data.estim.bound_relative = sim_data.estim.bound/sim_data.estim.normPN;
    sim_data.estim.bound_relative2 = 2*gamma*epsilon/sim_data.estim.normPN;
    time.estim = toc(t);
end
time.total = time.estim + time.solve + time.residual;

sim_data.time = time;
end