Infeasibilities with economic dispatch model

Specific Domains Optimization (Mathematical)
, , ,
1

Hello,

I am currently working on a network economic dispatch model. The basic transport model below works with only conventional generators. However, whenever I add renewable injection constraints with hourly capacity factors, the model cannot solve and I get an infeasibility message. I have tried various options to no avail. Any assistance would be highly appreciated.

#using DataFrame, JuMP

generators = DataFrame(index = ["conv_gas", "conv_nuc", "conv_oil", "ror_plant"], 
                      node = ["zone_1", "zone_1", "zone_2", "zone_2"],
                      R_ID = 1:4, node_id = [1,1,2,2], 
                                g_max = [100,250,150,60], cost = [50,10,20,0], conventional = [1,1,1,0], 
                                ROR = [0,0,0,1])

demand = DataFrame(load_z1 = [100,150,80], load_z2 = [150,165,170])

variability = DataFrame(conv_gas = [1,1,1], conv_nuc = [1,1,1], conv_oil = [1,1,1], ror_plant = [0.6,0.6,0.6])

network_zones = ["z1","z2"]

lines = DataFrame(network_lines = 1, z1 = 1, z2 = -1)

# Sets
G = generators.R_ID #all generators
T = convert(Array{Int64}, 1:3)
Z = convert(Array{Int64}, 1:2) #zones
L = lines.network_lines 

#SUBSETS

#C = generators.R_ID[generators.conventional.==1] #conventional generators
#ROR = intersect(generators.R_ID[generators.ROR.==1], G) #Run-Of-River

#Dispatch_Model =  Model(Clp.Optimizer)

#@variable(Dispatch_Model, vGEN[T, c in C] >=0) 
#@variable(Dispatch_Model, vROR[T, r in ROR] >=0)
#@variable(Dispatch_Model, vFLOW[T,L])

#Constraints

#@constraint(Dispatch_Model, cMaxPower[t in T, c in C], 
#        vGEN[t,c] <= generators[generators.R_ID .== c,:g_max][1])

@constraint(Dispatch_Model, cRORInjection[t in T, r in ROR], vROR[t,r] <= generators[r,:g_max] * variability[t,r])

#@expression(Dispatch_Model, eCosts, sum(vGEN[t,c]*generators.cost[c] for t in T, c in C))
#@objective(Dispatch_Model, Min, eCosts)

@constraint(Dispatch_Model, cEnergyBalance[t in T, z in Z],
            sum(vGEN[t,c] for c in intersect(generators[generators.node.==z, :R_ID],G))+
            sum(vRUNOFRIVER[t,r] for r in intersect(generators[generators.node.==z, :R_ID],G))
            -demand[t,z] -
            sum(lines[l,Symbol(string("z",z))] * vFLOW[t,l] for l in L) == 0
            )

@time optimize!(Dispatch_Model)

Model output message:

Coin0507I Presolve determined that the problem was infeasible with tolerance of 1e-08
Clp3003W Analysis indicates model infeasible or unbounded
Clp0006I 0  Obj 0 Primal inf 815 (5)
Clp0001I Primal infeasible - objective value 0
Clp0032I PrimalInfeasible objective 0 - 1 iterations time 0.012
  0.011694 seconds (2 allocations: 32 bytes)


Status
  Termination status : INFEASIBLE
  Primal status      : UNKNOWN_RESULT_STATUS
  Dual status        : INFEASIBILITY_CERTIFICATE
  Message from the solver:
  "1 - primal infeasible"

* Candidate solution
  Objective value      : 0.0

* Work counters
  Solve time (sec)   : 0.01200

Thanks

2

Do you expect there to always be a feasible solution? I can’t run your problem, because I don’t know which lines are meant to be commented or not.

Some tips to work on this:

  • Simplify the problem. Remove constraints until it gets feasible
  • Once you know which constraint causes the infeasibility, check the data. Is there a problem with the units? Did you get the sign of each term correct? It’s easy to have a - instead of a +.
  • Make the problem smaller. Simplify the data, use less generators and less time periods
  • Try changing the constraint. If it’s an equality, relax into >=, or add a slack variable that you penalize in the objective. What is the optimal value of the slack?
3

Thanks. I re-posted the uncommented codes below. Only one of the constraints is causing the infeasibility (cRORInjection). The model is able to generate a feasible solution with the generator capacity constraint (cMaxPower) and the energy balance (cEnergyBalance) activated.

#using DataFrame, JuMP

