Numerical solution of the P model
Benjamin D. Stocker
2019-01-11
library(rsofun)
library(dplyr)
lseq <- function(from=1, to=100000, length.out=6) {
# logarithmic spaced sequence
exp(seq(log(from), log(to), length.out = length.out))
}Get the analytical solution for one specific set of arguments.
## Run for different temperatures and plot
tc <- 20
ppfd <- 800 / 0.05
vpd <- 100
co2 <- 400
elv <- 0
fapar <- 1
out_analytical <- rpmodel( tc = tc, vpd = vpd, co2 = co2, elv = elv, kphio = 0.05, fapar = fapar, ppfd = ppfd, method_optci="prentice14", method_jmaxlim = "none", do_ftemp_kphio = FALSE )Is it internally consistent? \[ c_i = c_a - A / g_s \]
print( out_analytical$ci )## [1] 35.45723print( out_analytical$ca - (out_analytical$gpp / 12.0107) / out_analytical$gs )## [1] 35.45723Yes.
And what about … \[ A = V_{\text{cmax}} \frac{c_i-\Gamma^{\ast}}{c_i + K} \]
print( out_analytical$gpp / 12.0107 )## [1] 609.8003print( out_analytical$vcmax * (out_analytical$ci - out_analytical$gammastar) / (out_analytical$ci + out_analytical$kmm ))## [1] 609.8003print( out_analytical$gs * (out_analytical$ca - out_analytical$ci) )## [1] 609.8003Yes.
Maximise the following term numerically, subject to \(V_{\mathrm{cmax}}\) and \(g_s\): \[ A_n = A - aE -bV_{\text{cmax}} = max. \;\;\;\;\;\;\;\;\;\;\;\;(1) \]
Note, that this is (almost) equivalent to the optimality criterium in Prentice et al., 2014: \[ a \; \frac{\partial (E/A)}{\partial \chi} = -b \; \frac{\partial (V_{\mathrm{cmax}}/A)}{\partial \chi} \;\;\;\;\;\;\;\;\;\;\;\;(2) \] The twist is, that in above equation, \(V_{\mathrm{cmax}}\) cancels and we get “just” an optimal \(\chi\).
Instead of analytically taking the derivatives in Eq. 2 and solving for \(\chi\), we can apply a numerical optimisation algorithm to maximise the following expression (which is strictly equivalent to Eq. 2): \[ E/A + \beta V_{\mathrm{cmax}}/A = max. \] with \(\beta=b/a\). This can be written out as \[ \frac{1.6\;D}{c_a(1-\chi)} - \beta \frac{\chi c_a + K}{\chi c_a - \Gamma^{\ast}} = max. \] To find the maximum, let’s define the numerical optimization function:
calc_optimal_chi_num <- function( kmm, gammastar, ns_star, ca, vpd, beta ){
#-----------------------------------------------------------------------
# Input: - float, 'kmm' : Pa, Michaelis-Menten coeff.
# - float, 'ns_star' : (unitless) viscosity correction factor for water
# - float, 'vpd' : Pa, vapor pressure deficit
# Output: float, ratio of ci/ca (chi)
# Features: Returns an estimate of leaf internal to ambient CO2
# partial pressure following the "simple formulation".
# Depends: - kc
# - ns
# - vpd
#-----------------------------------------------------------------------
maximise_this <- function( chi, kmm, gammastar, ns_star, ca, vpd, beta ){
out <- 1.6 * ns_star * vpd / (ca * (1.0 - chi)) + beta * (chi * ca + kmm)/(chi * ca - gammastar)
return(out)
}
out_optim <- optimr::optimr(
par = 0.7,
lower = 0.01,
upper = 0.99,
fn = maximise_this,
kmm = kmm,
gammastar = gammastar,
ns_star = ns_star,
ca = ca,
vpd = vpd,
beta = beta,
method = "L-BFGS-B",
control = list( maxit = 100, maximize = TRUE )
)
return(out_optim$par)
}Now, let’s calculate the optimal \(\chi\) with the same parameters as used to run the P-model above.
chi_opt <- calc_optimal_chi_num( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146 )
print(chi_opt)## [1] 0.8748358As a check: is this the same as returned by the P-model?
print(out_analytical$chi)## [1] 0.8748392Practically, yes.
Let’s call this \(\chi^{\ast}\). For the P-model, it is then assumed that \(A=A_J\) with \[ A_J = \phi_0 \; I_{\mathrm{abs}}\;\frac{c_a \chi^{\ast} - \Gamma^{\ast}}{c_a \chi^{\ast} + 2\Gamma^{\ast}} \] A check again: is this the same as returned by the P-model?
print( 0.05 * ppfd * fapar * (out_analytical$ca * chi_opt - out_analytical$gammastar)/(out_analytical$ca * chi_opt + 2 * out_analytical$gammastar))## [1] 609.7997print( out_analytical$gpp / 12.0107 )## [1] 609.8003Yes, practically.
The analytically calculated \(A_n\) from the P-model (rpmodel() call above) is:
cost_scalar <- 0.001
cost_transp <- cost_scalar * 1.6 * out_analytical$ns_star * out_analytical$gs * vpd
cost_vcmax <- 146 * cost_scalar * out_analytical$vcmax
net_assim_analytical <- out_analytical$gpp / 12.0107 - cost_transp - cost_vcmax
print( net_assim_analytical )## [1] 362.0813Define the function which is to be maximised. This corresponds to Eq. 1.
calc_net_assim <- function( par, args, iabs, kphio, cost_scalar, maximize=FALSE, return_chi=FALSE, return_a_j_net=FALSE ){
kmm <- args[1]
gammastar <- args[2]
ns_star <- args[3]
ca <- args[4]
vpd <- args[5]
beta <- args[6]
vcmax <- par[1]
gs <- par[2]
## Get assimilation and ci, given Vcmax and gs
## by solving the equation system
## assim = vcmax * (ci - gammastar)/(ci + kmm)
## assim = gs * (ca - ci)
a_quad <- -1.0 * gs
b_quad <- gs * ca - gs * kmm - vcmax
c_quad <- gs * ca * kmm + vcmax * gammastar
root_ci <- try( polyroot( c(c_quad, b_quad, a_quad) ) )
if (class(root_ci)=="try-error"){
return( NA )
} else {
ci <- Re(root_ci)[which(Re(root_ci)>0)] # take positive root
# det <- sqrt( b_quad^2 - 4 * a_quad * c_quad )
# ci <- -b_quad - det / (2 * a_quad) # smaller root
## A_j
a_j <- iabs * kphio * (ci - gammastar) / (ci + 2 * gammastar)
## A_c
a_c <- vcmax * (ci - gammastar) / (ci + kmm)
# ## I don't understand if this is correct... It must be inconsistent with the Vcmax-related ci calculated above. But anyhow, this the logic of the P-model...
# assim <- a_j
## ... when it should actually be the minimum of A_j and A_c:
# assim <- min(a_j, a_c)
## All above is correct. That is, A and ci are correctly back-calculated from Vcmax and gs
## only cost ratio is defined. for this here we need absolute values. Set randomly
cost_transp <- cost_scalar * 1.6 * ns_star * gs * vpd
cost_vcmax <- cost_scalar * beta * vcmax
# ## Option A, Eq. 1. This does not give the same result as options B and C as the
# ## cost scalar (cost_scalar) does not cancel.
# net_assim <- assim - cost_transp - cost_vcmax
# ## Option B: This is equivalent to the P-model with its optimization of ci:ca.
# net_assim <- -(cost_transp + cost_vcmax) / a_c
## Option C: Equivalent (gives same result for ci:ca) as option B and the original P-model
## And: optimization target has ecological meaning.
## The trouble is that A_c is used to determine optimality here, but actual assimilation is
## later calculated following the light-limited function for A_j.
net_assim <- a_c - cost_transp - cost_vcmax
if (maximize) net_assim <- -net_assim
if (return_chi){
return(ci/ca)
} else if (return_a_j_net){
a_j_net <- a_j - cost_transp - cost_vcmax
return(a_j_net)
} else {
return(net_assim)
}
}
}Use the analytically optimal vcmax and ci as starting values and calculate net assimilation with these
## starting values for par = [vcmax, gs]. Use the ones returned by rpmodel()
par <- c( out_analytical$vcmax, out_analytical$gs )
args <- c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 )
out_numerical <- calc_net_assim( par, args, iabs=(fapar * ppfd), kphio = 0.05, cost_scalar=0.0003 )
print(out_numerical)## [1] 535.4846print( calc_net_assim( par, args, iabs=(fapar * ppfd), kphio = 0.05, cost_scalar=0.0003, return_chi = TRUE ) )## [1] 0.8748392print(out_analytical$chi)## [1] 0.8748392Ok, this is calculated correctly.
Optimise \(V_{\text{cmax}}\) and \(g_s\) simultaneously so that net assimilation is maximised (Implementing Eq. 1). This should return the same values for Vcmax and gs as the starting values (which were taken from the analytical solution), and the same \(\chi\) as the analytical solution of the P-model.
out_optim <- optimr::optimr(
par = c( out_analytical$vcmax, out_analytical$gs ),
lower = c( out_analytical$vcmax*0.0001, out_analytical$gs*0.001 ),
upper = c( out_analytical$vcmax*20, out_analytical$gs*30 ),
fn = calc_net_assim,
args = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ),
iabs = (fapar * ppfd),
kphio = 0.05,
cost_scalar = 0.0003,
method = "L-BFGS-B",
maximize = TRUE,
control = list( maxit = 1000 )
)
print( paste( "Vcmax, gs from gs/Vcmax optimization:", out_optim$par ) )## [1] "Vcmax, gs from gs/Vcmax optimization: 30969.0621926482"
## [2] "Vcmax, gs from gs/Vcmax optimization: 3606.31687358757"print( paste("A_n from gs/Vcmax optimization:", -out_optim$value) )## [1] "A_n from gs/Vcmax optimization: 11004.0401416738"print( paste( "Vcmax,gs from P-model:", c( out_analytical$vcmax, out_analytical$gs )) )## [1] "Vcmax,gs from P-model: 1548.45310963241"
## [2] "Vcmax,gs from P-model: 120.210562452919"print( paste( "chi from gs/Vcmax optimization:", calc_net_assim( par = out_optim$par, args = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ), iabs = (fapar * ppfd), kphio = 0.05, cost_scalar=0.0003, return_chi = TRUE ) ))## [1] "chi from gs/Vcmax optimization: 0.914101374016631"print(paste("... this should be identical to ", out_analytical$chi))## [1] "... this should be identical to 0.874839171881005"No. This doesn’t work.
The problem of all of this is that the cost factor that scales both the cost of \(V_{\text{cmax}}\) and \(g_s\) equally is missing from the optimization criterium when $ E/A + ; Vcmax/A = max.$ is used. Let’s term this cost scalar \(x\) and find it numerically.
setzero <- function( x_cost, chi_target, out_analytical, fapar, ppfd, vpd ){
out_optim <- optimr::optimr(
par = c( out_analytical$vcmax, out_analytical$gs ),
lower = c( out_analytical$vcmax*0.0001, out_analytical$gs*0.001 ),
upper = c( out_analytical$vcmax*20, out_analytical$gs*30 ),
fn = calc_net_assim,
args = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ),
iabs = (fapar * ppfd),
kphio = 0.05,
cost_scalar = x_cost,
maximize = TRUE,
method = "L-BFGS-B",
control = list( maxit = 100 )
)
args <- c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 )
chi_numerical <- calc_net_assim( out_optim$par, args, iabs=(fapar * ppfd), kphio = 0.05, cost_scalar=x_cost, return_chi = TRUE )
out <- chi_target - chi_numerical
return(abs(out))
}
sapply( lseq( from=0.001, to=0.003, length.out = 100), FUN = function(x) setzero( x, out_analytical$chi, out_analytical, fapar, ppfd, vpd ) ) %>%
plot( lseq( from=0.001, to=0.003, length.out = 100), ., type="l" )
abline(h=0, lty=2)
## Watch out: numerical noise is happening. For this example, chose starting point where visually determined in plot above
unitcost_root <- optimr::optimr(
par = 0.002,
# lower = 0.001,
# upper = 0.003,
fn = setzero,
chi_target = out_analytical$chi,
out_analytical = out_analytical,
fapar = fapar,
ppfd = ppfd,
vpd = vpd,
# method = "L-BFGS-B",
control = list( maxit = 100000 )
)
abline( v = unitcost_root$par, col="red" )
print( paste("My guess for the cost scalar is:", unitcost_root$par ) )## [1] "My guess for the cost scalar is: 0.00246166076660156"Now, try the same again and see if, with the determined fitted cost scalar, we get the same for the analytical P-model and the \(V_{\text{cmax}}\), \(g_s\) optimization: Optimise \(V_{\text{cmax}}\) and \(g_s\) simultaneously so that net assimilation is maximised (Implementing Eq. 1). This should now definitely return the same \(c_i\) as the analytical solution of the P-model.
out_optim <- optimr::optimr(
par = c( out_analytical$vcmax, out_analytical$gs ),
lower = c( out_analytical$vcmax*0.0001, out_analytical$gs*0.001 ),
upper = c( out_analytical$vcmax*5, out_analytical$gs*5 ),
fn = calc_net_assim,
args = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ),
iabs = (fapar * ppfd),
kphio = 0.05,
cost_scalar = unitcost_root$par,
method = "L-BFGS-B",
maximize = TRUE,
control = list( maxit = 1000 )
)
print( paste( "chi from gs/Vcmax optimization:", calc_net_assim( par = out_optim$par, args = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ), iabs = (fapar * ppfd), kphio = 0.05, cost_scalar=unitcost_root$par, return_chi = TRUE ) ))## [1] "chi from gs/Vcmax optimization: 0.874839171997369"print(paste("... this should be identical to ", out_analytical$chi))## [1] "... this should be identical to 0.874839171881005"Yes!!!
Plot a surface of \(A_n\) as a function of \(V_{\text{cmax}}\) and \(g_s\), with a red point where the numerical optimization found the maximum, and a green point, where the analytical solution lies.
library(dplyr)
len_vcmax <- 50
len_gs <- 50
vcmax_vec <- lseq(out_analytical$vcmax*0.0001, out_analytical$vcmax*5, length.out = len_vcmax)
gs_vec <- lseq(out_analytical$gs*0.001, out_analytical$gs*5, length.out = len_gs)
A_net <- expand.grid(vcmax_vec, gs_vec) %>%
setNames( c("vcmax", "gs") ) %>%
rowwise() %>%
do( net_assim = calc_net_assim(
c(.$vcmax, .$gs),
args=c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ),
iabs=(fapar * ppfd),
kphio = 0.05,
cost_scalar=unitcost_root$par,
maximize=FALSE
) ) %>%
tidyr::unnest(net_assim) %>%
as.matrix() %>%
array( dim=c(len_vcmax,len_gs) )
pmat <- persp( (vcmax_vec), (gs_vec), A_net,
phi = 10, # Azimut, 0 is along gs_vec axis
theta = 30,
col="lightblue"
)
points( trans3d( (out_optim$par[1]), (out_optim$par[2]), -out_optim$value, pmat ), pch=16, col="red", cex=1.8 )
points( trans3d( (out_analytical$vcmax), (out_analytical$gs), net_assim_analytical, pmat ), pch=16, col="green", cex=1 )
# numerical (of course, the value is not the same as analytical, because it's the optimization target and not A_j. Re-calculate A_j from optimized Vcmax and gs.):
print( c(out_optim$par[1], out_optim$par[2], -out_optim$value) )## [1] 1.548453e+03 1.202106e+02 1.798418e-04print( calc_net_assim(
par = c(out_optim$par[1], out_optim$par[2]),
args = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, 146.0 ),
iabs = (fapar * ppfd),
kphio = 0.05,
cost_scalar = unitcost_root$par,
return_a_j_net = TRUE
))## [1] 0.0001797671# analytical:
cost_transp <- unitcost_root$par * 1.6 * out_analytical$ns_star * out_analytical$gs * vpd
cost_vcmax <- 146 * cost_scalar * out_analytical$vcmax
net_assim_analytical <- out_analytical$gpp / 12.0107 - cost_transp - cost_vcmax
print( c(out_analytical$vcmax, out_analytical$gs, net_assim_analytical) )## [1] 1548.4531 120.2106 330.4439And in 2D. For fixed (optimal) \(g_s\).
library(dplyr)
## cost_scalar is searched by hand
sapply( as.list(vcmax_vec), FUN = function(x) calc_net_assim( c(x, out_optim$par[2]), args, iabs=(fapar * ppfd), kphio = 0.05, cost_scalar = unitcost_root$par ) ) %>% plot( vcmax_vec, ., type="l", xlab="Vcmax", ylab="net A", ylim=c(-500,500), xlim=c(0,5000))
abline(v=out_optim$par[1], lty=3)
points( (out_optim$par[1]), -out_optim$value, pch=16, col="red", cex=1.8 )
points( (out_optim$par[1]), net_assim_analytical, pch=16, col="green", cex=1 )
And for fixed (optimal) \(V_{\text{cmax}}\):
library(dplyr)
## cost_scalar is searched by hand
sapply( as.list(gs_vec), FUN = function(x) calc_net_assim( c(out_optim$par[1], x), args, iabs=(fapar * ppfd), kphio = 0.05, cost_scalar = unitcost_root$par ) ) %>% plot( gs_vec, ., type="l", xlab="gs", ylab="net A", ylim=c(-500,500), xlim=c(0,300))
abline( v=out_optim$par[2], lty=3 )
points( (out_optim$par[2]), -out_optim$value, pch=16, col="red", cex=1.8 )
points( (out_optim$par[2]), net_assim_analytical, pch=16, col="green", cex=1 )
The reason why the red and green dots are not exactly on top of each other is because the optimisation (red) returns \(A_{n}=A_c - \psi_{\text{transp}} - \psi_{\text{Vcmax}}\), while the analytical (green point) returns \(A_n = A_j - \psi_{\text{transp}} - \psi_{\text{Vcmax}}\).