MTK & linearization -- MTKParameters access out of bounds

[B-LIE]

Linearization of a semi-complex model fails -- is this a bug, or do I do something wrong❓

I have done linearization of a "simple" model as follows:

1. I use @mtkmodel to create a model (Tank), which I instantiate as tank via @mtkbuild. Model tank is reduced to an ODE with a single unknown. Simulation of tank works fine.
2. I then use @mtkmodel to create a new model (Tank_noi with no input), i.e., Tank_noi is identical to Tank except that I have commented out the equation of Tank where the input variable is given a value.
3. I instantiate the unbalanced model tank_noi via @named tank_noi = Tank_noi() followed by tank_noi = complete(tank_noi).
4. Finally, I linearize the model doing:

linearize(tank_noi, [tank_noi.md_i], [tank_noi.h]; op = Dict(tank_noi.m=>m_ss, tank_noi.md_i=>md(0)))

This works fine.

Next, I follow exactly the same strategy for a more complex model. This more complex model fails in the linearization... Here is the "identical" procedure I follow as above, but which gives an error message. I include the complex model at the bottom.

i. The more complex, balanced model is named HydropowerSystem (instead of Tank), and I instantiate it as hps (instead of tank) -- see point 1 above. For this more complex model, hps is reduced to a DAE with 4 differential equations and 3 algebraic equations with 7 unknowns, of which 2 of the unknowns have been introduced by ModelingToolkit and are time derivatives of two of my variables.
ii. The unbalanced model with no input is named HydropowerSystem_noi (instead of Tank_noi), and I instantiate it as hps_noi (instead of tank_noi), se points 2-3 above.
iii. When I try to linearize the model, I get an error message:

julia> linearize(hps_noi, [hps_noi.u_v], [hps_noi.w], op=Dict(hps_noi.u_v=>0.89))
BoundsError: attempt to access ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}} at index [2]

