Optimal control with free termination time (gekko)

user_16549177 0
Viewed 108

I am making a numerical problem to show as an example and am trying to find an optimal control using gekko for the following problem:

minimize the integral of a*x(t) from 0 to T, where T is the first time x(t) is 0, i.e., it is a random time. The constraints are such that x(t) follows some dynamic f(x(t),u(t)), x(t) >= 0, and u(t) is between 0 and 1.

I followed the tutorials on GEKKO website and youtube for fixed final time, but I could not find any information on a random final time. The following is the current code I have, but how would I be able to move from a fixed final time to a random final time? Any help would be appreciated! Thanks!

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO

# Initial conditions
xhh0 = 3; xhi0 = 0; 
xvh0 = 30; xvi0 = 0;
hin0 = 0; vin0 = 0;
tt0 = 0

# Parameters
a1 = 0.1; a2 = 0.1;
b1 = 0.01; b2 = 0.5;
delta1 = 0.1; delta2 = 0.5;
rho1 = 0.3; rho2 = 0.01
mu = 1

# Gekko
m = GEKKO()

# Control variable
u = m.MV(0.5, lb = 0, ub = 1)

# Final time <------------------------ currently a fixed final time
T = 10

# Initialize
xhh, xhi, xvh, xvi, Ah, Av = m.Array(m.Var, 6)
xhh.value = xhh0; xhi.value = xhi0;
xvh.value = xvh0; xvi.value = xvi0;
Ah.value = hin0; Av.value = vin0; 

# System dynamics
m.Equations([xhh.dt() == -a1*xhh - mu*u - b1*xhi*xhh,\
             xhi.dt() == a1*xhh + b1*xhi*xhh - delta1*xhi - rho1*xhi,\
             xvh.dt() == -a2*xvh - mu*(1-u) - b2*xvi*xvh,\
             xvi.dt() == a2*xvh + b2*xvi*xvh - delta2*xvi - rho2*xvi,\
             Ah.dt() == a1*xhh,\
             Av.dt() == a2*xvh])

# Time space
t = np.linspace(0, T, 101) 

m.time = t

# initialize with simulation
m.options.IMODE = 7
m.options.NODES = 3
m.solve(disp = False)

# optimization
m.options.IMODE = 6
xhh.LOWER = 0; xhi.LOWER = 0; xvh.LOWER = 0; xvi.LOWER = 0
u.STATUS = 1
m.options.SOLVER = 3
xhh.value = xhh.value.value
xhi.value = xhi.value.value
xvh.value = xvh.value.value
xvi.value = xvi.value.value
Ah.value = Ah.value.value
Av.value = Av.value.value

# Objective function
m.Minimize(Ah + Av)
m.solve()
gekko
1 Answers

The final time is adjustable with T = m.FV() and T.STATUS=1 when each differential is divided by T. This scales the problem to any arbitrary final time when t = np.linspace(0,1).

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from gekko import GEKKO

# Initial conditions
xhh0 = 3; xhi0 = 0; 
xvh0 = 30; xvi0 = 0;
hin0 = 0; vin0 = 0;
tt0 = 0

# Parameters
a1 = 0.1; a2 = 0.1;
b1 = 0.01; b2 = 0.5;
delta1 = 0.1; delta2 = 0.5;
rho1 = 0.3; rho2 = 0.01; mu = 1

# Gekko
m = GEKKO()

# Control variable
u = m.MV(0.5, lb = 0, ub = 1)

# Final time
T = m.FV(10,lb=1e-2,ub=100); T.STATUS = 1

# Initialize
xhh, xhi, xvh, xvi, Ah, Av = m.Array(m.Var, 6)
xhh.value = xhh0; xhi.value = xhi0;
xvh.value = xvh0; xvi.value = xvi0;
Ah.value = hin0; Av.value = vin0; 
xhh.LOWER = 0; xhi.LOWER = 0; xvh.LOWER = 0; xvi.LOWER = 0
u.STATUS = 1

# System dynamics
m.Equations([xhh.dt()/T == -a1*xhh - mu*u - b1*xhi*xhh,\
             xhi.dt()/T == a1*xhh + b1*xhi*xhh - delta1*xhi - rho1*xhi,\
             xvh.dt()/T == -a2*xvh - mu*(1-u) - b2*xvi*xvh,\
             xvi.dt()/T == a2*xvh + b2*xvi*xvh - delta2*xvi - rho2*xvi,\
             Ah.dt()/T == a1*xhh,\
             Av.dt()/T == a2*xvh])

# Time space
t = np.linspace(0, 1, 101) 
m.time = t

# optimization
m.options.IMODE = 6
m.options.SOLVER = 3

# Objective function
m.Minimize(Ah + Av)
m.solve()

print('Final time: ', T.value[0])

There may be a missing constraint or some other information because the optimal final time always goes to the lower bound. The Jennings problem is a related example with variable final time.

user_23669410
Related