[SciPy-User] maximising a function

mdekauwe mdekauwe@gmail....
Mon Mar 11 00:30:58 CDT 2013


Hi,

I am trying to use the scipy library to optimise 4 parameters so that they
return the maximum for a given function. To do this I am using the COBYLA
function and setting a few constraints, essentially that these parameter
values are positive and that when summed they are less than a given value
(N_p). Can I check what I have is the correct implementation (see below)? 

The second issue is that part of the function depends on the solution to a
quadratic. If I don't resolve positive roots I want to be able to return
some value to let the minimiser know that these are bad values. I can't work
out how to do this? I tried returning a very large value but this didn't see
to work, so then I tried to use NaNs. What is the correct way to do this?

thanks,

Code (this will run).

#!/usr/bin/env python

"""
Optimise the distribution of nitrogen within the photosynthetic system to 
maximise photosynthesis for a given PAR and CO2 concentration.

Reference:
=========
* Medlyn (1996) The Optimal Allocation of Nitrogen within the C3
Photsynthetic 
  System at Elevated CO2. Australian Journal of Plant Physiology, 23,
593-603.

"""

import sys
import numpy as np
import scipy.optimize

class PhotosynthesisModel(object):
    """
    Photosynthesis model based on Evans (1989)
    
    References:
    -----------
    * Evans, 1989
    """
    def __init__(self, alpha=0.425, beta=0.7, g_m=0.4, theta=0.7,
K_cat=24.0,
                 K_s=1.25E-04, MAXIMISE=False):
        """
        Parameters
        ----------
        alpha : float
            quantum yield of electron transport (mol mol-1) 
        beta : float
            constant of proportionality between atmospheric and
intercellular 
            CO2 concentrations (-)
        g_m : float
            Conductance to CO2 transfer between intercellular spaces and
sites 
            of carboxylation (mol m-2 s-1 bar-1) 
        theta : float
            curvature of the light response (-)    
        K_cat : float
            specific activitt of Rubisco (mol CO2 mol-1 Rubisco s-1)
        K_s : float
            constant of proportionality between N_s and Jmax (g N m2 s
umol-1)
        MAXIMISE : logical
            if we want to minimise a function we need to return f(x), if we
            want to maximise a function we return -f(x)
        """
        self.alpha = alpha
        self.beta = beta
        self.g_m = g_m
        self.theta = theta
        self.K_cat = K_cat
        self.K_s = K_s
        self.sign = 1.0
        if MAXIMISE:
            self.sign = -1.0
        else:
            self.sign = 1.0
        
    def calc_photosynthesis(self, N_pools=None, par=None, Km=None, Rd=None, 
                            gamma_star=None, Ca=None, N_p=None,
error=False):
        """
        Leaf temperature is assumed to be constant = 25 deg C.
        
        
        Parameters
        ----------
        N_pools : list of floats
            list containing the 5 N pools: 
            N_c - amount of N in chlorophyll (mol m-2)
            N_s - amount of N in soluble protein other than Rubisco (mol
m-2)
            N_e - amount of N in electron transport components (mol m-2)
            N_r - amount of N in Rubisco (mol m-2)
            N_other - amount of leaf N not involved in photosynthesis (mol
m-2)
            
        par : float
            incident PAR (umol m-2 s-1)
        Km : float
            effective Michalis-Menten coefficent of Rubisco (umol mol-1)   
        Rd : float
            rate of dark respiration (umol m2 s-1)
        gamma_star : float
            leaf CO2 compensation point (umol mol-1)
        Ca : float
            atmospheric CO2 concentration (umol mol-1)
        N_p : float
            total amount of leaf N to be distrubuted (mol m-2)
        
        Returns:
        --------
        An : float
            Net leaf assimilation rate [umol m-2 s-1] 
          
        """
        # unpack...we need to pass as a list for the optimisation step
        (N_e, N_c, N_r, N_other) = N_pools
        
        # Leaf absorptance depends on the chlorophyll protein complexes
