|
3432
|
1 ## Copyright (C) 1997 Jose Daniel Munoz Frias |
|
|
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} {} place (@var{sys}, @var{p}) |
|
3432
|
21 ## Computes the matrix K such that if the state |
|
|
22 ## is feedback with gain K, then the eigenvalues of the closed loop |
|
3502
|
23 ## system (i.e. A-BK) are those specified in the vector @var{p}. |
|
3432
|
24 ## |
|
|
25 ## Version: Beta (May-1997): If you have any comments, please let me know. |
|
|
26 ## (see the file place.m for my address) |
|
|
27 ## @end deftypefn |
|
|
28 |
|
|
29 ## Author: Jose Daniel Munoz Frias |
|
|
30 |
|
|
31 ## Universidad Pontificia Comillas |
|
|
32 ## ICAIdea |
|
|
33 ## Alberto Aguilera, 23 |
|
|
34 ## 28015 Madrid, Spain |
|
|
35 ## |
|
|
36 ## E-Mail: daniel@dea.icai.upco.es |
|
|
37 ## |
|
|
38 ## Phone: 34-1-5422800 Fax: 34-1-5596569 |
|
|
39 ## |
|
|
40 ## Algorithm taken from "The Control Handbook", IEEE press pp. 209-212 |
|
|
41 ## |
|
|
42 ## code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls |
|
|
43 ## toolbox |
|
|
44 |
|
|
45 function K = place (sys, P) |
|
|
46 |
|
|
47 sav_val = empty_list_elements_ok; |
|
|
48 empty_list_elements_ok = 1; |
|
|
49 |
|
|
50 ## check arguments |
|
|
51 |
|
4030
|
52 if(!isstruct(sys)) |
|
3432
|
53 error("sys must be in system data structure format (see ss2sys)"); |
|
|
54 endif |
|
|
55 sys = sysupdate(sys,"ss"); # make sure it has state space form up to date |
|
|
56 if(!is_controllable(sys)) |
|
|
57 error("sys is not controllable."); |
|
|
58 elseif( min(size(P)) != 1) |
|
|
59 error("P must be a vector") |
|
|
60 else |
|
|
61 P = reshape(P,length(P),1); # make P a column vector |
|
|
62 endif |
|
|
63 ## system must be purely continuous or discrete |
|
|
64 is_digital(sys); |
|
|
65 [n,nz,m,p] = sysdimensions(sys); |
|
|
66 nx = n+nz; # already checked that it's not a mixed system. |
|
|
67 if(m != 1) |
|
|
68 error(["sys has ", num2str(m)," inputs; need only 1"]); |
|
|
69 endif |
|
|
70 |
|
|
71 ## takes the A and B matrix from the system representation |
|
|
72 [A,B]=sys2ss(sys); |
|
|
73 sp = length(P); |
|
|
74 if(nx == 0) |
|
|
75 error("place: A matrix is empty (0x0)"); |
|
|
76 elseif(nx != length(P)) |
|
|
77 error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ... |
|
|
78 "entries."]) |
|
|
79 endif |
|
|
80 |
|
|
81 ## arguments appear to be compatible; let's give it a try! |
|
|
82 ## The second step is the calculation of the characteristic polynomial ofA |
|
|
83 PC=poly(A); |
|
|
84 |
|
|
85 ## Third step: Calculate the transformation matrix T that transforms the state |
|
|
86 ## equation in the controllable canonical form. |
|
|
87 |
|
|
88 ## first we must calculate the controllability matrix M: |
|
|
89 M=B; |
|
|
90 AA=A; |
|
|
91 for n = 2:nx |
|
|
92 M(:,n)=AA*B; |
|
|
93 AA=AA*A; |
|
|
94 endfor |
|
|
95 |
|
|
96 ## second, construct the matrix W |
|
|
97 PCO=PC(nx:-1:1); |
|
|
98 PC1=PCO; # Matrix to shift and create W row by row |
|
|
99 |
|
|
100 for n = 1:nx |
|
|
101 W(n,:) = PC1; |
|
|
102 PC1=[PCO(n+1:nx),zeros(1,n)]; |
|
|
103 endfor |
|
|
104 |
|
|
105 T=M*W; |
|
|
106 |
|
|
107 ## finaly the matrix K is calculated |
|
|
108 PD = poly(P); # The desired characteristic polynomial |
|
|
109 PD = PD(nx+1:-1:2); |
|
|
110 PC = PC(nx+1:-1:2); |
|
|
111 |
|
|
112 K = (PD-PC)/T; |
|
|
113 |
|
|
114 ## Check if the eigenvalues of (A-BK) are the same specified in P |
|
|
115 Pcalc = eig(A-B*K); |
|
|
116 |
|
|
117 Pcalc = sortcom(Pcalc); |
|
|
118 P = sortcom(P); |
|
|
119 |
|
|
120 if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1) |
|
|
121 disp("Place: Pole placed at more than 10% relative error from specified"); |
|
|
122 endif |
|
|
123 |
|
|
124 empty_list_elements_ok = sav_val; |
|
|
125 endfunction |
|
|
126 |
