Inventory Management Optimization problem

[Mini_Miudo]

Hey everyone. So, I don't know if this is the right place to post this, I'm starting my dissertation so this isn't really "homework" per se, but my first task for this was to solve a "simple" (not so simple for me it seems) Inventory Management problem, where I have to optimize the points at which I should order and how much I should order. However, I never used python before (thus the "simple" problem to start with), having only started 2 days ago, so I'm struggling. I'm trying to use the "minimize" function from SciPy. I'll try to give you a detailed explanation of the problem and then I'll post the code I have so far (which has no errors but isn't working as intended) if someone wants to try and help out, I've been trying to solve this pretty much non-stop during these 2 days to no avail.

The problem

Parameters (set by the user/me):
T - time span for the problem
d - demand during each day
RT - review time, which is the time between 2 placed orders
LT - lead time, which is the time that goes from when the order is placed until it arrives
h - on-hand inventory cost
I0 - initial stock
s - target stock

The variables to calculate: (NOTE: I used i as the counter in my code, but using j here so everything doesn't become italic)
Qt - quantity to order each day - main variable as it will influence all the others
It - on-hand inventory at the end of each day - I defined "today's" inventory as It[j] = It[j-1] - d[j] + Qt[j-LT] -> j-LT because the quantity I ordered arrives LT days later. It[0] = I0 - d[0]
NIt - net inventory at the start of each day - I defined it as NIt[j] = NIt[j-1] - d[j] + Qt[j-1] with NIt[0] = I0
I+ - on-hand inventory (It > 0), defined as I+ = It
I- - backorder (It < 0), defined as I- = -It

Constraints:
Qt >= 0
I+ >= 0
I- >= 0

Objective function:
minimize cost = h*I+ + b*I-

The code

Taking the above into account, I tried to create a program that optimizes the objective function by changing the value of Qt. However, it isn't working (Qt is always equal to the initial value I give it - if I say the initial value is 1 for every day, then the minimize function returns Qt = 1 for every day as well, same for every number) and it was brought to my attention that the while loop I'm doing in the "Objective" function I created doesn't make sense. So, how do I make the minimize function actually optimize the values for Qt, while also making sure these values for Qt affect It/NIt/I+/I- and, of course, the objective function? The code is in the spoiler:

Hide/Show

import random
import numpy as np
from scipy.optimize import minimize

#Parameters
T = 30 #time span (days)
d = [] #demand (random)
for j in range(T):
    d.append(random.randint(0, 2)) 
RT = 1  # review time (days) - time between 2 placed orders
LT = 1  # lead time (days) - time that goes from when an order is placed to the moment it arrives
h = 1   # on-hand inventory cost (unit cost)
b = 2   # backorder cost (unit cost)
I0 = 1  # initial stock (units)
s = 2   # target stock (units)

#Variables
Qt = [0] * T    # quantity to order each day (units)
It = [0] * T    # on-hand inventory at the end of each day (units)
NIt = [0] * T   # net inventory at the start of each day (units) = (on-hand + in transit) inventory
Ip = [0] * T    # on-hand inventory (units) (It > 0)
Im = [0] * T    # backorder (units) (It < 0)

#Initial Values
It[0] = I0 - d[0]
NIt[0] = I0
Ip[0] = It[0]

#Counter
i=1

def objective (Qt):
    global i
    sumIp=0
    sumIm=0
    while i < T:
        if (i - LT) >= 0:
            It[i] = It[i-1] - d[i] + Qt[i-LT]
        else:
            It[i] = It[i-1] - d[i]
        NIt[i] = NIt[i-1] - d[i-1] + Qt[i-1]
        if It[i] >= 0:
            sumIp = sumIp+It[i]
            sumIm = sumIm+0
        elif It[i] < 0:
            sumIp = sumIp+0
            sumIm = sumIm-It[i]
        i += 1
    return h*sumIp+b*sumIm

    
def constraint1 (Qt):
    return Qt[0]

x0 = np.zeros(T)
con1 = {'type': 'ineq', 'fun': constraint1}
cons = [con1]
sol = minimize(objective, x0, constraints=cons)
print(sol)

Apologies if this is too much or if I said some nonsense and thanks in advance.

[Mini_Miudo]

Been working on this all day and it seems I've been able to come up with a better code, hopefully easier to help now. I still can't get my head around the Objective function. It isn't done properly whatsoever so the result of it is totally wrong but at least the rest seems to be calculated properly, based on the (wrong) value the objective function gives me.

Still based on the description of the first post, any idea how I should do my Objective function? Here's the new code (can't seem to edit the opening post):

import random
import numpy as np
from scipy.optimize import minimize

#Parameters
T = 30  #time span (days)
d = []  #demand (random, units per day)
for j in range(T):
    d.append(random.randint(0, 3)) 
RT = 1  # review time (days) - time between 2 placed orders
LT = 1  # lead time (days) - time that goes from when an order is placed to the moment it arrives
h = 1   # on-hand inventory cost (unit cost)
b = 5   # backorder cost (unit cost)
I0 = 1  # initial stock (units)
s = []  # target stock (units)
for j in range(T):
    s.append((RT+LT)*d[j]) 

#Variable Initialization
Qt = [0] * T    # quantity to order each day (units)
It = [0] * T    # on-hand inventory at the end of each day (units)
NIt = [0] * T   # net inventory at the start of each day (units) = (on-hand + in transit) inventory
Ip = [0] * T    # on-hand inventory (units) (It > 0)
Im = [0] * T    # backorder (units) (It < 0)

#Initial Values
It[0] = I0 - d[0]
NIt[0] = I0 
Qt[0] = 0

if (It[0] >= 0):
    sumIp = It[0]
    sumIm = 0

elif (It[0] < 0):
    sumIp = 0
    sumIm = -It[0]

#Counter
i=1

def objective (x):
    global sumIp, sumIm
    if (It[i-1] - d[i] + Qt[i-LT]) >= 0:
        sumIp = sumIp + x[2]
        sumIm = sumIm + 0
    else:
        sumIp = sumIp + 0
        sumIm = sumIm - x[2]
    return h*sumIp*b*sumIm

    
def constraint1 (x):
    return x[0]

def constraint2 (x):
    if (NIt[i-1] < s[i-1]):
        return x[1] - NIt[i-1] + d[i-1] - x[0]
    else:
        return x[1] - NIt[i-1] + d[i-1]

def constraint3 (x):
    return x[2] - x[1] + d[i] + x[0] - Qt[i-LT]

def test(T):
    global i
    while i < T:
        x0 = [1, NIt[0], It[0]]
        con1 = {'type': 'ineq', 'fun': constraint1}
        con2 = {'type': 'eq', 'fun': constraint2}
        con3 = {'type': 'eq', 'fun': constraint3}
        cons = [con1, con2, con3]
        sol = minimize(objective, x0, constraints=cons)
        Qt[i-1] = round(sol.x[0], 0)
        NIt[i] = round(sol.x[1], 0)
        It[i] = round(sol.x[2], 0)
        i += 1
    return Qt, NIt, It

[Qt, NIt, It] = test(T)
print ('It = [%s]' % ', '.join(map(str, It)))
print ('NIt = [%s]' % ', '.join(map(str, NIt)))
print ('Qt = [%s]' % ', '.join(map(str, Qt)))
print ('d = [%s]' % ', '.join(map(str, d)))
print ('s = [%s]' % ', '.join(map(str, s)))
print(sumIp)
print(sumIm)