ExtraVariablesSystemException: The system is unbalanced. There are 47 highest order derivative variables and 44 equations

Hi,
I am solving a large nonlinear system of equations. It contains 44 variables and equations.
But I face the error :
ExtraVariablesSystemException: The system is unbalanced. There are 47 highest order derivative variables and 44 equations.

using Pkg
Pkg.add(["DataFrames", "XLSX"])
Pkg.add("Optim")
##############Reading data from excel #################
using DataFrames, XLSX

df = DataFrame(XLSX.readtable("??????\SolverJ\\psc_ft_scenario_new.xlsx", 1))

i = 1 # Julia indexing starts at 1, reading te first row
Qsl = df[i, "Qsl"]
Qdil = df[i, "Qdil"]
Qmo = df[i, "Qm"]
Qto = df[i, "Qt"]
Qft = df[i, "Qft"]
α_bsl_1 = df[i, "a_bsl_1"]
α_bsl_2 = df[i, "a_bsl_2"]
α_bsl_3 = df[i, "a_bsl_3"]
α_bsl_4 = df[i, "a_bsl_4"]
α_ssl_1 = df[i, "a_ssl_1"]
α_ssl_2 = df[i, "a_ssl_2"]
α_ssl_3 = df[i, "a_ssl_3"]

Qfd = Qsl + Qdil

# PSV FEED COMPOSITION
α_bfd_1 = α_bsl_1 * Qsl / Qfd
α_bfd_2 = α_bsl_2 * Qsl / Qfd
α_bfd_3 = α_bsl_3 * Qsl / Qfd
α_bfd_4 = α_bsl_4 * Qsl / Qfd
α_sfd_1 = α_ssl_1 * Qsl / Qfd
α_sfd_2 = α_ssl_2 * Qsl / Qfd
α_sfd_3 = α_ssl_3 * Qsl / Qfd
α_wfd = 1 - (α_bfd_1 + α_bfd_2 + α_bfd_3 + α_bfd_4 + α_sfd_1 + α_sfd_2 + α_sfd_3)

# ITERATIVELY RETRIEVE THE DATA FOR DIFFERENT SCENARIO
###########################################################
using NLsolve, NonlinearSolve, Plots, ModelingToolkit

@variables begin 
    v_bf_1 = -0.04660483043350957 
    v_bf_2 = -0.09345981917569582
    v_sf_1 = -0.08130625749545124 
    v_bmu_1 = -0.09613872404222351
    v_bmu_2 = -0.030677683191747712 
    v_smu_1 = -0.07997942776010468 
    v_bt_1 = 0.013060048582814041
    v_st_1 = 0.05644695757124584
    v_st_2 = 0.09751454192570927
    v_st_3 = 0.08123282460225396
    v_wf = -0.058046308842861466
    v_wmu = -0.00036375845116324483
    v_wt = 0.007696945147263263
    #####
    ρ_f = 960.8014204117135 
    ρ_mu = 975.0320325950338
    ρ_m = 1138.387734137675
    ρ_t= 1476.4611449382703
    #######
    α_bf = 0.7397525076061395  # 0.6397525076
    α_sf = 0.05714768930793848   # 0.57147689307938
    α_bmu = 0.6852046880649323   
    α_smu = 0.053434669298424966  # 0.43434669298424966
    α_bm = 0.07531046509265504  # makes sense
    α_sm = 0.21740611040100571
    α_bt = 0.09304537667828224    #  0.49008855975103
    α_st = 0.49008855975103  #  0.09304537667828224
    α_bf_1 = 0.018487075313053238  # 0.37705165301787147
    α_bf_2 = 0.711263232293  # 0.018487075313053238
    α_sf_1 = 0.040193096
    α_bmu_1 = 0.2058157585672766  # 0.4758001863282879
    α_bmu_2 = 0.4758001863282879  # 0.9258157585672766
    α_smu_1 = 0.05702341         # 0.8259463792344883
    α_bm_1 = 0.010464741819863431 # 0.10464741819863431 
    α_bm_2 = 0.06144341881817391  # 0.08144341881817391 
    α_sm_1 = 0.10082698              # 0.14524751024977683
    α_sm_2 = 0.010710206940301835  # 0.1871020694030183
    α_sm_3 = 0.10586891754887175  # 0.2586891754887175
    α_bt_1 = 0.0045             # 0.11805236842596298
    α_st_1 = 0.01384434783970221
    α_st_2 = 0.042886064930470547 #  0.42886064930470547
    α_st_3 = 0.60029809495899 # 0.6190129809495899
    ###
    α_wf = 0.20309981069   # 0.4277137177198411
    α_wmu = 0.2613606426   # 0.9023711547833565
    α_wm =  0.7072834245 # changed with the previous value and its interchanged
    α_wt =  0.3384714   # 0.4533852950765186
    ###
    VI = -9.541152104376794e-11
    ###
    w_bf = 0.6403761770580065 # 0.8403761770580065
    w_sf = 0.1299167398195812 # 0.3799167398195812
    w_wf = 0.255356814994774  # 0.735356814994774
    w_bm = 0.031998081785725  # 0.831998081785725
    w_sm = 0.1946105277644868 # 0.7946105277644868
    w_wm = 0.7622801692643188 # 0.5722801692643188
    w_bt = 0.0442609239538598 # 0.44426092395385985
    w_st = 0.5819373884863865 # 0.9819373884863865  
    w_wt = 0.2737807493701219 # 0.2737807493701219
    ###
    Qff = 0.399794587515803
    Qf = 0.5301731026307072
    Qm = 2.1561373280078153
    Qt = 0.12985582342677932
