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

[azev77]

Why shouldn’t it work?

redundant equations usually cause no problem in these types of packages

[YingboMa]

For what types of packages? Could you give me an example which this works? It shouldn't work because it is not a nonlinear system.

[azev77]

I can take a look tomorrow, but this should be true in general.

People overconstrain problems for numerical reasons.

You need at least as many equations as unknowns. You can have redundant equations as long as they are truly redundant restrictions (otherwise the problem may have no solution)...

[YingboMa]

No, that's not for numerical reasons. You probably just want to solve an optimization problem instead.

[YingboMa]

Alright, let's take a look at your system.

sin(x) - x ~ 0 => x = 0. 0.0 ~ y + log(x) - 5.2 is not well defined because of the log(0).

[YingboMa]

If you run structural_simplify on your original system, you will get

julia> sys = structural_simplify(ns);

julia> prob = NonlinearProblem(ns,guess,ps);

julia> sol = solve(prob,NewtonRaphson())
u: 11-element Vector{Float64}:
  8.105960280618913
  9.99999994216286
  8.193892851162932
  0.32775572634556627
  8.193892851162932
 15.894039719381087
  8.105960280618913
 10.0
  2.4346666406933113
  0.36872517697473867
  0.5155324225949323

which just works fine.

[YingboMa]

Ah, I don't know what happened last night. I must have copied the wrong equations. When I tried the first system in the morning, I got

julia> structural_simplify(ns)
ERROR: InvalidSystemException: The system is unbalanced. There are 8 highest order derivative variables and 9 equations.
Stacktrace:
 [1] check_consistency(s::SystemStructure)
   @ ModelingToolkit.StructuralTransformations ~/src/julia/ModelingToolkit/src/structural_transformation/utils.jl:59
 [2] structural_simplify(sys::NonlinearSystem)
   @ ModelingToolkit ~/src/julia/ModelingToolkit/src/systems/abstractsystem.jl:545
 [3] top-level scope
   @ REPL[55]:1

[azev77]

Sorry, I meant cos(x) = x

#Ex. 
using ModelingToolkit, NonlinearSolve
@variables  x, y 
@parameters θ
eqs   = [ 0.0 ~ cos(x) - x, 0.0 ~ log(x) + y - 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())
0.739085133385284
5.502342163171829

#Ex. 
using ModelingToolkit, NonlinearSolve
@variables  x, y 
@parameters θ
eqs   = [ 0.0 ~ cos(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]

Notice: eqs = [ 0.0 ~ cos(x) - x, 0.0 ~ log(x) + y - 5.2,] and eqs = [ 0.0 ~ cos(x) - x, 0.0 ~ log(x) + y - 5.2, 0.0 ~ 2.0log(x) + 2.0y - 2.0*5.2,] has the same set of solutions.

However Julia crashes when solving the second problem.