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

happy

1. clear;
   
2. A=[2,1,4,0;2,2,0,4];
   
3. B=[2,3];
   
4. C=[12,8,16,12];
   
5. f=y(A,B,C)

复制代码

1. %对偶单纯形法
   
2. function result=y(A,B,C)
   
3. a=-1*A;
   
4. [row line]=size(A);
   
5. a=[a eye(row,row)];
   
6. b=-1*B;
   
7. c=[C zeros(1,row)];
   
8. dex=c;                                                %dex为判别数                          
   
9. basic=line+1:1;line+row;                              %basic为基可行解
   
10. X=zeros(1,line);
    
11. X1=zeros(1,line+row);
    
12. b1=sort(b);
    
13. while(b1(1)<0)
    
14.     m=0;   
    
15.     for i=1:length(b)
    
16.         if(b(i)<m)
    
17.             m=b(i);   j=i;                            %求出最小的b
    
18.         end
    
19.     end;   
    
20.    n=a(j,:);
    
21.    count=0;
    
22.    for i=1:length(n)
    
23.        if(n(i)>=0)   
    
24.            count=count+1;
    
25.        end
    
26.    end;
    
27.    if(count==length(n))   
    
28.        error('无最优解');
    
29.    end
    
30.    seta=dex./n;
    
31.    seta1=-Inf;
    
32.    for i=1:length(seta)
    
33.       if(seta(i)>=seta1&n(i)<0)                    %求出最大且为负的seta
    
34.              seta1=seta(i);   k=i;
    
35.       end
    
36.    end;   
    
37.    basic(j)=k;
    
38.    b(j)=b(j)/a(j,k);
    
39.    a(j,:)=a(j,:)/a(j,k);
    
40.    for i=1:length(b)
    
41.        if(i~=j)
    
42.            b(i)=b(i)-b(j)*a(i,k);
    
43.            a(i,:)=a(i,:)-a(j,:)*a(i,k);
    
44.        end   
    
45.    end;
    
46.   q=dex(k);
    
47.   dex=dex-q*a(j,:);
    
48.   b1=sort(b);
    
49. end;
    
50. for i=1:length(basic)
    
51.     temp=basic(i);
    
52.     X1(temp)=b(i);
    
53. end;
    
54. for i=1:length(X)
    
55.     X(i)=X1(i);
    
56. end;
    
57. X=X'
    
58. X1=X1';
    
59. result=c*X1;

复制代码

happy

1. %不精确一维搜索法
   
2. function result=Usearch2(f,x1,x2,df,x0,p)
   
3. mu=0.4;sgma=0.8;a=0;b=inf;arf=1;m_count=0;s_count=0;
   
4. pk=p;
   
5. x3=x0;
   
6. x4=x3+arf*pk;m_count=m_count+2;s_count=s_count+2;
   
