riccati_residual_norm.m

function [res, nres] = riccati_residual_norm(model, reduced_data, ...
    sim_data)
%function res = riccati_residual_norm(model, reduced_data, sim_data, full)
%
%  This function calculates the residual norm in a fast way by applying
%  the online/offline Frobenius norm technique as explained in the paper.
%
%  Outputs:
%    res ................ The Frobenius norm of the Residual
%    nres ............... Normalized Residual norm
%    
%  Andreas Schmidt, 2015
%

estim = reduced_data.estim;
PN = sim_data.PN;

% Calculate the coefficients
model.decomp_mode = 2;
A_coeff = model.A(model, reduced_data);
B_coeff = model.B(model, reduced_data);
C_coeff = model.C(model, reduced_data);
E_coeff = model.E(model, reduced_data);

% Assemble all matrices:
R_ECCA = zeros(size(PN,1));
for i = 1:length(A_coeff)
    for j = 1:length(C_coeff)
        for k = 1:length(C_coeff)
            for l = 1:length(E_coeff)
                R_ECCA = R_ECCA + A_coeff(i)*C_coeff(j)*C_coeff(k)*...
                    E_coeff(l)*estim.M1{i,j,k,l};
            end
        end
    end
end

R_EE = zeros(size(PN,1));
for i = 1:length(E_coeff)
    for j = 1:length(E_coeff)
        R_EE = R_EE + E_coeff(i)*E_coeff(j)*estim.M8{i,j};
    end
end

R_AA = zeros(size(PN,1));
for i = 1:length(A_coeff)
    for j = 1:length(A_coeff)
        R_AA = R_AA + A_coeff(i)*A_coeff(j)*estim.M5{i,j};
    end
end

R_BB = zeros(size(PN,1));
for i = 1:length(B_coeff)
    for j = 1:length(B_coeff)
        R_BB = R_BB + B_coeff(i)*B_coeff(j)*estim.M2{i,j};
    end
end

R_EA = zeros(size(PN,1));
for i = 1:length(A_coeff)
    for j = 1:length(E_coeff)
        R_EA = R_EA + A_coeff(i)*E_coeff(j)*estim.M3{i,j};
    end
end

R_ECCE = zeros(size(PN,1));
for i = 1:length(C_coeff)
    for j = 1:length(C_coeff)
        for k = 1:length(E_coeff)
            for l = 1:length(E_coeff)
                R_ECCE = R_ECCE + C_coeff(i)*C_coeff(j)*E_coeff(k)*...
                    E_coeff(l)*estim.M4{i,j,k,l};
            end
        end
    end
end

d = 0;
res = 0;
for i = 1:length(C_coeff)
    for j = 1:length(C_coeff)
        for k = 1:length(C_coeff)
            for l = 1:length(C_coeff)
                tmp = estim.M7{i,j,k,l}*C_coeff(i)*C_coeff(j)*...
                    C_coeff(k)*C_coeff(l);
                res = res + tmp;
                d = d + tmp;
            end
        end
    end
end

res = res + 4*trace(R_ECCA*PN) - 4*trace(R_EE*PN*R_BB*PN*R_EA*PN) + ...
    -2*trace(R_ECCE*PN*R_BB*PN) + trace(R_EE*PN*R_BB*PN*R_EE*PN*R_BB*PN) ...
    + 2*trace(R_EE*PN*R_AA*PN) + 2*trace(R_EA'*PN*R_EA'*PN);

res = sqrt(abs(res));
nres = res/sqrt(d);