1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
|
<HTML>
<HEAD><TITLE>AB09MD - SLICOT Library Routine Documentation</TITLE>
</HEAD>
<BODY>
<H2><A Name="AB09MD">AB09MD</A></H2>
<H3>
Balance & Truncate model reduction for the stable part of a system
</H3>
<A HREF ="#Specification"><B>[Specification]</B></A>
<A HREF ="#Arguments"><B>[Arguments]</B></A>
<A HREF ="#Method"><B>[Method]</B></A>
<A HREF ="#References"><B>[References]</B></A>
<A HREF ="#Comments"><B>[Comments]</B></A>
<A HREF ="#Example"><B>[Example]</B></A>
<P>
<B><FONT SIZE="+1">Purpose</FONT></B>
<PRE>
To compute a reduced order model (Ar,Br,Cr) for an original
state-space representation (A,B,C) by using either the square-root
or the balancing-free square-root Balance & Truncate (B & T)
model reduction method for the ALPHA-stable part of the system.
</PRE>
<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
<PRE>
SUBROUTINE AB09MD( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR, ALPHA,
$ A, LDA, B, LDB, C, LDC, NS, HSV, TOL, IWORK,
$ DWORK, LDWORK, IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, EQUIL, JOB, ORDSEL
INTEGER INFO, IWARN, LDA, LDB, LDC, LDWORK, M, N, NR,
$ NS, P
DOUBLE PRECISION ALPHA, TOL
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), HSV(*)
</PRE>
<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
<P>
<B>Mode Parameters</B>
<PRE>
DICO CHARACTER*1
Specifies the type of the original system as follows:
= 'C': continuous-time system;
= 'D': discrete-time system.
JOB CHARACTER*1
Specifies the model reduction approach to be used
as follows:
= 'B': use the square-root Balance & Truncate method;
= 'N': use the balancing-free square-root
Balance & Truncate method.
EQUIL CHARACTER*1
Specifies whether the user wishes to preliminarily
equilibrate the triplet (A,B,C) as follows:
= 'S': perform equilibration (scaling);
= 'N': do not perform equilibration.
ORDSEL CHARACTER*1
Specifies the order selection method as follows:
= 'F': the resulting order NR is fixed;
= 'A': the resulting order NR is automatically determined
on basis of the given tolerance TOL.
</PRE>
<B>Input/Output Parameters</B>
<PRE>
N (input) INTEGER
The order of the original state-space representation, i.e.
the order of the matrix A. N >= 0.
M (input) INTEGER
The number of system inputs. M >= 0.
P (input) INTEGER
The number of system outputs. P >= 0.
NR (input/output) INTEGER
On entry with ORDSEL = 'F', NR is the desired order of the
resulting reduced order system. 0 <= NR <= N.
On exit, if INFO = 0, NR is the order of the resulting
reduced order model. For a system with NU ALPHA-unstable
eigenvalues and NS ALPHA-stable eigenvalues (NU+NS = N),
NR is set as follows: if ORDSEL = 'F', NR is equal to
NU+MIN(MAX(0,NR-NU),NMIN), where NR is the desired order
on entry, and NMIN is the order of a minimal realization
of the ALPHA-stable part of the given system; NMIN is
determined as the number of Hankel singular values greater
than NS*EPS*HNORM(As,Bs,Cs), where EPS is the machine
precision (see LAPACK Library Routine DLAMCH) and
HNORM(As,Bs,Cs) is the Hankel norm of the ALPHA-stable
part of the given system (computed in HSV(1));
if ORDSEL = 'A', NR is the sum of NU and the number of
Hankel singular values greater than
MAX(TOL,NS*EPS*HNORM(As,Bs,Cs)).
ALPHA (input) DOUBLE PRECISION
Specifies the ALPHA-stability boundary for the eigenvalues
of the state dynamics matrix A. For a continuous-time
system (DICO = 'C'), ALPHA <= 0 is the boundary value for
the real parts of eigenvalues, while for a discrete-time
system (DICO = 'D'), 0 <= ALPHA <= 1 represents the
boundary value for the moduli of eigenvalues.
The ALPHA-stability domain does not include the boundary.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the state dynamics matrix A.
On exit, if INFO = 0, the leading NR-by-NR part of this
array contains the state dynamics matrix Ar of the reduced
order system.
The resulting A has a block-diagonal form with two blocks.
For a system with NU ALPHA-unstable eigenvalues and
NS ALPHA-stable eigenvalues (NU+NS = N), the leading
NU-by-NU block contains the unreduced part of A
corresponding to ALPHA-unstable eigenvalues in an
upper real Schur form.
The trailing (NR+NS-N)-by-(NR+NS-N) block contains
the reduced part of A corresponding to ALPHA-stable
eigenvalues.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
On entry, the leading N-by-M part of this array must
contain the original input/state matrix B.
On exit, if INFO = 0, the leading NR-by-M part of this
array contains the input/state matrix Br of the reduced
order system.
LDB INTEGER
The leading dimension of array B. LDB >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading P-by-N part of this array must
contain the original state/output matrix C.
On exit, if INFO = 0, the leading P-by-NR part of this
array contains the state/output matrix Cr of the reduced
order system.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
NS (output) INTEGER
The dimension of the ALPHA-stable subsystem.
HSV (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, the leading NS elements of HSV contain the
Hankel singular values of the ALPHA-stable part of the
original system ordered decreasingly.
HSV(1) is the Hankel norm of the ALPHA-stable subsystem.
</PRE>
<B>Tolerances</B>
<PRE>
TOL DOUBLE PRECISION
If ORDSEL = 'A', TOL contains the tolerance for
determining the order of reduced system.
For model reduction, the recommended value is
TOL = c*HNORM(As,Bs,Cs), where c is a constant in the
interval [0.00001,0.001], and HNORM(As,Bs,Cs) is the
Hankel-norm of the ALPHA-stable part of the given system
(computed in HSV(1)).
If TOL <= 0 on entry, the used default value is
TOL = NS*EPS*HNORM(As,Bs,Cs), where NS is the number of
ALPHA-stable eigenvalues of A and EPS is the machine
precision (see LAPACK Library Routine DLAMCH).
This value is appropriate to compute a minimal realization
of the ALPHA-stable part.
If ORDSEL = 'F', the value of TOL is ignored.
</PRE>
<B>Workspace</B>
<PRE>
IWORK INTEGER array, dimension (LIWORK)
LIWORK = 0, if JOB = 'B';
LIWORK = N, if JOB = 'N'.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The length of the array DWORK.
LDWORK >= MAX(1,N*(2*N+MAX(N,M,P)+5) + N*(N+1)/2).
For optimum performance LDWORK should be larger.
</PRE>
<B>Warning Indicator</B>
<PRE>
IWARN INTEGER
= 0: no warning;
= 1: with ORDSEL = 'F', the selected order NR is greater
than NSMIN, the sum of the order of the
ALPHA-unstable part and the order of a minimal
realization of the ALPHA-stable part of the given
system. In this case, the resulting NR is set equal
to NSMIN.
= 2: with ORDSEL = 'F', the selected order NR is less
than the order of the ALPHA-unstable part of the
given system. In this case NR is set equal to the
order of the ALPHA-unstable part.
</PRE>
<B>Error Indicator</B>
<PRE>
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: the computation of the ordered real Schur form of A
failed;
= 2: the separation of the ALPHA-stable/unstable diagonal
blocks failed because of very close eigenvalues;
= 3: the computation of Hankel singular values failed.
</PRE>
<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
<PRE>
Let be the following linear system
d[x(t)] = Ax(t) + Bu(t)
y(t) = Cx(t) (1)
where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1)
for a discrete-time system. The subroutine AB09MD determines for
the given system (1), the matrices of a reduced order system
d[z(t)] = Ar*z(t) + Br*u(t)
yr(t) = Cr*z(t) (2)
such that
HSV(NR+NS-N) <= INFNORM(G-Gr) <= 2*[HSV(NR+NS-N+1)+...+HSV(NS)],
where G and Gr are transfer-function matrices of the systems
(A,B,C) and (Ar,Br,Cr), respectively, and INFNORM(G) is the
infinity-norm of G.
The following procedure is used to reduce a given G:
1) Decompose additively G as
G = G1 + G2
such that G1 = (As,Bs,Cs) has only ALPHA-stable poles and
G2 = (Au,Bu,Cu) has only ALPHA-unstable poles.
2) Determine G1r, a reduced order approximation of the
ALPHA-stable part G1.
3) Assemble the reduced model Gr as
Gr = G1r + G2.
To reduce the ALPHA-stable part G1, if JOB = 'B', the square-root
Balance & Truncate method of [1] is used, and for an ALPHA-stable
continuous-time system (DICO = 'C'), the resulting reduced model
is balanced. For ALPHA-stable systems, setting TOL < 0, the
routine can be used to compute balanced minimal state-space
realizations.
If JOB = 'N', the balancing-free square-root version of the
Balance & Truncate method [2] is used to reduce the ALPHA-stable
part G1.
</PRE>
<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
<PRE>
[1] Tombs M.S. and Postlethwaite I.
Truncated balanced realization of stable, non-minimal
state-space systems.
Int. J. Control, Vol. 46, pp. 1319-1330, 1987.
[2] Varga A.
Efficient minimal realization procedure based on balancing.
Proc. of IMACS/IFAC Symp. MCTS, Lille, France, May 1991,
A. El Moudui, P. Borne, S. G. Tzafestas (Eds.),
Vol. 2, pp. 42-46.
</PRE>
<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
<PRE>
The implemented methods rely on accuracy enhancing square-root or
balancing-free square-root techniques.
3
The algorithms require less than 30N floating point operations.
</PRE>
<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
<PRE>
None
</PRE>
<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
<P>
<B>Program Text</B>
<PRE>
* AB09MD EXAMPLE PROGRAM TEXT
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX, PMAX
PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER LDA, LDB, LDC
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = NMAX*( 2*NMAX +
$ MAX( NMAX, MMAX, PMAX ) + 5 ) +
$ ( NMAX*( NMAX + 1 ) )/2 )
* .. Local Scalars ..
DOUBLE PRECISION ALPHA, TOL
INTEGER I, INFO, IWARN, J, M, N, NR, NS, P
CHARACTER*1 DICO, EQUIL, JOB, ORDSEL
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ DWORK(LDWORK), HSV(NMAX)
INTEGER IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL AB09MD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M, P, NR, ALPHA, TOL, DICO, JOB, EQUIL,
$ ORDSEL
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99990 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1, N )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99988 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
* Find a reduced ssr for (A,B,C).
CALL AB09MD( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR,
$ ALPHA, A, LDA, B, LDB, C, LDC, NS, HSV,
$ TOL, IWORK, DWORK, LDWORK, IWARN, INFO)
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) NR
WRITE ( NOUT, FMT = 99987 )
WRITE ( NOUT, FMT = 99995 ) ( HSV(J), J = 1,NS )
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, NR
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 40 I = 1, NR
WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
40 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 60 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,NR )
60 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' AB09MD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB09MD = ',I2)
99997 FORMAT (' The order of reduced model = ',I2)
99996 FORMAT (/' The reduced state dynamics matrix Ar is ')
99995 FORMAT (20(1X,F8.4))
99993 FORMAT (/' The reduced input/state matrix Br is ')
99992 FORMAT (/' The reduced state/output matrix Cr is ')
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
99988 FORMAT (/' P is out of range.',/' P = ',I5)
99987 FORMAT (/' The Hankel singular values of ALPHA-stable part are')
END
</PRE>
<B>Program Data</B>
<PRE>
AB09MD EXAMPLE PROGRAM DATA (Continuous system)
7 2 3 0 -.6D0 1.D-1 C N N A
-0.04165 0.0000 4.9200 -4.9200 0.0000 0.0000 0.0000
-5.2100 -12.500 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 3.3300 -3.3300 0.0000 0.0000 0.0000 0.0000
0.5450 0.0000 0.0000 0.0000 -0.5450 0.0000 0.0000
0.0000 0.0000 0.0000 4.9200 -0.04165 0.0000 4.9200
0.0000 0.0000 0.0000 0.0000 -5.2100 -12.500 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 3.3300 -3.3300
0.0000 0.0000
12.500 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 12.500
0.0000 0.0000
1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000
</PRE>
<B>Program Results</B>
<PRE>
AB09MD EXAMPLE PROGRAM RESULTS
The order of reduced model = 5
The Hankel singular values of ALPHA-stable part are
1.9178 0.8621 0.7666 0.0336 0.0246
The reduced state dynamics matrix Ar is
-0.5181 -1.1084 0.0000 0.0000 0.0000
8.8157 -0.5181 0.0000 0.0000 0.0000
0.0000 0.0000 0.5124 0.0000 1.7910
0.0000 0.0000 0.0000 -1.4460 0.0000
0.0000 0.0000 -4.2167 0.0000 -2.9900
The reduced input/state matrix Br is
-1.2837 1.2837
-0.7522 0.7522
-0.7447 -0.7447
1.9275 -1.9275
-3.6872 -3.6872
The reduced state/output matrix Cr is
-0.1380 -0.6445 -0.6582 -0.5771 0.2222
0.6246 0.0196 0.0000 0.4131 0.0000
0.1380 0.6445 -0.6582 0.5771 0.2222
</PRE>
<HR>
<p>
<A HREF=..\libindex.html><B>Return to index</B></A></BODY>
</HTML>
|
