基于GEKKO的MHE与MPC库存模型：实测数据替代均值计算需求

Let's break down and fix the two key issues in your inventory model: integrating actual sales data directly and correctly tracking backlog when sales exceed available inventory.

Key Issues in Current Model #

1. Inventory Update Logic: Your current equation m.Inv.dt() == m.Order_delay - m.Inv_set uses a setpoint (Inv_set) instead of actual sales data, so it doesn't respond to real-time sales fluctuations.
2. Backlog Calculation: There's no explicit link between unmet sales and backlog growth—your current equations don't account for when sales exceed available inventory.

Step-by-Step Fixes #

1. Pass Actual Sales Data as a GEKKO Parameter

Instead of treating sales as a state variable (SV) with measurements, define it as a time-varying parameter (Param) in both MHE and MPC models. This lets you directly use real sales data in your dynamic equations.

2. Correct Inventory & Backlog Dynamics

Replace your existing inventory and backlog equations with logic that:

• Calculates how much of the sales demand can be fulfilled with current inventory + incoming shipments
• Updates inventory based on fulfilled sales
• Accumulates backlog for unmet demand

Here's the revised equation core:

# Define actual sales as a parameter
m.Sales = m.Param(value=np.zeros(len(m.time)))

# Calculate fulfilled sales (min of available stock + incoming, or sales demand)
m.fulfilled = m.min2(m.Inv + m.Order_delay, m.Sales)

# Inventory dynamics: incoming shipments minus fulfilled sales
m.Equation(m.Inv.dt() == m.Order_delay - m.fulfilled)

# Backlog dynamics: unfulfilled sales add to backlog
m.Equation(m.Backlog.dt() == m.Sales - m.fulfilled)

3. Adjust MHE/MPC Configuration

