home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / BOILDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  12.0 KB  |  302 lines
  1. echo on
  2. % This file demonstrates MATLAB's ability for LQG control system design by 
  3. % going through the design of a controller for a boiler system.
  4.  
  5. echo off
  6. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  7. echo on
  8.  
  9. pause % Press any key to continue ...
  10.  
  11. % Define boiler system
  12. A = [ -3.93 -3.15e-3 0 0 0 4.03e-5 0 0 0
  13.       3.68e2 -3.05 3.03 0 0 -3.77e-3 0 0 0
  14.       2.74e1  7.87e-2 -5.96e-2 0 0 -2.81e-4 0 0 0
  15.       -6.47e-2 -5.2e-5 0 -2.55e-1 -3.35e-6 3.6e-7 6.33e-5 1.94e-4 0
  16.       3.85e3 1.73e1 -1.28e1 -1.26e4 -2.91 -1.05e-1 1.27e1 4.31e1 0
  17.       2.24e4 1.8e1 0 -3.56e1 -1.04e-4 -4.14e-1 9e1 5.69e1 0
  18.       0 0 2.34e-3 0 0 2.22e-4 -2.03e-1 0 0
  19.       0 0 0 -1.27 1e-3 7.86e-5 0 -7.17e-2 0
  20.       -2.2 -1.77e-3 0 -8.44 -1.11e-4 1.38e-5 1.49e-3 6.02e-3 -1e-10];
  21. B = [ 0 0 0 -1.0e-2
  22.       0 0 0 0
  23.       1.56 0 0 0
  24.       0 -5.13e-6 0 0
  25.       8.28 -1.5 3.95e-2 5.2e1
  26.       0 1.78 0 0
  27.       2.33 0 0 0
  28.       0 -2.45e-2 2.84e-3 0
  29.       0 2.94e-5 0 0];
  30. C = [0 0 0 0 0 1 0 0 0 
  31.      0 0 0 0 0 0 0 0 1];
  32. D = [0 0 0 0
  33.      0 0 0 0];
  34.  
  35. states='x1 x2 x3 x4 x5 x6 x7 x8 x9';
  36. inputs='fuel-flow water-flow water-temp steam-flow';
  37. outputs='water-level pressure';
  38. printsys(A,B,C,D,inputs,outputs,states);
  39.  
  40. pause % Press any key to continue ...
  41.  
  42. % This is the state space model for a boiler system consisting of a superheater
  43. % and a riser in cascade.  It is desired to regulate the water-level and
  44. % pressure to desired set points.  The control inputs are fuel-flow,water-flow
  45. % and water-temp.  The steam-flow is a disturbance.  This system has one very
  46. % slow pole at -1.e-10 and is non-minimum phase.
  47.  
  48. % Let's look at some step responses to get a feel for how the system behaves.
  49. t = 0:5:500;        % Define time vector (automatic vector too long here)
  50. %subplot(221), step(A,B,C(1,:),D(1,:),1,t); title('Input 1 Output 1')
  51. %subplot(222), step(A,B,C(2,:),D(2,:),1,t); title('Input 1 Output 2')
  52. %subplot(223), step(A,B,C(1,:),D(1,:),2,t); title('Input 2 Output 1')
  53. %subplot(224), step(A,B,C(2,:),D(2,:),2,t); title('Input 2 Output 2'), pause
  54. %clg
  55. %subplot(221), step(A,B,C(1,:),D(1,:),3,t); title('Input 3 Output 1')
  56. %subplot(222), step(A,B,C(2,:),D(2,:),3,t); title('Input 3 Output 2')
  57. %subplot(223), step(A,B,C(1,:),D(1,:),4,t); title('Input 4 Output 1')
  58. %subplot(224), step(A,B,C(2,:),D(2,:),4,t); title('Input 4 Output 2')
  59. %pause
  60. echo off
  61. subplot(221), step(A,B,C(1,:),D(1,:),1,t); title('Input 1 Output 1')
  62. subplot(222), step(A,B,C(2,:),D(2,:),1,t); title('Input 1 Output 2')
  63. subplot(223), step(A,B,C(1,:),D(1,:),2,t); title('Input 2 Output 1')
  64. subplot(224), step(A,B,C(2,:),D(2,:),2,t); title('Input 2 Output 2'), pause
  65. clg
  66. subplot(221), step(A,B,C(1,:),D(1,:),3,t); title('Input 3 Output 1')
  67. subplot(222), step(A,B,C(2,:),D(2,:),3,t); title('Input 3 Output 2')
  68. subplot(223), step(A,B,C(1,:),D(1,:),4,t); title('Input 4 Output 1')
  69. subplot(224), step(A,B,C(2,:),D(2,:),4,t); title('Input 4 Output 2')
  70. pause
  71. echo on
  72.  
  73. % The time responses show the effect of the slow pole.  The water-level 
  74. % (output 1) responds in 500 seconds but the pressure shows ramp-like behavior
  75. % over this time frame.
  76.  
  77. pause % Press any key to continue ...
  78.  
  79. % Since we plan to use LQG control methods, let's look at the singular value
  80. % responses for this system to determine the dynamic controllability.
  81. subplot(111)
  82. sigma(A,B,C,D), pause % Press any key after plot
  83.  
  84. % The frequency responses show a lot of gain at low frequencies so that
  85. % achieving acceptable tracking of input commands should be easy. 
  86.  
  87. % To start lets design a full state feedback controller.  MATLAB contains two
  88. % commands to do full state feedback controller design, LQR and LQRY.
  89. % They differ only in the the specific cost function used.  LQR uses a cost
  90. % function that weights the states and LQRY uses a cost function that weights
  91. % the outputs.  We will use LQRY below.  This means that we must form a set
  92. % of criteria outputs to be weighted in the cost function.  We will use 
  93. % the plant outputs as criteria outputs.
  94.  
  95. pause % Press any key to continue ...
  96.  
  97. % First extract a plant containing the control inputs and plant outputs.
  98. [a,b,c,d] = ssselect(A,B,C,D,[1:3],[1:2]);
  99.  
  100. % Perform the LQR gain matrix calculation.  Use an "optimal root locus" to 
  101. % chose the "best" relative weighting between the Q and R matrices.  Note
  102. % the weights below were chosen after a few design iterations.
  103. Q = [.001 0
  104.      0   10];    % Output weighting matrix
  105. R = [1 0 0
  106.      0 .2 0
  107.      0 0  1];    % Input weighting matrix
  108.  
  109. [nx,na] = size(a);
  110. E = zeros(nx,11);
  111.  
  112. for i = 1:10;
  113.   rho = i*.1;
  114.   [K,S,E(:,i)] = lqry(a,b,c,d,rho*Q,R);
  115. end
  116. plot(real(E),imag(E),'x'), title('Optimal Root Locus'), pause % Press any key
  117.  
  118. % Using the last design, the closed loop eigenvalues are
  119. damp(esort(E(:,10)));
  120.  
  121. pause % Press any key to continue ...
  122.  
  123. % Now analyze the design.  To do this we must form the closed loop system.
  124. %      --->O---[Plant]--+--->
  125. %          |            | x
  126. %          +-----[K]----+
  127. % 
  128. % So append the states as outputs, and do a feedback connection. 
  129. [ac,bc,cc,dc] = augstate(a,b,c,d);
  130. inputs = [1:3]; outputs=[3:11];    % Define selection vectors
  131. [ac,bc,cc,dc] = feedback(ac,bc,cc,dc,[],[],[],K,-inputs,outputs);
  132. [ac,bc,cc,dc] = ssdelete(ac,bc,cc,dc,[],outputs); % Remove state outputs
  133.  
  134. % Look at the closed-loop time response.
  135. %subplot(221), step(ac,bc,cc(1,:),dc(1,:),1,t); title('Input 1 Output 1')
  136. %subplot(222), step(ac,bc,cc(2,:),dc(2,:),1,t); title('Input 1 Output 2')
  137. %subplot(223), step(ac,bc,cc(1,:),dc(1,:),2,t); title('Input 2 Output 1')
  138. %subplot(224), step(ac,bc,cc(2,:),dc(2,:),2,t); title('Input 2 Output 2'), pause
  139. echo off
  140. subplot(221), step(ac,bc,cc(1,:),dc(1,:),1,t); title('Input 1 Output 1')
  141. subplot(222), step(ac,bc,cc(2,:),dc(2,:),1,t); title('Input 1 Output 2')
  142. subplot(223), step(ac,bc,cc(1,:),dc(1,:),2,t); title('Input 2 Output 1')
  143. subplot(224), step(ac,bc,cc(2,:),dc(2,:),2,t); title('Input 2 Output 2'), pause
  144. echo on
  145.  
  146. % The time responses show that the water-level is controlled quickly, while
  147. % the pressure takes about 11 min (660 secs) to achieve steady state.
  148.  
  149. pause % Press any key to continue ...
  150.  
  151. % Let's plot the step response for the first 10 secs to see how the 
  152. % water-level behaves.
  153. tq = 0:.1:10;
  154. clg
  155. %subplot(221), step(ac,bc,cc(1,:),dc(1,:),1,tq); title('Input 1 Output 1')
  156. %subplot(222), step(ac,bc,cc(2,:),dc(2,:),1,tq); title('Input 1 Output 2')
  157. %subplot(223), step(ac,bc,cc(1,:),dc(1,:),2,tq); title('Input 2 Output 1')
  158. %subplot(224), step(ac,bc,cc(2,:),dc(2,:),2,tq); title('Input 2 Output 2')
  159. %pause % Press any key after plot ...
  160. echo off
  161. subplot(221), step(ac,bc,cc(1,:),dc(1,:),1,tq); title('Input 1 Output 1')
  162. subplot(222), step(ac,bc,cc(2,:),dc(2,:),1,tq); title('Input 1 Output 2')
  163. subplot(223), step(ac,bc,cc(1,:),dc(1,:),2,tq); title('Input 2 Output 1')
  164. subplot(224), step(ac,bc,cc(2,:),dc(2,:),2,tq); title('Input 2 Output 2')
  165. pause % Press any key to continue ...
  166. echo on
  167.  
  168. % The behavior looks good.  
  169.  
  170. % To recap, we've designed a state feedback gain matrix that requires all 
  171. % the states to be fedback.  However, we only have the two outputs so it will
  172. % be necessary to design a state estimator for our controller.  MATLAB
  173. % contains two commands to do this design, LQE and LQEW.  They only differ in
  174. % that LQEW allows direct feedthrough of process noise into the sensors.
  175. % Since our D matrix is zero we will use LQE.
  176.  
  177. % For estimator design we need to identify the sensor outputs and process 
  178. % noise inputs.  The sensors will be our two outputs, and we will use the
  179. % control inputs as well as the steam flow (output 4) as process noise inputs.
  180. % The control inputs are included as process noise inputs since to do so helps
  181. % to make the system robust to parameter changes at the control inputs.
  182.  
  183. pause % Press any key to continue ...
  184.  
  185. % Estimator design model
  186. a = A; b = B; c = C; d = D;
  187.  
  188. % Perform the LQE gain matrix calculation.  Use an "optimal root locus" to 
  189. % chose the "best" relative weighting between the Q and R matrices.
  190.  
  191. % Define process noise and sensor noise covariance.
  192. Q = [.1 0 0 0
  193.      0 100 0 0
  194.      0 0 1 0
  195.      0 0 0 100];    % Process noise covariance
  196. R = [1000 0
  197.      0 1];    % Sensor noise covariance
  198.  
  199. [nx,na] = size(a);
  200. E = zeros(nx,11);
  201.  
  202. for i = 1:10;
  203.   rho = i*.1;
  204.   [L,P,E(:,i)] = lqe(a,b,c,rho*Q,R);
  205. end
  206. subplot(111)
  207. plot(real(E),imag(E),'x'), title('Optimal Root Locus'), pause % Press any key
  208.  
  209. % Using the last design, the closed loop eigenvalues are
  210. damp(esort(E(:,10)));
  211.  
  212. pause % Press any key to continue ...
  213.  
  214. % Now form the controller so that we can form the closed-loop system.
  215. sensors= [1:2];  controls=[1:3];
  216. [ak,bk,ck,dk] = reg(A,B,C,D,K,L,sensors,[],controls);
  217.  
  218. % Form the closed loop system so that the inputs are output commands, yc:
  219. %  yc --->O--[Controller]---[Plant]--+---> y
  220. %         |                          |
  221. %         +--------------------------+
  222. % Use a combination of SERIES and CLOOP.
  223. [ac,bc,cc,dc] = series(ak,bk,ck,dk,A,B,C,D,[1:3],controls);
  224. [ac,bc,cc,dc] = cloop(ac,bc,cc,dc);
  225.  
  226. pause % Press any key to continue ...
  227.  
  228. % Closed loop eigenvalues.  Should be a combination of the closed-loop feedback
  229. % and estimator eigenvalues as guaranteed by the separation theorem.
  230. damp(ac)
  231.  
  232. pause % Press any key to continue ...
  233.  
  234. % Now plot the step response to water-level-command and pressure-command
  235. % (inputs 1 and 2 of the closed-loop system)
  236. step(ac,bc,cc,dc,1); title('Water-Level-Command'), pause % press any key ...
  237.  
  238. % The plot shows that the water is controlled quickly (in about 5 seconds)
  239. % and that after about 20 seconds the steady state error is small.  In addition
  240. % the pressure response is very small and is thus decoupled from this input.
  241.  
  242. pause % Press any key to continue ...
  243.  
  244. % Plot the response to pressure command
  245. step(ac,bc,cc,dc,2); title('Pressure-Command'), pause % Press any key ...
  246.  
  247. % This plot shows that the pressure achieves its steady state value in
  248. % about 50 minutes (3000 seconds) and that the water level is adversely
  249. % affected during this transient.  One way to remove this transient is to
  250. % use integral control on the water-level.  The integral control would 
  251. % help to decouple the pressure response from the water-level response.
  252.  
  253. pause % Press any key to continue ...
  254.  
  255. % Lets look at how a constant steam-flow disturbance can affect our steady 
  256. % state tracking.  Form a new closed loop system with the plant inputs as
  257. % inputs.   --->O----[Plant]---+--->
  258. %               |              |
  259. %               +-[Controller]-+
  260. [ac,bc,cc,dc] = feedback(A,B,C,D,ak,bk,ck,dk,-[1:3],[1:2]);
  261.  
  262. % Do step response from steam-flow (input 4)
  263. step(ac,bc,cc,dc,4); title('Steam-flow Disturbance'), pause % Press any key
  264.  
  265. % So the steam-flow affects the water-level substantially.  Again an integrator
  266. % on the water-level would provide better disturbance rejection.  
  267.  
  268. pause % Press any key to continue ...
  269.  
  270. % Although we will not redo the design, to include the integrator in the 
  271. % design we would:
  272.  
  273. % 1) Add integrator to the plant model by doing a series connection with the
  274. %    plant and the integrator:
  275. %    [ai,bi,ci,di] = zp2ss([],[0],[1]); % Integrator model
  276.  
  277. % 2) form error output by adding command input:
  278. %    B = [B,zeros(9,1)]; D = [D,[-1;0]];
  279.  
  280. % 3) Connect integrator in series using APPEND and CLOOP since SERIES removes
  281. % the outputs of system1 and the inputs of system2 from the resulting model.
  282. % The integral output will now be output 3.
  283. %   [a,b,c,d] = append(A,B,C,D,ai,bi,ci,di);
  284. %   [a,b,c,d] = cloop(a,b,c,d,[1],[5]);
  285.  
  286. % 4) Design the feedback gain matrix as above but also weight the integrator.
  287.  
  288. pause % Press any key to continue ...
  289.  
  290. % 5) Design the estimator exactly as above without the integrator.  The
  291. % estimator formed will be a reduced order estimator since the integrator 
  292. % state will not be estimated.  
  293.  
  294. % 6) So before forming the controller add an extra row of zeros to the Kalman 
  295. % gain matrix for the integral state.  
  296.  
  297. % 7) Finally when forming the controller use the plant containing the integral
  298. % and specify the integral command (input 4) as a known input.
  299.  
  300. echo off
  301.  
  302.