Parameter Estimate Problem

Specific Domains Optimization (Mathematical)
1

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)
'''
2

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
3

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

4

I don’t understand your response.

5

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

6

Use FMinBox for the optimizer to set the upper bound.

7

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

image
image699×460 28.5 KB

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))
'''