path: root/modules/cacsd/macros/linfn.sci

// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA -
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt

function [x,frequ]=linfn(G,PREC,RELTOL,options);
    //[x,frequ]=linfn(G,PREC,RELTOL,options);
    //  Computes the Linf (or Hinf) norm of a transfer function
    //                 -1
    //      G(s) = D + C (sI - A)  B
    //
    //  This norm is well-defined as soon as the realization
    //  (A,B,C) has no imaginary eval which is both controllable and observable.
    //
    //  The algorithm follows the paper by G. Robel (AC-34 pp. 882-884, 1989).
    //  The case D=0 is not treated separately due to superior accuracy of
    //  the general method when (A,B,C) is nearly nonminimal.
    //
    //  In the general case (A neither stable nor antistable), no upper bound is
    //  prespecified. If by contrast A is stable or antistable, lower
    //  and upper bounds are computed using the associated Lyapunov
    //  solutions (see Glover).


    //  On input:
    //  ---------
    //     *  G is a syslin list
    //     *  PREC is the desired relative accuracy on the norm
    //     *  RELTOL: relative threshold to decide when an eigenvalue can be
    //          considered on the imaginary axis.
    //     *  available options are
    //          - 'trace': traces each bisection step, i.e., displays the lower
    //               and upper bounds and the current test point.
    //          - 'cond':  estimates a confidence index on the computed value
    //               and issues a warning if computations are
    //               ill-conditioned
    //
    //  On output:
    //  ---------
    //     *  x is the computed norm.
    //     *  freq:  list of the frequencies for which ||G|| is attained, i.e.,
    //          such that ||G (j om)|| = ||G||. If -1 is in the list, the norm
    //          is attained at infinity. If -2 is in the list, G is all-pass in
    //          some direction so that ||G (j omega)|| = ||G|| for all
    //          frequencies omega.
    //!
    //
    //  Called macros:
    //  -------------
    //    heval_test, cond_test, list_set
    //
    //  History:
    //  -------
    //      author: P. Gahinet, INRIA
    //      last modification: Oct 3nd, 1991
    //****************************************************************************


    //constants
    //*********
    INIT_LOW=1.0e-4; INIT_UPP=1.0e5;  INFTY=10e10;
    frequ=[];


    //user interface. The default values are:
    //     PREC=1.0e-3; RELTOL=1.0e-10; options='nul';
    //************************************************
    [lhs,rhs]=argn(0);
    select rhs,
    case 0 then
        error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"linfn",1))
    case 1 then
        PREC=1.0e-3; RELTOL=1.0e-10; options="nul";
    case 2 then
        RELTOL=1.0e-10;
        if type(PREC)==10 then
            iopt=2
            options=PREC; PREC=1.0e-3;
        else
            options="nul";
        end,
    case 3 then
        if type(RELTOL)==10 then
            iopt=3
            options=RELTOL; RELTOL=1.0e-10;
        else
            options="nul";
        end,
    end

    if typeof(G)<>"state-space" then
        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space  expected.\n"),"linfn",1))
    end
    if G.dt<>"c" then
        error(msprintf(gettext("%s: Wrong type for argument #%d: In continuous time expected.\n"),"linfn",1))
    end

    if type(options)<>10|and(options<>["t","nul"]) then
        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"linfn",iopt,"""t"",""nul"""))
    end
    if type(PREC)<>1|size(PREC,"*")<>1|~isreal(PREC)|PREC<=0 then
        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be a positive scalar.\n"),"linfn",2))
    end
    if type(RELTOL)<>1|size(RELTOL,"*")<>1|~isreal(RELTOL)|RELTOL<=0 then
        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be a positive scalar.\n"),"linfn",3))
    end


    //recover realization
    //*******************
    [a,b,c,d]=abcd(G);


    //SCALING
    //*******                          ||B||.||C||
    //  Scale A,B,C so that ||A||=||B||=||C||=1. With scale:= -----------,
    //                                 ||A||
    //  and DD:=D/scale, AA:=A/||A||, BB:=B/||B||, CC:=C/||C||, we have
    //                           -1
    //       || G || = scale * || DD + CC (sI - AA)   BB ||
    //
    //  From now on, the parameters A,B,C,D are scaled

    noa=norm(a,"inf"); nob=norm(b,"inf"); noc=norm(c,"inf"); nobc=nob*noc;

    if nobc==0, x=norm(d); return, end

    scale=nobc/noa; a=a/noa; b=b/nob; c=c/noc; d=d/scale; nd=norm(d);


    //test the spectrum of A
    //**********************

    s=real(spec(a));

    if min(abs(s)) < RELTOL then
        mprintf(gettext("%s: WARNING: the A matrix has eigenvalues near the imaginary axis.\n"),"linfn");
    end


    //Search window initialization
    //****************************
    // Initialize the search window [lower,upper] where `lower' and `upper'
    // are lower and upper bounds on the Linf norm of G. When no such
    // bounds are available, the window is arbitrarily set to [INIT_LOW,INIT_UPP]
    // and the variables LOW and UPP keep record of these initial values so that
    // the window can be extended if necessary.

    if max(s)*min(s) > 0 then
        // A is stable or antistable: use associated Lyapunov equations
        // to derive lower and upper bounds.
        p=lyap(a',-b*b',"c");
        q=lyap(a,-c'*c,"c");
        s=sqrt(abs(spec(p*q)));

        lower=max(nd,max(s)); LOW=0;
        upper=nd+2*sum(diag(s));   UPP=100*upper;

    else
        if nd==0 then
            lower=INIT_LOW; LOW=INIT_LOW;
        else
            lower=nd;       LOW=0;
        end
        upper=INIT_UPP;     UPP=INIT_UPP;
    end



    //form the constant parts of the pencil (E,F) (see G. Robel).
    //***********************************************************
    [na,na]=size(a); twona=2*na;
    [p,m]=size(d);
    nf=twona+min(m,p); //size of e and f

    // to ensure that D'*D is of size min(m,p), replace (a,b,c,d) by
    // (a',c',b',d') if m>p
    if m>p then
        a=a'; d=d'; aux=b; b=c'; c=aux';
    end

    e=eye(2*na,2*na); e(nf,nf)=0;
    nul=0; nul(na,na)=0;
    f=[a,nul;-c'*c,-a']; f(nf,nf)=0;
    dd=d'*d; Cd=c'*d;


    //----------------------
    //     BISECTION STARTS
    //----------------------

    while %t,

        ga=sqrt(lower*upper);   //test point gamma = log middle of [lower,upper]

        if part(options,1)=="t" then
            write(%io(2),[scale*lower,scale*ga,scale*upper],..
            "(''lower,current,upper = '',3e20.10)");
        end

        bga=b/ga; Cdga=Cd/ga;
        f(1:na,twona+1:nf)=-bga;
        f(na+1:twona,twona+1:nf)=Cdga;
        f(twona+1:nf,1:nf)=[Cdga',bga',eye(dd)-dd/(ga**2)];


        // Test for generalized eigenvalues on the imaginary axis
        // ------------------------------------------------------

        [dist,frequ]=heval_test(e,f,RELTOL,"test");

        if dist < RELTOL then
            lower=ga;  LOW=0;
            //  eigenvalue on the imaginary axis: gamma < ||G||
        else
            upper=ga;  UPP=100*upper;
            //  gamma > ||G||
        end


        // Search window management:
        //--------------------------
        //   If the gamma-iteration runs into one of the initial arbitrary bounds
        //   LOW or UPP, extend the search window to allow for continuation

        if ga<10*LOW then
            lower=LOW/10; LOW=lower;
        end
        // expand search window toward gamma<<1
        if ga>UPP/10 then
            upper=UPP*10; UPP=upper;
        end
        // expand search window toward gamma>>1


        // Termination tests
        //------------------

        if lower > INFTY then
            mprintf(gettext("%s: Controllable & observable mode(s) of A near the imaginary axis"),"linfn");
            x=scale*lower;
            return;
        else
            if upper < 1.0e-10 then
                x=scale*upper;
                mprintf(gettext("%s: All modes of A are nearly nonminimal so that || G || is almost 0.\n"),"linfn");
                return;
            else
                if 1-lower/upper < PREC,
                    ga=sqrt(lower*upper);
                    x=scale*ga;

                    // Compute all the frequencies achieving ||G||
                    if lower<>0 then
                        bga=b/lower; Cdga=Cd/lower;
                        f(1:na,twona+1:nf)=-bga;
                        f(na+1:twona,twona+1:nf)=Cdga;
                        f(twona+1:nf,1:nf)=[Cdga',bga',eye(dd)-dd/(lower**2)];
                        [dist,frequ]=heval_test(e,f,RELTOL,"freq");
                    end
                    if frequ==[] then
                        mprintf(gettext("%s: The computed value of || G || may be inaccurate.\n"),"linfn");
                    end

                    // evaluate the condition of the eigenproblem of (e,f) near || G ||
                    if part(options,1)=="c" then
                        gt=1.1*ga;
                        f=[a,nul,-b/gt;cc,at,Cd/gt;dc/gt,bt/gt,eye(dd)-dd/(gt**2)]
                        co=cond_test(e,f,frequ,RELTOL);
                        if co < RELTOL then
                            mprintf(gettext("%s: The computed value of || G || may be inaccurate.\n"),"linfn");
                        end
                    end
                    //-----------

                    return;
                end,
            end,
        end

    end//end while


endfunction


function [dist,frequ]=heval_test(e,f,TOL,option);
    //[dist,frequ]=heval_test(e,f,TOL,option);
    // This procedure estimates the distance of the generalized spectrum
    // of the pencil  f - lambda e  to the imaginary axis. Here e is always
    // of the form  diag(I_(nf-nz),0_nz). The distance is 0 whenever there are
    // (nearly) infinite eigenvalues or eigenvalues of the form 0/0.
    //
    // The eigenvalues are computed via a generalized Schur decomposition
    // of  f - lambda e . Let (a(i),b(i)) : i=1..nf be the output of gspec.
    // Three cases must be distinguished:
    //    * both a(i) and b(i) are << 1 -> singularity of the pencil
    //    * b(i)<<1 and a(i) close to 1 -> infinite eigenvalue
    //    * both a(i) and b(i) are close to 1 -> finite eigenvalue.
    //
    // Let nz denote the rank deficiency of e which is also the size of D'*D.
    // For gamma > ||G||, the generalized spectrum of (e,f) consists of exactly
    // nz infinite modes and nf-nz finite ones. By contrast, there may be
    // additional singularities or infinite modes for ||D|| <= gamma <= ||G||,
    // depending on whether ||D|| = ||G|| or ||D|| < ||G||.
    //
    // If ||D|| < ||G||, there are still exactly nz infinite modes for gamma in
    // [ ||D|| , ||G|| ]. At gamma=||G||, some finite mode(s) hit the imaginary
    // axis and their imaginary part omega is such that ||G(j omega)|| = ||G||.
    //
    // If ||D|| == ||G|| now, we always have ||G (infinity)|| = ||G||
    // and if moreover some pair (a(i),b(i)) is nearly (0,0), then
    // || G (j omega) || = || G || for all omega's (direction along which G is
    // all-pass). Note that G is all-pass iff there are nz pairs
    // (a(i),b(i)) nearly equal to (0,0). Finally, finite modes which hit
    // the imaginary axis still yield frequencies for which ||G|| is attained.
    //
    // Two options are available in this function:
    // * When option='test', the function counts the number of finite modes.
    // If less than nf-nz (nf=order of f), it concludes gamma <= ||G|| and returns
    // dist=0. Otherwise, it estimates the distance of the finite spectrum
    //                  min | Re(l_i) |
    // to the imag. axis computed as  --------------- where the l_i's denote
    //                    max | l_i |
    // the pencil finite eigenvalues.
    // * With the option 'freq' (used for gamma = ||G||), it furthermore returns
    // all frequencies for which ||G|| is attained. Infinite frequencies are
    // denoted by -1 and if ||G(j omega)|| == ||G|| for all omega's, frequ=[-2];
    //
    //
    // Input:
    //    * (e,f): pencil
    //    * TOL: relative tolerance on the size of eigenvalue real parts.
    //    * option: 'test' or 'freq'.
    //
    // Output
    //    * dist: distance of the spectrum to the imaginary axis as defined above.
    //    * frequ: list of frequencies for which ||G|| is attained.
    //
    //!
    //balancing
    frequ=[]; evals=[];
    [f,xx]=balanc(f);
    [nf,nf]=size(f);
    nz=nf-sum(diag(e)); //rank deficiency of e


    //Generalized Schur decomposition of the pencil (f,e)
    //---------------------------------------------------
    [a,b]=spec(f,e);


    if option=="test" then
        //***********************************
        //Simple test and computation of dist
        //***********************************

        // Check that there are exactly nz infinite modes of (e,f) and compute dist

        nai=0; //nai: number of infinite or (0,0) modes (b(i) << 1)
        for i=1:nf,
            bi=abs(b(i));
            if bi < 100*TOL then
                nai=nai+1;
            else
                evals=[evals,a(i)/bi];
            end
        end

        if nai>nz then
            dist=0;
        else
            dist=min(abs(real(evals)))/max(abs(evals));
        end

    else
        //option = 'freq'
        //*************************************************
        //Compute the frequency for which ||G|| is attained
        //*************************************************

        // Here gamma is appx equal to ||G||. Distinguish two cases:
        // ||D|| < ||G|| and ||D|| = ||G||.

        if min(svd(f(nf-nz+1:nf,nf-nz+1:nf))) < TOL then
            //-----------------------------------------------
            // f(nf-nz+1:nf,nf-nz+1:nf)= I - (D'*D)/||G||**2 -> case ||D||=||G||

            noa=max(abs(a));  na=0;     //na -> # pairs (0,0)
            frequ=[-1];  //||G|| is always attained for s=infinity in this case
            for i=1:nf,
                bi=b(i);
                if abs(bi) < 100*TOL then
                    if abs(a(i)) < 100*TOL*noa, na=na+1; end
                else
                    evals=[evals,a(i)/bi];
                end
            end

            if na>0 then    // G is all-pass along some direction
                frequ=[-2];
                if na>=nz, mprintf(_("G is all-pass")); end
            else
                if evals<>[] then
                    maxabs=max(abs(evals));
                    for i=1:max(size(evals)),
                        if abs(real(evals(i))) <= TOL*maxabs,
                            frequ=[frequ,abs(imag(evals(i)))];
                        end
                    end,
                end,
            end


        else
            // case ||D|| < ||G||
            //------------------------

            for i=1:nf,
                bi=b(i);
                if abs(bi) > 100*TOL, evals=[evals,a(i)/bi]; end
            end

            maxabs=max(abs(evals));
            for i=1:max(size(evals)),
                if abs(real(evals(i))) <= TOL*maxabs then
                    frequ=[frequ,abs(imag(evals(i)))];
                end
            end

        end //endif

        frequ=list_set(frequ,1.0e-5); //eliminate redundancy in frequ

    end //endif

endfunction


function [c]=cond_test(e,f,frequ,TOL);
    //[c]=cond_test(e,f,frequ,TOL);
    //   This procedure returns a confidence index for the computed gamma = ||G||
    //   at which some generalized eigenvalue(s) of (e,f) meets the imaginary
    //   axis. Specifically, it considers gamma := 1.1*computed value of ||G||
    //   and computes how close (e,f) is to have imaginary eigenvalues
    //   for this gamma. If very close, this indicates (e,f) has generalized
    //   eigenvalues near the imaginary axis for all gamma's in an interval
    //   around ||G|| whence the computed value is likely to be inaccurate.
    //!
    [nf,nf]=size(f);
    c=1;

    for i=1:max(size(frequ)),
        s=svd(f-%i*frequ(i)*e);
        c=min(c,s(nf)/s(1));
        if c < TOL then return; end
    end

endfunction

function [l]=list_set(l,TOL);
    //[l]=list_set(l,TOL);
    //  eliminates redundant elements in a list. Two entries are considered
    //  identical when their difference is smaller then TOL (in relative terms)
    //!
    nl=max(size(l));
    i=1;

    while i < nl,
        entry=l(i);  TOLabs=TOL*entry;
        j=i+1;
        while j <= nl,
            if abs(l(j)-entry) <= TOLabs then
                l=[l(1:j-1),l(j+1:nl)];
                nl=nl-1;
            else
                j=j+1;
            end
        end
        i=i+1;
    end
endfunction