(N_c)
        absorptance = self.calc_absorptance(N_c)
        
        # Max rate of electron transport (umol m-2 s-1) 
        Jmax = self.calc_jmax(N_e, N_c)
        
        # Max rate of carboxylation velocity (umol m-2 s-1) 
        Vcmax = self.calc_vcmax(N_r)
        N_s = self.calc_n_alloc_soluble_protein(Jmax)
        
        # store values so I can get them all later
        self.N_pool_store = np.array([N_e, N_c, N_r, N_s, N_other])
        
        # rate of electron transport, a saturating function of absorbed PAR
        J, error = self.quadratic(a=self.theta, 
                                  b=-(self.alpha * absorptance * par +
Jmax), 
                                  c=self.alpha * absorptance * par * Jmax)
        if error:
            return None, error
        
        
        # CO2 concentration in the intercellular air spaces 
        Ci = self.beta * Ca
        
        # Photosynthesis when Rubisco is limiting
        a = 1. / self.g_m
        b = (Rd - Vcmax) / self.g_m - Ci - Km
        c = Vcmax * (Ci - gamma_star) - Rd * (Ci + Km)
        Wc, error = self.quadratic(a=a, b=b, c=c)
        if error:
            return None, error

        # Photosynthesis when electron transport is limiting
        a = -4. / self.g_m
        b = (Rd - J / 4.0) / self.g_m - Ci - 2.0 * gamma_star
        c = J / 4.0 * (Ci - gamma_star) - Rd * (Ci + 2.0 * gamma_star)
        Wj, error = self.quadratic(a=a, b=b, c=c)
        if error:
            return None, error
        
        Aleaf = np.minimum(Wc, Wj)
        
        Cc = Ci - (Aleaf / self.g_m);
        An = (1.0 - (gamma_star / Cc)) * Aleaf - Rd
        
        
        return An, error
 
        
    def assim(self, Cc, a1, a2):
        """Photosynthesis with the limitation defined by the variables
passed 
        as a1 and a2, i.e. if we are calculating vcmax or jmax limited.
        
        Parameters:
        ----------
        Cc : float
            CO2 concentration at the site of of carboxylation
        gamma_star : float
            CO2 compensation point in the abscence of mitochondrial
respiration
        a1 : float
            variable depends on whether the calculation is light or rubisco 
            limited.
        a2 : float
            variable depends on whether the calculation is light or rubisco 
            limited.

        Returns:
        -------
        assimilation_rate : float
            assimilation rate assuming either light or rubisco limitation.
        """
        return a1 * (Cc / (Cc + a2))
    
    def quadratic(self, a=None, b=None, c=None, error=False):
        """ minimilist quadratic solution as root for J solution should
always
        be positive, so I have excluded other quadratic solution steps. If 
        roots for J are negative or equal to zero then the data is bogus 
        
        Parameters:
        ----------
        a : float
            co-efficient 
        b : float
            co-efficient
        c : float
            co-efficient
        
        Returns:
        -------
        val : float
            positive root
        """
        d = b**2 - 4.0 * a * c # discriminant
        if np.any(d <= 0.0) or np.any(np.isnan(d)):
            error = True
            root1 = -999.9
            root2 = -999.9 
        else:
            root1 = (-b + np.sqrt(d)) / (2.0 * a)
            root2 = (-b - np.sqrt(d)) / (2.0 * a)
            
        return np.minimum(root1, root2), error
        
    
    def calc_jmax(self, N_e, N_c, a_j=15870.0, b_j=2775.0):
        """ Evans (1989( found a linear reln between the rate of electron
        transport and the total amount of nitrogen in the thylakoids per
unit
        chlorophyll 
        
        Parameters:
        ----------
        N_e : float
            amount of N in electron transport components (mol m-2)
        N_c : float
            amount of N in chlorophyll (mol m-2)
        a_j : float
            umol (mol N)-1 s-1
        b_j : float
            umol (mol N)-1 s-1
            
        Returns:
        --------
        Jmax : float
            Max rate of electron transport (umol m-2 s-1)
        """
        Jmax = a_j * N_e + b_j * N_c
        
        return Jmax
    
    def calc_vcmax(self, N_r):
        """ 
        Calculate Rubisco activity 
        
        Parameters:
        ----------
        N_r : float
            amount of N in Rubisco (mol m-2)
        
        Returns:
        --------
        Vcmax : float
            Max rate of carboxylation velocity (umol m-2 s-1) 
        """
        conv = 7000.0 / 44.0 # mol N to mol Rubiso
        
        return self.K_cat * conv * N_r
        
    
    def calc_absorptance(self, N_c):
        """ Calculate absorptance in the leaf based on N in the chlorophyll
        protein complexes
        
        Parameters:
        ----------
        N_c : float
            amount of N in chlorophyll (mol m-2)
            
        Returns:
        --------
        absorptance : float
            Leaf absorptance (-)
        """
        return (25.0 * N_c) / (25.0 * N_c + 0.076)
    
    def calc_n_alloc_soluble_protein(self, Jmax):
        """ Calculate N allocated to soluble protein
        
        Parameters:
        --------
        Jmax : float
            Max rate of electron transport (umol m-2 s-1)
        
        Returns:
        --------
        N_s : float
            amount of N in soluble protein other than Rubisco (mol m-2)
        """
        return self.K_s * Jmax 
  
    def optimise_func(self, x0, par=None, Km=None, Rd=None, 
                   gamma_star=None, Ca=None, N_p=None):
        """ Minimisation and maximisation are equivalent, to get the global
        maximium we just need to return -f(x) instead of f(x) """
        
        An, error = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star,
