您的位置:首页 > 汽车 > 时评 > 东莞倣网站_宁波最新消息今天_软文写作经验是什么_今日头条新闻在线看

东莞倣网站_宁波最新消息今天_软文写作经验是什么_今日头条新闻在线看

2024/12/22 22:02:24 来源:https://blog.csdn.net/mpt0816/article/details/144576997  浏览:    关键词:东莞倣网站_宁波最新消息今天_软文写作经验是什么_今日头条新闻在线看
东莞倣网站_宁波最新消息今天_软文写作经验是什么_今日头条新闻在线看

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[0N1Lk(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)+xVkxk(xxk)+21(xxk)Tx22Vk(xxk),可以得到: δ 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} Pk1中。
在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(xkx~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=0N121α2dkTPuudk+αdkTPu。(类似Armijo Conditon,保证每一步迭代充分下降。)
如果 z ∈ [ β 1 , β 2 ] z \in [\beta_1, \beta_2] z[β1,β2],一般为 [ 1 e − 4 , 10 ] [1e-4, 10] [1e4,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)tkk1[((k1)ϵgkϵ)k1]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的方法主要有以下几个问题:

  1. 对约束的惩罚函数有很大的权重,cost function在约束边界处急剧下降,会有weak condition的情况,甚至会产生奇异矩阵;
  2. Barrier function的方法,是在约束范围内迭代优化,需要轨迹的初始解(nominal trajectory)在约束范围内。现在的一些方法通过修改barrier function可以处理初始解轨迹不满足约束的情况;
  3. 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,iIotherwise 
求解过程为:

  1. 保持 λ \lambda λ μ \mu μ为常量,求解 min ⁡ x J ( x , λ , μ ) \min\limits_{x} J(x, \lambda, \mu) xminJ(x,λ,μ)
  2. 更新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))iEiI
  3. 更新惩罚系数: μ + = ϕ μ , ϕ > 1 \mu^+ = \phi \mu, \phi > 1 μ+=ϕμ,ϕ>1
  4. 检查约束是否满足;
  5. 判断是否收敛;
    其中, 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=0N1[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))iEiI
μ + = ϕ μ \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[0N1(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} xlowxxhigh

4.2 控制约束

u l o w ≤ u ≤ u h i g h u_{low} \le u \le u_{high} ulowuuhigh

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∣∣pkpo∣∣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} diduiWego/2dlidiWego/2dobs,j,kirobs+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(xkxo)TRPRT(xkxo)<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θosinθocosθo]
在这里插入图片描述

这么做优点:
● 在结构化道路场景中,障碍物的圆形近似并不能表达障碍物车辆的纵向运动;
● 椭圆是凸的,平滑的;
● 椭圆符合高斯分布;

4.3.4 Rectangle and Ellipse Representation

论文[9]提出将主车处理为rectangle,将障碍物处理为椭圆,公式描述同4.3.3,并没有介绍主车rectangle的处理方式。
推测为:

  1. 障碍物的椭圆不需要以车辆等效圆半径膨胀;
  2. 车辆Rectangle轮廓上的采样点的坐标 p c p_c pc和后轴非线性关系可以得到;
  3. 约束 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={xR2:φj(x)<0},Oj={xR2:φ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:={minyOj∣∣xy∣∣minyOj∣∣xy∣∣x/OjxOj

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:={minyO j[∣∣xy∣∣+∣∣yδ(y)∣∣]minyO j[∣∣xy∣∣∣∣yδ(y)∣∣]x/O jxO 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{xR2h:φj(lq(x))0}=xj=1NΓj=xj=1NFj
其中, h h h是规划时域,即需要规划的轨迹点的数目。 l q ( x ) : R 2 h → R 2 l_q(x) : \R^{2h} \rightarrow \R^2 lq(x):R2hR2,即选择规划轨迹上的某个点。

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)(xxr)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)(xxr)21(xxr)RHj(xxr)}
其中, 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{1p1,p11,p21,1p2}min{∣∣xa∣∣,∣∣xb∣∣,∣∣xc∣∣,∣∣xd∣∣}min{p11∣,p1+1∣}min{p21∣,p2+1∣}∣∣x1p1>1,p2>1p1>1,p2<1p1<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)(xxr)0}={x:22 (p1p2+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)2p2,p11,p21,1p2}min{∣∣xa∣∣,∣∣xb∣∣,∣∣xc∣∣,∣∣xd∣∣}min{p11∣,p1+1∣}p21(p1)2p2∣∣x1,p2(p1)2p1>1,p2>1p1>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)(xxr)21(xxr)RHj(xxr)}={x:p2(p1)20}

4.3.5.4 凸走廊的可行性

当存在的障碍物都是凸的,并且没有相交时,在 x r x^r xr处必然可以找到凸走廊。反之,则不然。

4.3.6 Minkowski Sum of Polygons

论文[12]基于Minkowski sum的理论(GJK)来处理碰撞约束:

  1. 将ego的box沿着obstacle的box移动,得到Minkowski sum的多边形 p M S p_{MS} pMS
    在这里插入图片描述

  2. 找到 p M S p_{MS} pMS距离主车最近的顶点 O O O
    在这里插入图片描述

  3. 在顶点 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

  1. Constrained iterative LQR for on-road autonomous driving motion planning
  2. Autonomous Driving Motion Planning With Constrained Iterative LQR
  3. Integrating Higher-Order Dynamics and Roadway-Compliance into Constrained ILQR-based Trajectory Planning for Autonomous Vehicles
  4. Iterative LQR with Discrete Barrier States for Efficient Collision Avoidance of Autonomous Vehicle
  5. Decentralized iLQR for Cooperative Trajectory Planning of Connected Autonomous Vehicles via Dual Consensus ADMM
  6. Model Predictive Control of Autonomous Vehicles With Integrated Barriers Using Occupancy Grid Maps
  7. Adaptive High-Order Control Barrier Function-Based Iterative LQR for Real Time Safety-Critical Motion Planning
  8. Safe Planning for Self-Driving Via Adaptive Constrained ILQR
  9. Alternating Direction Method of Multipliers for Constrained Iterative LQR in Autonomous Driving
  10. Constrained iterative lqg for real-time chance-constrained gaussian belief space planning
  11. The convex feasible set algo rithm for real time optimization in motion planning
  12. Real Time Motion Planning Using Constrained Iterative Linear Quadratic Regulator for On-Road Self-Driving
  13. AL-iLQR Tutorial
  14. Trajectory Planning for BERTHA - a Local, Continuous Method
  15. Improved Trajectory Planning for On-Road Self-Driving Vehicles Via Combined Graph Search, Optimization & Topology Analysis

版权声明:

本网仅为发布的内容提供存储空间,不对发表、转载的内容提供任何形式的保证。凡本网注明“来源:XXX网络”的作品,均转载自其它媒体,著作权归作者所有,商业转载请联系作者获得授权,非商业转载请注明出处。

我们尊重并感谢每一位作者,均已注明文章来源和作者。如因作品内容、版权或其它问题,请及时与我们联系,联系邮箱:809451989@qq.com,投稿邮箱:809451989@qq.com