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
|
## Copyright (C) 1996 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.
## Usage: x = dare (a, b, c, r {,opt})
##
## Solves discrete-time algebraic riccati equation
##
## a' x a - x + a' x b (r + b' x b)^{-1} b' x a + c = 0
##
## for
##
## a: nxn
## b: nxm
## c: nxn, symmetric positive semidefinite
## r: mxm, invertible
##
## If c is not square, then the function attempts to use c'*c instead.
##
## Solution method: Laub's Schur method (IEEE Trans Auto Contr, 1979) applied
## to the appropriate symplectic matrix.
##
## See also: Ran and Rodman, "Stable Hermitian Solutions of Discrete
## Algebraic Riccati Equations," Mathematics of Control, Signals and
## Systems, Vol 5, no 2 (1992) pp 165-194.
##
## opt is an option passed to the eigenvalue balancing routine default
## is "B".
##
## See also: balance, are
## Author: A. S. Hodel <scotte@eng.auburn.edu>
## Created: August 1993
## Adapted-By: jwe
function x = dare (a, b, c, r, opt)
if (nargin == 4 || nargin == 5)
if (nargin == 5)
if (opt != "N" || opt != "P" || opt != "S" || opt != "B")
warning ("dare: opt has an invalid value -- setting to B");
opt = "B";
endif
else
opt = "B";
endif
## Check a matrix dimensions
if ((n = is_square (a)) == 0)
error ("dare: a is not square");
endif
## Check a,b compatibility.
[n1, m] = size (b);
if (n1 != n)
warning ("dare: a,b are not conformable");
endif
if (is_controllable (a, b) == 0)
warning ("dare: a,b are not controllable");
endif
## Check a,c compatibility.
if (is_observable (a, c) == 0)
warning ("dare: a,c are not observable");
endif
if ((p = is_square (c)) == 0)
c = c'*c;
p = rows (c);
endif
if (n != p)
error ("dare: a,c are not conformable");
endif
## Check r dimensions.
if ((m1 = is_square (r)) == 0)
warning ("dare: r is not square");
elseif (m1 != m)
warning ("b,r are not conformable");
endif
brb = (b/r)*b';
atc = a'\c;
[d, sy] = balance ([a + brb*atc, -brb/(a'); -atc, (inv (a'))], opt);
[u, s] = schur (sy, 'D');
u = d*u;
n1 = n+1;
n2 = 2*n;
x = u (n1:n2, 1:n)/u(1:n, 1:n);
else
usage ("x = dare (a, b, c, r {,opt})");
endif
endfunction
|
