第六章 约束规划.doc

第六章 约束优化方法
工程实际优化问题绝大多数属于约束非线性规划问题，其一般数学表达式为：

求解方法可分为：直接法和间接法。两类方法的特点见表1

分类
举例
优点
缺点
应用场合
直接法
坐标轮换法
随机方向搜索法复合形法
算法简单、直观
计算量大，收敛慢
维数低，函数复杂，精度要求低的优化问题
间接法
简约梯度法
惩罚函数法
序列二次规划
收敛快，整体计算量小
公式较复杂
大型优化问题

6.1 约束随机方向搜索法
6.1.1 基本原理
在可行域内选一初始点X0，以一初始步长沿一随机方向S1，求得探索点X。S1应保证X在可行域内且使目标函数值下降，即新点应有可行性和下降性。改变步长，继续在S1方向上探索，得到S1方向上的最优点X1。
以X1为初始点，在另一随机方向S2上重复上述过程，得到S2方向上的最优点X2。
如此重复下去。
当某一成功点X*沿着N（N=50~500）个随机方向的探索均失败时，以X*为最优解。
6.1.2 初始点的选取
要求初始点可行点。当约束条件比较简单时，可人为确定；当约束条件复杂时，人为方法比较困难，可随机选取，即利用计算机产生的伪随机数来生成初始点。具体方法为：
设设计变量的分量xi在取值范围为区间[ai,bi]，qi为区间（0，1）内的随机数，则xi的随机数为

由此可得到X所有分量随机数，然后将X代入约束条件中检验，若满足所有约束条件，则X是可行点。否则应重新选取初始点。
6.1.3 随机搜索方向的产生
以二维问题为例，说明随机搜索方向的产生方法。若y1、y2为区间[-1，1]上两个随机数，则向量[y1,y2]可为平面内的任意方向。取向量[y1,y2]的单位向量e

则e的端点位于单位圆的圆周上。对于三维问题，单位向量e的端点位于以单位长度为半径的球面上。
上式中要求y1、y2在区间[-1，1]内随机取值，若计算机只能产生[0,1]区间内的伪随机数ri，可用下式将其转换为[-1,1]区间的伪随机数：

6.1.4 举例
例：用约束随机方向搜索法求解

解：人工选取一初始点X0=[5，5]T，初始点在可行域内。相应的目标函数值为F(X0)=50。
第一次迭代
1)产生两个伪随机数，求出第一个随机方向。
生成两个伪随机数r1=0.9, r2=0.2

由此得第一个随机方向为：
2)求第一个迭代点
第一个迭代点表达式为：
式中a为步长。
将X1的表达式代入目标函数中，进行一维搜索，令目标函数对步长a的一阶导数为0，即可求出沿e1方向的最优步长。

第一个迭代点为：
3）检验X（1）点是否满足约束条件
g(X(1))=2╳4.2+5.6=14>12，X（1）满足约束条件，是可行点。相应的目标函数值为:

第二次迭代
以X1为初始点，重新生成随机搜索方向e2。
（以下过程略）

6.2 复合形法
6.2.1 复合形法的基本原理
在n维空间的可行域中选取k个设计点（通常取n+1≤k≤2n）作为初始复合形的顶点，然后比较复合形各顶点目标函数值的大小，目标函数值最大的点为坏点，以坏点之外的点的中心为映射点为中心，求出坏点的映射点。一般映射点优于坏点，以映射点代替坏点，构成新的复形。如此循环，使复合形不断向最优点移动和收缩。当收缩到复合形的各个顶点与形心非常接近、满足迭代精度要求时为止。最后输出复合形顶点中目标函数值最小的顶点为近似最优点。

如图二维约束优化问题，n=2, 取4个顶点构成初始复合形，设这4个点中目标函数值最大的点为X1 ，记为XH，目标函数值最小的点为X3，记为XL。丢弃XH点，求另外3点的中心XC,连接XC和XH，在连线上确定出映射点XR。

式中α为映射系数，一般α>1，通常取α=1.3。

若XR在可行域内，且XR的目标函数值优于XH，则以XR替换XH形成新的复合形，重复上述过程。
若XR不满足上述两个条件，则将映射系数α减半，α可以多次减半。若α变得很小仍不能使XR满足条件，可以次坏点代替坏点进行映射。

