Gekko中受控变量(CV)无法跟踪设定值(SP)的问题排查求助

我最近在用Gekko做双水箱液位控制系统的试验，目前先尝试单独控制其中一个水箱的液位，但遇到了一个头疼的问题：受控变量h1完全不跟着设定值区间走，一直在下降。我试着调整了各类参数都没效果，想请大家帮我看看代码里是不是哪里的定义或者逻辑出问题了？

问题背景与代码说明 #

我定义了三个进水流量（其中flow2会分配到两个水箱，flow1只进水箱2，flow3只进水箱1），还有一个总出水流量。然后建立了Gekko模型，尝试通过调整flow2分配到水箱1的比例qt1_v来控制水箱1的液位h1，但运行后h1就是不跟踪设定区间（0.8*~0.95*水箱1高度），结果一直下降，和预期完全不符。

以下是完整的代码和数据：

数据与模型定义代码

import numpy as np
from gekko import GEKKO
from scipy.integrate import odeint

# Influent flows
flow1 = np.array([1.25282567, 1.18243383, 1.25445315, 1.1755483 , 1.17953059,
       1.164077  , 1.18473485, 1.18064562, 1.14043986, 1.23820012,
       1.24054883, 1.20570221, 1.21603648, 1.20606003, 1.15339559,
       1.18137919, 1.20935554, 1.20335999, 1.21147486, 1.23357391,
       1.21825142, 1.19723155, 1.21343981, 1.11376114, 1.19094124,
       1.20130132, 1.21304258, 1.18635852, 1.26465022, 1.18148956,
       1.20242284, 1.23000755, 1.26954079, 1.16466064, 1.23057836,
       1.23027471]) # goes only to tank 2

flow2 = np.array([2.59677412, 2.58811843, 2.544673  , 2.55571262, 2.57453146,
       2.59810488, 2.56841369, 2.45821758, 2.4281892 , 2.43671096,
       2.4474268 , 2.46243303, 2.43683383, 2.39242551, 2.4513853 ,
       2.45447782, 2.42802739, 2.4221822 , 2.42537811, 2.44381434,
       2.44872821, 2.44426348, 2.4158827 , 2.42840086, 2.47415572,
       2.50235597, 2.47706445, 2.44440066, 2.44883323, 2.46050179,
       2.47231121, 2.49074611, 2.48412978, 2.47831733, 2.47747716,
       2.48604111]) # goes to tank 1 and 2

flow3 = np.array([1.57011312, 1.50738452, 1.52807646, 1.57724425, 1.56187045,
       1.64383615, 1.57872569, 1.57535121, 1.5540041 , 1.51887366,
       1.58146572, 1.57476749, 1.60523332, 1.5716739 , 1.56179746,
       1.53903109, 1.58319486, 1.56445033, 1.58821278, 1.61287479,
       1.57210261, 1.6220108 , 1.5820968 , 1.55307617, 1.52589127,
       1.57522106, 1.55631243, 1.55931072, 1.48830357, 1.57045034,
       1.50157251, 1.55527039, 1.56454099, 1.56372395, 1.51546389,
       1.552174  ]) #goes only to tank 1

# Effluent
ts_out = np.array([6.53681895, 6.211645  , 6.12570737, 6.36272278, 6.47968763,
       6.88415985, 6.53169619, 6.44752413, 6.76828499, 7.24453761,
       7.3429685 , 7.51939125, 7.13905894, 7.45646591, 7.53236135,
       7.64722176, 7.36411756, 7.18427607, 7.24093132, 7.13198001,
       6.98640455, 6.64083468, 6.36015803, 6.41666819, 6.58571269,
       6.52822367, 6.79654961, 6.46079076, 6.4651886 , 6.68471489,
       7.02170868, 7.23842901, 7.3350867 , 7.32969907, 7.27846311,
       7.27119013]) # sum of both effluent tanks

# Tank t1
# Constants
t1_area = 116.89
t1_height = 15.2
#initial state
t1_h0 = 13.9

# Tank t2
# Constants
t2_area = 226.98
t2_height = 15
#initial state
t2_h0 = 13.38

#Define models
def tanks(x, t, qd_in1, q_out, qt1_v, qt2_v_out, 
          qt2_in1,  A2, 
          qt1_in2,  A1):
    
    qt1_in1 = qd_in1*qt1_v
    qt2_in2 = qd_in1*(1-qt1_v)
    
    qt2_out = qt2_v_out*q_out
    qt1_out = (1 - qt2_v_out)*q_out

    dh2dt = (qt2_in1 + qt2_in2 - qt2_out)/A2
    dh1dt = (qt1_in1 + qt1_in2 - qt1_out)/A1
    
    return dh1dt, dh2dt

# Gekko model
m = GEKKO(remote=False)
m.time = np.arange(0,12,2)

