# Request help solving a Mathematical Programming with Equilibrium Constraints problem in JuMP?

**URL:** https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802
**Category:** Optimization (Mathematical)
**Created:** [June 21, 2020, 1:53am UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802 "2020-06-21T01:53:40Z")
**Posts on this page:** 1
**Showing post:** 1

<div class="post-metadata">

### Author: ![jack\_rabbit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jack_rabbit/32/10914_2.png) [@jack\_rabbit](https://discourse.julialang.org/u/jack_rabbit)
#### Post date: [June 21, 2020, 1:53am UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/1 "2020-06-21T01:53:40Z")

</div>

Hi,

I am pretty new to both Julia and JuMP. Please feel free to point out any aspects that I didn’t understand well.  
I am trying to solve a Mathematical Programming with Equilibrium Constraints (MPEC) problem using the interface of JuMP. My problem can be formulated as

\max\_{\theta, EV} \mathcal{L}(\theta, EV) \\ s.t. EV = T(EV, \theta)

where \theta and EV are both vectors. And T is a contracting mapping that for a given \theta, there should exist a vector EV (a fixed-point in the jargon) satisfies this constraint. This problem is quite similar to a problem [here](https://discourse.julialang.org/t/a-hacky-guide-to-using-automatic-differentiation-in-nested-optimization-problems/39123/5). So I guess it is of great possibility that I may get help here.  
The following steps document how I was trying to solve this problem in JuMP. Since I have a user-defined likelihood function to be maximized, after reading [the relevant document](https://jump.dev/JuMP.jl/stable/nlp/#User-defined-Functions-1) for JuMP, I first wrapped a likelihood function like this

```nohighlight
const S = 90
# wrap a likelihood function for Optimizer
function wrap_ll(Θ...)
    vΘ = collect(Θ)
    θ₁ = vΘ[1]
    θ₃ = vΘ[2:4]
    RC = vΘ[5]
    EV = reshape(vΘ[6:end], 2, S)

    θ = BusParameters(θ₁=θ₁, θ₃=θ₃, RC=RC)
    return loglikelihood(xmt, dmt, θ, EV)
end

```

where `BusParameters` is a self-defined struct to hold parameters \theta in the problem and `xmt` and `dmt` are my data.  
Then I define my model in JuMP using the Ipopt solver

```nohighlight
model = Model(with_optimizer(Ipopt.Optimizer, max_cpu_time=300.0))
JuMP.register(model, :wrap_ll, 5+2*S, wrap_ll, autodiff=true)

```

Then add variables with start values

```nohighlight
@variable(model, Θ[1:5+2*S])
# set start value for all parameters
set_start_value(Θ[1], 2.2923)
set_start_value(Θ[2], 0.3010)
set_start_value(Θ[3], 0.6884)
set_start_value(Θ[4], 0.0106)
set_start_value(Θ[5], 10.07)
for i in 6:5+2*S
    set_start_value(Θ[i], 0.0)
end

```

And add some constraints for my variables

```nohighlight
# add constraints
@constraint(model, Θ[1]>=0)
@constraint(model, 0.00001 .<= Θ[2:4] .<= 1)
@constraint(model, Θ[6:end] .>= 0)
@constraint(model, sum(Θ[2:4]) == 1)

```

Up to this point, I think it’s quite smooth.  
The first difficulty I encounter is how can I add the Equilibrium Constraint EV = T(EV, \theta). I have a user-defined function T and after check [this question](https://discourse.julialang.org/t/how-to-add-user-defined-nonlinear-functions-as-constraints-with-jump/11436) in the forum, I wrapped a T function like this

```nohighlight
function constraint_T(Θ...)
    vΘ = collect(Θ)
    θ₁ = vΘ[1]
    θ₃ = vΘ[2:4]
    RC = vΘ[5]
    EV = reshape(vΘ[6:end], 2, S)

    θ = BusParameters(θ₁=θ₁, θ₃=θ₃, RC=RC)
    EVₙ = T(EV, θ)
    return maximum(abs.(EVₙ .- EV))
end

```

Here the `T` function returns a new `EV`, I was trying to return `EVₙ .- EV` instead of `maximum(abs.(EVₙ .- EV))` and add a constraint like this

```nohighlight
JuMP.register(model, :constraint_T, 5+2*S, constraint_T, autodiff=true )
@NLconstraint(model, constraint_T(Θ...) .== 0)

```

But this broadcasting form is not allowed in JuMP, so I used a trick that the `constraint_T` function return `maximum(abs.(EVₙ .- EV))` and add the constraint like this instead

```nohighlight
@NLconstraint(model, constraint_T(Θ...) == 0)

```

I think, conceptually, `maximum(abs.(EVₙ .- EV)) == 0` is equivalent to `EVₙ .- EV .== 0`. I don’t know if this has any problem with the solver. I know JuMP has a way to define [complementarity constraints](https://jump.dev/JuMP.jl/v0.21.1/constraints/index.html#Complementarity-constraints-1), but I am not sure how can I do it that way in my example.  
And finally, define my objective function and try to solve the model

```nohighlight
@NLobjective(model, Max, wrap_ll(Θ...))
optimize!(model)

```

Formulate my problem this way seems to work without errors. However, after running the solver for a certain period of time, the results I got seems not quite right (I have solved this problem using [Knitro solver](https://neos-server.org/neos/solvers/cp:Knitro/AMPL.html) in NEOS).  
So I wonder whether my way of defining the problem in JuMP is correct and `Ipopt` is a right solver for my problem and what may be the alternatives? (I plan to try [Knitro.jl](https://github.com/JuliaOpt/KNITRO.jl) later.) Thanks for any comments or suggestions.

---

_[View the full topic](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802)._