Linear programming : stock problem

See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html



In [ ]:

    
%matplotlib inline

import matplotlib.pyplot as plt

import scipy.optimize

Scipy's syntax

Example for a problem of 2 dimensions:

$$ \begin{align} \min_{x_1,x_2} & \quad \color{\red}{c_1} x_1 + \color{\red}{c_2} x_2 \\ \text{s.t.} & \quad \color{\orange}{A_{1,1}} x_1 + \color{\orange}{A_{1,2}} x_2 \leq \color{\green}{b_1} \\ & \quad \color{\orange}{A_{2,1}} x_1 + \color{\orange}{A_{2,2}} x_2 \leq \color{\green}{b_2} \\ & \quad \color{\purple}{B_{1,1}} \geq x_1 \geq \color{\purple}{B_{1,2}} \\ & \quad \color{\purple}{B_{2,1}} \geq x_2 \geq \color{\purple}{B_{2,2}} \\ \end{align} $$

$$ \color{\red}{ \boldsymbol{c} = \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} } \quad \color{\orange}{ \boldsymbol{A} = \begin{pmatrix} A_{1,1} & A_{1,2} \\ A_{2,1} & A_{2,2} \end{pmatrix} } \quad \color{\green}{ \boldsymbol{b} = \begin{pmatrix} b_1 \\ b_2 \end{pmatrix} } \quad \color{\purple}{ \boldsymbol{B} = \begin{pmatrix} B_{1,1} & B_{1,2} \\ B_{2,1} & B_{2,2} \end{pmatrix} } $$

$$ \text{scipy.optimize.linprog}(\color{\red}{\boldsymbol{c}}, ~ \color{\orange}{\boldsymbol{A}}, ~ \color{\green}{\boldsymbol{b}}, ~ \color{\purple}{\boldsymbol{B}}) $$

Stock problem



In [ ]:

    
plt.plot([1, 2, 3], [10, 30, 20], "o-")
plt.xlabel("Unit of time (t)")
plt.ylabel("Price of one unit of energy (c)")
plt.title("Cost of energy on the market")
plt.show();

TODO: put pictures and explanations from my PhD

$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 \\ \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 \\ \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 [ ]:

    
# Price of energy on the market
price = [10, 30, 20]

plt.plot(price);



In [ ]:

    
stock_max = 100      # battery capacity

Hand write the model



In [ ]:

    
# Coefficients of the linear objective function to be minimized
p = -np.array(price)

# 2-D array which, when matrix-multiplied by x, gives the values of the upper-bound inequality constraints at x.
A = [[-1,  0,  0],
     [ 1,  0,  0],
     [-1, -1,  0],
     [ 1,  1,  0],
     [-1, -1, -1],
     [ 1,  1,  1]]

# 1-D array of values representing the upper-bound of each inequality constraint (row) in A.
b = [stock_max, 0, stock_max, 0, stock_max, 0]

# Sequence of (min, max) pairs for each element in x, defining the bounds on that parameter.
# Use None for one of min or max when there is no bound in that direction.
# By default bounds are (0, None) (non-negative).
# If a sequence containing a single tuple is provided, then min and max will be applied to all variables in the problem.
x0_bounds = (None, None)
x1_bounds = (None, None)
x2_bounds = (None, None)
bounds = (x0_bounds, x1_bounds, x2_bounds)

scipy.optimize.linprog(p, A_ub=A, b_ub=b, bounds=bounds)

Automatically make the model



In [ ]:

    
# Cost of energy on the market

#price = [10, 30, 20]          # -> -100, 100, 0
#price = [10, 30, 10, 30]      # -> [-100.,  100., -100.,  100.]
#price = [10, 30, 10, 30, 30]  # -> [-100.,  100., -100.,  100.,    0.]
#price = [10, 20, 30, 40]      # -> [-100.,  0., 0.,  100.]
price = [10, 30, 20, 50]
price = [10, 30, 20, 50]

plt.plot(price);



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 [ ]:

    
bounds = tuple((None, None) for _p in price)
bounds



In [ ]:

    
%%time
res = scipy.optimize.linprog(p, A_ub=A, b_ub=b, bounds=bounds)   # , method='revised simplex'



In [ ]:

    
res



In [ ]:

    
res.x.round(decimals=2)