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.