MPC in gekko - add dynamics in process

Question

I'm working on an MPC using gekko and I have face a problem to add dynamics to process. I have equations to two control properties thus once the program is running I need to add dynamic on MV, one by one. This is a non-linear process. In the future I will use a equation to calculate the gain. The following is the original code.

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)
#m.open_folder()
t= np.linspace(0,500,100)
m.time=t

m.options.CV_TYPE = 2 # squared error
m.options.IMODE = 6  #MPC
m.options.SOLVER = 2

# Parameters
Tm = np.ones(len(t))*260
Tm[60:] = 265

Tr = m.Param(value=Tm)
Ti=m.Param(value=30)
HT=m.Param(value=230)

# Manipulated variable
F = m.Param(value=0.11, lb=0, ub=0.65)

j = m.MV(value=1, lb=0, ub=40)
j.STATUS = 1  
j.DCOST = 0.1 
j.DMAX = 5   

mv1 = m.MV(value=20, lb=18, ub=22)
mv1.STATUS = 1  
mv1.DCOST = 0.1 
mv1.DMAX = 0.1  
To = 8.6576*mv1 + 112.66

#MI
sp_m = np.ones(len(t))*4
sp_m[30:] = 10
sp_m[50:]=20
O_sp = m.Param(value=sp_m)

# Controlled Variable
P2 = m.Var(value=10) 

m.Obj((O_sp-O2)**2)

P1 = m.CV(value=94)
P1.STATUS = 1 
P1.SP = 95    
P1.TR_INIT = 1 
P1.TAU = 30   
P1.FDELAY=2  
P1.WSP=1   

# Process model
m.Equation(P1==(((To-(Ti+4))/(((16.15-0.019*(To+Ti)/2))*mv1))-0.06*F)*100)
m.Equation(P2==m.exp(0.792*m.log(j)-3.09*m.log(HT)+13.87*m.log(Tr)+3.06*m.log(1+F)-59.9))

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=(10,15))

plt.subplot(5,1,1)
plt.plot(m.time,mv1,'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('MV1')

plt.subplot(5,1,2)
plt.plot(m.time,results['v2.tr'],'k-',label='Setpoint')
plt.plot(m.time,P1,'r--',label='CV Response')
plt.ylim([93, 96])
plt.ylabel('Conversion(%)')
plt.legend(loc='best')

plt.subplot(5,1,3)
plt.plot(m.time,results['p6'],'b-',label='MV Optimized')
plt.legend(loc='best')
plt.ylabel('MV2')

plt.subplot(5,1,4)
plt.plot(m.time,results['p8'],'k-',label='Setpoint')
plt.plot(m.time,results['v1'],'r--',label='CV Response')
plt.ylabel('P1(mu)')

plt.subplot(5,1,5)
plt.plot(m.time,Tm,'k-',label='Reactor temp(°C)')
plt.ylabel('T (°C)')
plt.xlabel('Time (min)')
plt.legend(loc='best')
plt.show()

For example I need to change the MV in equation

P2==m.exp(0.792*m.log(j)-3.09*m.log(HT)+13.87 \ 
          *m.log(Tr)+3.06*m.log(1+F)-59.9))

I have been tried:

Option 1: Doesn't work

U1 =0.11/(1.8*(s**2)+2.88*s+1)*m.exp(-s) #transfer function
u1 = inverse_laplace_transform(U1,s,t)
P2==m.exp(0.792*m.log(j***u1**)-3.09*m.log(HT)\
          +13.87*m.log(Tr)+3.06*m.log(1+F)-59.9))

Option 2: Doesn't work

m.Equation(ft==k*(1-m.exp(-(t-tau_d)/tau))
m.Equation(P2==m.exp(0.792*m.log(j***ft**)-3.09*m.log(HT)\
            +13.87*m.log(Tr)+3.06*m.log(1+F)-59.9)))

Result: Equation without an equality (=) or inequality (>,<)
(((-(10.1010101010101-p1)))/(p2))
(((-(15.15151515151515-p1)))/(p2))
STOPPING...

Thanks in advance

Answer 1

There are several ways to add dynamics in Gekko including ARX (Time Series), Differential equations, state space, discrete state space, and others. The dynamic equation for the MV isn't clear from the question. In the first option, there is a second order system with a one time unit delay:

U1 =0.11/(1.8*(s**2)+2.88*s+1)*m.exp(-s) #transfer function
u1 = inverse_laplace_transform(U1,s,t)

Here is material to help with getting the second order system into a differential equation form. For the time delay, use a discrete state space model or ARX model that is an additional model element.

The second option that didn't work appears to be an explicit solution to a dynamic step response:

m.Equation(ft==k*(1-m.exp(-(t-tau_d)/tau))

Instead of the step response, use the original differential equation as show in the first order systems material.

The original code that you posted does not produce the chart. There are several missing variable declarations. For future posts, include a minimal and complete working example of the problem for more specific feedback.