5

I have been working on this model what simulates the inventory of a store. Only thing i can not get right is using the measurement data with in the calculation of the inventory.

Currently it only uses the average measurement data for the calculation. Is there a method where i can use directly the measurement data for calculations? I have been working of the examples provided for GEKKO example Lab H.

What i currently goes wrong is that the first steps of te simulation the "sales data" is higher than the inventory. This should result in Backlog to rise. But currently it does not do that. The inventory is only reacting on the average measurement data (Inv_set)

Thank you in advance, 

from random import random
from gekko import GEKKO
import matplotlib.pyplot as plt
import json


Loops = 50

# number of data points (every day)
n = Loops + 1

# model time
tm = np.linspace(0,Loops,(Loops+1))



safetystock = 5

# time MPC
t = np.linspace(0,40,41)

MPC_time = len(t)

Transport_time = 3


## Output variables

Inv     = np.ones(Loops)*0
Order   = np.ones(Loops)*0
Order_delay = np.ones(Loops)*0
Order_unfilled = np.ones(Loops)*0
Order_mhe = np.ones(Loops)*0

Setpoint        = np.ones(Loops)*0
Setpoint_mhe    = np.ones(Loops)*0
Error           = np.ones(Loops)*0
Error_backlog   = np.ones(Loops)*0

Backlog         = np.ones(Loops)*0
Store_inventory = np.ones(Loops)*0





Sales_sp = np.ones(len(t))*15
Sales_sp_full = np.ones(len(tm)+len(t))*15



## Sales data for real store
SalesData = np.ones(len(tm)+len(t))*15
#SalesData[0:Loops] = Sold_schiphol_globaltime

#SalesData[5:10] = 5
#SalesData[10:15] = 5
SalesData[20:30] = 30
SalesData[30:35] = 30
#SalesData[35:40] = 5
#SalesData[50:75] = 5
#SalesData[85:100] = 10

SalesData[0] = 0

SalesData[11] = 5






# constants
mass = 1

# Parameters

mhe = GEKKO(name='mhe',remote=True)
mpc = GEKKO(name='mpc',remote=True)



for m in [mhe,mpc]:

    Z = m.Param(value=0)

    m.tau = m.Param(value=1)
    m.Kp = m.Param(value=1)

    # Manipulated variable
    m.Test = m.MV(value=0, lb=0, ub=100)

    m.Order = m.MV(value=0, lb=0, ub=100,integer=True)
    m.Order_delay = m.MV(value=0, lb=0, ub=100,integer=True)
    Delay_Order = Transport_time            # leadtime of transport

    m.delay(m.Order,m.Order_delay,Delay_Order)

    # Controlled Variable
    m.Inv = m.SV(value=0,name='Inv1' ,integer=True)          # inventory of store
    m.Inv_set = m.SV(value=0 , name='Inv_set',integer=True)  # Demand of shirts

    m.Backlog = m.SV(value=0 , lb=0)  
    m.Inv_test = m.SV(value=0)        
    m.Order_unfilled = m.SV(value=0)
    m.Storage_Location = m.SV(value=0 ,lb=0)

    # Error
    m.e = m.CV(value=0,name='e')
    m.e_backlog = m.CV(value=0,name='e_backlog')
    m.e_Storage = m.CV(value=0)
    m.e_Inv = m.CV(value=0)




    ############################# New Equations ##############################

    m.Equation( Z ==   - m.Inv_set + m.Order_unfilled ) 

    m.Equation( Z ==  m.Inv - m.Inv_set + m.Backlog - m.Storage_Location )
    m.Equation(m.Inv.dt() ==   m.Order_delay - m.Inv_set )




    ############################# New Equations ##############################
    ## Objective
    m.Equation(m.e==m.Inv_set)

    m.Equation(m.e_Inv==m.Inv-m.Inv_set -safetystock) #
    m.Equation(m.e_backlog ==  ( m.Backlog) )
    m.Equation(m.e_Storage ==  ( m.Storage_Location) )

    m.Obj(m.e_backlog)
    m.Obj(m.e_Storage)

#################################################

# Configure MHE
# 120 second time horizon, steps of 4 sec

ntm = 20
mhe.time = np.linspace(0,ntm,(ntm+1))


# MV tuning
mhe.Order.STATUS = 0 #allow optimizer to change
mhe.Order_delay.STATUS = 0 #allow optimizer to change

# CV tuning
mhe.e.STATUS = 0 #add the CV to the objective
mhe.e.FSTATUS = 0 

mhe.e_backlog.STATUS = 0
mhe.e_backlog.FSTATUS = 0


mhe.e_Storage.STATUS = 0
mhe.e_Storage.FSTATUS = 0
#
mhe.e_Inv.STATUS = 0
mhe.e_Inv.FSTATUS = 0

#mhe.Inv_set.STATUS = 0
mhe.Inv_set.FSTATUS = 1


