home *** CD-ROM | disk | FTP | other *** search
- function [a,b,c,d]=drmodel(n,p,m)
- %DRMODEL Generates random stable discrete nth order test models.
- %
- % [NUM,DEN]=DRMODEL(N) forms an Nth order SISO transfer function
- % model.
- %
- % [NUM,DEN]=DRMODEL(N,P) generates an Nth order single input, P
- % output transfer function model.
- %
- % [A,B,C,D]=DRMODEL(N) generates an Nth order SISO state space model
- %
- % [A,B,C,D]=DRMODEL(N,P,M) generates an Nth order, P output, M
- % input, state space model.
- %
- % See also: RMODEL.
-
- % Clay M. Thompson 12-27-90
- % Copyright (c) 1986-93 by the MathWorks, Inc.
-
- error(nargchk(0,3,nargin));
-
- if nargin==0,
- n=max([1,round(abs(10*randn(1,1)))]);
- if nargout==4,
- p=max([1,round(4*randn(1,1))]);
- m=max([1,round(4*randn(1,1))]);
- else
- m=1;
- p=1;
- end
- end
- if nargin<2, p=1; end
- if nargin<3, m=1; end
-
- rand('uniform')
- % Prob of an integrator is 0.10 for the first and 0.01 for all others
- nint = (rand(1,1)<0.10)+sum(rand(n-1,1)<0.01);
- % Prob of repeated roots is 0.05
- nrepeated = floor(sum(rand(n-nint,1)<0.05)/2);
- % Prob of complex roots is 0.5
- ncomplex = floor(sum(rand(n-nint-2*nrepeated,1)<0.5)/2);
- nreal = n-nint-2*nrepeated-2*ncomplex;
-
- % Determine random poles
- rep = 2*rand(nrepeated,1)-1;
- mag = rand(ncomplex,1);
- ang = pi*rand(ncomplex,1);
- jay = sqrt(-1);
- complex = mag.*exp(jay*ang);
- re = real(complex); im = imag(complex);
-
- if nargout==2, % Random transfer function
- poles = [re+jay*im;re-jay*im;ones(nint,1);rep;rep;2*rand(nreal,1)-1];
-
- % Determine random zeros
- zer = inf*ones(n,p);
- jay = sqrt(-1);
- for i=1:p
- rand('uniform')
- % Prob of complex zeros is 0.35
- ncomplex = floor(sum(rand(n,1)<0.35)/2);
- % Prob of real zero is 0.35
- nreal = sum(rand(n-2*ncomplex,1)<0.35);
- nzeros = 2*ncomplex+nreal;
- if nzeros>0,
- re = randn(ncomplex,1); im = randn(ncomplex,1);
- zer(1:nzeros,i) = [re+jay*im;re-jay*im;randn(nreal,1)];
- end
- end
- [a,b] = zp2tf(zer,poles,randn(p,1));
-
- else % Random state space model
- a = zeros(n,n);
- for i=1:ncomplex,
- ndx = [2*i-1,2*i];
- a(ndx,ndx) = [re(i),im(i);-im(i),re(i)];
- end
- ndx = [2*ncomplex+1:n];
- if ~isempty(ndx),
- a(ndx,ndx) = diag([ones(nint,1);rep;rep;2*rand(nreal,1)-1]);
- end
- T = orth(rand(n,n));
- a = T\a*T;
- b = rand(n,m);
- c = rand(p,n);
- d = rand(p,m);
- rand('uniform')
- b = b .* (rand(n,m)<0.75);
- c = c .* (rand(p,n)<0.75);
- d = d .* (rand(p,m)<0.5);
- end
-