|
3430
|
1 ## Copyright (C) 1996, 1998 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 -*- |
|
3502
|
20 ## @deftypefn {Function File} {[@var{xinf}, @var{x_ha_err}] =} hinfsyn_ric (@var{a}, @var{bb}, @var{c1}, @var{d1dot}, @var{r}, @var{ptol}) |
|
3439
|
21 ## Forms |
|
|
22 ## @example |
|
5016
|
23 ## xx = ([bb; -c1'*d1dot]/r) * [d1dot'*c1 bb']; |
|
|
24 ## Ha = [a 0*a; -c1'*c1 - a'] - xx; |
|
3439
|
25 ## @end example |
|
|
26 ## and solves associated Riccati equation. |
|
|
27 ## The error code @var{x_ha_err} indicates one of the following |
|
|
28 ## conditions: |
|
|
29 ## @table @asis |
|
|
30 ## @item 0 |
|
|
31 ## successful |
|
|
32 ## @item 1 |
|
3502
|
33 ## @var{xinf} has imaginary eigenvalues |
|
3439
|
34 ## @item 2 |
|
3502
|
35 ## @var{hx} not Hamiltonian |
|
3439
|
36 ## @item 3 |
|
3502
|
37 ## @var{xinf} has infinite eigenvalues (numerical overflow) |
|
3439
|
38 ## @item 4 |
|
3502
|
39 ## @var{xinf} not symmetric |
|
3439
|
40 ## @item 5 |
|
3502
|
41 ## @var{xinf} not positive definite |
|
3439
|
42 ## @item 6 |
|
3502
|
43 ## @var{r} is singular |
|
3430
|
44 ## @end table |
|
|
45 ## @end deftypefn |
|
|
46 |
|
|
47 function [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol) |
|
|
48 |
|
|
49 x_ha_err = 0; # assume success |
|
|
50 Xinf = []; # default return value |
|
4030
|
51 n = issquare(A); |
|
|
52 nw = issquare(R); |
|
3430
|
53 if(rank(R) != nw) x_ha_err = 6; |
|
|
54 else # build hamiltonian Ha for X_inf |
|
|
55 xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB']; |
|
|
56 Ha = [A, 0*A; -C1'*C1, -A'] - xx; |
|
|
57 x_ha_err = 0; |
|
|
58 [d, Ha] = balance(Ha); |
|
|
59 [u, s] = schur(Ha, "A"); |
|
|
60 rev = real(eig(s)); |
|
|
61 |
|
|
62 if (any(abs(rev) <= ptol)) # eigenvalues near the imaginary axis |
|
|
63 x_ha_err = 1; |
|
|
64 elseif (sum(rev > 0) != sum(rev < 0)) |
|
|
65 ## unequal number of positive and negative eigenvalues |
|
|
66 x_ha_err = 2; |
|
|
67 else |
|
|
68 ## compute positive Riccati equation solution |
|
|
69 u = d * u; |
|
|
70 Xinf = u(n+1:2*n,1:n) / u(1:n,1:n); |
|
|
71 if (!all(all(finite(Xinf)))) |
|
|
72 x_ha_err = 3; |
|
|
73 elseif (norm(Xinf-Xinf') >= 10*ptol) |
|
|
74 ## solution not symmetric |
|
|
75 x_ha_err = 4; |
|
|
76 else |
|
|
77 ## positive semidefinite? |
|
|
78 ## force symmetry (faster, avoids some convergence problems) |
|
|
79 Xinf = (Xinf + Xinf')/2; |
|
|
80 rev = eig(Xinf); |
|
|
81 if (any(rev <= -ptol)) |
|
|
82 x_ha_err = 5; |
|
|
83 endif |
|
|
84 endif |
|
|
85 endif |
|
|
86 endif |
|
|
87 endfunction |
