http://www.ma-xy.com
目录
第一章 几何规划 1
1.1 问题的引入与分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 模型规范化及其基础理论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2.1 几何规划解法介绍 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 正定式几何规划 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3.1 一般形式及对偶规划 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3.2
最优化方法
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 广义几何规划 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.1 最优化算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5 MATLAB 应用实例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1
http://www.ma-xy.com
第一章 几何规划
1.1 问题的引入与分析
考虑如下优化问题
max x/y
s.t.
2 ⩽ x ⩽ 3
x
2
+ 3y/z ⩽
√
y
x/y = z
2
xyz > 0
其中：x, y, z ∈ R 为决策变量。将上述问题转化为 GP(几何规划) 的标准形
min x
−1
y
s.t.
2x
−1
⩽ 1
1
3
x ⩽ 1
x
2
y
−1/2
+ 3y
1/2
z
−1
⩽ 1
xy
−1
z
−2
= 1
1.2 模型规范化及其基础理论
几何规划 GP 的一般形式为
min f
0
(x)
s.t. f
l
(x) ⩽ 1 l ∈ I
x > 0
其中：x ∈ R
n
为决策变量，I = {1, 2, . . . , L}，f
l
(x) 如下
f
l
(x) =
m
l
j=1
σ
lj
c
lj
n
i=1
x
a
lij
i
1
http://www.ma-xy.com
1.3 正定式几何规划 第一章 几何规划
l ∈ I, σ
lj
= ±1, c
lj
> 0, a
lij
为任意给定的实数，a
lij
∈ R。
几何规划 GP 最初由 Zener 于 1961 年提出，是非线性规划的一类分支。Zener 毕业于哈佛
大学，主攻物理。在后来的实践中，他发现许多工程设计问题均可归纳为一种特殊的形式 GP：
目标函数和约束函数 f
l
(x) 为自变量 x 乘幂之代数和的形式。Richard Dun 主要从事于非线性
对偶理论的研究，在发现 Zener 的工作后，与他一起从事这项研究。这种方法最初的发展使用
了算数几何平均不等式，因此，Dun 称这种问题为几何规划。Dun 和 Zener 主要是讨论正系
数的问题，也即正定式几何规划。1967 年，Passy 和 Wilde 提出了广义几何规划的伪对偶理论。
再后来，Peterson 提出一般无约束几何规划的对称对偶定理。70 年代压缩方法在求解广义几何
规划中很受欢迎。
1.2.1 几何规划解法介绍
正定式几何规划的局部极小点即为全局极小点，只需要对正定式几何规划进行变量替换：
x
i
= e
i
(i = 1, 2, . . . , n)，替换后的原正定式几何规划就等价转化为凸规划问题，原问题就等于凸
区域上极小化一个凸函数，这种性质是正定式几何规划的一个非常重要的特征。正定式几何规划
的主要解法有：¬割平面法；对偶方法；®压缩法；¯在半无限线性规划基础上的另一种新的对
偶方法。广义几何规划的主要解法有：对偶方法和两种压缩法等。
1.3 正定式几何规划
1.3.1 一般形式及对偶规划
下面，我们给出正定式几何规划 (PGP) 的一般形式
min
x
f
0
(x)
s.t. f(x) ⩽ 1 l = {1, 2, . . . , L}
x ⩾ 0
其中：x ∈ R
n
，f
l
(x) =
m
l
j=1
c
lj
n
i=1
x
a
lij
i
，l = 0, 1, 2, . . . , L。
我们对 PGP 进行变量替换，令 y
i
= logx
i
，从而将其变为凸规划：
f
l
(x) =
m
l
i=1
c
lj
n
i=1
(e
y
i
)
a
lij
=
m
l
j=1
c
lj
e
a
T
lj
·y
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 几何规划 1.3 正定式几何规划
再令 β
lj
= logc
lj
，有
f
l
(y) =
m
l
j=1
(e
β
lj
) · e
a
T
lj
·y
=
m
l
j=1
e
a
T
lj
·y+β
lj
于是，原问题变为如下问题
min
y
f
0
(y)
s.t. f
l
(y) ⩽ 1
其中：y 为优化变量，f
l
(y) =
m
l
j=1
e
a
T
lj
·y +β
lj
。对上式两边取对数，有
min
˜
f
0
(y) = log
m
l
j=1
e
a
T
0j
·y +β
j
s.t.
˜
f
l
(y) = log
m
l
j=1
e
a
T
lj
y +β
lj
⩽ 0
j = 1, 2, . . . , L
上述规划问题的目标函数和约束条件
˜
f
l
(y) 均为凸函数，所以原 PGP 转化为凸形式的几何
规划。下面讨论 PGP 的对偶问题。
设一个 m
L
行 n 列的矩阵为正定式几何规划的指数矩阵。
A =
a
011
a
012
. . . a
01n
.
.
.
.
.
.
.
.
.
.
.
.
a
0m
0
1
a
0m
0
2
. . . a
0m
0
n
.
.
.
.
.
.
.
.
.
.
.
.
a
L11
a
L12
. . . a
L1n
.
.
.
.
.
.
.
.
.
.
.
.
a
Lm
L
1
a
Lm
L
2
. . . a
Lm
L
n
m
L
×n
其中：a
lij
，l = 0, 1, . . . , L; i = 1, 2, . . . , n; j = 1, 2, . . . , m
l
，n 为自变量 x 的维度，m 是所有约束
函数 f
l
(x) 中的项数总和，即 m = m
0
+ m
1
+ . . . + m
L
，此时，记 d(p) =
L
l=0
m
l
j=1
(
p
l
0
c
lj
p
lj
)
p
lj
。
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.3 正定式几何规划 第一章 几何规划
称下面规划为 PGP 的对偶规划 (SGP)
max lnd(p) =
L
l=0
m
l
j=1
p
lj
lnc
lj
+
m
L
j=1
p
lj
− lnp
lj
s.t.
m
0
j=1
p
0j
= p
00
= 1
L
l=0
m
l
j=1
a
lij
p
lj
= 0 i = 1, 2, . . . , n
p
lj
> 0 l = 0, 1, . . . , L; j = 1, 2, . . . , m
l
其中：p
l0
=
m
l
j=1
p
lj
, l = 0, 1, . . . , L。
定理 如果一个几何规划 (不包含等式约束) 存在极小值，则有
¬对偶问题的极大值存在；
d(p
∗
) = f
0
(x
∗
)，其中 p
∗
, x
∗
为最优点；
®在最优点处，原变量与对偶变量的关系为：
c
0j
n
i=1
(x
∗
i
)
a
0ij
= p
∗
0j
f
0
(x
∗
)x
∗
i
j = 1, 2, . . . , m
0
c
lj
n
i=1
(x
∗
i
)
a
lij
=
p
∗
lj
p
∗
l0
p
∗
lj
(1 − f
l
(x
∗
)) = 0
其中：l = 1, 2, . . . , L。
此定理不仅告诉我们原问题的最小值等于对偶问题的最大值，而且只要求出对偶问题的最优
点和最优值通过上面的关系¬ - ®即可求出原问题的最小点。从形式上看，上面的关系®为未知
数的非线性方程组，取对数之后，就得到以 lnx
i
为未知数的线性方程组，即
n
i=1
a
0ij
lnx
∗
i
= ln
p
∗
0j
d(p
∗
)
c
0j
j = 1, 2, . . . , m
0
n
i=1
a
lij
lnx
∗
i
= ln
p
∗
ij
p
∗
l0
c
lj
l = 0, 1, . . . , L; j = 1, 2, . . . , m
l
此方程的个数为
m
L
l=0
= m，未知量的个数为 n。根据线性方程组的理论，可以解出 lnx
∗
i
(i =
1, 2, . . . , n)，从而得到 x
∗
= (x
∗
1
, x
∗
2
, . . . , x
∗
n
)。
由于对偶规划 (SGD) 是一凹规划，因此，许多算法都是基于它来求解原 PGP 的，但是目标
函数 d(p) 在 p
lj
= 0 上是不可微的。尽管这样，我们仍能发现对偶规划 (SGD) 有非常重要的特
征，其约束函数是线性的，而且目标函数 d(p) 是一凹函数。Beck 提出一种改进的单纯形算法来
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 几何规划 1.3 正定式几何规划
求解对偶规划 (SGD)。Alegandre 充分利用了对偶规划的结构及其线性约束的特征，并为了克服
不可微的缺陷，对所有变量增加严格正的下界约束，从而提出一种基于对偶的有效求解算法。
在求解 SGD 时，需要定义一个困难度
d = m − n − 1 =
l
l=0
m
l
− n − 1
困难度用来衡量几何规划求解难度的一个重要指标。一般来说，困难度越大，求解越困难。
1.3.2 最优化方法
信赖域内点算法
我们将 SGD 写为一般规划形式
min
x∈R
n
f(x)
s.t. x ∈ G
其中：
f(x) = C
T
x +
n
i=1
x
i
lnx
i
−
L
l=0
e
T
l
xlne
T
l
x
∆
= −lnd(P )
G = {x|Ax = b, x > 0}
x = (x
1
, x
2
, ··· , x
n
)
T
≜ P = (P
01
, P
0m
0
, ··· , P
L
1
, ··· , P
Lm
L
)
T
C = (c
1
, c
2
, ··· , c
n
)
T
≜ −(lnC
01
, lnC
0m
0
, ··· , lnC
L
1
, ··· , lnC
Lm
L
)
T
b = (1, 0, 0, ··· , 0)
T
∈ R
m
, n =
L
l=0
m
l
e
l
为从第 1 +
l−1
k=0
m
k
位分量到第
l
k=0
m
k
位分量 (共 m
l
个分量) 为 1，其余分量为 0 的 n
维列向量。令 A 为
A =
1 ... 1 ... 0 ... 0 ... 0 ... 0
a
011
... a
01m
0
... a
111
... a
11m
1
... a
L11
... a
L1m
L
a
021
... a
02m
0
... a
121
... a
12m
1
... a
L21
... a
L2m
L
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
a
0m−11
... a
0m−1m
0
... a
1m−11
... a
1m−1m
1
... a
Lm−11
... a
Lm−1m
L
则
g(x) = −∇f(x) = −[C + (ln(x
1
/e
T
0
x), ··· , ln(x
m
0
/e
T
0
x), ln(x
m
0
+1
/e
T
1
x), ··· , ln(x
m
0
+m
1
/e
T
1
x),
··· , ln(x
m
0
+m
1
/e
T
1
x), ··· , ln(x
n−m
L
+1
/e
T
L
x))]
T
∇
2
f(x) =
1
x
1
1
x
2
.
.
.
1
x
n
−
D
0
e
T
0
x
D
1
e
T
1
x
.
.
.
D
L
e
T
L
x
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.3 正定式几何规划 第一章 几何规划
其中：D
l
为全 1 的 m
l
× m
l
方阵。
我们已经知道，在最优解 x
∗
= p
∗
之后，只需解方程组即可求解原 PDP 的解 (t
∗
)，这里为
了防止混淆，我们都用 t
∗
来作为原问题的决策变量。
n
i=1
a
0ij
lnt
∗
i
= ln
p
∗
0j
d(p
∗
)
c
0j
j = 1, 2, . . . , m
0
(1.1)
n
i=1
a
lij
lnt
∗
i
= ln
p
∗
1j
p
∗
0l
/c
lj
l = 0, 1, . . . , L; j = 1, 2, . . . , m
l
(1.2)
其中：p
∗
l0
=
m
l
j=1
p
∗
lj
l = 0, 1, . . . , L。解得 lnt
∗
i
即为 SGP 的最优解 t
∗
。
引入线性约束内部可微凸规划的信赖域解法。假设 1：SGP 有最优解 x
∗
= p
∗
；假设 2：
rank(A) = m，即 A 的行向量组线性无关。
首先，我们定义如下矩阵
H = H(x) =
∆
0
(x)
∆
1
(x)
.
.
.
∆
l
(x)
−
D
0
σ +e
T
0
x
D
1
σ+e
T
1
x
.
.
.
D
L
σ +e
T
L
x
其中：
∆
0
(x)
∆
1
(x)
.
.
.
∆
l
(x)
=
1
x
1
1
x
2
.
.
.
1
x
n
∆
l
(x) 为 m
l
阶对角矩阵，1 > σ > 0 为一常数，记 H
k
= H(x
k
)。
引理 (1) H = H(x) 为正定对称矩阵。
引理 (2) 原问题的一阶稳定条件 (KKT 条件)：∃x
∗
∈ R
n
, λ
∗
∈ R
m+1
, µ
∗
∈ R
n
满足：
∇f(x
∗
) + A
T
λ
∗
= 0
Ax
∗
= b
u
∗
i
x
∗
i
= 0 i = 1, 2, ··· , n
x
∗
⩾ 0, u
∗
⩾ 0
定义 F = {x ∈ R
n
|Ax = b, x ⩾ 0}，
ˆ
F = {x ∈ R
n
|Ax
∗
= b, x > 0}。设 F 有界，
ˆ
F 非空，用
x
k
表示第 k 次求出的点。求试探步 d
k
，若接受，则 x
k+1
= x
k
+ d
k
。记 X
k
= diag(x
k
)，所有范
数为 L
2
范数，f
k
= f(x
k
), ∇f
k
= ∇f(x
k
)。
算法步骤：
step1. 初始化。
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 几何规划 1.3 正定式几何规划
x
0
∈
ˆ
F , δ ∈ (0, 1), η ∈ (0, 1), 0 < σ
1
< σ
2
< 1，置 k := 0.
step2. 求出 H
0
，∆k ∈ (δ, 1/δ)
0
。
step3. 求下面问题的全局 d
k
。
min ψ
k
(d) = ∇f
T
k
d +
1
2
d
T
H
k
d
s.t.
Ad = 0
∥X
−1
k
d∥ ⩽ ∆
k
使 x
k
+ d
k
> 0。
step4. 如果 ψ
k
(d
k
) = ψ
k
(0) = 0，停止计算。
step5. 做下降性测试，计算
r
k
=
f(x
k
) − f(x
k
+ d
k
)
−ψ
k
(d
k
)
step6. 若 r
k
< η，取 ∆
k
∈ (σ
1
∥X
−1
k
d
k
∥, σ
2
∆
k
)，返回 step3，重新求 d
k
；否则 ∆k = ∆k; x
k−1
=
x
k
+ d
k
。
step7. 计算 f
k+1
，∇f
k+1
，置 k := k + 1，返回 step2。
下面给出上述算法的全局收敛性。
假设 1：∃x
0
∈
ˆ
F ，且水平集 L = {x ∈ F |f (x) ⩽ f (x
∗
)} 是紧集；
假设 2：{H
k
} 一致有界；
假设 3：∀I ⊂ N, λ ∈ R
m+1
, u ∈ R
n
∇f(x) + A
T
λ − u = 0
Ax = b
x
I
= 0
u
i
= 0
i /∈ I
若有界，则只有孤立解。
引理 ψ
k
(d
k
) = ψ(0) = 0，当且仅当 x
k
是问题 P 的稳定点。
引理 算法产生的 x
k
如果不是问题 P 的稳定点，则存在可行下降方向 d = 0，也即存在有界
数 λ，有 d
T
∇f
k
< 0, x
k
+ λd ⩾ 0, d
T
∇f
k
< 0。
引理 算法产生的序列 {x
k
} 有上界。
定理 算法是有限的。
定理 x 为问题 P 的 KKT 点，则 x 为问题 P 的最优解。
定理 {x
k
} 由算法产生，若有限步终止于 x
L
，则 x
k
为问题 P 的最优解；{x
k
} 的任一聚点
¯x 为 P 的 KKT 点即最优解。
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.4 广义几何规划 第一章 几何规划
广义梯度投影内点算法
引理 ∀x ∈ G, X = diag(x), A = (a
0
, a
1
, ··· , a
m−1
)
T
，其中 a
j
(j = 0, 1, . . . , m − 1) 为 A 的
行向量，则矩阵 AXXA
T
是一对正定矩阵。
算法步骤：
step1. 初始化。
x
0
∈ G, σ
1
> 1, < σ
2
> 0，置 k := 0.
step2. 求 g(x
k
) = −∇f(x
k
)，B
k
, P
k
, µ(x
k
)。
B
k
= B(x
k
) = (AX
k
X
k
A
T
)
−1
AX
k
P (x
k
) = X
k
X
k
− X
k
AB(x
k
)
µ(x
k
) = (µ
0
(x
L
), µ
1
(x
L
), . . . , µ
m−1
(x
L
)) = B(x
k
)X
k
g(x
k
)
step3. 求 s
k
= P (X
k
) = X
k
g(x
k
)。
ρ
k
=
|g(x
k
)
T
D
k
s
k
|
σ
1
|µ(x
l
)
T
w|+ σ
2
其中：w = (−1, −1, . . . , −1)
T
∈ R
m
。
step4. 求 d
k
= s
k
+ ρ
k
p
k
s
k
。
step5. 如果 g(x
k
)
T
X
k
d
k
= 0，则 x
k
为问题 P 的 KKT 点，否则转到 step6。
step6. 取 x
k+1
= X
k
(e + λ
k
d
k
) ∈ G，其中 e = (1, 1, ···)
T
∈ R
n
，λ
k
为搜索步长，置 k := k + 1
返回 step2。
下面给出上述算法的全局收敛性：
(1) 算法有限步终止于问题 P 的某一 KKT 点，x
k
也即最优点；
(2) 算法产生一无穷序列 {x
k
}，记 x
∗
为其聚点，若 x
∗
> 0，且 ∇f(x) 在 G 上关于 {x
k
} 一
致连续，则 x
∗
为问题 P 的 KKT 点，也即最优点。
1.4 广义几何规划
广义几何规划的一般形式为
min f
0
(x)
s.t.
f(x) ⩽ a
l
l = 1, 2, . . . , L
x > 0
其中：x ∈ R
n
，f
l
(x) =
m
L
j=1
σ
lj
c
lj
n
i=1
x
a
lij
i
, l = 0, 1, . . . , L。σ
lj
= ±1 至少有一个是负的，a
l
= ±1，
c
lj
> 0, l = 0, 1, . . . , m
L
，a
lij
为任一实数。
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 几何规划 1.4 广义几何规划
广义几何规划的对偶形式为
max d(P ) = a
0
L
l=0
m
l
j=1
p
l0
c
lj
p
lj
σ
lj
p
lj
a
0
(1.3)
s.t.
m
0
l=1
σ
0j
p
0j
= a
0
L
l=0
m
l
j=1
σ
lj
a
lij
p
lj
= 0
p
lj
⩾ 0, i = 1, 2, ··· , n
(1.4)
其中，规定
lim
p
lj
→0
p
l0
c
lj
p
lj
σ
lj
p
lj
= 1
L 个线性不等式约束
p
l0
= a
l
m
l
j=1
σ
lj
p
lj
⩾ 0 l = 1, 2, . . . , L (1.5)
其中：p
00
= 1，
a
0
=
1 f
0
(x
∗
) > 0
− 1 f
0
(x
∗
) < 0
定理 设有 m =
L
l=0
m
l
个对偶变量 p
lj
(l = 0, 1, . . . , L; j = 1, 2, . . . , m
L
) 满足上述约束条件，
x
∗
是原问题的最小点，则对 f
0
(x) 的每个局部极小点 x
0
，存在一个对偶变量 p
0
满足式 (1.3)(1.5)
使得 d(p
0
) = f
0
(x
∗
)，并且原有变量 x
0
与对偶变量 p
0
之间有如下关系：
c
0j
n
i=1
x
i
a
0ij
= p
0j
f
0
(x
0
) j = 1, 2, . . . , m
0
c
lj
n
i=1
x
i
a
0ij
=
p
lj
p
l0
l = 1, 2, . . . , L; j = 1, 2, . . . , m
L
以及
p
0
lj
(1 − f
l
(x
0
)) = 0 l = 1, 2, . . . , L; j = 1, 2, . . . , m
l
下面，给出此问题的几点说明：
¬f(x) ⩾ d(P ) 不一定成立。
对于每一个 f (x) 的局部极小点 x
0
，必存在 p
0
使得 f (x
0
) = d(p
0
)，x
0
并非总体极小。
®当困难度为 0 时，由对偶约束确定唯一 p
0
与 d(p
0
)，又可确定 x
0
，可以证明原问题有极小
点时，x
0
即为极小点，f(x
0
) 即为极小值。
¯由于事先不知道原问题 (P) 的极小点，a
0
的取值带的符号难定，不过定错了也没关系，若
a
0
取错，所有对偶变量全为负值，与要求不符，就可以纠正 a
0
的值，再重新计算。
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.4 广义几何规划 第一章 几何规划
1.4.1 最优化算法
考虑不等式约束下的广义几何规划 (GPE)
min f
0
(t)
s.t. f
l
(t) ⩽ 0 l = 1, 2, . . . , L
其中：t ∈ R
n
，
f
l
(t) =
m
L
j=1
c
lj
n
i=1
t
i
a
lij
l = 1, 2, . . . , L; j = 1, 2, . . . , m
l
, i = 1, 2, . . . , n
c
lj
和 a
lij
为任意实数。
令 x
i
= lnt
i
，则 GPE 转化为如下等价形式 GP E
1
：
min C
0
(x) =
m
0
j=1
c
0j
e
n
∑
i=1
a
0ij
x
i
s.t. C
i
(x) =
m
l
j=1
c
lj
e
n
∑
i=1
a
lij
x
i
⩽ 0 l = 1, 2, . . . , L
将上式写为向量形式 GP E
2
，记 e
x
= (e
x
1
, e
x
2
, ··· , e
x
n
)
T
，A
T
l
= (a
lij
)
n×m
l
，b
l
= (b
l
1
, b
l
2
, ··· , b
l
m
l
)
T
，
l = 1, 2, . . . , L; j = 1, 2, . . . , m
L
, i = 1, 2, . . . , n，有
min C
0
(x) = c
T
0
e
A
0
x
s.t. C
l
(x) = c
T
l
e
A
l
x
⩽ 0 l = 1, 2, . . . , L
构造无约束优化问题。设上面 GPE 的拉格朗日函数为
L(x, λ) = C
0
(x) +
L
l=1
λ
l
C
l
(x)
记 λ = (λ
l
, λ
2
, ··· , λ
L
)
T
，C(x) = (C
1
(x), C
2
(x), ··· , C
L
(x))
T
。
引入拉格朗日变量 Z
l
，将 GP E
2
中的不等式约束变为等式约束。
C
L
(x) + Z
2
l
= 0 l = 1, 2, ··· , L
考虑 C
0
(x) 在这组等式约束下的增广拉格朗日函数
¯
L(x, λ, z, σ) = C
0
(x) +
L
l=1
λ
l
(C
l
(x) + Z
2
l
) +
σ
l
+
L
l=1
(C
L
(x) + Z
2
l
)
2
其中：σ > 0 为罚权重。为消去 Z
l
将
¯
L 关于 z 求极小，即令 ∇
z
¯
L(x, λ, z, σ) = 0，得 Z
l
(λ
l
+
σ
l
(C
l
(x) + Z
2
l
)) = 0, l = 1, 2, ··· , L。若 σ
l
C
l
(x) + λ
l
⩽ 0，则 Z
2
l
= −
λ
k
σ
l
− C
l
(x)，否则 Z
l
= 0。
因此，Z
2
l
= −
1
σ
l
max(0, −λ
l
−σ
l
C
l
(x)), l = 1, 2, ··· , L。从而可得问题 GP E
2
的拉格朗日函数为
L(x, λ, σ) = C
0
(x) +
1
2σ
L
l=1
{[max(0, σC
l
(x) + λ
l
)]
2
− λ
l
}
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 几何规划 1.5 MATLAB 应用实例
之后，确定搜索方向 d
k
和步长 α
k
即可确定算法
¬
。
1.5 MATLAB 应用实例
悬梁臂的设计示意图如图 (1.1) 所示
图 1.1: 悬梁臂示意图
设悬梁臂共有 N 块组成，第 i 块的长为 l，宽为 w
i
，高为 h
i
(各块长均为 L)。在悬梁臂末
端悬挂重物 F ，要求设计各臂块 i 的 h
i
, w
i
，来使悬梁臂总体体积 V 最小。当然，悬梁臂要满足
一定的约束条件：
¬边界约束
w
min
⩽ w
i
⩽ w
max
, h
min
⩽ h
i
⩽ h
max
, i = 1, 2, ··· , N
长宽约束
s
min
⩽ h
i
/w
i
⩽ s
max
®最大弯曲应力约束。要求对每一臂块 i-th，在其弯曲应力应有最大上限。引入：弯曲应力
σ
b
= M(x) · y/I
其中：M (x) 为在 x 处的力矩，x 是终端到负载的距离，I 为悬梁转动的惯量面积 (惯性力矩)，
所以对具体的第 i-th 悬梁，F · D
i
为最大力矩，其中 D
i
为最大距离，所以弯曲应力为
σ
i
=
F · D
i
· (
h
i
2
)
I
i
惯性力矩为
I
i
= w
i
· h
3
i
/12
将 I
i
代入 σ
i
中，有
σ
i
=
F · D
i
w
i
h
i
¬
几何规划问题的算法研究-晋香莲 P16
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.5 MATLAB 应用实例 第一章 几何规划
要求
σ
i
⩽ σ
max
i = 1, 2, ··· , N
¯最后的约束是在悬梁末端 i-th= 1 时垂直偏转的限制。我们要求末端偏转有上限，即 y
i
⩽
y
max
，偏转 y
1
可以由递归关系式得到：
v
i
= 12
D
i
−
1
2
F
Ew
i
h
3
i
+ v
i+1
y
i
= 6
D
i
−
1
3
F
Ew
i
h
3
i
+ v
i+1
+ y
i+1
对应 i = N, N − 1, ··· , 1，我们从 y
N+1
= v
N+1
= 0 开始设计递归。
综上，悬梁臂的设计为
min V = l ·
N
i=1
w
i
h
i
s.t.
w
min
⩽ w
i
⩽ w
max
h
min
⩽ h
i
⩽ h
max
s
min
⩽ h
i
/w
i
⩽ s
max
6F D
i
/(w
i
h
2
i
) ⩽ σ
max
y
1
⩽ y
max
上面的设计是一个几何规划问题。当然，对于约束¯，在处理悬梁臂末端偏转 y
1
时，我们可
以用 Castigliano’s 第二定理进行设计。
y
1
=
∂u
∂F
其中：y 为梁臂的偏转，u 是由外力 F 产生的能量，u =
L
0
M
2
/2EIdx，M 为力 F 在 x 处的力
矩。给定 M = F x，我们可以写出 u
u = F
2
/2E
L
0
N
i=1
(x + (N − 1)L)
2
/I
i
dx
将
u
关于
F
求偏导
y
1
=
∂u
∂F
= P · E ·
L
0
N
i=1
((x + (N − 1)L)
2
/I
i
dx
要求 y
1
⩽ y
max
。
MATLAB 求解时的设置为：设置参数：F = 50000N，L = 500cm，l = 100cm，y
max
= 2.7cm，
σ
max
= 14000N /cm
2
，E = 2 × 10
T
N/cm
2
，s
max
= 20cm。并且我们额外要求：w
1
h
1
为整数，
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 几何规划 1.5 MATLAB 应用实例
w
i
h
i
的边界要求为
1 ⩽ w
1
⩽ 5
30 ⩽ h
1
⩽ 65
2.4 ⩽ w
2
, w
3
⩽ 3.1
45 ⩽ h
2
, h
3
⩽ 60
1 ⩽ w
4
, w
5
⩽ 5
30 ⩽ h
4
, h
5
⩽ 65
用 MATLAB 进行求解，求解程序如下
1 fun = @cantileverVolume ;
2 l b = [ 1 30 2.4 45 2. 4 45 1 30 1 3 0 ] ;
3 ub = [ 5 65 3 .1 60 3.1 60 5 65 5 6 5 ] ;
4 A= [ ] ; b = [ ] ;
5 Aeq= [ ] ; beq =[ ] ;
6 nonlcon = @ca nti lev er Con str ain ts ;
7 opts = gaoptimset ( . . .
8 ’ Popu l ationS ize ’ , 150 , . . .
9 ’ Generations ’ , 200 , . . .
10 ’ EliteCount ’ , 10 , . . .
11 ’TolFun ’ , 1e −8, . . .
12 ’ PlotFcns ’ , @ga plotbestf ) ;
13 rng (0 , ’ t w i st e r ’ ) ;
14 [ xbest , fbe st , e x i t f l a g ] = ga ( fun , 10 ,A, b , Aeq, beq , . . .
15 lb , ub , nonlcon , [1 2] , opts ) ;
16
http://www.ma-xy.com 13 http://www.ma-xy.com
