classdef BalancedTruncation < mor.GenericBTSpec & matlab.mixin.CustomDisplay
   % Specification for reducing LTI model using Balanced Truncation.
   %
   %   For non-sparse LTI models SYS, R = REDUCESPEC(SYS,'balanced') creates 
   %   a mor.BalancedTruncation spec R for reducing SYS with the Balanced
   %   Truncation method. Use R.Options to further configure this spec,
   %   type "help(R.Options)" for details.
   %
   %   R = PROCESS(R) processes the specification R and collects the
   %   following information:
   %     Sigma     Hankel singular values or equivalent. They measure the
   %               contribution of each state to the I/O map.
   %     Energy    Normalized energy of each state, measured as
   %               (Sigma(k)/Sigma(1))^2.
   %     Error     Bound on absolute or relative approximation error in the
   %               frequency domain. Error(k) is the error when setting the
   %               order to r=k-1.
   %     TL,TR     Left/right transformations
   %                 TL'*A*TR=[A1 0;0 A2],  TL'*B=[B1;B2],  C*TR=[C1 C2]
   %               where A2 (A1) collects all stable (unstable) modes.
   %     Lr,Lo     Cholesky factors of controllability and observability 
   %               gramians for the stable part (A2,B2,C2).
   %     Regularization  Regularization level when Goal='relative'.
   %
   %   Use VIEW to decide which states to keep with one of these plots:
   %     'sigma'    Bar chart of Hankel singular values and associated error
   %                bound.
   %     'energy'   Bar chart of normalized energies.
   %   You can customize these plots using Name/Value pairs, for example,
   %      view(R,'sigma',Axes=gca,Yscale='linear')
   %   See HSVOPTIONS for available options besides 'Axes'.
   %
   %   Use GETROM to obtain reduced-order models. Select the order with
   %   these options:
   %     Order      Desired order (scalar or vector for batch selection).
   %     MaxError   Maximum approximation error. This picks the lowest order
   %                for which the error stays below this value.
   %     MinEnergy  States with lower energy are discarded.
   %   and eliminate states with this option:
   %     Method     State elimination method ('matchDC' or 'truncate').
   %                Specifies how to eliminate weakly contributing states.
   %                The "matchDC" method (default) preserves the DC gain
   %                (steady-state value of step response) and tends to better
   %                match the time response. The "truncate" method tends to
   %                better match the frequency response and is recommended
   %                with relative-error control.
   %
   %   Example: Reduce a 30th-order state-space model using balanced 
   %   truncation with relative error control and a max error of 5%:
   %      sys = rss(30);
   %      R = reducespec(sys,'balanced');
   %      R.Options.Goal = 'relative'
   %      % Run algorithm once for efficiency (optional)
   %      R = process(R)
   %      % View Hankel singular values
   %      view(R)
   %      % Get reduced-order model
   %      rsys = getrom(R,MaxError=0.05)
   %      sigma(sys,sys-rsys)
   %
   %   See also REDUCESPEC, MOR.PROCESS, MOR.VIEW, MOR.GETROM, HSVOPTIONS.

   % Copyright 2020-2023 The MathWorks, Inc.

   % Note: In dense case, the cached Lr and Lo are not factors of the original
   % gramians because of preliminary STABSEP scaling and transformation. Since
   % they are nonessential, current design just leaves them out.
   properties (SetAccess=protected, Dependent)
      TL
      TR
      Regularization
   end

   properties (Access=protected)
      TL_
      TR_
   end

   methods

      function R = BalancedTruncation(sys)
         % Convert to @ss
         try
            R.Model_ = ss(sys);
         catch
            error(message('Control:transformation:ROM1'))
         end
         R.Options_ = mor.BalancedTruncationOptions;
      end

      function value = get.TL(this)
         % Note: Cache uses TL*A*TR but BalancedTruncation uses TL'*A*TR
         % for consistency with ModalTruncation.
         if isempty(this.Cache_)
            value = [];
         else
            value = this.TL_;
         end
      end

      function value = get.TR(this)
         if isempty(this.Cache_)
            value = [];
         else
            value = this.TR_;
         end
      end

      function value = get.Regularization(this)
         if isempty(this.Cache_) || ~strcmp(this.Options.Goal,'relative')
            value = [];
         else
            value = this.Cache_.Regularization;
         end
      end

      function R = process(R)
         % Recompute derived quantities
         if isempty(R.Cache_)
            try
               OPT = validate(R.Options_,R.Model_);
               R.Cache_ = runBT(R.Model_,OPT);
               [T,R.TR_] = getBlockDiagTransform(R.Cache_);
               % Note: Cache uses TL*A*TR but BalancedTruncation uses TL'*A*TR
               % for consistency with ModalTruncation.
               R.TL_ = T';
            catch ME
               throw(ME)
            end
         end
      end

      function [rsys,varargout] = getrom(R,arg2,optArgs)
         %GETROM  Obtain reduced-order model (ROM).
         %
         %   RSYS = GETROM(R) returns a simplified model RSYS where all states 
         %   associated with numerically zero Hankel singular values have been
         %   removed. This amounts to a minimal realization of the original 
         %   system SYS.
         %
         %   RSYS = GETROM(R,Name1=Value1,Name2=Value2,...) returns a reduced-order
         %   model RSYS based on the specified options. You can select the order
         %   using one of the following criteria:
         %     Order      Desired order (scalar or vector for batch selection).
         %     MaxError   Maximum approximation error. This picks the lowest order
         %                for which the error does not exceed this value.
         %     MinEnergy  This discards all states with lower normalized energy.
         %   You can also specify the following option:
         %     Method     State elimination method ('matchDC' or 'truncate').
         %                Specifies how to eliminate the weakly contributing states
         %                The "MatchDC" method (default) preserves the DC gain
         %                (steady-state value of step response) and tends to better
         %                match the time response. The "Truncate" method tends to
         %                better match the frequency response and is recommended
         %                with relative-error control.
         %
         %   [RSYS,INFO] = GETROM(R,...) also returns a struct INFO containing the
         %   left and right projectors PL and PR that map the original matrices
         %   (A,B,C,D) to the reduced-order matrices (PL'*A*PR,PL'*B,C*PR,D).
         %   when using the 'truncate' method.
         %
         %   See also REDUCESPEC.
         arguments
            R
            arg2 = '';
            optArgs.?mor.getrom.BalancedTruncationOptions;
         end

         % Handle help request getrom(R,'-help')
         if nargin==2
            if nargout>0 || ~strcmp(arg2,'-help')
               error(message('Control:transformation:getrom2'))
            end
            help mor.BalancedTruncation.getrom
            return
         end

         % Update cache, get options, and resolve order
         try
            R = process(R);
            ROMopt = initOptions(mor.getrom.BalancedTruncationOptions,namedargs2cell(optArgs));
            MatchDC = strcmp(ROMopt.Method,'matchDC');
         catch ME
            throw(ME)
         end
         orders = resolveOrder(ROMopt,R.Cache_);

         if MatchDC
             warnMatchDCIntervals(R.Options_)
         end
        
         % Get reduced model(s).
         [rsys,varargout{1:nargout-1}] = reduceBT(R.Model_,orders,R.Cache_,MatchDC);
         % Clear notes, userdata, etc
         rsys.Name = '';  rsys.Notes = {};  rsys.UserData = [];
      end

      % BALRED support
      function Info = getBalredInfo(R)
         % Get data structure formerly returned by HSVD/BALRED
         C = R.Cache_;
         if ~(isempty(R.Options.FreqIntervals) && isempty(R.Options.TimeIntervals))
            % BALRED returned NaN errors with these options
            C.Error(:) = NaN;
         end
         Info = ltipack.balredInfo(C);
      end

      function R = setBalredInfo(R,Info)
         % Info = ltipack.balredInfo
         R.Cache_ = strip(Info);
      end

   end

   methods (Hidden)

      function [okpole,gaps,dcgaps] = qeCheckROM(R,w,Sensitivity,ElimMethod,orders)
         % For QE only.
         sys = R.Model_;  opt = R.Options_;
         if nargin<5
            orders = (0:order(sys))';
         end
         % Determine number of unstable poles
         p = pole(sys);
         if isct(sys)
            nns = sum(real(p)>=-opt.Offset*max(1,abs(imag(p))));
         else
            nns = sum(abs(p)>=1-opt.Offset);
         end
         % Reduce and compare
         R = process(R);
         rsys = getrom(R,Order=orders,Method=ElimMethod);
         balInfo = ltipack.balredInfo(R.Cache_);
         [gaps,nnsr] = qeCheckModelAppx(sys,rsys,w,Sensitivity,opt,balInfo);
         nr = order(rsys);
         for ct=1:numel(nr)
            [~,~,~,~,e] = dssdata(rsys,ct);
            if ~isempty(e)
               % Watch for DAE form
               nr(ct) = nr(ct)-sum(all(e==0,1));
            end
         end
         okpole = [all(nr<=max(nns,orders(:))) ; (max(nnsr)<=nns)];
         % DC gaps for MatchDC
         if nns==0 && strncmp(ElimMethod,'m',1)
            dc = dcgain(sys);
            atol = 1e-12*norm(dc,1);
            dcgaps = qeCompareDC(rsys,sys,Sensitivity(1),atol);
         else 
            dcgaps = 0;
         end
      end
   end

   methods (Access = protected)
      % Reorder properties
      function propgrp = getPropertyGroups(~)
         proplist = {'Sigma','Energy','Error','TL','TR','Lr','Lo','Regularization','Options'};
         propgrp = matlab.mixin.util.PropertyGroup(proplist);
      end
   end

end