end
ρ_b = [1050, 700]  # % kg/m^3
d_b = [d_b_1, d_b_2] * 1e-6  # m
d_s = [d_s_1, d_s_2, d_s_3] * 1e-6 

@parameters begin
    # FIXED VARIABLE FOR PSC
    A = 7.065 * 10^2  # ; % m2
    A_t = 555
    μ = 0.000325  # ; % kg/m *s
    g = 9.81  # ; % m/s
    
    # FIXED VARIABLE FOR FT CELL
    Aft = 163  # ; % m2
    cd = 0.47  # drag coefficient
    
    # ADJUSTABLE 
    ρ_s = 2650 # [2650, 2650, 2650]  # % kg/m^3
    ρ_w = 971.8  # % kg/m^3
    d_b_1 = 21
    d_b_2 = 240
    d_s_1 = 7.5
    d_s_2 = 160
    d_s_3 = 380
    
    # Tunable parameters
    k_fines = 1.31 
    n_bf_1 = 6.05
    n_bf_2 = 5.6
    n_sf_1 = 7.3
    n_bmu_1 = 5.3
    n_bmu_2 = 4.2
    n_smu_1 = 6.4
    n_bt_1 = 6.5 
    n_st_1 = 7
    n_st_2 = 6
    n_st_3 = 5

    λ_bf = 1.0906647723893508
    λ_ff = 1.0909709890659598
    λ_bm = 0.9505347443893531
    λ_sm = 1.0616129443794053
    λ_bt = 0.7502931250710971
    λ_st = 1.0505261676728082
    λ_fm = 0.8900741387772585
    λ_ft = 1.060947106931769


end


