1D Plane flight optimal control problem in gekko

user_15692976 0
Viewed 97

I need to solve 1D plane flight optimal control problem. I have a plane that is 1000m high. I need it to travel 1000m forward along x-axis while minimizing fuel consumption. And when it achieves x=1000 I need program to stop.

My objective function looks like this: min J => -mass(tf). (By maximizing mass, the consumed fuel gets minimized). Optimizer controls the accelerator pedal.

The problem is subject to the following constraints: dx/dt = V; dV/dt = (Throttle - Drag)/mass; dm/dt = -Throttle * fuelflow.

I have developed a Gekko / Python script as follows. However, the optimization script does not achieve a solution, raising:

Exception: @error: Solution Not Found

I have tried to change stuff here and there, but it did not work and I'm not sure what the problem is. Here is 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 = 101
m.time = np.linspace(0,100,nt)

# Parameters
Tcontr = m.MV(value=0.5,lb=0.3,ub=1) #throttle pedal position
Tcontr.STATUS = 1
Tcontr.DCOST = 0

# Variables
Ro = 1.1 #air density
S = 122.6
Cd = value=0.1
FuelFlow = m.Var(value=0.7)
D = m.Var() #drag
Thrmax = 200000 #maximum theoretical throttle
Thr = m.Var() #real throttle
V = m.Var() #Velocity
x = m.Var(value=0)
mass = m.Var(value=60000)

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

# Equations
m.Equation(x.dt()==V)
m.Equation(Thr==Tcontr*Thrmax) #Throttle
m.Equation(V.dt()==(Thr-D)/mass)
m.Equation(mass.dt()==-Thr*FuelFlow)
m.Equation(D==0.5*Ro*(V**2)*Cd*S) #Drag 
m.Equation(final*x==1000) # to stop when x==1000 is achieved

# Objective Function
m.Obj(-mass*final)
m.options.IMODE = 6
m.options.NODES = 2
m.options.MV_TYPE = 1
m.options.SOLVER = 3
m.open_folder() 
m.solve()
gekko python
1 Answers

It seems like the major problem was that the upper bound of your MV (Tcontr) is little bit too low to meet all the constraints.

And, I guess the "FuelFlow" is need to be a m.Const as I assumed it is a fuel flowrate per throttle.

Lastly, when you use the final for designating the value at the last time point, you need to make sure the associated Equation is also feasible for the entire time steps. I suggest to use this instead.

m.Equation(final*(x-1000)==0)

That way, you can achive your goal without posing any infeasible equation as your previous equation did. Your previous equation is infeasible at the most of the time points except the last time point with non-zero value. It's been resulting "0 (LHS) = 1000 (RHS)" before the final time point.

Please see the modified code below.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
# import math

#Gekko model
m = GEKKO(remote=True)

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

# Parameters
Tcontr = m.MV(value=0.5,lb=0.3,ub=10) #throttle pedal position
Tcontr.STATUS = 1
Tcontr.DCOST = 1e-2

# Variables
Ro = 1.1 #air density
S = 122.6
Cd = value=0.1
FuelFlow = m.Const(value=0.7)
D = m.Var() #drag
Thrmax = 200000 #maximum theoretical throttle
Thr = m.Var() #real throttle
V = m.Var() #Velocity
x = m.Var(value=0)
mass = m.Var(value=60000)

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

# Equations
m.Equation(x.dt()==V)
m.Equation(Thr==Tcontr*Thrmax) #Throttle
m.Equation(mass*V.dt()==(Thr-D))
m.Equation(mass.dt()==-Thr*FuelFlow)
m.Equation(D==0.5*Ro*(V**2)*Cd*S) #Drag 
m.Equation(final*(x-1000)==0) # to stop when x==1000 is achieved

# Objective Function
m.Obj(-mass*final)
m.options.IMODE = 6
m.options.NODES = 3
m.options.MV_TYPE = 1
m.options.SOLVER = 3
# m.open_folder() 
m.solve()
user_113419940
Related