I am trying to use gekko to solve a MPC problem. My task is to develop a MPC to control three diferents variables (CV1, CV2 and CV3). I have equations (models) to Control variables, although I created a Neural Network like a Virtual Online Anlyser for each one. I would like to use the RNA as measurement and use it to adjust the MPC by MEAS (gekko). This process is no linear, so the gain is variable depend on CV1 and CV2 value.
I created an algorithm to solve this mathematical problem using gekko but I don't know how to use RNA as MEAS and when I add dynamics the solution is no found.
I do appreciate if someone can help me understand what the problem is and how to do the bias correction using RNA model. Thanks in advance.
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()