# Define the  equations
psv_eqs = [# WATER BALANCE
    0 ~ 1 - (α_wf + (α_bf_1 + α_bf_2) + α_sf_1), # 1
    0 ~ 1 - (α_wmu + (α_bmu_1 + α_bmu_2) + α_smu_1), # 2
    0 ~ α_wm  - (1 - (α_bm_1 + α_bm_2) - (α_sm_1 + α_sm_2 + α_sm_3)), # 3
    0 ~ 1 - (α_wt + α_bt_1 + (α_st_1 + α_st_2 + α_st_3)), # 4 
    # MIXTURE DENSITY
    0 ~ ρ_f - (ρ_w * α_wf + ρ_b[1] * α_bf_1  + ρ_s * α_sf_1),                                      # 5  # ρ_b[2] * α_bf_2 + ρ_s[2] * α_sf_2 + ρ_s[3] * α_sf_3 = 0
    0 ~ ρ_mu - (ρ_w * α_wmu + ρ_b[1] * α_bmu_1 + ρ_b[2] * α_bmu_2 + ρ_s * α_smu_1),                # 6 # ρ_s * (α_smu_1 + α_smu_2 +  α_smu_3)
    0 ~ ρ_m - (ρ_w * α_wm + ρ_b[1] * α_bm_1 + ρ_b[2] * α_bm_2 + ρ_s * (α_sm_1 + α_sm_2 +  α_sm_3)), # 7 ρ_s * (α_sm_1 + α_sm_2 +  α_sm_3)
    0 ~ ρ_t - (ρ_w * α_wt + ρ_b[1] * α_bt_1  + ρ_s * (α_st_1 + α_st_2 + α_st_3)),                    # 8      ρ_b[2] * α_bf_2
    # Water (Bulf flow) Velocity
    0 ~ -v_wf - (Qf + A * (α_bmu_1 * v_bf_1 + α_bmu_2 * v_bf_2 + α_smu_1 * v_sf_1)) / (α_wf * A) , # 9 + α_smu_2 * v_sf_2 + α_smu_3 * v_sf_3
    0 ~ -v_wmu - (Qf + A * (α_bm_1 * v_bmu_1 + α_bm_2 * v_bmu_2 + α_sm_1 * v_smu_1)) / (α_wmu * A), # 10 + α_sm_2 * v_smu_2 + α_sm_3 * v_smu_3
    0 ~ -v_wt + (Qt - A_t * (α_bm_1 * v_bt_1 + α_sm_1 * v_st_1 + α_sm_2 * v_st_2 + α_sm_3 * v_st_3)) / (α_wt * A_t), # 11 + α_bm_2 * v_bt_2
    ######### SETTLING VELOCITY ##########
    # Froth
    0 ~ v_bf_1 - (g * (d_b[1] ^ 2) * (ρ_b[1] - ρ_f) / (18 * μ) * max((α_wf - (k_fines - 1) * α_sf_1), 0) ^ n_bf_1 / α_wf ^ 2 + v_wf), # 12
    0 ~ v_bf_2 - (g * (d_b[2] ^ 2) * (ρ_b[2] - ρ_f) / (18 * μ) * max((α_wf - (k_fines - 1) * α_sf_1), 0) ^ n_bf_2 / α_wf ^ 2 + v_wf), # 13
    0 ~ v_sf_1 - (g * (d_s[1] ^ 2) * (ρ_s - ρ_f) / (18 * μ) * max((α_wf - (k_fines - 1) * α_sf_1), 0) ^ n_sf_1 / α_wf ^ 2 + v_wf), # 14
    # MIDDLING UPPER
    0 ~ v_bmu_1 - (g * (d_b[1] ^ 2) * (ρ_b[1] - ρ_mu) / (18 * μ) * max((α_wmu - (k_fines - 1) * α_smu_1), 0) ^ n_bmu_1 / α_wmu ^ 2 + v_wmu), # 15
    0 ~ v_bmu_2 - (g * (d_b[2] ^ 2) * (ρ_b[2] - ρ_mu) / (18 * μ) * max((α_wmu - (k_fines - 1) * α_smu_1), 0) ^ n_bmu_2 / α_wmu ^ 2 + v_wmu), # 16
    0 ~ v_smu_1 - (g * (d_s[1] ^ 2) * (ρ_s - ρ_mu) / (18 * μ) * max((α_wmu - (k_fines - 1) * α_smu_1), 0) ^ n_smu_1 / α_wmu ^ 2 + v_wmu), # 17
    # Tailing
    0 ~ v_bt_1 - (g * (d_b[1] ^ 2) * (ρ_b[1] - ρ_t) / (18 * μ) * max((α_wt - (k_fines - 1) * α_st_1), 0) ^ n_bt_1 / α_wt ^ 2 + v_wt), # 18
    0 ~ v_st_1 - (g * (d_s[1] ^ 2) * (ρ_s - ρ_t) / (18 * μ) * max((α_wt - (k_fines - 1) * α_st_1), 0) ^ n_st_1 / α_wt ^ 2 + v_wt), # 19
    0 ~ v_st_2 - (g * (d_s[2] ^ 2) * (ρ_s - ρ_t) / (18 * μ) * max((α_wt - (k_fines - 1) * α_st_1), 0) ^ n_st_2 / α_wt ^ 2 + v_wt), # 20
    0 ~ v_st_3 - (g * (d_s[3] ^ 2) * (ρ_s - ρ_t) / (18 * μ) * max((α_wt - (k_fines - 1) * α_st_1), 0) ^ n_st_3 / α_wt ^ 2 + v_wt), # 21
    ########## VOLUME FRACTION CHANGE ##########
    # % Froth
    0 ~ α_bf_1 + (α_bmu_1 * A * v_bf_1) / (λ_bf * Qf), # 22
    0 ~ α_bf_2 + (α_bmu_2 * A * v_bf_2) / (λ_bf * Qf), # 23
    0 ~ α_sf_1 + (α_smu_1 * A * v_sf_1) / (λ_ff * Qf), # 24
    # % MIDDLING UPPER
    0 ~ α_bmu_1 - (α_bm_1 * A * v_bmu_1) / (A * v_bf_1), # 25
    0 ~ α_bmu_2 - (α_bm_2 * A * v_bmu_2) / (A * v_bf_2), # 26
    0 ~ α_smu_1 - (α_sm_1 * A * v_smu_1) / (A * v_sf_1), # 27
    # % MIDDLING (SOURCE ZONE)
    0 ~ (Qfd * α_bfd_1 + Qff * 0 - λ_bm * Qm * α_bm_1 - α_bm_1 * A_t * v_bt_1 + α_bm_1 * A * v_bmu_1), # 28 α_bff_1
    0 ~ (Qfd * α_bfd_2 + Qff * 0 - λ_bm * Qm * α_bm_2 - α_bm_2 * A_t * v_bt_2 + α_bm_2 * A * v_bmu_2), # 29 α_bff_2
    0 ~ (Qfd * α_sfd_1 + Qff * 0 - λ_bm * Qm * α_sm_1 - α_sm_1 * A_t * v_st_1 + α_sm_1 * A * v_smu_1), # 30 α_sff_1
    0 ~ (Qfd * α_sfd_2 + Qff * 0 - λ_bm * Qm * α_sm_2 - α_sm_2 * A_t * v_st_2 + α_sm_2 * A * v_smu_2), # 31 α_sff_2
    0 ~ (Qfd * α_sfd_3 + Qff * 0 - λ_bm * Qm * α_sm_3 - α_sm_3 * A_t * v_st_3 + α_sm_3 * A * v_smu_3), # 32 α_sff_3
    # TAILING
    0 ~ α_st_1 - (α_sm_1 * A_t * v_st_1) / (λ_ft * Qt), # 33
    0 ~ α_st_2 - (α_sm_2 * A_t * v_st_2) / (λ_st * Qt), # 34
    0 ~ α_st_3 - (α_sm_3 * A_t * v_st_3) / (λ_st * Qt), # 35
    0 ~ α_bt_1 - (α_bm_1 * A_t * v_bt_1) / (λ_bt * Qt), # 36
    ########### % COMBINED VOLUME FRACTION ##########
    0 ~ α_bf + (α_bmu_1 * A * v_bf_1 + α_bmu_2 * A * v_bf_2) / (λ_bf * Qf), # 37
    0 ~ α_sf - α_sf_1, # 38
    0 ~ α_bmu - (α_bmu_1 + α_bmu_2), # 39
    0 ~ α_smu - α_smu_1, # 40 
    0 ~ α_bm - (α_bm_1 + α_bm_2), # 41 
    0 ~ α_sm - (α_sm_1 + α_sm_2 + α_sm_3), # 42
    0 ~ α_bt - α_bt_1, # 43 
    0 ~ α_st - (α_st_1 + α_st_2 + α_st_3), # 44
]