Update the models to use the new Sales parameter instead of Inv_set, and tweak the objective function to prioritize backlog reduction (since that's a key concern).

Modified Full Code #

Here's the updated code with all fixes applied, with comments marking changes:

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

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[20:30] = 30
SalesData[30:35] = 30
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.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: Add actual sales as a parameter
    m.Sales = m.Param(value=np.zeros(len(m.time)))

############################# Revised Equations ##############################
    # Calculate fulfilled sales (min of available stock + incoming, or sales demand)
    m.fulfilled = m.min2(m.Inv + m.Order_delay, m.Sales)
    # Inventory dynamics: incoming shipments minus fulfilled sales
    m.Equation(m.Inv.dt() == m.Order_delay - m.fulfilled)
    # Backlog dynamics: unfulfilled sales add to backlog
    m.Equation(m.Backlog.dt() == m.Sales - m.fulfilled)
    # Optional: Link storage location to safety stock target
    m.Equation(m.Storage_Location == m.Inv - safetystock)
############################# Revised Equations ##############################

## Objective
    m.Equation(m.e_Inv == m.Inv - safetystock) # Track safety stock target
    m.Equation(m.e_backlog == m.Backlog) # Prioritize backlog reduction
    m.Obj(100*m.e_backlog**2) # High weight on backlog to minimize it
    m.Obj(m.e_Inv**2) # Maintain safety stock

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

# MV tuning
mhe.Order.STATUS = 0
mhe.Order_delay.STATUS = 0

# CV tuning
mhe.e.STATUS = 0
mhe.e.FSTATUS = 0
mhe.e_backlog.STATUS = 0
mhe.e_backlog.FSTATUS = 1 # Use measured backlog
mhe.e_Storage.STATUS = 0
mhe.e_Storage.FSTATUS = 0
mhe.e_Inv.STATUS = 0
mhe.e_Inv.FSTATUS = 1 # Use measured inventory

# Solve settings
mhe.options.IMODE = 5 # MHE mode
mhe.options.NODES = 2
mhe.options.CV_TYPE = 3

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

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

# MV tuning
mpc.Order.STATUS = 1 # Allow optimizer to adjust orders
mpc.Order.DCOST = 0.1 # Smooth order changes
mpc.Order_delay.STATUS = 0 # Don't adjust delayed orders (fixed by transport time)

# CV tuning
mpc.e.STATUS = 0
mpc.e.FSTATUS = 0
mpc.e_backlog.STATUS = 1 # Add backlog to objective
mpc.e_backlog.COST = 100 # High priority to reduce backlog
mpc.e_Storage.STATUS = 0
mpc.e_Inv.STATUS = 1 # Track safety stock
mpc.e_Inv.FSTATUS = 0

db_inv = 2
mpc.e_Inv.SPHI = db_inv # Target range around safety stock
mpc.e_Inv.SPLO = -db_inv
mpc.e_Inv.TR_INIT = 1
mpc.e_Inv.TAU = 1

# Solve settings
mpc.options.IMODE = 6 # MPC mode
mpc.options.NODES = 2
mpc.options.CV_TYPE = 3

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

print(1)
for i in range(1,Loops):
    print(i)
    #################################
    ### Moving Horizon Estimation ###
    #################################
    # Update MHE sales parameter with historical data
    mhe.Sales.value = SalesData[i-ntm:i+1] if i>=ntm else SalesData[0:i+1]
    mhe.Order.MEAS = Order[i-1]
    mhe.Order_delay.MEAS = Order_delay[i-1]
    mhe.Inv.MEAS = Inv[i-1]
    mhe.Backlog.MEAS = Backlog[i-1]
    mhe.solve(disp=False)
    # Retrieve MHE results if successful
    if (mhe.options.APPSTATUS==1):
        Setpoint_mhe[i] = mhe.Sales.value[-1]
        Order_mhe[i] = mhe.Order.NEWVAL
        print('Mhe ', i)

    #################################
    ### Model Predictive Control ###
    #################################
    # Update MPC sales parameter with current + future data (simplified here)
    mpc.Sales.value = SalesData[i:i+MPC_time+1]
    # Update current inventory/backlog measurements
    mpc.Inv.MEAS = Inv[i-1]
    mpc.Backlog.MEAS = Backlog[i-1]
    mpc.solve(disp=False,GUI=False)

    # Capture output values
    Inv[i] = mpc.Inv.MODEL
    Order[i] = mpc.Order.NEWVAL
    Order_delay[i] = mpc.Order_delay.NEWVAL
    Error[i] = mpc.Inv.MODEL - safetystock
    Setpoint[i] = SalesData[i]
    Error_backlog[i] = mpc.Backlog.MODEL
    Backlog[i] = mpc.Backlog.MODEL
    Store_inventory[i] = mpc.Storage_Location.MODEL
    print('Mpc ', i)

    with open(mpc.path+'//results.json') as f:
        results = json.load(f)

    # Print performance metrics
    Total_sales = sum(SalesData[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 Shipped = %.2f" % Total_send)
    print("Current Inventory = %.2f" % Inv[i])
    print("Current Backlog = %.2f" % Backlog[i])

    # Update plot
    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='Actual Sales')
    plt.plot(tm[0:i+1],Inv[0:i+1],'g.-',label=r'Current Inventory',linewidth=2)
    plt.plot(tm[i]+mpc.time,mpc.Inv.value,'g--',label=r'Predicted Inventory',linewidth=2)
    plt.plot([tm[i],tm[i]],[-50,100],'k-')
    plt.axhline(y=safetystock, color='b', linestyle='--', label='Safety Stock')
    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(tm[0:i+1],Order[0:i+1],'--',label=r'Order History',linewidth=1)
    plt.plot(tm[i]+mpc.time,mpc.Order.value,'b--',label=r'Order Plan',linewidth=1)
    plt.plot(tm[0:i+1],Order_delay[0:i+1],'b.-',label=r'Shipments Received',linewidth=2)
    plt.plot(tm[i]+mpc.time,mpc.Order_delay.value,'b-',label=r'Predicted Shipments',linewidth=3)
    plt.ylabel('Replenishment')
    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='Excess Inventory (Above Safety)')
    plt.plot(tm[0:i],Backlog[0:i],'r--',MarkerSize=1,label='Backlog')
    plt.ylabel('Retail Store Metrics')
    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')

Key Changes Explained #

• Direct Sales Integration: We added m.Sales as a parameter to both models, updated with real SalesData in each loop, so the model reacts directly to actual sales fluctuations.
• Fulfillment Logic: The m.min2 function ensures we only deduct fulfilled sales from inventory, preventing negative stock values.
• Backlog Dynamics: Backlog now correctly accumulates when sales exceed available inventory, fixing your initial issue where backlog didn't rise as expected.
• Objective Prioritization: We increased the weight on backlog in the objective function to prioritize reducing unmet demand, while keeping inventory aligned with safety stock targets.
• MHE/MPC Tweaks: We added inventory and backlog measurements to MHE for better estimation accuracy, and locked the delayed order in MPC since it's fixed by transport lead time.

This should resolve both of your core issues: the model now uses actual sales data for calculations, and backlog correctly tracks unmet demand when sales exceed available inventory.

内容的提问来源于stack exchange，提问作者Duco Schaefer