Model.m

classdef Model < ARE.Model
% MODEL
%   Main class for the DARE model

    properties
        RB_closed_loop_norm = 'fro';
        RB_closed_loop_stable = false;
    end

    methods
        
        function pt = problem_type(this)
            pt = 'DARE';
        end
        
        function g = gamma(this, model_data, dsim)
            calculator = DARE.GammaCalculation.Lyapunov();
            g = calculator.calculate(this, model_data, dsim);
        end
        
        function [n] = closed_loop_norm(this, model_data, dsim)
            % CLOSED_LOOP_NORM Implementation of an efficient algorithm for
            % the approximation of the 2-norm of the closed-loop matrix.
            [E,A,B,~,~,R] = this.assemble(model_data);
            
            M1 = A'*B;
            Rxinv = inv(R + (B'*dsim.Z)*(dsim.Z'*B));
            M2 = Rxinv*(B'*dsim.Z)*(dsim.Z'*A);
            M3 = (A'*dsim.Z)*(dsim.Z'*B)*Rxinv;
            BtB = B'*B;
            n = sqrt(eigs( @(x) A'*(A*x) - M1*(M2*x) - M2'*(M1'*x) + ...
                M3*(BtB*(M3'*x)), this.n, 1));
        end
        
        function [stable,cleig] = closed_loop_stable(this,model_data,dsim)
            % CLOSED_LOOP_STABLE Implementation of an efficient algorithm
            % for checking whether the closed loop matrix gives a stable
            % system. This is done by calculting the maximum eigenvalue of
            % the closed loop system.
            [E,A,B,~,~,R] = this.assemble(model_data);
            Rxinv = inv(R + (B'*dsim.Z)*(dsim.Z'*B));
            K = Rxinv*(B'*dsim.Z)*(dsim.Z'*A);
            [EL,EU] = lu(E);
            cleig = eigs( @(x) EU\(EL\(A*x - B*(K*x))), this.n, 1 );
            stable = abs(cleig) < 1;   
        end
            
    end
    
end