OptimizationSystem from a set of multiple equations

I have a set of equations where there are some input variables and one output variable. I want to minimize the output variable, and thought to do this with an OptimizationProblem. In this example, distance is the only input, and I want to minimize total_acc. My question is: how can I find the equation of total_acc relative to distance programmatically, so I can pass it to the OptimizationSystem?

# Falling mass, attached to linear spring
using ModelingToolkit, Optimization, OptimizationOptimJL, LinearAlgebra
using ModelingToolkit: t_nounits as t, D_nounits as D

G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]

# defing the model, Z component upwards
@parameters mass=1.0 c_spring=50.0 damping=0.5 l0=L0
@variables   pos(t)[1:3]
@variables   vel(t)[1:3]
@variables   acc(t)[1:3]
@variables unit_vector(t)[1:3]
@variables spring_force(t)[1:3]
@variables norm1(t) spring_vel(t)
@variables total_acc(t)
@variables distance(t)

eqs = [
    pos          ~ distance * [0.0, 0.0, 1.0]
    vel          ~ [0.0, 0.0, 0.0]
    norm1        ~ norm(pos)
    unit_vector  ~ -pos / norm1         # direction from point mass to origin
    spring_vel   ~ -unit_vector ⋅ vel
    spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector
    acc          ~ G_EARTH + spring_force/mass
    total_acc ~ norm(acc)
]

eqs = Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs)))

@mtkbuild sys = OptimizationSystem(total_acc, [distance], [])
prob = OptimizationProblem(sys, [1.0], [], grad = true, hess = true)
u_opt = solve(prob, GradientDescent())
# solution should be the distance where the acceleration is the smallest

I’m trying to figure out what you’re trying to do. Are you trying to minimize total_acc as a function of distance? Distance is not defined.

Yes, I am trying to minimize total_acc as a function of distance. I forgot to define distance but changed it now… The variables are time-dependent because I want to use a very similar model in an ODE and I want to re-use the code, but the time-dependence is not used in this minimization problem.

So basically I want to programmatically be able to rewrite the system of equations eqs to one big equation, with total_acc on the lhs.

new_eq = total_acc ~ ...
@mtkbuild sys = OptimizationSystem(new_eq.rhs, [distance], [mass, c_spring, damping, l0])

@ChrisRackauckas I cant find any existing functions in Symbolics.jl that can do this. I will try to write my own function for this, that takes in a set of equations and tries to combine it into one big equation.

This works.

# Example one: Falling mass.
using Groebner, ModelingToolkit, Optimization, OptimizationOptimJL, LinearAlgebra
using ModelingToolkit: t_nounits as t, D_nounits as D

G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]

# defing the model, Z component upwards
@parameters mass=1.0 c_spring=50.0 damping=0.5 l0=L0
@variables   pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  0.0] 
@variables   acc(t)[1:3]
@variables unit_vector(t)[1:3]
@variables spring_force(t)[1:3]
@variables norm1(t) spring_vel(t)
@variables total_acc(t)
@variables distance(t)

eqs = [
    pos          ~ distance * [0.0, 0.0, 1.0]
    vel          ~ [0.0, 0.0, 0.0]
    norm1        ~ norm(pos)
    unit_vector  ~ -pos / norm1         # direction from point mass to origin
    spring_vel   ~ -unit_vector ⋅ vel
    spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector
    acc          ~ G_EARTH + spring_force/mass
    total_acc    ~ norm(acc)
]
eqs = Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs)))

function merge_eqs(eqs, opt_var)
    subs = Dict{Num, Num}()
    eqs = Vector{Any}(eqs)
    opt_i = 1
    for i in eachindex(eqs)
        if occursin(eqs[i].lhs, opt_var)
            opt_i = i
        end
        subs[eqs[i].lhs] = eqs[i].rhs
        eqs[i] = eqs[i].rhs
    end
    eqs[opt_i] = Symbolics.fixpoint_sub(eqs[opt_i], subs)
    return eqs[opt_i]
end

opt_eq = merge_eqs(eqs, total_acc)

@mtkbuild sys = OptimizationSystem(opt_eq, [distance], [mass, c_spring, l0])
prob = OptimizationProblem(sys, [-0.9], [mass => 1.0, c_spring => 50.0, l0 => L0], grad = true, hess = true)
@time sol = solve(prob, Optim.NelderMead(); maxiters=10_000)
1 Like

