riccati_gamma.m

Go to the documentation of this file.

function gamma = riccati_gamma(model, model_data, sim)
% This function calculates gamma.
%  Simdata must be a high dimensional solution.
%
%  This is done by solving the generalized Lyapunov equation
%    A' H E + E' H A = - I
%  and then gamma = || H ||
%
%  Andreas Schmidt, 2015

% Build the full dimensional system matrices:
[E,A,B,~] = riccati_assemble_system_matrices(model, model_data);

if ~isfield(model, 'gamma_mode') || strcmp(model.gamma_mode, 'lyap')
    Y = A-B*B'*sim.Z*sim.Z'*E;
    gamma = normest(lyap(Y', -eye(model.n), [], E));
elseif strcmp(model.gamma_mode, 'power_iteration')
    % This function is work in progress!!!!!
    % Initial guess for the "eigenvector"
    [x0,~,~] = eigs(A,E,1,'sm');
    x0 = x0/norm(x0);
    
    e0 = norm(x0);
    tol = 1e-8;
    
    Tend = 1e10;
    Z = sim.Z;
    
    % Outer iteration:
    for l = 1:20
        % First inner iteration. Integrate the system
        % E \dotx = (A-BB'PE)x, x(0) = E\x0
        t = 0;

        X = zeros(5000, model.n);
        dt0 = 1e-3;
        dt = dt0;
        X(1,:) = x0;
        norm0 = norm(X(1,:));
        k = 1;
        h = [];
        T = [];

        M = sparse(diag(diag(E)));
        while t<Tend
            rhs = E*X(k,:)';
            
            [x,~,~,iter,~] = bicgstab(afun(E,A,B,Z,dt), rhs, tol, 200, [], []);
            
            X(k+1, :) = x;
            
            t = t + dt;
            
            % Adjust dt:
            tau = x - X(k,:)';
            dt = dt*0.9*(1e-2/norm(tau))^0.3;

            T = [T t];
            k = k + 1;
            norm(X(k,:))
            if norm(X(k,:)) < 1e-4
                break
            end
        end
        
        % Integrate backwards to get gamma:
        
        % Perform the "backward" solution:
        K = k - 2;
        y = zeros(model.n, 1);

        for k = 1:K
            j = K+1-k;
            g = X(j,:)';
            dt = T(j+1) - T(j);
            rhs = E'*y + dt*g;
            [y,~,~,~,~] = bicgstab(afun2(E,A,B,Z,dt), rhs, tol, 200, [], [], y);
        end

        y = E'\y;
        gamma = norm(y);
        display(['Relative change ', num2str(abs(gamma-e0)/e0)])
        if abs(gamma-e0)/e0 < 1e-3
            break
        end
        
        %sim2 = riccati_simulate_system(model, model_data, 'use_feedback', sim, 'x0', zeros(model.n,1), 'use_rhs', flip(sim1.x), 'nt', 200, 'T', 10);
        %y = E'\sim2.x(end,:)';
        e0 = gamma;
        x0 = y/norm(y);
    end
end

end

function f = afun(E,A,B,Z,dt)
    f = @(x) E*x - dt*A*x + dt*B*(B'*(Z*(Z'*(E*x))));
end
function f = afun2(E,A,B,Z,dt)
    f = @(x) E'*x - dt*A'*x + dt*E'*(Z*(Z'*(B*(B'*x))));
end