Redundant equations causes a crash

[azev77]

The following nonlinear system causes a crash:

using ModelingToolkit, NonlinearSolve, Plots
if 1==1
      @variables LD, KD, Y, Π,   c, ℓ, LS, KS, v,   r, w
      @parameters αK, αL, γ, EH, EK     

      F(K,L,αK, αL) = (K^αK)*(L^αL)
      FK(K,L,αK,αL) = αK*(K^αK)*(L^αL)/K
      FL(K,L,αK,αL) = αL*(K^αK)*(L^αL)/L

      U(c,ℓ,γ) = γ*log(c) + (1.0 -γ)*log(ℓ)
      Uc(c,ℓ,γ) = γ/c 
      Uℓ(c,ℓ,γ) = (1.0 -γ)/ℓ

      eqs = [
            0.0 ~ r - FK(KD,LD,αK,αL),                    # r = F_K
            0.0 ~ w - FL(KD,LD,αK,αL),                    # w = F_L
            0.0 ~ Y - F(KD,LD,αK, αL),                    # Y = F(K,L)
            0.0 ~ Π - (Y- r*KD - w*LD),                   # Π = Y -rK -wL
            #
            0.0 ~ Uc(c,ℓ,γ)/Uℓ(c,ℓ,γ) - 1.0/w,            # Uc/Uℓ = 1/w   MRS = price ratio
            0.0 ~ c + w*ℓ - (w*EH + r*KS + Π),            # BC
            0.0 ~ EH - LS - ℓ,                            # ℓ = EH - L
            0.0 ~ EK - KS,                                # K ≤ EK  
            0.0 ~ v -  U(c,ℓ,γ),                          # value function 
            #
            0.0 ~ KD - KS,                                # Capital Market Clearing 
            0.0 ~ LD - LS,                                # Labor Market Clearing 
            0.0 ~ c - Y,                                  # Goods Market Clearing 
      ]
      ns = NonlinearSystem(eqs, 
            [LD, KD, Y, Π,   c, ℓ, LS, KS, v,   r, w,], 
            [αK, αL, γ, EH, EK,]
            )
      guess = [LD=>1.0, KD=>1.0, Y=>1.0, Π=>1.0,
            c=>1.0, ℓ=>1.0, LS=>1.0, KS=>1.0, v=>1.0,
            r=>1.0, w=>1.0]
      ps = [ αK=>0.45; αL=>0.51; γ=>0.5; EH=>24.0; EK=>10.0;]
      #ps = [ αK=>0.5; αL=>0.5; γ=>0.5; EH=>24.0; EK=>10.0;]
      prob = NonlinearProblem(ns,guess,ps)
      sol = solve(prob,NewtonRaphson())
      tuple(sol...)
end 

Message:
The terminal process "C:\Users\azevelev\AppData\Local\Programs\Julia-1.6.0\bin\julia.exe '-i', '--banner=no', '--project=C:\Users\azevelev.julia\environments\v1.6', 'c:\Users\azevelev.vscode\extensions\julialang.language-julia-1.1.37\scripts\terminalserver\terminalserver.jl', '\.\pipe\vsc-jl-repl-04d88e42-6073-49e7-a1c7-7d343e77fe5a', '\.\pipe\vsc-jl-cr-1ed06595-e8e0-49c4-9402-f2ce31c35611', 'USE_REVISE=true', 'USE_PLOTPANE=true', 'USE_PROGRESS=true', 'DEBUG_MODE=undefined'" terminated with exit code: 1.

It works if we comment out one (and only one) of the redundant market clearing equations:

using ModelingToolkit, NonlinearSolve, Plots
if 1==1
      @variables LD, KD, Y, Π,   c, ℓ, LS, KS, v,   r, w
      @parameters αK, αL, γ, EH, EK     

      F(K,L,αK, αL) = (K^αK)*(L^αL)
      FK(K,L,αK,αL) = αK*(K^αK)*(L^αL)/K
      FL(K,L,αK,αL) = αL*(K^αK)*(L^αL)/L

      U(c,ℓ,γ) = γ*log(c) + (1.0 -γ)*log(ℓ)
      Uc(c,ℓ,γ) = γ/c 
      Uℓ(c,ℓ,γ) = (1.0 -γ)/ℓ

      eqs = [
            0.0 ~ r - FK(KD,LD,αK,αL),                    # r = F_K
            0.0 ~ w - FL(KD,LD,αK,αL),                    # w = F_L
            0.0 ~ Y - F(KD,LD,αK, αL),                    # Y = F(K,L)
            0.0 ~ Π - (Y- r*KD - w*LD),                   # Π = Y -rK -wL
            #
            0.0 ~ Uc(c,ℓ,γ)/Uℓ(c,ℓ,γ) - 1.0/w,            # Uc/Uℓ = 1/w   MRS = price ratio
            0.0 ~ c + w*ℓ - (w*EH + r*KS + Π),            # BC
            0.0 ~ EH - LS - ℓ,                            # ℓ = EH - L
            0.0 ~ EK - KS,                                # K ≤ EK  
            0.0 ~ v -  U(c,ℓ,γ),                          # value function 
            #
            #0.0 ~ KD - KS,                                # Capital Market Clearing 
            0.0 ~ LD - LS,                                # Labor Market Clearing 
            0.0 ~ c - Y,                                  # Goods Market Clearing 
      ]
      ns = NonlinearSystem(eqs, 
            [LD, KD, Y, Π,   c, ℓ, LS, KS, v,   r, w,], 
            [αK, αL, γ, EH, EK,]
            )
      guess = [LD=>1.0, KD=>1.0, Y=>1.0, Π=>1.0,
            c=>1.0, ℓ=>1.0, LS=>1.0, KS=>1.0, v=>1.0,
            r=>1.0, w=>1.0]
      ps = [ αK=>0.45; αL=>0.51; γ=>0.5; EH=>24.0; EK=>10.0;]
      #ps = [ αK=>0.5; αL=>0.5; γ=>0.5; EH=>24.0; EK=>10.0;]
      prob = NonlinearProblem(ns,guess,ps)
      sol = solve(prob,NewtonRaphson())
      tuple(sol...)
