|
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 |
|
|
50 x_ha_err = 0; # assume success |
|
|
51 Xinf = []; # default return value |
|
4030
|
52 n = issquare(A); |
|
|
53 nw = issquare(R); |
|
3430
|
54 if(rank(R) != nw) x_ha_err = 6; |
|
|
55 else # build hamiltonian Ha for X_inf |
|
|
56 xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB']; |
|
|
57 Ha = [A, 0*A; -C1'*C1, -A'] - xx; |
|
|
58 x_ha_err = 0; |
|
|
59 [d, Ha] = balance(Ha); |
|
|
60 [u, s] = schur(Ha, "A"); |
|
|
61 rev = real(eig(s)); |
|
|
62 |
|
|
63 if (any(abs(rev) <= ptol)) # eigenvalues near the imaginary axis |
|
|
64 x_ha_err = 1; |
|
|
65 elseif (sum(rev > 0) != sum(rev < 0)) |
|
|
66 ## unequal number of positive and negative eigenvalues |
|
|
67 x_ha_err = 2; |
|
|
68 else |
|
|
69 ## compute positive Riccati equation solution |
|
|
70 u = d * u; |
|
|
71 Xinf = u(n+1:2*n,1:n) / u(1:n,1:n); |
|
|
72 if (!all(all(finite(Xinf)))) |
|
|
73 x_ha_err = 3; |
|
|
74 elseif (norm(Xinf-Xinf') >= 10*ptol) |
|
|
75 ## solution not symmetric |
|
|
76 x_ha_err = 4; |
|
|
77 else |
|
|
78 ## positive semidefinite? |
|
|
79 ## force symmetry (faster, avoids some convergence problems) |
|
|
80 Xinf = (Xinf + Xinf')/2; |
|
|
81 rev = eig(Xinf); |
|
|
82 if (any(rev <= -ptol)) |
|
|
83 x_ha_err = 5; |
|
|
84 endif |
|
|
85 endif |
|
|
86 endif |
|
|
87 endif |
|
|
88 endfunction |
