# [SciPy-User] Scipy.optimize.fmin_slsqp problem

Hendrik hendrik.venter1@gmail....
Wed Jan 11 11:01:24 CST 2012

```Hey Folks!! I'm new to the Forum world. I looked for a topic on
optimization, but didn't see one, so sorry if I'm duplicating...

I run Python 2.6.6.1

I'm using scipy.optimize.fmin_slsqp  to solve a non-linear function
but it only returns the function value with the initial guess

Below this message is my program... I'm struggling with this problem

I want to minimize the function "changevars(self.gibbs)" within the
constraints "changevars(self.atombalance)".

The optimization program only returns the function value with the
initial guess values... And i don't know why.
When I run the program, I can see the minimum value graphically, but
the optimization can't calculate it...

Regards

*************************************************************************************************************************************
#!/usr/bin/env python

# Calculate mixture equilibriums using minimisation of the Gibbs
energy

# 201110 Originally by Hendrik Venter
# 201111 Significantly reworked by Carl Sandrock

from __future__ import division
import scipy.optimize
import scipy.linalg
import math
import numpy
import atomparser
import sys
import matplotlib.pyplot as pl
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

smallvalue = 1e-10
fig = pl.figure()

T = 298.15
R = 8.314
RT = R*T
mustplot = 1    # 1 for Plotting, 0 for not plotting

class compound:
""" Basic container for compound properties """
def __init__(self, name, DGf):
self.name = name
self.DGf = DGf
self.parsed = atomparser.parseformula(name)

class mixture:
""" Container for mixture properties """
def __init__(self, charge):
""" Initialise mixture - charge contains tuples of
(initialcharge, compound) """
self.N, self.compounds = zip(*charge)
self.N = numpy.array(self.N)
self.compoundnames = [c.name for c in self.compounds]
self.DGf = numpy.array([c.DGf for c in self.compounds])
self.elements = reduce(set.union,
(c.parsed.distinctelements()
for c in self.compounds))
self.S = numpy.array([c.parsed.counts(self.elements)
for c in self.compounds]).T
Srank = numpy.count_nonzero(numpy.linalg.svd(self.S,
compute_uv=False) > smallvalue)

# Calculating the Degrees of freedom
self.DOF = self.S.shape[1] - Srank
print 'Degrees of Freedom =', self.DOF

# Coefficients of the Gibbs function
Ncomps = len(self.compounds)
# self.A = numpy.tile(numpy.repeat([-1, 1], [1, Ncomps-1]),
[Ncomps, 1])
# number of atoms of each element
self.atoms = numpy.dot(self.S, self.N)

def gibbs(self, N=None):
""" Gibbs energy of mixture with N of each compound """
if N is None: N = self.N
# TODO: There is every chance that this function is not
correct. It needs to be checked.
#return sum(N*(self.DGf/(R*T) + numpy.log(numpy.dot(self.A,
N))))
logs = numpy.log(sum(N))
return sum(N*(self.DGf/RT + numpy.log(N) - logs))

def atombalance(self, N):

""" Atom balance with N of each compound """
#if N is None: N = self.N
return numpy.dot(self.S, N) - self.atoms

def conversion(self, conversionspec):
""" Calculate the composition given a certain conversion.
conversionspec is a list of 2-tuples containing a component
index or name and a conversion """
#TODO: A and B should only be calculated once
#TODO: This does not take into account any existing products
in the mixture
#TODO: This works only for conversion specified in terms of
reagents
if len(conversionspec) < self.DOF:
raise Exception("Not enough conversions specified.")
C = numpy.zeros([len(conversionspec), self.S.shape[1]])
Ic = C.copy()
for i, (j, c) in enumerate(conversionspec):
if type(j) is str:
j = self.compoundnames.index(j)
C[i, j] = 1-c
Ic[i, j] = 1
A = numpy.vstack([self.S, C])
B = numpy.vstack([self.S, Ic])
# A ni = B nf
nf, _, _, _ = scipy.linalg.lstsq(B, numpy.dot(A, self.N))
# assert residuals are neglegable
nf[nf<0] = 0
return nf

def equilibrium(self, initialconversion):
""" Return equilibrium composition as minimum of Gibbs Energy
"""
# guess initial conditions
N0 = self.conversion(initialconversion)
logN0 = numpy.log(N0)

# This decorator modifies a function to be defined in terms of
new variables
def changevars(f):
def newf(newX):
#         print 'newX', newX
#         print 'X', numpy.exp(newX)
r = f(numpy.exp(newX))
#          print 'f(X)', r
return r
return newf

# Find optimal point in terms of a change of variables
logN = scipy.optimize.fmin_slsqp(changevars(self.gibbs),
logN0,

f_eqcons=changevars(self.atombalance),
acc = 1.0E-12)
scipy.optimize.fmin
N = numpy.exp(logN)
print ''
print 'Calculated Optimmum Values'
print N
print ''
print 'Calculated Optimum Function value'
print m.gibbs(N)

return N

# Here for a particular mixture:
m = mixture([[1.,  compound('MgO',                      -568343.0)],
[0.5, compound('Al2O3',                    -1152420.0)],
[4.,  compound('H2O',                      -237141.0)],
[0,   compound('Mg(OH)2',                  -833644.0)],
[0,   compound('Al(OH)3',                  -1835750.0)]])

# find Gibbs energy surface for many conversion possibilities
Nsteps = 10
convrange = numpy.linspace(0.001, 0.999, Nsteps)
gibbssurface = numpy.zeros((Nsteps, Nsteps))
for iMgO, xMgO in enumerate(convrange):
for iAl2O3, xAl2O3 in enumerate(convrange):
N = m.conversion([('MgO', xMgO),
('Al2O3', xAl2O3/50)])
gibbssurface[iMgO, iAl2O3] = m.gibbs(N)

# Try to find equilibrium
# Initial conversion guess
print ''
print "***************Optimizing******************"
X0 = [0.5, 0.5]
conversionspec = zip(['MgO', 'Al2O3'], X0)
# Solve to equilibrium
N = m.equilibrium(conversionspec)

print ''
print '*************Final Values***************'
for n, init, name in zip(N, m.N, m.compoundnames):
print name, 'initial =', init, 'final =', n

# Plot results
if mustplot:
Xaxis, Yaxis = numpy.meshgrid(convrange, convrange)
Zaxis = gibbssurface
surf = ax.plot_surface(Xaxis, Yaxis, Zaxis, rstride=1, cstride=1,
cmap=cm.jet, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_xlabel('Conversion of MgO')
ax.set_ylabel('Conversion of Al2O3')
ax.set_zlabel('Gibbs free energy (G/RT)')
pl.show()

```