簡介

自動(dòng)駕駛運(yùn)動(dòng)規(guī)劃中會(huì)用到各種曲線，主要用于生成車輛的軌跡，常見的軌跡生成算法，如貝塞爾曲線，樣條曲線，以及apollo EM Planner的五次多項(xiàng)式曲線，城市場景中使用的是分段多項(xiàng)式曲線，在EM Planner和Lattice Planner中思路是，都是先通過動(dòng)態(tài)規(guī)劃生成點(diǎn)，再用5次多項(xiàng)式生成曲線連接兩個(gè)點(diǎn)（雖然后面的版本改動(dòng)很大，至少lattice planner目前還是這個(gè)方法）。

在這里可以看出5次多項(xiàng)式的作用，就是生成軌跡，這里的軌跡不一定是車行駛的軌跡，比如S—T圖中的線，是用來做速度規(guī)劃的。 如下圖：

在apollo里面用到了，3-5次多項(xiàng)式，

• cubic_polynomial_curve1d，三次多項(xiàng)式曲線
• quartic_polynomial_curve1d，四次多項(xiàng)式曲線-
• quintic_polynomial_curve1d，五次多項(xiàng)式曲線

3次多項(xiàng)式

$$f(x)=C_0+C_{1}x+C_2{x}^2+C_3{x}^3????(1)$$
對(duì)x求導(dǎo)====>
$$f^{\prime }(x)=C_1+2C_{2}x+3C_3{x}^2$$
$$f(x)=2C_2+6C_{3}x$$
因?yàn)槭沁B接兩個(gè)已知點(diǎn)，所以兩個(gè)軌跡點(diǎn)的信息是已知的。
如已知，當(dāng)x=0，代如得：

所以

將終點(diǎn) (X_end，f(X_end))帶入得到C_3。

$$C_3=(f(x_p)?C_0?C_1{x}_p?C_2{x}_p^2)/x_p^3$$

（1）式的所有未知參數(shù)都求出，代入x可以得到兩點(diǎn)間的軌跡。

apollo中3次多項(xiàng)式的計(jì)算

// cubic_polynomial_curve1d.cc
CubicPolynomialCurve1d::CubicPolynomialCurve1d(const double x0,
                                               const double dx0,
                                               const double ddx0,
                                               const double x1,
                                               const double param) {
  ComputeCoefficients(x0, dx0, ddx0, x1, param);
  param_ = param;
  start_condition_[0] = x0;
  start_condition_[1] = dx0;
  start_condition_[2] = ddx0;
  end_condition_ = x1;
}

void CubicPolynomialCurve1d::ComputeCoefficients(const double x0,
                                                 const double dx0,
                                                 const double ddx0,
                                                 const double x1,
                                                 const double param) {
  DCHECK(param > 0.0);
  const double p2 = param * param;
  const double p3 = param * p2;
  coef_[0] = x0;
  coef_[1] = dx0;
  coef_[2] = 0.5 * ddx0;
  coef_[3] = (x1 - x0 - dx0 * param - coef_[2] * p2) / p3;
}

5次多項(xiàng)式

$$f(x)=C_0+C_{1}x+C_2{x}^2+C_3{x}^3+C_4{x}^4+C_5{x}^5????(3?1)$$
對(duì)x求導(dǎo)====>
$$f^{\prime }(x)=C_1+2C_{2}x+3C_3{x}^2+4C_4{x}^3+5C_5{x}^4????(3?2)$$

$$f^{\prime \prime }(x)=2C_2+6C_{3}x+12C_4{x}^2+20C_5{x}^3????(3?3)$$
因?yàn)槭沁B接兩個(gè)已知點(diǎn)，所以兩個(gè)軌跡點(diǎn)的信息是已知的。
如已知，當(dāng)x=0，代如得：

所以

將終點(diǎn) (x_e，f(x_e))帶入 (3-1) (3-2) (3-3)。

