Bilevel JuMP

You cannot have Integer (nor binary) variables in the lower model.
You have them in the upper model though.

2 Likes

You’re right. It’s OK.
In this way, the constraint of decision variables originally belonging to the lower model is put into the upper model, which is equivalent to writing it in the lower model?

No it is no equivalent.

The problem is:
The solution technique implemented in BilevelJuMP is the KKT method.
This method requires the lower level problem to be convex.

If the lower level problem has integer variables you need other methods, almost always, these other methods are either problem specific or do not guarantee global optimality.

You can add integrality in the upper level and create an equality with a variable in the lower level, however, this is NOT equivalent to solving a bilevel problem with the integrality constraints in the lower level. This might or might not not be a good approximation, will depend a lot on the problem and for some problems it is actually completely wrong.

2 Likes

I understand your meaning roughly as follows:


@variable(Lower(model), PES_BT_D[1:T,1:3])
@variable(Lower(model), PES_BT_C[1:T,1:3])
@variable(Upper(model), 0<=W1[1:T,1:3]<=1 ,Int)
@variable(Upper(model), 0<=W2[1:T,1:3]<=1 ,Int)
@constraints(Upper(model), begin
0 .<=PES_BT_D.*W1
      PES_BT_D.*W1 .<=40
  
      0 .<=PES_BT_C.*W2
      PES_BT_C.*W2 .<=80
      W1+W2 .<=1

I tried to use Matlab to push KKT by myself. After transforming it into a single-layer problem, I transformed the complementary constraint into a mixed integer problem by using the large M method. Unfortunately, it didn’t work out. Later, I found that Julia had your package, which filled me with hope

I dont get where you are going.
But I can say that the problem you just wrote wont work because of the product of variables.

1 Like

Yes, this is equivalent to adding complementary constraints. Model nonconvex

Would you mind writing with latex formatting a small version of the problem you want to solve?
In BilevelJuMP you should not write the complementarity constraints. They are written automatically. You just write the upper level primal and the lower level primal.

1 Like

Sorry, I don’t have latex installed. I only have MathType. I wrote it in word and found that it can’t be pasted here. Can you give me an email? Or I’ll send it here in the form of a screenshot

For the energy storage device, I don’t want it to be both charged and discharged in a certain period of time, so I added a 0-1 variable.

is fine.

Using binaries for those constraints is brute force. You should consider other options.

You seem to have a loss factor in charging and discharging so you should see the model deciding to to both, unless its beneficial to throw away energy, which indicates other modelling issues.

1 Like

OK,thank you so much!

@joaquimg @odow
My Bilevel optimization code did not report an error. When calling the gurobi solver solver, the gap has always been "- and I haven’t found the optimal solution after 9 hours of solving. So I would like to ask, is it possible to set the gap size and iteration time for the gurobi solver in jump? For example, set the iteration time to 10 minutes?

For your convenience, I attach the source code.

using JuMP, BilevelJuMP, Gurobi
using CSV ,DataFrames


model = BilevelModel(Gurobi.Optimizer, mode = BilevelJuMP.SOS1Mode())

#Known parameters
T=24
a_m1=0.0043;
b_m1=0.15;
c_m1=10;
a_m2=0.0065;
b_m2=0.12;
c_m2=10;
PGm1min=10
PGm1max=600
PGm2min=10
PGm2max=150
PGD_EXmin=-70
PGD_EXmax=150
Plpgmin=0;
Plpgmax=50;
PMTmin=0
PMTmax=60
#Initial value
WES_BT0=[100 100 100]
WES_CD0=[100 100 100]
WES_HC0=[100 100 100]
PMT0=[30 30 30]


Cb=DataFrame(CSV.File("E:\\Julia程序\\Cb.csv"))
Cs=DataFrame(CSV.File("E:\\Julia程序\\Cs.csv"))
Pload=DataFrame(CSV.File("E:\\Julia程序\\Pload.csv"))
PPV=DataFrame(CSV.File("E:\\Julia程序\\PPV.csv"))
PWT=DataFrame(CSV.File("E:\\Julia程序\\PWT.csv"))
PDL=DataFrame(CSV.File("E:\\Julia程序\\PDL.csv"))
PCL=DataFrame(CSV.File("E:\\Julia程序\\PCL.csv"))
PHL=DataFrame(CSV.File("E:\\Julia程序\\PHL.csv"))

Cb=convert(Array,Cb)
Cs=convert(Array,Cs)
Pload=convert(Array,Pload)
PPV=convert(Array,PPV)
PWT=convert(Array,PWT)
PDL=convert(Array,PDL)
PCL=convert(Array,PCL)
PHL=convert(Array,PHL)



#Set lower level variables
@variable(Lower(model), 0 <=PMT[1:T,1:3] <= 60)

@variable(Lower(model), 0<=PES_BT_D[1:T,1:3])
@variable(Lower(model), 0<=PES_BT_C[1:T,1:3])

@variable(Lower(model), 0 <=PES_CD_C[1:T,1:3])
@variable(Lower(model), 0 <=PES_CD_D[1:T,1:3])


@variable(Lower(model), 0 <=PES_HC_C[1:T,1:3])
@variable(Lower(model), 0 <=PES_HC_D[1:T,1:3])


@variable(Lower(model), PEC[1:T,1:3])
@variable(Lower(model), 0<=PEC_EX[1:T,1:3])


@variable(Lower(model), PEH[1:T,1:3])
@variable(Lower(model), 0<=PEH_EX[1:T,1:3])


@variable(Lower(model), PAC[1:T,1:3])
@variable(Lower(model), 0<=PAC_EX[1:T,1:3])


@variable(Lower(model), PHX[1:T,1:3])
@variable(Lower(model), 0<=PHX_EX[1:T,1:3])


@variable(Lower(model), 0<=PREC[1:T,1:3])
@variable(Lower(model), PREC_EX[1:T,1:3])


@variable(Lower(model), 0<=PGB[1:T,1:3])
@variable(Lower(model), PGB_EX[1:T,1:3])


@variable(Lower(model), 0<=WES_BT[t=1:T,i=1:3])
@variable(Lower(model), 0<=WES_CD[t=1:T,i=1:3])
@variable(Lower(model), 0<=WES_HC[t=1:T,i=1:3])



    
#Set Upper level variables

@variable(Upper(model), PGD_EX[1:T,1:3])
@variable(Upper(model), PGm1[1:T,1:1])
@variable(Upper(model), PGm2[1:T,1:2])
@variable(Upper(model), Plpg[1:T,1:1])

#Upper objective  and constraints
@objective(Upper(model), Min, a_m1*sum(PGm1.*PGm1)+b_m1*sum(PGm1)+c_m1+a_m2*sum(PGm2.*PGm2)+b_m2*sum(sum(PGm2))+c_m2-sum((Cb+Cs).*PGD_EX))
@constraints(Upper(model), begin
    PGm1+sum(PGm2,dims=2)+Plpg.==sum(PGD_EX,dims=2)+Pload
    10 .<=PGm1 .<= 600
    10 .<=PGm2 .<= 150
    0 .<=Plpg .<=50
    -70 .<=PGD_EX.<=150
end)


#Upper objective  and constraints
@objective(Lower(model), Min, sum(0.025*PPV+0.029*PWT+0.025*2*PMT+0.015*PEC_EX+0.021*PEH_EX+0.035*PAC_EX+0.012*PHX_EX+0.012*PGB_EX+0.03*PREC_EX+0.028*(PES_BT_D+PES_BT_C)+0.018*(PES_CD_C+PES_CD_D)+0.016*(PES_HC_C+PES_HC_D)+3.45/9.7*(2/0.3*PMT+PGB/0.73)))
@constraints(Lower(model), begin
  #Power balance
    PGD_EX+PWT+PPV+2*PMT+PES_BT_D .== PDL+PEC+PEH+PES_BT_C  

    PEC_EX+PAC_EX+PES_CD_D .== PCL+PES_CD_C 

    PEH_EX+PHX_EX+PES_HC_D .== PES_HC_C+PHL

    PREC_EX*(1-0.8)+PGB_EX.== PHX 

    #Power conversion

    PEC*3 .== PEC_EX

    PEH*3 .== PEH_EX

    PAC .== PREC_EX*0.8

    PAC*0.8 .==PAC_EX

    0.9*PHX .==PHX_EX

    PREC*0.73 .==PREC_EX

    PGB*0.73 .==PGB_EX
    
#The initial value of capacity is equal to the final value
   WES_BT0.==WES_BT[24,:] 

   WES_CD0.==WES_CD[24,:]

   WES_HC0.==WES_HC[24,:]
   
    end)


    #Climbing constraints
    for i =1:3
        @constraint(Lower(model), -30 <=PMT[1,i]-PMT0[1,i] )
        @constraint(Lower(model),   PMT[1,i]-PMT0[1,i]<=30)


    end

    for i=1:3
        for t=2:T
            @constraint(Lower(model),  -30 <=PMT[t,i]-PMT[t-1,i] )
            @constraint(Lower(model),  PMT[t,i]-PMT[t-1,i]<=30 )
        end
    end




    

#Energy storage device constraints
    for i in 1:3
            for t in 2:T
                @constraint(Lower(model),WES_BT[1,:] .== WES_BT0+PES_BT_C[1,:]*0.98-PES_BT_D[1,:]/0.98)
                @constraint(Lower(model), WES_BT[t,i] ==  WES_BT[t-1,i]+PES_BT_C[t,i]*0.98-PES_BT_D[t,i]/0.98)

            end
        end


        for i in 1:3
                for t in 2:T
                    @constraint(Lower(model),WES_CD[1,:] .== WES_CD0+PES_CD_C[1,:]*0.98-PES_CD_D[1,:]/0.98)
                    @constraint(Lower(model), WES_CD[t,i] ==  WES_CD[t-1,i]+PES_CD_C[t,i]*0.98-PES_CD_D[t,i]/0.98)

                end
            end

            for i in 1:3
                    for t in 2:T
                        @constraint(Lower(model),WES_HC[1,:] .== WES_HC0+PES_HC_C[1,:]*0.98-PES_HC_D[1,:]/0.98)
                        @constraint(Lower(model), WES_HC[t,i] ==  WES_HC[t-1,i]+PES_HC_C[t,i]*0.98-PES_HC_D[t,i]/0.98)
                    end
            end











optimize!(model)

objective_value(model)


PGm1=value.(PGm1)
PGm2=value.(PGm2)
PGD_EX=value.(PGD_EX)
Plpg=value.(Plpg)



PMT=value.(PMT)
PEC_EX=value.(PEC_EX)
PEH_EX=value.(PEH_EX)

PAC_EX=value.(PAC_EX)
PHX_EX=value.(PHX_EX)

PREC_EX=value.(PREC_EX)
PGB_EX=value.(PGB_EX)
PHX=value.(PHX)
PEH=value.(PEH)
PEC=value.(PEC)
PES_BT_C=value.(PES_BT_C)
PES_BT_D=value.(PES_BT_D)
WES_BT=value.(WES_BT)
PES_CD_C=value.(PES_CD_C);
PES_CD_D=value.(PES_CD_D);

PES_HC_C=value.(PES_HC_C);
PES_HC_D=value.(PES_HC_D);

WES_CD=value.(WES_CD)
WES_HC=value.(WES_HC)

@joaquimg @odow I don’t think the model is complicated. Why did it take so long to solve? Maybe there’s something wrong with my writing?
If the iteration time is set, will it output the solution when the iteration time is reached?

I added the following code

model = BilevelModel(Gurobi.Optimizer, mode = BilevelJuMP.SOS1Mode())
set_optimizer_attribute(model, "TimeLimit", 100)

The result is wrong. The same code can be run in JuMP

ERROR: LoadError: MethodError: no method matching set_optimizer_attribute(::BilevelModel, ::String, ::Int64)
Closest candidates are:
  set_optimizer_attribute(::Model, ::String, ::Any) at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\JuMP\YXK4e\src\JuMP.jl:520
  set_optimizer_attribute(::Model, ::MathOptInterface.AbstractOptimizerAttribute, ::Any) at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\JuMP\YXK4e\src\JuMP.jl:541
Stacktrace:
 [1] top-level scope at E:\Julia程序\加整数0-1测试.jl:6
in expression starting at E:\Julia程序\加整数0-1测试.jl:6

See https://jump.dev/JuMP.jl/v0.21.1/solvers/#JuMP.optimizer_with_attributes

model = BilevelModel(
    optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => 100)
)
1 Like

