1. % 3)
2.  
3. echo on,
4. clc,
5. % wersja 12.04.2013
6. % najpierw uruchom m-plik SScz01.m
7. % SScz01
8. %-------------------------------------
9. % model w przestrzeni stanu
10. %-------------------------------------
11. % dane - macierze - forma kanonicna sterowalna - zmiennych fazowych
12. A = [0,1,0; 0,0,1; -12,-18.4,-8.2]  % macierz stanu
13. B = [0; 0; 1]                       % macierz wejścia
14. C = [8 0 2]                         % macierz wyjścia
15. D = [0]                             % macierz przejścia
16. pause
17. %-------------------------------------
18. sys = ss(A,B,C,D)               % model przestrzeni stanu
19. pause
20. %-------------------------------------
21. % funkcja eig() - wartości własne
22. eigA = eig(A)
23. pause
24. %-------------------------------------
25. % funkcja ss2zp()
26. [zss,bss,wss] = ss2zp(A,B,C,D)
27. pause
28. %-------------------------------------
29. % funkcja tfdata()
30. [lG,mG] = tfdata(sysG,'v')
31. pause
32. %-------------------------------------
33. [zG,bG,wG] = tf2zp(lG,mG)
34. pause
35. %-------------------------------------
36.  
37. % 5)
38.  
39. % odpowiedzi na step()
40. pause
41. %t = [0:0.01:12];
42. [y,t,x] = step(sys); [yc,tc,xc] = step(sysc); yG = step(sysG,t);
43. %-------------------------------------
44. disp('!!! Przebiegi czasowej odpowiedzi y(t) i zmiennych stanu x(t) !!!')
45. disp('!!! Naciśnij ENTER !!!')
46. pause
47. %-------------------------------------
48. figure
49. subplot(2,2,1), plot(t,y), grid, title('y(t)-ABCD')  %, axis([0,max(t),0,0.4])
50. subplot(2,2,3), plot(t,x), grid, title('x(t)-ABCD')  %, axis([0,max(t),0,1.1])
51. subplot(2,2,2), plot(tc,yc), grid, title('y(t)-AcBcCcDc')%, axis([0,max(t),0,0.4])
52. subplot(2,2,4), plot(tc,xc), grid, title('x(t)-AcBcCcDc')%, axis([0,max(t),0,1.1])
53. figure
54. plot(t,y,tc,yc,'x',t,yG,'.'), grid              %, axis([0,max(t),0,0.4])
55. legend('step()-y-ABCD','step()-yc-AcBcCcDc','step()-yG-G(s)')
56. pause, close all
57. %-------------------------------------
58. % niezerowe warunki początkowe -> x0 = [0.5; 0; 1]
59. pause
60. tu=[0:0.1:max(t)]';
61. x0 = [0.5; 0; 1];
62. u = ones(length(tu),1);
63. %u = 2*ones(20,1); u(21:201,1) = 0.5*ones(181,1);
64. [yu,xu] = lsim(A,B,C,D,u',tu,x0);
65. %plot(t,u,t,xu)
66. figure
67. subplot(2,2,1), plot(t,y), grid, title('y(t) x0=[0;0;0]')  %, axis([0,max(t),0,0.4])
68. subplot(2,2,3), plot(t,x), grid, title('x(t) x0=[0;0;0]')  %, axis([0,max(t),0,1.1])
69. subplot(2,2,2), plot(tu,yu), grid, title('yu(t) x0=[0.5;0;1]')%, axis([0,max(t),0,0.4])
70. subplot(2,2,4), plot(tu,xu), grid, title('xu(t) x0=[0.5;0;1]')%, axis([0,max(t),0,1.1])
71. figure
72. plot(t,y,'b',tu,yu,'r'), grid %, axis([0,max(t),0,0.4])
73. legend('step() - y','step() - yu')
74. pause, close all
75. %-------------------------------------
76.  
77. % 6)
78.  
79. % postać kanoniczna: "modal" i "companion"
80. % "modal" - macierz Am zawiera tylko elementy rzeczywiste
81. % "companion" - forma kanoniczna obserwowalna
82. % funcja canon()
83. [Am,Bm,Cm,Dm] = canon(A,B,C,D,'modal')
84. %[Amx,Bmx,Cmx,Dmx] = canon(Ac,Bc,Cc,Dc,'modal')
85. pause
86. %-------------------------------------
87. [Amc,Bmc,Cmc,Dmc] = canon(A,B,C,D,'companion')
88. %[Amcx,Bmcx,Cmcx,Dmcx] = canon(Ac,Bc,Cc,Dc,'companion')
89. % Ac' - odpowiada Amc - tu tak jest, ale dla Acx=inv(P1)*Amc*P1
90. % Acx'=Amc, ale A'=Amc
91. pause
92. %-------------------------------------
93. % macierz transformacyjna T = inv(P)
94. % nowe współrzędne stanu z = inv(P)*x = T*x
95. [Am,Bm,Cm,Dm,Tm] = canon(A,B,C,D,'modal')
96. %[Amx,Bmx,Cmx,Dmx,Tmx] = canon(Ac,Bc,Cc,Dc,'modal')
97. %[sysm,Tm] = canon(sys,'modal')
98. [Amc,Bmc,Cmc,Dmc,Tmc] = canon(A,B,C,D,'companion')
99. %[Amcx,Bmcx,Cmcx,Dmcx,Tmcx] = canon(Ac,Bc,Cc,Dc,'companion')
100. pause
101. %-------------------------------------
102. % porównanie wartości własnych
103. [eig(Ac), eig(A), eig(Am), eig(Amc)]
104. pause
105. %-------------------------------------
106.  
107. % 7)
108.  
109. % transformacja do nowego układu współrzędnych stanu
110. P  = inv(Tm);
111. Az = inv(P)*A*P;
112. Bz = inv(P)*B;
113. Cz = C*P;
114. Dz = D;
115. pause
116. %-------------------------------------
117. Pmc  = inv(Tmc);
118. Azmc = inv(Pmc)*A*Pmc;
119. Bzmc = inv(Pmc)*B;
120. Czmc = C*Pmc;
121. Dzmc = D;
122. pause
123. %-------------------------------------
124. % porównanie wartości własnych
125. [eig(A), eig(Am), eig(Amc), eig(Az), eig(Azmc)]
126. pause
127. %-------------------------------------
128. % transformacja tożsamościowa (modalna) - diagonalizacja A
129. % Ad = macierz.wekt.wł.A^(-1)*A*mac.wek.wł.A -> macierz A będzie diagonalna
130. % Bd = inv(macierz.wekt.wł.A)*B
131. % Cd = C*macierz.wekt.wł.A
132. % Dd = D
133. pause
134. %-------------------------------------
135. % M macierz wektorów własnych
136. % Ad macierz wartosci własnych (diagonalna)
137. [M Ad] = eig(A)                                        
138. pause
139. %-------------------------------------
140. Bd = inv(M)*B    % transformowana macierz wejścia
141. Cd = C*M         % transformowana macierz wyjścia
142. Dd = D           % transformowana macierz przejścia
143. pause
144. %-------------------------------------
145. % eigAd = eig(Ad)
146. [eig(Az), eig(Azmc), eig(Ad), eig(inv(M)*A*M)]
147. pause
148. %-------------------------------------
149. % funkcja ss2tf()
150. [numGm,denGm]   = ss2tf(Am,Bm,Cm,Dm);
151. [numGmc,denGmc] = ss2tf(Amc,Bmc,Cmc,Dmc);
152. [numGd,denGd]   = ss2tf(Ad,Bd,Cd,Dd);
153. %-------------------------------------
154. printsys(numGm,denGm)
155. printsys(numGmc,denGmc)
156. printsys(numGd,denGd)
157. pause
158. %-------------------------------------
159. % funkcja roundaf() - zaokrąglenie
160. numGm = roundaf(numGm); denGm = roundaf(denGm);
161. numGmc = roundaf(numGmc); denGmc = roundaf(denGmc);
162. numGd = roundaf(numGd); denGd = roundaf(denGd);
163. printsys(numGm,denGm)
164. printsys(numGmc,denGmc)
165. printsys(numGd,denGd)
166. pause
167. %-------------------------------------
168. % odwrócenie zmiennych - przejście między Ac i A,
169. % Amc->AmcP1 - zamiana kolumn
170. % AmcP1=inv(P1)*Amc*P1
171. P1=[0 0 1; 0 1 0; 1 0 0]
172. AcP1=inv(P1)*Ac*P1
173. BcP1=inv(P1)*Bc
174. CcP1=Cc*P1
175. DcP1=Dc
176. eigAcP1=eig(AcP1)
177. pause
178. %-------------------------------------
179.  
180. % 8)
181.  
182. % sterowalność
183. Ster = ctrb(A,B)
184. rzad_Mc = rank(Ster)
185. pause
186. %-------------------------------------
187. % obserwowalność
188. Obser = obsv(A,C)
189. rzad_Mo = rank(Obser)
190. pause
191. %-------------------------------------