Implementing LP on Python

This notebook use the CPLEX optimization package

Libraries

from docplex.mp.model import Model
import numpy as np
import pandas as pd
import math
from matplotlib import pyplot as plt
plt.style.use('seaborn')

LP Function

####### LP-Optimization Function ########
def lp_program(t):
    ### create model object ##
    m = Model(name = "plattoon formation")
    
    ## add decision variables to model object ##
    a = m.continuous_var(name = "acceleration")
    
    ## add constraints to model object ##
    m.add_constraint( a_min <= a )
    m.add_constraint( a <= a_max)  

    
    ## define position functions ##
    x_lead = x_lead_arrival + v_lead_arrival*(t-t_lead_arrival)
    
    x_follower = x_follower_arrival + v_follower_arrival*(t-t_follower_arrival
                    ) + 0.5*a*(t**2-t_follower_arrival**2
                    ) - a*t_follower_arrival*(t-t_follower_arrival)
    ## define speed functions ##                         
    v_lead = v_lead_arrival
    v_follower = v_follower_arrival + a*(t - t_follower_arrival)
    
    
    
    m.add_constraint(x_lead - x_follower  - (v_lead)*sh >= 0.00)  
   
    # m.add_constraint(x_lead - x_follower - v_lead*sh <= 1) 
    m.add_constraint(v_follower <= v_lead)
    m.add_constraint(v_follower >= v_min)
                       
    m.minimize(x_lead - x_follower)
    
    # Solve
    s = m.solve()
    m.print_solution()
    variables = s.as_dict()
    # print(variables["acceleration"])
        
    return  round(variables["acceleration"],3)

Input Variables

########## Input Variables ###########
sim_time = 60
resolution = 1
t = float()
a_max = 1.4691
a_min = 0
dacc_min = 1.999
dacc_max = 0.899
v_max = 15.665
v_min = 0
sh = 1.5

####### Lead arrival_properties #######
x_lead_arrival = 1000
v_lead_arrival = 13.4
t_lead_arrival = 0

####### Follower arrival_properties #######
x_follower_arrival = 900
v_follower_arrival = 11
t_follower_arrival = 2

Simulation

data = {"time":[], "lead_pos": [], "foll_pos": [], "lead_speed": [], "foll_speed":[],
        "acc":[], "diff": []}
while t <= sim_time:
    print("\n",15*"#"," time is: ", t, 15*"#", "\n")
    try:
        a = lp_program(t)
        print(a)
    
        x_lead = x_lead_arrival + v_lead_arrival*(t-t_lead_arrival)
        print("leader_position ", x_lead)
        x_follower = x_follower_arrival + v_follower_arrival*(t-t_follower_arrival
                        ) + 0.5*a*(t**2-t_follower_arrival**2
                        ) - a*t_follower_arrival*(t-t_follower_arrival)
        print("follower_position ",x_follower)
        
        
        
        data["time"].append(t)
        data["lead_pos"].append(x_lead)
        data["foll_pos"].append(x_follower)
        data["lead_speed"].append(v_lead_arrival)
        data["foll_speed"].append(v_follower_arrival + a*(t - t_follower_arrival))
        data["acc"].append(a)
        data["diff"].append(x_lead - x_follower)  
        print(data)
        print(t, " ", v_follower_arrival, " ", x_follower, " ", x_follower_arrival, " ")
    except:
        print("No solution found")
        
    t += resolution

Visualization

df = pd.DataFrame.from_dict(data)

time = df.loc[:, "time"].values
lead_pos = df.loc[:, "lead_pos"].values
foll_pos = df.loc[:, "foll_pos"].values
lead_speed = df.loc[:, "lead_speed"].values
foll_speed = df.loc[:, "foll_speed"].values
acc = df.loc[:, "acc"].values

plt.style.use("seaborn")

fig, axes = plt.subplots(nrows = 1, ncols = 3)
custom_xlim = (0, 60)
plt.setp(axes, xlim=custom_xlim)

axes[0].plot(time, lead_pos, c = "k", label="Leader, $x(t)_{n-1}$")
axes[0].legend(loc = "lower right")
axes[0].plot(time, foll_pos, c = "b", label = "Follower, $x(t)_{n}$")
axes[0].legend(loc = "lower right")

axes[1].plot(time, lead_speed, c = "r", label="Des Speed")
axes[1].legend(loc = "lower right")
axes[1].plot(time, foll_speed, c = "c", label="$v(t)_{n-1}$")
axes[1].legend(loc = "lower right")

axes[2].bar(time, acc, label="$a(t)_{n-1}$")
axes[2].legend(loc = "upper right",)

axes[0].title.set_text('Position')
axes[1].title.set_text('Speed')
axes[2].title.set_text('Accelaration')
fig.suptitle('$x(t), v(t), a(t)$ vs time', fontsize=16)
plt.savefig("acc", dpi=300)