Battery Problem¶

We have a vehicle which needs to pass a known test track. You are given values $Preq(t)$ for each time $t=1,...,T$ which the vehicle's wheels needs to exert in order to successfully pass this test track.

The vehicle has a combustion engine, a motor/generator connected to a battery, and a friction break. Motor/generator can act either as a motor when it uses energy stored in the battery to power the wheels, or as a generator when it extracts the power from the wheels (regenerative breaking) or the engine to store it in the battery.

$Preq(t)$ is positive when the wheels need power, e.g. vehicle is ascending a hill or accelerating: power for the wheels needs to be provided by the combustion engine and/or the motor/generator which extracts the power stored in the battery. When $Preq(t)$ is negative (e.g. because vehicle is descending a hill), the power is extracted by the motor/generator which can store it in the battery and/or by the friction break.

Power is conserved, i.e., at each time $t$, we have $Preq(t) = Peng(t) + Pmg(t) - Pbr(t)$ where $0≤Peng(t)≤Peng_{max}$ is power produced by the combustion engine, $Pmg_{min} ≤ Pmg(t) ≤ Pmg_{max}$ is power produced by the motor/generator (can be also negative if motor/generator absorbs power to charge the battery) and $Pbr(t)≥0$ is the power absorbed by the friction break. See data bellow for the definition of the constants.

For every $t=1,...,T+1$, the energy $E(t)$ in the battery has to be between $0$ and $Ebatt_{max}$ representing the power stored in empty and full battery respectively. Moreover, we need to take into account the charging and discharging of the battery: we have $E(t+1) = E(t) - Pmg(t) - η|Pmg(t)|, for t=1,...,T,$ see variable eta in the data. The term with coefficient η represents the energy lost due to the inefficiency of the battery and motor/generator. We also require $E(T+1)=E(1)$ to make a fair comparison with a non-hybrid vehicle which has no battery.

The objective is to minimize the total fuel consumption of the vehicle over time, where the consumption at time $t$ is given by formula $Peng(t) + γ(Peng(t))^2$, i.e., a quadratic function. See data for the value of $γ$.

Task 1: formulate this minimization problem as a convex program (10 points)¶

It is almost convex, just one set of constraints is problematic. Hint: Try to relax it, only one inequality is important, but provide an explanation (e.g. in the comments of your code) why your formulation is equivalent.

Task 2: solve your convex program with cvxpy library (10 points)¶

Link to the library: https://www.cvxpy.org/. Present the solution using a clear and well explained plot generated by matplotlib which was also used in the previous assignment.

Task 3: comparison with a battery-less car (5 points)¶

Change $Ebatt_{max}$ to $0$ and solve your program again to see how does it affect the power consumption.

Task 4: handle glitches (5 points)¶

