How to maximise a function with non-differentiable constraints while integrating over a time series data series with linear programming?

Question

In this problem, I need to pick the optimal charging and discharging times and durations of an electric car-sharing battery to maximize income - not just for a single (dis)charging cycle, but a whole bunch, with some car-sharing bookings that must work in between.

The curve in the graph is the time series data of the energy price $P(t)$. The red box is a car-sharing booking, which pauses the (dis)charging. It discharges the battery by $\Delta \text{SoC}_{b_n}$. Left of the booking is a discharging period starting at $I_n$, lasting for $\Delta i_n$. On the right is a charging time at $S_n$ for the duration of $\Delta s_n$.

The red line signifies the state of the battery charge (state of charge, SoC), which follows the charging and discharging cycle and the usage of the car.

My objective function is

$$ f_{objective} = \sum_{n=0}^{m} \int\limits_{I_n}^{I_n + \Delta i_n} P(t) \,dt - \int\limits_{S_n}^{S_n + \Delta s_n} P(t) \,dt $$

Here, I directly integrate the price data $P(t)$ with my variables, setting the intervals and the times.

The constraints are primarily for SoC: At no time can the battery be fuller than 100% and emptier than 0%. I try to express that with the upper bound of the sum of $1 .. m$. ($c$ is the charge/time factor)

$$ 0 \leq \sum_{n=0}^{1..m} -\Delta i_n \cdot c - \Delta \text{SoC}_{b_n} + \Delta s_n \cdot c \leq 1 $$

and even the individual charges can never leave the SoC-interval of 0..1:

$$ 0 \leq \Delta i_n \cdot c \leq 1, \quad 0 \leq \Delta \text{SoC}_{b_n} \leq 1, \quad 0 \leq \Delta s_n \cdot c \leq 1 $$ The (dis)charging and booking times can not overlap. $$ I_n + \Delta i_n < B_n, \quad B_n + \Delta b_n < S_n, \quad S_n + \Delta s_n < I_{n+1} $$

I wrote some python code for this and ran into problems with the integration. I have difficulty understanding the pyomo Integral function and how to use it, especially with my time series data. I tried to integrate with NumPy functionality, but that fails. Can you please help me make this work? I am a pyomo noob and still struggle with concepts and terminology.

import numpy as np
from pyomo.environ import *
from pyomo.opt import SolverFactory

def make_demo_data():
    x = np.linspace(0, 2 * np.pi, 100)
    y = np.sin(x) + 1.5 -x*.025  # add offset to make sure all values are positive
    return np.array([x, y])

def integrate(data_a, x, delta_x):
    x = x.values
    delta_x = delta_x.values
    scaling_factor = 100 / (4*np.pi)
    x_start = x
    x_end = x + delta_x
    x_data = np.linspace(x_start, x_end, 5)
    y_interp = np.interp(x_data, data_a[0], data_a[1])
    earning = np.trapz(y_interp, x_data * scaling_factor)
    return earning

data = make_demo_data()

model = ConcreteModel(name="(DischargePlan)")

# define variables
model.i = Var(domain=NonNegativeReals, name="i")
model.delta_i = Var(domain=NonNegativeReals, name="delta_i")
model.s = Var(domain=NonNegativeReals, name="s")
model.delta_s = Var(domain=NonNegativeReals, name="delta_s")

model.objectives = Objective(expr=integrate(data, model.i, model.delta_i) - integrate(data, model.s, model.delta_s), sense=maximize)

model.constraints = ConstraintList()
model.constraints.add(model.i + model.delta_i <= np.pi)
model.constraints.add(model.s + model.delta_s <= 2*np.pi)
model.constraints.add(1 - model.delta_i + model.delta_s <= 1)
model.constraints.add(1 - model.delta_i  <= 1)
model.constraints.add(1 - model.delta_i + model.delta_s > 0)
model.constraints.add(1 - model.delta_i > 0)

solver = SolverFactory('CPLEX')
solver.solve(model)
print(model.objectives())
print(f"i = {model.i.values}, delta_i = {model.delta_i.values}, s = {model.s.values}, delta_s = {model.delta_s.values}")

Answer 1

This is an attempt at a solution. It has several shortcomings.

• it is a discrete solution
• it does not have bookings yet
• it should handle different buying and selling prices, but it does not.

On the upside

• it does buy and sell
• it handles a few thousand variables nicely

Some of this is most likely a z3 limitation. I will try with another optimizer next.

import numpy as np
from z3 import *

CHARGING_CONST = 1/3
START_VAL = 1.0
NUM = 20

def make_demo_data():
    x = np.linspace(0, (4 * np.pi), NUM)
    y = np.sin(x) +1.5 -.025 * x
    z = np.sin(x) * 3 -.027 * x
    data = np.zeros(NUM, dtype=([("time", float), ("price", float), ("price_co2", float)]))
    data["time"] = x
    data["price"] = y
    data["price_co2"] = z
    return data

data = make_demo_data()
soc = []

o = z3.Optimize()
soc.append(z3.Real("soc_0"))
o.add(soc[0] <= START_VAL + CHARGING_CONST)
o.add(soc[0] >= START_VAL - CHARGING_CONST)
o.add(soc[0] >= 0)
o.add(soc[0] <= 1)
f_obj = 0
f_obj += data["price"][0] * (START_VAL - soc[0])

for cnt in range(1, NUM):
    soc.append(z3.Real(f"soc_{cnt}"))
    o.add(soc[cnt] <= soc[cnt - 1] + CHARGING_CONST)
    o.add(soc[cnt] >= soc[cnt - 1] - CHARGING_CONST)
    o.add(soc[cnt] >= 0)
    o.add(soc[cnt] <= 1)
    f_obj += (soc[cnt - 1] - soc[cnt]) *\
        z3.If(soc[cnt - 1] >= soc[cnt],
           data["price"][cnt],
           data["price_co2"][cnt])
# data["price"][cnt]

o.add(soc[-1] == START_VAL)
o.maximize(f_obj)
print(o.check())
m=o.model()

for cnt in range(NUM):
    print(f"""{cnt}: {m.evaluate(soc[cnt]).as_decimal(2)}""")

print(f"Result: {m.evaluate(f_obj).as_decimal(2).strip('?')}")