# Solve
mhe.options.IMODE = 5 # control
mhe.options.NODES   = 2 # Collocation nodes
mhe.options.CV_TYPE = 3 #Dead-band
#mhe.options.SOLVER = 1 #Dead-band

##############################################
##################################################################
# Configure MPC



mpc.time = np.linspace(0,MPC_time,(MPC_time+1))


# MV tuning
mpc.Order.STATUS = 1 #allow optimizer to change
mpc.Order.DCOST = 0.1 #smooth out MV

mpc.Order_delay.STATUS = 1 #allow optimizer to change
mpc.Order_delay.DCOST = 0.1 #smooth out MV

# CV tuning
mpc.e.STATUS = 1 #add the CV to the objective
mpc.e.FSTATUS = 0 


mpc.e_backlog.STATUS = 0
mpc.e_backlog.COST = 100

mpc.e_Storage.STATUS = 0

mpc.e_Inv.STATUS = 1
mpc.e_Inv.FSTATUS = 0

mpc.Inv_set.FSTATUS = 1

db = 0
db_inv = 2
mpc.e.SPHI = db #set point
mpc.e.SPLO = -db #set point
mpc.e.TR_INIT = 1 #setpoint trajectory
mpc.e.TAU = 1 #time constant of setpoint trajectory

mpc.e_Inv.SPHI = db_inv  #set point #+ safetystock
mpc.e_Inv.SPLO = -db_inv #set point + safetystock
mpc.e_Inv.TR_INIT = 1 #setpoint trajectory
mpc.e_Inv.TAU = 1 #time constant of setpoint trajectory
#


# Solve
mpc.options.IMODE = 6 # control
mpc.options.NODES   = 2 # Collocation nodes
mpc.options.CV_TYPE = 3 #Dead-band
#mpc.options.SOLVER = 1 #Dead-band

##################################


print(1)

for i in range(1,Loops):



    print(i)


    #################################
    ### Moving Horizon Estimation ###
    #################################
    mhe.Inv_set.MEAS = SalesData[i]
    mhe.Order.MEAS = Order[i-1]
    mhe.Order_delay.MEAS = Order_delay[i-1]


    mhe.solve(disp=False)
    # Parameters from MHE to MPC (if successful)
    if (mhe.options.APPSTATUS==1):

        # CVs
        Setpoint_mhe[i] = mhe.Inv_set.MODEL    
        Order_mhe[i] = mhe.Order.NEWVAL
        print('Mhe ', i)

    #################################
    ### Model Predictive Control  ###
    #################################

    mpc.Inv_set.MEAS = SalesData[i]



    db = 0
    mpc.e.SPHI = SalesData[i] + db  #set point
    mpc.e.SPLO = SalesData[i] -db #set point


    mpc.solve(disp=False,GUI=False) #

    if (mpc.options.APPSTATUS==1):
        mpc.Order.MEAS = mpc.Order.NEWVAL #update production value
        mpc.Order_delay.MEAS = mpc.Order_delay.NEWVAL #update production value

        # output:

        Inv[i] = mpc.Inv.MODEL   

        Order[i] = mpc.Order.NEWVAL
        Order_delay[i] = mpc.Order_delay.NEWVAL
#        Order_unfilled[i] = mpc.Order_unfilled.MODEL
        Error[i] =  mpc.Inv_set.MODEL
        Setpoint[i] = Sales_sp_full[i]
#        Error[i] = mpc.e.MODEL
        Error_backlog[i] = mpc.e_backlog.MODEL

        Backlog[i]   = mpc.Backlog.MODEL
        Store_inventory[i] = mpc.Storage_Location.MODEL


        print('Mpc ', i)
        with open(m.path+'//results.json') as f:
                results = json.load(f)




    Total_sales = sum(SalesData[j] for j in range(i))
    Total_sales_mhe = sum(Setpoint_mhe[j] for j in range(i))
    Total_send = sum(Order_delay[j] for j in range(i))

    print("Total Sales        = %.2f" % Total_sales)
    print("Total sold mhe     = %.2f" % Total_sales_mhe)

    print("Total send         = %.2f" % Total_send)
    print("Shirt left in store          = %.2f"  % Inv[i])



    plt.clf()
    j = max(0,i-ntm-1)
    ax=plt.subplot(3,1,1)
    ax.grid()

    ax.axvspan(tm[j], tm[i], alpha=0.2, color='purple')
    ax.axvspan(tm[i], tm[i]+mpc.time[-1], alpha=0.2, color='orange')
    plt.text(tm[i]+10,40,'Future: MPC')
    plt.text(tm[j]+1,40,'Past: MHE')
    plt.plot(tm[0:i+1],SalesData[0:i+1],'ro',MarkerSize=2,label='Sales data')
    plt.plot(mhe.time-ntm+tm[i],mhe.Inv_set.value,'k.-',linewidth=1,alpha=0.7,label=r'Demand arrived history mhe')
    plt.plot(tm[0:i+1],Error[0:i+1],'b-',linewidth=1,alpha=0.7,label=r'Demand arrived history mpc')

