Dummy/phantom constraints

[DarioSlaifsteinSk]

Hi everyone,
I'm building an Energy Management System to control an electrified household's power flows (electricity and heat). Right now the system looks like this:

The system is quite simple given that the following measurements (in the physical sense) are given:
$P_{PV}, P_{load}^{e}, P_{load}^{th}$.
Then the Solar thermal heat is:
$P_{ST}=\eta_{ST}.P_{PV}$
The heat pump also has a very simple model:
$P_{hp}^{t}=\eta_{hp}.P_{hp}^{e}$
So in the end the problem can be summed as:

$$\min_{P^* } C_{grid}$$

$$C_{grid}=\int{\lambda_{sell} .P_{grid}^{+}-\lambda_{buy} .P^{-}_{grid}.dt}$$

subject to:
$$P_{hp}^{t}=\eta_{hp}.P_{hp}^{e}$$
$$P_{ST}=\eta_{ST}.P_{PV}$$
TESS model

$$\underline{SoC}{tess} \leq SoC{tess} \leq \overline{SoC}_{tess}$$

$$ \underline{P}{tess} \leq {P}{tess} \leq \overline{P}_{tess}$$

$$ \underline{P}{tess} ≤ P{tess}^{-} ≤ 0$$

$$ 0 ≤ P_{tess}^{+} ≤ \overline{P}_{tess}$$

$$ P_{tess}^{-} + P_{tess}^{+} = P_{tess}$$

$$ P_{tess}^{-} . P_{tess}^{+} = 0$$

$$ SoC_{tess}(0) = SoC_{tess,0}$$

$$ \frac{dSoC_{tess}}{dt} = \frac{1}{\eta_{tess}}.\frac{P_{tess}}{Q_{tess}}$$

Thermal balance

$$ P_{ST}+P_{HP}^{t}+P_{TESS}=P_{load}^{t}$$

Power balance

$$ P_{PV}+P_{g}=P_{load}^{e}+P_{HP}^{e}$$

This is implemented in the following way:

using JuMP, InfiniteOpt, HiGHS, Ipopt, Juniper, Plots, LinearAlgebra, Distributions, AmplNLWriter, GAMS
using CSV, XLSX, DataFrames, LaTeXStrings

function heatpump(model::InfiniteModel, data::Dict) # heat pump
    # The heat pump has a variable (electrical) and a subordinate finite_param (thermal)
    # Extract params
    t = model[:t];
    PhpRated = data["HP"].RatedPower; # rated power of the heat pump
    # ηHP = data["HP"].η; # electrothermal conv efficiency $/kW
    Phpe0=data["HP"].P0;
    
    @variable(model, - PhpRated ≤ Phpe ≤ PhpRated, Infinite(t), start = Phpe0);  # Electric power
    # @variable(model, - PhpRated*ηHP ≤ Phpt ≤ PhpRated*ηHP, Infinite(t));  # Thermal power
    
    # Dummy variables for bidirectional flow
    @variable(model, 1e-4 ≤ PhpPos ≤ PhpRated, Infinite(t), start = Phpe0);
    @variable(model, -PhpRated ≤ PhpNeg ≤ -1e-4, Infinite(t), start = 0);
    # Bidirectional power flow, ensuring only export or import
    @constraint(model, PhpPos + PhpNeg == Phpe);
    @constraint(model, PhpPos * PhpNeg == 0);
    
    # Thermal power
    # @constraint(model, Phpt == ηHP*Phpe);
    return model;
end;

