http://www.ma-xy.com
目录
第一章 约束非线性规划 1
1.1 问题的引入与分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 模型规范化及基本理论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 解的存在性 - 最优性条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 等式约束的最优化条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.2 不等式约束问题的最优性条件 . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.3
一般约束问题的最优性条件
. . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 对偶问题与鞍点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.1 弱对偶定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.2 强对偶定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.3 鞍点定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5 最优化算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.1 外罚函数法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.2 内罚函数法 - 内点法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5.3 乘子法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.4 SQP 方法 (序列二次规划方法) . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 MATLAB 应用实例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1
http://www.ma-xy.com
第一章 约束非线性规划
1.1 问题的引入与分析
前面介绍的是无约束非线性规划，如果我们在求最小化的过程中要求 x 满足一定的约束条
件，则形成约束非线性规划问题。考虑如下含不等式约束的非线性优化问题
min
x
f(x) = x
1
exp(−(x
2
1
+ x
2
2
)) + (x
2
1
+ x
2
2
)/20
s.t. xy/2 + (x + 2)
2
+ (y − 2)
2
/2 ⩽ 2
其中：x = (x
1
, x
2
) ∈ R
2
。
Optimization Toolbox 使用四种算法来求解约束非线性规划问题：¬内点算法：特别适用于
具有稀疏性或其他结构的大规模问题，它基于障碍函数，且在优化迭代过程中对于边界严格可行；
SQP 算法；®动态序列算法；¯信赖域反射算法：仅用于边界约束或线性等。在内点算法和信
赖域反射算法中，对 Hesse 矩阵的近似：
(1) 对内点法而言，可以通过以下方法来近似 Hesse 矩阵：
1) BFGS(稠密)；
2) 有限内存 BFGS(用于大规模问题)；
3) 海赛 - 乘函数；
4) 实际海赛矩阵 (稀疏式稠密)；
5) 有限差分法，但不要求预先知道稀疏性结构。
(2) 对信赖域反射算法，可以通过以下方法来近似 Hesse 矩阵：
1) 有限差分法；
2) 实际海赛矩阵 (稀疏式稠密)；
3) 海赛 - 乘函数；
用 MATLAB 求解上面的问题，程序如下
1 f = @(x , y) x .* exp(−x.^2−y . ^ 2) +(x.^2+y . ^2 ) / 2 0 ;
2 g = @(x , y) x .* y/2+(x+2).^2+(y−2).^2/2−2;
3 ezp l ot (g ,[ −6 ,0 , −1 ,7])
4 hold on
5 ezcontour ( f ,[ −6 ,0 , −1 ,7])
6 plot ( −.9727 ,.4685 , ’ ro ’ ) ;
7 legend ( ’ c on s tr ai nt ’ , ’ f contours ’ , ’minimum ’ ) ;
8 hold o f f
1
http://www.ma-xy.com
1.2 模型规范化及基本理论 第一章 约束非线性规划
9 x0 = [−2 1 ] ;
10 o p t i ons = optimoptions ( ’ fmincon ’ , ’ Algorithm ’ , ’ i n t e r i o r −point ’ , ’ Display ’ , ’ i t e r ’ ) ;% 求 解 器
fmincon 使 用 内 点 算 法 （ i n t e r i o r−point ali g orit h m ） ， 并 打 开 每 一 次 迭 代 的 结 果 。
11 gfun = @(x) deal ( g (x (1) , x( 2 ) ) , [ ] ) ;% 求 解 要 求 非 线 性 约 束 有 两 个 ： 一 个 是 非 线 性 不 等 式 ， 另 一 个 是 非
线 性 等 式 。 我 们 用 deal 函 数 写 这 个 约 束 。
12 [ x , fv al , e x i t f l a g , output ] = fmincon ( fun , x0 , [ ] , [ ] , [ ] , [ ] , [ ] , [ ] , gfun , option s ) ;
13
1.2 模型规范化及基本理论
我们将前面的引例规范化，写出约束非线性规划的一般形式
min
x∈R
n
f(x)
s.t.
h
i
(x) = 0 i = 1, 2, . . . , l
g
j
(x) ⩾ 0 j = 1, 2, . . . , m
其中：h
i
(x), g
j
(x) 为定义在 R
n
上的实值连续函数，h
i
(x) = 0 是等式约束，g
j
(x) ⩾ 0 是不等式
约束，f 是目标函数。如果 m = 0，即不存在不等式约束，则称规划问题为非线性等式约束规划；
如果 h
i
(x), g
j
(x) 为线性函数，则称为线性约束规划。进一步，如果一个线性约束规划的目标函
数 f 为二次函数，则称为二次规划。二次规划是最简单的约束非线性规划问题。
¬量的规定：x ∈ R
n
；f(x) : R
n
→ R；h
i
(x) : R
n
→ R；g
j
(x) : R
n
→ R；至于 f, h, g 的具
体函数类型则有待进一步讨论，例：f ∈ C
1
(D) 或者 f ∈ C
2
(D)。
解空间的定义 (可行域)：满足约束 h
i
(x) = 0 和 g
j
(x) ⩾ 0 的解 x 称为可行解，可行解集
合称为可行域，记为 D = {x|h
i
(x) = 0 & g
j
(x) ⩾ 0}。我们的目标是在 D 中求 x
∗
使得 f 最小。
®解的定义 (全局极小点，局部极小点)：
定义 (全局极小点) 设 x
∗
∈ D，如果对 ∀x ∈ D，有
f(x) ⩾ f(x
∗
)
则称 x
∗
为全局极小点。如果对 ∀x ∈ D, x = x
∗
，有
f(x) > f(x
∗
)
则称 x
∗
为全局严格极小点。
定义 (局部极小点) 设 x
∗
∈ D，如果 ∃δ > 0, ∀x ∈ D ∩ N (x
∗
, δ)，有
f(x) ⩾ f(x
∗
)
则称 x
∗
为 x 的 δ 局部极小点
¬
。其中：N (x
∗
, δ) = {x
∥x−x
∗
∥
2
⩽ δ}。如果对 ∀x ∈ D∩N(x
∗
, δ)/x
∗
，
有
f(x) > f(x
∗
)
则称 x
∗
为严格局部极小点。
¬
注：在全局最优章节之前，我们讨论的 x
∗
皆为局部极小点。
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.3 解的存在性 - 最优性条件
毫无疑问，全局极小值 f(x
∗
) ⩽ 局部极小值f (x
∗
)，某种特殊情况下等式成立，这是有待讨
论的。
1.3 解的存在性 - 最优性条件
假设 f, h
i
, g
j
都是连续可微的。仅依据局部极小点的定义来判断一个 x 是否为 x
∗
是困难的，
因此，我们有必要给出极小点存在的充分必要条件，以方便计算极小点。
我们先给出只有等式约束 h
i
(x) = 0 时解的最优性条件，然后再讨论不等式约束 g
j
(x) ⩾ 0，
最终将二者结合。这样做的理由是：等式约束是相对易于处理的，并且在数学分析中，我们也学
过条件极值和拉格朗日乘数法。
1.3.1 等式约束的最优化条件
考虑如下等式约束规划问题
min
x∈R
n
f(x) (1.1)
s.t. h
i
(x) = 0 i = 1, 2, . . . , l
我们用拉格朗日函数处理上述等式约束极值问题，做 Largrange 函数
L(x, λ) = f(x) −
l
i=1
λ
i
h
i
(x) = λ
T
h(x)
其中：λ = (λ
1
, λ
2
, . . . , λ
l
)
T
为拉格朗日乘子。
一阶必要条件 设 x
∗
是问题 (1.1) 的局部极小点，f, h
i
在 x
∗
的某邻域内连续可微，若向量组
∇h
i
(x
∗
) 线性无关，则 ∃λ
∗
= (λ
∗
1
, λ
∗
2
, . . . , λ
∗
l
)
T
使得
∇
x
L(x
∗
, λ
∗
) = 0
即
∇f(x
∗
) −
l
i=1
λ
∗
i
∇h
i
(x
∗
) = 0
二阶充分条件 定义 L(x, λ) 的梯度以及关于 x 的 Hesse 矩阵
∇L(x, λ) =
∇
x
L(x, λ)
∇
λ
L(x, λ)
=
∇f(x) −
l
i=1
λ
i
∇h
i
(x)
−h(x)
∇
2
xx
L(x, λ) = ∇
2
f(x) −
l
i=1
λ
i
∇
2
h
i
(x) ≜ G
设 f, h
i
是二阶连续可微，并 ∃(x
∗
, λ
∗
) ∈ R
n
×R
l
使 ∇L(x
∗
, λ
∗
) = 0。若 ∀d ∈ R
n
/0，∇h
i
(x
∗
)
T
d = 0，
有 d
T
∇
2
xx
L(x
∗
, λ)d > 0，则 x
∗
是优化问题的一个严格局部极小点。
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.3 解的存在性 - 最优性条件 第一章 约束非线性规划
1.3.2 不等式约束问题的最优性条件
考虑如下不等式约束优化问题
min
x∈R
n
f(x)
s.t. g
j
(x) ⩾ 0 j = 1, 2, . . . , m
记可行域为 D = {x|g
j
(x) ⩾ 0}，指标集 I = {1, 2, . . . , m}。对于某一个可行点 x ∈ D 而言，其
s.t. 会有两种情况：有些约束是 g
j
(x) = 0，另一些约束是 g
j
(x) > 0。对于后一种情形，在 x 的
某个邻域内仍然保持 g
j
(x) > 0 成立，而前者不具备这种性质。因此，我们将二者分开，定义积
极集：
定义 (积极集) 若某个可行点 x ∈ D，使得 g
j
(x) = 0, j ∈ I，则称不等式约束 g
j
(x) ⩾ 0 为
x 的有效约束，称集合 I(x) = {j : g
j
(x) = 0} 为 x 处的有效约束指标集。
给出上述问题的广义拉格朗日函数
L(x, µ) = f(x) − µ
T
g(x) = f (x) −
m
j=1
µ
j
g
j
(x)
其中：µ = (µ
1
, . . . , µ
m
) 为广义拉格朗日乘子。
一阶必要条件 (KKT 条件或 KT 条件) 设 x
∗
为极小点，有效约束指标集 I(x
∗
) = {j|g
j
(x
∗
) =
0, j ∈ I}，假设 f, g
j
在 x
∗
处可微，若向量组 ∇g
j
(x
∗
), j ∈ I(x
∗
) 线性无关，则 ∃µ
∗
= (µ
∗
1
, µ
∗
2
, . . . , mu
∗
m
)
T
使得
∇f(x
∗
) −
m
j=1
µ
∗
i
∇g
j
(x
∗
) = 0
g
j
(x
∗
) ⩾ 0
µ
∗
j
⩾ 0
µ
∗
j
g
j
(x
∗
) = 0
反过来，如果 ∃x
∗
∈ D, µ
∗
∈ R
m
使得上述 KKT 条件成立，则 x
∗
为最优化问题的 K-T 点。
由上述必要性易知，极小点 ⇒K-T 点，但 K-T 点 ⇏ 极小点。
由于 KKT 条件是 1951 年 Kuhn 和 Tucher 给出，故它常被称为 K-T 条件。同时，1939 年
Karush 也类似地考虑了约束的最优性条件，所以也被称为 K-K-T 条件 (Karush-Kuhn-Tucher
定理)。
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.3 解的存在性 - 最优性条件
1.3.3 一般约束问题的最优性条件
上面，我们将等式约束和不等式约束分开讨论，下面来讨论如下一般约束最优问题
min
x∈R
n
f(x)
s.t.
h
i
(x) = 0 i = 1, 2, . . . , l
g
j
(x) ⩾ 0 j = 1, 2, . . . , m
记指标集 E = {1, 2, . . . , l}, I = {1, 2, . . . , m}, i ∈ E, j ∈ I，可行域为 D = {x ∈ R
n
|h
i
(x) =
0, g
j
(x) ⩾ 0, ∀i ∈ E, j ∈ I}。
定义上述优化问题的广义拉格朗日函数为
L(x, λ, µ) = f(x) −
l
i=1
λ
i
h
i
(x) −
m
j=1
µ
j
g
j
(x)
其中：µ, λ 为广义拉格朗日乘子。
一阶必要条件 (KKT 条件) 设 x
∗
是局部极小点，在 x
∗
处的有效约束集为
s(x
∗
) = E ∪ I(x
∗
) = E ∪ {j|g
j
(x
∗
) = 0, i ∈ I}
并假设 f, h
i
, g
j
在 x 处可微。若向量组 ∇h
i
(x
∗
), ∇g
j
(x
∗
), j ∈ I(x
∗
) 线性无关，则存在向量
(λ
∗
, µ
∗
) ∈ R
l
× R
m
，得到
∇f(x
∗
) −
l
i=1
λ
∗
i
∇h
i
(x
∗
) −
m
j=1
µ
∗
j
∇g
j
(x
∗
) = 0
h
i
(x
∗
) = 0, i ∈ E
g
j
(x
∗
) ⩾ 0
µ
∗
j
⩾ 0
µ
∗
j
g
j
(x
∗
) = 0, j ∈ I
其中：µ, λ 为广义拉格朗日乘子。称 µ
∗
j
g
j
(x
∗
) = 0(j ∈ I(x
∗
)) 为互补松弛条件。这意味着 µ
∗
j
, g
j
(x
∗
)
中至少有一个必为 0。
定义广义拉格朗日函数 L 关于 x 的梯度和 Hesse 矩阵为
∇
x
L(x, λ, µ) = ∇f(x) −
l
i=1
λ
i
∇h
i
(x) −
m
j=1
µ
j
∇g
j
(x)
∇
xx
L(x, λ, µ) = ∇
2
f(x) −
l
i=1
λ
i
∇
2
h
i
(x) −
m
j=1
µ
j
∇g
2
j
(x)
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.4 对偶问题与鞍点 第一章 约束非线性规划
二阶充分条件 假设 f, g, h 是二阶连续可微的，设 (x
∗
, (µ
∗
, λ
∗
)) 是优化问题的一个 KKT 点 (即
(x
∗
, (µ
∗
, λ
∗
)) 满足 KKT 条件)，若 ∀d ∈ R
n
/0，∇g
j
(x
∗
)
T
d = 0(j ∈ I(x
∗
))，∇h
i
(x
∗
)
T
d = 0(i ∈ E)
均有 d
T
∇
xx
L(x
∗
, λ
∗
, µ
∗
)d > 0，则 x
∗
是一个严格局部极小点。
由一阶必要条件，我们知道，优化问题的 KKT 点不一定是局部极小点，但如果优化问题是
一个凸优化问题，则 KKT 点 ≡ 局部极小点 ≡ 全局极小点。凸优化是如下优化问题
min
x∈R
n
f(x)
s.t.
h
i
(x) = 0 i = 1, 2, . . . , l
g
j
(x) ⩾ 0 j = 1, 2, . . . , m
如果 f 是凸函数，h
i
(x) 是线性函数 (仿射函数)，g
j
(x) 是凹函数 (即 −g
j
(x) 是凸函数)，那么该
优化问题为凸优化。
1.4 对偶问题与鞍点
考虑如下约束非线性规划问题
min
x∈R
n
f(x)
s.t.
h
i
(x) = 0 i ∈ E
g
j
(x) ⩾ 0 j ∈ I
称上述问题为此非线性规划的原始问题 (PNLP)，相对的对偶问题 (DNLP) 定义如下
max θ(λ, µ) = inf L(x, λ, µ) = inf{f(x) − λ
T
h(x) − µ
T
g(x)|x ∈ D}
s.t. µ ⩾ 0
其中：θ(λ, µ) 称为原问题的拉格朗日对偶函数。此外，如果设原问题和对偶问题的可行域分别为
D 和 ∆，相应的拉格朗日函数为
L(x, λ, µ) = f(x) − λ
T
h(x) − µ
T
g(x), x ∈ D, µ ∈ R
m
+
, λ ∈ R
l
亦可记 m = |I|, l = |E|。值得指出的是，对 ∀x ∈ D，拉格朗日函数 L(x, λ, µ) 是 λ, µ 的线性函
数，于是，拉格朗日对偶函数 θ(λ, µ) 作为线性函数的逐点下确界，必然是一个凹函数。这说明
上述对偶问题是一个凸规划问题。
上面，我们给出了对偶问题的定义，需要明确的是：原始问题的最优值是否等于对偶问题的
最优值呢？即 min
x∈s
f(x) = max
(λ,µ)∈∆
θ(λ, µ)？
1.4.1 弱对偶定理
设 x ∈ D, (λ, µ) ∈ ∆ 分别为原问题和对偶问题的可行解，则
f(x) ⩾ θ(λ, µ)
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.4 对偶问题与鞍点
对于原问题 (PNLP) 和对偶问题 (DNLP)，下面 3 个结论成立：
1) 若 D = ϕ, ∆ = ϕ，则
f
min
= inf{f(x)|x ∈ D} ⩾ θ
max
= sup{θ(λ, µ)|(λ, µ) ∈ ∆}
2) 若 ∃x
∗
∈ D 和 (λ
∗
, µ
∗
) ∈ ∆，使得 f (x
∗
) ⩽ θ(λ
∗
, µ
∗
)，则 x
∗
和 (λ
∗
, µ
∗
) 的解为原问题和
对偶问题的最优解。
3) 若 f
min
= −∞，则 θ(λ, µ) ∈ ∆, θ(λ, µ) = −∞。若 θ
max
= +∞，则原问题没有可行解。
对于一般问题，我们知道 δ = f
min
− θ
max
= 0，下面，我们讨论在什么条件下 δ = 0？
1.4.2 强对偶定理
定义 (slater 约束规格) 对于凸优化来说，若相对内点的集合 riS = {x|h(x) = 0, g(x) >
0, x ∈ D} = ϕ，则称约束函数满足 slater 约束规格。
对于凸优化而言，假设 D 为非空的凸开集，f 为凸函数，g
j
(j ∈ I) 为凹函数，h
i
(i ∈ E) 为
线性函数。若函数 g, h 满足 slater 约束规格，则
f
min
:= inf{f(x)|x ∈ D} = sup{θ(λ, µ)|(λ, µ) ∈ ∆} =: θ
max
此外，若 f
min
> −∞，则 ∃(λ
∗
, µ
∗
) ∈ ∆，使得 θ(λ
∗
, µ
∗
) = θ
max
。特别地，若存在 x
∗
∈ s，使
f(x
∗
) = f
min
，则互补松弛条件 λ
∗T
g(x
∗
) = 0 成立。
1.4.3 鞍点定理
下面介绍鞍点的定义以及鞍点与 KKT 点、最优解之间的关系。
定义 (鞍点) 设 L(x, λ, µ) 是原问题的拉格朗日函数，如果 ∃(x
∗
, λ
∗
, µ
∗
) ∈ D × R
|E|
× R
|I|
+
，
使得对 ∀(x, λ, µ)，有
L(x
∗
, λ, µ) ⩽ L(x
∗
, λ
∗
, µ
∗
) ⩽ L(x, λ
∗
, µ
∗
)
则称 (x
∗
, λ
∗
, µ
∗
) 为函数 L(x, λ, µ) 的一个鞍点 (saddle point)。
鞍点与最优点
在一般情况下，鞍点的存在性并不是最优解存在的必要条件，但是，如果鞍点存在的话，它
一定是约束优化问题原问题的最优解。
1) 设 (x
∗
, λ
∗
, µ
∗
) 是原问题的拉格朗日函数 L(x, λ, µ) 的鞍点，则 x
∗
和 (λ
∗
, µ
∗
) 分别是原问
题和对偶问题的最优解。
2) 如果原问题为凸优化问题，x
∗
为最优解，并且 g
j
(x), h
i
(x) 满足 slater 约束规格，则
∃λ
∗
, µ
∗
⩾ 0，使得 (x
∗
, λ
∗
, µ
∗
) 是原问题的拉格朗日函数的鞍点。
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.5 最优化算法 第一章 约束非线性规划
鞍点与 KKT 点
假设原问题 (PNLP) 中的 f, g, h 在包含可行域 D 的开集 D 内连续可微，那么：
1) 若 (x
∗
, λ
∗
, µ
∗
) ∈ D×R
|I|
+
×R
|E|
是原问题的拉格朗日函数 L(x, λ, µ) 的鞍点，则 (x
∗
, λ
∗
, µ
∗
)
是原问题 (PNLP) 的一个 KKT 点对。
2) 若原问题是一个凸规划，并且 (x
∗
, λ
∗
, µ
∗
) 是原问题的一个 KKT 点，则 (x
∗
, λ
∗
, µ
∗
) 是原
问题的拉格朗日函数的鞍点。
1.5 最优化算法
1.5.1 外罚函数法
罚函数的基本思想：根据约束条件的特点，将其转化为某种罚函数加到目标函数中去，从而
将约束优化问题转化为无约束优化问题来求解。下面，我们来介绍一种罚函数方法 - 外罚函数法，
也称外点法。
考虑如下一般的约束优化问题
min
x∈R
n
f(x)
s.t.
h
i
(x) = 0 i ∈ E
g
j
(x) ⩾ 0 j ∈ I
记可行域为 D = {x|h
i
(x) = 0, g
j
(x) ⩾ 0}, |E| = l, |I| = m，构造罚函数
¯p(x) =
l
i=1
h
2
i
(x) +
m
j=1
[min{0, g
j
(x)}]
2
将罚函数 ¯p(x) 增加到目标中，构建新的目标函数
p(x, σ) = f (x) + σ ¯p(x)
其中：σ 为罚权重，σ > 0。不难发现，当 x ∈ D 时，min f(x) ≜ min p(x, σ)；当 x /∈ D 时，
p(x, σ) > f (x)，σ 越大，p 越大，当 σ 足够大时，要使 p 达到极小，罚函数 ¯p(x) 要充分小才可
以，从而 p(x, σ) 的极小点充分逼近可行域 D，因此，最优化问题等价于
min
x∈R
n
p(x|σ
k
) = f(x) + σ
k
¯p(x)
其中：σ
k
为第 k 次迭代的罚权重，σ
k
> 0，且 σ
k
→ +∞。
由上述思想可知：x(σ) 是从可行域 D 的外部来趋近于 x
∗
，因此上述罚函数也称为外罚函
数/外点法。(问：如何确定初点 x
0
/∈ D)。上述思路有一些缺点：¬当 σ
k
较大时，p(x, σ
k
) 的
Hesse 矩阵的条件数很大，在数值上求解不易。¯p(x) 一般不可微，因此，不易于直接采用导数
求解方法。
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.5 最优化算法
1.5.2 内罚函数法 - 内点法
(1) 对于只有不等式约束的优化问题
min
x∈R
n
f(x)
s.t. g
j
(x) ⩾ 0 j ∈ I
内点法的基本思想是：保持每一个迭代点 x
k
都在可行域 D 内，可行域的边界被筑起一道很高的
“围墙”作为障碍。当迭代点 x
k
靠近边界时，增广目标函数值骤然增大，以示“惩罚”，并阻止
迭代点窜越边界。因此，内点法也称为内罚函数法或障碍函数法。它只适用于可行域内点集非空
的情形。
类似于外罚函数法，我们构造增广目标函数
H(x, τ ) = f (x) + τ
¯
H(x)
其中；
¯
H(x) 为障碍函数。当 x ∈ D，至少有一个 g
j
(x) 趋近于 0，而
¯
H(x) 趋近于无穷大。因此，
可取约束函数倒数之和为障碍函数
¯
H(x) =
m
j=1
1
g
j
(
x
)
或者反对数障碍函数
¯
H(x) = −
m
j=1
ln(g
j
(x))
τ > 0 为罚权重，τ
k
→ 0(k → ∞)。于是，约束优化问题转化为无约束优化问题
min
x
H(x, τ ) = f (x) + τ
¯
H(x)
问：如何确定初始点 x ∈ D？
(2) 如果约束中存在等式约束，我们可以将外罚函数法和内罚函数法相结合，构建增广目标
函数
H(x, µ) = f (x) +
1
2µ
l
i=1
h
2
i
(x) + µ
m
j=1
1
g
j
(x)
或者
H(x, µ) = f (x) +
1
2µ
l
i=1
h
2
i
(x) − µ
m
j=1
ln[g
j
(x)]
另外，我们还可以引入松弛变量 ξ
j
, j = 1, 2, . . . , m，将原问题转化为
min f(x)
s.t.
h
i
(x) = 0 i = 1, 2, . . . , l
g
j
(x) − ξ
j
= 0 j = 1, 2, . . . , m
ξ
j
j = 1, 2, . . . , m
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.5 最优化算法 第一章 约束非线性规划
然后构建混合增广目标函数
φ(x, ξ, µ) = f (x) +
1
2µ
l
i=1
h
2
i
(x) +
1
2µ
m
j=1
[g
j
(x) − ξ
j
]
2
+ µ
m
j=1
1
ξ
j
1.5.3 乘子法
乘子法是 Povell 和 Hestenes 于 1969 年针对等式约束问题同时独立提出的。Rockfellar 于
1973 年将该方法推广到不等式约束优化中。其基本思想是：从原问题的拉格朗日函数出发，再加
上适应的罚函数，从而将原问题转化为一个无约束优化问题。
先只考虑等式约束优化问题
min f(x)
s.t. h
i
(x) = 0 i ∈ E
作上述问题的拉格朗日函数
L(x, λ) = f(x) − λ
T
h(x)
其中：λ
T
= (λ
1
, λ
2
, . . . , λ
l
)
T
为 Lagrange 乘子向量。h(x) = (h
1
(x), h
2
(x), . . . , h
l
(x))
T
。
设 (x
∗
, λ
∗
) 是原问题的 KKT 点，则由最优性条件，有
∇
x
L(x
∗
, λ
∗
) = 0
∇
λ
L(x
∗
, λ
∗
) = −h(x
∗
) = 0
此外，对 ∀x ∈ D，有
L(x
∗
, λ
∗
) = f(x
∗
) ⩽ f(x) − (λ
∗
)
T
h(x) = L(x, λ
∗
)
上式表明，如果已知 λ
∗
，则原问题等价于
min L(x, λ
∗
)
s.t. h(x) = 0
可以考虑用外罚函数法求解上述问题，其增广目标函数为
φ(x, λ
∗
, σ) = L(x, λ
∗
) +
σ
2
∥h(x)∥
2
但 λ
∗
事先并不知道，故可以考虑如下增广目标函数
φ(x, λ, σ) = L(x, λ) +
σ
2
∥h(x)∥
2
= f(x) − λ
T
h(x) +
σ
2
∥h(x)∥
2
我们求 min φ(x, λ, σ)。首先固定一个 λ =
¯
λ，求 φ(x,
¯
λ, σ) 的极小点 ¯x；然后再适当改变 λ 的取
值，再求新的 ¯x，直到求得满意的 x
∗
, λ
∗
为止。
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.5 最优化算法
具体而言，我们在第 k 次迭代求无约束的问题 min φ(x, λ
k
, σ) 的极小点 x
k
时，其取极值的
必要条件为
∇
x
φ(x
k
, λ
k
, σ) = ∇f (x
k
) − ∇h(x
k
)[λ
k
− σh(x
k
)] = 0
而在原问题的 KKT 点 (x
∗
, λ
∗
) 处，有
∇f(x
∗
) − ∇h(x
∗
)λ
∗
= 0 h(x
∗
) = 0
我们自然是希望 {x
k
} → x
∗
, {λ
k
} → λ
∗
，于是，可以取 λ
∗
的更新公式为
λ
k+1
=
λ
k
−
σh
(
x
k
)
且 {h(x
k
)} → 0 是 {λ
k
} 收敛的充要条件，也是 (x
k
, λ
k
) 为 KKT 对的充要条件。
上面讨论的是只含有等式约束的优化问题，如果含有不等式约束，我们可以引入松弛变量 ξ
j
以进行转化。一般约束最优化的增广 Lagrange 函数为
L(x, λ, µ) = f(x) +
1
2µ
j∈I
min
2
{µg
j
(x) − λ
j
, 0} − λ
2
j
−
i∈E
λ
i
h
i
(x) +
1
2
µ
i∈E
h
2
i
(x)
相应的 Lagrange 乘子迭代格式为
λ
+
i
=
λ
i
− µh
i
(x) i ∈ E
max{λ
j
− µg
j
(x), 0} j ∈ I
其中：x, µ, λ 为当前迭代值，λ
+
i
表示下一次迭代值。
给出一般的乘子法算法流程：
step1. 初始化。x
0
∈ R
n
，初始乘子向量 λ
0
，罚参数序列 {µ
k
}，容许误差 ε > 0, k := 0。
step2. 构建增广 Lagrange 函数 L(x, λ, µ)。
step3. 以 x
k−1
作为初始点 (k = 0 时，初始点任意)，求无约束优化
min
x∈R
n
L(x, µ
k
, λ
k
)
解得 x
k
。
step4. 若
∥h(x
k
)∥ + ∥min{g(x
k
), µ
−1
k
λ
k
}∥ ⩽ ε
则解得 x
k
，否则转到 step5。
step5. 更新 λ
k
。令 x := x
k
, λ := λ
k
，由 λ
+
i
更新公式得到 λ
k+1
，转至 step2。
1.5.4 SQP 方法 (序列二次规划方法)
SQP(sequential quadratic programing) 是求解约束优化问题最有效的算法之一。其基本思
想是：在每一次迭代步通过求解一个二次规划子问题来确定一个下降方向。然后，以减少价值函
数来取得步长，重复这些步骤直到收敛。
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.5 最优化算法 第一章 约束非线性规划
Newton - Lagrange 方法
先来考虑等式约束
min
x∈R
n
f(x)
s.t. h(x) = 0
其中：f : R
n
→ R, h
i
: R
n
→ R, h(x) = (h
1
(x), h
2
(x), . . . , h
l
(x))
T
, E = {1, 2, . . . , l}, |E| = l。设
f, h
i
为二阶连续可微函数。
写出原问题的 Lagrange 函数
L(x, µ) = f(x) − µ
T
h(x)
记约束 h 的梯度矩阵为
∇h(x) = (∇h
1
(x), ∇h
2
(x), . . . , ∇h
l
(x))
h 的 Jacobi 矩阵为 A(x) = ∇h(x)
T
。根据原问题的 KKT 条件，可得到如下方程组
∇L(x, µ) =
∇
x
L(x, µ)
∇
µ
L(x, µ)
=
∇f(x) − A(x)
T
µ
−h(x)
= 0 (1.2)
若 A(x
∗
) 行满秩，则原问题的最优解 (x
∗
, µ
∗
) 必然满足上述非线性方程 (1.2)。由于 KKT 条
件为 Lagrange 函数平稳点的条件，所以人们通常把基于求解上述非线性方程的优化方法称为
Lagrange 方法。特别地，如果使用 Newton 方法求解上述方程组，那么相应的优化方法称为
Newton-Lagrange 方法。
现在用牛顿法求解上述的非线性方程组 (1.2)。记函数 ∇L(x, µ) 的 Jacobi 矩阵为
N(x, µ) =
G(x, µ) −A(x)
T
−A(x) 0
其中：
G(x, µ) = ∇
2
xx
L(x, µ) = ∇
2
f(x) −
l
i=1
µ
i
∇
2
h
i
(x)
为 Lagrange 函数 L(x, µ) 关于 x 的 Hesse 矩阵。
对于给定的点 (x
∗
, µ
∗
)，牛顿法的迭代公式为
x
k+1
µ
k+1
=
x
k
k
µ
k
k
+ p
k
k
其中：p
k
k
= (p
k
x
, p
k
µ
)
T
为牛顿方向，满足方程
N(x
k
, µ
k
)p
k
= −∇L(x
k
, µ
k
)
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.5 最优化算法
即
G
k
−A
T
k
−A
k
0
p
k
x
p
k
µ
=
−∇f(x
k
) + A(x
k
)
T
µ
k
h(x
k
)
其中：G
k
= G(x
k
, µ
k
)。如果下面条件成立，那么上述方程的系数矩阵非奇异，并且该方程有唯
一解。
假设约束函数 h 在 x
k
的 Jacobi 矩阵 A
k
行满秩；假设矩阵 g
k
在约束函数 h 的切空间
N(A
k
) 上是正定的，即 ∀d ∈ N (A
k
)/{0}, d
T
w
k
d > 0，则 N-L 方法具有局部二次收敛性质。但由
于每一次迭代均求解非线性方程组，导致数值上的不稳定。鉴于这种不稳定性，所以将其转化为
一个严格凸二次规划问题。转化的条件是原问题的解点 x
∗
处最优化二阶充分条件成立，即对满
足 A(x
∗
)
T
d = 0 的任一非零向量 d，有
d
T
G(x
∗
, µ
∗
)d > 0
这时，当 τ > 0 充分小时，有
G
(
x
∗
, µ
∗
) +
1
2τ
A
(
x
∗
)
T
A
(
x
∗
)
正定。考虑将方程组 (1.2) 中的 G(x
k
, µ
k
) 用一个正定矩阵来代替，记
B(x
k
, µ
k
) = G(x
∗
, µ
∗
) +
1
2τ
A(x
k
)
T
A(x
k
)
则当 (x
k
, µ
k
) → (x
∗
, µ
∗
) 时，矩阵 B(x
∗
, µ
∗
) 正定。注意到 (1.2) 方程组的展开式为
G(x
k
, µ
k
)d
k
− A(x
k
)
T
v
k
= −∇f(x
k
) + A(x
k
)
T
µ
k
将上式变形为
[G(x
k
, µ
k
) +
1
2τ
A(x
k
)
T
A(x
k
)]d
k
− A(x
k
)
T
µ
k
+ v
k
+
1
2τ
A(x
k
)d
k
= −∇f(x
k
)
令
¯µ
k
:= µ
k
+ v
k
+
1
2τ
A(x
k
)d
k
即得
B(x
k
, µ
k
)d
k
− A(x
k
)
T
¯µ
k
= −∇f(x
k
)
因此，方程组 (1.2) 等价于
B(x
k
, µ
k
) −A(x
k
)
T
−A(x
k
) 0
d
k
¯µ
k
= −
∇f(x
k
)
h(x
k
)
进一步，可以将上述方程转化为严格凸二次规划
min
d
q
k
(d) =
1
2
d
T
B(x
k
, µ
k
)d + ∇f(x
k
)
T
d (1.3)
s.t. h(x
k
) + A(x
k
)d = 0
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.5 最优化算法 第一章 约束非线性规划
其中：B(x
k
, µ
k
) 是 n × n 正定矩阵，A(x
k
) 是 m × n 行满秩矩阵。
上述凸二次规划的全局极小点等价于方程中的 d
k
。下面，给出纯等式约束优化问题的 SQP
算法流程：
step1. 初始化。x
0
∈ R
n
, µ
0
∈ R
l
, ρ, r ∈ (0, 1), 0 ⩽ ε ≪ 1，置 k := 0。
step2. 计算 p(x
k
, µ
k
) 的值。若 p(x
k
, µ
k
) ⩽ ε 停止；否则转置 step3。
step3. 求解 (1.3) 式的凸二次规划，得到 ¯µ
k
, d
k
。并置
v
k
= ¯µ
k
− µ
k
−
1
2τ
A(x
k
)d
k
step4. 若 p(x
k
+ d
k
, µ
k
+ v
k
) ⩽ (1 − r)p(x
k
, µ
k
)，则置 α
k
:= 1，转到 step6，否则转到 step5。
step5. 令 m
k
是使下面的不等式成立的最小非负整数 m
p(x
k
+ ρ
m
d
k
, µ
k
+ ρ
m
v
k
) ⩽ (1 − rρ
m
)p(x
k
, µ
k
)
置 α
k
= ρ
m
k
。
step6. 令 x
k+1
= x
k
+ α
k
d
k
，µ
k+1
= µ
k
+ α
k
v
k
，置 k := k + 1，转到 step2。
不难发现，在上面的算法中，若 α
k
< 1，则必有
p(x
k
+ ρ
m
k−1
d
k
, µ
k
+ ρ
m
k−1
v
k
) > (1 − rρ
m
k−1
)p(x
k
, µ
k
)
并且该算法具有全局收敛性：
若 SQP 生成的序列 {(x
k
, µ
k
)} 使得 KKT 矩阵的逆矩阵 N (x
k
, µ
k
)
−1
一致有界，则 {(x
k
, µ
k
)}
的任何聚点 (x
∗
, µ
∗
) 都满足 p(x
∗
, µ
∗
) = 0。特别地，{x
k
} 的任一聚点是原问题的 KKT 点。下面
给出的是线性 SQP 的收敛速度。
设 SQP 产生的序列 {x
k
} 收敛到一个局部极小点 x
∗
，若 f, h 在 x
∗
附近二阶连续可微，Jacobi
矩阵 A(x
∗
) = ∇h(x)
T
行满秩，且二阶最优化充分条件成立，则有
(1) {µ
k
} → µ
∗
，其中，µ
∗
是等式约束问题的 Lagrange 乘子，且 {(x
k
, µ
k
)} 是二阶收敛的，
即
∥(x
k+1
− x
∗
, µ
k+1
− µ
∗
)∥ = O(∥(x
k
− x
∗
, µ
k
− µ
∗
)∥
2
)
(2) 序列 {x
k
} 超线性收敛到 x
∗
，且 t ∈ Z(正整数)
∥x
k+1
− x
∗
∥ = O(∥x
k
− x
∗
∥
t
Π
i=1
∥x
k−i
− x
∗
∥)
一般约束优化的 SQP 算法
将前面的等式约束问题的 SQP 思想推广到一般形式的约束优化问题。在给定 (x
k
, µ
k
, λ
k
) 之
后，将约束函数线性化，并且对 L(x, λ, µ) 进行二次多项式近似，得到下列形式的二次规划子问
题
min
1
2
d
T
G
k
d + ∇f(x
k
)d (1.4)
s.t.
h
i
(x
k
) + ∇h
i
(x
k
)
T
d = 0 i ∈ E
g
j
(x
k
) + ∇g
j
(x
k
)
T
d ⩾ 0 j ∈ I
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.5 最优化算法
其中：
G
k
= G(x
k
, µ
k
, λ
k
) = ∇
2
xx
L(x
k
, µ
k
, λ
k
)
L(x, µ, λ) = f(x) −
i
µ
i
h
i
(x) −
j
λ
j
g
j
(x)
于是 x
k
的校正步 d
k
以及新的乘子估计量 µ
k+1
, λ
k+1
可以分别定义为问题 (1.4) 的最优解 d
∗
和
相应的拉格朗日乘子 µ
∗
, λ
∗
。
上述的二次规划子问题 (1.4) 可能不存在可行点。为此，Powell 引进一辅助变量 ξ
min −ξ
s.t.
− ξh
i
(x
k
) + ∇h
i
(x
k
)
T
d = 0 i ∈ E
− ξg
k
(x
k
) + ∇g
j
(x
k
)
T
d ⩾ 0 j ∈ u
k
g
k
(x
k
) + ∇g
i
(x
k
)
T
d ⩾ 0 i ∈ v
k
− 1 ⩽ ξ ⩽ 0
其中：u
k
= {i|g
i
(x
k
) < 0, i ∈ I}，v
k
= {i|g
i
(x
k
) ⩾ 0, i ∈ I}。
注意到，在构建二次规划子问题 (
1.4) 时，需要计算 L(x, λ, µ) 在迭代点 x
k
处的 Hesse 矩阵
G
k
= G(x
k
, µ
k
, λ
k
)。但其计算量巨大，为了克服这一缺陷，1976 年，华裔数学家韩世平 (Han)
基于 N-L 方法提出了一种利用对称正定矩阵 B
k
替代 G
k
的 SQP。另外，Wilson 于 1963 年较
早考虑 N-L 方法。Powell 于 1977 年修正了 Ham 的方法，所以也称这种 SQP 为 WHP 方法。
在迭代点 (x
k
, µ
k
, λ
k
) 处，WHP 需要构造一个下列形式的二次规划子问题
min
1
2
d
T
B
k
d
+
∇
f
(
x
k
)
T
d (1.5)
s.t.
h
i
(x
k
) + ∇h
i
(x
k
)
T
d = 0 i ∈ E
g
j
(x
k
) + ∇g
j
(x
k
)
T
d ⩾ 0 j ∈ I
(1.6)
并且用此二次规划子问题的解 d
k
作为原问题变量 x
k
的搜索方向
。
下面，给出 WHP 的算法步骤：
step1. 初始化。x
0
∈ R
n
，初始对称矩阵 B
0
∈ R
n×n
，容许误差 0 < ξ ≪ 1，非负数列 {η
k
}，
∞
k=0
η
k
< +∞, σ > 0, δ > 0，置 k := 0。
step2. 求解二次规划子问题，解得 d
k
。
step3. 若 ∥d
k
∥ < ε 停止，输出 x
k
，否则转到 step4。
step4. 计算 ℓ
1
罚函数 p(x, σ)
p(x, σ) = f (x) +
1
σ
i
|h
i
(x)| +
j
|[g
j
(x)]
−
|
注：搜索方向 d
k
是许多罚函数的下降方向，比如 ℓ
1
罚函数。
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.5 最优化算法 第一章 约束非线性规划
其中：罚权重 σ > 0, [g
j
(x)]
−
= max{0, −g
j
(x)}。并按照某种线搜索规则确定步长 α
k
∈ (0, δ] 使
得
p(x
k
+ α
k
d
k
, σ) ⩽ min
α∈(0,δ ]
p(x
k
+ αd
k
, σ) + η
k
step5. 置 x
k+1
:= x
k
+ α
k
d
k
，更新 B
k
为 B
k+1
，置 k := k + 1，转到 step2。
下面，我们给出 WHP 的全局收敛性。设 f, h
i
, g
j
是连续可微的，且 ∃0 < m ⩽ M ，使对称
正定矩阵 B
k
满足
m∥d∥
2
⩽ d
T
B
k
d ⩽ M∥d∥
2
(∀d ∈ R
n
)
若罚权重 σ > 0 和二次规划子问题的拉格朗日乘子向量 µ
k+1
, λ
k+1
⩾ 0 满足
σ max{∥λ
k+1
∥
∞
, ∥µ
k+1
∥
∞
} ⩽ 1 (∀k)
则 WHP 产生的序列 {x
k
} 的任何聚点都是原问题的 KKT 点。在上述 WHP 中有两个存留的问
题：¬关于 B
k+1
的计算。关于 α
k
的求解。下面介绍求解 B
k+1
的 Powell 方法和增广拉格朗
日函数法，以及求解 α
k
的价值函数法。
Powell 的方法 B
k+1
的计算一般是用拟牛顿修正公式逐步迭代产生，我们希望 B
k+1
是 La-
grange 函数 Hesse 矩阵 G ≜ ∇
2
xx
L 的近似。我们令
s
k
= x
k+1
− x
k
y
k
= ∇
x
L(x
k+1
, µ
k+1
) − ∇
x
L(x
k
, µ
k+1
)
因为 BFGS 校正公式要求 s
k
, y
k
满足曲率条件。s
T
k
y
k
> 0，但上式确定的 s
k
, y
k
可能不满足这一
条件。为此，有必要对 y
k
进行修正。Powell 于 1978 年建议用下式对 y
k
进行修正
¯y
k
=
y
k
s
T
k
y
k
⩾ 0.2s
T
k
B
k
s
k
θ
k
y
k
+ (1 − θ
k
)B
k
s
k
其它
其中：θ
k
=
0.8s
T
k
B
k
s
k
s
T
k
B
k
s
k
−s
T
k
y
k
。
这种选取 ¯y
k
的基本思想是：利用 y
k
和 B
k
s
k
的凸组合构造一可以用来修正矩阵的向量。由
于 B
k
s
k
可理解为 y
k
的一种近似估计，且满足 (因为 B
k
正定)
s
T
k
(B
k
s
k
) > 0
故利用 y
k
和 B
k
s
k
的凸组合是一种很自然的选择。于是，矩阵 B
k
的约束 BFGS 校正公式为
B
k+1
= B
k
−
B
k
s
k
s
T
k
B
T
k
s
T
k
B
k
s
k
+
¯y
k
¯y
T
k
s
T
k
¯y
k
由 ¯y
k
的定义，不难验证 s
T
k
¯y
k
> 0。
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.5 最优化算法
基于增广拉格朗日函数的选择 另一种选择 B
k+1
的方法是基于增广拉格朗日函数。考虑下面的
增广拉格朗日函数
L
A
(x, λ, µ) = f(x) − λ
T
h(x) +
1
2µ
∥h∥
2
2
其中：罚权重 µ > 0。在局部极小点 (x
∗
, λ
∗
) 处，根据 h(x
∗
) = 0，可知增广拉格朗日函数的 Hesse
矩阵
∇
2
xx
L
A
(x
∗
, λ
∗
, µ) = ∇
2
xx
L(x
∗
, λ
∗
) + µ
−1
A(x
∗
)
T
A(x
∗
)
其中：等式右边第二项为曲率，对应着约束函数在 x
∗
点的法向量张成的空间 R(A
T
)。对于纯等
式约束的优化问题，若约束函数在 x
∗
处的 Jacobi 矩阵 A(x
∗
) 行满秩，并且二阶最优化充分条
件成立，则存在某个阈值
¯
µ
，使得
∀
µ
∈
(0
,
¯
µ
]
,
∇
2
xx
L
A
(x
k
, λ
k
, µ) 是正定的。于是，G(x
k
, λ
k
) 可以
取成正定矩阵 ∇
2
xx
L
A
(x
k
, λ
k
, µ) 或者对 ∇
2
xx
L
A
进行拟牛顿法近似的校正矩阵 B
k
。
步长 α
k
的选取 为了保证 SQP 的全局收敛性，通常借助某价值函数来确定搜索步长 α
k
。如：
目标函数 f、罚函数 p、增广拉格朗日函数 L
A
等都可以作为价值函数。
ℓ
1
价值函数：对于一般的约束优化问题，可以考虑相应的二次规划子问题 (1.5) 并且将罚函
数写为如下的 ℓ
1
罚函数 (ℓ
1
价值函数的形式)
p(x, σ) = f (x) +
1
σ
∥g(x)
−
∥
1
+ ∥h(x)∥
1
(1.7)
设 d
k
λ
k+1
⩾ 0, µ
k+1
为问题 (1.5) 的最优解和拉格朗日乘子向量，则 p(x, σ) 沿 d
k
的方向导数满
足
D(p(x
k
, σ), d
k
) ⩽ −(d
k
)
T
B
k
d
k
− (σ
−1
− ∥λ
k+1
∥
∞
)∥(g(x))
+
∥
1
− (σ
−1
− ∥µ
k+1
∥
∞
)∥h
k
∥
其中：(g
k
)
+
= max{0, −g(x
k
)}。
当然我们还有一些其它的价值函数可以选择，比如增广拉格朗日价值函数等
®
。
下面给出一般约束优化问题的 SQP 算法流程
step1. 初始化。给定初始点对 (x
0
, λ
0
, µ
0
)，对称正定矩阵 B
0
，计算 f
0
= f(x
0
)，∇f
0
= ∇f(x
0
)，
g
0
= g(x
0
)，h
0
= h(x
0
)，A
T
0
= (∇g(x
0
), ∇h(x
0
))，选择参数 η ∈ (0, 1/2)，τ ∈ (0, 1)，容许误差
ϵ
1
, ϵ
2
> 0，置 k := 0。
step2. 计算改进方向。求解二次规划自问题，得到变量 x
k
的改进方向 d
k
。
step3. 收敛性检验。若 ||d
k
||
1
⩽ ϵ
1
，并且 ||(g
k
)
−
||
1
+ ||h
k
||
1
⩽ ϵ
2
成立，则得到约束优化问题的
一个近似 KKT 点 (x
k
, λ
k
, µ
k
)，算法终止；否则，转到 step4。
step4. 确定罚因子。对于某种价值函数 ϕ(x, σ)，选择罚因子 σ
k
，使得 d
k
为该函数在 x
k
点的下
降方向。
step5. 步长选择。在序列 1, τ, τ
2
, . . . 中，选择最大的项作为 α
k
，使得
ϕ(x
k
+ α
k
d
k
, σ
k
) ⩽ ϕ(x
k
, σ
k
) + ηα
k
D(ϕ(x
k
, σ
k
), d
k
)
®
可以参考《数学规划》黄红选 P302 或者《最优化方法与 Matlab 程序设计》P231。
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.6 MATLAB 应用实例 第一章 约束非线性规划
step6. 改进迭代点。令 x
k+1
= x
k
+ α
k
d
k
，并且计算
f
k+1
= f(x
k+1
)
∇f
k+1
= ∇f(x
k+1
)
g
k+1
= g(x
k+1
)
h
k+1
= h(x
k+1
)
A
T
k+1
= (∇g(x
k+1
), ∇h(x
k+1
))
以及最小二乘乘子
λ
k+1
µ
k+1
= (A
k+1
A
T
k+1
)
−1
A
k+1
∇f
k+1
step7. 校正 Hesse 矩阵。令
s
k
= α
k
d
k
, y
k
= ∇
x
L(x
k+1
, λ
k+1
, µ
k+1
) − ∇
x
L(x
k
, λ
k+1
, µ
k+1
)
利用约束牛顿公式修正矩阵 B
k
，生成新的对称正定矩阵 B
k+1
，令 k := k + 1，返回 step2。
在上述算法中，我们隐设了矩阵 A
k
是行满秩的。如果这个条件不成立，那么在计算最小二
乘乘子时，就需要使用计算矩阵广义逆的技巧。最后，值得指出的是，对于无约束优化问题，如
果 x
∗
是驻点，并且目标函数 f(x) 在该点满足二阶最优性充分条件，即 Hesse 矩阵正定，那么只
要迭代点列 {x
k
} 收敛到 x
∗
，并且搜索方向列 {d
k
} 满足
lim
k→∞
||x
k
+ d
k
− x
∗
||
||x
k
− x
∗
||
= 0
对于充分大的 k，必然有
f(x
k
+ d
k
) < f(x
k
)
我们把这样的 d
k
称为超线性收敛步。
1.6 MATLAB 应用实例
MATLAB 通过 fmincon 函数来求解约束优化问题，其调用格式为
[x,fval,exitag,output,lambda,grad,hessian]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlan,options)
其中：fun 为目标函数；x0 为初始点；A,b 为线性不等式约束 Ax ⩽ b；Aeq,beq 为线性等式约束
Aeq x=beq；lb,ub 为自变量 x 的上下限，lb ⩽ x ⩽ ub；nonlcon 为非线性约束条件；options 为
结构体参数；lambda 为最优点 x 处的拉格朗日乘子；grad 为最优点 x 处的梯度；hessian 为最
优点 x 处的 Hesse 矩阵；
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 约束非线性规划 1.6 MATLAB 应用实例
我们用 fmincon 函数来求解如下非线性约束优化问题
min f(x) = x
4
1
− 4x
1
− 8x
2
+ 15
s.t.
9 − x
2
1
− x
2
2
⩽ 0
2x
1
+ 3x
2
⩽ 2
x
2
− x
1
⩽ 5
其求解程序为
1 fun = @(x) x(1 ) ^4−4x (1 )−8x (2) +15;
2 x0 = [ 1 , 2 ] ;
3 A = [ 2 ,3 ;1 , − 1] ;
4 b = [ 2 ; 5 ] ;
5 Aeq = [ ] ;
6 beq = [ ] ;
7 Lb = [ ] ;
8 Ub = [ ] ;
9 f unct i on [ c , ceq ] = ConFun(x)
10 c = 9−x( 1 )^2−x(2 ) ^2;
11 ceq = [ ] ;
12 end
13 nonlcon = @ConFun
14 o p t i ons = optimoptions ( ’ fmincon ’ , ’ Display ’ , ’ i t e r ’ , ’ Algorithm ’ , ’ sqp ’ ) ;
15 x = fmincon ( fun , x0 ,A, b , Aeq , beq , lb , ub , nonlcon , option )
16
http://www.ma-xy.com 19 http://www.ma-xy.com
