罚函数求解带约束的规划问题

罚函数求解带约束的规划问题（教案）

01级混合八班徐涛3013001231

01级混合八班王菁3013001215

01级混合六班赵晓楠3013001155

§1 求解带约束的非线性规划问题

罚函数法求解带约束的非线形规划问题的基本思想是：利用问题的目标函数和约束函数构造出带参数的所谓增广目标函数，把约束非线形规划问题转化为一系列无约束非线形规划问题来求解。增广目标函数由两个部分构成，一部分是原问题的目标函数，另一部分是由约束函数构造出的“惩罚”项，“惩罚”项的作用是对“违规”的点进行“惩罚”。罚函数法主要有两种形式。一种称为外部罚函数法，或称外点法，这种方法的迭代点一般在可行域的外部移动，随着迭代次数的增加，“惩罚”的力度也越来越大，从而迫使迭代点向可行域靠近；另一种成为内部罚函数法，或称内点法，它从满足约束条件的可行域的内点开始迭代，并对企图穿越可行域边界的点予以“惩罚”，当迭代点越接近边界，“惩罚”就越大，从而保证迭代点的可行性。

1. 1.外部罚函数法(外点法)

约束非线形规划问题

min f(x),

s.t. g(x)>=0,

其中g (x) = (g 1(x),…,gm(x)),

将带约束的规划问题转化为无约束非线形规划问题来求解的一个直观想法是：设法加大不可行点处对应的目标函数值，使不可行点不能成为相应无约束问题的最优解，于是对于可行域S= { x | g(x) >= 0} 作一惩罚函数

P(x) = 0, x∈S;

K, else

其中K是预先选定的很大的数。然后构造一个增广目标函数

F (x) = f (x) + P (x) ,

显然x∈S时，F（x）与f (x)相等，而x S时，相应的F值很大。因此以F(x)为目标函数的无约束问题

minF x) = f(x) + P (x) （1）的最优解也是原问题（NP）的最优解。

上述P(x)虽然简单，但因它的不连续性导致无约束问题（1）求解的困难。为此将P(x)修改为带正参数M（称为罚因子）的函数

P(x) =M ∑[min (0,gj(x))]2

则

min F(x,M) = f(x) + M∑[min (0,gj(x))]2

的最优解x(M) 为原问题的最优解或近似最优解。这时，若x (M) ∈S 则它必定是问题的最优解；若对于某一个罚因子M ，使得 x (M) -∈S ，则加大M 的值，罚函数的“惩罚”作用也将随之加大，因此当 M 是很大的数时，即使x (M) -∈S ，它与 S 的“距离”也不会太远，而且随M 的增大，“距离会越来越近，因此外部罚函数法就是选区一个丹增且趋于无穷的罚因子列

0 < M1 < M2 < … < Mk < …,

从而构成一系列无约束非线性规划问题

min F(x,Mk) = f(x) + Mk∑[min (0,gj(x))]2

2．内部罚函数（内点法）

对于仅带不等式约束的非线性规划问题，也可考虑使用另一种“惩罚”方式。引进的罚函数的作用相当于在可行域的边界上设置障碍，是求解的迭代过程始终在可行域内部进行。由于这种罚函数使得迭代点保持在可行域内部，故称为内部罚函数或障碍函数。

记可行域内部为

S0={ x | g(x) > 0 , j=1, 2, …, m}

且 S0≠? 我们可以仿照外部罚函数法的叠加办法来构造增广目标函数，使得该增广目标函数在可行域内部离边界较远处与原问题的目标函数f(x) 尽可能接近，而在靠近边界是函数之迅速增大

常取

B(x,r) = r ∑ 1/gj(x), (r>0)

或

B(x,r) = r ∑ ln (gj(x)), (r>0)

为障碍函数。在S 的边界上，B(x,r) 为正无穷大。

社选区一旦剪切区域0的“障碍”引子列{ rk} k=1, 2, …,，由每一 rk 作一对应的障碍函数B(x,rk) ，在利用它构造出定义在 S0 内的增广目标函数列