Stacktrace:
  [1] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\6UxiU\src\systems\parameter_buffer.jl:0 [inlined]
  [2] getindex
    @ C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\6UxiU\src\systems\parameter_buffer.jl:681 [inlined]
  [3] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\SymbolicUtils\EGhOJ\src\code.jl:375 [inlined]
  [4] macro expansion
    @ C:\Users\Bernt_Lie\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ .\none:0 [inlined]
  [6] generated_callfunc
    @ .\none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8959787b, 0xbf2e6982, 0x8e8ecb53, 0x1b78ba97, 0x5c827a52), Nothing})(::Vector{Float64}, ::ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, ::Float64)
    @ RuntimeGeneratedFunctions C:\Users\Bernt_Lie\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:150
  [8] (::ModelingToolkit.var"#339#347"{ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x037a3fd5, 0x6a56ca50, 0x92f48148, 0x0ddc0727, 0xe833d44a), Nothing}, SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}}}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8959787b, 0xbf2e6982, 0x8e8ecb53, 0x1b78ba97, 0x5c827a52), Nothing}})(u::Vector{Float64}, p::ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, t::Float64)
    @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\6UxiU\src\systems\abstractsystem.jl:2259
  [9] (::ModelingToolkit.var"#341#349"{UnitRange{Int64}, UnitRange{Int64}, Vector{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, ModelingToolkit.var"#339#347"{ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x037a3fd5, 0x6a56ca50, 0x92f48148, 0x0ddc0727, 0xe833d44a), Nothing}, SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}}}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8959787b, 0xbf2e6982, 0x8e8ecb53, 0x1b78ba97, 0x5c827a52), Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#775"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9e666d13, 0xf85a5fec, 0xdc6334ed, 0x5f3aad02, 0x2dbbfcab), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x92b43b5f, 0xd06be8e9, 0xc149ed4d, 0xdfa6e6f5, 0x2cfe6464), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#608"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6f213571, 0x4ad2b8da, 0x4f151bf6, 0x6ad63a30, 0x3d3983a8), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x94a28684, 0x6cefae0c, 0x7e19d41a, 0x9d81884f, 0x4de48cb1), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{NonlinearSystem}, Nothing, NonlinearSystem, Vector{Float64}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#11238017669661590846"), Symbol("##arg#5890982026656049690")), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xedcaabcc, 0xfa19692d, 0x2a6e82b4, 0x727a617b, 0x9a61a907), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf816b0af, 0xe66bcce8, 0x3ea1015e, 0x1a04f2ec, 0x02f4e157), Nothing}, ForwardDiff.Chunk{1}, ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, Bool, GeneralizedFirstOrderAlgorithm{nothing, :TrustRegion, Missing, NonlinearSolve.GenericTrustRegionScheme{NonlinearSolve.RadiusUpdateSchemes.__Simple, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing, Nothing, Nothing, Nothing}, Dogleg{NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, SteepestDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}}, Nothing, Nothing, Nothing}, ODESystem})(u::Vector{Float64}, p::Dict{Num, Float64}, t::Float64)
    @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\6UxiU\src\systems\abstractsystem.jl:2313
 [10] linearize(sys::ODESystem, lin_fun::ModelingToolkit.var"#341#349"{UnitRange{Int64}, UnitRange{Int64}, Vector{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, ModelingToolkit.var"#339#347"{ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x037a3fd5, 0x6a56ca50, 0x92f48148, 0x0ddc0727, 0xe833d44a), Nothing}, SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}}}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8959787b, 0xbf2e6982, 0x8e8ecb53, 0x1b78ba97, 0x5c827a52), Nothing}}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#775"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9e666d13, 0xf85a5fec, 0xdc6334ed, 0x5f3aad02, 0x2dbbfcab), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x92b43b5f, 0xd06be8e9, 0xc149ed4d, 0xdfa6e6f5, 0x2cfe6464), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#608"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6f213571, 0x4ad2b8da, 0x4f151bf6, 0x6ad63a30, 0x3d3983a8), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x94a28684, 0x6cefae0c, 0x7e19d41a, 0x9d81884f, 0x4de48cb1), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{NonlinearSystem}, Nothing, NonlinearSystem, Vector{Float64}}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#11238017669661590846"), Symbol("##arg#5890982026656049690")), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xedcaabcc, 0xfa19692d, 0x2a6e82b4, 0x727a617b, 0x9a61a907), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#9952546503938205633"), Symbol("##arg#17241324348203558476"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf816b0af, 0xe66bcce8, 0x3ea1015e, 0x1a04f2ec, 0x02f4e157), Nothing}, ForwardDiff.Chunk{1}, ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, Bool, GeneralizedFirstOrderAlgorithm{nothing, :TrustRegion, Missing, NonlinearSolve.GenericTrustRegionScheme{NonlinearSolve.RadiusUpdateSchemes.__Simple, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing, Nothing, Nothing, Nothing}, Dogleg{NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, SteepestDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}}, Nothing, Nothing, Nothing}, ODESystem}; t::Float64, op::Dict{Num, Float64}, allow_input_derivatives::Bool, p::SciMLBase.NullParameters)
    @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\6UxiU\src\systems\abstractsystem.jl:2594
 [11] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; op::Dict{Num, Float64}, t::Float64, allow_input_derivatives::Bool, zero_dummy_der::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\Bernt_Lie\.julia\packages\ModelingToolkit\6UxiU\src\systems\abstractsystem.jl:2645
 [12] top-level scope
    @ c:\Users\Bernt_Lie\OneDrive\Documents\researchUSN\Supervisor\MScThesis\Fall2023\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X40sZmlsZQ==.jl:2

A potential problem here may be that the time derivatives of the variables v_i and Vd_d have not been specified. If I do:

linearize(hps_noi, [hps_noi.u_v], [hps_noi.w], op=Dict(hps_noi.u_v=>0.89, Dt(hps_noi.v_i)=>0, Dt(hps_noi.Vd_d)=>0))

I get a different error message...

BoundsError: attempt to access Tuple{SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, NonlinearLeastSquaresProblem{Nothing, true, ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa148094c, 0xba54127a, 0xe432ae56, 0x70bc53cc, 0xa7efd9f2), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a3c5c5b, 0x93608ad1, 0xbcc33a9c, 0x1971196f, 0xe3aae881), Nothing}}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#608"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd837c3bf, 0x60809d7d, 0x8877b772, 0x1889088f, 0xe46171ee), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd8d967f, 0x4d7a5092, 0x22ac55bd, 0xa045a117, 0x7d12ef84), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{NonlinearSystem}, Nothing, NonlinearSystem, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing}} at index [2]
...

QUESTION:

• Is this a bug with linearization, or if not -- what is it that I do that is wrong?

---

CODE