Thank you very much!The code you gave is feasible.
But my idea of ending the solution process ahead of time seems to be wrong, because gap is still infinite, and there is no solution after the final arrival time.
I think this seems to have something to do with the length of my underlying objective function? I really don’t know how to shorten the solution time. The solver has not worked out for 10 hours, so I feel helpless

In theory Bilevel Optimization is very hard, it is at least as hard MIP.
So it might happen that moderate sized instances take forever to solve.
In the readme: https://github.com/joaquimg/BilevelJuMP.jl

There are some additional reformulations for bilevel programming that might be better in practice. No guarantees though.

I you just need a feasible solution and there is no need for optmality gurantees, you can try product mode and ipopt/knitro.

If you have good bounds for you variables you can try the fourtuny amat mode with MIP solvers, note that if your bounds are not carefully chosen you might reach suboptimal solutions or even have an infeasible problem.

2 Likes

@joaquimg
First of all, thank you for your reply.
I’ve gone https://github.com/joaquimg/BilevelJuMP.jl/blob/master/test/jump.jl
After reading your writing method of fourtuny Amat mode,but I have limited ability and I really can’t understand some places. (for example, how to select the bound value of dual variable and original variable).
May I trouble you for a few minutes, can you use the model of the readme file BilevelJuMP.FortunyAmatMcCarlMode () write it down,

