进退法、0.618法、Powell法、罚函数法的程序设计,急！！

alei

RT,《最优化方法》的题目
涉及进退法、0.618法、Powell法、罚函数法的程序设计
可用Matlab或C程序语言编程，需包含这四类算法。
实例已给出如下，编一例类似程序。实在不行，可在此实例上适当修改。
有感兴趣的，完成后可发到本人邮箱：aleiho@163.com
见回复确认后，立即加分结贴。
p.s.:谢谢，急需。19日16时之前。

一、进退法、0.618法、Powell法、罚函数法的Matlab程序设计
罚函数法（通用）
function y=ff(x,k)
y=-17.86*0.42*x(1)/(0.8+0.42*x(1))*(1-exp(-2*(0.8+0.42*x(1))/3))*exp(-1.6)*x(2)-22.99*x(1)/(0.8+x(1))*(1-exp(-2*(0.8+x(1))/3))*x(3)+k*(x(2)-(1.22*10^2*(9517.8*exp(-1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)))/(1.22*10^2+9517.8*exp(-1.6-2*0.42*x(1)/3)*x(2)+19035.6*exp(-2*x(1)/3)*x(3)))^2+k*(x(3)-exp(-0.8-2*x(1)/3)*x(3)-exp(-2.4-2*0.42*x(1)/3)*x(2))^2;
% 主函数，参数包括未知数的个数n，惩罚因子q，惩罚因子增长系数k，初值x0，以及允许的误差r
function G=FHS(x0,q,k,n,r,h,a)  
l=1;
while (l)
    x=powell(x0,n,q,r(1),h,a);  %调用powell函数
        g(1)=ff1(x),g(2)=ff2(x) . . . g(p)=ffp(x);  %调用不等式约束函数，将其值
            %存入数组g
        h(1)=hh1(x),h(2)=hh2(x) . . . h(t)=hht(x);  %调用等式约束函数，将其值
%存入数组h
        for i=1:p
            if g(i)<r
                for j=1:t
                    if abs(h(j))<r(2)
                        continue;
                    else
                        break;
                    end
                end
            else
                break;
            end
        end
        if (i==p)&(j==t)   %如果所有约束条件的值都在允许的范围内就结束
%搜索，否则加大惩罚力度继续搜索
            l=0;
        else
            x0=x;
            q=k*q;
    end
end
   G=x
%powell算法，用于寻找无约束最优值点
 function  powel=powell(x0,n,q,r,h,a)  
    d=eye(n);         %n个线性无关的初始搜索方向
    k=1;
    kk=1;
    xx(1,1:n)=x0;
    while (kk)
        y(1,1:n)=xx(k,1:n);
        for j=1:n
            s(j)=HJ(y(j,1:n),d(j,1:n),q);       %调用0.618算法
            y(j+1,1:n)=y(j,1:n)+s(j).*d(j,1:n);
        end
        d(n+1,1:n)=y(n+1,1:n)-y(1,1:n);
        if (norm(d(n+1,1:n),2)<r)
            kk=0;
            break;
            else
                ww=0;
                m=1;
                for i=1:n
                    gg=ff(y(i,1:n),q)-ff(y(i+1,1:n),q);
                    if (gg>=ww)
                        m=i;
                    end
                end
                cha=ff(y(1,1:n),q)-2*ff(y(n+1,1:n),q)+ff(2*y(n+1,1:n)-y(1,1:n),q);
                cha1=2*(ff(y(m,1:n),q)-ff(y(m+1,1:n),q));
                if (cha<cha1)
                    s(n+1)=HJ(y(n+1,1:n),h,a,d(n+1,1:n),q)
                    xx(k+1,1:2)=y(n+1,1:n)+s(n+1).*d(n+1,1:n)
                    for j=m+1:n
                        d(j,1:n)=d(j+1,1:n);
                    end
                    k=k+1;
                    else
                        xx(k+1,1:n)=y(n+1,1:n);
                        k=k+1;
                  end
              end
          end
      powel=y(n+1,1:n)