得到
$$f(x_e)=C_0+C_1{x}_e+C_2{x}_e^2+C_3{x}_e^3+C_4{x}_e^4+C_5{x}_e^5????(3?4)$$
對(duì)x求導(dǎo)====>
$$f^{\prime }(x)=C_1+2C_2{x}_e+3C_3{x}_e^2+4C_4{x}_e^3+5C_5{x}_e^4????(3?5)$$

$$f^{\prime \prime }(x)=2C_2+6C_3{x}_e+12C_4{x}_e^2+20C_5{x}_e^3????(3?6)$$

聯(lián)立可得
$$C_3=\frac{10(f(x_p)?1/2?f^{\prime \prime }(0)x_p^2?f^{\prime }(0)x_p?f(0))?4(f^{\prime }(x_p)?f^{\prime \prime }(0)x_p?f^{\prime }(0))x_p+1/2?(f^{\prime \prime }(x_p)?f^{\prime \prime }(0))x_p^2}{x_p^3}$$

$$C_4=\frac{?15(f(x_p)?1/2?f^{\prime \prime }(0)x_p^2?f^{\prime }(0)x_p?f(0))+7(f^{\prime }(x_p)?f^{\prime \prime }(0)x_p?f^{\prime }(0))x_p?1/2?(f^{\prime \prime }(x_p)?f^{\prime \prime }(0))x_p^2}{x_p^4}$$

$$C_5=\frac{6(f(x_p)?1/2?f^{\prime \prime }(0)x_p^2?f^{\prime }(0)x_p?f(0))?3(f^{\prime }(x_p)?f^{\prime \prime }(0)x_p?f^{\prime }(0))x_p+1/2?(f^{\prime \prime }(x_p)?f^{\prime \prime }(0))x_p^2}{x_p^5}$$

apollo中5次多項(xiàng)式的計(jì)算

void QuinticPolynomialCurve1d::ComputeCoefficients(
    const double x0, const double dx0, const double ddx0, const double x1,
    const double dx1, const double ddx1, const double p) {
  CHECK_GT(p, 0.0);

  coef_[0] = x0;
  coef_[1] = dx0;
  coef_[2] = ddx0 / 2.0;

  const double p2 = p * p;
  const double p3 = p * p2;

  // the direct analytical method is at least 6 times faster than using matrix
  // inversion.
  const double c0 = (x1 - 0.5 * p2 * ddx0 - dx0 * p - x0) / p3;
  const double c1 = (dx1 - ddx0 * p - dx0) / p2;
  const double c2 = (ddx1 - ddx0) / p;

  coef_[3] = 0.5 * (20.0 * c0 - 8.0 * c1 + c2);
  coef_[4] = (-15.0 * c0 + 7.0 * c1 - c2) / p;
  coef_[5] = (6.0 * c0 - 3.0 * c1 + 0.5 * c2) / p2;
}

5次多項(xiàng)式擬合

import math

import matplotlib.pyplot as plt
import numpy as np

# parameter
MAX_T = 100.0  # maximum time to the goal [s]
MIN_T = 5.0  # minimum time to the goal[s]

show_animation = True

class QuinticPolynomial:

    def __init__(self, xs, vxs, axs, xe, vxe, axe, time):
        # calc coefficient of quintic polynomial
        # See jupyter notebook document for derivation of this equation.
        # 根據(jù)初始狀態(tài)c0 c1 c2
        self.a0 = xs # x(t)
        self.a1 = vxs # v(t)
        self.a2 = axs / 2.0 # a(t)

        A = np.array([[time ** 3, time ** 4, time ** 5],
                      [3 * time ** 2, 4 * time ** 3, 5 * time ** 4],
                      [6 * time, 12 * time ** 2, 20 * time ** 3]])
        b = np.array([xe - self.a0 - self.a1 * time - self.a2 * time ** 2,
                      vxe - self.a1 - 2 * self.a2 * time,
                      axe - 2 * self.a2])
        # Ax=b 求解x
        x = np.linalg.solve(A, b)
        # 計(jì)算c3 c4 c5
        self.a3 = x[0]
        self.a4 = x[1]
        self.a5 = x[2]
    # 根據(jù)時(shí)間計(jì)算點(diǎn)x(t)
    def calc_point(self, t):
        xt = self.a0 + self.a1 * t + self.a2 * t ** 2 + \
             self.a3 * t ** 3 + self.a4 * t ** 4 + self.a5 * t ** 5
        return xt
    # 計(jì)算v(t)
    def calc_first_derivative(self, t):
        xt = self.a1 + 2 * self.a2 * t + \
             3 * self.a3 * t ** 2 + 4 * self.a4 * t ** 3 + 5 * self.a5 * t ** 4

        return xt
    # 計(jì)算a(t)  
    def calc_second_derivative(self, t):
        xt = 2 * self.a2 + 6 * self.a3 * t + 12 * self.a4 * t ** 2 + 20 * self.a5 * t ** 3

        return xt
    # 計(jì)算jerk(t)  
    def calc_third_derivative(self, t):
        xt = 6 * self.a3 + 24 * self.a4 * t + 60 * self.a5 * t ** 2

        return xt


