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

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. 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 for JuMP, I first wrapped a likelihood function like this

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

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

```
@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

```
# 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 in the forum, I wrapped a T function like this

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

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

```
@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, 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

```
@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 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 later.) Thanks for any comments or suggestions.