classdef (Hidden) SystemArray < DynamicSystem & ltipack.ModelArray
   % System Array Interface (abstract class).
   %
   % @SystemArray implements the notion of system arrays as arrays of
   % LTIPACK data containers and translates operations on system arrays
   % into operations on the underlying LTIPACK arrays.
   
%   Author(s): P. Gahinet
%   Copyright 1986-2014 The MathWorks, Inc.
   
   %% DATA ABSTRACTION INTERFACE
   methods (Access = protected)
      
      %% MODEL CHARACTERISTICS
      boo = isstatic_(sys,~)
      boo = hasdelay_(sys,~)
      boo = hasInternalDelay_(sys,~)
      boo = isstable_(sys,~)
      [boo,sys] = isproper_(sys,ExplicitFlag,varargin)
      
      function ns = order_(sys)
         % Returns order of each array entry
         D = sys.Data_;
         nD = numel(D);
         if nD==1
            ns = order(D);
         else
            ns = zeros(size(D));
            for ct=1:numel(D)
               ns(ct) = order(D(ct));
            end
         end
      end
      
      function n = normh2_(sys)
         % H2 norm
         D = sys.Data_;
         n = zeros(size(D));
         for ct=1:numel(D)
            n(ct) = normh2(D(ct));
         end
      end
      
      function [n,fpeak] = norminf_(sys,tol,fBand,NeedStable)
         % Linf and Hinf norms
         D = sys.Data_;
         n = zeros(size(D));
         fpeak = zeros(size(D));
         for ct=1:numel(D)
            [n(ct),fpeak(ct)] = norminf(D(ct),tol,fBand,NeedStable);
         end
      end
      
      function [Index,fIndex] = sectorbnd_(sys,M0,W1,W2,tol,fBand)
         % Conic sector bound sys' * M * sys < 0
         D = sys.Data_;
         Index = zeros(size(D));
         fIndex = zeros(size(D));
         for ct=1:numel(D)
            [Index(ct),fIndex(ct)] = sectorbnd(D(ct),M0,W1,W2,tol,fBand);
         end
      end

      function [Index,fIndex] = ifpofp_(sys,IndexType,tol,fBand)
         % Computes IFP or OFP passivity index
         D = sys.Data_;
         Index = zeros(size(D));
         fIndex = zeros(size(D));
         for ct=1:numel(D)
            [Index(ct),fIndex(ct)] = ifpofp(D(ct),IndexType,tol,fBand);
         end
      end
      
      %% MODEL CONVERSIONS
      function [a,b,c,d,Ts] = ssdata_(sys,varargin)
         % Quick access to explicit state-space data
         % Select single model
         Data = localSelectSingleModel(sys.Data_,varargin);
         % Convert to state-space and get data
         Data = Data.ss();
         [a,b,c,d] = getABCD(Data);
         Ts = Data.Ts;
      end
      
      function [a,b,c,d,e,Ts] = dssdata_(sys,varargin)
         % Quick access to descriptor state-space data
         Data = localSelectSingleModel(sys.Data_,varargin);
         Data = Data.ss();
         [a,b,c,d,e] = getABCDE(Data);
         Ts = Data.Ts;
      end
      
      function [a,b,c,d,e,Ts] = sparssdata_(sys,varargin)
         % Quick access to SPARSS data
         Data = localSelectSingleModel(sys.Data_,varargin);
         Data = Data.sparss();
         [a,b,c,d,e] = getABCDE(Data);
         Ts = Data.Ts;
      end
      
      function [m,c,k,b,f,g,d,Ts] = mechssdata_(sys,varargin)
         % Quick access to MECHSS data
         Data = localSelectSingleModel(sys.Data_,varargin);
         Data = Data.mechss();
         [m,c,k,b,f,g,d] = getMCKBFGD(Data);
         Ts = Data.Ts;
      end
      
      function [num,den,Ts,sdnum,sdden] = tfdata_(sys,varargin)
         % Quick access to transfer function data
         Data = localSelectSingleModel(sys.Data_,varargin);
         Data = Data.tf();
         num = Data.num;  den = Data.den;  Ts = Data.Ts;
         sdnum = []; sdden = []; % default for systems with no covariance info
      end
      
      function [z,p,k,Ts,covz,covp,covk] = zpkdata_(sys,varargin)
         % Quick access to ZPK data
         Data = localSelectSingleModel(sys.Data_,varargin);
         Data = Data.zpk();
         z = Data.z;  p = Data.p;  k = Data.k;   Ts = Data.Ts;
         covz = []; covp = []; covk = [];
      end
      
      function [resp,freq,Ts,covresp] = frdata_(sys,varargin)
         % Quick access to FRD data
         Data = localSelectSingleModel(sys.Data_,varargin);
         Data = Data.frd();  % for GENFRD/UFRD
         resp = Data.Response;  
         freq = funitconv('rad/TimeUnit',sys.FrequencyUnit,sys.TimeUnit) * ...
            Data.Frequency;
         Ts = Data.Ts;
         covresp = []; % default implementation for response covariance data
      end
      
      function [Kp,Ki,Kd,Tf,Ts] = piddata_(sys,varargin)
         % Quick access to PID data
         Data = localSelectSingleModel(sys.Data_,varargin);
         try
            Data = Data.pid();
         catch ME %#ok<*NASGU>
            error(message('Control:ltiobject:piddata1'))
         end
         Kp=Data.Kp; Ki=Data.Ki; Kd=Data.Kd; Tf=Data.Tf; Ts=Data.Ts;
      end
      
      function [Kp,Ti,Td,N,Ts] = pidstddata_(sys,varargin)
         % Quick access to PID data
         Data = localSelectSingleModel(sys.Data_,varargin);
         try
             Data = Data.pidstd();
         catch ME
             error(message('Control:ltiobject:pidstddata1'))
         end
         Kp=Data.Kp; Ti=Data.Ti; Td=Data.Td; N=Data.N; Ts=Data.Ts;
      end
      
      function [Kp,Ki,Kd,Tf,b,c,Ts] = piddata2_(sys,varargin)
          % Quick access to 2-dof PID data
          Data = localSelectSingleModel(sys.Data_,varargin);
          try
              Data = Data.pid2();
          catch ME %#ok<*NASGU>
              error(message('Control:ltiobject:piddata21'))
          end
          Kp=Data.Kp; Ki=Data.Ki; Kd=Data.Kd; Tf=Data.Tf; b=Data.b; c=Data.c; Ts=Data.Ts;
      end
      
      function [Kp,Ti,Td,N,b,c,Ts] = pidstddata2_(sys,varargin)
          % Quick access to 2-dof PID data
          Data = localSelectSingleModel(sys.Data_,varargin);
          try
              Data = Data.pidstd2();
          catch ME
              error(message('Control:ltiobject:pidstddata21'))
          end
          Kp=Data.Kp; Ti=Data.Ti; Td=Data.Td; N=Data.N; b=Data.b; c=Data.c; Ts=Data.Ts;
      end
      
      %% BINARY OPERATIONS
      function [sys,SingularFlag] = connect_(sys,k,feedin,feedout,iu,iy,Options)
         % CONNECT(SYS,K,...) for SYS of type *ss or *frd
         SingularFlag = false;
         Data = sys.Data_;
         for ct=1:numel(Data)
            [Data(ct),warnflag] = connect(Data(ct),k,feedin,feedout,iu,iy,Options);
            SingularFlag = SingularFlag || warnflag;
         end
         sys.Data_ = Data;
         sys.IOSize_ = [length(iy) length(iu)];
         % Metadata
         sw = ctrlMsgUtils.SuspendWarnings('Control:ltiobject:DuplicatedIONames');
         sys = parenReferenceMetaData(sys,iy,iu);
         sys.Name_ = '';
     end
      
      function sys = removeInputBranching_(sys,I,J)
         Data = sys.Data_;
         for ct=1:numel(Data)
            Data(ct) = elimBranching(Data(ct),I,J);
         end
         sys.Data_ = Data;
      end
      
      %-----------------------------------------------------------
      function sys1 = times_(sys1,sys2,ScalarFlags)
         % Generic implementation of SYS1.*SYS2 for two systems
         % of the same type. Can be overloaded by subclasses.
         % NOTE: This function should not modify SYS.IOSize_
         [sys1,sys2] = matchArraySize(sys1,sys2);   % must come first
         [sys1,sys2] = matchAttributes(sys1,sys2);  % overloadable
         % Combine data
         Data1 = sys1.Data_;  Data2 = sys2.Data_;
         for ct=1:numel(Data1)
            Data1(ct) = times(Data1(ct),Data2(ct),ScalarFlags);
         end
         sys1.Data_ = Data1;
         if ScalarFlags(1)
            sys1.IOSize_ = sys2.IOSize_;
         end
      end
      
      %% ANALYSIS
      [g,dceq] = dcgain_(sys)
      fresp = evalfr_(sys,s)
      [h,SingularWarn,covH] = freqresp_(sys,w)
      s = allmargin_(sys,opt)
      fb = bandwidth_(sys,drop)
      [m,p,w,FocusInfo,sdm,sdp] = magphaseresp_(sys,grade,wspec)
      [re,im,w,FocusInfo,sdre,sdim] = nyquistresp_(sys,wspec)
      varargout = loopmargin_(sys,istr)
                  
      function [sv,w,FocusInfo] = sigmaresp_(sys,type,wspec)
         % Computes singular value response
         [sv,w,FocusInfo] = sigmaresp(sys.Data_,type,wspec,false);
      end
      
      function [indx,w,FocusInfo,InfFlag] = sectorresp_(sys,M0,W1,W2,wspec)
         % Computes sector index vs frequency
         [indx,w,FocusInfo,InfFlag] = sectorresp(sys.Data_,M0,W1,W2,wspec,false);
      end
      
      function [indx,w,FocusInfo,InfFlag] = ifpofpresp_(sys,type,wspec)
         % Computes sector index vs frequency
         [indx,w,FocusInfo,InfFlag] = ifpofpresp(sys.Data_,type,wspec,false);
      end
      
      function wc = getGainCrossover_(sys,g)
         % Gain crossover frequencies
         wc = getGainCrossover(sys.Data_,g);
      end
      
      function wc = getSectorCrossover_(sys,varargin)
         % Sector crossover frequencies
         wc = getSectorCrossover(sys.Data_,varargin{:});
      end
      
      function p = pole_(sys,varargin)
         % Poles of single system (SYS or SYS(J1,J2,...))
         D = localSelectSingleModel(sys.Data_,varargin);
         p = pole(D);
      end
      
      function [z,g] = zero_(sys,varargin)
         % Zeros of single SISO system (SYS or SYS(J1,J2,...))
         D = localSelectSingleModel(sys.Data_,varargin);
         [z,g] = zero(D);
      end
      
      function [z,nrk] = tzero_(sys,varargin)
         % Transmission zeros of single MIMO system
         D = sys.Data_;
         if issparse(D)
            error(message('Control:analysis:NoSparseDynamics',class(sys)))
         end
         [z,nrk] = tzero(ss(D),varargin{:});
      end
      
      function varargout = timeresp_(sys,RespType,tspec,~,Config)
         % Simulate step/impulse/initial response
         % [y,t,Focus,x,ysd] = timeresp_(sys,...)
         [varargout{1:nargout}] = timeresp(sys.Data_,RespType,tspec,Config);
      end
      
      function varargout = lsim_(sys,u,t,xinit,~,InterpRule)
         % Linear simulation
         no = nargout;
         varargout = cell(1,no);
         [varargout{1:min(no,2)}] = lsim(sys.Data_,u,t,xinit,InterpRule);
      end
      
      function [r,k] = rlocus_(sys,k)
         % Root locus
         [r,k] = rlocus(sys.Data_,k);
      end
      
      function [op,SINGULAR] = findop_(sys,~,~,opspec)
         % Compute operating condition for each model
         % Note: SYS is a StateSpaceModel subclass
         D = sys.Data_;
         [op,SINGULAR] = findop(D(1),opspec);
         nsys = numel(D);
         if nsys>1
            op = repmat(op,size(D));
            SINGULAR = repmat(SINGULAR,size(D));
            for ct=2:nsys
               [op(ct),SINGULAR(ct)] = findop(D(ct),opspec);
            end
         end
      end

      function sys = fixInput_(sys,indices,values)
         D = sys.Data_;
         for ct=1:numel(D)
            D(ct) = fixInput(D(ct),indices,values);
         end
         sys.Data_ = D;
      end

      %% TRANSFORMATIONS
      sys = chgTimeUnit_(sys,newUnits)
      [sys,gic] = c2d_(sys,Ts,options)
      [sys,gic] = d2c_(sys,options)
      sys = d2d_(sys,Ts,options)
      sys = upsample_(sys,L)
      [sys,icmap] = absorbDelay_(sys,SCOPE)
      sys = pade_(sys,Ni,No,Nf)
      [sys,varargout] = minreal_(sys,tol,dispflag)
      [sys,g,varargout] = balreal_(sys,Options)

      function [sys,blks,TL,TR] = modalreal_(sys,Options)
         try
            sys = ss(sys);
         catch
            error(message('Control:general:NotSupportedModelsofClass','modalreal',class(sys)))
         end
         [sys.Data_,blks,TL,TR] = modreal(sys.Data_,Options);
      end
      
      function [sys,T] = compreal_(sys,type)
         try
            sys = ss(sys);
         catch
            error(message('Control:general:NotSupportedModelsofClass','compreal',class(sys)))
         end
         [sys.Data_,T] = compreal(sys.Data_,type);
      end
      
      function [H,tau] = getDelayModel_(H)
         try
            [H.Data_,tau] = getInternalDelayModel(H.Data_);
         catch
            error(message('Control:ltiobject:setDelayModel4'))
         end
      end
      
      function H = setDelayModel_(H,tau)
         D = H.Data_;
         try
            for ct=1:numel(D)
               D(ct) = setInternalDelayModel(D(ct),tau);
            end
         catch
            error(message('Control:ltiobject:setDelayModel4'))
         end
         H.Data_ = D;
         H.IOSize_ = H.IOSize_ - numel(tau);
      end

      function sys = interp_(sys,w)
         % FRD interpolation. W is assumed in rad/TimeUnits
         Data = sys.Data_;
         for ct=1:numel(Data)
            Data(ct) = interp(Data(ct),w);
         end
         sys.Data_ = Data;
      end
      
      %% DESIGN
      function [C,Info,Merit] = pidtune_(G,C,Options,fixBC)
         % PID tuning
         ArraySize = getArraySize(G);
         nsys = prod(ArraySize);
         Info = struct('Stable',false(ArraySize),'CrossoverFrequency',0,'PhaseMargin',0);
         Merit = zeros(ArraySize);
         GData = G.Data_;  % expects ltipack.ltidata
         Dc = repmat(C.Data_(1),ArraySize);
         NeedsInteg = false;
         for ct=1:nsys
            % Create tuning data structure
            try
               DS = ltipack.PIDTuningData(GData(ct),C,Options);
            catch ME
               error(message('Control:ltiobject:GetPIDTuningData',class(G)))
            end
            % Tune PID
            [Dc(ct),Tinfo] = tune(DS,Options,true,fixBC);
            Info(ct).Stable = Tinfo.Stable;
            Info(ct).CrossoverFrequency = Tinfo.wc;
            Info(ct).PhaseMargin = Tinfo.PM;
            Merit(ct) = Tinfo.F;
            NeedsInteg = NeedsInteg || Tinfo.NeedsIntegrator;
         end
         C.Data_ = Dc;
         if NeedsInteg
            warning(message('Control:design:pidtune11'))
         end
      end
      
      function TuningData = getPIDTuningData_(G,C,Options)
         TuningData = ltipack.PIDTuningData(G.Data_,C,Options);
      end
      
      function isStable = checkNyquistStability_(OL,LoopSign,nUnstable)
         % Reduce OL.Data_ from ltipack.ltidata to zpkdata or frddata
         Data = getPIDPlantData(OL.Data_);
         isStable = checkNyquistStability(Data,LoopSign,nUnstable);
      end

      
      %% STATE-SPACE MODELS
      sys = ss2ss_(sys,T)
      sys = dss2ss_(sys,varargin)
      sys = ssequiv_(sys,TL,eL,TR,eR)
      sys = xperm_(sys,perm)
      W = gram_(sys,type,options)
      [sys,varargout] = canon_(sys,Type,varargin)
      [p,q] = covar_(sys,w,rw)
      
      function [sys,nx] = augstate_(sys)
         D = sys.Data_;
         nx = order(D(1));
         for ct=1:numel(D)
            if order(D(ct))~=nx
               error(message('Control:general:RequiresUniformNumberOfStates','augstate'))
            end
            D(ct) = augstate(D(ct));
         end
         sys.Data_ = D;
      end

      function sys = augoffset_(sys)
         D = sys.Data_;
         for ct=1:numel(D)
            D(ct) = augoffset(D(ct));
         end
         sys.Data_ = D;
      end

      function sys = xelim_(sys,elim,method)
         % State elimination or reduction
         sys.Data_ = xelim(sys.Data_,elim,method);
      end
      
      % RCTB
      function S = diskmargin_(L,e,JOB)
         % Gateway to ssdata/diskmargin and frddata/diskmargin
         D = L.Data_;
         [nL,~] = iosize(L); % number loop channels
         S = struct('LowerBound',[],'UpperBound',[],'Frequency',[],...
             'WorstPerturbation',[]);

         % Gather disk margin data
         if strcmp(JOB,'siso') && nL>1
            % loop-at-a-time case
            S = repmat(S,[nL size(D)]);
            for ct=1:numel(D)
               S(:,ct) = diskmargin(D(ct),e,JOB);
            end
         else
            S = repmat(S,size(D));
            for ct=1:numel(D)
               S(ct) = diskmargin(D(ct),e,JOB);
            end
         end
      end
      
      function [SM,INFO,IRD] = wcdiskmargin_(L,e,JOB,opt)
         % Gateway to lftdatSS/diskmargin and lftdatSSFRD/diskmargin
         D = L.Data_;  % assumed non-empty
         [nL,~] = iosize(D(1)); % number loop channels
         if strcmp(JOB,'siso') && nL>1
            % loop-at-a-time case
            for ct=1:numel(D)
               [X,Y,Z] = wcdiskmargin(D(ct),e,JOB,opt);  % nL-by-1 structs
               for j=nL:-1:1  % min disk margin over models
                  if ct==1 || X(j).LowerBound<SM(j).LowerBound
                     SM(j,1) = X(j);  INFO(j,1) = Y(j);  IRD(j,1) = Z(j);
                     INFO(j).Model = ct;
                  end
               end
            end
         else
            for ct=1:numel(D)
               [X,Y,Z] = wcdiskmargin(D(ct),e,JOB,opt);
               if ct==1 || X.LowerBound<SM.LowerBound
                  SM = X;  INFO = Y;  IRD = Z;  INFO.Model = ct;
               end
            end
         end
      end
            
   end
         
   %% STATIC METHODS
   methods (Static, Hidden)
      
      function Sinfo = getStateInfo(Data,Prop)
         % Helper for get.StateName, get.StatePath, and get.StateUnit.
         % Returns string vector of state names/units when the state 
         % dimension is uniform. If the state names or units clash, 
         % returns vector of "" strings.
         if isempty(Data)
            Sinfo = strings(0,1);
         else
            Sinfo = getStateInfo(Data(1),Prop);
            nx = numel(Sinfo);
            for ct = 2:numel(Data)
               S = getStateInfo(Data(ct),Prop);
               if numel(S)~=nx
                  switch Prop
                     case 'StateName'
                        error(message('Control:ltiobject:get2'))
                     case 'StatePath'
                        error(message('Control:ltiobject:get8'))
                     case 'StateUnit'
                        error(message('Control:ltiobject:get6'))
                  end
               end
               [Sinfo,clash] = ltipack.mrgname(Sinfo,S);
               if clash
                  Sinfo = strings(nx,1);  break;
               end
            end
         end
      end
      
   end
   
   
end

%----------------------------------

function Data = localSelectSingleModel(Data,Indices)
% Single model selection in model array
if ~isempty(Indices)
   try
      Data = Data(Indices{:});
   catch %#ok<CTCH>
      error(message('Control:ltiobject:access2'))
   end
end
if ~isscalar(Data)
   error(message('Control:ltiobject:access1'))
end
end
