How to solve issue of maximum iteration getting exceeded in python gekko in this case (explained in the body)?

user_16462828 0
Viewed 308

I am trying to track multiple set-points in the case of interacting quadruple tank system process. Here, the upper limits of tanks are 25 and lower limits are 0. I want to track the set-point values of 5,12,7 and 5. Although, I am able to track the initial 3 set-points (i.e. 5,12 and 7), I am not able to track the last set-point due to solver exceeding max. iterations. I have attached the code below->

#MHE+MPC model

#to measure computational time of the code
start=time.time()

#Process Model
p = GEKKO(remote=False)
process=0
p.time = [0,0.5]
noise = 0.25

#Constants
g = 981
g1 = .9
g2 = .9
A1=32
A3=32
A2=32
A4=32
a1=0.057
a3=0.057
a2=0.057
a4=0.057
init_h=5

#Controlled process variables 
p.h1=p.SV(lb=0,ub=25)
p.h2=p.SV(lb=0,ub=25)
p.h3=p.SV(lb=0,ub=25)
p.h4=p.SV(lb=0,ub=25)

#Manipulated process variables
p.v1=p.MV(value=3.15,lb=0.1,ub=8)
p.v2=p.MV(value=3.15,lb=0.1,ub=8)

#Parameters of process
p.k1=p.Param(value=3.14,lb=0,ub=10)
p.k2=p.Param(value=3.14,lb=0,ub=10)

#Equations process
p.Equation(A1*p.h1.dt()==a3*((2*g*p.h3)**0.5)-(a1*((2*g*p.h1)**0.5))+(g1*p.k1*p.v1))
p.Equation(A2*p.h2.dt()==a4*((2*g*p.h4)**0.5)-(a2*((2*g*p.h2)**0.5))+(g2*p.k2*p.v2))
p.Equation(A3*p.h3.dt()==-a3*((2*g*p.h3)**0.5)+((1-g2)*p.k2*p.v2))
p.Equation(A4*p.h4.dt()==-a4*((2*g*p.h4)**0.5)+((1-g1)*p.k1*p.v1))

#options
p.options.IMODE = 4
#p.h1.TAU=-10^10
#p.h2.TAU=-10^10


#%% MHE Model
m = GEKKO(remote=False)

#prediction horizon
m.time = np.linspace(0,40,41) #0-20 by 0.5 -- discretization must match simulation

#MHE control, manipulated variables and parameters
m.h1=m.CV(lb=0,ub=25)
m.h2=m.CV(lb=0,ub=25)
m.h3=m.SV(lb=0,ub=25)
m.h4=m.SV(lb=0,ub=25)

m.v1=m.MV(value=3.15,lb=0.10,ub=8)
m.v2=m.MV(value=3.15,lb=0.10,ub=8)

m.k1=m.FV(value=3.14,lb=0,ub=10)
m.k2=m.FV(value=3.14,lb=0,ub=10)

#m.h1.TAU=0
#m.h2.TAU=0

#Equations
m.Equation(A1*m.h1.dt()==a3*((2*g*m.h3)**0.5)-(a1*((2*g*m.h1)**0.5))+(g1*m.k1*m.v1))
m.Equation(A2*m.h2.dt()==a4*((2*g*m.h4)**0.5)-(a2*((2*g*m.h2)**0.5))+(g2*m.k2*m.v2))
m.Equation(A3*m.h3.dt()==-a3*((2*g*m.h3)**0.5)+((1-g2)*m.k2*m.v2))
m.Equation(A4*m.h4.dt()==-a4*((2*g*m.h4)**0.5)+((1-g1)*m.k1*m.v1))

#Options
m.options.IMODE = 5 #MHE
m.options.EV_TYPE = 2

# STATUS = 0, optimizer doesn't adjust value
# STATUS = 1, optimizer can adjust
m.v1.STATUS = 0
m.v2.STATUS = 0
m.k1.STATUS=1
m.k2.STATUS=1

