From e1c435885c578ee44e3c0631afb077970f5b8a6b Mon Sep 17 00:00:00 2001 From: wclaerhout Date: Mon, 19 Oct 2026 16:48:24 +0200 Subject: [PATCH 1/2] merged ferementer sensitivity and uncertainty into one notebook and changed order of intro a bit to be more logical --- .../P7_sens_uncert/fermenter_mutant.jl | 711 ++++++++++++++++ .../sens_fermenter_monod_sol.jl | 781 ------------------ .../uncert_bitrophic_model_sol.jl | 2 +- .../uncert_fermenter_monod_sol.jl | 182 ---- .../P7_sens_uncert/uncert_intro_sol.jl | 34 +- 5 files changed, 729 insertions(+), 981 deletions(-) create mode 100644 exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant.jl delete mode 100644 exercises/solved_notebooks/P7_sens_uncert/sens_fermenter_monod_sol.jl delete mode 100644 exercises/solved_notebooks/P7_sens_uncert/uncert_fermenter_monod_sol.jl diff --git a/exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant.jl b/exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant.jl new file mode 100644 index 00000000..6903dc70 --- /dev/null +++ b/exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant.jl @@ -0,0 +1,711 @@ +### A Pluto.jl notebook ### +# v0.20.28 + +using Markdown +using InteractiveUtils + +# ╔═╡ 9d8acf40-e635-4dc4-9938-ec63ab68e3bd +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 + +# ╔═╡ 1ec7def2-0c43-11ef-0850-ddefa680d2a4 +using Markdown + +# ╔═╡ 2861bc17-7810-461e-8857-a55b6392977a +using InteractiveUtils + +# ╔═╡ 093679cc-122b-4b0a-84c3-15cf3a6f5c1b +using Catalyst + +# ╔═╡ 0720732e-105a-40c2-bf45-f8fbbc10e5bc +using OrdinaryDiffEq + +# ╔═╡ 28bb60c9-ab56-497f-b039-0d6577501be5 +using StatsPlots, PlutoUI; TableOfContents() + +# ╔═╡ 0a65afbc-643a-4fe4-b22d-c5fe54e3be1f +using Turing + +# ╔═╡ d1b6739f-1e0c-452d-8763-717dd650b2ec +using ForwardDiff + +# ╔═╡ 162b330b-02f7-4179-b950-d960eca68206 +using Random + +# ╔═╡ 437364b6-a953-428b-a891-e17ffd531656 +md""" +# Exercise: Fermenter - Switching to a mutant strain? +""" + +# ╔═╡ c44e10f3-8452-4a60-b840-f4a7e6ec3238 +md""" +As a process engineer, you are responsible for a continuously stirred 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$ [$L/h$] containing substrate at a concentration $S_{in}$ [$g/L$]. The yield $Y$ [$g$ biomass$/g$ substrate] tells how much biomass is formed per gram of substrate consumed. As in the previous practica, the growth is modelled with Monod kinetics: + +$$\begin{eqnarray*} +S \xrightarrow[\quad\quad]{r} Y \, X \quad\quad\quad\quad r = \mu \, X \quad \textrm{with} \quad \mu = \mu_{max} \, \cfrac{S}{S + K_s} +\end{eqnarray*}$$ +""" + +# ╔═╡ f3065171-5526-476a-bfd7-a2b95bb24d61 +md""" +## The question + +A colleague has found a mutant strain with a **higher yield** $Y$, and management asks you a simple question: *should we switch to the mutant?* + +The catch: the maximal growth rate $\mu_{max}$ of the mutant is not known. A higher yield often comes at the cost of slower growth, and a slower grower leaves more substrate unused in the outflow. So the mutant could be better or worse than the current strain. + +You will answer this question with the tools of this practical: +1. **Uncertainty propagation**: what does our uncertainty about $\mu_{max}$ mean for the biomass? How sure are we that the mutant is better? +2. **Sensitivity analysis**: when does the biomass depend most strongly on $\mu_{max}$? That is the moment where a measurement tells us the most. +3. **Bayesian inference**: update our belief about $\mu_{max}$ with a few measurements, and propagate the reduced uncertainty again. +""" + +# ╔═╡ 356822df-bad3-41f9-8556-de04234b0a09 +md""" +The *reaction network object* for this model is: +""" + +# ╔═╡ c551c33e-95b0-4ea8-85d4-6c68db72b8bc +fermenter_monod = @reaction_network begin + X * mm(S, μmax, Ks), S => Y*X + Q/V, (S, X) --> 0 + Q/V*Sin, ∅ --> S +end + +# ╔═╡ 69f8c265-fe91-4e4b-971b-395325d4474e +md""" +which results in the following differential equations: +""" + +# ╔═╡ 2806654f-7a24-4f8e-8b5a-3129b2431c95 +md""" +$$\begin{eqnarray*} +\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ +\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X +\end{eqnarray*}$$ +""" + +# ╔═╡ 5d6a881d-421e-4dea-9668-38b3f3378bf1 +# osys = missing +osys = convert(ODESystem, fermenter_monod) + +# ╔═╡ 5802078d-12bc-489e-8566-15f77a11e8ba +md""" +## The current strain + +For the strain currently in use, the parameters are well known: +- $\mu_{max} = 0.40\;h^{-1}$ +- $K_s = 0.015 \;g/L$ +- $Y = 0.67$ +- $S_{in} = 0.022 \;g/L$ +- $Q = 2.0\;L/h$ +- $V = 40.0\;L$ + +At $t=0$ no substrate is present in the reactor and the biomass concentration is $0.0005\;g/L$. We judge a strain by the biomass concentration in the outflow after $100\,h$. + +Initialize a vector `u0` with the initial conditions, set the timespan `tspan` to $[0, 100]\,h$, the time step `dt` to `0.5`, and initialize a vector `parms` with the parameter values. +""" + +# ╔═╡ 6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde +# u0 = missing +u0 = [:S => 0.0, :X => 0.0005] + +# ╔═╡ 1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b +# tspan = missing +tspan = (0.0, 100.0) + +# ╔═╡ 9e3c427b-8a56-4787-9b28-43ec3932ddbe +# dt = missing +dt = 0.5 + +# ╔═╡ 62f4049b-617f-4268-beba-b2c634b46642 +# parms = missing +parms = [:μmax => 0.40, :Ks => 0.015, :Y => 0.67, :Q => 2.0, :V => 40.0, :Sin => 0.022] + +# ╔═╡ 0e4c73d3-1630-463b-994d-073170c1a6c3 +md""" +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. +""" + +# ╔═╡ 3cad11d6-c2eb-4865-858a-d407aa25d780 +# oprob = missing +oprob = ODEProblem(fermenter_monod, u0, tspan, parms, combinatoric_ratelaws=false) + +# ╔═╡ 40168902-0780-4f92-8eca-23982c80511e +# osol = missing +osol = solve(oprob, Tsit5(), saveat=dt) + +# ╔═╡ 4d753fdb-089d-41dc-89ef-45a1733af786 +# missing +plot(osol, title="Current strain", xlabel="Time (hours)", ylabel="concentration (g/l)", linewidth=2) + +# ╔═╡ ca660e2c-9348-40b1-b92d-415068bea514 +md""" +Store the biomass concentration at the end of the simulation in `X_current`. This is the value the mutant has to beat. +""" + +# ╔═╡ ec042e61-04b8-4d51-8e7d-35e1646cc5ba +# X_current = missing +X_current = osol[:X][end] + +# ╔═╡ 9ee22028-7a37-4510-9b61-198d3bd99a73 +md""" +## The mutant strain +""" + +# ╔═╡ cc9ac43b-448a-4a67-ba0b-8fdf21a4236e +md""" +Some batch tests show that the mutant has a yield of $Y = 0.70$ instead of $0.67$. All other parameters are the same as for the current strain, except $\mu_{max}$, which is unknown. + +From the literature on mutants with an improved yield, you expect $\mu_{max}$ to be around $0.35\;h^{-1}$, give or take $15\%$. You express this **prior belief** with a log-normal distribution `LogNormal(log(0.35), 0.15)`. +""" + +# ╔═╡ f5482531-b23f-472b-968f-20f7a6471778 +begin + plot(LogNormal(log(0.35), 0.15), xlabel="μmax (1/h)", ylabel="density", label="Prior μmax mutant", linewidth=2) + vline!([0.40], label="μmax current strain", linestyle=:dash) +end + +# ╔═╡ ff21b3cd-64b1-4d09-8eb1-afc617585fc8 +md""" +Create an ODE problem `oprob_mut` for the mutant by remaking `oprob` with $Y = 0.70$ and the expected value $\mu_{max} = 0.35$. + +!!! note + Remember that you can remake an ODEProblem with new parameter values by using `remake(oprob, p = [...])`. +""" + +# ╔═╡ 225eaaff-42a2-4aca-94e1-afcd988de4eb +# oprob_mut = missing +oprob_mut = remake(oprob, p = [:μmax => 0.35, :Y => 0.70]) + +# ╔═╡ 18ffe877-c579-4b4c-80bd-b5d50baf36bb +md""" +## Monte Carlo error propagation + +We don't know $\mu_{max}$, but we do know which values are plausible. Using Monte Carlo simulation, we draw many values from the prior, simulate the mutant for each value, and look at the spread of the outputs. This shows how the uncertainty on the parameter **propagates** to the uncertainty on the output. + +Make a Turing model `mutant_deviation` in which you implement the prior for `μmax_dev` and return the solution of the remade ODE problem. +""" + +# ╔═╡ 7821dc21-3e6b-4ad1-829e-f29dcd020144 +# @model function mutant_deviation() +# μmax_dev ~ missing +# new_prob = missing +# return missing +# end +@model function mutant_deviation() + μmax_dev ~ LogNormal(log(0.35), 0.15) + new_prob = remake(oprob_mut, p = [:μmax => μmax_dev]) + return solve(new_prob, Tsit5(), saveat=dt) +end + +# ╔═╡ 3051d39a-fc2f-40e4-a688-4d6e75cfc2d1 +md""" +Now sample 500 solutions from your `mutant_deviation` model. +""" + +# ╔═╡ 33df05c9-08d3-4107-9240-931658cda2a4 +# missing +begin + dev_model = mutant_deviation() + chain = sample(dev_model, Prior(), 500) + solutions = generated_quantities(dev_model, chain) +end; + +# ╔═╡ 02c5cf65-ebd3-477b-9a3e-8855bd46f55c +md""" +Plot the biomass of all Monte Carlo simulations together with the biomass of the current strain. Use `idxs=:X`, `label=false`, `color=:gray`, `alpha=0.4` and `linestyle=:dot`. +""" + +# ╔═╡ aa833d92-567f-4971-a360-e60bce424c8e +# begin +# p1 = plot(title="Monte Carlo of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") +# for missing in missing +# missing +# end +# plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) +# p1 +# end +begin + p1 = plot(title="Monte Carlo of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") + for solution in solutions + plot!(solution, idxs=:X, label=false, color=:gray, alpha=0.4, linestyle=:dot) + end + plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) + p1 +end + +# ╔═╡ ffc66a94-5435-4ca6-af87-76c85a599906 +md""" +Create a vector `x1` that stores the biomass at the end of the simulations, make a histogram of it, and compute the probability `P_better_prior` that the mutant produces more biomass than the current strain. + +!!! note + The fraction of elements of a vector `x` that are larger than `a` can be computed with `mean(x .> a)`. +""" + +# ╔═╡ d3afd5b4-034f-415c-a64a-776deed74d5c +# x1 = missing +x1 = [solution[:X][end] for solution in solutions]; + +# ╔═╡ 1d95728e-079e-4e9b-952a-c4f692661813 +# missing +begin + histogram(x1, title="Prior histogram of X(100 h)", xlabel="X (g/l)", ylabel="counts", label="Mutant", legend=:topleft) + vline!([X_current], linewidth=3, color=:red, label="Current strain") +end + +# ╔═╡ 02213469-a233-4759-a18d-29a76a4c13c6 +# P_better_prior = missing +P_better_prior = mean(x1 .> X_current) + +# ╔═╡ 97d2065c-6234-4b13-8e4c-6782b8301416 +md""" +!!! questions + - When is the spread on $X$ largest, and when is it smallest? Can you explain this? + - How likely is it that the mutant beats the current strain? Would you advise management to switch based on this? +""" + +# ╔═╡ d2acbfb9-a7ff-48b0-9d9a-9492ecb9c1bf +md""" +Answers: +- missing +- missing +""" +#= +- The spread is largest during the start-up phase (around 40-60 h). There, the growth rate decides how far the population has grown. Towards the end the spread shrinks again, because almost all substrate is consumed and X tends to Y*Sin, regardless of μmax. +- About 65-70%. That is better than a coin flip, but switching a production process is costly, and a 1 in 3 chance of ending up worse is too risky. We need more information about μmax. The question is: which information? +=# + +# ╔═╡ eb613ead-6aea-41fe-acb7-f679b5a97c35 +md""" +## Local Sensitivity Analysis (LSA) + +The lab can run the mutant in a pilot reactor under the same conditions, but you can only afford **three measurements** of $X$, all at the same time point. Measurements are expensive, so we want them at the moment where $X$ depends most strongly on $\mu_{max}$. There, a small change in $\mu_{max}$ causes a large change in $X$. The other way around, a measurement of $X$ at that moment pins down $\mu_{max}$ best. + +So, we will first perform a normalized local sensitivity analysis to pin down the moment, $X$ is most sensitive to a change in $μmax$. +""" + +# ╔═╡ 67f572b0-3485-40bc-8f64-cc58471b38a9 +md""" +Write a solution function `mutant_sim` with as argument a vector with the parameter `μmax` of the mutant, and that returns the solution. Based on it, write a function `mutant_sim_X` that only returns the output $X$. +""" + +# ╔═╡ 565b6e18-012f-4b38-bf21-bdc683688107 +# function mutant_sim(parms) +# missing +# ... +# end +function mutant_sim(parms) + μmax_val = parms[1] + oprob = remake(oprob_mut, p = [:μmax => μmax_val]) + osol = solve(oprob, Tsit5(), saveat=dt) + return osol +end + +# ╔═╡ b7f14c48-ba19-4fbc-bded-829a34bd7ac8 +# mutant_sim_X(parms) = missing +mutant_sim_X(parms) = mutant_sim(parms)[:X] + +# ╔═╡ 73345fe4-6877-476d-8ce1-5a0056c73cc6 +md""" +Make the time vector `t_vals` with timestep `dt` equal to $0.5$ and compute the output `X_sim` for the expected value $\mu_{max} = 0.35$, stored in `μmax_val`. +""" + +# ╔═╡ 75ab68b2-1e88-421f-a089-4c0b3b72788d +# t_vals = missing +t_vals = tspan[1]:dt:tspan[2] + +# ╔═╡ 82ae802b-2a7a-4525-9cc6-dfaa552ba990 +# μmax_val = missing +μmax_val = 0.35 + +# ╔═╡ 8a655d11-0d45-4e7c-bcbf-8ea9659516aa +# X_sim = missing +X_sim = mutant_sim_X([μmax_val]) + +# ╔═╡ 9a97a7c7-bc3e-434e-9dce-7c859417aac2 +md""" +Use `ForwardDiff.jacobian` to compute the sensitivity `sens_X_on_μmax`, and compute the normalized sensitivity `sens_X_on_μmax_rel`. +""" + +# ╔═╡ d7d3028a-51c1-41e0-a865-b521ed235548 +# sens_X_on_μmax = missing +sens_X_on_μmax = ForwardDiff.jacobian(mutant_sim_X, [μmax_val])[:, 1] + +# ╔═╡ cb6047f2-3884-4fbd-96f2-f15cf684bea8 +# sens_X_on_μmax_rel = missing +sens_X_on_μmax_rel = sens_X_on_μmax .* μmax_val ./ X_sim + +# ╔═╡ f623b4b4-6010-45fa-b9a9-bc796f0fb106 +md""" +Identify the time `t_star` at which the biomass is most sensitive to variations in $\mu_{max}$, and plot the normalized sensitivity function together with `t_star`. + +!!! 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. +""" + +# ╔═╡ 9c3fc6d8-8e5e-425d-9caa-56a9e0b03b99 +# begin +# t_star_idx = missing +# t_star = t_vals[t_star_idx] +# +# plot(t_vals, abs.(sens_X_on_μmax_rel), +# title="Sensitivity of X w.r.t. μmax", xlabel="Time (hours)", +# ylabel="normalized sensitivity", label=false, linewidth=2) +# vline!([t_star], label="t* = $(t_star) h", linestyle=:dash, color=:red) +# end +begin + t_star_idx = argmax(abs.(sens_X_on_μmax_rel)) + t_star = t_vals[t_star_idx] + + plot(t_vals, abs.(sens_X_on_μmax_rel), + title="Sensitivity of X w.r.t. μmax", xlabel="Time (hours)", + ylabel="normalized sensitivity", label=false, linewidth=2) + vline!([t_star], label="t* = $(t_star) h", linestyle=:dash, color=:red) +end + +# ╔═╡ b1825a21-e711-4f99-9a4b-51f3fb73ff82 +md""" +!!! questions + - Compare the sensitivity function with the Monte Carlo simulations above. What do you notice? + - The sensitivity at $t = 100\,h$ is small, while that is exactly the time we care about. Why is it small? Is it a good idea to measure at $t = 100\,h$? +""" + +# ╔═╡ 44234af3-6015-49b6-8709-1f3a3cc0f69d +md""" +Answers: +- missing +- missing +""" +#= +- The spread of the Monte Carlo simulations is largest where the sensitivity is largest. To first order, the uncertainty on the output is σ_X(t) ≈ |∂X/∂μmax(t)| σ_μmax: the sensitivity function tells us *where* parameter uncertainty propagates into output uncertainty. +- At the end, the substrate is almost completely consumed, so X ≈ Y(Sin - S) with S small, and μmax hardly matters. A measurement at 100 h therefore tells us very little about μmax. Yet this small dependency still decides whether the mutant beats the current strain, so we need to know μmax precisely, and we get that by measuring at t*. +=# + +# ╔═╡ a00afabe-7bc1-4d8b-9b55-2953c7add9f5 +md""" +## Measuring campaign + +The lab ran the mutant in the pilot reactor and sampled it at your chosen time. We have defined a function `measure_X(t)` in this notebook that returns three hypothetical measurements of $X$ at the given timepoint `t`, with a relative measurement error of $5\%$. Use this function to get the measurements `X_meas` at the optimal time `t_star`. +""" + +# ╔═╡ 635ff400-639c-4c80-93e6-e8721a836c16 +begin + μmax_real = 0.37 + real_solution = solve(remake(oprob_mut, p = [:μmax => μmax_real]), Tsit5(), saveat=dt) + measure_X(t) = real_solution(t, idxs=:X) .* (1 .+ 0.05 .* randn(Xoshiro(round(Int, 10t)), 3)) +end; + +# ╔═╡ 93cc5a4f-f40a-4d28-abc3-8d6790fd60dc +# X_meas = missing +X_meas = measure_X(t_star) + +# ╔═╡ 3e93a001-9401-4a1f-b283-c79fd0c43434 +md""" +### Conditioning on the measurements +""" + +# ╔═╡ 04ae6cd7-2d59-4bc0-a894-d9d5129e243d +md""" +Based on these measurements we can update what we believed before (the prior) with what the data tell us (the likelihood) into an updated belief (the posterior): + +Create a Turing model `mutant_meas(t_meas)` with: +- the same prior for `μmax_dev` as before; +- a prior `Exponential()` for the measurement error `σ_X`; +- the simulated biomass `X_pred` at time `t_meas`, which you can get with `sol(t_meas, idxs=:X)`; +- the three measurements `X_meas`, each normally distributed around `X_pred` with standard deviation `σ_X`; +- the solution as return value. +""" + +# ╔═╡ 02096fef-224e-4896-9f8d-3f7aba7334e1 +# @model function mutant_meas(t_meas) +# μmax_dev ~ missing +# σ_X ~ missing +# +# sol = missing +# X_pred = missing +# +# X_meas = zeros(3) +# for i in 1:length(X_meas) +# X_meas[i] ~ missing +# end +# return sol +# end +@model function mutant_meas(t_meas) + μmax_dev ~ LogNormal(log(0.35), 0.15) + σ_X ~ Exponential() + + sol = solve(remake(oprob_mut, p = [:μmax => μmax_dev]), Tsit5(), saveat=dt) + X_pred = sol(t_meas, idxs=:X) + + X_meas = zeros(3) + for i in 1:length(X_meas) + X_meas[i] ~ Normal(X_pred, σ_X) + end + return sol +end + +# ╔═╡ 354f0768-efc0-4304-aa04-c6a68f961aa6 +md""" +Condition the model on the measurements and store it in `mutant_cond_model`. +""" + +# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 +# mutant_cond_model = missing +mutant_cond_model = mutant_meas(t_star) | (X_meas = X_meas,) + +# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 +md""" +### MCMC +""" + +# ╔═╡ f99ca771-8c14-4e19-a1ae-478a2103716a +# mutant_chain = missing +mutant_chain = sample(mutant_cond_model, NUTS(), 500) + +# ╔═╡ e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +md""" +Plot the chain to check if it properly converged to the posterior distribution. +""" + +# ╔═╡ 50b30e65-72d3-4d8d-a739-ae969e0f50a8 +# missing +plot(mutant_chain) + +# ╔═╡ a65710f6-70ce-4348-923d-ed2ed2b4a469 +md""" +Now, get the sampled posterior solutions `solutions_mutant` and plot them on top of the prior Monte Carlo simulations, together with the current strain and the measurements. +""" + +# ╔═╡ 12c9f810-ef5b-49a4-8fdd-6efb1a040440 +# solutions_mutant = missing +solutions_mutant = generated_quantities(mutant_cond_model, mutant_chain); + +# ╔═╡ e2e11092-4498-437b-8c67-8d77302893fe +md""" +The posterior chance that the mutant is better is: +""" + +# ╔═╡ cae4d9c4-3790-423e-8321-d44c0833aa40 +begin + sol_end = [solutions_mutant[i][:X][end] for i in eachindex(solutions_mutant)] + P_better_posterior = mean(sol_end .> X_current) + println("The chance the mutant is better is $(P_better_posterior*100)%") +end + +# ╔═╡ 49307c36-3518-4118-b453-91e634f92daf +# begin +# p2 = plot(title="Prior vs posterior of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") +# for missing in missing +# missing +# end +# for missing in missing +# missing +# end +# plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) +# scatter!(fill(t_star, 3), X_meas, label="Measurements", color=:black) +# p2 +# end +begin + p2 = plot(title="Prior vs posterior of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") + for solution in solutions + plot!(solution, idxs=:X, label=false, color=:gray, alpha=0.2, linestyle=:dot) + end + for solution_mutant in solutions_mutant + plot!(solution_mutant, idxs=:X, label=false, color=:orange, alpha=0.4) + end + plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) + scatter!(fill(t_star, 3), X_meas, label="Measurements", color=:black) + p2 +end + +# ╔═╡ aa479bcd-2314-47a5-83ba-b0d1a2654fa0 +md""" +Create a vector `x2` with the posterior biomass at the end of the simulations, compare its histogram with the prior one, and compute the probability `P_better_post` that the mutant beats the current strain. +""" + +# ╔═╡ 2f55a771-0d37-48e5-a9d9-24693bf1d177 +# x2 = missing +x2 = [solution_mutant[:X][end] for solution_mutant in solutions_mutant]; + +# ╔═╡ 18059d0c-cf37-44a1-86a8-859ac66a66f0 +# missing +begin + histogram(x1, normalize=:pdf, alpha=0.5, label="Prior", title="Histogram of X(100 h)", xlabel="X (g/l)", ylabel="density", legend=:topleft) + histogram!(x2, normalize=:pdf, alpha=0.5, label="Posterior") + vline!([X_current], linewidth=3, color=:red, label="Current strain") +end + +# ╔═╡ 40e73ddf-6c0d-4f83-860f-434d934193ed +# P_better_post = missing +P_better_post = mean(x2 .> X_current) + +# ╔═╡ e1c87260-25ee-4d71-8744-ed0d7cb75f6a +md""" +### What if we had measured at the end? + +Intuitively, you might have measured $X$ at $t = 100\,h$, since that is the quantity we care about. Repeat the inference with three measurements at $t = 100\,h$, and compute the probability `P_better_100` that the mutant is better. +""" + +# ╔═╡ caa5c04d-49af-4576-83a2-605bae1c2494 +# X_meas_100 = missing +X_meas_100 = measure_X(100.0) + +# ╔═╡ 6a53aff6-29d4-4d41-8236-f621ce519361 +# missing +begin + mutant_cond_model_100 = mutant_meas(100.0) | (X_meas = X_meas_100,) + mutant_chain_100 = sample(mutant_cond_model_100, NUTS(), 500) + solutions_mutant_100 = generated_quantities(mutant_cond_model_100, mutant_chain_100) + x3 = [solution_mutant[:X][end] for solution_mutant in solutions_mutant_100] +end; + +# ╔═╡ 9292c593-1f0f-47a6-a2ec-91a8aadbe6e6 +# P_better_100 = missing +P_better_100 = mean(x3 .> X_current) + +# ╔═╡ b1814ea6-0065-4ac2-a012-505d2bbd8e40 +md""" +Compare what we know about $\mu_{max}$ before and after both measuring campaigns: +""" + +# ╔═╡ 44cb7615-86b2-45f0-b87a-69fcbfe863bb +begin + density(vec(chain[:μmax_dev]), label="Prior", linewidth=2, title="What did we learn about μmax?", xlabel="μmax (1/h)", ylabel="density") + density!(vec(mutant_chain_100[:μmax_dev]), label="Posterior, measured at t = 100 h", linewidth=2) + density!(vec(mutant_chain[:μmax_dev]), label="Posterior, measured at t* = $(t_star) h", linewidth=2) +end + +# ╔═╡ f364fa1f-5c47-4ca6-99df-17e2b735d98b +md""" +!!! questions + - Compare `P_better_prior`, `P_better_100` and `P_better_post`. +""" + +# ╔═╡ 5f9873a5-af70-4ac2-be08-4f249b71eb76 +md""" +Answers: +- missing +- missing +""" +#= +- Measuring at t* shrinks the uncertainty on μmax from about 15% to about 3%, while measuring at 100 h only brings it down to about 10%. At t*, X is about ten times more sensitive to μmax than at 100 h, so the same 5% measurement error translates into a much smaller error on μmax. Sensitivity analysis tells you *where* to measure. +=# + +# ╔═╡ 3980da6b-7036-44d2-9d8f-3ddc1062d0ad +md""" +## Conclusions + +!!! questions + - What do you advise management? Support your advice with numbers. + - What did you learn/what did you use the information for from uncertainty propagation, sensitivity analysis and the Bayesian inference. +""" + +# ╔═╡ 8021123d-1b19-47bf-9db0-0d7c587bc00e +md""" +Answers: +- missing +- missing +- missing +""" +#= +- Switch to the mutant. With the measurements, μmax ≈ 0.37 1/h and the probability that the mutant produces more biomass after 100 h is practically 1, with about 3-4% more biomass than the current strain. +- Uncertainty propagation turned "we don't know μmax" into a probability for the decision (about 65-70%), and showed that this was not good enough. Sensitivity analysis showed when a measurement is informative (t* ≈ 45 h, not at the end). Bayesian inference combined the prior knowledge with only three measurements into a posterior that is narrow enough to decide, and that posterior was propagated again in the same way. +- E.g. the dose of a drug for a patient with an unknown clearance rate, deciding whether a new fertilizer increases crop yield, choosing when to sample a river to estimate a degradation rate, ... +=# + +# ╔═╡ Cell order: +# ╠═1ec7def2-0c43-11ef-0850-ddefa680d2a4 +# ╠═2861bc17-7810-461e-8857-a55b6392977a +# ╠═9d8acf40-e635-4dc4-9938-ec63ab68e3bd +# ╠═093679cc-122b-4b0a-84c3-15cf3a6f5c1b +# ╠═0720732e-105a-40c2-bf45-f8fbbc10e5bc +# ╠═28bb60c9-ab56-497f-b039-0d6577501be5 +# ╠═0a65afbc-643a-4fe4-b22d-c5fe54e3be1f +# ╠═d1b6739f-1e0c-452d-8763-717dd650b2ec +# ╠═162b330b-02f7-4179-b950-d960eca68206 +# ╟─437364b6-a953-428b-a891-e17ffd531656 +# ╟─c44e10f3-8452-4a60-b840-f4a7e6ec3238 +# ╟─f3065171-5526-476a-bfd7-a2b95bb24d61 +# ╟─356822df-bad3-41f9-8556-de04234b0a09 +# ╠═c551c33e-95b0-4ea8-85d4-6c68db72b8bc +# ╟─69f8c265-fe91-4e4b-971b-395325d4474e +# ╟─2806654f-7a24-4f8e-8b5a-3129b2431c95 +# ╠═5d6a881d-421e-4dea-9668-38b3f3378bf1 +# ╟─5802078d-12bc-489e-8566-15f77a11e8ba +# ╠═6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde +# ╠═1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b +# ╠═9e3c427b-8a56-4787-9b28-43ec3932ddbe +# ╠═62f4049b-617f-4268-beba-b2c634b46642 +# ╟─0e4c73d3-1630-463b-994d-073170c1a6c3 +# ╠═3cad11d6-c2eb-4865-858a-d407aa25d780 +# ╠═40168902-0780-4f92-8eca-23982c80511e +# ╠═4d753fdb-089d-41dc-89ef-45a1733af786 +# ╟─ca660e2c-9348-40b1-b92d-415068bea514 +# ╠═ec042e61-04b8-4d51-8e7d-35e1646cc5ba +# ╟─9ee22028-7a37-4510-9b61-198d3bd99a73 +# ╟─cc9ac43b-448a-4a67-ba0b-8fdf21a4236e +# ╟─f5482531-b23f-472b-968f-20f7a6471778 +# ╟─ff21b3cd-64b1-4d09-8eb1-afc617585fc8 +# ╠═225eaaff-42a2-4aca-94e1-afcd988de4eb +# ╟─18ffe877-c579-4b4c-80bd-b5d50baf36bb +# ╠═7821dc21-3e6b-4ad1-829e-f29dcd020144 +# ╟─3051d39a-fc2f-40e4-a688-4d6e75cfc2d1 +# ╠═33df05c9-08d3-4107-9240-931658cda2a4 +# ╟─02c5cf65-ebd3-477b-9a3e-8855bd46f55c +# ╠═aa833d92-567f-4971-a360-e60bce424c8e +# ╟─ffc66a94-5435-4ca6-af87-76c85a599906 +# ╠═d3afd5b4-034f-415c-a64a-776deed74d5c +# ╠═1d95728e-079e-4e9b-952a-c4f692661813 +# ╠═02213469-a233-4759-a18d-29a76a4c13c6 +# ╟─97d2065c-6234-4b13-8e4c-6782b8301416 +# ╟─d2acbfb9-a7ff-48b0-9d9a-9492ecb9c1bf +# ╟─eb613ead-6aea-41fe-acb7-f679b5a97c35 +# ╟─67f572b0-3485-40bc-8f64-cc58471b38a9 +# ╠═565b6e18-012f-4b38-bf21-bdc683688107 +# ╠═b7f14c48-ba19-4fbc-bded-829a34bd7ac8 +# ╟─73345fe4-6877-476d-8ce1-5a0056c73cc6 +# ╠═75ab68b2-1e88-421f-a089-4c0b3b72788d +# ╠═82ae802b-2a7a-4525-9cc6-dfaa552ba990 +# ╠═8a655d11-0d45-4e7c-bcbf-8ea9659516aa +# ╟─9a97a7c7-bc3e-434e-9dce-7c859417aac2 +# ╠═d7d3028a-51c1-41e0-a865-b521ed235548 +# ╠═cb6047f2-3884-4fbd-96f2-f15cf684bea8 +# ╟─f623b4b4-6010-45fa-b9a9-bc796f0fb106 +# ╠═9c3fc6d8-8e5e-425d-9caa-56a9e0b03b99 +# ╟─b1825a21-e711-4f99-9a4b-51f3fb73ff82 +# ╟─44234af3-6015-49b6-8709-1f3a3cc0f69d +# ╟─a00afabe-7bc1-4d8b-9b55-2953c7add9f5 +# ╟─635ff400-639c-4c80-93e6-e8721a836c16 +# ╠═93cc5a4f-f40a-4d28-abc3-8d6790fd60dc +# ╟─3e93a001-9401-4a1f-b283-c79fd0c43434 +# ╟─04ae6cd7-2d59-4bc0-a894-d9d5129e243d +# ╠═02096fef-224e-4896-9f8d-3f7aba7334e1 +# ╟─354f0768-efc0-4304-aa04-c6a68f961aa6 +# ╠═11e105de-270c-434a-bfb4-9118dfc9e461 +# ╟─e46afb98-16de-4d19-86b4-b920b6b929c5 +# ╠═f99ca771-8c14-4e19-a1ae-478a2103716a +# ╟─e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +# ╠═50b30e65-72d3-4d8d-a739-ae969e0f50a8 +# ╟─a65710f6-70ce-4348-923d-ed2ed2b4a469 +# ╠═12c9f810-ef5b-49a4-8fdd-6efb1a040440 +# ╟─e2e11092-4498-437b-8c67-8d77302893fe +# ╠═cae4d9c4-3790-423e-8321-d44c0833aa40 +# ╠═49307c36-3518-4118-b453-91e634f92daf +# ╟─aa479bcd-2314-47a5-83ba-b0d1a2654fa0 +# ╠═2f55a771-0d37-48e5-a9d9-24693bf1d177 +# ╠═18059d0c-cf37-44a1-86a8-859ac66a66f0 +# ╠═40e73ddf-6c0d-4f83-860f-434d934193ed +# ╟─e1c87260-25ee-4d71-8744-ed0d7cb75f6a +# ╠═caa5c04d-49af-4576-83a2-605bae1c2494 +# ╠═6a53aff6-29d4-4d41-8236-f621ce519361 +# ╠═9292c593-1f0f-47a6-a2ec-91a8aadbe6e6 +# ╟─b1814ea6-0065-4ac2-a012-505d2bbd8e40 +# ╟─44cb7615-86b2-45f0-b87a-69fcbfe863bb +# ╟─f364fa1f-5c47-4ca6-99df-17e2b735d98b +# ╟─5f9873a5-af70-4ac2-be08-4f249b71eb76 +# ╟─3980da6b-7036-44d2-9d8f-3ddc1062d0ad +# ╟─8021123d-1b19-47bf-9db0-0d7c587bc00e 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 deleted file mode 100644 index 0b56a62b..00000000 --- a/exercises/solved_notebooks/P7_sens_uncert/sens_fermenter_monod_sol.jl +++ /dev/null @@ -1,781 +0,0 @@ -### A Pluto.jl notebook ### -# v0.20.21 - -using Markdown -using InteractiveUtils - -# ╔═╡ 3ef93246-657d-4e77-9bf0-8380c64bfcfd -using Pkg; Pkg.activate("..") - -# ╔═╡ 55cdebd2-0881-11ef-2722-91de1447877a -using Markdown, InteractiveUtils - -# ╔═╡ a355b0ba-baaf-49f4-a5dc-965364a884f0 -using Catalyst, OrdinaryDiffEq - -# ╔═╡ 00fd6d49-f561-42e9-9413-d33af92f83dc -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 -""" - -# ╔═╡ 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 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} \, . -\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: -""" - -# ╔═╡ 935ca610-7a7a-4692-8908-fc26abb880b4 -fermenter_monod = @reaction_network begin - @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 -md""" -which resulted in the following differential equations: -""" - -# ╔═╡ fa93e2c3-8b43-418e-ba24-406645b2e397 -md""" -$$\begin{eqnarray*} -\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ -\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X -\end{eqnarray*}$$ -""" - -# ╔═╡ 06730f54-7293-43f2-b772-84eec3e5528a -md""" -Convert the system to a symbolic differential equation model and verify, by analyzing the differential equation, that your model is correctly implemented. - -Keep in mind that `mm(S, μmax, Ks)` stands for $\mu_{max} \, \cfrac{S}{S + K_s}$. -""" - -# ╔═╡ 7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd -# osys = missing -osys = convert(ODESystem, fermenter_monod) - -# ╔═╡ 1ee7359c-2757-44cd-8650-09b762cf30d0 -md""" -## Goals of this exercise -""" - -# ╔═╡ 55f1d688-0c53-481b-9965-5e92ca87ad83 -md""" -Compute the following in a timespan of `(0.0, 100.0)`$\,h$: - -- 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$ 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 `parms` with the parameter values: -""" - -# ╔═╡ 2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c -# 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 -tspan = (0.0, 100.0) - -# ╔═╡ 201dfb54-2056-4846-a64c-ff4951dc084d -md""" -For practical reasons, we will define the time step size as `dt = 0.5`. -""" - -# ╔═╡ 5627068f-6442-43bb-bacb-bcba917372f3 -# dt = missing -dt = 0.5 - -# ╔═╡ 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. -""" - -# ╔═╡ 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 -# 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=dt`, and store the solution in `osol`. Finally plot the results. -""" - -# ╔═╡ 8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 -# oprob = missing -oprob = ODEProblem(fermenter_monod, [], tspan, parms, combinatoric_ratelaws=false) - -# ╔═╡ 89a31c32-88a4-479f-a688-ffcb75ee8e91 -# osol = missing -osol = solve(oprob, Tsit5(), saveat=dt) - -# ╔═╡ 51a9b7e6-8ad9-477d-9596-ffd614df2c79 -# 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 -# function fermenter_monod_sim(parms) -# missing -# ... -# end -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 - -# ╔═╡ 4a5971b1-f4d0-43b6-805f-e17f5052ae92 -md""" -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(parms) = missing -fermenter_monod_sim_S(parms) = fermenter_monod_sim(parms)[:S] - -# ╔═╡ 3ae5bd00-2e06-4789-aab3-d897824d5e29 -# 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""" -Make the time vector. -""" - -# ╔═╡ fbd98975-aa32-46ae-8db0-0e65cdf48309 -# t_vals = missing -t_vals = tspan[1]:dt:tspan[2] - -# ╔═╡ 93791eb3-1eaa-4146-90b5-c4811fb3485b -md""" -Compute the two outputs $S$ and $X$ for the given parameter values. -""" - -# ╔═╡ dc0557d6-81b9-4759-8ed7-3129f60c6dc3 -# S_sim = missing -S_sim = fermenter_monod_sim_S([μmax_val, Ks_val, Sin_val]) - -# ╔═╡ 95bc683c-f6e6-4b42-b90b-b5a863edd4d5 -# 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 -sens_S = ForwardDiff.jacobian(fermenter_monod_sim_S, [μmax_val, Ks_val, Sin_val]) - -# ╔═╡ 49a94b9a-a543-495e-b4f1-c8579e59304d -# sens_X = missing -sens_X = ForwardDiff.jacobian(fermenter_monod_sim_X, [μmax_val, Ks_val, Sin_val]) - -# ╔═╡ 9cace6c1-e678-4dd7-8705-92a55eb32fa9 -md""" -Extract the (absolute) sensitivities of the outputs on the different parameters. -""" - -# ╔═╡ a6dc2b60-6a0a-4140-892e-02cde8dc79d3 -# begin -# sens_S_on_μmax = missing -# sens_S_on_Ks = missing -# sens_S_on_Sin = missing -# end -begin - sens_S_on_μmax = sens_S[:, 1] - sens_S_on_Ks = sens_S[:, 2] - sens_S_on_Sin = sens_S[:, 3] -end; - -# ╔═╡ f806c243-9032-46b7-add3-4714344691c7 -# begin -# sens_X_on_μmax = missing -# sens_X_on_Ks = missing -# sens_X_on_Sin = missing -# end -begin - sens_X_on_μmax = sens_X[:, 1] - sens_X_on_Ks = sens_X[:, 2] - sens_X_on_Sin = sens_X[:, 3] -end; - -# ╔═╡ f4d0aaaa-bb3f-4335-901a-dbbcde643c78 -md""" -#### Normalized sensitivities -""" - -# ╔═╡ 5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 -md""" -Compute the normalized sensitivities. -""" - -# ╔═╡ 76846731-929c-408f-a3de-970581c497e9 -# 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_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 -# 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_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""" -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 -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""" -!!! 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 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 -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=...`). -""" - -# ╔═╡ c0223da4-9959-48d0-b607-633b2e82986c -# 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""" -!!! 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 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 -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=...`). -""" - -# ╔═╡ 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 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""" -!!! 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 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 -# ╠═3ef93246-657d-4e77-9bf0-8380c64bfcfd -# ╠═a355b0ba-baaf-49f4-a5dc-965364a884f0 -# ╠═00fd6d49-f561-42e9-9413-d33af92f83dc -# ╠═7ae714c4-d25d-4f9f-ab3d-cc067db9c156 -# ╠═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 -# ╠═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_bitrophic_model_sol.jl b/exercises/solved_notebooks/P7_sens_uncert/uncert_bitrophic_model_sol.jl index c40b2a1c..4af5e5e7 100644 --- a/exercises/solved_notebooks/P7_sens_uncert/uncert_bitrophic_model_sol.jl +++ b/exercises/solved_notebooks/P7_sens_uncert/uncert_bitrophic_model_sol.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.46 +# v0.20.28 using Markdown using InteractiveUtils diff --git a/exercises/solved_notebooks/P7_sens_uncert/uncert_fermenter_monod_sol.jl b/exercises/solved_notebooks/P7_sens_uncert/uncert_fermenter_monod_sol.jl deleted file mode 100644 index fca64a97..00000000 --- a/exercises/solved_notebooks/P7_sens_uncert/uncert_fermenter_monod_sol.jl +++ /dev/null @@ -1,182 +0,0 @@ -### A Pluto.jl notebook ### -# v0.19.46 - -using Markdown -using InteractiveUtils - -# ╔═╡ 9d8acf40-e635-4dc4-9938-ec63ab68e3bd -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 - -# ╔═╡ 1ec7def2-0c43-11ef-0850-ddefa680d2a4 -using Markdown - -# ╔═╡ 2861bc17-7810-461e-8857-a55b6392977a -using InteractiveUtils - -# ╔═╡ 093679cc-122b-4b0a-84c3-15cf3a6f5c1b -using Catalyst - -# ╔═╡ 0720732e-105a-40c2-bf45-f8fbbc10e5bc -using OrdinaryDiffEq, StatsPlots - -# ╔═╡ 28bb60c9-ab56-497f-b039-0d6577501be5 -using Measurements - -# ╔═╡ 437364b6-a953-428b-a891-e17ffd531656 -md" -### Exercise: Fermenter - Monod kinetics - Uncertainty analysis -" - -# ╔═╡ c44e10f3-8452-4a60-b840-f4a7e6ec3238 -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: - -$$\begin{eqnarray*} -%S \xrightarrow[\quad\quad]{\beta} Y \, X -S \xrightarrow[\quad\quad]{r} Y \, X \quad\quad\quad\quad r = \mu \, X \quad \textrm{with} \quad \mu = \mu_{max} \, \cfrac{S}{S + K_s} -\end{eqnarray*}$$ -" - -# ╔═╡ 356822df-bad3-41f9-8556-de04234b0a09 -md" -The *reaction network object* for this model could be set-up as: -" - -# ╔═╡ c551c33e-95b0-4ea8-85d4-6c68db72b8bc -fermenter_monod = @reaction_network begin - X * mm(S, μmax, Ks), S => Y*X - Q/V, (S, X) --> 0 - Q/V*Sin, ∅ --> S -end - -# ╔═╡ 69f8c265-fe91-4e4b-971b-395325d4474e -md" -which resulted in the following differential equations: -" - -# ╔═╡ 2806654f-7a24-4f8e-8b5a-3129b2431c95 -md" -$$\begin{eqnarray*} -\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ -\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X -\end{eqnarray*}$$ -" - -# ╔═╡ 5d6a881d-421e-4dea-9668-38b3f3378bf1 -# osys = ... -osys = convert(ODESystem, fermenter_monod) - -# ╔═╡ 5802078d-12bc-489e-8566-15f77a11e8ba -md" -Assume the uncertainties in the following parameter values: - -- $\mu_{max} = 0.30 \pm 0.06\;h^{-1}$ -- $K_s = 0.15 \pm 0.03 \;g/L$ -- $S_{in} = 2.2 \pm 0.4\;g/L$ - -and that the uncertainty in the other parameters values: $Y = 0.80$, $Q = 2.0\;L/h$, $V = 40.0\;L$ are negligible. Suppose that at $t$ no substrate $S$ is present in the reactor but that there is initially some biomass with a concetration of $0.01\;g/L$. Perform an uncertainty analysis by plotting the uncertainty bands on the simulation results of $S$ and $X$ in a timespan of $[0, 100]\,h$. - -Interpret your results. -" - -# ╔═╡ 562e0006-bb79-4895-98b4-c376ea6ce856 -md" -Initialize a vector `u₀` with the initial conditions, and set the timespan: -" - -# ╔═╡ 6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde -# u₀ = missing # Uncomment and complete the instruction -# u0 = [:S => 0.0, :X => 0.01] -u0 = [:S => 0.0, :X => 0.0005] - -# ╔═╡ 1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b -# tspan = missing # Uncomment and complete the instruction -tspan = (0.0, 100.0) - -# ╔═╡ 40a8c697-56d1-4488-b089-00ddf41cdb5a -md" -We initialize a vector `params_uncert` with the parameter values and their corresponding uncertainty: -" - -# ╔═╡ 62f4049b-617f-4268-beba-b2c634b46642 -# params_uncert = missing # Uncomment and complete the instruction -# params_uncert = [:μmax => 0.30±0.06, :Ks => 0.15±0.03, :Y => 0.80, :Q => 2, :V => 40, :Sin => 2.2±0.4] -# params_uncert = [:μmax => 0.30±0.06, :Ks => 0.15, :Y => 0.80, :Q => 2, :V => 40, :Sin => 2.2] -# params_uncert = [:μmax => 0.30, :Ks => 0.15±0.03, :Y => 0.80, :Q => 2, :V => 40, :Sin => 2.2] -# params_uncert = [:μmax => 0.30, :Ks => 0.15, :Y => 0.80, :Q => 2, :V => 40, :Sin => 2.2±0.4] -# Use these: -params_uncert = [:μmax => 0.40±0.06, :Ks => 0.015±0.003, :Y => 0.67, :Q => 2, :V => 40, :Sin => 0.022±0.004] -# params_uncert = [:μmax => 0.40±0.06, :Ks => 0.015, :Y => 0.67, :Q => 2, :V => 40, :Sin => 0.022] -# params_uncert = [:μmax => 0.40, :Ks => 0.015±0.003, :Y => 0.67, :Q => 2, :V => 40, :Sin => 0.022] -# params_uncert = [:μmax => 0.40, :Ks => 0.015, :Y => 0.67, :Q => 2, :V => 40, :Sin => 0.022±0.004] - -# ╔═╡ 0e4c73d3-1630-463b-994d-073170c1a6c3 -md" -We create the corresponding ODE problem and store it in `oprob_uncert`: -" - -# ╔═╡ 3cad11d6-c2eb-4865-858a-d407aa25d780 -# oprob_uncert = missing # Uncomment and complete the instruction -oprob_uncert = ODEProblem(fermenter_monod, u0, tspan, params_uncert, combinatoric_ratelaws=false) - -# ╔═╡ a7b3b68f-d34e-4e9a-aa44-224919c44fec -md" -We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol_uncert`: -" - -# ╔═╡ 40168902-0780-4f92-8eca-23982c80511e -# osol_uncert = missing # Uncomment and complete the instruction -osol_uncert = solve(oprob_uncert, Tsit5(), saveat=2.0) - -# ╔═╡ 6b2abf2a-f1f3-4c2d-96e8-cc626848fa68 -md" -Plot the results (simulation of the output variables $S$ and $X$ together with their uncertainty band): -" - -# ╔═╡ 4d753fdb-089d-41dc-89ef-45a1733af786 -plot(osol_uncert) - -# ╔═╡ 4f3900d7-4ff8-4b4b-a59a-f7079ac6c457 -md" -Try to relate the local sensitivity analysis to the uncertainty analysis. Hence, study the effect of the individual parameter uncertainties on the output variables $S$ and $X$ and compare with your local sensitivity results of the corresponding parameter. - -In order to do that, analyse the effect on the uncertainty bands for $S$ and $X$ by taking one uncertainty on a parameter at a time. In other words, analyse the subsequent cases seperately: -- Assume uncertainty only in $\mu_{max}$ -- Assume uncertainty only in $K_s$ -- Assume uncertainty only in $S_{in}$ - -Draw your conclusions: -- Answer: missing -" - -# ╔═╡ Cell order: -# ╠═1ec7def2-0c43-11ef-0850-ddefa680d2a4 -# ╠═2861bc17-7810-461e-8857-a55b6392977a -# ╠═9d8acf40-e635-4dc4-9938-ec63ab68e3bd -# ╠═093679cc-122b-4b0a-84c3-15cf3a6f5c1b -# ╠═0720732e-105a-40c2-bf45-f8fbbc10e5bc -# ╠═28bb60c9-ab56-497f-b039-0d6577501be5 -# ╠═437364b6-a953-428b-a891-e17ffd531656 -# ╠═c44e10f3-8452-4a60-b840-f4a7e6ec3238 -# ╠═356822df-bad3-41f9-8556-de04234b0a09 -# ╠═c551c33e-95b0-4ea8-85d4-6c68db72b8bc -# ╠═69f8c265-fe91-4e4b-971b-395325d4474e -# ╠═2806654f-7a24-4f8e-8b5a-3129b2431c95 -# ╠═5d6a881d-421e-4dea-9668-38b3f3378bf1 -# ╠═5802078d-12bc-489e-8566-15f77a11e8ba -# ╠═562e0006-bb79-4895-98b4-c376ea6ce856 -# ╠═6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde -# ╠═1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b -# ╠═40a8c697-56d1-4488-b089-00ddf41cdb5a -# ╠═62f4049b-617f-4268-beba-b2c634b46642 -# ╠═0e4c73d3-1630-463b-994d-073170c1a6c3 -# ╠═3cad11d6-c2eb-4865-858a-d407aa25d780 -# ╠═a7b3b68f-d34e-4e9a-aa44-224919c44fec -# ╠═40168902-0780-4f92-8eca-23982c80511e -# ╠═6b2abf2a-f1f3-4c2d-96e8-cc626848fa68 -# ╠═4d753fdb-089d-41dc-89ef-45a1733af786 -# ╠═4f3900d7-4ff8-4b4b-a59a-f7079ac6c457 diff --git a/exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl b/exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl index 12b85542..4d7ee999 100644 --- a/exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl +++ b/exercises/solved_notebooks/P7_sens_uncert/uncert_intro_sol.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.21 +# v1.0.3 using Markdown using InteractiveUtils @@ -199,27 +199,11 @@ In uncertainty propagation, we are interested in the effect of uncertainties (or 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. @@ -233,6 +217,22 @@ Firstly, solve the original system (Wf = 10 and μ = 0.07) and afterwards we can # ╔═╡ 62afc495-29c3-4baa-b45a-94b84149d15d oprob_log = ODEProblem(growth_mod_log, u0_log,tspan, [:μ => 0.07, :Wf => 10]); +# ╔═╡ 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 + +# ╔═╡ 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); + # ╔═╡ 691a8b68-935e-4012-aa25-9cfdc323b1de osol_log = solve(oprob_log, Tsit5(), saveat = 0.5); From 049091e9db7b77ff14290a730e93aa5ee450cdee Mon Sep 17 00:00:00 2001 From: wclaerhout Date: Tue, 20 Oct 2026 14:39:47 +0200 Subject: [PATCH 2/2] Update P7 student notebooks to match the merged fermenter solution - Add student version of the fermenter mutant notebook (sensitivity + uncertainty merged) and remove the old sens/uncert fermenter notebooks - Apply the uncert_intro cell reordering to the student notebook - Rename solved fermenter_mutant.jl to fermenter_mutant_sol.jl - Update the website frontmatter table and regenerate src/exercises Co-Authored-By: Claude Opus 5.5 --- ...nter_mutant.jl => fermenter_mutant_sol.jl} | 0 .../P7_sens_uncert/fermenter_mutant.jl | 595 +++++++++++++++ .../P7_sens_uncert/sens_fermenter_monod.jl | 666 ----------------- .../P7_sens_uncert/uncert_fermenter_monod.jl | 174 ----- .../P7_sens_uncert/uncert_intro.jl | 34 +- notebook-checks/sync_exercises.jl | 3 +- src/exercises/fermenter_mutant.jl | 603 ++++++++++++++++ src/exercises/sens_fermenter_monod.jl | 678 ------------------ src/exercises/uncert_fermenter_monod.jl | 182 ----- src/exercises/uncert_intro.jl | 34 +- 10 files changed, 1233 insertions(+), 1736 deletions(-) rename exercises/solved_notebooks/P7_sens_uncert/{fermenter_mutant.jl => fermenter_mutant_sol.jl} (100%) create mode 100644 exercises/student_notebooks/P7_sens_uncert/fermenter_mutant.jl delete mode 100644 exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl delete mode 100644 exercises/student_notebooks/P7_sens_uncert/uncert_fermenter_monod.jl create mode 100644 src/exercises/fermenter_mutant.jl delete mode 100644 src/exercises/sens_fermenter_monod.jl delete mode 100644 src/exercises/uncert_fermenter_monod.jl diff --git a/exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant.jl b/exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant_sol.jl similarity index 100% rename from exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant.jl rename to exercises/solved_notebooks/P7_sens_uncert/fermenter_mutant_sol.jl diff --git a/exercises/student_notebooks/P7_sens_uncert/fermenter_mutant.jl b/exercises/student_notebooks/P7_sens_uncert/fermenter_mutant.jl new file mode 100644 index 00000000..3489209a --- /dev/null +++ b/exercises/student_notebooks/P7_sens_uncert/fermenter_mutant.jl @@ -0,0 +1,595 @@ +### A Pluto.jl notebook ### +# v0.20.28 + +using Markdown +using InteractiveUtils + +# ╔═╡ 9d8acf40-e635-4dc4-9938-ec63ab68e3bd +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 + +# ╔═╡ 1ec7def2-0c43-11ef-0850-ddefa680d2a4 +using Markdown + +# ╔═╡ 2861bc17-7810-461e-8857-a55b6392977a +using InteractiveUtils + +# ╔═╡ 093679cc-122b-4b0a-84c3-15cf3a6f5c1b +using Catalyst + +# ╔═╡ 0720732e-105a-40c2-bf45-f8fbbc10e5bc +using OrdinaryDiffEq + +# ╔═╡ 28bb60c9-ab56-497f-b039-0d6577501be5 +using StatsPlots, PlutoUI; TableOfContents() + +# ╔═╡ 0a65afbc-643a-4fe4-b22d-c5fe54e3be1f +using Turing + +# ╔═╡ d1b6739f-1e0c-452d-8763-717dd650b2ec +using ForwardDiff + +# ╔═╡ 162b330b-02f7-4179-b950-d960eca68206 +using Random + +# ╔═╡ 437364b6-a953-428b-a891-e17ffd531656 +md""" +# Exercise: Fermenter - Switching to a mutant strain? +""" + +# ╔═╡ c44e10f3-8452-4a60-b840-f4a7e6ec3238 +md""" +As a process engineer, you are responsible for a continuously stirred 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$ [$L/h$] containing substrate at a concentration $S_{in}$ [$g/L$]. The yield $Y$ [$g$ biomass$/g$ substrate] tells how much biomass is formed per gram of substrate consumed. As in the previous practica, the growth is modelled with Monod kinetics: + +$$\begin{eqnarray*} +S \xrightarrow[\quad\quad]{r} Y \, X \quad\quad\quad\quad r = \mu \, X \quad \textrm{with} \quad \mu = \mu_{max} \, \cfrac{S}{S + K_s} +\end{eqnarray*}$$ +""" + +# ╔═╡ f3065171-5526-476a-bfd7-a2b95bb24d61 +md""" +## The question + +A colleague has found a mutant strain with a **higher yield** $Y$, and management asks you a simple question: *should we switch to the mutant?* + +The catch: the maximal growth rate $\mu_{max}$ of the mutant is not known. A higher yield often comes at the cost of slower growth, and a slower grower leaves more substrate unused in the outflow. So the mutant could be better or worse than the current strain. + +You will answer this question with the tools of this practical: +1. **Uncertainty propagation**: what does our uncertainty about $\mu_{max}$ mean for the biomass? How sure are we that the mutant is better? +2. **Sensitivity analysis**: when does the biomass depend most strongly on $\mu_{max}$? That is the moment where a measurement tells us the most. +3. **Bayesian inference**: update our belief about $\mu_{max}$ with a few measurements, and propagate the reduced uncertainty again. +""" + +# ╔═╡ 356822df-bad3-41f9-8556-de04234b0a09 +md""" +The *reaction network object* for this model is: +""" + +# ╔═╡ c551c33e-95b0-4ea8-85d4-6c68db72b8bc +fermenter_monod = @reaction_network begin + X * mm(S, μmax, Ks), S => Y*X + Q/V, (S, X) --> 0 + Q/V*Sin, ∅ --> S +end + +# ╔═╡ 69f8c265-fe91-4e4b-971b-395325d4474e +md""" +which results in the following differential equations: +""" + +# ╔═╡ 2806654f-7a24-4f8e-8b5a-3129b2431c95 +md""" +$$\begin{eqnarray*} +\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ +\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X +\end{eqnarray*}$$ +""" + +# ╔═╡ 5d6a881d-421e-4dea-9668-38b3f3378bf1 +# osys = missing + +# ╔═╡ 5802078d-12bc-489e-8566-15f77a11e8ba +md""" +## The current strain + +For the strain currently in use, the parameters are well known: +- $\mu_{max} = 0.40\;h^{-1}$ +- $K_s = 0.015 \;g/L$ +- $Y = 0.67$ +- $S_{in} = 0.022 \;g/L$ +- $Q = 2.0\;L/h$ +- $V = 40.0\;L$ + +At $t=0$ no substrate is present in the reactor and the biomass concentration is $0.0005\;g/L$. We judge a strain by the biomass concentration in the outflow after $100\,h$. + +Initialize a vector `u0` with the initial conditions, set the timespan `tspan` to $[0, 100]\,h$, the time step `dt` to `0.5`, and initialize a vector `parms` with the parameter values. +""" + +# ╔═╡ 6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde +# u0 = missing + +# ╔═╡ 1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b +# tspan = missing + +# ╔═╡ 9e3c427b-8a56-4787-9b28-43ec3932ddbe +# dt = missing + +# ╔═╡ 62f4049b-617f-4268-beba-b2c634b46642 +# parms = missing + +# ╔═╡ 0e4c73d3-1630-463b-994d-073170c1a6c3 +md""" +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. +""" + +# ╔═╡ 3cad11d6-c2eb-4865-858a-d407aa25d780 +# oprob = missing + +# ╔═╡ 40168902-0780-4f92-8eca-23982c80511e +# osol = missing + +# ╔═╡ 4d753fdb-089d-41dc-89ef-45a1733af786 +# missing + +# ╔═╡ ca660e2c-9348-40b1-b92d-415068bea514 +md""" +Store the biomass concentration at the end of the simulation in `X_current`. This is the value the mutant has to beat. +""" + +# ╔═╡ ec042e61-04b8-4d51-8e7d-35e1646cc5ba +# X_current = missing + +# ╔═╡ 9ee22028-7a37-4510-9b61-198d3bd99a73 +md""" +## The mutant strain +""" + +# ╔═╡ cc9ac43b-448a-4a67-ba0b-8fdf21a4236e +md""" +Some batch tests show that the mutant has a yield of $Y = 0.70$ instead of $0.67$. All other parameters are the same as for the current strain, except $\mu_{max}$, which is unknown. + +From the literature on mutants with an improved yield, you expect $\mu_{max}$ to be around $0.35\;h^{-1}$, give or take $15\%$. You express this **prior belief** with a log-normal distribution `LogNormal(log(0.35), 0.15)`. +""" + +# ╔═╡ f5482531-b23f-472b-968f-20f7a6471778 +begin + plot(LogNormal(log(0.35), 0.15), xlabel="μmax (1/h)", ylabel="density", label="Prior μmax mutant", linewidth=2) + vline!([0.40], label="μmax current strain", linestyle=:dash) +end + +# ╔═╡ ff21b3cd-64b1-4d09-8eb1-afc617585fc8 +md""" +Create an ODE problem `oprob_mut` for the mutant by remaking `oprob` with $Y = 0.70$ and the expected value $\mu_{max} = 0.35$. + +!!! note + Remember that you can remake an ODEProblem with new parameter values by using `remake(oprob, p = [...])`. +""" + +# ╔═╡ 225eaaff-42a2-4aca-94e1-afcd988de4eb +# oprob_mut = missing + +# ╔═╡ 18ffe877-c579-4b4c-80bd-b5d50baf36bb +md""" +## Monte Carlo error propagation + +We don't know $\mu_{max}$, but we do know which values are plausible. Using Monte Carlo simulation, we draw many values from the prior, simulate the mutant for each value, and look at the spread of the outputs. This shows how the uncertainty on the parameter **propagates** to the uncertainty on the output. + +Make a Turing model `mutant_deviation` in which you implement the prior for `μmax_dev` and return the solution of the remade ODE problem. +""" + +# ╔═╡ 7821dc21-3e6b-4ad1-829e-f29dcd020144 +# @model function mutant_deviation() +# μmax_dev ~ missing +# new_prob = missing +# return missing +# end + +# ╔═╡ 3051d39a-fc2f-40e4-a688-4d6e75cfc2d1 +md""" +Now sample 500 solutions from your `mutant_deviation` model. +""" + +# ╔═╡ 33df05c9-08d3-4107-9240-931658cda2a4 +# missing + +# ╔═╡ 02c5cf65-ebd3-477b-9a3e-8855bd46f55c +md""" +Plot the biomass of all Monte Carlo simulations together with the biomass of the current strain. Use `idxs=:X`, `label=false`, `color=:gray`, `alpha=0.4` and `linestyle=:dot`. +""" + +# ╔═╡ aa833d92-567f-4971-a360-e60bce424c8e +# begin +# p1 = plot(title="Monte Carlo of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") +# for missing in missing +# missing +# end +# plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) +# p1 +# end + +# ╔═╡ ffc66a94-5435-4ca6-af87-76c85a599906 +md""" +Create a vector `x1` that stores the biomass at the end of the simulations, make a histogram of it, and compute the probability `P_better_prior` that the mutant produces more biomass than the current strain. + +!!! note + The fraction of elements of a vector `x` that are larger than `a` can be computed with `mean(x .> a)`. +""" + +# ╔═╡ d3afd5b4-034f-415c-a64a-776deed74d5c +# x1 = missing + +# ╔═╡ 1d95728e-079e-4e9b-952a-c4f692661813 +# missing + +# ╔═╡ 02213469-a233-4759-a18d-29a76a4c13c6 +# P_better_prior = missing + +# ╔═╡ 97d2065c-6234-4b13-8e4c-6782b8301416 +md""" +!!! questions + - When is the spread on $X$ largest, and when is it smallest? Can you explain this? + - How likely is it that the mutant beats the current strain? Would you advise management to switch based on this? +""" + +# ╔═╡ d2acbfb9-a7ff-48b0-9d9a-9492ecb9c1bf +md""" +Answers: +- missing +- missing +""" + +# ╔═╡ eb613ead-6aea-41fe-acb7-f679b5a97c35 +md""" +## Local Sensitivity Analysis (LSA) + +The lab can run the mutant in a pilot reactor under the same conditions, but you can only afford **three measurements** of $X$, all at the same time point. Measurements are expensive, so we want them at the moment where $X$ depends most strongly on $\mu_{max}$. There, a small change in $\mu_{max}$ causes a large change in $X$. The other way around, a measurement of $X$ at that moment pins down $\mu_{max}$ best. + +So, we will first perform a normalized local sensitivity analysis to pin down the moment, $X$ is most sensitive to a change in $μmax$. +""" + +# ╔═╡ 67f572b0-3485-40bc-8f64-cc58471b38a9 +md""" +Write a solution function `mutant_sim` with as argument a vector with the parameter `μmax` of the mutant, and that returns the solution. Based on it, write a function `mutant_sim_X` that only returns the output $X$. +""" + +# ╔═╡ 565b6e18-012f-4b38-bf21-bdc683688107 +# function mutant_sim(parms) +# missing +# ... +# end + +# ╔═╡ b7f14c48-ba19-4fbc-bded-829a34bd7ac8 +# mutant_sim_X(parms) = missing + +# ╔═╡ 73345fe4-6877-476d-8ce1-5a0056c73cc6 +md""" +Make the time vector `t_vals` with timestep `dt` equal to $0.5$ and compute the output `X_sim` for the expected value $\mu_{max} = 0.35$, stored in `μmax_val`. +""" + +# ╔═╡ 75ab68b2-1e88-421f-a089-4c0b3b72788d +# t_vals = missing + +# ╔═╡ 82ae802b-2a7a-4525-9cc6-dfaa552ba990 +# μmax_val = missing + +# ╔═╡ 8a655d11-0d45-4e7c-bcbf-8ea9659516aa +# X_sim = missing + +# ╔═╡ 9a97a7c7-bc3e-434e-9dce-7c859417aac2 +md""" +Use `ForwardDiff.jacobian` to compute the sensitivity `sens_X_on_μmax`, and compute the normalized sensitivity `sens_X_on_μmax_rel`. +""" + +# ╔═╡ d7d3028a-51c1-41e0-a865-b521ed235548 +# sens_X_on_μmax = missing + +# ╔═╡ cb6047f2-3884-4fbd-96f2-f15cf684bea8 +# sens_X_on_μmax_rel = missing + +# ╔═╡ f623b4b4-6010-45fa-b9a9-bc796f0fb106 +md""" +Identify the time `t_star` at which the biomass is most sensitive to variations in $\mu_{max}$, and plot the normalized sensitivity function together with `t_star`. + +!!! 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. +""" + +# ╔═╡ 9c3fc6d8-8e5e-425d-9caa-56a9e0b03b99 +# begin +# t_star_idx = missing +# t_star = t_vals[t_star_idx] +# +# plot(t_vals, abs.(sens_X_on_μmax_rel), +# title="Sensitivity of X w.r.t. μmax", xlabel="Time (hours)", +# ylabel="normalized sensitivity", label=false, linewidth=2) +# vline!([t_star], label="t* = $(t_star) h", linestyle=:dash, color=:red) +# end + +# ╔═╡ b1825a21-e711-4f99-9a4b-51f3fb73ff82 +md""" +!!! questions + - Compare the sensitivity function with the Monte Carlo simulations above. What do you notice? + - The sensitivity at $t = 100\,h$ is small, while that is exactly the time we care about. Why is it small? Is it a good idea to measure at $t = 100\,h$? +""" + +# ╔═╡ 44234af3-6015-49b6-8709-1f3a3cc0f69d +md""" +Answers: +- missing +- missing +""" + +# ╔═╡ a00afabe-7bc1-4d8b-9b55-2953c7add9f5 +md""" +## Measuring campaign + +The lab ran the mutant in the pilot reactor and sampled it at your chosen time. We have defined a function `measure_X(t)` in this notebook that returns three hypothetical measurements of $X$ at the given timepoint `t`, with a relative measurement error of $5\%$. Use this function to get the measurements `X_meas` at the optimal time `t_star`. +""" + +# ╔═╡ 635ff400-639c-4c80-93e6-e8721a836c16 +begin + μmax_real = 0.37 + real_solution = solve(remake(oprob_mut, p = [:μmax => μmax_real]), Tsit5(), saveat=dt) + measure_X(t) = real_solution(t, idxs=:X) .* (1 .+ 0.05 .* randn(Xoshiro(round(Int, 10t)), 3)) +end; + +# ╔═╡ 93cc5a4f-f40a-4d28-abc3-8d6790fd60dc +# X_meas = missing + +# ╔═╡ 3e93a001-9401-4a1f-b283-c79fd0c43434 +md""" +### Conditioning on the measurements +""" + +# ╔═╡ 04ae6cd7-2d59-4bc0-a894-d9d5129e243d +md""" +Based on these measurements we can update what we believed before (the prior) with what the data tell us (the likelihood) into an updated belief (the posterior): + +Create a Turing model `mutant_meas(t_meas)` with: +- the same prior for `μmax_dev` as before; +- a prior `Exponential()` for the measurement error `σ_X`; +- the simulated biomass `X_pred` at time `t_meas`, which you can get with `sol(t_meas, idxs=:X)`; +- the three measurements `X_meas`, each normally distributed around `X_pred` with standard deviation `σ_X`; +- the solution as return value. +""" + +# ╔═╡ 02096fef-224e-4896-9f8d-3f7aba7334e1 +# @model function mutant_meas(t_meas) +# μmax_dev ~ missing +# σ_X ~ missing +# +# sol = missing +# X_pred = missing +# +# X_meas = zeros(3) +# for i in 1:length(X_meas) +# X_meas[i] ~ missing +# end +# return sol +# end + +# ╔═╡ 354f0768-efc0-4304-aa04-c6a68f961aa6 +md""" +Condition the model on the measurements and store it in `mutant_cond_model`. +""" + +# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 +# mutant_cond_model = missing + +# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 +md""" +### MCMC +""" + +# ╔═╡ f99ca771-8c14-4e19-a1ae-478a2103716a +# mutant_chain = missing + +# ╔═╡ e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +md""" +Plot the chain to check if it properly converged to the posterior distribution. +""" + +# ╔═╡ 50b30e65-72d3-4d8d-a739-ae969e0f50a8 +# missing + +# ╔═╡ a65710f6-70ce-4348-923d-ed2ed2b4a469 +md""" +Now, get the sampled posterior solutions `solutions_mutant` and plot them on top of the prior Monte Carlo simulations, together with the current strain and the measurements. +""" + +# ╔═╡ 12c9f810-ef5b-49a4-8fdd-6efb1a040440 +# solutions_mutant = missing + +# ╔═╡ e2e11092-4498-437b-8c67-8d77302893fe +md""" +The posterior chance that the mutant is better is: +""" + +# ╔═╡ cae4d9c4-3790-423e-8321-d44c0833aa40 +begin + sol_end = [solutions_mutant[i][:X][end] for i in eachindex(solutions_mutant)] + P_better_posterior = mean(sol_end .> X_current) + println("The chance the mutant is better is $(P_better_posterior*100)%") +end + +# ╔═╡ 49307c36-3518-4118-b453-91e634f92daf +# begin +# p2 = plot(title="Prior vs posterior of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") +# for missing in missing +# missing +# end +# for missing in missing +# missing +# end +# plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) +# scatter!(fill(t_star, 3), X_meas, label="Measurements", color=:black) +# p2 +# end + +# ╔═╡ aa479bcd-2314-47a5-83ba-b0d1a2654fa0 +md""" +Create a vector `x2` with the posterior biomass at the end of the simulations, compare its histogram with the prior one, and compute the probability `P_better_post` that the mutant beats the current strain. +""" + +# ╔═╡ 2f55a771-0d37-48e5-a9d9-24693bf1d177 +# x2 = missing + +# ╔═╡ 18059d0c-cf37-44a1-86a8-859ac66a66f0 +# missing + +# ╔═╡ 40e73ddf-6c0d-4f83-860f-434d934193ed +# P_better_post = missing + +# ╔═╡ e1c87260-25ee-4d71-8744-ed0d7cb75f6a +md""" +### What if we had measured at the end? + +Intuitively, you might have measured $X$ at $t = 100\,h$, since that is the quantity we care about. Repeat the inference with three measurements at $t = 100\,h$, and compute the probability `P_better_100` that the mutant is better. +""" + +# ╔═╡ caa5c04d-49af-4576-83a2-605bae1c2494 +# X_meas_100 = missing + +# ╔═╡ 6a53aff6-29d4-4d41-8236-f621ce519361 +# missing + +# ╔═╡ 9292c593-1f0f-47a6-a2ec-91a8aadbe6e6 +# P_better_100 = missing + +# ╔═╡ b1814ea6-0065-4ac2-a012-505d2bbd8e40 +md""" +Compare what we know about $\mu_{max}$ before and after both measuring campaigns: +""" + +# ╔═╡ 44cb7615-86b2-45f0-b87a-69fcbfe863bb +begin + density(vec(chain[:μmax_dev]), label="Prior", linewidth=2, title="What did we learn about μmax?", xlabel="μmax (1/h)", ylabel="density") + density!(vec(mutant_chain_100[:μmax_dev]), label="Posterior, measured at t = 100 h", linewidth=2) + density!(vec(mutant_chain[:μmax_dev]), label="Posterior, measured at t* = $(t_star) h", linewidth=2) +end + +# ╔═╡ f364fa1f-5c47-4ca6-99df-17e2b735d98b +md""" +!!! questions + - Compare `P_better_prior`, `P_better_100` and `P_better_post`. +""" + +# ╔═╡ 5f9873a5-af70-4ac2-be08-4f249b71eb76 +md""" +Answers: +- missing +- missing +""" + +# ╔═╡ 3980da6b-7036-44d2-9d8f-3ddc1062d0ad +md""" +## Conclusions + +!!! questions + - What do you advise management? Support your advice with numbers. + - What did you learn/what did you use the information for from uncertainty propagation, sensitivity analysis and the Bayesian inference. +""" + +# ╔═╡ 8021123d-1b19-47bf-9db0-0d7c587bc00e +md""" +Answers: +- missing +- missing +- missing +""" + +# ╔═╡ Cell order: +# ╠═1ec7def2-0c43-11ef-0850-ddefa680d2a4 +# ╠═2861bc17-7810-461e-8857-a55b6392977a +# ╠═9d8acf40-e635-4dc4-9938-ec63ab68e3bd +# ╠═093679cc-122b-4b0a-84c3-15cf3a6f5c1b +# ╠═0720732e-105a-40c2-bf45-f8fbbc10e5bc +# ╠═28bb60c9-ab56-497f-b039-0d6577501be5 +# ╠═0a65afbc-643a-4fe4-b22d-c5fe54e3be1f +# ╠═d1b6739f-1e0c-452d-8763-717dd650b2ec +# ╠═162b330b-02f7-4179-b950-d960eca68206 +# ╟─437364b6-a953-428b-a891-e17ffd531656 +# ╟─c44e10f3-8452-4a60-b840-f4a7e6ec3238 +# ╟─f3065171-5526-476a-bfd7-a2b95bb24d61 +# ╟─356822df-bad3-41f9-8556-de04234b0a09 +# ╠═c551c33e-95b0-4ea8-85d4-6c68db72b8bc +# ╟─69f8c265-fe91-4e4b-971b-395325d4474e +# ╟─2806654f-7a24-4f8e-8b5a-3129b2431c95 +# ╠═5d6a881d-421e-4dea-9668-38b3f3378bf1 +# ╟─5802078d-12bc-489e-8566-15f77a11e8ba +# ╠═6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde +# ╠═1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b +# ╠═9e3c427b-8a56-4787-9b28-43ec3932ddbe +# ╠═62f4049b-617f-4268-beba-b2c634b46642 +# ╟─0e4c73d3-1630-463b-994d-073170c1a6c3 +# ╠═3cad11d6-c2eb-4865-858a-d407aa25d780 +# ╠═40168902-0780-4f92-8eca-23982c80511e +# ╠═4d753fdb-089d-41dc-89ef-45a1733af786 +# ╟─ca660e2c-9348-40b1-b92d-415068bea514 +# ╠═ec042e61-04b8-4d51-8e7d-35e1646cc5ba +# ╟─9ee22028-7a37-4510-9b61-198d3bd99a73 +# ╟─cc9ac43b-448a-4a67-ba0b-8fdf21a4236e +# ╟─f5482531-b23f-472b-968f-20f7a6471778 +# ╟─ff21b3cd-64b1-4d09-8eb1-afc617585fc8 +# ╠═225eaaff-42a2-4aca-94e1-afcd988de4eb +# ╟─18ffe877-c579-4b4c-80bd-b5d50baf36bb +# ╠═7821dc21-3e6b-4ad1-829e-f29dcd020144 +# ╟─3051d39a-fc2f-40e4-a688-4d6e75cfc2d1 +# ╠═33df05c9-08d3-4107-9240-931658cda2a4 +# ╟─02c5cf65-ebd3-477b-9a3e-8855bd46f55c +# ╠═aa833d92-567f-4971-a360-e60bce424c8e +# ╟─ffc66a94-5435-4ca6-af87-76c85a599906 +# ╠═d3afd5b4-034f-415c-a64a-776deed74d5c +# ╠═1d95728e-079e-4e9b-952a-c4f692661813 +# ╠═02213469-a233-4759-a18d-29a76a4c13c6 +# ╟─97d2065c-6234-4b13-8e4c-6782b8301416 +# ╟─d2acbfb9-a7ff-48b0-9d9a-9492ecb9c1bf +# ╟─eb613ead-6aea-41fe-acb7-f679b5a97c35 +# ╟─67f572b0-3485-40bc-8f64-cc58471b38a9 +# ╠═565b6e18-012f-4b38-bf21-bdc683688107 +# ╠═b7f14c48-ba19-4fbc-bded-829a34bd7ac8 +# ╟─73345fe4-6877-476d-8ce1-5a0056c73cc6 +# ╠═75ab68b2-1e88-421f-a089-4c0b3b72788d +# ╠═82ae802b-2a7a-4525-9cc6-dfaa552ba990 +# ╠═8a655d11-0d45-4e7c-bcbf-8ea9659516aa +# ╟─9a97a7c7-bc3e-434e-9dce-7c859417aac2 +# ╠═d7d3028a-51c1-41e0-a865-b521ed235548 +# ╠═cb6047f2-3884-4fbd-96f2-f15cf684bea8 +# ╟─f623b4b4-6010-45fa-b9a9-bc796f0fb106 +# ╠═9c3fc6d8-8e5e-425d-9caa-56a9e0b03b99 +# ╟─b1825a21-e711-4f99-9a4b-51f3fb73ff82 +# ╟─44234af3-6015-49b6-8709-1f3a3cc0f69d +# ╟─a00afabe-7bc1-4d8b-9b55-2953c7add9f5 +# ╟─635ff400-639c-4c80-93e6-e8721a836c16 +# ╠═93cc5a4f-f40a-4d28-abc3-8d6790fd60dc +# ╟─3e93a001-9401-4a1f-b283-c79fd0c43434 +# ╟─04ae6cd7-2d59-4bc0-a894-d9d5129e243d +# ╠═02096fef-224e-4896-9f8d-3f7aba7334e1 +# ╟─354f0768-efc0-4304-aa04-c6a68f961aa6 +# ╠═11e105de-270c-434a-bfb4-9118dfc9e461 +# ╟─e46afb98-16de-4d19-86b4-b920b6b929c5 +# ╠═f99ca771-8c14-4e19-a1ae-478a2103716a +# ╟─e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +# ╠═50b30e65-72d3-4d8d-a739-ae969e0f50a8 +# ╟─a65710f6-70ce-4348-923d-ed2ed2b4a469 +# ╠═12c9f810-ef5b-49a4-8fdd-6efb1a040440 +# ╟─e2e11092-4498-437b-8c67-8d77302893fe +# ╠═cae4d9c4-3790-423e-8321-d44c0833aa40 +# ╠═49307c36-3518-4118-b453-91e634f92daf +# ╟─aa479bcd-2314-47a5-83ba-b0d1a2654fa0 +# ╠═2f55a771-0d37-48e5-a9d9-24693bf1d177 +# ╠═18059d0c-cf37-44a1-86a8-859ac66a66f0 +# ╠═40e73ddf-6c0d-4f83-860f-434d934193ed +# ╟─e1c87260-25ee-4d71-8744-ed0d7cb75f6a +# ╠═caa5c04d-49af-4576-83a2-605bae1c2494 +# ╠═6a53aff6-29d4-4d41-8236-f621ce519361 +# ╠═9292c593-1f0f-47a6-a2ec-91a8aadbe6e6 +# ╟─b1814ea6-0065-4ac2-a012-505d2bbd8e40 +# ╟─44cb7615-86b2-45f0-b87a-69fcbfe863bb +# ╟─f364fa1f-5c47-4ca6-99df-17e2b735d98b +# ╟─5f9873a5-af70-4ac2-be08-4f249b71eb76 +# ╟─3980da6b-7036-44d2-9d8f-3ddc1062d0ad +# ╟─8021123d-1b19-47bf-9db0-0d7c587bc00e diff --git a/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl b/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl deleted file mode 100644 index 34d954b0..00000000 --- a/exercises/student_notebooks/P7_sens_uncert/sens_fermenter_monod.jl +++ /dev/null @@ -1,666 +0,0 @@ -### A Pluto.jl notebook ### -# v0.20.21 - -using Markdown -using InteractiveUtils - -# ╔═╡ 3ef93246-657d-4e77-9bf0-8380c64bfcfd -using Pkg; Pkg.activate("..") - -# ╔═╡ 55cdebd2-0881-11ef-2722-91de1447877a -using Markdown, InteractiveUtils - -# ╔═╡ a355b0ba-baaf-49f4-a5dc-965364a884f0 -using Catalyst, OrdinaryDiffEq - -# ╔═╡ 00fd6d49-f561-42e9-9413-d33af92f83dc -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 -""" - -# ╔═╡ 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 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} \, . -\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: -""" - -# ╔═╡ 935ca610-7a7a-4692-8908-fc26abb880b4 -fermenter_monod = @reaction_network begin - @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 -md""" -which resulted in the following differential equations: -""" - -# ╔═╡ fa93e2c3-8b43-418e-ba24-406645b2e397 -md""" -$$\begin{eqnarray*} -\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ -\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X -\end{eqnarray*}$$ -""" - -# ╔═╡ 06730f54-7293-43f2-b772-84eec3e5528a -md""" -Convert the system to a symbolic differential equation model and verify, by analyzing the differential equation, that your model is correctly implemented. - -Keep in mind that `mm(S, μmax, Ks)` stands for $\mu_{max} \, \cfrac{S}{S + K_s}$. -""" - -# ╔═╡ 7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd -# osys = missing - -# ╔═╡ 1ee7359c-2757-44cd-8650-09b762cf30d0 -md""" -## Goals of this exercise -""" - -# ╔═╡ 55f1d688-0c53-481b-9965-5e92ca87ad83 -md""" -Compute the following in a timespan of `(0.0, 100.0)`$\,h$: - -- 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$ 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 `parms` with the parameter values: -""" - -# ╔═╡ 2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c -# u0 = missing # in principle not necessary since we will use the default - -# ╔═╡ 3fdc6b17-cdeb-4dc5-8886-9d3a62caac8d -# tspan = missing - -# ╔═╡ 201dfb54-2056-4846-a64c-ff4951dc084d -md""" -For practical reasons, we will define the time step size as `dt = 0.5`. -""" - -# ╔═╡ 5627068f-6442-43bb-bacb-bcba917372f3 -# dt = missing - -# ╔═╡ 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. -""" - -# ╔═╡ 0f995929-4d2b-4a7a-8da1-04e4d501385f -begin - # μmax_val = missing - # Ks_val = missing - # Sin_val = missing -end; - -# ╔═╡ 79b0eb65-5a0f-40b3-aa97-4088421c562e -# 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=dt`, and store the solution in `osol`. Finally plot the results. -""" - -# ╔═╡ 8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 -# oprob = missing - -# ╔═╡ 89a31c32-88a4-479f-a688-ffcb75ee8e91 -# osol = missing - -# ╔═╡ 51a9b7e6-8ad9-477d-9596-ffd614df2c79 -# 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 (that you want the sensitivity on), and that returns the outputs. -""" - -# ╔═╡ 9622f7ca-f71a-4ad9-a309-d7d10a1c3e3b -# function fermenter_monod_sim(parms) -# missing -# ... -# end - -# ╔═╡ 4a5971b1-f4d0-43b6-805f-e17f5052ae92 -md""" -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(parms) = missing - -# ╔═╡ 3ae5bd00-2e06-4789-aab3-d897824d5e29 -# fermenter_monod_sim_X(params) = missing - -# ╔═╡ 6913fb1b-1986-4b37-819e-a6a7bd91f1a5 -md""" -### Compute the outputs -""" - -# ╔═╡ 4bd2bcca-9c42-4333-b062-2aaa9f7be3fe -md""" -Make the time vector. -""" - -# ╔═╡ fbd98975-aa32-46ae-8db0-0e65cdf48309 -# t_vals = missing - -# ╔═╡ 93791eb3-1eaa-4146-90b5-c4811fb3485b -md""" -Compute the two outputs $S$ and $X$ for the given parameter values. -""" - -# ╔═╡ dc0557d6-81b9-4759-8ed7-3129f60c6dc3 -# S_sim = missing - -# ╔═╡ 95bc683c-f6e6-4b42-b90b-b5a863edd4d5 -# 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 - -# ╔═╡ 49a94b9a-a543-495e-b4f1-c8579e59304d -# sens_X = missing - -# ╔═╡ 9cace6c1-e678-4dd7-8705-92a55eb32fa9 -md""" -Extract the (absolute) sensitivities of the outputs on the different parameters. -""" - -# ╔═╡ a6dc2b60-6a0a-4140-892e-02cde8dc79d3 -# begin -# sens_S_on_μmax = missing -# sens_S_on_Ks = missing -# sens_S_on_Sin = missing -# end - -# ╔═╡ f806c243-9032-46b7-add3-4714344691c7 -# begin -# sens_X_on_μmax = missing -# sens_X_on_Ks = missing -# sens_X_on_Sin = missing -# end - -# ╔═╡ f4d0aaaa-bb3f-4335-901a-dbbcde643c78 -md""" -#### Normalized sensitivities -""" - -# ╔═╡ 5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 -md""" -Compute the normalized sensitivities. -""" - -# ╔═╡ 76846731-929c-408f-a3de-970581c497e9 -# begin -# sens_S_on_μmax_rel = missing -# sens_S_on_Ks_rel = missing -# sens_S_on_Sin_rel = missing -# end - -# ╔═╡ b6c57444-547c-4e82-8526-6a30566e07c5 -# begin -# sens_X_on_μmax_rel = missing -# sens_X_on_Ks_rel = missing -# sens_X_on_Sin_rel = missing -# end - -# ╔═╡ 7f0b9e63-2856-4601-9f14-c20c6b1b3707 -md""" -### Plotting + questions -""" - -# ╔═╡ 5388c2a7-5a11-4da8-be09-46045cde8a4e -md""" -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 - -# ╔═╡ d41375ef-6958-4705-a417-4c6a491232ee -md""" -!!! 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 -""" - -# ╔═╡ 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=...`). -""" - -# ╔═╡ c0223da4-9959-48d0-b607-633b2e82986c -# missing - -# ╔═╡ ff86a29f-9308-473b-aa1c-dfd4af8179c7 -md""" -!!! 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 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 -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=...`). -""" - -# ╔═╡ 53134149-0bf7-41c1-9b35-e5037744211f -# missing - -# ╔═╡ 355ca6a7-466b-4969-ab48-28e2257f9810 -md""" -!!! 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 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. -""" - -# ╔═╡ ecc28b8e-0eaf-489f-8104-7c6383d4e01f -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 - -# ╔═╡ 62791c48-5e38-46fe-8336-3a21003f3fe3 -md""" -Now get sample 100 solutions from your monod_deviation model. -""" - -# ╔═╡ 9e049653-347a-4198-90fd-23f15c404d0c -# missing - -# ╔═╡ 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 - -# ╔═╡ 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""" -#### 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 - -# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 -# monod_cond_model = missing - -# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 -md""" -### MCMC -""" - -# ╔═╡ b8cd84de-57e4-4d8b-92da-e7c5b6d627bd -# monod_chain = missing - -# ╔═╡ 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 - -# ╔═╡ 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 -# ╠═3ef93246-657d-4e77-9bf0-8380c64bfcfd -# ╠═a355b0ba-baaf-49f4-a5dc-965364a884f0 -# ╠═00fd6d49-f561-42e9-9413-d33af92f83dc -# ╠═7ae714c4-d25d-4f9f-ab3d-cc067db9c156 -# ╠═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 -# ╠═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 -# ╟─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_fermenter_monod.jl b/exercises/student_notebooks/P7_sens_uncert/uncert_fermenter_monod.jl deleted file mode 100644 index efb0f57c..00000000 --- a/exercises/student_notebooks/P7_sens_uncert/uncert_fermenter_monod.jl +++ /dev/null @@ -1,174 +0,0 @@ -### A Pluto.jl notebook ### -# v0.20.4 - -using Markdown -using InteractiveUtils - -# ╔═╡ 9d8acf40-e635-4dc4-9938-ec63ab68e3bd -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 - -# ╔═╡ 1ec7def2-0c43-11ef-0850-ddefa680d2a4 -using Markdown - -# ╔═╡ 2861bc17-7810-461e-8857-a55b6392977a -using InteractiveUtils - -# ╔═╡ 093679cc-122b-4b0a-84c3-15cf3a6f5c1b -using Catalyst - -# ╔═╡ 0720732e-105a-40c2-bf45-f8fbbc10e5bc -using OrdinaryDiffEq, StatsPlots - -# ╔═╡ 28bb60c9-ab56-497f-b039-0d6577501be5 -using Measurements - -# ╔═╡ 437364b6-a953-428b-a891-e17ffd531656 -md""" -### Exercise: Fermenter - Monod kinetics - Uncertainty analysis -""" - -# ╔═╡ c44e10f3-8452-4a60-b840-f4a7e6ec3238 -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 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: - -$$\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} -\end{eqnarray*}$$ -""" - -# ╔═╡ 356822df-bad3-41f9-8556-de04234b0a09 -md""" -The *reaction network object* for this model could be set-up as: -""" - -# ╔═╡ c551c33e-95b0-4ea8-85d4-6c68db72b8bc -fermenter_monod = @reaction_network begin - μmax/(S+Ks), S + X --> (1 + Y)*X - Q/V, (S, X) --> 0 - Q/V*Sin, 0 --> S -end - -# ╔═╡ 69f8c265-fe91-4e4b-971b-395325d4474e -md""" -which resulted in the following differential equations: -""" - -# ╔═╡ 2806654f-7a24-4f8e-8b5a-3129b2431c95 -md""" -$$\begin{eqnarray*} -\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ -\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X -\end{eqnarray*}$$ -""" - -# ╔═╡ 5d6a881d-421e-4dea-9668-38b3f3378bf1 -# osys = missing - -# ╔═╡ 5802078d-12bc-489e-8566-15f77a11e8ba -md""" -Assume the uncertainties in the following parameter values: - -- $\mu_{max} = 0.40 \pm 0.06\;h^{-1}$ -- $K_s = 0.015 \pm 0.003 \;g/L$ -- $S_{in} = 0.022 \pm 0.004\;g/L$ - -and that the uncertainty in the other parameters values: $Y = 0.67$, $Q = 2.0\;L/h$, $V = 40.0\;L$ are negligible. 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$. Perform an uncertainty analysis by plotting the uncertainty bands on the simulation results of $S$ and $X$ in a timespan of $[0, 100]\,h$. - -Interpret your results. -""" - -# ╔═╡ 562e0006-bb79-4895-98b4-c376ea6ce856 -md""" -Initialize a vector `u0` with the initial conditions, and set the timespan: -""" - -# ╔═╡ 6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde -# u0 = missing # Uncomment and complete the instruction - -# ╔═╡ 1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b -# tspan = missing # Uncomment and complete the instruction - -# ╔═╡ 40a8c697-56d1-4488-b089-00ddf41cdb5a -md""" -We initialize a vector `params_uncert` with the parameter values and their corresponding uncertainty: -""" - -# ╔═╡ 62f4049b-617f-4268-beba-b2c634b46642 -# params_uncert = missing # Uncomment and complete the instruction - -# ╔═╡ 0e4c73d3-1630-463b-994d-073170c1a6c3 -md""" -We create the corresponding ODE problem and store it in `oprob_uncert`: -""" - -# ╔═╡ 3cad11d6-c2eb-4865-858a-d407aa25d780 -# oprob_uncert = missing # Uncomment and complete the instruction - -# ╔═╡ a7b3b68f-d34e-4e9a-aa44-224919c44fec -md""" -We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol_uncert`: -""" - -# ╔═╡ 40168902-0780-4f92-8eca-23982c80511e -# osol_uncert = missing # Uncomment and complete the instruction - -# ╔═╡ 6b2abf2a-f1f3-4c2d-96e8-cc626848fa68 -md""" -Plot the results (simulation of the output variables $S$ and $X$ together with their uncertainty band): -""" - -# ╔═╡ 4d753fdb-089d-41dc-89ef-45a1733af786 -# missing - -# ╔═╡ 4f3900d7-4ff8-4b4b-a59a-f7079ac6c457 -md""" -Try to relate the local sensitivity analysis to the uncertainty analysis. Hence, study the effect of the individual parameter uncertainties on the output variables $S$ and $X$ and compare with your local sensitivity results of the corresponding parameter. - -In order to do that, analyse the effect on the uncertainty bands for $S$ and $X$ by taking one uncertainty on a parameter at a time. In other words, analyse the subsequent cases separately: -- Assume uncertainty only in $\mu_{max}$ -- Assume uncertainty only in $K_s$ -- Assume uncertainty only in $S_{in}$ -""" - -# ╔═╡ ff3e6f75-c02b-4530-8b9a-6ff40915cb5c -md""" -!!! question - Draw your conclusions. -""" - -# ╔═╡ 33b4d1c8-5da2-4490-afe9-4778e11d3cc9 -md"- Answer: missing" - -# ╔═╡ Cell order: -# ╠═1ec7def2-0c43-11ef-0850-ddefa680d2a4 -# ╠═2861bc17-7810-461e-8857-a55b6392977a -# ╠═9d8acf40-e635-4dc4-9938-ec63ab68e3bd -# ╠═093679cc-122b-4b0a-84c3-15cf3a6f5c1b -# ╠═0720732e-105a-40c2-bf45-f8fbbc10e5bc -# ╠═28bb60c9-ab56-497f-b039-0d6577501be5 -# ╟─437364b6-a953-428b-a891-e17ffd531656 -# ╟─c44e10f3-8452-4a60-b840-f4a7e6ec3238 -# ╟─356822df-bad3-41f9-8556-de04234b0a09 -# ╠═c551c33e-95b0-4ea8-85d4-6c68db72b8bc -# ╟─69f8c265-fe91-4e4b-971b-395325d4474e -# ╟─2806654f-7a24-4f8e-8b5a-3129b2431c95 -# ╠═5d6a881d-421e-4dea-9668-38b3f3378bf1 -# ╟─5802078d-12bc-489e-8566-15f77a11e8ba -# ╟─562e0006-bb79-4895-98b4-c376ea6ce856 -# ╠═6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde -# ╠═1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b -# ╟─40a8c697-56d1-4488-b089-00ddf41cdb5a -# ╠═62f4049b-617f-4268-beba-b2c634b46642 -# ╟─0e4c73d3-1630-463b-994d-073170c1a6c3 -# ╠═3cad11d6-c2eb-4865-858a-d407aa25d780 -# ╟─a7b3b68f-d34e-4e9a-aa44-224919c44fec -# ╠═40168902-0780-4f92-8eca-23982c80511e -# ╟─6b2abf2a-f1f3-4c2d-96e8-cc626848fa68 -# ╠═4d753fdb-089d-41dc-89ef-45a1733af786 -# ╟─4f3900d7-4ff8-4b4b-a59a-f7079ac6c457 -# ╟─ff3e6f75-c02b-4530-8b9a-6ff40915cb5c -# ╠═33b4d1c8-5da2-4490-afe9-4778e11d3cc9 diff --git a/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl b/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl index b5f13f8c..b1064c14 100644 --- a/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl +++ b/exercises/student_notebooks/P7_sens_uncert/uncert_intro.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.21 +# v1.0.3 using Markdown using InteractiveUtils @@ -199,27 +199,11 @@ In uncertainty propagation, we are interested in the effect of uncertainties (or 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. @@ -233,6 +217,22 @@ Firstly, solve the original system (Wf = 10 and μ = 0.07) and afterwards we can # ╔═╡ 62afc495-29c3-4baa-b45a-94b84149d15d oprob_log = ODEProblem(growth_mod_log, u0_log,tspan, [:μ => 0.07, :Wf => 10]); +# ╔═╡ 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 + +# ╔═╡ 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); + # ╔═╡ 691a8b68-935e-4012-aa25-9cfdc323b1de osol_log = solve(oprob_log, Tsit5(), saveat = 0.5); diff --git a/notebook-checks/sync_exercises.jl b/notebook-checks/sync_exercises.jl index 973009e0..8ba1cba7 100644 --- a/notebook-checks/sync_exercises.jl +++ b/notebook-checks/sync_exercises.jl @@ -72,11 +72,10 @@ const FRONTMATTER = Dict{String,NamedTuple}( "optim_wastewater_treatment.jl" => (order = "34", title = "6. Optimisation wastewater treatment", description = "Optimisation wastewater treatment", author = "Gauthier Vanhaelewyn"), # P7_sens_uncert "sens_intro.jl" => (order = "35", title = "7. Sensitivity intro", description = "Sensitivity intro", author = "Gauthier Vanhaelewyn"), - "sens_fermenter_monod.jl" => (order = "36", title = "7. Sensitivity fermenter monod", description = "Sensitivity fermenter monod", author = "Gauthier Vanhaelewyn"), "sens_bitrophic_model.jl" => (order = "37", title = "7. Sensitivity bitrophic model", description = "Sensitivity bitrophic model", author = "Gauthier Vanhaelewyn"), "sens_insuline.jl" => (order = "38", title = "7. Sensitivity insuline", description = "Sensitivity insuline", author = "Gauthier Vanhaelewyn"), "uncert_intro.jl" => (order = "39", title = "7. Uncertainty intro", description = "Uncertainty intro", author = "Gauthier Vanhaelewyn"), - "uncert_fermenter_monod.jl" => (order = "40", title = "7. Uncertainty fermenter monod", description = "Uncertainty fermenter monod", author = "Gauthier Vanhaelewyn"), + "fermenter_mutant.jl" => (order = "40", title = "7. Fermenter mutant", description = "Sensitivity and uncertainty: switching to a mutant strain?", author = "Ward Claerhout"), "uncert_bitrophic_model.jl" => (order = "41", title = "7. Uncertainty bitrophic model", description = "Uncertainty bitrophic model", author = "Gauthier Vanhaelewyn"), # P8_modselect "model_selection_intro.jl" => (order = "42", title = "8. Model selection intro", description = "Model selection intro", author = nothing), diff --git a/src/exercises/fermenter_mutant.jl b/src/exercises/fermenter_mutant.jl new file mode 100644 index 00000000..ef3897db --- /dev/null +++ b/src/exercises/fermenter_mutant.jl @@ -0,0 +1,603 @@ +### A Pluto.jl notebook ### +# v0.20.28 + +#> [frontmatter] +#> order = "40" +#> title = "7. Fermenter mutant" +#> tags = ["exercises"] +#> layout = "layout.jlhtml" +#> description = "Sensitivity and uncertainty: switching to a mutant strain?" +#> +#> [[frontmatter.author]] +#> name = "Ward Claerhout" + +using Markdown +using InteractiveUtils + +# ╔═╡ 9d8acf40-e635-4dc4-9938-ec63ab68e3bd +# Running this yourself? Point this at your own environment — +# we advise one shared project in the parent folder: Pkg.activate("..") +using Pkg; Pkg.activate("../../pluto-deployment-environment") + +# ╔═╡ 1ec7def2-0c43-11ef-0850-ddefa680d2a4 +using Markdown + +# ╔═╡ 2861bc17-7810-461e-8857-a55b6392977a +using InteractiveUtils + +# ╔═╡ 093679cc-122b-4b0a-84c3-15cf3a6f5c1b +using Catalyst + +# ╔═╡ 0720732e-105a-40c2-bf45-f8fbbc10e5bc +using OrdinaryDiffEq + +# ╔═╡ 28bb60c9-ab56-497f-b039-0d6577501be5 +using StatsPlots, PlutoUI; TableOfContents() + +# ╔═╡ 0a65afbc-643a-4fe4-b22d-c5fe54e3be1f +using Turing + +# ╔═╡ d1b6739f-1e0c-452d-8763-717dd650b2ec +using ForwardDiff + +# ╔═╡ 162b330b-02f7-4179-b950-d960eca68206 +using Random + +# ╔═╡ 437364b6-a953-428b-a891-e17ffd531656 +md""" +# Exercise: Fermenter - Switching to a mutant strain? +""" + +# ╔═╡ c44e10f3-8452-4a60-b840-f4a7e6ec3238 +md""" +As a process engineer, you are responsible for a continuously stirred 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$ [$L/h$] containing substrate at a concentration $S_{in}$ [$g/L$]. The yield $Y$ [$g$ biomass$/g$ substrate] tells how much biomass is formed per gram of substrate consumed. As in the previous practica, the growth is modelled with Monod kinetics: + +$$\begin{eqnarray*} +S \xrightarrow[\quad\quad]{r} Y \, X \quad\quad\quad\quad r = \mu \, X \quad \textrm{with} \quad \mu = \mu_{max} \, \cfrac{S}{S + K_s} +\end{eqnarray*}$$ +""" + +# ╔═╡ f3065171-5526-476a-bfd7-a2b95bb24d61 +md""" +## The question + +A colleague has found a mutant strain with a **higher yield** $Y$, and management asks you a simple question: *should we switch to the mutant?* + +The catch: the maximal growth rate $\mu_{max}$ of the mutant is not known. A higher yield often comes at the cost of slower growth, and a slower grower leaves more substrate unused in the outflow. So the mutant could be better or worse than the current strain. + +You will answer this question with the tools of this practical: +1. **Uncertainty propagation**: what does our uncertainty about $\mu_{max}$ mean for the biomass? How sure are we that the mutant is better? +2. **Sensitivity analysis**: when does the biomass depend most strongly on $\mu_{max}$? That is the moment where a measurement tells us the most. +3. **Bayesian inference**: update our belief about $\mu_{max}$ with a few measurements, and propagate the reduced uncertainty again. +""" + +# ╔═╡ 356822df-bad3-41f9-8556-de04234b0a09 +md""" +The *reaction network object* for this model is: +""" + +# ╔═╡ c551c33e-95b0-4ea8-85d4-6c68db72b8bc +fermenter_monod = @reaction_network begin + X * mm(S, μmax, Ks), S => Y*X + Q/V, (S, X) --> 0 + Q/V*Sin, ∅ --> S +end + +# ╔═╡ 69f8c265-fe91-4e4b-971b-395325d4474e +md""" +which results in the following differential equations: +""" + +# ╔═╡ 2806654f-7a24-4f8e-8b5a-3129b2431c95 +md""" +$$\begin{eqnarray*} +\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ +\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X +\end{eqnarray*}$$ +""" + +# ╔═╡ 5d6a881d-421e-4dea-9668-38b3f3378bf1 +# osys = missing + +# ╔═╡ 5802078d-12bc-489e-8566-15f77a11e8ba +md""" +## The current strain + +For the strain currently in use, the parameters are well known: +- $\mu_{max} = 0.40\;h^{-1}$ +- $K_s = 0.015 \;g/L$ +- $Y = 0.67$ +- $S_{in} = 0.022 \;g/L$ +- $Q = 2.0\;L/h$ +- $V = 40.0\;L$ + +At $t=0$ no substrate is present in the reactor and the biomass concentration is $0.0005\;g/L$. We judge a strain by the biomass concentration in the outflow after $100\,h$. + +Initialize a vector `u0` with the initial conditions, set the timespan `tspan` to $[0, 100]\,h$, the time step `dt` to `0.5`, and initialize a vector `parms` with the parameter values. +""" + +# ╔═╡ 6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde +# u0 = missing + +# ╔═╡ 1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b +# tspan = missing + +# ╔═╡ 9e3c427b-8a56-4787-9b28-43ec3932ddbe +# dt = missing + +# ╔═╡ 62f4049b-617f-4268-beba-b2c634b46642 +# parms = missing + +# ╔═╡ 0e4c73d3-1630-463b-994d-073170c1a6c3 +md""" +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. +""" + +# ╔═╡ 3cad11d6-c2eb-4865-858a-d407aa25d780 +# oprob = missing + +# ╔═╡ 40168902-0780-4f92-8eca-23982c80511e +# osol = missing + +# ╔═╡ 4d753fdb-089d-41dc-89ef-45a1733af786 +# missing + +# ╔═╡ ca660e2c-9348-40b1-b92d-415068bea514 +md""" +Store the biomass concentration at the end of the simulation in `X_current`. This is the value the mutant has to beat. +""" + +# ╔═╡ ec042e61-04b8-4d51-8e7d-35e1646cc5ba +# X_current = missing + +# ╔═╡ 9ee22028-7a37-4510-9b61-198d3bd99a73 +md""" +## The mutant strain +""" + +# ╔═╡ cc9ac43b-448a-4a67-ba0b-8fdf21a4236e +md""" +Some batch tests show that the mutant has a yield of $Y = 0.70$ instead of $0.67$. All other parameters are the same as for the current strain, except $\mu_{max}$, which is unknown. + +From the literature on mutants with an improved yield, you expect $\mu_{max}$ to be around $0.35\;h^{-1}$, give or take $15\%$. You express this **prior belief** with a log-normal distribution `LogNormal(log(0.35), 0.15)`. +""" + +# ╔═╡ f5482531-b23f-472b-968f-20f7a6471778 +begin + plot(LogNormal(log(0.35), 0.15), xlabel="μmax (1/h)", ylabel="density", label="Prior μmax mutant", linewidth=2) + vline!([0.40], label="μmax current strain", linestyle=:dash) +end + +# ╔═╡ ff21b3cd-64b1-4d09-8eb1-afc617585fc8 +md""" +Create an ODE problem `oprob_mut` for the mutant by remaking `oprob` with $Y = 0.70$ and the expected value $\mu_{max} = 0.35$. + +!!! note + Remember that you can remake an ODEProblem with new parameter values by using `remake(oprob, p = [...])`. +""" + +# ╔═╡ 225eaaff-42a2-4aca-94e1-afcd988de4eb +# oprob_mut = missing + +# ╔═╡ 18ffe877-c579-4b4c-80bd-b5d50baf36bb +md""" +## Monte Carlo error propagation + +We don't know $\mu_{max}$, but we do know which values are plausible. Using Monte Carlo simulation, we draw many values from the prior, simulate the mutant for each value, and look at the spread of the outputs. This shows how the uncertainty on the parameter **propagates** to the uncertainty on the output. + +Make a Turing model `mutant_deviation` in which you implement the prior for `μmax_dev` and return the solution of the remade ODE problem. +""" + +# ╔═╡ 7821dc21-3e6b-4ad1-829e-f29dcd020144 +# @model function mutant_deviation() +# μmax_dev ~ missing +# new_prob = missing +# return missing +# end + +# ╔═╡ 3051d39a-fc2f-40e4-a688-4d6e75cfc2d1 +md""" +Now sample 500 solutions from your `mutant_deviation` model. +""" + +# ╔═╡ 33df05c9-08d3-4107-9240-931658cda2a4 +# missing + +# ╔═╡ 02c5cf65-ebd3-477b-9a3e-8855bd46f55c +md""" +Plot the biomass of all Monte Carlo simulations together with the biomass of the current strain. Use `idxs=:X`, `label=false`, `color=:gray`, `alpha=0.4` and `linestyle=:dot`. +""" + +# ╔═╡ aa833d92-567f-4971-a360-e60bce424c8e +# begin +# p1 = plot(title="Monte Carlo of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") +# for missing in missing +# missing +# end +# plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) +# p1 +# end + +# ╔═╡ ffc66a94-5435-4ca6-af87-76c85a599906 +md""" +Create a vector `x1` that stores the biomass at the end of the simulations, make a histogram of it, and compute the probability `P_better_prior` that the mutant produces more biomass than the current strain. + +!!! note + The fraction of elements of a vector `x` that are larger than `a` can be computed with `mean(x .> a)`. +""" + +# ╔═╡ d3afd5b4-034f-415c-a64a-776deed74d5c +# x1 = missing + +# ╔═╡ 1d95728e-079e-4e9b-952a-c4f692661813 +# missing + +# ╔═╡ 02213469-a233-4759-a18d-29a76a4c13c6 +# P_better_prior = missing + +# ╔═╡ 97d2065c-6234-4b13-8e4c-6782b8301416 +md""" +!!! questions + - When is the spread on $X$ largest, and when is it smallest? Can you explain this? + - How likely is it that the mutant beats the current strain? Would you advise management to switch based on this? +""" + +# ╔═╡ d2acbfb9-a7ff-48b0-9d9a-9492ecb9c1bf +md""" +Answers: +- missing +- missing +""" + +# ╔═╡ eb613ead-6aea-41fe-acb7-f679b5a97c35 +md""" +## Local Sensitivity Analysis (LSA) + +The lab can run the mutant in a pilot reactor under the same conditions, but you can only afford **three measurements** of $X$, all at the same time point. Measurements are expensive, so we want them at the moment where $X$ depends most strongly on $\mu_{max}$. There, a small change in $\mu_{max}$ causes a large change in $X$. The other way around, a measurement of $X$ at that moment pins down $\mu_{max}$ best. + +So, we will first perform a normalized local sensitivity analysis to pin down the moment, $X$ is most sensitive to a change in $μmax$. +""" + +# ╔═╡ 67f572b0-3485-40bc-8f64-cc58471b38a9 +md""" +Write a solution function `mutant_sim` with as argument a vector with the parameter `μmax` of the mutant, and that returns the solution. Based on it, write a function `mutant_sim_X` that only returns the output $X$. +""" + +# ╔═╡ 565b6e18-012f-4b38-bf21-bdc683688107 +# function mutant_sim(parms) +# missing +# ... +# end + +# ╔═╡ b7f14c48-ba19-4fbc-bded-829a34bd7ac8 +# mutant_sim_X(parms) = missing + +# ╔═╡ 73345fe4-6877-476d-8ce1-5a0056c73cc6 +md""" +Make the time vector `t_vals` with timestep `dt` equal to $0.5$ and compute the output `X_sim` for the expected value $\mu_{max} = 0.35$, stored in `μmax_val`. +""" + +# ╔═╡ 75ab68b2-1e88-421f-a089-4c0b3b72788d +# t_vals = missing + +# ╔═╡ 82ae802b-2a7a-4525-9cc6-dfaa552ba990 +# μmax_val = missing + +# ╔═╡ 8a655d11-0d45-4e7c-bcbf-8ea9659516aa +# X_sim = missing + +# ╔═╡ 9a97a7c7-bc3e-434e-9dce-7c859417aac2 +md""" +Use `ForwardDiff.jacobian` to compute the sensitivity `sens_X_on_μmax`, and compute the normalized sensitivity `sens_X_on_μmax_rel`. +""" + +# ╔═╡ d7d3028a-51c1-41e0-a865-b521ed235548 +# sens_X_on_μmax = missing + +# ╔═╡ cb6047f2-3884-4fbd-96f2-f15cf684bea8 +# sens_X_on_μmax_rel = missing + +# ╔═╡ f623b4b4-6010-45fa-b9a9-bc796f0fb106 +md""" +Identify the time `t_star` at which the biomass is most sensitive to variations in $\mu_{max}$, and plot the normalized sensitivity function together with `t_star`. + +!!! 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. +""" + +# ╔═╡ 9c3fc6d8-8e5e-425d-9caa-56a9e0b03b99 +# begin +# t_star_idx = missing +# t_star = t_vals[t_star_idx] +# +# plot(t_vals, abs.(sens_X_on_μmax_rel), +# title="Sensitivity of X w.r.t. μmax", xlabel="Time (hours)", +# ylabel="normalized sensitivity", label=false, linewidth=2) +# vline!([t_star], label="t* = $(t_star) h", linestyle=:dash, color=:red) +# end + +# ╔═╡ b1825a21-e711-4f99-9a4b-51f3fb73ff82 +md""" +!!! questions + - Compare the sensitivity function with the Monte Carlo simulations above. What do you notice? + - The sensitivity at $t = 100\,h$ is small, while that is exactly the time we care about. Why is it small? Is it a good idea to measure at $t = 100\,h$? +""" + +# ╔═╡ 44234af3-6015-49b6-8709-1f3a3cc0f69d +md""" +Answers: +- missing +- missing +""" + +# ╔═╡ a00afabe-7bc1-4d8b-9b55-2953c7add9f5 +md""" +## Measuring campaign + +The lab ran the mutant in the pilot reactor and sampled it at your chosen time. We have defined a function `measure_X(t)` in this notebook that returns three hypothetical measurements of $X$ at the given timepoint `t`, with a relative measurement error of $5\%$. Use this function to get the measurements `X_meas` at the optimal time `t_star`. +""" + +# ╔═╡ 635ff400-639c-4c80-93e6-e8721a836c16 +begin + μmax_real = 0.37 + real_solution = solve(remake(oprob_mut, p = [:μmax => μmax_real]), Tsit5(), saveat=dt) + measure_X(t) = real_solution(t, idxs=:X) .* (1 .+ 0.05 .* randn(Xoshiro(round(Int, 10t)), 3)) +end; + +# ╔═╡ 93cc5a4f-f40a-4d28-abc3-8d6790fd60dc +# X_meas = missing + +# ╔═╡ 3e93a001-9401-4a1f-b283-c79fd0c43434 +md""" +### Conditioning on the measurements +""" + +# ╔═╡ 04ae6cd7-2d59-4bc0-a894-d9d5129e243d +md""" +Based on these measurements we can update what we believed before (the prior) with what the data tell us (the likelihood) into an updated belief (the posterior): + +Create a Turing model `mutant_meas(t_meas)` with: +- the same prior for `μmax_dev` as before; +- a prior `Exponential()` for the measurement error `σ_X`; +- the simulated biomass `X_pred` at time `t_meas`, which you can get with `sol(t_meas, idxs=:X)`; +- the three measurements `X_meas`, each normally distributed around `X_pred` with standard deviation `σ_X`; +- the solution as return value. +""" + +# ╔═╡ 02096fef-224e-4896-9f8d-3f7aba7334e1 +# @model function mutant_meas(t_meas) +# μmax_dev ~ missing +# σ_X ~ missing +# +# sol = missing +# X_pred = missing +# +# X_meas = zeros(3) +# for i in 1:length(X_meas) +# X_meas[i] ~ missing +# end +# return sol +# end + +# ╔═╡ 354f0768-efc0-4304-aa04-c6a68f961aa6 +md""" +Condition the model on the measurements and store it in `mutant_cond_model`. +""" + +# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 +# mutant_cond_model = missing + +# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 +md""" +### MCMC +""" + +# ╔═╡ f99ca771-8c14-4e19-a1ae-478a2103716a +# mutant_chain = missing + +# ╔═╡ e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +md""" +Plot the chain to check if it properly converged to the posterior distribution. +""" + +# ╔═╡ 50b30e65-72d3-4d8d-a739-ae969e0f50a8 +# missing + +# ╔═╡ a65710f6-70ce-4348-923d-ed2ed2b4a469 +md""" +Now, get the sampled posterior solutions `solutions_mutant` and plot them on top of the prior Monte Carlo simulations, together with the current strain and the measurements. +""" + +# ╔═╡ 12c9f810-ef5b-49a4-8fdd-6efb1a040440 +# solutions_mutant = missing + +# ╔═╡ e2e11092-4498-437b-8c67-8d77302893fe +md""" +The posterior chance that the mutant is better is: +""" + +# ╔═╡ cae4d9c4-3790-423e-8321-d44c0833aa40 +begin + sol_end = [solutions_mutant[i][:X][end] for i in eachindex(solutions_mutant)] + P_better_posterior = mean(sol_end .> X_current) + println("The chance the mutant is better is $(P_better_posterior*100)%") +end + +# ╔═╡ 49307c36-3518-4118-b453-91e634f92daf +# begin +# p2 = plot(title="Prior vs posterior of the mutant", ylabel="X (g/l)", xlabel="Time (hours)") +# for missing in missing +# missing +# end +# for missing in missing +# missing +# end +# plot!(osol, idxs=:X, label="Current strain", color=:red, linewidth=3) +# scatter!(fill(t_star, 3), X_meas, label="Measurements", color=:black) +# p2 +# end + +# ╔═╡ aa479bcd-2314-47a5-83ba-b0d1a2654fa0 +md""" +Create a vector `x2` with the posterior biomass at the end of the simulations, compare its histogram with the prior one, and compute the probability `P_better_post` that the mutant beats the current strain. +""" + +# ╔═╡ 2f55a771-0d37-48e5-a9d9-24693bf1d177 +# x2 = missing + +# ╔═╡ 18059d0c-cf37-44a1-86a8-859ac66a66f0 +# missing + +# ╔═╡ 40e73ddf-6c0d-4f83-860f-434d934193ed +# P_better_post = missing + +# ╔═╡ e1c87260-25ee-4d71-8744-ed0d7cb75f6a +md""" +### What if we had measured at the end? + +Intuitively, you might have measured $X$ at $t = 100\,h$, since that is the quantity we care about. Repeat the inference with three measurements at $t = 100\,h$, and compute the probability `P_better_100` that the mutant is better. +""" + +# ╔═╡ caa5c04d-49af-4576-83a2-605bae1c2494 +# X_meas_100 = missing + +# ╔═╡ 6a53aff6-29d4-4d41-8236-f621ce519361 +# missing + +# ╔═╡ 9292c593-1f0f-47a6-a2ec-91a8aadbe6e6 +# P_better_100 = missing + +# ╔═╡ b1814ea6-0065-4ac2-a012-505d2bbd8e40 +md""" +Compare what we know about $\mu_{max}$ before and after both measuring campaigns: +""" + +# ╔═╡ 44cb7615-86b2-45f0-b87a-69fcbfe863bb +begin + density(vec(chain[:μmax_dev]), label="Prior", linewidth=2, title="What did we learn about μmax?", xlabel="μmax (1/h)", ylabel="density") + density!(vec(mutant_chain_100[:μmax_dev]), label="Posterior, measured at t = 100 h", linewidth=2) + density!(vec(mutant_chain[:μmax_dev]), label="Posterior, measured at t* = $(t_star) h", linewidth=2) +end + +# ╔═╡ f364fa1f-5c47-4ca6-99df-17e2b735d98b +md""" +!!! questions + - Compare `P_better_prior`, `P_better_100` and `P_better_post`. +""" + +# ╔═╡ 5f9873a5-af70-4ac2-be08-4f249b71eb76 +md""" +Answers: +- missing +- missing +""" + +# ╔═╡ 3980da6b-7036-44d2-9d8f-3ddc1062d0ad +md""" +## Conclusions + +!!! questions + - What do you advise management? Support your advice with numbers. + - What did you learn/what did you use the information for from uncertainty propagation, sensitivity analysis and the Bayesian inference. +""" + +# ╔═╡ 8021123d-1b19-47bf-9db0-0d7c587bc00e +md""" +Answers: +- missing +- missing +- missing +""" + +# ╔═╡ Cell order: +# ╠═1ec7def2-0c43-11ef-0850-ddefa680d2a4 +# ╠═2861bc17-7810-461e-8857-a55b6392977a +# ╠═9d8acf40-e635-4dc4-9938-ec63ab68e3bd +# ╠═093679cc-122b-4b0a-84c3-15cf3a6f5c1b +# ╠═0720732e-105a-40c2-bf45-f8fbbc10e5bc +# ╠═28bb60c9-ab56-497f-b039-0d6577501be5 +# ╠═0a65afbc-643a-4fe4-b22d-c5fe54e3be1f +# ╠═d1b6739f-1e0c-452d-8763-717dd650b2ec +# ╠═162b330b-02f7-4179-b950-d960eca68206 +# ╟─437364b6-a953-428b-a891-e17ffd531656 +# ╟─c44e10f3-8452-4a60-b840-f4a7e6ec3238 +# ╟─f3065171-5526-476a-bfd7-a2b95bb24d61 +# ╟─356822df-bad3-41f9-8556-de04234b0a09 +# ╠═c551c33e-95b0-4ea8-85d4-6c68db72b8bc +# ╟─69f8c265-fe91-4e4b-971b-395325d4474e +# ╟─2806654f-7a24-4f8e-8b5a-3129b2431c95 +# ╠═5d6a881d-421e-4dea-9668-38b3f3378bf1 +# ╟─5802078d-12bc-489e-8566-15f77a11e8ba +# ╠═6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde +# ╠═1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b +# ╠═9e3c427b-8a56-4787-9b28-43ec3932ddbe +# ╠═62f4049b-617f-4268-beba-b2c634b46642 +# ╟─0e4c73d3-1630-463b-994d-073170c1a6c3 +# ╠═3cad11d6-c2eb-4865-858a-d407aa25d780 +# ╠═40168902-0780-4f92-8eca-23982c80511e +# ╠═4d753fdb-089d-41dc-89ef-45a1733af786 +# ╟─ca660e2c-9348-40b1-b92d-415068bea514 +# ╠═ec042e61-04b8-4d51-8e7d-35e1646cc5ba +# ╟─9ee22028-7a37-4510-9b61-198d3bd99a73 +# ╟─cc9ac43b-448a-4a67-ba0b-8fdf21a4236e +# ╟─f5482531-b23f-472b-968f-20f7a6471778 +# ╟─ff21b3cd-64b1-4d09-8eb1-afc617585fc8 +# ╠═225eaaff-42a2-4aca-94e1-afcd988de4eb +# ╟─18ffe877-c579-4b4c-80bd-b5d50baf36bb +# ╠═7821dc21-3e6b-4ad1-829e-f29dcd020144 +# ╟─3051d39a-fc2f-40e4-a688-4d6e75cfc2d1 +# ╠═33df05c9-08d3-4107-9240-931658cda2a4 +# ╟─02c5cf65-ebd3-477b-9a3e-8855bd46f55c +# ╠═aa833d92-567f-4971-a360-e60bce424c8e +# ╟─ffc66a94-5435-4ca6-af87-76c85a599906 +# ╠═d3afd5b4-034f-415c-a64a-776deed74d5c +# ╠═1d95728e-079e-4e9b-952a-c4f692661813 +# ╠═02213469-a233-4759-a18d-29a76a4c13c6 +# ╟─97d2065c-6234-4b13-8e4c-6782b8301416 +# ╟─d2acbfb9-a7ff-48b0-9d9a-9492ecb9c1bf +# ╟─eb613ead-6aea-41fe-acb7-f679b5a97c35 +# ╟─67f572b0-3485-40bc-8f64-cc58471b38a9 +# ╠═565b6e18-012f-4b38-bf21-bdc683688107 +# ╠═b7f14c48-ba19-4fbc-bded-829a34bd7ac8 +# ╟─73345fe4-6877-476d-8ce1-5a0056c73cc6 +# ╠═75ab68b2-1e88-421f-a089-4c0b3b72788d +# ╠═82ae802b-2a7a-4525-9cc6-dfaa552ba990 +# ╠═8a655d11-0d45-4e7c-bcbf-8ea9659516aa +# ╟─9a97a7c7-bc3e-434e-9dce-7c859417aac2 +# ╠═d7d3028a-51c1-41e0-a865-b521ed235548 +# ╠═cb6047f2-3884-4fbd-96f2-f15cf684bea8 +# ╟─f623b4b4-6010-45fa-b9a9-bc796f0fb106 +# ╠═9c3fc6d8-8e5e-425d-9caa-56a9e0b03b99 +# ╟─b1825a21-e711-4f99-9a4b-51f3fb73ff82 +# ╟─44234af3-6015-49b6-8709-1f3a3cc0f69d +# ╟─a00afabe-7bc1-4d8b-9b55-2953c7add9f5 +# ╟─635ff400-639c-4c80-93e6-e8721a836c16 +# ╠═93cc5a4f-f40a-4d28-abc3-8d6790fd60dc +# ╟─3e93a001-9401-4a1f-b283-c79fd0c43434 +# ╟─04ae6cd7-2d59-4bc0-a894-d9d5129e243d +# ╠═02096fef-224e-4896-9f8d-3f7aba7334e1 +# ╟─354f0768-efc0-4304-aa04-c6a68f961aa6 +# ╠═11e105de-270c-434a-bfb4-9118dfc9e461 +# ╟─e46afb98-16de-4d19-86b4-b920b6b929c5 +# ╠═f99ca771-8c14-4e19-a1ae-478a2103716a +# ╟─e76a7bf7-4d8d-4f2f-8272-b63e49cd4ae9 +# ╠═50b30e65-72d3-4d8d-a739-ae969e0f50a8 +# ╟─a65710f6-70ce-4348-923d-ed2ed2b4a469 +# ╠═12c9f810-ef5b-49a4-8fdd-6efb1a040440 +# ╟─e2e11092-4498-437b-8c67-8d77302893fe +# ╠═cae4d9c4-3790-423e-8321-d44c0833aa40 +# ╠═49307c36-3518-4118-b453-91e634f92daf +# ╟─aa479bcd-2314-47a5-83ba-b0d1a2654fa0 +# ╠═2f55a771-0d37-48e5-a9d9-24693bf1d177 +# ╠═18059d0c-cf37-44a1-86a8-859ac66a66f0 +# ╠═40e73ddf-6c0d-4f83-860f-434d934193ed +# ╟─e1c87260-25ee-4d71-8744-ed0d7cb75f6a +# ╠═caa5c04d-49af-4576-83a2-605bae1c2494 +# ╠═6a53aff6-29d4-4d41-8236-f621ce519361 +# ╠═9292c593-1f0f-47a6-a2ec-91a8aadbe6e6 +# ╟─b1814ea6-0065-4ac2-a012-505d2bbd8e40 +# ╟─44cb7615-86b2-45f0-b87a-69fcbfe863bb +# ╟─f364fa1f-5c47-4ca6-99df-17e2b735d98b +# ╟─5f9873a5-af70-4ac2-be08-4f249b71eb76 +# ╟─3980da6b-7036-44d2-9d8f-3ddc1062d0ad +# ╟─8021123d-1b19-47bf-9db0-0d7c587bc00e diff --git a/src/exercises/sens_fermenter_monod.jl b/src/exercises/sens_fermenter_monod.jl deleted file mode 100644 index 55fb2d38..00000000 --- a/src/exercises/sens_fermenter_monod.jl +++ /dev/null @@ -1,678 +0,0 @@ -### A Pluto.jl notebook ### -# v0.20.21 - -#> [frontmatter] -#> order = "36" -#> title = "7. Sensitivity fermenter monod" -#> tags = ["exercises"] -#> layout = "layout.jlhtml" -#> description = "Sensitivity fermenter monod" -#> -#> [[frontmatter.author]] -#> name = "Gauthier Vanhaelewyn" - -using Markdown -using InteractiveUtils - -# ╔═╡ 3ef93246-657d-4e77-9bf0-8380c64bfcfd -# Running this yourself? Point this at your own environment — -# we advise one shared project in the parent folder: Pkg.activate("..") -using Pkg; Pkg.activate("../../pluto-deployment-environment") - -# ╔═╡ 55cdebd2-0881-11ef-2722-91de1447877a -using Markdown, InteractiveUtils - -# ╔═╡ a355b0ba-baaf-49f4-a5dc-965364a884f0 -using Catalyst, OrdinaryDiffEq - -# ╔═╡ 00fd6d49-f561-42e9-9413-d33af92f83dc -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 -""" - -# ╔═╡ 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 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} \, . -\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: -""" - -# ╔═╡ 935ca610-7a7a-4692-8908-fc26abb880b4 -fermenter_monod = @reaction_network begin - @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 -md""" -which resulted in the following differential equations: -""" - -# ╔═╡ fa93e2c3-8b43-418e-ba24-406645b2e397 -md""" -$$\begin{eqnarray*} -\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ -\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X -\end{eqnarray*}$$ -""" - -# ╔═╡ 06730f54-7293-43f2-b772-84eec3e5528a -md""" -Convert the system to a symbolic differential equation model and verify, by analyzing the differential equation, that your model is correctly implemented. - -Keep in mind that `mm(S, μmax, Ks)` stands for $\mu_{max} \, \cfrac{S}{S + K_s}$. -""" - -# ╔═╡ 7f8b7a2e-bc65-4b51-ad59-bd7ac98604dd -# osys = missing - -# ╔═╡ 1ee7359c-2757-44cd-8650-09b762cf30d0 -md""" -## Goals of this exercise -""" - -# ╔═╡ 55f1d688-0c53-481b-9965-5e92ca87ad83 -md""" -Compute the following in a timespan of `(0.0, 100.0)`$\,h$: - -- 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$ 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 `parms` with the parameter values: -""" - -# ╔═╡ 2ee277e5-ce4a-4ade-be0e-9bba7a4dc08c -# u0 = missing # in principle not necessary since we will use the default - -# ╔═╡ 3fdc6b17-cdeb-4dc5-8886-9d3a62caac8d -# tspan = missing - -# ╔═╡ 201dfb54-2056-4846-a64c-ff4951dc084d -md""" -For practical reasons, we will define the time step size as `dt = 0.5`. -""" - -# ╔═╡ 5627068f-6442-43bb-bacb-bcba917372f3 -# dt = missing - -# ╔═╡ 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. -""" - -# ╔═╡ 0f995929-4d2b-4a7a-8da1-04e4d501385f -begin - # μmax_val = missing - # Ks_val = missing - # Sin_val = missing -end; - -# ╔═╡ 79b0eb65-5a0f-40b3-aa97-4088421c562e -# 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=dt`, and store the solution in `osol`. Finally plot the results. -""" - -# ╔═╡ 8b2f23f6-80b2-4e63-942e-e5cd17d8ba72 -# oprob = missing - -# ╔═╡ 89a31c32-88a4-479f-a688-ffcb75ee8e91 -# osol = missing - -# ╔═╡ 51a9b7e6-8ad9-477d-9596-ffd614df2c79 -# 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 (that you want the sensitivity on), and that returns the outputs. -""" - -# ╔═╡ 9622f7ca-f71a-4ad9-a309-d7d10a1c3e3b -# function fermenter_monod_sim(parms) -# missing -# ... -# end - -# ╔═╡ 4a5971b1-f4d0-43b6-805f-e17f5052ae92 -md""" -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(parms) = missing - -# ╔═╡ 3ae5bd00-2e06-4789-aab3-d897824d5e29 -# fermenter_monod_sim_X(params) = missing - -# ╔═╡ 6913fb1b-1986-4b37-819e-a6a7bd91f1a5 -md""" -### Compute the outputs -""" - -# ╔═╡ 4bd2bcca-9c42-4333-b062-2aaa9f7be3fe -md""" -Make the time vector. -""" - -# ╔═╡ fbd98975-aa32-46ae-8db0-0e65cdf48309 -# t_vals = missing - -# ╔═╡ 93791eb3-1eaa-4146-90b5-c4811fb3485b -md""" -Compute the two outputs $S$ and $X$ for the given parameter values. -""" - -# ╔═╡ dc0557d6-81b9-4759-8ed7-3129f60c6dc3 -# S_sim = missing - -# ╔═╡ 95bc683c-f6e6-4b42-b90b-b5a863edd4d5 -# 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 - -# ╔═╡ 49a94b9a-a543-495e-b4f1-c8579e59304d -# sens_X = missing - -# ╔═╡ 9cace6c1-e678-4dd7-8705-92a55eb32fa9 -md""" -Extract the (absolute) sensitivities of the outputs on the different parameters. -""" - -# ╔═╡ a6dc2b60-6a0a-4140-892e-02cde8dc79d3 -# begin -# sens_S_on_μmax = missing -# sens_S_on_Ks = missing -# sens_S_on_Sin = missing -# end - -# ╔═╡ f806c243-9032-46b7-add3-4714344691c7 -# begin -# sens_X_on_μmax = missing -# sens_X_on_Ks = missing -# sens_X_on_Sin = missing -# end - -# ╔═╡ f4d0aaaa-bb3f-4335-901a-dbbcde643c78 -md""" -#### Normalized sensitivities -""" - -# ╔═╡ 5bf3a62d-d2aa-4653-8ee8-e90caa9504e8 -md""" -Compute the normalized sensitivities. -""" - -# ╔═╡ 76846731-929c-408f-a3de-970581c497e9 -# begin -# sens_S_on_μmax_rel = missing -# sens_S_on_Ks_rel = missing -# sens_S_on_Sin_rel = missing -# end - -# ╔═╡ b6c57444-547c-4e82-8526-6a30566e07c5 -# begin -# sens_X_on_μmax_rel = missing -# sens_X_on_Ks_rel = missing -# sens_X_on_Sin_rel = missing -# end - -# ╔═╡ 7f0b9e63-2856-4601-9f14-c20c6b1b3707 -md""" -### Plotting + questions -""" - -# ╔═╡ 5388c2a7-5a11-4da8-be09-46045cde8a4e -md""" -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 - -# ╔═╡ d41375ef-6958-4705-a417-4c6a491232ee -md""" -!!! 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 -""" - -# ╔═╡ 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=...`). -""" - -# ╔═╡ c0223da4-9959-48d0-b607-633b2e82986c -# missing - -# ╔═╡ ff86a29f-9308-473b-aa1c-dfd4af8179c7 -md""" -!!! 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 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 -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=...`). -""" - -# ╔═╡ 53134149-0bf7-41c1-9b35-e5037744211f -# missing - -# ╔═╡ 355ca6a7-466b-4969-ab48-28e2257f9810 -md""" -!!! 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 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. -""" - -# ╔═╡ ecc28b8e-0eaf-489f-8104-7c6383d4e01f -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 - -# ╔═╡ 62791c48-5e38-46fe-8336-3a21003f3fe3 -md""" -Now get sample 100 solutions from your monod_deviation model. -""" - -# ╔═╡ 9e049653-347a-4198-90fd-23f15c404d0c -# missing - -# ╔═╡ 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 - -# ╔═╡ 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""" -#### 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 - -# ╔═╡ 11e105de-270c-434a-bfb4-9118dfc9e461 -# monod_cond_model = missing - -# ╔═╡ e46afb98-16de-4d19-86b4-b920b6b929c5 -md""" -### MCMC -""" - -# ╔═╡ b8cd84de-57e4-4d8b-92da-e7c5b6d627bd -# monod_chain = missing - -# ╔═╡ 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 - -# ╔═╡ 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 -# ╠═3ef93246-657d-4e77-9bf0-8380c64bfcfd -# ╠═a355b0ba-baaf-49f4-a5dc-965364a884f0 -# ╠═00fd6d49-f561-42e9-9413-d33af92f83dc -# ╠═7ae714c4-d25d-4f9f-ab3d-cc067db9c156 -# ╠═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 -# ╠═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 -# ╟─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/src/exercises/uncert_fermenter_monod.jl b/src/exercises/uncert_fermenter_monod.jl deleted file mode 100644 index a18375a1..00000000 --- a/src/exercises/uncert_fermenter_monod.jl +++ /dev/null @@ -1,182 +0,0 @@ -### A Pluto.jl notebook ### -# v0.20.4 - -#> [frontmatter] -#> order = "40" -#> title = "7. Uncertainty fermenter monod" -#> tags = ["exercises"] -#> layout = "layout.jlhtml" -#> description = "Uncertainty fermenter monod" -#> -#> [[frontmatter.author]] -#> name = "Gauthier Vanhaelewyn" - -using Markdown -using InteractiveUtils - -# ╔═╡ 9d8acf40-e635-4dc4-9938-ec63ab68e3bd -# Running this yourself? Point this at your own environment — -# we advise one shared project in the parent folder: Pkg.activate("..") -using Pkg; Pkg.activate("../../pluto-deployment-environment") - -# ╔═╡ 1ec7def2-0c43-11ef-0850-ddefa680d2a4 -using Markdown - -# ╔═╡ 2861bc17-7810-461e-8857-a55b6392977a -using InteractiveUtils - -# ╔═╡ 093679cc-122b-4b0a-84c3-15cf3a6f5c1b -using Catalyst - -# ╔═╡ 0720732e-105a-40c2-bf45-f8fbbc10e5bc -using OrdinaryDiffEq, StatsPlots - -# ╔═╡ 28bb60c9-ab56-497f-b039-0d6577501be5 -using Measurements - -# ╔═╡ 437364b6-a953-428b-a891-e17ffd531656 -md""" -### Exercise: Fermenter - Monod kinetics - Uncertainty analysis -""" - -# ╔═╡ c44e10f3-8452-4a60-b840-f4a7e6ec3238 -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 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: - -$$\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} -\end{eqnarray*}$$ -""" - -# ╔═╡ 356822df-bad3-41f9-8556-de04234b0a09 -md""" -The *reaction network object* for this model could be set-up as: -""" - -# ╔═╡ c551c33e-95b0-4ea8-85d4-6c68db72b8bc -fermenter_monod = @reaction_network begin - μmax/(S+Ks), S + X --> (1 + Y)*X - Q/V, (S, X) --> 0 - Q/V*Sin, 0 --> S -end - -# ╔═╡ 69f8c265-fe91-4e4b-971b-395325d4474e -md""" -which resulted in the following differential equations: -""" - -# ╔═╡ 2806654f-7a24-4f8e-8b5a-3129b2431c95 -md""" -$$\begin{eqnarray*} -\cfrac{dS}{dt} &=& \cfrac{Q}{V} \left(S_{in} - S \right) - \mu_{max}\cfrac{S}{S + K_s} X\\ -\cfrac{dX}{dt} &=& -\cfrac{Q}{V} X + Y \mu_{max}\cfrac{S}{S + K_s} X -\end{eqnarray*}$$ -""" - -# ╔═╡ 5d6a881d-421e-4dea-9668-38b3f3378bf1 -# osys = missing - -# ╔═╡ 5802078d-12bc-489e-8566-15f77a11e8ba -md""" -Assume the uncertainties in the following parameter values: - -- $\mu_{max} = 0.40 \pm 0.06\;h^{-1}$ -- $K_s = 0.015 \pm 0.003 \;g/L$ -- $S_{in} = 0.022 \pm 0.004\;g/L$ - -and that the uncertainty in the other parameters values: $Y = 0.67$, $Q = 2.0\;L/h$, $V = 40.0\;L$ are negligible. 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$. Perform an uncertainty analysis by plotting the uncertainty bands on the simulation results of $S$ and $X$ in a timespan of $[0, 100]\,h$. - -Interpret your results. -""" - -# ╔═╡ 562e0006-bb79-4895-98b4-c376ea6ce856 -md""" -Initialize a vector `u0` with the initial conditions, and set the timespan: -""" - -# ╔═╡ 6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde -# u0 = missing # Uncomment and complete the instruction - -# ╔═╡ 1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b -# tspan = missing # Uncomment and complete the instruction - -# ╔═╡ 40a8c697-56d1-4488-b089-00ddf41cdb5a -md""" -We initialize a vector `params_uncert` with the parameter values and their corresponding uncertainty: -""" - -# ╔═╡ 62f4049b-617f-4268-beba-b2c634b46642 -# params_uncert = missing # Uncomment and complete the instruction - -# ╔═╡ 0e4c73d3-1630-463b-994d-073170c1a6c3 -md""" -We create the corresponding ODE problem and store it in `oprob_uncert`: -""" - -# ╔═╡ 3cad11d6-c2eb-4865-858a-d407aa25d780 -# oprob_uncert = missing # Uncomment and complete the instruction - -# ╔═╡ a7b3b68f-d34e-4e9a-aa44-224919c44fec -md""" -We solve the ODE problem. Use `Tsit5()` and `saveat=2.0`. Store the solution in `osol_uncert`: -""" - -# ╔═╡ 40168902-0780-4f92-8eca-23982c80511e -# osol_uncert = missing # Uncomment and complete the instruction - -# ╔═╡ 6b2abf2a-f1f3-4c2d-96e8-cc626848fa68 -md""" -Plot the results (simulation of the output variables $S$ and $X$ together with their uncertainty band): -""" - -# ╔═╡ 4d753fdb-089d-41dc-89ef-45a1733af786 -# missing - -# ╔═╡ 4f3900d7-4ff8-4b4b-a59a-f7079ac6c457 -md""" -Try to relate the local sensitivity analysis to the uncertainty analysis. Hence, study the effect of the individual parameter uncertainties on the output variables $S$ and $X$ and compare with your local sensitivity results of the corresponding parameter. - -In order to do that, analyse the effect on the uncertainty bands for $S$ and $X$ by taking one uncertainty on a parameter at a time. In other words, analyse the subsequent cases separately: -- Assume uncertainty only in $\mu_{max}$ -- Assume uncertainty only in $K_s$ -- Assume uncertainty only in $S_{in}$ -""" - -# ╔═╡ ff3e6f75-c02b-4530-8b9a-6ff40915cb5c -md""" -!!! question - Draw your conclusions. -""" - -# ╔═╡ 33b4d1c8-5da2-4490-afe9-4778e11d3cc9 -md"- Answer: missing" - -# ╔═╡ Cell order: -# ╠═1ec7def2-0c43-11ef-0850-ddefa680d2a4 -# ╠═2861bc17-7810-461e-8857-a55b6392977a -# ╠═9d8acf40-e635-4dc4-9938-ec63ab68e3bd -# ╠═093679cc-122b-4b0a-84c3-15cf3a6f5c1b -# ╠═0720732e-105a-40c2-bf45-f8fbbc10e5bc -# ╠═28bb60c9-ab56-497f-b039-0d6577501be5 -# ╟─437364b6-a953-428b-a891-e17ffd531656 -# ╟─c44e10f3-8452-4a60-b840-f4a7e6ec3238 -# ╟─356822df-bad3-41f9-8556-de04234b0a09 -# ╠═c551c33e-95b0-4ea8-85d4-6c68db72b8bc -# ╟─69f8c265-fe91-4e4b-971b-395325d4474e -# ╟─2806654f-7a24-4f8e-8b5a-3129b2431c95 -# ╠═5d6a881d-421e-4dea-9668-38b3f3378bf1 -# ╟─5802078d-12bc-489e-8566-15f77a11e8ba -# ╟─562e0006-bb79-4895-98b4-c376ea6ce856 -# ╠═6afb23f2-3da7-43e0-8a3b-fc4bb5eb9bde -# ╠═1b09c8cb-d3b7-4a15-b6ed-d03019c5f57b -# ╟─40a8c697-56d1-4488-b089-00ddf41cdb5a -# ╠═62f4049b-617f-4268-beba-b2c634b46642 -# ╟─0e4c73d3-1630-463b-994d-073170c1a6c3 -# ╠═3cad11d6-c2eb-4865-858a-d407aa25d780 -# ╟─a7b3b68f-d34e-4e9a-aa44-224919c44fec -# ╠═40168902-0780-4f92-8eca-23982c80511e -# ╟─6b2abf2a-f1f3-4c2d-96e8-cc626848fa68 -# ╠═4d753fdb-089d-41dc-89ef-45a1733af786 -# ╟─4f3900d7-4ff8-4b4b-a59a-f7079ac6c457 -# ╟─ff3e6f75-c02b-4530-8b9a-6ff40915cb5c -# ╠═33b4d1c8-5da2-4490-afe9-4778e11d3cc9 diff --git a/src/exercises/uncert_intro.jl b/src/exercises/uncert_intro.jl index f676d283..b50320af 100644 --- a/src/exercises/uncert_intro.jl +++ b/src/exercises/uncert_intro.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.21 +# v1.0.3 #> [frontmatter] #> order = "39" @@ -211,27 +211,11 @@ In uncertainty propagation, we are interested in the effect of uncertainties (or 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. @@ -245,6 +229,22 @@ Firstly, solve the original system (Wf = 10 and μ = 0.07) and afterwards we can # ╔═╡ 62afc495-29c3-4baa-b45a-94b84149d15d oprob_log = ODEProblem(growth_mod_log, u0_log,tspan, [:μ => 0.07, :Wf => 10]); +# ╔═╡ 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 + +# ╔═╡ 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); + # ╔═╡ 691a8b68-935e-4012-aa25-9cfdc323b1de osol_log = solve(oprob_log, Tsit5(), saveat = 0.5);