Pyomo: Minimize for Max Value in Vector

I am optimizing the behavior of battery storage combined with solar PV to generate the highest possible revenue stream. I now want to add one more revenue stream: Peak Shaving (or Demand Charge Reduction)

My approach is as follows:

• Next to the price per kWh, an industrial customer pays for the maximal amount of power (kW) he was drawing from the grid in one period (i=1:end), so called demand charges
• This maximum amount is found in the vector P_Grid = P_GridLoad (energy self-consumed from the grid) + P_GridBatt (energy used to charge the battery)
• There exists a price vector which tells the price per kW for all points in time
• I now want to generate a vector P_GridMax that is zero for all points in time but the moment when the maximal value of P_Grid occurs (then it equals max(P_Grid).
• Thus, the vector P_GridMax consists of zeros and one nonzero element (not more!)
• In doing so, I can now multiply this vector with the price vector, sum up over all points in time and receive the billed demand charges
• By including this vector into the objective of my model I can minimize these charges

Now, does anybody see a solution for how to formulate such a constraint (P_GridMax)? I already updated my objective function and defined P_Grid. Any other approach would also be welcome.

This is the relevant part of my model, with P_xxx = power flow vectors, C_xxx = price vectors, ...

m.P_Grid = Var(m.i_TIME, within = NonNegativeReals)
m.P_GridMax = Var(m.i_TIME, within = NonNegativeReals)


# Minimize electricity bill
def Total_cost(m):
    return ... + sum(m.P_GridMax[i] * m.C_PowerCosts[i] for i in m.i_TIME) - ...
m.Cost = Objective(rule=Total_cost)


## Peak Shaving constraints
def Grid_Def(m,i):
    return m.P_Grid[i] = m.P_GridLoad[i] + m.P_GridBatt[i]
m.Bound_Grid = Constraint(m.i_TIME,rule=Grid_Def)

def Peak_Rule(m,i):
    ????
    ????
    ????
    ????
m.Bound_Peak = Constraint(m.i_TIME,rule=Peak_Rule)

Thank you very much in advance! Please be aware that I have very little experience with python/pyomo coding, I would really appreciate you giving extensive explanations :)

Best, Mathias

Here is one way to do this:

introduce a binary helper variable ismax[i] for i in i_TIME. This variable is 1 if the maximum is obtained in period i and 0 otherwise. Then obviously you have a constraint sum(ismax[i] for i in i_TIME) == 1: the maximum must be attained in exactly one period.

Now you need two additional constraints:

1. if ismax[i] == 0 then P_GridMax[i] == 0.
2. if ismax[i] == 1 then for all j in i_TIME we must have P_GridMax[i] >= P_GridMax[j].

The best way to formulate this would be to use indicator constraints but I don't know Pyomo so I don't know whether it supports that (I suppose it does but I don't know how to write them). So I'll give instead a big-M formulation.

For this formulation you need to define a constant M so that P_Grid[i] can not exceed that value for any i. With that the first constraint becomes

P_GridMax[i] <= M * ismax[i]

That constraint forces P_GridMax[i] to 0 unless ismax[i] == 1. For ismax[i] == 1 it is redundant. The second constraint would be for all j in i_TIME

P_GridMax[i] + M * (1 - ismax[i]) >= P_Grid[j]

If ismax[i] == 0 then the left-hand side of this constraint is at least M, so by the definition of M it will be satisfied no matter what the value of P_GridMax[i] is (the first constraint forces P_Grid[i] == 0 in that case). For ismax[i] == 1 the left-hand side of the constraint becomes just P_GridMax[i], exactly what we want.