F(x,rk) =f(x) + B(x,rk)

则若点 x(k) 从S0 内向S 的边界趋近时，F(x,rk) 的值将无限增大，由此关于该增广目标函数的无约束问题

min F(x,rk) (1)

得最优解必落在可行域内部，且难以接近可行域边界。若原余额书问题的最优解在内部，则党渠道某一适当值时，无约束问题1的最优解可以达到它。若原问题的最优解在 S 的边界上，则随障碍因子rk 逐渐减小，相应的问题的最优解点烈将向S边界上的问题的最优解逼近。这就是内部罚函数的求解过程。很显然该方法的初始点 x(0) 必须在可行域内部。

§1 求解带约束的线性规划问题

罚函数法还可用于解0-1整数规划问题。0-1整数规划是NP完全的，问题复杂度较高，目前已知的分枝定界法和隐枚举法在需要处理的元素太多时效果并不理想，且不能保证一定能找到最优解，用计算机处理问题涉及矩阵运算时具有较大的空间和实践复杂度，算法还需要改进。在这里我们可以利用罚函数，将整系数多项式0-1规划问题转化为箱约束多项式规划问题，便于用各种数学软件进行算法处理。

考虑一般的0-1整规划形式：

min f(x);

s.t ai≤gi(x)≤bi, i=1,2,……m;

x∈X I

这里f(x),g i(x)都是整系数多项式函数,a i,b i均是整数,X I={0,1}n ;

我们可对目标函数做下述处理，将其表示成约束函数gi(x)的线性组合：

f(x)=

22

2

1

(())(())

[]

()0.5

m

k i

i

g x ai gi x bi

bi ai

=

-+-

-+

∑

记： α= 12||

l i i λ=≤∑ ；

22()1min{,1,2,...}()0.5bi ai i m bi ai β-+==-+；

2

2()max{,1,2,...}()0.5bi ai i m bi ai γ-==-+；

K=max{ln(1)/ln ,1,ln(1)/ln }m αβαγ++-+

; 定义罚函数： Pk(x)=22

21(())(())[]()0.5m

k i i g x ai gi x bi bi ai =-+--+∑; k>K 作如下整理，得到箱约束多项式规划问题：

min f k (x)=f(x)+Pk(x);

st. x ∈X n

其中X n ≡［0,1］n ；

可以证明，若x *=(x *1,x *2,…,x *n )T 是该问题的的最优解,不失一般性,设x i *∈(0,1),i=1,2,…,t,x i *∈{0,1},i=t+1,t+2,…,n.

设

=(1

,…,t

,

,…,)T ,其

中i ∈{0,1}(i=1,2,…,t).那么,也是此问题的最优解。

§3 应用举例

下面我们应用罚函数方法来解决一个实际问题。试考虑如下情形的飞行管理策略： 在约10000米高空的某边长为160公里的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘的时候，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角度，以避免碰撞。现假定条件如下：

（1） （1） 不碰撞的标准为任意两架飞机的距离达于8公里；

（2） （2） 飞机飞行方向角调整的幅度不应超过30度；

（3） （3） 所有飞机飞行速度均为每小时800公里；

（4） （4） 进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60公

里以上；

（5） （5） 最多需考虑6架飞机；

（6） （6） 不必考虑飞机离开次区域后的状况。

请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01度），要求飞机飞行方向角调整的幅度尽量小。

设该区域4个定点的坐标 (0,0)，(160,0)，(160,160)，(0,160)。记录数据为：

注：方向角指飞行方向与x轴正向的夹角。

这是一个有约束最优化问题。初步分析后可以发现，约束条件是非线性的。既要求较高的精度，又要求在尽可能短的时间内给出解答。在这里，综合精度的要求和技术的可实现性，我们将其转化为无约束的非线性规划，从中可以看到罚函数的重要作用。

首先，我们做出如下假设以简化问题：

（1）（1）飞机进入控制区域后完全服从地面控制台的调度，飞机未接到指令时保持飞行状态不变；

