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

I don’t understand your response.

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”,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

add these code

function condition(u,t,integrator)
        t-t_discontinue
    end
    
    function affect!(integrator)
        println("u1=",integrator.u[1],"t=",integrator.t)
    end

    cb = ContinuousCallback(condition,affect!,save_positions=(false,false))
'''