Just want to add that this is not a good solution, as it is very slow for bigger systems.

@ChrisRackauckas How difficult would it be to make OptimizationProblem compatible with a system of multiple equations?

Not so hard. We need to finalize the multi objective optimization interface in Optimization.jl though. But your case is only a single objective? You would want most of those equations as constraints instead

1 Like

Thanks Chris, using all equations as constraints was the solution. I somehow didn’t understand the constraints correctly… Anyways: here is the working code:

# Example one: Falling mass.
using Groebner, ModelingToolkit, Optimization, OptimizationOptimJL, LinearAlgebra
using ModelingToolkit: t_nounits as t, D_nounits as D
using ProfileView

G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]

# defing the model, Z component upwards
@parameters mass=1.0 c_spring=50.0 damping=0.5 l0=L0
@variables pos(t)[1:3] = [0.0, 0.0,  L0]
@variables vel(t)[1:3] = [0.0, 0.0,  0.0] 
@variables acc(t)[1:3]
@variables unit_vector(t)[1:3]
@variables spring_force(t)[1:3]
@variables norm1(t) spring_vel(t)
@variables total_acc(t)
@variables distance(t)

eqs = [
    pos          ~ distance * [0.0, 0.0, 1.0]
    vel          ~ [0.0, 0.0, 0.0]
    norm1        ~ norm(pos)
    unit_vector  ~ -pos / norm1         # direction from point mass to origin
    spring_vel   ~ -unit_vector ⋅ vel
    spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector
    acc          ~ G_EARTH + spring_force/mass
    total_acc    ~ norm(acc)
]

eqs = Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs)))

@mtkbuild sys = OptimizationSystem(total_acc, [distance], [mass, c_spring, l0]; constraints=eqs)
prob = OptimizationProblem(sys, [L0], [mass => 1.0, c_spring => 50.0, l0 => L0], grad = true, hess = true)
@time sol = solve(prob, Optim.NelderMead(); maxiters=10_000)

I tried to do this with a bigger system now. But the following code causes an UndefVarError, and I am not sure why.

Stacktrace of running Tether_11.jl:

ERROR: LoadError: UndefVarError: `total_force` not defined
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:430 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/Symbolics/YbNrd/src/build_function.jl:546 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:387 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ ./none:0 [inlined]
  [6] generated_callfunc
    @ ./none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::MTKParameters{…})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [8] (::OptimizationBase.var"#198#209"{OptimizationFunction{…}, MTKParameters{…}})(res::Vector{Float64}, x::Vector{Float64})
    @ OptimizationBase ~/.julia/packages/OptimizationBase/gvXsf/src/function.jl:180
  [9] __solve(cache::OptimizationCache{…})
    @ Optimization ~/.julia/packages/Optimization/cfp9i/src/lbfgsb.jl:116
 [10] solve!
    @ ~/.julia/packages/SciMLBase/ZyZAV/src/solve.jl:186 [inlined]
 [11] #solve#684
    @ ~/.julia/packages/SciMLBase/ZyZAV/src/solve.jl:94 [inlined]
 [12] macro expansion
    @ ./timing.jl:279 [inlined]
 [13] calc_initial_state(se::Settings3, ss::ODESystem, eqs::Vector{Equation}, p1::Vector{Int64}, elevation::Float64, tether_length::Float64)
    @ Main ~/Code/Tethers.jl/src/Tether_11.jl:104
 [14] main(; p1::Vector{Int64}, elevation::Float64, tether_length::Float64, fix_p1::Bool, fix_p2::Bool)
    @ Main ~/Code/Tethers.jl/src/Tether_11.jl:263
 [15] top-level scope
    @ ~/Code/Tethers.jl/src/Tether_11.jl:273
 [16] include(fname::String)
    @ Base.MainInclude ./client.jl:494
 [17] top-level scope
    @ REPL[17]:1
in expression starting at /home/bart/Code/Tethers.jl/src/Tether_11.jl:273
Some type information was truncated. Use `show(err)` to see complete types.

Maybe an issue with observed in the objective function codegen? Can you try to isolate that to an MWE? CC @cryptic.ax