m.h1.STATUS = 1
m.h2.STATUS = 1
#m.h3.STATUS = 0
#m.h4.STATUS = 0

# FSTATUS = 0, no measurement
# FSTATUS = 1, measurement used to update model
m.v1.FSTATUS = 1
m.v2.FSTATUS = 1
m.k1.FSTATUS=0
m.k2.FSTATUS=0

m.h1.FSTATUS = 1
m.h2.FSTATUS = 1
m.h3.FSTATUS = 1
m.h4.FSTATUS = 1

#m.options.MAX_ITER=1000
m.options.SOLVER=3
m.options.NODES=3

#%% MPC Model
c = GEKKO(remote=False)

c.time = np.linspace(0,10,11) #0-5 by 0.5 -- discretization must match simulation

c.v1=c.MV(value=3.15,lb=0.10,ub=8)
c.v2=c.MV(value=3.15,lb=0.10,ub=8)

c.k1=c.FV(value=3.14,lb=0,ub=10)
c.k2=c.FV(value=3.14,lb=0,ub=10)
#Variables
c.h1=c.CV(lb=0,ub=25)
c.h2=c.CV(lb=0,ub=25)
c.h3=c.SV(lb=0,ub=25)
c.h4=c.SV(lb=0,ub=25)
#Equations
c.Equation(A1*c.h1.dt()==a3*((2*g*c.h3)**0.5)-(a1*((2*g*c.h1)**0.5))+(g1*c.k1*c.v1))
c.Equation(A2*c.h2.dt()==a4*((2*g*c.h4)**0.5)-(a2*((2*g*c.h2)**0.5))+(g2*c.k2*c.v2))
c.Equation(A3*c.h3.dt()==-a3*((2*g*c.h3)**0.5)+((1-g2)*c.k2*c.v2))
c.Equation(A4*c.h4.dt()==-a4*((2*g*c.h4)**0.5)+((1-g1)*c.k1*c.v1))

#Options
c.options.IMODE = 6 #MPC
c.options.CV_TYPE = 2

# STATUS = 0, optimizer doesn't adjust value
# STATUS = 1, optimizer can adjust
c.v1.STATUS = 1
c.v2.STATUS = 1
c.k1.STATUS=0
c.k2.STATUS=0
c.h1.STATUS = 1
c.h2.STATUS = 1
#c.h3.STATUS = 0
#c.h4.STATUS = 0

# FSTATUS = 0, no measurement
# FSTATUS = 1, measurement used to update model
c.v1.FSTATUS = 0
c.v2.FSTATUS = 0
c.k1.FSTATUS=1
c.k2.FSTATUS=1
c.h1.FSTATUS = 1
c.h2.FSTATUS = 1
c.h3.FSTATUS = 1
c.h4.FSTATUS = 1


sp=5
c.h1.SP=sp
c.h2.SP=sp

p1 = GEKKO(remote=False)
p1.time = [0,0.5]


#Parameters

p1.h1=p1.CV(lb=0,ub=25)
p1.h2=p1.CV(lb=0,ub=25)
p1.h3=p1.CV(lb=0,ub=25)
p1.h4=p1.CV(lb=0,ub=25)

p1.v1=p1.MV(value=3.15,lb=0.1,ub=8)
p1.v2=p1.MV(value=3.15,lb=0.1,ub=8)

p1.k1=p1.Param(lb=0,ub=10,value=3.14)
p1.k2=p1.Param(lb=0,ub=10,value=3.14)

#Equations
p1.Equation(A1*p1.h1.dt()==a3*((2*g*p1.h3)**0.5)-a1*((2*g*p1.h1)**0.5)+g1*p1.k1*p1.v1)
p1.Equation(A2*p1.h2.dt()==a4*((2*g*p1.h4)**0.5)-a2*((2*g*p1.h2)**0.5)+g2*p1.k2*p1.v2)
p1.Equation(A3*p1.h3.dt()==-a3*((2*g*p1.h3)**0.5)+(1-g2)*p1.k2*p1.v2)
p1.Equation(A4*p1.h4.dt()==-a4*((2*g*p1.h4)**0.5)+(1-g1)*p1.k1*p1.v1)

