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Θ θ₃ = vΘ[2:4] RC = vΘ EV = reshape(vΘ[6:end], 2, S) θ = BusParameters(θ₁=θ₁, θ₃=θ₃, RC=RC) return loglikelihood(xmt, dmt, θ, EV) end
BusParameters is a self-defined struct to hold parameters \theta in the problem 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(Θ, 2.2923) set_start_value(Θ, 0.3010) set_start_value(Θ, 0.6884) set_start_value(Θ, 0.0106) set_start_value(Θ, 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, Θ>=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Θ θ₃ = vΘ[2:4] RC = vΘ EV = reshape(vΘ[6:end], 2, S) θ = BusParameters(θ₁=θ₁, θ₃=θ₃, RC=RC) EVₙ = T(EV, θ) return maximum(abs.(EVₙ .- EV)) end
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.