python-3.x - 为什么延迟参数会影响预测曲线？

标签 python-3.x gekko

我将Gekko的计算结果放入队列中，延迟一段时间后，将队列的第一个值写入TCLab Arduino。我用这种方法模拟了一个工厂大延时系统，然后优化Gekko参数以达到更好的控制效果。当我在模型中添加delay=1时，我得到了一个非常好的预测曲线: enter image description here

最终的控制效果也还不错: enter image description here

但是当我设置delay=80，没有修改其他参数时，预测曲线并不理想: enter image description here

最终的控制效果也不好: enter image description here

为什么延迟参数会影响预测曲线？我认为引用轨迹也应该随着时间延迟而移动。我该如何解决这个问题？

这是代码:

import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from gekko import GEKKO
import _thread
import queue
import threading
import json

# Connect to Arduino
a = tclab.TCLabModel()

R = threading.Lock()

# Turn LED on
print('LED On')
a.LED(100)

# Simulate a time delay
delay = 80

# Run time in minutes
run_time = 60.0

# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)

# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 40.0 # set point (degC)

# heater values
Q1s = np.ones(loops) * 0.0

#########################################################
# Initialize Model
#########################################################
# use remote=True for MacOS
m = GEKKO(name='tclab-mpc',remote=False)

# 240 second time horizon
m.time = np.linspace(0,240,121)

# Parameters
Q1_ss = m.Param(value=0)
TC1_ss = m.Param(value=a.T1)
Kp = m.Param(value=1.0)
tau = m.Param(value=120.0)

# Manipulated variable
Q1 = m.MV(value=Q1_ss.VALUE, name='q1')
Q1.STATUS = 1  # use to control temperature
Q1.FSTATUS = 0 # no feedback measurement
Q1.LOWER = 0.0
Q1.UPPER = 100.0
Q1.DMAX = 50
Q1.DCOST = 2.0

# Controlled variable
TC1 = m.CV(value=TC1_ss.VALUE, name='tc1')
TC1.STATUS = 1     # minimize error with setpoint range
TC1.FSTATUS = 1    # receive measurement
TC1.TR_INIT = 2    # reference trajectory
TC1.TAU = 30      # time constant for response

# 添加延时
Q1d=m.Var()
m.delay(Q1, Q1d, delay)
# Equation
m.Equation(tau * TC1.dt() + (TC1-TC1_ss) == Kp * Q1d)

# Steady-state initialization
m.options.IMODE = 1
m.solve()

# Global Options
m.options.IMODE   = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES   = 2 # Collocation nodes
m.options.SOLVER  = 1 # 1=APOPT, 3=IPOPT
##################################################################

# Create plot
plt.figure()
plt.ion()
plt.show()

filter_tc1 = []

q = queue.Queue()
flag = False
def setQ1():
    global flag
    global a
    last_time = time.time()
    while True:
        sleep_max = 2.0
        sleep = sleep_max - (time.time() - last_time)
        print("q.qsize()", q.qsize())
        if sleep >= 0.01:
            time.sleep(sleep)
        else:
            time.sleep(0.01)

        # Record time and change in time
        last_time = time.time()
        while q.qsize() >= delay:
            q1 = q.get()
            print(f'Q1: {q1} ')
            try:
                R.acquire()
                print("write Q1：", a.Q1(float(q1)))
                R.release()
            except:
                print("convertion error!")

_thread.start_new_thread(setQ1, ())

# Rolling average filtering
def movefilter(predata, new, n):
    if len(predata) < n:
        predata.append(new)
    else:
        predata.pop(0)
        predata.append(new)
    return np.average(predata)

