pythonoptimizationcontrollersimulationgekko
Controlled variable does not follow SP in gekko
I am trying to control a level system (of two tanks) by using gekko. In the trial related in this question, just one of the level is under control. However, my CV does not seem to follow the SP.
Code and data below:
# Influent flows
flow1 = np.array([1.25282567, 1.18243383, 1.25445315, 1.1755483 , 1.17953059,
1.164077 , 1.18473485, 1.18064562, 1.14043986, 1.23820012,
1.24054883, 1.20570221, 1.21603648, 1.20606003, 1.15339559,
1.18137919, 1.20935554, 1.20335999, 1.21147486, 1.23357391,
1.21825142, 1.19723155, 1.21343981, 1.11376114, 1.19094124,
1.20130132, 1.21304258, 1.18635852, 1.26465022, 1.18148956,
1.20242284, 1.23000755, 1.26954079, 1.16466064, 1.23057836,
1.23027471]) # goes only to tank 2
flow2 = np.array([2.59677412, 2.58811843, 2.544673 , 2.55571262, 2.57453146,
2.59810488, 2.56841369, 2.45821758, 2.4281892 , 2.43671096,
2.4474268 , 2.46243303, 2.43683383, 2.39242551, 2.4513853 ,
2.45447782, 2.42802739, 2.4221822 , 2.42537811, 2.44381434,
2.44872821, 2.44426348, 2.4158827 , 2.42840086, 2.47415572,
2.50235597, 2.47706445, 2.44440066, 2.44883323, 2.46050179,
2.47231121, 2.49074611, 2.48412978, 2.47831733, 2.47747716,
2.48604111]) # goes to tank 1 and 2
flow3 = np.array([1.57011312, 1.50738452, 1.52807646, 1.57724425, 1.56187045,
1.64383615, 1.57872569, 1.57535121, 1.5540041 , 1.51887366,
1.58146572, 1.57476749, 1.60523332, 1.5716739 , 1.56179746,
1.53903109, 1.58319486, 1.56445033, 1.58821278, 1.61287479,
1.57210261, 1.6220108 , 1.5820968 , 1.55307617, 1.52589127,
1.57522106, 1.55631243, 1.55931072, 1.48830357, 1.57045034,
1.50157251, 1.55527039, 1.56454099, 1.56372395, 1.51546389,
1.552174 ]) #goes only to tank 1
# Effluent
ts_out = np.array([6.53681895, 6.211645 , 6.12570737, 6.36272278, 6.47968763,
6.88415985, 6.53169619, 6.44752413, 6.76828499, 7.24453761,
7.3429685 , 7.51939125, 7.13905894, 7.45646591, 7.53236135,
7.64722176, 7.36411756, 7.18427607, 7.24093132, 7.13198001,
6.98640455, 6.64083468, 6.36015803, 6.41666819, 6.58571269,
6.52822367, 6.79654961, 6.46079076, 6.4651886 , 6.68471489,
7.02170868, 7.23842901, 7.3350867 , 7.32969907, 7.27846311,
7.27119013]) # sum of both effluent tanks
# Tank t1
# Constants
t1_area = 116.89
t1_height = 15.2
#initial state
t1_h0 = 13.9
# Tank t2
# Constants
t2_area = 226.98
t2_height = 15
#initial state
t2_h0 = 13.38
#Define models
def tanks(x, t, qd_in1, q_out, qt1_v, qt2_v_out,
qt2_in1, A2,
qt1_in2, A1):
qt1_in1 = qd_in1*qt1_v
qt2_in2 = qd_in1*(1-qt1_v)
qt2_out = qt2_v_out*q_out
qt1_out = (1 - qt2_v_out)*q_out
dh2dt = (qt2_in1 + qt2_in2 - qt2_out)/A2
dh1dt = (qt1_in1 + qt1_in2 - qt1_out)/A1
return dh1dt, dh2dt
# Gekko model
m = GEKKO(remote=False)
m.time = np.arange(0,12,2)
#Flow 2 goes to both tanks
qd_in1 = m.SV(value = flow2[0], name ='flow2') #flow2 to t1+t2
# Tank t1
A1 = m.Const(value = t1_area, name = 'A1') # m2
h1 = m.CV(t1_h0, name = 'h1', lb=0.6*t1_height, ub=0.95*t1_height) #m, controlled variable
qt1v = 0.8 # Initial distribuition to tank 1
qt1_v = m.MV(qt1v, lb=0, ub=1, name = 'flow2_t1_v') # manipulate flow of flow2 to t1
qt1_in1 = m.Intermediate(qd_in1*qt1_v, name = 'flow2_t1') # flow2 to t1
qt1_in2 = m.SV(flow3[0],name ='flow3')
# Tank t2
A2 = m.Const(value = t2_area, name = 'A2') # m2,Cross section Area
h2 = m.CV(t2_h0, name = 'h2', lb=0.6*t2_height, ub=0.95*t2_height) #m , the controlled variable
qt2_in1 = m.SV(flow1[0], name = 'flow1')
qt2_in2 = m.Intermediate(qd_in1*(1-qt1_v), name = 'flow2_t2') # flow2 to t2
#Out
q_out = m.SV(ts_out[0], name = 'q_out')
qt2_vout = 0.32 # 32% of outflow belongs to t2 (later will be optimized)
qt2_v_out = m.SV(qt2_vout, name = 't2_out_v')
qt2_out = m.Intermediate(qt2_v_out*q_out, name = 't2_out') #t2 to out
qt1_out = m.Intermediate((1-qt2_v_out)*q_out, name = 't1_out') #t1 to out
m.Equation(h2.dt() == ((qt2_in1 + qt2_in2 - qt2_out)/A2))
m.Equation(h1.dt() == ((qt1_in1 + qt1_in2 - qt1_out)/A1))
m.options.IMODE = 6
m.options.NODES = 3
m.CV_TYPE=1
#CV tuning t1
h1.STATUS = 1 #add the SP to the objective
h1.FSTATUS = 1
h1.SPHI = 0.90*t1_height
h1.SPLO = 0.80*t1_height
h1.TR_INIT = 1
h1.WSPHI = 800
h1.WSPLO = 800
# h1.TAU = 10 #time constant of setpoint trajectory
#MV tuning for t1
qt1_v.STATUS = 1 #allow optimizer to change
qt1_v.DCOST = 0.0000001 #smooth out movement
qt1_v.DMAXLO = - 0.3
qt1_v.DMAXHI = 0.3
Then I'm running both models in a for loop:
for i in range(len(t)-1):
# Simulate
inputs = (flow2[i], ts_out[i], qt1_v_store[i], qt2_vout_store[i],
flow1[i], t2_area,
flow3[i], t1_area)
ts = [t[i],t[i+1]]
y = odeint(tanks,x,ts,args=inputs)
# Record measurements
h1.MEAS = y[-1][0]
h2.MEAS = y[-1][1]
# Store results
t1_h_ode[i+1] = y[-1][0]
t2_h_ode[i+1] = y[-1][1]
# Adjust initial condition for next loop
x[0] = y[-1][0]
x[1] = y[-1][1]
m.solve(disp=False)
t1_h_store[i+1] = h1.model
t2_h_store[i+1] = h2.model
qd_in1.MEAS = flow2[i+1]
qt1_in2.MEAS = flow3[i+1]
qt2_in1.MEAS = flow1[i+1]
q_out.MEAS = ts_out[i+1]
qt1_v_store[i+1] = qt1_v.NEWVAL
The code runs but it doesn't matter what I try, the CV h1 won't follow the SP. It just decreases... I tried tuning the parameters but with no luck. Is there something wrong in my definition?
Solution
The m.SV() definitions should be m.FV() as input parameters to the model. Each m.Var(), m.SV(), and m.CV() needs a corresponding equation to define the value of that calculated value. Also set the qt1_v.FSTATUS=0 to not reset the value to the initially specified conditions each cycle.
Below are the results with the modified script. The controller adjusts the two MV values to stay within the upper (SPHI) and lower (SPLO) limits that are 85-90% of the tank 1 height.
import matplotlib.pyplot as plt
from gekko import GEKKO
import numpy as np
from scipy.integrate import odeint
# Influent flows
flow1 = np.array([1.25282567, 1.18243383, 1.25445315, 1.1755483 , 1.17953059,
1.164077 , 1.18473485, 1.18064562, 1.14043986, 1.23820012,
1.24054883, 1.20570221, 1.21603648, 1.20606003, 1.15339559,
1.18137919, 1.20935554, 1.20335999, 1.21147486, 1.23357391,
1.21825142, 1.19723155, 1.21343981, 1.11376114, 1.19094124,
1.20130132, 1.21304258, 1.18635852, 1.26465022, 1.18148956,
1.20242284, 1.23000755, 1.26954079, 1.16466064, 1.23057836,
1.23027471]) # goes only to tank 2
flow2 = np.array([2.59677412, 2.58811843, 2.544673 , 2.55571262, 2.57453146,
2.59810488, 2.56841369, 2.45821758, 2.4281892 , 2.43671096,
2.4474268 , 2.46243303, 2.43683383, 2.39242551, 2.4513853 ,
2.45447782, 2.42802739, 2.4221822 , 2.42537811, 2.44381434,
2.44872821, 2.44426348, 2.4158827 , 2.42840086, 2.47415572,
2.50235597, 2.47706445, 2.44440066, 2.44883323, 2.46050179,
2.47231121, 2.49074611, 2.48412978, 2.47831733, 2.47747716,
2.48604111]) # goes to tank 1 and 2
flow3 = np.array([1.57011312, 1.50738452, 1.52807646, 1.57724425, 1.56187045,
1.64383615, 1.57872569, 1.57535121, 1.5540041 , 1.51887366,
1.58146572, 1.57476749, 1.60523332, 1.5716739 , 1.56179746,
1.53903109, 1.58319486, 1.56445033, 1.58821278, 1.61287479,
1.57210261, 1.6220108 , 1.5820968 , 1.55307617, 1.52589127,
1.57522106, 1.55631243, 1.55931072, 1.48830357, 1.57045034,
1.50157251, 1.55527039, 1.56454099, 1.56372395, 1.51546389,
1.552174 ]) #goes only to tank 1
# Effluent
ts_out = np.array([6.53681895, 6.211645 , 6.12570737, 6.36272278, 6.47968763,
6.88415985, 6.53169619, 6.44752413, 6.76828499, 7.24453761,
7.3429685 , 7.51939125, 7.13905894, 7.45646591, 7.53236135,
7.64722176, 7.36411756, 7.18427607, 7.24093132, 7.13198001,
6.98640455, 6.64083468, 6.36015803, 6.41666819, 6.58571269,
6.52822367, 6.79654961, 6.46079076, 6.4651886 , 6.68471489,
7.02170868, 7.23842901, 7.3350867 , 7.32969907, 7.27846311,
7.27119013]) # sum of both effluent tanks
# Tank t1
# Constants
t1_area = 116.89
t1_height = 15.2
#initial state
t1_h0 = 13.9
# Tank t2
# Constants
t2_area = 226.98
t2_height = 15
#initial state
t2_h0 = 13.38
#Define models
def tanks(x, t, qd_in1, q_out, qt1_v, qt2_v_out,
qt2_in1, A2,
qt1_in2, A1):
qt1_in1 = qd_in1*qt1_v
qt2_in2 = qd_in1*(1-qt1_v)
qt2_out = qt2_v_out*q_out
qt1_out = (1 - qt2_v_out)*q_out
dh2dt = (qt2_in1 + qt2_in2 - qt2_out)/A2
dh1dt = (qt1_in1 + qt1_in2 - qt1_out)/A1
return dh1dt, dh2dt
# Gekko model
m = GEKKO(remote=False)
m.time = np.arange(0,12,2)
# Flow 2 goes to both tanks
qd_in1 = m.FV(value = flow2[0], name ='flow2') #flow2 to t1+t2
qd_in1.FSTATUS = 1
# Tank t1
A1 = m.Const(value = t1_area, name = 'A1') # m2
# Do not add hard constraints - can lead to infeasibility
# Removed: lb=0.6*t1_height, ub=0.95*t1_height
h1 = m.CV(t1_h0, name = 'h1') # controlled variable
qt1v = 0.8 # Initial distribuition to tank 1
qt1_v = m.MV(qt1v, lb=0, ub=1, name = 'flow2_t1_v') # manipulate flow of flow2 to t1
qt1_in1 = m.Intermediate(qd_in1*qt1_v, name = 'flow2_t1') # flow2 to t1
qt1_in2 = m.FV(flow3[0],name ='flow3')
qt1_in2.FSTATUS = 1
# Tank t2
A2 = m.Const(value = t2_area, name = 'A2') # m2,Cross section Area
h2 = m.CV(t2_h0, name = 'h2') #, lb=0.6*t2_height, ub=0.95*t2_height) #m , the controlled variable
qt2_in1 = m.FV(flow1[0], name = 'flow1')
qt1_in1.FSTATUS = 1
qt2_in2 = m.Intermediate(qd_in1*(1-qt1_v), name = 'flow2_t2') # flow2 to t2
#Out
q_out = m.FV(ts_out[0], name = 'q_out')
q_out.FSTATUS = 1
qt2_vout = 0.32 # 32% of outflow belongs to t2 (later will be optimized)
qt2_v_out = m.MV(qt2_vout, lb=0.05, ub=0.5, name = 't2_out_v')
qt2_v_out.STATUS = 1
qt2_v_out.FSTATUS = 0
qt2_v_out.DCOST = 0.0000001 #smooth out movement
qt2_v_out.DMAXLO = -0.3
qt2_v_out.DMAXHI = 0.3
qt2_out = m.Intermediate(qt2_v_out*q_out, name = 't2_out') #t2 to out
qt1_out = m.Intermediate((1-qt2_v_out)*q_out, name = 't1_out') #t1 to out
m.Equation(A2 * h2.dt() == qt2_in1 + qt2_in2 - qt2_out)
m.Equation(A1 * h1.dt() == qt1_in1 + qt1_in2 - qt1_out)
m.options.IMODE = 6
m.options.NODES = 3
m.CV_TYPE=1
#CV tuning t1
h1.STATUS = 1 #add the SP to the objective
h1.FSTATUS = 0
h1.SPHI = 0.90*t1_height
h1.SPLO = 0.85*t1_height
h1.TR_INIT = 0
h1.WSPHI = 800
h1.WSPLO = 800
# h1.TAU = 10 #time constant of setpoint trajectory
#MV tuning for t1
qt1_v.STATUS = 1 #allow optimizer to change
qt1_v.FSTATUS = 0
qt1_v.DCOST = 0.0000001 #smooth out movement
qt1_v.DMAXLO = -0.3
qt1_v.DMAXHI = 0.3
t = np.arange(0,72,2)
t1_h_ode = [t1_h0]
t2_h_ode = [t2_h0]
t1_h_store = [t1_h0]
t2_h_store = [t2_h0]
qt1_v_store = [qt1v]
qt2_vout_store = [qt2_vout]
x = [t1_h0,t2_h0]
for i in range(len(t)-1):
print(i)
# Simulate
inputs = (flow2[i], ts_out[i], qt1_v_store[i], qt2_vout_store[i],
flow1[i], t2_area,
flow3[i], t1_area)
ts = [t[i],t[i+1]]
y = odeint(tanks,x,ts,args=inputs)
# Record measurements (not needed, should be perfect match)
#h1.MEAS = y[-1][0]
#h2.MEAS = y[-1][1]
qd_in1.MEAS = flow2[i]
qt1_in2.MEAS = flow3[i]
qt2_in1.MEAS = flow1[i]
q_out.MEAS = ts_out[i]
# Store results
t1_h_ode.append(y[-1][0])
t2_h_ode.append(y[-1][1])
# Adjust initial condition for next loop
x[0] = y[-1][0]
x[1] = y[-1][1]
m.solve(disp=True)
t1_h_store.append(h1.PRED[1])
t2_h_store.append(h2.PRED[1])
qt1_v_store.append(qt1_v.NEWVAL)
qt2_vout_store.append(qt2_v_out.NEWVAL)
# Create a figure with three subplots (stacked vertically)
fig, ax = plt.subplots(3, 1, figsize=(8, 6), sharex=True)
# Top subplot: Tank 1 Level (h1)
ax[0].plot(t, t1_h_ode, 'ko-', label='Measured h1')
ax[0].plot(t, t1_h_store, 'b.-', label='Predicted h1')
ax[0].plot([t[0],t[-1]], [h1.SPHI,h1.SPHI], 'r--', label='SPHI')
ax[0].plot([t[0],t[-1]], [h1.SPLO,h1.SPLO], 'r--', label='SPLO')
ax[0].set_title('Tank 1 Level (h1)')
ax[0].set_ylabel('Level (m)')
ax[0].legend()
ax[0].grid(True)
# Middle subplot: Tank 2 Level (h2)
ax[1].plot(t, t2_h_ode, 'ko-', label='Measured h2')
ax[1].plot(t, t2_h_store, 'r.-', label='Predicted h2')
ax[1].set_title('Tank 2 Level (h2)')
ax[1].set_ylabel('Level (m)')
ax[1].legend()
ax[1].grid(True)
# Bottom subplot: Flow distribution
ax[2].plot(t, qt1_v_store, 'g.-', label='qt1_v')
ax[2].plot(t, qt2_vout_store, 'b--', label='qt2_vout')
ax[2].set_title('Flow Distribution')
ax[2].set_xlabel('Time (s)')
ax[2].set_ylabel('Fraction')
ax[2].legend()
ax[2].grid(True)
plt.tight_layout()
plt.savefig('results.png',dpi=300)
plt.show()
The plot and some of the missing information is added. A new m.MV() is defined for qt2_vout and this can be turned on or off with qt2_vout.STATUS. For future questions, please create a complete and minimal example of the problem that reproduces the issue.
- Unexpected list append
- Force matrix_world to be recalculated in Blender
- SQLAlchemy and empty columns
- ValueError: time data '24:00' does not match format '%H:%M'
- Convert RDD of LabeledPoint to DataFrame toDF() Error
- How to cancel trigonometric expressions in SymPy
- Get view used in Django tests
- Precompiled sasl python 3.9+ package for windows
- Regex: Substitute pattern in string multiple times without leftovers
- How to render raw html in the PyHTML library
- Why does my implementation of trilateration give wrong results?
- Django admin: how to sort by one of the custom list_display fields that has no database field
- TypeError: not all arguments converted during string formatting - psycopg2
- Is there a Python equivalent of the C# null-coalescing operator?
- Kraken API - Account balances request returning Invalid Nonce
- configparser without whitespace surrounding operator
- Pytorch tensor to numpy array
- Django: How to get a person whose birthday is today from a database?
- Performance impact of inheriting from many classes
- How can I do a line break (line continuation) in Python (split up a long line of source code)?
- Using pydantic to change int to string
- Breaking long method chains into multiple lines in Python
- What do ** (double star/asterisk) and * (star/asterisk) mean in a function call?
- How to install Pygame on Python 3.4?
- Rotating values in a list [Python]
- Launch default image viewer from pygtk program
- what's the inverse of the quantile function on a pandas Series?
- How can I install packages using pip according to the requirements.txt file from a local directory?
- Python generate all n-permutations of n lists
- FastAPI error when handling file together with form-data defined in a Pydantic model