有了之前的三篇知识基础,这一篇主要是工程实践。
无人车辆在惯性坐标系中,车辆必须从一个给定的初始状态出发,然后通过某种控制器控制车辆尽量按照规划好的轨迹行驶。如下图所示: 如图中所示,期望轨迹是一条几何曲线 f ( x t ( t ) ) f(x_t(t)) f(xt(t)),自变量 x t x_t xt是时间 t t t的函数。 为了进一步深入研究,做出了以下假设:
环境感知部分无问题,规划出的路径为可行驶区域定位部分可以精确获得车辆的位置信息,包括车辆后轴中点的坐标信息,偏航角可以获取到车辆底层信息,车辆前轮的转角、车身速度等相关信息换句话说,就是轨迹跟踪是在行驶环境以及车辆内部状态信息完全已知的情况下进行的。
由上一篇文章可知,目前主要使用的三种车辆模型,由于不考虑车辆高速行驶的情况,本文采用第二种车辆运动学模型。 [ x ˙ y ˙ φ ˙ ] = [ cos φ sin φ tan δ / l ] ∗ v (1) \begin{bmatrix} \dot{x} \\ \dot{y} \\ \dot{\varphi} \\ \end{bmatrix} = \begin{bmatrix} \cos\varphi \\ \sin\varphi \\ \tan\delta/l \end{bmatrix}*v \tag{1} ⎣⎡x˙y˙φ˙⎦⎤=⎣⎡cosφsinφtanδ/l⎦⎤∗v(1) 式中, ( x , y ) (x,y) (x,y)为车辆后轴中心的坐标, φ \varphi φ为车辆的偏航角, δ \delta δ为车辆前轮转角, v v v为车辆后轴中心速度, l l l为车辆前后轴之间的距离。
系统化为一般形式,输入为 u ( v , δ ) u(v,\delta) u(v,δ)和状态量 χ ( x , y , φ ) \chi(x,y,\varphi) χ(x,y,φ),此时有 χ ˙ = f ( χ , u ) (2) \dot{\chi} = f(\chi,u)\tag{2} χ˙=f(χ,u)(2) 对于给定的参考轨迹,可以用参考车辆的运动轨迹描述,用 r r r代表参考量,参考轨迹可以转化为: χ r ˙ = f ( χ r , u r ) (3) \dot{\chi_r}=f(\chi_r,u_r) \tag{3} χr˙=f(χr,ur)(3) 式中, χ r = [ x r y r φ r ] T \chi_r=[x_r \space y_r \space \varphi_r]^T χr=[xr yr φr]T, u r = [ v r δ r ] u_r=[v_r \space \delta_r] ur=[vr δr]。
将式(2)用泰勒展开并忽略高阶项,有: χ ˙ = f ( χ r , u r ) + α f ( χ , u ) α χ ∣ χ = χ r u = u r ∗ ( χ − χ r ) + α f ( χ , u ) α u ∣ χ = χ r u = u r ∗ ( u − u r ) (4) \dot{\chi}=f(\chi_r,u_r)+\frac{\alpha f(\chi,u)}{\alpha\chi}\vert_{{\chi=\chi_r}\atop{u=u_r}}*(\chi-\chi_r)\\+ \frac{\alpha f(\chi,u)}{\alpha u}\vert_{{\chi=\chi_r}\atop{u=u_r}}*(u-u_r) \tag{4} χ˙=f(χr,ur)+αχαf(χ,u)∣u=urχ=χr∗(χ−χr)+αuαf(χ,u)∣u=urχ=χr∗(u−ur)(4)
式(4)与式(3)相减,可以得到 χ ~ ˙ = [ x ˙ − x ˙ r y ˙ − y ˙ r φ ˙ − φ ˙ r ] = [ 0 0 − v r ∗ s i n φ r 0 0 v r ∗ c o s φ r 0 0 0 ] [ x − x r y − y r φ − φ r ] + [ c o s φ r 0 s i n φ r 0 t a n δ r l v r l ∗ c o s 2 δ r ] [ v − v r δ − δ r ] (5) \dot{\tilde{\chi}}= \begin{bmatrix} \dot{x}-\dot{x}_r \\ \dot{y}-\dot{y}_r \\ \dot{\varphi}-\dot{\varphi}_r \end{bmatrix} \\ = \begin{bmatrix} 0 && 0 && -v_r*sin\varphi_r \\ 0 && 0&& v_r*cos\varphi_r \\ 0 && 0 && 0 \end{bmatrix} \begin{bmatrix} x-x_r \\ y-y_r \\ \varphi-\varphi_r \\ \end{bmatrix} + \begin{bmatrix} cos\varphi_r && 0 \\ sin\varphi_r && 0 \\ \frac{tan\delta_r}{l} && \frac{v_r}{l*cos^2\delta_r} \end{bmatrix} \begin{bmatrix} v-v_r \\ \delta-\delta_r \end{bmatrix} \tag{5} χ~˙=⎣⎡x˙−x˙ry˙−y˙rφ˙−φ˙r⎦⎤=⎣⎡000000−vr∗sinφrvr∗cosφr0⎦⎤⎣⎡x−xry−yrφ−φr⎦⎤+⎣⎡cosφrsinφrltanδr00l∗cos2δrvr⎦⎤[v−vrδ−δr](5) 式(5)即线性化的无人驾驶车辆误差模型。 工程中,可以参考的一般只有规划好的路点坐标而已,与上述的误差模型还是有区别的。
为了便于控制器的设计,对式(5)离散化处理后,得到 χ ~ ( k + 1 ) = A k , t ∗ χ ~ ( k ) + B k , t ∗ u ~ ( k ) (6) \tilde{\chi}(k+1)=A_{k,t}*\tilde{\chi}(k)+B_{k,t}*\tilde{u}(k) \tag{6} χ~(k+1)=Ak,t∗χ~(k)+Bk,t∗u~(k)(6) 式中,
A k , t = [ 1 0 − v r ∗ s i n φ r ∗ T 0 1 v r ∗ c o s φ r ∗ T 0 0 1 ] A_{k,t}=\begin{bmatrix} 1 && 0 && -v_r*sin\varphi_r*T \\ 0 && 1 && v_r*cos\varphi_r*T \\ 0 && 0 && 1 \end{bmatrix} Ak,t=⎣⎡100010−vr∗sinφr∗Tvr∗cosφr∗T1⎦⎤, B k , t = [ c o s φ r ∗ T 0 s i n φ r ∗ T 0 t a n δ f ∗ T l v r ∗ T l ∗ c o s 2 δ r ] B_{k,t}=\begin{bmatrix} cos\varphi_r*T && 0 \\ sin\varphi_r*T && 0 \\ \frac{tan\delta_f*T}{l} && \frac{v_r*T}{l*cos^2\delta_r} \end{bmatrix} Bk,t=⎣⎡cosφr∗Tsinφr∗Tltanδf∗T00l∗cos2δrvr∗T⎦⎤,T为采样时间。
首先,需要将式(6)转化为线性误差模型,令 ξ ( k ∣ t ) = [ χ ~ ( k ∣ t ) u ~ ( k − 1 ∣ t ) ] (7) \xi(k|t)=\begin{bmatrix} \tilde\chi(k|t) \\ \tilde u(k-1|t) \\ \end{bmatrix} \tag{7} ξ(k∣t)=[χ~(k∣t)u~(k−1∣t)](7) 得到新的空间状态表达式: ξ ( k + 1 ∣ t ) = A ~ k , t ∗ ξ ( k ∣ t ) + B ~ k , t ∗ Δ U ( k ∣ t ) η ( k ∣ t ) = C ~ k , t ∗ ξ ( k ∣ t ) (8) \xi(k+1|t)=\tilde A_{k,t}*\xi(k|t)+\tilde B_{k,t}*\Delta U(k|t) \\ \eta(k|t)=\tilde C_{k,t}*\xi(k|t) \tag{8} ξ(k+1∣t)=A~k,t∗ξ(k∣t)+B~k,t∗ΔU(k∣t)η(k∣t)=C~k,t∗ξ(k∣t)(8) 式中, A ~ k , t = [ A k , t B k , t 0 m ∗ n I m ] , B ~ k , t = [ B k , t I m ] \tilde A_{k,t}=\begin{bmatrix} A_{k,t} && B_{k,t} \\ 0_{m*n} && I_m \end{bmatrix}, \tilde B_{k,t}=\begin{bmatrix} B_{k,t} \\ I_m \\ \end{bmatrix} A~k,t=[Ak,t0m∗nBk,tIm],B~k,t=[Bk,tIm]
n n n为状态量维度, m m m为控制量维度。
为了简化计算,做出如下假设 A k , t = A t , t , k = 1 , … , t + N − 1 B k , t = B t , t , k = 1 , … , t + N − 1 (9) A_{k,t}=A_{t,t},\space k=1,\dots,t+N-1 \\ B_{k,t}=B_{t,t},\space k=1,\dots,t+N-1 \tag{9} Ak,t=At,t, k=1,…,t+N−1Bk,t=Bt,t, k=1,…,t+N−1(9) 推导后,得到系统的预测输出表达式: Y ( t ) = Ψ t ∗ ξ ( t ∣ t ) + Θ ∗ Δ U ( t ) (10) Y(t)=\Psi_t*\xi(t|t)+\varTheta*\Delta U(t) \tag{10} Y(t)=Ψt∗ξ(t∣t)+Θ∗ΔU(t)(10) 式中, Y ( t ) = [ η ( t + 1 ∣ t ) η ( t + 2 ∣ t ) … η ( t + N c ∣ t ) … η ( t + N p ∣ t ) ] Y(t)=\begin{bmatrix} \eta(t+1|t) \\ \eta(t+2|t) \\ \dots \\ \eta(t+N_c|t) \\ \dots \\ \eta(t+N_p|t) \\ \end{bmatrix} Y(t)=⎣⎢⎢⎢⎢⎢⎢⎡η(t+1∣t)η(t+2∣t)…η(t+Nc∣t)…η(t+Np∣t)⎦⎥⎥⎥⎥⎥⎥⎤, Ψ t = [ C ~ t , t ∗ A ~ t , t C ~ t , t ∗ A ~ t , t 2 … C ~ t , t ∗ A ~ t , t N c … C ~ t , t ∗ A ~ t , t N p ] \Psi_t=\begin{bmatrix} \tilde C_{t,t}*\tilde A_{t,t} \\ \tilde C_{t,t}*\tilde A_{t,t}^2 \\ \dots \\ \tilde C_{t,t}*\tilde A_{t,t}^{N_c} \\ \dots \\ \tilde C_{t,t}*\tilde A_{t,t}^{N_p} \\ \end{bmatrix} Ψt=⎣⎢⎢⎢⎢⎢⎢⎡C~t,t∗A~t,tC~t,t∗A~t,t2…C~t,t∗A~t,tNc…C~t,t∗A~t,tNp⎦⎥⎥⎥⎥⎥⎥⎤, Θ t = [ C ~ t , t B ~ t , t 0 0 0 C ~ t , t A ~ t , t B ~ t , t C ~ t , t B ~ t , t 0 0 ⋯ ⋯ ⋱ ⋯ C ~ t , t A ~ t , t N c − 1 B ~ t , t C ~ t , t A ~ t , t N c − 2 B ~ t , t ⋯ C ~ t , t B ~ t , t C ~ t , t A ~ t , t N c B ~ t , t C ~ t , t A ~ t , t N c − 1 B ~ t , t ⋯ C ~ t , t A ~ t , t B ~ t , t ⋮ ⋮ ⋱ ⋮ C ~ t , t A ~ t , t N p − 1 B ~ t , t C ~ t , t A ~ t , t N p − 2 B ~ t , t ⋯ C ~ t , t A ~ t , t N p − N c − 1 B ~ t , t ] \varTheta_{t}= \begin{bmatrix} \tilde{C}_{t,t}\tilde{B}_{t,t} && 0 && 0 && 0 \\ \tilde{C}_{t,t}\tilde{A}_{t,t}\tilde{B}_{t,t} && \tilde{C}_{t,t}\tilde{B}_{t,t} && 0 && 0 \\ \cdots && \cdots && \ddots && \cdots \\ \tilde{C}_{t,t}\tilde{A}^{N_c - 1}_{t,t}\tilde{B}_{t,t} && \tilde{C}_{t,t}\tilde{A}^{N_c - 2}_{t,t}\tilde{B}_{t,t} && \cdots && \tilde{C}_{t,t}\tilde{B}_{t,t} \\ \tilde{C}_{t,t}\tilde{A}^{N_c}_{t,t}\tilde{B}_{t,t} && \tilde{C}_{t,t}\tilde{A}^{N_c - 1}_{t,t}\tilde{B}_{t,t} && \cdots && \tilde{C}_{t,t}\tilde{A}_{t,t}\tilde{B}_{t,t} \\ \vdots && \vdots && \ddots && \vdots \\ \tilde{C}_{t,t}\tilde{A}^{N_{p}-1}_{t,t}\tilde{B}_{t,t} && \tilde{C}_{t,t}\tilde{A}^{N_{p}-2}_{t,t}\tilde{B}_{t,t} && \cdots && \tilde{C}_{t,t}\tilde{A}^{N_{p} - N_{c} -1}_{t,t}\tilde{B}_{t,t} \end{bmatrix} Θt=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡C~t,tB~t,tC~t,tA~t,tB~t,t⋯C~t,tA~t,tNc−1B~t,tC~t,tA~t,tNcB~t,t⋮C~t,tA~t,tNp−1B~t,t0C~t,tB~t,t⋯C~t,tA~t,tNc−2B~t,tC~t,tA~t,tNc−1B~t,t⋮C~t,tA~t,tNp−2B~t,t00⋱⋯⋯⋱⋯00⋯C~t,tB~t,tC~t,tA~t,tB~t,t⋮C~t,tA~t,tNp−Nc−1B~t,t⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ Δ U ( t ) = [ Δ u ( t ∣ t ) Δ u ( t + 1 ∣ t ) ⋯ Δ u ( t + N c ∣ t ) ] \Delta U(t)=\begin{bmatrix} \Delta u(t|t) \\ \Delta u(t+1|t) \\ \cdots \\ \Delta u(t+N_c|t) \\ \end{bmatrix} ΔU(t)=⎣⎢⎢⎡Δu(t∣t)Δu(t+1∣t)⋯Δu(t+Nc∣t)⎦⎥⎥⎤ 。
目标函数要能够保证无人驾驶车辆快速且平稳地追踪期望轨迹。以下介绍基于无人驾驶车辆线性误差模型预测车辆未来时刻的输出。 一般采用如下形式的目标函数: J ( k ) = ∑ i = 1 N p ∥ η ( k + i ∣ t ) − η r e f ( k + i ∣ t ) ∥ Q 2 + ∑ i = 1 N c − 1 ∥ Δ U ( k + i ∣ t ) ∥ R 2 + ρ ∗ ε 2 (11) J(k) \space =\space \sum\limits_{i=1}^{N_p} \begin{Vmatrix} \eta(k+i|t) \space - \space \eta_{ref}(k+i|t) \end{Vmatrix}_Q^2 \space + \space \sum\limits_{i=1}^{N_c-1} \begin{Vmatrix} \Delta U(k+i|t) \end{Vmatrix}_R^2 \space + \space \rho*\varepsilon^2 \tag{11} J(k) = i=1∑Np∥∥η(k+i∣t) − ηref(k+i∣t)∥∥Q2 + i=1∑Nc−1∥∥ΔU(k+i∣t)∥∥R2 + ρ∗ε2(11) 式中, N p N_p Np为预测时域, N c N_c Nc为控制时域, Q , R , ρ Q,R,\rho Q,R,ρ为权重系数, ε \varepsilon ε为松弛因子,加入松弛因子是为了防止求解过程出现无可行解的情况, η ( k ∣ t ) \eta(k|t) η(k∣t)为系统状态量, Δ U ( k ∣ t ) \Delta U(k|t) ΔU(k∣t)为系统控制量。 第一项表示控制系统对期望轨迹的跟踪能力,第二项表示控制系统控制量的约束情况。
实际过程中,受到实际控制过程中,不同控制器之间的控制量的可控制范围不同,大部分情况下控制量可以表达为如下形式: u m i n ( t + k ) ≤ u t + k ≤ u m a x ( t + k ) k = 0 , 1 , … , N c − 1 (12) u_{min}(t+k) \space \leq \space u_{t+k} \space \leq \space u_{max}(t+k) \\ \space\space\space\space\space\space\space k \space = \space 0,1,\dots,N_c-1 \tag{12} umin(t+k) ≤ ut+k ≤ umax(t+k) k = 0,1,…,Nc−1(12) 控制增量约束表达式: Δ u m i n ( t + k ) ≤ Δ u t + k ≤ Δ u m a x ( t + k ) k = 0 , 1 , … , N c − 1 (13) \Delta u_{min}(t+k) \space \leq \space \Delta u_{t+k} \space \leq \space \Delta u_{max}(t+k) \\ \space\space\space\space\space\space\space k \space = \space 0,1,\dots,N_c-1 \tag{13} Δumin(t+k) ≤ Δut+k ≤ Δumax(t+k) k = 0,1,…,Nc−1(13)
将式(12),(13)转换为以下形式: U m i n ≤ A ∗ Δ U t + U t ≤ U m a x U_{min} \space \leq \space A*\Delta U_t \space + \space U_t \space \leq \space U_{max} Umin ≤ A∗ΔUt + Ut ≤ Umax 式中, U m i n , U m a x U_{min},U_{max} Umin,Umax分别为控制时域内的控制量最小值、最大值集合。 将目标函数转化为标准二次型并结合约束条件,可以解决如下形式的优化问题: J ( ξ ( t ) , u ( t − 1 ) , Δ U ( t ) ) = [ Δ U ( t ) T , ε ] T ∗ H t ∗ [ Δ U ( t ) T , ε ] + G t ∗ [ Δ U ( t ) T , ε ] (14) J(\xi(t),u(t-1),\Delta U(t)) \space = \space [\Delta U(t)^T,\varepsilon]^T*H_t*[\Delta U(t)^T,\varepsilon] \space + \space G_t*[\Delta U(t)^T,\varepsilon] \tag{14} J(ξ(t),u(t−1),ΔU(t)) = [ΔU(t)T,ε]T∗Ht∗[ΔU(t)T,ε] + Gt∗[ΔU(t)T,ε](14) 式中, H t = [ Θ t T ∗ Q ∗ Θ t + R 0 0 ρ ] H_t = \begin{bmatrix} \Theta_t^T*Q*\Theta_t+R && 0 \\ 0 && \rho \end{bmatrix} Ht=[ΘtT∗Q∗Θt+R00ρ], G t = [ 2 ∗ e t T ∗ Q ∗ Θ 0 ] G_t \space = \space \begin{bmatrix} 2*e_t^T*Q*\Theta && 0 \end{bmatrix} Gt = [2∗etT∗Q∗Θ0], e t e_t et为预测时域内的跟踪误差。 每个空置周期内对式(14)求解后,得到了一组控制时域内的一系列控制输入增量: Δ U t ∗ = [ Δ u t ∗ , Δ u t + 1 ∗ , … , Δ u t + N c − 1 ∗ ] T (15) \Delta U_t^* \space = \space [\Delta u_t^*,\Delta u_{t+1}^*,\dots,\Delta u_{t+N_c-1}^*]^T \tag{15} ΔUt∗ = [Δut∗,Δut+1∗,…,Δut+Nc−1∗]T(15) 取控制序列的第一个控制增量与上一时刻的控制量作用于系统, u ( t ) = u ( t − 1 ) + Δ u t ∗ (16) u(t) \space = \space u(t-1) \space + \space \Delta u_t^* \tag{16} u(t) = u(t−1) + Δut∗(16) 进入下一周期后,重复上述的计算,即可实现对车辆的轨迹跟踪控制。