二次规划是非线性规划中一种特殊情形,它的目标函数是二次实函数,约束是线性的。由于二次规划比较简单,便于求解,且一些非线性规划可以转化为求解一系列二次规划问题,因此二次规划算法较早引起人们的重视,成为求解非线性规划的一个重要通径。二次规划的算法较多,本章介绍其中几个典型的方法,它们是 Lagrange 方法,起作用集方法,Lemke 方法和路径路踪法。
一、Lagrange 方法
考虑二次规划问题
min
1
2
x
T
H
x
+
c
T
x
s
.
t
.
A
x
=
b
\begin{aligned} &\min\quad\quad \dfrac{1}{2} \pmb x^T\pmb H \pmb x + \pmb c^T \pmb x\\[2ex] &\mathrm{ \ s.t.}\quad\quad\ \ \pmb A\pmb x=\pmb b \end{aligned}
min21xTHx+cTx s.t. Ax=b
其中, H \pmb H H 是 n n n 阶对称矩阵, A \pmb A A 是 m × n m\times n m×n 矩阵, A \pmb A A 的秩为 m m m, b \pmb b b 是 m m m 维列向量。
通过 Lagrange 乘子法求解:首先定义 Language 函数
L
(
x
,
λ
)
=
1
2
x
T
H
x
+
c
T
x
−
λ
T
(
A
x
−
b
)
L(\pmb x,\pmb\lambda) = \frac{1}{2}\pmb x^T\pmb H \pmb x + \pmb c^T \pmb x - \pmb\lambda^T(\pmb A\pmb x-\pmb b)
L(x,λ)=21xTHx+cTx−λT(Ax−b)
令
∇
x
L
(
x
,
λ
)
=
0
,
∇
λ
L
(
x
,
λ
)
=
0
\nabla_xL(\pmb x,\pmb\lambda)=0,\quad\nabla_{\pmb\lambda}L(\pmb x,\pmb\lambda)=0
∇xL(x,λ)=0,∇λL(x,λ)=0
得到方程组
H
x
+
c
−
A
T
λ
=
0
−
A
+
b
=
0
\begin{aligned} &\pmb H\pmb x + \pmb c - \pmb A^T\pmb\lambda=\pmb 0\\[1ex] &-\pmb A\pmb +\pmb b = \pmb 0 \end{aligned}
Hx+c−ATλ=0−A+b=0
将此方程组写成
[
H
−
A
T
−
A
−
0
]
[
x
λ
]
=
[
−
c
−
b
]
\left[ \begin{matrix} \pmb H & - \pmb A^T\\[1ex] -\pmb A & - \pmb 0\\ \end{matrix} \right] \left[ \begin{matrix} \pmb x \\[2ex] \pmb \lambda\\ \end{matrix} \right]= \left[ \begin{matrix} -\pmb c \\[2ex] -\pmb b\\ \end{matrix} \right]
[H−A−AT−0][xλ]=[−c−b]
将系数矩阵称为 Lagrange 矩阵。
设上述 Lagrange 矩阵可逆,且为对称矩阵,则其逆矩阵也为对称矩阵,可表示为
[
H
−
A
T
−
A
−
0
]
−
1
=
[
Q
−
R
T
−
R
−
S
]
\left[ \begin{matrix} \pmb H & - \pmb A^T\\[1ex] -\pmb A & - \pmb 0\\ \end{matrix} \right]^{-1}= \left[ \begin{matrix} \pmb Q & - \pmb R^T\\[1ex] -\pmb R & - \pmb S\\ \end{matrix} \right]
[H−A−AT−0]−1=[Q−R−RT−S]
由可逆矩阵性质,即
[
H
−
A
T
−
A
−
0
]
[
Q
−
R
T
−
R
−
S
]
=
I
m
+
n
\left[ \begin{matrix} \pmb H & - \pmb A^T\\[1ex] -\pmb A & - \pmb 0\\ \end{matrix} \right] \left[ \begin{matrix} \pmb Q & - \pmb R^T\\[1ex] -\pmb R & - \pmb S\\ \end{matrix} \right]=\pmb I_{m+n}
[H−A−AT−0][Q−R−RT−S]=Im+n
推得
H
Q
+
A
T
R
=
I
n
H
R
T
+
(
A
T
S
)
=
0
n
×
m
A
Q
=
0
m
×
n
A
R
T
=
I
m
\begin{aligned} &\pmb{HQ}+\pmb{A^TR}=\pmb I_n\\[1ex] &\pmb{HR^T}+\pmb(A^TS)=\pmb0_{n\times m}\\[1ex] &\pmb{AQ}=\pmb 0_{m\times n}\\[1ex] &\pmb{AR^T} = \pmb I_m \end{aligned}
HQ+ATR=InHRT+(ATS)=0n×mAQ=0m×nART=Im
假设矩阵
H
\pmb H
H 可逆,则可以得到矩阵
Q
,
R
,
S
\pmb{Q,R,S}
Q,R,S 的表达式
Q
=
H
−
1
−
H
−
1
A
T
(
A
H
−
1
A
T
)
−
1
A
H
−
1
,
R
=
(
A
H
−
1
A
T
)
−
1
A
H
−
1
,
S
=
−
(
A
H
−
1
A
T
)
−
1
\begin{aligned} &\pmb{Q}=\pmb H^{-1} - \pmb H^{-1}\pmb A^T(\pmb A\pmb H^{-1}\pmb A^T)^{-1}\pmb A\pmb H^{-1},\\[1ex] &\pmb R = (\pmb A \pmb H^{-1}\pmb A^T)^{-1}\pmb A \pmb H^{-1},\\[1ex] &\pmb S = -(\pmb A \pmb H^{-1}\pmb A^T)^{-1} \end{aligned}
Q=H−1−H−1AT(AH−1AT)−1AH−1,R=(AH−1AT)−1AH−1,S=−(AH−1AT)−1
从而可得问题的解
x
∗
=
−
Q
c
+
R
T
b
λ
∗
=
R
c
−
S
b
\begin{aligned} &\pmb x^\ast = -\pmb{Qc} + \pmb R^T\pmb b\\[1ex] &\pmb \lambda^\ast = \pmb {Rc}-\pmb{Sb} \end{aligned}
x∗=−Qc+RTbλ∗=Rc−Sb
二、起作用集方法
1、起作用集方法的分析推导
考虑具有不等式约束的二次规划问题
min
f
(
x
)
=
1
2
x
T
H
x
+
c
T
x
s
.
t
.
A
x
≥
b
\begin{aligned} &\min\quad\quad f(\pmb x)=\dfrac{1}{2} \pmb x^T\pmb H \pmb x + \pmb c^T \pmb x\\[2ex] &\mathrm{ \ s.t.}\quad\quad\ \ \pmb A\pmb x\geq\pmb b \end{aligned}
minf(x)=21xTHx+cTx s.t. Ax≥b
其中, H \pmb H H 是 n n n 阶对称正定矩阵, A \pmb A A 是 m × n m\times n m×n 矩阵, A \pmb A A 的秩为 m m m, b \pmb b b 是 m m m 维列向量。
由于不等式约束的存在,不能直接用 Lagrange 方法求解,因此需将它转化为求解等式约束问题。运用起作用集方法,在每次追代中,以已知的可行点为起点,把在该点起作用约束作为等式约束,在此约束下极小化目标函数 f ( x ) f(\pmb x) f(x),而其余的约束暂且不管。求得新的比较好的可行点后,再重复以上做法,下面加以具体分析。
设在第
k
k
k 此迭代中,已知可行点
x
(
k
)
\pmb x^{(k)}
x(k),在该点起作用约束指标集用
I
(
k
)
I^{(k)}
I(k) 表示。这时需要求解等式约束
min
f
(
x
)
=
1
2
x
T
H
x
+
c
T
x
s
.
t
.
a
i
x
=
b
i
,
i
∈
I
(
k
)
\begin{aligned} &\min\quad\quad f(\pmb x)=\dfrac{1}{2} \pmb x^T\pmb H \pmb x + \pmb c^T \pmb x\\[2ex] &\mathrm{ \ s.t.}\quad\quad\ \ \pmb a^i\pmb x=b_i,\quad i\in I^{(k)} \end{aligned}
minf(x)=21xTHx+cTx s.t. aix=bi,i∈I(k)
其中 a i \pmb a^i ai 是矩阵 A \pmb A A 的第 i i i 行。
为方便起见,现将坐标原点移至
x
(
k
)
\pmb x^{(k)}
x(k),令
δ
=
x
−
x
(
k
)
\pmb\delta = \pmb x - \pmb x^{(k)}
δ=x−x(k)
则
f
(
x
)
=
1
2
(
δ
+
x
(
k
)
)
T
H
(
δ
+
x
(
k
)
)
+
c
T
(
δ
+
x
(
k
)
)
=
1
2
δ
T
H
δ
+
δ
T
H
x
(
k
)
+
1
2
x
(
k
)
T
H
x
(
k
)
+
c
T
δ
+
c
T
x
(
k
)
=
1
2
δ
T
H
δ
+
∇
f
(
x
(
k
)
)
T
δ
+
f
(
x
(
k
)
)
\begin{aligned} f(\pmb x) &=\dfrac{1}{2} (\pmb\delta + \pmb x^{(k)})^T\pmb H (\pmb\delta + \pmb x^{(k)}) + \pmb c^T (\pmb\delta + \pmb x^{(k)})\\[2ex] &=\dfrac{1}{2}\pmb\delta^T\pmb H \pmb\delta + \pmb\delta ^T\pmb H\pmb x^{(k)}+\frac{1}{2}{\pmb x^{(k)}}^T \pmb H\pmb x^{(k)} +\pmb c^T\pmb\delta +\pmb c^T\pmb x^{(k)} \\[2ex] &=\frac{1}{2}\pmb\delta^T\pmb H \pmb\delta +\nabla f(\pmb x^{(k)})^T\pmb\delta + f(\pmb x^{(k)}) \end{aligned}
f(x)=21(δ+x(k))TH(δ+x(k))+cT(δ+x(k))=21δTHδ+δTHx(k)+21x(k)THx(k)+cTδ+cTx(k)=21δTHδ+∇f(x(k))Tδ+f(x(k))
于是问题转化为求校正量
δ
(
k
)
\pmb\delta^{(k)}
δ(k) 的问题
min
1
2
δ
T
H
δ
+
∇
f
(
x
(
k
)
)
T
δ
s
.
t
.
a
i
δ
=
0
,
i
∈
I
(
k
)
\begin{aligned} &\min\quad\quad \frac{1}{2}\pmb\delta^T\pmb H \pmb\delta +\nabla f(\pmb x^{(k)})^T\pmb\delta\\[2ex] &\mathrm{ \ s.t.}\quad\quad\ \ \pmb a^i\pmb\delta=0,\quad i\in I^{(k)} \end{aligned}
min21δTHδ+∇f(x(k))Tδ s.t. aiδ=0,i∈I(k)
解二次规划,求出最优解 δ ( k ) \pmb\delta^{(k)} δ(k),然后区别不同情形,决定下面应采取的步骤。
- 如果 x ( k ) + δ ( k ) \pmb x^{(k)} + \pmb\delta^{(k)} x(k)+δ(k) 是可行点,且 δ ( k ) ≠ 0 \pmb\delta^{(k)}\neq\pmb0 δ(k)=0,则在第 k + 1 k+1 k+1 次迭代中,已知点取作 x ( k + 1 ) = x ( k ) + δ ( k ) \pmb x^{(k+1)}=\pmb x^{(k)}+\pmb\delta^{(k)} x(k+1)=x(k)+δ(k)。
- 如果
x
(
k
)
+
δ
(
k
)
\pmb x^{(k)} + \pmb\delta^{(k)}
x(k)+δ(k) 不是可行点,则沿方向
d
(
k
)
=
δ
(
k
)
\pmb d^{(k)}=\pmb\delta^{(k)}
d(k)=δ(k) 搜索,令
x ( k + 1 ) = x ( k ) + a k d ( k ) \pmb x^{(k+1)} = \pmb x^{(k)} + a_k\pmb d^{(k)} x(k+1)=x(k)+akd(k)
现在分析怎样确定步长
a
k
a_k
ak,根据保持可行性的要求,其应满足
a
i
(
x
(
k
)
+
a
k
d
(
k
)
)
≥
b
i
,
i
∉
I
(
k
)
\pmb a^i(\pmb x^{(k)} + a_k\pmb d^{(k)})\geq b_i,\quad i\notin I^{(k)}
ai(x(k)+akd(k))≥bi,i∈/I(k)
由于
x
(
k
)
\pmb x^{(k)}
x(k) 是可行点,即
a
i
x
(
k
)
≥
b
i
\pmb a^i\pmb x^{(k)}\geq b_i
aix(k)≥bi,因此
当
a
i
d
(
k
)
≥
0
\pmb a^i\pmb d^{(k)}\geq 0
aid(k)≥0 时,对于任意非负数
a
k
a_k
ak,上式总成立;
当
a
i
d
(
k
)
<
0
\pmb a^i\pmb d^{(k)}< 0
aid(k)<0 时,只要取正数
a
k
≤
min
{
b
i
−
a
i
x
(
k
)
a
i
d
(
k
)
∣
i
∉
I
(
k
)
,
a
i
d
(
k
)
<
0
}
⏟
a
^
k
a_k\leq\underbrace{\min\Bigg\lbrace\frac{b_i-\pmb a^i\pmb x^{(k)}}{\pmb a^i\pmb d^{(k)}}\bigg|i\notin I^{(k)},\ \pmb a^i\pmb d^{(k)}<0\Bigg\rbrace}_{\hat a_k}
ak≤a^k
min{aid(k)bi−aix(k)
i∈/I(k), aid(k)<0}
为了在第
k
k
k 次迭代中得到较好可行点,应取
a
k
=
min
{
1
,
a
^
k
}
,
a_k=\min\lbrace1,\hat a_k\rbrace,
ak=min{1,a^k},
并令
x
(
k
+
1
)
=
x
(
k
)
+
a
k
d
(
k
)
\pmb x^{(k+1)} = \pmb x^{(k)} + a_k\pmb d^{(k)}
x(k+1)=x(k)+akd(k)
如果
a
k
=
b
p
−
a
p
x
(
k
)
a
p
d
(
k
)
<
1
a_k=\frac{b_p-\pmb a^p\pmb x^{(k)}}{\pmb a^p\pmb d^{(k)}}<1
ak=apd(k)bp−apx(k)<1
则在点
x
(
k
+
1
)
\pmb x^{(k+1)}
x(k+1),有
a
p
x
(
k
+
1
)
=
a
p
(
x
(
k
)
+
a
k
d
(
k
)
)
=
b
p
\pmb a^p\pmb x^{(k+1)}=\pmb a^p(\pmb x^{(k)} + a_k\pmb d^{(k)})=b_p
apx(k+1)=ap(x(k)+akd(k))=bp
因此,在 x ( k + 1 ) \pmb x^{(k+1)} x(k+1) 处, a p x ( k ) ≥ b p \pmb a^p\pmb x^{(k)}\geq b_p apx(k)≥bp 为起作用约束。这时,把指标 p p p 加入 I ( k ) I^{(k)} I(k),得到在 x ( k + 1 ) \pmb x^{(k+1)} x(k+1) 处的起作用约束指标集 I ( k + 1 ) I^{(k+1)} I(k+1)。