You may find out that although your program is equivalent, the solution found by the solver does not fulfill all the relaxed constraints with equality. This might happen during long periods of breaking when there is a lot of opportunities to charge the battery and sometimes the solver may decide to waste part of the power available to charge the battery (because it may charge it fully in the following time steps). There are two options how to proceed:

  • Postprocess the solution to get another solution with the same objective value where the relaxed constraints are preserved with equality. Such solution must exist since your program is equivalent to the original one. Hint: just move all the power from motor/generator to the battery and, if the battery is full, move it to the friction break.

  • Add a small term to the objective which will discourage such situation: for every $t=1,...,T$, you can add a term $ε·max(0, -Pmg(t))$ for some small positive $ε$ in order to discourage absorbing power by the motor/generator if it is not going to be used for charging the battery (because it is cheaper to absorb the power using the friction break) In case you decide to modify the objective, compare the power consumption achieved by your program to the power consumption achieved by the program with the original objective (to show that they are very similar with your choice of $ε$.

Data¶

Include the following code in your solution which generates data for you. The array Preq contains power requirements in all time steps. Bellow, are the bounds specifying the parameters of the engine, motor/generator, capacity of the battery, coefficient eta of the inefficiency of charging/discharging and coefficient gamma in the objective function.

import numpy as np
a=[0.5, -0.5, 0.2, -0.7, 0.6, -0.2, 0.7, -0.5, 0.8, -0.4]
l=[40, 20, 40, 40, 20, 40, 30, 40, 30, 60]
Preq=np.arange(a[0],a[0]*(l[0]+0.5),a[0])

for i in range(1, len(l)):
    Preq=np.r_[ Preq, np.arange(Preq[-1]+a[i],Preq[-1]+a[i]*(l[i]+0.5),a[i]) ]

T = sum(l)

Peng_max = 20.0
Pmg_min = -6.0
Pmg_max = 6.0
Ebatt_max = 100.0
eta = 0.1
gamma = 0.1

A note on precision¶

LP solvers which we used so far always provided a precise feasible solutions. This is not always the case with the solvers for convex optimization. If you want better precision, see "eps" parameters in Solver options. However, default values are good enough for the purpose of this assignment.

In [ ]:
import numpy as np
import cvxpy as cp

#######
# DATA, do not change this part!
#######
a=[0.5, -0.5, 0.2, -0.7, 0.6, -0.2, 0.7, -0.5, 0.8, -0.4]
l=[40, 20, 40, 40, 20, 40, 30, 40, 30, 60]
Preq=np.arange(a[0],a[0]*(l[0]+0.5),a[0])
for i in range(1, len(l)):
    Preq=np.r_[ Preq, np.arange(Preq[-1]+a[i],Preq[-1]+a[i]*(l[i]+0.5),a[i]) ]

T = sum(l)

Peng_max = 20.0
Pmg_min = -6.0
Pmg_max = 6.0
eta = 0.1
gamma = 0.1
#####
# End of DATA part
#####

# Implement the following functions
# they should return a dictionary retval such that
# retval['Peng'] is a list of floats of length T such that retval['Peng'][t] = P_eng(t+1) for each t=0,...,T-1
# retval['Pmg'] is a list of floats of length T such that retval['Pmg'][t] = P_mg(t+1) for each t=0,...,T-1
# retval['Pbr'] is a list of floats of length T such that retval['Pbr'][t] = P_br(t+1) for each t=0,...,T-1
# retval['E'] is a list of floats of length T+1 such that retval['E'][t] = E(t+1) for each t=0,...,T


def car_with_battery():
    Ebatt_max = 100.0
    retval = car(Ebatt_max)
    return retval

def car_without_battery():
    Ebatt_max=0
    retval = car(Ebatt_max)
    return retval


def car(Ebatt_max):
    # Variables
    Peng = cp.Variable(T)
    Pmg = cp.Variable(T)
    Pbr = cp.Variable(T)
    E = cp.Variable(T+1)

    # Create a small term to add to the objective value
    epsilon = 0

    # Constraints
    constraints = []

    # Objective value
    y = 0

    # The energy in the battery at the end of the test track is equal to the initial energy
    constraints.append(E[T] == E[0])

    for t in range(T):


        # The energy of the battery is always non-negative and below the maximum battery capacity
        constraints.append(E[t] >= 0)
        constraints.append(E[t] <= Ebatt_max)

        # The power produced by the engine is always non-negative and below the maximum engine power
        constraints.append(Peng[t] >= 0)
        constraints.append(Peng[t] <= Peng_max)

        #The power produced or absorbed by the motor/generator is above the minimum and below the maximum motor/generator power
        constraints.append(Pmg[t] >= Pmg_min)
        constraints.append(Pmg[t] <= Pmg_max)

        # The power absorbed by the friction is always non-negative
        constraints.append(Pbr[t] >= 0)

        # Conservation of power at each time step
        constraints.append(Preq[t] == Peng[t] + Pmg[t] - Pbr[t])

        # The energy in the battery at the next time step is equal to the
        # energy in the battery at the current step minus the power absorbed or
        # produced by the motor/generator minus the energy loss due to the
        # inefficiency of the battery and moto/generator. This constraint is
        # problematic, since it is not convex. In order to find a way around
        # this problem, we can change the "==" into "<=", by discarding the
        # possibility that we have ">=". Indeed, this formulation should be
        # equivalent, because E[t+1] should always be chosen as large as
        # possible to keep the battery as charged as possible. This allows the
        # car to satisfy the power requirements through the motor/generator
        # instead of the the engine, which is better since we are trying to
        # minimize fuel consumption. So, we reach
        # 'E[t+1] <= E[t]-Pmg[t]-eta*abs(Pmg[t])'. Then, we split the
        # inequality by exploiting the fact that for any real number 'x',
        # '|x|<y <=> -y < x < y'. Hence, I accordingly split my constraint into
        # two.
        constraints.append(E[t+1] - E[t] + Pmg[t] <= -eta*Pmg[t])
        constraints.append(-E[t+1] + E[t] - Pmg[t] >= -eta*Pmg[t])

        # Add a small term that discourages the possibility to not satisfy
        # relaxed constraints.
        x = cp.maximum(0, -Pmg[t])

        # Objective value
        y += cp.sum(Peng[t] + gamma*(Peng[t]**2)+epsilon*x)

    # Problem
    prob = cp.Problem(cp.Minimize(y), constraints)

    # Solve
    prob.solve(solver=cp.OSQP, eps_abs=1e-10, eps_rel=1e-10)

    retval = {
        'Peng': [float(Peng[t].value) for t in range(T)],
        'Pmg': [float(Pmg[t].value) for t in range(T)],
        'Pbr': [float(Pbr[t].value) for t in range(T)],
        'E': [float(E[t].value) for t in range(T+1)],
# =============================================================================
        'Sol':  prob.value
# =============================================================================
        }

    return retval





# Print optimal solutions to the two problems
print('Sol:', car_with_battery()['Sol'])
print('Sol:', car_without_battery()['Sol'])

# check for errors
def check_error(res):
    for t in range(T):
        error = abs(res['E'][t]-res['E'][t+1]-res['Pmg'][t]-eta*abs(res['Pmg'][t]))
        if error > 0.0003:
            print(error)
error_with_battery = check_error(car_with_battery())
error_without_battery = check_error(car_without_battery() )