diff --git a/exercises/solved_notebooks/P7_sens_uncert/sens_fermenter_monod_sol.jl b/exercises/solved_notebooks/P7_sens_uncert/sens_fermenter_monod_sol.jl index 84d704a..0b56a62 100644 --- a/exercises/solved_notebooks/P7_sens_uncert/sens_fermenter_monod_sol.jl +++ b/exercises/solved_notebooks/P7_sens_uncert/sens_fermenter_monod_sol.jl @@ -1,48 +1,44 @@ ### A Pluto.jl notebook ### -# v0.19.46 +# v0.20.21 using Markdown using InteractiveUtils # ╔═╡ 3ef93246-657d-4e77-9bf0-8380c64bfcfd -begin - # add this cell if you want the notebook to use the environment from where the Pluto server is launched - using Pkg - Pkg.activate("..") -end +using Pkg; Pkg.activate("..") # ╔═╡ 55cdebd2-0881-11ef-2722-91de1447877a -using Markdown - -# ╔═╡ 03ae0690-06a0-4276-9f00-d07b206fe124 -using InteractiveUtils +using Markdown, InteractiveUtils # ╔═╡ a355b0ba-baaf-49f4-a5dc-965364a884f0 -using Catalyst +using Catalyst, OrdinaryDiffEq # ╔═╡ 00fd6d49-f561-42e9-9413-d33af92f83dc -using OrdinaryDiffEq, StatsPlots +using StatsPlots, PlutoUI; TableOfContents() # ╔═╡ 7ae714c4-d25d-4f9f-ab3d-cc067db9c156 using ForwardDiff +# ╔═╡ dfce0717-e33e-4b6b-bb5c-d8754a74aba4 +using Turing + # ╔═╡ 31d294d1-3a1f-41db-abff-54f2a67c7ed9 md""" -### Exercise: Fermenter - Monod kinetics - Sensitivity analysis +# Exercise: Fermenter - Monod kinetics - Sensitivity analysis """ # ╔═╡ 5ffe7dcb-620d-4f22-95fe-2f77cda6fbe7 md""" -In one of the previous practica we were introduced to a fermenter in which biomass $X$ [$g/L$] grows by breaking down substrate $S$ [$g/L$]. The reactor is fed with a inlet flow rate $Q_{in}$ [$L/h$], which consist of a (manipulable) input concentration of substrate $S_{in}$ [$g/L$]. This process was modelled using Monod kinetics, resulting in the model below: +In one of the previous practica, we were introduced to a fermenter in which biomass $X$ [$g/L$] grows by breaking down substrate $S$ [$g/L$]. The reactor is fed with an inlet flow rate $Q_{in}$ [$L/h$], which consists of a (manipulable) input concentration of substrate $S_{in}$ [$g/L$]. This process was modelled using Monod kinetics, resulting in the model below: $$\begin{eqnarray*} -%S \xrightarrow[\quad\quad]{\beta} Y \, X -S \xrightarrow[\quad\quad]{r} Y \, X \quad\quad\quad\quad r = \mu \, X +S + X \xrightarrow[\quad\quad]{k} (1 + Y) \, X \quad\quad\quad\quad \textrm{with} \quad k = \cfrac{\mu_{max}}{S + K_s} \, . \end{eqnarray*}$$ +""" -where - -$$\mu = \mu_{max} \, \cfrac{S}{S + K_s}$$ +# ╔═╡ d635d577-4d6d-40d7-a84c-3871981d59a4 +md""" +Suppose that at $t=0$ no substrate $S$ is present in the reactor but that there is initially some biomass with a concentration of $0.0005\;g/L$. The (default) parameter values are $Q = 2.0$, $V = 40.0$ and $Y = 0.67$. The parameters $\mu_{max}$, $K_s$ and $S_{in}$ will be assigned later. """ # ╔═╡ 6ec6da23-853b-4129-94cf-67b5cadb1f95 @@ -52,7 +48,11 @@ The *reaction network object* model for this problem could be defined as: # ╔═╡ 935ca610-7a7a-4692-8908-fc26abb880b4 fermenter_monod = @reaction_network begin - mm(S, μmax, Ks)*X, S => Y*X + @species S(t)=0.0 X(t)=0.0005 + @parameters Q=2.0 V=40.0 Y=0.67 Sin μmax Ks + μmax/(S + Ks), S + X --> (1 + Y)*X + # Alternative: + # mm(S, μmax, Ks)*X, S => Y*X Q/V, (S, X) --> ∅ Q/V*Sin, ∅ --> S end @@ -78,101 +78,118 @@ Keep in mind that `mm(S, μmax, Ks)` stands for $\mu_{max} \, \cfrac{S}{S + K_s} """ # ╔═╡ 7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd -# osys = missing # Uncomment and complete the instruction +# osys = missing osys = convert(ODESystem, fermenter_monod) -# ╔═╡ 911b22bf-4cec-455d-a7f7-967bb55afea9 +# ╔═╡ 1ee7359c-2757-44cd-8650-09b762cf30d0 md""" -Check out the order of the parameters: +## Goals of this exercise """ -# ╔═╡ be15ae00-163c-44e5-bc33-e939ec63ed05 -# missing # Uncomment and complete the instruction -parameters(fermenter_monod) - # ╔═╡ 55f1d688-0c53-481b-9965-5e92ca87ad83 md""" -The parameter values are $\mu_{max} = 0.40$, $K_s = 0.015$, $Y = 0.67$, $Q = 2.0$, $V = 40.0$ and $S_{in} = 0.022\;g/L$. Suppose that at $t$ no substrate $S$ is present in the reactor but that there is initially some biomass with a concetration of $0.0005\;g/L$. Compute the following in a timespan of $[0, 100]\,h$: +Compute the following in a timespan of `(0.0, 100.0)`$\,h$: -- The sensitivities of $S$ and $X$ on $\mu_{max}$, $K_s$ and $S_{in}$. +- The sensitivities of $S$ and $X$ w.r.t. $\mu_{max}$, $K_s$ and $S_{in}$. Plot the following: -- A figure with the sensitivity functions of $S$ and $X$ on $S_{in}$. -- A figure with the sensitivity functions of $S$ on $\mu_{max}$, $K_s$ and $S_{in}$. -- A figure with the sensitivity functions of $X$ on $\mu_{max}$, $K_s$ and $S_{in}$. +- A figure with the sensitivity functions of $S$ and $X$ w.r.t. $S_{in}$. +- A figure with the sensitivity functions of $S$ w.r.t. $\mu_{max}$, $K_s$ and $S_{in}$. +- A figure with the sensitivity functions of $X$ w.r.t. $\mu_{max}$, $K_s$ and $S_{in}$. Interpret your results. + +Use the following parameter values $\mu_{max}$ = `0.40`, $K_s$ = `0.015` and $S_{in}$ = `0.022` for the calculations of the sensitivities. """ # ╔═╡ 6c4e3c09-4b84-4f5c-8739-2ac18e6f2af6 md""" -Initialize a vector `u0` with the initial conditions, set the timespan and initialize a vector `param` with the parameter values: +Initialize a vector `u0` with the initial conditions, set the timespan and initialize a vector `parms` with the parameter values: """ # ╔═╡ 2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c -# u0 = missing # Uncomment and complete the instruction +# u0 = missing # in principle not necessary since we will use the default u0 = [:S => 0.0, :X => 0.0005] # ╔═╡ 3fdc6b17-cdeb-4dc5-8886-9d3a62caac8d -# tspan = missing # Uncomment and complete the instruction +# tspan = missing tspan = (0.0, 100.0) -# ╔═╡ 0139da85-02e3-4021-9b39-84af7e68d428 +# ╔═╡ 201dfb54-2056-4846-a64c-ff4951dc084d md""" -For the sake of clarity, we will use the variables `μmax`, `Ks` and `Sin` to store the parameter values that are used for the calculation of the sensitivity functions. +For practical reasons, we will define the time step size as `dt = 0.5`. """ -# ╔═╡ 0f995929-4d2b-4a7a-8da1-04e4d501385f -# μmax = missing # Uncomment and complete the instruction -μmax = 0.40 +# ╔═╡ 5627068f-6442-43bb-bacb-bcba917372f3 +# dt = missing +dt = 0.5 -# ╔═╡ 262e8346-df6d-49bf-9186-92f5afb421e0 -# Ks = missing # Uncomment and complete the instruction -Ks = 0.015 +# ╔═╡ 0139da85-02e3-4021-9b39-84af7e68d428 +md""" +For practical reasons, we will use the variables `μmax_val`, `Ks_val`, and `Sin_val` to store the parameter values that are used for the calculation of the sensitivity functions. +""" -# ╔═╡ baa777d2-abb8-45f8-87aa-b3b17c8dc07c -# Sin = missing # Uncomment and complete the instruction -Sin = 0.022 +# ╔═╡ 0f995929-4d2b-4a7a-8da1-04e4d501385f +begin + # μmax_val = missing + # Ks_val = missing + # Sin_val = missing + μmax_val = 0.40 + Ks_val = 0.015 + Sin_val = 0.022 +end; # ╔═╡ 79b0eb65-5a0f-40b3-aa97-4088421c562e -# params = missing # Uncomment and complete the instruction -params = [:μmax => μmax, :Ks => Ks, :Y => 0.67, :Q => 2, :V => 40, :Sin => Sin] +# parms = missing # no need to include Q, V, and Y since we will use the default +parms = [:μmax => μmax_val, :Ks => Ks_val, :Sin => Sin_val] + +# ╔═╡ af882cf4-51fd-45ff-9b05-cfd52d6467b1 +md""" +## Preliminary simulation +""" # ╔═╡ f0f4fa14-6f99-4f21-a743-be61e08444a7 md""" -Create the ODE problem and store it in `oprob`. Next, solve the ODE problem using `Tsit5()` and `saveat=0.5`, and store the solution in `osol`. Finally plot the results. +Create the ODE problem and store it in `oprob`. Next, solve the ODE problem using `Tsit5()` and `saveat=dt`, and store the solution in `osol`. Finally plot the results. """ # ╔═╡ 8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 -# oprob = missing # Uncomment and complete the instruction -oprob = ODEProblem(fermenter_monod, u0, tspan, params) +# oprob = missing +oprob = ODEProblem(fermenter_monod, [], tspan, parms, combinatoric_ratelaws=false) # ╔═╡ 89a31c32-88a4-479f-a688-ffcb75ee8e91 -# osol = missing # Uncomment and complete the instruction -osol = solve(oprob, Tsit5(), saveat=0.5) +# osol = missing +osol = solve(oprob, Tsit5(), saveat=dt) # ╔═╡ 51a9b7e6-8ad9-477d-9596-ffd614df2c79 -# missing # Uncomment and complete the instruction +# missing plot(osol) +# ╔═╡ 570ebda9-b187-42ab-a761-bed0fa3ce097 +md""" +## Local Sensitivity Analysis (LSA) +""" + +# ╔═╡ 0343674d-32d1-49e1-a7c1-47196ae760ed +md""" +### Setting up function +""" + # ╔═╡ 693844d0-3858-4861-bae0-b47e78809f17 md""" Write a solution function with as argument a vector of the parameters (that you want the sensitivity on), and that returns the outputs. """ # ╔═╡ 9622f7ca-f71a-4ad9-a309-d7d10a1c3e3b -# Uncomment and complete the instruction -# function fermenter_monod_sim(params) +# function fermenter_monod_sim(parms) # missing # ... # end -function fermenter_monod_sim(params) - μmax, Ks, Sin = params - u0 = [:S => 0.0, :X => 0.0005] - tspan = (0.0, 100.0) - params = [:μmax => μmax, :Ks => Ks, :Y => 0.67, :Q => 2, :V => 40, :Sin => Sin] - oprob = ODEProblem(fermenter_monod, u0, tspan, params) # , combinatoric_ratelaws=false - osol = solve(oprob, Tsit5(), saveat=0.5) +function fermenter_monod_sim(parms) + μmax_val, Ks_val, Sin_val = parms + parms_new = [:μmax => μmax_val, :Ks => Ks_val, :Sin => Sin_val] + oprob = ODEProblem(fermenter_monod, [], tspan, parms_new, combinatoric_ratelaws=false) + osol = solve(oprob, Tsit5(), saveat=dt) return osol end @@ -182,12 +199,17 @@ Make two functions based on the solution function that each returns a single out """ # ╔═╡ f40c6402-3c28-4a7d-b629-83507a9f29bd -# fermenter_monod_sim_S(params) = missing # Uncomment and complete the instruction -fermenter_monod_sim_S(params) = fermenter_monod_sim(params)[:S] +# fermenter_monod_sim_S(parms) = missing +fermenter_monod_sim_S(parms) = fermenter_monod_sim(parms)[:S] # ╔═╡ 3ae5bd00-2e06-4789-aab3-d897824d5e29 -# fermenter_monod_sim_X(params) = missing # Uncomment and complete the instruction -fermenter_monod_sim_X(params) = fermenter_monod_sim(params)[:X] +# fermenter_monod_sim_X(params) = missing +fermenter_monod_sim_X(parms) = fermenter_monod_sim(parms)[:X] + +# ╔═╡ 6913fb1b-1986-4b37-819e-a6a7bd91f1a5 +md""" +### Compute the outputs +""" # ╔═╡ 4bd2bcca-9c42-4333-b062-2aaa9f7be3fe md""" @@ -195,8 +217,8 @@ Make the time vector. """ # ╔═╡ fbd98975-aa32-46ae-8db0-0e65cdf48309 -# t_vals = missing # Uncomment and complete the instruction -t_vals = 0:0.5:100.0 +# t_vals = missing +t_vals = tspan[1]:dt:tspan[2] # ╔═╡ 93791eb3-1eaa-4146-90b5-c4811fb3485b md""" @@ -204,25 +226,35 @@ Compute the two outputs $S$ and $X$ for the given parameter values. """ # ╔═╡ dc0557d6-81b9-4759-8ed7-3129f60c6dc3 -# S_sim = missing # Uncomment and complete the instruction -S_sim = fermenter_monod_sim_S([μmax, Ks, Sin]) +# S_sim = missing +S_sim = fermenter_monod_sim_S([μmax_val, Ks_val, Sin_val]) # ╔═╡ 95bc683c-f6e6-4b42-b90b-b5a863edd4d5 -# X_sim = missing # Uncomment and complete the instruction -X_sim = fermenter_monod_sim_X([μmax, Ks, Sin]) +# X_sim = missing +X_sim = fermenter_monod_sim_X([μmax_val, Ks_val, Sin_val]) + +# ╔═╡ 35717c5a-9895-4e62-a890-03ef04aa07b9 +md""" +### Compute the sensitivities +""" # ╔═╡ fa970c0e-fb3b-486f-bbc1-345d44f8f0da md""" Using `ForwardDiff.jacobian` to compute the sensitivities for the single ouputs $S$ and $X$. Hence, you need to call `ForwardDiff.jacobian` twice. """ +# ╔═╡ 07f2b8f8-5162-4077-a8a5-68fdec646644 +md""" +#### Absolute sensitivities +""" + # ╔═╡ 64354302-f4cc-4592-9302-5db0f5bccb2e -# sens_S = missing # Uncomment and complete the instruction -sens_S = ForwardDiff.jacobian(fermenter_monod_sim_S, [μmax, Ks, Sin]) +# sens_S = missing +sens_S = ForwardDiff.jacobian(fermenter_monod_sim_S, [μmax_val, Ks_val, Sin_val]) # ╔═╡ 49a94b9a-a543-495e-b4f1-c8579e59304d -# sens_X = missing # Uncomment and complete the instruction -sens_X = ForwardDiff.jacobian(fermenter_monod_sim_X, [μmax, Ks, Sin]) +# sens_X = missing +sens_X = ForwardDiff.jacobian(fermenter_monod_sim_X, [μmax_val, Ks_val, Sin_val]) # ╔═╡ 9cace6c1-e678-4dd7-8705-92a55eb32fa9 md""" @@ -230,7 +262,6 @@ Extract the (absolute) sensitivities of the outputs on the different parameters. """ # ╔═╡ a6dc2b60-6a0a-4140-892e-02cde8dc79d3 -# Uncomment and complete the instruction # begin # sens_S_on_μmax = missing # sens_S_on_Ks = missing @@ -240,10 +271,9 @@ begin sens_S_on_μmax = sens_S[:, 1] sens_S_on_Ks = sens_S[:, 2] sens_S_on_Sin = sens_S[:, 3] -end +end; # ╔═╡ f806c243-9032-46b7-add3-4714344691c7 -# Uncomment and complete the instruction # begin # sens_X_on_μmax = missing # sens_X_on_Ks = missing @@ -253,7 +283,12 @@ begin sens_X_on_μmax = sens_X[:, 1] sens_X_on_Ks = sens_X[:, 2] sens_X_on_Sin = sens_X[:, 3] -end +end; + +# ╔═╡ f4d0aaaa-bb3f-4335-901a-dbbcde643c78 +md""" +#### Normalized sensitivities +""" # ╔═╡ 5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 md""" @@ -261,30 +296,33 @@ Compute the normalized sensitivities. """ # ╔═╡ 76846731-929c-408f-a3de-970581c497e9 -# Uncomment and complete the instruction # begin # sens_S_on_μmax_rel = missing # sens_S_on_Ks_rel = missing # sens_S_on_Sin_rel = missing # end begin - sens_S_on_μmax_rel = sens_S_on_μmax .* μmax ./ S_sim - sens_S_on_Ks_rel = sens_S_on_Ks .* Ks ./ S_sim - sens_S_on_Sin_rel = sens_S_on_Sin .* Sin ./ S_sim -end + sens_S_on_μmax_rel = sens_S_on_μmax .* μmax_val ./ S_sim + sens_S_on_Ks_rel = sens_S_on_Ks .* Ks_val ./ S_sim + sens_S_on_Sin_rel = sens_S_on_Sin .* Sin_val ./ S_sim +end; # ╔═╡ b6c57444-547c-4e82-8526-6a30566e07c5 -# Uncomment and complete the instruction # begin # sens_X_on_μmax_rel = missing # sens_X_on_Ks_rel = missing # sens_X_on_Sin_rel = missing # end begin - sens_X_on_μmax_rel = sens_X_on_μmax .* μmax ./ X_sim - sens_X_on_Ks_rel = sens_X_on_Ks .* Ks ./ X_sim - sens_X_on_Sin_rel = sens_X_on_Sin .* Sin ./ X_sim -end + sens_X_on_μmax_rel = sens_X_on_μmax .* μmax_val ./ X_sim + sens_X_on_Ks_rel = sens_X_on_Ks .* Ks_val ./ X_sim + sens_X_on_Sin_rel = sens_X_on_Sin .* Sin_val ./ X_sim +end; + +# ╔═╡ 7f0b9e63-2856-4601-9f14-c20c6b1b3707 +md""" +### Plotting + questions +""" # ╔═╡ 5388c2a7-5a11-4da8-be09-46045cde8a4e md""" @@ -292,20 +330,26 @@ Plot the sensitivity functions of $S$ and $X$ on $S_{in}$. Provide a suitable ti """ # ╔═╡ db840c76-a6c6-49fb-a0bb-d9149f947bc0 -# missing # Uncomment and complete the instruction -plot(t_vals, [sens_S_on_Sin_rel, sens_X_on_Sin_rel], title="Normalized sensitivities", label=["S on Sin" "X on Sin"], xlabel="Time (hours)", linewidth=2) +# missing +plot(t_vals, [sens_S_on_Sin_rel, sens_X_on_Sin_rel], title="Normalized sensitivities", label=["S w.r.t. Sin" "X w.r.t. Sin"], xlabel="Time (hours)", linewidth=2) # ╔═╡ d41375ef-6958-4705-a417-4c6a491232ee md""" -Interpret your results. Try to answer the following question(s): -- Which output variable $S$ or $X$ is most sensitive on $S_{in}$ in steady state? - - Answer: missing -- Why is the sensitivity function of $S$ on $S_{in}$ at first positive but then becomes zero? - - Answer: missing +!!! questions + Interpret your results. Try to answer the following question(s): + - Which output variable $S$ or $X$ is most sensitive to $S_{in}$ in steady state? + - Why is the sensitivity function of $S$ on $S_{in}$ at first positive but then becomes zero? +""" + +# ╔═╡ 39c36eb4-5f15-4577-9b84-116e29d7891d +md""" +Answers: +- missing +- missing """ #= -- X is most sensitive on Sin because the higher Sin, the more X can be produced. -- In the beginning the substrate S is positively affected by Sin because S enters the tank through Sin and only little biomass X is present, so the biomass cannot consume the substrate very fast. +- X is most sensitive to Sin because the higher Sin, the more X can be produced. +- In the beginning, the substrate S is positively affected by Sin because S enters the tank through Sin, and only a little biomass X is present, so the biomass cannot consume the substrate very fast. =# # ╔═╡ be89600a-4927-4afc-9813-d8a70adb2852 @@ -314,23 +358,29 @@ Plot the sensitivity functions of $S$ on $\mu_{max}$, $K_s$ and $S_{in}$. Provid """ # ╔═╡ c0223da4-9959-48d0-b607-633b2e82986c -# missing # Uncomment and complete the instruction -plot(t_vals, [sens_S_on_μmax_rel, sens_S_on_Ks_rel, sens_S_on_Sin_rel], title="Normalized sensitivities", label=["S on μmax" "S on Ks" "S on Sin"], xlabel="Time (hours)", linewidth=2) +# missing +plot(t_vals, [sens_S_on_μmax_rel, sens_S_on_Ks_rel, sens_S_on_Sin_rel], title="Normalized sensitivities", label=["S w.r.t. μmax" "S w.r.t. Ks" "S w.r.t. Sin"], xlabel="Time (hours)", linewidth=2) # ╔═╡ ff86a29f-9308-473b-aa1c-dfd4af8179c7 md""" -Interpret your results. Try to answer the following question(s): -- Which parameter $\mu_{max}$, $K_s$ or $S_{in}$ affects the output $S$ the most in steady state? - - Answer: missing -- Why is the sensitivity function of $S$ on $K_s$ positive? - - Answer: missing -- Why is the sensitivity function of $S$ on $\mu_{max}$ negative? - - Answer: missing +!!! questions + Interpret your results. Try to answer the following question(s): + - Which parameter $\mu_{max}$, $K_s$ or $S_{in}$ affects the output $S$ the most in steady state? + - Why is the sensitivity function of $S$ w.r.t. $K_s$ positive? + - Why is the sensitivity function of $S$ w.r.t. $\mu_{max}$ negative? +""" + + +# ╔═╡ 430494da-402e-409d-85fc-c87cdfa3b427 +md""" +Answers: +- missing +- missing """ #= -- It seems like μmax is affecting S the most in steady state and its influence is negative, hence, the larger μmax, the smaller S. -- The sensitivity function of S on Ks is positive, because the larger Ks, the smaller the reaction rate r (S => Y*X). Hence, less X will be produced, so less S will be consumed. Remember that r = μmax*S*X/(S + Ks) and Ks is in the denominator. -- The sensitivity function of S on μmax is negative, because the larger μmax, the larger the reaction rate r (S => Y*X). Hence, more X will be produced, so more S will be consumed. Remember that r = μmax*S*X/(S + Ks) and μmax is in the numerator. +- It seems like μmax is affecting S the most in steady state, and its influence is negative; hence, the larger μmax, the smaller S. +- The sensitivity function of S w.r.t. Ks is positive, because the larger Ks, the smaller the reaction rate r (S => Y*X). Hence, less X will be produced; so less S will be consumed. Remember that r = μmax*S*X/(S + Ks) and Ks is in the denominator. +- The sensitivity function of S w.r.t. μmax is negative, because the larger μmax, the larger the reaction rate r (S => Y*X). Hence, more X will be produced; so more S will be consumed. Remember that r = μmax*S*X/(S + Ks) and μmax is in the numerator. =# # ╔═╡ 16a84fdb-8ce2-45b9-bfb7-7f4e1284a1d7 @@ -340,79 +390,392 @@ Plot the sensitivity functions of $X$ on $\mu_{max}$, $K_s$ and $S_{in}$. Provid # ╔═╡ 53134149-0bf7-41c1-9b35-e5037744211f # missing -plot(t_vals, [sens_X_on_μmax_rel, sens_X_on_Ks_rel, sens_X_on_Sin_rel], title="Normalized sensitivities", label=["X on μmax" "X on Ks" "X on Sin"], xlabel="Time (hours)", linewidth=2) +plot(t_vals, [sens_X_on_μmax_rel, sens_X_on_Ks_rel, sens_X_on_Sin_rel], title="Normalized sensitivities", label=["X w.r.t. μmax" "X w.r.t. Ks" "X w.r.t. Sin"], xlabel="Time (hours)", linewidth=2) # ╔═╡ 355ca6a7-466b-4969-ab48-28e2257f9810 md""" -Interpret your results. Try to answer the following question(s): -- Which parameter $\mu_{max}$, $K_s$ or $S_{in}$ affects the output $X$ the most in steady state? - - Answer: missing -- Why is the sensitivity function of $X$ on $K_s$ negative? - - Answer: missing -- Why is the sensitivity function of $X$ on $\mu_{max}$ positive? - - Answer: missing +!!! questions + Interpret your results. Try to answer the following question(s): + - Which parameter $\mu_{max}$, $K_s$ or $S_{in}$ affects the output $X$ the most in steady state? + - Why is the sensitivity function of $X$ w.r.t. $K_s$ negative? + - Why is the sensitivity function of $X$ w.r.t. $\mu_{max}$ positive? +""" + +# ╔═╡ 9a18a8bd-27d3-4771-a844-f6b65c7c8918 +md""" +Answers: +- missing +- missing +- missing """ #= -- It seems like Sin is affecting X the most in steady state and its influence is positive, hence, the larger Sin, the larger X. -- The sensitivity function of X on Ks is negative, because the larger Ks, the smaller the reaction rate r (S => Y*X). Hence, less X will be produced. Remember that r = μmax*S*X/(S + Ks) and Ks is in the denominator. -- The sensitivity function of X on μmax is positive, because the larger μmax, the larger the reaction rate r (S => Y*X). Hence, more X will be produced. Remember that r = μmax*S*X/(S + Ks) and μmax is in the numerator. +- It seems like Sin is affecting X the most in steady state, and its influence is positive; hence, the larger Sin, the larger X. +- The sensitivity function of X w.r.t. Ks is negative, because the larger Ks, the smaller the reaction rate r (S => Y*X). Hence, less X will be produced. Remember that r = μmax*S*X/(S + Ks) and Ks is in the denominator. +- The sensitivity function of X w.r.t. μmax is positive, because the larger μmax, the larger the reaction rate r (S => Y*X). Hence, more X will be produced. Remember that r = μmax*S*X/(S + Ks) and μmax is in the numerator. =# +# ╔═╡ 6536f76a-0f6d-4155-ab3b-38e2e3189886 +md""" +## Monte Carlo error propagation +""" + +# ╔═╡ 9fd59148-6993-4598-af95-ac85b6a4bb45 +md""" +Let's take a look at how the uncertainty propagates in your model by using Monte Carlo simulations. Based on some literature search you can assume that the parameter `Ks` is normally distributed around the original value with a standard deviation of 20%. +""" + +# ╔═╡ 5a559bc2-3d55-499c-93cf-4b0845ffc4b0 +md""" +Start with making a Turing model in which you implement the prior and return the solution of the solved problem. +""" + +# ╔═╡ 1ee0fbfb-0e64-4db5-bed9-caa2aff7cf0b +md""" +!!! note + Instead of creating a new ODEProblem everytime, you can simply remake an old ODEProblem by using `new_prob = remake(oprob, p = [:Ks => Ks])`. +""" + +# ╔═╡ 26ae811d-11fd-44d8-b328-79381fa7472c +# @model function monod_deviation() +# Ks_dev ~ missing # 20% standard deviation +# new_prob = missing +# return missing +# end +@model function monod_deviation() + Ks_dev ~ Normal(0.015, 0.0030) # 20% standard deviation + new_prob = remake(oprob, p = [:Ks => Ks_dev]) + return solve(new_prob, AutoTsit5(Rosenbrock23()), saveat=0.5); +end + +# ╔═╡ 62791c48-5e38-46fe-8336-3a21003f3fe3 +md""" +Now get sample 100 solutions from your monod_deviation model. +""" + +# ╔═╡ 9e049653-347a-4198-90fd-23f15c404d0c +# missing + +begin + dev_model = monod_deviation() + chain = sample(dev_model, Prior(), 100) + solutions = generated_quantities(dev_model, chain) +end; + +# ╔═╡ f11827e5-de69-4a0c-b385-2c6ea031bb79 +md""" +We can now visualise the results of our Monte Carlo simulation by looping through our solutions and plotting them. It is advised to put `label = false` (avoids 100 labels), `color=:gray` (neutral color) and `alpha=0.4` (makes it more transparent) and `linestyle=:dot`. +""" + +# ╔═╡ 9d0f8147-3c82-4ace-a408-20e32b4d8ea0 +# begin +# p1 = plot(title="Monte Carlo of logistic growth") # make an empty plot +# for missing in missing +# missing +# end +# plot!(osol, label = "Original simulation", linewidth = 3) # Plots original result +# p1 +# end + +begin + p1 = plot(title="Monte Carlo of logistic growth", ylabel = "concentration (g/l)", xlabel = "time (hours)") # make an empty plot + for solution in solutions + plot!(solution, label = false, color =:gray, alpha = 0.4, + linestyle = :dot) + end + plot!(osol, label = "Original simulation", linewidth = 3) + p1 +end + +# ╔═╡ 30c270f1-cbd5-401b-92db-7123efb01db6 +md""" +As can be seen from the graph, the uncertainty is rather high, especially if we want to carefully monitor the concentration of biomass through time. Note that, since we have a sample of values for $S$ and $X$ at every time step, it is also perfectly possible here to quantify the uncertainty by calculating the standard deviation at every time step, but the graph already gives a clear indication. +""" + +# ╔═╡ 16098df3-35c0-4d06-a482-35137fd53b97 +md""" +Now create a vector x1 that stores the values of the biomass at the end of the simulations and make a histogram of it. Put `xlims = [0.010, 0.015]`. +""" + +# ╔═╡ de5736d4-1c2f-4b0c-a84b-6d8b07d29a0a +# missing +x1 = [solution[:X][end] for solution in solutions]; + +# ╔═╡ 1b1d0997-b854-4445-8c19-e73db6784a9b +# missing +histogram(x1, xlims = [0.010, 0.015]) + +# ╔═╡ 43aaf563-f799-43e4-9f4c-8a19be1586dd +md""" +### Determininig optimal measurements +""" + +# ╔═╡ f7f1f035-b48e-4108-bb84-0c00ab68a93a +md""" +As a process operator, you may want to reduce the uncertainty in biomass simulations. However, since taking measurements can be costly, it is important to perform them at the most informative time points. + +Local sensitivity analysis is a useful tool in this context, as it indicates when the model output is most sensitive to changes in parameter values. At these time points, accurate measurements provide the greatest insight into the true parameter values. + +Start off by identifying the time at which the biomass concentration is most sensitive to variations in `Ks`. +""" + +# ╔═╡ 6178c0ca-f8ab-49cf-aa9c-82c1dea9b8cf +md""" +!!! note + You can find the index of the maximum value in a vector with the `argmax()` function. For example `argmax([1, 2, 5, 4])` will give you 3. +""" + +# ╔═╡ ddfb530a-ba79-43c0-973a-4deb8cad1d33 +# begin +# t_star_idx = missing +# t_star = t_vals[t_star_idx] + +# plot(t_vals, abs.(sens_X_on_Ks_rel), +# title="Sensitivity of Ks wrt. X", xlabel="Time (h)", +# ylabel="normalized sensitivity") +# vline!([t_star], label="t* = $(round(t_star, digits=2)) h", linestyle=:dash, color=:red) +# end +begin + t_star_idx = argmax(abs.(sens_X_on_Ks_rel)) + t_star = t_vals[t_star_idx] + + plot(t_vals, abs.(sens_X_on_Ks_rel), + title="Sensitivity of Ks wrt. X", xlabel="Time (h)", + ylabel="normalized sensitivity") + vline!([t_star], label="t* = $(round(t_star, digits=2)) h", linestyle=:dash, color=:red) +end + +# ╔═╡ 8359ea0e-a916-48f4-a07a-c905c3efaf15 +begin + μmax_real = 0.42 + Ks_real = 0.0165 + Sin_real = 0.0207 + real_solution = solve(remake(oprob, p = [:μmax => μmax_real, :Ks => Ks_real, + :Sin => Sin_real])); +end; + +# ╔═╡ 9654dea2-e13e-4aa9-8a5e-c3c368c26e24 +md""" +We have defined a function `real_solution(t)` in this notebook that returns a hypothetical measurement (`S_meas, X_meas`) at the given timepoint `t`. Use this function to get a measurement at the optimal time (`t_star`). +""" + +# ╔═╡ 77eab4ec-5012-40dd-b702-9fb48982ff38 +# S_meas, X_meas = real_solution(t_star) +S_meas, X_meas = real_solution(t_star) + +# ╔═╡ 3e93a001-9401-4a1f-b283-c79fd0c43434 +md""" +#### Conditioning on the measurements +""" + +# ╔═╡ 684f9163-ac84-49be-beb1-66d7bdc21da8 +md""" +Now we can use this measurement to calibrate our model. In this part we will use Monte Carlo Markov chain to show the decrease in uncertainty. + +We start by creating a Turing model, for which you can use the priors from the first part of this exercise. For the standard deviation `σ_X` you can assume a value of 0.0002 and a value of 0.00055 for `σ_S`. + +Make sure to truncate the domain of your parameters to a reasonable domain or you might run into issues with the solver. +""" + +# ╔═╡ 4fe68607-d67e-4e71-865a-b10a06e5908d +# @model function monod_meas() +# σ_X = missing +# σ_S = missing + +# Ks_dev ~ missing + +# sol = missing + +# S_pred, X_pred = sol(t_meas) + +# X_meas ~ missing +# S_meas ~ missing +# end +@model function monod_meas(t_meas) + σ_X = 0.05 * 0.004 # ~5% relative noise + σ_S = 0.05 * 0.011 + + Ks_dev ~ truncated(Normal(0.015, 0.0030), lower = 0.0, upper = 1.0) + + sol = solve(remake(oprob, p = [:Ks => Ks_dev]), AutoTsit5(Rosenbrock23()), saveat=0.1) + + S_pred, X_pred = sol(t_meas) + + X_meas ~ Normal(X_pred, σ_X) + S_meas ~ Normal(S_pred, σ_S) + return solve(remake(oprob, p = [:Ks => Ks_dev]), saveat=0.5) +end + +# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 +# monod_cond_model = missing +monod_cond_model = monod_meas(t_star) | (X_meas = X_meas, + S_meas = S_meas) + +# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 +md""" +### MCMC +""" + +# ╔═╡ b8cd84de-57e4-4d8b-92da-e7c5b6d627bd +# monod_chain = missing +monod_chain = sample(monod_cond_model, NUTS(), 100) + +# ╔═╡ e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +md""" +Plot the chain to check if it properly converged to the posterior distribution. +""" + +# ╔═╡ 834c3500-7556-4999-b465-08a8a07c7f28 +# missing +plot(monod_chain) + +# ╔═╡ eeaf3255-95c8-4813-bdd7-38635e546606 +md""" +Now, get the sampled posterior solutions and afterwards plot the resulting simulations. +""" + +# ╔═╡ a66b6d91-e529-4c20-b5d6-91830eae1437 +# missing +solutions_monod = generated_quantities(monod_cond_model, monod_chain); + +# ╔═╡ c4192ff2-fca0-4ab1-8ac1-ce9771bc69d3 +md""" +Plot the resulting Monte Carlo simulation with the original simulations. +""" + +# ╔═╡ 54f20cff-21fc-4169-9f47-5e670b7ec621 +# begin +# p = plot(ylabel = "concentration (g/l)", xlabel = "Time (hours)") +# for missing in missing +# missing +# end +# plot!(osol, label = ["Original simulation S" "Original simulation X"], linewidth = 3) + +# p +# end +begin + p = plot(ylabel = "concentration (g/l)", xlabel = "Time (hours)") + for solution in solutions_monod + plot!(solution, label=false, color=:grey, linestyle=:dot, alpha=0.4) + end + plot!(osol, label = ["Original simulation S" "Original simulation X"], linewidth = 3) + + p +end + +# ╔═╡ 988aab6d-7193-4f3d-bd08-c3c36c920f81 +md""" +Lastly, create a histogram of the value of X at the end of the simulations based on your posterior distribution of Ks. Put the `xlims = [0.010, 0.015]` to compare with the prior distribution. +""" + +# ╔═╡ 2f55a771-0d37-48e5-a9d9-24693bf1d177 +# missing +x2 = [solution_monod[:X][end] for solution_monod in solutions_monod]; + +# ╔═╡ 18059d0c-cf37-44a1-86a8-859ac66a66f0 +# missing +histogram(x2, title="Posterior histogram of X", ylabel = "counts", xlabel="X", xlims = [0.010, 0.015]) + # ╔═╡ Cell order: +# ╟─31d294d1-3a1f-41db-abff-54f2a67c7ed9 # ╠═55cdebd2-0881-11ef-2722-91de1447877a -# ╠═03ae0690-06a0-4276-9f00-d07b206fe124 # ╠═3ef93246-657d-4e77-9bf0-8380c64bfcfd # ╠═a355b0ba-baaf-49f4-a5dc-965364a884f0 # ╠═00fd6d49-f561-42e9-9413-d33af92f83dc # ╠═7ae714c4-d25d-4f9f-ab3d-cc067db9c156 -# ╟─31d294d1-3a1f-41db-abff-54f2a67c7ed9 +# ╠═dfce0717-e33e-4b6b-bb5c-d8754a74aba4 # ╟─5ffe7dcb-620d-4f22-95fe-2f77cda6fbe7 +# ╟─d635d577-4d6d-40d7-a84c-3871981d59a4 # ╟─6ec6da23-853b-4129-94cf-67b5cadb1f95 # ╠═935ca610-7a7a-4692-8908-fc26abb880b4 # ╟─79e6056a-881c-442f-8989-5bc284d3d777 # ╟─fa93e2c3-8b43-418e-ba24-406645b2e397 # ╟─06730f54-7293-43f2-b772-84eec3e5528a # ╠═7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd -# ╟─911b22bf-4cec-455d-a7f7-967bb55afea9 -# ╠═be15ae00-163c-44e5-bc33-e939ec63ed05 +# ╟─1ee7359c-2757-44cd-8650-09b762cf30d0 # ╟─55f1d688-0c53-481b-9965-5e92ca87ad83 # ╟─6c4e3c09-4b84-4f5c-8739-2ac18e6f2af6 # ╠═2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c # ╠═3fdc6b17-cdeb-4dc5-8886-9d3a62caac8d +# ╟─201dfb54-2056-4846-a64c-ff4951dc084d +# ╠═5627068f-6442-43bb-bacb-bcba917372f3 # ╟─0139da85-02e3-4021-9b39-84af7e68d428 # ╠═0f995929-4d2b-4a7a-8da1-04e4d501385f -# ╠═262e8346-df6d-49bf-9186-92f5afb421e0 -# ╠═baa777d2-abb8-45f8-87aa-b3b17c8dc07c # ╠═79b0eb65-5a0f-40b3-aa97-4088421c562e +# ╟─af882cf4-51fd-45ff-9b05-cfd52d6467b1 # ╟─f0f4fa14-6f99-4f21-a743-be61e08444a7 # ╠═8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 # ╠═89a31c32-88a4-479f-a688-ffcb75ee8e91 # ╠═51a9b7e6-8ad9-477d-9596-ffd614df2c79 +# ╟─570ebda9-b187-42ab-a761-bed0fa3ce097 +# ╟─0343674d-32d1-49e1-a7c1-47196ae760ed # ╟─693844d0-3858-4861-bae0-b47e78809f17 # ╠═9622f7ca-f71a-4ad9-a309-d7d10a1c3e3b # ╟─4a5971b1-f4d0-43b6-805f-e17f5052ae92 # ╠═f40c6402-3c28-4a7d-b629-83507a9f29bd # ╠═3ae5bd00-2e06-4789-aab3-d897824d5e29 +# ╟─6913fb1b-1986-4b37-819e-a6a7bd91f1a5 # ╟─4bd2bcca-9c42-4333-b062-2aaa9f7be3fe # ╠═fbd98975-aa32-46ae-8db0-0e65cdf48309 # ╟─93791eb3-1eaa-4146-90b5-c4811fb3485b # ╠═dc0557d6-81b9-4759-8ed7-3129f60c6dc3 # ╠═95bc683c-f6e6-4b42-b90b-b5a863edd4d5 +# ╟─35717c5a-9895-4e62-a890-03ef04aa07b9 # ╟─fa970c0e-fb3b-486f-bbc1-345d44f8f0da +# ╟─07f2b8f8-5162-4077-a8a5-68fdec646644 # ╠═64354302-f4cc-4592-9302-5db0f5bccb2e # ╠═49a94b9a-a543-495e-b4f1-c8579e59304d # ╟─9cace6c1-e678-4dd7-8705-92a55eb32fa9 # ╠═a6dc2b60-6a0a-4140-892e-02cde8dc79d3 # ╠═f806c243-9032-46b7-add3-4714344691c7 +# ╟─f4d0aaaa-bb3f-4335-901a-dbbcde643c78 # ╟─5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 # ╠═76846731-929c-408f-a3de-970581c497e9 # ╠═b6c57444-547c-4e82-8526-6a30566e07c5 +# ╟─7f0b9e63-2856-4601-9f14-c20c6b1b3707 # ╟─5388c2a7-5a11-4da8-be09-46045cde8a4e # ╠═db840c76-a6c6-49fb-a0bb-d9149f947bc0 # ╟─d41375ef-6958-4705-a417-4c6a491232ee +# ╠═39c36eb4-5f15-4577-9b84-116e29d7891d # ╟─be89600a-4927-4afc-9813-d8a70adb2852 # ╠═c0223da4-9959-48d0-b607-633b2e82986c # ╟─ff86a29f-9308-473b-aa1c-dfd4af8179c7 +# ╟─430494da-402e-409d-85fc-c87cdfa3b427 # ╟─16a84fdb-8ce2-45b9-bfb7-7f4e1284a1d7 # ╠═53134149-0bf7-41c1-9b35-e5037744211f # ╟─355ca6a7-466b-4969-ab48-28e2257f9810 +# ╟─9a18a8bd-27d3-4771-a844-f6b65c7c8918 +# ╟─6536f76a-0f6d-4155-ab3b-38e2e3189886 +# ╟─9fd59148-6993-4598-af95-ac85b6a4bb45 +# ╟─5a559bc2-3d55-499c-93cf-4b0845ffc4b0 +# ╟─1ee0fbfb-0e64-4db5-bed9-caa2aff7cf0b +# ╠═26ae811d-11fd-44d8-b328-79381fa7472c +# ╟─62791c48-5e38-46fe-8336-3a21003f3fe3 +# ╠═9e049653-347a-4198-90fd-23f15c404d0c +# ╟─f11827e5-de69-4a0c-b385-2c6ea031bb79 +# ╠═9d0f8147-3c82-4ace-a408-20e32b4d8ea0 +# ╟─30c270f1-cbd5-401b-92db-7123efb01db6 +# ╟─16098df3-35c0-4d06-a482-35137fd53b97 +# ╠═de5736d4-1c2f-4b0c-a84b-6d8b07d29a0a +# ╠═1b1d0997-b854-4445-8c19-e73db6784a9b +# ╟─43aaf563-f799-43e4-9f4c-8a19be1586dd +# ╟─f7f1f035-b48e-4108-bb84-0c00ab68a93a +# ╟─6178c0ca-f8ab-49cf-aa9c-82c1dea9b8cf +# ╠═ddfb530a-ba79-43c0-973a-4deb8cad1d33 +# ╟─8359ea0e-a916-48f4-a07a-c905c3efaf15 +# ╟─9654dea2-e13e-4aa9-8a5e-c3c368c26e24 +# ╠═77eab4ec-5012-40dd-b702-9fb48982ff38 +# ╟─3e93a001-9401-4a1f-b283-c79fd0c43434 +# ╟─684f9163-ac84-49be-beb1-66d7bdc21da8 +# ╠═4fe68607-d67e-4e71-865a-b10a06e5908d +# ╠═11e105de-270c-434a-bfb4-9118dfc9e461 +# ╟─e46afb98-16de-4d19-86b4-b920b6b929c5 +# ╠═b8cd84de-57e4-4d8b-92da-e7c5b6d627bd +# ╟─e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +# ╠═834c3500-7556-4999-b465-08a8a07c7f28 +# ╟─eeaf3255-95c8-4813-bdd7-38635e546606 +# ╠═a66b6d91-e529-4c20-b5d6-91830eae1437 +# ╟─c4192ff2-fca0-4ab1-8ac1-ce9771bc69d3 +# ╠═54f20cff-21fc-4169-9f47-5e670b7ec621 +# ╟─988aab6d-7193-4f3d-bd08-c3c36c920f81 +# ╠═2f55a771-0d37-48e5-a9d9-24693bf1d177 +# ╠═18059d0c-cf37-44a1-86a8-859ac66a66f0 diff --git a/exercises/solved_notebooks/P7_sens_uncert/uncert_intro.jl b/exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl similarity index 59% rename from exercises/solved_notebooks/P7_sens_uncert/uncert_intro.jl rename to exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl index caeda3d..12b8554 100644 --- a/exercises/solved_notebooks/P7_sens_uncert/uncert_intro.jl +++ b/exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl @@ -5,45 +5,41 @@ using Markdown using InteractiveUtils # ╔═╡ 8349306c-e98c-4221-9a1b-322fde3e18cb -begin - # add this cell if you want the notebook to use the environment from where the Pluto server is launched - using Pkg - Pkg.activate("..") -end +using Pkg; Pkg.activate("..") # ╔═╡ d2c4d230-0943-11ef-3aad-5719e74bb20e -using Markdown - -# ╔═╡ 47790920-74d2-4f21-b0f2-45dc19d25156 -using InteractiveUtils +using Markdown, InteractiveUtils # ╔═╡ 83f3d978-aefa-40c1-a647-b26e837aeed6 -using Catalyst - -# ╔═╡ 73d5f921-82e4-4160-a3c4-fc1935f4c58d -using OrdinaryDiffEq, StatsPlots +using OrdinaryDiffEq, Catalyst # ╔═╡ 5185d0eb-7392-4775-9336-3e0f9e1449ce using Measurements -# ╔═╡ 1f5e389e-e003-4e73-8f18-b5ad6340a912 -Catalyst.ModelingToolkit.NaNMath.log(x::Measurement) = Measurements.log(x) - # using `log` on a Measurement-type variable in the reaction rate of a Catalyst model will otherwise error in the current Catalyst version +# ╔═╡ 1c8ea6fe-650b-4c52-ba14-47efe3bd3e39 +using StatsPlots, PlutoUI; TableOfContents() + +# ╔═╡ 3eb3e651-afcf-41f8-a652-154a4f7ae07f +using Turing # ╔═╡ 241a8a65-c59f-44f1-be39-d5edd1321b49 -md" +md""" # Introduction to uncertainty analysis -" +""" + +# ╔═╡ 1f5e389e-e003-4e73-8f18-b5ad6340a912 +Catalyst.ModelingToolkit.NaNMath.log(x::Measurement) = Measurements.log(x) + # using `log` on a Measurement-type variable in the reaction rate of a Catalyst model will otherwise error in the current Catalyst version # ╔═╡ 0f518e36-bc96-4599-807e-728504e5ca7b -md" +md""" ## Goal of this practicum -" +""" # ╔═╡ da3cdc89-b911-4af9-9d4a-526301cba581 -md" -Parameter uncertainty plays a crucial role in shaping the behavior of output variables. Model equations describe how systems evolve over time, often incorporating parameters representing various aspects of the system's characteristics. However, these parameters are rarely known with absolute certainty and they often come with inherent uncertainty due to measurement errors, variability in real-world conditions, or incomplete knowledge about the system. This uncertainty can propagate through the model equations, leading to uncertainties in the predicted outcomes. Consequently, understanding the influence of parameter uncertainty becomes essential for assessing the reliability and robustness of the model predictions, as well as for making informed decisions based on these predictions. Techniques such as sensitivity analysis and uncertainty quantification are employed to explore and quantify the impact of parameter uncertainty on the output variables, providing insights into the system's behavior and guiding the refinement of models for improved accuracy and reliability. -" +md""" +Parameter uncertainty plays a crucial role in shaping the behavior of output variables. Model equations describe how systems evolve, often incorporating parameters representing various aspects of the system's characteristics. However, these parameters are rarely known with absolute certainty and often carry inherent uncertainty due to measurement errors, variability in real-world conditions, or incomplete knowledge of the system. This uncertainty can propagate through the model equations, leading to uncertainties in the predicted outcomes. Consequently, understanding the influence of parameter uncertainty becomes essential for assessing the reliability and robustness of the model predictions, as well as for making informed decisions based on these predictions. Techniques such as sensitivity analysis and uncertainty quantification are employed to explore and quantify the impact of parameter uncertainty on the output variables, providing insights into the system's behavior and guiding the refinement of models for improved accuracy and reliability. +""" # ╔═╡ 78afaded-5a19-4386-aa76-7974977ea354 md" @@ -56,108 +52,111 @@ We will now compute the variability in the output variables replected as error b " # ╔═╡ e5ab4490-6fd4-4b51-bfa7-1366438efefc -md" +md""" ## Grass growth models -" +""" # ╔═╡ 3e2750e7-220a-47b2-b445-eb3315734dec -md" +md""" In this notebook, three different models will be used, each modelling the yield of grass in a grassland: - Logistic growth model: $\cfrac{dW}{dt} = \mu \left( 1 - \cfrac{W}{W_f} \right) W$ - Exponential growth model: $\cfrac{dW}{dt} = \mu \left( W_f - W \right)$ -- Gompertz growth model: $\cfrac{dW}{dt} = \left( \mu - D \ln(W) \right) W$ +- Gompertz growth model: $\cfrac{dW}{dt} = \left( \mu - d \ln(W) \right) W$ with output $W$ the grass yield, and $W_f$, $\mu$ and $D$ parameters. The table below show some typical values for the parameters together with their uncertainties: -| | $\mu$ | $W_f$ | $D$ | +| | $\mu$ | $W_f$ | $d$ | |:----------- |:----------:|:-----------:|:------------:| | Logistic | 0.07$\pm$0.02 | 10.0$\pm$0.15 | | | Exponential | 0.02$\pm$0.01 | 10.0$\pm$0.15 | | | Gompertz | 0.09$\pm$0.01 | | 0.040$\pm$0.002 | We will use an initial condition of $W_0 = 2.0$ for each and a simulation time of $100$ days. -" +""" # ╔═╡ c874a08c-7d82-4631-b611-c598dcaded09 -md" +md""" We will illustrate how to compute the error bars, due to model parameter uncertainty, in conjunction with the output variable simulation for the logistic model. The same will be left as exercises below for the exponential and gompertz models. **Important:** - We will use consequently `_log`, `_exp` and `_gom` appended to relevant variables names in order to indicate their model origin **and** to prevent cell-disabling that occurs when using the same variables names in these Notebooks. -" +""" -# ╔═╡ 9d3192bc-ee69-44a1-9273-06b8a55c60ef -md" +# ╔═╡ 87272d30-d95a-4b36-84ef-84bb9363c19e +md""" ### Uncertainty analysis of the logistic growth model +""" -We will start by modelling our system and simulating using the aforementioned parameters values and uncertainties, initial condition and timespan. -" +# ╔═╡ 9d3192bc-ee69-44a1-9273-06b8a55c60ef +md""" +We will start by modelling our system and simulating using the aforementioned parameter values and uncertainties, initial condition, and timespan. +""" # ╔═╡ 8acb7ea2-54ef-428b-bb4f-364377d2c38d -md" +md""" Implementation of the system: -" +""" # ╔═╡ 460d98ef-2177-49bf-87a4-35412f4183ed growth_mod_log = @reaction_network begin - #μ*W, 0 --> W - #μ/Wf*W, W --> 0 + @species W(t)=2.0 + @parameters μ Wf μ*(1-W/Wf), W --> 2W end # ╔═╡ 38e0f746-1d04-46ba-81cc-21522b9ce6a4 -md" +md""" Convert the *reaction model* to check that we work with the correct differential equation: -" +""" # ╔═╡ 2538bd6e-027f-47fe-813c-db0a02586d2c osys_log = convert(ODESystem, growth_mod_log) # ╔═╡ 72b53ba6-dd05-4f2d-a65e-52783e76a1e9 -md" -In order to use uncertainties in the parameter values, we will need to load the package `Measurements`: -" +md""" +In order to use uncertainties in the parameter values, we will need to load the package `Measurements` (see above). +""" # ╔═╡ 1a67db96-6973-446d-aa5b-253dd0008a73 -md" +md""" Setting initial conditions: -" +""" # ╔═╡ 048112b0-10f1-47c5-834f-c3851d3078e8 u0_log = [:W => 2.0] # ╔═╡ 037cbf2b-dc0c-4b62-879b-aa702c964877 -md" +md""" Set the timespan for the simulation: -" +""" # ╔═╡ 40c1192c-bc76-4d3a-8144-dd76ba72fba2 tspan = (0.0, 100.0) # ╔═╡ 0b5493eb-6ce6-43c1-9289-4c2e6f20d014 -md" -In order to see the effect of the individual parameter uncertainties on the output variable, we will analyse this (only for didactical reasons) assuming different study cases: +md""" +To see the effect of the individual parameter uncertainties on the output variable, we will analyse this (only for didactical reasons), assuming different study cases: -1. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ while $W_f$ is exactly known and equal to $0.15$. -2. $K_s$ has a nominal value of $10.0$ and a standard deviation of $0.15$ while $\mu$ is exactly known and equal to $0.07$. -3. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ and similarily $K_s$ has a nominal value of $10.0$ and a standard deviation of $0.15$. +1. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ while $W_f$ is exactly known and equal to $10.0$. +2. $W_f$ has a nominal value of $10.0$ and a standard deviation of $0.15$ while $\mu$ is exactly known and equal to $0.07$. +3. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ and similarily $W_f$ has a nominal value of $10.0$ and a standard deviation of $0.15$. -The last study case is the true case because in reality all encertainties will contribute at once. -" +The last case study is the true case because, in reality, all uncertainties will contribute simultaneously. +""" # ╔═╡ b4c76f21-9839-4db5-b952-43869d3a1efe -md" -#### Study case 1 -" +md""" +### Case study 1 +""" # ╔═╡ 6ab73fc5-7315-4ec6-a3ee-9b5c689ce82e -md" -We initialize a vector `params1_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the first parameter**: -" +md""" +We initialize a vector `parms1_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the first parameter**: +""" # ╔═╡ b0f32f69-b6a2-4050-97e2-8018d07463d8 -params1_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0] +parms1_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0] # ╔═╡ c1a0702c-6fa5-47d7-bd72-6bfcbdeba93f md" @@ -165,20 +164,20 @@ We create the corresponding ODE problem and store it in `oprob1_uncert_log`: " # ╔═╡ 2fc91add-8468-4402-bcd8-42da5544f612 -oprob1_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, params1_uncert_log) +oprob1_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, parms1_uncert_log) # ╔═╡ 5fa951d3-23ee-4510-873c-8158df4f2faf -md" +md""" We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol1_uncert_log`: -" +""" # ╔═╡ d54635e8-1693-42bc-91f4-7e504a8661c7 osol1_uncert_log = solve(oprob1_uncert_log, Tsit5(), saveat=2.0) # ╔═╡ 707c6023-49ae-4a13-a220-d827a76a352c -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ 48e92497-f6b3-484c-bc5f-3ce69f4f1a91 plot(osol1_uncert_log) @@ -188,101 +187,204 @@ md" When thinking back of the sensitivity of $W$ to the parameter $\mu$, we saw that the corresponding sensitivity function had a maximum around $t=33\;s$. Looking at the above plot with error bars, we can see that de largest error bars (largest uncentrainty in the output variable) occurs at the timepoints where the sensitivity is strongest. " +# ╔═╡ 5dcb1ef0-d2a8-41fa-9592-2a86cf8364e6 +md""" +### Alternative: Monte Carlo uncertainty propagation +""" + +# ╔═╡ 0d883d09-1c92-4bb5-9bbf-3c87c856a5fa +md""" +In uncertainty propagation, we are interested in the effect of uncertainties (or errors) on the final output of the model. Monte Carlo simulations are a straightforward way for testing different parameters and checking the results of the outcomes. + +Below a Turing model is defined that samples values from a normal distribution with a 5% standard deviation. This is our Prior belief of how accurate the estimate of our intial parameters is. Remember that this is similar to what we did in practical 4. +""" + +# ╔═╡ 0e0dc268-faaf-462e-a980-9398c6015e02 +@model function logistic_deviation() + μ_dev ~ Normal(0.07, 0.0035) # 5% standard deviation + Wf_dev ~ Normal(10, 0.5) # 5% standard deviation + return solve(remake(oprob_log, p = [:Wf => Wf_dev, :μ => μ_dev]), saveat=0.5); +end + +# ╔═╡ e34f4b23-ac1b-4b88-b959-65085cef4b4f +md""" +Using the `sample()` function we can sample values from our distribution. +""" + +# ╔═╡ 738448dd-d1c9-4fdc-af5d-fad4d91adf38 +μ_model = logistic_deviation(); + +# ╔═╡ 51579212-dd33-42b9-a19e-221ab7630c1a +chain = sample(μ_model, Prior(), 500); + +# ╔═╡ 22bf3f50-1909-4a2c-820f-45b87f6bce14 +solutions = generated_quantities(μ_model, chain); + +# ╔═╡ 616ee1d4-1f56-4c4a-a658-9f9351949f59 +md""" +Finally, we can use theses perturbed values to model the effects of these small perturbations on the output of the model. +""" + +# ╔═╡ 6b4a862e-6fb6-40ba-9e3f-43d57f169994 +md""" +Firstly, solve the original system (Wf = 10 and μ = 0.07) and afterwards we can compare this with the outputs of the Monte Carlo simulations. +""" + +# ╔═╡ 62afc495-29c3-4baa-b45a-94b84149d15d +oprob_log = ODEProblem(growth_mod_log, u0_log,tspan, [:μ => 0.07, :Wf => 10]); + +# ╔═╡ 691a8b68-935e-4012-aa25-9cfdc323b1de +osol_log = solve(oprob_log, Tsit5(), saveat = 0.5); + +# ╔═╡ f49a0085-396f-469e-afb8-2b7884f9bb1b +md""" +Let's now plot the solution of the Monte Carlo experiments. +""" + +# ╔═╡ ff5fcd0a-a37a-46d5-8b4b-f3f7a9a457a3 +begin + p1 = plot(title="Monte Carlo of logistic growth") + for solution in solutions + plot!(solution, label = false, color =:gray, alpha = 0.4, + linestyle = :dot) + end + plot!(osol_log, label = "Original simulation", linewidth = 3) + p1 +end + +# ╔═╡ c0e72c96-23b6-4903-bd31-2978f6aebb47 +md""" +An advantage of using Monte Carlo over the error bars, is that the resulting graphs create a distribution from which we can calculate probabilities. We could now for example calculate what the average yield is after 50 days: +""" + +# ╔═╡ 988e9cd3-29a8-4fb1-a706-82d8356b9f20 +md""" +It is worth noting that you can obtain the solution at a specific time point by indexing it using parentheses `()`, as shown below. However, the result is still returned as a vector (indicated by the square brackets `[]`). To extract the actual value, you need to index the first element of that vector. +""" + +# ╔═╡ 980596b3-52e3-4f15-9678-f15f13f22688 +osol_log(50) # gives a vector with inside the value of W at time = 50 + +# ╔═╡ 851a7da1-ab7f-4a7a-ae1e-07a05c1e8052 +osol_log(50)[1] # This gives the value of W at time = 50 + +# ╔═╡ 8a9ec64f-2492-4f86-98d5-0e7c9512f5bc +x = [solution(50)[1] for solution in solutions]; + +# ╔═╡ 67f1348d-979c-42c4-a9d4-3031dbb5179a +md""" +Looking at the histogram, we observe that the yield is approximately normally distributed at 50 days, near the inflection point where logistic growth begins to plateau. This is intuitive, as both the carrying capacity $W_f$ and the growth rate μ are normally distributed, and at this point in the growth curve their combined influence on the yield remains approximately linear, preserving normality. Note, however, that this need not hold at other time points or for other models +""" + +# ╔═╡ d4b354f7-369a-497b-bc9d-68ab6cb5f2b4 +histogram(x, xlabel="W", ylabel="Count", title="Histogram of W at t=50") + +# ╔═╡ 6859b3f8-bd1c-47b7-a4e8-0c38501f5342 +p = mean(x) + +# ╔═╡ 63975327-5391-4de0-a538-4dcb503211b0 +md""" +What can you see on the resulting graph? What does this tell us about the uncertainty? +- missing +- missing +- missing +""" + # ╔═╡ 86f8d937-c42a-47ff-a0f6-6cae9d637fcf -md" -#### Study case 2 -" +md""" +### Case study 2 +""" # ╔═╡ 0755ea1c-50a0-4c8a-9adc-be2d63948a10 -md" -We initialize a vector `params2_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the second parameter**: -" +md""" +We initialize a vector `parms2_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the second parameter**: +""" # ╔═╡ 7befd387-f3a6-4b34-90b0-5f04b5429252 -params2_uncert_log = [:μ => 0.07, :Wf => 10.0±1.5] +parms2_uncert_log = [:μ => 0.07, :Wf => 10.0±1.5] # ╔═╡ 9e969176-f197-482f-9e2e-8b6d89de8f03 -md" +md""" We create the corresponding ODE problem and store it in `oprob2_uncert_log`: -" +""" # ╔═╡ 50cc4ee1-6601-442c-bac7-4d23b07db87e -oprob2_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, params2_uncert_log) +oprob2_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, parms2_uncert_log) # ╔═╡ 6d69f64e-8e3a-46e5-bb2d-bfb6f55e129e -md" +md""" We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol2_uncert_log`: -" +""" # ╔═╡ ee548df5-1334-47dc-9258-a25062de3f67 osol2_uncert_log = solve(oprob2_uncert_log, Tsit5(), saveat=2.0) # ╔═╡ afdb7450-e883-477b-b9a2-8e8585d16c40 -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ 25bfbf46-69b1-47eb-a473-788aad066d4a plot(osol2_uncert_log) # ╔═╡ b8896997-ddec-4d4a-a7c4-0e1008089d1b -md" -When thinking back of the sensitivity of $W$ to the parameter $W_f$, we saw that the corresponding sensitivity function was strongest in the tail of the curve around the steady state value. Looking at the above plot with error bars, we can see that de largest error bars (largest uncentrainty in the output variable) occurs at the tail of the curve where the sensitivity is strongest. -" +md""" +When thinking back to the sensitivity of $W$ to the parameter $W_f$, we saw that the corresponding sensitivity function was strongest in the tail of the curve around the steady state value. Looking at the above plot with error bars, we can see that the largest error bars (the largest uncertainty in the output variable) occur at the tail of the curve, where the sensitivity is strongest. +""" # ╔═╡ 5522a5c2-ff7e-4093-86d0-ffb879fcc6e8 md" -#### Study case 3 +### Case study 3 " # ╔═╡ 5e3da6ca-656e-41ad-ba04-bf3d6b51939d -md" -We initialize a vector `params_uncert_log` with the parameter values and uncertainty (standard deviation) **in all parameters**: -" +md""" +We initialize a vector `parms_uncert_log` with the parameter values and uncertainty (standard deviation) **in all parameters**: +""" # ╔═╡ 7be29c69-7fa0-4f41-80c9-f69aa01e8ea4 -params_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0±1.5] +parms_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0±1.5] # ╔═╡ d48c378e-df1c-45f7-9c3a-c5e9640e25a5 -md" +md""" We create the corresponding ODE problem and store it in `oprob_uncert_log`: -" +""" # ╔═╡ f209f4fb-985c-4a9b-b4dd-d50cc3fae1cf -oprob_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, params_uncert_log) +oprob_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, parms_uncert_log) # ╔═╡ de96646b-b752-4f72-ac78-2a503478bf50 -md" +md""" We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol_uncert_log`: -" +""" # ╔═╡ 6b1a3c2a-0a88-41ad-8866-41a33da5cd90 osol_uncert_log = solve(oprob_uncert_log, Tsit5(), saveat=2.0) # ╔═╡ 127ead82-2e2d-4c2e-9417-0490b6839514 -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ a5c254ad-b629-4a9e-89e9-90460e84620c plot(osol_uncert_log) # ╔═╡ c181d973-4fcc-4fbe-b37c-5f71dd33b0d2 -md" +md""" Now you have the combined effect of the uncertainty in the parameter $\mu$ as well as in the parameter $W_f$. -" +""" # ╔═╡ f4406d4e-d9f9-45fe-865c-405dbbb3127f -md" +md""" ## Exercises -" +""" # ╔═╡ e9f49fbc-d877-4db7-b320-a2db964ecaaa -md" +md""" ### Exercise 1 - Uncertainty analysis of the exponential growth model Perform an uncertainty analysis of the exponential growth model. Use the parameter uncertainties mentioned in the *Table* in the *Grass growth models* sections. -" +""" # ╔═╡ 5201995b-3298-4984-8839-dd25dc73a20f md" @@ -291,36 +393,38 @@ A possible *reaction network object* for the exponential growth model can be imp # ╔═╡ 97bdfce1-a670-469a-b8a9-34e1d1964409 growth_exp = @reaction_network begin + @species W(t)=2.0 + @parameters μ Wf μ*Wf, 0 --> W μ, W --> 0 end # ╔═╡ 52991cca-0365-48fb-b157-62a71c5573ee md" -The vector `u₀_exp` with the initial condition is: +The vector `u0_exp` with the initial condition is: " # ╔═╡ f7241453-755e-466c-a963-a504bc9ea446 u0_exp = [:W => 2.0] # ╔═╡ 11b81610-90b1-4d71-aeba-0bf33709dbda -md" -Initialize a vector `params_uncert_exp` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ -**Remark**: you can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. -" +md""" +Initialize a vector `parms_uncert_exp` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ +**Remark**: You can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. +""" # ╔═╡ 1af01af7-f5ee-4029-94ff-418c2f9748f2 -# params_uncert_exp = missing # Uncomment and complete the instruction -params_uncert_exp = [:μ => 0.02±0.01, :Wf => 10.0±1.5] +# parms_uncert_exp = missing +parms_uncert_exp = [:μ => 0.02±0.01, :Wf => 10.0±1.5] # ╔═╡ c25160d9-471a-414c-be99-10aa2efbb6d4 -md" +md""" Create the corresponding ODE problem and store it in `oprob_uncert_exp`: -" +""" # ╔═╡ 7fcec514-192b-4694-b149-14f5f87bae17 -# oprob_uncert_exp = missing # Uncomment and complete the instruction -oprob_uncert_exp = ODEProblem(growth_exp, u0_exp, tspan, params_uncert_exp) +# oprob_uncert_exp = missing +oprob_uncert_exp = ODEProblem(growth_exp, u0_exp, tspan, parms_uncert_exp) # ╔═╡ bd120dbc-827a-499a-a11c-25423ec5f397 md" @@ -328,31 +432,31 @@ Solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `os " # ╔═╡ 2dd58666-80ae-4346-ae15-34792eea2d25 -# osol_uncert_exp = missing # Uncomment and complete the instruction +# osol_uncert_exp = missing osol_uncert_exp = solve(oprob_uncert_exp, Tsit5(), saveat=2.0) # ╔═╡ 8059e00b-8e9c-4c11-a484-bea2a0bf155e -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ a89677f3-602b-4ef4-8f1a-ab203584f8df -# missing # Uncomment and complete the instruction +# missing plot(osol_uncert_exp) # ╔═╡ 2a470290-bca7-459f-b4f8-a5ac9fc031ed -md" +md""" Draw your conclusions: - missing -" +""" # ╔═╡ ded3aa76-ba09-499e-8c2d-52ff81ccfe02 -md" +md""" ### Exercise 2 - Uncertainty analysis of the exponential Gompertz model Perform an uncertainty analysis of the Gompertz growth model. Use the parameter uncertainties mentioned in the *Table* in the *Grass growth models* sections. -" +""" # ╔═╡ e245e5dd-bb11-455b-8de3-ae0ec1763399 md" @@ -361,37 +465,37 @@ A possible *reaction network object* for the Gompertz growth model can be implem # ╔═╡ 3e16bd3a-898a-49f9-9ed1-8998be7dc045 growth_gom = @reaction_network begin - # -μ, W --> ∅ - # D*log(W), W --> ∅ - μ-D*log(W), W --> 2W + @species W(t)=2.0 + @parameters μ d + μ-d*log(W), W --> 2W end # ╔═╡ db2f002f-20b7-4b5c-8d59-85b089fb3e59 md" -The vector `u₀_gom` with the initial condition is: +The vector `u0_gom` with the initial condition is: " # ╔═╡ 22a0013b-04f9-4148-9344-10bee1def0cf u0_gom = [:W => 2.0] # ╔═╡ f4fff573-6533-4ac3-897d-840a9aaf55f1 -md" -Initialize a vector `params_uncert_gom` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ -**Remark**: you can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. -" +md""" +Initialize a vector `parms_uncert_gom` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ +**Remark**: You can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. +""" # ╔═╡ 5ab3aca2-e570-4802-815b-5614680b427d -# params_uncert_gom = missing # Uncomment and complete the instruction -params_uncert_gom = [:μ => 0.09±0.01, :D => 0.040±0.002] +# parms_uncert_gom = missing +parms_uncert_gom = [:μ => 0.09±0.01, :d => 0.040±0.002] # ╔═╡ 890fe099-28ae-4e4f-aa1e-8d9ab6240130 -md" +md""" Create the corresponding ODE problem and store it in `oprob_uncert_log`: -" +""" # ╔═╡ 1aed2273-7b3e-4a49-a8d0-91251c58b700 -# oprob_uncert_gom = missing # Uncomment and complete the instruction -oprob_uncert_gom = ODEProblem(growth_gom, u0_gom, tspan, params_uncert_gom) +# oprob_uncert_gom = missing +oprob_uncert_gom = ODEProblem(growth_gom, u0_gom, tspan, parms_uncert_gom) # ╔═╡ 76da6edd-9e9f-4f5b-8673-b737c8f9c5f8 md" @@ -399,33 +503,34 @@ Solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `os " # ╔═╡ 8b10af49-e3d0-426e-9147-7666262d1c54 -# osol_uncert_gom = missing # Uncomment and complete the instruction +# osol_uncert_gom = missing osol_uncert_gom = solve(oprob_uncert_gom, Tsit5(), saveat=2.0) # ╔═╡ 6e8318bb-c884-4726-bb3d-9e0cfb7f40c1 -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ 92785455-14d2-4adb-9731-8575527c5217 -# missing # Uncomment and complete the instruction +# missing plot(osol_uncert_gom) # ╔═╡ 9d11df31-57b2-4825-a9db-34b422d5f084 -md" +md""" Draw your conclusions: - missing -" +""" # ╔═╡ Cell order: +# ╟─241a8a65-c59f-44f1-be39-d5edd1321b49 # ╠═d2c4d230-0943-11ef-3aad-5719e74bb20e -# ╠═47790920-74d2-4f21-b0f2-45dc19d25156 # ╠═8349306c-e98c-4221-9a1b-322fde3e18cb # ╠═83f3d978-aefa-40c1-a647-b26e837aeed6 -# ╠═73d5f921-82e4-4160-a3c4-fc1935f4c58d +# ╠═5185d0eb-7392-4775-9336-3e0f9e1449ce +# ╠═1c8ea6fe-650b-4c52-ba14-47efe3bd3e39 +# ╠═3eb3e651-afcf-41f8-a652-154a4f7ae07f # ╟─1f5e389e-e003-4e73-8f18-b5ad6340a912 -# ╟─241a8a65-c59f-44f1-be39-d5edd1321b49 # ╟─0f518e36-bc96-4599-807e-728504e5ca7b # ╟─da3cdc89-b911-4af9-9d4a-526301cba581 # ╟─78afaded-5a19-4386-aa76-7974977ea354 @@ -433,13 +538,13 @@ Draw your conclusions: # ╟─e5ab4490-6fd4-4b51-bfa7-1366438efefc # ╟─3e2750e7-220a-47b2-b445-eb3315734dec # ╟─c874a08c-7d82-4631-b611-c598dcaded09 +# ╟─87272d30-d95a-4b36-84ef-84bb9363c19e # ╟─9d3192bc-ee69-44a1-9273-06b8a55c60ef # ╟─8acb7ea2-54ef-428b-bb4f-364377d2c38d -# ╟─460d98ef-2177-49bf-87a4-35412f4183ed +# ╠═460d98ef-2177-49bf-87a4-35412f4183ed # ╟─38e0f746-1d04-46ba-81cc-21522b9ce6a4 # ╠═2538bd6e-027f-47fe-813c-db0a02586d2c # ╟─72b53ba6-dd05-4f2d-a65e-52783e76a1e9 -# ╠═5185d0eb-7392-4775-9336-3e0f9e1449ce # ╟─1a67db96-6973-446d-aa5b-253dd0008a73 # ╠═048112b0-10f1-47c5-834f-c3851d3078e8 # ╟─037cbf2b-dc0c-4b62-879b-aa702c964877 @@ -452,9 +557,31 @@ Draw your conclusions: # ╠═2fc91add-8468-4402-bcd8-42da5544f612 # ╟─5fa951d3-23ee-4510-873c-8158df4f2faf # ╠═d54635e8-1693-42bc-91f4-7e504a8661c7 -# ╠═707c6023-49ae-4a13-a220-d827a76a352c +# ╟─707c6023-49ae-4a13-a220-d827a76a352c # ╠═48e92497-f6b3-484c-bc5f-3ce69f4f1a91 # ╟─9b30b23b-6e92-49cf-bf93-e12a187fa2ba +# ╟─5dcb1ef0-d2a8-41fa-9592-2a86cf8364e6 +# ╟─0d883d09-1c92-4bb5-9bbf-3c87c856a5fa +# ╠═0e0dc268-faaf-462e-a980-9398c6015e02 +# ╟─e34f4b23-ac1b-4b88-b959-65085cef4b4f +# ╠═738448dd-d1c9-4fdc-af5d-fad4d91adf38 +# ╠═51579212-dd33-42b9-a19e-221ab7630c1a +# ╠═22bf3f50-1909-4a2c-820f-45b87f6bce14 +# ╟─616ee1d4-1f56-4c4a-a658-9f9351949f59 +# ╟─6b4a862e-6fb6-40ba-9e3f-43d57f169994 +# ╠═62afc495-29c3-4baa-b45a-94b84149d15d +# ╠═691a8b68-935e-4012-aa25-9cfdc323b1de +# ╟─f49a0085-396f-469e-afb8-2b7884f9bb1b +# ╠═ff5fcd0a-a37a-46d5-8b4b-f3f7a9a457a3 +# ╟─c0e72c96-23b6-4903-bd31-2978f6aebb47 +# ╟─988e9cd3-29a8-4fb1-a706-82d8356b9f20 +# ╠═980596b3-52e3-4f15-9678-f15f13f22688 +# ╠═851a7da1-ab7f-4a7a-ae1e-07a05c1e8052 +# ╠═8a9ec64f-2492-4f86-98d5-0e7c9512f5bc +# ╟─67f1348d-979c-42c4-a9d4-3031dbb5179a +# ╠═d4b354f7-369a-497b-bc9d-68ab6cb5f2b4 +# ╠═6859b3f8-bd1c-47b7-a4e8-0c38501f5342 +# ╟─63975327-5391-4de0-a538-4dcb503211b0 # ╟─86f8d937-c42a-47ff-a0f6-6cae9d637fcf # ╟─0755ea1c-50a0-4c8a-9adc-be2d63948a10 # ╠═7befd387-f3a6-4b34-90b0-5f04b5429252 @@ -482,14 +609,14 @@ Draw your conclusions: # ╟─52991cca-0365-48fb-b157-62a71c5573ee # ╠═f7241453-755e-466c-a963-a504bc9ea446 # ╟─11b81610-90b1-4d71-aeba-0bf33709dbda -# ╟─1af01af7-f5ee-4029-94ff-418c2f9748f2 +# ╠═1af01af7-f5ee-4029-94ff-418c2f9748f2 # ╟─c25160d9-471a-414c-be99-10aa2efbb6d4 # ╠═7fcec514-192b-4694-b149-14f5f87bae17 # ╟─bd120dbc-827a-499a-a11c-25423ec5f397 # ╠═2dd58666-80ae-4346-ae15-34792eea2d25 -# ╠═8059e00b-8e9c-4c11-a484-bea2a0bf155e +# ╟─8059e00b-8e9c-4c11-a484-bea2a0bf155e # ╠═a89677f3-602b-4ef4-8f1a-ab203584f8df -# ╟─2a470290-bca7-459f-b4f8-a5ac9fc031ed +# ╠═2a470290-bca7-459f-b4f8-a5ac9fc031ed # ╟─ded3aa76-ba09-499e-8c2d-52ff81ccfe02 # ╟─e245e5dd-bb11-455b-8de3-ae0ec1763399 # ╠═3e16bd3a-898a-49f9-9ed1-8998be7dc045 @@ -497,10 +624,10 @@ Draw your conclusions: # ╠═22a0013b-04f9-4148-9344-10bee1def0cf # ╟─f4fff573-6533-4ac3-897d-840a9aaf55f1 # ╠═5ab3aca2-e570-4802-815b-5614680b427d -# ╠═890fe099-28ae-4e4f-aa1e-8d9ab6240130 +# ╟─890fe099-28ae-4e4f-aa1e-8d9ab6240130 # ╠═1aed2273-7b3e-4a49-a8d0-91251c58b700 # ╟─76da6edd-9e9f-4f5b-8673-b737c8f9c5f8 # ╠═8b10af49-e3d0-426e-9147-7666262d1c54 # ╟─6e8318bb-c884-4726-bb3d-9e0cfb7f40c1 # ╠═92785455-14d2-4adb-9731-8575527c5217 -# ╟─9d11df31-57b2-4825-a9db-34b422d5f084 +# ╠═9d11df31-57b2-4825-a9db-34b422d5f084 diff --git a/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl b/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl index a61a520..34d954b 100644 --- a/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl +++ b/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl @@ -1,45 +1,46 @@ ### A Pluto.jl notebook ### -# v0.20.4 +# v0.20.21 using Markdown using InteractiveUtils # ╔═╡ 3ef93246-657d-4e77-9bf0-8380c64bfcfd -begin - # add this cell if you want the notebook to use the environment from where the Pluto server is launched - using Pkg - Pkg.activate("..") -end +using Pkg; Pkg.activate("..") # ╔═╡ 55cdebd2-0881-11ef-2722-91de1447877a -using Markdown - -# ╔═╡ 03ae0690-06a0-4276-9f00-d07b206fe124 -using InteractiveUtils +using Markdown, InteractiveUtils # ╔═╡ a355b0ba-baaf-49f4-a5dc-965364a884f0 -using Catalyst +using Catalyst, OrdinaryDiffEq # ╔═╡ 00fd6d49-f561-42e9-9413-d33af92f83dc -using OrdinaryDiffEq, StatsPlots +using StatsPlots, PlutoUI; TableOfContents() # ╔═╡ 7ae714c4-d25d-4f9f-ab3d-cc067db9c156 using ForwardDiff +# ╔═╡ dfce0717-e33e-4b6b-bb5c-d8754a74aba4 +using Turing + # ╔═╡ 31d294d1-3a1f-41db-abff-54f2a67c7ed9 md""" -### Exercise: Fermenter - Monod kinetics - Sensitivity analysis +# Exercise: Fermenter - Monod kinetics - Sensitivity analysis """ # ╔═╡ 5ffe7dcb-620d-4f22-95fe-2f77cda6fbe7 md""" -In one of the previous practicals we were introduced to a fermenter in which biomass $X$ [$g/L$] grows by breaking down substrate $S$ [$g/L$]. The reactor is fed with an inlet flow rate $Q_{in}$ [$L/h$], which consists of a (manipulable) input concentration of substrate $S_{in}$ [$g/L$]. This process was modelled using Monod kinetics, resulting in the model below: +In one of the previous practica, we were introduced to a fermenter in which biomass $X$ [$g/L$] grows by breaking down substrate $S$ [$g/L$]. The reactor is fed with an inlet flow rate $Q_{in}$ [$L/h$], which consists of a (manipulable) input concentration of substrate $S_{in}$ [$g/L$]. This process was modelled using Monod kinetics, resulting in the model below: $$\begin{eqnarray*} -S + X \xrightarrow[\quad\quad]{k} (1 + Y) \, X \quad\quad\quad\quad \textrm{with} \quad k = \cfrac{\mu_{max}}{S + K_s} +S + X \xrightarrow[\quad\quad]{k} (1 + Y) \, X \quad\quad\quad\quad \textrm{with} \quad k = \cfrac{\mu_{max}}{S + K_s} \, . \end{eqnarray*}$$ """ +# ╔═╡ d635d577-4d6d-40d7-a84c-3871981d59a4 +md""" +Suppose that at $t=0$ no substrate $S$ is present in the reactor but that there is initially some biomass with a concentration of $0.0005\;g/L$. The (default) parameter values are $Q = 2.0$, $V = 40.0$ and $Y = 0.67$. The parameters $\mu_{max}$, $K_s$ and $S_{in}$ will be assigned later. +""" + # ╔═╡ 6ec6da23-853b-4129-94cf-67b5cadb1f95 md""" The *reaction network object* model for this problem could be defined as: @@ -47,12 +48,13 @@ The *reaction network object* model for this problem could be defined as: # ╔═╡ 935ca610-7a7a-4692-8908-fc26abb880b4 fermenter_monod = @reaction_network begin - μmax/(S+Ks), S + X --> (1 + Y)*X - # Alternatives: - # mm(S, μmax, Ks)*X, S => Y*X - # mm(S, μmax, Ks)*X, S + X => (1 + Y)*X - Q/V, (S, X) --> 0 - Q/V*Sin, 0 --> S + @species S(t)=0.0 X(t)=0.0005 + @parameters Q=2.0 V=40.0 Y=0.67 Sin μmax Ks + μmax/(S + Ks), S + X --> (1 + Y)*X + # Alternative: + # mm(S, μmax, Ks)*X, S => Y*X + Q/V, (S, X) --> ∅ + Q/V*Sin, ∅ --> S end # ╔═╡ 79e6056a-881c-442f-8989-5bc284d3d777 @@ -72,96 +74,122 @@ $$\begin{eqnarray*} md""" Convert the system to a symbolic differential equation model and verify, by analyzing the differential equation, that your model is correctly implemented. -In case you want to use the `mm` function, keep in mind that `mm(S, μmax, Ks)` stands for $\mu_{max} \, \cfrac{S}{S + K_s}$. +Keep in mind that `mm(S, μmax, Ks)` stands for $\mu_{max} \, \cfrac{S}{S + K_s}$. """ # ╔═╡ 7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd -# osys = missing # Uncomment and complete the instruction +# osys = missing + +# ╔═╡ 1ee7359c-2757-44cd-8650-09b762cf30d0 +md""" +## Goals of this exercise +""" # ╔═╡ 55f1d688-0c53-481b-9965-5e92ca87ad83 md""" -The parameter values are $\mu_{max} = 0.40$, $K_s = 0.015$, $Y = 0.67$, $Q = 2.0$, $V = 40.0$ and $S_{in} = 0.022\;g/L$. Suppose that at $t=0$ no substrate $S$ is present in the reactor but that there is initially some biomass with a concetration of $0.0005\;g/L$.\ -Compute the following in a timespan of $[0, 100]\,h$: +Compute the following in a timespan of `(0.0, 100.0)`$\,h$: -- The sensitivities of $S$ and $X$ wrt. $\mu_{max}$, $K_s$ and $S_{in}$. +- The sensitivities of $S$ and $X$ w.r.t. $\mu_{max}$, $K_s$ and $S_{in}$. Plot the following: -- A figure with the sensitivity functions of $S$ and $X$ wrt. $S_{in}$. -- A figure with the sensitivity functions of $S$ wrt. $\mu_{max}$, $K_s$ and $S_{in}$. -- A figure with the sensitivity functions of $X$ wrt. $\mu_{max}$, $K_s$ and $S_{in}$. +- A figure with the sensitivity functions of $S$ and $X$ w.r.t. $S_{in}$. +- A figure with the sensitivity functions of $S$ w.r.t. $\mu_{max}$, $K_s$ and $S_{in}$. +- A figure with the sensitivity functions of $X$ w.r.t. $\mu_{max}$, $K_s$ and $S_{in}$. Interpret your results. + +Use the following parameter values $\mu_{max}$ = `0.40`, $K_s$ = `0.015` and $S_{in}$ = `0.022` for the calculations of the sensitivities. """ # ╔═╡ 6c4e3c09-4b84-4f5c-8739-2ac18e6f2af6 md""" -Initialize a vector `u0` with the initial conditions, set the timespan and initialize a vector `param` with the parameter values: +Initialize a vector `u0` with the initial conditions, set the timespan and initialize a vector `parms` with the parameter values: """ # ╔═╡ 2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c -# u0 = missing # Uncomment and complete the instruction +# u0 = missing # in principle not necessary since we will use the default # ╔═╡ 3fdc6b17-cdeb-4dc5-8886-9d3a62caac8d -# tspan = missing # Uncomment and complete the instruction +# tspan = missing -# ╔═╡ 0139da85-02e3-4021-9b39-84af7e68d428 +# ╔═╡ 201dfb54-2056-4846-a64c-ff4951dc084d md""" -For the sake of clarity, we will use the variables `μmax`, `Ks` and `Sin` to store the parameter values that are used for the calculation of the sensitivity functions. +For practical reasons, we will define the time step size as `dt = 0.5`. """ -# ╔═╡ 0f995929-4d2b-4a7a-8da1-04e4d501385f -# μmax = missing # Uncomment and complete the instruction +# ╔═╡ 5627068f-6442-43bb-bacb-bcba917372f3 +# dt = missing -# ╔═╡ 262e8346-df6d-49bf-9186-92f5afb421e0 -# Ks = missing # Uncomment and complete the instruction +# ╔═╡ 0139da85-02e3-4021-9b39-84af7e68d428 +md""" +For practical reasons, we will use the variables `μmax_val`, `Ks_val`, and `Sin_val` to store the parameter values that are used for the calculation of the sensitivity functions. +""" -# ╔═╡ baa777d2-abb8-45f8-87aa-b3b17c8dc07c -# Sin = missing # Uncomment and complete the instruction +# ╔═╡ 0f995929-4d2b-4a7a-8da1-04e4d501385f +begin + # μmax_val = missing + # Ks_val = missing + # Sin_val = missing +end; # ╔═╡ 79b0eb65-5a0f-40b3-aa97-4088421c562e -# params = missing # Uncomment and complete the instruction +# parms = missing # no need to include Q, V, and Y since we will use the default + +# ╔═╡ af882cf4-51fd-45ff-9b05-cfd52d6467b1 +md""" +## Preliminary simulation +""" # ╔═╡ f0f4fa14-6f99-4f21-a743-be61e08444a7 md""" -Create the ODE problem and store it in `oprob`. Next, solve the ODE problem using `Tsit5()` and `saveat=0.5`, and store the solution in `osol`. Finally plot the results. +Create the ODE problem and store it in `oprob`. Next, solve the ODE problem using `Tsit5()` and `saveat=dt`, and store the solution in `osol`. Finally plot the results. """ # ╔═╡ 8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 -# oprob = missing; # Uncomment and complete the instruction +# oprob = missing # ╔═╡ 89a31c32-88a4-479f-a688-ffcb75ee8e91 -# osol = missing # Uncomment and complete the instruction +# osol = missing # ╔═╡ 51a9b7e6-8ad9-477d-9596-ffd614df2c79 -# missing # Uncomment and complete the instruction +# missing + +# ╔═╡ 570ebda9-b187-42ab-a761-bed0fa3ce097 +md""" +## Local Sensitivity Analysis (LSA) +""" + +# ╔═╡ 0343674d-32d1-49e1-a7c1-47196ae760ed +md""" +### Setting up function +""" # ╔═╡ 693844d0-3858-4861-bae0-b47e78809f17 md""" -Write a solution function with as argument a vector of the parameters (with the values for which we want to calculate the sensitivity), and that returns the outputs. +Write a solution function with as argument a vector of the parameters (that you want the sensitivity on), and that returns the outputs. """ # ╔═╡ 9622f7ca-f71a-4ad9-a309-d7d10a1c3e3b -# Uncomment and complete the instruction -# function fermenter_monod_sim(params) - # μmax, Ks, Sin = missing - # u0 = missing - # tspan = missing - # params = missing - # oprob = missing; - # osol = missing - # return missing +# function fermenter_monod_sim(parms) +# missing +# ... # end # ╔═╡ 4a5971b1-f4d0-43b6-805f-e17f5052ae92 md""" -Make two functions based on the solution function, where each returns a single output; hence, one function that returns the output $S$, and another function that returns the output $X$. +Make two functions based on the solution function that each returns a single output, hence, one function that returns the output $S$, and another function that returns the output $X$. """ # ╔═╡ f40c6402-3c28-4a7d-b629-83507a9f29bd -# fermenter_monod_sim_S(params) = missing # Uncomment and complete the instruction +# fermenter_monod_sim_S(parms) = missing # ╔═╡ 3ae5bd00-2e06-4789-aab3-d897824d5e29 -# fermenter_monod_sim_X(params) = missing # Uncomment and complete the instruction +# fermenter_monod_sim_X(params) = missing + +# ╔═╡ 6913fb1b-1986-4b37-819e-a6a7bd91f1a5 +md""" +### Compute the outputs +""" # ╔═╡ 4bd2bcca-9c42-4333-b062-2aaa9f7be3fe md""" @@ -169,7 +197,7 @@ Make the time vector. """ # ╔═╡ fbd98975-aa32-46ae-8db0-0e65cdf48309 -# t_vals = missing # Uncomment and complete the instruction +# t_vals = missing # ╔═╡ 93791eb3-1eaa-4146-90b5-c4811fb3485b md""" @@ -177,21 +205,31 @@ Compute the two outputs $S$ and $X$ for the given parameter values. """ # ╔═╡ dc0557d6-81b9-4759-8ed7-3129f60c6dc3 -# S_sim = missing # Uncomment and complete the instruction +# S_sim = missing # ╔═╡ 95bc683c-f6e6-4b42-b90b-b5a863edd4d5 -# X_sim = missing # Uncomment and complete the instruction +# X_sim = missing + +# ╔═╡ 35717c5a-9895-4e62-a890-03ef04aa07b9 +md""" +### Compute the sensitivities +""" # ╔═╡ fa970c0e-fb3b-486f-bbc1-345d44f8f0da md""" Using `ForwardDiff.jacobian` to compute the sensitivities for the single ouputs $S$ and $X$. Hence, you need to call `ForwardDiff.jacobian` twice. """ +# ╔═╡ 07f2b8f8-5162-4077-a8a5-68fdec646644 +md""" +#### Absolute sensitivities +""" + # ╔═╡ 64354302-f4cc-4592-9302-5db0f5bccb2e -# sens_S = missing # Uncomment and complete the instruction +# sens_S = missing # ╔═╡ 49a94b9a-a543-495e-b4f1-c8579e59304d -# sens_X = missing # Uncomment and complete the instruction +# sens_X = missing # ╔═╡ 9cace6c1-e678-4dd7-8705-92a55eb32fa9 md""" @@ -199,20 +237,23 @@ Extract the (absolute) sensitivities of the outputs on the different parameters. """ # ╔═╡ a6dc2b60-6a0a-4140-892e-02cde8dc79d3 -# Uncomment and complete the instruction # begin # sens_S_on_μmax = missing # sens_S_on_Ks = missing # sens_S_on_Sin = missing -# end; +# end # ╔═╡ f806c243-9032-46b7-add3-4714344691c7 -# Uncomment and complete the instruction # begin # sens_X_on_μmax = missing # sens_X_on_Ks = missing # sens_X_on_Sin = missing -# end; +# end + +# ╔═╡ f4d0aaaa-bb3f-4335-901a-dbbcde643c78 +md""" +#### Normalized sensitivities +""" # ╔═╡ 5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 md""" @@ -220,7 +261,6 @@ Compute the normalized sensitivities. """ # ╔═╡ 76846731-929c-408f-a3de-970581c497e9 -# Uncomment and complete the instruction # begin # sens_S_on_μmax_rel = missing # sens_S_on_Ks_rel = missing @@ -228,197 +268,399 @@ Compute the normalized sensitivities. # end # ╔═╡ b6c57444-547c-4e82-8526-6a30566e07c5 -# Uncomment and complete the instruction # begin # sens_X_on_μmax_rel = missing # sens_X_on_Ks_rel = missing # sens_X_on_Sin_rel = missing -# end; +# end + +# ╔═╡ 7f0b9e63-2856-4601-9f14-c20c6b1b3707 +md""" +### Plotting + questions +""" # ╔═╡ 5388c2a7-5a11-4da8-be09-46045cde8a4e md""" -Plot the sensitivity functions of $S$ and $X$ wrt. $S_{in}$. Provide a suitable title (`title="..."`), labels (`label=["..." "..."]`) and an x-label (`xlabel="..."`), and set the line width to 2 (`linewidth=...`). +Plot the sensitivity functions of $S$ and $X$ on $S_{in}$. Provide a suitable title (`title="..."`), labels (`label=["..." "..."]`) and an x-label (`xlabel="..."`), and set the line width to 2 (`linewidth=...`). """ # ╔═╡ db840c76-a6c6-49fb-a0bb-d9149f947bc0 -# missing # Uncomment and complete the instruction - -# ╔═╡ 02314e9f-6186-4c97-9b28-03a4a1a15668 -maximum(sens_X_on_Sin_rel) +# missing # ╔═╡ d41375ef-6958-4705-a417-4c6a491232ee md""" -Interpret your results. Try to answer the following question(s): +!!! questions + Interpret your results. Try to answer the following question(s): + - Which output variable $S$ or $X$ is most sensitive to $S_{in}$ in steady state? + - Why is the sensitivity function of $S$ on $S_{in}$ at first positive but then becomes zero? """ -# ╔═╡ f6be14f3-e0d0-41dc-bc7b-54174b1ea5ea +# ╔═╡ 39c36eb4-5f15-4577-9b84-116e29d7891d md""" -!!! question - Which output variable, $S$ or $X$, is most sensitive to $S_{in}$ in steady state? +Answers: +- missing +- missing """ -# ╔═╡ b2583ee2-412b-4a0a-a4d0-fe476e0510e5 -md"- Answer: missing" +# ╔═╡ be89600a-4927-4afc-9813-d8a70adb2852 +md""" +Plot the sensitivity functions of $S$ on $\mu_{max}$, $K_s$ and $S_{in}$. Provide a suitable title (`title="..."`), labels (`label=["..." "..." "..."]`) and an x-label (`xlabel="..."`), and set the line width to 2 (`linewidth=...`). +""" -# ╔═╡ dbbf694b-942a-4757-a255-f67b208ba03b +# ╔═╡ c0223da4-9959-48d0-b607-633b2e82986c +# missing + +# ╔═╡ ff86a29f-9308-473b-aa1c-dfd4af8179c7 md""" -!!! question - Why is the sensitivity function of $S$ wrt. $S_{in}$ at first positive but then becomes zero? +!!! questions + Interpret your results. Try to answer the following question(s): + - Which parameter $\mu_{max}$, $K_s$ or $S_{in}$ affects the output $S$ the most in steady state? + - Why is the sensitivity function of $S$ w.r.t. $K_s$ positive? + - Why is the sensitivity function of $S$ w.r.t. $\mu_{max}$ negative? """ -# ╔═╡ f10ac0f2-6b14-4805-b649-56c63f5b523a -md"- Answer: missing" -# ╔═╡ be89600a-4927-4afc-9813-d8a70adb2852 +# ╔═╡ 430494da-402e-409d-85fc-c87cdfa3b427 md""" -Plot the sensitivity functions of $S$ wrt. $\mu_{max}$, $K_s$ and $S_{in}$. Provide a suitable title (`title="..."`), labels (`label=["..." "..." "..."]`) and an x-label (`xlabel="..."`), and set the line width to 2 (`linewidth=...`). +Answers: +- missing +- missing """ +#= +- It seems like μmax is affecting S the most in steady state, and its influence is negative; hence, the larger μmax, the smaller S. +- The sensitivity function of S w.r.t. Ks is positive, because the larger Ks, the smaller the reaction rate r (S => Y*X). Hence, less X will be produced; so less S will be consumed. Remember that r = μmax*S*X/(S + Ks) and Ks is in the denominator. +- The sensitivity function of S w.r.t. μmax is negative, because the larger μmax, the larger the reaction rate r (S => Y*X). Hence, more X will be produced; so more S will be consumed. Remember that r = μmax*S*X/(S + Ks) and μmax is in the numerator. +=# -# ╔═╡ c0223da4-9959-48d0-b607-633b2e82986c -# missing # Uncomment and complete the instruction +# ╔═╡ 16a84fdb-8ce2-45b9-bfb7-7f4e1284a1d7 +md""" +Plot the sensitivity functions of $X$ on $\mu_{max}$, $K_s$ and $S_{in}$. Provide a suitable title (`title="..."`), labels (`label=["..." "..." "..."]`) and an x-label (`xlabel="..."`), and set the line width to 2 (`linewidth=...`). +""" -# ╔═╡ 260ede0e-2584-4fe9-ae92-69990ca8f854 +# ╔═╡ 53134149-0bf7-41c1-9b35-e5037744211f +# missing + +# ╔═╡ 355ca6a7-466b-4969-ab48-28e2257f9810 md""" -Interpret your results. Try to answer the following question(s): +!!! questions + Interpret your results. Try to answer the following question(s): + - Which parameter $\mu_{max}$, $K_s$ or $S_{in}$ affects the output $X$ the most in steady state? + - Why is the sensitivity function of $X$ w.r.t. $K_s$ negative? + - Why is the sensitivity function of $X$ w.r.t. $\mu_{max}$ positive? """ -# ╔═╡ ff86a29f-9308-473b-aa1c-dfd4af8179c7 +# ╔═╡ 9a18a8bd-27d3-4771-a844-f6b65c7c8918 md""" -!!! question - Which parameter, $\mu_{max}$, $K_s$ or $S_{in}$, affects the output $S$ the most in steady state? Why is this? +Answers: +- missing +- missing +- missing """ +#= +- It seems like Sin is affecting X the most in steady state, and its influence is positive; hence, the larger Sin, the larger X. +- The sensitivity function of X w.r.t. Ks is negative, because the larger Ks, the smaller the reaction rate r (S => Y*X). Hence, less X will be produced. Remember that r = μmax*S*X/(S + Ks) and Ks is in the denominator. +- The sensitivity function of X w.r.t. μmax is positive, because the larger μmax, the larger the reaction rate r (S => Y*X). Hence, more X will be produced. Remember that r = μmax*S*X/(S + Ks) and μmax is in the numerator. +=# -# ╔═╡ 555228b2-8075-49ac-a7ef-f51aa51dd95d -md"- Answer: missing" +# ╔═╡ 6536f76a-0f6d-4155-ab3b-38e2e3189886 +md""" +## Monte Carlo error propagation +""" -# ╔═╡ 6704e43b-30bc-4168-93db-fb853c9761fa +# ╔═╡ 9fd59148-6993-4598-af95-ac85b6a4bb45 md""" -!!! question - Why is the sensitivity function of $S$ wrt. $\mu_{max}$ negative? +Let's take a look at how the uncertainty propagates in your model by using Monte Carlo simulations. Based on some literature search you can assume that the parameter `Ks` is normally distributed around the original value with a standard deviation of 20%. """ -# ╔═╡ 4d2094d8-98b3-4ac9-9612-942b5cd18ce0 -md"- Answer: missing" +# ╔═╡ 5a559bc2-3d55-499c-93cf-4b0845ffc4b0 +md""" +Start with making a Turing model in which you implement the prior and return the solution of the solved problem. +""" -# ╔═╡ dc8ef7cb-ae97-401c-9a25-12439e4363c6 +# ╔═╡ ecc28b8e-0eaf-489f-8104-7c6383d4e01f md""" -!!! question - Why is the sensitivity function of $S$ wrt. $K_s$ positive? Does this correspond to the meaning of the half-saturation constant $K_s$ or substrate concentration that allows to achieve half the maximum growth rate? *In other words: do higher values of $K_s$ support growth at low substrate concentrations?* +!!! note + Instead of creating a new ODEProblem everytime, you can simply remake an old ODEProblem by using `new_prob = remake(oprob, p = [:Ks => Ks])`. """ -# ╔═╡ 66712fb2-50ea-41be-bcca-305319d158e9 -md"- Answer: missing" +# ╔═╡ 26ae811d-11fd-44d8-b328-79381fa7472c +# @model function monod_deviation() +# Ks_dev ~ missing # 20% standard deviation +# new_prob = missing +# return missing +# end -# ╔═╡ 16a84fdb-8ce2-45b9-bfb7-7f4e1284a1d7 +# ╔═╡ 62791c48-5e38-46fe-8336-3a21003f3fe3 md""" -Plot the sensitivity functions of $X$ on $\mu_{max}$, $K_s$ and $S_{in}$. Provide a suitable title (`title="..."`), labels (`label=["..." "..." "..."]`) and an x-label (`xlabel="..."`), and set the line width to 2 (`linewidth=...`). +Now get sample 100 solutions from your monod_deviation model. """ -# ╔═╡ 53134149-0bf7-41c1-9b35-e5037744211f +# ╔═╡ 9e049653-347a-4198-90fd-23f15c404d0c # missing -# ╔═╡ 0355bbf6-853d-4bd8-b32f-98fdbd2b761c +# ╔═╡ f11827e5-de69-4a0c-b385-2c6ea031bb79 md""" -Interpret your results. Try to answer the following question(s): +We can now visualise the results of our Monte Carlo simulation by looping through our solutions and plotting them. It is advised to put `label = false` (avoids 100 labels), `color=:gray` (neutral color) and `alpha=0.4` (makes it more transparent) and `linestyle=:dot`. """ -# ╔═╡ 355ca6a7-466b-4969-ab48-28e2257f9810 +# ╔═╡ 9d0f8147-3c82-4ace-a408-20e32b4d8ea0 +# begin +# p1 = plot(title="Monte Carlo of logistic growth") # make an empty plot +# for missing in missing +# missing +# end +# plot!(osol, label = "Original simulation", linewidth = 3) # Plots original result +# p1 +# end + +# ╔═╡ 30c270f1-cbd5-401b-92db-7123efb01db6 +md""" +As can be seen from the graph, the uncertainty is rather high, especially if we want to carefully monitor the concentration of biomass through time. Note that, since we have a sample of values for $S$ and $X$ at every time step, it is also perfectly possible here to quantify the uncertainty by calculating the standard deviation at every time step, but the graph already gives a clear indication. +""" + +# ╔═╡ 16098df3-35c0-4d06-a482-35137fd53b97 +md""" +Now create a vector x1 that stores the values of the biomass at the end of the simulations and make a histogram of it. Put `xlims = [0.010, 0.015]`. +""" + +# ╔═╡ de5736d4-1c2f-4b0c-a84b-6d8b07d29a0a +# missing + +# ╔═╡ 1b1d0997-b854-4445-8c19-e73db6784a9b +# missing + +# ╔═╡ 43aaf563-f799-43e4-9f4c-8a19be1586dd +md""" +### Determininig optimal measurements +""" + +# ╔═╡ f7f1f035-b48e-4108-bb84-0c00ab68a93a +md""" +As a process operator, you may want to reduce the uncertainty in biomass simulations. However, since taking measurements can be costly, it is important to perform them at the most informative time points. + +Local sensitivity analysis is a useful tool in this context, as it indicates when the model output is most sensitive to changes in parameter values. At these time points, accurate measurements provide the greatest insight into the true parameter values. + +Start off by identifying the time at which the biomass concentration is most sensitive to variations in `Ks`. +""" + +# ╔═╡ 8dbadbc6-acd9-42f5-abf4-81fe6aef7fb7 +md""" +!!! note + You can find the index of the maximum value in a vector with the `argmax()` function. For example `argmax([1, 2, 5, 4])` will give you 3. +""" + +# ╔═╡ ddfb530a-ba79-43c0-973a-4deb8cad1d33 +# begin +# t_star_idx = missing +# t_star = t_vals[t_star_idx] + +# plot(t_vals, abs.(sens_X_on_Ks_rel), +# title="Sensitivity of Ks wrt. X", xlabel="Time (h)", +# ylabel="normalized sensitivity") +# vline!([t_star], label="t* = $(round(t_star, digits=2)) h", linestyle=:dash, color=:red) +# end + +# ╔═╡ 8359ea0e-a916-48f4-a07a-c905c3efaf15 +begin + μmax_real = 0.42 + Ks_real = 0.0165 + Sin_real = 0.0207 + real_solution = solve(remake(oprob, p = [:μmax => μmax_real, :Ks => Ks_real, + :Sin => Sin_real])); +end; + +# ╔═╡ 9654dea2-e13e-4aa9-8a5e-c3c368c26e24 +md""" +We have defined a function `real_solution(t)` in this notebook that returns a hypothetical measurement (`S_meas, X_meas`) at the given timepoint `t`. Use this function to get a measurement at the optimal time (`t_star`). +""" + +# ╔═╡ 77eab4ec-5012-40dd-b702-9fb48982ff38 +# S_meas, X_meas = real_solution(t_star) + +# ╔═╡ 3e93a001-9401-4a1f-b283-c79fd0c43434 md""" -!!! question - Which parameter, $\mu_{max}$, $K_s$ or $S_{in}$, affects the output $X$ the most in steady state? +#### Conditioning on the measurements """ -# ╔═╡ 4e08baf1-eeb6-49be-9751-204dd9ae1103 -md"- Answer: missing" +# ╔═╡ 684f9163-ac84-49be-beb1-66d7bdc21da8 +md""" +Now we can use this measurement to calibrate our model. In this part we will use Monte Carlo Markov chain to show the decrease in uncertainty. + +We start by creating a Turing model, for which you can use the priors from the first part of this exercise. For the standard deviation `σ_X` you can assume a value of 0.0002 and a value of 0.00055 for `σ_S`. + +Make sure to truncate the domain of your parameters to a reasonable domain or you might run into issues with the solver. +""" + +# ╔═╡ 4fe68607-d67e-4e71-865a-b10a06e5908d +# @model function monod_meas() +# σ_X = missing +# σ_S = missing + +# Ks_dev ~ missing -# ╔═╡ 182e0c2b-9aa4-4514-96e8-b6a9f53b371b +# sol = missing + +# S_pred, X_pred = sol(t_meas) + +# X_meas ~ missing +# S_meas ~ missing +# end + +# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 +# monod_cond_model = missing + +# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 md""" -!!! question - Why is the sensitivity function of $X$ wrt. $\mu_{max}$ positive? How does this compare to substrate $S$? +### MCMC """ -# ╔═╡ 05ec320e-ac87-45af-904e-5d69639e1f27 -md"- Answer: missing" +# ╔═╡ b8cd84de-57e4-4d8b-92da-e7c5b6d627bd +# monod_chain = missing -# ╔═╡ 22e98c77-0e9f-444a-a2f0-b940940996b7 +# ╔═╡ e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 md""" -!!! question - Why is the sensitivity function of $X$ wrt. $K_s$ negative? Compare it to that of substrate $S$. +Plot the chain to check if it properly converged to the posterior distribution. """ -# ╔═╡ 91795413-c9a2-43cb-a7b0-d9e055e1cf94 -md"- Answer: missing" +# ╔═╡ 834c3500-7556-4999-b465-08a8a07c7f28 +# missing + +# ╔═╡ eeaf3255-95c8-4813-bdd7-38635e546606 +md""" +Now, get the sampled posterior solutions and afterwards plot the resulting simulations. +""" + +# ╔═╡ a66b6d91-e529-4c20-b5d6-91830eae1437 +# missing + +# ╔═╡ c4192ff2-fca0-4ab1-8ac1-ce9771bc69d3 +md""" +Plot the resulting Monte Carlo simulation with the original simulations. +""" + +# ╔═╡ 54f20cff-21fc-4169-9f47-5e670b7ec621 +# begin +# p = plot(ylabel = "concentration (g/l)", xlabel = "Time (hours)") +# for missing in missing +# missing +# end +# plot!(osol, label = ["Original simulation S" "Original simulation X"], linewidth = 3) + +# p +# end + +# ╔═╡ 988aab6d-7193-4f3d-bd08-c3c36c920f81 +md""" +Lastly, create a histogram of the value of X at the end of the simulations based on your posterior distribution of Ks. Put the `xlims = [0.010, 0.015]` to compare with the prior distribution. +""" + +# ╔═╡ 2f55a771-0d37-48e5-a9d9-24693bf1d177 +# missing + +# ╔═╡ 18059d0c-cf37-44a1-86a8-859ac66a66f0 +# missing # ╔═╡ Cell order: +# ╟─31d294d1-3a1f-41db-abff-54f2a67c7ed9 # ╠═55cdebd2-0881-11ef-2722-91de1447877a -# ╠═03ae0690-06a0-4276-9f00-d07b206fe124 # ╠═3ef93246-657d-4e77-9bf0-8380c64bfcfd # ╠═a355b0ba-baaf-49f4-a5dc-965364a884f0 # ╠═00fd6d49-f561-42e9-9413-d33af92f83dc # ╠═7ae714c4-d25d-4f9f-ab3d-cc067db9c156 -# ╟─31d294d1-3a1f-41db-abff-54f2a67c7ed9 +# ╠═dfce0717-e33e-4b6b-bb5c-d8754a74aba4 # ╟─5ffe7dcb-620d-4f22-95fe-2f77cda6fbe7 +# ╟─d635d577-4d6d-40d7-a84c-3871981d59a4 # ╟─6ec6da23-853b-4129-94cf-67b5cadb1f95 # ╠═935ca610-7a7a-4692-8908-fc26abb880b4 # ╟─79e6056a-881c-442f-8989-5bc284d3d777 # ╟─fa93e2c3-8b43-418e-ba24-406645b2e397 # ╟─06730f54-7293-43f2-b772-84eec3e5528a # ╠═7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd +# ╟─1ee7359c-2757-44cd-8650-09b762cf30d0 # ╟─55f1d688-0c53-481b-9965-5e92ca87ad83 # ╟─6c4e3c09-4b84-4f5c-8739-2ac18e6f2af6 # ╠═2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c # ╠═3fdc6b17-cdeb-4dc5-8886-9d3a62caac8d +# ╟─201dfb54-2056-4846-a64c-ff4951dc084d +# ╠═5627068f-6442-43bb-bacb-bcba917372f3 # ╟─0139da85-02e3-4021-9b39-84af7e68d428 # ╠═0f995929-4d2b-4a7a-8da1-04e4d501385f -# ╠═262e8346-df6d-49bf-9186-92f5afb421e0 -# ╠═baa777d2-abb8-45f8-87aa-b3b17c8dc07c # ╠═79b0eb65-5a0f-40b3-aa97-4088421c562e +# ╟─af882cf4-51fd-45ff-9b05-cfd52d6467b1 # ╟─f0f4fa14-6f99-4f21-a743-be61e08444a7 # ╠═8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 # ╠═89a31c32-88a4-479f-a688-ffcb75ee8e91 # ╠═51a9b7e6-8ad9-477d-9596-ffd614df2c79 +# ╟─570ebda9-b187-42ab-a761-bed0fa3ce097 +# ╟─0343674d-32d1-49e1-a7c1-47196ae760ed # ╟─693844d0-3858-4861-bae0-b47e78809f17 # ╠═9622f7ca-f71a-4ad9-a309-d7d10a1c3e3b # ╟─4a5971b1-f4d0-43b6-805f-e17f5052ae92 # ╠═f40c6402-3c28-4a7d-b629-83507a9f29bd # ╠═3ae5bd00-2e06-4789-aab3-d897824d5e29 +# ╟─6913fb1b-1986-4b37-819e-a6a7bd91f1a5 # ╟─4bd2bcca-9c42-4333-b062-2aaa9f7be3fe # ╠═fbd98975-aa32-46ae-8db0-0e65cdf48309 # ╟─93791eb3-1eaa-4146-90b5-c4811fb3485b # ╠═dc0557d6-81b9-4759-8ed7-3129f60c6dc3 # ╠═95bc683c-f6e6-4b42-b90b-b5a863edd4d5 +# ╟─35717c5a-9895-4e62-a890-03ef04aa07b9 # ╟─fa970c0e-fb3b-486f-bbc1-345d44f8f0da +# ╟─07f2b8f8-5162-4077-a8a5-68fdec646644 # ╠═64354302-f4cc-4592-9302-5db0f5bccb2e # ╠═49a94b9a-a543-495e-b4f1-c8579e59304d # ╟─9cace6c1-e678-4dd7-8705-92a55eb32fa9 # ╠═a6dc2b60-6a0a-4140-892e-02cde8dc79d3 # ╠═f806c243-9032-46b7-add3-4714344691c7 +# ╟─f4d0aaaa-bb3f-4335-901a-dbbcde643c78 # ╟─5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 # ╠═76846731-929c-408f-a3de-970581c497e9 # ╠═b6c57444-547c-4e82-8526-6a30566e07c5 +# ╟─7f0b9e63-2856-4601-9f14-c20c6b1b3707 # ╟─5388c2a7-5a11-4da8-be09-46045cde8a4e # ╠═db840c76-a6c6-49fb-a0bb-d9149f947bc0 -# ╠═02314e9f-6186-4c97-9b28-03a4a1a15668 # ╟─d41375ef-6958-4705-a417-4c6a491232ee -# ╟─f6be14f3-e0d0-41dc-bc7b-54174b1ea5ea -# ╠═b2583ee2-412b-4a0a-a4d0-fe476e0510e5 -# ╟─dbbf694b-942a-4757-a255-f67b208ba03b -# ╠═f10ac0f2-6b14-4805-b649-56c63f5b523a +# ╠═39c36eb4-5f15-4577-9b84-116e29d7891d # ╟─be89600a-4927-4afc-9813-d8a70adb2852 # ╠═c0223da4-9959-48d0-b607-633b2e82986c -# ╟─260ede0e-2584-4fe9-ae92-69990ca8f854 # ╟─ff86a29f-9308-473b-aa1c-dfd4af8179c7 -# ╠═555228b2-8075-49ac-a7ef-f51aa51dd95d -# ╟─6704e43b-30bc-4168-93db-fb853c9761fa -# ╠═4d2094d8-98b3-4ac9-9612-942b5cd18ce0 -# ╟─dc8ef7cb-ae97-401c-9a25-12439e4363c6 -# ╠═66712fb2-50ea-41be-bcca-305319d158e9 +# ╟─430494da-402e-409d-85fc-c87cdfa3b427 # ╟─16a84fdb-8ce2-45b9-bfb7-7f4e1284a1d7 # ╠═53134149-0bf7-41c1-9b35-e5037744211f -# ╟─0355bbf6-853d-4bd8-b32f-98fdbd2b761c # ╟─355ca6a7-466b-4969-ab48-28e2257f9810 -# ╠═4e08baf1-eeb6-49be-9751-204dd9ae1103 -# ╟─182e0c2b-9aa4-4514-96e8-b6a9f53b371b -# ╠═05ec320e-ac87-45af-904e-5d69639e1f27 -# ╟─22e98c77-0e9f-444a-a2f0-b940940996b7 -# ╠═91795413-c9a2-43cb-a7b0-d9e055e1cf94 +# ╟─9a18a8bd-27d3-4771-a844-f6b65c7c8918 +# ╟─6536f76a-0f6d-4155-ab3b-38e2e3189886 +# ╟─9fd59148-6993-4598-af95-ac85b6a4bb45 +# ╟─5a559bc2-3d55-499c-93cf-4b0845ffc4b0 +# ╟─ecc28b8e-0eaf-489f-8104-7c6383d4e01f +# ╠═26ae811d-11fd-44d8-b328-79381fa7472c +# ╟─62791c48-5e38-46fe-8336-3a21003f3fe3 +# ╠═9e049653-347a-4198-90fd-23f15c404d0c +# ╟─f11827e5-de69-4a0c-b385-2c6ea031bb79 +# ╠═9d0f8147-3c82-4ace-a408-20e32b4d8ea0 +# ╟─30c270f1-cbd5-401b-92db-7123efb01db6 +# ╟─16098df3-35c0-4d06-a482-35137fd53b97 +# ╠═de5736d4-1c2f-4b0c-a84b-6d8b07d29a0a +# ╠═1b1d0997-b854-4445-8c19-e73db6784a9b +# ╟─43aaf563-f799-43e4-9f4c-8a19be1586dd +# ╟─f7f1f035-b48e-4108-bb84-0c00ab68a93a +# ╟─8dbadbc6-acd9-42f5-abf4-81fe6aef7fb7 +# ╠═ddfb530a-ba79-43c0-973a-4deb8cad1d33 +# ╟─8359ea0e-a916-48f4-a07a-c905c3efaf15 +# ╟─9654dea2-e13e-4aa9-8a5e-c3c368c26e24 +# ╠═77eab4ec-5012-40dd-b702-9fb48982ff38 +# ╟─3e93a001-9401-4a1f-b283-c79fd0c43434 +# ╟─684f9163-ac84-49be-beb1-66d7bdc21da8 +# ╠═4fe68607-d67e-4e71-865a-b10a06e5908d +# ╠═11e105de-270c-434a-bfb4-9118dfc9e461 +# ╟─e46afb98-16de-4d19-86b4-b920b6b929c5 +# ╠═b8cd84de-57e4-4d8b-92da-e7c5b6d627bd +# ╟─e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +# ╠═834c3500-7556-4999-b465-08a8a07c7f28 +# ╟─eeaf3255-95c8-4813-bdd7-38635e546606 +# ╠═a66b6d91-e529-4c20-b5d6-91830eae1437 +# ╟─c4192ff2-fca0-4ab1-8ac1-ce9771bc69d3 +# ╠═54f20cff-21fc-4169-9f47-5e670b7ec621 +# ╟─988aab6d-7193-4f3d-bd08-c3c36c920f81 +# ╠═2f55a771-0d37-48e5-a9d9-24693bf1d177 +# ╠═18059d0c-cf37-44a1-86a8-859ac66a66f0 diff --git a/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl b/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl index b84645f..b5f13f8 100644 --- a/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl +++ b/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl @@ -5,45 +5,41 @@ using Markdown using InteractiveUtils # ╔═╡ 8349306c-e98c-4221-9a1b-322fde3e18cb -begin - # add this cell if you want the notebook to use the environment from where the Pluto server is launched - using Pkg - Pkg.activate("..") -end +using Pkg; Pkg.activate("..") # ╔═╡ d2c4d230-0943-11ef-3aad-5719e74bb20e -using Markdown - -# ╔═╡ 47790920-74d2-4f21-b0f2-45dc19d25156 -using InteractiveUtils +using Markdown, InteractiveUtils # ╔═╡ 83f3d978-aefa-40c1-a647-b26e837aeed6 -using Catalyst - -# ╔═╡ 73d5f921-82e4-4160-a3c4-fc1935f4c58d -using OrdinaryDiffEq, StatsPlots +using OrdinaryDiffEq, Catalyst # ╔═╡ 5185d0eb-7392-4775-9336-3e0f9e1449ce using Measurements -# ╔═╡ 1f5e389e-e003-4e73-8f18-b5ad6340a912 -Catalyst.ModelingToolkit.NaNMath.log(x::Measurement) = Measurements.log(x) - # using `log` on a Measurement-type variable in the reaction rate of a Catalyst model will otherwise error in the current Catalyst version +# ╔═╡ 1c8ea6fe-650b-4c52-ba14-47efe3bd3e39 +using StatsPlots, PlutoUI; TableOfContents() + +# ╔═╡ 3eb3e651-afcf-41f8-a652-154a4f7ae07f +using Turing # ╔═╡ 241a8a65-c59f-44f1-be39-d5edd1321b49 -md" +md""" # Introduction to uncertainty analysis -" +""" + +# ╔═╡ 1f5e389e-e003-4e73-8f18-b5ad6340a912 +Catalyst.ModelingToolkit.NaNMath.log(x::Measurement) = Measurements.log(x) + # using `log` on a Measurement-type variable in the reaction rate of a Catalyst model will otherwise error in the current Catalyst version # ╔═╡ 0f518e36-bc96-4599-807e-728504e5ca7b -md" +md""" ## Goal of this practicum -" +""" # ╔═╡ da3cdc89-b911-4af9-9d4a-526301cba581 -md" -Parameter uncertainty plays a crucial role in shaping the behavior of output variables. Model equations describe how systems evolve over time, often incorporating parameters representing various aspects of the system's characteristics. However, these parameters are rarely known with absolute certainty and they often come with inherent uncertainty due to measurement errors, variability in real-world conditions, or incomplete knowledge about the system. This uncertainty can propagate through the model equations, leading to uncertainties in the predicted outcomes. Consequently, understanding the influence of parameter uncertainty becomes essential for assessing the reliability and robustness of the model predictions, as well as for making informed decisions based on these predictions. Techniques such as sensitivity analysis and uncertainty quantification are employed to explore and quantify the impact of parameter uncertainty on the output variables, providing insights into the system's behavior and guiding the refinement of models for improved accuracy and reliability. -" +md""" +Parameter uncertainty plays a crucial role in shaping the behavior of output variables. Model equations describe how systems evolve, often incorporating parameters representing various aspects of the system's characteristics. However, these parameters are rarely known with absolute certainty and often carry inherent uncertainty due to measurement errors, variability in real-world conditions, or incomplete knowledge of the system. This uncertainty can propagate through the model equations, leading to uncertainties in the predicted outcomes. Consequently, understanding the influence of parameter uncertainty becomes essential for assessing the reliability and robustness of the model predictions, as well as for making informed decisions based on these predictions. Techniques such as sensitivity analysis and uncertainty quantification are employed to explore and quantify the impact of parameter uncertainty on the output variables, providing insights into the system's behavior and guiding the refinement of models for improved accuracy and reliability. +""" # ╔═╡ 78afaded-5a19-4386-aa76-7974977ea354 md" @@ -56,108 +52,111 @@ We will now compute the variability in the output variables replected as error b " # ╔═╡ e5ab4490-6fd4-4b51-bfa7-1366438efefc -md" +md""" ## Grass growth models -" +""" # ╔═╡ 3e2750e7-220a-47b2-b445-eb3315734dec -md" +md""" In this notebook, three different models will be used, each modelling the yield of grass in a grassland: - Logistic growth model: $\cfrac{dW}{dt} = \mu \left( 1 - \cfrac{W}{W_f} \right) W$ - Exponential growth model: $\cfrac{dW}{dt} = \mu \left( W_f - W \right)$ -- Gompertz growth model: $\cfrac{dW}{dt} = \left( \mu - D \ln(W) \right) W$ +- Gompertz growth model: $\cfrac{dW}{dt} = \left( \mu - d \ln(W) \right) W$ with output $W$ the grass yield, and $W_f$, $\mu$ and $D$ parameters. The table below show some typical values for the parameters together with their uncertainties: -| | $\mu$ | $W_f$ | $D$ | +| | $\mu$ | $W_f$ | $d$ | |:----------- |:----------:|:-----------:|:------------:| | Logistic | 0.07$\pm$0.02 | 10.0$\pm$0.15 | | | Exponential | 0.02$\pm$0.01 | 10.0$\pm$0.15 | | | Gompertz | 0.09$\pm$0.01 | | 0.040$\pm$0.002 | We will use an initial condition of $W_0 = 2.0$ for each and a simulation time of $100$ days. -" +""" # ╔═╡ c874a08c-7d82-4631-b611-c598dcaded09 -md" +md""" We will illustrate how to compute the error bars, due to model parameter uncertainty, in conjunction with the output variable simulation for the logistic model. The same will be left as exercises below for the exponential and gompertz models. **Important:** - We will use consequently `_log`, `_exp` and `_gom` appended to relevant variables names in order to indicate their model origin **and** to prevent cell-disabling that occurs when using the same variables names in these Notebooks. -" +""" -# ╔═╡ 9d3192bc-ee69-44a1-9273-06b8a55c60ef -md" +# ╔═╡ 87272d30-d95a-4b36-84ef-84bb9363c19e +md""" ### Uncertainty analysis of the logistic growth model +""" -We will start by modelling our system and simulating using the aforementioned parameters values and uncertainties, initial condition and timespan. -" +# ╔═╡ 9d3192bc-ee69-44a1-9273-06b8a55c60ef +md""" +We will start by modelling our system and simulating using the aforementioned parameter values and uncertainties, initial condition, and timespan. +""" # ╔═╡ 8acb7ea2-54ef-428b-bb4f-364377d2c38d -md" +md""" Implementation of the system: -" +""" # ╔═╡ 460d98ef-2177-49bf-87a4-35412f4183ed growth_mod_log = @reaction_network begin - #μ*W, 0 --> W - #μ/Wf*W, W --> 0 + @species W(t)=2.0 + @parameters μ Wf μ*(1-W/Wf), W --> 2W end # ╔═╡ 38e0f746-1d04-46ba-81cc-21522b9ce6a4 -md" +md""" Convert the *reaction model* to check that we work with the correct differential equation: -" +""" # ╔═╡ 2538bd6e-027f-47fe-813c-db0a02586d2c osys_log = convert(ODESystem, growth_mod_log) # ╔═╡ 72b53ba6-dd05-4f2d-a65e-52783e76a1e9 -md" -In order to use uncertainties in the parameter values, we will need to load the package `Measurements`: -" +md""" +In order to use uncertainties in the parameter values, we will need to load the package `Measurements` (see above). +""" # ╔═╡ 1a67db96-6973-446d-aa5b-253dd0008a73 -md" +md""" Setting initial conditions: -" +""" # ╔═╡ 048112b0-10f1-47c5-834f-c3851d3078e8 u0_log = [:W => 2.0] # ╔═╡ 037cbf2b-dc0c-4b62-879b-aa702c964877 -md" +md""" Set the timespan for the simulation: -" +""" # ╔═╡ 40c1192c-bc76-4d3a-8144-dd76ba72fba2 tspan = (0.0, 100.0) # ╔═╡ 0b5493eb-6ce6-43c1-9289-4c2e6f20d014 -md" -In order to see the effect of the individual parameter uncertainties on the output variable, we will analyse this (only for didactical reasons) assuming different study cases: +md""" +To see the effect of the individual parameter uncertainties on the output variable, we will analyse this (only for didactical reasons), assuming different study cases: -1. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ while $W_f$ is exactly known and equal to $0.15$. -2. $K_s$ has a nominal value of $10.0$ and a standard deviation of $0.15$ while $\mu$ is exactly known and equal to $0.07$. -3. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ and similarily $K_s$ has a nominal value of $10.0$ and a standard deviation of $0.15$. +1. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ while $W_f$ is exactly known and equal to $10.0$. +2. $W_f$ has a nominal value of $10.0$ and a standard deviation of $0.15$ while $\mu$ is exactly known and equal to $0.07$. +3. $\mu$ has a nominal value of $0.07$ and a standard deviation of $0.02$ and similarily $W_f$ has a nominal value of $10.0$ and a standard deviation of $0.15$. -The last study case is the true case because in reality all encertainties will contribute at once. -" +The last case study is the true case because, in reality, all uncertainties will contribute simultaneously. +""" # ╔═╡ b4c76f21-9839-4db5-b952-43869d3a1efe -md" -#### Study case 1 -" +md""" +### Case study 1 +""" # ╔═╡ 6ab73fc5-7315-4ec6-a3ee-9b5c689ce82e -md" -We initialize a vector `params1_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the first parameter**: -" +md""" +We initialize a vector `parms1_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the first parameter**: +""" # ╔═╡ b0f32f69-b6a2-4050-97e2-8018d07463d8 -params1_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0] +parms1_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0] # ╔═╡ c1a0702c-6fa5-47d7-bd72-6bfcbdeba93f md" @@ -165,20 +164,20 @@ We create the corresponding ODE problem and store it in `oprob1_uncert_log`: " # ╔═╡ 2fc91add-8468-4402-bcd8-42da5544f612 -oprob1_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, params1_uncert_log) +oprob1_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, parms1_uncert_log) # ╔═╡ 5fa951d3-23ee-4510-873c-8158df4f2faf -md" +md""" We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol1_uncert_log`: -" +""" # ╔═╡ d54635e8-1693-42bc-91f4-7e504a8661c7 osol1_uncert_log = solve(oprob1_uncert_log, Tsit5(), saveat=2.0) # ╔═╡ 707c6023-49ae-4a13-a220-d827a76a352c -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ 48e92497-f6b3-484c-bc5f-3ce69f4f1a91 plot(osol1_uncert_log) @@ -188,101 +187,204 @@ md" When thinking back of the sensitivity of $W$ to the parameter $\mu$, we saw that the corresponding sensitivity function had a maximum around $t=33\;s$. Looking at the above plot with error bars, we can see that de largest error bars (largest uncentrainty in the output variable) occurs at the timepoints where the sensitivity is strongest. " +# ╔═╡ 5dcb1ef0-d2a8-41fa-9592-2a86cf8364e6 +md""" +### Alternative: Monte Carlo uncertainty propagation +""" + +# ╔═╡ 0d883d09-1c92-4bb5-9bbf-3c87c856a5fa +md""" +In uncertainty propagation, we are interested in the effect of uncertainties (or errors) on the final output of the model. Monte Carlo simulations are a straightforward way for testing different parameters and checking the results of the outcomes. + +Below a Turing model is defined that samples values from a normal distribution with a 5% standard deviation. This is our Prior belief of how accurate the estimate of our intial parameters is. Remember that this is similar to what we did in practical 4. +""" + +# ╔═╡ 0e0dc268-faaf-462e-a980-9398c6015e02 +@model function logistic_deviation() + μ_dev ~ Normal(0.07, 0.0035) # 5% standard deviation + Wf_dev ~ Normal(10, 0.5) # 5% standard deviation + return solve(remake(oprob_log, p = [:Wf => Wf_dev, :μ => μ_dev]), saveat=0.5); +end + +# ╔═╡ e34f4b23-ac1b-4b88-b959-65085cef4b4f +md""" +Using the `sample()` function we can sample values from our distribution. +""" + +# ╔═╡ 738448dd-d1c9-4fdc-af5d-fad4d91adf38 +μ_model = logistic_deviation(); + +# ╔═╡ 51579212-dd33-42b9-a19e-221ab7630c1a +chain = sample(μ_model, Prior(), 500); + +# ╔═╡ 22bf3f50-1909-4a2c-820f-45b87f6bce14 +solutions = generated_quantities(μ_model, chain); + +# ╔═╡ 616ee1d4-1f56-4c4a-a658-9f9351949f59 +md""" +Finally, we can use theses perturbed values to model the effects of these small perturbations on the output of the model. +""" + +# ╔═╡ 6b4a862e-6fb6-40ba-9e3f-43d57f169994 +md""" +Firstly, solve the original system (Wf = 10 and μ = 0.07) and afterwards we can compare this with the outputs of the Monte Carlo simulations. +""" + +# ╔═╡ 62afc495-29c3-4baa-b45a-94b84149d15d +oprob_log = ODEProblem(growth_mod_log, u0_log,tspan, [:μ => 0.07, :Wf => 10]); + +# ╔═╡ 691a8b68-935e-4012-aa25-9cfdc323b1de +osol_log = solve(oprob_log, Tsit5(), saveat = 0.5); + +# ╔═╡ f49a0085-396f-469e-afb8-2b7884f9bb1b +md""" +Let's now plot the solution of the Monte Carlo experiments. +""" + +# ╔═╡ ff5fcd0a-a37a-46d5-8b4b-f3f7a9a457a3 +begin + p1 = plot(title="Monte Carlo of logistic growth") + for solution in solutions + plot!(solution, label = false, color =:gray, alpha = 0.4, + linestyle = :dot) + end + plot!(osol_log, label = "Original simulation", linewidth = 3) + p1 +end + +# ╔═╡ c0e72c96-23b6-4903-bd31-2978f6aebb47 +md""" +An advantage of using Monte Carlo over the error bars, is that the resulting graphs create a distribution from which we can calculate probabilities. We could now for example calculate what the average yield is after 50 days: +""" + +# ╔═╡ 988e9cd3-29a8-4fb1-a706-82d8356b9f20 +md""" +It is worth noting that you can obtain the solution at a specific time point by indexing it using parentheses `()`, as shown below. However, the result is still returned as a vector (indicated by the square brackets `[]`). To extract the actual value, you need to index the first element of that vector. +""" + +# ╔═╡ 980596b3-52e3-4f15-9678-f15f13f22688 +osol_log(50) # gives a vector with inside the value of W at time = 50 + +# ╔═╡ 851a7da1-ab7f-4a7a-ae1e-07a05c1e8052 +osol_log(50)[1] # This gives the value of W at time = 50 + +# ╔═╡ 8a9ec64f-2492-4f86-98d5-0e7c9512f5bc +x = [solution(50)[1] for solution in solutions]; + +# ╔═╡ 67f1348d-979c-42c4-a9d4-3031dbb5179a +md""" +Looking at the histogram, we observe that the yield is approximately normally distributed at 50 days, near the inflection point where logistic growth begins to plateau. This is intuitive, as both the carrying capacity $W_f$ and the growth rate μ are normally distributed, and at this point in the growth curve their combined influence on the yield remains approximately linear, preserving normality. Note, however, that this need not hold at other time points or for other models +""" + +# ╔═╡ d4b354f7-369a-497b-bc9d-68ab6cb5f2b4 +histogram(x, xlabel="W", ylabel="Count", title="Histogram of W at t=50") + +# ╔═╡ 6859b3f8-bd1c-47b7-a4e8-0c38501f5342 +p = mean(x) + +# ╔═╡ 63975327-5391-4de0-a538-4dcb503211b0 +md""" +What can you see on the resulting graph? What does this tell us about the uncertainty? +- missing +- missing +- missing +""" + # ╔═╡ 86f8d937-c42a-47ff-a0f6-6cae9d637fcf -md" -#### Study case 2 -" +md""" +### Case study 2 +""" # ╔═╡ 0755ea1c-50a0-4c8a-9adc-be2d63948a10 -md" -We initialize a vector `params2_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the second parameter**: -" +md""" +We initialize a vector `parms2_uncert_log` with the parameter values and uncertainty (standard deviation) **only in the second parameter**: +""" # ╔═╡ 7befd387-f3a6-4b34-90b0-5f04b5429252 -params2_uncert_log = [:μ => 0.07, :Wf => 10.0±1.5] +parms2_uncert_log = [:μ => 0.07, :Wf => 10.0±1.5] # ╔═╡ 9e969176-f197-482f-9e2e-8b6d89de8f03 -md" +md""" We create the corresponding ODE problem and store it in `oprob2_uncert_log`: -" +""" # ╔═╡ 50cc4ee1-6601-442c-bac7-4d23b07db87e -oprob2_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, params2_uncert_log) +oprob2_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, parms2_uncert_log) # ╔═╡ 6d69f64e-8e3a-46e5-bb2d-bfb6f55e129e -md" +md""" We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol2_uncert_log`: -" +""" # ╔═╡ ee548df5-1334-47dc-9258-a25062de3f67 osol2_uncert_log = solve(oprob2_uncert_log, Tsit5(), saveat=2.0) # ╔═╡ afdb7450-e883-477b-b9a2-8e8585d16c40 -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ 25bfbf46-69b1-47eb-a473-788aad066d4a plot(osol2_uncert_log) # ╔═╡ b8896997-ddec-4d4a-a7c4-0e1008089d1b -md" -When thinking back of the sensitivity of $W$ to the parameter $W_f$, we saw that the corresponding sensitivity function was strongest in the tail of the curve around the steady state value. Looking at the above plot with error bars, we can see that de largest error bars (largest uncentrainty in the output variable) occurs at the tail of the curve where the sensitivity is strongest. -" +md""" +When thinking back to the sensitivity of $W$ to the parameter $W_f$, we saw that the corresponding sensitivity function was strongest in the tail of the curve around the steady state value. Looking at the above plot with error bars, we can see that the largest error bars (the largest uncertainty in the output variable) occur at the tail of the curve, where the sensitivity is strongest. +""" # ╔═╡ 5522a5c2-ff7e-4093-86d0-ffb879fcc6e8 md" -#### Study case 3 +### Case study 3 " # ╔═╡ 5e3da6ca-656e-41ad-ba04-bf3d6b51939d -md" -We initialize a vector `params_uncert_log` with the parameter values and uncertainty (standard deviation) **in all parameters**: -" +md""" +We initialize a vector `parms_uncert_log` with the parameter values and uncertainty (standard deviation) **in all parameters**: +""" # ╔═╡ 7be29c69-7fa0-4f41-80c9-f69aa01e8ea4 -params_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0±1.5] +parms_uncert_log = [:μ => 0.07±0.02, :Wf => 10.0±1.5] # ╔═╡ d48c378e-df1c-45f7-9c3a-c5e9640e25a5 -md" +md""" We create the corresponding ODE problem and store it in `oprob_uncert_log`: -" +""" # ╔═╡ f209f4fb-985c-4a9b-b4dd-d50cc3fae1cf -oprob_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, params_uncert_log) +oprob_uncert_log = ODEProblem(growth_mod_log, u0_log, tspan, parms_uncert_log) # ╔═╡ de96646b-b752-4f72-ac78-2a503478bf50 -md" +md""" We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol_uncert_log`: -" +""" # ╔═╡ 6b1a3c2a-0a88-41ad-8866-41a33da5cd90 osol_uncert_log = solve(oprob_uncert_log, Tsit5(), saveat=2.0) # ╔═╡ 127ead82-2e2d-4c2e-9417-0490b6839514 -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ a5c254ad-b629-4a9e-89e9-90460e84620c plot(osol_uncert_log) # ╔═╡ c181d973-4fcc-4fbe-b37c-5f71dd33b0d2 -md" +md""" Now you have the combined effect of the uncertainty in the parameter $\mu$ as well as in the parameter $W_f$. -" +""" # ╔═╡ f4406d4e-d9f9-45fe-865c-405dbbb3127f -md" +md""" ## Exercises -" +""" # ╔═╡ e9f49fbc-d877-4db7-b320-a2db964ecaaa -md" +md""" ### Exercise 1 - Uncertainty analysis of the exponential growth model Perform an uncertainty analysis of the exponential growth model. Use the parameter uncertainties mentioned in the *Table* in the *Grass growth models* sections. -" +""" # ╔═╡ 5201995b-3298-4984-8839-dd25dc73a20f md" @@ -291,34 +393,36 @@ A possible *reaction network object* for the exponential growth model can be imp # ╔═╡ 97bdfce1-a670-469a-b8a9-34e1d1964409 growth_exp = @reaction_network begin + @species W(t)=2.0 + @parameters μ Wf μ*Wf, 0 --> W μ, W --> 0 end # ╔═╡ 52991cca-0365-48fb-b157-62a71c5573ee md" -The vector `u₀_exp` with the initial condition is: +The vector `u0_exp` with the initial condition is: " # ╔═╡ f7241453-755e-466c-a963-a504bc9ea446 u0_exp = [:W => 2.0] # ╔═╡ 11b81610-90b1-4d71-aeba-0bf33709dbda -md" -Initialize a vector `params_uncert_exp` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ -**Remark**: you can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. -" +md""" +Initialize a vector `parms_uncert_exp` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ +**Remark**: You can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. +""" # ╔═╡ 1af01af7-f5ee-4029-94ff-418c2f9748f2 -# params_uncert_exp = missing # Uncomment and complete the instruction +# parms_uncert_exp = missing # ╔═╡ c25160d9-471a-414c-be99-10aa2efbb6d4 -md" +md""" Create the corresponding ODE problem and store it in `oprob_uncert_exp`: -" +""" # ╔═╡ 7fcec514-192b-4694-b149-14f5f87bae17 -# oprob_uncert_exp = missing # Uncomment and complete the instruction +# oprob_uncert_exp = missing # ╔═╡ bd120dbc-827a-499a-a11c-25423ec5f397 md" @@ -326,29 +430,29 @@ Solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `os " # ╔═╡ 2dd58666-80ae-4346-ae15-34792eea2d25 -# osol_uncert_exp = missing # Uncomment and complete the instruction +# osol_uncert_exp = missing # ╔═╡ 8059e00b-8e9c-4c11-a484-bea2a0bf155e -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ a89677f3-602b-4ef4-8f1a-ab203584f8df -# missing # Uncomment and complete the instruction +# missing # ╔═╡ 2a470290-bca7-459f-b4f8-a5ac9fc031ed -md" +md""" Draw your conclusions: - missing -" +""" # ╔═╡ ded3aa76-ba09-499e-8c2d-52ff81ccfe02 -md" +md""" ### Exercise 2 - Uncertainty analysis of the exponential Gompertz model Perform an uncertainty analysis of the Gompertz growth model. Use the parameter uncertainties mentioned in the *Table* in the *Grass growth models* sections. -" +""" # ╔═╡ e245e5dd-bb11-455b-8de3-ae0ec1763399 md" @@ -357,35 +461,35 @@ A possible *reaction network object* for the Gompertz growth model can be implem # ╔═╡ 3e16bd3a-898a-49f9-9ed1-8998be7dc045 growth_gom = @reaction_network begin - # -μ, W --> ∅ - # D*log(W), W --> ∅ - μ-D*log(W), W --> 2W + @species W(t)=2.0 + @parameters μ d + μ-d*log(W), W --> 2W end # ╔═╡ db2f002f-20b7-4b5c-8d59-85b089fb3e59 md" -The vector `u₀_gom` with the initial condition is: +The vector `u0_gom` with the initial condition is: " # ╔═╡ 22a0013b-04f9-4148-9344-10bee1def0cf u0_gom = [:W => 2.0] # ╔═╡ f4fff573-6533-4ac3-897d-840a9aaf55f1 -md" -Initialize a vector `params_uncert_gom` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ -**Remark**: you can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. -" +md""" +Initialize a vector `parms_uncert_gom` with the parameter values and their uncertainties (standard deviation) **in all parameters**.\ +**Remark**: You can use the same variable and leave a single uncertainty if you want to see the effect of the uncertainty in only one parameter later on. +""" # ╔═╡ 5ab3aca2-e570-4802-815b-5614680b427d -# params_uncert_gom = missing # Uncomment and complete the instruction +# parms_uncert_gom = missing # ╔═╡ 890fe099-28ae-4e4f-aa1e-8d9ab6240130 -md" +md""" Create the corresponding ODE problem and store it in `oprob_uncert_log`: -" +""" # ╔═╡ 1aed2273-7b3e-4a49-a8d0-91251c58b700 -# oprob_uncert_gom = missing # Uncomment and complete the instruction +# oprob_uncert_gom = missing # ╔═╡ 76da6edd-9e9f-4f5b-8673-b737c8f9c5f8 md" @@ -393,31 +497,32 @@ Solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `os " # ╔═╡ 8b10af49-e3d0-426e-9147-7666262d1c54 -# osol_uncert_gom = missing # Uncomment and complete the instruction +# osol_uncert_gom = missing # ╔═╡ 6e8318bb-c884-4726-bb3d-9e0cfb7f40c1 -md" +md""" Plot the results (simulation of the output variable $W$ and uncertainty band): -" +""" # ╔═╡ 92785455-14d2-4adb-9731-8575527c5217 -# missing # Uncomment and complete the instruction +# missing # ╔═╡ 9d11df31-57b2-4825-a9db-34b422d5f084 -md" +md""" Draw your conclusions: - missing -" +""" # ╔═╡ Cell order: +# ╟─241a8a65-c59f-44f1-be39-d5edd1321b49 # ╠═d2c4d230-0943-11ef-3aad-5719e74bb20e -# ╠═47790920-74d2-4f21-b0f2-45dc19d25156 # ╠═8349306c-e98c-4221-9a1b-322fde3e18cb # ╠═83f3d978-aefa-40c1-a647-b26e837aeed6 -# ╠═73d5f921-82e4-4160-a3c4-fc1935f4c58d +# ╠═5185d0eb-7392-4775-9336-3e0f9e1449ce +# ╠═1c8ea6fe-650b-4c52-ba14-47efe3bd3e39 +# ╠═3eb3e651-afcf-41f8-a652-154a4f7ae07f # ╟─1f5e389e-e003-4e73-8f18-b5ad6340a912 -# ╟─241a8a65-c59f-44f1-be39-d5edd1321b49 # ╟─0f518e36-bc96-4599-807e-728504e5ca7b # ╟─da3cdc89-b911-4af9-9d4a-526301cba581 # ╟─78afaded-5a19-4386-aa76-7974977ea354 @@ -425,13 +530,13 @@ Draw your conclusions: # ╟─e5ab4490-6fd4-4b51-bfa7-1366438efefc # ╟─3e2750e7-220a-47b2-b445-eb3315734dec # ╟─c874a08c-7d82-4631-b611-c598dcaded09 +# ╟─87272d30-d95a-4b36-84ef-84bb9363c19e # ╟─9d3192bc-ee69-44a1-9273-06b8a55c60ef # ╟─8acb7ea2-54ef-428b-bb4f-364377d2c38d -# ╟─460d98ef-2177-49bf-87a4-35412f4183ed +# ╠═460d98ef-2177-49bf-87a4-35412f4183ed # ╟─38e0f746-1d04-46ba-81cc-21522b9ce6a4 # ╠═2538bd6e-027f-47fe-813c-db0a02586d2c # ╟─72b53ba6-dd05-4f2d-a65e-52783e76a1e9 -# ╠═5185d0eb-7392-4775-9336-3e0f9e1449ce # ╟─1a67db96-6973-446d-aa5b-253dd0008a73 # ╠═048112b0-10f1-47c5-834f-c3851d3078e8 # ╟─037cbf2b-dc0c-4b62-879b-aa702c964877 @@ -444,9 +549,31 @@ Draw your conclusions: # ╠═2fc91add-8468-4402-bcd8-42da5544f612 # ╟─5fa951d3-23ee-4510-873c-8158df4f2faf # ╠═d54635e8-1693-42bc-91f4-7e504a8661c7 -# ╠═707c6023-49ae-4a13-a220-d827a76a352c +# ╟─707c6023-49ae-4a13-a220-d827a76a352c # ╠═48e92497-f6b3-484c-bc5f-3ce69f4f1a91 # ╟─9b30b23b-6e92-49cf-bf93-e12a187fa2ba +# ╟─5dcb1ef0-d2a8-41fa-9592-2a86cf8364e6 +# ╟─0d883d09-1c92-4bb5-9bbf-3c87c856a5fa +# ╠═0e0dc268-faaf-462e-a980-9398c6015e02 +# ╟─e34f4b23-ac1b-4b88-b959-65085cef4b4f +# ╠═738448dd-d1c9-4fdc-af5d-fad4d91adf38 +# ╠═51579212-dd33-42b9-a19e-221ab7630c1a +# ╠═22bf3f50-1909-4a2c-820f-45b87f6bce14 +# ╟─616ee1d4-1f56-4c4a-a658-9f9351949f59 +# ╟─6b4a862e-6fb6-40ba-9e3f-43d57f169994 +# ╠═62afc495-29c3-4baa-b45a-94b84149d15d +# ╠═691a8b68-935e-4012-aa25-9cfdc323b1de +# ╟─f49a0085-396f-469e-afb8-2b7884f9bb1b +# ╠═ff5fcd0a-a37a-46d5-8b4b-f3f7a9a457a3 +# ╟─c0e72c96-23b6-4903-bd31-2978f6aebb47 +# ╟─988e9cd3-29a8-4fb1-a706-82d8356b9f20 +# ╠═980596b3-52e3-4f15-9678-f15f13f22688 +# ╠═851a7da1-ab7f-4a7a-ae1e-07a05c1e8052 +# ╠═8a9ec64f-2492-4f86-98d5-0e7c9512f5bc +# ╟─67f1348d-979c-42c4-a9d4-3031dbb5179a +# ╠═d4b354f7-369a-497b-bc9d-68ab6cb5f2b4 +# ╠═6859b3f8-bd1c-47b7-a4e8-0c38501f5342 +# ╟─63975327-5391-4de0-a538-4dcb503211b0 # ╟─86f8d937-c42a-47ff-a0f6-6cae9d637fcf # ╟─0755ea1c-50a0-4c8a-9adc-be2d63948a10 # ╠═7befd387-f3a6-4b34-90b0-5f04b5429252 @@ -474,14 +601,14 @@ Draw your conclusions: # ╟─52991cca-0365-48fb-b157-62a71c5573ee # ╠═f7241453-755e-466c-a963-a504bc9ea446 # ╟─11b81610-90b1-4d71-aeba-0bf33709dbda -# ╟─1af01af7-f5ee-4029-94ff-418c2f9748f2 +# ╠═1af01af7-f5ee-4029-94ff-418c2f9748f2 # ╟─c25160d9-471a-414c-be99-10aa2efbb6d4 # ╠═7fcec514-192b-4694-b149-14f5f87bae17 # ╟─bd120dbc-827a-499a-a11c-25423ec5f397 # ╠═2dd58666-80ae-4346-ae15-34792eea2d25 -# ╠═8059e00b-8e9c-4c11-a484-bea2a0bf155e +# ╟─8059e00b-8e9c-4c11-a484-bea2a0bf155e # ╠═a89677f3-602b-4ef4-8f1a-ab203584f8df -# ╟─2a470290-bca7-459f-b4f8-a5ac9fc031ed +# ╠═2a470290-bca7-459f-b4f8-a5ac9fc031ed # ╟─ded3aa76-ba09-499e-8c2d-52ff81ccfe02 # ╟─e245e5dd-bb11-455b-8de3-ae0ec1763399 # ╠═3e16bd3a-898a-49f9-9ed1-8998be7dc045 @@ -489,10 +616,10 @@ Draw your conclusions: # ╠═22a0013b-04f9-4148-9344-10bee1def0cf # ╟─f4fff573-6533-4ac3-897d-840a9aaf55f1 # ╠═5ab3aca2-e570-4802-815b-5614680b427d -# ╠═890fe099-28ae-4e4f-aa1e-8d9ab6240130 +# ╟─890fe099-28ae-4e4f-aa1e-8d9ab6240130 # ╠═1aed2273-7b3e-4a49-a8d0-91251c58b700 # ╟─76da6edd-9e9f-4f5b-8673-b737c8f9c5f8 # ╠═8b10af49-e3d0-426e-9147-7666262d1c54 # ╟─6e8318bb-c884-4726-bb3d-9e0cfb7f40c1 # ╠═92785455-14d2-4adb-9731-8575527c5217 -# ╟─9d11df31-57b2-4825-a9db-34b422d5f084 +# ╠═9d11df31-57b2-4825-a9db-34b422d5f084