是否可以为GEKKO地平线不同区域的简历指定不同的SPHI或SPLO限制？

Question

我想为系统的动态调度优化建立一个GEKKO模型。我正在探索一个玩具问题的GEKKO功能(下面的代码)。我预见到，它将被要求为某些CV为地平线的不同部分指定不同的目标，如下所示。

​

我尝试为级别SPLO的CV参数分配一个数组，但它只是将SPHI和SPLO折叠为CV的起始值。

我喜欢使用目标函数来驱动解决方案的灵活性，而不是“硬”约束。这能在一个非迭代的实现中实现吗?如果是的话，如何实现？

from gekko import GEKKO
import numpy as np
import json
import pandas as pd
from matplotlib import pyplot as plt


    
def G1_offline(timespace=100):
    
    tk_lowlimit=[37]*100       #init low limit
    tk_lowlimit[40:70]=[38]*30 #increase low limit for portion of horizon  
    
    m=GEKKO(remote=False)
    #tk_lowlimit_hard=m.Param(tk_lowlimit)  
    
    rundown_schedule=[100]*timespace  #init rundown schedule
    rundown_schedule[40:45]=[95]*5   #adjust schedile for few points
    
    
    m.time=np.linspace(0,timespace-1,timespace)
    
    m.Unit1_Feed=m.MV(value=25,lb=0,ub=60,name='Unit1 Feed')
    m.Unit2_Feed=m.MV(value=27,lb=0,ub=60,name='Unit2 Feed')
    m.Fuel=m.MV(value=10,lb=0,ub=100,name='Fuel')
    m.Rundown=m.MV(name='Rundown')  #This is a DV
    m.Efficiency=m.FV(value=0.99,lb=0.95,ub=1,name='Efficiency')                 
    m.Rundown.value=rundown_schedule
    
        
    m.Flare=m.SV(value=30,lb=0,ub=100,name='Flare')
    m.TankLevel=m.CV(value=25, lb=0,ub=300,name='tklevel')
    
    
    m.Consumers=m.MV(value=30,lb=0,ub=130,name='Consumers')
    m.Product=m.Intermediate((m.Unit1_Feed+m.Unit2_Feed)*m.Efficiency,name='Product')
    m.Balance=m.Intermediate(m.Product-m.Consumers,name='Balance')
    
    m.Equation(m.TankLevel.dt()==m.Balance)
    
    m.Equation(m.Flare==m.Rundown-(m.Unit1_Feed+m.Unit2_Feed+m.Fuel))
    #m.Equation(m.Flare>=1)
     
    #GLOBAL OPTIONS
    m.options.IMODE=6  #control mode,dynamic control, simultaneous
    m.options.NODES=2  #collocation nodes
    m.options.SOLVER=1 # 1=APOPT, 2=BPOPT, 3=IPOPT
    m.options.CV_TYPE=1  #2 = squared error from reference trajectory
    m.options.CTRL_UNITS=3   #control time steps units (3= HOURS)
    m.options.CTRL_TIME=1    #1=1 hour per time step
    m.options.REQCTRLMODE=3  #3= CONTROL
    #m.options.SCALING=2
    m.options.RTOL=1e-6
    m.options.OTOL=1e-6
    #m.options.CV_WGT_START=5
    m.options.CSV_WRITE=2
   
    
    #MV/DV modes
    m.Unit1_Feed.STATUS=1   #1 = can change
    m.Unit2_Feed.STATUS=1   #1 = can change
    m.Fuel.STATUS=1         #1 = can change
    m.Consumers.STATUS=1    #1 = can change
    m.Rundown.STATUS=0      #0  = cannot change, this is a DV
    m.Efficiency.STATUS=0
    m.Efficiency.FSTATUS=1
    
    #CV Modes
    m.TankLevel.STATUS=1    #1 = Control this CV
    #m.Flare.STATUS=0        #0 = Do Not Control this CV
   
        
    m.TankLevel.FSTATUS=1  #Allow Feedback
    m.TankLevel.STATUS=1  #Control this CV
    m.TankLevel.TAU=12     #Time constant for trajectory
    m.TankLevel.SPHI=40 #Upper limit for trajectory
    m.TankLevel.SPLO=37  #Lower limit for trajectory
    m.TankLevel.WSPLO=20   #Penalty for crossing LO limit
    m.TankLevel.WSPHI=20   #Penalty for crossing HI limit
    m.TankLevel.TR_INIT=0  #0 -Do not re-center.
    m.TankLevel.TR_OPEN=1    #Openi#ng shape of trajectory
    
    
    m.Consumers.COST=-40
    m.Unit1_Feed.COST=5
    m.Unit2_Feed.COST=4
    m.Fuel.COST=-2
    #m.Flare.COST=0
    
    m.Consumers.DCOST=15
    m.Unit1_Feed.DCOST=5
    m.Unit2_Feed.DCOST=5
    m.Fuel.DCOST=1
    
    
    m.Consumers.DMAX=10
    m.Unit1_Feed.DMAX=10
    m.Unit2_Feed.DMAX=8
    m.Fuel.DMAX=10
    
    m.Consumers.MV_STEP_HOR=1
    m.Unit1_Feed.MV_STEP_HOR=1
    m.Unit2_Feed.MV_STEP_HOR=1
    m.Fuel.MV_STEP_HOR=1
 
    m.solve(GUI=False)   
           
    with open(m.path+'//results.json') as f:
       results = json.load(f)
       #print(results)           
       results_df=pd.DataFrame(results)
    print(results_df)
    #results_df.to_excel(r'c:\data\toyproblem.xlsx')
    
    fig = plt.figure(figsize=(14,6))
    plt.plot(results_df['time'],results_df['tklevel'],color='red',label='Level')
    plt.fill_between(x=results_df['time'],y1=results_df['tklevel.tr_lo'], y2=results_df['tklevel.tr_hi'],color='green',alpha=0.2, label='Tklevel CV bounds')
    
    plt.xlabel('TIME')
    plt.title('Controlled solution')
    plt.ylabel('TankLevel')
    plt.legend(bbox_to_anchor=(0.0, 1), loc='upper left', borderaxespad=0.5)
    plt.minorticks_on()
    plt.grid(color = 'b', linestyle = '--', linewidth = 0.5, axis='y')
    plt.show()
    
    fig = plt.figure(figsize=(14,6))
    plt.plot(results_df['time'],results_df['unit1_feed'],color='red',label='Unit1')
    plt.plot(results_df['time'],results_df['unit2_feed'],color='green',label='Unit2')
    plt.plot(results_df['time'],results_df['consumers'],color='black',label='Consumers')
    plt.plot(results_df['time'],results_df['flare'],color='orange',label='Flare')
    plt.plot(results_df['time'],results_df['fuel'],color='blue',label='Fuel')
    plt.plot(results_df['time'],results_df['rundown'],color='purple',label='Rundown')
    plt.xlabel('TIME'), plt.ylabel('knm3/h'), plt.title('Independent variables'),  
    plt.legend(bbox_to_anchor=(0.0, 1), loc='upper left', borderaxespad=0.5)
    plt.minorticks_on()
    plt.grid(color = 'b', linestyle = '--', linewidth = 0.5, axis='y')
    trj_hi=results_df['tklevel.tr_hi']
    trj_lo=results_df['tklevel.tr_lo']
        
    return m,results_df
    


