DFP变尺度法

格式：docx
大小：16.14 KB
文档页数：6

下载文档原格式

/ 6

约束变尺度法

约束变尺度法Newton 法最突出的优点是收敛速度快，在这一点上其它算法无法比拟的。

因此，建议凡是Hesse 矩阵比较容易求出的问题，尽可能使用Newton 法求解。

但是，Newton 法也有一个严重缺陷，就是每次迭代都要计算目标函数的Hesse 矩阵和它的逆矩阵，当问题的维数较大时，计算量迅速增加，从而就抵消了Newton 法的优点。

为此，人们开始寻找一种算法既可以保持Newton 法收敛速度快的优点，又可以摆脱关于Hesse 矩阵的计算，这就是变尺度算法。

变尺度法是一种非常好的方法，其中DFP 算法和BFGS 算法。

可以说，直到目前为止，在不用Hesse 矩阵的方法中是最好的算法。

一、拟Newton 法为了吸收Newton 法收敛速度快的优点，同时避免Newton 法每次迭代都要计算目标函数的Hesse 矩阵和它的逆矩阵，人们提出了具有超线性收敛的拟Newton 法。

（一）拟Newton 法的基本原理在Newton 法中的基本迭代公式kk k k P t X X +=+1，其中1=k t ，)()]([12kkk Xf Xf P ∇∇-=-令)()(2kkkkXf gXf G∇=∇=，于是有，，，，21011=-=-+k g G X X k k k k其中X0是初始点, gk 和 Gk 分别是目标函数f (X )在点 Xk 的梯度和Hesse 矩阵.为了消除这个迭代公式中的Hesse 逆矩阵G-1k ,可用某种近似矩阵Hk=Hk(Xk)来替换它,即构造一矩阵序列{Hk}去逼近Hesse 逆矩阵序列{G-1k}，此时kk k k g H X X -=+1事实上,式中 Pk= -Hk gk 无非是确定了第k 次迭代的搜索方向.为了取得更大的灵活性,考虑更一般的迭代公式kk k k k g H t X X -=+1其中步长tk 通过从Xk 出发沿Pk= -Hk gk 作直线搜索来确定.此式代表很广的一类迭代公式.例如,当Hk=I (单位矩阵)时,它变为最速下降法的迭代公式。

第七节变尺度法

5）检查迭代次数，若k=n(完成一轮迭代仍未能满足精度要求)，则令 = 进行下一轮迭代计算；若k n,则进行下一步；
6）构造新的迭代方向
五、DFP法的特点
1）DFP法具有牛顿法的二阶收敛速度，因变尺度矩阵最终逼近。
2）对任意给定的初始点，都具有最速下降方向 =-▽f( )。这是由于初始变尺度矩阵 =I(单位阵)，所以搜索方向 =- ▽f( )=-▽f( )。
3）每次搜产生的方向都是共轭的。
4）计算公式具有递推性，便于迭代。只要给出，，即可往下计算。
5）计算方便，无需计算二阶导数矩阵及其逆矩阵。
六、BFGS法
本章学习要点
1）首先要抓住并理解各种优化方法的基本思想、具体的迭代步骤，读懂计算程序框图；其次是较复杂的数学推导。对于工程技术人员来说，主要是应用。因此，对每种优化方法书中都有较简单的数学模型的例题，用手算的迭代过程，读者将例题与计算程序框图结合起来，认真地阅读理解。有条件最好再上机运算，在实践中进一步理解每种优化方法的迭代逻辑。
一、基本思想
二、构造变尺度矩阵的基本要求
三、校正矩阵 de推导
四、变尺度法的迭代步骤
1）给定初始点、收敛精度ε和维数n;
2）置k=0, = =I（n×n阶单位阵），计算梯度 =▽f( ),搜索方向 =- ▽f( )=- ；
3）进行一维搜索求最优步长因子，迭代点：
minf( + = +
4）收敛精度判断
第七节
变尺度法是在克服了梯度法收敛速度慢和牛顿法计算量、存储量大的特点基础上而发展起来的，被公认为是求解无约束优化问题最有效的算法之一，现已在工程优化设计中得到了广泛的应用。变尺度法的种类较多，这里只介绍其中最常用的两种方法：DFP法和BFGS法。DFP变尺度法先由戴维顿（Davidon）于1959年提出，后经费莱彻（Fletcher）和鲍威尔（Powell）于1963年改进而来，故又称DFP法。

第四节无约束--DFP变尺度法6

2.6、DFP变尺度法：变尺度法是Davidon 由于1959年提出后又经Fletcher和Powell加以发展和完善了后的一种变尺度法，故称DFP变尺度法
1、基本思想：变尺度法是克服了梯度法收敛慢和牛顿法计算量大的缺点而发展起来的，是求解无约束问题最有效的算法，在工程优化设计中得到了广泛的应用。利用牛顿法的迭代公式，然后并不是直接计算 [ H ( x )] 而是用一个对称正定矩阵 G( x )近似 (k ) 地代替 [ H ( x )] , G( x ) 在迭代过程中不断改进。
⑥检查迭代次数若若 k n 则转(7)
( k 1)
k n
置
( k 1)
x
x
则转(2)
(k )
⑦计算
f ( x
( k 1)
)g
(k )
( k 1)
g
g
(k )
g
(k )
(k )
x
x
x
E
(k )
G
(k )
E
(k )
G
( k 1)
然后置 k 1 k转(3) G
说明: 已有文献证明：对于二次函数 f ( x)DFP 变尺度法所构成的搜索方向 S S S S 为一组关于海色矩阵共轭的矢量。所以DFP变尺度法属于共轭方向法，具有二次收敛性,在任意情况下，这样方法对于二次目标函数都将在几步内搜索到目标函数的最优点，而是最后的构造矩阵必等于海色矩阵.
0 1
计算 x
点的导数值
S
( 0)
( 0)
g
( 0)
f ( x
( 0)
8(8 5) 24 ) 2(9 6) 6

DPF变尺度法理论证明

( j)
H k +1 Ap
( j)
( j)
= H k Ap
(k )
( j)
+
p p p p( j )
(k )
( k )T
Ap
(k )
( j)
( k )T
q
−
Hkq q q
( k ) ( k )T ( k )T
H k q( k )
=p
−
Hkq
( Ap
(k ) T (k )
)
q
( k )T
Hkq
= p( j )
(k ) ( k )T (k )
∆x ∆x H k ∆g ∆g H k H k +1 = H k + ( k ) T ( k ) − ( k )T (k ) ∆x ∆g ∆g H k ∆g
(k ) ( k )T (k ) ( k )T
p p H k q q H k (k ) H k +1q = H k + ( k )T ( k ) − ( k )T q (k ) p q q Hkq ( k ) ( k )T ( k ) ( k ) ( k )T (k ) p p q Hkq q Hkq (k ) =H k q + − =p ( k ) ( k )T ( k ) ( k )T (k ) p q q Hkq
( i )T
p
( k )T
Ap
( j)
= λk d
( k )T
Ap
( j)
= −λk ( H k g k ) Ap
T
( j)
T T T = −λk g k H k Ap ( j ) = −λk g k p ( j ) = −λk λ j g k d ( j ) = 0.

DFP变尺度法例题终结版

(0)
) = (4,100) , P
T
(0)
= − H ( X (0) )∇f ( X (0) ) = (−4, −100)T
_
1.求迭代点X (1) , 令ϕ0 α）=f ( X (0) + α 0 P (0) )得ϕ0 α 0）的最小值点α 0 = 0.02 （（
= X
(1) (0)
+ α 0 P ( 0 ) = (2, 2) T + 0.02( − 4, − 100) T = (1.92, 0) T
表2 DFP变尺度法结果表
K αK
0 1 2 0.02 0.5
X1
2 1.92 0
X2 2 0 0.0015 5
∂ X K ∂X 1
_
∇f (X
_ K
_ K
f (X )
4 3.84 0
100 0
100.2 3.84
104 3.686
0.0775 0.07 75
6.006×10−5
下一个搜索方向：P = − H( X )∇f ( X ) = [ −3.84, −0.0031]
(1) (1) (1)
_
T
2.求迭代点 X ，令 ϕ1 = f ( X (1) + ∂1 P (1) )
(2)
得 ϕ(∂1) 的最小值点∂ 1
1
=0.5
T
[ 0, − 0 .0 0 1 5 5 ] T (2) ∇f ( X ) = [ 0, −0.0775]
∆X(0) (∆X(0) )T H(X(0) )∆g(0) (∆g(0) )T H(X(0) ) H(X(0) ) = H(X(0) ) + − − (0) (0) (∆g )∆X (∆g(0) )T H(X(0) )∆g(0)

【MATLAB与机械设计】多维优化之DFP变尺度法

【MATLAB与机械设计】多维优化之DFP变尺度法1.程序框图其中E矩阵的求解公式为：2.MATLAB可执⾏程序function [k,x_min,f_min]= variable_metric_algorithm(f,x,exp)%% 程序说明%{该程序应⽤于多维⽆条件优化问题中的变尺度法变量说明：输⼊值部分f为⽬标函数，对于⽬标函数中⾃由变量的个数没有要求x为初值矩阵，要求在调⽤函数是必须为⼀⾏n列的形式，否则会导致后期维度出现错误exp为精度返回值部分k为迭代次数x_min为函数的局部极⼩值数列f_min为函数的局部极⼩值调⽤⽅法：clcclearsyms x1 x2%所有的⾃由变量且必须由x1,x2,x3……表⽰x=[0,0];f=x1^2+x2^2-x1*x2-10*x1-4*x2+60;exp=0.01;[k,x_min,f_min]=variable_metric_algorithm(f,x,exp)%}%% 确定维度x_size=size(x,2);x_i= sym(zeros(1,x_size));%class(x_i)for i=1:x_sizesyms(['x' num2str(i)]);x_i(1,i)=['x' num2str(i)];end%% 迭代主题grad_f=gradient(f,x_i);H=eye(x_size);fz_d=-subs(grad_f,x_i,x);fz_d=double(fz_d);g0=fz_d;k=0;while 1s=-H*g0;syms as=s';f_b=subs(f,x_i,x+a.*s);f_bd=diff(f_b,a);a=solve(f_bd==0,a);a=double(a);x_1=x+a*s;fz_d1=-subs(grad_f,x_i,x_1);fz_d1=double(fz_d1);fz_dm1=norm(fz_d1);g1=fz_d1;if fz_dm1<=expx_min=x_1;f_min=subs(f,x_i,x_min);k=k;break;elseif k==x_sizefz_d=-subs(grad_f,x_i,x_1);fz_d=double(fz_d);k=0;H=eye(x_size);elsedeta_x=x_1-x;deta_f_d=fz_d1-fz_d;A=(deta_x')*(deta_x);sub_A=(deta_x)*deta_f_d;B=H*deta_f_d*(deta_f_d')*H; sub_B=(deta_f_d')*H*deta_f_d; E=A./sub_A-B./sub_B;H=H+E;k=k+1;endx=x_1;g0=g1;endendend。

变尺度法ppt课件

满足上述方程的解很多，我们可以如下确定一组解：
k uk uTk gk xk kvkvTk gk Hk gk
这样，我们可以取：
uk xk ,
k uTk gk 1,
vk Hk gk , kvTk gk 1。
10
即：
uk xk , vk Hk gk ,
H I时梯度法, k
最速下降方向 d k f ( xk ) , 度量为 x xT I x
Hk Gk1时 Newton法最速下降方向d k Gk1f ( xk ), 度量为 x xTGk1 x
也称拟Newton法为变尺度算法。
4
3. 如何对H 附加某些条件使得： k （1）迭代公式具有下降性质 Hk （ 0 正定）
牛顿法特点
|| ||
xk1 x xk x*
* || ||2
,
0.
▪ 收敛速度快，为二阶收敛。 ▪ 初始点要选在最优点附近。
存在缺点及修正 (1) f ( xk1 ) f ( xk ) ?
xk 1
xk
Gk1
g
，
k
g
为梯度，
k
Gk1海塞阵的逆
(2) 初始点的选取困难，甚至无法实施。 (3) G1的存在性和计算量问题。
k
Step3. 令 xk 1 xk kd k , 其中tk : f ( xk kd k ) min f ( xk d k )。
7
Step 4. 判断 xk 1 是否满足终止准则： yes: 计算 stop, 则x* : xk1
No : 转 step 5 。
Step5. 令 gk1 f ( xk1), gk f ( xk ), gk f ( xk1) f ( xk ) gk1 gk , xk xk1 xk 。

张飞---变尺度法—

H k 1 H k
T H k g k g T x k (x k ) kHk k T (x ) g k g T k H k g k
DFP算法的性质：（1）校正公式的分母总大于零，校正公式总有意义。（2） H k都是正定的，所以每次迭代方向都是下降方向。（3） H n A 1 ，说明H k 逐次逼近 A 1 。
2
k 1
k
H kBFGS Hk 1
k T T x k g T H H g ( x ) x k (x k ) k k k k (1 ) T T (x k ) g k (x k ) g k (x k )T g k
g T k H k g k
变尺度法的计算步骤：
一般函数举例： f(x ) e DFP：初值
X0=[0.1,0.1,0.1] X0=[0.5,0.5,0.5] X0=[1.0,1.0,1.0]
2 2 2 0.01x1 0.2x 2 x3
e
2 2 2 x1 0.1x 2 x3
2 2 (sin( x 12 1.1x 2 x3 ))2
变尺度法，也称拟牛顿法。它既吸收了共轭梯度法计算量小，适合求解大规模优化问题的优点，又吸收了牛顿法收敛速度快的优点。它被认为是无约束优化问题中最有效的方法之一。
变尺度法的思想： 2 k 1 p f ( x ) gk 牛顿方向： k
近似为： p k H k g k
H k 称为尺度矩阵，且迭代过程中 H k是变化的，该算法为变尺度法；
g k ：梯度（一阶导数） Ak ：Hesse矩阵（二阶导数）
x2
X0 SN
-S
x1
变尺度法
牛顿法的最大优点是靠近极小点时收敛速度快，主要缺点是计算比较麻烦，对于n维问题，要计算1/2n(n+1)个二维偏导数，还要计算Hesse矩阵的逆矩阵，计算量大。

变尺度法

变尺度法 (拟牛顿法)
一.尺度矩阵的引入
尺度变化可以降低函数的偏心程度（通过放大或者缩小各个坐标）
若f (x) 1 xT Gx bT x c，G是对称正定阵，则G与单位阵E合同，故存在一个 2
可逆矩阵Q, 使得QT GQ=E,即函数偏心程度为0 QT Q=G 1 xTGx 1 xTQTGQx
满足上述方程的解很多，我们可以如下确定一组解：
kukukT yk sk k vk vkT yk H k yk
这样，我们可以取：
uk sk ,
kukukT yk sk ,
vk H k yk , k vk vkT yk H k yk ;
8
即：
uk sk , vk Hk yk ,
计算H
，令
k
k

k
1，转step
3.
10
四、变尺度法算法框图：
k k1
给定x0 , H 0
, k 0
d k -H kf (xk )
一维搜索得 k xk 1 xk k d k
修正Hk
N
产生H k 1
xk1 xk
Stop y xk1 解
11
4
(1)d k

H
k
gk为下降方向，则d
k
与梯度方向成钝角，即gT k
d
k

0

gT k
H
k
gk

0,
gT k
H
k
gk

0
H k 对称正定
（3）
f
(x)

1 2
xT Gx

bT
x

第2章优化设计-2

（2-45）
上式中的搜索方向 S (k)
称为牛顿方向，
可见原始牛顿法的步长因子恒取： (k) 1 ，因此，原始牛顿法是一
种定步长的迭代过程。
牛顿算法对于二次函数是非常有效的，迭代一步就可达到极值点，而这一步根本不需要进行一维搜索。
对于高次函数，只有当迭代靠近极值点附近，目标函数近似二次函数时，才会保证很快收敛，否则也可能导致算法失败。
(2-57)
但其中的校正矩阵的计算公式为
A(k )

X
1
(
k
)
T

g (k )
{X
(k)

X
(
k
)
T

X (k)
X
(
k
)
T

g
(k
)
T

A(k )g (k )
X
(k)
T

g (k )

A(k )g (k )

X
(k
)
T

X (k )
g
BFGS法的优点:
在于计算中它的数值稳定性强，所以它是目前变尺度法中最受欢迎的一种算法。
2.5 约束优化方法
工程中的大量优化设计问题，都是约束优化问题，这类问题的一般数学模型为：
min f ( X ), X Rn s.t.gu ( X ) 0u 1, 2, hv (X ) 0v 1, 2, ,
[X (k ) ]T g (k )
[g (k ) ]T A(k )g (k )
（2-54）（2-55）
式中: X (k) X (k1) X (k)
第 k 次迭代中前后迭代点的向量差；

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
g1=g2; k=-1;
}else{
int j;
double a1=0,a2=0;
for(i=0;i<n;i++){
dg[i]=g2[i]-g1[i];
dx[i]=xk1[i]-xk[i];
}
for(i=0;i<n;i++){
if(ia==1) { aa=a[1]; *ft=f[1]; break; }
d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]);
if(fabs(d)==0) break;
c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d;
ib=0;
if(ia==1) { xx=xk; break; }
ib=0;
for(i=0;i<n;i++) s[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
s[i]+= -h[i][j]*g1[j];
aa=lagrange(xk,ft,s);
xk1=iterate(xk,aa,s);
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}else{
if(*ft>f[1]) {a[0]=aa;f[0]=*ft;}
else if(*ft<f[1]) {a[2]=a[1];a[1]=aa;f[2]=f[1];f[1]=*ft;}
else if(*ft==f[1])
printf("\n\tf*=%f",f);
getch();
}
程序运行结果为：
x[1]*=4.000000
x[2]*=2.000000
f*=-8.000000
x1=iterate(x,a,s);
f=fny(x1);
return f;
}
void finding(double a[3],double f[3],double *xk,double *s){
double t=tt;
int i;
double a1,f1;
a[0]=0;f[0]=func(xk,a[0],s);
}
void main ()
{
int k;
double *xx,f;
double xk[n]={1,1};
xx=bfgs(xk);
f=fny(xx);
printf("\n\nThe Optimal Design Result Is:\n");
for(k=0;k<n;k++)
{printf("\n\tx[%d]*=%f",k+1,xx[k]);}
double *g,f1,f2,q;
int i;
g=(double*)malloc(n*sizeof(double));
f1=fny(xk);
for(i=0;i<n;i++){
q=ff; xk[i]=xk[i]+q; f2=fny(xk); g[i]=(f2-f1)/q; xk[i]=xk[i]-q;}
for(i=0;;i++){
a[1]=a[0]+t;
f[1]=func(xk,a[1],s);
if(f[1]<f[0]) break;
if(fabs(f[1]-f[0])>=ad) {
t=-t;
a[0]=a[1];f[0]=f[1];
}else{
if(ia==1) return; //break
t=t/2;ia=1;
}
}
for(i=0;;i++){
a[2]=a[1]+t;
f[2]=func(xk,a[2],s);
if(f[2]>f[1]) break;
t=2*t;
a[0]=a[1];f[0]=f[1];
a[1]=a[2];f[1]=f[2];
}
if(a[0]>a[2]){
a1=a[0];
f1=f[0];
a[0]=a[2];
return g;
}
double * bfgs(double *xk){
double u[n],v[n],h[n][n],dx[n],dg[n],s[n];
double aa,ib;
double *ft,*xk1,*g1,*g2,*xx,*x0=xk;
double fi;
int i,j,k;
ft=(double *)malloc(sizeof(double));
f[0]=f[2];
a[2]=a1;
f[2]=f1;
}
return;
}
double lagrange(double *xk,double *ft,double *s){
int i;
double a[3],f[3];
double b,c,d,aa;
finding(a,f,xk,s);
for(i=0;;i++){
int j;
u[i]=0;v[i]=0;
for(j=0;j<n;j++){
u[i]=u[i]+dg[j]*h[j][i];
v[i]=v[i]+dg[j]*h[i][j];
}
}
for(j=0;j<n;j++){
a1+=dx[j]*dg[j];
a2+=v[j]*dg[j];
}
if(fabs(a1)!=0){
if(fabs(c)==0) break;
b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]);
aa=-b/(2*c);
*ft=func(xk,aa,s);
if(fabs(aa-a[1])<=ad) {if(*ft>f[1]) aa=a[1];break;}
double fny(double *x){
double x1=x[0],x2=x[1];
double f;
f=x1*x1+2*x2*x2-4*x1-2*x1*x2;
return f;
}
double * iterate(double *x,double a,double *s){
double *x1;
xk1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++){
s[i]=0;
for(j=0;j<n;j++){
h[i][j]=0;
if(j==i) h[i][j]=1;
}
}
g1=gradient(xk);
fi=fny(xk);
x0=xk;
for(k=0;k<n;k++){
int i;
x1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
x1[i]=x[i]+a*s[i];
return x1;
}
double func(double *x,double a,double *s){
double *x1;
double f;
a2=1+a2/a1;
for(i=0;i<n;i++)
for(j=0;j<;j++)
h[i][j]+=(a2*dx[i]*dx[j]-v[i]*dx[j]-dx[i]*u[j])/a1;
}
xk=xk1; g1=g2;
}
}
if(*ft>fi) { *ft=fi; xx=xk;}
xk=x0;
return xx;
{a[0]=aa;a[2]=a[1]; f[0]=*ft;f[2]=f[1]; a[1]=(a[0]+a[2])/2; f[1]=func(xk,a[1],s); }
}
}
if(*ft>f[1]) {*ft=f[1];aa=a[1];}
return aa;
}
double *gradient(double *xk){
g2=gradient(xk1);
for(i=0;i<n;i++)
if((fabs(g2[i])>=ac)&&(fabs(g2[i]-g1[i])>=ac))
{ib=ib+1;}
if(ib==0) { xx=xk1; break; }
fi=*ft;
if(k==n-1){
int j;
xk=xk1;
for(i=0;i<n;i++)
if(aa>a[1]){
if(*ft>f[1]) {a[2]=aa;f[2]=*ft;}
else if(*ft<f[1]) {a[0]=a[1];a[1]=aa;f[0]=f[1];f[1]=*ft;}
else if(*ft==f[1]){
a[2]=aa;a[0]=a[1];