Follow up question: GEKKO optimization in matrix form

Question

this is a follow up question to the one I posted earlier:

GEKKO - optimization in matrix form

I need to add one more constraint that tracks "inventory" ("Inv"), which tracks the sum(q[i,:] - q[:,i]). "Inv" will be a 4X1 column vector. I tried the following:

 m = GEKKO(remote=False)
    q = m.Array(m.Var,(4,4),lb=0,ub=10)

    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

    Inv[0]=0
    for i in range(4):
        m.Equation(np.sum(q[i,:])<=10)
        m.Equation(np.sum(q[:,i])<=8)
        m.Equation(I[i] = I[i-1] + (np.sum(q[i,:]) - np.sum(q[:,i]))) # New Line 1a inventory 
        Inv[i] = Inv[i-1] + (np.sum(q[i,:]) - np.sum(q[:,i])) # New Line 1b inventory. Keep either 1a or 1b 
        m.Equation(Inv[i] <= 15) # New Line 2 inventory constraint
        m.Equation(Inv[4] = 0) # New Line 3 ending inventory should be equal to starting inventory

    m.Maximize(profit(q))

    m.solve()
    print(q)
qr = np.array([[q[i,j].value[0] for j in range(4)] for i in range(4)])
Ir = np.array([Inv[i].value[0] for i in range(4)]) #New Line 4

Errors :

1a. Adding New Line 1a: "keyword can't be an expression"

1b. Replacing New Line 1a with 1b: no issues (but, I'm not sure if GEKKO will keep track of I or not.Also, I need to define "I", the way "q" was done...not sure how). Replacing = comment out 1a, then run the code with 1b.

3. New Line 2: Error = "object of type 'int' has no len()"; but type(I) shows as ndarray. (Kept New Lines 1a and 1b, and then added New Line 2)

4. New Line 3:Error = "keyword can't be an expression" (Kept Line 32and then added Line 3)

5. New Line 4: Error "'numpy.ndarray' object has no attribute 'value'" [removed Lines 3 and 4. This makes sense as if I can't capture "Inv" in the model, then it won't have the value attribute)

Questions: 1. Am I defining inventory correctly?

2. If yes, can it be done in the current model specification, or will it need an entirely different formulation? If so, could you please guide on what that would be?

3. Of the various videos posted on the GEKKO website, is there a specific one I should look at for more information? I was thinking the DO video, but I don't believe this is quite a dynamic optimization problem (as I'm not trying to optimize the best path)....

Thanks again for all your help,

----UPDATE 5/10 Also tried:

Inv = m.SV()
for i in range(4):
    m.Equation(np.sum(q[i,:])<=10)
    m.Equation(np.sum(q[:,i])<=8)
    #m.Equation(I[i] = I[i-1] + (np.sum(q[i,:]) - np.sum(q[:,i])))
    m.Equation(m.Inv.dt() == m.Inv + (np.sum(q[i,:]) - np.sum(q[:,i])))
    #I[i] = I[i-1] + (np.sum(q[i,:]) - np.sum(q[:,i])
    m.Equation(m.Inv <= 15)
    #m.Equation(I[4] = 0)
m.Maximize(profit(q))

New Error: 'GEKKO' object has no attribute 'Inv'

Answer 1

One way to do this is to start with zero inventory with Inv[0]=0 and then track the inventory amount with gekko variables Inv[1:4]. A couple tips on building the model:

• Use double equal signs for equality constraints
• You can define a temporary variable such as I = Inv[2]+Inv[3] but it won't be a gekko variable
• You may also want to look at Intermediate variables for those that are explicitly calculated. This can speed up the calculation.
• I recommend this tutorial in the Dynamic Optimization course

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
Inv = m.Array(m.Var,5,lb=0,ub=15)
Inv[0].upper = 0 # start with 0 inventory
for i in range(4):
    m.Equation(np.sum(q[i,:])<=10)
    m.Equation(np.sum(q[:,i])<=8)
    # track inventory    
    m.Equation(Inv[i+1]==Inv[i] + (m.sum(q[i,:])-m.sum(q[:,i]))) 
m.Equation(Inv[4] == 0) # Use double == sign, not needed in loop
m.Maximize(profit(q))
m.solve()
print(q)
# convert to matrix form
qr = np.array([[q[i,j].value[0] for j in range(4)] for i in range(4)])
for i in range(4):
    rs = qr[i,:].sum()
    print('Row sum ' + str(i) + ' = ' + str(rs))
    cs = qr[:,i].sum()
    print('Col sum ' + str(i) + ' = ' + str(cs))    
Ir = np.array([Inv[i].value[0] for i in range(4)])
print(Ir)