home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / covar < prev    next >
Encoding:
Text File  |  1995-11-15  |  3.7 KB  |  139 lines
  1. //-----------------------------------------------------------------------
  2. //
  3. // covar
  4. //
  5. // Syntax: R=covar(A,B,C,D,W)  or  R=covar(NUM,DEN,W)
  6. //
  7. // This routine computes the covariance response of a continuous
  8. // system to white noise input. Calling the routine as,
  9. //
  10. //    R=covar(A,B,C,D,W)
  11. //
  12. // computes the covariance response of the continuous state-space
  13. // system defined by A,B,C,D to Gaussian white noise input with
  14. // intensity W given by the following,
  15. //  
  16. //    E[w(t)w(tau)'] = W delta(t-tau) 
  17. //  
  18. // where delta(t) is the dirac delta function and E is the expectation
  19. // operator.
  20. //
  21. // The routine returns both the state covariance matrix X and the
  22. // output covariance matrix Y given by the following,
  23. //
  24. //     X=E[xx']
  25. //     Y=E[yy']
  26. //
  27. // where x is the state vector and y is the output vector.
  28. //
  29. // This routine can also be used to compute the covariance response
  30. // of a transfer function type system by calling it as,
  31. //
  32. //     R=covar(NUM,DEN,W)
  33. //
  34. // where NUM is a vector containing the coefficients of the numerator
  35. // transfer function polynomial, DEN contains the coefficients of the
  36. // denominator transfer function polynomial, and W is the intensity
  37. // of the white noise input.
  38. //
  39. //  ==================================================================
  40. //  Warning: Unstable systems or systems with a non-zero D matrix have
  41. //  an infinite covariance response.
  42. //  ==================================================================
  43. //
  44. // The results are returned in a list. For example,  R=covar(A,B,C,D,W)
  45. // produces,
  46. //
  47. //  R.Y = output covariance matrix (Square of Output RMS)
  48. //  R.X = state covariance matrix
  49. //
  50. // Ref: Skelton, R. "Dynamic Systems Control Linear System Analysis and
  51. //      Synthesis," John Wiley and Sons, 1988.
  52. //
  53. // Copyright (C), by Jeffrey B. Layton, 1994
  54. // Version JBL 940917
  55. //-----------------------------------------------------------------------
  56.  
  57. rfile tfchk
  58. rfile abcdchk
  59. rfile tf2ss
  60. rfile lyap
  61. rfile isempty
  62.  
  63. covar = function(a,b,c,d,w)
  64. {
  65.    local(narg,D,W,Dum,A,B,C,D,num,den,msg,estr,X,Y,pinf)
  66.  
  67. // count number of input arguments
  68.    narg=0;
  69.    if (exist(a)) { narg=narg+1; }
  70.    if (exist(b)) { narg=narg+1; }
  71.    if (exist(c)) { narg=narg+1; }
  72.    if (exist(d)) { narg=narg+1; }
  73.    if (exist(w)) { narg=narg+1; }
  74.  
  75.    if ( narg != 3) {
  76.         if ( narg != 5) {
  77.             error("covar: Wrong number of input arguments.");
  78.         }
  79.    }
  80.  
  81.    if (narg == 3) {   // T.F. syntax
  82.        Dum=tfchk(a,b);
  83.        num=Dum.numc;
  84.        den=Dum.denc;
  85.        W=c;
  86. // Convert to state-space form
  87.        Dum=tf2ss(num,den,2);
  88.        A=Dum.a;
  89.        B=Dum.b;
  90.        C=Dum.c;
  91.        D=Dum.d;
  92.    else
  93.        A=a;
  94.        B=b;
  95.        C=c;
  96.        D=d;
  97.        W=w;
  98.    }
  99.  
  100.    if (narg == 5) {
  101.        msg="";
  102.        msg=abcdchk(A,B,C,D);
  103.        if (msg != "") {
  104.            estr="COVAR: "+msg;
  105.            error(estr);
  106.        }
  107.    }
  108.  
  109.    X=lyap(A,B*W*B');
  110.    Y=C*X*C';
  111.  
  112. // Systems with non-zero D matrix have infinite output covariance.
  113.    if ( any( any( abs(d) > epsilon()) ) ) {
  114.    printf("%s","Warning: Systems with non-zero D matrix have infinite output covariance.\n");
  115.        pinf=inf()*ones(Y.nr,Y.nc);
  116.        if (!isempty(Y)) {
  117.            Y=inf()*ones(Y.nr,Y.nc);
  118.        }
  119.    }
  120.  
  121. // A valid covariance must be positive semi-definite.
  122.    if (min(real(eig(X).val)) < -epsilon()) {
  123.    printf("%s","Warning: Invalid covaraince - not positive semi-definite. Returning infinity.\n");
  124.        X=inf()*ones(X.nr,X.nc);
  125.        Y=inf()*ones(Y.nr,Y.nc);
  126.        return << X=X; Y=Y >>
  127.    }
  128.  
  129. // The system must be stable for a valid covariance
  130.    if (any(real(eig(a).val) > 0) ) {
  131.        printf("%s","Warning: Unstable system. Returning Infinity.\n");
  132.        X=inf()*ones(length(X));
  133.        Y=inf()*ones(length(Y));
  134.        return << X=X; Y=Y >>
  135.    }
  136.  
  137.    return << X=X; Y=Y >>
  138. };
  139.