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

Domains Optimization (Mathematical)
#1

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

1 Like
#2

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?

#3

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?

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?

#4

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.

1 Like
#5

I see. Thanks for the note. I was not accustomed to thinking about the underlying math of programs.

#6

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

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

JuMP.register(model, :constraint_T, 1+5+2*S, constraint_T, autodiff=true )

And add the separate constraints

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?

#8

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

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.]