# Parameter optimization using Multiple shooting

**URL:** <https://discourse.julialang.org/t/parameter-optimization-using-multiple-shooting/104104>\
**Category:** Machine Learning\
**Created:** [September 21, 2023, 12:16pm UTC](https://discourse.julialang.org/t/parameter-optimization-using-multiple-shooting/104104 "2023-09-21T12:16:07Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marco\_Nesta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_nesta/32/50318_2.png) [@Marco\_Nesta](https://discourse.julialang.org/u/Marco_Nesta)\
**Post date:** [September 21, 2023, 12:16pm UTC](https://discourse.julialang.org/t/parameter-optimization-using-multiple-shooting/104104/1 "2023-09-21T12:16:07Z")

</div>

> **[Multiple Shooting · DiffEqFlux.jl](https://docs.sciml.ai/DiffEqFlux/stable/examples/multiple_shooting/)**
>
> Documentation for DiffEqFlux.jl.

'm using the same structure as in the link, but how can I view the updated parameters?

this is my code:

```julia-auto

```

using DifferentialEquations, Plots, Flux,Optim, DiffEqFlux, DataInterpolations,Random, ComponentArrays, Lux  
using Optimization, OptimizationOptimisers, OptimizationOptimJL,OptimizationNLopt  
rng = Random.default\_rng()  
using CSV  
using DataFrames  
using Plots  
using Flux  
using Statistics: mean  
using DiffEqFlux, Optimization, OptimizationOptimJL,Plots  
using ComponentArrays, Lux, DiffEqFlux, Optimization, OptimizationPolyalgorithms, DifferentialEquations, Plots  
using DiffEqFlux: group\_ranges

# Load the data

df = CSV.read(“/Users/marco/OneDrive/Desktop/Tesi/Julia\_files/measured\_data\_ex\_00.csv”, DataFrame)  
df=repeat(df, outer=10)

registered\_gas\_flow = df[:, :2]  
Registered\_Temperature = df[:, 3]  
sec\_Temp = [mean(Registered\_Temperature[i:i+59]) for i in 1:60:length(Registered\_Temperature)-59]  
sec\_gas = [mean(registered\_gas\_flow[i:i+59]) for i in 1:60:length(registered\_gas\_flow)-59]  
#create a 3600 time vector  
function RC!(du,u,p,t)  
Rv, Ci,Rv2,Ci2 = p  
Text = ext\_Temp(t)  
P= ext\_power(t)  
P2= ext\_power2(t)

```
du[1] = 1/(Rv*Ci) .* (Text .- u[1]) .+ P/Ci
end

```

u0= [20.0]

tspan= (0.0f0,3000.0f0)  
datasize = 3000  
tsteps= range(tspan[1], tspan[2], length = datasize)  
p= [10000, 10, 10000, 10]

Ext\_temperature = LinearInterpolation(sec\_Temp,tsteps);  
power = LinearInterpolation(sec\_gas,tsteps);  
power2= LinearInterpolation(sec\_gas,tsteps);  
function ext\_Temp(tsteps)  
return Ext\_temperature(tsteps)  
end

function ext\_power(tsteps)  
return power(tsteps)  
end

function ext\_power2(tsteps)  
return power2(tsteps)  
end

prob= ODEProblem(RC!, u0, tspan, p)  
sol= solve(prob, Tsit5(), saveat=tsteps)[1,:]

# Verify ODE solution

ode\_data =Array(solve(prob, Tsit5(), saveat=tsteps))

plot(tsteps, ode\_data[1,:], label = “data”)

#resize ode data  
ode\_data=ode\_data[:,1:20:3000]  
tsteps=tsteps[1:20:3000]

plot(tsteps, ode\_data[1,:], label = “data”)

ode\_data\_train=ode\_data[:,1:75]  
ode\_data\_test=ode\_data[:,76:150]  
tspan= (0.0f0, 75.0f0)  
tsteps= range(tspan[1], tspan[2], length = 75)  
datasize=75

plot(ode\_data[1,1:75], label = “datatrain”)  
plot!(ode\_data[1,76:150], label = “datatest”)

anim = Plots.Animation()

# Define the Neural Network

nn = Lux.Chain(  
Lux.Dense(1, 84, tanh),  
Lux.Dense(84, 84, tanh),  
Lux.Dense(84, 44, tanh),  
Lux.Dense(44,1))  
p\_init, st = Lux.setup(rng, nn)

neuralode = NeuralODE(nn, tspan, Tsit5(), saveat = tsteps)  
prob\_node = ODEProblem((u,p,t)-\>nn(u,p,st)[1], u0, tspan, ComponentArray(p\_init))

function plot\_multiple\_shoot(plt, preds, group\_size)  
step = group\_size-1  
ranges = group\_ranges(datasize, group\_size)

```
for (i, rg) in enumerate(ranges)
	plot!(plt, tsteps[rg], preds[i][1,:], markershape=:circle)
end

```

end

# Animate training, cannot make animation on CI server

# anim = Plots.Animation()

iter = 0  
callback = function (p, l, preds; doplot = true)  
display(l)  
global iter  
iter += 1  
if doplot && iter%1 == 0  
# plot the original data  
plt = scatter(tsteps, ode\_data\_train[1,:], label = “Data”)

```
# plot the different predictions for individual shoot
plot_multiple_shoot(plt, preds, group_size)

frame(anim)
display(plot(plt))

```

end  
return false  
end

# Define parameters for Multiple Shooting

group\_size = 3  
continuity\_term = 100

function loss\_function(data, pred)  
return sum(abs2, data - pred)  
end

function loss\_multiple\_shooting(p)  
return multiple\_shoot(p, ode\_data\_train, tsteps, prob\_node, loss\_function, Tsit5(),  
group\_size; continuity\_term)  
end

adtype = Optimization.AutoZygote()  
optf = Optimization.OptimizationFunction((x,p) → loss\_multiple\_shooting(x), adtype)  
optprob = Optimization.OptimizationProblem(optf, ComponentArray(p\_init))  
res\_ms = Optimization.solve(optprob, PolyOpt(), callback = callback)

gif(anim, “multiple\_shooting.gif”, fps=15)

---

<div class="post-metadata">

**Author:** ![Marco\_Nesta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marco_nesta/32/50318_2.png) [@Marco\_Nesta](https://discourse.julialang.org/u/Marco_Nesta)\
**Post date:** [September 21, 2023, 12:18pm UTC](https://discourse.julialang.org/t/parameter-optimization-using-multiple-shooting/104104/2 "2023-09-21T12:18:19Z")

</div>

[https://docs.sciml.ai/DiffEqFlux/stable/examples/multiple\_shooting/](https://docs.sciml.ai/DiffEqFlux/stable/examples/multiple_shooting/)  
I’m using the same structure as in the link, but how can I view the updated parameters?

```julia
using DifferentialEquations, Plots, Flux,Optim, DiffEqFlux, DataInterpolations,Random, ComponentArrays, Lux
using Optimization, OptimizationOptimisers, OptimizationOptimJL,OptimizationNLopt
rng = Random.default_rng()
using CSV
using DataFrames
using Plots
using Flux
using Statistics: mean
using DiffEqFlux, Optimization, OptimizationOptimJL,Plots
using ComponentArrays, Lux, DiffEqFlux, Optimization, OptimizationPolyalgorithms, DifferentialEquations, Plots
using DiffEqFlux: group_ranges
# Load the data
df = CSV.read("/Users/marco/OneDrive/Desktop/Tesi/Julia_files/measured_data_ex_00.csv", DataFrame)
df=repeat(df, outer=10)

registered_gas_flow = df[:, :2]
Registered_Temperature = df[:, 3]
sec_Temp = [mean(Registered_Temperature[i:i+59]) for i in 1:60:length(Registered_Temperature)-59]
sec_gas = [mean(registered_gas_flow[i:i+59]) for i in 1:60:length(registered_gas_flow)-59]
#create a 3600 time vector
function RC!(du,u,p,t)
    Rv, Ci,Rv2,Ci2 = p
    Text = ext_Temp(t)
    P= ext_power(t)
    P2= ext_power2(t)
    
    du[1] = 1/(Rv*Ci) .* (Text .- u[1]) .+ P/Ci
    end

u0= [20.0]

tspan= (0.0f0,3000.0f0)
datasize = 3000
tsteps= range(tspan[1], tspan[2], length = datasize)
p= [10000, 10, 10000, 10]

Ext_temperature = LinearInterpolation(sec_Temp,tsteps);
power = LinearInterpolation(sec_gas,tsteps);
power2= LinearInterpolation(sec_gas,tsteps);
function ext_Temp(tsteps)
return Ext_temperature(tsteps)
end
    
function ext_power(tsteps)
return power(tsteps)
end   

function ext_power2(tsteps)
return power2(tsteps)
end

prob= ODEProblem(RC!, u0, tspan, p)
sol= solve(prob, Tsit5(), saveat=tsteps)[1,:]

# Verify ODE solution
ode_data =Array(solve(prob, Tsit5(), saveat=tsteps))

plot(tsteps, ode_data[1,:], label = "data")

#resize ode data 
ode_data=ode_data[:,1:20:3000]
tsteps=tsteps[1:20:3000]

plot(tsteps, ode_data[1,:], label = "data")

ode_data_train=ode_data[:,1:75]
ode_data_test=ode_data[:,76:150]
tspan= (0.0f0, 75.0f0)
tsteps= range(tspan[1], tspan[2], length = 75)
datasize=75

plot(ode_data[1,1:75], label = "datatrain")
plot!(ode_data[1,76:150], label = "datatest")

anim = Plots.Animation()

# Define the Neural Network
nn = Lux.Chain(
                Lux.Dense(1, 84, tanh),
                Lux.Dense(84, 84, tanh),
                Lux.Dense(84, 44, tanh),
                Lux.Dense(44,1))
p_init, st = Lux.setup(rng, nn)

neuralode = NeuralODE(nn, tspan, Tsit5(), saveat = tsteps)
prob_node = ODEProblem((u,p,t)->nn(u,p,st)[1], u0, tspan, ComponentArray(p_init))

function plot_multiple_shoot(plt, preds, group_size)
	step = group_size-1
	ranges = group_ranges(datasize, group_size)

	for (i, rg) in enumerate(ranges)
		plot!(plt, tsteps[rg], preds[i][1,:], markershape=:circle)
	end
end

# Animate training, cannot make animation on CI server
# anim = Plots.Animation()
iter = 0
callback = function (p, l, preds; doplot = true)
  display(l)
  global iter
  iter += 1
  if doplot && iter%1 == 0
    # plot the original data
    plt = scatter(tsteps, ode_data_train[1,:], label = "Data")

    # plot the different predictions for individual shoot
    plot_multiple_shoot(plt, preds, group_size)

    frame(anim)
    display(plot(plt))
  end
  return false
end

# Define parameters for Multiple Shooting
group_size = 3
continuity_term = 100

function loss_function(data, pred)
	return sum(abs2, data - pred)
end

function loss_multiple_shooting(p)
    return multiple_shoot(p, ode_data_train, tsteps, prob_node, loss_function, Tsit5(),
                          group_size; continuity_term)
end

adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x,p) -> loss_multiple_shooting(x), adtype)
optprob = Optimization.OptimizationProblem(optf, ComponentArray(p_init))
res_ms = Optimization.solve(optprob, PolyOpt(), callback = callback)

gif(anim, "multiple_shooting.gif", fps=15)

```