# OptimizationSystem from a set of multiple equations
**URL:** https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809
**Category:** Modelling & Simulations
**Created:** [November 19, 2024, 12:09pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809 "2024-11-19T12:09:17Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 19, 2024, 12:09pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/1 "2024-11-19T12:09:17Z")
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?
```julia
# 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
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 19, 2024, 6:37pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/2 "2024-11-19T18:37:20Z")
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.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 20, 2024, 5:09am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/3 "2024-11-20T05:09:52Z")
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.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 20, 2024, 7:49am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/4 "2024-11-20T07:49:48Z")
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.
```julia
new_eq = total_acc ~ ...
@mtkbuild sys = OptimizationSystem(new_eq.rhs, [distance], [mass, c_spring, damping, l0])
```
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 22, 2024, 6:05am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/5 "2024-11-22T06:05:09Z")
@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.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 22, 2024, 7:08am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/6 "2024-11-22T07:08:36Z")
This works.
```julia
# 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)
```
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 22, 2024, 3:52pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/7 "2024-11-22T15:52:40Z")
Just want to add that this is not a good solution, as it is very slow for bigger systems.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 22, 2024, 4:32pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/8 "2024-11-22T16:32:20Z")
@ChrisRackauckas How difficult would it be to make `OptimizationProblem` compatible with a system of multiple equations?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 22, 2024, 8:23pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/9 "2024-11-22T20:23:24Z")
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
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 23, 2024, 8:23am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/10 "2024-11-23T08:23:22Z")
Thanks Chris, using all equations as constraints was the solution. I somehow didn’t understand the constraints correctly… Anyways: here is the working code:
```julia
# 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)
```
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 25, 2024, 4:06pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/11 "2024-11-25T16:06:58Z")
> [@ChrisRackauckas](#):
>
> You would want most of those equations as constraints instead
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:
```julia
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.
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 26, 2024, 12:39pm UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/12 "2024-11-26T12:39:36Z")
Maybe an issue with observed in the objective function codegen? Can you try to isolate that to an MWE? CC @cryptic.ax
---
### Author:  [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)
#### Post date: [November 27, 2024, 5:55am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/13 "2024-11-27T05:55:18Z")
Yeah that looks like a codegen issue. An MWE is greatly appreciated.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 27, 2024, 10:28am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/14 "2024-11-27T10:28:49Z")
I tried to minimize the example as much as possible, without the error disappearing:
> **MWE with bug**
>
> ```julia
> 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.**
>
> ```julia
> 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
>
> ```
---
### Author:  [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)
#### Post date: [November 28, 2024, 7:32am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/15 "2024-11-28T07:32:35Z")
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:
```julia
@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.
---
### Author:  [@cryptic.ax](https://discourse.julialang.org/u/cryptic.ax)
#### Post date: [November 28, 2024, 7:51am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/16 "2024-11-28T07:51:06Z")
[feat: add automatic variable discovery for `OptimizationSystem` by AayushSabharwal · Pull Request #3243 · SciML/ModelingToolkit.jl · GitHub](https://github.com/SciML/ModelingToolkit.jl/pull/3243) adds the required functionality.
---
### Author:  [@Bart\_van\_de\_Lint](https://discourse.julialang.org/u/Bart_van_de_Lint)
#### Post date: [November 28, 2024, 8:53am UTC](https://discourse.julialang.org/t/optimizationsystem-from-a-set-of-multiple-equations/122809/17 "2024-11-28T08:53:30Z")
Thanks!