# Parameter Estimate Problem

I find the IDA() can’t converge,the program keeping runing ,can’t not stop,how to select DAE algorithm to solve the prob.

``````using Plots,DifferentialEquations,SmoothingSplines,Optim,Sundials
using DiffEqParamEstim
t=[ 0,
16.78100,35.70300,53.78100,	71.92200,90.09400,108.2810,	126.4530,144.6410,162.8440,	181.0160,
199.2190,217.3910,235.5940,	253.7500,271.9840,290.1410,	308.3590,326.5000,344.7500,	362.9060,
381.0470,399.1870,417.2810,	435.4840,453.6720,471.8120,	490.2030,508.3440,526.6250,	544.9220,
563.1250,581.4060,599.6560,	617.8750,636.0470,654.4220,	672.5940,690.8120,708.9370,	727.2660,
745.5780,763.7970,782.1090,	800.4840,818.7970,837.0160,	855.2500,873.5780,891.8120,	910.1410,
928.3910,946.7190,965.0470,	983.3910,1001.766,1020.172,	1038.594,1056.922,1075.203,	1093.687, 1112,
1130.328,1148.750,1167.125,	1185.641,1203.984,1222.406,	1240.797,1259.156,1277.547,	1295.906,
1314.297,1332.687,1351.203,	1369.687,1388.125,1406.578,	1425.062,1443.609,1462.047,	1480.547,
1499.031,1517.672,1536.109,	1554.781,1573.281,1591.859,	1610.406,1628.953,1647.562,	1666.062,
1684.656,1703.219,1721.891,	1740.437,1759.016,1777.641,	1796.328,1815.047,1833.750,	1852.391,
1871.078,1889.844,1908.578,	1927.281,1946.031,1964.891,	1983.703,2002.484,2021.187,	2039.906,
2058.641,2077.406,2096.312,	2115.047,2133.969,2152.812,	2171.625,2190.547,2209.422,	2228.297,
2247.219,2265.969,2284.937,	2303.859,2322.797,2341.766,	2360.594,2379.578,2398.672,	2417.797,
2436.906,2455.969,2474.984,	2494.031,2513.187,2532.391,	2551.547,2570.578,2589.703,	2608.797,
2627.953,2647.234,2666.516,	2685.812,2704.969,2724.203,	2743.453,2762.625,2781.922,	2801.250,
2820.625,2839.937,2859.234,	2878.578,2897.937,2917.281,	2936.687,2956.078,2975.453,	2994.969,
3014.453,3033.891,3053.375,	3072.859,3092.328,3111.828,	3131.234,3150.672,3170.187,	3189.734,
3209.297,3229.016,3248.625,	3268.328,3287.969,3307.625,	3327.234,3346.937,3366.781,	3386.641,
3406.656,3426.625,3446.875,	3466.984,3487.078,3507.328,	3527.562,3547.781,3568.062,	3588.328,
3608.594,3628.953,3649.375,	3669.875,3690.234]
cur=-[
-0.00490158920746269,	-0.00147800555164251,	-2.01252832408604,	-2.01397936209874,	-2.01114359599115,
-2.01300662781294,-2.01439965705264,-2.01160288530172,-2.01801468763330,-2.01313475351151,-2.01316220224595,
-2.01303023418461,-2.01371049851857,-2.01415013745466,-2.01346398324100,-2.01264583638754,-2.01340930058713,
-2.01314549982481,-2.01204128943031,-2.01446325299940,-2.01314002971595,-2.00922643619629,-2.01466618038391,
-2.01321832884032,-2.01151422011174,-2.01220227981274,-2.01385334477125,-2.01213146011162,-2.01070173140752,
-2.01208603752232,-2.01229391932774,-2.01346597836118,-2.01358796093004,-2.01284705723582,-2.01242628276702,
-2.01245565050621,-2.01338480049785,-2.01262885734551,-2.01550536193371,-2.01132965757290,-2.01087904153185,
-2.01299613028587,-2.01075794856787,-2.01370126925639,-2.01114443646185,-2.01217407576836,-2.01231592461865,
-2.01201237565114,-2.01208037472717,-2.01071108108784,-2.01334411961785,-2.01103073973055,-2.01268369370295,
-2.01144124181049,-2.01279047704850,-2.01156919175811,-2.01266319129602,-2.01390538807608,-2.01071577554447,
-2.01202154557025,-2.00921432744165,-2.01370984965265,-2.01281710553717,-2.01406265698721,-2.01107376465680,
-2.01089290844406,-2.01240394459316,-2.01323590772332,-2.01419384651964,-2.01263667547701,-2.01221031969746,
-2.01155223845763,-2.01406420553574,-2.01189939235907,-2.01541865942174,-2.01274965002258,-2.01385129522284,
-2.01355408336219,-2.01630990968773,-2.01277099624401,-2.01282271630919,-2.01251995347493,-2.01158502785775,
-2.01622851328796,-2.01222751351201,-2.01382425980576,-2.01388780465926,-2.01084308459760,-2.01067819804807,
-2.01347742984303,-2.01328656682240,-2.01174038327055,-2.01302220679742,-2.01341143203016,-2.01369895987642,
-2.01306414177861,-2.01105141466717,-2.01281524432994,-2.01484013680141,-2.01296867524278,-2.01375543740492,
-2.01199770761214,-2.01126505496849,-2.01018698097930,-2.01252046603508,-2.01359547876020,-2.01306955725642,
-2.01212855894781,-2.01165083343349,-2.01092615557540,-2.01508274597528,-2.01250901677532,-2.01533995414952,
-2.01178460382067,-2.01414060527905,-2.01195543862615,-2.01170081554690,-2.01229053985690,-2.01129035042103,
-2.01162470119891,-2.01361557545006,-2.01193090976318,-2.01518164543474,-2.01404367333819,-2.01321776062893,
-2.01166963894327,-2.01289334202513,-2.01194438315430,-2.01242748450159,-2.01029962998185,-2.01265850112307,
-2.01396821918466,-2.01445353476959,-2.01075963984222,-2.01217799840737,-2.00870934653250,-2.01667912857894,
-2.01385867085400,-2.01157581274602,-2.01111904367465,-2.01352387230768,-2.01129222159290,-2.01389014130306,
-2.01266874049439,-2.01284143757402,-2.01109682920084,-2.01176990089421,-2.01056085398907,-2.01244788793771,
-2.00991454145122,-2.01524819475563,-2.01509534741295,-2.01287592478801,-2.00971985804389,-2.01252560977573,
-2.01247177160775,-2.01239144498529,-2.01177114049426,-2.01305721063971,-2.01312977446971,-2.01050450679999,
-2.01344675453145,-2.01199821133029,-2.01407644805371,-2.01006854225820,-2.01386694513532,-2.01101273322598,
-2.01231137821328,-2.01109019781861,-2.01123506748447,-2.01171702739523,-2.01106145318504,-2.01168998518774,
-2.01339607190516,-2.01480817601455,-2.01109469569165,-2.01421697020321,-2.01060914921418,-2.01411159372913,
-2.01263909909732,-0.00420062513414368,	-0.00292880440280358,-0.000231200171749473,-0.00123136610093889,-0.00101438395288505,
-0.00165459735373620,	-0.00306905502911012,	-0.00155671726282023,	-0.00345773355377072,	 0.000231405806106993,
0.000728585272845447,	-0.00265918084987476,	-0.00151257983033383,	-0.00153169125712880,	-0.000405677789241947,
-0.000387595838485931,	-0.00652835125160699]
vol=[
4.19149180750530,4.19074906777610,3.97487091222999,	3.95171670751089,3.93435248943260,3.92005844166791,
3.90790350991309,3.89703572169109,3.88747658054044,	3.87895859247514,3.87101598480779,3.86347967749340,
3.85666050703146,3.85001357916834,3.84384190541381,	3.83768797571906,3.83193825298352,3.82641444300992,
3.82109573438256,3.81589666466863,3.81069682190636,	3.80583738096969,3.80112163346330,3.79629313055025,
3.79144890639008,3.78662468670613,3.78206292303341,	3.77717540186865,3.77240629397936,3.76793745919496,
3.76322930364035,3.75851149202913,3.75421321007525,	3.74962920289956,3.74516885963638,3.74100500876119,
3.73669982885741,3.73248270531015,3.72835255440004,	3.72453323547867,3.72065102608771,3.71640900822778,
3.71279183783938,3.70842761443677,3.70479177165460,	3.70068257779630,3.69675771477012,3.69303933221240,
3.68917679669778,3.68521448154197,3.68149374141573,	3.67779221378976,3.67406181318209,3.67030176499603,
3.66673391153023,3.66299804950626,3.65966438680144,	3.65589413612439,3.65232651300318,3.64876348229365,
3.64554685927165,3.64207625483732,3.63846406305335,	3.63507606938654,3.63166482425027,3.62857853646882,
3.62501259928352,3.62158222681123,3.61846125166182,	3.61506725193750,3.61188180875812,3.60857867726732,
3.60553494110992,3.60218589935641,3.59885714511492,	3.59576709379380,3.59301666356893,3.58971009399727,
3.58665220327661,3.58344203823537,3.58048534434714,	3.57725047813027,3.57436883391113,3.57128512768717,
3.56872066550615,3.56555467329979,3.56259540878202,	3.55978447041058,3.55689375383135,3.55397358754084,
3.55119782742195,3.54851559957747,3.54580583110077,	3.54311072155012,3.54029355759787,3.53740100559044,
3.53518704275422,3.53250894116367,3.52990333773161,	3.52754753455626,3.52501550598722,3.52260335738846,
3.52011277923015,3.51796894567226,3.51538624843045,	3.51328369043578,3.51113986209176,3.50880484091507,
3.50689114847090,3.50471159524183,3.50278636905226,	3.50064119243639,3.49884032472450,3.49651009173077,
3.49476208394121,3.49264904176421,3.49090381147511,	3.48884863255138,3.48686928955369,3.48499587233813,
3.48317381201077,3.48099635282812,3.47932882883719,	3.47721288139096,3.47531054138315,3.47344329038601,
3.47118790190817,3.46915803990222,3.46715799430632,	3.46509799303156,3.46319748544490,3.46076570511258,
3.45852954324684,3.45637594635014,3.45382105374173,	3.45159071222752,3.44934341523932,3.44691177545508,
3.44433708002687,3.44186104809679,3.43894026608122,	3.43636103495106,3.43367080904850,3.43089658037290,
3.42782089757862,3.42484405913391,3.42175726911613,	3.41874632404291,3.41545780861659,3.41177600871450,
3.40829396937228,3.40435779393306,3.40023624064031,	3.39613071787179,3.39181337842552,3.38704829198639,
3.38203092034683,3.37654809410201,3.37045168080920,	3.36378702895450,3.35685219381610,3.34867749480244,
3.33959089746460,3.32963512743319,3.31820790626283,	3.30494144277741,3.29037736566309,3.27371569393002,
3.25451646129461,3.23310264309898,3.20827869027007,	3.17996176471995,3.14722299323229,3.10933058143290,
3.06538157935910,3.01306122356690,2.94920549054709,	2.86601801412344,2.75725191265609,2.61246734790709,
2.99812519730241,3.07035107039746,3.11626935891479,	3.14865270656010,3.17290605528670,3.19102610813801,
3.20603613243278,3.21844940269493,3.22875341147208,	3.23809513428383,3.24584399920062,3.25261433963303,
3.25870451344492,3.26412058058012,3.26902995262555,	3.27320673224271,3.27716997682520]

