[SciPy-User] Scipy.optimize.fmin_slsqp problem

[josef.pktd at gmail.com]

On Wed, Jan 11, 2012 at 12:01 PM, Hendrik <hendrik.venter1 at gmail.com> wrote:
> 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
> 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...

I have no idea what might be wrong, and I don't have an atomparser
available, so I cannot run the script.

What I would do to debug and see whether it's a coding problem or a
problem with slsqp is to convert the atombalance constraint into a
quadratic penalty (with a penalization factor that can be increased)
and see whether an approximate "unconstraint" solution can be found.

Josef

>
> Can you please help??
>
> 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()
> 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,
> (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()
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user