happy 发表于 2006-11-17 10:07

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

clear;
A=;
B=;
C=;
f=y(A,B,C)

%对偶单纯形法
function result=y(A,B,C)
a=-1*A;
=size(A);
a=;
b=-1*B;
c=;
dex=c;                                                %dex为判别数                        
basic=line+1:1;line+row;                              %basic为基可行解
X=zeros(1,line);
X1=zeros(1,line+row);
b1=sort(b);
while(b1(1)<0)
    m=0;   
    for i=1:length(b)
      if(b(i)<m)
            m=b(i);   j=i;                            %求出最小的b
      end
    end;   
   n=a(j,:);
   count=0;
   for i=1:length(n)
       if(n(i)>=0)   
         count=count+1;
       end
   end;
   if(count==length(n))   
       error('无最优解');
   end
   seta=dex./n;
   seta1=-Inf;
   for i=1:length(seta)
      if(seta(i)>=seta1&n(i)<0)                  %求出最大且为负的seta
             seta1=seta(i);   k=i;
      end
   end;   
   basic(j)=k;
   b(j)=b(j)/a(j,k);
   a(j,:)=a(j,:)/a(j,k);
   for i=1:length(b)
       if(i~=j)
         b(i)=b(i)-b(j)*a(i,k);
         a(i,:)=a(i,:)-a(j,:)*a(i,k);
       end   
   end;
q=dex(k);
dex=dex-q*a(j,:);
b1=sort(b);
end;
for i=1:length(basic)
    temp=basic(i);
    X1(temp)=b(i);
end;
for i=1:length(X)
    X(i)=X1(i);
end;
X=X'
X1=X1';
result=c*X1;

happy 发表于 2006-11-17 10:08

%不精确一维搜索法
function result=Usearch2(f,x1,x2,df,x0,p)
mu=0.4;sgma=0.8;a=0;b=inf;arf=1;m_count=0;s_count=0;
pk=p;
x3=x0;
x4=x3+arf*pk;m_count=m_count+2;s_count=s_count+2;
f1=subs(f,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+4;s_count=s_count+3;
f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+4;s_count=s_count+3;
gk1=subs(df,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+7;s_count=s_count+4;
gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+7;s_count=s_count+4;
while (f1-f2<=-mu*arf*gk1'*pk)
   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)});
   m_count=m_count+19;s_count=s_count+12;
end;
m_count=m_count+5;s_count=s_count+2;
while (1>0)
    if(gk2'*pk<sgma*gk1'*pk)
      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)});
      m_count=m_count+21;s_count=s_count+12;
    while (f1-f2<=-mu*arf*gk1'*pk)
      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)});
      m_count=m_count+19;s_count=s_count+12;
    end;
    m_count=m_count+5;s_count=s_count+2;
    else
      break;
    end
end;
result=;

happy 发表于 2006-11-17 10:08

%不精确一维搜索法
function result=Usearch1(f,x1,x2,df,x0,p)
mu=0.001;sgma=0.99;a=0;b=inf;arf=1;m_count=0;s_count=0;
pk=p;
x3=x0;
x4=x3+arf*pk;m_count=m_count+2;s_count=s_count+2;
f1=subs(f,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+4;s_count=s_count+3;
f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+4;s_count=s_count+3;
gk1=subs(df,{x1,x2},{x3(1,1),x3(2,1)});m_count=m_count+7;s_count=s_count+4;
gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});m_count=m_count+7;s_count=s_count+4;
while (f1-f2<=-mu*arf*gk1'*pk)
   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)});
   m_count=m_count+19;s_count=s_count+12;
end;
m_count=m_count+5;s_count=s_count+2;
while (1>0)
    if(gk2'*pk<sgma*gk1'*pk)
      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)});
      m_count=m_count+21;s_count=s_count+12;
    while (f1-f2<=-mu*arf*gk1'*pk)
      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)});
      m_count=m_count+19;s_count=s_count+12;
    end;
    m_count=m_count+5;s_count=s_count+2;
    else
      break;
    end
end;
result=;

happy 发表于 2006-11-17 10:09

%精确一维搜索
function result=Asearch(y,ar,b1)
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;
M_count=M_count+30;S_count=S_count+12;
while abs(b-a)>=e
          if f1<f2
             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;
          else
             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;
          end
          S_count=S_count+1;
end;
answer=0.5*(a+b);
M_count=M_count+1;S_count=S_count+1;
result=;

happy 发表于 2006-11-17 10:09

