% 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

eiga = 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);
sizev =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

% compute the eigenvalues and eigenvectors of a

[T,L] = eig(a)

% compute b tilda

btil = inv(T)*b;

% compute H tilda

htil = btil*inv(r)*btil';

% set the dynamic response of the conc with eigenvalue of 0.07


qtil = zeros(2);
qtil(1,1) = (alpha1^2-L(1,1)^2)/htil(1,1);
Tp = T';
q1 = inv(Tp)*qtil*inv(T);
[k1,p1] = lqr(a,b,q1,r)
a1 = a-b*k1;

% set the dynamic response of the temp with eigenvalue of .1

[T1,L1] = eig(a1)

% note that matlab always puts the smallest eigenvalue first which
% in this case is the temperature value

btil1 = inv(T1)*b;
htil1 = btil1*inv(r)*btil1';
qtil1 = zeros(2);
qtil1(1,1) =(.1^2-L1(1,1)^2)/htil1(1,1);
T1p= T1';
q2 = inv(T1p)*qtil1*inv(T1);
[k2,p2] = lqr(a1,b,q2,r)
a2 = a1-b*k2;
eiga2 = eig(a2)
y2 =  impulse(a2,bb,c,d,1,t);
u2 = -k2*y2';
plot(t,y2,t,u2')
title('Modal Control')
xlabel('T(time')
ylabel('State and Control Variables')
text(10,-10,'Control')
text(10,1.5,'Conc')
text(20,4,'Temp')
pause
plot(t,y(:,1),t,y2(:,1))
title('Concentration Response')
xlabel('time')
ylabel('Concentration')
text(20,.1,'Modal')
text(30,.5,'Regular')
pause
plot(t,y(:,2),t,y2(:,2))
title('Temperature Response')
xlabel('time')
ylabel('Temperature')
text(10,1,'Regular')
text(40,2,'Modal')


