I’m converting a model from Modelica to ModelingToolkit [MTK]. The Modelica-version works fine using OpenModelica, but the ModelingToolkit has problems wrt. initialization.
Here is the Julia+MTK code that is equivalent to the Modelica code (including initial conditions).
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using DifferentialEquations
#
# Friction force computation
#
@register_symbolic fForce(v,rho,mu,D,epsilon,Aw)
function fForce(v,rho,mu,D,epsilon,Aw)
# Friction force
#= Function input arguments
v : velocity, m/s
rho : density, kg/m3
mu : viscositiy, Pa.s
D : pipe diameter, m
epsilon : pipe roughness heigth, m
Aw : pipe wetting surface, m2
=#
#= Function return (response) value
F_f : Friction force, N
=#
#= Local (protected) quantities
N_Re : Reynolds number
arg : temporary variable
fD : Darcy friction factor
=#
N_Re = rho*abs(v)*D/mu
if N_Re < 2.3e3
F_f = 8mu*Aw/D*v
else
arg = epsilon/3.7/D + 5.74/(N_Re + 1e-3)^0.9
if arg <= 0
fD = 0
else
fD = 1/(2 * log10(arg))^2
end
F_f = rho*v*abs(v)*Aw*fD/8
end
return F_f
end
#
# Registering inputs
#
@register_symbolic u_u_v(t) # Turbine valve signal function
@register_symbolic u_Wd_e(t) # Electrical power usage, W
@register_symbolic u_H_t(t) # Tailwater level, m
#
u_u_v(t) = 0.89062946
u_Wd_e(t) = 111.729e6
u_H_t(t) = 5
#
# MTK model
#
@mtkmodel HydropowerSystem begin
# Model parameters
# -- structure parameters
@structural_parameters begin
# :io = standard input-output, :zero = zero dynamics
# :lin_io = linearizing input-output, :lin_zero = linearizing zero dynamics
case = :io
end
#=
# -- constants
@constants begin
# Constant to keep or kill surge tank
keep_surge_tank = 1, [description = "Parameter to keep surge tank, -"]
g = 9.81, [description = "Acceleration of gravity, m/s2"]
end
=#
#
# -- parameters
@parameters begin
# Constant to keep or kill surge tank
g = 9.81, [description = "Acceleration of gravity, m/s2"]
# General parameters
gamma_sw = 1, [description = "Swirl constant,"]
k_g = 0.675
gamma_f = 2, [description = "Friction constant, -"]
gamma_sh = 0.88, [description = "shock constant, -"]
rho = 997, [description = "Water density, kg/m3"]
mu = 8.9e-4, [description = "Viscosity of water, Pa.s"]
epsilon = 1.5e-5, [description = "Pipe roughness height, m"]
H_r = 48, [description = "Reservoir level above intake, m"]
p_a = 1.01325e5, [description = "Atmospheric pressure, Pa"]
#
# Rated values
Vd_n = 24.3, [description = "Nominal flow rate, m3/s"]
Np = 6, [description = "Number of pole pairs, -"]
f_0 = 50, [description = "Electric freqency, Hz"]
w_n = f_0*2pi/Np, [description = "Nominal angular velocity, rad/s"]
W_nom = 104.4e6, [description = "Nominal power output, MW"]
#
# Intake race
H_i = 23, [description = "Vertical drop of intake race, m"]
L_i = 6600, [description = "Length of intake race, m"]
D_i = 5.8, [description = "Diameter of intake race, m"]
A_i = pi*D_i^2/4, [description = "Cross section area of intake race, m2"]
V_i = A_i*L_i, [description = "Volume of intake race, m3"]
m_i = rho*V_i, [description = "Mass of intake race, kg"]
A_wi = pi*D_i*L_i, [description = "Wetting surface of intake race, m2"]
p_i1 = p_a + rho*g*H_r, [description = "Intake race pressure, Pa"]
#
# Penstock
H_p = 428.5, [description = "Vertical drop of penstock, m"]
L_p = 600, [description = "Length of penstock, m"]
D_p = 3, [description = "Diameter of penstock, m"]
A_p = pi*D_p^2/4, [description = "Cross section area of penstock, m2"]
V_p = A_p*L_p, [description = "Volume of penstock, m3"]
m_p = rho*V_p, [description = "Mass of penstock, kg"]
A_wp = pi*D_p*L_p, [description = "Wetting surface of penstock, m2"]
#
# Discharge
H_d = 0.5, [description = "Vertical drop of discharge race, m"]
L_d = 600, [description = "Length of discharge race, m"]
D_d = 5.8, [description = "Diameter of discharge race, m"]
A_d = pi*D_d^2/4, [description = "Cross section area of discharge race, m2"]
V_d = A_d*L_d, [description = "Volume of discharge race, m3"]
m_d = rho*V_d, [description = "Mass of discharge race, kg"]
A_wd = pi*D_d*L_d, [description = "Wetting surface of discharge race, m2"]
#
# Aggregate
J_a = 2e5, [description = "Moment of inertia of aggregate, kg.m2"]
k_ba = 1e3, [description = "Friction factor in the aggregate bearing box, W.s3/rad3"]
eta_e = 0.99, [description = "Electricity generator efficiency, -"]
#
# Turbine
w_1 = 0.25, [description = "Width of the turbine/blades inlet, m"]
beta_1 = 62.85, [description = "Turbine inlet angle, degree"]
beta_2 = 17.5, [description = "Turbine outlet angles, degree"]
r_1 = 1.32, [description = "Radius of the turbine outer, m"]
r_2 = 0.777, [description = "Radius of the turbine inner, m"]
r_v = 1.445, [description = "Radius of the outer guide vanes, m"]
A_1 = 2pi*r_1*w_1, [description = "Maximal area of opening, m2"]
end
#
# Model variables, with initial values needed
@variables begin
# Variables
#
md_i(t), [description = "Inlet race mass flow rate, kg/s"]
md_p(t), [description = "Penstock mass flow rate, kg/s"]
md_d(t), [description = "Discharge mass flow rate, kg/s"]
#
Vd_i(t)=Vd_n, [description = "Inlet race volumetric flow rate, m3/s"]
Vd_p(t), [description = "Penstock volumetric flow rate, m3/s"]
Vd_d(t), [description = "Discharge volumetric flow rate, m3/s"]
#
v_i(t), [description = "Inlet race linear velocity, m/s"]
v_p(t), [description = "Penstock linear velocity, m/s"]
v_d(t), [description = "Discharge linear velocity, m/s"]
#
N_Re_i(t), [description = "Inlet race Reynolds number, -"]
N_Re_p(t), [description = "Penstock Reynolds number, -"]
N_Re_d(t), [description = "Discharge Reynolds number, -"]
#
F_fi(t), [description = "Inlet race tank friction force, N"]
F_fp(t), [description = "Penstock tank friction force, N"]
F_fd(t), [description = "Discharge tank friction force, N"]
#
F_gi(t), [description = "Inlet race tank gravity force, N"]
F_gp(t), [description = "Penstock tank gravity force, N"]
F_gd(t), [description = "Discharge tank gravity force, N"]
#
M_i(t), [description = "Inlet race tank linear momentum, kg.m/s"]
M_p(t), [description = "Penstock tank linear momentum, kg.m/s"]
M_d(t), [description = "Discharge tank linear momentum, kg.m/s"]
#
p_m(t), [description = "Pressure in manifold, Pa"]
p_i2(t), [description = "Inlet race exit pressure, Pa"]
p_p1(t), [description = "Penstock inlet pressure, Pa"]
p_p2(t), [description = "Penstock outlet pressure, Pa"]
p_tr2(t), [description = "Turbine outlet pressure, Pa??"]
p_d1(t), [description = "Discharge inlet pressure, Pa"]
p_d2(t), [description = "Discharge outlet pressure, Pa"]
#
Dp_tt(t), [description = "Turbine total pressure drop, Pa"]
Dp_tf(t), [description = "Turbine friction pressure drop, Pa"]
Dp_sw(t), [description = "Turbine swirl pressure drop, Pa"]
Dp_sh(t), [description = "Turbine shock pressure drop, Pa"]
Dp_ts(t), [description = "Turbine shaft pressure drop, Pa"]
#
alpha_1(t), [description = "Guide vane opening angle, deg"]
gamma(t), [description = "Inlet flow angle, deg"]
#
v_1(t), [description = "XX linear velocity, m/s"]
v_w1(t), [description = "XX linear velocity, m/s"]
v_w2(t), [description = "XX linear velocity, m/s"]
# Variables related to height loss ex turbine
Dp_i(t), [description = "Intake race pressure drop, Pa"]
Dp_p(t), [description = "Penstock pressure drop, Pa"]
Dp_d(t), [description = "Discharge pressure drop, Pa"]
Dp_xtr(t), [description = "XX pressure drop, Pa"]
H_xtr(t), [description = "XX height drop, m"]
H_tr(t), [description = "XX heigth drop, m"]
H_n(t), [description = "XX nominal heigth drop, m"]
K_a(t), [description = "Aggregate kinetic energy, J"]
Wd_ts(t), [description = "Turbine shaft power, W"]
Wd_af(t), [description = "Aggregate friction power, W"]
# Real Wd_e; // Work powers
w(t)=w_n, [description = "Aggregate angular velocity, rad/s"]
Wd_tt(t), [description = "Turbine total power, W"]
# Inputs
u_v(t), [description = "Input guide vane signal, -"]
Wd_e(t), [description = "Electrical power usage, W"]
H_t(t), [description = "Tail water level, m"]
end
#
# Model equations
@equations begin
# Discharge counter pressure
p_d2 ~ p_a + rho*g*H_t
# Mass flow rates
md_i ~ rho*Vd_i
md_p ~ rho*Vd_p
md_d ~ rho*Vd_d
# Linear velocities
Vd_i ~ A_i*v_i
Vd_p ~ A_p*v_p
Vd_d ~ A_d*v_d
# Manifold
p_i2 ~ p_m
p_p1 ~ p_m
Vd_i ~ Vd_p
# Penstock-Discharge string
Vd_d ~ Vd_p
# Momentums
M_i ~ m_i*v_i
M_p ~ m_p*v_p
M_d ~ m_d*v_d
# Gravitational forces
F_gi ~ m_i*g*H_i/L_i
F_gp ~ m_p*g*H_p/L_p
F_gd ~ m_d*g*H_d/L_d
# Reynold's numbers
N_Re_i ~ rho*abs(v_i)*D_i/mu
N_Re_p ~ rho*abs(v_p)*D_p/mu
N_Re_d ~ rho*abs(v_d)*D_d/mu
# Friction forces
F_fi ~ fForce(v_i,rho,mu,D_i,epsilon,A_wi)
F_fp ~ fForce(v_p,rho,mu,D_p,epsilon,A_wp)
F_fd ~ fForce(v_d,rho,mu,D_d,epsilon,A_wd)
# Turbine
#alpha versus guide vane signal
alpha_1 ~ (90-90*u_v) # "relation between alpha and u_v"
v_1 ~ Vd_p/(A_1*sind(alpha_1)) # "linear velocity"
gamma ~ atand(1/((v_w1/(sind(alpha_1)*v_1)) - (1/tand(alpha_1) )))
v_w1 ~ w*r_1
v_w2 ~ w*r_2
# Turbine shaft power
Wd_ts ~ rho*Vd_p*w*( r_1*Vd_p/(A_1*tand(alpha_1))-gamma_sw*k_g*r_2^2*w*(1- (w_n*Vd_p/(w*Vd_n))) )
# Turbine pressure losses
Dp_tf ~ gamma_f*0.5*rho*(Vd_p/(pi*r_2^2))^2
Dp_sh ~ gamma_sh*0.5*rho*((sind(gamma-beta_1)^2)/(sind(beta_1)^2))*(Vd_p/(A_1*sind(gamma)))^2
Dp_sw ~ rho*(gamma_sw*k_g*v_w2)^2*(1-(w_n*Vd_p/(w*Vd_n)))^2
Dp_ts ~ Wd_ts/Vd_p
Dp_tt ~ Dp_tf + Dp_sh + Dp_sw + Dp_ts
Wd_tt ~ Vd_p*Dp_tt
Dp_tt ~ p_p2-p_tr2
p_d1 ~ p_tr2
#Pressure friction losses in pipes
Dp_i ~ F_fi/A_i
Dp_p ~ F_fp/A_p
Dp_d ~ F_fd/A_d
Dp_xtr ~ Dp_i + Dp_p + Dp_d
H_xtr ~ Dp_xtr/rho/g
H_n ~ H_i + H_p
H_tr ~ H_n + H_r - H_t - H_xtr
# Aggregate
K_a ~ J_a*w^2/2
Wd_af ~ k_ba*w*abs(w)/2 # k_ba*w^2/2
# Balance laws
Dt(M_i) ~ (p_i1 - p_i2)*A_i + F_gi - F_fi
Dt(M_p) ~ (p_p1 - p_p2)*A_p + F_gp - F_fp
Dt(M_d) ~ (p_d1 - p_d2)*A_d + F_gd - F_fd
Dt(K_a) ~ Wd_ts - Wd_af - Wd_e
# Injecting inputs
#if case == :io
# Standard input-output simulation
u_v ~ u_u_v(t)
#elseif case == :lin_io
# Linearization of input-output model
# u_v ~ u_u_v(t)
#end
Wd_e ~ u_Wd_e(t)
H_t ~ u_H_t(t)
end
end
#
# Instantiation, etc.
#
@mtkbuild hps = HydropowerSystem()
tspan = (0,15)
prob_sys = ODEProblem(hps, [], tspan)
The model looks as follows:
The dynamics seems ok – one ode for rotational (kinetic) energy, and one for velocity + some algebraic equations.
In OpenModelica, it suffices to provide initial conditions for Vd_p which is a volumetric flowrate, directly proportional to v_d, and w which is the angular velocity, directly related to K_a. In other words, with Vd_p and w specified, it is trivial to find initial conditions for K_a and v_d.
Julia/MTK complains that
ArgumentError: SymbolicUtils.BasicSymbolic{Real}[p_d1(0)] are either missing from the variable map or missing from the system's unknowns/parameters list.
OK… if I try to specify an initial value for p_d1, I get another error message, related to over specified system …
Warning: Initialization system is overdetermined. 2 equations for 1 unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false
If I try to solve this overdetermined system (default solver or Rodas5()), there is no solution – sol.t contains the single element 0.0.
Question:
- What is a good procedure to understand which initial conditions are needed?