Constrained minimization to find equilibrium compositions
adapated from Chemical Reactor analysis and design fundamentals, Rawlings and Ekerdt, appendix A.2.3.
The equilibrium composition of a reaction is the one that minimizes the total Gibbs free energy. The Gibbs free energy of a reacting ideal gas mixture depends on the mole fractions of each species, which are determined by the initial mole fractions of each species, the extent of reactions that convert each species, and the equilibrium constants.
The equilibrium constants for these reactions are known, and we seek to find the equilibrium reaction extents so we can determine equilibrium compositions.
clc, close all, clear all
we have the following constraints, written in standard less than or equal to form:
We express this in matrix form as Ax=b where and
A = [[-1 0];[0 -1];[1 1]]; b = [0;0;0.5]; % this is our initial guess X0 = [0.1 0.1];
X = fmincon(@gibbs,X0,A,b)
Warning: Trust-region-reflective algorithm does not solve this type of problem, using active-set algorithm. You could also try the interior-point or sqp algorithms: set the Algorithm option to 'interior-point' or 'sqp' and rerun. For more help, see Choosing the Algorithm in the documentation. Solver stopped prematurely. fmincon stopped because it exceeded the function evaluation limit, options.MaxFunEvals = 200 (the default value). X = NaN + NaNi NaN + NaNi
the error message suggests we try another solver. Here is how you do that.
options = optimset('Algorithm','interior-point'); X = fmincon(@gibbs,X0,A,b,,,,,,options) % note the options go at the end, so we have to put a lot of empty  as % placeholders for other inputs to fmincon. % we can see what the gibbs energy is at our solution gibbs(X)
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the default value of the function tolerance, and constraints were satisfied to within the default value of the constraint tolerance. X = 0.1334 0.3507 ans = -2.5594
these are simply mole balances to account for the consumption and production of each species
e1 = X(1); e2 = X(2); yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0; %initial mole fractions d = 1 - e1 - e2; yI = (yI0 - e1 - e2)/d; yB = (yB0 - e1 - e2)/d; yP1 = (yP10 + e1)/d; yP2 = (yP20 + e2)/d; sprintf('y_I = %1.3f y_B = %1.3f y_P1 = %1.3f y_P2 = %1.3f',yI,yB,yP1,yP2)
ans = y_I = 0.031 y_B = 0.031 y_P1 = 0.258 y_P2 = 0.680
One way we can verify our solution is to plot the gibbs function and see where the minimum is, and whether there is more than one minimum. We start by making grids over the range of 0 to 0.5. Note we actually start slightly above zero because at zero there are some numerical imaginary elements of the gibbs function or it is numerically not defined since there are logs of zero there.
[E1, E2] = meshgrid(linspace(0.001,0.5),linspace(0.001,0.5)); % we have to remember that when e1 + e2 > 0.5, there is no solution to our % gibbs functions, so we find the indices of the matrix where the sum is % greater than 0.5, then we set those indices to 0 so no calculation gets % done for them. sumE = E1 + E2; ind = sumE > 0.5; E1(ind) = 0; E2(ind) = 0; % we are going to use the arrayfun function to compute gibbs at each e1 and % e2, but our gibbs function takes a vector input not two scalar inputs. So % we create a function handle to call our gibbs function with two scalar % values as a vector. f = @(x,y) gibbs([x y]); % finally, we compute gibbs at each reaction extent G = arrayfun(f, E1, E2); % and make a contour plot at the 100 levels figure; hold all; contour(E1,E2,G,100) xlabel('\epsilon_1') ylabel('\epsilon_2') % let's add our solved solution to the contour plot plot([X(1)],[X(2)],'ko','markerfacecolor','k') colorbar % you can see from this contour plot that there is only one minimum.
ans = done
function retval = gibbs(X) % we do not derive the form of this function here. See chapter 3.5 in the % referenced text book. e1 = X(1); e2 = X(2); K1 = 108; K2 = 284; P = 2.5; yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0; d = 1 - e1 - e2; yI = (yI0 - e1 - e2)/d; yB = (yB0 - e1 - e2)/d; yP1 = (yP10 + e1)/d; yP2 = (yP20 + e2)/d; retval = -(e1*log(K1) + e2*log(K2)) + ... d*log(P) + yI*d*log(yI) + ... yB*d*log(yB) + yP1*d*log(yP1) + yP2*d*log(yP2); % categories: Optimization % tags: thermodynamics, reaction engineering