Pyomo - Stock problem v1

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

pip install pyomo



In [ ]:

    
%matplotlib inline

import matplotlib.pyplot as plt



In [ ]:

    
from pyomo.environ import *

Version 1 (P1)

$x_t < 0$ = buy $|x|$ at time $t$
$x_t > 0$ = sell $|x|$ at time $t$
$s_t$ = the battery level at time $t$



In [ ]:

    
# Cost of energy on the market

#cost = [10, 30, 20]          # -> -100, 100, 0
#cost = [10, 30, 10, 30]      # -> [-100.,  100., -100.,  100.]
#cost = [10, 30, 10, 30, 30]  # -> [-100.,  100., -100.,  100.,    0.]
#cost = [10, 20, 30, 40]      # -> [-100.,  0., 0.,  100.]
#cost = [10, 30, 20, 50]
price = [10, 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 [ ]:

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

model.x = Var(T)
model.s = Var(T)


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

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

########

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[t]
    else:
        return Constraint.Skip

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


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

model.constraint_xt_inf = Constraint(T, rule=constraint_decision_inf)


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

model.constraint_xt_sup = Constraint(T, rule=constraint_decision_sup)

########

model.pprint()

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

opt = SolverFactory('glpk')

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

model.display()

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

Version 2 (P2)

$$ \begin{align} \max_{x_1,x_2,x_3} & \quad p_1 x_1 + p_2 x_2 + p_3 x_3 \quad\quad\quad\quad \text{(P1)} \\ \text{s.t.} & \quad s_0 = 0 \\ & \quad s_t = s_{t-1} - x_t \\ & \quad x_t \geq -(s_{\max} - s_{t-1}) \quad \quad ~ \text{can't buy more than the "free space" of the battery} \\ & \quad x_t \leq s_{t-1} \quad \quad \quad \quad \quad \quad \text{can't sell more than what have been stored in the battery} \\ & \\ & \\ \Leftrightarrow \min_{x_1,x_2,x_3} & \quad - p_1 x_1 - p_2 x_2 - p_3 x_3 \\ \text{s.t.} & \quad s_0 = 0 \\ & \quad s_t = s_{t-1} - x_t \\ & \quad x_t \geq -(s_{\max} - s_{t-1}) \Leftrightarrow x_t \geq -s_{\max} + s_{t-1} \\ & \quad x_t \leq s_{t-1} \\ & \\ & \\ \Leftrightarrow \min_{x_1,x_2,x_3} & \quad - p_1 x_1 - p_2 x_2 - p_3 x_3 \\ \text{s.t.} & \quad s_1 = -x_1 \\ & \quad x_1 \geq -s_{\max} \\ & \quad x_1 \leq 0 \\ & \quad s_2 = s_1 - x_2 = -x_1 - x_2 \\ & \quad x_2 \geq -s_{\max} + s_1 \Leftrightarrow x_2 \geq -s_{\max} - x_1 \\ & \quad x_2 \leq s_1 \Leftrightarrow x_2 \leq -x_1 \\ & \quad s_3 = s_2 - x_3 = -x_1 - x_2 - x_3 \\ & \quad x_3 \geq -s_{\max} + s_2 \Leftrightarrow x_3 \geq -s_{\max} - x_1 - x_2 \\ & \quad x_3 \leq s_2 \Leftrightarrow x_3 \leq -x_1 - x_2 \\ & \\ & \\ \Leftrightarrow \min_{x_1,x_2,x_3} & \quad - p_1 x_1 - p_2 x_2 - p_3 x_3 \\ \text{s.t.} & \quad x_1 \geq -s_{\max} \\ & \quad x_1 \leq 0 \\ & \quad x_2 \geq -s_{\max} - x_1 \\ & \quad x_2 \leq -x_1 \\ & \quad x_3 \geq -s_{\max} - x_1 - x_2 \\ & \quad x_3 \leq -x_1 - x_2 \\ & \\ & \\ \Leftrightarrow \min_{x_1,x_2,x_3} & \quad - p_1 x_1 - p_2 x_2 - p_3 x_3 \\ \text{s.t.} & \quad x_1 \geq -s_{\max} \\ & \quad x_1 \leq 0 \\ & \quad x_1 + x_2 \geq -s_{\max} \\ & \quad x_1 + x_2 \leq 0 \\ & \quad x_1 + x_2 + x_3 \geq -s_{\max} \\ & \quad x_1 + x_2 + x_3 \leq 0 \\ & \\ & \\ \Leftrightarrow \min_{x_1,x_2,x_3} & \quad - p_1 x_1 - p_2 x_2 - p_3 x_3 \quad\quad\quad\quad \text{(P2)} \\ \text{s.t.} & \quad -x_1 \leq s_{\max} \\ & \quad x_1 \leq 0 \\ & \quad -x_1 - x_2 \leq s_{\max} \\ & \quad x_1 + x_2 \leq 0 \\ & \quad -x_1 - x_2 - x_3 \leq s_{\max} \\ & \quad x_1 + x_2 + x_3 \leq 0 \\ & \\ & \\ \end{align} $$$$ \color{\red}{ \boldsymbol{c} = \begin{pmatrix} -p_1 \\ -p_2 \\ -p_3 \end{pmatrix} } \quad \color{\orange}{ \boldsymbol{A} = \begin{pmatrix} -1 & 0 & 0 \\ 1 & 0 & 0 \\ -1 & -1 & 0 \\ 1 & 1 & 0 \\ -1 & -1 & -1 \\ 1 & 1 & 1 \end{pmatrix} } \quad \color{\green}{ \boldsymbol{b} = \begin{pmatrix} s_{\max} \\ 0 \\ s_{\max} \\ 0 \\ s_{\max} \\ 0 \end{pmatrix} } \quad \color{\purple}{ \boldsymbol{B} = \begin{pmatrix} - & - \\ - & - \\ - & - \end{pmatrix} } $$



In [ ]:

    
T = list(range(len(price)))       # num decision variables
M = list(range(len(price) * 2))       # num decision variables



In [ ]:

    
plt.plot(T, 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 [ ]:

    
p = -np.array(price)



In [ ]:

    
A = np.repeat(np.tril(np.ones(len(price))), 2, axis=0)
A[::2, :] *= -1
A



In [ ]:

    
b = np.zeros(A.shape[0])
b[::2] = stock_max
b



In [ ]:

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

model.x = Var(T)


def objective_fn(model):
    return sum(p[t] * model.x[t] for t in T)

model.obj = Objective(rule=objective_fn)


def constraint_fn(model, m):
    return sum(A[m][t] * model.x[t] for t in T) <= b[m]

model.constraint = Constraint(M, rule=constraint_fn)  

model.pprint()

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

opt = SolverFactory('glpk')

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

model.display()

print()
print("Optimal solution: ", [value(model.x[t]) for t in T])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail