Julia JuMP optimization help

To partly answer my own question. See this discussion.

Using “complementarity.jl” works with nonlinear expressions. Here is the working piece of code:

using PATH, JuMP, Complementarity

function geological_constraints_MCP()

    Α = 50.0
    Β = 0.1
    δ = 0.05
    γ = 0.1
    ψ₁ = 0.0
    ψ₂ = 2.0
    λ₁ = 5.0
    λ₂ = 4.0
    R₀ = 0.0
    S₀ = 10.0
    T = 4
    t = 1:T

    gxt =MCPModel()
    @variable(gxt, q[i in t]>=0)
    @mapping(gxt, FOC[i in t], -(Α - Β*q[i] - ψ₁*q[i]^(ψ₂-1)))
    @complementarity(gxt,FOC,q)

    print(gxt)
    status = solveMCP(gxt)
    #value.(q)
    @show result_value.(q)
end


geological_constraints_MCP()