% 10.34 - Fall 2006
% Hw Set #7 - Problem #2
% Rob Ashcraft - Oct. 20, 2006

%This function solve the reactor optimization problem in Beers' text 5.B.4.
% This implementation uses a scaled variable LN(V_rxtr) as well as the
% other normal variables.  

function p2_main_logV
clear; clc;

% rate constant data at 298K and 310K for k1 --> k5
T_data = [298, 310];  % kelvin
k_T = [0.01, 0.02;  0.01, 0.02;  0.001, 0.005;  0.001, 0.005;  0.001, 0.005];  % L/mol-s

T_plot = linspace(273,335,20);
for i=1:length(k_T)
    EaR(i,1) = -log(k_T(i,1)/k_T(i,2)) / (1/T_data(1) - 1/T_data(2));  % solve for Ea/R for each reaction
    Af(i,1) = k_T(i,1)*exp(EaR(i)/T_data(1));  % solve for the A-factor using k(298K) 
    k_plot(:,i) = Af(i).*exp(-EaR(i)./T_plot);  % make a vector k(T) to plot
end

% plot the rate constants to make sure they fit the data
figure;
semilogy(T_data,k_T(1,:),'ok',T_plot,k_plot(:,1),'-',T_data,k_T(3,:),'*k',T_plot,k_plot(:,3),'--'); 
xlabel('Temperature (K)'); ylabel('Rate Constant (L/mole-s)');
legend('k_1 and k_2 data','k_1 and k_2 fit','k_3, k_4, k_5 data','k_3, k_4, k_5 fit','location','SouthEast');


% now deal with the reactor problem
V_flow = 1;  % L/s

% initial guesses for the optimization:
C_Ao_guess = 1;  % molar
C_Bo_guess = 1;  % molar
T_guess = 315;  % Kelvin
V_rxtr_guess = 5000;  % liters

% Choose the Volume variable to be LN(V)
var0 = [C_Ao_guess; C_Bo_guess; T_guess; log(V_rxtr_guess)];

% note that in the scaled case, you do not very high tolerances.  The
% tolerance in this case only serves to increase the precision of the
% solution (even default tolerance will get you very close)
options = optimset('Display','iter','TolFun',1e-10,'TolX',1e-10);

% Define the constraint that Ca + Cb < 2 M
A = [1, 1, 0, 0];
b = 2;

% Include the given upper and lower bounds on the system varaible
lb = [0; 0; 298; log(10)];
ub = [2; 2; 335; log(10000)];

% call fmincon to solve the problem
[var, fval] = fmincon(@obj_fun,var0,A,b,[],[],lb,ub,[],options,V_flow,Af,EaR);

disp(['Inlet [A] = ',num2str(var(1)),' M']);
disp(['Inlet [B] = ',num2str(var(2)),' M']);
disp(['Reactor Temperature = ',num2str(var(3)),' K']);
disp(['Reactor Volume = ',num2str(exp(var(4))),' L']);
disp(['Outlet [C] = ',num2str(-fval),' M']);

C_Ao = var(1);
C_Bo = var(2);
T = var(3);
V_rxtr = exp(var(4));

% Since we solved for [C]final by integrating to the SS condition, it is
% useful to make sure that the reactor actually reach steady-state.  We can
% perform a non-rigorous check by integrating based on the optimal
% conditions, and then visually inspect the plot.  This does not mean that
% all other parameters sets encountered during the optimization also
% reached SS, but it is a useful check none-the-less
t_span = [0, 1e5];
[time, C_all] = ode45(@rxtr_dCdt_ode,t_span,zeros(8,1),[],C_Ao,C_Bo,T,V_rxtr,V_flow,Af,EaR);

figure; plot(time,C_all(:,1),'--',time,C_all(:,2),':',time,C_all(:,3),'-',time,C_all(:,4),'-.');
xlabel('Time (s)'); ylabel('Concentration (M)'); title('Did Reactor Reach Steady-State?');
legend('[A]','[B]','[C]','[D]');

% Since we found out early that the V_rxtr term was causing problems, it is
% also useful to look at the cost function as a function of V_rxtr, and
% make sure our fmincon solution appears to be at the minimum.
test = logspace(1,4,50);
for i=1:length(test)
cost_i(i) = obj_fun([var(1:3);log(test(i))],V_flow,Af,EaR);
end

figure;
semilogx(test,cost_i,'-',V_rxtr,fval,'o')
xlabel('Reactor Volume (L)'); ylabel('Cost Function (-[C]_o_u_t  M)');
legend('Cost Function','fmincon Solution');
return;


% Other functions used in the program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the function called by fmincon with the cost function
function cost = obj_fun(var,V_flow,Af,EaR)

    C_Ao = var(1);
    C_Bo = var(2);
    T = var(3);
    LN_V_rxtr = var(4);
    
    V_rxtr = exp(var(4)); % since we are using LN(V) as the variable

    C_guess = [1; 1; 0; 0; 0; 0; 0; 0];

    % we solve for [C]final by integrating to a SS condition.  You must
    % make sure you integrate for a long enough time to achieve SS
    t_span = [0, 1e5];
    [time, C_all] = ode15s(@rxtr_dCdt_ode,t_span,C_guess,odeset('RelTol',1e-8),C_Ao,C_Bo,T,V_rxtr,V_flow,Af,EaR);

    % the cost function in the negative of [C]final
    cost = -(C_all(end,3));

return;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the reactor governing equations as a function.  It will return
% dC/dt for each of the species in the reactor
function dAlldt = rxtr_dCdt_ode(t,C_all,C_Ao,C_Bo,T,V_rxtr,V_flow,Af,EaR)
%  all of the C variable are concentrations in moles/L

    A = C_all(1);
    B = C_all(2);
    C = C_all(3);
    D = C_all(4);
    S1 = C_all(5);
    S2 = C_all(6);
    S3 = C_all(7);
    S4 = C_all(8);

    k = Af.*exp(-EaR./T);  % calculates the rate constants

    dAdt = -k(1)*A*B - k(3)*A*D - k(4)*A*B + V_flow/V_rxtr*(C_Ao - A);
    dBdt = -k(1)*A*B - k(2)*B*C - k(4)*A*B - k(5)*B*C + V_flow/V_rxtr*(C_Bo - B);
    dCdt = k(1)*A*B - k(2)*B*C - k(5)*B*C - V_flow/V_rxtr*C;
    dDdt = k(1)*A*B + k(2)*B*C - k(3)*A*D - V_flow/V_rxtr*D;
    dS1dt = k(2)*B*C - V_flow/V_rxtr*S1;
    dS2dt = k(3)*A*D - V_flow/V_rxtr*S2;
    dS3dt = k(4)*A*B - V_flow/V_rxtr*S3;
    dS4dt = k(5)*B*C - V_flow/V_rxtr*S4;

    dAlldt = [dAdt; dBdt; dCdt; dDdt; dS1dt; dS2dt; dS3dt; dS4dt];

return;