def quintic_polynomials_planner(sx, sy, syaw, sv, sa, gx, gy, gyaw, gv, ga, max_accel, max_jerk, dt):
    """
    quintic polynomial planner

    input
        s_x: start x position [m]
        s_y: start y position [m]
        s_yaw: start yaw angle [rad]
        sa: start accel [m/ss]
        gx: goal x position [m]
        gy: goal y position [m]
        gyaw: goal yaw angle [rad]
        ga: goal accel [m/ss]
        max_accel: maximum accel [m/ss]
        max_jerk: maximum jerk [m/sss]
        dt: time tick [s]

    return
        time: time result
        rx: x position result list
        ry: y position result list
        ryaw: yaw angle result list
        rv: velocity result list
        ra: accel result list

    """

    vxs = sv * math.cos(syaw)# 起點(diǎn)速度在x方向分量
    vys = sv * math.sin(syaw)# 起點(diǎn)速度在y方向分量
    vxg = gv * math.cos(gyaw)# 終點(diǎn)速度在x方向分量
    vyg = gv * math.sin(gyaw)# 終點(diǎn)速度在y方向分量

    axs = sa * math.cos(syaw)# 起點(diǎn)加速度在x方向分量
    ays = sa * math.sin(syaw)# 終點(diǎn)加速度在y方向分量
    axg = ga * math.cos(gyaw)# 起點(diǎn)加速度在x方向分量
    ayg = ga * math.sin(gyaw)# 終點(diǎn)加速度在y方向分量

    time, rx, ry, ryaw, rv, ra, rj = [], [], [], [], [], [], []

    for T in np.arange(MIN_T, MAX_T, MIN_T):
        xqp = QuinticPolynomial(sx, vxs, axs, gx, vxg, axg, T)
        yqp = QuinticPolynomial(sy, vys, ays, gy, vyg, ayg, T)

        time, rx, ry, ryaw, rv, ra, rj = [], [], [], [], [], [], []

        for t in np.arange(0.0, T + dt, dt):
            time.append(t)
            rx.append(xqp.calc_point(t))
            ry.append(yqp.calc_point(t))

            vx = xqp.calc_first_derivative(t)
            vy = yqp.calc_first_derivative(t)
            v = np.hypot(vx, vy)
            yaw = math.atan2(vy, vx)
            rv.append(v)
            ryaw.append(yaw)

            ax = xqp.calc_second_derivative(t)
            ay = yqp.calc_second_derivative(t)
            a = np.hypot(ax, ay)
            if len(rv) >= 2 and rv[-1] - rv[-2] < 0.0:
                a *= -1
            ra.append(a)

            jx = xqp.calc_third_derivative(t)
            jy = yqp.calc_third_derivative(t)
            j = np.hypot(jx, jy)
            if len(ra) >= 2 and ra[-1] - ra[-2] < 0.0:
                j *= -1
            rj.append(j)

        if max([abs(i) for i in ra]) <= max_accel and max([abs(i) for i in rj]) <= max_jerk:
            print("find path!!")
            break

    if show_animation:  # pragma: no cover
        for i, _ in enumerate(time):
            plt.cla()
            # for stopping simulation with the esc key.
            plt.gcf().canvas.mpl_connect('key_release_event',
                                         lambda event: [exit(0) if event.key == 'escape' else None])
            plt.grid(True)
            plt.axis("equal")
            plot_arrow(sx, sy, syaw)
            plot_arrow(gx, gy, gyaw)
            plot_arrow(rx[i], ry[i], ryaw[i])
            plt.title("Time[s]:" + str(time[i])[0:4] +
                      " v[m/s]:" + str(rv[i])[0:4] +
                      " a[m/ss]:" + str(ra[i])[0:4] +
                      " jerk[m/sss]:" + str(rj[i])[0:4],
                      )
            plt.pause(0.001)

    return time, rx, ry, ryaw, rv, ra, rj


