剑桥系统辨识工具箱(system identification toolbox)

function [model,extra] = bilinid(data,n,i,Dzero,algorithm,lsmethod) %BILINID Identify a bilinear model from input-output data. % % [M,EXTRA] = bilinid(DATA) % [M,EXTRA] = bilinid(DATA,n) % [M,EXTRA] = bilinid(DATA,n,k) % [M,EXTRA] = bilinid(DATA,n,k,DMAT) % [M,EXTRA] = bilinid(DATA,n,k,DMAT,ALG) % [M,EXTRA] = bilinid(DATA,n,k,DMAT,ALG,LS) % % DATA is the input-output data given as an IDDATA object. M is the % estimated model returned as a discrete-time BILIN state-space system. % % n is the system order (optional): % (1) If n = [] then the algorithm will automatically select the order % such that n <= 10 (default). % (2) If n is a scalar, then n is the system order. % (3) If entered as a vector (e.g. [1 2 3 4 5]), a plot of singular % values will be given and the user will be prompted to select an order. % % k is the block size given as a positive integer (optional). If not % specified, k is chosen to be as large as possible while still trying to % be compatible with the size of DATA and choice of ALG. It is recommended % that k >= max(n). % % DMAT determines whether the D matrix is estimated or set to zero (optional). % 'Estimate': Estimate D (default). % 'Zero' : Set D = 0. % % ALG determines the choice of system identification algorithm and can be % one of {'general' | 'fast' | 'accurate'} (optional). ALG='fast' is only % allowed if the number of outputs < min(n), and ALG='accurate' is only % allowed if the number of outputs >= max(n). If ALG is not specified, % then the most appropriate algorithm is automatically chosen based on the % number of outputs, i.e. if number of outputs < max(n), then % ALG='general'. If number of outputs >= max(n), then ALG='accurate'. % % LS determines the choice of least squares method (optional) used in % estimating the system matrices and can be one of 'cls' (default) or % 'ols'. LS='cls' will result in using a constrained least squares method % and EXTRA.ls=0. LS='ols' will result in the use of an ordinary least % squares method. If LS='ols', then EXTRA.ls is a matrix containing % information about the accuracy of the estimation and should be close to % 0 for a good estimate. % % EXTRA is a structure containing additional information about the model % and data. EXTRA.SV is a vector containing the singular values resulting % from the SVD decomposition. EXTRA.Q, EXTRA.R and EXTRA.S are the % matrices that form the noise covariance matrix EXTRA.COV = [Q S; S' R]. % % See also BILIN, IDDATA, BILIN/SIM, BILIN/COMPARE, BILIN/PE. % % CUED System Identification Toolbox. % Cambridge University Engineering Department. % Copyright (C) 1998-2002. All Rights Reserved. % Version 1.00, Date: 01/06/2002 % Created by H. Chen and E.C. Kerrigan. % Check the arguments if nargin < 1 error('BILINID needs at least one input argument.') end % Turn the data into row vectors % Y=[y(1) y(2) ... y(N)]; U=[u(1) u(2) ... u(N)] if ~isa(data,'iddata') error('DATA should be an IDDATA object.') end Y = get(data,'OutputData'); U = get(data,'InputData'); Ts = get(data,'Ts'); %clear data Y = Y'; [p,cy] = size(Y); if p < 1 error('Need a non-empty output vector in DATA.') end U = U'; m = size(U,1); if m < 1 error('Need a non-empty input vector in DATA.') end if nargin > 1 if ~isempty(n) n = sort(n(:)); if any(n < 1 | round(n) ~= n) error('Only positive integers are allowed for the choice of system order.') end % if max(n) > i*p % warning('The block size might be too small to support this choice of order.') % end end else n = []; end if nargin > 3 if ~isempty(Dzero) switch lower(Dzero) case {'e','estimate','estimated'} isDzero=0; case {'z','zero','zeros'} isDzero=1; otherwise error('Invalid choice for DMAT. Select ''Estimate'' or ''Zero''.') end clear Dzero; else isDzero=0; end else isDzero=0; end if nargin > 4 if isempty(algorithm) if isempty(n) if p >= 10 % Maximum n will be 10 algorithm = 42; else algorithm = 3; end else if p < max(n) algorithm = 3; else algorithm = 42; end end else algorithm = lower(algorithm); switch algorithm case {'general','g'} algorithm = 3; case {'fast','f'} if ~isempty(n) if p < min(n) algorithm = 41; else error('ALG = ''fast'' is not allowed if the number of outputs is not less than min(n).') end else error('ALG = ''fast'' is not allowed if n is empty.') end case {'accurate','a'} if ~isempty(n) if p >= max(n) algorithm = 42; else error('ALG = ''accurate'' is not allowed if the number of outputs is less than max(n).') end else error('ALG = ''accurate'' is not allowed if n is empty.') end otherwise error('Invalid choice for ALG - choose ''general'', ''fast'' or ''accurate''.') end % case {'3','iii'} % algorithm = 3; % case {'4.1','iv.i','iv.1'} % algorithm = 41; % case {'4.2','iv.ii','iv.2'} % algorithm = 42; % otherwise % error('Invalid choice for ALG - choose ''III'', ''IV.I'' or ''IV.II''.') % end end else if isempty(n) if p >= 10 % maximum n is 10 algorithm = 42; else algorithm = 3; end else if p < max(n) algorithm = 3; else algorithm = 42; end end end switch algorithm case 3 disp(' ') disp('Using general four-block, deterministic-stochastic algorithm.') case 41 disp(' ') disp('Number of outputs < system order.') disp('Using fast four-block, deterministic-stochastic algorithm.') case 42 disp(' ') disp('Number of outputs >= system order.') disp('Using accurate four-block, deterministic-stochastic algorithm.') end if nargin > 5 if isempty(lsmethod) lsmethod = 'cls'; else lsmethod = lower(lsmethod); switch lsmethod case {'ols','o','ordinaryleastsquares'} lsmethod = 'ols'; case {'cls','c','constrainedleastsquares'} lsmethod = 'cls'; otherwise error('Invalid choice for LSMETHOD - choose ''ols'' or ''cls''.') end end else lsmethod = 'cls'; end switch lsmethod case 'ols' disp('Ordinary least squares will be used when estimating system matrices.') case 'cls' disp('Constrained least squares will be used when estimating system matrices.') end disp(' ') if nargin > 2 if isempty(i) i=1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end if algorithm ~= 3 di=0; for j=0:i-1 di=di+p*(m+1)^(j); end end switch algorithm case 3 while cy >= ei^(2)+ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1) i=i+1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end end case 41 while cy >= ei^(2)+ei+di i=i+1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end di=0; for j=0:i-1 di=di+p*(m+1)^(j); end end case 42 while cy >= di+3*(ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1)) i=i+1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end di=0; for j=0:i-1 di=di+p*(m+1)^(j); end end end i=i-1; disp(sprintf('Block size k = %d.',i)) else if i < 1 | round(i) ~= i error('Block size k should be a positive integer.') end ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end if algorithm ~= 3 di=0; for j=0:i-1 di=di+p*(m+1)^(j); end end switch algorithm case 3 if (cy < ei^(2)+ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1)) error('Not enough data points for this choice of block size.') end case 41 if (cy < ei^(2)+ei+di ) error('Not enough data points for this choice of block size.') end case 42 if (cy < di+3*(ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1))) error('Not enough data points for this choice of block size.') end end end else i=1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end if algorithm ~= 3 di=0; for j=0:i-1 di=di+p*(m+1)^(j); end end switch algorithm case 3 while cy >= ei^(2)+ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1) i=i+1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end end case 41 while cy >= ei^(2)+ei+di i=i+1; ei=0; for j=0:i-1 ei=ei+m*(m+1)^(j); end di=0