7. f1=subs(f,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+4;s_count=s_count+3;
   
8. f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+4;s_count=s_count+3;
   
9. gk1=subs(df,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+7;s_count=s_count+4;
   
10. gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+7;s_count=s_count+4;
    
11. while (f1-f2<=-mu*arf*gk1'*pk)
    
12.    b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
    
13.    m_count=m_count+19;s_count=s_count+12;
    
14. end;
    
15. m_count=m_count+5;s_count=s_count+2;
    
16. while (1>0)
    
17.     if(gk2'*pk<sgma*gk1'*pk)
    
18.         a=arf;a=min(2*arf,(a+b)/2);x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
    
19.         m_count=m_count+21;s_count=s_count+12;
    
20.     while (f1-f2<=-mu*arf*gk1'*pk)
    
21.         b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
    
22.         m_count=m_count+19;s_count=s_count+12;
    
23.     end;
    
24.     m_count=m_count+5;s_count=s_count+2;
    
25.     else
    
26.         break;
    
27.     end
    
28. end;
    
29. result=[arf,m_count,s_count];

复制代码

happy

1. %不精确一维搜索法
   
2. function result=Usearch1(f,x1,x2,df,x0,p)
   
3. mu=0.001;sgma=0.99;a=0;b=inf;arf=1;m_count=0;s_count=0;
   
4. pk=p;
   
5. x3=x0;
   
6. x4=x3+arf*pk;m_count=m_count+2;s_count=s_count+2;
   
7. f1=subs(f,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+4;s_count=s_count+3;
   
8. f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+4;s_count=s_count+3;
   
9. gk1=subs(df,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+7;s_count=s_count+4;
   
10. gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+7;s_count=s_count+4;
    
11. while (f1-f2<=-mu*arf*gk1'*pk)
    
12.    b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
    
13.    m_count=m_count+19;s_count=s_count+12;
    
14. end;
    
15. m_count=m_count+5;s_count=s_count+2;
    
16. while (1>0)
    
17.     if(gk2'*pk<sgma*gk1'*pk)
    
18.         a=arf;a=min(2*arf,(a+b)/2);x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
    
19.         m_count=m_count+21;s_count=s_count+12;
    
20.     while (f1-f2<=-mu*arf*gk1'*pk)
    
21.         b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
    
22.         m_count=m_count+19;s_count=s_count+12;
    
23.     end;
    
24.     m_count=m_count+5;s_count=s_count+2;
    
25.     else
    
26.         break;
    
27.     end
    
28. end;
    
29. result=[arf,m_count,s_count];

复制代码

happy

1. %精确一维搜索
   
2. function result=Asearch(y,ar,b1)
   
3. a=0;b=b1;e=1e-5;a1=a+0.382*(b-a);f1=subs(y,{ar},{a1});a2=a+0.618*(b-a);f2=subs(y,{ar},{a2});M_count=0;S_count=0;
   
4. M_count=M_count+30;S_count=S_count+12;
   
5. while abs(b-a)>=e
   
6.           if f1<f2
   
7.              b=a2;a2=a1;f2=f1;a1=a+0.382*(b-a);f1=subs(y,{ar},{a1});M_count=M_count+15;S_count=S_count+6;
   
8.           else
   
9.              a=a1;a1=a2;f1=f2;a2=a+0.618*(b-a);f2=subs(y,{ar},{a2});M_count=M_count+15;S_count=S_count+6;
   
10.           end
    
11.           S_count=S_count+1;
    
12. end;
    
13. answer=0.5*(a+b);
    
14. M_count=M_count+1;S_count=S_count+1;
    
15. result=[answer,M_count,S_count];

复制代码

happy

1. %BFGS计算算法(不精确一维搜索,k=21,mul_count=2380,sum_count=1448,花费时间很少)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson&k<300)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-H0*g1;mul_count=mul_count+4;sum_count=sum_count+2;
    
13.     else
    
14.         w1=sqrt(yk'*H0*yk)*((sk/(yk'*sk)-H0*yk/(yk'*H0*yk)));mul_count=mul_count+22;sum_count=sum_count+10;
    
15.         H1=H0-H0*yk*yk'*H0/(yk'*H0*yk)+sk*sk'/(yk'*sk)+w1*w1';mul_count=mul_count+40;sum_count=sum_count+22;
    
16.         p=-H1*g1;mul_count=mul_count+4;sum_count=sum_count+2;
    
17.         H0=H1;
    
18.     end
    
19.     x00=x0;
    
20.     result=Usearch1(f,x1,x2,df,x0,p);
    
21.     arf=result(1);Mul=result(2);Sum=result(3);
    
22.     mul_count=mul_count+Mul;sum_count=sum_count+Sum;
    
23.     x0=x0+arf*p;
    
24.     g0=g1;
    
25.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
26.     p0=p;
    
27.     yk=g1-g0; sk=x0-x00;
    
28.     k=k+1;
    
29.     mul_count=mul_count+9;sum_count=sum_count+11;
    
30. end;
    
31. k
    
32. x0
    
33. mul_count
    
34. sum_count

复制代码

happy

1. %BFGS计算算法(精确一维搜索,k=16,mul_count=8888,sum_count=4227,花费时间很少)
   
2. syms x1 x2 ar;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-H0*g1;mul_count=mul_count+4;sum_count=sum_count+2;
    
13.     else
    
14.         w1=sqrt(yk'*H0*yk)*((sk/(yk'*sk)-H0*yk/(yk'*H0*yk)));mul_count=mul_count+22;sum_count=sum_count+10;
    
15.         H1=H0-H0*yk*yk'*H0/(yk'*H0*yk)+sk*sk'/(yk'*sk)+w1*w1';mul_count=mul_count+40;sum_count=sum_count+22;
    
16.         p=-H1*g1;mul_count=mul_count+4;sum_count=sum_count+2;
    
17.         H0=H1;
    
18.     end
    
19.     x00=x0;
    
20.     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;
    
21.     result1=qujian(y,ar);
    
22.     b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
    
23.     result2=Asearch(y,ar,b1);
    
24.     arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
    
25.     x0=x0+arf*p;
    
26.     g0=g1;
    
27.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
28.     p0=p;
    
29.     yk=g1-g0; sk=x0-x00;
    
30.     k=k+1;
    
31.     mul_count=mul_count+9;sum_count=sum_count+11;
    
32. end;
    
33. k
    
34. x0
    
35. mul_count
    
36. sum_count

复制代码

happy

1. %最速下降法(不精确一维搜索,k＝1000,花费时间太多,精度远不能达到要求)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson&k<1000)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     p=-g1;
    
12.     result=Usearch1(f,x1,x2,df,x0,p);
    
13.     arf=result(1);Mul=result(2);Sum=result(3);
    
14.     mul_count=mul_count+Mul;sum_count=sum_count+Sum;
    
15.     x0=x0+arf*p;
    
16.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
17.     k=k+1;
    
18.     mul_count=mul_count+9;sum_count=sum_count+7;
    
19. end;
    
20. k
    
21. x0

复制代码

happy

1. %最速下降法(精确一维搜索,k=1000,花费时间太多,精度远不能达到要求)
   
2. syms x1 x2 ar;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson&k<1000)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     p=-g1;
    
12.     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;
    
13.     result1=qujian(y,ar);
    
14.     b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
    
15.     result2=Asearch(y,ar,b1);
    
16.     arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
    
17.     x0=x0+arf*p;
    
18.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
19.     k=k+1;
    
20.     mul_count=mul_count+9;sum_count=sum_count+7;
    
21. end;
    
22. k

复制代码

happy

1. %阻尼Newton算法(k=2,mul_count=44,sum_count=28,最少迭代)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. G=jacobian(df,v);
   
8. epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;
   
9. mul_count=mul_count+12;sum_count=sum_count+6;
   
10. while(norm(g1)>epson)
    
11.     p=-G1\g1;
    
12.     x0=x0+p;
    
13.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
14.     G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
    
15.     k=k+1;
    
16.     mul_count=mul_count+16;sum_count=sum_count+11;
    
17. end;
    
18. k
    
19. x0
    
20. mul_count
    
21. sum_count

复制代码

happy

1. %进退法
   
2. function result=qujian(y,ar)
   
3. t0=0; dex=1;M_count=0;S_count=0;
   
4. t1=t0+dex;S_count=S_count+1;
   
5. ft0=subs(y,{ar},{t0});M_count=M_count+14;S_count=S_count+4;
   
6. ft1=subs(y,{ar},{t1});M_count=M_count+14;S_count=S_count+4;
   
7. if(ft1<=ft0)
   
8.     dex=2*dex;t2=t1+dex;ft2=subs(y,{ar},{t2});M_count=M_count+15;S_count=S_count+5;
   
9.     while(ft1>ft2)
   
10.         t1=t2;dex=2*dex;t2=t1+dex;ft1=subs(y,{ar},{t1});ft2=subs(y,{ar},{t2});M_count=M_count+29;S_count=S_count+9;
    
11.     end
    
12. else
    
13.         dex=dex/2;t2=t1;t1=t2-dex;ft1=subs(y,{ar},{t1});M_count=M_count+15;S_count=S_count+5;
    
14.     while(ft1>ft0)
    
15.         dex=dex/2;t2=t1;t1=t2-dex;ft1=subs(y,{ar},{t1});M_count=M_count+15;S_count=S_count+5;
    
16.     end
    
17. end
    
18. result=[t2,M_count,S_count];

复制代码

happy

1. %PRP直接计算算法(不精确一维搜索,k=83,mul_count=19842,sum_count=12722时间略长)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson&k<1000)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-g1;
    
13.     else
    
14.         bta=g1'*(g1-g0)/(g0'*g0);
    
15.         p=-g1+bta*p0;
    
16.         mul_count=mul_count+7;sum_count=sum_count+6;
    
17.     end
    
18.     result=Usearch2(f,x1,x2,df,x0,p);
    
19.     arf=result(1);Mul=result(2);Sum=result(3);
    
20.     mul_count=mul_count+Mul;sum_count=sum_count+Sum;
    
21.     x0=x0+arf*p;
    
22.     g0=g1;
    
23.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
24.     p0=p;
    
25.     k=k+1;
    
26.     mul_count=mul_count+9;sum_count=sum_count+7;
    
27. end;
    
28. k
    
29. x0
    
30. mul_count
    
31. sum_count

复制代码

happy

1. %PRP重新开始计算算法(不精确一维搜索,k=133,mul_count=18421,sum_count=11850时间略长)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson&k<1000)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-g1;
    
