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

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

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

``````

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

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