function  w=HJ(x0,h,d,dd,q)   %0.618算法
 [a,b]=JTF(x0,h,d,dd,q);        %调用进退法算法，确定范围
    r=0.618;
    r1=a+(1-r)*(b-a);
    r2=a+r*(b-a);
    y1=ff(x0+r1.*dd,q);
    y2=ff(x0+r2.*dd,q);
    k=1;
    while (abs(r1-r2)>=0.1)
        if y1<y2
            b=r2;
            r2=r1;
            y2=y1;
            r1=a+(1-r)*(b-a);
            y1=ff(x0+r1.*dd,q);
        else
            a=r1;
            r1=r2;
            y1=y2;
            r2=a+r*(b-a);
            y2=ff(x0+r2.*dd,q);
        end
    end
    w=(r1+r2)/2
%进退法
function [a,b]=JTF(x0,h,d,dd,q)
 r0=0;
    y0=ff(x0+r0.*dd,q);
    k=0;
    l=1;
    while (l)
        r1=r0+h;
        y1=ff(x0+r1.*dd,q);
        if y1<y0
            h=d*h;
            r=r0;
            r0=r1;
            y0=y1;
        else
            if k==0;
                h=-h;
                r=r0;
            else
                l=0;
                break;
            end
        end
        k=k+1;
    end
        a=min(r,r1);
        b=max(r,r1);
二、进退法、0.618法的C语言程序设计
例1  设f(x)=x^2-2*x+4 ，试确定初始搜索区间，并输出迭代次数

main（）
{int n=0；
float t，h，a，b，c，d，f1，f2；
printf（“t，h：”）；
scanf（“%f，%f”，&t，&h）；
f1=t*t-2*t+4；
if（（t+h）*（t+h）-2*（t+h）+4<=f1）
{c=t+h;f2=c*c-2*c+4;}
else
{h=-h;c=t+h; f2=c*c-2*c+4;}
while(f2<=f1)
{h=2*h;d=t;t=c;f1=f2;c=t+h; f2=c*c-2*c+4;n=n+1;}
if(c<=d)
{a=c;b=d;}
else
{a=d;b=c;}
printf(“[a,b]=[%.2f,%f]\n”,a,b);
printf(“n=%d\n”,n);
}
运行情况如下：
（1）t,h:0,0.05
     [a,b]=[0.35,1.55]
     n=4
(2) t,h:-1.0,0.02
     [a,b]=[0.26,4.10]
     n=7
(3)t,h:-100,0.01
     [a,b]=[-59.05,63.83]
     n=13
例2  用0.618法求f(x)=8x^3-2*x^2-7*x+3 的局部最优解．允许误差ε=0.0001 ，初始点设为x0=0 ，初始步长h=1 ．

#include“math.h”
float f(x)
float x；
float y;
{y=((8.0*x-2.0)*x-7.0)*x+3.0;
return(y);
}
main()
{int n=0;
float t,h,a,b,c,d,x0,f0,f1,f2,fc,fd;
printf(“t,h:”);
scanf(“%f,%f”,&t,&h);
f1=f(t);
if (f(t+h)<=f1)c=t+h;
else{h=-h;c=t+h;}
f2=f(c);
while(f2<=f1)
{d=t;t=c;f1=f2;c=t+h;f2=f(c);}
if(c<=d){a=c;b=d;}
else{a=d;b=c;}
c=a+0.382*(b-a);fc=f(c);d=a+0.618*(b-a);fd=f(d);
do
{if (fc>=fd) {a=c;c=d;fc=fd;d=a+0.618*(b-a); fd=f(d);}
else {b=d;d=c;fd=fc;c=a+0.382*(b-a);fc=f(c);}
n++;
}
while (b-a>0.0001);
x0=(a+b)/2.0;f0=f(x0);
pintf(“n=%d,x0=%7.4f,f0=%7.4f\n”,n,x0,f0);
}
运行情况如下：
t,h:0,1.0
n=21,xo=0.6298,f0=-0.2034．

 

zhcf

呵呵~~好长...

powell

函数里面的参数的顺序是不不对啊，需要从新调整

ammon

不太会啊，等待牛人

ammon

确实感觉有些参数不对

zhoyi

有没有c++ 版？