我一直在研究这个模拟商店库存的模型。我唯一不能做对的是使用测量数据来计算库存。
目前，它仅使用平均测量数据进行计算。 有没有可以直接使用测量数据进行计算的方法？ 我一直在研究为GEKKO示例实验室H提供的示例。
我目前出错的是，模拟的第一步是"销售数据"高于库存。这应该会导致积压增加。但目前它没有这样做。库存仅对平均测量数据做出反应 (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')
以下是当前应用程序的一些问题：
- 它使用默认的 IPOPT 求解器，该求解器将为您提供分数库存，而不是整数解。您需要将
m.options.SOLVER=1用于混合整数求解器。 - 您需要让求解器有机会更改 MHE 中的
backlog或其他变量。它们应定义为具有STATUS=1的m.MV()类型。 - MHE和MPC应用程序都试图通过
m.Obj(m.e_backlog)和m.Obj(m.e_Storage)最小化。这是对的吗？您可能需要将它们分别定义为mhe.Obj()或mpc.Obj()。
我建议您仅从MHE应用程序开始，看看它是否可以通过估计未测量的干扰来使您的模型与测量值保持一致。然后，可以将这些未测量的干扰传输到您的 MPC 应用程序，以提供用于前向预测和优化的更新模型。如果要更新MPC中的库存之类的内容，则可以使用FSTATUS=1来更新内部偏差并与测量值保持一致。机器学习和动态优化课程中有关于 MHE 和 MPC 的其他教程。