（2）（2）飞机接到指令后可立即转到所需的角度，即不考虑转弯半径的影响，调整在瞬时完成；

（3）（3）飞机在区域内至多调整一次方向；

（4）（4）已在区域内的飞机已经调整完全，不会相撞；

设0i a为第i架飞机的初始方向角，i a为第i架飞机的方位角，(xi,yi)是飞机的坐标，可以用i a的正弦函数来表示，d ij（t）表示时刻t时i,j两架飞机的距离，则问题的目标函数和约束条件可以直接的表述如下：

min f=

2

1

()

n

i i

i

a a

=

-

∑

(n≤6)

s.t. min d ij2≥64 (i,j=1,2,……6,i≠j) t>0

│i a│≤π/6

其中，d ij2=

实际上，由于每架飞机飞过该区域的事件不可能超过0.28小时（即飞正方形对角线时间），所以我们可以将t的上限设定在0.28小时，以最大限度的减少计算。

由于题目所涉及数据变量不是太多，可以先考虑用逐步求精的直接搜索法来求解，MATLAB软件也提供了相应的函数可以方便的调用。这种方法每次用一定的步长以较少的循环次数进行“粗选”，在“粗选”出的几个解的附近一间小的步长进行“精选”，逐次推进，直到达到所需精度。为了控制计算时间，还可以采用以下的优化方法：

（1）（1）将底层循环内判别相撞的函数分散设在每层循环下，使在高层发现相撞后可提前结束循环；

（2）（2）进入新一层循环后以积累的偏差平房与已得最小偏差平方和进行比较，若大则结束该层循环。

这些措施可以大大平均搜索次数，节省运算时间。就题中特例，该算法用MA TLAB解决，耗时约为6-7秒。所的结果为（为了与后面的方法作比较，保留了较高的精度）：

容易检验，这里给出的的确是质量很高的最优解。可是算法的耗费时间却并不令人满意。即便如此，这种方法也有其可用之处。它能在较短的时间内给出一个较为接近最优解的可行解，从这个可行解出发，我们可以构造相应的罚函数较快地得出满足精度要求的最优解。

为了使求导等计算更加简便，我们对目标函数和约束条件作一些形式上的修正：

min F （X ）=22,0,01()n i i i i i C S =+∑（n ≤6）

s.t. g ij （X ）=min d ij 2-64≥0 (i,j=1,2,……6,i ≠00j) t>0

│i a │≤π/6

其中，,0i i C =0cos cos i i a a -, ,0i i S =0sin sin i i a a -

X= （1,10C ，2,20C ，……6,60C ，1,10S ，2,20S ，……6,60S ）

此时dij 也用,0i i C 和,0i i S 来表示.

构造罚函数

P （X k ，M k ）= F （X k ）+M k 2

()1[min(0,)]n ij X i g =∑ (n ≤6)

其中,权因子M k 是一单调递增趋向于∞的序列。对每个M k 值求取相应的X k ，使P 取得极小值。设精度要求为ε，当1k k X X --<ε时结束运算。X k 即为所求最优解。

考虑到精确性要求和运算的便利，我们取M 1=1，M k = k-1，ε=0.5*10-2。

我们使用求解无约束规划问题的经典算法SUMT 来具体处理题中所给的数据记录，初始值由短时间的直接搜索所得的近似解带入，可得结果：

可见，只有第二架飞机的调整角度有了较为显著的变化，但是这种显著是由于我们保留了较高的精确度引起的。绝对变化0.0045弧度远远小于我们的精度要求0.01度。而所需时间则降至2秒以下了。

通过这个事例我们可以看出，在规划中，罚函数方法是一种富有效率的求解手段，在保证精确度的同时，能够大幅度的降低运算的时间复杂性，因此，它在实际操作中得到了广泛的应用。

参考： 高峰 张连生 《多项式0－1整规划的两个连续化途径》

万中 曾金平 《数学实验》

上海大学学报1999年第2期