f(t+1) in JuMP

Specific Domains Optimization (Mathematical)
,
1

Hi everyone!

I am just starting with JuMP and I’m trying to solve a little battery optimization problem that I used to solve easy in GAMS. Though, I am really stuck in the SoC code because I don’t know how to write f(t+1) in Julia and I cannot find the answer in the forums.

Can someone help me with that?

thank you!.

2

You have not posted enough information for us to answer your question. For example, the way to write f(t+1) in Julia is usually as a function call, i.e. f(t+1). However, in a JuMP model with discrete variables you may want an array reference like variablename[t+1]. Or it could be other things depending on your code.

So if you want help, please include a small but relevant section of your code. And please format your code using backticks, like this.

3

Thank you Niclas,

I’m going to try to clarify everything.

This is the ecuation I want to implement:

image

and this is the implementation in GAMS:

image

The idea of this ecuation is that the energy in the following period is going to be de energy in the current period plus Pc ( energy consume in the current period ) or minus Ph ( energy generated in the current period ).

and finally here is my code in Julia, sorry for all the mistakes but is my first full code and I still don’t understand some aspects of this language


#reading demand data.

csvfilename = "demand.csv"
csvdata = readcsv(csvfilename,  header=true)
data = csvdata[1]
header = csvdata[2]

Pc= data[:,1]

#Prmax                                                                 max grid potency
#Pr                                                                       grid potency
#Pb                                                                      Battery storage potency ( + - ).
#Pc                                                                      Direct consume potency from grid.
#Pbmin y Pbmax                                                 Minimum and maximum potency battery
#Cpot                                                                  Cost kW
#Cen                                                                   Cost kWh
#Enom                                                                Battery nominal energy

using JuMP, Cbc
m = Model(solver=CbcSolver())
Prmax = 3
Pbmax = 1.7
Pbmin = -1.7
Cp = 0.15
Enom = 1.7
Cen = [0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.012,0.012,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.12,0.12,0.12,0.12,0.12]

# Declaring variables

@variable(m, 0<= Pr[t=1:24] <= Prmax)
@variable(m, Pbmin <= Pb[t=1:24] <= Pbmax )
@variable(m, 0.1*Enom <= E[t=1:24] <= 0.9*Enom)

# Setting the objective

@objective(m, Min, (Prmax*Cpot) + ( Pr[t=1:24]*Cen[t=1:24] ) 

# Adding constraints
@constraint(m, constraint1, Pr[t=1:24] .== Pc[t=1:24] - Pb[t=1:24] )

# This is the problematic constraint

@constraint(m, constraint2, E[(t=1:24)+1] .== E[t=1:24]- Pr[t=1:24] )




@constraint(m, constraint3, Pbmin .<= Pb[t=1:24] .<= Pbmax)
@constraint(m, constraint4, (Enom*0.1) .<= E[t=1:24] .<= (Enom*0.9) )

# Printing the prepared optimization model

print(m)

# Solving the optimization problem

solve(m)

# Printing the optimal solutions obtained

println("Optimal Solutions:")
println("Pb = ", getvalue(Pb))
println("Pr = ", getvalue(Pr))
println("E = ", getvalue(E))

I hope this could help you to understand the problem.

thank you very much!!!!
4

Good, now it’s quite clear. You’re actually quite close to a working model, just a few problems:

  • Cpot is not defined. I set it to 0 for now.
  • Your objective function is a vector. You probably want to sum its elements to make a scalar. You also need elementwise multiplication .* here.
  • And now your question, how to do E[t+1] in constraint2. You currently make a range E[(1:24)+1] which is the same as writing E[2:25]. That errors since your vectors only have 24 elements. What you probably want is 2:24 on the left and 1:23 on the right. Depending on your problem you may need special explicit conditions for the first element, possibly the last as well.

So the minimum changes you need to do to get a model that runs without errors (note that I haven’t checked the logic of your equations) is this:

@objective(m, Min, Prmax*0 + sum(Pr[t=1:24].*Cen[t=1:24])) 

@constraint(m, constraint2, E[t=2:24] .== E[t=1:23] - Pr[t=1:23])

Also, you aren’t using any of the summation indexes t, so you can drop all the occurances of t= in your model. (In many cases you can drop the ranges 1:24 too, e.g. just use sum(Pr.*Cen) in the objective.) Finally, there are many ways to rewrite your model more cleanly, but I’m short on time at the moment so I’ll stop here. Maybe someone else will jump in or I can do it some other day if you want tips.

5

Hi Niclas!

I’m going to try all your suggestions! and when it works I may try to rewrite it cleanly.

Thank you very much for your time!!