Closed
Description
The following attempt to get the initial values gives the error: ERROR: ArgumentError: drho_1(t) is neither an observed nor an unknown variable.
. The variable drho_1(t)
is a part of the system, I should be able to retrieve the initial value.
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using NonlinearSolve
using Plots
reg_pow(x,a) = sign(x)*abs(x)^a
reg_pow_safe(x, a, delta = 0.1) = ifelse(abs(x/delta) >= 1, sign(x)*abs(x/delta)^a * delta^a, (delta^a*x)/delta)
pars = @parameters begin
ρ₀ = 1000
β = 2e9
A = 0.1
m = 100
L = 1
p_s = 100e5
p_r = 10e5
C1 = 1.35
C2 = 1.35
c = 1000
A_p = 0.00094
transition=0.1
end
vars = @variables begin
p_1(t) = p_s
p_2(t) = p_r
x(t)=0
dx(t)=0
ddx(t), [guess=0]
rho_1(t), [guess=ρ₀]
rho_2(t), [guess=ρ₀]
drho_1(t), [guess=0]
drho_2(t), [guess=0]
m_1(t), [guess=0]
m_2(t), [guess=0]
u_1(t), [guess=0]
u_2(t), [guess=0]
end
eqs = [
D(x) ~ dx
D(dx) ~ ddx
D(m_1) ~ u_1*ρ₀*A_p
D(m_2) ~ u_2*ρ₀*A_p
# volume & mass
m_1 ~ rho_1*(L+x)*A
m_2 ~ rho_2*(L-x)*A
# density
rho_1 ~ ρ₀*(1 + p_1/β)
rho_2 ~ ρ₀*(1 + p_2/β)
# actuator force balance
m*ddx ~ (p_1 - p_2)*A - c*dx
# orifices
(p_s - p_1) ~ C1*ρ₀*reg_pow(u_1, 2) # orifice 1
(p_r - p_2) ~ C2*ρ₀*reg_pow(u_2, 2) # orifice 2
]
@mtkbuild sys = ODESystem(eqs, t, vars, pars)
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob, FastShortcutNonlinearPolyalg()) # Success
initsol[sys.drho_1] #ERROR: ArgumentError: drho_1(t) is neither an observed nor an unknown variable.
Versions...
[961ee093] ModelingToolkit v9.61.0
[8913a72c] NonlinearSolve v4.3.0
[1dea7af3] OrdinaryDiffEq v6.90.1