The Kitchin Research Group

The problem we solve here is that we have several compounds containing Ni and Al, and a bulk mixture of a particular composition of Ni and Al. We want to know which mixture of phases will minimize the total energy. The tricky part is that the optimization is constrained because the mixture of phases must have the overall stoichiometry we want. We formulate the problem like this.

Basically, we want to minimize the function \(E = \sum w_i E_i\), where \(w_i\) is the mass of phase \(i\), and \(E_i\) is the energy per unit mass of phase \(i\). There are some constraints to ensure conservation of mass. Let us consider the following compounds: Al, NiAl, Ni3Al, and Ni, and consider a case where the bulk composition of our alloy is 93.8% Ni and balance Al. We want to know which phases are present, and in what proportions. There are some subtleties in considering the formula and molecular weight of an alloy. We consider the formula with each species amount normalized so the fractions all add up to one. For example, Ni_3Al is represented as Ni_{0.75}Al_{0.25}, and the molecular weight is computed as 0.75*MW_{Ni} + 0.25*MW_{Al}.

We use scipy.optimize.fmin_slsqp to solve this problem, and define two equality constraint functions, and the bounds on each weight fraction.

Note: the energies in this example were computed by density functional theory at 0K.

import numpy as np
from scipy.optimize import fmin_slsqp

# these are atomic masses of each species
Ni = 58.693
Al = 26.982

COMPOSITIONS = ['Al', 'NiAl',              'Ni3Al',  'Ni']
MW = np.array(  [Al,  (Ni + Al)/2.0, (3*Ni + Al)/4.0, Ni])

xNi = np.array([0.0, 0.5, 0.75, 1.0])  # mole fraction of nickel in each compd
WNi = xNi*Ni / MW                      # weight fraction of Ni in each cmpd

ENERGIES = np.array([0.0, -0.7, -0.5, 0.0])

BNi = 0.938

def G(w):
    'function to minimize. w is a vector of weight fractions, ENERGIES is defined above.'
    return np.dot(w, ENERGIES)

def ec1(w):
    'conservation of Ni constraint'
    return BNi - np.dot(w, WNi)

def ec2(w):
    'weight fractions sum to one constraint'
    return 1 - np.sum(w)

w0 = np.array([0.0, 0.0, 0.5, 0.5]) # guess weight fractions

y = fmin_slsqp(G,   
               w0,
               eqcons=[ec1, ec2], 
               bounds=[(0,1)]*len(w0))

for ci, wi in zip(COMPOSITIONS, y):
    print '{0:8s} {1:+8.2%}'.format(ci, wi)

Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.233299644373
            Iterations: 2
            Function evaluations: 12
            Gradient evaluations: 2
Al         -0.00%
NiAl       +0.00%
Ni3Al     +46.66%
Ni        +53.34%

So, the sample will be about 47% by weight of Ni3Al, and 53% by weight of pure Ni.

It may be convenient to formulate this in terms of moles.

import numpy as np
from scipy.optimize import fmin_slsqp

COMPOSITIONS = ['Al', 'NiAl', 'Ni3Al',  'Ni']
xNi = np.array([0.0, 0.5, 0.75, 1.0])   # define this in mole fractions

ENERGIES = np.array([0.0, -0.7, -0.5, 0.0]) 

xNiB = 0.875  # bulk Ni composition

def G(n):
    'function to minimize'
    return np.dot(n, ENERGIES)

def ec1(n):
    'conservation of Ni'
    Ntot = np.sum(n)
    return (Ntot * xNiB) - np.dot(n,  xNi)

def ec2(n):
    'mole fractions sum to one'
    return 1 - np.sum(n)

n0 = np.array([0.0, 0.0, 0.45, 0.55]) # initial guess of mole fractions

y = fmin_slsqp(G,   
               n0,
               eqcons=[ec1, ec2], 
               bounds=[(0, 1)]*(len(n0)))

for ci, xi in zip(COMPOSITIONS, y):
    print '{0:8s} {1:+8.2%}'.format(ci, xi)

Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.25
            Iterations: 2
            Function evaluations: 12
            Gradient evaluations: 2
Al         +0.00%
NiAl       -0.00%
Ni3Al     +50.00%
Ni        +50.00%

This means we have a 1:1 molar ratio of Ni and Ni_{0.75}Al_{0.25}. That works out to the overall bulk composition in this particular problem.

Let us verify that these two approaches really lead to the same conclusions. On a weight basis we estimate 53.3%wt Ni and 46.7%wt Ni3Al, whereas we predict an equimolar mixture of the two phases. Below we compute the mole fraction of Ni in each case.

# these are atomic masses of each species
Ni = 58.693
Al = 26.982

# Molar case
# 1 mol Ni + 1 mol Ni_{0.75}Al_{0.25}
N1 = 1.0; N2 = 1.0
mol_Ni = 1.0 * N1 + 0.75 * N2
xNi = mol_Ni / (N1 + N2)
print xNi

# Mass case
M1 = 0.533; M2 = 0.467
MW1 = Ni; MW2 = 0.75*Ni + 0.25*Al

