Why MHE after 12 seconds runs skeptically over my model?

Question

I am new to control systems and have tried to use GEKKO. For my model simulation geos correct and run however at MHE and also at MPC It runs correctly when t = 15 seconds after that if t = 16 seconds at the 12th it goes skeptical as shown on the attached figure. Is there a reason for that?

m = GEKKO(remote = False)

g  = m.Const(value = 9.81)
Cc = m.Const(value = 2*10**-5)
D1 = m.Const(value = 0.1016)
D2 = m.Const(value = 0.1016)
h1 = m.Const(value = 100) 
hv = m.Const(value = 1000)
L1 = m.Const(value = 500)
L2 = m.Const(value = 1100)
V1 = m.Const(value = 4.054)
V2 = m.Const(value = 9.729)
A1 = m.Const(0.008107)
A2 = m.Const(0.008107)

beta1 = m.Const(value = 1.5 * 10**9)
beta2 = m.Const(value = 1.5* 10**9)
Pres = m.Const(value = 1.261 * 10**7)
M = m.Const(value = 1.992 * 10**8) 
rho = m.Const(value = 932)
PI = m.Const(value = 2.32 * 10**(-9))
visc = m.Const(value = 0.024)

fo = m.Const(value = 60)
Inp = m.Const(value = 65)
Pnp = m.Const(value = 1.625 * 10**5)

z = m.MV(1, lb = 0.1, ub = 1)

f = m.MV(35, lb = 35, ub = 65)
f.STATUS = 1
f.DMAX = 30   
f.DCOST = 0.0002 # slow down change of frequency
f.UPPER = 65

Ppin = m.CV()
Ppin.STATUS = 1 # add the SP to the objective
Ppin.SP = 6e6 # set point
Ppin.TR_INIT = 1 # set point trajectory
Ppin.TAU = 5 # time constant of trajectory

Pwf = m.Var(ub = Pres)
Ppout = m.Var()
Pwh = m.Var()
Pm = m.Var()

P_fric_drop = m.Var()
F1 = m.Var()
F2 = m.Var()

q = m.Var(lb = 0)
qc = m.Var(lb = 0)
qr = m.Var(lb = 0)

I = m.Var()
H = m.Var()
BHP = m.Var(lb = 0)
BHPo = m.Var()
qo = m.Var()
Ho = m.Var()

m.Equation([Pwh.dt() == beta2*(q-qc)/V2,
            Pwf.dt() == beta1*(qr-q)/V1,
            q.dt() == (1/M)*(Pwf - Pwh - rho*g*hv - P_fric_drop + (Ppout-Ppin))])

m.Equation([qr == PI*(Pres - Pwf),
            qc == Cc*z*(m.sqrt((Pwh - Pm))),
            Pm == Pwh/2,
            P_fric_drop == F1 + F2,
            F1 == 0.158 * rho * L1* q**2 * visc**0.5 /(D1 * A1**2 * (rho*D1*q)**0.5),
            F2 == 0.158 * rho * L2* q**2 * visc**0.5 /(D2 * A2**2 * (rho*D2*q)**0.5)])

m.Equation([
              qo == q * fo/f,
              H == Ho * (f/fo)**2,
              BHPo == -2.3599e9*qo**3 - 1.8082e7*qo**2 + 4.3346e6*qo + 9.4355e4,
              BHP == BHPo * (f/fo)**3,
              Ppin == Pwf - rho*g*h1 - F1,
              Ho == 9.5970e2+7.4959e3*qo+-1.2454e6*qo**2,
              I == Inp * BHP / Pnp,
              Ppout == H*rho*g + Ppin
             ])

m.options.solver = 1
m.options.MAX_ITER = 250   

m.options.IMODE = 1
m.solve(debug=True)

tf = 30  # final time (sec)
st = 2.0   # simulation time intervals
nt = int(tf/st)+1 # simulation time points
m.time = np.linspace(0,tf,nt)

m.options.CV_TYPE = 2 # squared error
m.options.solver = 3
m.options.NODES=2
m.options.IMODE = 5
m.solve(debug = False)

Answer 1

I wasn't able to reproduce the problem you encountered so I'll give a few specific suggestions for troubleshooting your application and improving the ability of the solver to find a solution.

Help the Solver Find a Solution

• Avoid sqrt(-negative)

#qc == Cc*z*(m.sqrt((Pwh - Pm)))
qc**2 == Cc**2 * z**2 * (Pwh - Pm) # avoid negative in sqrt

• Avoid potential for divide by zero. The value of q is constrained away from zero but look for other equations that could be improved by re-arrangement.

#qo == q * fo/f
qo * f == q * fo

#H == Ho * (f/fo)**2
H * fo**2 == Ho * f**2

#I == Inp * BHP / Pnp
I * Pnp == Inp * BHP

• Constrain the FVs and MVs

It can often help to put small DMAX or limit the upper and lower bound for FVs and MVs. Constraints on the Var types can lead to infeasible solutions.

• Reduce the degrees of freedom

For MHE applications, I recommend that you use FVs instead of MVs for your unknown parameters. The FVs calculate a single value over all of the measurements. The MVs calculate a new parameter value for every time point. Using an FV instead of an MV can help the solver because there are fewer decisions for the optimizer and the parameter values won't be adjusted to track signal noise.

