Pyomo - Stock problem v2

Pyomo installation: see http://www.pyomo.org/installation

pip install pyomo



In [ ]:

    
%matplotlib inline

import matplotlib.pyplot as plt



In [ ]:

    
from pyomo.environ import *

$s_t$ = the battery level at time $t$

$$ \begin{align} \min_{x_{inj}, x_{ext}, x_{grid}} & \quad \sum_t p_t x_{inj,t} + \sum_t p_t x_{grid,t} \quad\quad\quad\quad \text{(P1)} \\ \text{s.t.} & \quad s_0 = 0 \\ & \quad s_t = s_{t-1} + x_{inj,t} - x_{ext,t} \\ & \quad x_{inj,t} \leq s_{\max} - s_{t-1} \quad\quad\quad \quad ~~ (L4) \\ & \quad x_{inj,t} \geq 0 \\ & \quad x_{ext,t} \leq s_{t-1} \quad\quad \quad \quad \quad \quad \quad (L6) \\ & \quad x_{ext,t} \geq 0 \\ & \quad x_{ext,t} + x_{grid,t} = \text{needs}_{t} \quad \quad \quad ~ (L8) \end{align} $$

L4: can't inject more than the "free space" of the battery
L6: can't extract more than what have been stored in the battery
L8: fulfill the needs at time t



In [ ]:

    
# Cost of energy on the market

price = [10, 20, 10, 20]
needs = [ 0, 30, 20, 50]

stock_max = 100           # battery capacity



In [ ]:

    
T = list(range(len(price) + 1))       # num decision variables



In [ ]:

    
plt.plot(T[:-1], price, "o-")
plt.xlabel("Unit of time (t)")
plt.ylabel("Price of one unit of energy (p)")
plt.title("Price of energy on the market")
plt.show();



In [ ]:

    
plt.plot(T[:-1], needs, "o-")
plt.xlabel("Unit of time (t)")
plt.ylabel("Needs")
plt.title("Needs")
plt.show();



In [ ]:

    
model = ConcreteModel(name="Stock problem v1")

model.x_inj  = Var(T, within=NonNegativeReals)
model.x_ext  = Var(T, within=NonNegativeReals)
model.x_grid = Var(T, within=NonNegativeReals)
model.s = Var(T, within=NonNegativeReals)


def objective_fn(model):
    return sum(price[t-1] * model.x_inj[t] + price[t-1] * model.x_grid[t] for t in T if t != 0)

model.obj = Objective(rule=objective_fn, sense=minimize)

########

model.constraint_s0 = Constraint(expr = model.s[0] == 0)

def constraint_stock_level(model, t):
    if t > 0:
        return model.s[t] == model.s[t-1] + model.x_inj[t] - model.x_ext[t]
    else:
        return Constraint.Skip

model.constraint_st = Constraint(T, rule=constraint_stock_level)


def constraint_x_inj_sup(model, t):
    if t > 0:
        return model.x_inj[t] <= stock_max - model.s[t-1]
    else:
        return Constraint.Skip

model.constraint_x_inj_sup = Constraint(T, rule=constraint_x_inj_sup)


def constraint_x_ext_sup(model, t):
    if t > 0:
        return model.x_ext[t] <= model.s[t-1]
    else:
        return Constraint.Skip

model.constraint_x_ext_sup = Constraint(T, rule=constraint_x_ext_sup)


def constraint_needs(model, t):
    if t > 0:
        return model.x_ext[t] + model.x_grid[t] == needs[t-1]
    else:
        return Constraint.Skip

model.constraint_needs = Constraint(T, rule=constraint_needs)

########

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)    # solves and updates instance

model.display()

print()
print("Optimal solution (x_inj):  ", [value(model.x_inj[t]) for t in T if t != 0])
print("Optimal solution (x_ext):  ", [value(model.x_ext[t]) for t in T if t != 0])
print("Optimal solution (x_grid): ", [value(model.x_grid[t]) for t in T if t != 0])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail