GEKKO Optimal control. How to force values to be above certain thresholds?

user_15692976 0
Viewed 98

So I am simulating plane flight. The plane flies certain distance (path variable) and then the simulation stops. The solver is trying to minimize the fuel consumption (m.Maximize(mass*tf*final), by maximizing the mass value. The solver controls the accelerator pedal position (Tcontr). I need it to try and keep Velocity value V above certain threshold(for example >=100). What kind of equation I need to add for that to happen?

My code:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
import math
#Gekko model
m = GEKKO(remote=False)

#Time points
nt = 11
tm =  np.linspace(0,100,nt)
m.time = tm

# Variables
Ro = m.Const(value=1.1)#
g = m.Const(value=9.80665)
pressure = m.Var(value=101325)#
T = m.Var()
T0 = m.Const(value=273)
S = m.Const(value=122.6)
Cd = m.Const(value=0.1)
FuelFlow = m.Var()
D = m.Var()
Thrmax = m.Const(value=200000)
Thr = m.Var()
V = m.Var(value=100,lb=0,ub=240)
Vmin = m.Var(value=100)
nu = m.Var(value=0)
nuu = nu.value
x = m.Var(value=0,lb=0)
h = m.Var(value=1000)
mass = m.Var(value=60000)
path = m.Const(value=5000000) #intended distance  travelled
L = m.Var()

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

m.options.MAX_ITER=50000 # iteration number

#Fixed Variable
tf = m.FV(value=1,lb=0.0001,ub=1000.0)#
tf.STATUS = 1

# Parameters
Tcontr = m.MV(value=0.2,lb=0.2,ub=1)
Tcontr.STATUS = 1
Tcontr.DCOST = 0

# Equations
m.Equation(x.dt()==tf*(V*(math.cos(nuu.value))))#
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(V.dt()==tf*((Thr-D)/mass))#
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60000)))#
m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(x*final<=path)

m.Equation(T==T0-(0.0065*h))    

# Objective Function
m.Minimize(final*(x-path)**2)
m.Maximize(mass*tf*final)
m.options.IMODE = 6
m.options.NODES = 2 # it was 3 before
m.options.MV_TYPE = 1
m.options.SOLVER = 3
#m.open_folder() # to search for infeasibilities
m.solve()

tm = tm * tf.value[0]

print('Final Time: ' + str((tf.value[-1])*100))
print('Final Speed: ' + str(V.value[-1]))
print('Final X: ' + str(x.value[-1]))
print('Final Mass: ' + str(mass.value[-1]))

fig, axs = plt.subplots(5)
fig.suptitle('Results')
axs[0].plot(tm,Tcontr,'r-',LineWidth=2,label=r'$Tcontr$')
axs[0].legend(loc='best')
axs[1].plot(tm,V.value,'b-',LineWidth=2,label=r'$V$')
axs[1].legend(loc='best')
axs[2].plot(tm,x.value,'r--',LineWidth=2,label=r'$x$')
axs[2].legend(loc='best')
axs[3].plot(tm,D.value,'g-',LineWidth=2,label=r'$D$')
axs[3].legend(loc='best')
axs[4].plot(tm,mass.value,'g:',LineWidth=2,label=r'$mass$')
axs[4].legend(loc='best')
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()
gekko python
1 Answers

There are two ways to implement that in Gekko.

Hard Constraint

The first way is to set the lower bound for the Velocity to 100 with:

V = m.Var(value=100,lb=100,ub=240)

This is a hard constraint - the solver does not return a successful solution unless this is satisfied. This can lead to an infeasible solution, however.

Soft Constraint

If it is desirable, but not required for the Velocity to stay above 100 or there are problems with infeasible solutions then a soft constraint may be a better option. Set the following to implement a dead-band between 100 and 240:

V = m.CV(value=100,lb=0,ub=240)
V.SPHI = 240; V.SPLO = 100
V.WSPHI = 1; V.WSPLO = 10  # penalty for violating bound
V.STATUS = 1; V.TR_INIT = 0

More information on this type of variable is given here in the Machine Learning and Dynamic Optimization Course with a detailed description of the l1-norm objective.

user_23669410
Related