Pyomo: Minimize for Max Value in Vector

Question

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

score 3 · Answer 1 · answered Oct 16 '19 at 22:37

Another idea is that you don't actually need to index your P_GridMax variable with time.

If you're dealing with demand costs they tend to be fixed over some period, or in your case it seems that they are fixed over the entire problem horizon (since you're only looking for one max value).

In that case you would just need to do:

m.P_GridMax = pyo.Var(domain=pyo.NonNegativeReals)

def Peak_Rule(m, i):
    return m.P_GridMax >= m.P_Grid[i]

m.Bound_Peak = pyo.Constraint(m.i_TIME,rule=Peak_Rule)

if you're really set on multiplying your vectors element-wise you can also just create a new variable that represents that indexed product and apply the same principle to extract the max value.

score 1 · Accepted Answer · answered May 08 '19 at 19:02

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:

if ismax[i] == 0 then P_GridMax[i] == 0.
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.

Pyomo: Minimize for Max Value in Vector

2 Answers2