Modeling Toolkit crash

Specific Domains Modelling & Simulations
1

The following causes a crash in Julia v 1.6

using ModelingToolkit, NonlinearSolve, Plots

@variables c, ℓ, L, K, Y, Π, 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-1.0))*(L^αL)
FL(K,L,αK,αL) = αL*(K^αK)*(L^(αL-1.0))

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

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

Crash occurs at solve()

Crash 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-95aefdc3-19e2-4af6-b729-7a88dcfed920’, ‘\.\pipe\vsc-jl-cr-0dc6ed39-f5cb-40b5-bbb6-83e0abb7fcb3’, ‘USE_REVISE=true’, ‘USE_PLOTPANE=true’, ‘USE_PROGRESS=true’, ‘DEBUG_MODE=undefined’”
terminated with exit code: 1.

2

Open an issue.

3

Thanks I’ll open an issue.
I think I got to the bottom of it: it crashes whenever I post a redundant equation

Basically, I have 3 market clearing equations, and only need two (any two) of the three to pin down the r the price of capital, and w the price of labor. It’s a shame that over-constraining a non-linear system w/ redundant equations causes a crash. I’m sure it’s a fixable bug…

For whoever wants to see, the following code works if you comment out any one (and only one) of the three 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
4

@ChrisRackauckas is there a good rule of thumb for when I should open an issue in the repo vs when I should post my question on Discourse?

5

When is a bug report and not a discussion. Sometimes it can be hard to tell if your not a dev so don’t worry too much

2 Likes