someoptimiummethod.txt

一些简单最优化方法的matlab实现和非线性规划算法

振动论坛 ? Matlab讨论区 ? 一些简单最优化方法的matlab实现
查看完整版本: 一些简单最优化方法的matlab实现
happy 2006-11-17 10:07一些简单最优化方法的matlab实现
[code]clear;
A=[2,1,4,0;2,2,0,4];
B=[2,3];
C=[12,8,16,12];
f=y(A,B,C)[/code]

[code]%对偶单纯形法
function result=y(A,B,C)
a=-1*A;
[row line]=size(A);
a=[a eye(row,row)];
b=-1*B;
c=[C zeros(1,row)];
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;[/code]
happy 2006-11-17 10:08[code]%不精确一维搜索法
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=[arf,m_count,s_count];[/code]
happy 2006-11-17 10:08[code]%不精确一维搜索法
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=[arf,m_count,s_count];[/code]
happy 2006-11-17 10:09[code]%精确一维搜索
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=[answer,M_count,S_count];[/code]
happy 2006-11-17 10:09[code]%BFGS计算算法(不精确一维搜索,k=21,mul_count=2380,sum_count=1448,花费时间很少)
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;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<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[/code]
happy 2006-11-17 10:10[code]%BFGS计算算法(精确一维搜索,k=16,mul_count=8888,sum_count=4227,花费时间很少)
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;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+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[/code]
happy 2006-11-17 10:10[code]%最速下降法(不精确一维搜索,k＝1000,花费时间太多,精度远不能达到要求)
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;
    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[/code]
happy 2006-11-17 10:10[code]%最速下降法(精确一维搜索,k=1000,花费时间太多,精度远不能达到要求)
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&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[/code]
happy 2006-11-17 10:10[code]%阻尼Newton算法(k=2,mul_count=44,sum_count=28,最少迭代)
syms x1 x2;
f=50*(x2-x1^2)^2+(1-x1)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
G=jacobian(df,v);
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;
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[/code]
happy 2006-11-17 10:11[code]%进退法
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=[t2,M_count,S_count];[/code]
happy 2006-11-17 10:11[code]%PRP直接计算算法(不精确一维搜索,k=83,mul_count=19842,sum_count=12722时间略长)
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'*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[/code]
happy 2006-11-17 10:11[code]%PRP重新开始计算算法(不精确一维搜索,k=133,mul_count=18421,sum_count=11850时间略长)
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'*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[/code]
happy 2006-11-17 10:11[code]%PRP直接计算降法(精确一维搜索,k=19,mul_count=8390,sum_count=3745,花费时间很少)
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'*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