Solving Mixed Complementarity Problem in Julia

Question

I am new to Julia. Trying to implement a toy economic equilibrium model as an MCP using Julia, JuMP and PATHSolver.

The problem assumes two regions trading a single commodity. Each region has linear demand and a fixed exogenous supply, A. Problem is to find the market clearing equilibrium prices P, subject to trade upper and lower bounds. Here is the code:

using JuMP, PATHSolver

β_0 = [10.0, 5.0]
β_1 = [-0.5, -0.25]
t_min = 0.0
t_max = 2.0
A = [1, 1]

model = Model(PATHSolver.Optimizer)

@variables(model, begin
    t_min <= trade <= t_max
    p1 >= 0
    p2 >= 0
end)

@constraints(model, begin
    p2 - p1 ⟂ trade 
    β_0[1] + β_1[1] * p1 - A[1] + trade ⟂ p1
    β_0[2] + β_1[2] * p2 - A[2] - trade ⟂ p2
end)

optimize!(model)

solution_summary(model; verbose=true)

Here is the output:

* Solver : Path 5.0.03

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "The problem was solved"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 0.00000e+00
  Primal solution :
    p1 : 0.00000e+00
    p2 : 0.00000e+00
    trade : 0.00000e+00

* Work counters
  Solve time (sec)   : 2.99100e-03

Answer 1

Was there a particular question? It looks like that solution is a valid equilibrium.

You probably want to set a different starting solution:

@variables(model, begin
    t_min <= trade <= t_max, (start = 1)
    p1 >= 0, (start = 1)
    p2 >= 0, (start = 1)
end)

Answer 2

Changed the order of the constraints (excess supply rather than excess demand) and this seems to work now:

using JuMP, PATHSolver

β_0 = [10.0, 5.0]
β_1 = [-0.5, -0.25]
t_min = 0.0
t_max = 2.0
A = [0.5, 1]

model = Model(PATHSolver.Optimizer)

@variables(model, begin
    t_min <= trade <= t_max 
    p1 >= 0, (start = 1)
    p2 >= 0, (start = 1)
end)

@constraints(model, begin
    p2 - p1 ⟂ trade 
    A[1] + trade - (β_0[1] + β_1[1] * p1) ⟂ p1
    A[2] - trade - (β_0[2] + β_1[2] * p2) ⟂ p2
end)

optimize!(model)

solution_summary(model; verbose=true)

* Solver : Path 5.0.03

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "The problem was solved"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 0.00000e+00
  Primal solution :
    p1 : 1.80000e+01
    p2 : 1.80000e+01
    trade : 5.00000e-01

* Work counters
  Solve time (sec)   : 2.53900e-03