using JuMP, BilevelJuMP, Cbc

model = BilevelModel(Cbc.Optimizer, mode = BilevelJuMP.SOS1Mode())

@variable(Lower(model), x)
@variable(Upper(model), y)

@objective(Upper(model), Min, 3x + y)
@constraints(Upper(model), begin
    x <= 5
    y <= 8
    y >= 0
end)

@objective(Lower(model), Min, -x)
@constraints(Lower(model), begin
     x +  y <= 8
    4x +  y >= 8
    2x +  y <= 13
    2x - 7y <= 0
end)

optimize!(model)

objective_value(model) # = 3 * (3.5 * 8/15) + 8/15 # = 6.13...
value(x) # = 3.5 * 8/15 # = 1.86...
value(y) # = 8/15 # = 0.53...
using JuMP, BilevelJuMP, Gurobi

model = BilevelModel(Gu.Optimizer, mode = BilevelJuMP.ProductMode())

@variable(Lower(model), x)
@variable(Upper(model), y)

@objective(Upper(model), Min, 3x + y)
@constraints(Upper(model), begin
    x <= 5
    y <= 8
    y >= 0
end)

@objective(Lower(model), Min, -x)
@constraints(Lower(model), begin
     x +  y <= 8
    4x +  y >= 8
    2x +  y <= 13
    2x - 7y <= 0
end)