Ca, N_p)
        
        # I want only positive values for xopt
        if np.any(x0<0.0):
            return np.nan
        
        (N_e, N_c, N_r, N_other) = x0
        
        # Leaf absorptance depends on the chlorophyll protein complexes
(N_c)
        absorptance = self.calc_absorptance(N_c)
        
        # Max rate of electron transport (umol m-2 s-1) 
        Jmax = self.calc_jmax(N_e, N_c)
        
        # Max rate of carboxylation velocity (umol m-2 s-1) 
        Vcmax = self.calc_vcmax(N_r)
        N_s = self.calc_n_alloc_soluble_protein(Jmax)
        
        # xopt can't be > N_p
        if (N_s + N_e + N_r + N_c + N_other) > N_p:
           return np.nan
        
        An, err = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star, Ca,
N_p)
        if err:
            return np.nan
        else:
            return self.sign * np.mean(An)
       
    
    
    def constraint1(self, x, *args):
        """ N_s + N_e + N_r + N_c < N_p 
        
        returns a positive number if within bound and 0.0 it is exactly on
the 
        edge of the bound
        """
        
        (par, Km, Rd, gamma_star, Ca, N_p) = args
        
        # unpack...we need to pass as a list for the optimisation step
        (N_e, N_c, N_r, N_other) = x
        
        # Leaf absorptance depends on the chlorophyll protein complexes
(N_c)
        absorptance = self.calc_absorptance(N_c)
        
        # Max rate of electron transport (umol m-2 s-1) 
        Jmax = self.calc_jmax(N_e, N_c)
        
        # Max rate of carboxylation velocity (umol m-2 s-1) 
        Vcmax = self.calc_vcmax(N_r)
        N_s = self.calc_n_alloc_soluble_protein(Jmax)
        
        return N_p - (N_s + N_e + N_r + N_c + N_other)

    
    def constraint2(self, x, *args):
        """ N_e > 0.0 
        
        returns a positive number if within bound and 0.0 it is exactly on
the 
        edge of the bound
        """
        return x[0]
        
    def constraint3(self, x, *args):
        """ N_c > 0.0 
        
        returns a positive number if within bound and 0.0 it is exactly on
the 
        edge of the bound
        """
        return x[1]
    
    def constraint4(self, x, *args):
        """ N_r > 0.0 
        
        returns a positive number if within bound and 0.0 it is exactly on
the 
        edge of the bound
        """
        return x[2]
    
    def constraint5(self, x, *args):
        """ N_other > 0.0 
        
        returns a positive number if within bound and 0.0 it is exactly on
the 
        edge of the bound
        """
        return x[3]

 
    

def optimise_test():
    # default values
    Km = 544.
    Rd = 0.5
    gamma_star = 38.6
    
    par = np.arange(10.0, 1500.0, 20.0)
    Ca = np.ones(len(par)) * 350.0
    
    P = PhotosynthesisModel(MAXIMISE=True)
    func = P.optimise_func # short name for func
    
    # initial guess
    # Testing, just equally divide up the total N available
    N_p = 0.09
    x0 = np.array([N_p/5.0, N_p/5.0, N_p/5.0, N_p/5.0])
    args = (par, Km, Rd, gamma_star, Ca, N_p)

    cons = ({'type': 'ineq', 'fun': P.constraint1, 'args': args},
            {'type': 'ineq', 'fun': P.constraint2, 'args': args},
            {'type': 'ineq', 'fun': P.constraint3, 'args': args},
            {'type': 'ineq', 'fun': P.constraint4, 'args': args},
            {'type': 'ineq', 'fun': P.constraint5, 'args': args})
    result = scipy.optimize.minimize(func, x0, method="COBYLA", 
                                     constraints=cons, args=args)
    if result["success"]:
        print result.x
      

if __name__ == "__main__":
    
    optimise_test()



--
View this message in context: http://scipy-user.10969.n7.nabble.com/maximising-a-function-tp17980.html
Sent from the Scipy-User mailing list archive at Nabble.com.


More information about the SciPy-User mailing list