function tess(model::InfiniteModel, data::Dict) # stationary battery pack
    # Thermal Energy Storage System
    t = model[:t];
    # Bucket model
    # Extract data
    Qtess = data["TESS"].Q; # Capacity [kWh]
    PtessMin = data["TESS"].PowerLim[1]; # Min power [kW]
    PtessMax = data["TESS"].PowerLim[2]; # Max power [kW]
    Ptess0 = data["TESS"].P0; # Initial State of Charge [p.u.]
    SoCtessMin = data["TESS"].SoCLim[1]; # Min State of Charge [p.u.]
    SoCtessMax = data["TESS"].SoCLim[2]; # Max State of Charge [p.u.]
    SoCtess0 = data["TESS"].SoC0; # Initial State of Charge [p.u.]
    ηtess = data["TESS"].η; # thermal efficiency
       
    # Add variables
    @variable(model, SoCtessMin ≤ SoCtess ≤ SoCtessMax, Infinite(t), start=SoCtess0); # State of Charge
    @variable(model, PtessMin ≤ Ptess ≤ PtessMax, Infinite(t), start=Ptess0); # Thermal power
    # Dummy variables for balance
    @variable(model, PtessMin ≤ PtessNeg ≤ -1e-4, Infinite(t), start=0); # Ptess^- in power
    @variable(model, 1e-4 ≤ PtessPos ≤ PtessMax, Infinite(t), start=Ptess0); # Ptess^+ out power

    # Bidirectional power flow, ensuring only export or import
    @constraint(model, PtessNeg + PtessPos == Ptess);
    @constraint(model, PtessNeg * PtessPos == 0);
    
    # Bucket model
    @constraint(model, ∂.(SoCtess, t) .== 1/ηtess*Ptess/Qtess);
    return model;
end;

function gridConn(model::InfiniteModel, sets, data::Dict) # power electronic interface
    Dtw = sets.dTime;
    t = model[:t];
    # Grid limits
    Pgmax = data["grid"].PowerLim[1]; # max power going out
    Pgmin = data["grid"].PowerLim[2]; # max power coming in
    Pg0 = data["grid"].P0; # initial condition
    
    # Add variables
    @variable(model, Pgmin ≤ Pg ≤ Pgmax, Infinite(t), start=Pg0);  # grid power
    
    # Dummy variables for balance
    @variable(model, Pgmin ≤ Pgneg ≤ 0, Infinite(t), start=0); # Pg^- in power
    @variable(model, 0 ≤ Pgpos ≤ Pgmax, Infinite(t), start=Pg0);  # Pg^+ out power
    
    # Bidirectional power flow, ensuring only export or import
    @constraint(model, Pgneg + Pgpos == Pg);
    @constraint(model, Pgneg * Pgpos == 0);
    
    return model;
end;

function gridThermal(model::InfiniteModel, sets, data::Dict) # thermal grid
    # Load data
    Dtw = sets.dTime;
    t = model[:t];
    loadTh = data["grid"].loadTh[1:length(Dtw)]; 
    Plt = JuMP.Containers.DenseAxisArray(loadTh, Dtw);
    
    # Extract params
    ηST = data["ST"].η; # conversion factor from Electric PV to thermal
    ηHP=data["HP"].η; # conversion factor from Electric to thermal heat pump
    MPPTmeas = data["SPV"].MPPT[1:length(Dtw)];
    Pst = JuMP.Containers.DenseAxisArray(ηST*MPPTmeas, Dtw);
        
    # Thermal Power balance
    model[:thBalance]=@constraint(model, [tw ∈ Dtw], Pst[tw]+model[:Phpe](tw).*ηHP+model[:Ptess](tw) .==  Plt[tw]);
    return model;
end;

function pei(model::InfiniteModel, sets, data::Dict) # power electronic interface
    # Extract data
    Dtw = sets.dTime;
    t = model[:t];
    # nEV = sets.nEV;
    
    # Load data
    loadElec = data["grid"].loadE[1:length(Dtw)];
    Ple = JuMP.Containers.DenseAxisArray(loadElec, Dtw);
       
    # MPPT measurement
    MPPTmeas = data["SPV"].MPPT[1:length(Dtw)]; 
    PpvMPPT = JuMP.Containers.DenseAxisArray(MPPTmeas, Dtw);
    
    # Power balance DC busbar
    model[:powerBalance]=@constraint(model, [tw ∈ Dtw], PpvMPPT[tw] + model[:Pg](tw) .== Ple[tw] + model[:Phpe](tw));
    return model;
end;