optimize!(model)

objective_value(model) # = 3 * (3.5 * 8/15) + 8/15 # = 6.13...
value(x) # = 3.5 * 8/15 # = 1.86...
value(y) # = 8/15 # = 0.53...

@odow @joaquimg
In response to your suggestion, I also tried to modify the model given in the readme file.However, errors are reported as follows. I don’t quite understand why gurobi reports errors. According to the rules, the gurobi solver can solve MIP and QP problems.

Presolve removed 3 rows and 0 columns
ERROR: LoadError: Gurobi.GurobiError(10020, "Q matrix is not positive semi-definite (PSD). Set NonConvex parameter to 2 to solve model.")
Stacktrace:
 [1] optimize at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\Gurobi\7YNJV\src\grb_solve.jl:7 [inlined]
 [2] optimize!(::Gurobi.Optimizer) at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\Gurobi\7YNJV\src\MOI_wrapper.jl:1920
 [3] optimize!(::MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}) at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\MathOptInterface\bygN7\src\Bridges\bridge_optimizer.jl:239
 [4] optimize!(::BilevelModel; lower_prob::String, upper_prob::String, bilevel_prob::String, solver_prob::String) at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\BilevelJuMP\Ll9sk\src\jump.jl:993
 [5] optimize!(::BilevelModel) at C:\Users\61940\.juliapro\JuliaPro_v1.4.2-1\packages\BilevelJuMP\Ll9sk\src\jump.jl:919
 [6] top-level scope at untitled-bc0232a4617cfc2f29e46dfcf3e660ff:23
in expression starting at untitled-bc0232a4617cfc2f29e46dfcf3e660ff:23