# Main Loop
start_time = time.time()
prev_time = start_time
last_Q1 = 0
try:
    for i in range(1,loops):
        # Sleep time
        sleep_max = 2.0
        sleep = sleep_max - (time.time() - prev_time)
        print(time.time() - prev_time)
        if sleep>=0.01:
            time.sleep(sleep)
        else:
            time.sleep(0.01)

        # Record time and change in time
        t = time.time()
        dt = t - prev_time
        prev_time = t
        tm[i] = t - start_time

        # Read temperatures in Kelvin
        R.acquire()
        curr_T1 = a.T1
        R.release()
        last_T1 = curr_T1
        avg_T1 = movefilter(filter_tc1, last_T1, 3)
        T1[i] = curr_T1
        # T1[i] = a.T1
        # T2[i] = a.T2

        ###############################
        ### MPC CONTROLLER          ###
        ###############################
        TC1.MEAS = avg_T1
        # input setpoint with deadband +/- DT
        DT = 1.0
        TC1.SPHI = Tsp1[i] + DT
        TC1.SPLO = Tsp1[i] - DT
        try:
            # solve MPC
            m.solve(disp=False)
        except:
            pass
        # test for successful solution
        if (m.options.APPSTATUS==1):
            # retrieve the first Q value
            Q1s[i] = Q1.NEWVAL
            with open(m.path+'//results.json') as f:
                results = json.load(f)
        else:
            # not successful, set heater to zero
            Q1s[i] = 0

        # Write output (0-100)
        q.put(Q1s[i])

        # Plot
        plt.clf()
        ax=plt.subplot(2,1,1)
        ax.grid()
        plt.plot(tm[0:i],T1[0:i],'ro',markersize=3,label=r'$T_1$')
        plt.plot(tm[0:i],Tsp1[0:i],'b-',markersize=3,label=r'$T_1 Setpoint$')
        plt.plot(tm[i]+m.time,results['tc1.bcv'],'k-.',\
                 label=r'$T_1$ predicted',linewidth=3)
        plt.plot(tm[i]+m.time,results['tc1.tr_hi'],'k--',\
                 label=r'$T_1$ trajectory')
        plt.plot(tm[i]+m.time,results['tc1.tr_lo'],'k--')
        plt.ylabel('Temperature (degC)')
        plt.legend(loc='best')
        ax=plt.subplot(2,1,2)
        ax.grid()
        plt.plot(tm[0:i],Q1s[0:i],'r-',linewidth=3,label=r'$Q_1$')
        plt.plot(tm[i]+m.time,Q1.value,'k-.',\
                 label=r'$Q_1$ plan',linewidth=3)
        # plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
        plt.ylabel('Heaters')
        plt.xlabel('Time (sec)')
        plt.legend(loc='best')
        plt.draw()
        plt.pause(0.05)

    # Turn off heaters
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()

# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()

# Make sure serial connection still closes when there's an error
except:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Error: Shutting down')
    a.close()
    raise

最佳答案

Controller 正在按预期工作。该问题是众所周知的模型不匹配和延迟问题。经典的控制方法是使用史密斯预测器来补偿延迟。使用 MPC，延迟被内置到模型中，类似于史密斯预测器的功能。如果将模型更改为:

Kp = m.Param(value=0.7)
tau = m.Param(value=180.0)

那么性能就提高了:

improved control performance

对引用轨迹的修改也可以提高性能，因为 Controller 专注于长期目标而不是它无法控制的近期误差。

方法 1 - 在开始时打开引用轨迹

TC1.TR_INIT = 2    # reference trajectory
TC1.TR_OPEN = 20   # more focus on final destination
TC1.TAU = 60       # time constant for response

方法 2 - 稍后开始惩罚

使用CV_WGT_START仅在一定时间延迟后进行处罚。

TC1.TR_INIT = 0    # reference trajectory
m.options.CV_WGT_START = 80

参数需要再调整一次，以更好地匹配工艺。

Kp = m.Param(value=0.8)
tau = m.Param(value=160.0)

improved control performance

方法三:定义自己的引用轨迹

如果您确实想要自定义引用轨迹，则可以定义一个新变量 SP，该变量延迟为 SPd。定义一个新的 CV error，它是 TC1 与延迟引用轨迹之间的差异。