在求映射点之前，若XC在可行域之外，则可行域可能是一个非凸集，这时可由XC与XL构成一个超立方体（二维情况下为长方形），在该区域内重新利用伪随机数产生k个顶点，构成新的复合形。

6.2.2 初始复合形法的产生
对于维数较低的简单优化问题，可以人为给出初始复合形顶点。
对于复杂优化问题，多采用随机方法产生初始复合形。过程如下：
（1） 利用随机函数，生成第一个顶点，确保该顶点在可行域内。

（2） 产生其他（k-1）个随机点

若X2、X3、…、Xq都在可行域内，则它们可作为初始复合形的顶点。若Xq+1点不在可行域，如图，

先求出已在可行域中的q个顶点的中心XD，然后在XD 与Xq+1连线上移动Xq+1点，使Xq+1到XD的距离减半，即：

若推移后Xq+1点进入了可行域，则将Xq+1点作为初始复合形的第q+1个顶点，否则继续缩短Xq+1到XD的距离。
6.2.3复合形法举例
例：用复合形法求解约束优化问题，迭代精度为0.01。

解：n=2,取复合形顶点数k=2n=4。
（1） 人为给出4个复合形顶点。

经检验4个顶点均在可行域内。
（2） 进行迭代，获得新的复合形。
求4个点的目标函数值：

知.
计算除坏点外所有顶点的中心XC：

经检验XC在可行域内。
求XH的映射点XR，取映射系数为1.3：

经检验，映射点XR在可行域内。
因F(XR)=5.1092< F(XH)，故用映射点替换X3，构成新的复合形：

4个点的目标函数值为：
最坏点为X4，最好点为X3。
（3） 检验迭代终止条件
常用各顶点与好点的目标函数之差的均方根小于误差作为终止迭代条件，即

若满足以上条件，XL为最优解。对于本题，

不满足迭代终止条件，需继续迭代。在新的复合形中重复上述过程……
6.3 惩罚函数法
惩罚函数法的基本思想是：根据约束条件构造惩罚函数，并将其加到目标函数中，将约束问题转化为无约束问题。
惩罚函数法可分为外罚函数法、内罚函数法和混合函数法。
6.3.1 外罚函数法
6.3.1.1 外罚函数法基本原理
该法的主要特点是惩罚函数定义在可行域外部，求解过程中从可行域的外部逐渐逼近原来有约束优化问题的最优解。它既可用于求解不等式约束优化问题，也可用于求解等式约束问题。先介绍这种方法对不等约束问题的处理。
设原约束优化问题为：

构造一般形式的外罚函数为

式中第二项为惩罚项，其中γk是惩罚因子，在迭代过程中γk的取值会越来越大取值，即γk是一个递增的正数列：

由此可知，当X在可行域内和约束边界上时，无论γk取何值，惩罚项均为0。
当X违反约束条件时，无论γk取何值，惩罚项不为0，表明X在可行域外面，且X离约束边界越远，惩罚项的值越大。
当γk由小变大时，迭代点必需由边界外向可行域靠拢才能使目标函数有最优解。
例：用外罚函数法求解下列约束优化问题

解：构造惩罚函数，将原问题转化为无约束优化问题，即：

选择一个适当的初始外点X0和一个初始惩罚因子γ0，用无约束优化问题求解方法（如牛顿法、变尺度法等）解出与之对应的最优值，然后令γk不断增大，当γk趋于无穷大时所得X即为最优解。
但本题比较简单，可用解析法求解,即令
得到：

当γk趋于无穷大时，X=1，P=1.此二值即为最优解和最小值。

若优化问题中约束条件为等式约束，

也可作类似处理，得到一个罚函数

综合不等式约束和等式约束，得到一般约束优化问题的外罚函数表达式：

6.3.1.2 约束裕量问题

实际计算中惩罚因子γk不可能达到无穷大，因此最后所得的最优点也就不可能收敛到原问题的最优点，而是落在可行域的外部，这就不能严格满足约束条件。为克服外罚函数法的这一缺点，引入约束裕量，将可行域边界向内收缩一个微量，即重新定义约束条件：

式中即为约束裕量，不能取值过大，以免所得结果与最优点相差太远，一般取。

