Modelling availability of assets

[DarioSlaifsteinSk]

Continuing with @parameter_function
I´m adding Electric Vehicles to the problem but I´m having trouble with one particular constraint.
The house has EV chargers with power $P_{EV,n}(t)$ with $n=1, ..., N_{EV}$ the set of EVs.
While the car is not connected its driving and consuming $P_{drive}$
Then the total power in the EV is the power charged and the power for driving following

$P_{tot, n}=\gamma_n P_{EV, n} - (1-\gamma_n) P_{drive} $

$\gamma=1$ when the EV is plugged and $0$ when it´s not.
The implementation I did is first creating a matrix containing the discrete values along the simulation and then passing that to the continuous time domain

nEV=1:2;  ndays=1
tend=24*ndays
Dtw=0:Δt:tend

# EV availability
    # First, we create availability vectors for each EV in the disc t-domain.
    tDep=zeros(length(nEV), ndays);
    tArr=zeros(length(nEV), ndays);
    γ=ones(length(nEV),1); # start with dummy col
    for day in 1:ndays # loop through the days
        tDep[:,day] = rand(6:12, nEV[end]); # departure time
        tArr[:,day] = tDep[:,day]+rand(1:8, nEV[end]); # arrival time
        av = ones(length(nEV), 24);
        for n in nEV
            td=Integer(tDep[n,day]); ta=Integer(tArr[n,day]);
            av[n, td:ta] .= 0;
        end
        γ=[γ av]
        # correction for t-axis
        tDep[:,day] .= tDep[:,day].+24*(day-1);
        tArr[:,day] .= tArr[:,day].+24*(day-1);
    end
    γ=repeat(γ, inner=(1,Integer(1/Δt))) # upsample
    γ=γ[:,Integer(1/Δt):end] # eliminate dummy cols
    # Now we need to project it into the cont t-domain.
    γ_interp = linear_interpolation((nEV,Dtw), γ)
    @parameter_function(model, γ_cont[nEV] == (t) -> γ_interp(nEV, t)) # make InfiniteOpt compatible

After that I sample the driving power from a distribution and create the variables

    # Driving consumption
    μDrive=5; σDrive=2
    Pdrive = rand(Normal(μDrive,σDrive), length(nEV)); # Gaussian distribution
    # Add variables
    @variable(model, Pev[nEV], Infinite(t))  # EV charger power
    @variable(model, PevTot[nEV], Infinite(t))  # total power of each EV, driving+V2G
    # Dummy variables for bidirectional flow
    @variable(model, 0 ≤ PevPos[n ∈ nEV] ≤ PevMax[n], Infinite(t)) # Pev^+ out power
    @variable(model, PevMin[n] ≤ PevNeg[n ∈ nEV] ≤ 0, Infinite(t)) # Pev^- in power
    @constraint(model, [n ∈ nEV], PevNeg[n] + PevPos[n] == Pev[n])
    @constraint(model, [n ∈ nEV], PevNeg[n] * PevPos[n] == 0)

The constraint that is giving me problems is the following:

@constraint(model, [n ∈ nEV], PevTot[n] == Pev[n]*γ_cont[n] - (1-γ_cont[n])*Pdrive[n]) # power balance

I think it is correctly generated before transcribing the model. But when I do:

build_optimizer_model!(model)

I get

MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number at twiceprecision.jl:273
  convert(::Type{T}, ::AbstractChar) where T<:Number at char.jl:185
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number at multidimensional.jl:130
Stacktrace:
  [1] call_function(fref::ParameterFunctionRef, supps::Float64)
    @ InfiniteOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\expressions.jl:310
  [2] call_function(fref::GeneralVariableRef, support::Float64)
    @ InfiniteOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\general_variables.jl:696
  [3] lookup_by_support(model::Model, fref::GeneralVariableRef, index_type::Type{ParameterFunctionIndex}, support::Vector{Float64})
    @ InfiniteOpt.TranscriptionOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\TranscriptionOpt\model.jl:483
  [4] transcription_expression(trans_model::Model, vref::GeneralVariableRef, index_type::Type{ParameterFunctionIndex}, support::Vector{Float64})
    @ InfiniteOpt.TranscriptionOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\TranscriptionOpt\transcribe.jl:431
  [5] transcription_expression(trans_model::Model, vref::GeneralVariableRef, support::Vector{Float64})
    @ InfiniteOpt.TranscriptionOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\TranscriptionOpt\transcribe.jl:419
  [6] #88
    @ C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\TranscriptionOpt\transcribe.jl:509 [inlined]
  [7] _ast_process_node(map_func::InfiniteOpt.TranscriptionOpt.var"#88#89"{Model, Vector{Float64}}, quad::GenericQuadExpr{Float64, GeneralVariableRef})
    @ InfiniteOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\nlp.jl:316
  [8] _tree_map_to_ast(map_func::Function, node::LeftChildRightSiblingTrees.Node{NodeData})
    @ InfiniteOpt C:\Users\daros\.julia\packages\InfiniteOpt\DB7Ch\src\nlp.jl:326
  [9] map_nlp_to_ast
