一些简单最优化方法的matlab实现

happy · 发表于 2006-11-17 10:12

%DFP计算算法(精确一维搜索,k=16,mul_count=9048,sum_count=4251花费时间很少)
syms x1 x2 ad;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;H0=[1,0;0,1];mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson)
mul_count=mul_count+2;sum_count=sum_count+1;
if k==0
p=-H0*g1;mul_count=mul_count+7;sum_count=sum_count+4;
else
H1=H0-H0*yk*yk'*H0/(yk'*H0*yk)+sk*sk'/(yk'*sk);mul_count=mul_count+36;sum_count=sum_count+18;
p=-H1*g1;mul_count=mul_count+7;sum_count=sum_count+4;
H0=H1;
end
x00=x0;
y=subs(f,{x1,x2},{x0(1,1)+ad*p(1,1),x0(2,1)+ad*p(2,1)});mul_count=mul_count+18;sum_count=sum_count+12;
result1=interzone(y,ad);
b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
result2=Goldsearch(y,ad,b1);
arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
yk=g1-g0; sk=x0-x00;
k=k+1;
mul_count=mul_count+9;sum_count=sum_count+11;
end;
k
x0
mul_count
sum_count

复制代码

happy · 发表于 2006-11-17 10:13

%DFP计算算法(精确一维搜索,k=16,mul_count=9048,sum_count=4251花费时间很少)
syms x1 x2 ar;
f=x1^2-x1*x2+x2^2+2*x1-4*x2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=[2,2]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;H0=[1,0;0,1];mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson&k<=2)
mul_count=mul_count+2;sum_count=sum_count+1;
if k==0
p=-H0*g1;mul_count=mul_count+7;sum_count=sum_count+4;
else
H1=H0-H0*yk*yk'*H0/(yk'*H0*yk)+sk*sk'/(yk'*sk);mul_count=mul_count+36;sum_count=sum_count+18;
p=-H1*g1;mul_count=mul_count+7;sum_count=sum_count+4;
H0=H1;
end
x00=x0;
y=subs(f,{x1,x2},{x0(1,1)+ar*p(1,1),x0(2,1)+ar*p(2,1)});mul_count=mul_count+18;sum_count=sum_count+12;
result1=qujian(y,ar);
b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
result2=Asearch(y,ar,b1);
arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
yk=g1-g0; sk=x0-x00;
k=k+1;
mul_count=mul_count+9;sum_count=sum_count+11;
end;
k
x0
mul_count
sum_count

复制代码

happy · 发表于 2006-11-17 10:13

%FR直接计算算法(不精确一维搜索,k=163,mul_count=32380,sum_count=20502时间略长)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson&k<=1000)
mul_count=mul_count+2;sum_count=sum_count+1;
if k==0
p=-g1;
else
bta=g1'*g1/(g0'*g0);
p=-g1+bta*p0;
mul_count=mul_count+7;sum_count=sum_count+4;
end
result=Usearch1(f,x1,x2,df,x0,p);
arf=result(1);Mul=result(2);Sum=result(3);
mul_count=mul_count+Mul;sum_count=sum_count+Sum;
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
k=k+1;
mul_count=mul_count+9;sum_count=sum_count+7;
end;
k
x0
mul_count
sum_count

复制代码

happy · 发表于 2006-11-17 10:13

%FR重新开始计算算法(不精确一维搜索,k=496,mul_count=76209,sum_count=47950花费时间最长)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson&k<1000)
mul_count=mul_count+2;sum_count=sum_count+1;
if k==0
p=-g1;
else
bta=g1'*g1/(g0'*g0); mul_count=mul_count+5;sum_count=sum_count+2;
if(mod(k,2))
p=-g1+bta*p0; mul_count=mul_count+3;sum_count=sum_count+2;
else
p=-g1;
end
end
result=Usearch1(f,x1,x2,df,x0,p);
arf=result(1);Mul=result(2);Sum=result(3);
mul_count=mul_count+Mul;sum_count=sum_count+Sum;
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
k=k+1;
mul_count=mul_count+9;sum_count=sum_count+7;
end;
k
x0
mul_count
sum_count

复制代码

happy · 发表于 2006-11-17 10:13

%FR重新开始计算算法(精确一维搜索,k=39,mul_count=17815,sum=7965,花费时间较少)
syms x1 x2 ar;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson)
mul_count=mul_count+2;sum_count=sum_count+1;
if k==0
p=-g1;
else
bta=g1'*g1/(g0'*g0); mul_count=mul_count+5;sum_count=sum_count+2;
if(mod(k,2))
p=-g1+bta*p0; mul_count=mul_count+3;sum_count=sum_count+2;
else
p=-g1;
end
end
y=subs(f,{x1,x2},{x0(1,1)+ar*p(1,1),x0(2,1)+ar*p(2,1)}); mul_count=mul_count+18;sum_count=sum_count+12;
result1=qujian(y,ar);
b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
result2=Asearch(y,ar,b1);
arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
k=k+1;
mul_count=mul_count+9;sum_count=sum_count+7;
end;
k
x0
mul_count
sum_count

复制代码

happy · 发表于 2006-11-17 10:13

