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 Solvers · JuMP
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