最优控制小结
本文主要参考刘豹、唐万生主编的《现代控制理论》第三版第6章,非常感谢这些大师将这么复杂的内容讲得这么透彻、通俗易懂。
最优控制做为一种控制理论,通过构造目标函数可以非常直观地对系统的状态、输入量进行控制,在系统状态方程准确的情况下能有非常好的控制效果。本文主要针对工程应用,不进行公式推导、数学证明。
参考练习题.
最优控制小结
两端固定的最优控制问题
问题描述:
J
=
∫
t
0
t
f
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
x
(
t
f
)
=
x
f
(1)
\begin{array}{l} J= \int _{t_0}^{t_f} L[x(t),\dot x(t),t]dt \\ x(t_0) = x_0 \\ x(t_f)=x_{f} \end{array} \tag1
J=∫t0tfL[x(t),x˙(t),t]dtx(t0)=x0x(tf)=xf(1)
极值条件为:
∂
L
∂
x
−
d
d
t
∂
L
∂
x
˙
=
0
x
(
t
0
)
=
x
0
x
(
t
f
)
=
x
f
(2)
\begin{array}{l} \frac{\partial L}{\partial x} - \frac{\mathrm{d} }{\mathrm{d} t} \frac{\partial L}{\partial \dot x} = 0 \\ x(t_0) = x_0 \\ x(t_f) = x_f \end{array} \tag2
∂x∂L−dtd∂x˙∂L=0x(t0)=x0x(tf)=xf(2)
端点自由的最优控制问题
问题描述:
J
=
∫
t
0
t
f
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
(3)
\begin{array}{l} J= \int _{t_0}^{t_f} L[x(t),\dot x(t),t]dt \\ x(t_0) = x_0 \end{array} \tag3
J=∫t0tfL[x(t),x˙(t),t]dtx(t0)=x0(3)
极值条件:
∂
L
∂
x
−
d
d
t
∂
L
∂
x
˙
=
0
x
(
t
0
)
=
x
0
∂
L
∂
x
˙
∣
t
f
=
0
(4)
\begin{array}{l} \frac{\partial L}{\partial x} - \frac{\mathrm{d} }{\mathrm{d} t} \frac{\partial L}{\partial \dot x} = 0 \\ x(t_0) = x_0 \\ \frac{\partial L}{\partial \dot x}|_{t_f} = 0 \end{array} \tag4
∂x∂L−dtd∂x˙∂L=0x(t0)=x0∂x˙∂L∣tf=0(4)
端点可变的最优控制问题
问题描述:
J
=
∫
t
0
t
f
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
x
(
t
f
)
=
C
(
t
f
)
(5)
\begin{array}{l} J= \int _{t_0}^{t_f} L[x(t),\dot x(t),t]dt \\ x(t_0) = x_0 \\ x(t_f) = C(t_f) \end{array} \tag5
J=∫t0tfL[x(t),x˙(t),t]dtx(t0)=x0x(tf)=C(tf)(5)
极值条件:
∂
L
∂
x
−
d
d
t
∂
L
∂
x
˙
=
0
x
(
t
0
)
=
x
0
{
L
−
[
C
˙
(
t
)
−
x
˙
(
t
)
]
∂
L
∂
x
˙
}
∣
t
f
(6)
\begin{array}{l} \frac{\partial L}{\partial x} - \frac{\mathrm{d} }{\mathrm{d} t} \frac{\partial L}{\partial \dot x} = 0 \\ x(t_0) = x_0 \\ \left\{ L-[\dot C(t) -\dot x(t)] \frac{\partial L}{\partial \dot x} \right\}|_{t_f} \end{array} \tag6
∂x∂L−dtd∂x˙∂L=0x(t0)=x0{L−[C˙(t)−x˙(t)]∂x˙∂L}∣tf(6)
终端性能要求的最优控制问题
问题描述:
J
=
ϕ
(
t
f
)
+
∫
t
0
t
f
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
(7)
\begin{array}{l} J= \phi(t_f)+ \int _{t_0}^{t_f} L[x(t),\dot x(t),t]dt \\ x(t_0) = x_0 \end{array} \tag7
J=ϕ(tf)+∫t0tfL[x(t),x˙(t),t]dtx(t0)=x0(7)
极值条件:
∂
L
∂
x
−
d
d
t
∂
L
∂
x
˙
=
0
x
(
t
0
)
=
x
0
∂
L
∂
x
˙
∣
t
f
=
−
∂
ϕ
[
x
(
t
f
)
]
∂
x
(
t
f
)
(8)
\begin{array}{l} \frac{\partial L}{\partial x} - \frac{\mathrm{d} }{\mathrm{d} t} \frac{\partial L}{\partial \dot x} = 0 \\ x(t_0) = x_0 \\ \frac{\partial L}{\partial \dot x}|_{t_f} = - \frac{\partial \phi[x(t_f)]}{\partial x(t_f)} \end{array} \tag8
∂x∂L−dtd∂x˙∂L=0x(t0)=x0∂x˙∂L∣tf=−∂x(tf)∂ϕ[x(tf)](8)
等式约束的最优控制问题
拉格朗日问题
问题描述:
J
=
ϕ
(
t
f
)
+
∫
t
0
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
x
(
t
f
)
=
x
f
x
˙
(
t
)
=
f
[
x
(
t
)
,
u
(
t
)
,
t
]
(9)
\begin{array}{l} J= \phi(t_f) + \int _{t_0}^{t_f} L[x(t),u(t),t]dt \\ x(t_0) = x_0 \\ x(t_f) = x_{f} \\ \dot x(t) = f[x(t),u(t),t] \end{array} \tag9
J=ϕ(tf)+∫t0tfL[x(t),u(t),t]dtx(t0)=x0x(tf)=xfx˙(t)=f[x(t),u(t),t](9)
重新构造目标函数:
J
=
ϕ
(
t
f
)
+
∫
t
0
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
+
λ
T
(
t
)
[
f
[
x
(
t
)
,
u
(
t
)
,
t
]
−
x
˙
(
t
)
]
d
t
J = \phi(t_f) + \int _{t_0}^{t_f} L[x(t),u(t),t] + \lambda ^T (t) [f[x(t),u(t),t] - \dot x(t)] dt
J=ϕ(tf)+∫t0tfL[x(t),u(t),t]+λT(t)[f[x(t),u(t),t]−x˙(t)]dt
构造哈密顿函数:
H
[
x
(
t
)
,
u
(
t
)
,
t
]
=
L
[
x
(
t
)
,
u
(
t
)
,
t
]
+
λ
T
(
t
)
f
[
x
(
t
)
,
u
(
t
)
,
t
]
H[x(t),u(t),t] = L[x(t),u(t),t] + \lambda ^T (t) f[x(t),u(t),t]
H[x(t),u(t),t]=L[x(t),u(t),t]+λT(t)f[x(t),u(t),t]
极值条件:
∂
H
∂
x
=
−
λ
˙
∂
H
∂
λ
=
x
˙
∂
H
∂
u
=
0
x
(
0
)
=
x
0
λ
∣
t
0
t
f
=
0
(10)
\begin{array}{l} \frac{\partial H}{\partial x} = -\dot \lambda \\ \frac{\partial H}{\partial \lambda} = \dot x \\ \frac{\partial H}{\partial u} = 0 \\ x(0) = x_0 \\ \lambda |_{t_0}^{t_f} = 0 \end{array} \tag{10}
∂x∂H=−λ˙∂λ∂H=x˙∂u∂H=0x(0)=x0λ∣t0tf=0(10)
波尔扎问题
问题描述:
J
=
ϕ
(
t
f
)
+
∫
t
0
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
N
[
x
(
t
f
)
,
t
f
]
=
0
x
˙
(
t
)
=
f
[
x
(
t
)
,
u
(
t
)
,
t
]
(11)
\begin{array}{l} J= \phi(t_f) + \int _{t_0}^{t_f} L[x(t),u(t),t]dt \\ x(t_0) = x_0 \\ N[x(t_f), t_f] = 0 \\ \dot x(t) = f[x(t),u(t),t] \end{array} \tag{11}
J=ϕ(tf)+∫t0tfL[x(t),u(t),t]dtx(t0)=x0N[x(tf),tf]=0x˙(t)=f[x(t),u(t),t](11)
重新构造目标函数:
J
=
ϕ
(
t
f
)
+
μ
T
N
[
x
(
t
)
,
u
(
t
)
,
t
]
+
∫
t
0
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
+
λ
T
(
t
)
[
f
[
x
(
t
)
,
u
(
t
)
,
t
]
−
x
˙
(
t
)
]
d
t
J = \phi(t_f) + \mu ^T N[x(t), u(t), t] + \int _{t_0}^{t_f} L[x(t),u(t),t] + \lambda ^T (t) [f[x(t),u(t),t] - \dot x(t)] dt
J=ϕ(tf)+μTN[x(t),u(t),t]+∫t0tfL[x(t),u(t),t]+λT(t)[f[x(t),u(t),t]−x˙(t)]dt
构造哈密顿函数:
H
[
x
(
t
)
,
u
(
t
)
,
t
]
=
L
[
x
(
t
)
,
u
(
t
)
,
t
]
+
λ
T
(
t
)
f
[
x
(
t
)
,
u
(
t
)
,
t
]
H[x(t),u(t),t] = L[x(t),u(t),t] + \lambda ^T (t) f[x(t),u(t),t]
H[x(t),u(t),t]=L[x(t),u(t),t]+λT(t)f[x(t),u(t),t]
极值条件:
∂
H
∂
x
=
−
λ
˙
∂
H
∂
λ
=
x
˙
∂
H
∂
u
=
0
x
(
t
0
)
=
x
0
λ
(
t
f
)
=
∂
ϕ
[
x
(
t
f
)
]
∂
x
(
t
f
)
+
∂
N
T
[
x
(
t
f
)
,
t
f
]
∂
x
(
t
f
)
μ
N
[
x
(
t
f
)
,
t
f
]
=
0
(12)
\begin{array}{l} \frac{\partial H}{\partial x} = -\dot \lambda \\ \frac{\partial H}{\partial \lambda} = \dot x \\ \frac{\partial H}{\partial u} = 0 \\ x(t_0) = x_0 \\ \lambda (t_f) = \frac{\partial \phi [x(t_f)]}{\partial x(t_f)} + \frac{\partial N^T [x(t_f), t_f]}{\partial x(t_f)} \mu \\ N[x(t_f), t_f] = 0 \end{array} \tag{12}
∂x∂H=−λ˙∂λ∂H=x˙∂u∂H=0x(t0)=x0λ(tf)=∂x(tf)∂ϕ[x(tf)]+∂x(tf)∂NT[x(tf),tf]μN[x(tf),tf]=0(12)
不等式约束的最优控制问题
需要用到庞特里亚金的极小值原理进行求解。
问题描述:
J
=
ϕ
(
t
f
)
+
∫
t
0
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
d
t
x
(
t
0
)
=
x
0
x
(
t
f
)
=
x
f
x
˙
(
t
)
=
f
[
x
(
t
)
,
u
(
t
)
,
t
]
g
[
x
(
t
)
,
u
(
t
)
,
t
]
≥
0
(13)
\begin{array}{l} J= \phi(t_f) + \int _{t_0}^{t_f} L[x(t),u(t),t]dt \\ x(t_0) = x_0 \\ x(t_f) = x_{f} \\ \dot x(t) = f[x(t),u(t),t] \\ g[x(t),u(t),t] \geq 0 \end{array} \tag{13}
J=ϕ(tf)+∫t0tfL[x(t),u(t),t]dtx(t0)=x0x(tf)=xfx˙(t)=f[x(t),u(t),t]g[x(t),u(t),t]≥0(13)
重新构造目标函数:
J
=
ϕ
(
t
f
)
+
∫
t
0
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
+
λ
T
(
t
)
[
f
[
x
(
t
)
,
u
(
t
)
,
t
]
−
x
˙
(
t
)
+
μ
T
(
t
)
g
[
x
(
t
)
,
u
(
t
)
,
t
]
]
d
t
J = \phi(t_f) + \int _{t_0}^{t_f} L[x(t),u(t),t] + \lambda ^T (t) [f[x(t),u(t),t] - \dot x(t) + \mu ^T (t) g[x(t),u(t),t]] dt
J=ϕ(tf)+∫t0tfL[x(t),u(t),t]+λT(t)[f[x(t),u(t),t]−x˙(t)+μT(t)g[x(t),u(t),t]]dt
构造哈密顿函数:
H
[
x
(
t
)
,
u
(
t
)
,
t
]
=
L
[
x
(
t
)
,
u
(
t
)
,
t
]
+
λ
T
(
t
)
f
[
x
(
t
)
,
u
(
t
)
,
t
]
H[x(t),u(t),t] = L[x(t),u(t),t] + \lambda ^T (t) f[x(t),u(t),t]
H[x(t),u(t),t]=L[x(t),u(t),t]+λT(t)f[x(t),u(t),t]
极值条件:
∂
H
∂
x
+
∂
g
T
∂
x
γ
=
−
λ
˙
∂
H
∂
λ
=
x
˙
∂
H
∂
u
+
∂
g
T
∂
u
γ
=
0
H
[
x
∗
,
λ
∗
,
u
∗
,
t
]
≤
H
[
x
∗
,
λ
∗
,
u
,
t
]
[
H
+
∂
ϕ
∂
t
f
+
μ
T
∂
N
∂
t
f
]
∣
t
=
t
f
λ
(
t
f
)
=
[
∂
ϕ
∂
x
(
t
f
)
+
∂
N
T
∂
x
(
t
f
)
μ
]
∣
t
f
x
(
0
)
=
x
0
N
[
x
(
t
f
)
,
t
f
]
=
0
(14)
\begin{array}{l} \frac{\partial H}{\partial x} + \frac{\partial g^T}{\partial x} \gamma = -\dot \lambda \\ \frac{\partial H}{\partial \lambda} = \dot x \\ \frac{\partial H}{\partial u} + \frac{\partial g^T}{\partial u} \gamma = 0 \\ H[x^*, \lambda ^*, u^*, t] \leq H[x^*, \lambda ^*, u, t] \\ [H + \frac{\partial \phi}{\partial t_f} + \mu ^T \frac{\partial N}{\partial t_f}]|_{t = t_f} \\ \lambda(t_f) = [\frac{\partial \phi}{\partial x(t_f)} + \frac{\partial N^T}{\partial x(t_f)} \mu]|_{t_f} \\ x(0) = x_0 \\ N[x(t_f), t_f] = 0 \end{array} \tag{14}
∂x∂H+∂x∂gTγ=−λ˙∂λ∂H=x˙∂u∂H+∂u∂gTγ=0H[x∗,λ∗,u∗,t]≤H[x∗,λ∗,u,t][H+∂tf∂ϕ+μT∂tf∂N]∣t=tfλ(tf)=[∂x(tf)∂ϕ+∂x(tf)∂NTμ]∣tfx(0)=x0N[x(tf),tf]=0(14)
连续系统动态规划
贝尔曼方程:
−
∂
J
∗
[
x
,
t
]
∂
t
=
min
u
∈
U
{
L
[
x
,
u
,
t
]
+
[
∂
J
∗
[
x
,
t
]
∂
x
]
T
f
[
x
,
u
,
t
]
}
-\frac{\partial J^*[x, t]}{\partial t} = \min \limits_{u \in U} \left\{ L[x, u, t] + [\frac{\partial J^*[x, t]}{\partial x}]^Tf[x, u, t] \right\}
−∂t∂J∗[x,t]=u∈Umin{L[x,u,t]+[∂x∂J∗[x,t]]Tf[x,u,t]}
解题步骤:
- 构造哈密顿函数:
H [ x , u , t ] = L [ x , u , t ] + ( ∂ J ∗ ∂ x ) T f [ x , u , t ] H[x, u, t] = L[x, u, t] + (\frac{\partial J^*}{\partial x})^T f[x, u, t] H[x,u,t]=L[x,u,t]+(∂x∂J∗)Tf[x,u,t] - 以 H [ x , u , t ] H[x, u, t] H[x,u,t]取极值为条件求 u ∗ u^* u∗;
- 将 u ∗ u^* u∗代入哈密顿-贝尔曼方程,并根据边界条件,解出 J ∗ [ x ( t ) , t ] J^*[x(t), t] J∗[x(t),t]
- 将 J ∗ [ x ( t ) , t ] J^*[x(t), t] J∗[x(t),t]代回 u ∗ u^* u∗,即得最优控制 u ∗ [ x ( t ) , t ] u^*[x(t), t] u∗[x(t),t];
- 将 u ∗ [ x ( t ) , t ] u^*[x(t), t] u∗[x(t),t]代入状态方程,可进一步解出最优轨线 x ∗ ( t ) x^*(t) x∗(t);
- 再将 x ∗ ( t ) x*(t) x∗(t)代入求得最优性能泛函 J ∗ [ x ( t ) ] J*[x(t)] J∗[x(t)]。
线性二次型最优控制
有限时间状态调节器
问题描述:
J
=
1
2
∫
t
0
t
f
[
x
T
Q
1
(
t
)
x
+
u
T
Q
2
(
t
)
u
]
d
t
+
1
2
x
T
(
t
f
)
Q
0
x
(
t
f
)
x
˙
(
t
)
=
A
(
t
)
x
(
t
)
+
B
(
t
)
u
(
t
)
y
=
C
(
t
)
x
(
t
)
x
(
t
0
)
=
x
0
(15)
\begin{array}{l} J=\frac{1}{2} \int_{t_0}^{t_f}[x^TQ_1(t) x + u^TQ_2(t)u] \mathrm{d}t + \frac{1}{2}x^T(t_f)Q_0x(t_f) \\ \dot x(t) = A(t) x(t) + B(t) u(t) \\ y = C(t) x(t) \\ x(t_0) = x_0 \end{array} \tag{15}
J=21∫t0tf[xTQ1(t)x+uTQ2(t)u]dt+21xT(tf)Q0x(tf)x˙(t)=A(t)x(t)+B(t)u(t)y=C(t)x(t)x(t0)=x0(15)
极值条件:
P
˙
(
t
)
=
−
P
(
t
)
A
(
t
)
−
A
T
(
t
)
P
(
t
)
+
P
(
t
)
B
(
t
)
Q
2
−
1
(
t
)
B
T
(
t
)
P
(
t
)
−
Q
1
(
t
)
λ
(
t
)
=
P
(
t
)
x
(
t
)
u
∗
=
−
Q
2
−
1
(
t
)
B
T
(
t
)
λ
(
t
)
J
∗
=
1
2
x
T
(
t
0
)
P
(
t
0
)
x
(
t
0
)
(16)
\begin{array}{l} \dot P(t) = -P(t) A(t) - A^T(t) P(t) + P(t) B(t) Q_2^{-1}(t) B^T(t) P(t) - Q_1(t) \\ \lambda (t) = P(t) x(t) \\ u^* = -Q^{-1}_2(t) B^T(t) \lambda (t) \\ J^* = \frac{1}{2} x^T(t_0) P(t_0) x(t_0) \end{array} \tag{16}
P˙(t)=−P(t)A(t)−AT(t)P(t)+P(t)B(t)Q2−1(t)BT(t)P(t)−Q1(t)λ(t)=P(t)x(t)u∗=−Q2−1(t)BT(t)λ(t)J∗=21xT(t0)P(t0)x(t0)(16)
无限时间状态调节器
问题描述:
J
=
1
2
∫
t
0
∞
[
x
T
Q
1
(
t
)
x
+
u
T
Q
2
(
t
)
u
]
d
t
x
˙
(
t
)
=
A
(
t
)
x
(
t
)
+
B
(
t
)
u
(
t
)
y
=
C
(
t
)
x
(
t
)
x
(
t
0
)
=
x
0
(17)
\begin{array}{l} J=\frac{1}{2} \int_{t_0}^{\infty}[x^TQ_1(t) x + u^TQ_2(t)u] \mathrm{d}t \\ \dot x(t) = A(t) x(t) + B(t) u(t) \\ y = C(t) x(t) \\ x(t_0) = x_0 \end{array} \tag{17}
J=21∫t0∞[xTQ1(t)x+uTQ2(t)u]dtx˙(t)=A(t)x(t)+B(t)u(t)y=C(t)x(t)x(t0)=x0(17)
极值条件:
−
P
A
−
A
T
P
+
P
B
Q
2
−
1
B
T
P
−
Q
1
=
0
u
∗
=
−
Q
2
−
1
B
T
P
x
(
t
)
J
∗
=
1
2
x
T
(
t
0
)
P
x
(
t
0
)
(18)
\begin{array}{l} -P A - A^T P + P BQ_2^{-1} B^T P - Q_1 = 0 \\ u^* = -Q^{-1}_2 B^T P x (t) \\ J^* = \frac{1}{2} x^T(t_0) P x(t_0) \end{array} \tag{18}
−PA−ATP+PBQ2−1BTP−Q1=0u∗=−Q2−1BTPx(t)J∗=21xT(t0)Px(t0)(18)
跟踪控制器问题
问题描述:
e
(
t
)
=
z
(
t
)
−
C
(
t
)
y
(
t
)
J
=
1
2
∫
t
0
t
f
[
e
T
Q
1
(
t
)
e
+
u
T
Q
2
(
t
)
u
]
d
t
+
1
2
e
T
(
t
)
f
)
Q
0
e
(
t
f
)
x
˙
(
t
)
=
A
(
t
)
x
(
t
)
+
B
(
t
)
u
(
t
)
y
=
C
(
t
)
x
(
t
)
x
(
t
0
)
=
x
0
(19)
\begin{array}{l} e(t) = z(t) - C(t) y(t) \\ J=\frac{1}{2} \int_{t_0}^{t_f}[e^TQ_1(t) e + u^TQ_2(t)u] \mathrm{d}t + \frac{1}{2}e^T(t)f) Q_0 e(t_f) \\ \dot x(t) = A(t) x(t) + B(t) u(t) \\ y = C(t) x(t) \\ x(t_0) = x_0 \end{array} \tag{19}
e(t)=z(t)−C(t)y(t)J=21∫t0tf[eTQ1(t)e+uTQ2(t)u]dt+21eT(t)f)Q0e(tf)x˙(t)=A(t)x(t)+B(t)u(t)y=C(t)x(t)x(t0)=x0(19)
极值条件:
P
˙
(
t
)
=
−
P
(
t
)
A
(
t
)
−
A
T
(
t
)
P
(
t
)
+
P
(
t
)
B
(
t
)
Q
2
−
1
(
t
)
B
T
(
t
)
P
(
t
)
−
C
T
(
t
)
Q
1
(
t
)
C
(
t
)
g
˙
(
t
)
=
[
P
(
t
)
B
(
t
)
Q
2
−
1
(
t
)
B
T
(
t
)
−
A
T
(
t
)
]
g
(
t
)
−
C
T
Q
1
(
t
)
z
(
t
)
u
∗
=
−
Q
2
−
1
(
t
)
B
T
(
t
)
P
(
t
)
x
(
t
)
+
Q
2
−
1
(
t
)
B
T
(
t
)
g
(
t
)
(20)
\begin{array}{l} \dot P(t) = -P(t) A(t) - A^T(t) P(t) + P(t) B(t) Q_2^{-1}(t) B^T(t) P(t) - C^T(t) Q_1(t) C(t) \\ \dot g (t) = [P(t) B(t)Q^{-1}_2(t) B^T(t) - A^T(t)] g(t) - C^T Q_1(t)z(t) \\ u^* = -Q^{-1}_2(t) B^T(t) P(t) x (t) + Q^{-1}_2(t)B^T(t)g(t) \end{array} \tag{20}
P˙(t)=−P(t)A(t)−AT(t)P(t)+P(t)B(t)Q2−1(t)BT(t)P(t)−CT(t)Q1(t)C(t)g˙(t)=[P(t)B(t)Q2−1(t)BT(t)−AT(t)]g(t)−CTQ1(t)z(t)u∗=−Q2−1(t)BT(t)P(t)x(t)+Q2−1(t)BT(t)g(t)(20)
更多推荐
所有评论(0)