MPC多目标MEAS校正

Question

我正在尝试使用gekko来解决MPC问题。我的任务是开发一个MPC来控制三个不同的变量(CV1、CV2和CV3(。我有控制变量的方程(模型(，尽管我为每个变量创建了一个类似虚拟在线分析器的神经网络。我想用RNA作为测量，并用它通过MEAS(gekko(调整MPC。这个过程不是线性的，所以增益是可变的，取决于CV1和CV2的值。

我用gekko创建了一个算法来解决这个数学问题，但我不知道如何使用RNA作为MEAS，当我添加动力学时，没有找到解决方案。

如果有人能帮助我了解问题所在以及如何使用RNA模型进行偏差校正，我将不胜感激。提前谢谢。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt  
import math 
import pandas as pd
import control as ct
import numpy as np
m = GEKKO(remote=False)
t= np.linspace(0,500,100)
m.time=t
m.options.CV_TYPE = 2 
m.options.IMODE = 6  #MPC
m.options.SOLVER = 2
# Parameters
#Temperature middle reactor
T1 = np.ones(len(t))*260
T1[90:] = 265
tau_m_d= m.Param(value=1) #p1
tau_m= m.Param(value=5) #p2
tau_d_d= m.Param(value=5) #p1
tau_d= m.Param(value=15) #p2
T = m.Param(value=T1) #p3
Ti=m.Param(value=30) #p5
H1=m.Param(value=230) #p6
# Manipulated variable
#B use to control CV1
B = m.MV(value=0.12, lb=0, ub=0.8) 
B.STATUS = 1  
B.DCOST = 0.01 
B.DMAX = 0.2
Bi=B #I am using this because B is MV only for CV1
#ti to control CV2 by FODPT
ti_m = m.MV(value=1, lb=0, ub=300) #p8
ti_m.STATUS = 1  # allow optimizer to change
ti_m.DCOST = 0.1 # smooth  movement - Delta cost penalty for MV movement
ti_m.DMAX = 2   # maximum movement each cycle
ti_d = m.MV(value=1, lb=0, ub=300) #p8
ti_d.STATUS = 1  # allow optimizer to change
ti_d.DCOST = 0.01 # smooth  movement - Delta cost penalty for MV movement
ti_d.DMAX = 5   # maximum movement each cycle
# E to control CV3
E = m.MV(value=20, lb=18, ub=22) #p10
E.STATUS = 1  
E.DCOST = 0.1 
E.DMAX = 0.2
Ei=E # I am using this because E is MV only to CV3
To = 8.6576*E + 112.66 #T is function of E
#setpoint CV1
sp_d = np.ones(len(t))*940 #
sp_d[20:] =945
sp_d[60:]=930
d = m.Param(value=sp_d)
# setpoint CV2
sp_m = np.ones(len(t))*4 #
sp_m[40:] = 40
sp_m[80:]= 10
mi = m.Param(value=sp_m)
#setpoint CV3
sp_q = np.ones(len(t))*94 
sp_q[15:] =94.5
sp_q[70:]=95
q = m.Param(value=sp_q)
# Controlled Variable
h2 = m.Var(value=1, lb=0, ub=300)
CV1 = m.Var(value=940) 
CV2 = m.Var() 
DT = m.Var() 
k_2 = m.Var(value=0.4) #K by CV2
k_1 = m.Var(value=1) #K by CV1
CV3=m.Var(value=94)
m.Obj((d-CV1)**2)
m.Obj((mi-CV2)**2)
m.Obj((q-CV3)**2)
# Process model
m.Equation(h2==40*k_mi*(1-m.exp(-(ti_m-tau_m_d)/tau_m)))
m.Equation(B==0.8*k_d*(1-m.exp(-(ti_d-tau_d_d)/tau_d)))
m.Equation(CV1==(1000*(-0.149776*B+0.300058*B**2-0.237433*B**3+0.0001952*(h2)+0.000149*(T)+0.91439)))
m.Equation(CV2==m.exp(2.05433906*Bi+0.08761687*(h2)+0.19555572*Ei-0.02127295*H1+0.02838329*DT-5.59437532)) 
m.Equation(CV3==(((To-(Ti+4))/(((16.15-0.019*(To+Ti)/2))*E))-0.06*Bi)*100)
m.Equation(k_mi==CV2*0.1)
m.Equation(k_d==0.0909/B) #
m.Equation(DT==To-Ti)
m.options.IMODE = 6 # MPC
m.solve(disp=False)
import json
with open(m.path+'//results.json') as f:
results = json.load(f)
plt.figure(figsize=(15,25))
plt.subplot(6,1,1)
plt.plot(m.time,B.value,'b-',label='MV Optimized')
plt.ylabel('B')
plt.legend(loc='best')
plt.subplot(7,1,2)
plt.plot(m.time,results['v2'],'k-',label='CV Response')
plt.plot(m.time,results['p12'],'r--',label='Setpoint')
plt.legend(loc='best')
plt.ylabel('CV1')
plt.subplot(7,1,3)
plt.plot(m.time,E,'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('E')
plt.subplot(7,1,4)
plt.plot(m.time,q,'r--',label='CV3 Setpoint')
plt.plot(m.time,results['v7'],'b-',label='CV3 Response')
plt.ylim([93, 96])
plt.ylabel('CV3')
plt.legend(loc='best')
plt.subplot(7,1,5)
plt.plot(m.time,results['v1'],'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('H2')
plt.subplot(7,1,6)
plt.plot(m.time,results['v3'],'k-',label='CV Response')
plt.plot(m.time,results['p13'],'r--',label='Setpoint')
plt.plot(m.time,pred,'b-',label='RNA')
plt.ylabel('CV2')
plt.legend(loc='best')
plt.subplot(7,1,7)
plt.plot(m.time,T,'k-',label='Tmean')
plt.ylabel('T (°C)')
plt.xlabel('Tempo (min)')
plt.legend(loc='best')
plt.show()

Answer 1

需要通过创建新变量或删除方程来调整自由度。以下是指定k_mi = 1和k_d = m.Var()(以前未定义(。

----------------------------------------------------------------
APMonitor, Version 1.0.0
APMonitor Optimization Suite
----------------------------------------------------------------


--------- APM Model Size ------------
Each time step contains
Objects      :  0
Constants    :  0
Variables    :  22
Intermediates:  0
Connections  :  0
Equations    :  11
Residuals    :  11

Number of state variables:    1584
Number of total equations: -  1584
Number of slack variables: -  0
---------------------------------------
Degrees of freedom       :    0

----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------

Iter    Objective  Convergence
0  5.49691E+05  2.55812E+02
1  1.93664E+05  1.57614E-01
2  1.93664E+05  1.11022E-16
No feasible solution

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :  0.09949999999999999 sec
Objective      :  62217.828791846914
Unsuccessful with error code  0
---------------------------------------------------

Creating file: infeasibilities.txt
Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
@error: Solution Not Found

使用m.options.COLDSTART=1时会出现错误消息。此模式关闭MV，并确定方程/变量不平衡的位置。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt  
import math 
import pandas as pd
import numpy as np
m = GEKKO(remote=False)
t= np.linspace(0,500,100)
m.time=t
m.options.CV_TYPE = 2 
m.options.IMODE = 6  #MPC
m.options.SOLVER = 1
# Parameters
#Temperature middle reactor
T1 = np.ones(len(t))*260
T1[90:] = 265
tau_m_d= m.Param(value=1) #p1
tau_m= m.Param(value=5) #p2
tau_d_d= m.Param(value=5) #p1
tau_d= m.Param(value=15) #p2
T = m.Param(value=T1) #p3
Ti=m.Param(value=30) #p5
H1=m.Param(value=230) #p6
k_mi = 1
k_d = m.Var()
# Manipulated variable
#B use to control CV1
B = m.MV(value=0.12, lb=0, ub=0.8) 
B.STATUS = 1  
B.DCOST = 0.01 
B.DMAX = 0.2
Bi=B #I am using this because B is MV only for CV1
#ti to control CV2 by FODPT
ti_m = m.MV(value=1, lb=0, ub=300) #p8
ti_m.STATUS = 1  # allow optimizer to change
ti_m.DCOST = 0.1 # smooth  movement - Delta cost penalty for MV movement
ti_m.DMAX = 2   # maximum movement each cycle
ti_d = m.MV(value=1, lb=0, ub=300) #p8
ti_d.STATUS = 1  # allow optimizer to change
ti_d.DCOST = 0.01 # smooth  movement - Delta cost penalty for MV movement
ti_d.DMAX = 5   # maximum movement each cycle
# E to control CV3
E = m.MV(value=20, lb=18, ub=22) #p10
E.STATUS = 1  
E.DCOST = 0.1 
E.DMAX = 0.2
Ei=E # I am using this because E is MV only to CV3
To = 8.6576*E + 112.66 #T is function of E
#setpoint CV1
sp_d = np.ones(len(t))*940 #
sp_d[20:] =945
sp_d[60:]=930
d = m.Param(value=sp_d)
# setpoint CV2
sp_m = np.ones(len(t))*4 #
sp_m[40:] = 40
sp_m[80:]= 10
mi = m.Param(value=sp_m)
#setpoint CV3
sp_q = np.ones(len(t))*94 
sp_q[15:] =94.5
sp_q[70:]=95
q = m.Param(value=sp_q)
# Controlled Variable
h2 = m.Var(value=1, lb=0, ub=300)
CV1 = m.Var(value=940) 
CV2 = m.Var() 
DT = m.Var() 
k_2 = m.Var(value=0.4) #K by CV2
k_1 = m.Var(value=1) #K by CV1
CV3=m.Var(value=94)
m.Obj((d-CV1)**2)
m.Obj((mi-CV2)**2)
m.Obj((q-CV3)**2)
# Process model
m.Equation(h2==40*k_mi*(1-m.exp(-(ti_m-tau_m_d)/tau_m)))
m.Equation(B==0.8*k_d*(1-m.exp(-(ti_d-tau_d_d)/tau_d)))
m.Equation(CV1==(1000*(-0.149776*B+0.300058*B**2-0.237433*B**3+0.0001952*(h2)+0.000149*(T)+0.91439)))
m.Equation(CV2==m.exp(2.05433906*Bi+0.08761687*(h2)+0.19555572*Ei-0.02127295*H1+0.02838329*DT-5.59437532)) 
m.Equation(CV3==(((To-(Ti+4))/(((16.15-0.019*(To+Ti)/2))*E))-0.06*Bi)*100)
m.Equation(k_mi==CV2*0.1)
m.Equation(k_d==0.0909/B) #
m.Equation(DT==To-Ti)
m.options.IMODE = 6 # MPC
m.options.COLDSTART = 1
m.solve(disp=True)
import json
with open(m.path+'//results.json') as f:
results = json.load(f)
plt.figure(figsize=(15,25))
plt.subplot(6,1,1)
plt.plot(m.time,B.value,'b-',label='MV Optimized')
plt.ylabel('B')
plt.legend(loc='best')
plt.subplot(7,1,2)
plt.plot(m.time,results['v2'],'k-',label='CV Response')
plt.plot(m.time,results['p12'],'r--',label='Setpoint')
plt.legend(loc='best')
plt.ylabel('CV1')
plt.subplot(7,1,3)
plt.plot(m.time,E,'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('E')
plt.subplot(7,1,4)
plt.plot(m.time,q,'r--',label='CV3 Setpoint')
plt.plot(m.time,results['v7'],'b-',label='CV3 Response')
plt.ylim([93, 96])
plt.ylabel('CV3')
plt.legend(loc='best')
plt.subplot(7,1,5)
plt.plot(m.time,results['v1'],'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('H2')
plt.subplot(7,1,6)
plt.plot(m.time,results['v3'],'k-',label='CV Response')
plt.plot(m.time,results['p13'],'r--',label='Setpoint')
plt.plot(m.time,pred,'b-',label='RNA')
plt.ylabel('CV2')
plt.legend(loc='best')
plt.subplot(7,1,7)
plt.plot(m.time,T,'k-',label='Tmean')
plt.ylabel('T (°C)')
plt.xlabel('Tempo (min)')
plt.legend(loc='best')
plt.show()

Gekko需要一个模型的方程来建立自动微分。它目前不允许使用黑盒模型。要将模型预测用作偏差更新，请创建一个新变量，例如：

bias = m.Param()
h2_bias = m.Var()
m.Equation(h2_bias==h2+bias)

这样，外部模型可以更新bias值，例如：

bias.value = h2_meas - h2.value[0]

这发生在优化例程之外，因此bias是一个固定值。