# [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.
```