Check the Solutions

• Configure the objective function.

For MHE applications, you need to have some CV values that you are trying to match that have FSTATUS=1. In the example you posted there are no CVs with FSTATUS=1 or any inserted measurements. MHE is to align a model with measured CVs values by adjusting the unknown parameters as FVs or MVs.

• Check the solver status to make sure that it found a successful solution with APPSTAUTS==1

if m.options.APPSTATUS==1:
    print('Successful Solution)
else:
    print('Solver Failed to Find a Solution')

• Create plots of the variables, especially if there is a problematic period. You can often see where there is an integrating variable such as qr or qc that are approaching a variable constraint. Here is an example where I converted your MHE application into an MPC application. I'm adjusting the friction factor f (not realistic) and the valve opening z to show how the MPC drives Ppin to the target set point.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote = False)

g  = m.Const(value = 9.81)
Cc = m.Const(value = 2*10**-5)
D1 = m.Const(value = 0.1016)
D2 = m.Const(value = 0.1016)
h1 = m.Const(value = 100) 
hv = m.Const(value = 1000)
L1 = m.Const(value = 500)
L2 = m.Const(value = 1100)
V1 = m.Const(value = 4.054)
V2 = m.Const(value = 9.729)
A1 = m.Const(0.008107)
A2 = m.Const(0.008107)

beta1 = m.Const(value = 1.5 * 10**9)
beta2 = m.Const(value = 1.5* 10**9)
Pres = m.Const(value = 1.261 * 10**7)
M = m.Const(value = 1.992 * 10**8) 
rho = m.Const(value = 932)
PI = m.Const(value = 2.32 * 10**(-9))
visc = m.Const(value = 0.024)

fo = m.Const(value = 60)
Inp = m.Const(value = 65)
Pnp = m.Const(value = 1.625 * 10**5)

z = m.MV(1, lb = 0.1, ub = 1)
z.STATUS = 1
z.DMAX = 0.01

f = m.MV(35, lb = 35, ub = 65)
f.STATUS = 1
f.DMAX = 30   
f.DCOST = 0.0002 # slow down change of frequency
f.UPPER = 65

Ppin = m.CV()
Ppin.STATUS = 1 # add the SP to the objective
Ppin.SP = 6e6 # set point
Ppin.TR_INIT = 1 # set point trajectory
Ppin.TAU = 5 # time constant of trajectory

Pwf = m.Var(ub = Pres)
Ppout = m.Var()
Pwh = m.Var()
Pm = m.Var()

P_fric_drop = m.Var()
F1 = m.Var()
F2 = m.Var()

q = m.Var(lb = 0)
qc = m.Var(lb = 0)
qr = m.Var(lb = 0)

I = m.Var()
H = m.Var()
BHP = m.Var(lb = 0)
BHPo = m.Var()
qo = m.Var()
Ho = m.Var()

m.Equation([Pwh.dt() == beta2*(q-qc)/V2,
            Pwf.dt() == beta1*(qr-q)/V1,
            M * q.dt() == Pwf - Pwh - rho*g*hv - P_fric_drop + (Ppout-Ppin)])

m.Equation([qr == PI*(Pres - Pwf),
            qc == Cc*z*(m.sqrt((Pwh - Pm))),
            Pm == Pwh/2,
            P_fric_drop == F1 + F2,
            F1 == 0.158 * rho * L1* q**2 * visc**0.5 /(D1 * A1**2 * (rho*D1*q)**0.5),
            F2 == 0.158 * rho * L2* q**2 * visc**0.5 /(D2 * A2**2 * (rho*D2*q)**0.5)])

m.Equation([
              qo * f == q * fo,
              H * fo**2 == Ho * f**2,
              BHPo == -2.3599e9*qo**3 - 1.8082e7*qo**2 + 4.3346e6*qo + 9.4355e4,
              BHP == BHPo * (f/fo)**3,
              Ppin == Pwf - rho*g*h1 - F1,
              Ho == 9.5970e2+7.4959e3*qo+-1.2454e6*qo**2,
              I * Pnp == Inp * BHP,
              Ppout == H*rho*g + Ppin
             ])

m.options.solver = 1
m.options.MAX_ITER = 250   

m.options.IMODE = 1
m.solve()

tf = 30  # final time (sec)
st = 2.0   # simulation time intervals
nt = int(tf/st)+1 # simulation time points
m.time = np.linspace(0,tf,nt)

m.options.CV_TYPE = 2 # squared error
m.options.solver = 3
m.options.NODES=2
m.options.IMODE = 6

plt.figure(figsize=(10,5))
for i in range(10):
    m.solve(disp=False)
    print(i,m.options.APPSTATUS)
    plt.subplot(10,1,i+1)
    plt.plot(m.time+i*st,Ppin.value)
    plt.xlim([0,50])
    plt.ylim([6e6,7e6])
plt.show()

Your application looks like it is for hydraulic pressure control in drilling. There are additional models available on GitHub. I hope you'll consider contributing your model or case studies to the open-source repository.