Trying to setup a Parameter Estimation in Julia with DifferentialEquations, I started reading the documentation and some topics here on Discourse. But it seems I am actually facing three issues at once, which is why I am opening this topic.
I have a simple ODE with two states. Input parameters change at certain times, which I am handling through PresetTimeCallback, with additional suggestions from here. (Using Interpolations.jl instead didn’t work, since previous interpolation is yet to be implemented.)
Consider this MWE:
using DifferentialEquations
using DelimitedFiles
using Plots
function odm2prod(dx, x, adjParams, t, fixedParams)
V_liq, X_in, Y_in = fixedParams
k_1, f_1, q_in = adjParams
rho_1 = k_1*x[1]
# q_prod = 0.52*f_1*x[1]
# Differential Equations
dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
dx[2] = q_in/V_liq*(Y_in - x[2])
end
x0 = [11.5, 11.55*0.2]
tspan = (0.0, 7.0)
fixedParams = [159, 249, 58]
adjParams = [0.19, 29.0, 0]
prob = ODEProblem((dx, x, adjParams, t) -> odm2prod(dx, x, adjParams,
t, fixedParams), x0, tspan, adjParams)
inputArray = readdlm("rawDataFiles/Inputarray_Odm2prod_2h_7d.csv", ',',
Float64)
dosetimes = inputArray[:,1]/86400
function affect!(integrator)
ind_t = findall(t -> t==integrator.t, dosetimes)
integrator.p[3] = inputArray[ind_t[1], 2]
end
cb = PresetTimeCallback(dosetimes, affect!)
sol = solve(prob, Tsit5(), callback=cb)
(Actually, X_in and Y_in should also change at every callback, but I did not implement this due to problem #3, see below)
Content of input CSV file
3300,430,32.64,32.64
3521.0232558,0,32.64,32.64
3660,60,249.83539555,58.001892497
3828.0000005,0,249.83539555,58.001892497
10860,60,249.83539555,58.001892497
11028,0,249.83539555,58.001892497
18060.005,60,249.83539555,58.001892497
18228.005,0,249.83539555,58.001892497
25260.01,60,249.83539555,58.001892497
25428.01,0,249.83539555,58.001892497
32460.015,60,249.83539555,58.001892497
32628.015,0,249.83539555,58.001892497
39660.02,60,249.83539555,58.001892497
39828.02,0,249.83539555,58.001892497
46860.025,60,249.83539555,58.001892497
47028.025,0,249.83539555,58.001892497
54060.03,60,249.83539555,58.001892497
54228.03,0,249.83539555,58.001892497
61260.035,60,249.83539555,58.001892497
61428.035,0,249.83539555,58.001892497
68460.04,60,249.83539555,58.001892497
68628.04,0,249.83539555,58.001892497
75660.045,60,249.83539555,58.001892497
75828.045,0,249.83539555,58.001892497
82860.05,60,249.83539555,58.001892497
83028.05,0,249.83539555,58.001892497
89700,430,32.64,32.64
89921.023256,0,32.64,32.64
90060,60,249.83539555,58.001892497
90228,0,249.83539555,58.001892497
97260,60,249.83539555,58.001892497
97428,0,249.83539555,58.001892497
104460.005,60,249.83539555,58.001892497
104628.005,0,249.83539555,58.001892497
111660.01,60,249.83539555,58.001892497
111828.01,0,249.83539555,58.001892497
118860.015,60,249.83539555,58.001892497
119028.015,0,249.83539555,58.001892497
126060.02,60,249.83539555,58.001892497
126228.02,0,249.83539555,58.001892497
133260.025,60,249.83539555,58.001892497
133428.025,0,249.83539555,58.001892497
140460.03,60,249.83539555,58.001892497
140628.03,0,249.83539555,58.001892497
147660.035,60,249.83539555,58.001892497
147828.035,0,249.83539555,58.001892497
154860.04,60,249.83539555,58.001892497
155028.04,0,249.83539555,58.001892497
162060.045,60,249.83539555,58.001892497
162228.045,0,249.83539555,58.001892497
169260.05,60,249.83539555,58.001892497
169428.05,0,249.83539555,58.001892497
176100,430,32.64,32.64
176321.02326,0,32.64,32.64
176460,60,249.83539555,58.001892497
176628,0,249.83539555,58.001892497
183660,60,249.83539555,58.001892497
183828,0,249.83539555,58.001892497
190860.005,60,249.83539555,58.001892497
191028.005,0,249.83539555,58.001892497
198060.01,60,249.83539555,58.001892497
198228.01,0,249.83539555,58.001892497
205260.015,60,249.83539555,58.001892497
205428.015,0,249.83539555,58.001892497
212460.02,60,249.83539555,58.001892497
212628.02,0,249.83539555,58.001892497
219660.025,60,249.83539555,58.001892497
219828.025,0,249.83539555,58.001892497
226860.03,60,249.83539555,58.001892497
227028.03,0,249.83539555,58.001892497
234060.035,60,249.83539555,58.001892497
234228.035,0,249.83539555,58.001892497
241260.04,60,249.83539555,58.001892497
241428.04,0,249.83539555,58.001892497
248460.045,60,249.83539555,58.001892497
248628.045,0,249.83539555,58.001892497
255660.05,60,249.83539555,58.001892497
255828.05,0,249.83539555,58.001892497
262500,430,32.64,32.64
262721.02326,0,32.64,32.64
262860,60,249.83539555,58.001892497
263028,0,249.83539555,58.001892497
270060,60,249.83539555,58.001892497
270228,0,249.83539555,58.001892497
277260.005,60,249.83539555,58.001892497
277428.005,0,249.83539555,58.001892497
284460.01,60,249.83539555,58.001892497
284628.01,0,249.83539555,58.001892497
291660.015,60,249.83539555,58.001892497
291828.015,0,249.83539555,58.001892497
298860.02,60,249.83539555,58.001892497
299028.02,0,249.83539555,58.001892497
306060.025,60,249.83539555,58.001892497
306228.025,0,249.83539555,58.001892497
313260.03,60,249.83539555,58.001892497
313428.03,0,249.83539555,58.001892497
320460.035,60,249.83539555,58.001892497
320628.035,0,249.83539555,58.001892497
327660.04,60,249.83539555,58.001892497
327828.04,0,249.83539555,58.001892497
334860.045,60,249.83539555,58.001892497
335028.045,0,249.83539555,58.001892497
342060.05,60,249.83539555,58.001892497
342228.05,0,249.83539555,58.001892497
348900,430,32.64,32.64
349121.02326,0,32.64,32.64
349260,60,249.83539555,58.001892497
349428,0,249.83539555,58.001892497
356460,60,249.83539555,58.001892497
356628,0,249.83539555,58.001892497
363660.005,60,249.83539555,58.001892497
363828.005,0,249.83539555,58.001892497
370860.01,60,249.83539555,58.001892497
371028.01,0,249.83539555,58.001892497
378060.015,60,249.83539555,58.001892497
378228.015,0,249.83539555,58.001892497
385260.02,60,249.83539555,58.001892497
385428.02,0,249.83539555,58.001892497
392460.025,60,249.83539555,58.001892497
392628.025,0,249.83539555,58.001892497
399660.03,60,249.83539555,58.001892497
399828.03,0,249.83539555,58.001892497
406860.035,60,249.83539555,58.001892497
407028.035,0,249.83539555,58.001892497
414060.04,60,249.83539555,58.001892497
414228.04,0,249.83539555,58.001892497
421260.045,60,249.83539555,58.001892497
421428.045,0,249.83539555,58.001892497
428460.05,60,249.83539555,58.001892497
428628.05,0,249.83539555,58.001892497
435300,430,32.64,32.64
435521.02326,0,32.64,32.64
435660,60,249.83539555,58.001892497
435828,0,249.83539555,58.001892497
442860,60,249.83539555,58.001892497
443028,0,249.83539555,58.001892497
450060.005,60,249.83539555,58.001892497
450228.005,0,249.83539555,58.001892497
457260.01,60,249.83539555,58.001892497
457428.01,0,249.83539555,58.001892497
464460.015,60,249.83539555,58.001892497
464628.015,0,249.83539555,58.001892497
471660.02,60,249.83539555,58.001892497
471828.02,0,249.83539555,58.001892497
478860.025,60,249.83539555,58.001892497
479028.025,0,249.83539555,58.001892497
486060.03,60,249.83539555,58.001892497
486228.03,0,249.83539555,58.001892497
493260.035,60,249.83539555,58.001892497
493428.035,0,249.83539555,58.001892497
500460.04,60,249.83539555,58.001892497
500628.04,0,249.83539555,58.001892497
507660.045,60,249.83539555,58.001892497
507828.045,0,249.83539555,58.001892497
514860.05,60,249.83539555,58.001892497
515028.05,0,249.83539555,58.001892497
521700,430,32.64,32.64
521921.02326,0,32.64,32.64
522060,60,249.83539555,58.001892497
522228,0,249.83539555,58.001892497
529260,60,249.83539555,58.001892497
529428,0,249.83539555,58.001892497
536460.005,60,249.83539555,58.001892497
536628.005,0,249.83539555,58.001892497
543660.01,60,249.83539555,58.001892497
543828.01,0,249.83539555,58.001892497
550860.015,60,249.83539555,58.001892497
551028.015,0,249.83539555,58.001892497
558060.02,60,249.83539555,58.001892497
558228.02,0,249.83539555,58.001892497
565260.025,60,249.83539555,58.001892497
565428.025,0,249.83539555,58.001892497
572460.03,60,249.83539555,58.001892497
572628.03,0,249.83539555,58.001892497
579660.035,60,249.83539555,58.001892497
579828.035,0,249.83539555,58.001892497
586860.04,60,249.83539555,58.001892497
587028.04,0,249.83539555,58.001892497
594060.045,60,249.83539555,58.001892497
594228.045,0,249.83539555,58.001892497
601260.05,60,249.83539555,58.001892497
601428.05,0,249.83539555,58.001892497
For completeness, here are the observations, with the columns t (in seconds; needs to be divided by 86400, in order to match unit of tspan), x[1], q_prod (needs to be multiplied by 24, in order to match unit of q_prod in code):
Content of measurements CSV file
0,11.3526,14.02611
3600,11.41204,5.990548
7200,11.48661,14.13363
10800,11.39558,13.81184
14400,11.47054,14.30515
18000,11.37988,13.85038
21600,11.45519,14.30656
25200,11.36489,13.83604
28800,11.44058,14.28658
32400,11.35062,13.81046
36000,11.42665,14.25907
39600,11.33701,13.77983
43200,11.41337,14.22781
46800,11.32404,13.74663
50400,11.40072,14.19458
54000,11.31169,13.71217
57600,11.38867,14.16042
61200,11.29993,13.67726
64800,11.37718,14.12593
68400,11.28871,13.64237
72000,11.36625,14.09166
75600,11.27804,13.6079
79200,11.35585,14.05792
82800,11.26789,13.57412
86400,11.34595,14.02487
90000,11.40542,5.988436
93600,11.48,14.13261
97200,11.38898,13.81054
100800,11.46395,14.30391
104400,11.3733,13.84899
108000,11.44865,14.30534
111600,11.35835,13.8347
115200,11.43404,14.28531
118800,11.34408,13.8091
122400,11.42012,14.25784
126000,11.33049,13.77849
129600,11.40687,14.22666
133200,11.31755,13.74535
136800,11.39422,14.19333
140400,11.30519,13.71084
144000,11.38219,14.15921
147600,11.29345,13.67594
151200,11.37071,14.12471
154800,11.28225,13.64105
158400,11.3598,14.09051
162000,11.2716,13.60663
165600,11.34943,14.05678
169200,11.26147,13.57286
172800,11.33954,14.02368
176400,11.39906,5.985696
180000,11.47366,14.13173
183600,11.38264,13.80937
187200,11.45763,14.3028
190800,11.36698,13.84773
194400,11.44234,14.30423
198000,11.35205,13.83347
201600,11.42775,14.28423
205200,11.3378,13.80791
208800,11.41384,14.25672
212400,11.32421,13.77728
216000,11.40058,14.22547
219600,11.31127,13.74408
223200,11.38796,14.19223
226800,11.29894,13.70962
230400,11.37593,14.15805
234000,11.28719,13.67469
237600,11.36448,14.12362
241200,11.27602,13.63984
244800,11.35358,14.08939
248400,11.26539,13.60541
252000,11.34322,14.05565
255600,11.25527,13.57164
259200,11.33335,14.02261
262800,11.39295,5.979839
266400,11.46755,14.13088
270000,11.37654,13.8083
273600,11.45153,14.30173
277200,11.36089,13.84653
280800,11.43624,14.30311
284400,11.34595,13.83227
288000,11.42167,14.28316
291600,11.33172,13.80673
295200,11.40778,14.25573
298800,11.31816,13.77616
302400,11.39455,14.2245
306000,11.30523,13.74299
309600,11.38193,14.19126
313200,11.29292,13.70855
316800,11.36992,14.1571
320400,11.28119,13.67363
324000,11.35848,14.12262
327600,11.27002,13.63875
331200,11.34759,14.08839
334800,11.2594,13.60431
338400,11.33724,14.05465
342000,11.24929,13.57054
345600,11.32737,14.02153
349200,11.38695,5.98036
352800,11.46156,14.12997
356400,11.37055,13.80708
360000,11.44554,14.30051
363600,11.3549,13.84516
367200,11.43028,14.3019
370800,11.34,13.83092
374400,11.41572,14.28194
378000,11.32578,13.8054
381600,11.40186,14.25456
385200,11.31224,13.77487
388800,11.38864,14.22333
392400,11.29933,13.74171
396000,11.37604,14.19009
399600,11.28703,13.70726
403200,11.36404,14.15591
406800,11.27531,13.67234
410400,11.35261,14.12148
414000,11.26416,13.63749
417600,11.34174,14.08725
421200,11.25356,13.60306
424800,11.33139,14.05345
428400,11.24345,13.56925
432000,11.32154,14.02036
435600,11.38121,5.97516
439200,11.45584,14.12915
442800,11.36483,13.80602
446400,11.43984,14.2995
450000,11.3492,13.84401
453600,11.42458,14.30084
457200,11.3343,13.82977
460800,11.41004,14.28096
464400,11.3201,13.80431
468000,11.39617,14.25352
471600,11.30656,13.77375
475200,11.38297,14.22233
478800,11.29366,13.7406
482400,11.37037,14.18904
486000,11.28136,13.70612
489600,11.35837,14.15484
493200,11.26965,13.67118
496800,11.34695,14.12039
500400,11.2585,13.63631
504000,11.33609,14.08616
507600,11.24791,13.60188
511200,11.32576,14.05241
514800,11.23783,13.5681
518400,11.31593,14.01936
522000,11.37565,5.971109
525600,11.45027,14.12832
529200,11.35927,13.80494
532800,11.43428,14.29842
536400,11.34365,13.84281
540000,11.41904,14.29982
543600,11.32877,13.82862
547200,11.4045,14.27986
550800,11.31457,13.80312
554400,11.39064,14.25243
558000,11.30103,13.77256
561600,11.37744,14.22124
565200,11.28814,13.73942
568800,11.36487,14.18803
572400,11.27587,13.70499
576000,11.3529,14.15386
579600,11.26418,13.67009
583200,11.3415,14.11943
586800,11.25306,13.63524
590400,11.33065,14.08518
594000,11.24247,13.6008
597600,11.32032,14.05139
601200,11.23239,13.567
604800,11.3105,14.01833
Now, since my observations are for x[1] and q_prod, I need q_prod as output variable from the ODEProblem. I managed to write this function, which I am not sure will be of any help for this endeavor:
# Out-function for odm2prod to get algebraic states
function odm2prod_out(sol, adjParams, fixedParams)
V_liq, X_in, Y_in = fixedParams
k_1, f_1, q_in = adjParams
x = Array(sol)
# Algebraic equation
q_prod = 0.52*f_1*x[1,:]
return q_prod
end
y = odm2prod_out(sol, adjParams, fixedParams)
plot(sol.t, y/24)
So these are my problems:
- How to get an algebraic output variable, which is dependent on the states and parameters, as part of
sol? - While the callbacks may occur at any time, the observations are equally spaced, at fixed interval, so the number and location of points I get in
solare not the same as my measured times, even with thesaveat=option enabled (I checked). How would I reconcile both timescales? - If I understood correctly, all parameters in
adjParamswould be adjusted by Parameter Estimation routine, but some must be adjusted by the callback. How would I separate them?
Thanks in advance for any help, it is very much appreciated.