struct MyLoss <: DiffEqBase.DECostFunction
t
data
pos::Integer
weight
diff_weight
end

MyLoss(t,data,pos;weight=nothing,diff_weight=nothing)=MyLoss(t,data,pos,weight,diff_weight)

function (f::MyLoss)(sol::DiffEqBase.DESolution)

if sol.retcode != :Success
return Inf
end

sumsq = 0.0
t=f.t
data=f.data
pos=f.pos
if f.weight==nothing
for i in 2:length(t)
sumsq += (sol(t[i])[f.pos] - data[i])^2
if f.diff_weight!=nothing
sumsq+=f.diff_weight*(data[i]-data[i-1]-sol(t[i])[pos]+sol(t[i-1])[pos])^2
end
end
else
for i in 2:length(t)
sumsq += f.weight*(sol(t[i])[pos] - data[i])^2
if f.diff_weight!=nothing
sumsq+=f.diff_weight*(data[i]-data[i-1]-sol(t[i])[pos]+sol(t[i-1])[pos])^2
end
end
end
return sumsq
end
capacity(t, cur) = abs(sum((t[i + 1] - t[i]) * cur[i] for i = 1:length(t) - 1))

function optim_cycle(t, cur, vol, p,CC)
cur_spl = SmoothingSplines.fit(SmoothingSpline, t, cur, 1.0)

function f(out, du, u, p, t)
oc1, oc2, oc3, oc4, oc5, oc6, rs1, rs2, rs3, rts1, rts2, rts3, cts1, cts2, cts3, rtl1, rtl2, rtl3, ctl1, ctl2, ctl3 = p

cur = SmoothingSplines.predict(cur_spl, t)
Voc = oc1 * exp(oc2 * u[1]) + oc3 + oc4 * u[1] + oc5 * u[1]^2 + oc6 * u[1]^3
R = rs1 * exp(rs2 * u[1]) + rs3
Rts = rts1 * exp(rts2 * u[1]) + rts3
Cts = cts1 * exp(cts2 * u[1]) + cts3
Rtl = rtl1 * exp(rtl2 * u[1]) + rtl3
Ctl = ctl1 * exp(ctl2 * u[1]) + ctl3

out[1] = du[1] + cur / CC
out[2] = du[2] - cur / Cts + u[2] / Rts / Cts
out[3] = du[3] - cur / Ctl + u[3] / Rtl / Ctl
out[4] = u[4] - (Voc - cur * R - u[2] - u[3])
end

du = [0.0,0.0,0.0,0.0]
u = [1.0,0.0,0.0,vol[1]]

loss = MyLoss(t,vol,4;weight=1.0,diff_weight=1.0)
tspan = (0.0, t[end])
differential_vars = [true,true,true,false]
prob = DAEProblem(f, du, u, tspan, p, differential_vars = differential_vars)
cost_function = build_loss_objective(prob,IDA(), loss,verbose = true,saveat=t)
result=optimize(cost_function, p)
prob = DAEProblem(f, du, u, tspan, result.minimizer, differential_vars = differential_vars)

sol = solve(prob,IDA(),saveat=t)
plot(sol.t,sol[4,:])
plot!(t,vol)
end

function findreverse(t, vol)
d = []
c = diff(vol)
for i = 1:length(c) - 1
if c[i] * c[i + 1] < 0
push!(d, i)
m=i
end
end
return d
end

CC=capacity(t,cur)

p = [-1.031,-35,3.685,0.2156,-0.1178,0.3201,0.1562,-24.37,0.07446,0.3208,-29.14,0.04669,-752.9,-13.5,703.6,6.603,-155.2,0.04984,-6056,-27.12,4475]
#d=findreverse(t,vol)
#t=t[1:d[1]-6]
#cur=cur[1:d[1]-6]
#vol=vol[1:d[1]-6]
optim_cycle(t,cur,vol,p,CC)
'''``````

Are there parameters that lead to unstable ODEs? In that case, you might want to reparameterize or constrain the optimization so that way those parameters can’t be chosen.

1 Like

uncomment the four sentence，the result is correct，maybe need special adaptive step control ，but I do not know how to do

if init parameter is neg，upbound should be 0，if init parameter is positive，low bound should be 0。it is li battery model

Use FMinBox for the optimizer to set the upper bound.

I ref this “Handling Instability When Solving ODE Problems - #5 by ChrisRackauckas”,Thanks.
But I found in solver option ’ * `d_discontinuities:` Denotes locations of discontinuities in low order derivatives. This will force FSAL algorithms which assume derivative continuity to re-evaluate the derivatives at the point of discontinuity. The default is `[]` .',can you give some example how to use this d_discontinuities

``````function condition(u,t,integrator)