[SciPy-User] maximising a function

mdekauwe mdekauwe@gmail....
Mon Mar 11 01:21:21 CDT 2013


Apologies I just realised np.minimum did not do what I thought, which was
messing up the quadratic calculation. I think this seems like it is working.
I am still not sure if there is a better way to set up the constraints? But
I guess this seems OK.

thanks.




#!/usr/bin/env python

"""
Belinda's optimal N allocation model.

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.

"""
__author__ = "Martin De Kauwe"
__version__ = "1.0 (05.03.2013)"
__email__ = "mdekauwe@gmail.com"

import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import scipy.optimize
import warnings

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)
       
        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 J, 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 Wc, 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 Wj, 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. I am 
        only returning the smallest of the two roots 
        
        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 = 0.0
            root2 = 0.0
            warnings.warn("Imaginary root found")
        else:
            root1 = (-b - np.sqrt(d)) / (2.0 * a)
            #root2 = (-b + np.sqrt(d)) / (2.0 * a) # don't want to return
this
            
        return root1, 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, err = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star, Ca,
N_p)

        return self.sign * np.mean(An)
    
    def penalty(self, x0, par=None, Km=None, Rd=None, gamma_star=None,
Ca=None, 
                N_p=None):
    
        
        An, error = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star,
Ca, N_p)
        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)
        
        if (N_s + N_e + N_r + N_c + N_other) > N_p:
           
            return np.nan
        
        An, error = self.calc_photosynthesis(x0, par, Km, Rd, gamma_star,
Ca, N_p)
        if error:
            
            return np.nan
        else:
            return self.sign * np.mean(An)
        
        
        #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.calc_photosynthesis
    #F = P.penalty # short name for func
    F = P.optimise_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(F, x0, method="COBYLA",
constraints=cons, 
                                     args=args)
    if result["success"]:
        print result
        (N_c, N_e, N_r, N_s, N_other) = P.N_pool_store
        fitted_x0 = np.array([N_c, N_e, N_r, N_other])
    
        # bit of a hack as there can be small rounding errors,
        # not sure of a better solution?
        diff = np.sum(P.N_pool_store) - N_p
        print "Diff", diff
        print fitted_x0
    
        An, err = func(fitted_x0, par, Km, Rd, gamma_star, Ca, N_p)
        print err
        print np.mean(An)
        plt.plot(par, An, label="Cobyla")
        plt.legend(loc="best", numpoints=1)
        plt.show()


if __name__ == "__main__":
    
    
    optimise_test()
    
    



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


More information about the SciPy-User mailing list