Overview
Comment:Decided that integrating the IR is a bad idea.
So now explicitly computes the SR using matrix exponential -- much
more accurate.
Downloads: Tarball | ZIP archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 5fd859876619d06b4377a089c1e40ce0aaff92be229b6ae066d13dcf64dd790b
User & Date: gawthrop@users.sourceforge.net on 1996-12-05 10:38:36.000
Other Links: branch diff | manifest | tags
Context
1996-12-05
10:41:39
Saved many changes over past 2 weeks. check-in: fe65307e73 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
10:38:36
Decided that integrating the IR is a bad idea.
So now explicitly computes the SR using matrix exponential -- much
more accurate.
check-in: 5fd8598766 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
10:18:52
Saved (quite old) version which integrates IR to give SR check-in: e67a236326 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes
︙︙
12
13
14
15
16
17
18



19
20
21
22
23
24
25
% u0 input gain vector: u = u0*unit step.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$



% %% Revision 1.3  1996/10/27  10:39:04  peterg
% %% Only compute matrix exponential once.
% %%
% %% Revision 1.2  1996/09/10  16:48:21  peter
% %% Changed ar counts in default settings.
% %%
% %% Revision 1.1  1996/08/19 15:34:29  peter







>
>
>







12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
% u0 input gain vector: u = u0*unit step.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.4  1996/12/05 10:18:52  peterg
% %% Saved (quite old) version which integrates IR to give SR
% %%
% %% Revision 1.3  1996/10/27  10:39:04  peterg
% %% Only compute matrix exponential once.
% %%
% %% Revision 1.2  1996/09/10  16:48:21  peter
% %% Changed ar counts in default settings.
% %%
% %% Revision 1.1  1996/08/19 15:34:29  peter
︙︙
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88





89













90
91
92
if M>N
  T = T';
  N = M;
end;



% $$$ one = eye(Nx);
% $$$ 
% $$$ Y = zeros(N,Ny);
% $$$ X = zeros(N,Nx);
% $$$ 
% $$$ dt = T(2)-T(1);% Assumes fixed interval
% $$$ expAdt = expm(A*dt); % Compute matrix exponential
% $$$ i = 0;
% $$$ expAt = one;
% $$$ 
% $$$ for t = T'
% $$$   i=i+1;
% $$$   if Nx>0
% $$$     x = ( A\(expAt-one) )*B*u0 + expAt*x0;
% $$$     expAt = expAt+expAdt;
% $$$     X(i,:) = x';
% $$$     if Ny>0
% $$$       y = C*x + D*u0;
% $$$       Y(i,:) = y';
% $$$     end;
% $$$   elseif Ny>0
% $$$     y = D*u0;
% $$$     Y(i,:) = y';
% $$$   end;
% $$$ end;

% Compute the impulse response
[Y,X] = sm2ir(A,B,C,D,T,u0,x0);

% Assume fixed sample interval
dT = T(2)-T(1);

% Do an Euler integration on it
 Y = mtt_euler(Y,dT);
 
 if nargout>1
   X = mtt_euler(X,dT);
 end;





 























|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
<
<
<
<
|
<
<
<
<
<
<
<
<
<
<
<
<
<
<
|
>
>
>
>
>
|
>
>
>
>
>
>
>
>
>
>
>
>
>



47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71




72














73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
if M>N
  T = T';
  N = M;
end;



one = eye(Nx);

Y = zeros(N,Ny);
X = zeros(N,Nx);

dt = T(2)-T(1);% Assumes fixed interval
expAdt = expm(A*dt); % Compute matrix exponential
i = 0;
expAt = one;

for t = T'
  i=i+1;
  if Nx>0
    x = ( A\(expAt-one) )*B*u0 + expAt*x0;
    expAt = expAt*expAdt;
    X(i,:) = x';
    if Ny>0
      y = C*x + D*u0;




      Y(i,:) = y';














    end;
  elseif Ny>0
    y = D*u0;
    Y(i,:) = y';
  end;
end;

% $$$ % Compute the impulse response
% $$$ [Y,X] = sm2ir(A,B,C,D,T,u0,x0);
% $$$ 
% $$$ % Assume fixed sample interval
% $$$ dT = T(2)-T(1);
% $$$ 
% $$$ % Do an Euler integration on it
% $$$  Y = mtt_euler(Y,dT);
% $$$  
% $$$  if nargout>1
% $$$    X = mtt_euler(X,dT);
% $$$  end;
% $$$  




MTT: Model Transformation Tools
[ GitHub | SourceHut | Sourceforge | Fossil RSS ]