13.     else
    
14.         bta=g1'*(g1-g0)/(g0'*g0); mul_count=mul_count+5;sum_count=sum_count+4;
    
15.         if(mod(k,2))
    
16.                p=-g1+bta*p0; mul_count=mul_count+3;sum_count=sum_count+2;
    
17.         else
    
18.                p=-g1;
    
19.         end
    
20.     end
    
21.     result=Usearch2(f,x1,x2,df,x0,p);
    
22.     arf=result(1);Mul=result(2);Sum=result(3);
    
23.     mul_count=mul_count+Mul;sum_count=sum_count+Sum;
    
24.     x0=x0+arf*p;
    
25.     g0=g1;
    
26.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
27.     p0=p;
    
28.     k=k+1;
    
29.     mul_count=mul_count+9;sum_count=sum_count+7;
    
30. end;
    
31. k
    
32. x0
    
33. mul_count
    
34. sum_count

复制代码

happy

1. %PRP直接计算降法(精确一维搜索,k=19,mul_count=8390,sum_count=3745,花费时间很少)
   
2. syms x1 x2 ar;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-g1;
    
13.     else
    
14.         bta=(g1'*(g1-g0))/(g0'*g0);
    
15.         p=-g1+bta*p0;
    
16.         mul_count=mul_count+7;sum_count=sum_count+6;
    