I use the latest versions of Julia (v10.4) and ModelingToolkit (v9.30.1). The model below is not finished (i.e., some descriptions are missing, and some changes needs to be included), but it should work.

# Import packages
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as Dt
using DifferentialEquations
# using Plots
# 
# Function to compute friction force
@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
# 
# Input function for model
@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
# 
# Model instantiator
@mtkmodel HydropowerSystem begin
    # Model 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
    =#
    #
    # Model parameters
    @parameters 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"]
        # 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"]
        #
        # Surge tank
        H_s = 120,              [description = "Vertical drop of surge tank, m"]
        L_s = 140,              [description = "Length of surge tank, m"]
        D_s = 3.4,              [description = "Diameter of the surge tank, m"]
        A_s = pi*D_s^2/4,       [description = "Cross section area of surge tank, m2"]
        p_s2 = p_a,             [description = "Outlet pressure of surge tank, 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
        m_s(t)=rho*A_s*L_s/H_s*(H_r+H_i),   [description = "Mass in surge tank, kg"]
        V_s(t),             [description = "Volume in surge tank, m3"]
        ls_s(t),            [description = "Length of water string in surge tank, m"]
        # h_s(t)=H_r+H_i,     [description = "Level in surge tank, m"]
        h_s(t),             [description = "Level in surge tank, m"]
        A_ws(t),            [description = "Wetting surface of surge tank string, m2"]
        #
        md_i(t),            [description = "Inlet race mass flow rate, kg/s"]
        md_s(t),            [description = "Surge tank 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_i(t),            [description = "Inlet race volumetric flow rate, m3/s"]
        Vd_s(t),            [description = "Surge tank volumetric flow rate, m3/s"]
        # Vd_p(t)=Vd_n,       [description = "Penstock 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_s(t)=0,             [description = "Surge tank linear velocity, m/s"]
        v_p(t),             [description = "Penstock linear velocity, m/s"]
        v_d(t)=Vd_n/A_d,    [description = "Discharge linear velocity, m/s"]
        #
        N_Re_i(t),          [description = "Inlet race Reynolds number, -"]
        N_Re_s(t),          [description = "Surge tank 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_fs(t),            [description = "Surge 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_gs(t),            [description = "Surge 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_s(t),             [description = "Surge 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"]
        #
        Md_s(t),            [description = "Surge tank linear momentum flow, kg.m/s2"]
        #
        p_m(t),             [description = "Pressure in manifold, Pa"]
        p_i2(t),            [description = "Inlet race exit pressure, Pa"]
        p_s1(t),            [description = "Surge tank inlet 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_tr(t),           [description = "Turbine pressure drop, Pa"]
        Dp_ft(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"]
        #
        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_s(t),            [description = "Surge tank 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)=J_a*w_n^2/2, [description = "Aggregate kinetic energy, J"]
        Wd_s(t),            [description = "XX work, W"]
        Wd_fa(t),           [description = "XX work, W"]
        # Real Wd_e; // Work powers
        w(t)=w_n,           [description = "Aggregate angular velocity, rad/s"]
        Wd_t(t),            [description = "XX work, 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
        #w ~ w_n
        # Surge tank liquid string
        m_s ~ rho*V_s
        V_s ~ A_s*ls_s
        h_s ~ ls_s*H_s/L_s 
        A_ws ~ pi*D_s*ls_s # "Wetting surface of surge tank string, m2"
        # Discharge counter pressure
        p_d2 ~ p_a + rho*g*H_t
        # Mass flow rates
        md_i ~ rho*Vd_i
        md_s ~ rho*Vd_s
        md_p ~ rho*Vd_p
        md_d ~ rho*Vd_d
        # Linear velocities
        Vd_i ~ A_i*v_i
        Vd_s ~ A_s*v_s
        Vd_p ~ A_p*v_p
        Vd_d ~ A_d*v_d
        # Manifold
        p_i2 ~ p_m
        p_s1 ~ p_m
        p_p1 ~ p_m
        Vd_i ~ Vd_s + Vd_p
        # Penstock-Discharge string
        Vd_d ~ Vd_p
        # Momentums
        M_i ~ m_i*v_i
        M_s ~ m_s*v_s
        M_p ~ m_p*v_p
        M_d ~ m_d*v_d
        # Momentum flow rate into/out of surge tank
        Md_s ~ rho*Vd_s*abs(Vd_s)/A_s
        # Gravitational forces
        F_gi ~ m_i*g*H_i/L_i
        F_gs ~ m_s*g*H_s/L_s
        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_s ~ rho*abs(v_s)*D_s/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_fs ~ fForce(v_s,rho,mu,D_s,epsilon,A_ws)*5e3
        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_s ~ 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_ft ~ 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_tr ~ Dp_ft+Dp_sh+Dp_sw+(Wd_s/Vd_p)
        Wd_t ~ Vd_p*(Dp_ft + Dp_sh + Dp_sw) + Wd_s
        Dp_tr ~ p_p2-p_tr2
        p_d1 ~ p_tr2
        #Pressure friction losses in pipes
        Dp_i ~ F_fi/A_i
        Dp_s ~ F_fs/A_s
        Dp_p ~ F_fp/A_p
        Dp_d ~ F_fd/A_d
        Dp_xtr ~ Dp_i+Dp_p + Dp_s+Dp_d
        H_xtr ~ Dp_xtr/rho/g
        H_n ~ H_i+H_p+H_s
        H_tr ~ H_n + H_r - H_t - H_xtr
        # Aggregate
        K_a ~ J_a*w^2/2
        Wd_fa ~ k_ba*w^2/2
        # Balance laws
        Dt(M_i) ~ (p_i1 - p_i2)*A_i + F_gi - F_fi
        Dt(m_s) ~ md_s
        Dt(M_s) ~ (Md_s + (p_s1 - p_s2)*A_s - F_gs - F_fs)*keep_surge_tank
        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_s - Wd_fa - Wd_e
        # Injecting inputs 
        u_v ~ u_u_v(t)
        Wd_e ~ u_Wd_e(t)
        H_t ~ u_H_t(t)
    end
end
# 
# Building model
@mtkbuild hps = HydropowerSystem(; keep_surge_tank=1e-10)
# 
# Setting up and solving the model
#= 
tspan = (0,1_000)

prob_sys = ODEProblem(hps, [Dt(hps.Vd_d)=>0, Dt(hps.v_i)=>0], tspan)

sol = solve(prob_sys, Rodas5())

plot(sol, idxs=(hps.h_s))
=#
# 
# Instantiator for *unbalanced* model used for linearization -- one of the input equations is commented out
@mtkmodel HydropowerSystem_noi begin
    # Model 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
    =#
    #
    # Model parameters
    @parameters 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"]
        # 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"]
        #
        # Surge tank
        H_s = 120,              [description = "Vertical drop of surge tank, m"]
        L_s = 140,              [description = "Length of surge tank, m"]
        D_s = 3.4,              [description = "Diameter of the surge tank, m"]
        A_s = pi*D_s^2/4,       [description = "Cross section area of surge tank, m2"]
        p_s2 = p_a,             [description = "Outlet pressure of surge tank, 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
        m_s(t)=rho*A_s*L_s/H_s*(H_r+H_i),   [description = "Mass in surge tank, kg"]
        V_s(t),             [description = "Volume in surge tank, m3"]
        ls_s(t),            [description = "Length of water string in surge tank, m"]
        # h_s(t)=H_r+H_i,     [description = "Level in surge tank, m"]
        h_s(t),             [description = "Level in surge tank, m"]
        A_ws(t),            [description = "Wetting surface of surge tank string, m2"]
        #
        md_i(t),            [description = "Inlet race mass flow rate, kg/s"]
        md_s(t),            [description = "Surge tank 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_i(t),            [description = "Inlet race volumetric flow rate, m3/s"]
        Vd_s(t),            [description = "Surge tank volumetric flow rate, m3/s"]
        # Vd_p(t)=Vd_n,       [description = "Penstock 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_s(t)=0,             [description = "Surge tank linear velocity, m/s"]
        v_p(t),             [description = "Penstock linear velocity, m/s"]
        v_d(t)=Vd_n/A_d,    [description = "Discharge linear velocity, m/s"]
        #
        N_Re_i(t),          [description = "Inlet race Reynolds number, -"]
        N_Re_s(t),          [description = "Surge tank 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_fs(t),            [description = "Surge 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_gs(t),            [description = "Surge 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_s(t),             [description = "Surge 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"]
        #
        Md_s(t),            [description = "Surge tank linear momentum flow, kg.m/s2"]
        #
        p_m(t),             [description = "Pressure in manifold, Pa"]
        p_i2(t),            [description = "Inlet race exit pressure, Pa"]
        p_s1(t),            [description = "Surge tank inlet 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_tr(t),           [description = "Turbine pressure drop, Pa"]
        Dp_ft(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"]
        #
        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_s(t),            [description = "Surge tank 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)=J_a*w_n^2/2, [description = "Aggregate kinetic energy, J"]
        Wd_s(t),            [description = "XX work, W"]
        Wd_fa(t),           [description = "XX work, W"]
        # Real Wd_e; // Work powers
        w(t)=w_n,           [description = "Aggregate angular velocity, rad/s"]
        Wd_t(t),            [description = "XX work, 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
        #w ~ w_n
        # Surge tank liquid string
        m_s ~ rho*V_s
        V_s ~ A_s*ls_s
        h_s ~ ls_s*H_s/L_s 
        A_ws ~ pi*D_s*ls_s # "Wetting surface of surge tank string, m2"
        # Discharge counter pressure
        p_d2 ~ p_a + rho*g*H_t
        # Mass flow rates
        md_i ~ rho*Vd_i
        md_s ~ rho*Vd_s
        md_p ~ rho*Vd_p
        md_d ~ rho*Vd_d
        # Linear velocities
        Vd_i ~ A_i*v_i
        Vd_s ~ A_s*v_s
        Vd_p ~ A_p*v_p
        Vd_d ~ A_d*v_d
        # Manifold
        p_i2 ~ p_m
        p_s1 ~ p_m
        p_p1 ~ p_m
        Vd_i ~ Vd_s + Vd_p
        # Penstock-Discharge string
        Vd_d ~ Vd_p
        # Momentums
        M_i ~ m_i*v_i
        M_s ~ m_s*v_s
        M_p ~ m_p*v_p
        M_d ~ m_d*v_d
        # Momentum flow rate into/out of surge tank
        Md_s ~ rho*Vd_s*abs(Vd_s)/A_s
        # Gravitational forces
        F_gi ~ m_i*g*H_i/L_i
        F_gs ~ m_s*g*H_s/L_s
        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_s ~ rho*abs(v_s)*D_s/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_fs ~ fForce(v_s,rho,mu,D_s,epsilon,A_ws)*5e3
        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_s ~ 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_ft ~ 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_tr ~ Dp_ft+Dp_sh+Dp_sw+(Wd_s/Vd_p)
        Wd_t ~ Vd_p*(Dp_ft + Dp_sh + Dp_sw) + Wd_s
        Dp_tr ~ p_p2-p_tr2
        p_d1 ~ p_tr2
        #Pressure friction losses in pipes
        Dp_i ~ F_fi/A_i
        Dp_s ~ F_fs/A_s
        Dp_p ~ F_fp/A_p
        Dp_d ~ F_fd/A_d
        Dp_xtr ~ Dp_i+Dp_p + Dp_s+Dp_d
        H_xtr ~ Dp_xtr/rho/g
        H_n ~ H_i+H_p+H_s
        H_tr ~ H_n + H_r - H_t - H_xtr
        # Aggregate
        K_a ~ J_a*w^2/2
        Wd_fa ~ k_ba*w^2/2
        # Balance laws
        Dt(M_i) ~ (p_i1 - p_i2)*A_i + F_gi - F_fi
        Dt(m_s) ~ md_s
        Dt(M_s) ~ (Md_s + (p_s1 - p_s2)*A_s - F_gs - F_fs)*keep_surge_tank
        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_s - Wd_fa - Wd_e
        # Injecting inputs 
        # u_v ~ u_u_v(t)
        Wd_e ~ u_Wd_e(t)
        H_t ~ u_H_t(t)
    end
end
# 
# Instantiating unbalanced model
@named hps_noi = HydropowerSystem_noi()
hps_noi = complete(hps_noi)
# 
# Trying to linearize model... either of the following two attempts fail:
mats, hps_ = linearize(hps_noi, [hps_noi.u_v], [hps_noi.w], op=Dict(hps_noi.u_v=>0.89))

mats, hps_ = linearize(hps_noi, [hps_noi.u_v], [hps_noi.w], op=Dict(hps_noi.u_v=>0.89, Dt(hps_noi.v_i)=>0, Dt(hps_noi.Vd_d)=>0))