classdef (CaseInsensitiveProperties = true, TruncatedProperties = true) ...
        Poles < TuningGoal.BandLimited & TuningGoal.GenericLoop & ...
        TuningGoal.SystemLevel
    % Constraint on closed-loop dynamics.
    %
    %   R = TuningGoal.Poles(MINDECAY,MINDAMPING,MAXFREQ) creates a tuning
    %   requirement R for constraining the closed-loop system dynamics. The
    %   scalars MINDECAY, MINDAMPING, and MAXFREQ specify the minimum decay
    %   rate, minimum damping, and maximum natural frequency for the
    %   closed-loop poles. Set MINDECAY=0, MINDAMPING=0, or MAXFREQ=Inf to
    %   skip any of these three constraints. You can use this requirement to
    %   increase damping and prevent slow or fast dynamics.
    %
    %   Use the "Openings" property of R to open some feedback loops and
    %   constrain the dynamics of the corresponding open-loop configuration.
    %   For example, for a control system with cascaded loops,
    %      R = TuningGoal.Poles(0,0.3,Inf)
    %      R.Openings = 'OuterLoop'
    %   requires that the system remain stable with a minimum damping of 0.3
    %   when the outer loop is open. In MATLAB, use an AnalysisPoint block
    %   to mark the loop opening location 'OuterLoop'. In Simulink, use a
    %   Linear Analysis point or the addPoint method of the slTuner interface
    %   to mark this location.
    %
    %   R = TuningGoal.Poles(LOC,MINDECAY,MINDAMPING,MAXFREQ) constrains the
    %   poles of the sensitivity function measured at the location LOC (see
    %   TuningGoal.Sensitivity for details). Use this syntax to narrow the
    %   scope of the requirement to a particular feedback loop. For the
    %   cascaded loop example above,
    %      R = TuningGoal.Poles('InnerLoop',0,0.5,50))
    %      R.Openings = 'OuterLoop'
    %   constrains the inner-loop dynamics when the outer loop is open. The
    %   dynamics of blocks that do not participate to the inner loop are
    %   ignored. The location 'InnerLoop' can be anywhere along the inner
    %   loop.
    %
    %   Use VIEWSPEC(R) to visualize this requirement and use SYSTUNE and
    %   related commands to tune the control system parameters subject to
    %   this and other requirements.
    %
    %   Note: This requirement implicitly enforces closed-loop stability for
    %   the specified Openings configuration.
    %
    %   See also pole, evalGoal, viewGoal, AnalysisPoint, getPoints,
    %   TuningGoal.ControllerPoles, TuningGoal.Margins, systune, looptune,
    %   slTuner.

    % Author: P. Gahinet
    % Copyright 2009-2016 The MathWorks, Inc.

    % Note: For the default settings, this requirement contributes nothing
    % to f(x) or g(x) but still enforces stability of the closed-loop
    % configuration specified by Openings. In other words, this requirement
    % always contributes to the fStab term.

    properties
        % Minimum decay rate of closed-loop poles (default = 0).
        %
        % The closed-loop poles must satisfy
        %       Re(s) < -alpha       (continuous)
        %    log(|z|) < -alpha*Ts    (discrete)
        % where alpha is the specified minimum decay rate (a nonnegative value).
        % This constraint only applies to poles with natural frequencies in the
        % range specified by the "Focus" property.
        MinDecay = 0;

        % Minimum damping of closed-loop poles (default = 0).
        %
        % The closed-loop poles must satisfy
        %    Re(s) < -zeta * |s|
        % where zeta is the specified minimum damping ratio (a value between 0
        % and 1). This constraint only applies to poles with natural frequencies
        % in the range specified by the "Focus" property. In discrete time, the
        % damping ratio is computed using s=log(z)/Ts.
        MinDamping = 0;

        % Maximum natural frequency of closed-loop poles (default = Inf).
        %
        % The closed-loop poles must satisfy
        %         |s| < rho       (continuous)
        %    |log(z)| < rho*Ts    (discrete)
        % where rho is the specified maximum natural frequency (also known
        % as the spectral radius). This is useful to prevent fast closed-loop
        % dynamics. This constraint only applies to poles with natural
        % frequencies in the range specified by the "Focus" property.
        MaxFrequency = Inf;
    end

    methods

        function this = Poles(varargin)
            narginchk(0,4)
            ni = nargin;
            if ni>0 && isnumeric(varargin{1})
                % Set LOCATION=0x1 cell
                varargin = [{cell(0,1)},varargin];  ni = ni+1;
            end
            try
                if ni>0
                    this.Location = varargin{1};
                end
                if ni>1
                    this.MinDecay = varargin{2};
                end
                if ni>2
                    this.MinDamping = varargin{3};
                end
                if ni>3
                    this.MaxFrequency = varargin{4};
                end
            catch ME
                throw(ME)
            end
        end

        function this = set.MinDecay(this,Value)
            % SET function for MinDecay
            if ~(isnumeric(Value) && isscalar(Value) && ...
                    isreal(Value) && isfinite(Value) && Value>=0)
                error(message('Control:tuning:PoleReq1'))
            end
            this.MinDecay = double(Value);
        end

        function this = set.MinDamping(this,Value)
            % SET function for MinDamping
            if ~(isnumeric(Value) && isscalar(Value) && ...
                    isreal(Value) && isfinite(Value) && Value>=0 && Value<=1)
                error(message('Control:tuning:PoleReq5'))
            end
            this.MinDamping = double(Value);
        end

        function this = set.MaxFrequency(this,Value)
            % SET function for MaxFrequency
            if ~(isnumeric(Value) && isscalar(Value) && ...
                    isreal(Value) && Value>0)
                error(message('Control:tuning:PoleReq2'))
            end
            this.MaxFrequency = double(Value);
        end

    end

    methods (Access = protected)

        function H = getClosedLoopTransfer_(this,CL,varargin)
            % Computes closed-loop transfer from inputs to outputs.
            if isempty(this.Location)
                % Return model(s) with no I/Os and same dynamics as closed loop
                % Note: Same as getDynamics, but supports uncertainty sampling
                CL = genss(CL); % slTuner does not handle empty I/Os
                H = getIOTransfer(CL,cell(0,1),cell(0,1),this.Openings,this.Models);
                H = getValue(H,varargin{:});
            else
                H = getSensitivity(CL,this.Location,this.Openings,this.Models);
                H = sminreal(getValue(H,varargin{:}));
            end
        end


        function [H,fObj] = evalSpec_(this,CL)
            % Evaluates requirement for given closed-loop model
            % Note: Requirement is normalized so that in the fObj range [0.25,4],
            % a violation of x% (fObj=1+0.01*x) corresponds to a x% violation
            % of the spectral requirement.
            H = getClosedLoopTransfer_(this,CL);
            if nargout>1
                Focus = this.Focus;
                sH = size(H);
                fObj = NaN([sH(3:end) 1 1]);
                Ts = abs(H.Ts);
                MINDECAY = 1e-7;
                WarnFlag = false;
                FiniteFlag = isfinite(H,'elem');
                for ct=1:numel(fObj)
                    if FiniteFlag(ct)
                        E = pole(H,ct);
                        if Ts>0
                            E = log(E(E~=0))/Ts;
                        end
                        absE = abs(E);
                        WarnFlag = WarnFlag || any(absE>=Focus(1) & absE<MINDECAY);
                        inFocus = (absE>=max(Focus(1),MINDECAY) & absE<=Focus(2));
                        if any(inFocus)
                            absE = absE(inFocus);
                            fSpec = zeros(1,3);
                            if this.MinDecay>0
                                fSpec(1) = NSOptUtil.SpectralPenalty(...
                                    min(-real(E(inFocus)))/this.MinDecay);
                            end
                            if this.MinDamping>0
                                fSpec(2) = NSOptUtil.SpectralPenalty(...
                                    min(-real(E(inFocus))./absE)/this.MinDamping);
                            end
                            fSpec(3) = max(absE)/this.MaxFrequency;
                            fObj(ct) = max(fSpec);
                        else
                            fObj(ct) = 0;
                        end
                    end
                end
                if WarnFlag
                    % Warn that quasi-integrators are ignored when computing FOBJ
                    % (for consistency with the way SYSTUNE ignored fixed integrators)
                    warning(message('Control:tuning:PoleReq6'))
                end
            end
        end

        function delta = getViewExtent_(this)
            % Compute extent for view
            delta = 8;
            if this.MaxFrequency<Inf
                delta = max(delta,1.1*this.MaxFrequency);
            end
            alpha = this.MinDecay;
            if alpha>0
                delta = max(delta,alpha);
                if this.MinDamping>1e-2
                    delta = max(delta,2*alpha*sqrt(1/this.MinDamping^2-1));
                end
            end
        end

    end

    methods (Hidden) %% Tuning goal API

        % Compute data for Goal wave forms
        function [Ts,TU] = viewSpecGoalData(~,CL)
            if isequal(CL,[])
                Ts = 0;
                TU = 'seconds';
            else
                Ts = CL.Ts;
                TU = CL.TimeUnit;
            end
        end

        % Compute data for Design wave forms
        function Hp = viewSpecDesignData(this,CL)
            H = getClosedLoopTransfer_(this,CL,'usample');
            nsys = nmodels(H);  % array support
            Hp = zpk(ones(1,1,nsys),'Ts',H.Ts,'TimeUnit',H.TimeUnit);
            for ct=1:nsys
                p = pole(H,ct);
                wn = damp(p,H.Ts);
                p = p(wn>=this.Focus(1) & wn<=this.Focus(2));
                Hp(:,:,ct).p = p;
            end
            Hp.Name = getString(message('Control:systunegui:TGPlotPoles'));
        end

        % Validate if data is compatible for plotting
        function b = isCompatible(this,CL,CurrentDesignWaveforms)
            % Check for consistent size and sample time
            Hp = viewSpecDesignData(this,CL);
            NewValue = Hp;
            if controllib.chart.internal.utils.isChartResponse(CurrentDesignWaveforms.Hp)
                OldValue = CurrentDesignWaveforms.Hp.Model;
            else
                OldValue = CurrentDesignWaveforms.Hp.DataSrc.Model;
            end
            b = isequal(size(NewValue),size(OldValue)) && isequal(NewValue.Ts,OldValue.Ts);
        end

        % Update limits
        function updateLimits(~,Ts,plotHandle,XLimitsFocus,YLimitsFocus)
                % Pick square Y limits that will show constraints if nearly active
            if controllib.chart.internal.utils.isChart(plotHandle)
                view = qeGetView(plotHandle);
                axGrid = qeGetAxesGrid(view);
                Xlim = XLimitsFocus{1};
                Ylim = YLimitsFocus{1};
                if Ts==0
                    xc = (Xlim(1)+Xlim(2))/2;
                    delta = max(Xlim(2)-Xlim(1),Ylim(2)-Ylim(1));
                    Xlim = [xc-delta,xc+delta];
                    Ylim = [-delta,delta];
                else
                    xc = 0;
                    delta = max([2 Xlim(2)-Xlim(1) Ylim(2)-Ylim(1)]);
                    Xlim = [xc-0.6*delta,xc+0.6*delta];
                    Ylim = 0.6*[-delta,delta];
                end
                axGrid.XLimitsFocus = {Xlim};
                axGrid.YLimitsFocus = {Ylim};
            else
                ax = plotHandle;
                XLimMode = XLimitsFocus;
                YLimMode = YLimitsFocus;
                if strcmp(XLimMode,'auto') || strcmp(YLimMode,'auto')
                    Xlim = get(ax,'Xlim');
                    Ylim = get(ax,'Ylim');
                    if Ts==0
                        xc = (Xlim(1)+Xlim(2))/2;
                        delta = max(Xlim(2)-Xlim(1),Ylim(2)-Ylim(1));
                        set(ax,'XLim',[xc-delta,xc+delta],'YLim',[-delta,delta])
                    else
                        xc = 0;
                        delta = max([2 Xlim(2)-Xlim(1) Ylim(2)-Ylim(1)]);
                        set(ax,'XLim',[xc-0.6*delta,xc+0.6*delta],'YLim',0.6*[-delta,delta])
                    end
                end
            end
        end

        function [h, GoalWaveforms, CurrentDesignWaveforms,Ts] = createPlot(this,CL,ax)
            % Create plot
            [Ts,TU] = viewSpecGoalData(this,CL);
            GoalWaveforms = struct('Bounds',[]);
            if isequal(CL,[])
                % Just show bounds and target shape
                if controllibutils.CSTCustomSettings.getCSTPlotsVersion == 2
                    h = controllib.chart.PZPlot(Axes=ax,Options=pzoptions('cstprefs'));
                    delta = getViewExtent_(this);
                    xc = -0.75*delta;
                    h.XLimits = {[xc-delta,xc+delta]};
                    h.YLimits = {[-delta,delta]};
                else
                    h = ltiplot(ax,'pzmap',[],[],pzoptions('cstprefs'),cstprefs.tbxprefs);
                    ltiplotmenu(h,'pzmap');
                    delta = getViewExtent_(this);
                    xc = -0.75*delta;
                    set(getaxes(h.AxesGrid),'XLim',[xc-delta,xc+delta],'YLim',[-delta,delta])
                    h.Visible = 'on';  % new plot created with Visible='off'
                end
                CurrentDesignWaveforms = struct('Hp',[]);
            else
                Hp = viewSpecDesignData(this,CL);
                h = pzplot(ax,Hp,'b');
                if controllib.chart.internal.utils.isChart(h)
                    CurrentDesignWaveforms = struct('Hp',h.Systems(1));
                    h.Systems(1).Style.SemanticColor = controllib.plot.internal.utils.GraphicsColor(1).SemanticName;
                else
                    CurrentDesignWaveforms = struct('Hp',h.Response(1));
                    h.Response(1).setstyle(SemanticColor=controllib.plot.internal.utils.GraphicsColor(1).SemanticName);
                end
            end
            % Thicker lines
            if controllib.chart.internal.utils.isChart(h)
                for ct=1:numel(h.Systems)
                    h.Systems(ct).Style.LineWidth = 1.5;
                end
            else
                for ct=1:numel(h.Response)
                    h.Response(ct).setstyle('LineWidth',1.5)
                end
            end
            % Title
            if controllib.chart.internal.utils.isChart(h)
                h.Title.String = getString(message('Control:tuning:strPole1',this.Name));
                h.Title.Interpreter = 'none';
                h.Grid = 'on';
            else
                h.AxesGrid.Title = getString(message('Control:tuning:strPole1',this.Name));
                h.AxesGrid.TitleStyle.Interpreter = 'none';
                h.AxesGrid.Grid = 'on';
            end
            % Plot bounds
            if controllib.chart.internal.utils.isChart(h)
                tau = funitconv('rad/TimeUnit','rad/s',TU);
                h.addBoundSystem(MinDecay=tau*this.MinDecay,MinDamping=this.MinDamping,...
                    MaxFrequency=tau*this.MaxFrequency,Ts=Ts,...
                    SystemName=getString(message('Control:systunegui:TGPlotBounds')));
                h.Systems(end).Style.SemanticFaceColor = controllib.plot.internal.utils.GraphicsColor(8).SemanticName;
                h.Systems(end).Style.SemanticEdgeColor = controllib.plot.internal.utils.GraphicsColor(8).SemanticName;
                h.Systems(end).Style.FaceAlpha = 0.3;
                h.Systems(end).Style.EdgeAlpha = 0.3;
                GoalWaveforms.Bounds = h.Systems(end);
            else
                tau = funitconv('rad/TimeUnit','rad/s',TU);
                r = h.addPoleBound(tau*this.MinDecay,this.MinDamping,tau*this.MaxFrequency,Ts);
                r.Name = getString(message('Control:systunegui:TGPlotBounds'));
                r.setstyle(SemanticColor=controllib.plot.internal.utils.GraphicsColor(8,"tertiary").SemanticName);
                % Set Goal wave forms for TBound
                GoalWaveforms.Bounds = r;
            end
        end

        % Update Goal waveforms
        function [Ts,TU] = updateGoal(this,CL,hPlot,GoalWaveforms)
            % Compute data for goals
            [Ts,TU] = viewSpecGoalData(this,CL);
            h = getPlotHandle(hPlot);
            % Goal Focus
            tau = funitconv('rad/TimeUnit','rad/s',TU);
            % Update plotted goals
            if controllib.chart.internal.utils.isChart(h)
                GoalWaveforms.Bounds.MaxFrequency = tau*this.MaxFrequency;
                GoalWaveforms.Bounds.MinDecay = tau*this.MinDecay;
                GoalWaveforms.Bounds.MinDamping = this.MinDamping;
            else
                GoalWaveforms.Bounds.Data.MaxFrequency = tau*this.MaxFrequency;
                GoalWaveforms.Bounds.Data.MinDecay = tau*this.MinDecay;
                GoalWaveforms.Bounds.Data.MinDamping = this.MinDamping;
                GoalWaveforms.Bounds.draw
            end
            % Set Title for the plot
            setTitle(hPlot,getString(message('Control:tuning:strPole1',this.Name)));
        end

        % Update Design waveforms
        function updateDesign(this,DesignWaveforms,CL)
            Hp = viewSpecDesignData(this,CL);
            if controllib.chart.internal.utils.isChartResponse(DesignWaveforms.Hp)
                DesignWaveforms.Hp.Model = Hp;
            else
                DesignWaveforms.Hp.DataSrc.Model = Hp;
            end
        end

        % Add design
        function DesignWaveforms = addDesign(this,PlotHandle,Design,Name,Color)
            % Compute data for design
            Hp = viewSpecDesignData(this,Design);
            % Add response
            if controllib.chart.internal.utils.isChart(PlotHandle)
                addSystem(PlotHandle,Hp,SystemName=sprintf('%s: %s',getString(message('Control:systunegui:TGPlotPoles')),Name));
                PlotHandle.Systems(end).Style.SemanticColor = Color;
                PlotHandle.Systems(end).Style.LineWidth = 1.75;
                DesignWaveforms = struct('Hp',PlotHandle.Systems(end));
            else
                src = resppack.ltisource(Hp, 'Name', sprintf('%s: %s',getString(message('Control:systunegui:TGPlotPoles')),Name));
                rX = PlotHandle.addresponse(src);
                % Define characteristics
                chars = src.getCharacteristics('sigma');
                rX.setCharacteristics(chars);
                rX.DataFcn = {'pzmap' src rX};
                rX.setstyle('Color',Color,'LineWidth',1.75)
                rX.draw;
                DesignWaveforms = struct('Hp',rX);
            end
        end

        function Type = getComparisonStyleType(~)
            Type = 'Color';
        end

    end

    methods (Hidden)

        function validateGoal(this,CL)
            % Note: Needed to validate Models selection
            if this.MinDecay>=this.MaxFrequency
                error(message('Control:systunegui:PoleMinMax'))
            end
            if isempty(this.Location)
                getDynamics(CL,this.Openings,this.Models);
            else
                getSensitivity(CL,this.Location,this.Openings,this.Models);
            end
        end

        function [SPEC,LoopConfigs] = getSpecData(this,SPEC,LoopConfigs,uNames,yNames,sNames,~)
            % Builds standardized requirement description for SYSTUNE
            if this.MinDecay>=this.MaxFrequency
                error(message('Control:tuning:PoleReq3',getID(this)))
            end
            SPEC.Type = 0;
            SPEC.Band = this.Focus;
            SPEC.Spectral = struct(...
                'MinDecay',this.MinDecay,...
                'MinDamping',this.MinDamping,...
                'MaxFrequency',this.MaxFrequency);
            if ~isempty(this.Location)
                % Locate channels where sensitivity is measured
                [iLoop,MisMatch] = ltipack.resolveSignalID(this.Location,sNames,true);
                error(resolveSignalError(this,'Control:tuning:TuningReq8',MisMatch,sNames))
                SPEC.Input = numel(uNames) + iLoop;
                SPEC.Output = numel(yNames) + iLoop;
            end
            % Locate openings
            if isempty(this.Openings)
                iOpen = [];
            else
                [iOpen,MisMatch] = ltipack.resolveSignalID(this.Openings,sNames,true);
                error(resolveSignalError(this,'Control:tuning:TuningReq8',MisMatch,sNames))
            end
            [SPEC.Config,LoopConfigs] = ...
                TuningGoal.SystemLevel.getSwitchConfig(iOpen,LoopConfigs);
        end

    end

end
