mpc_python_learn/mpc_pybullet_demo/cvxpy_mpc.py

import numpy as np
np.seterr(divide='ignore', invalid='ignore')

from scipy.integrate import odeint
from scipy.interpolate import interp1d
import cvxpy as cp

from utils import *

from mpc_config import Params
P=Params()

def get_linear_model(x_bar,u_bar):
    """
    """
    L=0.3

    x = x_bar[0]
    y = x_bar[1]
    v = x_bar[2]
    theta = x_bar[3]

    a = u_bar[0]
    delta = u_bar[1]

    A = np.zeros((P.N,P.N))
    A[0,2]=np.cos(theta)
    A[0,3]=-v*np.sin(theta)
    A[1,2]=np.sin(theta)
    A[1,3]=v*np.cos(theta)
    A[3,2]=v*np.tan(delta)/L
    A_lin=np.eye(P.N)+P.dt*A

    B = np.zeros((P.N,P.M))
    B[2,0]=1
    B[3,1]=v/(L*np.cos(delta)**2)
    B_lin=P.dt*B

    f_xu=np.array([v*np.cos(theta), v*np.sin(theta), a,v*np.tan(delta)/L]).reshape(P.N,1)
    C_lin = P.dt*(f_xu - np.dot(A,x_bar.reshape(P.N,1)) - np.dot(B,u_bar.reshape(P.M,1)))

    return np.round(A_lin,4), np.round(B_lin,4), np.round(C_lin,4)


def optimize(state,u_bar,track,ref_vel=1.):
    '''
    :param state:
    :param u_bar:
    :param track:
    :returns:
    '''

    MAX_SPEED = 1.5 #m/s
    MAX_ACC = 1.0 #m/ss
    MAX_D_ACC = 1.0 #m/sss
    MAX_STEER = np.radians(30) #rad
    MAX_D_STEER = np.radians(30) #rad/s

    REF_VEL = ref_vel #m/s
    
    # dynamics starting state
    x_bar = np.zeros((P.N,P.T+1))
    x_bar[:,0] = state

    #prediction for linearization of costrains
    for t in range (1,P.T+1):
        xt=x_bar[:,t-1].reshape(P.N,1)
        ut=u_bar[:,t-1].reshape(P.M,1)
        A,B,C=get_linear_model(xt,ut)
        xt_plus_one = np.squeeze(np.dot(A,xt)+np.dot(B,ut)+C)
        x_bar[:,t]= xt_plus_one

    #CVXPY Linear MPC problem statement
    cost = 0
    constr = []
    x = cp.Variable((P.N, P.T+1))
    u = cp.Variable((P.M, P.T))
    
    # Cost Matrices
    Q = np.diag([20,20,10,0]) #state error cost
    Qf = np.diag([10,10,10,10]) #state final error cost
    R = np.diag([10,10])       #input cost
    R_ = np.diag([10,10])      #input rate of change cost

    #Get Reference_traj
    x_ref, d_ref = get_ref_trajectory(x_bar[:,0] ,track, REF_VEL)

    #Prediction Horizon
    for t in range(P.T):

        # Tracking Error
        cost += cp.quad_form(x[:,t] - x_ref[:,t], Q)

        # Actuation effort
        cost += cp.quad_form(u[:,t], R)

        # Actuation rate of change
        if t < (P.T - 1):
            cost += cp.quad_form(u[:,t+1] - u[:,t], R_)
            constr+= [cp.abs(u[0, t + 1] - u[0, t])/P.dt <= MAX_D_ACC]    #max acc rate of change
            constr += [cp.abs(u[1, t + 1] - u[1, t])/P.dt <= MAX_D_STEER] #max steer rate of change

        # Kinrmatics Constrains (Linearized model)
        A,B,C = get_linear_model(x_bar[:,t], u_bar[:,t])
        constr += [x[:,t+1] == A@x[:,t] + B@u[:,t] + C.flatten()]

    # sums problem objectives and concatenates constraints.
    constr += [x[:,0] == x_bar[:,0]]           #starting condition
    constr += [x[2,:] <= MAX_SPEED]           #max speed
    constr += [x[2,:] >= 0.0]                 #min_speed (not really needed)
    constr += [cp.abs(u[0,:]) <= MAX_ACC]     #max acc
    constr += [cp.abs(u[1,:]) <= MAX_STEER]   #max steer

    # Solve
    prob = cp.Problem(cp.Minimize(cost), constr)
    prob.solve(solver=cp.OSQP, verbose=False)

    if "optimal" not in prob.status:
        print("WARN: No optimal solution")
        return u_bar

    #retrieved optimized U and assign to u_bar to linearize in next step
    u_opt=np.vstack((np.array(u.value[0, :]).flatten(),
                    (np.array(u.value[1, :]).flatten())))

    return u_opt