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()