Hello, I have a (soil respiration) model, with soil temperature and soil moisture as input variables, and a bunch of parameters, see the code below.
I am trying to figure out how to fit the model when one (or more) parameter is known in advance (so we don’t want to fit them).
Below is a MWE, with first what works (sometimes), and then what does not work.
I hope you guys have a solution for me!
function DAMM(x,p)
# Constants
R = 8.314472e-3 # Universal gas constant, kJ K⁻¹ mol⁻¹
O₂ₐ = 0.209 # Volume of O₂ in the air, L L⁻¹
Dₗᵢ = 3.17 # Diffusion coefficient of substrate in liquid phase, dimensionless
Dₒₐ = 1.67 # Diffusion coefficient of oxygen in air, dimensionless
pₛₓ = 0.024 # Fraction of soil carbon that is considered soluble
# Independent variables
Tₛ = x[:, 1] # Soil temperature, °C
θ = x[:, 2] # Soil moisture, m³ m⁻³
# Parameters
αₛₓ = 1e8*p[1] # Pre-exponential factor, mgC cm⁻³ h⁻¹
Eaₛₓ = 1e1*p[2] # Activation energy, kJ mol⁻¹
kMₛₓ = 1e-8*p[3] # Michaelis constant, gC cm⁻³
kMₒ₂ = 1e-3*p[4] # Michaelis constant for O₂, L L⁻¹
porosity = p[5] # 1 - soil buld density / soil particle density
Sxₜₒₜ = 1e-2*p[6] # Total soil carbon, gC cm⁻³
# DAMM model
Vmax = @. (αₛₓ * exp(-Eaₛₓ/(R * (273.15 + Tₛ)))) # Maximum potential rate of respiration
Sₓ = @. pₛₓ * Sxₜₒₜ * Dₗᵢ * θ^3 # All soluble substrate, gC cm⁻³
MMₛₓ = @. Sₓ / (kMₛₓ + Sₓ) # Availability of substrate factor, 0-1
O₂ = @. Dₒₐ * O₂ₐ * ((porosity - θ)^(4/3)) # Oxygen concentration
MMₒ₂ = @. O₂ / (kMₒ₂ + O₂) # Oxygen limitation factor, 0-1
Resp = @. Vmax * MMₛₓ * MMₒ₂ # Respiration, mg C cm⁻³ hr⁻¹
Respc = Resp .* 2314.8148 # Respiration, μmol CO₂ m⁻² s⁻¹
end
x_fake = [18.0 0.35; 22.0 0.22; 30.0 0.30; 35.0 0.39] # fake input data
p_test = [1, 6, 3.46, 2, 0.4, 1.25] # Parameters test
R_test = DAMM(x_fake, p_test) # function test. output = [3.0, 5.3, 9.4, 4.3]
R_fake = [1.2, 2.5, 5, 1.5] # fake output data. A bit lower than R_test
using LsqFit # Let's try to fit param to better match Resp
function fitDAMM(x, R) #
lb = [0.0, 0.0, 0.0, 0.0, 0.4, 0.0] # Porosity can't be below max(M)
ub = [Inf, Inf, Inf, Inf, 1.0, Inf] # Porosity can't be above 1
p_ini = [1, 6, 3.46, 2, 0.5, 1.25] # same as p
fit = curve_fit(DAMM, x, R, p_ini, lower=lb, upper=ub)
params = coef(fit)
return params
end
p_out = fitDAMM(x_fake, R_fake) # works
R_out = DAMM(x_fake, p_out) # good fit! very close to R_fake.
# attempt with forced p[5]
poro_val = 0.4 # Porosity parameter (p[5]) is known
function fitDAMM_f(x, R, poro_val)
lb = [0.0, 0.0, 0.0, 0.0, poro_val, 0.0] # force p[5] to 0.4
ub = [Inf, Inf, Inf, Inf, poro_val, Inf] # force p[5] to 0.4
p_ini = [1, 6, 3.46, 2, poro_val, 1.25] # same as p
fit = curve_fit(DAMM, x, R, p_ini, lower=lb, upper=ub)
params = coef(fit)
return params
end
p_out_f = fitDAMM_f(x_fake, R_fake, poro_val) # bug
I say the first “solution” works sometime, because sometime it doesn’t, as the algorithm of curve_fit still attempt parameters below the lower bound lb , which crashes the model (the model equations crash if porosity is lower than moisture).
What I need is either:
A way to force a parameter value
A way for the algorithm not to attempt parameter value below the lower bound
