约束条件下多变量非线性函数的区间算法

一种求解约束条件下多变量非线性函数所有全局最优解的区间算法*
李爽，许才军， 王新洲
（武汉大学 测绘学院，湖北，武汉，430079）
摘 要：研究和实践中经常会遇到附有约束条件的非线性优化问题，对这类问题，通常采用随机搜索的方法来解决，但是，随机搜索法不能证明所得到的解就是全局最优解。

本文给出了一种求解约束条件下非线性优化问题所有全局最优点和最优值的区间算法，认为该算法非常宜于解决优化问题，它能求出问题的所有全局最优解，给出解的包含区间，并很容易获得解的逼近误差，这是随即搜索等其它方法做不到的。

理论分析和数值结果均表明，区间算法是稳定而可靠的。

关键字：区间算法；全局约束最优；非线性函数
0 概述
研究和实践中经常遇到某些问题的模型参数，由于要考虑到很多具体设计或参数的物理意义等方面的要求，会受到某种约束或限制的情况。

比如在地球物理反演当中，断层的长度和宽度不能为负值，甚至约束可能是非线性的，这就使问题成为附有约束条件的非线性优化问题。

对这种问题的求解，目前用的较多的是遗传算法和模拟退火算法，但这两种方法由于采用随机搜索，不能保证每次运算都能得到全局最优解，对这种随机搜索的改进也只是提高了运算获得全局最优解的概率，不能达到准确无误的效果。

本文提出一种求解该问题的扩展的区间算法，采用穷尽搜索原则，能够同时求出目标函数的所有全局最优解，解决了其他方法难以做到的问题。

