5

我一直在研究这个模拟商店库存的模型。唯一我不能正确的是在计算库存时使用测量数据。

目前它只使用平均测量数据进行计算。有没有一种方法可以直接使用测量数据进行计算？我一直在研究为 GEKKO 示例 Lab H 提供的示例。

我目前出错的是模拟“销售数据”的第一步高于库存。这应该会导致 Backlog 增加。但目前它并没有这样做。库存仅对平均测量数据（Inv_set）做出反应

先感谢您，

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')
4

1 回答 1

0

以下是您当前应用程序的一些问题：

  • 它使用默认的 IPOPT 求解器，它将为您提供分数库存，而不是整数解决方案。您需要使用m.options.SOLVER=1混合整数求解器。
  • 您需要让求解器有机会更改backlogMHE 中的其他变量或其他变量。它们应该被定义为m.MV()带有STATUS=1.
  • MHE 和 MPC 应用程序都在尝试使用 和 最小m.Obj(m.e_backlog)化m.Obj(m.e_Storage)。这个对吗？您可能需要将它们分别定义为mhe.Obj()或mpc.Obj()。

我建议您从一个 MHE 应用程序开始，并通过估计未测量的干扰来查看它是否可以将您的模型与测量值对齐。然后可以将这些未测量的干扰传输到您的 MPC 应用程序，以提供用于前向预测和优化的更新模型。如果您想更新 MPC 中的库存等内容，则可以使用FSTATUS=1更新内部偏差并与测量结果保持一致。机器学习和动态优化课程中还有关于 MHE 和 MPC 的附加教程。

于 2020-01-23T04:57:03.163 回答