ADM1F SRT: Synthetic Model-Data Calibration¶
Authors: Wenjuan Zhang and Elchin Jafarov
Before starting model clibration some additional softwares needs to be installed
Install Julia
Install Mads
“import Pkg; Pkg.add(”Mads“)” -> this command should be enough
Mads github page is https://github.com/madsjulia/Mads.jl
Mads documentation documentation page is http://madsjulia.github.io/Mads.jl/
To use jupyter notebook to display Julia, add JuliaHub, check this page https://juliahub.com/ui/Packages/IJulia/nfu7T/1.23.1
1. Sensitivity of the model output to parameter changes¶
Before starting model clibration we need to make sure that results are senstive to the parameter changes.
[2]:
import Mads
import DelimitedFiles
import OrderedCollections
[3]:
# check if you are in the right folder
pwd()
[3]:
"/Users/elchin/project/ADM1F_WM/calibration"
[4]:
# navigate to the calibration folder
cd("/Users/elchin/project/ADM1F_WM/calibration")
[5]:
# passing "tim.mads" to the Mads. "tim.mads" includes the "tim.jl" that executes the
# ADM1F model and saves seven outputs of interest.
filename = "tim.mads"
md = Mads.loadmadsfile(filename)
[5]:
Dict{String,Any} with 4 entries:
"Parameters" => OrderedCollections.OrderedDict{String,OrderedCollections.O…
"Observations" => OrderedCollections.OrderedDict{String,OrderedCollections.O…
"Julia command" => "./tim.jl"
"Filename" => "tim.mads"
[6]:
# run the model with the default pamaters from "tim.mads"
# we denote the synthetic truth/observations
output_truth=Mads.forward(md)
[6]:
OrderedCollections.OrderedDict{Any,Float64} with 7 entries:
"o2" => 145.362
"o3" => 48.578
"o1" => 363.333
"o4" => 35.1163
"o5" => 11.6628
"o6" => 6.51573
"o7" => 3.04116e-312
[7]:
# the variables in the dictinary that have "type" equal to "opt" will participate in the
# model-data calibration
md["Parameters"]["p48"]#["init"]
[7]:
OrderedCollections.OrderedDict{String,Any} with 4 entries:
"init" => 0.5
"max" => 0.6
"min" => 0.4
"type" => "opt"
[8]:
Mads.showobservations(md)
o1 target = 43.524 weight = 1
o2 target = 19.1588 weight = 1
o3 target = 6.23771 weight = 1
o4 target = 0.00779 weight = 1
o5 target = 0 weight = 1
o6 target = 6.87726 weight = 1
o7 target = 0.63904 weight = 1
Number of observations is 7
[10]:
# set observations equal to model runs with md["Parameters"]["p48"][init]=0.5
md["Observations"]["o1"]["target"]=output_truth["o1"]
md["Observations"]["o2"]["target"]=output_truth["o2"]
md["Observations"]["o3"]["target"]=output_truth["o3"]
md["Observations"]["o4"]["target"]=output_truth["o4"]
md["Observations"]["o5"]["target"]=output_truth["o5"]
md["Observations"]["o6"]["target"]=output_truth["o6"]
md["Observations"]["o7"]["target"]=output_truth["o7"]
[10]:
3.04116075262e-312
[11]:
# new synthetic observations
Mads.showobservations(md)
o1 target = 363.333 weight = 1
o2 target = 145.362 weight = 1
o3 target = 48.578 weight = 1
o4 target = 35.1163 weight = 1
o5 target = 11.6628 weight = 1
o6 target = 6.51573 weight = 1
o7 target = 3.04116e-312 weight = 1
Number of observations is 7
[27]:
#find out how many parameters are opt-in
for i in 1:100
s1=string.(i)
s2="p"
if md["Parameters"][s2*s1]["type"]=="opt"
println(i)
println(md["Parameters"][s2*s1])
end
end
48
OrderedCollections.OrderedDict{String,Any}("init" => 0.41,"max" => 0.6,"min" => 0.4,"type" => "opt")
54
OrderedCollections.OrderedDict{String,Any}("init" => 0.031,"max" => 0.045599999999999995,"min" => 0.0304,"type" => "opt")
60
OrderedCollections.OrderedDict{String,Any}("init" => 0.013,"max" => 0.018,"min" => 0.012,"type" => "opt")
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.
[20]:
#let's opt-in a new parameter
md["Parameters"]["p54"]["type"]="opt" #K_S_fa
md["Parameters"]["p54"]
[20]:
OrderedCollections.OrderedDict{String,Any} with 4 entries:
"init" => 0.032
"max" => 0.0456
"min" => 0.0304
"type" => "opt"
[21]:
println(md["Parameters"]["p48"])
println(md["Parameters"]["p60"])
println(md["Parameters"]["p54"])
OrderedCollections.OrderedDict{String,Any}("init" => 0.55,"max" => 0.6,"min" => 0.4,"type" => "opt")
OrderedCollections.OrderedDict{String,Any}("init" => 0.017,"max" => 0.018,"min" => 0.012,"type" => "opt")
OrderedCollections.OrderedDict{String,Any}("init" => 0.032,"max" => 0.045599999999999995,"min" => 0.0304,"type" => "opt")
[31]:
#change to new initial values
md["Parameters"]["p48"]["init"]=0.58
md["Parameters"]["p60"]["init"]=0.0165 #K_S_pro
md["Parameters"]["p54"]["init"]=0.045
println(md["Parameters"]["p48"])
println(md["Parameters"]["p60"])
println(md["Parameters"]["p54"])
OrderedCollections.OrderedDict{String,Any}("init" => 0.58,"max" => 0.6,"min" => 0.4,"type" => "opt")
OrderedCollections.OrderedDict{String,Any}("init" => 0.0165,"max" => 0.018,"min" => 0.012,"type" => "opt")
OrderedCollections.OrderedDict{String,Any}("init" => 0.045,"max" => 0.045599999999999995,"min" => 0.0304,"type" => "opt")
[9]:
#we can also change the boundary of the allowed interval
#md["Parameters"]["p48"]["init"]=0.59
#md["Parameters"]["p48"]["max"]=0.7
#md["Parameters"]["p48"]["min"]=0.35
[9]:
0.35
[32]:
#rerun the model and compare the results with the previous forward run
output1=Mads.forward(md)
[32]:
OrderedCollections.OrderedDict{Any,Float64} with 7 entries:
"o2" => 145.06
"o3" => 48.3842
"o1" => 363.488
"o4" => 35.0174
"o5" => 11.6574
"o6" => 6.51656
"o7" => -1.14602e-312
“o2” => 145.362 “o3” => 48.578 “o1” => 363.333 “o4” => 35.1163 “o5” => 11.6628 “o6” => 6.51573 “o7” => 3.04116e-312
2. Model Calibration¶
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.
[34]:
p, r = Mads.calibraterandom(md)
[34]:
(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
------ -------------- --------------
, 1101, 10, 0))
[35]:
Mads.showobservations(md)
o1 target = 363.333 weight = 1
o2 target = 145.362 weight = 1
o3 target = 48.578 weight = 1
o4 target = 35.1163 weight = 1
o5 target = 11.6628 weight = 1
o6 target = 6.51573 weight = 1
o7 target = 3.04116e-312 weight = 1
Number of observations is 7
[52]:
println("initial guess:")
println("p48=",md["Parameters"]["p48"]["init"])
println("p60=",md["Parameters"]["p60"]["init"])
println("p54=",md["Parameters"]["p54"]["init"])
println(" ")
println("calibrated parameters:")
println("p48=",round(p["p48"],digits=3))
println("p60=",round(p["p60"],digits=3))
println("p54=",round(p["p54"],digits=3))
println(" ")
println("true parameters:")
println("p48=",0.41)
println("p60=",0.013)
println("p54=",0.031)
initial guess:
p48=0.58
p60=0.0165
p54=0.045
calibrated parameters:
p48=0.486
p60=0.013
p54=0.039
true parameters:
p48=0.41
p60=0.013
p54=0.031
[51]:
#Only plot the outputs which participate in this match
est = Mads.forward(md, p);
obs = Mads.getobstarget(md);
idx_out = collect(1:1:7)
key_out = [string('o',i) for i in idx_out]
est_filt = [est[i] for i in key_out];
obs_filt = [obs[i] for i in idx_out];
Mads.plotseries([obs_filt est_filt]; names=["Truth", "Est."])
[ ]: