7017
|
1 ## Copyright (C) 1996, 1998, 2000, 2002, 2004, 2005, 2007 |
|
2 ## Auburn University. All rights reserved. |
3430
|
3 ## |
|
4 ## This file is part of Octave. |
|
5 ## |
|
6 ## Octave is free software; you can redistribute it and/or modify it |
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
9 ## your option) any later version. |
3430
|
10 ## |
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
14 ## General Public License for more details. |
3430
|
15 ## |
|
16 ## You should have received a copy of the GNU General Public License |
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
18 ## <http://www.gnu.org/licenses/>. |
3430
|
19 |
|
20 ## -*- texinfo -*- |
3502
|
21 ## @deftypefn {Function File} {[@var{xinf}, @var{x_ha_err}] =} hinfsyn_ric (@var{a}, @var{bb}, @var{c1}, @var{d1dot}, @var{r}, @var{ptol}) |
3439
|
22 ## Forms |
|
23 ## @example |
5016
|
24 ## xx = ([bb; -c1'*d1dot]/r) * [d1dot'*c1 bb']; |
|
25 ## Ha = [a 0*a; -c1'*c1 - a'] - xx; |
3439
|
26 ## @end example |
|
27 ## and solves associated Riccati equation. |
|
28 ## The error code @var{x_ha_err} indicates one of the following |
|
29 ## conditions: |
|
30 ## @table @asis |
|
31 ## @item 0 |
|
32 ## successful |
|
33 ## @item 1 |
3502
|
34 ## @var{xinf} has imaginary eigenvalues |
3439
|
35 ## @item 2 |
3502
|
36 ## @var{hx} not Hamiltonian |
3439
|
37 ## @item 3 |
3502
|
38 ## @var{xinf} has infinite eigenvalues (numerical overflow) |
3439
|
39 ## @item 4 |
3502
|
40 ## @var{xinf} not symmetric |
3439
|
41 ## @item 5 |
3502
|
42 ## @var{xinf} not positive definite |
3439
|
43 ## @item 6 |
3502
|
44 ## @var{r} is singular |
3430
|
45 ## @end table |
|
46 ## @end deftypefn |
|
47 |
|
48 function [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol) |
|
49 |
7125
|
50 if (nargin != 6) |
|
51 print_usage (); |
|
52 endif |
|
53 |
3430
|
54 x_ha_err = 0; # assume success |
|
55 Xinf = []; # default return value |
7125
|
56 n = issquare (A); |
|
57 nw = issquare (R); |
|
58 if (rank (R) != nw) |
|
59 x_ha_err = 6; |
3430
|
60 else # build hamiltonian Ha for X_inf |
|
61 xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB']; |
|
62 Ha = [A, 0*A; -C1'*C1, -A'] - xx; |
|
63 x_ha_err = 0; |
7125
|
64 [d, Ha] = balance (Ha); |
|
65 [u, s] = schur (Ha, "A"); |
|
66 rev = real (eig (s)); |
3430
|
67 |
7125
|
68 if (any (abs (rev) <= ptol)) # eigenvalues near the imaginary axis |
3430
|
69 x_ha_err = 1; |
7125
|
70 elseif (sum (rev > 0) != sum (rev < 0)) |
3430
|
71 ## unequal number of positive and negative eigenvalues |
|
72 x_ha_err = 2; |
|
73 else |
|
74 ## compute positive Riccati equation solution |
|
75 u = d * u; |
|
76 Xinf = u(n+1:2*n,1:n) / u(1:n,1:n); |
7125
|
77 if (! all (all (finite (Xinf)))) |
3430
|
78 x_ha_err = 3; |
7125
|
79 elseif (norm (Xinf-Xinf') >= 10*ptol) |
3430
|
80 ## solution not symmetric |
|
81 x_ha_err = 4; |
|
82 else |
|
83 ## positive semidefinite? |
|
84 ## force symmetry (faster, avoids some convergence problems) |
|
85 Xinf = (Xinf + Xinf')/2; |
7125
|
86 rev = eig (Xinf); |
|
87 if (any (rev <= -ptol)) |
3430
|
88 x_ha_err = 5; |
|
89 endif |
|
90 endif |
|
91 endif |
|
92 endif |
|
93 endfunction |