% input the a, b, and w matrices
a = [-0.033181 -0.002511;
   0.05237 0.010216];

b = [0 -0.0988]';

w = [0.030864 0;
   0 0.03086];

% check the eigenvalues of a

eig(a)

% the positive eigenvalue shows that the system is unstable

% let's look at the response of the uncontrolled system

% use a unit impulse on each state

bb = [1 1]';
c = eye(2);
t = [0:1:100];
d = [0 0]';
y = impulse(a,bb,c,d,1,t);
size (y)

x1 = y(1:101,1);
x2 = y(1:101,2);

plot (t,y) 
title ('Uncontrolled Response')
xlabel('Time')
ylabel('State Vector Components')
text(10,0,'Conc')
text(20,3,'Temp')
pause
print uncontrol

% compute the Kalman feedback gain matrix for the sysem

q = [1 0;
   0 1];
r = 1;
[k,p,e] = lqr(a,b,q,r)

% these results show that a stable optimal proportional feedback sysem
% has been designed

% now let's display the results of the controlled system to an unit
% impule disturbance to each state variable

% define the closed loop response matrix

aa = a -b*k;

y = impulse(aa,bb,c,d,1,t);
u = -k*y';
plot(t,y,t,u')
title('Optimal State and Control Variables')
xlabel('Time')
ylabel('State Vector Components and Control Variables')
text(30,.6,'Conc')
text(10,.2,'Temp')
text(15,1,'Control')
print optialfb
pause

% implement optimal PI controller

aa = [a(1,1) a(1,2) 0;
   a(2,1) a(2,2) 0;
   1 0 0];
bb = [b;0];

% define the weighting matrices

qp = [1 0;
   0 0];
qi = 1;
q = [qp(1,1) qp(1,2) 0;
   qp(2,1) qp(2,2)0;
   0 0 qi];
r = 1;

% compute the steady state feedback gains

[k,p,e] = lqr(aa,bb,q,r)

% plot the response to a step input to each state

wi = [w(1,1) w(1,2) ;
   w(2,1) w(2,2) ;
   0 0 ];
wid = wi*[1;1];
c = eye(3);
d = zeros(3);
api = aa-bb*k;
xpi = step(api,wid,c,d,1,t);
plot (t,xpi(:,1:2))
title('PI Control State Response')
xlabel('Time')
ylabel('Normalized Composition and Temperature')
text(20,1,'Conc')
text(30,9,'Temp')
