% implementation of a Gibbs reactor (optimization problem: minimzing Gibbs free enthalpy) without
% considering chemical reactions
% system of reaction basis 1 solved with fmincon for non-linear function
% with eqality and inequality constraint
% components: H2O, CO2, CO and H2, p=1bar and T=1300 K
clear all
clc

% reaction basis = number of species - number of elements
R = 8.3145; % in J/mol/K
T = 1300;   % in K
p = 1e5;    % in Pa
p0 = 1e5;   % in Pa

% at the moment works only for T=1300 K
[CO2, H2O, CO, H2] = read_component_data(T);

% molar mass of "elements" in kg/kmol
C.M = 12;
O.M = 16;
O2.M = 32;

CO2.m = CO2.n0 * CO2.M; 
CO.m = CO.n0 * CO.M;
H2O.m = H2O.n0 * H2O.M; 
H2.m = H2.n0 * H2.M; 

% initial values in kg (for m0) and kg/kmol (for M0_tot)
CO2.m0 = 90;
CO.m0 = 23;
H2.m0 = 35;
H2O.m0 = 15;
M0_tot = 30;

% overall mass
m = CO2.n0 * CO2.M + H2.n0 * H2.M;

% array m_ of initial values
m_ = [CO2.m0, CO.m0, H2O.m0, H2.m0, M0_tot];
[c, ceq] = constraints(m_, C, O, O2, CO, CO2, H2O, H2, m);

% objective function G_min in function
G_min = objective(m_, CO, CO2, H2O, H2, m, T, p, p0, R);

A = []; b = []; Aeq = []; beq = []; lb = []; ub = [];

% solver fmincon for non-linear function with eqality and inequality
% constraints 
nonlincon = @(m_) constraints(m_, C, O, O2, CO, CO2, H2O, H2, m);
[x, fval, ef, output, lambda] = fmincon(@(m_)G_min, m_, A, b, Aeq, beq, lb, ub, nonlincon);

m_CO2_sol = x(1);
m_CO_sol = x(2);
m_H2O_sol = x(3);
m_H2_sol = x(4);
M_tot_sol = x(5);

disp(['solver message: ', num2str(output.message)])
disp(['No of iterations: ', num2str(output.iterations)])
disp(['Solution for m(CO2): ', num2str(m_CO2_sol), ' kg, solution for m(CO): ', num2str(m_CO_sol), ' kg, solution for m(H2O): ', num2str(m_H2O_sol), ' kg, solution for m_H2: ', num2str(m_H2_sol), ' kg'])