home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / RMODEL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.5 KB  |  92 lines
  1. function [a,b,c,d]=rmodel(n,p,m)
  2. %RMODEL Generates random stable continuous nth order test models.
  3. %
  4. %    [NUM,DEN]=RMODEL(N) generates an Nth order SISO transfer function
  5. %         model.
  6. %    [NUM,DEN]=RMODEL(N,P) generates an Nth order single input, P 
  7. %        output transfer function model.
  8. %    [A,B,C,D]=RMODEL(N) generates an Nth order SISO state space model.
  9. %    [A,B,C,D]=RMODEL(N,P,M) generates an Nth order, P output, M input,
  10. %        state space model.
  11. %
  12. %    See also: DRMODEL.
  13.  
  14. %    Clay M. Thompson 1-22-91
  15. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  16.  
  17. if nargin==0,
  18.   n=max([1,round(abs(10*randn(1,1)))]);
  19.   if nargout==4,
  20.     p=max([1,round(4*randn(1,1))]);
  21.     m=max([1,round(4*randn(1,1))]);
  22.   else
  23.     m=1;
  24.     p=1;
  25.   end
  26. end
  27. if nargin<2, p=1; end
  28. if nargin<3, m=1; end
  29.  
  30. rand('uniform')
  31. % Prob of an integrator is 0.10 for the first and 0.01 for all others
  32.  nint = (rand(1,1)<0.10)+sum(rand(n-1,1)<0.01);
  33. % Prob of repeated roots is 0.05
  34.  nrepeated = floor(sum(rand(n-nint,1)<0.05)/2);
  35. % Prob of complex roots is 0.5
  36.  ncomplex = floor(sum(rand(n-nint-2*nrepeated,1)<0.5)/2);
  37.  nreal = n-nint-2*nrepeated-2*ncomplex;
  38.  
  39. if nargout==2, nzeros = sum(rand(n,p)<0.7); end
  40.  
  41. % Determine random poles
  42. rep = -exp(randn(nrepeated,1));
  43. re = -exp(randn(ncomplex,1));
  44. % Make imaginary part bigger for more oscillatory (and interesting) models
  45. im = 3 * exp(randn(ncomplex,1));
  46. jay = sqrt(-1);
  47. complex = -exp(randn(ncomplex,1))+sqrt(-1)*exp(randn(ncomplex,1));
  48.  
  49. if nargout==2, % Random transfer function
  50.   poles = [re+jay*im;re-jay*im;zeros(nint,1);rep;rep;-exp(randn(nreal,1))];
  51.  
  52.   % Determine random zeros
  53.   zer = inf*ones(n,p);
  54.   jay = sqrt(-1);
  55.   for i=1:p
  56.     rand('uniform')
  57.     % Prob of complex zeros is 0.35
  58.     ncomplex = floor(sum(rand(n,1)<0.35)/2);
  59.     % Prob of real zero is 0.35
  60.     nreal = sum(rand(n-2*ncomplex,1)<0.35);
  61.     nzeros = 2*ncomplex+nreal;
  62.     if nzeros>0,
  63.       re = randn(ncomplex,1); im = randn(ncomplex,1);
  64.       zer(1:nzeros,i) = [re+jay*im;re-jay*im;randn(nreal,1)];
  65.     end
  66.   end
  67.   [a,b] = zp2tf(zer,poles,randn(p,1));
  68.  
  69. else % Random state space model
  70.   a = zeros(n,n);
  71.   for i=1:ncomplex,
  72.     ndx = [2*i-1,2*i];
  73.     a(ndx,ndx) = [re(i),im(i);-im(i),re(i)];
  74.   end
  75.   ndx = [2*ncomplex+1:n];
  76.   if ~isempty(ndx),
  77.     a(ndx,ndx) = diag([zeros(nint,1);rep;rep;-exp(randn(nreal,1))]);
  78.   end
  79.   T = orth(randn(n,n));
  80.   a = T\a*T;
  81.   b = randn(n,m);
  82.   c = randn(p,n);
  83.   d = randn(p,m);
  84.   rand('uniform')
  85.   b = b .* (rand(n,m)<0.75);
  86.   c = c .* (rand(p,n)<0.75);
  87.   d = d .* (rand(p,m)<0.5);
  88. end
  89.  
  90.   
  91.  
  92.