1

PSA: I am very new to gekko, thus I might be missing something very obvious here.

I have been trying to find the solution to an optimal control problem, namely trajectory optimization of a regular vehicle, under certain speed constraints at certain distances along their trip. In order to do this, I tried using a pwl function based on the distance and speed constraint data and using v_max as a constraint to v. As a objective function, I use a Vehicle Specific Power (VSP) approximation.

The computation keeps going until the maximum no. of iterations is reached and cancels. Is there maybe a way to discretize the search space of this problem to make it solvable in acceptable time trading off computation time for accuracy?

  • goal_dist = The distance that needs to be covered
  • max_accel = maximum possible acceleration of the vehicle
  • max_decel = maximum possible deceleration of the vehicle
  • max_velocity = maximum possible velocity of the vehicle
  • min_velocity = minimum possible velocity of the vehicle
  • trip_time = No. of discrete data points (1s apart)
  • distances = array of length trip_time of discrete distance values based on GPS data points of the desired trip
  • speed_limits = array of length trip_time of discrete speed limits based on GPS data points of the desired trip
  • slope = array of length trip_time of discrete angle values
def optimal_trip(goal_dist, max_accel, max_decel, max_velocity, min_velocity, trip_time, distances ,speed_limits, slope):

    model = GEKKO(remote=True)
    model.time = [i for i in range(trip_time)]

    x = model.Var(value=0.0)
    v = model.Var(value=0.0, lb = min_velocity, ub = max_velocity)

    v_max = model.Var()
    slope_var = model.Var()
    
    a = model.MV(value=0, lb=max_decel ,ub=max_accel)
    a.STATUS = 1
    
    #define vehicle movement
    model.Equation(x.dt()==v)
    model.Equation(v.dt()==a)
    # path constraint
    model.Equation(x >= 0)

    #aggregated velocity constraint
    model.pwl(x, v_max, distances, speed_limits)
    model.Equation(v_max>=v)

#slope is modeled as a piecewise linear function
    model.pwl(x, slope_var, distances, slope)

    #End state constraints
    model.fix(x, pos=trip_time-1,val=goal_dist) # vehicle must arrive at destination
    model.fix(v, pos=trip_time-1,val=0) # vehicle must be fully stopped
    #VSPI Objective function
    obj = (v * (1.1 * a + 9.81 * slope_var + 0.132) +0.0003002*pow(v, 3))
    model.Obj(obj)
    # solve
    model.options.IMODE = 6
    model.options.REDUCE = 3
    model.solve(disp=True)
    return x.value, v.value, obj.value

Could someone shed some light onto this?

1 Answers1

0

Here is a version of the model with sample values that solves successfully:

from gekko import GEKKO
import numpy as np

min_velocity =  0
max_velocity = 10
max_decel    = -1
max_accel    =  1
distances    = np.linspace(0,20,21)
goal_dist    = 200

trip_time    = 100

# set up PWL functions
distances    = np.linspace(0,1000,10)
speed_limits = np.ones(10)*5
speed_limits[5:]=7
slope        = np.zeros(10)
slope[3:5]=1; slope[7:9]=-1

model = GEKKO(remote=True)
model.time = [i for i in range(trip_time)]

x = model.Var(value=0.0)
v = model.Var(value=0.0, lb = min_velocity, ub = max_velocity)

v_max = model.Var()
slope_var = model.Var()

a = model.MV(value=0, lb=max_decel ,ub=max_accel)
a.STATUS = 1

#define vehicle movement
model.Equation(x.dt()==v)
model.Equation(v.dt()==a)
# path constraint
model.Equation(x >= 0)

#aggregated velocity constraint
model.pwl(x, v_max, distances, speed_limits)
model.Equation(v_max>=v)

#slope is modeled as a piecewise linear function
model.pwl(x, slope_var, distances, slope)

#End state constraints
model.fix(x, pos=trip_time-1,val=goal_dist) # vehicle must arrive at destination
model.fix(v, pos=trip_time-1,val=0) # vehicle must be fully stopped
#VSPI Objective function
obj = (v * (1.1 * a + 9.81 * slope_var + 0.132) +0.0003002*pow(v, 3))
model.Obj(obj)
# solve
model.options.IMODE = 6
model.options.REDUCE = 3
model.solve(disp=True)