%BFGS计算算法(不精确一维搜索,k=21,mul_count=2380,sum_count=1448,花费时间很少)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;H0=;mul_count=0;sum_count=0;
mul_count=mul_count+7;sum_count=sum_count+4;
while(norm(g1)>epson&k<300)
    mul_count=mul_count+2;sum_count=sum_count+1;
    if k==0
         p=-H0*g1;mul_count=mul_count+4;sum_count=sum_count+2;
    else
      w1=sqrt(yk'*H0*yk)*((sk/(yk'*sk)-H0*yk/(yk'*H0*yk)));mul_count=mul_count+22;sum_count=sum_count+10;
      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;
      p=-H1*g1;mul_count=mul_count+4;sum_count=sum_count+2;
      H0=H1;
    end
    x00=x0;
    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;
    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:10

%BFGS计算算法(精确一维搜索,k=16,mul_count=8888,sum_count=4227,花费时间很少)
syms x1 x2 ar;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;H0=;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+4;sum_count=sum_count+2;
    else
      w1=sqrt(yk'*H0*yk)*((sk/(yk'*sk)-H0*yk/(yk'*H0*yk)));mul_count=mul_count+22;sum_count=sum_count+10;
      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;
      p=-H1*g1;mul_count=mul_count+4;sum_count=sum_count+2;
      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:10

%最速下降法(不精确一维搜索,k＝1000,花费时间太多,精度远不能达到要求)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';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;
    p=-g1;
    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;
    g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    k=k+1;
    mul_count=mul_count+9;sum_count=sum_count+7;
end;
k
x0

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

%最速下降法(精确一维搜索,k=1000,花费时间太多,精度远不能达到要求)
syms x1 x2 ar;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';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;
    p=-g1;
    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;
    g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    k=k+1;
    mul_count=mul_count+9;sum_count=sum_count+7;
end;
k

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

%阻尼Newton算法(k=2,mul_count=44,sum_count=28,最少迭代)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
G=jacobian(df,v);
epson=1e-12;x0=';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;
mul_count=mul_count+12;sum_count=sum_count+6;
while(norm(g1)>epson)
    p=-G1\g1;
    x0=x0+p;
    g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
    G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
    k=k+1;
    mul_count=mul_count+16;sum_count=sum_count+11;
end;
k
x0
mul_count
sum_count

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

%进退法
function result=qujian(y,ar)
t0=0; dex=1;M_count=0;S_count=0;
t1=t0+dex;S_count=S_count+1;
ft0=subs(y,{ar},{t0});M_count=M_count+14;S_count=S_count+4;
ft1=subs(y,{ar},{t1});M_count=M_count+14;S_count=S_count+4;
if(ft1<=ft0)
    dex=2*dex;t2=t1+dex;ft2=subs(y,{ar},{t2});M_count=M_count+15;S_count=S_count+5;
    while(ft1>ft2)
      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;
    end
else
      dex=dex/2;t2=t1;t1=t2-dex;ft1=subs(y,{ar},{t1});M_count=M_count+15;S_count=S_count+5;
    while(ft1>ft0)
      dex=dex/2;t2=t1;t1=t2-dex;ft1=subs(y,{ar},{t1});M_count=M_count+15;S_count=S_count+5;
    end
end
result=;

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

%PRP直接计算算法(不精确一维搜索,k=83,mul_count=19842,sum_count=12722时间略长)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';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'*g0);
      p=-g1+bta*p0;
      mul_count=mul_count+7;sum_count=sum_count+6;
    end
    result=Usearch2(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:11

%PRP重新开始计算算法(不精确一维搜索,k=133,mul_count=18421,sum_count=11850时间略长)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';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'*g0); mul_count=mul_count+5;sum_count=sum_count+4;
      if(mod(k,2))
               p=-g1+bta*p0; mul_count=mul_count+3;sum_count=sum_count+2;
      else
               p=-g1;
      end
    end
    result=Usearch2(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:11

%PRP直接计算降法(精确一维搜索,k=19,mul_count=8390,sum_count=3745,花费时间很少)
syms x1 x2 ar;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';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'*g0);
      p=-g1+bta*p0;
      mul_count=mul_count+7;sum_count=sum_count+6;
    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:12

%PRP重新开始计算算法(精确一维搜索,k=21,mul_count=9901,sum_count=4470,花费时间很少)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';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'*g0); mul_count=mul_count+5;sum_count=sum_count+4;
      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:12

%DFP计算算法(不精确一维搜索,k=27,mul_count=2326,sum_count=1462,花费时间很少)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=;
df=jacobian(f,v);
df=df.';
epson=1e-12;x0=';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;H0=;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=-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;
    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;
    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
页: [1] 2 3
查看完整版本: 一些简单最优化方法的matlab实现