在区间运算中约定，I 是实数R 上所有闭区间R b a b a ∈,],[，的集合，
m I 是I 上的m 维区间矢量，m I X ∈是优化问题的定义域，那么优化问题的一般表示为
⎪⎪⎩⎪
⎪⎨⎧+===≤∈r
k j x h k i x g x f j
i X x ,,10)(,,10)()(min
定义1 函数→m R f :R ，设m I X ∈在f 的定义域内，若对任意X x ∈，都有
)(}:)({X F X x x f ⊆∈
则称I I F m →:为函数f 的包含函数[1]。

定义 2 如果I A ∈，],[b a A =，则A 的中点为2/)()(b a A m i d +=，宽度为a b A w -=)(；如果m I A ∈，T 1),,(m A A A =，I A i ∈，则A 的中点为
T
1))(,),(()(m A mid A mid A mid =，宽度为))((max )(1
i m
i A w A w ==。

定义3 如果I A ∈，],[b a A =，则定义)sup(A 表示区间A 的上界，定义inf (A)表示区间A 的下界，即sup (A)=b ，inf (A)=a 。

其它相关的定义可参见文献[2]。

1 算法
设初始域为X ，H G F ,,分别是h g f ,,的包含函数，L 是以),,(R y Y 为元素的区间列，其中Y 是分割得到的子域、y 是子域上目标函数的极小值，),,,(21r R R R R =是子域的标志矢量，R 的分量取值0或1，若0=i R ，表明到目前为止，对任意Y x ∈，还不能判断是否满足0)(≤x g i 或0)(=x h i ，称Y 是未定的；若1=i R ，则表示，对任意Y x ∈，都满足0)(≤x g i 和0)(=x h i 。

f ~
表示到目前为止f 在可行点的最小值。

把多维区间向量称为盒子（box ）。

* 本文得到国家自然科学基金项目（编号：49904001）以及高等学校骨干教师资助计划资助(2123) 李爽（1971- ），女，博士生。

现主要从事大地测量地球物理反演和优化理论的研究，email ：
geo_lish@
1.1 穷尽准则
区间算法是使穷尽准则得以实现的主要工具，其它算法迫于无法解决该准则的计算量才转而求助随机搜索。

穷尽准则在区间算法中有两次体现。

体现1 初始化域X 被分割，测试每个子域是否满足等式和不等式约束。

假设Y 是分割产生的一个子域，计算k i Y G i ,,2,1),( =和r k k i Y H i ,,2,1),( ++=，判断Y 是否可行。

若Y 可行，则保留至L 中；不可行，则舍弃。

若不能做出Y 是否可行的判断，则Y 被分割，以至进一步判断Y 的子域是否可行。

把L 中的子域按照y 不递减的次序进行排列，以使得算法总分割第一个未定（无法判断是否可行）的子域，而且，在这个子域上，目标函数具有最小的下界，这利于优先求取全局最小点。

体现2 用),,(R z Z 表示L 中的所有元素，对它们进行中点测试，舍弃L 中满足z f <~
的不可能包含最优解的子域。

1.2 计算步骤
假设目标函数的参数个数为m ，用x ~表示一个已知的可行点，给定松弛参数01≥ε，如果01=ε，则原约束优化问题被解决，如果01>ε，则表示通过把原问题的等式约束
0)(=x h j )(X x ∈用不等式约束],[)(11εε-∈x h j ),,1(X x r k j ∈+= 所代替，解决原
问题的松弛问题。

计算步骤为
步骤1 设X Y =；
步骤2 如果x ~存在，取)~(~x f f =，否则∞=f ~
；
步骤3 设)0,,0(),,(1 ==r R R R ，)(inf Y F y =；
步骤4 初始化区间列)),,((R y Y L =；
步骤 5 在方向},,1),()(:{m i Y w Y w i v i ==∈分割Y ，得到21V V 、，并使
21V V Y =，然后从L 中删除),,(R y Y 的纪录，但保留R 作为21V V 、的标志；
步骤6 对2,1=j
1）设R R R R j
r j
j
==),,(1 ； 2）计算)(inf j j V F v =； 3）如果j v f <~
，转至⑺；
4）对于k i ,,1 =，如果0=i R ，计算)(j i V G 。

若0)(>j i V G ，转至⑺；若0)(1≤V G i ，则1=j i R ；
5）对于r k k i ,,2,1 ++=，如果0=i R ，计算)(j i V H 。

若)(0j i V H ∉，转至⑺；若],[)(11εε-⊆j i V H ，则1=j i R ；
6）若)1,,1( =j
R ，则)(j j V mid c =，))(,~
min(~j c f f f =，并对j V 进行单调性测试，若测试后未被删除，则把),,(j
j j R v V 按第二个元素不递减的顺序存入L ；
7）结束j 循环。

步骤7 若L 为空，则步骤13，终止计算，表明初始域内没有可行点； 步骤8 指定L 的第一个元素为)~
,~,~(R y Y ；
步骤9 中点测试：删除L 中满足z f <~
的所有元素； 步骤10 把L 的第一个元素指定为),,(R y Y ；
步骤11 如果终止条件成立，转至步骤12，否则，转至步骤5； 步骤12 终止时对L 中的l Z Z ,,1 ，有两种可能：
1）如果没有标志矢量)1,,1( =R ，则表明，或者不能判定可行域是否为空，或者如
果已知可行域D 不为空，能够得出i l i Z D X 1*=⊆⊆ ，*~f y ≤；
2）如果存在标志矢量)1,,1( =R ，则能够得出0/≠D ，i l i Z X 1*=⊆ ，f f y ~
*~≤≤；
步骤13 计算结束。

1.3 收敛条件
设n L 是上述算法第n 次迭代时的区间列，n l 是其长度，即元素个数，n nl n Z Z ,,1 是n
L 中的所有区间，ni l
i n Z U n 1== ，在n L 中，y ~的值被表示为n y ~，f ~
的值被表示为n f 。

对任何)(X I Y ∈，对于h g f ,,的包含H G F ,,，假设当0)(→Y w 时，
0))((→Y F w ，0))((→Y G w ，0))((→Y H w （1） 成立，即h g f ,,都是连续函数。

还假设，在可行域D 的内部存在一个点列n x ，满足当∞
→n 时有
**X x x n ∈→ （2）
引理1 设∞=1)(n n Y 是算法产生的区间序列，)(n n Y w w =，那么当∞→n 时，0→n w 。

证明：假设初始区间是X ，其边长数是m ，w 是所有边长中的最大值，则满足w Y w w ≤<)(2/的所有区间Y 的数量可由数列,11=N 121+=-m m mN N 算出，即至多有m N 个。

满足2/)(4/w Y w w ≤<的区间Y 的数量至多m N m !)!2(个，也就是有限个，由于∞=1)2/(k k w 是零序列，所以，对任意的0>ε，只有有限个盒子的宽度大于等于ε。

因此，可以得出，任给一0>ε，存在0n ，使得对任意0n n ≥，总有ε<)(n Y w 。

下面用反证法，假设对0>ε，存在0n ，使得对任意0n n ≥，ε≥)(n Y w ，这说明存在∞=1)(n n Y 的一个子列∞=0)(v k v Y ，并且对 ,1,0=v ，有ε≥)(v k Y w 。

这与前面只有有限个盒子的宽度大于等于ε相矛盾。

所以引理得证，当∞→n 时，0→n w 。

引理2：设X 为上述算法的初始区间，松弛参数01=ε，（1）、（2）两式成立，则)
(n U 是嵌套序列，并且n n U D ∞==1 ，即D U n →；当∞→n ，*~f y n ≤，*~f y n →，n f ↘
*f 。

证明：从算法本身可直接得出)(n U 是嵌套序列，D U n ⊇，*~f y n ≤，n f f ≤*和
n n f f ≤+1。

顾及对于n l i ,,1 =，)()(ni n Z w Y w ≥和引理1，能够得出
当∞→n ，0)(→ni Z w （3）
首先证明*~f y n →。

对所有n l i ,,1 =，有ni n n z Y F y ≤=))~
(inf(~成立。

根据（3）
，当∞→n 时，有0)~
(→n Y w ，再根据（1），有0))~((→n Y F w 。

由于)~(*n Y F f ∈，所以
*~f y n →。

用假设法证明)~(*n Y F f ∈，假设)~(*n Y F f ∉，根据*~
f y n ≤，则有*~0f f y n <→，由此得出0)~(f Y F n →和*))~(sup(f Y F n <，即对足够大的n ，n Y ~将不
是可行域，应该已经被算法在中间的处理中舍弃掉。

所以前面假设不成立，*~f y n
→得证。

下面证明n f ↘*f 。

由n f f ≤*和n n f f ≤+1足够得出*f f n →。

设)(n x 是满足（2）
的点列，∞=1)(n k n f 是∞=1)(n n f 的子列，
给定n x ，根据上述算法有)(n k x f f n ≤，由于*x x n →和（1），所以*)(f x f n →，)(*n k x f f f n ≤≤，顾及n n f f ≤+1，得到n f ↘*f 。

用反证法证明n n U X ∞==1* 。

假设存在*X x ∉，但对所有的n 有n U x ∈，意即n Z x ∈，n n L Z ∈。

根据（1）
，0)(→n Z w ，因此有x Z n →和)()(x f Z f n →，由于n f ↘*f 和*)(f x f >，所以n Z 将在中点测试中被舍弃，这与假设矛盾，得证。

2 算例
约定ε作为终止参数，即当区间列中所有区间的宽度都小于ε时，停止计算；t 表示迭代次数；n 表示终止时区间列中盒子的个数；L 表示终止时区间列中的盒子；tf 表示对F 的
估计次数；tg 表示对G 的估计次数。

算例1 最小化DeJong 函数F2[3]
⎪⎩⎪⎨
⎧=≤≤--+-=)2,1(048
.2048.2)1()(100),(2
12221212i x x x x x x f i 这是一个二维函数，它具有一个全局极小点,1(2f 0)1=，虽然是单峰函数，但它却是病态的。

计算中，取710-=ε，得到233=t ，657=tf ，6=n ，它们是
⎝
⎛=,7929680.99999987[896484,0.99999993
[1X ⎪⎪⎭⎫]8964850.99999993]0000011.00000000 ⎝⎛=,8964840.99999993[896484,0.99999993[2
X ⎪⎪⎭
⎫]0000011.00000000]0000011.00000000 ⎝
⎛=,0000001.00000000[896484,0.99999993
[3X ⎪⎪⎭⎫]1035161.00000006]0000011.00000000 ⎝⎛=,1035151.00000006[000000,1.00000000[4
X ⎪⎪⎭
⎫]2070321.00000012]1035161.00000006 ⎝
⎛=,0000001.00000000[000000,1.00000000
[5X ⎪⎪⎭⎫]1035161.00000006]1035161.00000006 ⎝⎛=,2070311.00000012[000000,1.00000000[6
X ⎪⎪⎭
⎫]3105471.00000018]1035161.00000006 030-4207637E 4.97968446~=f ，f 落在区间0,[]~f 内。

算例2 目标函数及其约束条件为
⎩⎨
⎧≤-+=0
)π4sin()π2sin(2)
(1.0),(122
22121x x x x x x f 在其定义域T ])1,1[],1,1([--内求目标函数的极小值。

该函数至少有24个局部最小值，但只有一个全局最小点T x )0,0(*=和全局最小值0)(*=x f 。

取410-=ε，得出115=t ，294=tf ，179=tg ，4=n ，这4个区间的并集为,46104.0[--E ]46104.0-E ，10-9230957E 1.86264514~=f ，函数值落在区间,0[]~f 内。

算例3 最小化Shubert 函数：
⎪⎩
⎪⎨⎧
=≤≤-++⋅++=∑∑==)2,1(1010]
)1cos[(])1cos[(),(5
1
251121i x i x i i i x i i x x f i i i 这是一个多峰值函数，在其定义域内总共有760个局部最小点，其中的18个点是全局最小点，全局最小值为731.186-=f 。

取610-=ε，计算得到，1450=t ，4559=tf ，18=n ，即恰好L 中的每个元素都是一个全局最小点的包含区间（这里没必要把18个包含区间全部
列出来），函数的全局最小值落在区间]8997309088309.186,6537309088332.186[--内。

算例4 最小化函数]10,10[],
)1sin[()(5
1
-∈++=
∑=x i x i i x f i
该函数在文献[4-5]中有给出，但据本文的计算，其给出的结果是不正确的，本文中取
610-=ε，得到函数的最大值44216712,[12.031249
max ∈f 4227943]12.0312494，101=t ，141=tf ，3=n ，函数的最大值点分别是
770722]5.79179465 6166076,[5.7917940029663]
6.77457571- 30634308,[-6.774576431793]
0.49139082- 2036438,-0.4913914[ 函数的最小值002584705,[-14.83795min ∈f 2570984]14.8379500-，79=t ，102=tf ，3=n ，函数的最小值点分别是
025239]5.16908586 6420593,[5.1690852775146]
7.39728450- 10379792,[-7.397285177276]
1.11409962- 1781922,-1.1141002[ 文献[5]给出该函数的442-1
2.031249*=f ，这个结果不对，我们得出
4212.0312494max =f ，026-14.837950
min =f ，看来应该是作者疏忽没把符号去掉，从下面该函数的图形可以知道本文给出的结果是对的，从表1还能看出本文的算法在所需的迭
算例5 最小化函数),18sin()4.13()(x x x f -=,0[∈x ]2.1 该函数在文献[4-6]都有运算。

本文取
610-=ε，得到函数的最小值
538732[-1.48
9072m i n ∈f 868161.48907253-，23=t ，29=tf ，1=n ，即函数的最小值点是689942]0.96608619 469482,0.96608562[。

通过表2的比较，可以看出，
本文的算法优于其它算法。

算例6 带有复杂约束条件的函数
⎪⎪⎪
⎪⎪⎪⎩
⎪
⎪⎪
⎪⎪⎪⎨
⎧
=≤-≤++-≤++≤++-≤-+-≤-+-+-+-++--+--=)3,2,1(01.461.29712128.341251211372428.023),,(min 321321321321321321321321321321i x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x f i 该函数的全局最小点为471428.2)0,0,1(-=f 。

计算中，初始区间取],100,0([0=X ],100,0[]),100,0[610-=ε,得到316=t ，184=n ，3435174
-2.4714256~
=f ，6521=tf ，8432=tg ，这里没必要把L 中这184个盒子都打印出来，但计算结果说明，这184个盒子，与全局最小点)0,0,1(的差距都小于ε，并且限
定了f 的取值落在区间]~
3080221,-2.4714321[f 内。

4 结论
通过运算发现，对于像算例1、2的目标函数，宜采用函数自然区间扩展，不宜采用中点或泰勒型包含，因为原目标函数是非负数，采用其它的包含会导致函数值出现负数的情况，反而破坏了包含的精确性。

由于区间算法是对全域的搜索，所以运算会比点估计的算法慢，为了加快算法的收敛速度，文献[7]给出了一种快速区间算法。

在具体应用方面，M. Tarvaine n 等采用区间算法进
行地震定位[8]
，取得了很好的反演结果。

但对于目标函数是复杂的分式的情况，目前国内外的成果很少，本文采用扩展的区间算法来进行运算，取得了很好的结果。

参考文献：
[1] FERNANDEZ, JOSE and BLAS PELEGRIN, Using Interval Analysis for Solving Planar Single-Facility Location Problems: New Discarding Tests. Journal of Global Optimization 19:61-81, 2001.
[2] Ratschek, H. and Rokne, J. (1988), New Computer Methods for Global Optimization . Ellis Horwood Ltd., Chichester.
[3] Zhou Ming, Sun Shudong, Theory and application of genetic algorithm [M]. Beijing: Publishing house of national defence and industry. 2001. (in Chinese)
[4] Shen Z, Zhu Y . An interval version of Shubert ’s iterative method for the localization of the global maximum[J]. Computing, 1987, 38: 275-280
[5] Shen, Peiping, Shouzhi, Yang. An Interval Algorithm for Finding all Global Minimizers of a Nonsmooth Function of Several Variable. MA THEMATICA APPLICATA. 2001,14(1):15-21 [6] Cao D X, Shen Z H. Interval algorithms for a class of nonsmooth optimization problems[J]. Numerical Mathematics. A Journal of Chinese Universities, 1998, 1: 23-33
[7] CASADO, L. G., J. A. MARTiNEZ and I. GARCiA, Experiments with a new selection criterion in a fast interval optimization algorithm. Journal of Global Optimization 19: 247-264, 2001.
[8] Tarvainen, M., Timo Tiira and Eystein S. Husebye, Locating regional seismic events with global optimization based on interval arithmetic, Geophys. J.Int , 1999, 138:879-885
An Interval Algorithm for Finding all Global Minimizers of a Constrained
Nonlinear Function with Several Variables
LI Shuang1Xu Caijun1Wang XinZhou1
(1 School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, Hubei, China)
(2 School of Remote-Sensing Information, Wuhan University, Wuhan 430079, Hubei, China) Abstract: It is frequently the case that a mathematical model of a system contains limitations on the acceptable values of the parameters. This leads to the most general global optimization topic, namely the global nonlinear constrained optimization problem. General, the random search is used to solve this kind of problem. But the random algorithms can’t test if it’s result is the optimum of problem. An interval algorithm is established in this paper. This algorithm can solve the global nonlinear constrained optimization problem and is very suitable to solve optimization problem. It can find all global optimums and all global optimizers of the problem, and the inclusion interval of optimizers is output. By these intervals, the approaching errors can be gotten easily. This can’t be done by the others algorithms, such as the random search method. Both theory and numerical analysis show that this algorithm is stable and reliable.
Key Word: Interval algorithm; global constrained optimization; nonlinear function。