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