Yeah that looks like a codegen issue. An MWE is greatly appreciated.

I tried to minimize the example as much as possible, without the error disappearing:

MWE with bug
using ModelingToolkit, LinearAlgebra, Optimization, OptimizationOptimJL
using ModelingToolkit: t_nounits as t, D_nounits as D

v_wind_tether::Vector{Float64} = [10.0, 0.0, 0.0]
segments::Int64 = 2                          # number of tether segments         [-]

function model()
    @variables begin 
        pos(t)[1:3, 1:segments+1]
        vel(t)[1:3, 1:segments+1]
        acc(t)[1:3, 1:segments+1]
        segment(t)[1:3, 1:segments]
        unit_vector(t)[1:3, 1:segments]
        len(t)[1:segments]
        v_apparent(t)[1:3, 1:segments]
        v_app_perp(t)[1:3, 1:segments]
        norm_v_app(t)[1:segments]
        half_drag_force(t)[1:3, 1:segments]
        total_force(t)[1:3, 1:segments+1]
    end

    # basic differential equations
    eqs = vcat(D.(pos) .~ vel,
                D.(vel) .~ acc)
    eqs = vcat(eqs...)
    # loop over all segments to calculate the spring and drag forces
    for i in 1:segments
        eqs = [
            eqs    
            segment[:, i]      ~ pos[:, i+1] - pos[:, i]
            len[i]             ~ norm(segment[:, i])
            unit_vector[:, i]  ~ -segment[:, i]/len[i]
            v_apparent[:, i]   ~ v_wind_tether .- (vel[:, i] + vel[:, i+1])/2
            v_app_perp[:, i]   ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i]
            norm_v_app[i]      ~ norm(v_app_perp[:, i])
            half_drag_force[:, i] ~ norm_v_app[i] * v_app_perp[:, i]
            ]
    end
    # loop over all tether particles to apply the forces and calculate the accelerations
    for i in 1:(segments+1)
        if i == segments+1
            push!(eqs, total_force[:, i] ~ half_drag_force[:, i-1])
        elseif i == 1
            push!(eqs, total_force[:, i] ~ half_drag_force[:, i])
        else
            push!(eqs, total_force[:, i] ~ half_drag_force[:, i-1])
        end
        push!(eqs, acc[:, i]         ~ total_force[:, i])
    end
    eqs = reduce(vcat, Symbolics.scalarize.(eqs))

    @mtkbuild sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs)), t)
    return sys, eqs
end

function calc_initial_state(ss, eqs, tether_length)
    # remove the differential equations
    diff_eqs = equations(ss)
    rm_idxs = []
    for i in eachindex(eqs)
        for d in diff_eqs
            if isequal(d.lhs, eqs[i].lhs)
                push!(rm_idxs, i)
                break
            end
        end
    end
    deleteat!(eqs, rm_idxs)

    # draw tethers
    @variables begin
        total_acc
        kite_distance
    end
    eqs = [
        eqs
        ss.pos[:, end] ~ kite_distance * [1, 0, 0]
        ss.vel[:, end] ~ [0, 0, 0]
        total_acc ~ norm(ss.acc)
    ]
    for i in 1:segments
        eqs = [
            eqs
            ss.pos[:, i] ~ [0, 0, 0]
            ss.vel[:, i] ~ zeros(3)
        ]
    end
    eqs = reduce(vcat, Symbolics.scalarize.(eqs))
    @mtkbuild sys = OptimizationSystem(total_acc, [kite_distance], []; constraints=eqs) simplify=false
    
    prob = OptimizationProblem(sys, [tether_length], []; grad = true, hess = true)
    sol = solve(prob, Optimization.LBFGS(); maxiters=1000)
    return sol
end


function main(; tether_length=10.0)
    simple_sys, eqs = model()
    sol = calc_initial_state(simple_sys, eqs, tether_length)
    nothing
end

main();

nothing

Making the MWE any more minimal than this, causes the error to disappear. For example:

MWE without bug, removing line 31-37 and changing total_force equations.
using ModelingToolkit, LinearAlgebra, Optimization, OptimizationOptimJL
using ModelingToolkit: t_nounits as t, D_nounits as D

v_wind_tether::Vector{Float64} = [10.0, 0.0, 0.0]
segments::Int64 = 2                          # number of tether segments         [-]

