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.

sketch of the price data and the naming of my variables

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')
print(f"i = {model.i.values}, delta_i = {model.delta_i.values}, s = {model.s.values}, delta_s = {model.delta_s.values}")
  • $\begingroup$ In the integrate func, you use values for $\delta_i$. $ \delta_i$ is your variable to be used in the obj. Does Pyomo allow to use $ \delta_i$.values() before solving? $\endgroup$
    – Sutanu
    Feb 1 at 21:10
  • $\begingroup$ @Sutanu I don't know. The problem requires $\Delta i$ and $\Delta s$ to be used in some shape or form. Can you suggest a workaround? $\endgroup$ Feb 1 at 21:28
  • $\begingroup$ try writing the objective function without using delta_i instead of delta_i .values. Also I think, obj needs to be an expression, not returning value of a function, like obj = <<expr>>. Then in the Objective(expr=obj,..) $\endgroup$
    – Sutanu
    Feb 1 at 21:31
  • $\begingroup$ @Sutanu is it ok to use $I_n$ and $S_n$ in the objective function? $\endgroup$ Feb 1 at 21:36
  • 1
    $\begingroup$ try this one: math.stackexchange.com/questions/1116267/…. You need to use the derivative of the integral function, basically the func and solve it. Repeat for the values within integral bounds. $\endgroup$
    – Sutanu
    Feb 1 at 22:11

1 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 *

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()
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):
    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]

o.add(soc[-1] == START_VAL)

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('?')}")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.