Solving dynamic complementarity problem using Optim.jl or Complementarity.jl

Specific Domains Optimization (Mathematical)
, ,
1

I am interested in solving a dynamic complementarity problem using Complementarity.jl or Optim.jl. In particular I would appreciate some assistance in solving or how I could solve the following.

A - B*q[t] = λ[t]            if q[t] > 0
λ[t]*(1 + δ) = λ[t+1]        if R[t+1] > 0
R[t] = R[t-1] - q[t]         if λ[t] > 0

where t in [0, T] and R[0] is given.

using Complementarity.jl, I provide below code that can be extended. I would like to add to this code, the dynamic aspects but so far my attempts have fallen flat (yet I can easily solve the dynamic model in GAMS).

using PATH, JuMP, Complementarity

function geological_constraints_MCP()

    Α = 50.0
    Β = 0.1
    δ = 0.05
    R₀ = 10.0
    T = 4
    t = 1:T

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

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

geological_constraints_MCP()
2

Let me provide a dynamic extension of the above code—but unfortunately it does not return a solution. Hopefully I get some ideas:

using PATH, JuMP, Complementarity

function geological_constraints_MCP()

    Α = 50.0
    Β = 0.1
    δ = 0.05
    R₀ = 0.5
    T = 10
    t = 1:T

    gxt =MCPModel()

    @variable(gxt, q[i in 1:T]>=0)
    @variable(gxt, λ[i in 1:T]>=0)
    @variable(gxt, R[i in 1:T]>=0)
    #@variable(gxt, R[i in 1:T])

    @mapping(gxt, FOCq[i in t], - (Α - Β*q[i] - λ[i]))
    @complementarity(gxt,FOCq,q)
    @mapping(gxt, FOCR[i in 1:T], - ((i<T ? λ[i+1] : Α) - λ[i]*(1 + δ)))
    @complementarity(gxt,FOCR,λ)
    @mapping(gxt, FOCλ[i in 1:T], - (R[i] + q[i] - (i>1 ? R[i-1] :  R₀)))
    @complementarity(gxt,FOCλ,R)
    #@constraint(gxt, [i in 1:T], (R[i] + q[i] - (i>1 ? R[i-1] :  R₀)) == 0)

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


geological_constraints_MCP()
3

You’re looking for something like

using JuMP
import PATH

function geological_constraints_MCP()
    A = 50.0
    B = 0.1
    δ = 0.05
    R₀ = 0.5
    T = 10
    model = Model(PATH.Optimizer)
    @variables(model, begin
        q[1:T] >= 0
        λ[1:T] >= 0
        R[1:T] >= 0
    end)
    @constraints(model, begin
        [t = 1:T], -(A - B * q[t] - λ[t])                     ⟂ q[t] 
        [t = 1:T], -((t < T ? λ[t + 1] : A) - (1 + δ) * λ[t]) ⟂ R[t]
        [t = 1:T], -(R[t] + q[t] - (t > 1 ? R[t - 1] :  R₀))  ⟂ λ[t]
    end)
    optimize!(model)
    @show termination_status(model)
    @show primal_status(model)
    @show value.(q)
    @show value.(λ)
    @show value.(R)
    return model
end

geological_constraints_MCP()

I don’t fully understand your model. For example, @complementarity(gxt,FOCλ,R) seemed like it should be @complementarity(gxt,FOCλ,λ) based on your original post?

Nevertheless, hope the syntax helps.

1 Like
4

Hi @odow . You deserve an Empire state building size like. Thanks for spotting @complementarity(gxt,FOCλ,R). This is where my frustrations started … and have apparently ended.

1 Like