SP = m.Var()
m.Equation(30 * SP.dt() == -SP + 40)
SPd = m.Var()
m.delay(SP, SPd, delay)

error = m.CV(value=0)
m.Equation(error==SPd - TC1)
error.STATUS = 1     # minimize error with setpoint range
error.FSTATUS = 0    # receive measurement
error.TR_INIT = 0    # reference trajectory
error.SPHI    = 0.5  # drive error close to 0
error.SPLO    = -0.5 # drive error close to 0

使用这种方法，您需要找出一种方法来更新测量反馈的误差测量值。

也有很多时候 MPC 计算会稍微超出循环时间。这会造成模型不匹配，您可以通过将 MPC 循环时间增加到 3 秒来解决这一问题。

修改后的脚本

import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from gekko import GEKKO
import _thread
import queue
import threading
import json

# Connect to Arduino
a = tclab.TCLabModel()

R = threading.Lock()

# Turn LED on
print('LED On')
a.LED(100)

# Simulate a time delay
delay = 80

# Run time in minutes
run_time = 60.0

# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)

# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 40.0 # set point (degC)

# heater values
Q1s = np.ones(loops) * 0.0

#########################################################
# Initialize Model
#########################################################
# use remote=True for MacOS
m = GEKKO(name='tclab-mpc',remote=False)

# 240 second time horizon
m.time = np.linspace(0,240,121)

# Parameters
Q1_ss = m.Param(value=0)
TC1_ss = m.Param(value=a.T1)
Kp = m.Param(value=0.7)
tau = m.Param(value=180.0)

# Manipulated variable
Q1 = m.MV(value=Q1_ss.VALUE, name='q1')
Q1.STATUS = 1  # use to control temperature
Q1.FSTATUS = 0 # no feedback measurement
Q1.LOWER = 0.0
Q1.UPPER = 100.0
Q1.DMAX = 50
Q1.DCOST = 2.0

# Controlled variable
TC1 = m.CV(value=TC1_ss.VALUE, name='tc1')
TC1.STATUS = 1     # minimize error with setpoint range
TC1.FSTATUS = 1    # receive measurement
TC1.TR_INIT = 0    # reference trajectory
m.options.CV_WGT_START = 80

# 添加延时
Q1d=m.Var()
m.delay(Q1, Q1d, delay)
# Equation
m.Equation(tau * TC1.dt() + (TC1-TC1_ss) == Kp * Q1d)

# Steady-state initialization
m.options.IMODE = 1
m.solve()

# Global Options
m.options.IMODE   = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES   = 2 # Collocation nodes
m.options.SOLVER  = 1 # 1=APOPT, 3=IPOPT
##################################################################

# Create plot
plt.figure()
plt.ion()
plt.show()

filter_tc1 = []

q = queue.Queue()
flag = False
def setQ1():
    global flag
    global a
    last_time = time.time()
    while True:
        sleep_max = 2.0
        sleep = sleep_max - (time.time() - last_time)
        print("q.qsize()", q.qsize())
        if sleep >= 0.01:
            time.sleep(sleep)
        else:
            time.sleep(0.01)

        # Record time and change in time
        last_time = time.time()
        while q.qsize() >= delay:
            q1 = q.get()
            print(f'Q1: {q1} ')
            try:
                R.acquire()
                print("write Q1：", a.Q1(float(q1)))
                R.release()
            except:
                print("convertion error!")

_thread.start_new_thread(setQ1, ())

# Rolling average filtering
def movefilter(predata, new, n):
    if len(predata) < n:
        predata.append(new)
    else:
        predata.pop(0)
        predata.append(new)
    return np.average(predata)

