From 8f3dd599ae4aecac2fc881cffbe0b671083796b9 Mon Sep 17 00:00:00 2001 From: ale Date: Thu, 13 Jun 2019 16:19:28 -0700 Subject: [PATCH 1/2] Create ODEwithFixParams.jl --- examples/ODEwithFixParams.jl | 118 +++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) create mode 100644 examples/ODEwithFixParams.jl diff --git a/examples/ODEwithFixParams.jl b/examples/ODEwithFixParams.jl new file mode 100644 index 00000000..dc7faa85 --- /dev/null +++ b/examples/ODEwithFixParams.jl @@ -0,0 +1,118 @@ +using Distributed +using Plots +addprocs(8) + + +@everywhere using DifferentialEquations +@everywhere using DiffEqParamEstim +@everywhere using BlackBoxOptim + +# Define function to fix parameters + +function fixvalues(vs, fixspec, dims) + newvs = zeros(Float64, dims) + vsidx = fixidx = 1 + for i in 1:dims + if (fixidx <= length(fixspec) && i == fixspec[fixidx][1]) + newvs[i] = fixspec[fixidx][2] + fixidx += 1 + else + newvs[i] = vs[vsidx] + vsidx += 1 + end + end + newvs +end + +# Define DAE System: Lorenz Attractor + Extra Variable to illustarte DAE syntax +function g(du,u,p,t) + σ,ρ,β = p + x,y,z,w = u + du[1] = σ*(y-x) + du[2] = x*(ρ-z) - y + du[3] = x*y - β*z + du[4] = w - (x + y + z) +end + +# Define Limits +lower=Float64.([0, 0, 0]) +upper=Float64.([20, 30, 10]) + +# Define Initial Conditions as well as time span +u0 = [1.0;0.0;0.0;1.0] +t = 0.0:0.05:1 +tspan = (0.0,1.0) + +# Generate Mass Matrox MM to transfrom ODE into DAE +MM=zeros(4,4) +MM[1,1]=MM[2,2]=MM[3,3]=1 + +# Generate Equation system +ODESystem = ODEFunction(g;mass_matrix=MM) + +# Generate Function that creates ODE problem as function of parameters (optmization argument) +model_ode(p_) = ODEProblem(ODESystem, u0, tspan,p_) + +# Generate Function that solves system as function of parameters, saves every s interval and save only vars we want. +solve_model(mp_,s, i) = DifferentialEquations.solve(model_ode(mp_), Rodas5(),saveat=s, save_idxs=i) + +# Generate Clean mock data / will add nisy mock data later +mock_data = Array(solve_model([10.0,28.0,8/3],0.05,[1,2,3,4]))#,4 + + +# Create Loss function +function L(t,dat,w=nothing) + return L2Loss(t,dat;data_weight=w) +end + +# Create objective function +loss_objective(mp_, dat)=build_loss_objective(model_ode(mp_), Rodas5(), L(t,dat)) + +# Simplify fitness/objective function to depend only on parameter we want to vary +function origfitness(args) + return loss_objective(args, mock_data)(args) +end + +dims = 3 + +# In specific opt problem A there is one var (at position 1) with fixed value of 10. + +problemfixed = [(1, 10.0)] + +# Create objective function with fix value =10 at position 1 +fitnessFix(vs) = origfitness(fixvalues(vs, problemfixed, dims)) + +# Generate options for optimization for problem with free parameters +opt_free = bbsetup(origfitness; Method=:separable_nes, SearchRange = (collect(zip(lower,upper))), + NumDimensions = 3, MaxFuncEvals = 100000, workers=workers(), TargetFitness=100.0 ) + +# Generate options for optimization for problem with fix parameters (note dimensions are 2 instead of 3) +opt_fix = bbsetup(fitnessFix; Method=:separable_nes, SearchRange = (collect(zip(lower[2:3],upper[2:3]))), + NumDimensions = 2, MaxFuncEvals = 100000, workers=workers()) + +# Call bboptim as usual and optimizing for fitnessFix and origfitness. +bbfix=bboptimize(opt_fix); +bbfree=bboptimize(opt_free); +bsfix=fixvalues(best_candidate(bbfix), problemfixed,dims) # then is our best solution to the general DAE problem with a fixed param. +bsfree=best_candidate(bbfree) # then is our best solution to the general DAE problem with full freedom. + +Note that within error bsfix and bs free are equal + +(bsfree, origfitness(bsfree)), (bsfix, origfitness(bsfix)) + + +# Plot results +pyplot() +# Shot Times to see Fitness to Data +Sol1=solve(ODEProblem(ODESystem, u0, (0,1.0), bsfix),Rodas5()) +plot(Sol1, layout=(2,2), xlabel="Time", ylabel=[:x :y :z :w], label=["" "" "" "w=x+y+z"]) +scatter!(t,[mock_data[i,:] for i in 1:4],layout=(2,2), label="") + +# Long times to see attractor behavior +Sol2=solve(ODEProblem(ODESystem, u0, (0,100.0), bsfix),Rodas5()) +#plot(Sol,layout=(2,2), labels=[:w :x :y "z=w+x+y"]) +P1=plot(Sol2.t, [Sol2[i,:] for i in 1:4], layout=(2,2), xlabel="Time", ylabel=[:x :y :z :w], labels=["" "" "" "w=x+y+z"]) +P2=plot(Sol2,vars=(1,2,3),xlabel="x", ylabel="y", zlabel="z", labels="") +l = @layout([a;b{0.5h}]) +PF=plot(P1,P2, layout=l, size=(900,600)) + From 6f70133c516c147bb94b056047ff7b41270a112f Mon Sep 17 00:00:00 2001 From: ale Date: Thu, 13 Jun 2019 16:25:02 -0700 Subject: [PATCH 2/2] Update ODEwithFixParams.jl commented line --- examples/ODEwithFixParams.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ODEwithFixParams.jl b/examples/ODEwithFixParams.jl index dc7faa85..20a3b6f7 100644 --- a/examples/ODEwithFixParams.jl +++ b/examples/ODEwithFixParams.jl @@ -96,7 +96,7 @@ bbfree=bboptimize(opt_free); bsfix=fixvalues(best_candidate(bbfix), problemfixed,dims) # then is our best solution to the general DAE problem with a fixed param. bsfree=best_candidate(bbfree) # then is our best solution to the general DAE problem with full freedom. -Note that within error bsfix and bs free are equal +# Note that within error bsfix and bs free are equal (bsfree, origfitness(bsfree)), (bsfix, origfitness(bsfix))