Optimization problem inside `ODESystem`

Specific Domains Modelling & Simulations
1

Hello,
I am trying to solve an optimization problem with the variables of an ODESystem. I have a toy code below:

using ModelingToolkit, DifferentialEquations, Optimization, OptimizationOptimJL,ForwardDiff, Optim
using ModelingToolkit: t_nounits as t, D_nounits as D

function SolveAlpha(x0,p)
    f(x,p) =  1  
    f_opt =  OptimizationFunction(f,AutoForwardDiff())
    prob = OptimizationProblem(f_opt,x0,p)
    sol = solve(prob,Optim.BFGS())
    return sol.u
end
@register_symbolic SolveAlpha(x0,p)

@component function system(;name)
    vars =@variables begin
        x(t),
        y(t),
        α(t)
    end
    para = @parameters begin
        x0 = [1.0]
    end
    eqs = [
        D(x) ~ -y
        D(y) ~ x
        α ~ SolveAlpha(x0,[x,y])
    ]
    System(eqs, t, vars, para;name=name)
end

@named mysys = system()
sys = structural_simplify(mysys)
u0 = [sys.x => 1.0, sys.y => 0.0]
tspan = (0,10)
para = []
prob = ODEProblem(sys,u0,tspan,para)
sol = solve(prob)

The error i get is the following:

ERROR: LoadError: MethodError: no method matching fill!(::Num, ::Num)
The function `fill!` exists, but no method is defined for this combination of argument types

This is internal to Optimization.jl → NLSolversBase.jl.
So my question would be is it possible to pass the local variables of a @component to optimizers? Or is there something wrong with the way iI have used @register_symbolic which fails here?

Thank you

2

You are getting into problems with packing your variables in vectors. Here is the working solution.

using ModelingToolkit, DifferentialEquations, Optimization, OptimizationOptimJL,ForwardDiff, Optim
using ModelingToolkit: t_nounits as t, D_nounits as D

function solve_alpha(x0,x,y)
    x0 = [x0]
    p = [x, y]
    f(x,p) =  1  
    f_opt =  OptimizationFunction(f,AutoForwardDiff())
    prob = OptimizationProblem(f_opt,x0,p)
    sol = solve(prob,Optim.BFGS())
    return sol.u
end
@register_symbolic solve_alpha(x0,x,y)

@component function system(;name)
    vars =@variables begin
        x(t),
        y(t),
        α(t)
    end
    para = @parameters begin
        x0 = 1.0
    end
    eqs = [
        D(x) ~ -y
        D(y) ~ x
        α ~ solve_alpha(x0,x,y)
    ]
    System(eqs, t, vars, para;name=name)
end

@named mysys = system()
sys = structural_simplify(mysys)
u0 = [sys.x => 1.0, sys.y => 0.0]
tspan = (0,10)
para = []
prob = ODEProblem(sys,u0,tspan,para)
sol = solve(prob)
2 Likes
3

Thanks!!

1 Like