...
    @ In[25]:1
 [20] eval
    @ .\boot.jl:368 [inlined]
 [21] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1428

Here are the constraints generated in latex before transcribing:

Can I change something on how its transcribed?
This should be a linear expression since $\gamma$ is a function, not a variable. If not, how else should I model this behaviour?

[pulsipher]

Here the problem is that each particular parameter function needs to return a single value, not a vector of values. Here it appears that γ_cont[nEV] at each index of nEV (which indexes individual parameter functions) returns a 2 element vector. This can be fixed by account for this when you declare the parameter function:

 @parameter_function(model, γ_cont[i in 1:2, nEV] == (t) -> γ_interp(nEV, t)[i])

You can place the additional index dimension either before or after nEV and then you can use γ_cont appropriately in your formulations.

P.S. For future posts, please provide a minimal working example if at all possible to make easier to help you. Please see https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757 for more info.

[DarioSlaifsteinSk]

Yes! Sorry I forgot a few lines. Here's the min. working example:

using JuMP, InfiniteOpt, Ipopt, Juniper, Plots, LinearAlgebra, Distributions, Statistics
using LaTeXStrings

nEV=1:2;  ndays=1; Δt=0.25;
tend=24*ndays
Dtw=0:Δt:tend

PevMax = [20, 20]; PevMin = [-20, -20]

model = InfiniteModel(Ipopt.Optimizer) # create model
@infinite_parameter(model, t in [0, tend], supports = Dtw, derivative_method = FiniteDifference(Backward()) )

# EV availability
# First, we create availability vectors for each EV in the disc t-domain.
tDep=zeros(length(nEV), ndays);
tArr=zeros(length(nEV), ndays);
γ=ones(length(nEV),1); # start with dummy col
for day in 1:ndays # loop through the days
    tDep[:,day] = rand(6:12, nEV[end]); # departure time
    tArr[:,day] = tDep[:,day]+rand(1:8, nEV[end]); # arrival time
    av = ones(length(nEV), 24);
    for n in nEV
        td=Integer(tDep[n,day]); ta=Integer(tArr[n,day]);
        av[n, td:ta] .= 0;
    end
    γ=[γ av]
    # correction for t-axis
    tDep[:,day] .= tDep[:,day].+24*(day-1);
    tArr[:,day] .= tArr[:,day].+24*(day-1);
end
γ=repeat(γ, inner=(1,Integer(1/Δt))) # upsample
γ=γ[:,Integer(1/Δt):end] # eliminate dummy cols
# Now we need to project it into the cont t-domain.
γ_interp = linear_interpolation((nEV,Dtw), γ);
@parameter_function(model, γ_cont[nEV] == (t) -> γ_interp(nEV, t)) # original
@parameter_function(model, γ_cont1[n in nEV] == (t) -> γ_interp(nEV, t)[n]) # I think its the right one
@parameter_function(model, γ_cont2[i in 1:2, nEV] == (t) -> γ_interp(nEV, t)[i]) # original from Pulsipher

Pdrive = rand(Normal(5,2), length(nEV)); # Gaussian distribution

@variable(model, Pev[nEV], Infinite(t))  # EV charger power
@variable(model, PevTot[nEV], Infinite(t))  # total power of each EV, driving+V2G
# Dummy variables for bidirectional flow
@variable(model, 0 ≤ PevPos[n ∈ nEV] ≤ PevMax[n], Infinite(t)) # Pev^+ out power
@variable(model, PevMin[n] ≤ PevNeg[n ∈ nEV] ≤ 0, Infinite(t)) # Pev^- in power
@constraint(model, [n ∈ nEV], PevNeg[n] + PevPos[n] == Pev[n])
@constraint(model, [n ∈ nEV], PevNeg[n] * PevPos[n] == 0)
@constraint(model, [n ∈ nEV], PevTot[n] == Pev[n]*γ_cont1[n] - (1-γ_cont1[n])*Pdrive[n]) # power balance
build_optimizer_model!(model)

A few things after playing a bit with the 3 definitions given above

