1. ilqr
ILQR算法是基于nominal trajectory ( x ~ , u ~ ) (\tilde{x}, \tilde{u}) (x~,u~)来优化求解的。ILQR是求解状态变量和控制变量的增量序列 ( δ x ∗ , δ u ∗ ) (\delta x^*, \delta u^*) (δx∗,δu∗)求解轨迹的局部最优值。
1.1 无约束轨迹优化问题形式
x ∗ , u ∗ = arg min x , u [ ∑ 0 N − 1 L k ( x k , u k ) + L F ( x N ) ] s . t . { x k + 1 = f ( x k , u k ) x 0 = x s t a r t \begin{split} & x^* ,u^* = \mathop{\arg \min}\limits_{x, u} [\sum_0^{N-1}L^k(x_k, u_k) + L_F(x_N) ] \\ & s.t. \quad \left\{\begin{aligned} x_{k+1} &= f(x_k, u_k) \\ x_0 &= x_{start} \end{aligned}\right. \end{split} x∗,u∗=x,uargmin[0∑N−1Lk(xk,uk)+LF(xN)]s.t.{xk+1x0=f(xk,uk)=xstart
1.2 Backward Pass
V k V^k Vk是第 k ( k ∈ [ 0 , N ] ) k(k \in [0,N]) k(k∈[0,N])的最优 cost-to-go,根据Bellman方程有:
{ V k ( x N ) = L F ( x N ) V k ( x k ) = min u k [ L k ( x k , u k ) + V k + 1 ( f ( x k , u k ) ) ] \left\{\begin{aligned} V^k(x_N) &= L_F(x_N) \\ V^k(x_k) &= \min\limits_{u_k}[L^k(x_k, u_k) + V^{k+1}(f(x_k, u_k))] \end{aligned}\right. ⎩ ⎨ ⎧Vk(xN)Vk(xk)=LF(xN)=ukmin[Lk(xk,uk)+Vk+1(f(xk,uk))]
根据Bellman方程的等式右侧公式,定义Perturbation如下:
P k ( δ x , δ u ) ≜ L k ( x ~ k + δ x k , u ~ k + δ u k ) − L k ( x ~ k , u ~ k ) + V k + 1 ( f ( x ~ k + δ x k , u ~ k + δ u k ) ) − V k + 1 ( f ( x ~ k , u ~ k ) ) P^k(\delta x, \delta u) \triangleq L^k(\tilde{x}_k + \delta x_k, \tilde{u}_k + \delta u_k) - L^k(\tilde{x}_k, \tilde{u}_k) + V^{k+1}(f(\tilde{x}_k + \delta x_k, \tilde{u}_k + \delta u_k)) - V^{k+1}(f(\tilde{x}_k, \tilde{u}_k)) Pk(δx,δu)≜Lk(x~k+δxk,u~k+δuk)−Lk(x~k,u~k)+Vk+1(f(x~k+δxk,u~k+δuk))−Vk+1(f(x~k,u~k))
其中, P k ( 0 , 0 ) = 0 P^k(0,0) = 0 Pk(0,0)=0,使用二阶泰勒展开:
P k ( δ x , δ u ) ≈ 1 2 [ 1 δ x δ u ] T [ 0 ( P x k ) T ( P u k ) T P x k P x x k P x u k P u k P u x k P u u k ] [ 1 δ x δ u ] = 1 2 ( δ x P x k + δ u P u k + ( P x k ) T δ x + δ x P x x k δ x + δ u P u x k δ x + ( P u k ) T δ u + δ x P x u k δ u + δ u P u u k δ u ) P^k(\delta x, \delta u) \approx \frac{1}{2} \begin{bmatrix} 1 \\ \delta x \\ \delta u \end{bmatrix} ^T \begin{bmatrix} 0 & (P_x^k)^T & (P_u^k)^T \\ P_x^k & P_{xx}^k & P_{xu}^k \\ P_u^k & P_{ux}^k & P_{uu}^k \\ \end{bmatrix} \begin{bmatrix} 1 \\ \delta x \\ \delta u \end{bmatrix} =\frac{1}{2} (\delta x P_x^k + \delta u P_u^k + (P_x^k)^T \delta x + \delta x P_{xx}^k \delta x + \delta u P_{ux}^k \delta x + (P_u^k)^T \delta u + \delta x P_{xu}^k \delta u + \delta u P_{uu}^k \delta u) Pk(δx,δu)≈21 1δxδu T 0PxkPuk(Pxk)TPxxkPuxk(Puk)TPxukPuuk 1δxδu =21(δxPxk+δuPuk+(Pxk)Tδx+δxPxxkδx+δuPuxkδx+(Puk)Tδu+δxPxukδu+δuPuukδu)
∂ P ∂ ( δ u ) = 1 2 ( P u k + P u x k δ x + P u k + P x u k δ x + 2 P u u k δ u ) = P u k + P u x k δ x + P u u k δ u \frac{\partial{P}}{\partial(\delta u)} = \frac{1}{2} (P_u^k + P_{ux}^k \delta x + P_u^k + P_{xu}^k \delta x + 2 P_{uu}^k \delta u) = P_u^k + P_{ux}^k \delta x + P_{uu}^k \delta u ∂(δu)∂P=21(Puk+Puxkδx+Puk+Pxukδx+2Puukδu)=Puk+Puxkδx+Puukδu
P k P^k Pk是标准二次型,因此满足一阶条件下, P k P^k Pk取到最小值,因此令 ∂ P ∂ ( δ u ) = P u k + P u x k δ x + P u u k δ u = 0 \frac{\partial{P}}{\partial(\delta u)} = P_u^k + P_{ux}^k \delta x + P_{uu}^k \delta u = 0 ∂(δu)∂P=Puk+Puxkδx+Puukδu=0,
可得 δ u = − ( P u u k ) − 1 ( P u k + P u x k δ x ) \delta u = -(P_{uu}^k)^{-1} (P_u^k + P_{ux}^k \delta x ) δu=−(Puuk)−1(Puk+Puxkδx)
其中:
{ P x k = L x k + f x T V x k + 1 P u k = L u k + f u T V x k + 1 P x x k = L x x k + f x T V x x k + 1 f x + V x k + 1 f x x P u u k = L u u k + f u T V x x k + 1 f u + V x k + 1 f u u P u x k = L u x k + f u T V x x k + 1 f x + V x k + 1 f u x \begin{cases} \begin{aligned} P_{x}^k &= L_{x}^k + f_{x}^T V_{x}^{k+1} \\ P_{u}^k &= L_{u}^k + f_{u}^T V_{x}^{k+1} \\ P_{xx}^k &= L_{xx}^k + f_{x}^T V_{xx}^{k+1} f_{x} + V_{x}^{k+1} f_{xx} \\ P_{uu}^k &= L_{uu}^k + f_{u}^T V_{xx}^{k+1} f_{u} + V_{x}^{k+1} f_{uu} \\ P_{ux}^k &= L_{ux}^k + f_{u}^T V_{xx}^{k+1} f_{x} + V_{x}^{k+1} f_{ux} \\ \end{aligned} \end{cases} ⎩ ⎨ ⎧PxkPukPxxkPuukPuxk=Lxk+fxTVxk+1=Luk+fuTVxk+1=Lxxk+fxTVxxk+1fx+Vxk+1fxx=Luuk+fuTVxxk+1fu+Vxk+1fuu=Luxk+fuTVxxk+1fx+Vxk+1fux
由于 V k + 1 ( x k + 1 ) = V k + 1 ( f k ( x , u ) ) V^{k+1}(x_{k+1}) = V^{k+1}(f^{k}(x, u)) Vk+1(xk+1)=Vk+1(fk(x,u)),因此链式求导 V V V都是对 x x x求偏导。对于ILQR来说,系统的二阶导数为 0 0 0,即忽略 f x x , f u u , f u x f_{xx}, f_{uu},f_{ux} fxx,fuu,fux。在DDP中,二阶导数不可忽略,这也是ILQR和DDP唯一的区别。
另外,在多数问题下, L u x = 0 L_{ux} = 0 Lux=0。
令 ∇ δ u k P k = 0 \nabla _{\delta u_k} P^k = 0 ∇δukPk=0,可得: δ u k ∗ = − ( P u u k ) − 1 ( P u k + P u x k δ x k ∗ ) = K k δ x k + d k \delta u_k^* = -(P^{k}_{uu})^{-1} (P^k_{u} + P^k_{ux} \delta x_k^*) = K_k \delta x_k + d_k δuk∗=−(Puuk)−1(Puk+Puxkδxk∗)=Kkδxk+dk。
将 δ u k ∗ \delta u_k^* δuk∗的结果带入 P k P^k Pk可得:
min δ u k P k ( δ x , δ u ) = − 1 2 ( P u k ) T ( P u u k ) − 1 P u k + ( ( P x k ) T − ( P u k ) T ( P u u k ) − 1 P u x k ) δ x + 1 2 δ x k T ( ( P x x k ) T − ( P u x k ) T ( P u u k ) − 1 P u x k ) δ x = δ V k \min\limits_{\delta u_k} P^k(\delta x, \delta u) = -\frac{1}{2} (P^k_u)^T (P^k_{uu})^{-1} P^k_u + ((P^k_x)^T - (P^k_u)^T (P^k_{uu})^{-1} P^k_{ux}) \delta x + \frac{1}{2} \delta x^T_k ((P^k_{xx})^T - (P^k_{ux})^T (P^k_{uu})^{-1} P^k_{ux}) \delta x = \delta V^k δukminPk(δx,δu)=−21(Puk)T(Puuk)−1Puk+((Pxk)T−(Puk)T(Puuk)−1Puxk)δx+21δxkT((Pxxk)T−(Puxk)T(Puuk)−1Puxk)δx=δVk
因为: V k ( x ) + δ V k ( x ) = V k ( x + δ x ) ≈ V k ( x k ) + ∂ V k ∂ x ∣ x k ( x − x k ) + 1 2 ( x − x k ) T ∂ 2 V k ∂ x 2 ( x − x k ) V^k(x) + \delta V^k(x) = V^k(x + \delta x) \approx V^k(x_k) + \frac{\partial V^k}{\partial x} |_{x_k} (x-x_k) + \frac{1}{2} (x-x_k)^T \frac{\partial^2 V^k}{\partial x^2} (x-x_k) Vk(x)+δVk(x)=Vk(x+δx)≈Vk(xk)+∂x∂Vk∣xk(x−xk)+21(x−xk)T∂x2∂2Vk(x−xk),可以得到: δ V k ( x k ) ≈ V x k δ x k + 1 2 δ x k T V x k x δ x k \delta V^k(x_k) \approx V^k_x \delta x_k + \frac{1}{2} \delta x^T_k V^k_xx \delta x_k δVk(xk)≈Vxkδxk+21δxkTVxkxδxk,因此:
{ Δ V k = − 1 2 ( P u k ) T ( P u u k ) − 1 P u k = 1 2 d k T P u u d k + d k T P u V x k = ( P x k ) T − ( P u k ) T ( P u u k ) − 1 P u x k = P x + K k T P u u d k + K k T P u + P u x T d k V x x k = ( P x x k ) T − ( P u x k ) T ( P u u k ) − 1 P u x k = P x x + K k T P u u K k + K k T P u x + P u x T K k \begin{cases} \begin{aligned} \Delta V^k &= -\frac{1}{2} (P^k_u)^T (P^k_{uu})^{-1} P^k_u = \frac{1}{2} d^T_k P_{uu} d_k + d^T_k P_u \\ V^k_x &= (P^k_x)^T - (P^k_u)^T (P^k_{uu})^{-1} P^k_{ux} = P_x + K^T_k P_{uu} d_k + K^T_k P_u + P^T_{ux} d_k \\ V^k_{xx} &= (P^k_{xx})^T - (P^k_{ux})^T (P^k_{uu})^{-1} P^k_{ux} = P_{xx} + K^T_k P_{uu} K_k + K^T_k P_{ux} + P^T_{ux} K_k\\ \end{aligned} \end{cases} ⎩ ⎨ ⎧ΔVkVxkVxxk=−21(Puk)T(Puuk)−1Puk=21dkTPuudk+dkTPu=(Pxk)T−(Puk)T(Puuk)−1Puxk=Px+KkTPuudk+KkTPu+PuxTdk=(Pxxk)T−(Puxk)T(Puuk)−1Puxk=Pxx+KkTPuuKk+KkTPux+PuxTKk
V x k , V x x k V^k_x, V^k_{xx} Vxk,Vxxk用在 P k − 1 P^{k-1} Pk−1中。
在Backward Pass中计算得到不是 δ u ∗ \delta u^* δu∗,而是反馈增益和开环增益序列。
1.3 Forward Pass
{ x 0 = x s t a r t u k = u ~ k + d k + K k ( x k − x ~ k ) x k + 1 = f ( x k , u k ) \begin{cases} \begin{aligned} x_0 &= x_{start} \\ u_k &= \tilde{u}_k + d_k + K_k(x_k - \tilde{x}_k) \\ x_{k+1} &= f(x_k, u_k) \\ \end{aligned} \end{cases} ⎩ ⎨ ⎧x0ukxk+1=xstart=u~k+dk+Kk(xk−x~k)=f(xk,uk)
其中, x k , x ~ k x_k, \tilde{x}_k xk,x~k分别是两次迭代计算得到的轨迹。
Line Search
论文[13]提出,和所有的非线性优化问题一样,在梯度下降方向进行线性搜索来保证cost有足够的下降。使用backtracking line search计算参数 α \alpha α应用在反馈中。
z = ( J ( X , U ) − J ( X ~ , U ~ ) ) − Δ V ( α ) z = \frac{(J(X,U) - J(\tilde{X}, \tilde{U}) )}{-\Delta V(\alpha)} z=−ΔV(α)(J(X,U)−J(X~,U~))
其中, Δ V ( α ) = ∑ k = 0 N − 1 1 2 α 2 d k T P u u d k + α d k T P u \Delta V(\alpha) = \sum_{k=0}^{N-1} \frac{1}{2} {\alpha}^2 d^T_k P_{uu} d_k +{\alpha} d^T_k P_u ΔV(α)=∑k=0N−121α2dkTPuudk+αdkTPu。(类似Armijo Conditon,保证每一步迭代充分下降。)
如果 z ∈ [ β 1 , β 2 ] z \in [\beta_1, \beta_2] z∈[β1,β2],一般为 [ 1 e − 4 , 10 ] [1e-4, 10] [1e−4,10]范围内,则认为此次迭代求解结果可以接受。否则,需要更新 α = γ α , 0 < γ < 1 \alpha = \gamma \alpha, 0 < \gamma < 1 α=γα,0<γ<1,一般 γ = 0.5 \gamma = 0.5 γ=0.5,重新进行Forward pass。
1.4 Regularization
论文[13]同时提出,当Line search超出最大迭代次数,或者cost超出最大阈值,应该放弃Forward pass,而增大backward pass中的正则项 ρ \rho ρ。
δ u k ∗ = − ( P u u k + ρ I ) − 1 ( P u k + P u x k δ x k ∗ ) = K k δ x k + d k \delta u_k^* = -(P^{k}_{uu} + \rho I)^{-1} (P^k_{u} + P^k_{ux} \delta x_k^*) = K_k \delta x_k + d_k δuk∗=−(Puuk+ρI)−1(Puk+Puxkδxk∗)=Kkδxk+dk
增大 ρ \rho ρ使 K k K_k Kk更加相似单位矩阵,使梯度下降方向更加相似Newton方向。
当 P u u k P^{k}_{uu} Puuk奇异矩阵时,也可以增大正则项 ρ \rho ρ,在这种情况下,Backward pass需要从头开始计算,当Backward pass成功后减小 ρ \rho ρ,一般 ρ ∈ [ 1.5 , 2.0 ] \rho \in [1.5, 2.0] ρ∈[1.5,2.0]。
2. cilqr
CILQR是将约束转换为cost的方式,然后进行ILQR计算。这里约束形式表示为:
{ g k ( x k , u k ) < 0 g N ( x N ) < 0 \begin{cases} \begin{aligned} g^k(x_k, u_k) &< 0 \\ g^N(x_N) &< 0 \\ \end{aligned} \end{cases} {gk(xk,uk)gN(xN)<0<0
2.1 Barrier function
Barrier function是一种被广泛采用的处理方式。将barrier function加到cost function中。
2.1.1 Exponential function
论文[1]采用了指数形式的barrier function处理约束。
b k = q 1 e q 2 g k ( x k , u k ) b^k = q_1 e^{q_2 g^k(x_k, u_k)} bk=q1eq2gk(xk,uk)
其中, q 1 , q 2 q_1,q_2 q1,q2是控制barrier function的参数,用来近似indicator function。
2.1.2 Logarithmic function
指数形式的barrier function有几个问题:
● 不能保证 hard constraint;
● q 1 , q 2 q_1, q_2 q1,q2难以调节,太小会违反约束,太大会使梯度和Hessian在边界处剧烈变化,造成ill-conditions。
论文[2]提出了对数形式的barrier function:
b ( g ( x , u ) ) = − 1 t l o g ( − g ( x , u ) ) b(g(x,u)) = -\frac{1}{t} log (-g(x, u)) b(g(x,u))=−t1log(−g(x,u))
其中, t > 0 t>0 t>0是调节参数。增大 t t t可以逼近indicator function。随着迭代进行, t t t逐渐增大。
2.1.3 ill-conditions
Penalty method的惩罚权重 σ \sigma σ越来越大,或者Barrier method的惩罚权重 σ \sigma σ越来越小时,无约束最优解越来越逼近有约束最优解。但是当无约束最优越来越逼近有约束最优解时,无约束优化问题的函数越来越病态。Penalty Method的曲率(条件数)从约束外部趋近无穷大,Barrier Method的曲率(条件数)从约束内部趋近无穷大。Hessian的最大奇异值变成无界的了。那么,使用梯度的优化算法会面临收敛越来越慢的情况。使用Sequential策略可以缓解这个问题,但是不能永远解决这个问题。
但是,在工程上,约束条件是有物理意义的,精度要求不需要异常的高,Penalty method和Barrier method可以很快的收敛。
2.1.4 Relaxation
但是,上述方法要求初始轨迹必须满足所有的不等式约束。论文[3]修改了barrier function,可以降低对初始解的要求:
b r e l a x ( g ) = { − 1 t l o g ( − g ) ( g < − ϵ ) k − 1 t k [ ( − g − k ϵ ( k − 1 ) ϵ ) k − 1 ] − 1 t l o g ( ϵ ) ( g ≥ − ϵ ) b_{relax}(g) = \begin{cases} -\frac{1}{t} log(-g) & (g < -\epsilon) \\ \frac{k-1}{tk}[(\frac{-g - k \epsilon}{(k-1)\epsilon})^k - 1] - \frac{1}{t} log(\epsilon) & (g \ge -\epsilon) \end{cases} brelax(g)={−t1log(−g)tkk−1[((k−1)ϵ−g−kϵ)k−1]−t1log(ϵ)(g<−ϵ)(g≥−ϵ)
论文中描述, k = 2 k=2 k=2保证 C 2 C^2 C2连续, ϵ , t \epsilon, t ϵ,t是大于 0 0 0的任意数,选择 t = 5 t=5 t=5, ϵ \epsilon ϵ和约束的类型有关。
论文[7]提出了另一种松弛的对数形式的barrier function。写法不同,其实是一样的。都是用的06年的一篇论文。
2.2 Augmented Lagrangian
使用Barrier function的方法主要有以下几个问题:
- 对约束的惩罚函数有很大的权重,cost function在约束边界处急剧下降,会有weak condition的情况,甚至会产生奇异矩阵;
- Barrier function的方法,是在约束范围内迭代优化,需要轨迹的初始解(nominal trajectory)在约束范围内。现在的一些方法通过修改barrier function可以处理初始解轨迹不满足约束的情况;
- Barrier function无法处理等式约束。虽然,一般的行车轨迹规划没有等式约束,但是个别情况下,仍会出现等式约束;
论文[13]详细描述了AL求解ILQR/DDP的过程,论文[6]使用AL-ilqr方法进行轨迹规划。
也是一种Relaxation方法。
2.2.1 Augmented Lagrangian process
J = f ( x ) + λ T c ( x ) + 1 2 c ( x ) T I μ c ( x ) J = f(x) + \lambda^T c(x) + \frac{1}{2} c(x)^T I_{\mu} c(x) J=f(x)+λTc(x)+21c(x)TIμc(x)
其中, λ \lambda λ是拉格朗日乘子, μ \mu μ是惩罚系数,并且有:
I μ = { 0 if c i ( x ) < 0 ∧ λ i = 0 , i ∈ I μ i otherwise I_{\mu} = \begin{cases} 0 & \text{if } c_i(x) < 0 \land \lambda_i =0, i \in \mathcal{I} \\ \mu_i & \text{otherwise } \end{cases} Iμ={0μiif ci(x)<0∧λi=0,i∈Iotherwise
求解过程为:
- 保持 λ \lambda λ和 μ \mu μ为常量,求解 min x J ( x , λ , μ ) \min\limits_{x} J(x, \lambda, \mu) xminJ(x,λ,μ);
- 更新Lagrange multipliers:
λ i + = { λ i + μ i c i ( x ∗ ) i ∈ E max ( 0 , λ i + μ i c i ( x ∗ ) ) i ∈ I \lambda^+_i = \begin{cases} \lambda_i + \mu_i c_i(x^*) & i \in \mathcal{E} \\ \max (0, \lambda_i + \mu_i c_i(x^*)) & i \in \mathcal{I} \\ \end{cases} λi+={λi+μici(x∗)max(0,λi+μici(x∗))i∈Ei∈I - 更新惩罚系数: μ + = ϕ μ , ϕ > 1 \mu^+ = \phi \mu, \phi > 1 μ+=ϕμ,ϕ>1;
- 检查约束是否满足;
- 判断是否收敛;
其中, E \mathcal{E} E和 I \mathcal{I} I分别是等式约束和不等式约束,一般 2 ≤ ϕ ≤ 10 2 \le \phi \le 10 2≤ϕ≤10。
2.2.2 AL-ilqr
2.2.2.1 Backward Pass
优化目标函数为:
J = L N ( x N ) + ( λ N + 1 2 c N ( x N ) I μ , N ) c N ( x N ) + ∑ k = 0 N − 1 [ L k ( x k , u k ) + ( λ k + 1 2 c k ( x k , u k ) I μ , k ) c k ( x k , u k ) ] J = L_N(x_N) + (\lambda_N + \frac{1}{2} c_N(x_N) I_{\mu,N}) c_N(x_N) + \sum^{N-1}_{k=0}[ L_k(x_k, u_k) + (\lambda_k + \frac{1}{2} c_k(x_k, u_k) I_{\mu,k}) c_k(x_k, u_k)] J=LN(xN)+(λN+21cN(xN)Iμ,N)cN(xN)+∑k=0N−1[Lk(xk,uk)+(λk+21ck(xk,uk)Iμ,k)ck(xk,uk)]
同样的,根据Bellman方程有:
{ V k ( x N ) ∣ λ , μ = L F ( x N ) V k ( x k ) ∣ λ , μ = min u k [ L k ( x k , u k ) + V k + 1 ( f ( x k , u k ) ) ∣ λ , μ ] \left\{\begin{aligned} V^k(x_N) | _{\lambda, \mu} &= L_F(x_N) \\ V^k(x_k) | _{\lambda, \mu} &= \min\limits_{u_k}[L^k(x_k, u_k) + V^{k+1}(f(x_k, u_k))| _{\lambda, \mu}] \end{aligned}\right. ⎩ ⎨ ⎧Vk(xN)∣λ,μVk(xk)∣λ,μ=LF(xN)=ukmin[Lk(xk,uk)+Vk+1(f(xk,uk))∣λ,μ]
同样的,使用泰勒二阶展开近似 V k ( x ) V^k(x) Vk(x)的Perturbation,有:
P k ( δ x , δ u ) ≈ 1 2 [ 1 δ x δ u ] T [ 0 ( P x k ) T ( P u k ) T P x k P x x k P x u k P u k P u x k P u u k ] [ 1 δ x δ u ] P^k(\delta x, \delta u) \approx \frac{1}{2} \begin{bmatrix} 1 \\ \delta x \\ \delta u \end{bmatrix} ^T \begin{bmatrix} 0 & (P_x^k)^T & (P_u^k)^T \\ P_x^k & P_{xx}^k & P_{xu}^k \\ P_u^k & P_{ux}^k & P_{uu}^k \\ \end{bmatrix} \begin{bmatrix} 1 \\ \delta x \\ \delta u \end{bmatrix} Pk(δx,δu)≈21 1δxδu T 0PxkPuk(Pxk)TPxxkPuxk(Puk)TPxukPuuk 1δxδu
其中:
{ P x k = L x k + f x T V x k + 1 + c x T ( λ + I μ c ) P u k = L u k + f u T V x k + 1 + c u T ( λ + I μ c ) P x x k = L x x k + f x T V x x k + 1 f x + V x k + 1 f x x + c x T I μ c x P u u k = L u u k + f u T V x x k + 1 f u + V x k + 1 f u u + c u T I μ c u P u x k = L x x k + f x T V x x k + 1 f x + V x k + 1 f u x + c u T I μ c x \begin{cases} \begin{aligned} P_{x}^k &= L_{x}^k + f_{x}^T V_{x}^{k+1} + c^T_x (\lambda + I_{\mu} c) \\ P_{u}^k &= L_{u}^k + f_{u}^T V_{x}^{k+1} + c^T_u (\lambda + I_{\mu} c) \\ P_{xx}^k &= L_{xx}^k + f_{x}^T V_{xx}^{k+1} f_{x} + V_{x}^{k+1} f_{xx} + c^T_x I_{\mu} c_x \\ P_{uu}^k &= L_{uu}^k + f_{u}^T V_{xx}^{k+1} f_{u} + V_{x}^{k+1} f_{uu} + c^T_u I_{\mu} c_u \\ P_{ux}^k &= L_{xx}^k + f_{x}^T V_{xx}^{k+1} f_{x} + V_{x}^{k+1} f_{ux} + c^T_u I_{\mu} c_x \\ \end{aligned} \end{cases} ⎩ ⎨ ⎧PxkPukPxxkPuukPuxk=Lxk+fxTVxk+1+cxT(λ+Iμc)=Luk+fuTVxk+1+cuT(λ+Iμc)=Lxxk+fxTVxxk+1fx+Vxk+1fxx+cxTIμcx=Luuk+fuTVxxk+1fu+Vxk+1fuu+cuTIμcu=Lxxk+fxTVxxk+1fx+Vxk+1fux+cuTIμcx
可以看出,AL-ilqr和常规的ilqr是相似的,只是在backward pass公式中明确的写出了 c ( x ) c(x) c(x)相关的参数,常规ilqr是隐式表达了。
同样的,可以求出最优控制增量: δ u k ∗ = − ( P u u k + ρ I ) − 1 ( P u k + P u x k δ x k ∗ ) = K k δ x k + d k \delta u_k^* = -(P^{k}_{uu} + \rho I)^{-1} (P^k_{u} + P^k_{ux} \delta x_k^*) = K_k \delta x_k + d_k δuk∗=−(Puuk+ρI)−1(Puk+Puxkδxk∗)=Kkδxk+dk
其中:
{ Δ V k = 1 2 d k T P u u d k + d k T P u V x k = P x + K k T P u u d k + K k T P u + P u x T d k V x x k = P x x + K k T P u u K k + K k T P u s + P u x T K k V x N = L x N + ( c N ) x T ( λ + I μ , N c N ) V x x N = L x x N + ( c N ) x T I μ , N ( c N ) x \begin{cases} \begin{aligned} \Delta V^k &= \frac{1}{2} d^T_k P_{uu} d_k + d^T_k P_u \\ V^k_x &= P_x + K^T_k P_{uu} d_k + K^T_k P_u + P^T_{ux} d_k \\ V^k_{xx} &= P_{xx} + K^T_k P_{uu} K_k + K^T_k P_{us} + P^T_{ux} K_k \\ V^N_x &= L_x^N + (c_N)^T_x (\lambda + I_{\mu,N} c_N) \\ V^N_{xx} &= L_{xx}^N + (c_N)^T_x I_{\mu,N} (c_N)_x \\ \end{aligned} \end{cases} ⎩ ⎨ ⎧ΔVkVxkVxxkVxNVxxN=21dkTPuudk+dkTPu=Px+KkTPuudk+KkTPu+PuxTdk=Pxx+KkTPuuKk+KkTPus+PuxTKk=LxN+(cN)xT(λ+Iμ,NcN)=LxxN+(cN)xTIμ,N(cN)x
2.2.2.2 Forward Pass
和前面所述一致。
2.2.2.3 Augmented Lagrangian update
λ i + = { λ i + μ i c i ( x ∗ ) i ∈ E max ( 0 , λ i + μ i c i ( x ∗ ) ) i ∈ I \lambda^+_i = \begin{cases} \lambda_i + \mu_i c_i(x^*) & i \in \mathcal{E} \\ \max (0, \lambda_i + \mu_i c_i(x^*)) & i \in \mathcal{I} \\ \end{cases} λi+={λi+μici(x∗)max(0,λi+μici(x∗))i∈Ei∈I
μ + = ϕ μ \mu^+ = \phi \mu μ+=ϕμ
论文[13]指出,通过大量的不同的启发式方法更新拉格朗日乘子和惩罚系数对比,在外部循环中更新的效果最好。在一些问题中,只更新拉格朗日乘子,或者更新active constraints对应的惩罚系数,会得到更好的性能,但是基本的方法对大范围的问题可以得到最好的性能。
3. cost function
CILQR的目标函数是ilqr问题的优化函数和不等式约束转化的cost部分的和。
x ∗ , u ∗ = arg min x , u [ ∑ 0 N − 1 ( L k ( x k , u k ) + b ( g k ( x k , u k ) ) ) + L F ( x N ) + b ( g N ( x N ) ) ] s . t . { x k + 1 = f ( x k , u k ) x 0 = x s t a r t \begin{split} & x^* ,u^* = \mathop{\arg \min}\limits_{x, u} [\sum_0^{N-1}(L^k(x_k, u_k) + b(g^k(x_k, u_k)))+ L_F(x_N) + b(g^N(x_N)) ] \\ & s.t. \quad \left\{\begin{aligned} x_{k+1} &= f(x_k, u_k) \\ x_0 &= x_{start} \end{aligned}\right. \end{split} x∗,u∗=x,uargmin[0∑N−1(Lk(xk,uk)+b(gk(xk,uk)))+LF(xN)+b(gN(xN))]s.t.{xk+1x0=f(xk,uk)=xstart
3.1 control cost
控制量的cost,使轨迹更加平滑和平顺。
3.2 tracking cost
对全局求解的粗轨迹的cost。
3.3 adjusting cost
论文[11]提出对上一次迭代的轨迹做cost,避免相邻两次迭代得到轨迹差距太大。
4. constraints
这里的约束是指不等式约束,包括状态变量的约束、控制变量的约束和避障约束。
4.1 状态约束
x l o w ≤ x ≤ x h i g h x_{low} \le x \le x_{high} xlow≤x≤xhigh
4.2 控制约束
u l o w ≤ u ≤ u h i g h u_{low} \le u \le u_{high} ulow≤u≤uhigh
4.3 避障约束
对于自动驾驶轨迹规划问题,避障约束是主要的问题,也有多种处理形式。
4.3.1 Mass point
论文[4,5]将主车和障碍物车辆当作质点,是两个点之间的欧式距离大于安全阈值。
g k ( x k , u k ) = d s a f e − ∣ ∣ p k − p o ∣ ∣ ≤ 0 g^k(x_k, u_k) = d_{safe} - || p_k - p_o || \leq 0 gk(xk,uk)=dsafe−∣∣pk−po∣∣≤0
其中, p k , p o p_k, p_o pk,po分别是主车和障碍物的位置。
4.3.2 Circles Representation
论文[6,7]采用等效圆的方式近似主车和障碍物车辆的形状。不同的是论文[6]采用的两个等效圆近似,论文[7]的研究对象是小型移动机器人,因此采用了一个等效圆。
{ d i − d u i ≥ W e g o / 2 d l i − d i ≥ W e g o / 2 d o b s , j , k i ≥ r o b s + r e g o \begin{cases} d^i - d^i_u \ge W_{ego} / 2 \\ d^i_l - d^i \ge W_{ego} / 2 \\ d^i_{obs,j,k} \ge r_{obs} + r_{ego} \\ \end{cases} ⎩ ⎨ ⎧di−dui≥Wego/2dli−di≥Wego/2dobs,j,ki≥robs+rego
其中, d u , d l d_u, d_l du,dl是道路的上下边界, d i d^i di是主车到道路边界的距离,但是论文中没有解释怎么计算 d i d^i di。 x e g o , f i x^i_{ego, f} xego,fi的位置和 x e g o , r i x^i_{ego, r} xego,ri满足几何约束。
4.3.3 Circles and Ellipse Representation
论文[1,8]将主车近似为两个等半径 r r r的圆形。将障碍物无车辆近似为一个椭圆,并将椭圆按照主车等效圆半径拓展,长轴是其速度方向, a = l + v × t s a f e + s s a f e + r a = l + v \times t_{safe} + s_{safe} + r a=l+v×tsafe+ssafe+r,短轴 b = w + s s a f e + r b = w + s_{safe} + r b=w+ssafe+r。
g = 1 − ( x k − x o ) T R P R T ( x k − x o ) < 0 g = 1 - (x_k - x_o)^T RPR^T (x_k - x_o) < 0 g=1−(xk−xo)TRPRT(xk−xo)<0
其中, P = d i a g ( 1 a , 1 b ) , R = [ cos θ o − sin θ o sin θ o cos θ o ] P = diag(\frac{1}{a}, \frac{1}{b}), R= \begin{bmatrix} \cos{\theta_o} & -\sin{\theta_o} \\ \sin{\theta_o} & \cos{\theta_o} \\ \end{bmatrix} P=diag(a1,b1),R=[cosθosinθo−sinθocosθo]。
这么做优点:
● 在结构化道路场景中,障碍物的圆形近似并不能表达障碍物车辆的纵向运动;
● 椭圆是凸的,平滑的;
● 椭圆符合高斯分布;
4.3.4 Rectangle and Ellipse Representation
论文[9]提出将主车处理为rectangle,将障碍物处理为椭圆,公式描述同4.3.3,并没有介绍主车rectangle的处理方式。
推测为:
- 障碍物的椭圆不需要以车辆等效圆半径膨胀;
- 车辆Rectangle轮廓上的采样点的坐标 p c p_c pc和后轴非线性关系可以得到;
- 约束 p c p_c pc在障碍物的椭圆外面;
4.3.5 Rectangle and Polygon Representation
论文[2,10]提出将主车处理为rectangle,将障碍物车辆处理为Polygon,显然Ploygon比圆形或者椭圆的精度更高。两个凸多边形的距离计算参考论文[11]。
论文[11]中没有描述主车的处理方式,推测和4.3.4中处理方式一样。
论文[11]提出了CFS算法,是一种迭代算法,根据参考轨迹 x k x^k xk计算凸的走廊优化求解得到 x k + 1 x^{k+1} xk+1,然后根据 x k + 1 x^{k+1} xk+1计算凸走廊求解直至收敛。
4.3.5.1 障碍物表达式
O j = { x ∈ R 2 : φ j ( x ) < 0 } , ∂ O j = { x ∈ R 2 : φ j ( x ) = 0 } {\mathcal{O}}_j = \{ x \in \R^2: \varphi_j(x) < 0 \}, \partial {\mathcal{O}}_j = \{ x \in \R^2: \varphi_j(x) =0 \} Oj={x∈R2:φj(x)<0},∂Oj={x∈R2:φj(x)=0}
其中, O j {\mathcal{O}}_j Oj障碍物 j j j的内部空间, ∂ O j \partial {\mathcal{O}}_j ∂Oj表示障碍物 j j j的边界。将障碍物分为三种形式,并给出了各自的表达式。
4.3.5.1.1 Convex obstacle
φ j : = { min y ∈ ∂ O j ∣ ∣ x − y ∣ ∣ x ∉ O j − min y ∈ ∂ O j ∣ ∣ x − y ∣ ∣ x ∈ O j {\varphi}_j := \begin{cases} {\min}_{y \in \partial {\mathcal{O}}_j} || x -y || & x \notin {\mathcal{O}}_j \\ -{\min}_{y \in \partial {\mathcal{O}}_j} || x -y || & x \in {\mathcal{O}}_j \\ \end{cases} φj:={miny∈∂Oj∣∣x−y∣∣−miny∈∂Oj∣∣x−y∣∣x∈/Ojx∈Oj
4.3.5.1.2 Boundary obstacle
φ j : = f ( p 1 ) − p 2 {\varphi}_j := f(p_1) - p_2 φj:=f(p1)−p2
4.3.5.1.3 Non-convex obstacle
φ j : = { min y ∈ ∂ O ^ j [ ∣ ∣ x − y ∣ ∣ + ∣ ∣ y − δ ( y ) ∣ ∣ ] x ∉ O ^ j − min y ∈ ∂ O ^ j [ ∣ ∣ x − y ∣ ∣ − ∣ ∣ y − δ ( y ) ∣ ∣ ] x ∈ O ^ j {\varphi}_j := \begin{cases} {\min}_{y \in \partial {\widehat{\mathcal{O}}}_j} [|| x -y || + || y - \delta(y) ||] & x \notin {\widehat{\mathcal{O}}}_j \\ -{\min}_{y \in \partial {\widehat{\mathcal{O}}}_j} [|| x -y || - || y - \delta(y) ||]& x \in \widehat{{\mathcal{O}}}_j \\ \end{cases} φj:={miny∈∂O j[∣∣x−y∣∣+∣∣y−δ(y)∣∣]−miny∈∂O j[∣∣x−y∣∣−∣∣y−δ(y)∣∣]x∈/O jx∈O j
其中, O ^ j \widehat{\mathcal{O}}_j O j是最小凸包络。
4.3.5.2 避障约束
x ∈ Γ = ⋂ j , q { x ∈ R 2 h : φ j ( l q ( x ) ) ≥ 0 } = x ∈ ⋂ j = 1 N Γ j = x ∈ ⋂ j = 1 N F j x \in \varGamma = \mathop{\bigcap}\limits_{j,q} \{ x \in \R^{2h} : \varphi_j(l_q(x)) \ge 0 \} = x \in \mathop{\bigcap}\limits_{j = 1}^{N} {\varGamma}_j = x \in \mathop{\bigcap}\limits_{j = 1}^{N} {\mathcal{F}}_j x∈Γ=j,q⋂{x∈R2h:φj(lq(x))≥0}=x∈j=1⋂NΓj=x∈j=1⋂NFj
其中, h h h是规划时域,即需要规划的轨迹点的数目。 l q ( x ) : R 2 h → R 2 l_q(x) : \R^{2h} \rightarrow \R^2 lq(x):R2h→R2,即选择规划轨迹上的某个点。
4.3.5.2.1 Convex obstacle
F j ( x r ) : = { x : φ j ( x ) + ∇ ^ φ j ( x r ) ( x − x r ) ≥ 0 } {\mathcal{F}}_j (x^r) := \{ x : \varphi_j(x) + \hat{\nabla} \varphi_j(x^r) (x - x^r) \ge 0 \} Fj(xr):={x:φj(x)+∇^φj(xr)(x−xr)≥0}
其中, ∇ ^ φ j ( x r ) \hat{\nabla} \varphi_j(x^r) ∇^φj(xr)是 φ j \varphi_j φj的梯度方向。
4.3.5.2.2 Non-convex obstalce
F j ( x r ) : = { x : φ j ( x ) + ∇ ^ φ j ( x r ) ( x − x r ) ≥ 1 2 ( x − x r ) R H j ∗ ( x − x r ) } {\mathcal{F}}_j (x^r) := \{ x : \varphi_j(x) + \hat{\nabla} \varphi_j(x^r) (x - x^r) \ge \frac{1}{2} (x - x^r)^R H^*_j (x - x^r) \} Fj(xr):={x:φj(x)+∇^φj(xr)(x−xr)≥21(x−xr)RHj∗(x−xr)}
其中, H j ∗ H^*_j Hj∗是 φ j \varphi_j φj的Hessian下边界。
4.3.5.3 Examples
4.3.5.3.1 Convex obstacle
多边形障碍物的四个顶点分别为 a = ( 1 , 1 ) , b = ( − 1 , 1 ) , c = ( − 1 , − 1 ) , d = ( 1 , − 1 ) a=(1, 1),b=(-1,1),c=(-1,-1),d=(1,-1) a=(1,1),b=(−1,1),c=(−1,−1),d=(1,−1),有:
φ ( x ) = { max { − 1 − p 1 , p 1 − 1 , p 2 − 1 , − 1 − p 2 } ∣ ∣ x ∣ ∣ ∞ ≤ 1 min { ∣ ∣ x − a ∣ ∣ , ∣ ∣ x − b ∣ ∣ , ∣ ∣ x − c ∣ ∣ , ∣ ∣ x − d ∣ ∣ } ∣ p 1 ∣ > 1 , ∣ p 2 ∣ > 1 min { ∣ p 1 − 1 ∣ , ∣ p 1 + 1 ∣ } ∣ p 1 ∣ > 1 , ∣ p 2 ∣ < 1 min { ∣ p 2 − 1 ∣ , ∣ p 2 + 1 ∣ } ∣ p 1 ∣ < 1 , ∣ p 2 ∣ > 1 \varphi(x) = \begin{cases} \max \{ -1 - p_1, p_1 - 1, p_2 - 1, -1 -p_2 \} & ||x||_{\infin} \le 1 \\ \min \{ || x - a||, || x - b||, || x - c||, || x - d|| \} & |p_1| > 1, |p_2| > 1 \\ \min \{ |p_1 - 1|, |p_1 + 1| \} & |p_1| > 1, |p_2| < 1 \\ \min \{ |p_2 - 1|, |p_2 + 1| \} & |p_1| < 1, |p_2| > 1 \\ \end{cases} φ(x)=⎩ ⎨ ⎧max{−1−p1,p1−1,p2−1,−1−p2}min{∣∣x−a∣∣,∣∣x−b∣∣,∣∣x−c∣∣,∣∣x−d∣∣}min{∣p1−1∣,∣p1+1∣}min{∣p2−1∣,∣p2+1∣}∣∣x∣∣∞≤1∣p1∣>1,∣p2∣>1∣p1∣>1,∣p2∣<1∣p1∣<1,∣p2∣>1
参考点 x r = ( − 1.5 , − 1.5 ) x^r = (-1.5, -1.5) xr=(−1.5,−1.5),并且 φ ( x r ) = 2 2 \varphi(x^r) = \frac{\sqrt{2}}{2} φ(xr)=22和 ∇ φ ( x r ) = 2 2 [ − 1 , − 1 ] \nabla \varphi(x^r) = \frac{\sqrt{2}}{2} [-1, -1] ∇φ(xr)=22[−1,−1],因此有:
F ( x r ) : = { x : φ j ( x ) + ∇ ^ φ j ( x r ) ( x − x r ) ≥ 0 } = { x : 2 2 ( − p 1 − p 2 + 2 ) ≥ 0 } \mathcal{F} (x^r) := \{ x : \varphi_j(x) + \hat{\nabla} \varphi_j(x^r) (x - x^r) \ge 0 \} = \{ x : \frac{\sqrt{2}}{2} (-p_1 - p_2 + 2) \ge 0 \} F(xr):={x:φj(x)+∇^φj(xr)(x−xr)≥0}={x:22(−p1−p2+2)≥0}
则,在 x r x^r xr处的凸的可行域如蓝色区域。
4.3.5.3.2 Non-convex obstalce
非凸障碍物的下边界为 p 2 = − ( p 1 ) 2 p_2 = -(p_1)^2 p2=−(p1)2。
φ ( x ) = { max { − ( p 1 ) 2 − p 2 , p 1 − 1 , p 2 − 1 , − 1 − p 2 } ∣ ∣ x ∣ ∣ ∞ ≤ 1 , p 2 ≥ − ( p 1 ) 2 min { ∣ ∣ x − a ∣ ∣ , ∣ ∣ x − b ∣ ∣ , ∣ ∣ x − c ∣ ∣ , ∣ ∣ x − d ∣ ∣ } ∣ p 1 ∣ > 1 , ∣ p 2 ∣ > 1 min { ∣ p 1 − 1 ∣ , ∣ p 1 + 1 ∣ } ∣ p 1 ∣ > 1 , ∣ p 2 ∣ < 1 p 2 − 1 p 2 > 1 − ( p 1 ) 2 − p 2 p 2 < − ( p 1 ) 2 \varphi(x) = \begin{cases} \max \{ -(p_1)^2 - p_2, p_1 - 1, p_2 - 1, -1 -p_2 \} & ||x||_{\infin} \le 1, p_2 \ge -(p_1)^2 \\ \min \{ || x - a||, || x - b||, || x - c||, || x - d|| \} & |p_1| > 1, |p_2| > 1 \\ \min \{ |p_1 - 1|, |p_1 + 1| \} & |p_1| > 1, |p_2| < 1 \\ p_2 - 1 & p_2 > 1 \\ -(p_1)^2 - p_2 & p_2 < -(p_1)^2 \\ \end{cases} φ(x)=⎩ ⎨ ⎧max{−(p1)2−p2,p1−1,p2−1,−1−p2}min{∣∣x−a∣∣,∣∣x−b∣∣,∣∣x−c∣∣,∣∣x−d∣∣}min{∣p1−1∣,∣p1+1∣}p2−1−(p1)2−p2∣∣x∣∣∞≤1,p2≥−(p1)2∣p1∣>1,∣p2∣>1∣p1∣>1,∣p2∣<1p2>1p2<−(p1)2
参考点 x r = ( 0 , − 1 ) x^r = (0, -1) xr=(0,−1),并且 φ ( x r ) = 1 \varphi(x^r) = 1 φ(xr)=1和 ∇ φ ( x r ) = [ 0 , − 1 ] , − H ∗ = [ − 2 , 0 ; 0 , 0 ] \nabla \varphi(x^r) = [0, -1], -H^* = [-2,0;0,0] ∇φ(xr)=[0,−1],−H∗=[−2,0;0,0],因此有:
F ( x r ) : = { x : φ j ( x ) + ∇ ^ φ j ( x r ) ( x − x r ) ≥ 1 2 ( x − x r ) R H j ∗ ( x − x r ) } = { x : − p 2 − ( p 1 ) 2 ≥ 0 } \mathcal{F} (x^r) := \{ x : \varphi_j(x) + \hat{\nabla} \varphi_j(x^r) (x - x^r) \ge \frac{1}{2} (x - x^r)^R H^*_j (x - x^r) \} = \{ x : -p_2 - (p_1)^2 \ge 0 \} F(xr):={x:φj(x)+∇^φj(xr)(x−xr)≥21(x−xr)RHj∗(x−xr)}={x:−p2−(p1)2≥0}
4.3.5.4 凸走廊的可行性
当存在的障碍物都是凸的,并且没有相交时,在 x r x^r xr处必然可以找到凸走廊。反之,则不然。
4.3.6 Minkowski Sum of Polygons
论文[12]基于Minkowski sum的理论(GJK)来处理碰撞约束:
-
将ego的box沿着obstacle的box移动,得到Minkowski sum的多边形 p M S p_{MS} pMS;
-
找到 p M S p_{MS} pMS距离主车最近的顶点 O O O;
-
在顶点 O O O分为五个区域,计算每个区域到ego的距离;
文章没有对计算细节有太多的描述。
4.4 车道线约束
对于结构化道路场景下的自动驾驶轨迹规划问题,车道线约束也是一个重要的约束条件。一般在大曲率道路下,左右车道线形成的行驶区域是非凸的。
4.4.1 polylines
论文[4]分别将左右车道线离散为 20 20 20段长度为 5 m 5m 5m的直线段,将车辆约束在左右对应的直线段内。论文[3]也是类似的处理。
论文[15]同样将车道线近似为 polylines (多条线段)。
4.4.2 polynomial function
论文[11]将左右车道线分别拟合为一条多项式,将轨迹约束在两个多项式之间。
但是,在大曲率,例如环岛,无法使用一条多项式曲线表达车道线。
5. Reference
- Constrained iterative LQR for on-road autonomous driving motion planning
- Autonomous Driving Motion Planning With Constrained Iterative LQR
- Integrating Higher-Order Dynamics and Roadway-Compliance into Constrained ILQR-based Trajectory Planning for Autonomous Vehicles
- Iterative LQR with Discrete Barrier States for Efficient Collision Avoidance of Autonomous Vehicle
- Decentralized iLQR for Cooperative Trajectory Planning of Connected Autonomous Vehicles via Dual Consensus ADMM
- Model Predictive Control of Autonomous Vehicles With Integrated Barriers Using Occupancy Grid Maps
- Adaptive High-Order Control Barrier Function-Based Iterative LQR for Real Time Safety-Critical Motion Planning
- Safe Planning for Self-Driving Via Adaptive Constrained ILQR
- Alternating Direction Method of Multipliers for Constrained Iterative LQR in Autonomous Driving
- Constrained iterative lqg for real-time chance-constrained gaussian belief space planning
- The convex feasible set algo rithm for real time optimization in motion planning
- Real Time Motion Planning Using Constrained Iterative Linear Quadratic Regulator for On-Road Self-Driving
- AL-iLQR Tutorial
- Trajectory Planning for BERTHA - a Local, Continuous Method
- Improved Trajectory Planning for On-Road Self-Driving Vehicles Via Combined Graph Search, Optimization & Topology Analysis