generators = DataFrame(index = ["conv_gas", "conv_nuc", "conv_oil", "ror_plant"], 
                      node = ["zone_1", "zone_1", "zone_2", "zone_2"],
                      R_ID = 1:4, node_id = [1,1,2,2], 
                                g_max = [100,250,150,60], cost = [50,10,20,0], conventional = [1,1,1,0], 
                                ROR = [0,0,0,1])

demand = DataFrame(load_z1 = [100,150,80], load_z2 = [150,165,170])

variability = DataFrame(conv_gas = [1,1,1], conv_nuc = [1,1,1], conv_oil = [1,1,1], ror_plant = [0.6,0.6,0.6])

network_zones = ["z1","z2"]

lines = DataFrame(network_lines = 1, z1 = 1, z2 = -1)

# Sets
G = generators.R_ID #all generators
T = convert(Array{Int64}, 1:3)
Z = convert(Array{Int64}, 1:2) #zones
L = lines.network_lines 

#SUBSETS

C = generators.R_ID[generators.conventional.==1] #conventional generators
ROR = intersect(generators.R_ID[generators.ROR.==1], G) #Run-Of-River

Dispatch_Model =  Model()

@variable(Dispatch_Model, vGEN[T, c in C] >=0) 
@variable(Dispatch_Model, vROR[T, r in ROR] >=0)
@variable(Dispatch_Model, vFLOW[T,L])

#Constraints

@constraint(Dispatch_Model, cMaxPower[t in T, c in C], 
       vGEN[t,c] <= generators[generators.R_ID .== c,:g_max][1])

@constraint(Dispatch_Model, cRORInjection[t in T, r in ROR], vROR[t,r] <= generators[r,:g_max] * variability[t,r])

@expression(Dispatch_Model, eCosts, sum(vGEN[t,c]*generators.cost[c] for t in T, c in C))
@objective(Dispatch_Model, Min, eCosts)

@constraint(Dispatch_Model, cEnergyBalance[t in T, z in Z],
            sum(vGEN[t,c] for c in intersect(generators[generators.node.==z, :R_ID],G))+
            sum(vRUNOFRIVER[t,r] for r in intersect(generators[generators.node.==z, :R_ID],G))
            -demand[t,z] -
            sum(lines[l,Symbol(string("z",z))] * vFLOW[t,l] for l in L) == 0
            )

@time optimize!(Dispatch_Model)


#Infeasibility results
Coin0507I Presolve determined that the problem was infeasible with tolerance of 1e-08
Clp3003W Analysis indicates model infeasible or unbounded
Clp0006I 0  Obj 0 Primal inf 815 (6)
Clp0001I Primal infeasible - objective value 0
Clp0032I PrimalInfeasible objective 0 - 1 iterations time 0.012
  7.858292 seconds (18.65 M allocations: 1.102 GiB, 4.93% gc time, 1.36% compilation time)

With regards to the slack variables, I have a couple of quick questions.

  1. When adding the slack variable into the constraint, do i convert the constraint to an equality, like this?
@variable(Dispatch_Model, a)

@constraint(Dispatch_Model, cRORInjection[t in T, r in ROR], vROR[t,r] + a == generators[r,:g_max] * variability[t,r])

  1. Should i use a unique slack variable for each constraint?

  2. Could you direct me to a MWE for penalizing the slack variables in my objective function. I’ve tried doing it, but getting error messages.

Thanks!

4

If I print your energy balance constraint:

julia> cEnergyBalance
2-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape},2,...} with index sets:
    Dimension 1, [1, 2, 3]
    Dimension 2, [1, 2]
And data, a 3×2 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 cEnergyBalance[1,1] : -vFLOW[1,1] = 100.0  cEnergyBalance[1,2] : vFLOW[1,1] = 150.0
 cEnergyBalance[2,1] : -vFLOW[2,1] = 150.0  cEnergyBalance[2,2] : vFLOW[2,1] = 165.0
 cEnergyBalance[3,1] : -vFLOW[3,1] = 80.0   cEnergyBalance[3,2] : vFLOW[3,1] = 170.0

cEnergyBalance[1,1] says that vFLOW[1,1] has to be -100, but cEnergyBalance[1,2] says that vFLOW[1,1] has to be 150.0. They can’t both be true, so that’s the source of your infeasibility.

I suggest you take another look at your formulation.

1 Like
5

Thanks @odow for the insight. The main problem was correctly mapping plants to nodes and then slightly reformulating the energy balance. Model now generates a feasible solution without the need for slack variables.

1 Like