1. First of all with your definition using both i and nEV I get the following:

print(γ_cont2)

2-dimensional DenseAxisArray{GeneralVariableRef,2,...} with index sets:
    Dimension 1, Base.OneTo(2)
    Dimension 2, 1:2
And data, a 2×2 Matrix{GeneralVariableRef}:
 γ_cont2[1,1](t)  γ_cont2[1,2](t)
 γ_cont2[2,1](t)  γ_cont2[2,2](t)

which is not what I need I think. What I need that every car in nEV has its own $\gamma(t)$ function varying over time.
This gives me 4 functions instead of 2. Additionally I don't know what's being broadcasted into the different elements [1,1], [1,2] and so on.

2. How do I check the dimension of the output of γ_cont?
   Because for what I see, when I print it:

print(γ_cont)

1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, [1, 2]
And data, a 2-element Vector{GeneralVariableRef}:
 γ_cont[1](t)
 γ_cont[2](t)

doesn't that mean that each γ_cont[n](t) for n in nEV I get scalar for every value of t?
furthermore when handling the γ_interp I plot them as:

plot(Dtw,γ_interp(1,Dtw), label=L"EV_{1}\ availability")
plot!(Dtw,γ_interp(2,Dtw), label=L"EV_{2}\ availability")

which means that for every value of Dtw I have an scalar output of γ_interp(n,t).
how is this syntax

@parameter_function(model, γ_cont[nEV] == (t) -> γ_interp(nEV, t)) # original

giving me a vector output for γ_cont

3. Finally I think that correct syntax might be with γ_cont1 since the build_optmizer_model!(model) doesn't crash.
   But yet again, what is actually being broadcasted into the n elements?

Thanks again for the help.
I'm gonna keep playing with it and see the full result if it works.

[pulsipher]

Sorry for the confusion in my last response, it greatly helps to see the full model. I know understand what is going on and can explain.

First, you should rewrite the time definition to:

@infinite_parameter(model, t in [0, tend], supports = collect(Dtw), derivative_method = FiniteDifference(Backward()) )

since a vector is expected for supports. I probably should add direct support for ranges to avoid this....

Now let's go through the varies parameter function definitions above. First we have:

@parameter_function(model, γ_cont[nEV] == (t) -> γ_interp(nEV, t))

This is problematic since nEV = 1:2 so this returns a vector for all the nEV values at time t. Since, InfiniteOpt expects scalar output, this throws an error:

julia> call_function(γ_cont[1], 0)
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number at Julia-1.7.1\share\julia\base\twiceprecision.jl:262
  convert(::Type{T}, ::AbstractChar) where T<:Number at Julia-1.7.1\share\julia\base\char.jl:185
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number at Julia-1.7.1\share\julia\base\multidimensional.jl:136
  ...
Stacktrace:
 [1] call_function(fref::ParameterFunctionRef, supps::Int64)
   @ InfiniteOpt ~InfiniteOpt.jl\src\expressions.jl:310
 [2] call_function(fref::GeneralVariableRef, support::Int64)
   @ InfiniteOpt ~InfiniteOpt.jl\src\general_variables.jl:696
 [3] top-level scope
   @ REPL[23]:1

Now for the 2nd definition we have:

@parameter_function(model, γ_cont1[n in nEV] == (t) -> γ_interp(nEV, t)[n])

This will work since it returns a scalar output:

julia> call_function(γ_cont1[1], 0)
1.0

However, this has some unnecessary complexity. γ_interp(nEV, t) returns a vector for all n that we then index at the desired n. We could instead index at the appropriate n in the first place:

julia> @parameter_function(model, γ_cont_direct[n in nEV] == (t) -> γ_interp(n, t))
1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{GeneralVariableRef}:
 γ_cont_direct[1](t)
 γ_cont_direct[2](t)

This creates one parameter function for each n as expected and the function calls are what we expect:

julia> call_function(γ_cont_direct[1], 0)
1.0

julia> γ_interp(1, 0)
1.0

The reason why the wrong syntax printed as you expected, is that the parameter functions where created as you expected, but the underlying function didn't return a scalar as they expect which triggers an error when the functions are called during the discretization model build. Ideally it would be nice to check the functions upfront when defining the parameter functions, but this is nontrivial to accomplish in a general way.

[DarioSlaifsteinSk]

Oh! I get it now!
The syntax:

@parameter_function(model, γ_cont[nEV] == (t) -> γ_interp(nEV, t))

Is broadcasting nEV elements into each nEV element of γ_cont!
Yes, the correct one is the last one! I'll get back to you with the results!