#    plt.plot(mhe.time-ntm+tm[i],mhe.Inv.value,'r-',linewidth=2,alpha=0.7,label=r'Inventory At retail store')

    plt.plot(tm[0:i+1],Inv[0:i+1],'g.-',label=r'Total inventory',linewidth=2)

    plt.plot(tm[i]+mpc.time,results['e.bcv'],'r-',label=r'Demand predicted',linewidth=2)             
    plt.plot(tm[i]+mpc.time,results['e.tr_hi'],'k--',label=r'Demand estimation bandwidth ')
    plt.plot(tm[i]+mpc.time,results['e.tr_lo'],'k--')          
    plt.plot([tm[i],tm[i]],[-50,100],'k-')
    plt.ylabel('Retail store inventory')
    plt.legend(loc=3)
    plt.xlim(0, tm[i]+mpc.time[-1])
    plt.ylim(-25, 50)

    ax=plt.subplot(3,1,2)
    ax.grid()
    ax.axvspan(tm[j], tm[i], alpha=0.2, color='purple')
    ax.axvspan(tm[i], tm[i]+mpc.time[-1], alpha=0.2, color='orange')
    plt.text(tm[i]-2,10,'Current Time',rotation=90)
    plt.plot([tm[i],tm[i]],[-20,70],'k-',label='Current Time',linewidth=1)

    plt.plot(mhe.time-ntm+tm[i],mhe.Inv_set.value,'k-',linewidth=1,alpha=0.7,label=r'Demand arrived history')

    plt.plot(tm[0:i+1],Order[0:i+1],'--',label=r'Ordering shirts history',linewidth=1)
    plt.plot(tm[i]+mpc.time,mpc.Order.value,'b--',label=r'Order shirts plan',linewidth=1)
    plt.plot(tm[i]+mpc.time[0],mpc.Order_delay.value[1],color='blue', marker='.',markersize=15)

    plt.plot(tm[0:i+1],Order_delay[0:i+1],'b.-',label=r'Shirt arrived history',linewidth=2)
    plt.plot(tm[i]+mpc.time,mpc.Order_delay.value,'b-',label=r'Shirt arrived plan',linewidth=3)


    plt.ylabel('Replenishment of the store')
    plt.legend(loc=3)    
    plt.xlim(0, tm[i]+mpc.time[-1])
    plt.ylim(-10, 50)


    ax=plt.subplot(3,1,3)
    ax.grid()
    ax.axvspan(tm[j], tm[i], alpha=0.2, color='purple')
    ax.axvspan(tm[i], tm[i]+mpc.time[-1], alpha=0.2, color='orange')
    plt.plot([tm[i],tm[i]],[-20,80],'k-',label='Current Time',linewidth=1)

    plt.plot(tm[0:i],Store_inventory[0:i],'b--',MarkerSize=1,label='Inventory at the store ')
    plt.plot(tm[0:i],Backlog[0:i],'r--',MarkerSize=1,label='Backlog of the store')

    plt.ylabel('Retail store')
    plt.xlabel('Time (Days)')
    plt.legend(loc=3)
    plt.xlim(0, tm[i]+mpc.time[-1])
    plt.ylim(-1, 40)



    plt.draw()
    plt.pause(0.09)

plt.savefig('tclab_mhe_mpc.png')
Duco Schaefer
  • 51
  • 1
  • I like the dynamic plot that you've modified for your project. It looks like it is off by one index value at the current time. Please post an updated solution when you resolve the problems with your application. MHE and MPC for inventory management is a nice application. There are probably some applications in refining and chemicals as well that could be inspired by your application in retail for supply chain management. – TexasEngineer Jan 25 '20 at 16:40

1 Answers1

0

Here are a few issues with your current application:

  • It is using the default IPOPT solver that will give you fractional inventory, not integer solutions. You need to use m.options.SOLVER=1 for a Mixed Integer solver.
  • You need to give the solver the opportunity to change backlog or other variables in the MHE. They should be defined as m.MV() types with STATUS=1.
  • The MHE and MPC applications are both trying to minimize with m.Obj(m.e_backlog) and m.Obj(m.e_Storage). Is this correct? You may need to define these separately as mhe.Obj() or mpc.Obj().

I recommend that you start with just an MHE application and see if it can align your model to the measurements by estimating you unmeasured disturbances. These unmeasured disturbances can then be transferred to your MPC application to give an updated model for forward prediction and optimization. If you want to update something like inventory in your MPC then you can use FSTATUS=1 to update the internal bias and align with the measurement. There are additional tutorials on MHE and MPC in the Machine Learning and Dynamic Optimization course.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25