[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

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

Below this message is my program... I'm struggling with this problem
for over 3 months now, asking for advice with no success!!

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...

Can you please help??


#!/usr/bin/env python

# Calculate mixture equilibriums using minimisation of the Gibbs

# 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()
ax = fig.add_subplot(111, projection='3d')

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,
                                            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,
        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
        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),
                                        acc = 1.0E-12)
        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)')

More information about the SciPy-User mailing list