 help lqr

 % input the a and b matrices
 a = [-0.05156  -0.05877;
         0.01503  -0.005887];
 b = [0  -0.03384]';
 a
 b
 pause
 help eye

 % input the q and r matrices
 
 q = eye(2)

 r = 1
 pause
 % compute the Kalman gain and Riccati Matrix
 
 [k,p] = lqr(a,b,q,r)
 pause
 % use IMPULSE to  get the response of the linear system
 
 help impulse

 pause
 % redifine some matrices to put in proper form so that the impulse is
 % applied to both state variables with magnitude unity
 aa = a - b*k


 bb = [1 1]'
 c = eye(2)
 d = [0 0]'
 pause
 sys = ss(aa,bb,c,d)
 pause
 [y,t] = impulse(sys,100);
 u = -k*y';
 uu = u';
 plot(t,y,t,uu)
 pause
 title('Plot of optimal state and control vs time')
 xlabel ('Time, min')
 ylabel ('State vector components and control variable')
 text (20,0,'Concentration')
 text (20,.6,'Temperature')
 text (15,.4,'Control')
 pause
 % readjusting the R weighting
 r = .1
 [k,p] = lqr(a,b,q,r)
 pause

 aa = a-b*k

 sys = ss(aa,bb,c,d);
 clear y,t;
 [y,t] = impulse(sys,100);
 u = -k*y';
 uu = u';
 plot(t,y,t,uu)
 pause
 title('Plot of optimal state and control vs time')
 xlabel ('Time, min')
 ylabel ('State vector components and control variable')
 text(1,.1,'Conc')
 text(8,.7,'Temp')
 text(20,1.4,'Control')
 pause
 % redefining the Q and R weighting matrices
 q(1,1) = 10
 r = 1
 [k,p] = lqr(a,b,q,r)
 pause
 aa = a-b*k
 sys = ss(aa,bb,c,d);
 u = -k*y';
 uu = u';
 clear y,t;
 clear u uu;
 [y,t] = impulse(sys,100);
 u = -k*y';
 uu = u';
 plot(t,y,t,uu)
 pause
 title('Plot of optimal state and control vs time')
 xlabel ('Time, min')
 ylabel ('State vector components and control variable')
 text (1,-.1,'Conc')
 text(20,.15,'Temp')
 text(40,.3,'Control')
 pause


