0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % DOBODE makes a bode computation from the given inputs to outputs.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % DESCRIPTION: makes a bode computation from all the inputs to all outputs.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % CALL: varargout = doBode(a, b, c, d, w, Ts)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % INPUTS:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % 'sys' - ssm object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 % 'inputnames, statenames, outputnames' - 3 cellstr
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % OUTPUTS:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 % 'as' - array of output TFs containing the requested responses.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 function varargout = doBode(a, b, c, d, w, Ts)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 f = w/(2*pi);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 %% dealing with timestep
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 if(Ts ~= 0)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 % Compute discrete freq:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 f_disc = f*Ts;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % Compute z = e^jw (with w discrete)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 z = exp(1i*2*pi*f_disc);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 %% getting matrices properly ordered
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 %% stop warnings
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 s = warning;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 warning('off','MATLAB:nearlySingularMatrix'); % turn all warnings off
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 %% looping over SISO systems
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 Ni = size(b,2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 No = size(c,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 G = zeros(No,Ni,length(f));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 G2 = zeros(No,Ni,length(f));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 %% compute responses
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 % To conserve accuracy, the system is transformed to Hessenberg form [1]
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 % Laub, A.J., "Efficient Multivariable Frequency Response Computations",
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 % IEEE Transactions on Automatic Control, AC-26 (1981), pp. 407-408
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 % This transformation leads to a faster resolution and at the
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 % same time to a more accurate and precise answer of C*(jw - A)^-1 *B + D
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 reduce_a = a;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 reduce_c = c;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 %% loop on inputs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 if (Ts == 0)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 for ii=1:Ni
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 reduce_b = b(:,ii);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 reduce_d = d(:,ii);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 [num,den] = ss2tf(reduce_a,reduce_b,reduce_c,reduce_d);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 denr = polyval(den, 2*pi*1i*f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 for jj = 1:No
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 numr = polyval(num(jj,:), 2*pi*1i*f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 G(jj,ii,:) = numr./denr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 [T,H]=hess(reduce_a);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 % Step 1:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 P = reduce_c*T;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 reduce_b = b(:,1:Ni);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 reduce_d = d(1:No,1:Ni);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 Q = T'*reduce_b;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 I = eye(size(reduce_a));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 %% time discrete
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 for ff = 1:length(f)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 G(1:No,1:Ni,ff) = (P/(z(ff)*I - H))*(Q) + reduce_d;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 %% execute warnings
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 [msg, msgid] = lastwarn;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 if(strcmp(msgid,'MATLAB:nearlySingularMatrix') == 1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 %warning('This frequency response may be unaccurate because the system matrix is close to singular or badly scaled.') %#ok<WNTAG>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 warning(s) % restore the warning state
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 %% output
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 varargout = {G};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90
|