end 

[YingboMa]

I cannot reproduce the error.

[YingboMa]

Ah, I see. You have 11 states and 12 equations. That's not a valid NonlinearSystem. The number of unknowns and the number of equations must match.

[ChrisRackauckas]

Are you adding consistency_check and a nice error message to the SciMLProblem constructors?

[ChrisRackauckas]

This has more states than equations, but only because one is redundant. @azev77 did you try prob = NonlinearProblem(structural_simplify(ns),guess,ps) to see if it was eliminated?

[azev77]

1. if I rerun the original problem (12 equations/11 variables):
   -It works the first time I run it.
   -It crashes the second time.
   -If I eliminate any one (and only one) of the three redundant eqs (11 eqs/11 vars) it always works!
2. structural_simplify(ns) gives this error:
   julia> prob = NonlinearProblem(structural_simplify(ns),guess,ps)
   ERROR: BoundsError: attempt to access 11-element BitVector at index [12]

[YingboMa]

Yeah, we need better checks.

[azev77]

BTW, including a redundant linear equation does not cause a crash:

@variables  x, y 
@parameters θ

#eqs   = [ 0.0 ~ 2.0x - 15.2, 0.0 ~ x + y - 5.2,]
eqs   = [ 0.0 ~ 2.0x - 15.2, 0.0 ~ x + y - 5.2, 0.0 ~ 2.0x + 2.0y - 2.0*5.2,]

ns    = NonlinearSystem(eqs, [x, y], [θ,])
guess = [x=>1.0, y=>1.0]
ps    = [ θ=>0.45; ]
prob  = NonlinearProblem(ns,guess,ps)
sol   = solve(prob,NewtonRaphson())
tuple(sol...)

[YingboMa]

Nah, that's a fluke.

[azev77]

Here is my minimal example which causes a crash:

using ModelingToolkit, NonlinearSolve
@variables  x, y 
@parameters θ
eqs   = [ 0.0 ~ sin(x) - x, 0.0 ~ log(x) + y - 5.2, 0.0 ~ 2.0log(x) + 2.0y - 2.0*5.2,]
ns    = NonlinearSystem(eqs, [x, y], [θ,])
guess = [x=>1.0, y=>1.0]
ps    = [ θ=>0.45; ]
prob  = NonlinearProblem(ns,guess,ps)
sol   = solve(prob,NewtonRaphson())

[azev77]

Still doesn't work on my machine:

julia> using ModelingToolkit, NonlinearSolve
[ Info: Precompiling ModelingToolkit [961ee093-0014-501f-94e3-6117800e7a78]

julia> @variables  x, y 
(x, y)

julia> @parameters θ
(θ,)

julia> eqs   = [ 0.0 ~ sin(x) - x, 0.0 ~ log(x) + y - 5.2, 0.0 ~ 2.0log(x) + 2.0y - 2.0*5.2,]
3-element Vector{Equation}:
 0.0 ~ sin(x) - x
 0.0 ~ y + log(x) - 5.2
 0.0 ~ 2.0y + 2.0log(x) - 10.4

julia> ns    = NonlinearSystem(eqs, [x, y], [θ,])
Model ##NonlinearSystem#257 with 3 equations
States (2):
  x
  y
Parameters (1):
  θ

julia> guess = [x=>1.0, y=>1.0]
2-element Vector{Pair{Num, Float64}}:
 x => 1.0
 y => 1.0

julia> ps    = [ θ=>0.45; ]
1-element Vector{Pair{Num, Float64}}:
 θ => 0.45

julia> prob  = NonlinearProblem(ns,guess,ps)
ERROR: ArgumentError: Equations (3), states (2), and initial conditions (2) are of different lengths.
Stacktrace:
 [1] check_eqs_u0(eqs::Vector{Equation}, dvs::Vector{Sym{Real, Nothing}}, u0::Vector{Float64})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/abstractsystem.jl:573
 [2] process_NonlinearProblem(constructor::Type, sys::NonlinearSystem, u0map::Vector{Pair{Num, Float64}}, parammap::Vector{Pair{Num, Float64}}; version::Nothing, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/nonlinear/nonlinearsystem.jl:223
 [3] process_NonlinearProblem
   @ ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/nonlinear/nonlinearsystem.jl:216 [inlined]
 [4] (NonlinearProblem{true, isinplace, P, F, K} where {isinplace, P, F, K})(sys::NonlinearSystem, u0map::Vector{Pair{Num, Float64}}, parammap::Vector{Pair{Num, Float64}}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/nonlinear/nonlinearsystem.jl:250
 [5] (NonlinearProblem{true, isinplace, P, F, K} where {isinplace, P, F, K})(sys::NonlinearSystem, u0map::Vector{Pair{Num, Float64}}, parammap::Vector{Pair{Num, Float64}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/nonlinear/nonlinearsystem.jl:250
 [6] NonlinearProblem(::NonlinearSystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Vector{Pair{Num, Float64}}, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/nonlinear/nonlinearsystem.jl:232
 [7] NonlinearProblem(::NonlinearSystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Vector{Pair{Num, Float64}}, N} where N)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/5x2CA/src/systems/nonlinear/nonlinearsystem.jl:232
 [8] top-level scope
   @ REPL[8]:1

[YingboMa]

Exactly. It shouldn't work