<jupyter_start><jupyter_text>Pyomo - Stock problem v1 Pyomo installation: see http://www.pyomo.org/installation```pip install pyomo```<jupyter_code>%matplotlib inline

import matplotlib.pyplot as plt
from pyomo.environ import *<jupyter_output><empty_output><jupyter_text>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$ $$\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}\end{align}$$<jupyter_code># 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
T = list(range(len(price) + 1))       # num decision variables
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();
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<jupyter_output><empty_output><jupyter_text>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}}$$<jupyter_code>T = list(range(len(price)))       # num decision variables
M = list(range(len(price) * 2))       # num decision variables
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();
p = -np.array(price)
A = np.repeat(np.tril(np.ones(len(price))), 2, axis=0)
A[::2, :] *= -1
A
b = np.zeros(A.shape[0])
b[::2] = stock_max
b
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<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1773040254237288, 'fraction_numerical': 0.04078389830508475, 'mean_word_length': 1.6271304347826088, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '90722', 'license': 'MIT', 'n_tokens_mistral': 3017, 'n_tokens_neox': 2803, 'n_words': 711, 'path': 'nb_dev_python/python_pyomo_stock_problem_1.ipynb', 'repo_name': 'jdhp-docs/python-notebooks', 'size': '7552'}