xNi2 = (1.0 * M1/MW1 + 0.75 * M2 / MW2) / (M1/MW1 + M2/MW2)
print xNi2

0.875
0.874192746385

You can see the overall mole fraction of Ni is practically the same in each case.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

adapated from Chemical Reactor analysis and design fundamentals, Rawlings and Ekerdt, appendix A.2.3.

Matlab post

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.

Reaction 1: \(I + B \rightleftharpoons P1\)

Reaction 2: \(I + B \rightleftharpoons P2\)

Here we define the Gibbs free energy of the mixture as a function of the reaction extents.

import numpy as np

def gibbs(E):
    'function defining Gibbs free energy as a function of reaction extents'
    e1 = E[0]
    e2 = E[1]
    # known equilibrium constants and initial amounts
    K1 = 108; K2 = 284; P = 2.5
    yI0 = 0.5; yB0 = 0.5; yP10 = 0.0; yP20 = 0.0
    # compute mole fractions
    d = 1 - e1 - e2
    yI = (yI0 - e1 - e2) / d
    yB = (yB0 - e1 - e2) / d
    yP1 = (yP10 + e1) / d
    yP2 = (yP20 + e2) / d
    G = (-(e1 * np.log(K1) + e2 * np.log(K2)) +
         d * np.log(P) + yI * d * np.log(yI) +
         yB * d * np.log(yB) + yP1 * d * np.log(yP1) + yP2 * d * np.log(yP2))
    return G

The equilibrium constants for these reactions are known, and we seek to find the equilibrium reaction extents so we can determine equilibrium compositions. The equilibrium reaction extents are those that minimize the Gibbs free energy. We have the following constraints, written in standard less than or equal to form:

\(-\epsilon_1 \le 0\)

\(-\epsilon_2 \le 0\)

\(\epsilon_1 + \epsilon_2 \le 0.5\)

In Matlab we express this in matrix form as Ax=b where

and

Unlike in Matlab, in python we construct the inequality constraints as functions that are greater than or equal to zero when the constraint is met.

def constraint1(E):
    e1 = E[0]
    return e1


def constraint2(E):
    e2 = E[1]
    return e2


def constraint3(E):
    e1 = E[0]
    e2 = E[1]
    return 0.5 - (e1 + e2)

Now, we minimize.

from scipy.optimize import fmin_slsqp

X0 = [0.2, 0.2]
X = fmin_slsqp(gibbs, X0, ieqcons=[constraint1, constraint2, constraint3],
               bounds=((0.001, 0.499),
                       (0.001, 0.499)))
print(X)

print(gibbs(X))

Optimization terminated successfully.    (Exit mode 0)
            Current function value: -2.55942394906
            Iterations: 7
            Function evaluations: 31
            Gradient evaluations: 7
[ 0.13336503  0.35066486]
-2.55942394906

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. We also set all elements where the sum of the two extents is greater than 0.5 to near zero, since those regions violate the constraints.

import numpy as np
import matplotlib.pyplot as plt

def gibbs(E):
    'function defining Gibbs free energy as a function of reaction extents'
    e1 = E[0]
    e2 = E[1]
    # known equilibrium constants and initial amounts
    K1 = 108; K2 = 284; P = 2.5;
    yI0 = 0.5; yB0 = 0.5; yP10 = 0.0; yP20 = 0.0;
    # compute mole fractions
    d = 1 - e1 - e2;
    yI = (yI0 - e1 - e2)/d;
    yB = (yB0 - e1 - e2)/d;
    yP1 = (yP10 + e1)/d;
    yP2 = (yP20 + e2)/d;
    G = (-(e1 * np.log(K1) + e2 * np.log(K2)) +
         d * np.log(P) + yI * d * np.log(yI) +
         yB * d * np.log(yB) + yP1 * d * np.log(yP1) + yP2 * d * np.log(yP2))
    return G


a = np.linspace(0.001, 0.5, 100)
E1, E2 = np.meshgrid(a,a)

sumE = E1 + E2
E1[sumE >= 0.5] = 0.00001
E2[sumE >= 0.5] = 0.00001

# now evaluate gibbs
G = np.zeros(E1.shape)
m,n = E1.shape

G = gibbs([E1, E2])

CS = plt.contour(E1, E2, G, levels=np.linspace(G.min(),G.max(),100))
plt.xlabel('$\epsilon_1$')
plt.ylabel('$\epsilon_2$')
plt.colorbar()

plt.plot([0.13336503],  [0.35066486], 'ro')

plt.savefig('images/gibbs-minimization-1.png')
plt.savefig('images/gibbs-minimization-1.svg')
plt.show()

You can see we found the minimum. We can compute the mole fractions pretty easily.

e1 = X[0];
e2 = X[1];

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

print('y_I = {0:1.3f} y_B = {1:1.3f} y_P1 = {2:1.3f} y_P2 = {3:1.3f}'.format(yI,yB,yP1,yP2))

y_I = 0.031 y_B = 0.031 y_P1 = 0.258 y_P2 = 0.680

Copyright (C) 2016 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 8.2.10

Matlab post (adapted from http://en.wikipedia.org/wiki/Lagrange_multipliers.)

Suppose we seek to maximize the function \(f(x,y)=x+y\) subject to the constraint that \(x^2 + y^2 = 1\). The function we seek to maximize is an unbounded plane, while the constraint is a unit circle. We want the maximum value of the circle, on the plane. We plot these two functions here.

import numpy as np

x = np.linspace(-1.5, 1.5)

[X, Y] = np.meshgrid(x, x)

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot_surface(X, Y, X + Y)

theta = np.linspace(0,2*np.pi);
R = 1.0
x1 = R * np.cos(theta)
y1 = R * np.sin(theta)

ax.plot(x1, y1, x1 + y1, 'r-')
plt.savefig('images/lagrange-1.png')

import numpy as np

def func(X):
    x = X[0]
    y = X[1]
    L = X[2] # this is the multiplier. lambda is a reserved keyword in python
    return x + y + L * (x**2 + y**2 - 1)

def dfunc(X):
    dLambda = np.zeros(len(X))
    h = 1e-3 # this is the step size used in the finite difference.
    for i in range(len(X)):
        dX = np.zeros(len(X))
        dX[i] = h
        dLambda[i] = (func(X+dX)-func(X-dX))/(2*h);
    return dLambda

from scipy.optimize import fsolve

# this is the max
X1 = fsolve(dfunc, [1, 1, 0])
print X1, func(X1)

# this is the min
X2 = fsolve(dfunc, [-1, -1, 0])
print X2, func(X2)

>>> ... >>> [ 0.70710678  0.70710678 -0.70710678] 1.41421356237
>>> ... >>> [-0.70710678 -0.70710678  0.70710678] -1.41421356237

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post

adapted from http://www.matrixlab-examples.com/linear-programming.html which solves this problem with fminsearch.

Let us suppose that a merry farmer has 75 roods (4 roods = 1 acre) on which to plant two crops: wheat and corn. To produce these crops, it costs the farmer (for seed, water, fertilizer, etc. ) $120 per rood for the wheat, and $210 per rood for the corn. The farmer has $15,000 available for expenses, but after the harvest the farmer must store the crops while awaiting favorable or good market conditions. The farmer has storage space for 4,000 bushels. Each rood yields an average of 110 bushels of wheat or 30 bushels of corn. If the net profit per bushel of wheat (after all the expenses) is $1.30 and for corn is $2.00, how should the merry farmer plant the 75 roods to maximize profit?

Let \(x\) be the number of roods of wheat planted, and \(y\) be the number of roods of corn planted. The profit function is: \( P = (110)($1.3)x + (30)($2)y = 143x + 60y \)

There are some constraint inequalities, specified by the limits on expenses, storage and roodage. They are:

\(\$120x + \$210y <= \$15000\) (The total amount spent cannot exceed the amount the farm has)

\(110x + 30y <= 4000\) (The amount generated should not exceed storage space.)

\(x + y <= 75\) (We cannot plant more space than we have.)

\(0 <= x and 0 <= y \) (all amounts of planted land must be positive.)

To solve this problem, we cast it as a linear programming problem, which minimizes a function f(X) subject to some constraints. We create a proxy function for the negative of profit, which we seek to minimize.

f = -(143*x + 60*y)

from scipy.optimize import fmin_cobyla

def objective(X):
    '''objective function to minimize. It is the negative of profit,
    which we seek to maximize.'''
    x, y = X
    return -(143*x + 60*y)

def c1(X):
    'Ensure 120x + 210y <= 15000'
    x,y = X
    return 15000 - 120 * x - 210*y

def c2(X):
    'ensure 110x + 30y <= 4000'
    x,y = X
    return 4000 - 110*x - 30 * y

def c3(X):
    'Ensure x + y is less than or equal to 75'
    x,y = X
    return 75 - x - y

def c4(X):
    'Ensure x >= 0'
    return X[0]

def c5(X):
    'Ensure y >= 0'
    return X[1]

X = fmin_cobyla(objective, x0=[20, 30], cons=[c1, c2, c3, c4, c5])

print 'We should plant {0:1.2f} roods of wheat.'.format(X[0])
print 'We should plant {0:1.2f} roods of corn'.format(X[1])
print 'The maximum profit we can earn is ${0:1.2f}.'.format(-objective(X))

   Normal return from subroutine COBYLA

   NFVALS =   40   F =-6.315625E+03    MAXCV = 4.547474E-13
   X = 2.187500E+01   5.312500E+01
We should plant 21.88 roods of wheat.
We should plant 53.12 roods of corn
The maximum profit we can earn is $6315.62.

This code is not exactly the same as the original post, but we get to the same answer. The linear programming capability in scipy is currently somewhat limited in 0.10. It is a little better in 0.11, but probably not as advanced as Matlab. There are some external libraries available:

1. http://abel.ee.ucla.edu/cvxopt/
2. http://openopt.org/LP

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source