Algorithm in gem

The Algorithm used in gem for calculating the minimum point of G,
Gibbs energy of the system.

A unique algorithm has been developed for gem in order to solve chemical equilibrium problems in general.
This algorithm is mainly based on the gradient projection method with aids of some concepts in linearprogramming. This algorithm neither needs any estimated initial values other than derived values fromthe input quantities for reactant mole numbers, temperature and pressure (or volume), nor passes
unreasonable states in which mass balance equations will not be satisfied. Therefore, this leads to excellent
stability in finding the solution.
1. The problem
Let the system consist of n compounds and m elements. Then, the problem can be described as:
Minimize the Gibbs energy of the systemnnn
G= ∑(gi/RT)Xi + ∑Xi ln (Xi P/ ∑Xj )
i=1
i=k
j=K

(1)

where xi: the quantity of the i-th compound (mol)the 1st to (k-1)th compounds are condensed,
and the rest of compounds are gaseous.
P: Pressure (atm)
T: Temperature(K)
R: Gas constantgi:Standard Gibbs formation energy of the i-th compound at temperature T
ln : natural logarithm
subject toa11x1+a12x2+…+a1nxn=b1
a21x1+a22x2+…+a2nxn=b2…

(2)…

am1x1+am2x2+…+amnxn=bm

andxi ≥0 i=1, ..., n(3)whereaji: stoichiometry coefficient of the i-th compound for the j-th element

We will use vector notation, and use primes (') to indicate the transposes of vectors and matrices.
Equations (2) can be written as

AX'=B
where
X = [x1,x2, …, xn]
A=A1
A2 =.

Am(2)

a11, a12, …, a1na21, a22, …, a2n…

…

am1, am2, …, amn

B is constant vector which is given by initial quantities in the system.

The steepest descent vector of G is-gradG = − [∂G/∂xi]
We will use this vector and find the minimum point of G under the constraints of(2) and (3). Here, we assume the G is convex in all area defined by (2) and (3).
2. The gradient projection method
If we project the vector -gradG into the space that satisfies (2), the moving point along that direction isalways satisfying the mass-balance law. Let {Ai} be the space spanned by the vectors A1, ..., Am, and S
be the space which satisfies (2), then {Ai} and S are the orthogonal complement of each other. In otherwords, if we subtract the components which belong to {Ai} from -gradG, the resultant vector belongs to S.
Hence we first transform the vector series Ai into the orthogonal vectors Li by the Gram-Schmidtorthogonalization.
Then, we obtain the vector Q as a result of the subtraction of the Li components from -gradG.
3. One stage of calculation for convergence
The initial point of calculation x0 is given as the quantities of input materials, and we derive the vector Q
as above mentioned. The next problem is how far the point should move along that direction. We canderivec = QQ '/QHQ ' where H=[∂2G/∂xi∂xj](4)and adopt X= X0+cQ as the next point based on the ordinary principle of the steepest descent method.
This point is, however, either over-going or short-going in general. In order to examine which situationis the case, we make the inner product of -gradG at X and Q. If the inner product is negative, the point isover-going, so that we will adopt some inner point between X and X0. If positive, we will go farther.
In the present algorithm, an inter(extra)polating coefficient for each case is given as functions of c, Q andr, which will be described later. By successively applying this method, we finally obtain the point on thisdirection as the next point where the norm of the inner product is less than an appropriate thresholdvalue.
4. Linear programming
In principle, the gradient projection method is the core of this algorithm. However, we cannot always
move along the direction of Q because of non-negative constraints (3). Moreover, if we handle all variables
equally at a calculation stage, there would be some problems arisen from the fact that each variable has
the value of quite different order.
In order to avoid this problem, our program has employed some concepts utilized in the linearprogramming. In the linear programming, the movement of point on the convex polyhedron is restricted
such as the point should go to the adjacent point only at a calculation stage. That is, the space in which
the point is movable in a calculation stage is restricted within one dimensional space including the point.
In our program, the space is being restricted in one dimensional subspace of S. Actually, we make one
calculation stage under the condition that (n-m-1) variables out of n variables are regarded to be fixed asconstant.

In the next stage, the variable which has the smallest value in the resultant variables is fixed and we take
into account the oldest variable which has been fixed earlier.
Initially, the bigger variable will be served earlier.
For the non-negative constraints, let

r=min{xi/|qi| : qi<0} where [qi]=Q(5)

and when r is smaller than c in (4), we use this r value instead of c in order to avoid negative variables.
In the case of r=0, no calculation is made at the stage and the variable xi which is 0 and qi <0 will be fixedin the next stage.
This calculation is iterated until the difference of each variable over n-m successive stages becomes lessthan 1/10000 of the value.
To confirm that the point satisfies the condition of minimum, the program outputs the norm of theprojected gradG (for i : xi >0 ) on S. The smaller value indicates that the solution can be regarded as moreprecise.