# Main Loop
start_time = time.time()
prev_time = start_time
last_Q1 = 0
try:
    for i in range(1,loops):
        # Sleep time
        sleep_max = 2.0
        sleep = sleep_max - (time.time() - prev_time)
        print(time.time() - prev_time)
        if sleep>=0.01:
            time.sleep(sleep)
        else:
            print('Warning: Dynamic mismatch from excess MPC time')
            print('         Consider increasing the cycle time to 3+ sec')
            time.sleep(0.01)

        # Record time and change in time
        t = time.time()
        dt = t - prev_time
        prev_time = t
        tm[i] = t - start_time

        # Read temperatures in Kelvin
        R.acquire()
        curr_T1 = a.T1
        R.release()
        last_T1 = curr_T1
        avg_T1 = movefilter(filter_tc1, last_T1, 3)
        T1[i] = curr_T1
        # T1[i] = a.T1
        # T2[i] = a.T2

        ###############################
        ### MPC CONTROLLER          ###
        ###############################
        TC1.MEAS = avg_T1
        # input setpoint with deadband +/- DT
        DT = 1.0
        TC1.SPHI = Tsp1[i] + DT
        TC1.SPLO = Tsp1[i] - DT
        try:
            # solve MPC
            m.solve(disp=False)
        except:
            pass
        # test for successful solution
        if (m.options.APPSTATUS==1):
            # retrieve the first Q value
            Q1s[i] = Q1.NEWVAL
            with open(m.path+'//results.json') as f:
                results = json.load(f)
        else:
            # not successful, set heater to zero
            Q1s[i] = 0

        # Write output (0-100)
        q.put(Q1s[i])

        # Plot
        plt.clf()
        ax=plt.subplot(2,1,1)
        ax.grid()
        plt.plot(tm[0:i],T1[0:i],'ro',markersize=3,label=r'$T_1$')
        plt.plot(tm[0:i],Tsp1[0:i],'b-',markersize=3,label=r'$T_1 Setpoint$')
        plt.plot(tm[i]+m.time,results['tc1.bcv'],'k-.',\
                 label=r'$T_1$ predicted',linewidth=3)
        plt.plot(tm[i]+m.time,results['tc1.tr_hi'],'k--',\
                 label=r'$T_1$ trajectory')
        plt.plot(tm[i]+m.time,results['tc1.tr_lo'],'k--')
        plt.ylabel('Temperature (degC)')
        plt.legend(loc='best')
        ax=plt.subplot(2,1,2)
        ax.grid()
        plt.plot(tm[0:i],Q1s[0:i],'r-',linewidth=3,label=r'$Q_1$')
        plt.plot(tm[i]+m.time,Q1.value,'k-.',\
                 label=r'$Q_1$ plan',linewidth=3)
        # plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
        plt.ylabel('Heaters')
        plt.xlabel('Time (sec)')
        plt.legend(loc='best')
        plt.draw()
        plt.pause(0.05)

    # Turn off heaters
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()

# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()

# Make sure serial connection still closes when there's an error
except:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Error: Shutting down')
    a.close()
    raise

关于python-3.x - 为什么延迟参数会影响预测曲线？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/66326361/

上一篇：javascript - 倒计时它不适用于单击按钮

下一篇：google-cloud-platform - 从 Kafka 读取数据流而不会丢失数据？

相关文章：

python-3.x - botocore.exceptions.ProfileNotFound : The config profile (default) could not be found

python - 如何在没有已知目标函数(比如一些随机函数)和已知变量和约束的情况下使用 gekko 优化器？

python - 如何在gekko中的最大迭代限制后获取决策变量的值

python - 评估随机长度的表达式

python - 如何在 Python3 中将字符串数字列表转换为 int 数字？

python - 使用 r2pipe 进行多处理

python - 在 Python 中格式化 JSON

python - 如何在 Python Gekko 中将变量值设置为 x[3]=6(不是初始条件)？

python - 在模式 6 python gekko 下运行时没有求解器显示

python - 从 MATLAB 调用 Python 运算符错误 x**2

©2024 IT工具网  联系我们