# 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:** 11
**Page:** 1
### Author:  [@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")
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [June 21, 2020, 2:58am UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/2 "2020-06-21T02:58:08Z")
Unfortunately, JuMP doesn’t support vector nonlinear functions at present. Correcting this is on our (long) list of things to do.
How complicated is the function `T`? Can you write out each component algebraically?
---
### Author:  [@jack\_rabbit](https://discourse.julialang.org/u/jack_rabbit)
#### Post date: [June 21, 2020, 2:20pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/3 "2020-06-21T14:20:01Z")
Thanks for your reply @odow.
> Unfortunately, JuMP doesn’t support vector nonlinear functions at present. Correcting this is on our (long) list of things to do.
What if I write the vector nonlinear function as a scalar function and use a loop to add a constraint for each EV[i], like in [this example](https://discourse.julialang.org/t/jump-rewriting-an-eval-nlexpression-with-a-user-defined-function-in-local-scope-w-expr/8757/5)?
> How complicated is the function `T` ? Can you write out each component algebraically?
Right now the `T` function is not very complex but also not trivial, I guess I may write out each component algebraically but the expression would be relatively long. The problem I would imagine is that in the future my problem will be more complex than this and to work out an algebra expression be hardly possible. What’s the question you had in mind? Do you think that my adding of the constraint using `maximum(abs.(EVₙ .- EV)) == 0` be a problem?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [June 21, 2020, 2:59pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/4 "2020-06-21T14:59:52Z")
Yes, writing out each component as a separate constraint will work.
> Do you think that my adding of the constraint using `maximum(abs.(EVₙ .- EV)) == 0` be a problem
Note that Ipopt assumes certain properties of the nonlinear functions (smooth, twice differentiable, etc.). It will still run, but you may have convergence issues, etc.
---
### Author:  [@jack\_rabbit](https://discourse.julialang.org/u/jack_rabbit)
#### Post date: [June 21, 2020, 6:31pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/5 "2020-06-21T18:31:33Z")
I see. Thanks for the note. I was not accustomed to thinking about the underlying math of programs.
---
### Author:  [@jack\_rabbit](https://discourse.julialang.org/u/jack_rabbit)
#### Post date: [June 21, 2020, 7:28pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/6 "2020-06-21T19:28:02Z")
When trying to write out each component as a separate constraint, I need to add an index `i` to the original `constraint_T` function like this
```nohighlight
function constraint_T(i, Θ...)
vΘ = collect(Θ)
θ₁ = vΘ[1]
θ₃ = vΘ[2:4]
RC = vΘ[5]
EV = reshape(vΘ[6:end], 2, S)
θ = BusParameters(θ₁=θ₁, θ₃=θ₃, RC=RC)
EVᵢ = T(EV, θ, i)
return EVᵢ - EV[i]
end
```
Then register the function with one additional parameter
```nohighlight
JuMP.register(model, :constraint_T, 1+5+2*S, constraint_T, autodiff=true )
```
And add the separate constraints
```nohighlight
for i in 1:2*S
@NLconstraint(model, constraint_T(i, Θ...) == 0)
end
```
But this introduces the additional parameter `i` that Autodiff is trying to differentiate, which of course will cause problems. How should let Autodiff skops to differentiate with `i`? Any suggestions @odow?
---
### Author:  [@PaoloBova](https://discourse.julialang.org/u/PaoloBova)
#### Post date: [December 2, 2020, 4:00pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/8 "2020-12-02T16:00:05Z")
Hi,
An alternative would be to define, register, and add each constraint in your for loop. This avoids needing to autodiff with respect to parameter i.
Solution:
for i in 1:number\_of\_constraints:
(i) Define your scalar function. Your function does not need to accept i as an argument.
(ii) Register your new function. Use Symbol() to create a unique keyword for each i.
(iii) Use JuMP.add\_NL\_constraint() and use [raw expression input](https://jump.dev/JuMP.jl/stable/nlp/#Raw-expression-input-1) to define your constraint. Note that we need Symbol() again to refer to the function created in (ii).
The relevant code in your case would be:
```julia
for i in 1:2*S
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, θ, i)
return EVᵢ - EV[i]
end
register(model, Symbol("constraint_$i"), 5+2*S, constraint_T, autodiff=true)
add_NL_constraint(model, :($(Expr(:call, Symbol("constraint_$i"), Θ...)) == 0))
end
```
Note: I do not know what JuMP.register does behind the scenes. If JuMP or Julia can already handle the extra parameter ‘i’ intelligently when performing AutoDiff, then there might be no advantage to my solution above. In that case, I’d recommend sticking to your current solution, so that you can use the JuMP macros. Can anyone confirm whether that extra parameter would trip up AutoDiff in Julia as jack\_rabbit suspects? [Quick Update: does not appear that JuMP can handle integer parameters like ‘i’ well at all. Happy to be proven wrong though.]
---
### Author:  [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)
#### Post date: [April 22, 2024, 5:28pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/9 "2024-04-22T17:28:17Z")
Hi everyone. Sorry for being late to the party, but is it solved now? And if yes: @jack_rabbit could you provide a full solution?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [April 24, 2024, 4:56am UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/10 "2024-04-24T04:56:24Z")
@CeterisPartybus can you start a new question with the problem you’re facing?
JuMP still doesn’t support vector-valued user-defined functions, but if your functions are simple enough, you can use function tracing: [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Function-tracing)
---
### Author:  [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)
#### Post date: [April 24, 2024, 12:53pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/11 "2024-04-24T12:53:54Z")
Thanks for the quick reply. I just think about doing exactly the same as @jack_rabbit. Thanks for the resource. I will have a look at it.
Is there any wish list for JuMP? It would probably be a simple but useful feature to add.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [April 24, 2024, 7:58pm UTC](https://discourse.julialang.org/t/request-help-solving-a-mathematical-programming-with-equilibrium-constraints-problem-in-jump/41802/12 "2024-04-24T19:58:36Z")
See [Roadmap · JuMP](https://jump.dev/JuMP.jl/stable/developers/roadmap/)
If you still need help, please start a new thread with a reproducible example of your problem.