6.3.1.3 外罚函数法的迭代步骤
（1） 给定初始条件：初始点，初始惩罚因子，迭代精度，惩罚因子递增系数。
（2） 用无约束优化方法求罚函数的极小点
（3） 用点距准则或函数值下降准则检验相邻两次迭代中的极小点是否满足迭代终止条件。若不满足，则增加惩罚因子，用当前极值点作为初始点，重新迭代。
6.3.1.4 外罚函数法应用中的问题
（1） 初始点可以在可行域内，也可以在可行域外，对实际计算很方便。
（2） 初始惩罚因子和递增系数的选取对方法的有效性和收敛速度都有影响。取值过小，会增加迭代次数，取值过大，罚函数的性态会变差，导致求解困难。
6.3.2 内罚函数法
特点为：只适用于不等式约束，惩罚函数定义在可行域的内部，每一迭代点都在可行域内部移动，并逐渐逼近原约束优化问题的最优解。（迭代过程的每个一设计点都对应一个可行方案）
内罚函数可构造为：

惩罚因子是一个递减的正数数值，即：
例：用内罚函数法求解下列约束优化问题

解:构造内罚函数为

用解析法求极值

6.3.3 混合罚函数法
混合罚函数法是内罚函数法和外罚函数法的综合。对于p个等式约束，构造外罚函数，对m个不等式约束，构造内罚函数，得到混合罚函数为

式中，使用统一的惩罚因子，它是一个递增数列：

例：求下列优化问题的混合罚函数

罚函数为

式中，惩罚因子是一个递增数列：。
罚函数法原理简单，算法易行，适用面广，可以和各种有效的无约束最优化方法结合。但也存在问题：从理论上讲，只有当惩罚因子趋于无穷时或0时，算法才收敛，因此收敛较慢。另外，当惩罚因子的初始值选取不当时，惩罚函数可能变得病态，使计算发生困难。
6.4 增广乘子法
6.4.1 拉格朗日乘子法
拉格朗日乘子法是一种古典的求约束极值的间接法。
设有具有等式约束的优化问题：(不等式约束也可以通过引入松弛变量转换为等式约束)

得到拉格朗日函数:
式中λ为拉格朗日乘子，为变量，故拉格朗日函数共有（n+p）个变量。为求L(X, λ)的极值，令▽L(X, λ)=0，得到(n+p)个方程：

拉格朗日乘子法看似简单，但实际存在许多问题：
其一，当L(X, λ)为非凸函数时，求得的最优点可能是鞍点，因为海塞矩阵H不一定为正定；
其二，对于大型非线性优化问题，(n+p)个方程为高次方程组，用数值方法求解(n+p)个方程的难度几乎与求解优化问题本身的难度相当；
其三，还必须分离出程组的重根。
因此，拉格朗日乘子法不是一种有效的求解约束优化问题的方法。
6.4.2 增广乘子法
增广乘子法来源于拉格朗日乘子法，是拉氏乘子法与罚函数法的综合。其搜索策略与外罚函数相仿。具体做法是，用外罚函数替代原目标函数，原约束条件不变，得到一个增广极值问题，然后构造该问题拉格朗日函数。
原问题为：

增广极值问题为：

可以证明，增广极值问题不仅与原问题同解，而且与原问题有相同的拉格朗日乘子。
增广拉格朗日函数：

当r足够大时，L的海塞矩阵为正定，L必有极小值。但确定r的临界值较困难，实际求解时是采取交替修改惩罚因子r和拉格朗日乘子λ，经反复迭代得到足够大的r和最优的λ。迭代公式的来历如下：
对增广拉格朗日函数求极值：

若对原问题采用拉格朗日函数求极值，得：

比较以上两式，得到迭代公式：

式中，惩罚因子是一个递增数列：,一般取为等比数列，比例系数为c’=10。
令,Powell等人证明，当r足够大时，必能满足

这时r不需要再增大。r一般取0.25。为方便起见，以后将此比值固定为0.25。
迭代收敛准则为
例：用增广乘子法求下列规划问题的最优解，取精度为0.01。

（本题最优解为 ，拉格朗日乘子为3。）
解：构造增广拉格朗日函数

可知，其海塞矩阵为，故当r>2时，H为正定矩阵，故L有最优解。
（1） 取r=8, λ=0,代入增广拉格朗日函数

求其极值,令：

得：

检验是否满足条件,取的初始值为一较大值，如10。则有