def plot_arrow(x, y, yaw, length=1.0, width=0.5, fc="r", ec="k"):  # pragma: no cover
    """
    Plot arrow
    """

    if not isinstance(x, float):
        for (ix, iy, iyaw) in zip(x, y, yaw):
            plot_arrow(ix, iy, iyaw)
    else:
        plt.arrow(x, y, length * math.cos(yaw), length * math.sin(yaw),
                  fc=fc, ec=ec, head_width=width, head_length=width)
        plt.plot(x, y)


def main():

    sx = 10.0  # start x position [m]
    sy = 10.0  # start y position [m]
    syaw = np.deg2rad(10.0)  # start yaw angle [rad]
    sv = 1.0  # start speed [m/s]
    sa = 0.1  # start accel [m/ss]
    gx = 30.0  # goal x position [m]
    gy = -10.0  # goal y position [m]
    gyaw = np.deg2rad(20.0)  # goal yaw angle [rad]
    gv = 1.0  # goal speed [m/s]
    ga = 0.1  # goal accel [m/ss]
    max_accel = 1.0  # max accel [m/ss]
    max_jerk = 0.5  # max jerk [m/sss]
    dt = 0.1  # time tick [s]

    time, x, y, yaw, v, a, j = quintic_polynomials_planner(
        sx, sy, syaw, sv, sa, gx, gy, gyaw, gv, ga, max_accel, max_jerk, dt)

    if show_animation:  # pragma: no cover
        plt.plot(x, y, "-r")

        plt.subplots()
        plt.plot(time, [np.rad2deg(i) for i in yaw], "-r")
        plt.xlabel("Time[s]")
        plt.ylabel("Yaw[deg]")
        plt.grid(True)

        plt.subplots()
        plt.plot(time, v, "-r")
        plt.xlabel("Time[s]")
        plt.ylabel("Speed[m/s]")
        plt.grid(True)

        plt.subplots()
        plt.plot(time, a, "-r")
        plt.xlabel("Time[s]")
        plt.ylabel("accel[m/ss]")
        plt.grid(True)

        plt.subplots()
        plt.plot(time, j, "-r")
        plt.xlabel("Time[s]")
        plt.ylabel("jerk[m/sss]")
        plt.grid(True)

        plt.show()


if __name__ == '__main__':
    main()

上圖模擬了車輛曲線的擬合過程的軌跡，箭頭代表方向。

過程的速度

感覺5次多項(xiàng)式和貝塞爾曲線還是挺像的，有時(shí)間寫一篇五次多項(xiàng)式、貝塞爾曲線、樣條曲線的對(duì)比。

還有一個(gè)問題就是為啥要用5次多項(xiàng)式，而不是4次，6次更高階或者低階呢？下回復(fù)習(xí)一下老王的講解。
參考：文章來源地址http://www.zghlxwxcb.cn/news/detail-436810.html

• https://blog.csdn.net/qq_36458461/article/details/110007656
• https://zhuanlan.zhihu.com/p/567375631

到了這里，關(guān)于自動(dòng)駕駛規(guī)劃 - 5次多項(xiàng)式擬合的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！