red_slv.m

% Script for computing the reduced order models (for examples generated by 
% AREBENCH.f for the continuous-time and/or discrete-time systems from the 
% AREs benchmark collection) using SLICOT hna and MATLAB ohklmr.
%
% Switches to activate some solvers:  
% UseSLhna, UseMThna.  These should be set to 1, for activation, 
% and to 0, otherwise.  If all solvers should be used, set UseAll = 1.
% Option for plotting the Hankel singular values, setnr, for choosing
% the reduced system order, nr.  The variable setnr should be set 
% (outside this script) to 1, for plotting the Hankel singular values.

% RELEASE 2.0 of SLICOT Model and Controller Reduction Toolbox.
% Based on SLICOT RELEASE 5.7, Copyright (c) 2002-2020 NICONET e.V.
%
% Contributor:
% V. Sima, Research Institute for Informatics, Bucharest, Nov. 2004,
%          Feb. 2009.
%
% Revisions:
% -
% 

% Set the number of solvers.

if ~exist( 'nsolv', 'var' ) || isempty( nsolv ),  nsolv = 2;  end

n = size( A, 1 );  m = size( B, 2 );  p = size( C, 1 );

sys = ss( A, B, C, 0, Ts );

space = [space; '   '];  l = l + 1;
prob(l,1) = nr1;  prob(l,2) = nr2;
dims(l,1) = n;    dims(l,2) = m;    dims(l,3) = p;
dims(l,4) = NaN;  dims(l,5) = Ts;

timing(l,1:nsolv) = NaN;  
err(l,1:nsolv)    = NaN;
if UseAll == 1, err(l,1:nsolv+1) = NaN;  end

if setnr,
    
    % Set the reduced system order, nr.
    
    [ sysr, slhsv ] = hna( sys );
    close( gcf )
    set( axes, 'FontSize', 16 )
    bar( slhsv )
    title( 'Hankel singular values (of the stable part)' )
    shg
    disp( ' ' )
    disp( 'Press any key to continue' )
    pause
    nr = size( sysr.a, 1 );  first = 1;  nus = n - length( slhsv );
    while first || nr <= 0 || nr > n,
        if first,
            disp( ' ' )
            disp( [ 'The computed reduced system order is nr = ', ...
                    num2str( nr ) ] );
            disp( ' ' )
            first = 0;
        else
            disp( ' ' )
            nr = input( 'Set the reduced system order, nr = ' );
            disp( ' ' )
            if nr < nus,  disp( [ 'nr must be larger than ',   num2str( nus ) ] ),   continue,  end
            if nr > n,    disp( [ 'nr must be not larger than ', num2str( n ) ]  ),  continue,  end
        end
        disp( [ 'H-infty error norm for the ', num2str( nr ), ...
                '-th order approximation = ', ...
                num2str( 2*sum( slhsv(min(max(1,nr+1-nus),n):end) ) ) ] )
        answer = input( 'Do you want to change nr? Type y or yes (default), or n or no: ', 's' );
        if ( ~isempty( answer ) && ...
           ( answer(1) == 'n' || answer(1) == 'N' ) ),  
            break;  
        else
            nr = -1;
        end
    end
    close( gcf )
end

% SLICOT calculations.

if UseSLhna == 1 || UseAll == 1, 
    time = cputime;    
    [ sysr, slhsv ] = hna( sys, [], nr );
    nr = size( sysr.a, 1 );
    boundS = 2*sum( slhsv(min(nr+1,n):end) );
    timing(l,1) = cputime - time;  
    err(l,1)    = boundS;
    dims(l,4)   = nr;    
end 
 
% MATLAB calculations.

if UseMThna == 1 || UseAll == 1, 
    try
        time = cputime;    
        if vers7,
            [ sysm, redinfo ] = hankelmr( sys, nr );
        else
            [ sysm, boundM, mhsv ] = ohklmr( sys, 1, nr );
        end
        timing(l,2) = cputime - time;  
        if vers7,
            err(l,2) = redinfo.ErrorBound;
        else
            err(l,2) = boundM;
        end
        
        if UseAll == 1, 
            if vers7,  mhsv = redinfo.StabSV;  end
            err(l,3) = norm( slhsv - mhsv(1:length(slhsv)) )/norm(slhsv);
        end
    catch rd_slv
        disp( [ 'For example ', num2str( nr1 ), '.', num2str( nr2 ), ' ', nmfct, ...
                ' returned with error:', sprintf( '\n' ), rd_slv.message ] )
    end
end 

disp(' ')