故不需要增大罚因子，仍保持为r=8。
现求λ，由迭代公式

得：
(2)将λ=4，r=8代入L函数中

求L的极值得：

条件不满足，故增大罚因子，即r=8*10=80
（3） 将λ=4，r=80代入L函数中

求L的极值得：

条件满足，故不需要增大罚因子
现求λ，由迭代公式

得
将新的λ和r代入L函数中

求L的极值得：
,满足收敛条件，迭代终止。
最优解为，f(X)=-0.000986。

6.5 二次规划
二次规划是指目标函数为二次型，约束条件为线性式的规划问题，即：

设式中变量数为n，约束数为m。

回顾库恩－塔克条件：
min f(x), x∈Rn
s.t. gu(x) ≤0 (u=1,2,...,m)
hv(x) =0 (v=1,2,...,p)

如果x*是一个局部极小值点，则该点的目标函数梯度可表示成该点诸约束面梯度的线性组合。即：
(2-13)
式中q为X*处不等式约束面的数目；j为X*处等式约束面的数目。
为非负乘子，也称拉格朗日乘子。若某个约束不起作用（即X*不在约束面上），则对应的为0.

对于二次规划，K-T条件为：

式中λ、μ分别为约束条件、变量的拉格朗日乘子，y为松弛变量。
上式中的一、二式构成了一个线性方程组，可由单纯形法求解。但没有目标函数，可以在方程中引入伪变量，并构造一个伪目标函数。使伪目标函数为0，并将伪变量变成非基本变量，即得到方程的解。三、四两式表明，不能同时进入基变量，也不能同时进入基变量。
例：用单纯形换基法求二次规划：

解：将原问题改写成

则C=[-8,-10], D=[2 0;0 2],A=[3 2],b=6

由

得：

即K-T条件为：

前三个方程为线性方程，变量个数为5，可用单纯形法求解。后三式表示相乘的两个变量不能同时进基。
引入人工变量和伪目标函数

建立单纯形表：

将of行中与z1、z2对应的系数1、1变换为0：

将of行中与r对应的系数-5变换为0，在该列生成一个基变量，即r列进基，z1列出基：

x2列进基，z2列出基

读出最优解
clc;
syms x1 x2 r u1 u2 y z1 z2 f
head=[x1 x2 r u1 u2 y z1 z2 f];
of=[0 0 0 0 0 0 1 1 0];
l1=[2 0 3 -1 0 0 1 0 8];
l2=[0 2 2 0 -1 0 0 1 10];
l3=[3 2 0 0 0 1 0 0 6];
[head;of;l1;l2;l3]
of=of-l1-l2;
[head;of;l1;l2;l3]
l1=l1/3;of=of+l1*5;l2=l2-l1*2;
[head;of;l1;l2;l3]
of=of+l2;l3=l3-l2;l2=l2/2;
[head;of;l1;l2;l3]
l3=l3*3/13;l1=l1-l3*2/3;l2=l2+l3*2/3;
[head;of;l1;l2;l3]

用matlab优化工具箱函数检验：
H=[2 0;0 2];C=[-8;-10];A=[3 2;-1 0;0 -1];b=[6;0;0];
x=quadprog(H,C,A,b)

序列二次规划：对于目标函数及约束条件均为非线性的规划问题，可将其按泰勒级数在迭代点展开成二次规划问题，解此二次规划，得到下一个迭代点。由于涉及到海塞矩阵，可采用变尺度法逼近海塞矩阵的逆。
序列线性规划：同样对于上述非线性的规划问题，可以在迭代点将其展开成线性规划问题，利用单纯性求出新的迭代点。
用以上两种方法求解非线性规划，存在两个问题：
1）由于线性化带来的误差，近似最优解可能落在原问题的可行域外，甚至与真正的最优解相距很远。当初始点离真实的最优解较远时，计算可能不收敛。
2）非线性规划的最优解不一定在可行域的顶点上，有时可能在约束曲面的某一点A上。在此情况下，采用序列二次规划或序列线性规划得到的解可能在与A点相邻的两个顶点之间来回振荡。
为克服上述现象，可以采用限步法。即定出最大移动量，人为限制近似最优解的移动范围，使每次计算都向有利方向前进一步，但步子比较小。此举相当于增加了约束条件，且约束条件为线性的。