psv_state = [v_bf_1, v_bf_2, v_sf_1, v_bmu_1, v_bmu_2, v_smu_1, v_bt_1, v_st_1, v_st_2, v_st_3, v_wf, v_wmu, v_wt,
             ρ_f, ρ_mu, ρ_m, ρ_t, α_bf, α_sf, α_bmu, α_smu, α_bm, α_sm, α_bt, α_st, α_bf_1, α_bf_2, α_sf_1,
             α_bmu_1, α_bmu_2, α_smu_1, α_bm_1, α_bm_2, α_sm_1, α_sm_2, α_sm_3, α_bt_1, α_st_1, α_st_2, α_st_3, α_wf, α_wmu, α_wm, α_wt]

psv_param = [A, A_t, μ, g, ρ_s, ρ_w, d_b_1, d_b_2, d_s_1, d_s_2, d_s_3, k_fines, n_bf_1, n_bf_2, n_sf_1, n_bmu_1, n_bmu_2, n_smu_1, n_bt_1, n_st_1, n_st_2, n_st_3, λ_bf, λ_ff, λ_bm, λ_sm, λ_bt, λ_st, λ_fm, λ_ft]

@named psv_sys = NonlinearSystem(psv_eqs, psv_state, psv_param)

simplified_psv_sys = structural_simplify(psv_sys,; verbose=true)

But the number of equations and variables are equal in practice and everyting is defined properly, I guess.
What can I do to resolve the issue?
modelingtoolkit
nonlinearsolve