# Pyomo - Stock problem v1

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

```
pip install pyomo
```


```python
%matplotlib inline

import matplotlib.pyplot as plt
```


```python
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$

$$
\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}
$$


```python
# 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
```


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


```python
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();
```


```python
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}
}
$$


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


```python
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();
```


```python
p = -np.array(price)
```


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


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


```python
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
```

5745ec4d7918c33193d2e1b54dfabced3b8c81fc
ipynb
Jupyter Notebook
nb_dev_python/python_pyomo_stock_problem_1.ipynb
jdhp-docs/python-notebooks
91a97ea5cf374337efa7409e4992ea3f26b99179
MIT
2017-05-03T12:23:36.000Z
2020-10-26T17:30:56.000Z
nb_dev_python/python_pyomo_stock_problem_1.ipynb
jdhp-docs/python-notebooks
91a97ea5cf374337efa7409e4992ea3f26b99179
MIT
nb_dev_python/python_pyomo_stock_problem_1.ipynb
jdhp-docs/python-notebooks
91a97ea5cf374337efa7409e4992ea3f26b99179
MIT
2020-10-26T17:30:57.000Z
2020-10-26T17:30:57.000Z
Qwen/Qwen-72B
1. YES
2. YES