{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ADM1F SRT: Synthetic Model-Data Calibration \n",
"\n",
"Authors: Wenjuan Zhang and Elchin Jafarov"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before starting model clibration some additional softwares needs to be installed\n",
"\n",
"1. Install Julia \n",
"\n",
"2. Install Mads \n",
"\n",
"- \"import Pkg; Pkg.add(\"Mads\")\" -> this command should be enough\n",
"\n",
"- Mads github page is https://github.com/madsjulia/Mads.jl\n",
"\n",
"- Mads documentation documentation page is http://madsjulia.github.io/Mads.jl/\n",
"\n",
"3. To use jupyter notebook to display Julia, add JuliaHub, check this page https://juliahub.com/ui/Packages/IJulia/nfu7T/1.23.1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 1. Sensitivity of the model output to parameter changes \n",
"\n",
"Before starting model clibration we need to make sure that results are senstive to the parameter changes."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import Mads\n",
"import DelimitedFiles\n",
"import OrderedCollections"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\"/Users/elchin/project/ADM1F_WM/calibration\""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# check if you are in the right folder\n",
"pwd()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# navigate to the calibration folder\n",
"cd(\"/Users/elchin/project/ADM1F_WM/calibration\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{String,Any} with 4 entries:\n",
" \"Parameters\" => OrderedCollections.OrderedDict{String,OrderedCollections.O…\n",
" \"Observations\" => OrderedCollections.OrderedDict{String,OrderedCollections.O…\n",
" \"Julia command\" => \"./tim.jl\"\n",
" \"Filename\" => \"tim.mads\""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# passing \"tim.mads\" to the Mads. \"tim.mads\" includes the \"tim.jl\" that executes the \n",
"# ADM1F model and saves seven outputs of interest.\n",
"filename = \"tim.mads\" \n",
"md = Mads.loadmadsfile(filename)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OrderedCollections.OrderedDict{Any,Float64} with 7 entries:\n",
" \"o2\" => 145.362\n",
" \"o3\" => 48.578\n",
" \"o1\" => 363.333\n",
" \"o4\" => 35.1163\n",
" \"o5\" => 11.6628\n",
" \"o6\" => 6.51573\n",
" \"o7\" => 3.04116e-312"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# run the model with the default pamaters from \"tim.mads\"\n",
"# we denote the synthetic truth/observations \n",
"output_truth=Mads.forward(md) "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OrderedCollections.OrderedDict{String,Any} with 4 entries:\n",
" \"init\" => 0.5\n",
" \"max\" => 0.6\n",
" \"min\" => 0.4\n",
" \"type\" => \"opt\""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the variables in the dictinary that have \"type\" equal to \"opt\" will participate in the \n",
"# model-data calibration\n",
"md[\"Parameters\"][\"p48\"]#[\"init\"]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"o1 target = 43.524 weight = 1\n",
"o2 target = 19.1588 weight = 1\n",
"o3 target = 6.23771 weight = 1\n",
"o4 target = 0.00779 weight = 1\n",
"o5 target = 0 weight = 1\n",
"o6 target = 6.87726 weight = 1\n",
"o7 target = 0.63904 weight = 1\n",
"Number of observations is 7\n"
]
}
],
"source": [
"Mads.showobservations(md)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.04116075262e-312"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# set observations equal to model runs with md[\"Parameters\"][\"p48\"][init]=0.5\n",
"md[\"Observations\"][\"o1\"][\"target\"]=output_truth[\"o1\"]\n",
"md[\"Observations\"][\"o2\"][\"target\"]=output_truth[\"o2\"]\n",
"md[\"Observations\"][\"o3\"][\"target\"]=output_truth[\"o3\"]\n",
"md[\"Observations\"][\"o4\"][\"target\"]=output_truth[\"o4\"]\n",
"md[\"Observations\"][\"o5\"][\"target\"]=output_truth[\"o5\"]\n",
"md[\"Observations\"][\"o6\"][\"target\"]=output_truth[\"o6\"]\n",
"md[\"Observations\"][\"o7\"][\"target\"]=output_truth[\"o7\"]"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"o1 target = 363.333 weight = 1\n",
"o2 target = 145.362 weight = 1\n",
"o3 target = 48.578 weight = 1\n",
"o4 target = 35.1163 weight = 1\n",
"o5 target = 11.6628 weight = 1\n",
"o6 target = 6.51573 weight = 1\n",
"o7 target = 3.04116e-312 weight = 1\n",
"Number of observations is 7\n"
]
}
],
"source": [
"# new synthetic observations\n",
"Mads.showobservations(md)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"48\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.41,\"max\" => 0.6,\"min\" => 0.4,\"type\" => \"opt\")\n",
"54\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.031,\"max\" => 0.045599999999999995,\"min\" => 0.0304,\"type\" => \"opt\")\n",
"60\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.013,\"max\" => 0.018,\"min\" => 0.012,\"type\" => \"opt\")\n"
]
}
],
"source": [
"#find out how many parameters are opt-in\n",
"for i in 1:100\n",
" s1=string.(i)\n",
" s2=\"p\"\n",
" if md[\"Parameters\"][s2*s1][\"type\"]==\"opt\"\n",
" println(i)\n",
" println(md[\"Parameters\"][s2*s1]) \n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Two paramaters are opt-in from the tim.mads. This file can be changed manually or here by chaning the initial value (`init`) to a slightly different one. "
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OrderedCollections.OrderedDict{String,Any} with 4 entries:\n",
" \"init\" => 0.032\n",
" \"max\" => 0.0456\n",
" \"min\" => 0.0304\n",
" \"type\" => \"opt\""
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#let's opt-in a new parameter\n",
"md[\"Parameters\"][\"p54\"][\"type\"]=\"opt\" #K_S_fa\n",
"md[\"Parameters\"][\"p54\"]"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.55,\"max\" => 0.6,\"min\" => 0.4,\"type\" => \"opt\")\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.017,\"max\" => 0.018,\"min\" => 0.012,\"type\" => \"opt\")\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.032,\"max\" => 0.045599999999999995,\"min\" => 0.0304,\"type\" => \"opt\")\n"
]
}
],
"source": [
"println(md[\"Parameters\"][\"p48\"])\n",
"println(md[\"Parameters\"][\"p60\"])\n",
"println(md[\"Parameters\"][\"p54\"])"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.58,\"max\" => 0.6,\"min\" => 0.4,\"type\" => \"opt\")\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.0165,\"max\" => 0.018,\"min\" => 0.012,\"type\" => \"opt\")\n",
"OrderedCollections.OrderedDict{String,Any}(\"init\" => 0.045,\"max\" => 0.045599999999999995,\"min\" => 0.0304,\"type\" => \"opt\")\n"
]
}
],
"source": [
"#change to new initial values\n",
"md[\"Parameters\"][\"p48\"][\"init\"]=0.58\n",
"md[\"Parameters\"][\"p60\"][\"init\"]=0.0165 #K_S_pro\n",
"md[\"Parameters\"][\"p54\"][\"init\"]=0.045\n",
"println(md[\"Parameters\"][\"p48\"])\n",
"println(md[\"Parameters\"][\"p60\"])\n",
"println(md[\"Parameters\"][\"p54\"])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.35"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#we can also change the boundary of the allowed interval \n",
"#md[\"Parameters\"][\"p48\"][\"init\"]=0.59\n",
"#md[\"Parameters\"][\"p48\"][\"max\"]=0.7\n",
"#md[\"Parameters\"][\"p48\"][\"min\"]=0.35"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OrderedCollections.OrderedDict{Any,Float64} with 7 entries:\n",
" \"o2\" => 145.06\n",
" \"o3\" => 48.3842\n",
" \"o1\" => 363.488\n",
" \"o4\" => 35.0174\n",
" \"o5\" => 11.6574\n",
" \"o6\" => 6.51656\n",
" \"o7\" => -1.14602e-312"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#rerun the model and compare the results with the previous forward run\n",
"output1=Mads.forward(md) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\"o2\" => 145.362\n",
"\"o3\" => 48.578\n",
"\"o1\" => 363.333\n",
"\"o4\" => 35.1163\n",
"\"o5\" => 11.6628\n",
"\"o6\" => 6.51573\n",
"\"o7\" => 3.04116e-312"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 2. Model Calibration \n",
"\n",
"Once the sensitivity of the model outputs to changes in model parameters are established we can start the calibration. The elaborated sensitivy analysis is required, for an in-depth model calibration. "
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(OrderedCollections.OrderedDict(\"p1\" => 0.09060657861205099,\"p2\" => 0.2964692632944874,\"p3\" => 0.18538822521329718,\"p4\" => 0.21021796952211014,\"p5\" => 0.20993757892890663,\"p6\" => 0.0025810430070856297,\"p7\" => 0.004153144361053321,\"p8\" => 0.006812701647314224,\"p9\" => 0.029624259152454625,\"p10\" => 0.028993549342676592…), OptimBase.MultivariateOptimizationResults{LsqFit.LevenbergMarquardt,Float64,1}(LsqFit.LevenbergMarquardt(), [-0.4889181593485484, 1.1927438925052574, -0.3739490066652115, 0.2583123288342893, -0.9293788161320726, -0.8970848356758778, -0.6686810831789546, -0.5646712759835035, 0.3221739693428631, -0.16853854003002233 … -0.515933804055272, -0.9755102180590636, -1.1189154560307426, 1.0593944723620046, 0.1384639308865269, -0.5026261086108612, -0.012128162836436718, 0.28788576425004025, -0.0029143278533535384, -1.1138366569543616], [-0.4889181593485484, 1.1927438925052574, -0.3739490066652115, 0.2583123288342893, -0.9293788161320726, -0.8970848356758778, -0.6686810831789546, -0.5646712759835035, 0.3221739693428631, -0.16853854003002233 … -0.515933804055272, -0.9755102180590636, -1.1189154560307426, 1.0593944723620046, 0.1384639308865269, -0.5026261086108612, -0.012128162836436718, 0.28788576425004025, -0.0029143278533535384, -1.1138366569543616], 150.99483823310115, 10, true, false, 0.0001, 0.0, false, 0.001, 0.0, false, 1.0e-6, 0.0, false, Iter Function value Gradient norm \n",
"------ -------------- --------------\n",
", 1101, 10, 0))"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p, r = Mads.calibraterandom(md) "
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"o1 target = 363.333 weight = 1\n",
"o2 target = 145.362 weight = 1\n",
"o3 target = 48.578 weight = 1\n",
"o4 target = 35.1163 weight = 1\n",
"o5 target = 11.6628 weight = 1\n",
"o6 target = 6.51573 weight = 1\n",
"o7 target = 3.04116e-312 weight = 1\n",
"Number of observations is 7\n"
]
}
],
"source": [
"Mads.showobservations(md)"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"initial guess:\n",
"p48=0.58\n",
"p60=0.0165\n",
"p54=0.045\n",
" \n",
"calibrated parameters:\n",
"p48=0.486\n",
"p60=0.013\n",
"p54=0.039\n",
" \n",
"true parameters:\n",
"p48=0.41\n",
"p60=0.013\n",
"p54=0.031\n"
]
}
],
"source": [
"println(\"initial guess:\")\n",
"println(\"p48=\",md[\"Parameters\"][\"p48\"][\"init\"])\n",
"println(\"p60=\",md[\"Parameters\"][\"p60\"][\"init\"])\n",
"println(\"p54=\",md[\"Parameters\"][\"p54\"][\"init\"])\n",
"println(\" \")\n",
"println(\"calibrated parameters:\")\n",
"println(\"p48=\",round(p[\"p48\"],digits=3))\n",
"println(\"p60=\",round(p[\"p60\"],digits=3))\n",
"println(\"p54=\",round(p[\"p54\"],digits=3))\n",
"println(\" \")\n",
"println(\"true parameters:\")\n",
"println(\"p48=\",0.41)\n",
"println(\"p60=\",0.013)\n",
"println(\"p54=\",0.031)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\r"
]
}
],
"source": [
"#Only plot the outputs which participate in this match\n",
"\n",
"est = Mads.forward(md, p);\n",
"obs = Mads.getobstarget(md);\n",
"\n",
"idx_out = collect(1:1:7)\n",
"key_out = [string('o',i) for i in idx_out]\n",
"\n",
"est_filt = [est[i] for i in key_out];\n",
"obs_filt = [obs[i] for i in idx_out];\n",
"Mads.plotseries([obs_filt est_filt]; names=[\"Truth\", \"Est.\"])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}