function costFunction(model, sets, data) # Objective function
    t = model[:t];
    Dtw = sets.dTime; Δt=Dtw[2]-Dtw[1];
    Pg0 = data["grid"].P0;
    # Energy prices
    @variable(model, -1e19 ≤ cgrid ≤ 1e19, Infinite(t), start = Pg0*Δt*data["grid"].λ[1,1])
    # Test with fixed prices
    @objective(model, Min, -∫(model[:Pgpos]*Δt*4 - model[:Pgneg]*Δt*2,t));
    return model;
end;

function solvePolicies(optimizer,
        sets,
        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
    tend=sets.dTime[end]
    
    model = InfiniteModel(optimizer) # create model
    # Ipopt attributes
    set_optimizer_attribute(model, "tol", 1e-4)
    set_optimizer_attribute(model, "max_iter", 1000)
    
    # define cont t
    @infinite_parameter(model, t in [0, tend], num_supports = length(sets.dTime), derivative_method = FiniteDifference(Central()) )
    
    # Add devices
    # thermal
    model = heatpump(model, data);
    model = tess(model, data);
    
    # # grids and balances
    model = gridConn(model, sets, data) # grid connection
    model = gridThermal(model, sets, data) # thermal balance
    model = pei(model, sets, data) # electrical balance
    # # Add objective function
    model = costFunction(model, sets, data)
    
    optimize!(model)
    return model;
end;

model=solvePolicies(Ipopt.Optimizer, s, EMSData)

An InfiniteOpt Model
Minimization problem with:
Finite Parameters: 386
Infinite Parameter: 1
Variables: 590
Derivative: 1
Measure: 1
Objective function type: GenericAffExpr{Float64, GeneralVariableRef}
`GeneralVariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 590 constraints
`GenericQuadExpr{Float64, GeneralVariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 3 constraints
`GeneralVariableRef`-in-`MathOptInterface.LessThan{Float64}`: 590 constraints
`GenericAffExpr{Float64, GeneralVariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 390 constraints
Names registered in the model: Pg, Pgneg, Pgpos, PhpNeg, PhpPos, Phpe, Ptess, PtessNeg, PtessPos, SoCtess, cgrid, powerBalance, t, thBalance, λbuy, λsell
Optimizer model backend information: 
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

The problem is that:

1. optimize!(model) never works. It runs for circa 2hs and goes nowhere.
2. Both gridThermal() and pei() introduce a weird behaviour in the constraints.
   When querying the containers as model[:thBalance] or model[:powerBalance] they return

model[:thBalance]
1-dimensional DenseAxisArray{InfOptConstraintRef,1,...} with index sets:
    Dimension 1, [0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25  …  45.75, 46.0, 46.25, 46.5, 46.75, 47.0, 47.25, 47.5, 47.75, 48.0]
And data, a 193-element Vector{InfOptConstraintRef}:
 3 Phpe(0) + Ptess(0) == 0.0
 3 Phpe(0.25) + Ptess(0.25) == 0.0
 3 Phpe(0.5) + Ptess(0.5) == 0.0
 3 Phpe(0.75) + Ptess(0.75) == 0.0
 3 Phpe(1) + Ptess(1) == 0.0
 3 Phpe(1.25) + Ptess(1.25) == 0.0
 3 Phpe(1.5) + Ptess(1.5) == 0.0

But when I

print(model)

I get those dummy/phantom constrains, that are accounted in the model, that repeat the information inside the first general constraint of the @variable.

[pulsipher]

Welcome @DarioSlaifsteinSk!

Interesting model, let me do my best to answer your questions.

With regard to the "dummy" constraints, this has to do with defining point variables which you are currently doing implicitly in your balance constraints. Point variables inherit all the properties of their parent infinite variable (e.g., bounds), but these can be changed to differ from their corresponding infinite variable. The good news here is that although the simple printing routine shows these point variable bounds as extra constraints (due to printing all the info of each point variable), the relationship between point variables and infinite variables is rigorously considered when the model is transcribed such that the model actually solved by Ipopt doesn't actually include these "redundant" constraints. However, the argument can be made to make this more transparent when printing which relates to #148.

Your balance constraints don't exactly do what I think you intend. Consider for instance:

model[:thBalance]=@constraint(model, [tw ∈ Dtw], Pst[tw]+model[:Phpe](tw).*ηHP+model[:Ptess](tw) .==  Plt[tw])

Here, I believe Dtw is the discrete support set to be used over the domain of the infinite parameter t. The problem with the current syntax is that you implicitly create many point variables with model[:Phpe](tw) and model[:Ptess](tw). This creates a lot of unnecessary place holder point variables to be kept track of and will be inefficiently discretized since this makes a lot of point constraints and instead of a single infinite constraint. I imagine the motivation for writing it this way is because of the discrete-time valued parameters. These should be lifted to continuous time using @parameter_function:

using Interpolations
Pst_interp = linear_interpolation(Dtw, Pst)
@parameter_function(model, Pst_cont == (t) -> Pst_interp(t)) # make InfiniteOpt compatible Pst(t) work for arbitrary times
Plt_interp = linear_interpolation(Dtw, Plt)
@parameter_function(model, Plt_cont == (t) -> Plt_interp(t)) # make InfiniteOpt compatible Plt(t) work for arbitrary times

Then we can rewrite the balance constraint in continuous time:

@constraint(model, thBalance, Pst_cont + model[:Phpe] .* ηHP + model[:Ptess] .==  Plt_cont)

I would also redefine your infinite parameter t to:

@infinite_parameter(model, t in [0, tend], supports = sets.dTime, derivative_method = FiniteDifference(Backward()) )

So you ensure you are using the supports you include in sets.dTime and use a more stable Euler scheme.

I am also not quite sure why you have Δt inside your integral.

As for why it is hanging at optimize!, I cannot run the model myself since I don't have the csv files you use. Either it is hanging when Ipopt was called or it is getting stuck in the discretization transformation. This is a relatively small model, so I wouldn't expect it to get stuck at the discretization step, but fixing the balance constraints should definitely help. Could you run:

build_optimizer_model!(model)

and see if that runs with the suggested changes. If that works without issue then, please share the Ipopt output when you optimize!.

[DarioSlaifsteinSk]

Worked as a charm! Some plots

As you said the whole problem came from the use of the discrete time data. A couple of hours later after the post I found @parameter_function but I was struggling with the implementation (similar to #207).
I agree that we can use Interpolations.jl but I was trying to make a Zero Order Hold object or something like that. The closest I could get was the following:

yp=hcat(data["grid"].λ[1:length(Dtw), 1])
λbuy = parameter_function(t, name = "buyPrice") do t_supp
                    if t_supp <= Dtw[1]
                        return yp[1]
                    end
                    if t_supp >= Dtw[end]
                        return yp[end]
                    end
                    k = 2
                    while t_supp > Dtw[k]
                        k += 1
                    end
                    return yp[k-1]
                 end

But this is "too" local and I´d prefer the syntax you did with Interpolations.jl.
I´ll check the rest of the package and see what´s the difference in results (which is actually cool thing in itself)

PS: the timestep $\Delta t$ was a typo

[pulsipher]

I am glad it works! I opened #302 to track adding to the documentation to make this more obvious in the future.

For a "zero order hold object", the most efficient way would be something like:

function zero_order_time_index(ts, t)
    inds = searchsorted(ts, t)
    if !isempty(inds) || inds.start == 1
        return inds.start
    elseif inds.start > length(ts)
        return length(ts)
    else # we don't have an exact match so let's choose the earlier time index
        return inds.stop
    end
end

@parameter_function(model, λbuy == (t) -> yp[zero_order_time_index(Dtw, t)])

Note, I haven't run the above so it may have a minor typo/bug.

[DarioSlaifsteinSk]

Thanks for the function!
The upgrade on the docs would be great! Just a simple implementation like this would be great. Probably the easiest place to put it is in the InfiniteOpt short course in the examples or something like that.
Using the measurement of a given plant to balance it out or something (which is what I´m doing)

[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,\ n}$.
Then the total power in the EV is the power charged and the power for driving following

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

$\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 continuos 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:

I have literally no idea on how to solve it. Can I change something on how its transcribed? Or is it the definition of $\gamma$??

[pulsipher]

See #304