%FR直接计算算法(精确一维搜索,k=71,mul_count=31360,sum_count=13875,时间略长)
syms x1 x2 ar;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson)
mul_count=mul_count+2;sum_count=sum_count+1;
if k==0
p=-g1;
else
bta=g1'*g1/(g0'*g0);
p=-g1+bta*p0;
mul_count=mul_count+7;sum_count=sum_count+4;
end
y=subs(f,{x1,x2},{x0(1,1)+ar*p(1,1),x0(2,1)+ar*p(2,1)});mul_count=mul_count+18;sum_count=sum_count+12;
result1=qujian(y,ar);
b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
result2=Asearch(y,ar,b1);
arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
k=k+1;
mul_count=mul_count+9;sum_count=sum_count+7;
end;
k
x0
mul_count
sum_count

复制代码

happy · 发表于 2006-11-17 10:29

%求min(x1^2+x2^2) 条件：1-x1-x2<=0
%0.618法，f=x1^2+x2^2+ck*(max(1-x1-x2),0)^2
c_jingdu=[0.001,0.001];
a=[0,0];
b=[20,20]; %单谷搜索区间及其精度
t1=[0.01,0.01];
x11=t1(1,1);
x12=t1(1,2);
t2=[9.90,9.90]; %最初的两个探索点 %给出最初两个探索点
x21=t2(1,1);
x22=t2(1,2);
ck=20; %给出罚参数
f1=x11^2+x12^2+ck*(max(1-x11-x12,0))^2;
f2=x21^2+x22^2+ck*(max(1-x21-x22,0))^2; %给出f的初值，进而迭代
t1=a+0.382*(b-a);
t2=a+0.618*(b-a);
N=10000; %给出迭代次数
for i=1:N
if f1<=f2
if x21-a(1,1)<=c_jingdu(1,1)
break
else
b=t2;
t2=t1;
t1=a+0.382*(b-a);
f2=f1;
f1=x11^2+x22^2+ck*(max(1-x11-x12,0))^2;
end
else
if b(1,1)-x11<=c_jingdu(1,1)
break
else
a=t1;
t1=t2;
t2=a+0.618*(b-a);
f1=f2;
f2=x21^2+x22^2+ck*(max(1-x21-x22,0))^2;
end
end
ck=ck*2;
end
(t1+t2)/2
t1
t2

复制代码

happy · 发表于 2006-11-17 10:29

% mg0523076 潘向昱
%求min(x1^2+x2^2),约束条件1-x1-x2<=0
%0.618法
a=[0,0];
b=[10,10]; %确定单谷搜索[a,b]
c_wucha=[1.0e-7,1.0e-7]; %给定最后精度c_wucha
t1=[0.1,0.1];
t2=[9.9,9.9];
ck=50; %给出罚参数
f1=t1(1,1)^2+t1(1,2)^2+ck*(max(1-t1(1,1)-t1(1,2),0))^2;
f2=t2(1,1)^2+t2(1,2)^2+ck*(max(1-t2(1,1)-t2(1,2),0))^2; %计算最初的两个探索点
t1=a+0.382*(b-a);
t2=a+0.618*(b-a);
while(1)
if f1<=f2 %如果f1<f2
if t2(1,1)-a(1,1)<=c_wucha(1,1)
break
else
b=t2;
t2=t1;
t1=a+0.382*(b-a);
f2=f1;
f1=t1(1,1)^2+t2(1,2)^2+ck*(max(1-t1(1,1)-t1(1,2),0))^2;
end
else
if b(1,1)-t1(1,1)<=c_wucha(1,1)
break
else
a=t1;
t1=t2;
t2=a+0.618*(b-a);
f1=f2;
f2=t2(1,1)^2+t2(1,2)^2+ck*(max(1-t2(1,1)-t2(1,2),0))^2;
end
end
ck=ck*2; %调整罚参数
end
(t1+t2)/2
t1
t2

复制代码

liweidlut · 发表于 2006-11-22 20:24

好贴哪能不顶
谢了啊

coldspring · 发表于 2006-11-23 19:26

15楼 result1=interzone(y,ad);
result2=Goldsearch(y,ad,b1);[/color]在那？
楼主能不能发一下

[ 本帖最后由 coldspring 于 2006-11-23 19:37 编辑 ]

coldspring · 发表于 2006-11-23 20:02

%DFP计算算法(精确一维搜索,k=16,mul_count=9048,sum_count=4251花费时间很少)
syms x1 x2 ar;
这里面的 ar 怎么解释！

suffer · 发表于 2006-11-25 10:46

原帖由 coldspring 于 2006-11-23 19:26 发表
15楼 result1=interzone(y,ad);
result2=Goldsearch(y,ad,b1);color]在那？
楼主能不能发一下

找了一下没有找到，Goldsearch估计是黄金分割搜索

suffer · 发表于 2006-11-25 10:48

原帖由 coldspring 于 2006-11-23 20:02 发表
%DFP计算算法(精确一维搜索,k=16,mul_count=9048,sum_count=4251花费时间很少)
syms x1 x2 ar;
这里面的 ar 怎么解释！

你是指在算法中的含义？

renxiaoyao007 · 发表于 2007-4-20 21:59

FR算法里面的这个函数怎么没有定义啊？
result1=qujian(y,ar);
qujian？？？？？？

qq8675332 · 发表于 2007-6-1 03:19

有没有K均值的代码呢

账号		自动登录	找回密码
密码			我要加入

[共享资源] 一些简单最优化方法的matlab实现