3431
|
1 ## Copyright (C) 1996 Auburn University. All rights reserved. |
|
2 ## |
|
3 ## This file is part of Octave. |
|
4 ## |
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
6 ## under the terms of the GNU General Public License as published by the |
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
8 ## later version. |
|
9 ## |
|
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT |
|
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
13 ## for more details. |
|
14 ## |
|
15 ## You should have received a copy of the GNU General Public License |
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
18 |
|
19 ## -*- texinfo -*- |
|
20 ## @deftypefn {Function File} {} analdemo () |
|
21 ## Octave Controls toolbox demo: State Space analysis demo |
|
22 ## @end deftypefn |
|
23 |
|
24 ## Author: David Clem |
|
25 ## Created: August 15, 1994 |
|
26 ## Updated by John Ingram December 1996 |
|
27 |
|
28 function analdemo () |
|
29 |
|
30 while (1) |
|
31 clc |
|
32 k=0; |
|
33 while(k > 8 || k < 1) |
|
34 k = menu("Octave State Space Analysis Demo", ... |
|
35 "System grammians (gram, dgram)", ... |
|
36 "System zeros (tzero)", ... |
|
37 "Continuous => Discrete and Discrete => Continuous conversions (c2d,d2c)", ... |
|
38 "Algebraic Riccati Equation (are, dare)", ... |
|
39 "Balanced realizations (balreal, dbalreal)", ... |
|
40 "Open loop truncations via Hankel singular values (balreal, dbalreal)", ... |
|
41 "SISO pole placement", ... |
|
42 "Return to main demo menu"); |
|
43 endwhile |
|
44 if (k == 1) |
|
45 clc |
|
46 help dgram |
|
47 prompt |
|
48 |
|
49 clc |
|
50 disp("System Grammians: (see Moore, IEEE T-AC, 1981) \n"); |
|
51 disp("Example #1, consider the discrete time state space system:\n"); |
|
52 a=[1, 5, -8.4; 1.2, -3, 5; 1, 7, 9] |
|
53 b=[1, 5; 2, 6; -4.4, 5] |
|
54 c=[1, -1.5, 2; 6, -9.8, 1] |
|
55 d=0 |
|
56 prompt |
|
57 disp("\nThe discrete controllability grammian is computed as follows:"); |
|
58 cmd = "grammian = dgram(a, b);"; |
|
59 run_cmd; |
|
60 disp("Results:\n"); |
|
61 grammian = dgram(a,b) |
|
62 disp("Variable Description:\n"); |
|
63 disp("grammian => discrete controllability grammian"); |
|
64 disp("a, b => a and b matrices of discrete time system\n"); |
|
65 disp("A dual approach may be used to compute the observability grammian."); |
|
66 prompt |
|
67 clc |
|
68 |
|
69 help gram |
|
70 prompt |
|
71 clc |
|
72 |
|
73 disp("Example #2, consider the continuous state space system:\n"); |
4139
|
74 a=[1, 3, -10.2; 3.7, -2, 9; 1, 3, 7] |
3431
|
75 b=[1, 12; 6, 2; -3.8, 7] |
|
76 c=[1, -1.1, 7; 3, -9.8, 2] |
|
77 d=0 |
|
78 prompt |
|
79 disp("\nThe continuous controllability grammian is computed as follows:"); |
|
80 cmd = "grammian = gram(a, b);"; |
|
81 run_cmd; |
|
82 disp("Results:\n"); |
|
83 grammian = gram(a,b) |
|
84 disp("Variable Description:\n"); |
|
85 disp("grammian => continuous controllability grammian"); |
|
86 disp("a, b => a and b matrices of continuous time system\n"); |
|
87 disp("A dual approach may be used to compute the observability grammian."); |
|
88 prompt |
|
89 clc |
|
90 |
|
91 |
|
92 elseif (k == 2) |
|
93 clc |
|
94 help tzero |
|
95 prompt |
|
96 |
|
97 disp("System zeros (tzero) example\n"); |
|
98 disp("Example #1, consider the state space system:\n"); |
|
99 a=[0, 1, 0; -10, -2, 0; -10, 0, -8] |
|
100 b=[0; 1; 9] |
|
101 c=[-10, 0, -4] |
|
102 d=1 |
|
103 prompt |
|
104 disp("\nTo compute the zeros of this system, enter the following command:\n"); |
|
105 cmd = "zer = tzero(a,b,c,d);"; |
|
106 run_cmd; |
|
107 disp("Results:\n"); |
|
108 zer = tzero(a,b,c,d) |
|
109 disp("Variable Description:\n"); |
|
110 disp("zer => zeros of state space system"); |
|
111 disp("a, b, c, d => state space system used as input argument"); |
|
112 prompt |
|
113 clc |
|
114 |
|
115 disp("Example #2, consider the state space system from example 1 again:"); |
4771
|
116 cmd = "sys = ss(a,b,c,d);"; |
3431
|
117 disp(cmd); |
|
118 eval(cmd); |
|
119 sysout(sys); |
|
120 disp("\nThe zeros of this system can also be calculated directly from the"); |
|
121 disp("system variable:"); |
|
122 cmd = "zer = tzero(sys);"; |
|
123 run_cmd; |
|
124 disp("Results:\n") |
|
125 zer = tzero(sys) |
|
126 disp("Variable Description:\n"); |
|
127 disp("zer => zeros of state space system"); |
|
128 disp("sys => state space system used as input argument"); |
|
129 prompt |
|
130 clc |
|
131 |
|
132 elseif (k == 3) |
|
133 clc |
|
134 help c2d |
|
135 prompt |
|
136 |
|
137 clc |
|
138 disp("Continuous => Discrete and Discrete => Continuous conversions (c2d,d2c)"); |
|
139 disp("\nExample #1, consider the following continuous state space system"); |
4771
|
140 cmd = "sys_cont = ss([-11, 6; -15, 8], [1; 2], [2, -1], 0);"; |
3431
|
141 eval(cmd); |
|
142 disp(cmd); |
|
143 disp("Examine the poles and zeros of the continuous system:"); |
|
144 sysout(sys_cont,"all"); |
|
145 disp("\nTo convert this to a discrete system, a sampling time is needed:"); |
|
146 cmd = "Tsam = 0.5;"; |
|
147 run_cmd; |
|
148 disp("\nNow convert to a discrete system with the command:"); |
|
149 cmd = "sys_disc = c2d(sys_cont,Tsam);"; |
|
150 run_cmd; |
|
151 disp("Examine the poles and zeros of the discrete system:"); |
|
152 sysout(sys_disc,"all"); |
|
153 prompt |
|
154 clc |
|
155 |
|
156 disp("\nNow we will convert the discrete system back to a continuous system"); |
|
157 disp("using the d2c command:"); |
|
158 help d2c |
|
159 prompt |
|
160 cmd = "new_sys_cont = d2c(sys_disc);"; |
|
161 run_cmd; |
|
162 disp("\nExamine the poles and zeros of the discrete system:"); |
|
163 sysout(new_sys_cont,"all"); |
|
164 prompt |
|
165 |
|
166 elseif (k == 4) |
|
167 clc |
|
168 help are |
|
169 prompt |
|
170 clc |
|
171 |
|
172 disp("Algebraic Riccati Equation (are, dare)"); |
|
173 |
|
174 disp("\nExample #1, consider the continuous state space system:\n"); |
|
175 a=[1, 3, -10.2; 3.7, -2, 9; 1, 3, 7] |
|
176 b=[1, 12; 6, 2; -3.8, 7] |
|
177 c=[1, -1.1, 7; 3, -9.8, 2] |
|
178 d=0 |
|
179 prompt |
|
180 disp("\nThe solution to the continuous algebraic riccati equation"); |
|
181 disp("is computed as follows:"); |
|
182 cmd = "x_cont = are(a, b, c);"; |
|
183 run_cmd; |
|
184 disp("Results:\n") |
|
185 x_cont = are(a,b,c) |
|
186 disp("Variable Description:\n") |
|
187 disp("x_cont => solution to the continuous algebraic riccati equation"); |
|
188 disp("a, b, c => a, b, and c matrices of continuous time system\n"); |
|
189 prompt |
|
190 |
|
191 clc |
|
192 help dare |
|
193 prompt |
|
194 clc |
|
195 |
|
196 disp("Example #2, consider the discrete time state space system:\n"); |
|
197 a=[1, 5, -8.4; 1.2, -3, 5; 1, 7, 9] |
|
198 b=[1, 5; 2, 6; -4.4, 5] |
|
199 c=[1, -1.5, 2; 6, -9.8, 1] |
|
200 d=0 |
|
201 r=eye(columns(b)) |
|
202 prompt |
|
203 disp("\nThe solution to the continuous algebraic riccati equation"); |
|
204 disp("is computed as follows:"); |
|
205 cmd = "x_disc = dare(a, b, c, r);"; |
|
206 run_cmd; |
|
207 disp("Results:\n") |
|
208 x_disc = dare(a,b,c,r) |
|
209 disp("Variable Description:\n"); |
|
210 disp("x_disc => solution to the discrete algebraic riccati equation"); |
|
211 disp("a, b, c => a, b and c matrices of discrete time system\n"); |
|
212 prompt |
|
213 clc |
|
214 |
|
215 elseif (k == 5) |
|
216 disp("--- Balanced realization: not yet implemented") |
|
217 elseif (k == 6) |
|
218 disp("--- Open loop balanced truncation: not yet implemented") |
|
219 elseif (k == 7) |
|
220 disp("SISO pole placement example:") |
4771
|
221 cmd = "sys=tf(1, [1, -2, 1]);"; |
3431
|
222 run_cmd |
|
223 disp("System in zero-pole form is:") |
|
224 cmd = "sysout(sys,\"zp\");"; |
|
225 run_cmd |
|
226 disp("and in state space form:") |
|
227 cmd = "sysout(sys,\"ss\");"; |
|
228 run_cmd |
|
229 disp("Desired poles at -1, -1"); |
|
230 cmd = "K=place(sys, [-1, -1])"; |
|
231 run_cmd |
|
232 disp("Check results:") |
|
233 cmd = "[A,B] = sys2ss(sys);"; |
|
234 run_cmd |
|
235 cmd = "poles=eig(A-B*K)"; |
|
236 run_cmd |
|
237 prompt |
|
238 elseif (k == 8) |
|
239 return |
|
240 endif |
|
241 endwhile |
|
242 endfunction |
|
243 |