3

I am using SCIPY to optimize a storage facility that uses forward prices for a deal term of 1 year. Gas can be injected and withdrawn from this facility, based on monthly spreads (e.g. March 21 vs May 20 spread) being high enough to cover the variable cost of operation. The attached picture represents the problem (the values here are arbitrary, don't match the values in the code; pic is just for concept)

enter image description here

The cells in blue are the "changing cells", volumes that SCIPY will adjust to maximize profits. The constraints need to be set up for each month separately. I get errors when I attempt to set up these constraints in SCIPY. Here's a reproducible version of the problem:

 import numpy as np
import scipy.optimize as opt

p= np.array([4, 5, 6.65, 12]) #p = prices
pmx = np.triu(p - p[:, np.newaxis]) #pmx = price matrix, upper triangular

q =np.triu(np.ones((4,4))) # q = quantity, upper triangular

def profit(q):
    profit = -np.sum(q.flatten() * pmx.flatten())
    return profit

bnds = (0,10)
bnds = [bnds for i in q.flatten()]

def cons1(q):
    np.sum(q,axis=1) -  10

#def cons2(q):
#    np.sum(q,axis=0) -  8

#con1 = {'type':'ineq','fun':cons1}
#con2 = {'type':'ineq','fun':cons2}
cons = [con1]    # using only  1 constraint (con1) to test the model

#sol = opt.minimize(profit,q,method='SLSQP', bounds= bnds,constraints = cons)
sol = opt.minimize(profit,q,method='SLSQP', bounds= bnds)
sol

The model runs fine when I exclude the constraints. When I add one of the constraints, the error I get is:

AxisError: axis 1 is out of bounds for array of dimension 1

I think this has to do with the the way I'm specifying the constraints....I'm not sure though. For the constraints, I do need to identify injections and withdrawals and set the constraints as shown in the picture. Help would be appreciated. Thanks!

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
user3654852
  • 753
  • 1
  • 5
  • 10
  • 1) The first argument `int` is not known in your program. Not really *reproducible*. 2) The first argument also has the name of some python-specific int-constructor. Not recommended 3) Your second argument is a 2d-array, where scipy [needs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) something with shape `(n,)` and not `(n, m)` – sascha May 06 '20 at 07:16
  • Hi - my error, I just updated the code and it runs fine (when I exclude the constraints). The price and quantity arrays are specified as shown in the picture., but are flattened inside the profit calculation. Constraints are the issue - I need to identity injection and withdrawals separately, and sum across different axis in the quantity matrix to set the constraints. Not sure how to do that. – user3654852 May 06 '20 at 13:16
  • No... `q` still needs to be 1d or every kind of strange stuff is allowed to happen. `np.triu` is not 1d. Running or not, it does not matter. There are probably > 20 questions here related to not working on *1d* optimization-vectors (and inconsistent hidden flattenings inside initial-vec calc, numerical-gradient calc or other components) – sascha May 06 '20 at 13:25
  • Thanks. I can use flatten in the q and p matrix directly, instead of in the profit function def.That does work. Any advice on the constraint specification? I guess I'll have to list each constraint; one for every injection, withdrawal pair in the q matrix. Depending on the deal term, this could be a long list, i was hoping to avoid that, by summing across the axis in the q matrix. – user3654852 May 06 '20 at 14:09
  • I guess I could build the constraint sequentially as a loop. I'm struggling with how to identify injection/wdral pairs without the matrix; summing down a month in the matrix identifies withdrawals, and summing across a row identifies an injection. That distinction vanishes in a 1 dimensional array. – user3654852 May 06 '20 at 14:27
  • 1
    You can do the inverse-flattening (`mat = my_1d_vec.reshape(my_2d_shape)`) in the first line (or even as part of a one-liner) to reuse your matrix-based original code. – sascha May 06 '20 at 14:34
  • AH! thanks for that tip...let me investigate it.Really appreciate your help here. – user3654852 May 06 '20 at 14:49

1 Answers1

1

As an alternative to Scipy.minimize.optimize, here is a solution with Python gekko.

import numpy as np
import scipy.optimize as opt
from gekko import GEKKO

p= np.array([4, 5, 6.65, 12]) #p = prices
pmx = np.triu(p - p[:, np.newaxis]) #pmx = price matrix, upper triangular

m = GEKKO(remote=False)
q = m.Array(m.Var,(4,4),lb=0,ub=10)
# only upper triangular can change
for i in range(4):
    for j in range(4):
        if j<=i:
            q[i,j].upper=0 # set upper bound = 0

def profit(q):
    profit = np.sum(q.flatten() * pmx.flatten())
    return profit

for i in range(4):
    m.Equation(np.sum(q[i,:])<=10)
    m.Equation(np.sum(q[:,i])<=8)
m.Maximize(profit(q))

m.solve()

print(q)

This gives the solution:

[[[0.0] [2.5432017412] [3.7228765674] [3.7339217013]]
 [[0.0] [0.0] [4.2771234426] [4.2660783187]]
 [[0.0] [0.0] [0.0] [0.0]]
 [[0.0] [0.0] [0.0] [0.0]]]
John Hedengren
  • 12,068
  • 1
  • 21
  • 25
  • 1
    Author created a new Gekko question that references this post: https://stackoverflow.com/questions/61646792/gekko-optimization-in-matrix-form I'll just leave the answer in both places with a little more explanation on the gekko question. – John Hedengren May 06 '20 at 23:26
  • 1
    I think you just saved my sanity. I'm going to test it and let you know. THANKS! – Chet May 06 '20 at 23:29