#Flow 2 goes to both tanks
qd_in1 = m.SV(value = flow2[0], name ='flow2') #flow2 to t1+t2

# Tank t1
A1 = m.Const(value = t1_area, name = 'A1') # m2
h1 = m.CV(t1_h0, name = 'h1', lb=0.6*t1_height, ub=0.95*t1_height) #m, controlled variable

qt1v = 0.8 # Initial distribuition to tank 1
qt1_v = m.MV(qt1v, lb=0, ub=1, name = 'flow2_t1_v') # manipulate flow of flow2 to t1
qt1_in1 = m.Intermediate(qd_in1*qt1_v, name = 'flow2_t1') # flow2 to t1
qt1_in2 = m.SV(flow3[0],name ='flow3')

# Tank t2
A2 = m.Const(value = t2_area, name = 'A2') # m2,Cross section Area
h2 = m.CV(t2_h0, name = 'h2', lb=0.6*t2_height, ub=0.95*t2_height) #m , the controlled variable

qt2_in1 = m.SV(flow1[0], name = 'flow1')
qt2_in2 = m.Intermediate(qd_in1*(1-qt1_v), name = 'flow2_t2') # flow2 to t2

#Out 
q_out = m.SV(ts_out[0], name = 'q_out')
qt2_vout = 0.32 # 32% of outflow belongs to t2 (later will be optimized)
qt2_v_out = m.SV(qt2_vout, name = 't2_out_v')

qt2_out = m.Intermediate(qt2_v_out*q_out, name = 't2_out') #t2 to out
qt1_out = m.Intermediate((1-qt2_v_out)*q_out, name = 't1_out') #t1 to out

m.Equation(h2.dt() == ((qt2_in1 + qt2_in2 - qt2_out)/A2))
m.Equation(h1.dt() == ((qt1_in1 + qt1_in2 - qt1_out)/A1))

m.options.IMODE = 6
m.options.NODES = 3
m.CV_TYPE=1 

#CV tuning t1
h1.STATUS = 1 #add the SP to the objective
h1.FSTATUS = 1
h1.SPHI = 0.90*t1_height
h1.SPLO = 0.80*t1_height
h1.TR_INIT = 1
h1.WSPHI   = 800
h1.WSPLO   = 800
# h1.TAU = 10 #time constant of setpoint trajectory

#MV tuning for t1
qt1_v.STATUS = 1 #allow optimizer to change
qt1_v.DCOST = 0.0000001  #smooth out movement
qt1_v.DMAXLO = - 0.3
qt1_v.DMAXHI = 0.3 

循环运行代码

# 初始化存储变量
t = np.arange(0, len(flow2)*2, 2)
t1_h_ode = np.zeros(len(t))
t2_h_ode = np.zeros(len(t))
t1_h_store = np.zeros(len(t))
t2_h_store = np.zeros(len(t))
qt1_v_store = np.zeros(len(t))
qt2_vout_store = np.full(len(t), qt2_vout)

t1_h_ode[0] = t1_h0
t2_h_ode[0] = t2_h0
t1_h_store[0] = t1_h0
t2_h_store[0] = t2_h0
qt1_v_store[0] = qt1v

x = [t1_h0, t2_h0]

for i in range(len(t)-1):
    # Simulate
    inputs = (flow2[i], ts_out[i], qt1_v_store[i], qt2_vout_store[i], 
              flow1[i],  t2_area, 
              flow3[i], t1_area)
    
    ts = [t[i],t[i+1]]
    y = odeint(tanks,x,ts,args=inputs)
    
    # Record measurements
    h1.MEAS = y[-1][0]
    h2.MEAS = y[-1][1]
    
    # Store results
    t1_h_ode[i+1] = y[-1][0]
    t2_h_ode[i+1] = y[-1][1]
    
    # Adjust initial condition for next loop
    x[0] = y[-1][0]
    x[1] = y[-1][1]
    
    m.solve(disp=False)
    
    t1_h_store[i+1] = h1.model
    t2_h_store[i+1] = h2.model
    
    qd_in1.MEAS = flow2[i+1]
        
    qt1_in2.MEAS = flow3[i+1]
    
    qt2_in1.MEAS = flow1[i+1]
    
    q_out.MEAS = ts_out[i+1]

    qt1_v_store[i+1] = qt1_v.NEWVAL

运行结果 #

运行后得到的h1变化曲线如下，完全没有跟踪设定区间，一直在下降：

我已经尝试调整了h1的权重WSPHI/WSPLO，还有MV的DCOST、DMAX等参数，但都没有改善。想请教大家，我的代码里是不是哪里的定义或者逻辑有问题？比如CV的配置、MV的设置，或者是模型方程的问题？

备注：内容来源于stack exchange，提问作者iarima