dlinmod.m

数字通信第四版原书的例程

function [a,b,c,d]=dlinmod(fun,st,x,u,para,xpert,upert) 
%DLINMOD Obtains linear models from systems of ODEs and discrete-time systems. 
%	[A,B,C,D]=DLINMOD('SFUNC',TS) obtains the state-space linear model of
%	the system of mixed continuous and discrete system described in the 
%	S-function 'SFUNC' when  the state variables and inputs are set to zero
%	at the sample time TS.  SFUNC may be a SIMULINK or other model.
%
%	[A,B,C,D]=DLINMOD('SFUNC',TS,X,U) allows the state vector, X, and 
%	input, U, to be specified. A linear model will then be obtained
%	at this operating point. 
%
%	[A,B,C,D]=DLINMOD('SFUNC',TS,X,U,PARA) allows a vector of parameters
%	to be set.  PARA(1) sets the perturbation level for obtaining the
%	linear model (default PARA(1)=1e-5).  For systems that are functions
%	of time PARA(2) may be set with the value of t at which the linear 
%	model is to be obtained (default PARA(2)=0).
%
%	To see more help, enter TYPE DLINMOD.
%	See also LINMOD, LINMOD2, SFUNC, DSFUNC, TRIM, and the file DLINMOD.M.

%
%	[A,B,C,D]=DLINMOD('SFUNC',TS,X,U,PARA,XPERT,UPERT) allows the
%	perturbation levels for all of the elements of X and U to be set. 
%	The default is otherwise  XPERT=PARA(1)+1e-3*PARA(1)*abs(X),
%	UPERT=PARA(1)+1e-3*PARA(1)*abs(U).
%
%	Any or all of  PARA, XPERT, UPERT may be empty matrices in which case
%	these parameters will be assumed to be undefined and the default
%	option will be used.

%	Copyright (c) 1990-94 by The MathWorks, Inc.
%	Andrew Grace 11-12-90.

% ---------------Options--------------------
eval(['sizes =',fun, ';']);
eval(['[sizes,x0,str,ts, tsx] =',fun,'([],[],[],0);']);
if ~isempty(tsx), tsx = tsx(:,1); end
sizes=[sizes(:); zeros(6-length(sizes),1)];
nxz=sizes(1)+sizes(2); nu=sizes(4); ny=sizes(3);
nx=sizes(1);

if nargin<2, st = max(ts(:,1)); end
if nargin<3, x=[]; end
if nargin<4, u=[]; end
if nargin<5, para=[]; end
if nargin<6, xpert=[]; end
if nargin<7, upert=[]; end

if isempty(u), u=zeros(nu,1); end
if isempty(x), x=zeros(nxz,1); end
if isempty(para) , para=[0;0]; end
if para(1)==0, para(1)=1e-5; end
if isempty(upert), upert=para(1)+1e-3*para(1)*abs(u); end
if isempty(xpert), xpert=para(1)+1e-3*para(1)*abs(x); end
if length(para)>1, t=para(2); else t=0; end



ts = [0 0; ts];

% Eliminate sample times that are the same 
% with different offsets.

tsnew = [];
[nts, dum] = size(ts);
for i=1:nts
	clash = 0;
	for j = 1:length(tsnew)
		if ts(i,1) == tsnew(j)
			clash = 1;
		end
	end
	if  ~clash 
		tsnew = [tsnew; ts(i,1)];
	end
end
		

[nts] = length(tsnew);
a = zeros(nxz,nxz);
b = zeros(nxz, nu);
Acd = a;
Bcd = b;
Aeye = eye(nxz,nxz);

for m = 1:nts
	[list,ind]  = sort(tsnew);
% Start of with next smallest sampling period
	if length(list) > 1
		stnext = min(st, list(2));
	else
		stnext = st;
	end
	storig = list(1);
	index = find(tsx == storig);
	nindex = find(tsx ~= storig);
	oldA = Acd;
	oldB = Bcd;
	[Acd,Bcd,c,d]=linall(fun,storig,x,u,para,xpert,upert);
	a = a + Aeye * (Acd - oldA);
	b = b + Aeye * (Bcd - oldB);
	n = length(index);
	if n & storig ~= stnext
		if storig ~=  0
			if stnext ~= 0
				[ad2,bd2] = d2d(a(index, index),eye(n,n),storig, stnext);
			else 
				[ad2,bd2] = d2ci(a(index, index),eye(n,n),storig);
			end
		else 
			[ad2,bd2] = c2d(a(index, index),eye(n,n),stnext);
		end	
               	a(index, index)  =  ad2;

		if length(nindex)
               	 	a(index, nindex) = bd2*a(index,nindex);
		end
		if nu
               		b(index,:) = bd2*b(index,:);
		end
		Aeye(index,index) = bd2*Aeye(index,index);
		tsx(index) =  stnext(ones(length(index),1));
	end
	tsnew(ind(1)) = [];
end
if norm(imag(a), 'inf') < sqrt(eps), a = real(a); end
if norm(imag(b), 'inf') < sqrt(eps), b = real(b); end

% Return transfer function model
if nargout == 2
	disp('Returning transfer function model')
	% Eval it in case its not on the path
	[a, b] = feval('ss2tf',a,b,c,d,1);
end