#options
p1.options.IMODE = 4

#%% problem configuration
# number of cycles
cycles = 480
# noise level


#%% run process, estimator and control for cycles
h1_meas = np.empty(cycles)
h2_meas =np.empty(cycles)
h3_meas =np.empty(cycles)
h4_meas=np.empty(cycles)
h1_est = np.empty(cycles)
h2_est = np.empty(cycles)
h3_est = np.empty(cycles)
h4_est = np.empty(cycles)
h1_plant=np.empty(cycles)
h2_plant=np.empty(cycles)
h3_plant=np.empty(cycles)
h4_plant=np.empty(cycles)
h1_measured=np.empty(cycles)
h2_measured=np.empty(cycles)
h3_measured=np.empty(cycles)
h4_measured=np.empty(cycles)
v1_est = np.empty(cycles)
v2_est = np.empty(cycles)
k1_est = np.empty(cycles)
k2_est = np.empty(cycles)
u_cont_k1 = np.empty(cycles)
u_cont_k2 = np.empty(cycles)
sp_store = np.empty(cycles)
sum_est=np.empty(cycles)
sum_model=np.empty(cycles)

# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()

p.MAX_ITER=20
c.MAX_ITER=20
m.MAX_ITER=20
p1.MAX_ITER=20

for i in range(cycles):
    print(i)
    # set point changes
    if i==cycles/4:
        sp = 12
    elif i==2*cycles/4:
        sp = 7
    elif i==3*cycles/4:
        sp = 5
    
    
    sp_store[i] = sp
    c.h1.SP=sp
    c.h2.SP=sp
   
    c.k1.MEAS = m.k1.NEWVAL
    c.k2.MEAS = m.k2.NEWVAL

    if p.options.SOLVESTATUS == 1:
       # print("going:",i)
        c.h1.MEAS = p.h1.MODEL
        c.h2.MEAS = p.h2.MODEL
        c.h3.MEAS = p.h3.MODEL
        c.h4.MEAS = p.h4.MODEL
   
    print(i,'Plant Model:',p.h1.MODEL,p.h2.MODEL,p.h3.MODEL,p.h4.MODEL)
 
    c.solve(disp=False,debug=0)
        
    #print("NEWVAL:",i,c.u,c.u.NEWVAL)
    u_cont_k1[i] = c.v1.NEWVAL
    u_cont_k2[i] = c.v2.NEWVAL
    
    #print("Horizon:",i,c.h1[0:],c.h2[0:])
    #print("Move:",i,c.v1.NEWVAL,c.v2.NEWVAL)

    ## process simulator
    #load control move
    p.v1.MEAS = u_cont_k1[i]
    p.v2.MEAS = u_cont_k2[i]

    #simulate
 
    p.solve(disp=False,debug=0)
    
    #plant model
    p1.k1=3.14
    p1.k2=3.14
    p1.v1.MEAS = u_cont_k1[i]
    p1.v2.MEAS = u_cont_k2[i]
    p1.solve(disp=False,debug=0)
    
    h1_plant[i]=p1.h1.MODEL
    h2_plant[i]=p1.h2.MODEL
    h3_plant[i]=p1.h3.MODEL
    h4_plant[i]=p1.h4.MODEL
    
    h1_measured[i]=p1.h1.MODEL+(random()*2)*noise
    h2_measured[i]=p1.h2.MODEL+(random()*2)*noise
    h3_measured[i]=p1.h3.MODEL+(random()*2)*noise
    h4_measured[i]=p1.h4.MODEL+(random()*2)*noise
    
    
    #print("Model process output:",i,p.h1.MODEL,p.h2.MODEL,p.h3.MODEL,p.h4.MODEL)
    #load output with white noise
    h1_meas[i] = p.h1.MODEL+(random()-0.5)*noise
    h2_meas[i] = p.h2.MODEL+(random()-0.5)*noise
    h3_meas[i] = p.h3.MODEL+(random()-0.5)*noise
    h4_meas[i] = p.h4.MODEL+(random()-0.5)*noise
    
    #Only MPC
    

    ## estimator
    #load input and measured output
    m.v1.MEAS = u_cont_k1[i]
    m.v2.MEAS = u_cont_k2[i]
    #m.h1.MEAS = h1_meas[i]+(random()*2)*noise
    #m.h2.MEAS = h2_meas[i]+(random()*2)*noise
    #m.h3.MEAS = h3_meas[i]+(random()*2)*noise
    #m.h4.MEAS = h4_meas[i]+(random()*2)*noise
    m.h1.MEAS = h1_meas[i]
    m.h2.MEAS = h2_meas[i]
    m.h3.MEAS = h3_meas[i]
    m.h4.MEAS = h4_meas[i]
    #m.COLDSTART=2
    #optimize parameters
    
       
    m.solve(disp=False,debug=0)
    #store results
    if i>=process:
        h1_est[i] = m.h1.MODEL
        h2_est[i] = m.h2.MODEL
        h3_est[i] = m.h3.MODEL
        h4_est[i] = m.h4.MODEL
        v1_est[i] = m.v1.NEWVAL
        v2_est[i] = m.v2.NEWVAL
        k1_est[i]= m.k1.NEWVAL
        k2_est[i] = m.k2.NEWVAL
        
    print("Estimated h:",i,h1_est[i],h2_est[i],h3_est[i],h4_est[i])
    print("Estimated k:",i,k1_est[i],k2_est[i],p.k1[0],p.k2[0])
    print("Estimated v:",i,v1_est[i],v2_est[i])
    print("dh1/dt:",(a3*((2*g*h3_est[i])**0.5)-(a1*((2*g*h3_est[i])**0.5))+(g1*k1_est[i]*v1_est[i]))/A3)
    print("dh2/dt:",(a4*((2*g*h4_est[i])**0.5)-(a2*((2*g*h2_est[i])**0.5))+(g2*k2_est[i]*v2_est[i]))/A2)
    print("dh3/dt:",(-a3*((2*g*h3_est[i])**0.5)+((1-g2)*k2_est[i]*v2_est[i]))/A3)
    print("dh4/dt:",(-a4*((2*g*h4_est[i])**0.5)+((1-g1)*k1_est[i]*v1_est[i]))/A4)
    
    if i%1==0:
        plt.clf()
        plt.subplot(4,1,1)
        #plt.plot(h1_meas[0:i])
        #plt.plot(h2_meas[0:i])
        #plt.plot(h3_meas[0:i])
        #plt.plot(h4_meas[0:i])
        plt.plot(h1_est[0:i])
        plt.plot(h2_est[0:i])
        plt.plot(sp_store[0:i])
        
        
        plt.subplot(4,1,2)
        
        plt.plot(h3_est[0:i])
        plt.plot(h4_est[0:i])
        
        #plt.legend(('h1_pred','h2_pred','h3_pred','h4_pred'))


        plt.subplot(4,1,3)
        plt.plot(k1_est[0:i])
        plt.plot(k2_est[0:i])


        plt.subplot(4,1,4)
        plt.plot(v1_est[0:i])
        plt.plot(v2_est[0:i])
        
        plt.draw()
        plt.pause(0.05)
end=time.time()
print("total time:",end-start)

I feel there is some issue with my MHE+MPC code. However, I am not able to realize the mistake?

jupyter gekko python-3.x
1 Answers

Nice application. I needed a few imports to make the script work. These may be loaded automatically for you.

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

successful solution

The script solves successfully if a lower bound is included on all the level variables (1e-6). There is a problem when the level goes below zero or is at zero when using m.sqrt(). This small adjustment helps it solve successfully so it doesn't get into a region where it is undefined. Gekko solvers can't deal with imaginary numbers.

Although the solution is successful, it appears that the control performance oscillates. There may need to be some tuning of the application.

user_23669410
Related