私は、店舗の在庫をシミュレートするこのモデルに取り組んできました。私が正しくできないのは、在庫の計算に測定データを使用することだけです。
現在、計算には平均測定データのみが使用されます。測定データを直接計算に使用する方法はありますか? GEKKO example 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')