GEKKO python, how to run a calculation until I get a desired value of a variable

user_19518038 0
Viewed 59

Does anyone know how to run a dynamic simulation in GEKKO Python (IMODE 7) until a certain variable takes on the desired value? For example, I have a differential equation and I need to do the calculation until the pressure is 1 bar. I understand how to make a calculation for 100 seconds, 200 seconds, etc., but I do not understand how to let it run until I have a variable of a specific desired value.

The code that I have is kind of large, so in principle this code that I want to modify:

I want to have it without a while statement because every while statement creates a new model that takes time. And I have much bigger code and this time makes a negative impact.

#Start calculation 
pressure_final = 0 
t_guess = 10
while pressure_final < 1:
    m = GEKKO(remote = False)  # Create GEKKO model           
    k = 0.1
    P = m.Var(0.1)
    #Time discretization    
    tf = t_guess
    nt = int(tf/1) + 1
    m.time = np.linspace(0,tf,nt)

    #Left boundary pressure
    m.Equation(P.dt() == k*P)

    print("Start Pressurization")
    m.options.IMODE = 7
    m.solve(disp = False)
    print("Finished pressurization")
    t_guess = t_guess + 1
    pressure_final = P[-1]

I will be very grateful for the answer!

gekko python
1 Answers

There are two ways to accomplish this.

The first method is with step-by-step simulation with IMODE=4. There is no need to recreate the model each cycle. A solve command time shifts the final conditions to the initial conditions automatically because m.options.TIME_SHIFT=1. The pressure may exceed the final target value using this method.

Simulation result 1

import numpy as np
from gekko import GEKKO

m = GEKKO(remote = False)  # Create GEKKO model once          
k = 0.1
P = m.Var(0.1)
m.time = [0,1]
m.Equation(P.dt() == k*P)
m.options.IMODE = 4
m.options.NODES = 3

t=[0]; Pr=[0.1] # store time / Pressure
for i in range(100):
    print(i)
    m.solve(disp = False)
    t.append(i+1)
    Pr.append(P[-1])
    if Pr[-1]>=1:
        break

import matplotlib.pyplot as plt
plt.plot(t,Pr,'b-o')
plt.plot([0,t[-1]],[1,1],'r--')
plt.ylabel('Pressure'), plt.xlabel('Time')
plt.savefig('pressure.png',dpi=300)
plt.show()

The second method is to make the final time an optimization variable an FV with tf.STATUS=1.

tf = m.FV(1,lb=1,ub=100)
tf.STATUS=1

Divide all differential terms by the tf value and use m.time between 0 and 1.

m.time = np.linspace(0,1,20)
m.Equation(P.dt()/tf == k*P)

Use IMODE=6 to minimize the deviation of the final pressure from 1 with

final = np.zeros(20); final[-1]=1
f = m.Param(final)
m.Minimize(f*(P-1)**2)

This gives the exact time of Final Time: 23.041309678 when the pressure limit is reached.

Exact solution

import numpy as np
from gekko import GEKKO

m = GEKKO(remote = False)  # Create GEKKO model once          
k = 0.1
P = m.Var(0.1)
m.time = np.linspace(0,1,20)
tf = m.FV(1,lb=1,ub=100); tf.STATUS=1
m.Equation(P.dt()/tf == k*P)
final = np.zeros(20); final[-1]=1
f = m.Param(final)
m.Minimize(f*(P-1)**2)
m.options.IMODE = 6
m.options.NODES = 3
m.solve(disp = False)

print('Final Time: ', tf[-1])

import matplotlib.pyplot as plt
plt.plot(m.time*tf[-1],P,'b-o')
plt.plot([0,tf[-1]],[1,1],'r--')
plt.ylabel('Pressure'), plt.xlabel('Time')
plt.savefig('pressure.png',dpi=300)
plt.show()
user_23669410
Related