It may be that the values you are using cause an infeasible solution. Here are some suggestions to help the model solve more reliably:

  1. Use upper bounds on variables instead of general inequality constraints when possible.
# remove these lines
#model.Equation(x >= 0)
#x = model.Var(value=0.0)

# put lower bound on x
x = model.Var(value=0,lb=0)
  1. Use soft terminal constraints instead of hard terminal constraints.
#End state constraints
# vehicle must arrive at destination
#model.fix(x, pos=trip_time-1,val=goal_dist)
# vehicle must be fully stopped
#model.fix(v, pos=trip_time-1,val=0) 
p = np.zeros_like(model.time); p[-1]=1
final = model.Param(p)
model.Minimize(1e4*final*(v**2))
model.Minimize(1e4*final*((x-goal_dist)**2))
  1. Increase maximum iterations. Sometimes the solver needs more iterations to find a solution.
model.options.MAX_ITER=1000

The final version of the model has these changes. I may help converge to a solution and avoid maximum iterations or an infeasible solution.

from gekko import GEKKO
import numpy as np

min_velocity =  0
max_velocity = 10
max_decel    = -1
max_accel    =  1
distances    = np.linspace(0,20,21)
goal_dist    = 200

trip_time    = 100

# set up PWL functions
distances    = np.linspace(0,1000,10)
speed_limits = np.ones(10)*5
speed_limits[5:]=7
slope        = np.zeros(10)
slope[3:5]=1; slope[7:9]=-1

model = GEKKO(remote=True)
model.time = [i for i in range(trip_time)]

x = model.Var(value=0.0, lb=0)
v = model.Var(value=0.0, lb = min_velocity, ub = max_velocity)

v_max = model.Var()
slope_var = model.Var()

a = model.MV(value=0, lb=max_decel ,ub=max_accel)
a.STATUS = 1

#define vehicle movement
model.Equation(x.dt()==v)
model.Equation(v.dt()==a)

#aggregated velocity constraint
model.pwl(x, v_max, distances, speed_limits)
model.Equation(v_max>=v)

#slope is modeled as a piecewise linear function
model.pwl(x, slope_var, distances, slope)

#End state constraints
# vehicle must arrive at destination
#model.fix(x, pos=trip_time-1,val=goal_dist)
# vehicle must be fully stopped
#model.fix(v, pos=trip_time-1,val=0) 
p = np.zeros_like(model.time); p[-1]=1
final = model.Param(p)
model.Minimize(1e4*final*(v**2))
model.Minimize(1e4*final*((x-goal_dist)**2))

#VSPI Objective function
obj = (v * (1.1 * a + 9.81 * slope_var + 0.132) +0.0003002*pow(v, 3))
model.Minimize(obj)
# solve
model.options.IMODE = 6
model.options.REDUCE = 3
model.options.MAX_ITER=1000
model.solve(disp=True)

If you ask another question on StackOverflow, don't forget to include a minimal and complete working example that replicates the problem.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
  • 2
    First off, I would like to thank you for your help! Your answers on here have been extremely helpful in getting to know gekko and how it can help solve interesting questions. I will keep the executable code in mind next time I ask something on here as well. Regarding the solution to the problem, with your changes, a solution was found, however, it took quite some time, which I will now try to work on. – Salty Spatula Aug 03 '22 at 06:32
  • I'm glad it helped. If you need help on the solution time, please post another question or look at answers such as https://stackoverflow.com/questions/71935405/how-to-speed-up-integer-nonlinear-programming-with-1446-variables-in-python-gekk/71975512#71975512 or https://stackoverflow.com/questions/59791146/using-gekko-optimization-why-is-my-model-builder-so-much-slower-than-my-solver/59826356#59826356 for MINLP problems. – John Hedengren Aug 03 '22 at 11:32
  • 1
    thanks a lot. Unfortunately the pwl functions used as constraints don't behave the way I expected and i posted another question: https://stackoverflow.com/questions/73331561/gekko-python-adjusting-defined-piecewise-linear-function-during-solving I would greatly appreciate it if you could have a look at this as well – Salty Spatula Aug 12 '22 at 08:51