控制horizon短于预测horizon时目标函数异常与求解器无界问题咨询

问题背景 #

我正在构建一个模拟“提纯-存储-使用”工艺流程的调度测试问题，优先级设定如下：

• 首要目标：最大化消耗存储中的纯产品
• 次要目标：若纯产品使用达到约束，则存储产品
• 最后目标：若存储达到上限（由存储高限定义），则使用未提纯的原料

我观察到在部分场景（特定初始条件和目标函数系数）下，解会出现最后一步操纵变量（MV）变动时，未考虑预测horizon外下一步的存储限制违反情况。我知道在积分过程中固定自变量最后两个点，并不意味着预测horizon外的因变量预测值会在限制内，除非最后一步积分自变量的导数设为0。但连接m.mv_Pure_Use的最后两个节点时，仍遇到了问题。我通过以下语句权衡库存累积与生产：

1. 第88行：m.mv_Pure_Use.COST=-10
2. 第95行：m.Minimize(store_acc_cost*m.Store_acc)

具体疑问 #

1. 未添加第65-66行的Connections时，目标函数最终值为-8117；添加后变为-1.E+15，但解仍合理。我预期无论是否添加连接，目标函数值应处于同一数量级。
2. 将store_acc_cost从-8改为-7时，求解器报告Unbounded solution（无界解）；改为-6时又能得到解。

对应的GEKKO代码 #

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 11 09:40:06 2023

@author: strydohj
"""

from gekko import GEKKO
import numpy as np
import json
import pandas as pd

m=GEKKO(remote=False)
duration=10   #schedule duration
m.time=np.linspace(0,9,10)  #10 day schedule

#set up vectors to instantiate Parameters with.
rundown_schedule=np.full(duration,100)  #init an upstream rundown schedule, tons/h, but schedule in days
raw_process_max  =np.full(duration,100)  #init raw process max sched, tons/h but schedule in days
pure_use_max  =np.full(duration,150)  #init  purification max sched, tons/h,but schedule in days 

pure_use_max[3]=89

#model Constants
m.c_Yield=m.Const(value=0.90,name='yield')

#define model Parameters
m.p_Raw_Rundown        =m.Param(value=rundown_schedule,name='rundown_schedule')
m.p_Raw_Process_Max    =m.Param(value=raw_process_max,name='raw_process_max')
m.p_Pure_Use_Max       =m.Param(value=pure_use_max,name='Purification Maximum')

#define MV and CV's
m.mv_Raw_Flare             =m.MV(value=2,lb=0,ub=100,name='raw_flare')     #tons/h
m.mv_Raw_Purify            =m.MV(value=90,lb=0, ub=100,name='raw_purify')   #tons/h
m.mv_Raw_Use               =m.MV(value=8,lb=0, ub=150,name='raw_use')      #tons/h
m.cv_Store                 =m.CV(value=100,name='store')         #tons
m.mv_Pure_Use              =m.MV(value=81,lb=0,ub=110,name='pure_use')      #tons/h

m.Store_acc                =m.SV(0,name='Store_acc')

#Intermediate calculations
m.Pure_Product=m.Intermediate(m.mv_Raw_Purify*m.c_Yield,name='Pure_Product')

#Setup Equations
m.Equation(m.mv_Raw_Purify+m.mv_Raw_Use+m.mv_Raw_Flare==m.p_Raw_Rundown)
m.Equation(m.cv_Store.dt()==(m.mv_Raw_Purify*m.c_Yield-m.mv_Pure_Use)*24)
m.Equation(m.mv_Raw_Purify<=m.p_Raw_Process_Max)
m.Equation(m.mv_Pure_Use<=m.p_Pure_Use_Max)
m.Equation(m.Store_acc==m.mv_Raw_Purify*m.c_Yield-m.mv_Pure_Use)
m.fix_final(m.Store_acc,0)

#--------SETTING UP GLOBAL OPTIONS
m.options.NODES=3  #collocation nodes
m.options.SOLVER=1 # 1=APOPT, 2=BPOPT, 3=IPOPT
m.options.CV_TYPE=1  #1 = linear from dead band trajectory

m.options.REQCTRLMODE=3  #3= CONTROL
m.options.IMODE=6    #dynamic control
m.options.MAX_ITER=50
m.options.MV_DCOST_SLOPE=1  #Increase MV Movement Cost over time
m.options.CV_WGT_SLOPE=1    #Increase CV Error Penalty over time

#Create prediction horizon - keep mv fixed for laste two time points.
#m.Connection(m.mv_Pure_Use,m.mv_Pure_Use,8,9)  #connect end point nodes
#m.Connection(m.mv_Pure_Use,m.mv_Pure_Use,8,9,1,2) #connect internal nodes


mv_list=[m.mv_Pure_Use,m.mv_Raw_Purify,m.mv_Raw_Flare,m.mv_Raw_Use]

for mv in mv_list:
    mv.STATUS=1   #Use this CV
    mv.DCOST=0.01    #Movement cost

#CV Options
m.cv_Store.STATUS=1  #control this cv
m.cv_Store.TAU=1     #Time constant for trajectory
m.cv_Store.TR_INIT=0  #1=Recenter at cold start.
m.cv_Store.TR_OPEN=1    #Opening shape of trajectory
m.cv_Store.WSPLO=200    #penalty for exceeding limits
m.cv_Store.WSPHI=200    

m.cv_Store.SPLO=20
m.cv_Store.SPHI=140

#set up specific COST variables
m.mv_Pure_Use.COST=-10   #First priority, up to constraints
m.mv_Raw_Purify.COST=0
m.cv_Store.COST=0
m.mv_Raw_Use.COST=0
m.mv_Raw_Flare.COST=0

store_acc_cost=-8
m.Minimize(store_acc_cost*m.Store_acc)  #additional objective on accumulation

m.solve()

with open(m.path+'//results.json') as f:
    results = json.load(f)
result_df=pd.DataFrame(results)

result_df[['store','store.tr_hi','store.tr_lo']].plot()
result_df[['raw_flare','raw_use','raw_purify','pure_use']].plot()
print('Total Accumulated:\t',result_df['store_acc'].sum())
print('Total Used:\t',result_df['pure_use'].sum()) 
#print(m.Store)     

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