使用GEKKO移动视界估计和模型预测控制，对测量数据而非平均测量做出反应 - python

我一直在研究模拟商店库存的模型。我唯一不能理解的是在计算清单时使用测量数据。

当前，它仅使用平均测量数据进行计算。
有没有一种方法可以直接将测量数据用于计算？
我一直在为GEKKO示例Lab 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或其他变量。应将它们定义为带有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更新内部偏差并与测量值保持一致。 Machine Learning and Dynamic Optimization course中还有有关MHE和MPC的其他教程。

在返回'Response'(Python)中传递多个参数 - python

我在Angular工作，正在使用Http请求和响应。是否可以在“响应”中发送多个参数。角度文件：this.http.get("api/agent/applicationaware").subscribe((data:any)... python文件：def get(request): ... return Response(seriali…

R'relaimpo'软件包的Python端口 - python

我需要计算Lindeman-Merenda-Gold（LMG）分数，以进行回归分析。我发现R语言的relaimpo包下有该文件。不幸的是，我对R没有任何经验。我检查了互联网，但找不到。这个程序包有python端口吗？如果不存在，是否可以通过python使用该包？ python参考方案 最近，我遇到了pingouin库。

如何用'-'解析字符串到节点js本地脚本？ - python

我正在使用本地节点js脚本来处理字符串。我陷入了将'-'字符串解析为本地节点js脚本的问题。render.js：#! /usr/bin/env -S node -r esm let argv = require('yargs') .usage('$0 [string]') .argv; console.log(argv…

TypeError:'str'对象不支持项目分配，带有json文件的python - python

以下是我的代码import json with open('johns.json', 'r') as q: l = q.read() data = json.loads(l) data['john'] = '{}' data['john']['use…

Python ThreadPoolExecutor抑制异常 - python

from concurrent.futures import ThreadPoolExecutor, wait, ALL_COMPLETED def div_zero(x): print('In div_zero') return x / 0 with ThreadPoolExecutor(max_workers=4) as execut…