#----main----
c1,results_df=G1_offline(100)

Answer 1

可以定制SPHI和SPLO，而不是固定的目标值。这是通过将CV重新定义为当前值和目标值之间的差异来实现的。目标值可以是前馈traj=m.Param()，并且控制器的每个周期的值都更新为类似于traj.value = [custom_setpoint]的值。在动态优化过程中有一个这种方法的例子(参见页面底部)。

# Error
e = m.CV(value=0,name='e')
m.Equation(e==v-traj)

# CV tuning
e.STATUS = 1 #add the CV to the objective
m.options.CV_TYPE = 1 #Dead-band
db = 2
e.SPHI = db #set point
e.SPLO = -db #set point
e.TR_INIT = 0 #dead-band

有些应用程序需要一个不符合标准格式的自定义参考轨迹。自定义参考轨迹是通过创建一个新的误差(e)变量来指定的，即指定轨迹(正弦、锯齿、随机等)与模型输出之间的差异。此误差被指定为受控变量(CV)，其上、下死区表示为SPHI和SPLO。CV值也可以是带平方误差目标(e.SP=0，m.options.CV_TYPE=2)的零值，以驱动目标，而不是死区范围。

​

​

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

# initialize GEKKO model
m = GEKKO()

# time
m.time = np.linspace(0,20,41)

# constants
mass = 500

# Parameters
b = m.Param(value=50)
K = m.Param(value=0.8)
# Manipulated variable
p = m.MV(value=0, lb=-100, ub=100)

# Reference trajectory
sine = 10*np.sin(m.time/20*4*np.pi)
traj = m.Param(value=sine)

# Controlled Variable
v = m.SV(value=0,name='v')

# Error
e = m.CV(value=0,name='e')

# Equations
m.Equation(mass*v.dt() == -v*b + K*b*p)
m.Equation(e==v-traj)

m.options.IMODE = 6 # control

# MV tuning
p.STATUS = 1 #allow optimizer to change
p.DCOST = 0.1 #smooth out MV
p.DMAX = 50 #slow down change of MV

# CV tuning
e.STATUS = 1 #add the CV to the objective
m.options.CV_TYPE = 1 #Dead-band
db = 2
e.SPHI = db #set point
e.SPLO = -db #set point
e.TR_INIT = 0 #dead-band

# Solve
m.solve()

# get additional solution information
import json
with open(m.path+'//results.json') as f:
    results = json.load(f)

# Plot solution
plt.figure()
plt.subplot(3,1,1)
plt.plot(m.time,p.value,'b-',lw=2,label='MV')
plt.legend(loc='best')
plt.ylabel('MV')
plt.subplot(3,1,2)
plt.plot(m.time,sine+db,'k-',label='SPHI')
plt.plot(m.time,sine-db,'k-',label='SPLO')
plt.plot(m.time,v.value,'r--',lw=2,label='CV')
plt.legend(loc='best')
plt.ylabel('CV')
plt.subplot(3,1,3)
plt.plot(m.time,results['e.tr_hi'],'k-',label='SPHI')
plt.plot(m.time,results['e.tr_lo'],'k-',label='SPLO')
plt.plot(m.time,e.value,'r--',lw=2,label='Error')
plt.legend(loc='best')
plt.ylabel('Error')
plt.xlabel('time')
plt.show()