function model()
    @variables begin 
        pos(t)[1:3, 1:segments+1]
        vel(t)[1:3, 1:segments+1]
        acc(t)[1:3, 1:segments+1]
        segment(t)[1:3, 1:segments]
        unit_vector(t)[1:3, 1:segments]
        len(t)[1:segments]
        v_apparent(t)[1:3, 1:segments]
        v_app_perp(t)[1:3, 1:segments]
        norm_v_app(t)[1:segments]
        half_drag_force(t)[1:3, 1:segments]
        total_force(t)[1:3, 1:segments+1]
    end

    # basic differential equations
    eqs = vcat(D.(pos) .~ vel,
                D.(vel) .~ acc)
    eqs = vcat(eqs...)
    # loop over all segments to calculate the spring and drag forces
    for i in 1:segments
        eqs = [
            eqs    
            segment[:, i]      ~ pos[:, i+1] - pos[:, i]
            # len[i]             ~ norm(segment[:, i])
            # unit_vector[:, i]  ~ -segment[:, i]/len[i]
            # v_apparent[:, i]   ~ v_wind_tether .- (vel[:, i] + vel[:, i+1])/2
            # v_app_perp[:, i]   ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i]
            # norm_v_app[i]      ~ norm(v_app_perp[:, i])
            # half_drag_force[:, i] ~ norm_v_app[i] * v_app_perp[:, i]
            ]
    end
    # loop over all tether particles to apply the forces and calculate the accelerations
    for i in 1:(segments+1)
        if i == segments+1
            push!(eqs, total_force[:, i] ~ segment[:, i-1])
        elseif i == 1
            push!(eqs, total_force[:, i] ~ segment[:, i])
        else
            push!(eqs, total_force[:, i] ~ segment[:, i-1])
        end
        push!(eqs, acc[:, i]         ~ total_force[:, i])
    end
    eqs = reduce(vcat, Symbolics.scalarize.(eqs))

    @mtkbuild sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs)), t)
    return sys, eqs
end

function calc_initial_state(ss, eqs, tether_length)
    # remove the differential equations
    diff_eqs = equations(ss)
    rm_idxs = []
    for i in eachindex(eqs)
        for d in diff_eqs
            if isequal(d.lhs, eqs[i].lhs)
                push!(rm_idxs, i)
                break
            end
        end
    end
    deleteat!(eqs, rm_idxs)

    # draw tethers
    @variables begin
        total_acc
        kite_distance
    end
    eqs = [
        eqs
        ss.pos[:, end] ~ kite_distance * [1, 0, 0]
        ss.vel[:, end] ~ [0, 0, 0]
        total_acc ~ norm(ss.acc)
    ]
    for i in 1:segments
        eqs = [
            eqs
            ss.pos[:, i] ~ [0, 0, 0]
            ss.vel[:, i] ~ zeros(3)
        ]
    end
    eqs = reduce(vcat, Symbolics.scalarize.(eqs))
    @mtkbuild sys = OptimizationSystem(total_acc, [kite_distance], []; constraints=eqs) simplify=false
    
    prob = OptimizationProblem(sys, [tether_length], []; grad = true, hess = true)
    sol = solve(prob, Optimization.LBFGS(); maxiters=1000)
    return sol
end


function main(; tether_length=10.0)
    simple_sys, eqs = model()
    sol = calc_initial_state(simple_sys, eqs, tether_length)
    nothing
end

main();

nothing

The reason for this is a bit weird. The constraints of the OptimizationSystem involve all the variables of the ODESystem. However, the unknowns of the OptimizationSystem only contains tether_length, which results in an incorrect system and it bugs out. The solution (right now) is to do:

 @mtkbuild sys = OptimizationSystem(total_acc, [unknowns(ss); getproperty.(observed(ss), (:lhs,))], []; constraints=eqs) simplify = false

In calc_initial_state. The better solution would be for OptimizationSystem to have a constructor which doesn’t require specifying unknowns and parameters, and instead discovers them from the objective and constraints.

feat: add automatic variable discovery for `OptimizationSystem` by AayushSabharwal · Pull Request #3243 · SciML/ModelingToolkit.jl · GitHub adds the required functionality.

1 Like

Thanks!