17.     end
    
18.     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;
    
19.     result1=qujian(y,ar);
    
20.     b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
    
21.     result2=Asearch(y,ar,b1);
    
22.     arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
    
23.     x0=x0+arf*p;
    
24.     g0=g1;
    
25.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
26.     p0=p;
    
27.     k=k+1;
    
28.     mul_count=mul_count+9;sum_count=sum_count+7;
    
29. end;
    
30. k
    
31. x0
    
32. mul_count
    
33. sum_count

复制代码

happy

1. %PRP重新开始计算算法(精确一维搜索,k=21,mul_count=9901,sum_count=4470,花费时间很少)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-g1;
    
13.     else
    
14.         bta=g1'*(g1-g0)/(g0'*g0); mul_count=mul_count+5;sum_count=sum_count+4;
    
15.         if(mod(k,2))
    
16.                p=-g1+bta*p0; mul_count=mul_count+3;sum_count=sum_count+2;
    
17.         else
    
18.                p=-g1;
    
19.         end
    
20.     end
    
21.     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;
    
22.     result1=qujian(y,ar);
    
23.     b1=result1(1);mul_count=mul_count+result1(2);sum_count=sum_count+result1(3);
    
24.     result2=Asearch(y,ar,b1);
    
25.     arf=result2(1);mul_count=mul_count+result2(2);sum_count=sum_count+result2(3);
    
26.     x0=x0+arf*p;
    
27.     g0=g1;
    
28.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
29.     p0=p;
    
30.     k=k+1;
    
31.     mul_count=mul_count+9;sum_count=sum_count+7;
    
32. end;
    
33. k
    
34. x0
    
35. mul_count
    
36. sum_count

复制代码

happy

1. %DFP计算算法(不精确一维搜索,k=27,mul_count=2326,sum_count=1462,花费时间很少)
   
2. syms x1 x2;
   
3. f=50*(x2-x1^2)^2+(1-x1)^2;
   
4. v=[x1,x2];
   
5. df=jacobian(f,v);
   
6. df=df.';
   
7. 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;
   
8. mul_count=mul_count+7;sum_count=sum_count+4;
   
9. while(norm(g1)>epson&k<1000)
   
10.     mul_count=mul_count+2;sum_count=sum_count+1;
    
11.     if k==0
    
12.          p=-H0*g1;mul_count=mul_count+7;sum_count=sum_count+4;
    
13.     else
    
14.         H1=H0-H0*yk*yk'*H0/(yk'*H0*yk)+sk*sk'/(yk'*sk);mul_count=mul_count+36;sum_count=sum_count+18;
    
15.         p=-H1*g1;mul_count=mul_count+7;sum_count=sum_count+4;
    
16.         H0=H1;
    
17.     end
    
18.     x00=x0;
    
19.     result=Usearch1(f,x1,x2,df,x0,p);
    
20.     arf=result(1);Mul=result(2);Sum=result(3);
    
21.     mul_count=mul_count+Mul;sum_count=sum_count+Sum;
    
22.     x0=x0+arf*p;
    
23.     g0=g1;
    
24.     g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    
25.     p0=p;
    
26.     yk=g1-g0; sk=x0-x00;
    
27.     k=k+1;
    
28.     mul_count=mul_count+9;sum_count=sum_count+11;
    
29. end;
    
30. k
    
31. x0
    
32. mul_count
    
33. sum_count

复制代码