Hi!
I’m working on an MINLP model and solving it with KNITRO.
Since it’s an MINLP, the solver solution heavily depends on solver settings and I’m looking for new ways to validate it.
I was trying to take the same model and solve it with Juniper + Ipopt but it doesn’t seem to work…
Any other solver recommendations? Or how to improve the Juniper performance? I don’t know if SCIP works with my problem formulation.
MINLP is highly dependent on the model formulation. If you post it here, people may have suggestions for improvements.
Well it took me a while but this is a first draft of a MWE:
using JuMP, InfiniteOpt
function solvePolicies(optimizer, # model optimizer
data::Dict, # information for the model (mainly parameters), these are the devices (EV, BESS, PV, etc.), the costs (interests, capex, etc) and the revenues
)
# Sets
Dt = 0:0.25:24; # time discretization
tend=Dt[end]
t0=Dt[1]; # initial time, can´t divide by 0
model = InfiniteModel(optimizer) # create model
# define cont t
@infinite_parameter(model, t ∈ [t0, tend], supports = collect(Dt),
derivative_method = FiniteDifference(Backward()))
# Electrical
# BESS
PbessMin = -12.5; # Min power [kW]
PbessMax = 12.5; # Max power [kW]
SoCbessMin = 0.2; # Min State of Charge [p.u.]
SoCbessMax = 1.0; # Max State of Charge [p.u.]
vmin = 2.8; vmax = 4.2; # voltage limits
Nsbess = 100; # number of cells in series
Npbess = 7; # number of cells in parallel
# EV
SoCevMin = 0.2; # Min State of Charge [p.u.]
SoCevMax = 1.0; # Max State of Charge [p.u.]
aOCV = 3.2; bOCV = 0.7;
SoCev0 = 0.5;
Qev0 = 5.2; # [Ah/cell]
Npev = 25; # number of cells in parallel
Nsev = 100; # number of cells in series
# Thermal
# ST
ηST = 0.6; # conversion factor from Electric PV to thermal
PhpRated = 5; # rated power of the heat pump
ηHP = 3.; # conversion factor from Electric to thermal heat pump
# Extract data
Qtess = 200.; # Capacity [kWh]
PtessMin = -5.; # Min power [kW]
PtessMax = 5.; # Max power [kW]
SoCtessMin = 0.2; # Min State of Charge [p.u.]
SoCtessMax = 0.9; # Max State of Charge [p.u.]
SoCtess0 = 0.5; # Initial State of Charge [p.u.]
ηtess = 0.95; # thermal efficiency
# grids and balances
PgMin = -17.5; # Min power going out
PgMax = 17.5; # Max power coming in
ηg = 0.9; # converter efficiency
# Electrical
# PV
MPPTmeas = data["SPV"].MPPTData[it0:itend];
loadTh = data["grid"].loadTh[it0:itend];
loadElec = data["grid"].loadE[it0:itend];
[γ[n,:] = data["EV"][n].driveInfo.γ[it0:itend] for n ∈ 1:sets.nEV];
γ_interp = linear_interpolation((1:nEV, Dt), γ)
priceBuy=hcat(data["grid"].λ[it0:itend, 1].*1e-3/3600) # convert from €/MWh to €/kWs
priceSell=hcat(data["grid"].λ[it0:itend, 2].*1e-3/3600) # convert from €/MWh to €/kWs
# simulated forecast
@parameter_function(model, PpvMPPT == (t) -> MPPTmeas[zero_order_time_index(Dt, t)])
@parameter_function(model, Pst == (t) -> ηST*MPPTmeas[zero_order_time_index(Dt, t)])
@parameter_function(model, Plt == (t) -> loadTh[zero_order_time_index(Dt, t)])
@parameter_function(model, Ple == (t) -> loadElec[zero_order_time_index(Dt, t)])
@parameter_function(model, γ_cont[n in 1:nEV] == (t) -> γ_interp(n, t))
@parameter_function(model, λbuy == (t) -> priceBuy[zero_order_time_index(Dt, t)])
@parameter_function(model, λsell == (t) -> priceSell[zero_order_time_index(Dt, t)])
OCVev0 = [aOCV[n]+bOCV[n]*SoCev0[n] for n in 1:nEV];
Pdrive = [rand(truncated(Normal(3.5, 1); lower = 0.01)) for n in 1:nEV]; # Gaussian distribution
Ereq=[sum(Pdrive[n].*(1 .-γ[n,:])*Δt)./3600 for n in 1:nEV];
# check if the driving power is greater than the energy in the battery pack.
while any(Ereq .> Qev0.*Npev.*Nsev.* vmax./1000*0.8)
Pdrive = [rand(truncated(Normal(3.5, 1); lower = 0.01)) for n in 1:nEV]; # Gaussian distribution
Ereq=[sum(Pdrive[n].*(1 .-γ[n,:])*Δt)./3600 for n in 1:nEV];
end
@variables(model, begin
# Grid connection
Pg, Infinite(t) # grid power
bPg, Infinite(t), Bin # Binary variable for Grid power
0 ≤ PgPos, Infinite(t) # Pg^+ out/buy power
PgNeg ≤ 0, Infinite(t), (start=0) # Pg^- in/sell power
# BESS
Pbess, Infinite(t) # output power
bPbess, Infinite(t), Bin # Binary variable for output power
SoCbessMin ≤ SoCbess ≤ SoCbessMax, Infinite(t) # State of Charge
0 ≤ PbessPos, Infinite(t) # Pbess^+ out power
PbessNeg ≤ 0, Infinite(t) # Pbess^- in power
vmin ≤ OCVbess ≤ vmax, Infinite(t), (start=OCVbess0) # open circuit voltage of the cell
ibess, Infinite(t), (start=1e3*Pbess0/Npbess/Nsbess/OCVbess0) # current per branch
# EV
Pev[n ∈ 1:nEV], Infinite(t) # EV charger power
bPev[n ∈ 1:nEV], Infinite(t), Bin # EV charger power
PevTot[n ∈ 1:nEV], Infinite(t) # total power of each EV, driving+V2G
SoCevMin[n] .≤ SoCev[n ∈ 1:nEV] .≤ SoCevMax[n], Infinite(t) # State of Charge
0 ≤ PevPos[n ∈ 1:nEV], Infinite(t) # Pev^+ out power
PevNeg[n ∈ 1:nEV] ≤ 0, Infinite(t) # Pev^- in power
vmin[n] ≤ OCVev[n ∈ 1:nEV] ≤ vmax[n], Infinite(t),(start=OCVev0[n]) # open circuit voltage of the cell
iev[1:nEV], Infinite(t) # current per branch
# HP
0 ≤ Phpe ≤ PhpRated, Infinite(t) # Heat pump power
# TESS
SoCtess ≥ SoCtessMin, Infinite(t) # State of Charge
auxTess ≥ 0, Infinite(t) # slack variable for SoCtess
Ptess, Infinite(t) # Thermal power
bPtess, Infinite(t), Bin # Binary variable for TESS power
0 ≤ PtessPos, Infinite(t) # Ptess^+ out power
PtessNeg ≤ 0, Infinite(t) # Ptess^- in power
end)
t1 = t0 + 6*3600; OCVbess0 = aOCV+bOCV*SoCbess0;
@constraints(model, begin
# Grid connection
PgPos ≤ (1-bPg)*PgMax
bPg*PgMin ≤ PgNeg
PgNeg * (1/ηg) + PgPos * ηg .== Pg
# BESS
PbessNeg * (1/ηbess) + PbessPos * ηbess .== Pbess
bPbess*PbessMin ≤ PbessNeg
PbessPos ≤ (1-bPbess)*PbessMax
OCVbess == aOCV+bOCV*SoCbess # Linear voltage model
ibess == 1e3*Pbess/Npbess/Nsbess/OCVbess # current per branch. 1e3 to convert kW->W
SoCbess(t0) == SoCbess0 # Initial conditions
SoCbess(t1) == SoCbess(t1+24*3600) # periodic condition
∂.(SoCbess, t) .== -(ηbess * bPbess + (1-bPbess))*ibess/Qbess0/3600 # Static Qbess
# EV
[n ∈ 1:nEV], PevNeg[n] * (1/ηev[n]) + PevPos[n] * ηev[n] == Pev[n]
[n ∈ 1:nEV], bPev[n]*PevMin[n] ≤ PevNeg[n]
[n ∈ 1:nEV], PevPos[n] ≤ (1-bPev[n])*PevMax[n]
[n ∈ 1:nEV], SoCev[n](t0) == SoCev0[n] # Initial conditions
availability[n ∈ 1:nEV], γ_cont[n].*Pev[n] + (1-γ_cont[n]).*Pdrive[n] - PevTot[n] .== 0 # power balance
[n ∈ 1:nEV], OCVev[n] .== aOCV[n]+bOCV[n]*model[:SoCev][n] # linear voltage model
[n ∈ 1:nEV], iev[n] .== 1e3*model[:PevTot][n]/Npev[n]/Nsev[n]/OCVev[n] # current per branch
[n ∈ 1:nEV], ∂.(SoCev[n], t) .== -(ηev[n] * bPev[n] + (1-bPev[n]))*iev[n]/Qev0[n]/3600
# TESS
bPtess*PtessMin ≤ PtessNeg
PtessNeg * (1/ηtess) + PtessPos * ηtess .== Ptess
PtessPos ≤ (1-bPtess)*PtessMax
auxTess ≥ SoCtess - SoCtessMax # slack variable soft constraint
SoCtess(t0) .== SoCtess0 # Initial condition
∂.(SoCtess, t) .== -Ptess/Qtess/3600 # Bucket model
end)
# CHECK ALL OF THIS tDep are the departure times from γ
model[:ϵSoC]=[]; # assign name in the model
for day in eachindex(tDep[1,:]) # if there's more than one day loop over them
for n in 1:nEV # loop over the EVs
td = tDep[n, day] # pick value
td = td * 3600 + (24 * 3600 * (day - 1)) # change to secs and add days
# find the nearest td inside supports(t)
td = findmin(abs.(td .- supports(t)))[1]
# Check if the time index is within the bounds of SoCev
if td >= t0 && td <= tend
ϵSoCexpr = SoCev[n](td) - SoCdep[n]
push!(model[:ϵSoC], ϵSoCexpr) # push to the expression vector
end
end
end
# Thermal balance
model[:thBalance]=@constraint(model, Pst+Phpe.*ηHP+Ptess .== Plt);
# Power balance DC busbar
Pev_sum = nEV != 1 ? sum(γ_cont[n].*Pev[n] for n in 1:nEV) : sum(γ_cont.*Pev[n] for n in 1:nEV)
model[:powerBalance]=@constraint(model, PpvMPPT + Pbess + Pev_sum + Pg .== Ple + Phpe)
set_start_value_function(Ptess, t -> loadTh[zero_order_time_index(Dt, t)])
set_start_value_function(Pg, t -> loadElec[zero_order_time_index(Dt, t)] .- MPPTmeas[zero_order_time_index(Dt, t)])
# Add objective function
costFunction!(model, sets, data)
# Grid costs
Wgrid = 1; # regularization factor for grid cost. max(λ)*max(P)
Wlims = 1000; # penalty for TESS overcharging
WSoCDep = 1000; # penalty for EVs not reaching the desired SoC
cgrid = Wgrid .* (PgPos*λbuy + PgNeg*λsell);
pDep = WSoCDep*sum(ϵSoC[n]^2 for n ∈ eachindex(ϵSoC))
@objective(model, Min, ∫(cgrid,t)/tend + pDep + Wlims*∫(auxTess,t)/tend)
return model
end;
it should be MIQCP/MINLP.
So what’s the question, specifically? We can’t run your problem because you haven’t provided the input data.
Also, since this is really an InfiniteOpt model, I don’t know if the usual JuMP suggestions directly apply. @pulsipher may have suggestions.
I am happy to help. Could you explain more of the problem you are running into? Does it throw an error? MINLPs like this are very difficult to solve, so some solvers might not be able to converge in any reasonable amount of time.
One potential issue with models that use a lot of parameter functions is the final constraint structure might not be what you expect. For instance, suppose p is a parameter function and y is a variable which are used in the following constraint:
@constraint(model, p*y^2 <= 42)
Even though we know that this a a quadratic constraint, it will be stored as a general nonlinear constraint since it has to store the references for p and y^2. So, if you call a solver like Gurobi that supports MIQCPs, it will throw an error since it will think you have a general nonlinear constraint.
JuMP has a similar limitation with the use of Parameter sets.
Developing an intelligent way of automatically detecting and converting such expressions down to quadratic or linear expressions during transcription is on the TODO list, but I haven’t had time to look into it much yet.