[SciPy-User] maximising a function
mdekauwe
mdekauwe at gmail.com
Mon Mar 11 02:21:21 EDT 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 at 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