Mercurial > octave-nkf

comparison scripts/control/base/place.m @ 7125:f084ba47812b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-11-08 02:29:23 by jwe]
author jwe
date Thu, 08 Nov 2007 02:29:24 +0000
parents a1dbe9d80eee
children 4a375de63f66
comparison
equal deleted inserted replaced
7124:d07cb867891b 7125:f084ba47812b
43 ## code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls 43 ## code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls
44 ## toolbox 44 ## toolbox
45 45
46 function K = place (sys, P) 46 function K = place (sys, P)
47 47
48 if (nargin != 2)
49 print_usage ();
50 endif
51
48 ## check arguments 52 ## check arguments
49 53
50 if(!isstruct(sys)) 54 if(! isstruct (sys))
51 error("sys must be in system data structure format (see ss)"); 55 error ("sys must be in system data structure format (see ss)");
52 endif 56 endif
53 sys = sysupdate(sys,"ss"); # make sure it has state space form up to date 57 sys = sysupdate (sys, "ss"); # make sure it has state space form up to date
54 if(!is_controllable(sys)) 58 if (! is_controllable (sys))
55 error("sys is not controllable."); 59 error ("sys is not controllable");
56 elseif( min(size(P)) != 1) 60 elseif (min (size (P)) != 1)
57 error("P must be a vector") 61 error ("P must be a vector")
58 else 62 else
59 P = reshape(P,length(P),1); # make P a column vector 63 P = P(:); # make P a column vector
60 endif 64 endif
61 ## system must be purely continuous or discrete 65 ## system must be purely continuous or discrete
62 is_digital(sys); 66 is_digital (sys);
63 [n,nz,m,p] = sysdimensions(sys); 67 [n, nz, m, p] = sysdimensions (sys);
64 nx = n+nz; # already checked that it's not a mixed system. 68 nx = n+nz; # already checked that it's not a mixed system.
65 if(m != 1) 69 if(m != 1)
66 error(["sys has ", num2str(m)," inputs; need only 1"]); 70 error ("sys has %d inputs; need only 1", m);
67 endif 71 endif
68 72
69 ## takes the A and B matrix from the system representation 73 ## takes the A and B matrix from the system representation
70 [A,B]=sys2ss(sys); 74 [A, B] = sys2ss (sys);
71 sp = length(P); 75 sp = length (P);
72 if(nx == 0) 76 if (nx == 0)
73 error("place: A matrix is empty (0x0)"); 77 error ("place: A matrix is empty (0x0)");
74 elseif(nx != length(P)) 78 elseif (nx != length (P))
75 error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ... 79 error ("A=(%dx%d), P has %d entries", nx, nx, length (P))
76 "entries."])
77 endif 80 endif
78 81
79 ## arguments appear to be compatible; let's give it a try! 82 ## arguments appear to be compatible; let's give it a try!
80 ## The second step is the calculation of the characteristic polynomial ofA 83 ## The second step is the calculation of the characteristic polynomial ofA
81 PC=poly(A); 84 PC = poly (A);
82 85
83 ## Third step: Calculate the transformation matrix T that transforms the state 86 ## Third step: Calculate the transformation matrix T that transforms the state
84 ## equation in the controllable canonical form. 87 ## equation in the controllable canonical form.
85 88
86 ## first we must calculate the controllability matrix M: 89 ## first we must calculate the controllability matrix M:
87 M=B; 90 M = B;
88 AA=A; 91 AA = A;
89 for n = 2:nx 92 for n = 2:nx
90 M(:,n)=AA*B; 93 M(:,n) = AA*B;
91 AA=AA*A; 94 AA = AA*A;
92 endfor 95 endfor
93 96
94 ## second, construct the matrix W 97 ## second, construct the matrix W
95 PCO=PC(nx:-1:1); 98 PCO = PC(nx:-1:1);
96 PC1=PCO; # Matrix to shift and create W row by row 99 PC1 = PCO; # Matrix to shift and create W row by row
97 100
98 for n = 1:nx 101 for n = 1:nx
99 W(n,:) = PC1; 102 W(n,:) = PC1;
100 PC1=[PCO(n+1:nx),zeros(1,n)]; 103 PC1 = [PCO(n+1:nx), zeros(1,n)];
101 endfor 104 endfor
102 105
103 T=M*W; 106 T = M*W;
104 107
105 ## finaly the matrix K is calculated 108 ## finaly the matrix K is calculated
106 PD = poly(P); # The desired characteristic polynomial 109 PD = poly (P); # The desired characteristic polynomial
107 PD = PD(nx+1:-1:2); 110 PD = PD(nx+1:-1:2);
108 PC = PC(nx+1:-1:2); 111 PC = PC(nx+1:-1:2);
109 112
110 K = (PD-PC)/T; 113 K = (PD-PC)/T;
111 114
112 ## Check if the eigenvalues of (A-BK) are the same specified in P 115 ## Check if the eigenvalues of (A-BK) are the same specified in P
113 Pcalc = eig(A-B*K); 116 Pcalc = eig (A-B*K);
114 117
115 Pcalc = sortcom(Pcalc); 118 Pcalc = sortcom (Pcalc);
116 P = sortcom(P); 119 P = sortcom (P);
117 120
118 if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1) 121 if (max ((abs(Pcalc)-abs(P))./abs(P) ) > 0.1)
119 disp("Place: Pole placed at more than 10% relative error from specified"); 122 warning ("place: Pole placed at more than 10% relative error from specified");
120 endif 123 endif
121 124
122 endfunction 125 endfunction
123 126