Gekko Python adjusting defined piecewise linear function during solving

user_13312159 0
Viewed 36

I am trying to solve a dynamic optimization problem using gekko. The goal is to minimize a form of energy consumption represented by VSP over a set distance under speed constraints. I define a piece-wise linear function as a the speed constraint and to model the slope of the road at different distances:

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,200,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=False)
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<=v_max)
#slope is modeled as a piecewise linear function
model.pwl(x, slope_var, distances, slope)
#End state constraints
p = np.zeros_like(model.time); p[-1]=1
final = model.Param(p)

model.Minimize(1e4*final*(v**2))# vehicle must be fully stopped
model.Minimize(1e4*final*((x-goal_dist)**2))# vehicle must arrive at destination

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

plt.plot(x.value, v_max.value, 'b-', label = r'$vmaxvals$')
plt.plot(x.value , v.value,'g-',label=r'$vopt$')
plt.plot(x.value, a.value, 'c-', label=r'$accel$')
plt.plot(x.value, slope_var.value, 'r-', label=r'$slope$')
plt.plot([i*20 for i in range(10)], slope, 'mx', label=r'$orig_slope$')
plt.plot([i*20 for i in range(10)], speed_limits, 'kx', label=r'$orig_spd_limit$')
plt.legend(loc='best')
plt.xlabel('Distance Covered')
plt.show()
print(model.options.APPSTATUS)

Unfortunately, however, the values of slope_var and v_max get adjusted in the process of solving the problem. I am sure this is intended in this case, so is there a way to fix these PWL functions in place similar to a Parameter? If I use a cspline object to apprximate the speed limits and slope, the values dont't change since it is pre-built as far as I understand, however, the accuracy of a cubic spline is limited to a few data points and few changes in slope, which is why I would like to model it using a piecewise linear function.

The Plot of my pwl values and the measured values of slope and speed limits

gekko python
1 Answers

The pwl function does give a linear interpolation but it relies on a Mathematical Program with Complementarity Constraints (MPCCs) that are challenging to solve with many local minima at saddle points. You mentioned that you don't want to use the cspline function, but it may be your best option. There are some slight errors at the transition points, but it can be fixed by adding additional points during transitions or by increasing the resolution.

cspline results

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

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,200,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=False)
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.cspline(x,v_max,distances,speed_limits,True)
#model.pwl(x, v_max, distances, speed_limits)
model.Equation(v<=v_max)
#slope is modeled as a piecewise linear function
#model.pwl(x, slope_var, distances, slope)
model.cspline(x,slope_var,distances,slope,True)
#End state constraints
p = np.zeros_like(model.time); p[-1]=1
final = model.Param(p)

model.Minimize(1e4*final*(v**2))# vehicle must be fully stopped
model.Minimize(1e4*final*((x-goal_dist)**2))# vehicle must arrive at destination

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

plt.plot(x.value, v_max.value, 'b-', label = 'vmaxvals')
plt.plot(x.value , v.value,'g-',label='vopt')
plt.plot(x.value, a.value, 'c-', label='accel')
plt.plot(x.value, slope_var.value, 'r-', label='slope')
plt.plot(distances, slope, 'mx', label='orig_slope')
plt.plot(distances, speed_limits, 'kx', label='orig_spd_limit')
plt.legend(loc='best')
plt.xlabel('Distance Covered')
plt.show()
print(model.options.APPSTATUS)

There was an error in the plotting with:

plt.plot([i*20 for i in range(10)], slope, 'mx', label=r'$orig_slope$')
plt.plot([i*20 for i in range(10)], speed_limits, 'kx', label=r'$orig_spd_limit$')
p

Use this instead:

plt.plot(distances, slope, 'mx', label='orig_slope')
plt.plot(distances, speed_limits, 'kx', label='orig_spd_limit')
user_23669410
Related