使用IRM法求解一般QCQP问题
参考文献[1]使用秩极小化的迭代法(The iterative rank minimization algorithm, IRM)求解一般二次约束二次规划问题，并将该方法用到无人机路径规划中。本文中涉及的程序(irma.m等程序)源代码都可以在代码托管平台码云上找到，我的用户名：olupengo
QCQP问题
二次约束二次规划问题(Quadratical Constraint Quadratic Programming, QCQP)的目标函数，不等式约束都是二次的，具有如下形式
(1)
J
=
m
i
n
1
2
x
T
Q
0
x
+
a
0
T
x
s
.
t
.
1
2
x
T
Q
j
x
+
a
j
T
x
+
c
j
⩽
0
f
o
r
j
=
1
,
2
,
.
.
.
,
m
A
x
=
b
\begin{aligned} J=\mathrm{min}\quad \frac{1}{2}& x^{T}Q_{0}x+a_{0}^{T}x \\ \mathrm{s.t.}\quad \frac{1}{2} x^{T}Q_{j}x+a_{j}^{T}x &+c_{j}\leqslant0 \quad \mathrm{for} \space j=1,2,...,m \\ Ax&=b \end{aligned} \tag{1}
J=min21s.t.21xTQjx+ajTxAxxTQ0x+a0Tx+cj⩽0for j=1,2,...,m=b(1)
其中，
Q
0
,
.
.
.
,
Q
m
Q_{0},...,Q_{m}
Q0,...,Qm是
n
×
n
n \times n
n×n矩阵，
x
∈
R
n
x \in \mathbb{R}^{n}
x∈Rn是优化变量。如果
Q
0
,
.
.
.
,
Q
m
Q_{0},...,Q_{m}
Q0,...,Qm都是半正定矩阵，这一问题就是凸的。
文献[1]将路径规划问题转化为了非凸QCQP问题，形式如下
(2)
J
=
m
i
n
x
x
T
Q
0
x
+
a
0
T
x
s
.
t
.
x
T
Q
j
x
+
a
j
T
x
+
c
j
⩽
0
∀
j
=
1
,
2
,
.
.
.
,
m
l
x
⩽
x
⩽
u
x
\begin{aligned} J=\mathrm{min}_{x} & x^{T}Q_{0}x+a_{0}^{T}x \\ \mathrm{s.t.}\quad x^{T}Q_{j}x+a_{j}^{T}x &+c_{j}\leqslant0 \quad \forall \space j=1,2,...,m \\ l_{x}\leqslant& x \leqslant u_{x} \end{aligned} \tag{2}
J=minxs.t.xTQjx+ajTxlx⩽xTQ0x+a0Tx+cj⩽0∀ j=1,2,...,mx⩽ux(2)
其中，
x
∈
R
n
x \in \mathbb{R}^{n}
x∈Rn是优化变量，
Q
j
∈
S
n
×
n
,
j
=
1
,
2
,
.
.
.
,
m
Q_{j}\in \mathbb{S^{n\times n}},j=1,2,...,m
Qj∈Sn×n,j=1,2,...,m是对称矩阵，
c
j
∈
R
,
j
=
1
,
.
.
.
,
m
c_{j}\in \mathbb{R},j=1,...,m
cj∈R,j=1,...,m，
a
j
∈
R
n
,
j
=
1
,
.
.
.
,
m
a_{j}\in \mathbb{R}^{n},j=1,...,m
aj∈Rn,j=1,...,m，
l
x
∈
R
n
l_{x}\in \mathbb{R}^{n}
lx∈Rn和
u
x
∈
R
n
u_{x}\in \mathbb{R}^{n}
ux∈Rn分别是
x
x
x的下界和上界。
IRM法介绍
如果遇到像问题(2)那样的非齐次QCQP问题，可以通过引入一个新变量
α
∈
R
\alpha \in \mathbb{R}
α∈R和一个新约束
α
2
=
1
\alpha^{2}=1
α2=1，将非齐次QCQP问题转化为齐次QCQP问题。例如问题(2)可以转化为如下式(3)的齐次形式
(3)
J
=
m
i
n
[
x
T
α
]
[
Q
0
a
0
/
2
a
0
T
/
2
0
]
[
x
α
]
s
.
t
.
[
x
T
α
]
[
Q
j
a
j
/
2
a
j
T
/
2
0
]
[
x
α
]
+
c
j
⩽
0
l
x
⩽
x
⩽
u
x
α
2
=
1
\begin{aligned} J=\mathrm{min}[x^{T}\quad \alpha] \begin{bmatrix} Q_{0} & a_{0}/2 \\ a_{0}^{T}/2 & 0 \end{bmatrix}& \begin{bmatrix} x\\ \alpha \end{bmatrix} \\ \mathrm{s.t.}\quad [x^{T}\quad \alpha] \begin{bmatrix} Q_{j} & a_{j}/2 \\ a_{j}^{T}/2 & 0 \end{bmatrix}& \begin{bmatrix} x\\ \alpha \end{bmatrix}+c_{j}\leqslant 0 \\ l_{x}\leqslant x \leqslant& u_{x} \\ \alpha^{2}=1& \end{aligned} \tag{3}
J=min[xTα][Q0a0T/2a0/20]s.t.[xTα][QjajT/2aj/20]lx⩽x⩽α2=1[xα][xα]+cj⩽0ux(3)
(4)
[
x
T
α
]
[
Q
j
a
j
/
2
a
j
T
/
2
0
]
[
x
α
]
+
c
j
⩽
0
\begin{aligned} [x^{T}\quad \alpha] \begin{bmatrix} Q_{j} & a_{j}/2 \\ a_{j}^{T}/2 & 0 \end{bmatrix}& \begin{bmatrix} x\\ \alpha \end{bmatrix}+c_{j}\leqslant 0 \\ \end{aligned} \tag{4}
[xTα][QjajT/2aj/20][xα]+cj⩽0(4)由(4)式可得
(5)
x
T
Q
j
x
+
α
a
j
T
x
+
c
j
⩽
0
\begin{aligned} x^{T}Q_{j}x+\alpha a_{j}^{T}x+c_{j}\leqslant0 \end{aligned} \tag{5}
xTQjx+αajTx+cj⩽0(5)
将(5)式两边同时除以
α
2
\alpha^{2}
α2可得
(6)
(
x
/
α
)
T
Q
j
(
x
/
α
)
+
a
j
T
(
x
/
α
)
+
c
j
⩽
0
\begin{aligned} (x/\alpha)^{T}Q_{j}(x/\alpha)+a_{j}^{T}(x/\alpha) +c_{j}\leqslant0 \end{aligned} \tag{6}
(x/α)TQj(x/α)+ajT(x/α)+cj⩽0(6)对比(2)式和(6)式可知，如果通过求解问题(3)得到最优解
(
x
∗
,
α
∗
)
(x^{*},\alpha^{*})
(x∗,α∗)，那么
x
∗
/
α
∗
x^{*}/\alpha^{*}
x∗/α∗一定是问题(2)的最优解。
齐次QCQP问题可以写成以下更简洁的形式
(7)
J
=
m
i
n
x
T
Q
0
x
s
.
t
.
x
T
Q
j
x
⩽
c
j
,
∀
j
=
1
,
.
.
.
,
m
\begin{aligned} J=\mathrm{min} &\space x^{T}Q_{0}x\\ \mathrm{s.t.} x^{T}Q_{j}x \leqslant c_{j},& \space \forall j=1,...,m \\ \end{aligned} \tag{7}
J=mins.t.xTQjx⩽cj, xTQ0x ∀j=1,...,m(7)
引入秩1半正定矩阵
X
=
x
x
T
X=xx^{T}
X=xxT，将齐次QCQP问题(7)松弛为式(8)所示半定规划问题(semidefinite programming，SDP)
(8)
J
=
m
i
n
X
⟨
X
,
Q
0
⟩
s
.
t
.
⟨
X
,
Q
j
⟩
⩽
c
j
,
∀
j
=
1
,
.
.
.
,
m
X
≽
0
\begin{aligned} J=\mathrm{min}_{X} & \left \langle X,Q_{0} \right \rangle \\ \mathrm{s.t.} \left \langle X,Q_{j} \right \rangle \leqslant c_{j}, &\space \forall j=1,...,m \\ X \succcurlyeq & 0 \\ \end{aligned} \tag{8}
J=minXs.t.⟨X,Qj⟩⩽cj,X≽⟨X,Q0⟩ ∀j=1,...,m0(8)
此时的问题(8)不能等价于原本的问题，因为问题(8)中的
X
X
X不满足秩1约束。已知，当
X
X
X是非零正定矩阵时，
X
X
X是秩1矩阵的充要条件为
r
I
n
−
1
−
V
T
X
V
≽
0
rI_{n-1}-V^{T} XV \succcurlyeq 0
rIn−1−VTXV≽0，其中
V
∈
R
n
×
(
n
−
1
)
V \in \mathbb{R}^{n \times (n-1)}
V∈Rn×(n−1)是
X
X
X的
n
−
1
n-1
n−1个较小特征值所对应的特征向量组成的矩阵，
r
r
r是趋于0的正数。IRM算法通过迭代的方法，逐渐减小
X
X
X的秩。所以可以把问题(8)转化为下面的凸优化问题
(9)
J
=
m
i
n
X
k
,
r
k
⟨
X
k
,
Q
0
⟩
+
w
k
r
k
s
.
t
.
⟨
X
,
Q
j
⟩
⩽
c
j
,
∀
j
=
1
,
.
.
.
,
m
X
k
≽
0
r
k
I
n
−
1
−
V
k
−
1
T
X
k
V
k
−
1
≽
0
\begin{aligned} J=\mathrm{min}_{X_{k},r_{k}} & \left \langle X_{k}, Q_{0} \right \rangle + w^{k}r_{k} \\ \mathrm{s.t.} \left \langle X,Q_{j} \right \rangle \leqslant c_{j}, &\space \forall j=1,...,m \\ X_{k} \succcurlyeq & 0 \\ r_{k}I_{n-1}-V_{k-1}^{T} & X_{k}V_{k-1} \succcurlyeq 0 \\ \end{aligned} \tag{9}
J=minXk,rks.t.⟨X,Qj⟩⩽cj,Xk≽rkIn−1−Vk−1T⟨Xk,Q0⟩+wkrk ∀j=1,...,m0XkVk−1≽0(9)
其中，
w
>
1
w>1
w>1为
r
k
r_k
rk的权重系数。
IRM法通过求解上面的SDP问题(8)获得
X
0
X_{0}
X0，通过
X
0
X_{0}
X0求得
V
0
V_{0}
V0，然后通过求解问题(9)获得
X
k
X_{k}
Xk，求得
V
k
V_{k}
Vk，不断迭代直到
r
k
r_{k}
rk足够小。IRM法的详细推导过程见参考文献[2]，如果只是想使用这个IRM算法，可以直接看下面的例子。
IRM法使用实例
- 求单位圆上 x x x最小的点的坐标。
(10)
J
=
m
i
n
(
x
)
s
.
t
.
x
2
+
y
2
=
1
\begin{aligned} J=&\mathrm{min}(x) \\ \mathrm{s.t.}\quad x^2&+y^2=1 \end{aligned} \tag{10}
J=s.t.x2min(x)+y2=1(10)使用irma.m子函数求解问题(10)的MATLAB程序untitled1.m已上传至码云。
- 求单位圆上 x 2 x^{2} x2最小的点的坐标，注意：代价函数和约束都是二次就不要加 α 2 = 1 \alpha^2=1 α2=1这个约束了。
(11)
J
=
m
i
n
(
x
2
)
s
.
t
.
x
2
+
y
2
=
1
\begin{aligned} J=&\mathrm{min}(x^{2}) \\ \mathrm{s.t.}\quad &x^2+y^2=1 \end{aligned} \tag{11}
J=s.t.min(x2)x2+y2=1(11)使用irma.m子函数求解问题(11)的MATLAB程序untitled2.m以上传至码云。
3. 求单位圆上距离
(
1
,
1
)
(1,1)
(1,1)最近的点的坐标。
(12)
J
=
m
i
n
(
(
x
−
1
)
2
+
(
y
−
1
)
2
)
s
.
t
.
x
2
+
y
2
=
1
\begin{aligned} J = \mathrm{min}((x-1&)^{2}+(y-1)^{2}) \\ \mathrm{s.t.}\quad x^2+&y^2=1 \end{aligned} \tag{12}
J=min((x−1s.t.x2+)2+(y−1)2)y2=1(12)
使用irma.m求解问题(12)的源代码untitled3.m以上传至码云。
- 求到两个固定端点 ( 0 , 0 ) (0,0) (0,0)和 ( 10 , 10 ) (10,10) (10,10)距离相等，并且到两个固定端点距离相等的点的坐标。
(13)
J
=
m
i
n
(
(
x
2
−
x
1
)
2
+
(
y
2
−
y
1
)
2
)
s
.
t
.
(
x
2
−
x
1
)
2
+
(
y
2
−
y
1
)
2
=
(
x
3
−
x
2
)
2
+
(
y
3
−
y
2
)
2
0
⩽
x
1
⩽
0
0
⩽
y
1
⩽
0
0
⩽
x
2
⩽
10
0
⩽
y
2
⩽
12
10
⩽
x
3
⩽
10
10
⩽
y
3
⩽
10
\begin{aligned} J=\mathrm{min}( (x_2-x_1&)^2+(y_2-y_1)^2 ) \\ \mathrm{s.t.}\quad (x_2-x_1)^2+(y_2-y_1)^2 &= (x_3-x_2)^2+(y_3-y_2)^2 \\ 0 \leqslant x_1 &\leqslant 0 \\ 0 \leqslant y_1 &\leqslant 0 \\ 0 \leqslant x_2 &\leqslant 10 \\ 0 \leqslant y_2 &\leqslant 12 \\ 10 \leqslant x_3 &\leqslant 10 \\ 10 \leqslant y_3 & \leqslant 10 \\ \end{aligned} \tag{13}
J=min((x2−x1s.t.(x2−x1)2+(y2−y1)20⩽x10⩽y10⩽x20⩽y210⩽x310⩽y3)2+(y2−y1)2)=(x3−x2)2+(y3−y2)2⩽0⩽0⩽10⩽12⩽10⩽10(13)
使用irma.m求解问题(13)的源代码untitled4.m以上传至码云。
参考
[1] Sun C , Liu Y C , Dai R , et al. Two Approaches for Path Planning of Unmanned Aerial Vehicles with Avoidance Zones[J]. Journal of Guidance Control & Dynamics, 2017, 40(8).
[2] Sun C , Dai R . An iterative approach to Rank Minimization Problems[C]// Decision & Control. IEEE, 2016.
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
